仓库源文站点原文


layout: post title: GROMACS质心牵引的几点说明 categories:


GROMACS的质心牵引功能也称为牵引代码(或拉伸代码, Pull Code), 可用于伞形采样, 并计算平均力势PMF, 在分子动力学模拟特别是涉及生物分子的模拟中经常用到. GROMACS手册中涉及到这个功能的说明主要有两处:

  1. 6.4 牵引代码
  2. 7.3.21 质心牵引

在使用本功能前, 请先仔细阅读并试着理解手册中的说明. 下面是我参考Sobereva的博文Gromacs的pull code介绍(3.3.x+4.0.x)并根据自己的一些经验整理的几点说明, 希望能帮助大家更好地理解并使用这个功能.

在进行质心牵引时, 首先要明确使用哪种牵引方式, 然后再明确沿什么方向进行牵引, 最后再是其他的一些设置, 如牵引速率, 力常数等. 这些牵引参数都可以直接在.mdp文件中进行定义. 下面对其中重要的几个关键词进行说明.

牵引方式: pull

牵引方式的设置使用pull关键词, 有下面几个选项.

伞形牵引的详细说明

使用伞形牵引时, 会在牵引点(牵引组质心)和弹簧点(用于定义牵引距离的虚拟位置)之间施加一个简谐势以维持牵引组相对于参考组的位置, 相当于牵引点与弹簧点之间存在一个虚拟的弹簧, 弹簧的一端处于牵引点, 另一端处于弹簧点. 牵引过程中, 弹簧点会匀速远离, 这样就会给牵引组原子一个外力, 拉动牵引组原子. 弹簧点的远离速度由pull-coord*-rate设定, 弹簧的弹性系数(简谐势的力常数)由pull-coord*-k设定. 因此, 牵引组原子所受的外力取决于弹簧点与牵引点之间的距离和弹簧的弹性常数.

举例来说吧, 对下面的设置

A(参考点)-------B(牵引点)====C(弹簧点)

A为参考组, B为牵引组, C为弹簧点.

t=0时刻AC间的距离(对应于C的初始位置)由pull-coord*-init设定. 牵引开始后AC间的距离会以pull-coord*-rate设定的速率均匀增加, C的移动会增大BC之间的距离, 产生弹力进而拉动B. B所受力的大小等于BC的距离乘以由pull-coord*-k设定的力常数.

从初始位置开始, C的移动方向有牵引几何设置决定: 当pull-geometry=distance时, C始终沿AB方向移动; 当pull-geometry=direction时, C沿由pull-coord*-vec设定的方向移动.

牵引时, 如果没有定义参考组, 牵引组就处于绝对坐标中, 参考点被设为(0 0 0). 如果参考点不是固定坐标而是一组可移动的原子, 牵引时, 如果参考组发生了移动, 牵引组也会随之移动, 将参考点的位置变化累加到弹簧点上, 以保持AC间距均匀增加.

此外, 弹簧点的移动速率可设置为0, 这样相当于在BC之间添加了一个固定的弹簧, 可用于伞形采样.

牵引几何: pull-geometry

牵引几何主要用于控制牵引的方向, 使用pull-geometry关键词设定. 有下面几个选项.

总的来说, 使用pull-geometry=distance时, 牵引方向是牵引组质心与参考组质心连线的瞬时方向; 使用pull-geometry=direction时则是直接定义牵引组的拉伸方向; pull-geometry=cylinder则适合用于生物膜体系.

牵引方式与牵引几何的组合

使用pull=umbrella

上面几种情况都是假定参考组是绝对坐标或者被固定住, 如果参考组在牵引过程中位置会发生变化, 弹簧点位置也会相应变化, 以维持间距均匀增长.

使用pull=constraint

pull=umbrella一样, 区别在于使用约束算法将牵引组位置固定在弹簧点位置上, 而不会被弹簧点拉着走.

比较常用的是pull-geometry=distance, 可使参考组与牵引组之间的距离保持固定不变(即pull-coord*-init)或者刻意让其逐渐改变(即pull-coord*-init+time`pull-coord-rate`), 而它们之间的相对角度方位可随意变化.

使用pull=constant-force

注意, 这种情况下拉力大小始终不变, 由pull-coord*-k设定. 无论使用哪种几何都没有弹簧点的概念, 因此pull-coord*-rate, pull-coord*-init对结果无影响.

总结

上述几种组合中, 假定参考组位置固定, 若不进行位置固定, 则参考组与牵引组的运动是对称的, 质心位置保持不变(包括pull=constant-force+pull-geometry=direction这种不直接涉及到参考组的模式也是如此), 两个组相当于等价.

常用的组合是pull=constraint+pull-geometry=distance可约束牵引组与参考组的距离; pull=constant-force+pull-geometry=direction方便地直接将牵引组往某个方向拉.

其他参数

输出结果

mdrun时可利用-pf选项指定输出xvg文件的名称, 其中包含时间和对应时刻牵引组的受力(外加的牵引力, 而非牵引组的净受力); -px指定的xvg输出文件中包含时间, 坐标和对应时刻的受力, 其中第一列为时间, 后面几列为参考组的绝对坐标, 牵引组相对于参考组的坐标. 只有pull-dim设为Y的维度才输出. 这些量的输出频率由pull-nstfoutpull-nstxout选项控制.

输出文件中的单位默认是GROMACS的标准单位, 时间ps, 位置nm, 力kJ/mol-nm.

如果运行不久就出现一堆.pdb文件, 提示write pdb之类, 说明体系崩溃了, 可能因为牵引速度过快, 或力常数过大, 可以试着把它们改小些. 可利用VMD等可视化软件查看牵引轨迹, 检查牵引过程是否合理, 比如牵引过程中是否出现过度碰撞, 挤压等.

不同版本关键词的差异

上文中所用的关键词都是基于GROMACS 5.x, 它们与GROMACS 4.x的存在一定差异.

<table id='tab-0'><caption>GROMACS不同版本中牵关键词的比较</caption> <tr> <th rowspan="1" colspan="1" style="text-align:center;">GMX 4.x</th> <th rowspan="1" colspan="1" style="text-align:center;">GMX 5.x</th> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">pull-geometry = position</td> <td rowspan="1" colspan="1" style="text-align:center;">废弃</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">无</td> <td rowspan="1" colspan="1" style="text-align:center;">pull-print-reference</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">无</td> <td rowspan="1" colspan="1" style="text-align:center;">pull-ncoords</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">pull-group0</td> <td rowspan="1" colspan="1" style="text-align:center;">废弃</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">pull-weights0</td> <td rowspan="1" colspan="1" style="text-align:center;">废弃</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">pull-pbcatom0</td> <td rowspan="1" colspan="1" style="text-align:center;">废弃</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">pull-group1</td> <td rowspan="1" colspan="1" style="text-align:center;">pull-group1-name</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">pull-weights1</td> <td rowspan="1" colspan="1" style="text-align:center;">pull-group1-weights</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">pull-pbcatom1</td> <td rowspan="1" colspan="1" style="text-align:center;">pull-group1-pbcatom</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">无</td> <td rowspan="1" colspan="1" style="text-align:center;">pull-coord1-groups</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">无</td> <td rowspan="1" colspan="1" style="text-align:center;">pull-coord1-origin</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">pull-vec1</td> <td rowspan="1" colspan="1" style="text-align:center;">pull-coord1-vec</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">pull-init1</td> <td rowspan="1" colspan="1" style="text-align:center;">pull-coord1-init</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">pull-rate1</td> <td rowspan="1" colspan="1" style="text-align:center;">pull-coord1-rate</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">pull-k1</td> <td rowspan="1" colspan="1" style="text-align:center;">pull-coord1-k</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">pull-kB1</td> <td rowspan="1" colspan="1" style="text-align:center;">pull-coord1-kB</td> </tr> </table>

示例说明

下面是一个运行示例, 将甲烷分子从水层的一侧牵引至另一侧. 运行文件基于GROMACS 4.x版本, 但稍加修改就应该可以用于GROMACS 5.x版本. 详细建模过程就不细说了, 大致如下:

  1. 创建1x1x1 nm^3^的水盒子
  2. 将盒子z方向长度改为4 nm, 并使体系在盒子内居中
  3. 将甲烷置于盒子z方向中轴线上, 其中心距盒子中心1 nm
  4. 设定运行参数和牵引参数, 运行1 ps, 沿z方向牵引, 牵引速度 2 nm/ps, 力常数 1000 kJ/mol-nm^2^
  5. 运行模拟

注意, 使用的各种数据和参数, 在用于实际体系时要根据情况做相应的更改.

模拟轨迹如下

看见在上面的模拟条件下, 1 ps时间内甲烷分子没有完全穿过水层, 需要增加模拟时间或是增大弹簧的力常数.

分析下输出文件pullx.xvgpullf.xvg. 初始时牵引点和弹簧点同时处于z~0~=1 nm处, 牵引开始后弹簧点以v=2 nm/ps的速度匀速远离, 同时带动牵引点远离. pullx.xvg中给出了牵引点每个时刻的位置, 再结合知道的力常数k=1000 kJ/mol-nm^2^, 我们就可以计算出每个时间的牵引力

$$F(t)=(z_0+vt-z)k$$

与文件中给出的数据完全一致

评论