stress/atom compute inconsistent with pressure in xz direction

Lammps developers,

I’ve started using the stress/atom compute to do some analysis of my simulations and have noticed a strange inconsistency in the values I get for the summed stress per atom/volume. All of the values in the stress tensor are exactly the same as the pressure except for the 5th value (stress in the xz direction) which can be several orders of magnitude larger than the pressure in that direction. Am I missing something in my input script or is lammps calculating this incorrectly?

I’ve modified (and attached) the input script for the CHO example in the reax examples directory to quickly reproduce the error.

(I’m using the latest version of lammps I compiled moments ago via SVN)

Ben Jensen
PhD Candidate
Department of Mechanical Engineering - Engineering Mechanics
Michigan Technological University
Houghton, MI

Added to the input script:
compute spa all stress/atom
compute sta all reduce sum c_spa[1] c_spa[2] c_spa[3] c_spa[4] c_spa[5] c_spa[6]
variable staxx equal -c_sta[1]/vol
variable stayy equal -c_sta[2]/vol
variable stazz equal -c_sta[3]/vol
variable staxy equal -c_sta[4]/vol
variable staxz equal -c_sta[5]/vol
variable stayz equal -c_sta[6]/vol

thermo_style custom step v_staxx pxx v_stayy pyy v_stazz pzz v_staxy pxy v_staxz pxz v_stayz pyz
thermo 1
thermo_modify flush yes norm yes

Log file:
Step staxx Pxx stayy Pyy stazz Pzz staxy Pxy staxz Pxz stayz Pyz
0 -93.792302 -93.792302 -163.23365 -163.23365 -62.793941 -62.793941 -156.96152 -156.96152 -3499.2057 -17.167756 23.825025 23.825025
1 -91.595083 -91.595083 -160.92794 -160.92794 -61.241049 -61.241049 -156.35643 -156.35643 -3499.1619 -17.264862 23.656523 23.656523
2 -54.096651 -54.096651 -120.39142 -120.39142 -28.420953 -28.420953 -147.34606 -147.34606 -3493.308 -17.140318 20.845109 20.845109
3 -1.1194428 -1.1194428 -63.185179 -63.185179 17.904336 17.904336 -134.70012 -134.70012 -3484.9421 -16.966068 16.860407 16.860407
4 63.298664 63.298664 6.2105432 6.2105432 73.905638 73.905638 -119.3871 -119.3871 -3474.7611 -16.836703 11.983748 11.983748
5 136.62205 136.62205 84.978475 84.978475 137.06027 137.06027 -102.0017 -102.0017 -3463.1947 -16.829739 6.3801103 6.3801103
6 216.79854 216.79854 170.81644 170.81644 205.24412 205.24412 -83.008251 -83.008251 -3450.5883 -17.030073 0.17572053 0.17572053
7 301.99324 301.99324 261.69454 261.69454 276.55134 276.55134 -62.77614 -62.77614 -3437.2267 -17.505977 -6.5198666 -6.5198666
8 390.4377 390.4377 355.66485 355.66485 349.13163 349.13163 -41.65545 -41.65545 -3423.3728 -18.328384 -13.600946 -13.600946
9 480.40706 480.40706 450.88004 450.88004 421.22083 421.22083 -19.932255 -19.932255 -3409.2509 -19.545214 -20.963824 -20.963824
10 570.1618 570.1618 545.50001 545.50001 491.07148 491.07148 2.1140262 2.1140262 -3395.0714 -21.200833 -28.506049 -28.506049
Loop time of 0.256282 on 1 procs for 10 steps with 105 atoms

in.CHO (1.01 KB)

Looks like a bookkeeping error, I'll ask Ray to look into it.

Steve

Change this line in src/pair.cpp in the v_tally() method:
from:
v[3] = x[i][0]*fi[1];
v[4] = x[i][0]*fi[1];
to:
v[3] = x[i][0]*fi[1];
v[4] = x[i][0]*fi[2];

Will be in the next patch ...

Thanks,
Steve

Steve,

Thanks! That fixed it for me.

Ben Jensen
PhD Candidate
Department of Mechanical Engineering - Engineering Mechanics
Michigan Technological University
Houghton, MI