本手册已过时, 不再更新. 如果需要最新手册, 请加入下方QQ群.
在本章中, 我们将介绍一些实现细节. 本章涉及的内容远远称不上完备, 但我们认为有必要澄清一些事情, 否则将很难理解.
B.1 GROMACS中的单个和维里
维里 \(\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].
B.1.1 维里
对一个周期性体系, 计算维里时必须考虑周期性:
\[\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个可能的映象.
B.1.2 非键力的维里
这里给出非键力子程序中单个和维里的推导. 在下面的所有公式中 \(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\) 的项加倍.
B.1.3 分子内移位(mol-shift)
对键合力和SHAKE, 可以创建一个 mol-shift 列表, 其中存储了周期性. 我们简单地使用一个数组mshift, 其中存储了每个原子在shiftvec数组中的索引.
生成此列表的算法可以从图论得到, 将分子中的每个粒子视为图中的节点, 将键视为图的边.
- 用双向图代表键和原子
- 令所有原子为白色
- 使白色原子中的一个原子(原子 \(i\))变为黑色, 并把它放在盒子中心
- 对 \(i\) 的每个邻居, 如果它目前为白色, 则将其变为灰色
- 选择灰色原子中的一个(原子 \(j\)), 相对于它的所有黑色邻居, 给它正确的周期性, 并将其变黑
- 对 \(j\) 的每个邻居, 如果它目前为白色, 则将其变为灰色
- 如果仍然存在任何一个灰色原子, 转到5
- 如果仍然存在任何一个白色原子, 转到3
使用这种算法, 我们可以
- 优化键合力计算以及SHAKE
- 使用单个和方法从键合力计算维里
- 获得键的双向图表示.
B.1.4 共价键的维里
由于共价键力对维里有贡献, 我们有:
\(\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\)
B.1.5 SHAKE的维里
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\) 为约束原子的数目.
B.2 优化
在这里, 我们将描述GROMACS使用的一些算法优化, 不包括并行化. 对其中的一个, 1.0/sqrt(x)函数的实现, 我们将在B.3节分开处理. 其他最重要优化的论述如下.
B.2.1 水的内部循环
GROMACS使用特殊的内部循环来计算水分子与其它原子的非键相互作用, 使用另一组循环计算水分子之间的相互作用. 这两组循环针对两种类型的水模型进行了高度优化. 对于类似于SPC[81]的三位点模型, 即:
- 分子中有三个原子.
- 整个分子属于单个电荷组.
- 第一个原子具有Lennard-Jones(4.1.1节)和库仑(4.1.3节)相互作用.
- 第二和第三个原子只具有库仑相互作用, 且电荷相等.
这些循环也适用于SPC/E[167]和TIP3P[125]水模型. 对类似于TIP4P[125]的四位点水模型:
- 分子中有四个原子.
- 整个分子属于单个电荷组.
- 第一个原子只具有Lennard-Jones(4.1.1节)相互作用.
- 第二和第三个原子具有库仑相互作用, 且电荷相等.
- 第四个原子只具有库仑相互作用.
这些实现方式的好处是, 在单个循环中有更多的浮点运算, 这意味着一些编译器可以更好地调度代码. 然而, 事实证明, 甚至一些最先进的编译器也存在调度问题, 这意味着需要进行手动调整以获得最佳性能. 这可能包括消去相同的子表达, 或移动代码到各处.
B.2.2 Fortran代码
不幸的是, 在一些平台上Fortran编译器仍好于C编译器. 对于一些机器(例如SGI Power Challenge)差异可高达3, 对矢量计算机差异可能更大. 因此, 针对英特尔和AMD的x86处理器, 有些占用大量计算时间的子程序被改写为Fortran甚至汇编代码. 在大多数情况下, 当适用时, Fortran或汇编循环会通过configure脚本自动选择, 但你也可以通过设置configure脚本的选项对此进行调整.
B.3 1.0/sqrt函数的计算
B.3.1 简介
GROMACS项目开始于开发一个 \(1/\sqrt{x}\) 的处理器, 用以计算:
\[Y(x)= {1\over \sqrt x} \tag{B.25}\]
随着项目的继续, 英特尔i860处理器被用于实现GROMACS, 现在几乎已经变成了一个完整的软件项目. \(1/\sqrt x\) 处理器的实现采用了一步的Newton-Raphson迭代方案, 为此, 需要查表以提供初始近似值. \(1/\sqrt x\) 函数使用了两个几乎独立的表格, 分别用于IEEE浮点表示的指数种子和分数种子.
B.3.2 通用
根据[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}\]
这是绝对误差的估计.
B.3.3 用于于浮点数
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\) 可以计算到指定的精度.
B.3.4 备查表格的要求
为了使用前面提到的迭代方案计算函数 \(1/\sqrt x\) , 很明显, 解的第一次估计必须足够精确, 这样才能获得精确的结果. 计算的要求是
- IEEE格式的最大可能精度
- 只使用一次迭代以达到最高速度
第一个要求指出, \(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比特.
B.3.5 指数和分数的独立计算
所使用的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导致问题. 通过将此值映射到<2的最近表示, 可以解决这个问题. 这种近似引入的小误差在允许范围内.
指数的整数表示是下一个问题. 计算 \((2^{ {127-E\over2} -1})\) 会引入分数的结果, 如果 \((127 - E)\) 为奇数. 通过将计算分为奇数部分和偶数部分, 很容易处理这个问题. 对 \((127 - E)\) 为偶数的情况, 方程(方程B.49)中的 \(E'\) 作为 \(E\) 的函数, 可以使用整数运算精确地计算出来
\[E' = {127 -E \over 2} + 126 \tag{B.51}\]
对 \((127 - E)\) 为奇数的情况, 方程(方程B.45)可以改写为:
\[y(x)=(2^{ {127-E-1\over2} }) ({1.F \over 2})^{-1/2} \tag{B.52}\]
因此:
\[E' = {126 - E \over 2} + 127 \tag{B.53}\]
这也可以使用整数运算进行精确计算. 注意, 对前面提到的范围分数部分是自动校正的, 因此不需要对指数部分进行额外的校正.
由此得到的结论是:
- 分数和指数的备查表是独立的. 存在两个分数备查表(奇指数和偶指数), 因此必须使用指数的奇/偶信息(1sb比特)选择适当的表格.
- 指数备查表是一个256 x 8比特的表, 对 奇数 和 偶数 进行了初始化.
B.3.6 实现
可以使用一个小的C程序来产生备查表, 它使用了符合IEEE 32位单精度格式的浮点数和操作. 需要注意的是, 因为需要 奇/偶 信息, 分数表是前面指定的两倍大(13比特 i.s.o 12比特).
必须实现方程B.29的函数. 应用到 \(1/\sqrt x\) 函数, 方程B.28给出:
\[f = a -{1\over y^2} \tag{B.54}\]
因此:
\[f' = {2\over y^3} \tag{B.55}\]
因此:
\[y_{n+1} = y_n -{a-{1\over y_n^2} \over {2\over y_n^3} } \tag{B.56}\]
或:
\[y_{n+1} = {y_n \over 2} (3 -a y_n^2) \tag{B.57}\]
其中 \(y_0\) 可以在备查表中找到, \(y_1\) 给出了具有最大精度的结果. 显然, 对双精度结果只需要一步额外的迭代(以双精度进行).
