What I would like to understand is the difference between compute 1 all pressure NULL pair, and compute pair/local fx*dx fy*dy fz*dz (where fx*dx etc. are simple additional pair/local computations I added to my local LAMMPS installation).
Generally, I would expect that using e.g.
compute mv all pair/local fx*dx fy*dy fz*dz
compute vas all reduce sum c_mv[1] c_mv[2] c_mv[3]
variable myPress2 equal (c_vas[1]+c_vas[2]+c_vas[3])/(3*vol)
would result in the same result as compute myPress all pressure NULL pair.
For an exemplaric system, such as the DPD one I recently posted here, this is not the case. The values are close, but differ by some ca. 1.5%. According to the documentation of compute pressure, the dof “necessarily includes atoms from neighboring subdomains (so-called ghost atoms) and the position and force vectors of ghost atoms are thus included in the summation”. In the documentation of the pair/local, no such statement is made. This could be a valid reason for the discrepancy, yet must not necessarily be the one – may I kindly as for some insights into how the pair/local compute handles ghost atoms, or what other reasons there could be (e.g., I misunderstand the computation of the pressure) for the discrepancy?
Well, you can read what the source code says also read the Thompson paper listed in the documentation of compute pressure which goes into the theoretical details of how this calculation is done with PBCs.
My understanding is that the (main) reason for the difference is that you are using pair style dpd.
When using compute pair/local all force contributions are recomputed using the Pair::single() function. With DPD this is a problem, though for two reasons: 1) its force includes a velocity dependent component, but when the virial is tallied from the regular force computation, it is done before the second velocity update of the velocity verlet algorithm while the pair/local computation is run after. 2) its force contains a random noise component and those random numbers are not the same.
The statement in the compute pressure docs has to be seen more in reference to using the “shortcut” computation of F dot r for the (global) pressure. This cannot be done from the position and force information available to atom style variables since you would need to use the position and force contribution to ghost atoms separately, but that information is not available since there is a reverse communication and the forces on ghosts atoms are added to their original local atoms. In LAMMPS this global virial computation is done before the reverse communication at the end of the pair style compute function in case the pair style supports it (not all do).
If you want to see, if your method is sound, you have to use a different, non-DPD pair style.
Indeed, the force being recalculated seems to be what I was missing. With the input below, for example, I can successfully reproduce myPress2 == c_pressVirial == c_pressPair (myPress above). Thank you very much!
# LAMMPS input file: test the pressure computation
units lj
atom_style dpd
variable kb equal 1.0
comm_modify vel yes
neigh_modify every 1 delay 0
variable seed equal 2570894
variable x0 equal 10
variable y0 equal 10
variable z0 equal 10
variable Natoms equal 1000
region thebox block 0 ${x0} 0 ${y0} 0 ${z0}
create_box 2 thebox
create_atoms 1 random ${Natoms} ${seed} thebox
mass * 1
variable T equal 1.
variable cutoff equal 1.0
pair_style lj/cut 1.12246
pair_coeff * * 1.0 1.0 1.12246 # type i, type j, eps, sigma, LJ cutoff
run_style verlet
fix 1 all langevin ${T} ${T} 0.5 ${seed}
fix 2 all nve
minimize 1.0e-8 1.0e-10 10000 100000
thermo_style custom step temp vol dt press pe ke density atoms
thermo 10
timestep 0.005
run 1000
# output we are interested in
compute mv all pair/local virx viry virz
compute vas all reduce sum c_mv[1] c_mv[2] c_mv[3]
variable myPress2 equal (c_vas[1]+c_vas[2]+c_vas[3])/(3*vol)
compute va all reduce ave c_mv[1] c_mv[2] c_mv[3]
variable myPress3 equal (c_va[1]+c_va[2]+c_va[3])/(3)
compute pressPair all pressure NULL pair
compute pressVirial all pressure NULL virial
compute press2 all pressure thermo_temp
thermo_style custom step temp vol dt press v_myPress2 v_myPress3 c_pressPair c_pressVirial c_press2
thermo 1
# actual simulation
# adjust the time-step
timestep 0.01
# run the simulation
run 20000
also read the Thompson paper listed in the documentation of