摘自：http://dxli75.blog.163.com/blog/static/106768289201142893419587/
lammps計算的應力有兩種：
一是體系整體的應力狀態，通過在thermo_style custom里加上pxx pyy pzz pxy pxz pyz字段可將給定時間步（由thermo N命令所指定）的體系應力值輸出，再求時間平均即可（實際上求出的是壓強張量，即負的應力值）； 該應力的計算用的是統計力學里的Virial定理（參見<<Computer simulation of liquid>> by Allen & Tildesley），所算出來的應力與宏觀應力是一致的（強調一下，用於平衡態）；
compute 1 all pressure thermo_temp virial
二是單個原子的應力，也就是樓主所討論的，通過compute stress/atom命令所計算出來的六個值；這個應力的計算則是通過對上述Virial定理所定義的體系應力按原子分解，即將公式中的按原子求和的算符拿掉（同時應將體系的體積換為單個原子的體積），不過，正如樓主所指出的，由於單個原子的體積計算太麻煩，lammps在計算時干脆去掉了體積項，這就是為什么用compute stress/atom命令所算出來的“應力”具有能量的單位的原因。
可以看出，在lammps里，如果要計算體系中某個區域（由region定義，可以是整個模擬盒）所圍成的“塊”的應力，只需將該區域里的所有原子的單原子應力值加起來，再除以這個區域的體積即可，無須進行單個原子體積的計算。
對原子應力進行Voronoi體積加權平均即可得到系統瞬時應力；系統瞬時應力的系綜平均值為宏觀測量的系統應力值。
compute s all stress/atom
compute p all reduce sum c_s[1] c_s[2] c_s[3]
variable Press equal -(c_p[1]+c_p[2]+c_p[3])/(3*vol)
thermo_style custom step temp etotal press v_Press
如果要求X方向的應力,就把所有原子用compute stress/atom得到的sigma(xx)加起來,去除掉"所有原子佔有的體積",就可以得到系統宏觀的sigma(xx).得到的值與別人發表的計算paper相差無幾(與實驗值相差大約10%)。