[lammps-users] verifying the pressure tensor from the trajectory

Dear LAMMPS users and developers,

I am trying to understand how the pressure tensor from “compute pressure” (https://lammps.sandia.gov/doc/compute_pressure.html) is calculated under PBC. According to this page, “the position and force vectors of ghost atoms are thus included in the summation”. It’s not clear whether the position and force vectors of original atoms are also included (force vectors of original and ghost atoms should be the same?).

In theory, the trajectory can output position (x, y, z), unwrapped position (xu,yu,zu), image flags(ix,iy,iz), force (fx,fy,fz). So whether we are only including ghost atoms (we get that information from image flags) or both original and ghost atoms, we should be able to calculate the pressure tensors that should be the same to pxx, pyy, pzz, pxy, pxz, pyz from the thermo_style, right?

I am using the following toy problem for this verification.

units real
variable T equal 86.4956
variable V equal vol
variable dt equal 4.0
variable p equal 400 # correlation length
variable s equal 5 # sample interval
variable d equal $p*$s # dump interval

convert from LAMMPS real units to SI

variable kB equal 1.3806504e-23 # [J/K] Boltzmann
variable atm2Pa equal 101325.0
variable A2m equal 1.0e-10
variable fs2s equal 1.0e-15
variable convert equal {atm2Pa}*{atm2Pa}{fs2s}*{A2m}{A2m}*{A2m}

setup problem

dimension 3
boundary p p p
lattice fcc 5.376 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
region box block 0 4 0 4 0 4
create_box 1 box
create_atoms 1 box
mass 1 39.948
pair_style lj/cut 13.0
pair_coeff * * 0.2381 3.405
timestep ${dt}
thermo $d

equilibration and thermalization

velocity all create $T 102486 mom yes rot yes dist gaussian
fix NVT all nvt temp $T $T 10 drag 0.2
run 8000

reset_timestep 0
compute 1 all temp

thermo_style custom step temp pxx pyy pzz pxy pxz pyz
dump dumpTrj all custom 10 traj.trj id x y z ix iy iz xu yu zu vx vy vz fx fy fz
dump_modify dumpTrj sort id

run 0

please note that unwrapped coordinates != ghost atoms.
the way how LAMMPS computes the virial and from it the pressure is outlined in the paper cited on the compute pressure documentation.
I suggest you study that paper. retracing the pressure computation after the fact is difficult without adding debug code to the source code and recompile.


Dear Axel,

Thanks for pointing to the paper. The equation on the webpage refers to the atom form in the paper. It makes sense to include all periodic images. In addition, since the forces of the ghost atoms are not saved in the trajectory, we can’t use the trajectory only to recalculate the pressure tensor.

Thanks, Wei