分子动力学模拟概述
四、经典力学
分子动力学模拟方法是建立在牛顿第二定律之上的,即牛顿运动方程:
F=ma
其中F为作用在粒子上的力,m表示粒子的质量,a表示粒子的运动加速度。只要知道作用在每一个粒子上的力,即可知道系统中每一个粒子的加速度。利用这个方程,模拟过程记录系统中各个粒子的位置、速度和加速度随时间变化的系统轨迹文件;而从这个轨迹,即可以求出系统各个性质的平均值。这是一个确定性的方法,系统中的每一个粒子的任一时刻位置跟速度一旦知道,系统过去或者未来的性质都可以推算得到。分子动力学模拟方法虽然要求计算时间长,计算资源多;但是计算机已经变得越来越快,价格也越来越便宜。现在的分子动力学一般在纳秒级别,毫秒级别还鲜有报道。
--------------------------------
牛顿运动方程:

其中,Fi表示作用在质点i上的力,mi表示质点i的质量,而ai表示质点i的加速度。作用力可以表达成势能的梯度函数,即

将上面两个式子进行整合,可得
其中V表示系统的势能。牛顿运动方程即表达成系统势能与随时间变化的质子位置之间的表达式。
--------------------------------
牛顿第二定律的一个简单应用:

如果加速度恒定

有:
因为速度是位移随时间的变化

有:
将上列式子进行整个,可以得到某一时刻t与位置x的关系:

其中x0和v0分别表示初始位置和初始速度。加速度可以表示为势能对位置的导数:

所以,为了求取一个运动轨迹,必须知道系统中粒子的初始位置,初始速度分布和加速度,这些都可以从势能导数方程推导得到。比如,开始时间的粒子的位置和时间可以准确得到某一时刻t的位置和速度。初始位置可以系统实验结构数据得到,如X光衍射或者NMR实验。系统初始速度分布主要来自计算机随机分布,同时根据系统温度进行修正以确保系统的整体动量为零:

其中,速度vi是根据温度在一个波尔兹曼分布或者高斯分布中随机抽取,所以原子i可以得到在温度T在x方向上的速度vx。

反过来,系统温度可以从速度求得:

其中N表示系统中粒子数。
-------------------------------
算法综合
由于系统地势能是系统内所有原子的位置的函数,方程十分复杂而没有解析解,只能用数值方法求解。几种常用的数值算法有:
1、Verlet算法
2、Leap-frog算法
3、Velocity Verlet算法
4、Beeman's算法
在选择不用算法是,要注意算法是否能保持系统内部动量的恒定,并可以高效地解决问题,同时允许进行比较长的步长模拟。
-------------------------------
所有以上的算法都假定系统中原子的位置、速度和加速度可以利用泰勒张开近似求得:

其中,r表示位置,v表示速度(随时间变化的第一个因变量),a为加速度(随时间变化的第二个因变量)。
Verlet算法可以表达成以下形式:

将两式合并可得到:

Verlet算法利用在t和t - &delta t时刻的位置和加速度以计算在t + &delta t时刻的新位置。这个算法并没有计算原子的精确速度。这个算法主要优点有:i )可以直接求解;ii )计算中加饥结果少,可以节省内存。但是这个的缺点就是精确度比较小。
-------------------------------
Leap-frog算法

在这个算法中,先计算在t+1/2&delta t时刻的速度,然后再用其计算在t+&delta t时刻的位置。这样,速度根据速度进行计算,而位置又根据得到的速度进行计算,一步步跳跃向前。这个算法得优势是速度可以精确计算,但是缺点是不能在同一时刻求得速度和位置。t时刻的速度可以用下列关系大约求得:

-------------------------------
Velocity Verlet算法
这个算法在同一时刻求取位置、速度和加速度,但是精确度不高,主要表达为:

-------------------------------
Beeman算法
这个算法从Verlet算法演变而来,表达为:

这个算法比较精确的表达了速度,同时也确保能量守恒。但是由于算式复杂,计算量较大。
-------------------------------