分子动力学模拟（Amber20）

分子动力学
分子动力学方案
发布日期:   2023-07-03
更新日期:   2023-09-14
文章字数:   4.4k
阅读时长:   16 分
阅读次数:  

1. 用途

使用Amber 20程序执行有机生物体系的常规平衡态分子动力学模拟，可采用显式溶剂模型和隐式溶剂模型（GB模型）。

2. 预备知识

2.1 Amber mask

Amber mask是Amber所用的选择原子或残基的记号，常用于能量最小化或分子动力学中约束原子或残基。

Amber mask表达式用法如下：

  • 以英文冒号:开始，后接残基编号或名称来指定残基。数字可用逗号分隔，也可用短横线指定范围，名称只能用逗号分隔。

    :1-10 表示残基1到10

    :1,3,5 表示残基1、3和5

    :1-3,5,7-9 表示残基1到3、5和7到9

    :LYS 表示所有赖氨酸残基

    :ARG,ALA,GLY 表示所有精氨酸、丙氨酸和甘氨酸

  • 以@开始，后接原子序号或名称来指定原子，使用原则与残基一致。

    @12,54-85,90 表示原子12、54到85和90

    @CA,C,N 表示所有名为CA、C和N的原子（即蛋白骨架）

  • 通配符用等号=表示（而非星号*，因为*在Amber的原子命名里有特殊含义）。

    @H= 表示所有氢原子

  • 逻辑组合：与&、或|、非!。必要时，用圆括号把不同表达式分隔开。

    :1-234&!@H= 表示残基1到234中所有重原子（非氢原子）

    :1-234@CA,C,N 等价于 :1-234&@CA,C,N，表示残基1到234中的骨架原子（CA、C、N）

    (:1-234@CA,C,N)|(:235&!@H=) 表示残基1到234中的骨架原子和残基235的重原子组成的集合

2.2 系综（ensemble）

在一定的宏观条件下，大量性质和结构完全相同的、处于各种运动状态的、各自独立的系统的集合，称为统计系综。

根据宏观约束条件，系综可分为：

  • 微正则系综（micro-canonical ensemble，NVE）

    具有确定的粒子数（N）、体积（V）和总能量（E）。平衡体系为孤立系统，与外界既无能量交换，也无粒子交换。通过调整原子的速度来调整能量，但有可能使系统失去平衡，可通过迭代迟豫来达到平衡。

  • 正则系综（canonical ensemble，NVT）

    全称宏观正则系综，具有确定的粒子数（N）、体积（V）和温度（T）。平衡体系为封闭系统，与大热源热接触平衡的恒温系统。通过调整原子的速度来保持系统动能恒定。

  • 等温等压（constant-pressure，constant-temperature，NPT）

    具有确定的粒子数（N）、压强（P）和温度（T）。通过调整系统的体积来保持压强恒定。

  • 等压等焓（constant-pressure，constant-enthalpy，NPH）

    具有确定的粒子数（N）、压强（P）和焓（H）。由于H=E+PVH=E+PV，在该系综下模拟需要保持压力和焓值固定，其调节技术实现起来有一定难度，在分子动力学模拟中并不常见（但LAMMPS支持）。

  • 巨正则系综（grand canonical ensemble，VTμ\mu）

    具有确定的体积（N）、温度（P）和化学势（μ\mu）。体系是开放系统，与大热源大粒子源热接触平衡而具有恒定的温度。

那在动力学模拟中如何选择系综？

本方案仅提供NVT和NPT两种系综，下面只谈论它们的使用。

常规体系的模拟是在等温等压（NPT）或等温等体积（NVT）条件下进行的。NPT是最接近常规生物实验条件的，但NPT的计算性能不如NVT，因此，NVT更有利于长时间模拟。另外，NPT会改变原胞的形状（调整体积以使压强恒定），因此，适用于形状变化较大的情况（如，打开又折叠的蛋白体系），而形状变化不大的，使用NVT更为合适。

由于大多数情况是体系形状不会有太大改变的，因此，常规思路是：先在NPT系综下进行预平衡，让体系迟豫到稳态，然后再转到NVT系综进行长时间动力学采样。

对于发生构象折叠的蛋白体系，可采用这样的思路：先在NPT系综下进行预平衡，让体系达到预期条件，然后再转到NVT系综进行长时间动力学采样。

下面列举一些特殊体系的处理方法：

  • 两相体系

    在构建处是构型前，需要对两相分别进行平衡处理，再组成新的结构。一种做法是：再NVT系综下限制固相位置，让液体充分适应固相位置。然后放开限制，转到NPT系综下进行预平衡。

  • 不同分子的尺寸相差较大

    例如，蛋白质在溶液中的模拟，可先在NVT系综下限制蛋白的位置，让溶液迟豫以适应蛋白的位置，然后放开限制，转到NPT系综下进行预平衡。

  • 平衡时间长的体系

    例如，离子液体。该类体系所需平衡时间较长，可考虑在NVT下高温跑一段时间，使体系空间分布变得均匀，然后退火到目标温度，再转到NPT系综下进行预平衡。

2.3 恒温器（thermostat）

在分子模拟中，一般需要通过恒温器使系统温度维持在给定值附近。系统温度和粒子的速度直接相关，通过调整粒子的速度可使系统温度维持在目标值。下面列举常用的恒温器：

  • 速度标度（scaling velocities）类方法

    这类方法不严格遵循正则系综，虽然在实践中，表现出来的偏离程度不是很大。因此，适用于加热阶段，理论上不太建议用于平衡采样，实际上看情况使用（使用Berendsen的文献仍在增长）。

    • 简单速度标度（Simple velocity scaling）

    方法简单，但速度不按照波尔兹曼分布，无法对应任何一个统计力学的系综。容易引起系统能量突变，使系统和真实的平衡态相距甚远。一般不采用。

    • Berendsen热浴和Bussi热浴

    弱耦合热浴。基本思想是将系统和一个恒温的外部热浴耦合，通过热浴吸收或释放能量来调节系统的温度。

    在系统远离平衡态时，对温度的调节效率较高，但动能不严格遵循波尔兹曼分布，可能产生“飞冰块现象”（Flying Ice Cube Effect）。因此，**它适用于加热阶段，不适用于预平衡阶段。**但当体系已经充分平衡，需要校正的程度很小，同时使用较弱的耦合常数（比如，10 ps），使总能量震荡较小，可用于生产（采样）阶段。

    Bussi热浴是Berendsen热浴的随机版，它将速度调整为从正则分布探得的随机温度，从而对正则分布进行采样。

    另外，Berendsen热浴不适合用于隐式溶剂模拟，因为无法通过与溶剂的碰撞来帮助维持恒温。

    Flying Ice Cube Effect：当分子振动越来越弱，而平动或转动越来越强的时候，分子像冻僵的方块飞来飞去的现象。造成这一现象的原因是这类热浴算法使分子的振动能被“砍掉”（标度）得多、补充得少，逐渐“冻僵”，而平动能和转动能受影响较小。

  • 随机力或随机速度

    这类方法的优点是能够从正则系综中正确采样，可以放心使用比NVE更长的时步。缺点是因引入了随机性而变得非确定性、非时间反演，破坏了动量传递，无法用于考擦系统的动力学性质（如扩散系数）。

    • Langevin热浴和Andersen热浴

    这两种方法的速度和精确度中等。两者采用虚拟随机碰撞的方式来调节速度，干扰了系统的正常演化，削弱了粒子间的速度相关性。但严格遵循正则系综，不影响各态遍历性（ergodicity）。因此，不适合计算动力学性质，可计算热力学性质（如，结合自由能）。

    在隐式溶剂环境下，Langevin热浴正好通过虚拟碰撞补回粘性效应，因此适合使用。

  • 扩展拉格朗日方法

    这类方法克服了上述方法的缺点，可时间反演。代表方法是Nosé-Hoover热浴。

    • Nosé-Hoover热浴

    该算法较为复杂，计算速度较慢。在系统远离平衡态时，温度振荡较大，不易收敛，不宜用于加热阶段。但它严格产生正则系综热力学，并近似获得真实的动力学，最适合用于平衡采样。

    值得一提的是，Nosé-Hoover热浴在特定体系中会表现出病态行为，近年来发展了一些改进方法，例如：Optimized Isokinetic Nose-Hoover chain（OIN）和Stochastic Isokinetic Nose-Hoover RESPA integrator（详见Amber 20的手册第346页及相关文献）。

Amber支持上述斜体标示的恒温器。

2.4 恒压器（barostat）

在NPT系综里，需要采用恒压器调控体系压力。常用的技术有：

  • Berendsen方法

    在系统远离平衡态时，该方法对压力的调节效率较高，适用于最初的压力迟豫，但它不按照正则分布来采样，一般不适用于平衡采样。

  • Nosé-Hoover方法和Parrinello-Rahman方法

    与Berendsen方法相反，适用于平衡态控压。

  • 蒙特卡洛（Monte Carlo，MC）方法

    方法简单，但效率相对较低。该方法适用于采样阶段。

Amber支持上述斜体标示的恒压器。

2.5 参数配置举例

分子动力学模拟是一项艺术。根据不同目的、针对不同情况，不同人有不同做法。但只要设置合理、达到目的，即为成功的模拟。下面列举一些例子以供参考：

  • 研究生物大分子与小分子结合模式，计算结合自由能（热力学性质）
阶段 系综 恒温器 恒压器 时长 备注
加热（Heat） NVT Langevin/Berendsen Berendsen ~20 ps 不宜过长，否则容易产生真空泡
平衡（Equilibration） NPT Langevin Berendsen 100~ 500 ps 使体积（密度）达到平衡
平衡（Equilibration） NVT Langevin Berendsen 1~5 ns 使能量、RMSD平稳
生产（Production） NVT Langevin/Berendsen/Nosé-Hoover Monte Carlo 按需
  • 研究蛋白折叠运动（动力学性质）
阶段 系综 恒温器 恒压器 时长
加热（Heat） NVT Langevin/Berendsen Berendsen ~20 ps
平衡（Equilibration） NPT Langevin Berendsen 100~ 500 ps
平衡（Equilibration） NVT Langevin Berendsen 1~5 ns
生产（Production） NPT Berendsen/Nosé-Hoover Monte Carlo 按需
  • 研究水溶液中有机小分子聚集行为，计算回转半径、径向分布函数等
阶段 系综 恒温器 恒压器 时长
加热（Heat） NVT Langevin/Berendsen Berendsen ~20 ps
平衡（Equilibration） NPT Langevin Berendsen 100~ 500 ps
平衡（Equilibration） NVT Langevin Berendsen 1~5 ns
生产（Production） NPT Berendsen/Nosé-Hoover Monte Carlo 按需
  • 隐式溶剂环境（GB模型）中研究DNA-小分子的相互作用
阶段 恒温器 时长
加热（Heat） Langevin 100 ps
平衡（Equilibration） Langevin 500 ps
生产（Production） Langevin 按需

3. 入口

平台地址：https://cloud.yinfotek.com/

功能入口：左侧菜单栏【计算方案】->【大方案】->【分子动力学】->【分子动力学（Amber 20）】

4. 平台操作步骤

4.1 提交任务

  1. 上传分子结构

    上传Amber的拓扑文件和坐标文件。

    坐标文件可以是：

    • 由【准备Amber文件】生成的inpcrd文件；
    • 上一次计算产生的rst文件；
    • 轨迹分析生成的pdb文件（只要原子坐标数目与拓扑文件对应）。

  2. 设置模拟环境

    显式溶剂的使用涉及温度（在计算步骤中控制）和压强（恒压）的控制，还需要设置非键截断值。

    注意，非键截断值绝对不能超过盒子边长！通常设为与水盒子大小一样即可（使用【准备Amber文件】工具添加水盒子时，盒子大小默认为10.0 Å）。

  3. 设置计算步骤

    设置计算步骤，选择（能量最小化）迭代步数或（分子动力学）模拟时长，选择系综（显式溶剂环境）、恒温器和设置目标温度。

    实际操作中，建议将计算步骤分成不同任务来计算，根据计算结果调整下一步设置。

    例如：先计算3步Minimization，看看最终能量是否不再下降而趋于平稳，否则可增设Minimization；再计算Heat和Equilibration，看看能量、RMSD是否平稳，否则可增设Equilibration；确认体系达到平衡状态后，再计算Production。具体例子见本公众号后续案例。

    这种做法还有个好处是，可以根据Equilibration所花时间和费用，决定Production的时长。

    此外，模拟过程中某些阶段通常还需要对体系施加一定约束，以免出现意外情况、偏离预期状态。点击约束条件下的【添加】，填写Amber mask（详见预备知识）和选择约束力，点击【确定】。

    如下图所示，在Minimization步骤，对1-267号残基（在本例中，1-266为蛋白质，267为小分子）施加10 kcal·mol−1^{-1}·Å−2^{-2}的约束力。

    注：目前仅支持对原子设置位置约束，将来会增加距离、角度、二面角度等约束方式。

  4. 点击【提交】。

4.2 分析结果

4.2.1 热力学性质

每一个阶段均给出热力学数据，包括：能量（总能量、动能、势能）、温度、压强、体积和密度。分析这些性质，判断模拟是否达到预期效果。

  • 能量

    在能量最小化（Minimization）阶段，能量一般会先下降后平稳；在加热（Heat）阶段，能量应在极短时间内达到平稳（范德华能量平稳下降，其他能量平稳上升）；在平衡（Equilibration）或生产（Production）阶段，能量应保持平稳。

  • 温度

    Minimization阶段无温度；在Heat阶段，温度应从初始温度几乎沿着直线上升（或下降）至目标温度；而当模拟进入Equilibration或Production阶段，温度应在目标温度附近波动（下图红色移动平均线显示平均值在300 K附近波动）。

  • 压强

    只有（显式溶剂环境）NPT系综才有压强，其曲线应在极短时间里达到期望值并保持平稳（下图显示，红色移动平均线在1 atm附近波动）。

  • 体积

    体积与压强负相关，只有（显式溶剂环境）NPT系综才有，其曲线与压强大致相反。

  • 密度

    密度与体积是倒数关系，也是只有（显式溶剂环境）NPT系综才有，其曲线与体积刚好相反。

4.2.2 RMSD

该处给出了Equilibration和Production阶段的RMSD数据。通过RMSD分析，了解体系的构象变化，是否达到平衡。

很多后续分析工作都建立在平衡状态的基础上，例如：进一步采样（production）、氢键分析、聚类分析、结合自由能分析；等等。

4.2.3 氢键

该处给出了Production阶段的氢键分析数据文件，将来会提供更多数据可视化功能。

  • 仅含有机小分子的体系才有氢键分析数据，纯生物大分子体系无此数据。若有分析需要，请使用专门的氢键分析方案。

数据文件hbond.xlsx中提供了5项内容，分别是：总体（Total）统计数据、溶质-溶质（solute-solute）氢键统计数据、溶质-溶剂（solute-solvent）氢键统计数据、溶质-溶质（solute-solute）氢键时间序列、溶质-溶剂（solute-solvent）氢键时间序列以及水桥统计数据。

  • 总体统计数据

该表给出各帧（Frame）的溶质-溶质（Solute-solute）氢键、溶质-溶剂（Solute-solvent）氢键和水桥（Bridge）个数。

  • 氢键统计数据

该表给出各个氢键的统计数据：Name为形成氢键的残基对（如，Glu191的OE1原子与Ser132的OG原子和HG原子形成氢键），Frames为轨迹中该氢键存在的帧数，Fraction为比例（等于Frames除以总采样帧数），Avg. Distance和Avg. Angle分别为平均键长和键角，Lifetimes为生命周期数，Max. Lifetime和Avg. Lifetime分别为最长寿命和平均寿命。

OE1表示名为OE1的氧原子，同理，HG为氢原子。

  • 氢键时间序列

时间序列给出了各氢键随时间（帧号）的生存情况。如下图所示，Frame为帧号，每个氢键的残基对都自成一列，其中1表示存在，0表示不存在。

当溶质和溶剂形成氢键时，溶剂分子用solvent表示。

  • 水桥统计数据

如下图所示，Bridged Residues为形成水桥的溶质残基（例如，Cyx95、Trp114和DE3126之间存在水介导的氢键），Frames为该水桥出现的帧数。

4.2.4 分子结构

该处给出各个阶段最后一帧的结构，其中，文件名带有-dry的是删除水和离子后的结构。

4.2.5 文件列表

该处给出分子动力学模拟的输出文件，其中rst文件（坐标文件）可作为后续模拟的起点。

5 参考资料

[1] http://www.pages.drexel.edu/~cfa22/msim/node32.html

[2] https://sites.engineering.ucsb.edu/~shell/che210d/Advanced_molecular_dynamics.pdf

[3] Efrem Braun, Seyed Mohamad Moosavi, and Berend Smit. Anomalous Effects of Velocity Rescaling Algorithms: The Flying Ice Cube Effect Revisited. Journal of Chemical Theory and Computation 2018 14 (10), 5262-5272. DOI: 10.1021/acs.jctc.8b00446

文章作者: 殷赋量子氢
文章链接: http://blog.yinfotek.com/posts/10401/
版权声明: 本博客所有文章除特別声明外，均采用 CC BY 4.0 许可协议。转载请注明来源 殷赋量子氢 !
分子动力学
  目录