Hello,
I am trying to understand how the virial is calculated for the tip4p water model in lammps, and was wondering if anyone could help me understand. Specifically, if I look in pair_lj_tip4p_long.cpp (line 339 for me) I see the following
if(vflag) {
domain->closest_image(x[i],x[iH1],xH1);
domain->closest_image(x[i],x[iH2],xH2);
v[0] = x[i][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
v[1] = x[i][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
v[2] = x[i][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
v[3] = x[i][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1];
v[4] = x[i][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2];
v[5] = x[i][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2];
}
Now I understand that the domain call returns the distances of the hydrogens from the oxygen wrapped using the minimum image convention, but why are we multiplying f0, the force added to the oxygen atom, by the position of the oxygen atom, x[i] here? Obviously, this would make the virial calculation depend on the frame of reference, which I don’t believe is correct. I assume I am missing something here. Could anyone point me to where I am incorrect please? I can’t seem to find anywhere in the .cpp file where x[i] is modified.
If the coords of N atoms (3 in this case) are close together
then the virial is V = sum (1 to N) Fi dot Ri.
That is translationally invariant if the Fi sum to 0. So it
doesn’t depend on the frame of reference (which image).
Steve
Hello,
I am trying to understand how the virial is calculated for the tip4p water
model in lammps, and was wondering if anyone could help me understand.
Specifically, if I look in pair_lj_tip4p_long.cpp (line 339 for me) I see
the following
if(vflag) {
domain->closest_image(x[i],x[iH1],xH1);
domain->closest_image(x[i],x[iH2],xH2);
v[0] = x[i][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
v[1] = x[i][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
v[2] = x[i][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
v[3] = x[i][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1];
v[4] = x[i][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2];
v[5] = x[i][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2];
}
Now I understand that the domain call returns the distances of the
hydrogens from the oxygen wrapped using the minimum image convention, but
why are we multiplying f0, the force added to the oxygen atom, by the
position of the oxygen atom, x[i] here? Obviously, this would make the
virial calculation depend on the frame of reference, which I don't believe
is correct. I assume I am missing something here. Could anyone point me to
where I am incorrect please? I can't seem to find anywhere in the .cpp file
where x[i] is modified.
please note that the code you are quoting is only adding the contributions
to the virial for the coulomb interactions of the point M, and only for one
of the two involved atoms. since the forces are projected back from M to
the sites of the oxygen and hydrogen atoms, it has to be treated like this.
axel.