分子动力学模拟

分子动力学模拟
分子动力学就是一门结合物理,数学与化学的综合技术。

分子动力学就是一套分子模拟方法,该方法主要就是依靠牛顿力学来模拟分子体系的运动,以在由分子体系的不同状态构成的系
统中抽取样本,从而计算体系的构型积分,并以构型积分的结果为基础进一步计算体系的热力
学量与其她宏观性质。

这门技术的发展进程就是:
1980年:恒压条件下的动力学方法(Andersenの方法、Parrinello-Rahman法)
1983年:非平衡态动力学方法(Gillan and Dixon)
1984年:恒温条件下的动力学方法(能势‐フーバーの方法)
1985年:第一原理分子动力学法(→カー?パリネロ法)
1991年:巨正则系综的分子动力学方法(Cagin and Pettit)、
最新的巨正则系综,即为组成系综的系统与一温度为T、化学势为μ的很大的热源、粒子源相接触,此时系统不仅同热源有能量交换,而且可以同粒子源有粒子的交换,最后达到平衡,这种系综称巨正则系综。

进行分子动力学模拟的第一步就是确定起始构型,一个能量较低的起始构型就是进行分子模
拟的基础,一般分子的其实构型主要就是来自实验数据或量子化学计算。

在确定起始构型之
后要赋予构成分子的各个原子速度,这一速度就是根据玻尔兹曼分布随机生成,由于速度的分布符合玻尔兹曼统计,因此在这个阶段,体系的温度就是恒定的。

另外,在随机生成各个原子的
运动速度之后须进行调整,使得体系总体在各个方向上的动量之与为零,即保证体系没有平动位移。

由上一步确定的分子组建平衡相,在构建平衡相的时候会对构型、温度等参数加以监控。

进入生产相之后体系中的分子与分子中的原子开始根据初始速度运动,可以想象其间会发生吸引、排斥乃至碰撞,这时就根据牛顿力学与预先给定的粒子间相互作用势来对各个例子的
运动轨迹进行计算,在这个过程中,体系总能量不变,但分子内部势能与动能不断相互转化,从而体系的温度也不断变化,在整个过程中,体系会遍历势能面上的各个点,计算的样本正就是在这个过程中抽取的。

用抽样所得体系的各个状态计算当时体系的势能,进而计算构型积分。

作用势的选择与动力学计算的关系极为密切,选择不同的作用势,体系的势能面会有不同的形状,动力学计算所得的分子运动与分子内部运动的轨迹也会不同,进而影响到抽样的结果与抽
样结果的势能计算,在计算宏观体积与微观成分关系的时候主要采用刚球模型的二体势,计算系统能量,熵等关系时早期多采用Lennard-Jones、morse势等双体势模型,对于金属计算,主要采用morse势,但就是由于通过实验拟合的对势容易导致柯西关系,与实验不符,因此在后来的模拟中有人提出采用EAM等多体势模型,或者采用第一性原理计算结果通过一定的物理方法
来拟合二体势函数。

但就是对于二体势模型,多体势往往缺乏明确的表达式,参量很多,模拟收敛速度很慢,给应用带来很大困难,因此在一般应用中,通过第一性原理计算结果拟合势函数
的L-J,morse等势模型的应用仍非常广泛。

分子动力学计算的基本思想就是赋予分子体系初始运动状态之后,利用分子的自然运动在相
空间中抽取样本进行统计计算,时间步长就就是抽样的间隔,因而时间步长的选取对动力学模
拟非常重要。

太长的时间步长会造成分子间的激烈碰撞,体系数据溢出;太短的时间步长会降低模拟过程搜索相空间的能力,因此一般选取的时间步长为体系各个自由度中最短运动周期
的十分之一。

但就是通常情况下,体系各自由度中运动周期最短的就是各个化学键的振动,而这种运动对计算某些宏观性质并不产生影响,因此就产生了屏蔽分子内部振动或其她无关运
动的约束动力学,约束动力学可以有效地增长分子动力学模拟时间步长,提高搜索相空间的能
力。

以下就是做模拟的一般性步骤,具体的步骤与过程依赖于确定的系统或者软件,但这不影响我们把它当做一个入门指南:
1）首先我们需要对我们所要模拟的系统做一个简单的评估,三个问题就是我们必须要明确的:
做什么(what to do)为什么做(why to do)怎么做(how to do)
2）选择合适的模拟工具,大前提就是它能够实现您所感兴趣的目标,这需要您非常谨慎的查阅文献,瞧瞧别人用这个工具都做了些什么,有没有与您相关的,千万不要做到一半才发现原来这个工具根本就不能实现您所感兴趣的idea,切记！
考虑1:软件的选择,这通常与软件主流使用的力场有关,而软件本身就具有一定的偏向性,比如说,做蛋白体系,Gromacs,Amber,Namd均可;做DNA,RNA体系,首选肯定就是Amber;做界面体系,DI_POLY比较强大,另外做材料体系,Lammps会就是一个不错的选择。

考虑2:力场的选择。

力场就是来描述体系中最小单元间的相互作用的,就是用量化等方法计算拟合后生成的经验式,有人会嫌它粗糙,但就是它确确实实给我们模拟大系统提供了可能,只能说关注的切入点不同罢了。

常见的有三类力场:全原子力场,联合力场,粗粒化力场;当然还有所谓的第一代,第二代,第三代力场的说法,这里就不一一列举了。

再次提醒注意:必须选择适合我们所关注体系与我们所感兴趣的性质及现象的力场。

3）通过实验数据或者就是某些工具得到体系内的每一个分子的初始结构坐标文件,之后,我们需要按我们的想法把这些分子按照一定的规则或就是随机的排列在一起,从而得到整个系统的初始结果,这也就是我们模拟的输入文件。

4）结构输入文件得到了,我们还需要力场参数输入文件,也就就是针对我们系统的力场文件,这通常就是由所选用的力场决定,比如键参数与非键参数等势能函数的输入参数。

5）体系的大小通常由您所选用的box大小决定,我们必须对可行性与合理性做出评估,从而确定体系的大小,这依赖于具体的体系,这里不细说了。

6）由于初始构象可能会存在两个原子挨的太近的情况(称之为bad connect),所以需要在正式模拟开始的第一步进行体系能量最小化,比较常用的能量最小化有两种,最速下降法与共轭梯度法,最速下降法就是快速移除体系内应力的好方法,但就是接近能量极小点时收
敛比较慢,而共轭梯度法在能量极小点附近收敛相对效率高一些,所以我们一般做能量最小化都就是在最速下降法优化完之后再用共轭梯度法优化,这样做能有效的保证后续模拟的进行。

7）以平衡态模拟为例,您需要设置适当的模拟参数,并且保证这些参数设置与力场的产生相一致, 举个简单的例子,gromos力场就是用的范德华势双截断来定范德华参数的,若您也用gromos力场的话也应该用双截断来处理范德华相互作用。

常见的模拟思路就是,先在NVT下约束住您的溶质(剂)做限制性模拟,这就是一个升温的过程,当温度达到您的设定后,接着做NPT模拟,此过程将调整体系的压强进而使体系密度收敛。

经过一段时间的平衡模拟,在确定系统弛豫已经完全消除之后,就可以开始取数据了,如何判断体系达到平衡,这个问题就是比较技术性的问题,简单的讲可以通过以下几种方式,一,瞧能量(势能,动能与总能)就是否收敛;二,瞧系统的压强,密度等等就是否收敛;三瞧系统的RMSD就是否达到您能接受的范围,等等。

8）运行足够长时间的模拟以确定我们所感兴趣的现象或就是性质能够被观测到,并且务必确保此现象的可重复性。

9）数据拿到手后,很容易通过一些可视化软件得到轨迹动画,但这并不能拿来发文章。

真正的工作才刚刚开始——分析数据,您所感兴趣的现象或性质只就是表面,隐含在它们之中的机理才就是文章中的主题。

分子动力学可以用于NPT,NVE,NVT等系综的计算,就是一种基于牛顿力学确定论的热力学计
算方法,与蒙特卡洛法相比在宏观性质计算上具有更高的准确度与有效性,可以广泛应用于物理,化学,生物,材料,医学等各个领域。

另外,在实际应用中,经常把分子动力学方法与蒙特卡罗法联合使用。