GROMACS中文手册：附录B　一些实现细节

本手册已过时, 不再更新. 如果需要最新手册, 请加入下方QQ群.

在本章中, 我们将介绍一些实现细节. 本章涉及的内容远远称不上完备, 但我们认为有必要澄清一些事情, 否则将很难理解.

维里 \(\X\) 可以写成全张量形式:

\[\X = -{1\over2}\Sum_{i < j}^N \bi r_{ij} \otimes \bi F_{ij} \tag{B.1}\]

其中 \(\otimes\) 表示两个向量的 直积.[1] 当此在MD程序的内部循环中计算直积时, 需要进行9次乘法运算和9加法运算.[2]

这里将展示如何才能将维里计算从内部循环中抽取出来[166].

对一个周期性体系, 计算维里时必须考虑周期性:

\[\X = -{1\over2}\Sum_{i < j}^N \bi r_{ij}^n \otimes \bi F_{ij} \tag{B.2}\]

其中 \(\bi r_{ij}^n\) 表示从原子 \(j\) 到原子 \(i\) 最近映象 的距离向量. 在这个定义中我们为原子 \(i\) 的位置向量 \(\bi r_i\) 添加一个 移位向量 \(\d_i\). 差向量 \(\bi r_{ij}^n\) 因此等于:

\[\bi r_{ij}^n = \bi r_i +\d_i - \bi r_j \tag{B.3}\]

或简写为:

\[\bi r_{ij}^n = \bi r_i^n - \bi r_j \tag{B.4}\]

对一个三斜体系, \(i\) 有27个可能的映象; 当使用截面八面体时, 有15个可能的映象.

这里给出非键力子程序中单个和维里的推导. 在下面的所有公式中 \(i \ne j\).

\(\alg \X &= -{1\over2}\Sum_{i < j}^N \bi r_{ij}^n \otimes \bi F_{ij} \tag{B.5} \\ &= -{1\over4} \Sum_{i=1}^N \Sum_{j=1}^N (\bi r_i +\d_i -\bi r_j) \otimes \bi F_{ij} \tag{B.6} \\ &= -{1\over4} \Sum_{i=1}^N \Sum_{j=1}^N (\bi r_i +\d_i ) \otimes \bi F_{ij} - \bi r_j \otimes \bi F_{ij} \tag{B.7} \\ &= -{1\over4} \left( \Sum_{i=1}^N \Sum_{j=1}^N (\bi r_i + \d_i) \otimes \bi F_{ij} - \Sum_{i=1}^N \Sum_{j=1}^N \bi r_j \otimes \bi F_{ij} \right) \tag{B.8} \\ &= -{1\over4} \left( \Sum_{i=1}^N (\bi r_i +\d_i) \otimes \Sum_{j=1}^N \bi F_{ij} - \Sum_{j=1}^N \bi r_j \otimes \Sum_{i=1}^N \bi F_{ij} \right) \tag{B.9} \\ &= -{1\over4} \left( \Sum_{i=1}^N (\bi r_i +\d_i) \otimes \bi F_i + \Sum_{j=1}^N \bi r_j \otimes \bi F_j \right) \tag{B.10} \\ &= -{1\over4} \left( 2\Sum_{i=1} \bi r_i \otimes \bi F_i + \Sum_{i=1} \d_i \otimes \bi F_i \right) \tag{B.11} \ealg\)

在上面这些公式中, 我们引入了

\[\bi F_i = \Sum_{j=1}^N \bi F_{ij} \tag{B.12}\]

\[\bi F_j = \Sum_{i=1}^N \bi F_{ji} \tag{B.13}\]

为 \(j\) 对 \(i\) 总的作用力. 因为我们使用了牛顿第三定律:

\[\bi F_{ij} = -\bi F_{ji} \tag{B.14}\]

在实现时我们必须将含移位 \(\d_i\) 的项加倍.

对键合力和SHAKE, 可以创建一个 mol-shift 列表, 其中存储了周期性. 我们简单地使用一个数组mshift, 其中存储了每个原子在shiftvec数组中的索引.

生成此列表的算法可以从图论得到, 将分子中的每个粒子视为图中的节点, 将键视为图的边.

使用这种算法, 我们可以

由于共价键力对维里有贡献, 我们有:

\(\alg b &= \| \bi r_{ij}^n \| \tag{B.15} \\ V_b &= {1\over2} k_b (b-b_0)^2 \tag{B.16} \\ \bi F_i &= - \nabla V_b \tag{B.17} \\ &= k_b (b-b_0) {\bi r_{ij}^n \over b} \tag{B.18} \\ \bi F_j &= - \bi F_i \tag{B.19} \ealg\)

来源于键的维里为:

\(\alg \X_b &= - {1\over2} (\bi r_i^n \otimes \bi F_i + \bi r_j \otimes \bi F_j ) \tag{B.20} \\ &= - {1\over2} (\bi r_{ij}^n \otimes \bi F_i) \tag{B.21} \ealg\)

SHAKE对维里有着重要贡献. 为满足约束, 力 \(\bi G\) 作用到“摇动”的粒子上. 如果此力不是来自算法(如在标准SHAKE中), 它可以在后面计算(当使用 蛙跳算法 时):

\(\alg \D \bi r_i &= \bi r_i (t+\D t) - [\bi r_i(t)+ \bi v_i(t-{\D t \over 2}) \D t+{\bi F_i \over m_i} \D t^2] \tag{B.22} \\ \bi G_i &= {m_i \D \bi r_i \over \D t^2} \tag{B.23} \ealg\)

在一般情况下, 上式对我们没有帮助. 只有当不需要使用周期性时(如刚性水), 才可以使用上面的公式, 否则我们必须在SHAKE的内部循环中增加维里的计算.

当 适用 时, 可以使用单个和方式计算维里:

\[\X = -{1\over2} \Sum_i^{N_c} \bi r_i \otimes \bi F_i \tag{B.24}\]

其中 \(N_c\) 为约束原子的数目.

在这里, 我们将描述GROMACS使用的一些算法优化, 不包括并行化. 对其中的一个, 1.0/sqrt(x)函数的实现, 我们将在B.3节分开处理. 其他最重要优化的论述如下.

GROMACS使用特殊的内部循环来计算水分子与其它原子的非键相互作用, 使用另一组循环计算水分子之间的相互作用. 这两组循环针对两种类型的水模型进行了高度优化. 对于类似于SPC[81]的三位点模型, 即:

这些循环也适用于SPC/E[167]和TIP3P[125]水模型. 对类似于TIP4P[125]的四位点水模型:

这些实现方式的好处是, 在单个循环中有更多的浮点运算, 这意味着一些编译器可以更好地调度代码. 然而, 事实证明, 甚至一些最先进的编译器也存在调度问题, 这意味着需要进行手动调整以获得最佳性能. 这可能包括消去相同的子表达, 或移动代码到各处.

不幸的是, 在一些平台上Fortran编译器仍好于C编译器. 对于一些机器(例如SGI Power Challenge)差异可高达3, 对矢量计算机差异可能更大. 因此, 针对英特尔和AMD的x86处理器, 有些占用大量计算时间的子程序被改写为Fortran甚至汇编代码. 在大多数情况下, 当适用时, Fortran或汇编循环会通过configure脚本自动选择, 但你也可以通过设置configure脚本的选项对此进行调整.

GROMACS项目开始于开发一个 \(1/\sqrt{x}\) 的处理器, 用以计算:

\[Y(x)= {1\over \sqrt x} \tag{B.25}\]

随着项目的继续, 英特尔i860处理器被用于实现GROMACS, 现在几乎已经变成了一个完整的软件项目. \(1/\sqrt x\) 处理器的实现采用了一步的Newton-Raphson迭代方案, 为此, 需要查表以提供初始近似值. \(1/\sqrt x\) 函数使用了两个几乎独立的表格, 分别用于IEEE浮点表示的指数种子和分数种子.

根据[168] \(1/\sqrt x\) 函数可以使用Newton-Raphson迭代方案进行计算. 反函数为:

\[X(y)= {1\over y^2} \tag{B.26}\]

因此不直接计算

\[Y(a)=q \tag{B.27}\]

而是使用Newton-Raphson方法求解方程

\[X(q) - a = 0 \tag{B.28}\]

通过计算:

\[y_{n+1}= y_n-{f(y_n) \over f'(y_n)} \tag{B.29}\]

进行迭代. 在这种近似下, 绝对误差 \(\ve\) 被定义为:

\[\ve \equiv y_n-q \tag{B.30}\]

利用Taylor级数展开来估计误差, 根据[168]的方程(3.2)得到:

\[\ve_{n+1} = - {\ve_n^2 \over 2}{f''(y_n) \over f'(y_n)} \tag{B.31}\]

这是绝对误差的估计.

IEEE 32位单精度格式的浮点数具有几乎恒定的相对误差 \(\D x/x = 2^{-24}\). 从前面所述的Taylor级数展开公式(方程B.31)可以看到, 每步迭代的误差是绝对的, 且一般与 \(y\) 有关. 如果将误差表示为相对误差, 有下面的式子

\[\ve_{r_{n+1} } \equiv {\ve_{n+1}\over y} \tag{B.32}\]

并且

\[\ve_{r_{n+1} } = - ({\ve_n \over y})^2 y {f'' \over 2f'} \tag{B.33}\]

对函数 \(f(y)= y^{-2}\), \(yf''/2f'\) 项为常数(等于–3/2), 因此相对误差 \(\ve_{r_n}\) 与 \(y\) 无关.

\[\ve_{r_{n+1} } = {3\over2} (\ve_{r_n})^2 \tag{B.34}\]

由此得到的结论是, 函数 \(1/\sqrt x\) 可以计算到指定的精度.

为了使用前面提到的迭代方案计算函数 \(1/\sqrt x\) , 很明显, 解的第一次估计必须足够精确, 这样才能获得精确的结果. 计算的要求是

第一个要求指出, \(1/\sqrt x\) 的结果可能具有的相对误差 \(\ve_r\) 等于IEEE 32位单精度浮点数的 \(\ve_r\). 由此可以得出初始近似的 \(1/\sqrt x\), 对后续步骤重写相对误差的定义(方程B.34):

\[{\ve_n \over y} = \sqrt{ \ve_{r_{n+1} }{2f' \over yf''} } \tag{B.35}\]

因此对于查表所需的精度为:

\[{\D Y\over Y} = \sqrt{\ {2\over3} 2^{-24}\ } \tag{B.36}\]

这定义了备查表的宽度必须≥13比特.

这样, 备查表的相对误差 \(\ve_{r_n}\) 是已知的, 由此可以计算参数的最大相对误差. 绝对误差 \(\D x\) 定义为:

\[\D x \equiv {\D Y \over Y'} \tag{B.37}\]

因此:

\[{\D x \over Y} = {\D Y \over Y} (Y')^{-1} \tag{B.38}\]

因此:

$$ \D x = \text{constant} {Y \over Y’} \tag{B.39}$$

对 \(1/\sqrt x\) 函数, 满足 \(Y/Y'~x\), 所以 \(\D x/x=\text{constant}\). 前面提到过, 这是所用浮点表示的性质. 备查表参数需要的精度符合:

\[{\D x \over x} = -2 {\D Y \over Y} \tag{B.40}\]

因此, 使用浮点精度(方程B.36):

\[{\D x \over x} = -2 \sqrt{ {2\over3}2^{-24} } \tag{B.41}\]

这定义了备查表的宽度必须≥12比特.

所使用的IEEE 32位单精度浮点格式规定, 一个数字由一个指数和一个分数表示. 上一节对每个可能的浮点数规定了备查表的长度和宽度. 精度仅由浮点数的分数部分的大小决定. 由此得到的结论是, 备查表的大小为其长度, 前面已规定, 再乘上指数的大小(21228, 1Mb). \(1/\sqrt x\) 函数具有指数部分与分数部分无关的性质. 如果使用浮点表示, 这很明显. 定义:

\[x \equiv (-1)^S (2^{E-127})(1.F) \tag{B.42}\]

参看图 B.1, 其中 \(0 \le S \le 1, 0 \le E \le 255, 1 \le 1.F \lt 2\), \(S, E, F\) 为整数(规范化条件). 符号位(\(S\))可以省略, 因为 \(1/\sqrt x\) 只对 \(x>0\) 有定义. \(1/\sqrt x\) 函数作用于 \(x\) 得到:

\[y(x)= {1 \over \sqrt x} \tag{B.43}\]

或:

\[y(x)= {1 \over \sqrt{(2^{E-127})(1.F)} } \tag{B.44}\]

这可以改写为:

\[y(x)=(2^{E-127})^{-1/2} (1.F)^{-1/2} \tag{B.45}\]

定义:

\[(2^{E'-127}) \equiv (2^{E-127})^{-1/2} \tag{B.46}\]

\[1.F' \equiv (1.F)^{-1/2} \tag{B.47}\]

这样 \({1\over\sqrt 2} \lt 1.F' \le 1\) 成立, 因此对规范的实数表示非常重要的条件 \(1\le 1.F' \lt 2\), 不再成立. 通过引入一个额外的项, 可以对此进行校正. 将 \(1/\sqrt x\) 函数应用于浮点数(方程B.45)改写为:

\[y(x)=(2^{ {127-E\over2}-1})(2(1.F)^{-1/2}) \tag{B.48}\]

和:

\[(2^{E'-127}) \equiv (2^{ {127-E\over2} -1} ) \tag{B.49}\]

\[1.F' \equiv 2(1.F)^{-1/2} \tag{B.50}\]

这样 \(\sqrt 2 \lt 1.F \le 2\) 成立. 这和方程B.42中定义的规范化浮点数的精确有效范围不同. 数值2导致问题. 通过将此值映射到