仓库源文，站点原文

layout: post title: GROMACS副本交换分子动力学 categories:

相对于标准的分子动力学模拟, 副本交换分子动力学(REMD)是一种增强采样的技术, 采用的方法是在不同温度下对具有相似势能的体系进行采样. 通过这种方法, 体系可能会越过势能面上的能垒, 从而探索新的构象空间.

实际上, REMD模拟的设置非常简单. 下面给出对指定体系运行REMD模拟的步骤. 至于模拟是否"成功"完全取决于要解决的问题和判断成功的标准. 虽然REMD模拟可以增加采样量, 但并不能提供最终答案. 这一点应该牢记在心.

一旦确定了多肽及其周围环境, 就需要确定待采样温度空间的范围, 使用的处理器数目以及模拟的时间长度. 体系, 副本数目, 温度空间的范围和温度分布决定了副本之间的平均交换概率. 这些值对所有副本都应该是相同的. 为此, 如果假定自由能空间中没有已知的瓶颈部分, 预期体系的势能会随温度的增加而增加, 因此副本的温度分布应服从指数分布

有充分的证据表明, 对多肽/蛋白质体系, 副本尝试交换的时间间隔不应小于1 ps. 这可以决定在两次交换之间, 一个副本探索相应构象空间所用的平均时间, 这个时间比实际的交换概率还重要. 在GROMACS的实现中, 特定一对副本的交换尝试会隔次进行, 因为奇数对和偶数对尝试交换会交替进行.

对低于4.0版本的GROMACS, 只允许每个副本使用一个处理器. 要运行REMD, 任何版本的GROMACS都必须使用MPI进行编译mdrun(即不能使用线程并行), 并且处理器的数目应该是副本数目的倍数.

一般步骤

  1. 定义体系, 例如: 多肽+溶剂(隐式或显式). GROMACS 4.5及以上版本才支持隐式溶剂.

  2. 根据可用的处理器数目和要采样的温度范围(实际上它们的相关性非常强), 选择温度分布. 使用指数分布: T_i=T_0 e^k*i, 其中k和T_0可以微调以获得合理的温度间隔保证能够进行交换. 指数形式可以保证温度间隔随温度升高而增大. 这是必要的, 因为总能量的分布随着温度的增加而增加, 因此交换概率也随之增加. 保持交换概率不随温度变化.

  3. 获得温度分布后, 在N个温度下分别对体系进行平衡(每个温度使用单独的.mdp文件). 然后根据平衡结构创建一系列(N个)运行输入文件(.tpr), 仍要使用不同的.mdp文件生成不同的.tpr文件. 如果有N个温度应该创建N个.tpr文件, 名称依次为prefix_0.tpr, ... prefix_N-1.tpr.
  4. 然后, 运行短时间的REMD模拟以获得交换概率的估计值(大约100 ps就可以获得良好的估计值), 如果得到的结果与所需值偏离较大, 试着修改温度. 典型情况下, 交换概率取0.2到0.3是合适的. 实际上, 比交换概率本身的值更重要的是副本在给定温度下的停留时间. 它的值为交换概率与交换频率的乘积. 例如, 交换概率为0.2, 每隔2 ps进行一次交换尝试, 可以得出相应温度下的平均时间为10 ps. 另外还需要考虑的是交换方式. 每次尝试时, 是尝试交换所有的副本对, 还是随机选择一对副本进行交换. 对后一种情况, 还应考虑一点, 必须将副本的停留时间乘以N-1, 即交换对的数目.
  5. 每个温度下的初始构象可以相同, 也可以不同. 如何选择取决于运行REMD的原因.

网上有一个可以选择T-REMD温度的在线工具, 它是基于溶剂化蛋白已知的能量分布计算的. 输入蛋白质的原子数, 水的原子数, 温度范围和所需的交换概率, 这个工具就会生成所需的温度.

具体步骤

  1. 创建一组.mdp文件, 每个文件指定要使用的不同温度. 根据索引对得到的.tpr文件进行编号, 从0到任意数字(例如, 如果有10个不同的温度, 则名称为prefix_0.tpr, prefix_1.tpr ... prefix_9.tpr). 对GROMACS 4.0之前的版本, 每个副本只能使用一个处理器, 因此要么省略gmx grompp的-np选项, 要么使用-np 1. 对于GROMACS 4.0, gmx grompp时不需要使用-np选项.
  2. 副本交换的间隔由mdrun选项-replex指定. 计算所用的核心数必须是副本数的倍数(由-multi给出, 必须等于.tpr文件的数目, 对于上面的示例, 使用prefix_0.tpr至prefix_9.tpr, 此数目为10). 如上所述, 对于4.0版之前的GROMACS, 该倍数必须为1. 输入文件的命名对于mdrun的正确运行至关重要.
  3. 命令行示例参见下文. 运行步数 是要达到所需的交换尝试周期所需的步数. 可能还需要其他命令行选项, 例如-o prefix

GROMACS 3.x: mpirun -np 10 mdrun -s prefix.tpr -np 10 -multi 10 -replex (步数) (后接输出选项) GROMACS 4.x: mpirun mdrun -s prefix.tpr -multi 10 -replex (步数) (后跟输出选项)

理解REMD相关的输出

mdrun会将重要信息输出到md.log文件, 你可以像这样提取它:

% grep Repl md0.log gives:

Initializing Replica Exchange
Repl  There are 6 replicas:
Repl      0     1     2     3     4     5
Repl  T 300.0 350.0 410.0 480.0 560.0 650.0
Repl
Repl  exchange interval: 1000
Repl  random seed: 525106
Repl  below: x=exchange, pr=probability
Replica exchange at step 1000 time 2
Repl 0 <-> 1  dE =  2.492e+00
Repl ex  0    1    2    3    4    5
Repl pr   .08       .00       .13
Replica exchange at step 2000 time 4
Repl ex  0    1    2    3 x  4    5
Repl pr        .00       1.0

这些解释确实很简略, 希望它能够提供足够的信息来帮助你理解得到的结果.

后处理

GROMACS输出的REMD轨迹是系综连续的, 但相对于模拟时间却不连续. 如果需要后者, 可以使用脚本scripts/demux.pl读取md0.log文件(必要时可以合并多个文件), 生成一些输出文件. 其中之一是.xvg文件(replica_ndx.xvg), trjcat可以使用该文件和原始轨迹文件得到连续轨迹. 另一个文件(replica_temp.xvg)包含了每个副本从原始温度开始的温度. 因此, 如果感兴趣的副本开始于, 比方说, 300 K, 你可以在温度空间中跟踪它的轨迹. 如果能给出每个副本的温度分布直方图可能更好, 根据大多数作者的说法, 这种直方图应该是平坦的. 使用这种方法处理过的轨迹已经用在论文中了, 可以根据REMD轨迹获得蛋白质的折叠动力学Rev. Lett. 96, 238102 (2006).(在g_kinetics中实现-在GROMACS 3.3.1或更早版本中不可用)