书城自然科学分子模拟力场方法与应用
8392100000012

第12章 分子动力学模拟概述(2)

,图4-3显示利用开关函数所调整的非键相互作用势能图,由此图可知,利用开关函数可将由RsRc区间的非键势能函数逐渐调整为零,形成连续且符合计算要求的势能函数。通常Rs不宜过小,约取0.9Rc。

4.5积分步长的选取

在分子动力学计算中,最重要的工作是如何选取适当的积分步长t(integrationstep),以节省计算的时间而又不失去计算的准确性。通常的原则是积分步长应小于系统中最快运动周期的1/10。以水分子的动力学计算为例,水分子内的运动为键长和键角的变化,分子间的运动为质心的平动与分子的转动。分子间的运动源自于vanderWaals作用,一般较慢,但分子内的运动则较快。由红外光谱的实验可知,水分子的最大振动频率为O-H的伸缩振动,

1=3660cm-1则振动频率可以表示为:

因此,分子动力学模拟的积分步长最大为t0.910-16sec。(目前,可能最快的就是O-H键的伸缩振动)。再以氩原子的分子动力学计算为例,氩原子间的势能函数只考虑vdw相互作用,利用Lennard-Jones(12-6)表示,则势能函数对坐标的二次微商可写为:

其中,为势阱深度,为碰撞直径,对氩原子体系通常取为0.24kcal/mol和3.504。

当势能有极小值时,6minr2,此时的rmin表示最低能量距离。将此代入到能量的二阶导数中得到:2minU"(r)57.14/=1.117kcal/mol/。因为势能函数最低点的二次微分相当于力的常数k,即U"(r)kmin,由简谐振动的关系式:

k4π22(4-28)式中,为折合质量(reducedmass)。氩原子的折合质量为:

故在氩原子体系中最快的振动周期为T=2.90610-12sec,积分步长t则可取约为2.910-13sec。比较以上的结果可知,在以分子动力学模拟过程中,氩原子的积分步长可取为计算水分子体系的积分步长的300倍。分子动力学计算中,积分步长愈大则可研究的时间范围就愈长。例如,研究时间范围为10-9sec的氩原子系统,需要3500步的分子动力学模拟计算,但研究相同时间范围内的水分子系统则至少需1.05106步的分子动力学模拟计算。以上利用Lennard-Jones(12-6)的势能函数估计原子系统积分步长的方法也可应用于不同原子的混合系统或多原子分子系统。在多原子分子体系中通常取(12-6)势能函数中最大的值原子来估计运动最快的原子间的振动频率,以此获得对积分步长的估算。例如水分子间的vdw作用包括氧原子-氧原子、氧原子-氢原子和氢原子-氢原子原子对间的作用,其中,氧原子-氧原子的值最大,故应以此来估计积分步长的标准。

时间步长大概是最快运动周期1/10的需求是一个非常严格的要求,但是当这些高频振动并不是我们所感兴趣的,并且和整个体系的行为并没有多大的关系时,如何来选取积分步长呢?有效地解决这个问题的方法就是“冻结”这样的高频运动,也就是通过限制键长为其平衡键长(刚体结构),同时在维持分子间和分子内力的存在条件下,保持其它自由度的变化,就可以允许一个较大的时间步长,这就是限制性的动力学计算方法。表4-1给出在分子动力学模拟计算中各种模拟系统一般采取的积分步长。

4.6分子动力学计算流程

执行分子动力学计算的起点就是将一定数目的分子置于立方体盒子中,使其密度与实验上获得的密度相一致,选定计算的温度,即可以开始计算。计算时,首先必须知道系统中分子的初始位置与速度,通常可将分子随机置于盒子中,或者是取其结晶形态的位置排列作为体系的初始位置。系统中所有原子运动的动能必需满足:

式中,KE为系统的总动能,N为总原子数,kB为波尔兹曼常数,T为绝对温度。

原子的初速度可利用上式产生,取一半的原子向右运动,另一半原子向左运动,或是令原子的初速度呈高斯分布(GaussianDistribution),则体系的总动能为3/2NkBT。

产生原子的初始位置与初速度后,就可以进行分子动力学模拟计算。

由初始位置与速度开始,分子动力学模拟计算的每一步都产生新的速度与位置,由新产生的速度可计算系统的温度Tcal:

若系统的计算温度Tcal与所定的温度T相比过高或过低,则需要校正速度,一般容许的温度范围为:

若计算的温度超过此容许的温度范围,则将所有原子的速度乘以一个校正因子。

这样系统的温度可表示为:

也就是将计算的温度校正回系统的设定温度。实际执行分子动力学模拟的过程中,在计算开始时每隔数步就需要校正温度;随后校正的时间间隔增长,每隔数百步或数千步才需校正,直至原子的速度不需要再校正,而系统的总动能在3/2NkBT上下10涨落,此时认为系统已达到热平衡(thermalequilibrium)状态,在达到热平衡状态前的轨迹与速度不需要保存,因其物理意义不够严密,只有当系统达到热力学平衡状态后,才开始存储计算的轨迹与速度。图4-4为Verlet跳蛙法的分子动力学模拟计算流程图,这种计算方法称为一般性动力学计算(conventionalmoleculardynamic),计算中系统的原子数N,体积V和总能量E维持不变,相当于统计力学中的微正则系综(micro-canonicalensemble,NVE)。

4.7分子动力学模拟的初始设定和平衡态

执行分子动力学计算必需选取适当的初始条件,如起始位置、速度、执行温度、积分步长等。若初始条件选择不当,则往往需要浪费相当长的计算时间才能达到系统的热平衡,甚至于根本无法实现热平衡。初始结构的选取,愈接近模拟系统的结构愈好。如果模拟面心立方的固体系统而选择了体心立方作为初始结构就是不明智的。通常在模拟前最好执行一些实验性的预测,以确定所设的系统中没有太高能量的作用,否则会导致模拟过程的不稳定。也可利用能量最小化的方法,以最低能量的结构作为模拟的起点,从而避免高能作用的产生。初始构象可以从实验数据中获得,也可以从理论模型方法中获得,或是从实验和理论相结合的方法中得到。另一个必不可少的过程就是分配给原子一定的初速度,它可以从定温下的Maxwell-Boltzmann分布中随机获取。

Maxwell-Boltzmann方程给出了质量为mi的i原子在温度T时有沿着x轴方向ixv速度的可能性。Maxwell-Boltzmann分布是Gaussian分布,可以通过随机数种子获取。

大多数的随机数种子产生的随机数在01范围内,把这样的随机数种子转变成Gaussian分布样本是很直接的,从Gaussian正则分布中获得平均值为(x>,变数波动为22xx2的可能性是:

从这样的分布中获得的数值(x)和高斯分布中的(x)相对应。平均值(x>和波动用下式表示:

xxx(4-37)初始速度也可以从单位分布或是简单Gaussian分布中得到,在这些方法中,Maxwell-Boltzmann获得速度的方法是最快的。

建立了体系并且分配了初始速度,就可以按照上一节介绍的分子动力学流程开始计算体系的运动轨迹。平衡态的到达是执行分子动力学模拟过程中最重要的一个步骤,这个过程的目的就是使系统从初始状态达到平衡状态,以使体系彻底地“忘记”初始的各种性质。在实现平衡态的过程中,调节体系的各种参数和各种组态,当利用这些参数所计算的动力学性质达到稳定的数值时,就达到了体系的产物空间。

在产物空间中可以计算各种动力学性质以及其它相关数据。用来描述体系是否达到平衡态的各种参数在某种程度上要依赖于所研究的体系,但是体系的动能、势能、总能量、速度、温度以及压力应该是保持不变的。正像我们所期望的,在巨正则系综下,势能是在不断波动的,而体系的总能量应该是保持不变的,速度的各个部分应该能够描述Maxwell-Boltzmann的各个分量(x,y,z轴),动能则应该在x、y、z轴三个方向上平均分配。通常我们要计算在给定温度下的分子动力学模拟,因而就要利用速度在平衡过程中来调节温度。

4.8等温计算方法

以上所介绍的分子一般性动力学模拟计算适用于微正则系综(NVE),能量恒定,温度在其平衡值附近波动。除了微正则系综外,分子动力学模拟计算也可以处理其它类型的系综,如正则系综(NVT)、等温等压系综(NPT)和巨正则系综等。使用何种系综进行计算要根据实际体系的需要来确定。例如研究材质的相变化,多采用等温等压系综,在定压下计算各温度时系统的能量或比热,由此来判断相变化的温度。

而研究一般溶液的行为,如水合能(hydrationenergy),解离效应等多采用正则系综(NVT)。简单介绍两种定温计算法,Andersen和Nosé-Hoove方法,Nosé-Hoover的定温计算方法较Andersen的计算方法稳定,目前普遍被采用。

在Anderson[196]提出的恒温方法中,体系与一个强加了指定温度的热浴相耦合。

与热浴的耦合利用随机选取的粒子上的随机脉冲力表示。这些与热浴的随机碰撞可以看做MonteCarlo移动,将体系从一个等能面转移到另一个等能面。在随机碰撞之间,体系遵从标准的牛顿运动定律,在恒定能量的情况下演变。随机碰撞保证可根据它们的Bolzmann权重来遍历所有可能的等势面。在Anderson的NVT系综下,粒子的位置与速度分别为:

在Anderson的等温分子动力学模拟计算方法中,恒定的温度通过与严格热浴的随机碰撞获得。Nosé则向人们展示了也可进行恒定温度下的确定的MD模拟[197~198]。

Nosé的这一方法基于扩展的Lagrange函数的巧妙使用,即Lagrange函数包括附加的、人为指定的坐标和速度。目前对恒温动力学模拟更多地使用在Hoover[199~200]公式中的Nosé算法。这一方法是以变量s调节过高或过低的动量p,使体系的温度保持不变。调节后的物理量为:

其中,上标()表示调节后的物理量。由此可以得运动方程为:

其中,Q是有效质量(effectivemass),是人为设定值,通常选取5或10,为热力摩擦系数(thermodymanicfrictioncoefficient)。执行分子动力学模拟计算时可将以上的运动方程积分来计算时刻t的物理量。