约束处理

约束算法在可以接受的前提下增大 MD模拟的时间步长以减少其计算量的一种方法。

在 MD 中，时间步长通常选取分子体系中最短振动周期的 1/10 以保证能够在振动周期范围内正确取样，如果约束模拟体系中振动频率最快的键，则时间步长便可以选取次短振动周期的 1/10。

如 C-H键和 C-O 键的振动周期分别是 10 fs 和 20 fs，如果约束C-H 键，则 MD 模拟的时间步长可选取 2 fs，显然对同样的模拟体系，时间步长增大一倍相当于 MD 计算时间减少一半，即计算效率提升一倍。

SHAKE 算法是一种迭代算法，通常迭代次数较多，同时约束之间具有关联性，每次迭代中约束之间需要相互通信，导致通信次数较多，不太适合并行。

LINCS 算法是一种投射方法，迭代次数很少，一般迭代一次即可满足精度需求，相应的通信次数少，适合并行。

约束动力学是指固定分子中的键长、键角或二面角不变，但是粒子仍然具有平动或转动等运动，约束动力学既符合牛顿定律又符合几何约束。

最常用的是键长约束，即在沿键的方向连续施加大小相等方向相反的力以保持键长不变，如图 1 所示。

键长约束示意图

SHAKE

SHAKE算法将一组不受约束的坐标 r′ 转变为一组坐标 r″, 使其满足一系列距离约束条件, 使用一组 r 作为参考,

算法所起的作用与求解一组约束运动方程的Lagrange乘子相同. SHAKE需要指定 相对容差. 算法会继续直到所有的约束都处于相对容差之内. 若由于偏差太大, SHAKE不能重置坐标, 或迭代次数超过指定数目, 会给出错误信息.

假定运动方程必须满足 K 个完整可积分约束, 表述为

例如 (r1−r2)2−b2=0. 力的定义为:

其中 λk 为Lagrange乘子, 必须求解它们以满足约束方程. 上面求和的第二部分决定了 约束力 Gi, 定义为:

在蛙跳式或Verlet算法中, 由约束力引起的位移为 (Gi/mi)(Δt)2. 求解Lagrange乘子(进而位移)需要求解一组二度耦合的方程, 这可使用SHAKE方法迭代求解

讲解LINCS

LINCS 约束算法在 GROMACS 中是一种键约束的标准算法，该算法不是直接沿键的方向连续施加约束力，而是沿变化率最快的位移梯度方向连续施加约束力，对大部分模拟体系，LINCS 约束算法收敛很快，一般一次迭代即可提供足够的精度，关联约束之间的通信很少，适合并行。LINCS 约束算法中的约束方程使用矩阵形式表示，采用蛙跳(Leap-Frog)有限差分法积分后可得式(1)和式(2)，式中变量均为矩阵形式：r 是粒子的位移；r*是非约束距离；T 是旋转矩阵；B 是约束方程对位移的梯度;d 是约束距离；p 是约束在投影方向上的距离；I 是单位向量。

LINCS 算法是一种投射法(梯度法)，基本执行步骤是初始键先在非约束条件下得到非约束键长，再施加约束力对非约束键长不断纠正直至满足约束条件。如图 2所示，图中虚线代表初始键长（老键），实线代表新键。图 2a 是没有使用约束时得到的新键；图 2b 是图 2a 中的新键在沿老键方向的投射力 (约束力)作用下得到的新键(对应式(1))；图 2c 是对图 2b 中得到的新键继续施加沿老键方向的投射力作用得到的新键(对应式(1))，依次类推，每次迭代都是投射力沿老键方向上不断进行“投射作用”得到新的键长直至经过 z 次迭代满足给定的约束精度。

投射法示意图

constraints

控制拓扑中的哪些键（bonds）将转换为刚性完整约束（constraints）。请注意，典型的刚性水模型没有键，而是一个专门的 [settles] 指令，因此不受此选项的影响。

none: 除在拓扑中明确定义的键外, 不使用任何约束, 即, 利用简谐(或其他)势或Morse势(取决于morse的设置)描述键, 利用简谐(或其他)势描述键角.

h-bonds: 将含有氢原子的键转换为约束.

all-bonds: 将所有的键都转换为约束.

h-angles: 将所有涉及氢原子的键和键角都转换为键约束.

all-angles: 将所有的键和键角都转换为键约束.

constraint-algorithm

能够满足所有非 SETTLE 完整约束的求解算法。

LINCS: 线性约束求解器（LINear Constraint Solver）。通过区域分解，使用并行版本的 P-LINCS。

精度由 lincs-order 设置，该值设置了矩阵逆展开中的矩阵数。在矩阵求逆校正之后，该算法进行迭代校正，以补偿由于旋转引起的延长。

这种迭代的次数可以用 lincs-iter 控制。每 nstlog 步都会将均方根相对约束偏差输出到日志文件中。如果一个键在一个步骤中旋转超过 lincs-warnangle，则会在日志文件和 stderr 中输出一条警告。

LINC 不应与耦合角度约束一起使用。

SHAKE: SHAKE 比 LINCs 稍慢且不太稳定，但可用于角度约束。相对公差用 shake-tol 设置，0.0001 适用于一般 MD。

SHAKE 不支持不同分解域上原子之间的约束，因此只有在使用所谓的更新组（update-groups）时，它才能用于区域分解，这通常是仅约束涉及氢的键的情况。SHAKE 不能与能量最小化一起使用。

continuation

 此选项以前称为 unconstrained-start。

no: 将约束应用于启动配置并重置 shell。

yes: 不要将约束应用于开始配置，也不要重置 shell，这对于精确的二分法和重新运行非常有用。

shake-tol (0.0001)

SHAKE 的相对公差

Shake-SOR

yes or no.（不）使用使用连续超松弛减少抖动迭代次数。

lincs-order (4)

约束耦合矩阵展开中的最高阶。当约束形成三角形时，相同阶次的附加展开仅应用于此类三角形内的耦合的法向展开之上。对于一般MD模拟，4阶通常就够了，对于具有虚拟位点或 BD 的大时间步长，需要6阶。对于双精度的精确能量最小化，可能需要8阶或更多。请注意，在单精度中，由于舍入误差的放大，高于6的阶数通常会导致精度降低。通过区域分解，单元大小受到 lincs-order+1 个约束的跨越距离限制。如果想要扩展到这个极限之外，可以降低 lincs-order 并增加 lincs-iter，因为当 (1+lincs-iter)*lincs-order 保持不变时，精度不会下降。

lincs-iter (1)

在 LINCs 中校正旋转延长的迭代次数。对于正常运行，一次就够了，但是对于想要精确节约能量的 NVE 模拟，或者对于双精度下的精确能量最小化，您可能希望将其增加到2。请注意，在单精度下，使用多于1次迭代通常会因舍入误差的放大而导致更差的精度。

lincs-warnangle (30) [deg]

LINCS失效前键能够旋转的最大角度.

morse

no: 键由谐波势表示。

yes: 键由 Mores 势表示。