Lammps例子：弹性常数的计算（转摘）
转摘自知乎：https://zhuanlan.zhihu.com/p/359014805?utm_medium=social&utm_oi=869273158199443456。这个解释很清晰，帮助理解。
一、广义胡克定律
材料在线弹性范围内，应力 与应变
成正。由于一些剪切条件的等价性和对称性，对于对称性最差的三斜晶系，独立的弹性常数可减至21个。可用下面的矩阵表示：
由于对称性中只需要考虑一个半三角就行。
二、Voigt标记
三、应力应变法求解弹性常数
Lammps例子中的文件有 in.elastic.lmp、init.mod、potential.mod、displace.mod 以及对应的势文件。下面挑点有用的讲。主要是跟着 in.elastic.lmp 的流程理一下。
（1） init.mod
应变量up（最好测试一下微应变与最后得到的弹性常数的无关性，换几个值试一下）
单位转换不多说
定义盒子和填充原子可以用Lammps的建模方法，也可以用read_data，但是对于data文件里面是正交盒子的，需要在data文件里面添加 0 0 0 xy xz yz 的tilt倾角值，然后在read_data前面加个命令 box tilt large 不然在后面施加切应变时Lammps无法对你的正交盒子操作。
Lammps的定义盒子的方式挺好，后面可以讲讲。
（2） potential.mod
这个没啥好说的 （不知道反应力场行不行，没试过）
（3） in.elastic.lmp
主流程，用来不断调用 displace.mod
先让整体静水压为零，然后弛豫，记录平衡体系的压强和盒子尺寸；然后定义计算弹性常数的variable（公式里面是去掉残余压强项，压强/应变*压强和弹性常数的系数）。
开始施加应变模式，此处直接转到下面一小段，后面没啥好说的了。
（4）displace.mod
dir决定了想到得到的弹性常数的列或行（对称性，故无区别），根据施加的是正应变还是切应变，确定了对应的len0，在前面定义的弹性常数variable中需要，如下：
if "${dir} == 1" then &
"variable len0 equal ${lx0}"
if "${dir} == 2" then &
"variable len0 equal ${ly0}"
if "${dir} == 3" then &
"variable len0 equal ${lz0}"
if "${dir} == 4" then &
"variable len0 equal ${lz0}"
if "${dir} == 5" then &
"variable len0 equal ${lz0}"
if "${dir} == 6" then &
"variable len0 equal ${ly0}" 由于Lammps定义盒子是固定x方向（在我看来是这样的），然后用 xy、xz、yz 三个倾角来定义另外两个方向与x方向的关系，很直观。所以在施加切应变的时候，以dir=4为例，对应原长是lz0（x方向不动，z方向转动）。
重新读取平衡状态，定义potential，然后分别做正、负应变（正负重复实验减少误差）
最关键的是下面对施加的应变模式的理解，使用了change_box命令，改变盒子尺寸，然后remap原子坐标，然后弛豫获取压强计算弹性常数，思路就是这样。
对于deformation那段命令的理解：
variable delta equal -${up}*${len0}
variable deltaxy equal -${up}*xy
variable deltaxz equal -${up}*xz
variable deltayz equal -${up}*yz
if "${dir} == 1" then &
"change_box all x delta 0 ${delta} xy delta ${deltaxy} xz delta ${deltaxz} remap units box"
if "${dir} == 2" then &
"change_box all y delta 0 ${delta} yz delta ${deltayz} remap units box"
if "${dir} == 3" then &
"change_box all z delta 0 ${delta} remap units box"
if "${dir} == 4" then &
"change_box all yz delta ${delta} remap units box"
if "${dir} == 5" then &
"change_box all xz delta ${delta} remap units box"
if "${dir} == 6" then &
"change_box all xy delta ${delta} remap units box"所有的delta就是施加的应变模式，这里是把每次取一个应变，如x方向：
应力应变比值就是对应的弹性常数，然后正负各1次，求平均值，6个应变量各取1次，共12次计算就得到全部的弹性常数。当然知道晶体的对称性也可以稍微改下直接求需要的元素，不过没啥意义哈。
之前困扰我的是x和y方面施加应变的时候，为啥有个 deltaxy 这些。后来请教了人特别好的戴老师，思考后才发现是因为考虑到晶胞如果不是正交的，那正应变也会导致切应变，就是得把正应变分解下。
