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.