【分子动力学教程】分子动力学模拟(Amber 20)

【分子动力学教程】分子动力学模拟(Amber 20)

用途

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

预备知识

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的重原子组成的集合

系综(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 + PV  ,在该系综下模拟需要保持压力和焓值固定,其调节技术实现起来有一定难度,在分子动力学模拟中并不常见(但LAMMPS支持)。

  • 巨正则系综(grand canonical ensemble,VTμ

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

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

本方案提供NVTNPTNVE三种系综,前两种最为常用,下面详细讨论。

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

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

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

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

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

不同分子的尺寸相差较大例如,蛋白质在溶液中的模拟,可先在NVT系综下限制蛋白的位置,让溶液弛豫以适应蛋白的位置,然后放开限制,转到NPT系综下进行预平衡

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

恒温器(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支持上述斜体标示的恒温器。

恒压器(barostat)

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

  • Berendsen方法

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

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

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

  • 蒙特卡洛(Monte Carlo,MC)方法

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

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

参数配置举例

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

  • 研究生物大分子与小分子结合模式,计算结合自由能(热力学性质)

  • 研究蛋白折叠运动(动力学性质)

  • 研究水溶液中有机小分子聚集行为,计算回转半径、径向分布函数等

  • 隐式溶剂环境(GB模型)中研究DNA-小分子的相互作用

入口

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

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

步骤

提交任务

1. 上传分子结构

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

坐标文件可以是:

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

2. 设置模拟环境

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

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

3. 设置计算步骤

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

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

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

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

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

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

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

4. 点击【提交】。

分析结果

热力学性质

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

  • 能量

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

  • 温度

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

  • 压强

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

  • 体积

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

  • 密度

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

RMSD

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

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

氢键

该处给出了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为该水桥出现的帧数。

分子结构

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

文件列表

该处给出分子动力学模拟的输出文件,其中*.rst文件(坐标文件)可作为后续模拟的起点(新建任务,上传到坐标文件处)。

参考资料

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

[2]https://sites.engineering.ucsb.edu/~shell/che210d/Advancedmoleculardynamics.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

殷赋科技创建学术交流群,请关注公众号“殷赋科技”,首页回复“加群”,验证后即可入群。