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

类别:    标签: gmx   阅读次数:   版权: (CC) BY-NC-SA

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

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数组中的索引.

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

  1. 用双向图代表键和原子
  2. 令所有原子为白色
  3. 使白色原子中的一个原子(原子 \(i\))变为黑色, 并把它放在盒子中心
  4. \(i\) 的每个邻居, 如果它目前为白色, 则将其变为灰色
  5. 选择灰色原子中的一个(原子 \(j\)), 相对于它的所有黑色邻居, 给它正确的周期性, 并将其变黑
  6. \(j\) 的每个邻居, 如果它目前为白色, 则将其变为灰色
  7. 如果仍然存在任何一个灰色原子, 转到5
  8. 如果仍然存在任何一个白色原子, 转到3

使用这种算法, 我们可以

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]的三位点模型, 即:

  1. 分子中有三个原子.
  2. 整个分子属于单个电荷组.
  3. 第一个原子具有Lennard-Jones(4.1.1节)和库仑(4.1.3节)相互作用.
  4. 第二和第三个原子只具有库仑相互作用, 且电荷相等.

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

  1. 分子中有四个原子.
  2. 整个分子属于单个电荷组.
  3. 第一个原子只具有Lennard-Jones(4.1.1节)相互作用.
  4. 第二和第三个原子具有库仑相互作用, 且电荷相等.
  5. 第四个原子只具有库仑相互作用.

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

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}\]

图 B.1: IEEE单精度浮点数格式
图 B.1: IEEE单精度浮点数格式

通过计算:

\[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\) , 很明显, 解的第一次估计必须足够精确, 这样才能获得精确的结果. 计算的要求是

第一个要求指出, \(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}\]

这也可以使用整数运算进行精确计算. 注意, 对前面提到的范围分数部分是自动校正的, 因此不需要对指数部分进行额外的校正.

由此得到的结论是:

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\) 给出了具有最大精度的结果. 显然, 对双精度结果只需要一步额外的迭代(以双精度进行).


  1. \((\bi u \otimes \bi v)^{\a\b} = \bi u_\a \bi v_\b\)  ↩

  2. Lennard-Jones和库仑力时计算大约需要50次浮点运算.  ↩

随意赞赏

微信

支付宝
◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: 
后一篇: 

访问人次(2015年7月 9日起): | 最后更新: 2017-06-23 03:12:17 UTC | 版权所有 © 2008 - 2017 Jerkwin