博文

ROMACS中文手册第4章：相互作用函数和力场 1-2节

已有 10479 次阅读 2015-1-28 08:16 |个人分类:成神颇难|系统分类:科研笔记

第四章　　相互作用函数和力场

感谢网友阮洋帮忙录入公式

为了能够使用一些流行力场中采用的势能函数(见4.10节)，GROMACS提供了一些势能函数，包括非键相互作用以及二面角的相互作用。我们将在适当的小节中对它们予以讨论$。$

势能函数可以被划分成三个部分：

非键项: Lennard-Jones或Buckingham，以及库仑或修正的库仑。非键相互作用的计算基于已移除了排除项的邻区列表(一定半径内非键原子的列表)。
键合项: 共价键伸缩，键角弯曲，反常二面角，和正常二面角。这些项的计算是基于固定的列表。
约束项: 位置约束，角度约束，距离约束，方向约束和二面角约束，均基于固定的列表进行计算。

4.1 非键相互作用

GROMACS的非键相互作用是对势累加的(pair-additive)，并满足中心对称：

V (r 1, \dots r N) = \sum i < j V i j (r i j) (4.1)

F i = - \sum j d V i j ( r i j ) d r i j r i j r i j = - F j (4.2)

非键相互作用包含排斥项，色散项和库仑项。组合起来的排斥项和色散项可以取Lennard-Jones(或6-12相互作用)或Buckingham(或exp-6势)形式。此外，(部分)带电的原子之间的非键相互作用通过库仑项表达。

图4.1：Lennard-Jones相互作用4.1.1 Lennard-Jones相互作用

两个原子之间Lennard-Jones势 VLJ(rij) 等于：

V L J (r i j) = C ( 12 ) i j r 12 i j - C ( 6 ) i j r 6 i j (4.3)

其函数曲线如图4.1所示。参数 C(12)ij 和 C(6)ij 取决于配对的原子类型; 因此来源于LJ-参数矩阵。在Verlet截断方案中，通过将势能移动一个恒量，可使其在截断距离处为零。

此势能函数的力是：

F i (r i j) = ⎛ ⎝ 12 C ( 12 ) i j r 13 i j - 6 C ( 6 ) i j r 7 i j ⎞ ⎠ r i j r i j (4.4)

LJ势也可以写成以下形式：

V L J (r i j) = 4 ε i j ⎛ ⎝ (σ i j r i j) 12 - (σ i j r i j) 6 ⎞ ⎠ (4.5)

在构造非键LJ-参数的参数矩阵时，GROMACS可使用两种不同的组合规则，只采用几何平均值（力场文件输入部分的类型1）：

C (6) i j C (12) i j = (C (6) i i C (6) j j) 1 / 2 = (C (12) i i C (12) j j) 1 / 2 (4.6)

或采用Lorentz-Berthelot规则。此规则利用算术平均计算 σij, 利用几何平均计算 εij （类型2）：

σ i j ε i j = 1 2 (σ i i + σ j j) = (ε i i ε j j) 1 / 2 (4.7)

最后，可以对两个参数都采用几何平均值（类型3）：

σ i j = (σ i i σ j j) 1 / 2 ε i j = (ε i i ε j j) 1 / 2 (4.8)

OPLS力场就采用这个规则。

图4.2：Buckingham相互作用4.1.2 Buckingham势

与Lennard-Jones势相比，Buckingham势的排斥项更加灵活，贴近实际，但其计算也更昂贵。Buckingham势的形式是：

V b h (r i j) = A i j exp (- B i j r i j) - C i j r 6 i j (4.9)

请参见图4.2。此势能函数的力为：

F i (r i j) = [A i j B i j exp (- B i j r i j) - 6 C i j r 7 i j] r i j r i j (4.10)

4.1.3 库仑相互作用

两个带电粒子之间的库仑相互作用由下式给出：

V c (r i j) = f q i q j ε r r i j (4.11)

参见图4.3，其中 f=14πε0=138.935485 （见第2章）

此势对应的力为：

F i (r i j) = f q i q j ε r r 2 i j r i j r i j (4.12)

$图4.3：反应场与非反应场的库仑相互作用（带同种电荷的粒子之间）。在非反应场的情况下 $.ve_r$ 为1， $.ve_{rf}$ 为78， $r_c$ 为0.9纳米。除相差一个常数外，点划线与虚线完全相同。$ 图4.3：反应场与非反应场的库仑相互作用（带同种电荷的粒子之间）。在非反应场的情况下 εr 为1， εrf 为78， rc 为0.9纳米。除相差一个常数外，点划线与虚线完全相同。

使用普通的库仑相互作用时要么不能截断，要么所有粒子对都处于截断半径范围内，因为力在截断处存在跃变，而且变化很大。如果你确实需要使用截断，可以将势能进行移动以使势能等于力的积分。在基于组的截断方案中，这种移动只适用于非排除对。在Verlet截断方案中，移动也适用于排除对和自相互作用，这使得势能等同于 εrf=1 的反应场（见下文）。

GROMACS中相对介电常数 εr 可在grompp的输入中设置。

4.1.4 反应场库仑作用

对均相体系，可通过假定截断半径 rc 以外的粒子周围是介电常数为 εrf 的恒定介电环境，对库仑相互作用进行校正。此时，库仑相互作用可写为：

V c r f = f q i q j ε r r i j [1 + ε r f - ε r 2 ε r f + ε r r 3 i j r 3 c] - f q i q j ε r r c 3 ε r f 2 ε r f + ε r (4.13)

上式右边的常量表达式使得在截断处势能为零。对带电的截断球体，上式对应于均匀背景电荷的中性化。我们可以将方程4.13简写为

V c r f = f q i q j ε r [1 r i j + k r f r 2 i j - c r f] (4.14)

其中,

k r f c r f = 1 r 3 c ε r f - ε r ( 2 ε r f + ε r ) = 1 r c + k r f r 2 c = 1 r c 3 ε r f ( 2 ε r f + ε r ) (4.15) (4.16)

对大的 εrf, krf 趋近 r−3c/2, 而当 εrf=εr 时校正消失。图4.3中画出了校正的相互作用，很显然，对 rij的导数(即负的力)在截断距离处为零。这种势函数对应的力为：

F i (r i j) = f q i q j ε r [1 r 2 i j - 2 k r f r i j] r i j r i j (4.17)

反应场校正也可以用于所有的排除原子对，包括自身对，在这种情况下，方程4.13和4.17中不出现正常的库仑项。

Tironi等引入了更一般的反应场理论，其截断半径 rc 以外的电介质也具有离子强度 I[75]。这种情况下，我们可以利用逆Debye屏蔽长度 κ, 将常数 krf 和 crf 重写为

κ 2 k r f c r f = 2 I F 2 ε 0 ε r f R T = F 2 ε 0 ε r f R T \sum i = 1 K c i z 2 i = 1 r 3 c ( ε r f - ε r ) ( 1 + κ r c ) + 1 2 ε r f ( κ r c ) 2 ( 2 ε r f + ε r ) ( 1 + κ r c ) + ε r f ( κ r c ) 2 = 1 r c 3 ε r f ( 1 + κ r c + 1 2 ( κ r c ) 2 ) ( 2 ε r f + ε r ) ( 1 + κ r c ) + ε r f ( κ r c ) 2 (4.18) (4.19) (4.20)

其中, F 是法拉第常数, R 是理想气体常数, T 是绝对温度, ci 是物种 i 的摩尔浓度, zi 是物种 i 的电荷数，共有 K 种不同的物种。在离子强度为零的极限情况下，方程4.19和4.20分别退化为简单形式的方程4.15和4.16。

4.1.5 修正的非键相互作用

在GROMACS中，可以利用移位函数对非键势能函数进行修改，其目的在于用连续的力取代截断的力，且使其在截断半径处具有连续的导数。使用这种力时，每步积分引入的误差更小，且避免了一些复杂情况，比如在截断过程中由偶极产生的电荷。实际上，使用移位力时，不需要对电荷组创建邻近列表。然而，移位函数使库仑势产生了一个相当大的变形。除非对“丢失”的远程势进行了正确的计算以及加和（通过使用PPPM，Ewald或PME），否则必须仔细评估这些修改的效果。对Lennard-Jones色散项和排斥项的修改较小，但它确实消除了截断效应引起的噪声。

切换函数（将势能乘上一个函数）和移位函数（将势能或力加上一个函数）之间没有本质区别[76]。切换函数是移位函数的一种特殊情况。将切换函数应用于与粒子 j 作用在粒子 i 上的静电力或范德华力有关的力函数 F(r) ：

F i = c F (r i j) r i j r i j (4.21)

对于纯的库仑或Lennard-Jones相互作用, F(r)=Fα(r)=r−α+1. 移位力 Fs(r) 一般可写为：

F s (r) F s (r) F s (r) = F α (r) = F α (r) + S (r) = 0 r < r 1 r 1 \leq r \leq r c r c \leq r (4.22)

当 r1 =0 时，这是一个通常的移动函数，否则它就是一个切换函数。相应的移位库仑势为：

V s (r i j) = f Φ s (r i j) q i q j (4.23)

其中, Φ(r) 是势函数

Φ s (r) = \int \infty r F s (x) d x (4.24)

GROMACS的移位函数在边界处应该是平滑的，因此移位函数应当满足以下边界条件：

S (r 1) S' (r 1) S (r c) S' (r c) = 0 = 0 = - F α (r c) = - F' α (r c) (4.25)

三阶多项式形式

S (r) = A (r - r 1) 2 + B (r - r 1) 3 (4.26)

即满足这些要求。常数 A 和 B 可利用 rc 处的边界条件求出

A B = - ( α + 4 ) r c - ( α + 1 ) r 1 r α + 2 c ( r c - r 1 ) 2 = ( α + 3 ) r c - ( α + 1 ) r 1 r α + 2 c ( r c - r 1 ) 3 (4.27)

因此，整个力函数为：

F s (r) = α r α + 1 + A (r - r 1) 2 + B (r - r 1) 3 (4.28)

相应的势函数为：

Φ (r) = 1 r α - A 3 (r - r 1) 3 - B 4 (r - r 1) 4 - C (4.29)

其中

C = 1 r α c - A 3 (r c - r 1) 3 - B 4 (r c - r 1) 4 (4.30)

图4.4：库仑力，移位力与移位函数 $S(r)$, 其中 $r_1=2, r_c=4$. 图4.4：库仑力，移位力与移位函数 S(r), 其中 r1=2,rc=4.

当 r1 =0 时，修正的库仑力函数

F s (r) = 1 r 2 - 5 r 2 r 4 c + 4 r 3 r 5 c (4.31)

与推荐使用的短程函数加上长程部分的泊松解是相同的[77]。修正的库仑势函数为：

Φ (r) = 1 r - 5 3 r c + 5 r 3 3 r 4 c - r 4 r 5 c (4.32)

参见图 4.4。

4.1.6 修改的Ewald求和短程相互作用

当利用Ewald求和方法或粒子网格Ewald方法计算长程相互作用时，短程库仑势必须进行修改，这与上面切换函数的情形类似。在这种情况下，短程势能由下式给出：

V (r) = f erfc ( β r i j ) r i j q i q j (4.33)

其中, β 参数用于决定直接空间和与倒易空间和之间的相对比重, erfc(x) 是余误差函数。关于长程静电相互作用的详细讨论见4.8节。

图4.5：键伸缩的原理（左）以及键伸缩的势能（右）。4.2 键合相互作用

键合相互作用基于固定的原子列表。它们不仅包括对相互作用(pair interactions)，也包括三体和四体相互作用。键合相互作用包括键伸缩(双体)，键角(三体)和二面角(四体)相互作用。一种特殊类型的二面角相互作用(称为异常二面角)可用于将原子保持在同一平面上，或是避免将分子转变为相反手性（镜像）的结构。

4.2.1 键伸缩

简谐势

两个共价键合原子 i 和 j 之间的键伸缩可利用简谐势进行描述

V b (r i j) = 1 2 k b i j (r i j - b i j) 2 (4.34)

参见图4.5，相应的力如下：

F i (r i j) = k b i j (r i j - b i j) r i j r i j (4.35)

四次幂势

在GROMOS-96力场[78]中，为提高计算效率，共价键势记为以下形式：

V b (r i j) = 1 4 k b i j (r 2 i j - b 2 i j) 2 (4.36)

相应的力为：

F i (r i j) = k b i j (r 2 i j - b 2 i j) r i j (4.37)

对这种形式的势函数，其力常数与通常的简谐力常数 kb,harm (4.2.1节) 存在如下关系

2 k b b 2 i j = k b, harm (4.38)

力常数大多是由GROMOS-87使用的简谐力常数导出[79]。尽管计算时这种形式效率更高（因为不需要计算平方根），但它在概念上比较复杂。它还存在一个特别的缺点，由于不是简谐势形式，单个键的平均能量不等于 12kT ，而在正常的简谐势情况下是 12kT.

4.2.2 Morse势键伸缩

对需要非简谐键伸缩势的一些体系，GROMACS提供了两个原子 i 和 j 之间的Morse势[80]。这个势能函数与简谐势的不同在于，它具有一个不对称的势阱，且无限远处的力为零。函数形式是：

V m o r s e (r i j) = D i j [1 - exp (- β i j (r i j - b i j))] 2 (4.39)

参见图4.6，相应的力是：

F m o r s e (r i j) = 2 D i j β i j r i j exp (- β i j (r i j - b i j)) * [1 - exp (- β i j (r i j - b i j))] r i j r i j (4.40)

其中 Dij 是势阱深度，以 kJ/mol 为单位， βij 定义了势阱的陡度（以 nm^{-1} 为单位）， bij 为平衡距离(以 nm 为单位)。陡度参数 βij 可用原子 i 和 j 的约化质量、基本振动频率 μij 和势阱深度 Dij 表示：

β i j = ω i j μ i j 2 D i j - - - - - \sqrt (4.41)

且因为 ω=kμ−−√ ，可将 βij 改用简谐力常数 kij 表示：

β i j = k i j 2 D i j - - - - - \sqrt (4.42)

对于小的偏差 (rij−bij) ，可利用泰勒展开，取指数项的一阶近似：

exp (- x) \approx 1 - x (4.43)

将方程 4.42和4.43代入到函数中：

V m o r s e (r i j) = D i j [1 - exp (- β i j (r i j - b i j))] 2 = D i j [1 - (1 - k i j 2 D i j - - - - - \sqrt (r i j - b i j))] 2 = 1 2 k i j (r i j - b i j) 2 (4.44)

就得到了简谐键伸缩势。

图4.6：Morse势阱，键长为0.15nm。4.2.3 立方键伸缩势

另一种比Morse势稍简单的非简谐键伸缩势是在简单的简谐形式上增加了距离的一个立方项：

V b (r i j) = k b i j (r i j - b i j) 2 + k b i j k c u b i j (r i j - b i j) 3 (4.45)

在Ferguson[82]发展的一个柔性水模型（基于SPC水模型[81]）中，OH键使用了立方键伸缩势。该模型给出了合理的红外光谱。GROMACS库中提供了Ferguson水模型(flexwat-ferguson.itp)。应当指出，这种势函数是不对称的：过分拉伸导致能量无限降低。因此积分的时间步长被限制为1 fs。

该势对应的力为：

F i (r i j) = 2 k b i j (r i j - b i j) r i j r i j + 3 k b i j k c u b i j (r i j - b i j) 2 r i j r i j (4.46)

4.2.4 FENE键伸缩势

在聚合物的粗粒化模拟中，珠子之间通常由FENE(finitely extensible nonlinear elastic, 有限扩展非线弹性)势[83]连接：

V FENE (r i j) = - 1 2 k b i j b 2 i j log (1 - r 2 i j b 2 i j) (4.47)

这种势能函数看起来很复杂，但力的表达式更简单：

F FENE (r i j) = - k b i j (1 - r 2 i j b 2 i j) - 1 r i j (4.48)

在短距离处此势不对称地接近力常数为 kb 的简谐势，而在距离为 b 时发散。

图4.7：角振动原理（左）和键角势（右）。4.2.5 简谐角势

三个原子 i−j−k 之间的键角振动也可以用键角 θijk 的简谐势描述

V a (θ i j k) = 1 2 k θ i j k (θ i j k - θ 0 i j k) 2 (4.49)

由于键角的振动由简谐势描述，其形式与键伸缩是一样的（图4.5）。

利用链式规则，可得到力的方程：

F i F k F j = - d V a ( θ i j k ) d r i = - d V a ( θ i j k ) d r k 其 中 θ i j k = arccos ( r i j \cdot r k j ) r i j r k j = - F i - F k (4.50)

编号 i,j,k 是共价键原子的序列。原子 j 处于中间；原子 i 和 k 处于顶端（见图4.7）。注意，在拓扑文件的输入中，键角以度为单位，力常数以 kJ/mol/{rad}2 为单位。

4.2.6 余弦键角势

GROMOS-96力场利用一个简化的函数表示角振动：

V a (θ i j k) = 1 2 k θ i j k (cos (θ i j k) - cos (θ 0 i j k)) 2 (4.51)

其中

cos (θ i j k) = r i j \cdot r k j r i j r k j (4.52)

相应的力可通过对原子位置的偏微分推导出来。这个函数中的力常数与简谐势的力常数 kθ,harm (4.2.5) 存在如下关系

k θ sin 2 (θ 0 i j k) = k θ, harm (4.53)

$图4.8：弯曲键角势：余弦简谐（黑实线），简谐角（虚黑线）和限制弯曲（红色），三者具有相同的弯曲常数 $k_{.theta} = 85 .text{kJ mol}^{-1}$ 和平衡角 $.theta_0 = 130^{.circ}$ 。橙色线表示余弦简谐($k = 50 .text{kJ mol}^{-1}$)与限制弯曲($k = 25 .text{kJ mol}^{-1}$)的和，二者的 $.theta_0 = 130^{.circ}$ 。$ 图4.8：弯曲键角势：余弦简谐（黑实线），简谐角（虚黑线）和限制弯曲（红色），三者具有相同的弯曲常数 kθ=85kJ mol−1 和平衡角 θ0=130° 。橙色线表示余弦简谐(k=50kJ mol−1)与限制弯曲(k=25kJ mol−1)的和，二者的 θ0=130° 。

在GROMOS-96手册中有一个更加复杂的，与温度相关的转换公式。这些公式在0 K时是等效的，300 K时的差异大约为0.1%到0.2%。注意在拓扑文件的输入中，角度的单位为度，力常数的单位为kJ/mol。

4.2.7 受限弯曲势

受限弯曲（ReB）势[84]可防止弯曲角度 θ 达到 180° 。这样，在进行粗粒分子动力学模拟时，避免了计算扭转角和势能带来的数值不稳定性。

为了系统地防止弯曲角度达到 180° ，在弯曲势4.51的分母项中引入了 sin2θ 因子

V ReB (θ i) = 1 2 k θ ( cos θ i - cos θ 0 ) 2 sin 2 θ i (4.54)

图4.8显示了ReB势4.54和标准势4.51之间的比较。ReB势的壁垒在靠近 180° 的区域非常排斥，因此，弯曲角能够保持在一定的安全区间内，远离不稳定性。之所以在分母中使用 sinθ 的二次方项是为了保证这种行为，并且能给出优雅的微分表达式：

F ReB (θ i) = 2 k θ sin 4 θ i (cos θ i - cos θ 0) (1 - cos θ i cos θ 0) \partial cos θ i \partial r ⃗ k (4.55)

因其构造，限制弯曲势不能用于平衡的 θ0 值太接近 0° 或 180° 的情形（根据经验，推荐至少有 10° 的差异）。重要的是，在初始构型中，所有的弯曲角度都必须处在安全的区间内，以避免初始的不稳定性。这种弯曲势可以与任何形式的扭转势联合使用。使用它将总是能够避免三个连续的粒子变成共线结构，因此，任何扭转势都不会出现奇点。也可以将它加入到标准的弯曲势中，以处理约 180° 的角，但在极小点附近保持它的原始形式（参见图4.8中的橙色曲线）。

4.2.8 Urey-Bradley势

三个原子 i−j−k 之间的Urey-Bradley键角振动用角 θijk 的简谐势和原子 i 和 k 之间的距离的简谐校正项来描述。尽管可以很容易地将它写成两项的简单加和，但方便的做法是把它作为拓扑文件中的单独一项，并在在输出文件中作为单独的能量项。这种势能主要用在CHARMM力场[85]中。其能量由下式给出：

V a (θ i j k) = 1 2 k θ i j k (θ i j k - θ 0 i j k) 2 + 1 2 k U B i j k (r i k - r 0 i k) 2 (4.56)

力的方程可从4.2.1和4.2.5节推导出来.

4.2.9 键键交叉项

形成键 i−j 和 k−j 的三个粒子 i,j,k 之间的键键交叉项由下式给出[86]：

V r r' = k r r' (∣ r i - r j ∣ - r 1 e) (∣ r k - r j ∣ - r 2 e) (4.57)

其中 $k{rr’}是力常数，r{1e}和r_{2e}分别是键i - j和j - k的平衡键长。此势对应的粒子i$ 受到的力是：

F i = - k r r' (∣ r k - r j ∣ - r 2 e) r i - r j ∣ r i - r j ∣ (4.58)

原子 k 受到的力可以通过交换上述方程中 i 和 k 得到。最后，原子 j 受到的力遵循内力之和应该为零的事实： Fj=−Fi−Fk 。

4.2.10 键-角交叉项

为形成键 i−j 和 k−j 的三个粒子 i,j,k 之间的键-角交叉项由下式给出[86]：

V k θ = k r θ (∣ r i - r k ∣ - r 3 e) (∣ r i - r j ∣ - r 1 e + ∣ r k - r j ∣ - r 2 e) (4.59)

其中 $k{r\theta}为力常数，r{3e}为i - k之间的距离，其他常数的含义与方程4.57中相同。原子i$ 受到的力为：

F i = - k r θ [(∣ r i - r k ∣ - r 3 e) r i - r j ∣ r i - r j ∣ + (∣ r i - r j ∣ - r 1 e + ∣ r k - r j ∣ - r 2 e) r i - r k ∣ r i - r k ∣] (4.60)

4.2.11 四次键角势

键角势可以使用四阶多项式以满足特殊目的

V q (θ i j k) = \sum n = 0 5 C n (θ i j k - θ 0 i j k) n (4.61)

图4.9：异常二面角原理。环的面外弯曲（左），环的取代（中），四面体偏离（右）。在所有情况下，异常二面角 ξ 被定义为平面 (i,j,k) 和 (j,k,l) 之间的夹角。4.2.12 异常二面角

异常二面角是为了保持平面基团（例如芳环）处于平面内，或者为了防止分子翻转成其镜像分子，见图4.9。

异常二面角：简谐型

最简单的异常二面角势是简谐势，曲线如图4.10。

V i d (ξ i j k l) = 1 2 k ξ (ξ i j k l - ξ 0) 2 (4.62)

简谐势函数不连续，但由于不连续点选择在距离 ξ0 为 180° 的位置，这种不连续性永远不会引出麻烦。注意，在拓扑输入文件中，角度的单位为度，而力常数的单位为 kJ/mol/rad2 。

异常二面角：周期型

这种势函数等同于周期性的正常二面角（见下文）。将此种二面角类型（类型4）单独分类只是为了在参数部分和输出中将异常二面角和正常二面角区分开来。

图4.10：异常二面角的势能4.2.13 正常二面角

对于正常的二面角相互作用，可以选择GROMOS的周期性函数或基于 cosϕ 的幂的多项式展开函数（被称为Ryckaert-Bellemans势）。这种选择使得二面角四个原子中的第一和第四个原子之间的特殊相互作用也被包含在内。

对周期性的GROMOS势，必须包含特殊的1-4 LJ-相互作用；对烷烃的Ryckaert-Bellemans势，1-4相互作用必须从非键列表中排除。注意: Ryckaert-Bellemans势也可在OPLS等力场中与1-4相互作用联合使用。因此，在这种情况下你不应该修改pdb2gmx生成的拓扑文件。

正常二面角：周期型

正常二面角的定义遵从IUPAC/IUB约定，其中, ϕ 是面 ijk 和面 jkl 之间的夹角，零相应于顺式(cis) 构型(i 和 l 处于面的同一侧)。GROMACS拓扑文件中有两类二面角函数类型。标准类型1的行为和其他的键相互作用类似。对某些力场，类型9很有用。当在[dihedraltypes]部分对同一原子类型定义了多个参数时，类型9可使多个势能函数自动用于在[dihedral]部分的一个单一的二面角上。

这种类型的势函数为

V d (ϕ i j k l) = k ϕ (1 + cos (n ϕ - ϕ s)) (4.63)

正常二面角：Ryckaert-Bellemans势函数

对烷烃，经常使用以下的正常二面角势（见图4.12）：

V r b (ϕ i j k l) = \sum n = 0 5 C n (cos (ψ)) n (4.64)

其中， ψ=ϕ−180° 。

注意：从一个约定到另一个约定的转换可以通过将每个系数 Cn 乘上 −1n 。

常数C的一个示例在表4.1中给出。

图4.11：正常二面角的原理（左，反式trans构型）和二面角势（右）

表4.1：Ryckaert-Bellemans势的常数值(kJ mol−1)
C0 9.28	C2 -13.12	C4 26.24
C1 12.16	C3 -3.06	C5 -31.5

图4.12：Ryckaert-Bellemans二面角势

(注意: 使用这种势函数意味着排除了组成二面角的第一个和最后一个原子之间的LJ相互作用，并且 ψ 的定义遵从“聚合物约定”(ψtrans=0)。

RB二面角函数也可用于包含傅立叶二面角（见下文）：

V r b (ϕ i j k l) = 1 2 [F 1 (1 + cos (ϕ)) + F 2 (1 - cos (2 ϕ)) + F 3 (1 + cos (3 ϕ)) + F 4 (1 - cos (4 ϕ))] (4.65)

由于存在等式 cos(2ϕ)=2cos2(ϕ)−1,cos(3ϕ)=4cos3(ϕ)−3cos(ϕ) 和 cos(4ϕ)=8cos4(ϕ)−8cos2(ϕ)+1 ，可以利用下面的方程将OPLS参数变换为Ryckaert-Bellemans参数

C 0 C 1 C 2 C 3 C 4 C 5 = F 2 + 1 2 (F 1 + F 3) = 1 2 (- F 1 + 3 F 3) = - F 2 + 4 F 4 = - 2 F 3 = - 4 F 4 = 0 (4.66)

OPLS参数遵从蛋白质约定，而RB参数遵从聚合物约定（这导致 cos(ϕ) 奇数次项前面的负号）。

注意：记得将文献中OPLS和RB参数的单位 kcal mol−1 转化为GROMACS中的单位 kcal mol−1 。

正常二面角：傅立叶函数

OPLS势能函数为傅立叶级数的前三或四次余弦项[87]。GROMACS中，四项函数为：

V F (ϕ i j k l) = 1 2 [C 1 (1 + cos (ϕ)) + C 2 (1 - cos (2 ϕ)) + C 3 (1 + cos (3 ϕ)) + C 4 (1 + cos (4 ϕ))] (4.67)

GROMACS在内部使用Ryckaert-Bellemans代码来计算傅立叶二面角（见上文），因为计算效率更高。

注意：记得将文献中OPLS参数的单位 kcal mol−1 转化为GROMACS中的单位 kcal mol−1 。

正常二面角：限制扭转势

非常类似于限制弯曲势（见4.2.7），限制扭转/二面角势的函数为：

V ReT (ϕ i) = 1 2 k ϕ ( cos ϕ i - cos ϕ 0 ) 2 sin 2 ϕ i (4.68)

这类函数具有 cosϕ 函数的优点（对 sinϕ 求导时不存在问题），并且能保持扭转角处于仅有的一个最小值。在这种情况下，因子 sin2ϕ 不会让二面角从[ −180°:0 ]区间移到[ 0:180° ]区间内，即，函数不能同时在 −ϕ0 和 +ϕ0 处取极大值，而只能在其中之一取极大值。由于这个原因，初始构型中所有的二面角的角度值都必须处于所需的角度区间内，并且平衡的 ϕ0 值不应太靠近区间的端点（与4.2.7中所讲的限制弯曲势类似，推荐至少 10° 的差异）.

正常二面角：组合弯曲-扭转势

当形成二面角的四个粒子共线时（这种情况永远不会发生在原子尺度的模拟中，但可能出现在粗粒模拟中），扭转角及其势函数的计算会引起数值不稳定性。避免这一问题的一种方法是使用限制弯曲势（见4.2.7），以防止二面角达到 180° 。

另一种方式是忽略二面角不合理定义时的任何效应，通过将扭转势（余弦形式）与相邻弯曲角的弯曲势耦合为一个特殊形式，保持二面角的力和势函数的计算在整个角度范围内连续：

V CBT (θ i - 1, θ i, ϕ i) = k ϕ sin 3 θ i - 1 sin 3 θ i \sum n = 0 4 a n cos n ϕ i (4.69)

这种组合的弯曲扭转（combined bending-torsion，CBT）势函数由文献[88]提出，用于聚合物熔融的模拟。文献[84]对此有详细的描述。

这种势函数有两个主要的优点：

它不仅取决于(第 i−2,i−1,i 和 i+1 号珠子之间的)二面角 ϕi ，而且也取决于由三个相邻的珠子(分别是 i−2,i−1,i 和 i−1,i,i+1)形成的弯曲角 θi−1 和 θi 。当两个弯曲角中的任意一个接近 180° 时，由文献[89]尝试性地提出，并由文献[90]从理论上进行了讨论的两个 sin3θ 前因子抵消了扭转势和力。
通过表达为 cosϕi 的多项式形式，它也取决于 ϕi ，这避免了 ϕ=0° 或 ϕ=180° 时计算扭转力的奇异性。

这两种特性使得CBT势用于弯曲角存在弱约束的MD模拟时表现良好，甚至对于弯曲角和扭转角变化剧烈的拉伸/非平衡MD模拟也是如此。当使用CBT势时，相邻的 θi−1 和 θi 的弯曲势可以取任意形式，也可完全省略这两个角弯曲项(θi−1 和 θi)。图4.13显示了包含和不含 sin3θ 因子（分别为蓝色和灰色的曲线）的扭转势之间的区别。此外， VCBT 对笛卡尔变量的导数直接易懂：

\partial V CBT ( θ i - 1 , θ i ϕ i ) \partial r ⃗ l = \partial V CBT \partial θ i - 1 \partial θ i - 1 \partial r ⃗ l + \partial V CBT \partial θ i \partial θ i \partial r ⃗ l + \partial V CBT \partial ϕ i \partial ϕ i \partial r ⃗ l (4.70)

CBT基于无多重性的余弦形式，因此它只在 0° 左右对称。为获得一个不对称的二面角分布（如在[−180°:180°]区间内只有一个最大值），应该使用标准的扭转势，如简谐角势或周期性的余弦势来代替CBT势。然而，这两种形式对于力的微分（ 1/sinϕ ）和珠子的直线排列（ θi 或 θi−1=0°,180° ）都不方便。非 cosϕ 势与 sin3θ 因子的耦合并不能提高模拟的稳定性，因为在某些情况下 θ 和 ϕ 同时为 180° 。在此步的积分是可能的（由于扭转势的抵消），但下一步积分将是病态的（ θ 不是 180° 而 ϕ 非常接近 180° ）。

$图4.13：蓝色：组合弯曲扭转势的表面图(方程4.69，其中 $k=10 .text{kJ mol}^{-1},a_0 = 2.41, a_1 = -2.95, a_2 = 0.36, a_3 = 1.33$)，为简单起见，弯曲角相同($.theta_1=.theta_2=.theta$)。灰色：同样的扭转势函数但不含 $.sin^3.theta$ 项（Ryckaert-Bellemans型）。 $.phi$ 为二面角。$ 图4.13：蓝色：组合弯曲扭转势的表面图(方程4.69，其中 k=10kJ mol−1,a0=2.41,a1=−2.95,a2=0.36,a3=1.33)，为简单起见，弯曲角相同(θ1=θ2=θ)。灰色：同样的扭转势函数但不含 sin3θ 项（Ryckaert-Bellemans型）。 ϕ 为二面角。4.2.14 表格式键合相互作用

为了有充分的灵活性，可通过用户提供的表格函数形式，对键，键角和二面角使用任何函数。它们包括：

V b (r i j) V a (θ i j k) V d (ϕ i j k l) = k f b n (r i j) = k f a n (θ i j k) = k f d n (ϕ i j k l) (4.71) (4.72) (4.73)

其中, k 是以能量为单位的力常数, f 是三次样条函数；详细信息见6.10.1。对每种相互作用，在拓扑信息中指定力常数 k 和表格编号 n. 有两种不同类型的键，一种产生排除（类型8），一种不产生排除（类型9）。详情见表5.5。表格文件会被提供给mdrun程序。表格文件名和下划线后，加入字母“b”表示键，“a”表示角或“d”表示二面角和表格编号。例如，对一个 n=0 的键（并使用默认的表格文件名），将从文件table_b0.xvg中读取表格。如拓扑信息（表5.5）中具体介绍的那样，可简单地通过使用不同的 n 值提供多个表格，并应用到恰当的键。表格文件的格式是分别为 x,f(x),−f′(x) 的三列，其中 x 必须是等间隔的。拓扑信息的输入要求列于表5.5。这些表格的设置如下：

键: x 是以纳米(nm)为单位的距离。若模拟时距离超出表格长度，mdrun将显示错误消息并退出。

角度: x 是以度为单位的角度。表格应从0开始直到并包括180度；导数以度为单位。

二面角: x 是以度为单位的二面角。表格应从-180直到并包括180度，使用IUPAC/IUB约定，即零为顺式，导数以度为单位.

◆本文地址: http://jerkwin.github.io/9999/04/01/GROMACS中文手册第4章：相互作用函数和力场_1-2节, 转载请注明◆

转载本文请联系原作者获取授权，同时请注明本文来自李继存科学网博客。
链接地址：https://blog.sciencenet.cn/blog-548663-863366.html

上一篇：gnuplot动画
下一篇：轨道查看器OribitalViewer汉化版

收藏 IP: 130.184.197.*| 热度|

当前推荐数：0

该博文允许注册用户评论请点击登录评论 (0 个评论)

数据加载中...

返回顶部

李继存

扫一扫，分享此博文

全部作者的精选博文

• 1927索尔维会议(SolvayConference)合影

Jerkwin分享 http://blog.sciencenet.cn/u/Jerkwin

博文

ROMACS中文手册第4章：相互作用函数和力场 1-2节

当前推荐数：0

该博文允许注册用户评论请点击登录评论 (0 个评论)

李继存

全部作者的精选博文

全部作者的其他最新博文

全部精选博文导读

相关博文

Jerkwin分享 http://blog.sciencenet.cn/u/Jerkwin

博文

ROMACS中文手册第4章：相互作用函数和力场 1-2节

当前推荐数：0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

李继存

全部作者的精选博文

全部作者的其他最新博文

全部精选博文导读

相关博文

该博文允许注册用户评论请点击登录评论 (0 个评论)