# stress calculation

Hi Steve,

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/vol
160.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.

Regards,
Wenting

Do you get the same answer when you let compute
And let LAMMPS take care of the units instead

I.e. show some thermo output results from a simple-as-possible
script (run 0) where the 2 LAMMPS computes calculate the
pressure.

Steve

Dear Steve,

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/vol
1602176.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

Step Temp Press v_pressure_command v_formula
0 100 660.67522 -28.135347 -1295.616
1 99.863336 661.11167 -26.783308 -433.81436
2 99.347107 661.67228 -22.742482 28280.698

Regards,
Wenting

I cut and paste your variable definitions into a script for Stillinger-Weber silicon, . I got the result you expected, not the result you got.

Step Temp Press v_pressure_command v_formula
0 1800 14991.815 -28.118473 -28.118475
10 1112.0424 18653.196 9484.1417 9484.1425
20 751.11539 14770.576 8505.5515 8505.5522

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/vol
1602176.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

thermo 10
run 100

Hi Aidan,

Thank you. I have turned off the periodic boundary condition. If you change “ fix 1 all nvt temp 1800 1800 0.1” to
fix 1 all langevin 1800 1800 1 7777

fix 2 all nve

The result is not the same. The thermo output is below:

Step Temp Press v_pressure_command v_formula
0 1800 14991.815 -28.118473 -14192.49
10 1103.3537 18527.535 9465.2731 57755.593
20 753.2934 14825.636 8381.3926 14700.003

Regards,
Wenting

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.

Aidan