【分子动力学教程】分子动力学模拟（Amber 20）
用途使用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: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）和化学势（）。体系是开放系统，与大热源大粒子源热接触平衡而具有恒定的温度。
那在动力学模拟中如何选择系综？本方案提供NVT、NPT和NVE三种系综，前两种最为常用，下面详细讨论。
常规体系的模拟是在等温等压（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·Å的约束力。注：目前仅支持对原子设置位置约束，将来会增加距离、角度、二面角度等约束方式。
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文件（坐标文件）可作为后续模拟的起点（新建任务，上传到坐标文件处）。
页:
[1]