[Python知识库]LAMMPS学习总结4

1、聚合物驰豫平衡的技巧

如果系统没有充分压缩和退火，有时会意外地产生人造气泡或空隙。因此，在初始阶段，应用具有较小ε值的高温和高压以确保聚合物良好缠结。

2、分子动力学时间尺度理解

MD时间步大部分时候设定在1fs这个尺度上，这是因为原子之间作用力的量级决定的，大了的话结构会无法收敛到势能面上，使得计算无法继续。（https://zhuanlan.zhihu.com/p/260828094）

3、minimize或run命令就会出现错误：Atom sorting has bin size = 0.0原因

1、力场没设置或者有问题
2、没有给cut off，可能是mini写在了力场设置之前

4、MS的car文件转出的data自带的力场参数对应lammps什么单位？

real

5、已经转换为标准格式的XYZ文件导入ISAACS中报错，提示格式错误，怎么回事？

答：标准XYZ文件最后空行只有一行，python处理dump转为标准xyz容易最后多空出来几行，这个引起的，手动删掉即可。正确如下

6、fix deform拉伸命令理解

fix 5 all deform 1000x erate 0.002 units box #设置沿x方向拉伸，应变率为0.002（1/ps）
fix ID group-ID deform N parameter args … keyword value …
N是每N步实施一次变形，如果这个数越大，那么每N步变形就越大，如果这个数较小，比如说是1，那么每跑一步就变形一次，这样模拟的更细，精度更高。总变形是由后面的参数确定的，N影响不了总变形。
fix 1 all deform 100 x erate 0.0001 units box remap x
run 40000
deform后面的N为1或100的区别在于，N为1时，是每运行1步，就变形1次；而N为100时，就是第一步变形了，第2-100步没有变形，第101步变形，第102-200步没有变形。
0.0001：应变速率为0.0001.
总应变是由erate=0.0001/ps、timestep=0.001ps、run 300000步三者决定，总应变为0.03，与deform后的N值无关！
upto关键字表示从当前时间步长开始直到指定的时间步长执行运行。例如。如果当前时间步长为10,000并且使用“运行100000 upto”，则将运行额外的90,000个时间步长。
（1）N=1时，意味着每运行一步，就变形一次，总共运行300000步要达到变形量0.03，也就是说每次运行的应变为10e-7（0.03/300000），即每运行一步，应变增加10e-7，运行到10000步时，总应变为10e-3;则在step=1~10000的过程中，总应变依次为10e-7，2 * 10e-7,3*10e-7…。
(2)N =10000时，意味着每运行10000步，就变形一次，总共运行300000步要达到变形量0.03，亦即在每个10000th变形一次，且应变为10^-3;;则在step=1~10000的过程中，总应变依次为0,0 0,…，0,10^{-3，step达到10000时，应变一下子就由0增加大10}-3。

上述分析说明，实际上N=1和N=10000时，模型的实际变形速率确实是不一致的，N=1是每一步变形，变形速率为10^-4/ps; N=10000时在10000步之前都保持0应变，而在每个10000 th时，直接变形达到规定值,变形速率为1/ps。N=10000时的实际变形速率是N=1的10000倍。

剪切应变是无单位的，定义为偏移量/长度，其中长度是垂直于剪切方向的盒子长度（例如，用于xy变形的y盒子长度），偏移量是剪切方向上的位移距离（例如，用于xy变形的x方向） 从不受约束的方向。

7、报错ERROR: Thermo keyword in variable requires thermo to use/init press (…/thermo.cpp:1399) Last command: run 25000

错误出现代码区域：

#计算应变值
variable tmp equal "lx"
variable L0 equal ${tmp}
variable strain equal "(lx - v_L0)/v_L0"
# 将拉伸强度改为Gpa
variable stressx equal "-pxx/10000"
variable stressy equal "-pyy/10000"
variable stressz equal "-pzz/10000"
#deform拉伸 x方向
fix 8 all npt temp 800 300 0.2 y 1 1 2 z 1 1 2
dump 6 all custom 1000 lashen.xyz id type x y z 
fix 9 all deform 1 x erate 0.005 units box remap x
thermo_style  custom step temp v_strain v_stressx v_stressy v_stressz
thermo          500
run 25000

出现原因及解决：

通过查询源码可以看出，问题出现在压力的变量设置上，这里有个判断，当pressure为0时才报这个错误，由于我是X方向拉伸，将variable stressx equal "-pxx/10000"命令去掉，换成pxx输出就没问题了，但是我看了一下pxx的输出值，从一开始就不是0啊，这我也不知道为啥了。

8、这个样子切片，compute计算的是每个原子的能量，那我的结果中最后一列是每个区域中原子的能量还是区域中所有原子的能量？

回答：Lammps手册中fix ave/chunk默认是norm，norm是这么说的，如下。输出是所有N个重复样本的原子之和，除以原子数。因此，最后一列代表的动能计算：每个原子的动能相加/原子数N。

9、ERROR: Expected floating point parameter in input script or data file (…/fix_nh.cpp:144)

答：fix 4 all npt temp 300 300 0.2 iso x 1000.0 1000.0 2多了个iso

10、MSD获得扩散系数的正确套路

分子动力学的一个重要的优势是能获得分子的运动行为，其中一个我们关心的物理量是分子或者原子的扩散行为。然后用LAMMPS给的计算命令(compute msd)得到的结果很难拟合出比较合理的扩散系数。有如下几个原因导致失败。

1、 compute msd命令输出的msd曲线平均次数太少，导致曲线光滑程度不够，计算的数据本身统计不够，数据可靠性比较低。
2、 拟合时间区间不合理：选取比较合理的拟合区间非常重要。我们计算自扩散系数使用的公式为MSD(t)=2dD*t.(d为模拟的维度)。这个公式表示 MSD(t)是t的1次方。而不是其他次方。一般我们算出MSD曲线后，要画log-log曲线来观察。你会发现，这个曲线并不一定都是 t的1次方。

一般来说，我们计算出来的 MSD(t)和t的关系是如下关系：

当alpha<1的时候，被称为sub-diffusive行为，这时候，粒子受到系统中其他粒子的约束，而且不容易逃出约束的势井。当1<alpha<2的时候，被称为 super-diffusive行为，一般体系里有能量的输入,导致粒子运动被加强。当alpha=2时，被称为弹道区，粒子行为类似于匀速运动。

为了确认alpha有没有出现1的区域，我们一般画log-log曲线来观察。如上图，然后确认 1 次方区的时间范围，取对应的时间区间进行拟合。而不是拿全线段来进行简单拟合。