用途 使用Amber 20程序执行有机生物体系的常规平衡态分子动力学模拟,可采用显式溶剂模型和隐式溶剂模型(GB模型)。 预备知识Amber maskAmber mask是Amber所用的选择原子或残基的记号,常用于能量最小化或分子动力学中约束原子或残基。 Amber mask表达式用法如下: :1-10 表示残基1到10 :1,3,5 表示残基1、3和5 :1-3,5,7-9 表示残基1到3、5和7到9 YS 表示所有赖氨酸残基 :ARG,ALA,GLY 表示所有精氨酸、丙氨酸和甘氨酸 @12,54-85,90 表示原子12、54到85和90 @CA,C,N 表示所有名为CA、C和N的原子(即蛋白骨架)
@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)和化学势( )。体系是开放系统,与大热源大粒子源热接触平衡而具有恒定的温度。
那在动力学模拟中如何选择系综? 本方案提供NVT、NPT和NVE三种系综,前两种最为常用,下面详细讨论。
常规体系的模拟是在等温等压(NPT)或等温等体积(NVT)条件下进行的。NPT是最接近常规生物实验条件的,但NPT的计算性能不如NVT,因此,NVT更有利于长时间模拟。另外,NPT会改变原胞的形状(调整体积以使压强恒定),因此,适用于形状变化较大的情况(如,打开又折叠的蛋白体系),而形状变化不大的,使用NVT更为合适。 由于大多数情况是体系形状不会有太大改变的,因此,常规思路是:先在NPT系综下进行预平衡,让体系弛豫到稳态,然后再转到NVT系综进行长时间动力学采样。 对于发生构象折叠的蛋白体系,可采用这样的思路:先在NPT系综下进行预平衡,让体系达到预期条件,然后再转到NVT系综进行长时间动力学采样。 下面列举一些特殊体系的处理方法: 两相体系 在构建系统构型前,需要对两相分别进行平衡处理,再组成新的结构。一种做法是:在NVT系综下限制固相位置,让液体充分适应固相位置。然后放开限制,转到NPT系综下进行预平衡。 不同分子的尺寸相差较大 例如,蛋白质在溶液中的模拟,可先在NVT系综下限制蛋白的位置,让溶液弛豫以适应蛋白的位置,然后放开限制,转到NPT系综下进行预平衡。 平衡时间长的体系 例如,离子液体。该类体系所需平衡时间较长,可考虑在NVT下高温跑一段时间,使体系空间分布变得均匀,然后退火到目标温度,再转到NPT系综下进行预平衡。
恒温器(thermostat)在分子模拟中,一般需要通过恒温器使系统温度维持在给定值附近。系统温度和粒子的速度直接相关,通过调整粒子的速度可使系统温度维持在目标值。下面列举常用的恒温器: 简单速度标度(Simple velocity scaling) 方法简单,但速度不按照波尔兹曼分布,无法对应任何一个统计力学的系综。容易引起系统能量突变,使系统和真实的平衡态相距甚远。一般不采用。 Berendsen热浴和Bussi热浴 弱耦合热浴。基本思想是将系统和一个恒温的外部热浴耦合,通过热浴吸收或释放能量来调节系统的温度。 在系统远离平衡态时,对温度的调节效率较高,但动能不严格遵循波尔兹曼分布,可能产生“飞冰块现象”(Flying Ice Cube Effect)。因此,它适用于加热阶段,不适用于预平衡阶段。但当体系已经充分平衡,需要校正的程度很小,同时使用较弱的耦合常数(比如,10 ps),使总能量震荡较小,可用于生产(采样)阶段。 Bussi热浴是Berendsen热浴的随机版,它根据正则分布随机调整温度。
另外,Berendsen热浴不适合用于隐式溶剂模拟,因为无法通过与溶剂的碰撞来帮助维持恒温。
Flying Ice Cube Effect:当分子振动越来越弱,而平动或转动越来越强的时候,分子像冻僵的方块飞来飞去的现象。造成这一现象的原因是这类热浴算法使分子的振动能被“砍掉”(标度)得多、补充得少,逐渐“冻僵”,而平动能和转动能受影响较小。
Nosé-Hoover热浴 该算法较为复杂,计算速度较慢。在系统远离平衡态时,温度振荡较大,不易收敛,不宜用于加热阶段。但它严格产生正则系综热力学,并近似获得真实的动力学,最适合用于平衡采样。 值得一提的是,Nosé-Hoover热浴在特定体系中会表现出病态行为,近年来发展了一些改进方法,例如:Optimized Isokinetic Nose-Hoover chain(OIN)和Stochastic Isokinetic Nose-Hoover RESPA integrator(详见Amber 20的手册第346页及相关文献)。
Amber支持上述斜体标示的恒温器。 恒压器(barostat)在NPT系综里,需要采用恒压器调控体系压力。常用的技术有: Amber支持上述斜体标示的恒压器。 参数配置举例分子动力学模拟是一项艺术。根据不同目的、针对不同情况,不同人有不同做法。但只要设置合理、达到目的,即为成功的模拟。下面列举一些例子以供参考:
研究蛋白折叠运动(动力学性质)
研究水溶液中有机小分子聚集行为,计算回转半径、径向分布函数等
隐式溶剂环境(GB模型)中研究DNA-小分子的相互作用
入口
功能入口:左侧菜单栏【计算方案】->【大方案】->【分子动力学】->【分子动力学(Amber 20)】 步骤提交任务1. 上传分子结构 上传Amber的拓扑文件和坐标文件。 坐标文件可以是:
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 ·Å 的约束力。注:目前仅支持对原子设置位置约束,将来会增加距离、角度、二面角度等约束方式。
4. 点击【提交】。 分析结果热力学性质每一个阶段均给出热力学数据,包括:能量(总能量、动能、势能)、温度、压强、体积和密度。分析这些性质,判断模拟是否达到预期效果。 能量 在能量最小化(Minimization)阶段,能量一般会先下降后平稳;在加热(Heat)阶段,能量应在极短时间内达到平稳(范德华能量平稳下降,其他能量平稳上升);在平衡(Equilibration)或生产(Production)阶段,能量应保持平稳。 温度 Minimization阶段无温度;在Heat阶段,温度应从初始温度几乎沿着直线上升(或下降)至目标温度;而当模拟进入Equilibration或Production阶段,温度应在目标温度附近波动(下图红色移动平均线显示平均值在300 K附近波动)。
压强 只有(显式溶剂环境)NPT系综才有压强,其曲线应在极短时间里达到期望值并保持平稳(下图显示,红色移动平均线在1 atm附近波动)。
体积 体积与压强负相关,只有(显式溶剂环境)NPT系综才有,其曲线与压强大致相反。
密度 密度与体积是倒数关系,也是只有(显式溶剂环境)NPT系综才有,其曲线与体积刚好相反。 RMSD该处给出了Equilibration和Production阶段的RMSD数据。通过RMSD分析,了解体系的构象变化,是否达到平衡。 很多后续分析工作都建立在平衡状态的基础上,例如:进一步采样(production)、氢键分析、聚类分析、结合自由能分析;等等。
RMSD该处给出了Equilibration和Production阶段的RMSD数据。通过RMSD分析,了解体系的构象变化,是否达到平衡。 很多后续分析工作都建立在平衡状态的基础上,例如:进一步采样(production)、氢键分析、聚类分析、结合自由能分析;等等。
该表给出各帧(Frame)的溶质-溶质(Solute-solute)氢键、溶质-溶剂(Solute-solvent)氢键和水桥(Bridge)个数。
该表给出各帧(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为氢原子。
该表给出各个氢键的统计数据:Name为形成氢键的残基对(如,Glu191的OE1原子与Ser132的OG原子和HG原子形成氢键),Frames为轨迹中该氢键存在的帧数,Fraction为比例(等于Frames除以总采样帧数),Avg. Distance和Avg. Angle分别为平均键长和键角,Lifetimes为生命周期数,Max. Lifetime和Avg. Lifetime分别为最长寿命和平均寿命。 OE1表示名为OE1的氧原子,同理,HG为氢原子。
如下图所示,Bridged Residues为形成水桥的溶质残基(例如,Cyx95、Trp114和DE3126之间存在水介导的氢键),Frames为该水桥出现的帧数。 分子结构该处给出各个阶段最后一帧的结构,其中,文件名带有 -dry的是删除水和离子后的结构。 文件列表该处给出分子动力学模拟的输出文件,其中*.rst文件(坐标文件)可作为后续模拟的起点(新建任务,上传到坐标文件处)。
|