一是体系整体的应力状态,通过在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%)。
联系客服