Does LAMMPS calculate pressure correctly?

LAMMPS LAMMPS Mailing List Mirror
1

Dear all

I output all atom’s unscaled atom coordinates, velocities and forces on atoms (‘x,y,z,vx,vy,vz,fx,fy,fz’ in dump command), then calculate the temperature and pressure through MATLAB. The simulation boundary is ‘pp pp pp’ and my system is polymer system. I use formula “xfx+yfy+z*fz” to calculate virial. When calculate pressure of the system, we should make sure all atom in the same box.

The pressure IS NOT equal the result calculate by LAMMPS. My result is 1.1364 and the LAMMPS gives me 1.252805. After checking the source codes, the problem may cause by the way of LAMMPS calculating virial. The LAMMPS calculates virial when calculate force. If one atom is out of boundary and one is in the box, the way LAMMPS use will not give the correct virial! Situation will be worse if system is small.

If the pressure isnot accuracy, other command using pressure will give the correct answer?

Best wishes!

2

I think the way you calculate virial is wrong. At least for pairwise force, it only depend on relative distance in the book of “computer simulation of simple liquid” formula 2.57.

Yi

3

The reason you get the result different from LAMMPS is that LAMMPS does not include forces from the “ghost” atoms into the virial but does include them in the force output.

The LAMMPS’ result is correct while yours is not. To understand this, consider the simple example of a perfect crystal: in the dump output, all forces would be zero, and so the virial pressure calculated by the naive method would also be zero, whichever lattice parameter you take. Obviously, this should not be the case. And it is not, because actually you should not include the forces acting through PBC when calculating virial pressure.

Vasily