An issue with rigid/nvt local stress calculation

LAMMPS LAMMPS Development
1

Dear LAMMP Developers,

We recently encountered an issue using rigid/nvt in local stress calculations. We did quite a lot tests on “rigid/nvt”. It seems to be a fundamental issue causing this problem.

In the example below, we used the TraPPE nitrogen model, a triatomic linear molecule with a ghost atom under NVT conditions. We used the 1D binning method (slices) and spherical binning method (shells) to calculate the stress tensor. The 1D binning method shows an anisotropic stress tensor where stress_xx is bigger than stress_yy and zz. However, for this homogenous system, stress components should be isotropic. In the spherical binning method, however, stress_xx, yy and zz are the same, as expected; however, they are larger than the results when using the 1D binning method. We are not sure what causes this issue. The input file is a sample script that we downloaded from the internet with further modifications. In the result file, the stress_xx, yy, zz are the first 3 columns. We have also tried with other N2 models as a diatomic molecule with no ghost atom, and the problem remains. Interestingly, when we use “SHAKE” for this rigid diatomic molecule, there is no such problem, and the local pressure is isotropic. However, we can not use “SHAKE” for a triatomic model of N2 with a ghost atom, which we need for our simulations.

If anyone knows how to fix this issue, please let us know. We would greatly appreciate it!
Thanks in advance!

Best regards,
Darren

Input file:

units		real	  
atom_style	full  
timestep	1.0	
boundary        p p p				   
read_data       Pure_Charged_N2.data
replicate	1 1 1
				 
variable        dt equal dt
variable        ext_temp equal 200
variable        Tdamp equal 100.0*${dt}		   

pair_style	lj/cut/coul/long 15.0 15.0          
pair_modify     tail yes			
bond_style      harmonic
angle_style     harmonic
kspace_style	pppm 1.0e-4                   

pair_coeff      * 1 0.0000 0.0000               
pair_coeff	2 2 0.07153934924399673 3.31000 

bond_style	harmonic   		
angle_style	harmonic		
dihedral_style	none				
improper_style	none			

bond_coeff	1 5000.00 0.550          
angle_coeff	1 500.0 180.00			

group ghost    type 1
group nitrogen type 2				

velocity nitrogen create ${ext_temp} 432567 dist uniform
reset_timestep  0                                

fix             2 all rigid/nvt molecule temp ${ext_temp} ${ext_temp} ${Tdamp}
variable        PotentialEnergy equal epair  
variable        Pressure equal press               
variable        Temperature equal temp        
variable        Density equal density           

fix             3 all ave/time 20 25000 500000 v_Temperature v_Density v_PotentialEnergy v_Pressure file ave.dens_22.0molL.out format %.8g
thermo          1000                           
thermo_style    custom step temp density epair press #Format for screen output of thermodynamics

###########################################Local Stress##################################
compute temp_thermo all temp
compute 33333 all chunk/atom bin/sphere 21.13061467 21.13061467 21.13061467 0.0 20.0 20 discard yes units box
compute stress_atom all stress/atom NULL
fix pressstrace_nitrogen all ave/chunk 20 25000 500000 33333 c_stress_atom[1] c_stress_atom[2] c_stress_atom[3] c_stress_atom[4] c_stress_atom[5] c_stress_atom[6] ave one file N2_Sphere.dat

compute 1D all chunk/atom bin/1d x lower 1 units box
fix pressstrace_bin1d all ave/chunk 20 25000 500000 1D _stress_atom[1] c_stress_atom[2] c_stress_atom[3] c_stress_atom[4] c_stress_atom[5] c_stress_atom[6] v_press_atom ave one file N2_1D.dat

#########################################################################################

run             1000000
write_data N2_Test.data nocoeff

Local stress spherical binning results:
c_stress_atom[1] c_stress_atom[2] c_stress_atom[3] c_stress_atom[4] c_stress_atom[5] c_stress_atom[6]
1000000 20 1333.2040000000002
-17633.2 -18276.1 -13162.8 1914.41 4169.87 3701.46
-20358.6 -21695.8 -21426.7 -3250.28 -4408.58 -4390.1
-12776.5 -13240.1 -15902.5 5872.22 5410.6 5998.98
-26185.3 -24278.9 -22579.4 -4545.38 -2727.44 -3235.05
-16111.8 -20611.7 -20110.8 -1253.68 244.896 -255.442
-20036.3 -18842.1 -19440.7 717.253 -75.9576 -67.3112
-22295.3 -21308.3 -20166 -1993.89 -706.762 -756.187
16692.4 -15995.7 -18999.1 2211.18 187.112 410.856
-17100.8 -18113.3 -18436.3 1200.35 842.616 711.992
-21065 -20513 -20670.5 -1929.48 -1242.32 -1876.48
-17290.4 -19276.9 -17602.6 19.0735 1112.23 1636.44
-18565.7 -18486.1 -19329.4 902.696 200.472 -244.609
-20114.7 -19827.7 -19773.9 -713.509 -967.597 -1029.73
-18362.9 -17209.8 -17497.9 1498 1737.4 1783.97
-18416.5 -19562.8 -20807.5 -887.017 -1526.39 -1572.01
-19985.7 -18304.9 -19300.1 821.711 -534.355 -646.35
-17998.2 -18854.5 -17454.5 -698.163 1406.9 1546.08
-19965.3 -19433.6 -19011.2 -289.413 143.979 502.981
18218.3 -18427.9 -18660.6 503.664 -249.583 -808.416

Local stress 1D binning results:
c_stress_atom[1] c_stress_atom[2] c_stress_atom[3] c_stress_atom[4] c_stress_atom[5] c_stress_atom[6]
1000000 43 3000.0000000000005
-18975.3 -16738 -18911.4 12.1673 -134.45 -421.765
-19045.6 -18322.6 -14983.8 -123.347 296.424 1200.33
-18801.4 -16605.9 -17463.9 81.6195 -123.494 -997.824
-19520.3 -14830.1 -17236.5 236.575 -302.237 -943.179
-18821.8 -18633.4 -15081.4 -462.574 451.501 1762.79
-17738.4 -16012.6 -17861.3 263.669 -215.727 -1167.35
-19726.7 -17710.7 -15413.5 -416.366 -169.879 1007.67
18674.9 -16540.1 -15607.6 -187.901 633.202 486.198
-19426 -14906.5 -17826.3 1029.56 -413.444 -679.495
-18912.7 -17923.1 -15637 -719.523 394.353 1065.08
-18151.7 -17134.7 -16490.1 -140.712 395.706 -224.941
-18968.6 -16505 -17728.8 288.631 -554.019 -840.03
-19603.5 -16444.3 -17547.4 186.959 -32.5546 -335.325
-18949.2 -15753.8 -16553 140.942 -745.093 -1369.78
-19274 -16649.9 -15052.9 -224.317 1018.45 1321.07
-18746.5 -16918.4 -16467.5 245.941 125.261 989.201
-19046.2 -15999.2 -16344.1 178.94 124.452 -103.387
-19536.9 -17228.4 -17802.1 352.23 -520.72 -1154.8
-18268.4 -17114.4 -16686.8 -470.838 439.22 484.356
-19664.4 -16210.4 -18298.5 141.125 -1506.72 -1453.68
-19978.6 -16234.8 -14021.1 506.765 1796.98 1972.85
-19160.2 -18595 -17421.8 -2029.89 -1014.55 -1054.93
-18074.2 -15049 -16938.7 2137.79 202.742 108.748
-17080.8 -17681.5 -15323.1 -1388.55 402.988 1535.94
-20088.6 -16461 -18204.1 483.316 -851.119 -2389.05
-17767.7 -16020.5 -16134.1 519.562 281.312 775.682
-21959.3 -18384.4 -16106.5 -1893.78 761.052 827.846
-16514.4 -15142.6 -17301.4 2083.16 44.2289 -106.988
-19248 -18290.1 -18290.3 -1959.32 -1683.53 -1078.77
-16986.8 -17516.7 -16684.8 -205.237 998.071 717.658
-19837.8 -15479.8 -15805.8 359.679 411.114 147.705
-19045.4 -14976.5 -15407.5 2441.7 800.705 887.143
-18323.8 -17097 -17742.3 -1734.05 -1424.16 -859.861
-19500.9 -16639.3 -16352.3 -578.791 1020.1 449.424
-18437.8 -14607.5 -15447.2 2339.29 445.652 -187.664
-23260 -18199.1 -15335.8 -1286.76 2217.86 1274.32
-16663.5 -17110.7 -18100.1 -1086.55 -2198.04 -1409.4
-19336.4 -15179.2 -16001.5 2863.89 336.3 1298.82
-21619.3 -16821.5 -17333.7 -2249.76 -1729.88 -1578.86

2

I have bumped into a problem in the past when doing a simulation using fix rigid/nvt and fix rigid/npt commands and trying to extract properties related to stress. The conclusion had been that if you use these commands without turning off the interactions inside the rigid body, the computation of stress is bogus.

To see if this is the case for you, you can try including the following lines in your input script:

delete_bonds	N2 multi
neigh_modify	exclude molecule/intra N2

PS: I am assuming a group-ID N2 that encapsulates all your N2 molecules and that you have set proper molecule-IDs for the atoms composing each of these N2 molecules in the LAMMPS data file.

3

Thanks for sharing this. I did some trials. It seems the results from the spherical binning method are still larger than the 1D binning method. Also, I noticed, if I delete bonds, sometimes it will lead to “Domain too large for neighbor bins” or “bonds missing” errors. Just wondering if you have more suggestions on these errors? Thanks very much!

4

A general remark: without knowing what your LAMMPS version is and what platform you are running on and which compilers and settings you use, it is hard to give any specific advice.
It can make a big difference if you are using some accelerator package (you should always confirm that you get the same results with and without). And some compilers have a non-zero risk to miscompile LAMMPS when using high optimization levels.

This typically only happens, if you have a very short cutoff. LAMMPS typically uses half the largest pairwise cutoff. But with short cutoffs this will overflow the list of bins. You can set the bin size to a larger value using the neigh_modify command.

I am confused. How can you have “bonds missing” errors, when you delete bonds?

While it is a good idea to post questions about possible issues in LAMMPS here on MatSci.org, so that experienced LAMMPS users can give their comments (they know the science often better than the developers), the ultimate location for a confirmed issue is the LAMMPS issue tracker on github: GitHub · Where software is built