分子动力学模拟之SETTLE约束算法

Python微信订餐小程序课程视频

https://edu.csdn.net/course/detail/36074

Python实战量化交易理财系统

https://edu.csdn.net/course/detail/35475目录

• 技术背景
• 算法流程
  • 约束前后
  • 约束算法中的约束条件
  • 建立新的坐标系
  • 三维旋转
• 算法解析表达式
• 三维空间坐标变换
• 坐标变换代码实现
• SETTLE代码实现
• 总结概要
• 版权声明
• 参考文章

---

回到顶部# 技术背景

在上一篇文章中，我们讨论了在分子动力学里面使用LINCS约束算法及其在具备自动微分能力的Jax框架下的代码实现。约束算法，在分子动力学模拟的过程中时常会使用到，用于固定一些既定的成键关系。例如LINCS算法一般用于固定分子体系中的键长关系，而本文将要提到的SETTLE算法，常用于固定一个构成三角形的体系，最常见的就是水分子体系。对于一个水分子而言，O-H键的键长在模拟的过程中可以固定，H-H的长度，或者我们更常见的作为一个H-O-H的夹角出现的参量，也需要固定。纯粹从计算量来考虑的话，RATTLE约束算法需要迭代计算，LINCS算法需要求矩阵逆（虽然已经给出了截断优化的算法），而SETTLE只涉及到坐标变换，显然SETTLE在约束大规模的水盒子时，性能会更加优秀。

回到顶部# 算法流程

类似于LINCS算法的，我们先引入一个核心的约束条件：

|rij|2−d2ij=0|rij|2−d2ij=0\left|r_{ij}\right|^2-d_{ij}2=0
意思就是说，原子i和原子j之间的距离在分子动力学模拟的迭代过程中保持不变。dijd_{ij}可以是初始值，也可以是一个给定值。在满足这个约束条件的同时，其实隐藏着另外一个约束条件：

rij⋅vij=0r_{ij}\cdot v_{ij}=0
换而言之，如果所有原子的运动方向都与其直接的键连方向垂直，那么在模拟迭代的过程中键长距离就不会发生改变了。

约束前后

我们还是需要从一个简单的三角形A0B0C0A_0B_0C_0模型出发来具体讲解这个算法的流程：

假定三角形A0B0C0A_0B_0C_0为原始位置，三角形A1B1C1A_1B_1C_1为没有加约束的坐标更新，而我们的目标是得到三角形A3B3C3A_3B_3C_3这个加了约束的三角形。这里先解释一下为什么不是三角形A2B2C2A_2B_2C_2而是三角形A3B3C3A_3B_3C_3，因为这里三角形的变换只是涉及到加了约束条件和没加约束条件的结果，但是实际上加约束条件是一个步骤比较复杂的过程，这里写成三角形A3B3C3A_3B_3C_3是为了方便跟后续的SETTLE约束所得到的结果下标对齐。

约束算法中的约束条件

在上一步的算法过程中，其实我们已经初步的把施加SETTLE约束算法的过程划分成了两个部分：先不加约束的演化，再考虑施加约束的演化。之所以能够这么划分，是因为我们把约束条件考虑成一个约束作用力，这样有两个好处：一是约束力也变成连续可微的参量，二是矢量的作用是可以叠加和拆分的，这里两步的拆分，就是用到了其矢量作用力的效果。

将三角形初始所在的平面定义为π0\pi_0，在经过未加约束的偏移之后得到了三角形A1B1C1A_1B_1C_1，此时可能已经偏离了原始的平面π0\pi_0，这里将三个顶点所在的位置都构造一个与π0\pi_0相互平行的平面，称为πA\pi_A、πB\pi_B和πC\pi_C。这里可能会有两个需要着重注意的点，第一，我们之所以要定义πA\pi_A、πB\pi_B和πC\pi_C这三个平面，是因为约束力总是在径向的，也就是说，不可能产生与平面π0\pi_0相垂直的约束力，因此约束力的作用只是让三个顶点在与π0\pi_0平面相平行的平面上运动，也就是这里的πA\pi_A、πB\pi_B和πC\pi_C三个平面。换而言之，三角形A3B3C3A_3B_3C_3的三个顶点必然在πA\pi_A、πB\pi_B和πC\pi_C三个平面内。第二，因为约束力是分子内部的作用力，也就是对重心的合力为0，不会导致整体分子的漂移，因此，在施加约束力前后，重心的位置不变。如果分别用D1D_1和D3D_3来标记三角形A1B1C1A_1B_1C_1和A3B3C3A_3B_3C_3的重心的话，那么有D1=D3D_1=D_3。并且，假如原始位置三角形A0B0C0A_0B_0C_0的重心d0d_0跟D1D_1和D3D_3如果是在同一个位置的话，那么原始位置的三角形A0B0C0A_0B_0C_0就可以通过绕重心的旋转跟三角形A3B3C3A_3B_3C_3重合。

建立新的坐标系

基于d0=D1=D3d_0=D_1=D_3的假定，我们可以基于三角形A0B0C0A_0B_0C_0建立这样一个新的坐标系X′Y′Z′X’Y’Z’，原点为d0d_0，X′−Y′X’-Y’平面与π0\pi_0平行，而Y′−Z′Y’-Z’平面过A0A_0点：

如果从Z′Z’轴向Z′Z’轴负方向看（常规理解就是从上往下看的俯视图），就是这样的一个平面：

三维旋转

在前面的章节中我们提到，通过旋转操作，可以将初始的坐标旋转到更施加约束之后的坐标完全重合的位置，因此我们先假设三个旋转角，对原始坐标进行旋转操作。最后再根据约束条件来计算对应的旋转角，进而得到施加约束之后的新的坐标，也就是我们最终所需要的结果。在新的坐标系下，我们把三角形A0B0C0A_0B_0C_0重新标记为三角形a0b0c0a_0b_0c_0，接下来开始对三角形a0b0c0a_0b_0c_0的一系列三维旋转操作：

在初始条件下，因为三角形跟X′−Y′X’-Y’处在同一个平面上，因此从Y′Y’轴向Y′Y’轴正方向看时，会看到一条直线，此时我们绕Y′Y’轴旋转一个ψ\psi的角度，得到了三角形a1b1c1a_1b_1c_1。

按照同样的操作，绕X′X’轴旋转ϕ\phi以及绕Z′Z’轴旋转θ\theta的角度，经过三次的旋转之后，得到一个新的三角形a3b3c3a_3b_3c_3。而此时的三角形a3b3c3a_3b_3c_3正是处于跟三角形A3B3C3A_3B_3C_3完全重合的状态。因此，我们就可以根据约束条件计算出来三个欧拉角ψ\psi、ϕ\phi和θ\theta，进而得到我们所需要的约束后的结果：三角形A3B3C3A_3B_3C_3。

回到顶部# 算法解析表