起源

写这篇博客主要是总结自己的这周和上周的一些经验，这两周被一个小程序折腾了好久。

其实这就是一个很简单的分子动力学(molecular dynamics,MD)的计算，但是程序始终在爆炸，能量等不收敛。

首先必须要说明，在科研中我们大量会使用 Charmm, Lammps,这些程序不需要自己编写程序，而且MD程序需要大量复杂的调优技术，比如如何安排计算减少计算量。

这个博文我假设读者知道全部的 MD 知识，至少上过基本的课程或者文献，这里不讨论那些在大部分课本上能学习到的细节。

重要细节

PBC

周期性边界是整个 MD 模拟的核心，重要的一点就是用小部分的模拟大的体系，提高计算效率。

其实核心就是这个文献讨论的这幅图。

在两个原子之间距离大于截断距离，一般是体系尺寸的一半，需要进行调节。因为需要把两个原子之间的作用距离调节为原子和另一个原子的镜像原子之间的作用。

这里的一个问题是，如果单单计算能量，其实这个距离的方向没有关系，因为在 LJ 能量计算中都是偶次项，向量的方向没有关系。

简单来说，对于 $R_{ij}$, 因为位移在物理学中是矢量，所以必须非常仔细如何求解。

同时 $R_{ij}$ 到底是 i - j 还是 j - i？

是 i -j ,因为这里还要涉及到矢量在数学中是如何定义的，依照正确的数学矢量定义，必须是 i - j 才能得到正确方向的 $R_{ij}$。

加速计算

其实这些都是些小细节，但是细节是魔鬼还是需要好好考虑:

(1) 计算能量和力的时候，虽然是一个双重循环，但是如果 i 从 0 到 N-1，那么内层循环 j 就只需要从 i + 1 到 N - 1即可，因为能量是共同分享的，力可以用牛顿第二定律求出来。

(2) 计算力的时候不需要带入系数，最好的，比如计算动能的时候，不需要每一个原子动能都乘以 0.5，可以在速度平方加和后直接在总和上乘以 0.5 即可，这样可以省下 N -1 次乘法操作。如果你的体系是溶液中的蛋白质，10万多的原子体系，这样的小细节还是能省下很多计算资源的。

(3) 周期性边界条件有两个影响，一个是计算距离需要考虑镜像原子，二是在原子超出体系边界的时候必须把原子拉回体系内。

约化单位

约化单位是必须在所有计算中保持的，只有在输入输出提供可以给人参考的数据时才需要考虑约化单位和真实单位的转化，否则都统一采用约化单位。

约化单位其实很简单，核心就是保持一个约化系数，这个约化系数在所有的真实单位中都统一转向约化单位。

在晶体结构的体系还有一个要点需要注意，就是晶体的原子有一个原胞坐标系，比如 FCC 是(0,0,0),(0,0.5,0.5),(0.5,0,0.5),(0.5,0.5,0)。但是这并不是 MD 希望的实验坐标系，所以还需要把原胞坐标系转化到实验坐标系。这儿只需要把原胞边长转化为约化单位即可。

总结

MD还是需要好好学习下的，这儿的很多细节是非常值得仔细思考思考。