The pressure calculated by compute pressure command and compute stress/atom should be the same. However, I found that the following result is not the same:
From compute pressure:
variable temp atom xfx+yfy+zfz #eV/A * a =eV
compute c9 all reduce sum v_temp
variable part2 equal c_c9/3/vol160.21766208 # eV/A3 > Gpa
From compute stress/atom:
compute c1 all stress/atom NULL virial
compute c11 all reduce sum c_c1[1]
compute c12 all reduce sum c_c1[2]
compute c13 all reduce sum c_c1[3] # bars
variable s10 equal -(c_c11+c_c12+c_c13)/3/vol/10000 # bars > Gpa
v_part2 ≠ v_s10.
The system is non-periodic system. So N’= N = total number of atoms in the equation in compute pressure command page. The unit is metal.
Could you explain why they are not equal to each other please? Or am I doing something wrong ? Thank you a lot in advance.
Do you get the same answer when you let compute
pressure do the calculation? Instead of your formula?
And let LAMMPS take care of the units instead
or your conversion factors?
I.e. show some thermo output results from a simple-as-possible
script (run 0) where the 2 LAMMPS computes calculate the
pressure.
Thank you for your reply. When I used compute pressure, the result is consistent with compute stress/atom. I want to check if compute pressure command really used the formula provided on the page. The pressure from kinetic energy of both formula NKbT/V and ( compute pressure ke command) are the same. The virial pressure from formula Σ r*f/3V and (compute pressure viral command) is not the same. The thermo output is below (unit metal, non-periodic system).
variable temp1 atom xfx # eV/A * A =eV, unit=eV
compute c1 all reduce sum v_temp1 # eV
variable formula equal c_c1/vol1602176.6208 # bar, eV/A3 --> bar , 1eV/A3=1602176.6208 bar
compute pressure_virial all pressure NULL virial
variable pressure_command equal c_pressure_virial[1] # bar
There must be something extra going on in your simulation. The large changes in v_formula suggest that you have not turned off the periodic boundary conditions. I append my script below.
Aidan
[athomps@…5834… threebody]$ more in.test
testing virial formula for non-periodic system
units metal
atom_style atomic
boundary s s s
lattice diamond 5.431
region myreg block 0 4 0 4 0 4
create_box 1 myreg
create_atoms 1 region myreg
mass * 28.06
velocity all create 1800 5287287
pair_style sw
pair_coeff * * Si.sw Si
fix 1 all nvt temp 1800 1800 0.1
variable temp1 atom xfx
compute c1 all reduce sum v_temp1
variable formula equal c_c1/vol1602176.6208
compute pressure_virial all pressure NULL virial
variable pressure_command equal c_pressure_virial[1]
thermo_style custom step temp press v_pressure_command v_formula
Langevin is a good example of the ‘something extra’ I was referring to. It applies extra contributions to the force on each atom. LAMMPS is clever about not including these extra force contributions in the per-atom virial calculation. Your v_formula is not. It is hard for me to understand why you would go to the trouble of checking the LAMMPS virial calculation for a non-periodic cluster that is also coupled to a Langevin heat bath and then assume that LAMMPS is somehow wrong. If you are going to test the rigor of LAMMPS, you need to do it rigorously. Otherwise, it is best to just assume the doc page is telling you the truth.