分子动力学模拟简介在线视频

同学们

大家好

这一节我们来学习

分子动力学简介

在前面的课程中

我们已经学习和强调过

分子势能面这一概念

对于计算化学的重要性

分子的稳定结构

一般可以通过寻找

势能面上的局域极小点

也就是通过优化计算来得到

基于分子结构

我们可以进一步求算

物质的物理化学性质

这一方法对于原子数很少的分子

特别是势能面极小点附近

形状比较陡峭的体系

是非常适用的

通常情况下

这样的分子基本上会保持在其

能量极小结构附近

所以它的性质

也基本由能量极小结构所决定

但是我们知道

由于存在着量子效应

每个化学键

每个分子

都存在零点振动能

这意味着即使在绝对零度下

对于氢分子这样的简单物质

它的键长也会在

一定范围内变化

而不会固定在

优化得到的最优键长上

随着温度的升高

分子中的原子运动

会更加剧烈

运动范围也相应增加

因此在通常情况下

因为热运动的存在

分子存在一定的构象变化区间

特别对于分子内

原子数比较多

结构比较复杂的分子来说

由于分子中的作用

也相应变得很复杂

分子会表现出相当的柔性

比如一些高聚物分子

生物分子等等

分子内存在着大量

很容易发生扭转的化学键

从而导致分子构象

极易发生变化

这样的分子实际上存在着

一个范围比较大的构象区间

而且随着温度的变化

分子的构象区间

还会发生变化

此时当我们提到

分子结构这个概念的时候

就不能仅仅考虑它的

势能面上的能量极小结构了

举一个极端的粒子

自然界中存在着一些

天然无序的蛋白质

通常情况下

它们根本就不存在一个

确定的三维结构

如果我们仍然依据

势能面极小点的结构

来对这样的分子进行研究

就会产生比较大的误差

另外一方面

对于蛋白质这样的生物大分子

它能够正常发挥其生物功能

也往往有赖于

它的构象的柔性

因此要对这样的分子的

性质进行计算

就必须考虑它们在

势能面上的构象分布

以及构象变化情况

而不能仅仅考察它们的

能量极小结构

而且实验测量得到的物质性质

往往也是针对大量分子

在时间平均和系综平均的

基础上进行的

因此为了更好的与实验进行对比

计算中也需要考虑

大量分子的构象分布

对体系性质的影响

如果我们把分子中的原子

当成经典粒子来处理

那么分子的变化

可以用其中所有原子的位置

和动量来进行完整的描述

在三维空间中

每个原子对应了六个自变量

包括XYZ三个位置坐标

和XYZ方向的三个动量分量

而整个分子体系的性质

则由6N个变量来决定

我们将这6N个变量

构成的空间

叫做对应分子的相空间

分子在任意时刻

根据其位置和运动

对应了相空间中的一点

X‘可以用这个表达式来表示

为了简化

我们可以把所有原子的位置

合并简写成q

而所有原子的动量

合并简写为P

所以相空间中的一点X’

可以用p，q来表示

它代表了分子的某一特定状态

而且随着时间的变化

分子在相空间中的位置

也会发生变化

从而在相空间中

生成对应的轨迹

相空间坐标p，q

决定了分子体系的性质

因此体系的性质A

可以表达成p，q的函数A（p，q）

由于p，q会随着时间而变化

因此A也会相应发生变化

在包含大量分子的体系中

不同分子对应相空间中的位置p，q

和相应的性质A（p，q）也会不同

在温度体积和粒子数

都一定的情况下

一个体系在平衡状态下的

某一性质A的平均值

应该由分子在相空间中的分布决定

表达成计算式就是A的平均值

由A（p，q）与体系中

分子处于相应点p，q的概率

相乘进行加权平均得到

这个表达式中的A（p，q）

代表了性质A

随体系中所有粒子的坐标

和动量的变化情况

而P（p，q）则代表了体系

在某一个相点对应坐标为q

动量为p所出现的概率

根据统计力学

我们知道这一几率

与该相点所具有的能量有关

根据体系的配分函数q

和该相点p，q所具有的能量E(p，q)

可以求得相应的几率P(p，q)

在这个几率表达式中

KB是波尔兹曼常数

T是温度

虽然原理非常简单

但是对于复杂体系来说

由于它们的配分函数

往往没有解析表达式

只能通过数值方法来拟合

这就要求我们

要将相空间中所有的点p，q

和对应的性质A（p，q）都计算出来

然后进行统计平均

然而由于相空间

是一个极大的空间

我们考虑一个

含有一百个粒子的体系

仅仅让其中每个粒子的动量

和坐标可以取正负2个值

相空间中就会有

2的600次方个点

这个数大概是

4.15×10的180次方

如果要把这些点

都找到显然是不可能的

但是分析一下配分函数

我们可以看到

由于相空间中某一点出现的概率

与该点对应的能量

呈负的指数相关

这意味着对于那些能量很高的点

它们出现的几率将会极小

从而对体系的性质的影响

也会微乎其微

甚至可以忽略不计

因此对我们来说

最关键的还是怎样找出

相空间中最重要的

贡献最大的那些点来

相对于完整地相空间来说

它们只占了比较小的一部分

在讨论怎样确定相空间中

哪些点更重要之前

我们首先来考察一下

如何来求算一个体系的平均性质

如果我们从分子的某一相点出发

该相点应该对应了

分子的一个合理的结构

分子中的原子有一定的动能

这样分子就会随着时间

而在相空间中

一定的区域运动形成轨迹

并对相空间进行取样

因此如果我们要求算

该分子所组成的体系的平均性质

可以一定的时间间隔

从相空间中找到对应的相点

并基于它求算对应的性质

最后对所有取出的点

进行平均就可以了

那么在这个计算式中

M就是取样的次数

A（ti）表示某一时刻ti

取得的相点所对应的性质A的值

如果我们能对完整的轨迹

进行连续的取样

那么平均值的求算

就转变成一个积分式

基于统计力学的各态历经假设

如果我们取样的时间足够长

那么从任意的起点出发

都可以使轨迹遍历完整的相空间

而且对相点取样的几率

满足波尔兹曼分布

由于一个宏观体系中

各个分子所表现出来的平均值

就取决于其在一定温度下

对应的玻尔兹曼分布

因此基于各态历经假设

我们可以认为对于处于一定系综

中的大量粒子的统计平均性质

可以通过对一个粒子

在相空间中的轨迹

进行求时间平均来得到

从而将宏观多粒子的统计平均问题

转变成了一个微观单粒子的

时间平均问题

现在的问题就是

我们应该怎样来得到

一个粒子在相空间中的轨迹了

对于服从经典力学的物体来说

因为一个相点已经规定了

粒子的位置和动量

因此只要给定了初始相点

可以认为整个轨迹都是确定的

例如我们来考察一个

经典的一维谐振子

它的相空间是二维的

仅有一个位置

和一个动量

相空间中的原点

显然对应了弹簧处于平衡长度

而振子动量为零的静息状态

如果起点是原点

那么由于这一点是驻点

所以谐振子在相空间中的轨迹

就仅仅是这样一个点

现在我们来考察另外一种情况

对于力常数为K的

弹簧上的谐振子

考虑在t0时刻

从长度为B的偏离平衡点的位置

将振子释放

此时振子的动量为零

但是在弹簧恢复力的作用下

振子开始朝着位置坐标

而减小的方向运动

此时坐标仍然为正值

同时由于弹簧的势能转换成

振子的动能

振子开始具有指向坐标轴

负方向的动量

所以动量是负值

随着振子运动经过平衡位置

坐标为零

它的动量也达到最大的负值

随着弹簧开始缩短

坐标进一步减少

变成负值

而振子也开始减速

动量相应增加

最终随着振子动量达到零

弹簧被压缩到最短

此时的位置应该是-B

此时振子又开始朝

坐标轴正方向运动

动量继续增加

坐标也开始增加

直到振子再一次达到B

点相应的动量

又开始由正转为负值

从而在相空间中

形成这样的椭圆形的

往复的轨迹

我们现在可以来考察一下

振子运动过程中的动力学过程

如果在时间t1振子的位置在qt1

经过任意的时间到t2时

根据牛顿动力学方程可以求的

qt2的表达式

式中的Pt为t时刻振子的动量

m为振子的质量

同样我们可以由t1时刻的动量pt1

求得t2时刻的动量pt2

表达式中的at表示t时刻

粒子的加速度

由牛顿第二定律可以得到

a等于F除以m

F为振子的受力

可以由胡克定律得到

因为谐振子的势能表达式非常简单

因此我们可以借此得到

任意时刻t的位置qt

和动量pt的解析表达式

这就是我们所看到的

相空间中的谐振子的

椭圆形轨迹的方程

然而由于分子的势能函数

远比谐振子要复杂

我们对分子的相空间轨迹

就无法得到解析形式的方程了

此时我们利用

欧拉近似将相空间中的轨迹

对无穷小时间的积分

改写为对有限时间△t的

差分计算式

根据坐标和动量的差分表达式

一旦给定了体系中

所有粒子的初始位置和动量

并且能够计算粒子

在任意时刻的受力的话

我们就可以据此来近似模拟

真实相空间中的轨迹了

通常来说

对于化学体系

其中原子的初始位置q

其实代表了

所要模拟的分子的初始结构

一个合理的初始结构

可以通过对分子

进行能量极小化计算

也就是优化计算来得到

而粒子的初始动量

则可以根据要模拟的温度

来对体系中的原子进行随机指定

经过指定后

粒子的平均动量与温度之间

应该满足这样一个对应的关系式

这也是温度的统计

力学定义式

在这个关系式中

N是原子总数

n代表的是

在模拟时被固定的坐标数量

这些坐标上没有动能

所以需要先从温度计算中扣除

Pi(t)表示t时刻

第i个原子的动量

mi代表这个原子的质量

为了能够进行动力学计算

我们需要计算原子的受力

也就是需要势能函数

显然对于原子数

很多的体系来说

分子力学方法比量子化学方法

要经济很多

对于化学体系

由于我们无法得到

动力学轨迹的精确解析表达式

因此能否精确重现

体系在相空间中的真实轨迹

将取决于我们对

动力学过程的计算的具体设置

特别是积分步长△t的设置

一方面我们把时间步长△t

设置的越小

计算得到的轨迹

就越接近真实的轨迹

当△t达到无穷小

就还原得到精确的轨迹

但是这样会导致

计算量大到无法接受

当然我们也可以通过

把时间步长△t设置的更大

来减少计算量

但是这样可能

又会导致我们得到的轨迹

偏离真实轨迹

甚至得到错误的结果

我们还是以刚才的一维谐振子为例

来进行说明

振子在弹簧的作用下

以一定的周期往复运动

假设振子的运动周期是τ

这意味着

在二分之τ的时间内

振子可以从离墙最远处

运动到离墙最近的地方

但是如果我们设定的△t步长

如果大于这个值的话

那么得到的轨迹中

振子的振动幅度

就会超过实际的振动幅度

从而产生误差

如果△t设置的太大

导致振子的振幅过大

我们甚至可以得到

振子从最远端

回复运动后进入墙体的荒谬结果

因此

为了避免动力学计算过程中

产生过大的误差

甚至带来不稳定

我们需要对

模拟的步长进行谨慎的选择

一般而言

计算所使用的步长

至少需要比体系中

最快的周期运动的周期

小一个数量级

分子体系中

最快的运动

是化学键的振动

而其中最快的

又是含有氢原子的化学键的振动

周期一般在10的-14次方秒

也就是十飞秒左右

因此如果需要模拟一个分子

在相空间中的轨迹

我们通常需要将步长

设置为小于一飞秒

而通常我们所感兴趣的

分子的动力学变化过程

例如蛋白质的折叠过程

往往发生在微秒

甚至毫秒的时间尺度上

这意味着动力学模拟

需要进行10的9次方步以上的计算

才能得到我们感兴趣的结果

而每一步还都需要计算出

分子的能量和梯度

这显然是无法承受的

针对这一问题

目前已经提出了很多解决方案

最简单的方案是

在计算中将所有含有氢原子的化学键固定住

从而可以将步长相应增大

而减少计算量

下面我们来对

动力学积分过程做一下简单的讨论

依据欧拉近似后的动力学方程

我们从粒子初始时刻t0的位置qt0

动量pt0和受力出发

指定△t

我们就可以模拟得到

完整的动力学过程中的轨迹

这样的方法又叫初始值法

根据我们刚才的分析

初始值法非常简单

但是误差较大

动力学模拟过程容易不稳定

如果我们仅仅讨论

轨迹中的位置变化情况

更好的近似方法是将

qt展开成泰勒级数的形式

表达式中的

vt和at分别代表了

粒子在t时刻的速度和加速度

它们分别是位置qt的

一阶和二阶导数

1967年Verlet基于泰勒展开

提出了一种动力学轨迹积分办法

它的基本思路是

先将展开式中的△t

替换成负的△t

可以得到q（t-△t)的表达式

然后将q(t+△t)

与q（t-△t)进行加和

可以消去其中的奇次项

例如含有速度vt的这一项

将在加和中消失

进一步将级数

在二阶处截断就可以得到一个

仅仅通过q（t-△t)

和qt以及at

来获得q(t+△t)的表达式

也对于一个粒子

其下一时刻的位置

将由它上一时刻的位置

和当前的位置

以及现在的受力所决定

而完全不需要粒子的速度信息

仅仅是计算第一步

需要用到初始值法

来设定并启动计算

由于计算中不需要粒子的动量信息

当我们需要求算的物理量

或者性质仅与结构有关时

Verlet方案就显得非常的方便

然而由于粒子的速度

在模拟过程中完全没有涉及

导致我们无法控制体系的温度

因此为了能在动力学模拟过程中

同时考虑所有粒子的动量与位置

有人提出了改进的Verlet方案

首先可以将qt表示成q(t+1/2△t)

再减1/2△t

而将q(t+△t)表示成

q(t+1/2△t)

再加1/2的△t

并且呢将这两个式子以

1/2△t为变量进行泰勒级数展开

将两个展开式加和后可以得到

q(t+△t)等于qt加上v(t+1/2△t)

乘以△t

这表明△t时刻的位置

需要依据t时刻的位置

和1/2△t时刻的速度

也就是动量来求算

而对于1/2△t时刻

粒子的速度

也可以进行同样的处理

得到v（t＋1/2△t）

等于v（t－1/2△t）

加上at乘以△t

这表明1/2△t时刻的速度

依赖于t时刻的加速度

而加速度需要通过势能函数

也就是t时刻的

粒子的位置来进行求算

从这两个表达式我们可以看到

在动力学计算过程中

位置与动量的计算是不同步的

因此呢这个算法也被称为蛙跳算法

蛙跳算法能够得到粒子的速度

并且能够很方便的

对温度进行控制

它的精度也Verlet方法更高

但是值得注意的是

蛙跳算法在计算过程中

速度与位置是不同步的

它们之间相差了半个时间布长

另外由于计算中

对泰勒级数中的高阶相机性的阶段

也会一定程度上

给动力学过程带来一些不稳定

针对这些缺陷

有一些新的更稳定的改进的数值算法

可以将高阶项引入进来

比如龙格库塔方法

和预估－校正方法等等

我们在这里就不再进行介绍了

一般来说

通过分子动力学模拟

我们可以对大量的结构

和性质变化进行分析

例如我们可以通过

动力学模拟计算得到体系的相应性质

并与实验结果进行对比

从而对模拟的精度作出评估

我们也可以计算得到

实验上难以得到的热力学性质

通过动力学模拟

我们还可以观察分子构象的变化过程

根据动力学模拟中分子能量的变化

我们可以求得体系的内能

基于模拟过程中对能量和温度的统计

我们可以计算得到体系的热容

我们还可以计算体系的压力和温度

在压力表达式中

第一项是理想气体项

基于维里定律得到

而第二项考察了分子间作用

对压力的贡献

是实际气体校正项

温度则简单的通过

它的统计热力学定义式

针对体系中所有的粒子求算动能平均

并结合能量均衡原理就可以得到

下面我们来简单介绍一下

分子动力学模拟的基本过程

首先我们需要选择一个

合适的分子力场

在前面的章节中我们介绍过

由于分子力场是经验势能函数

其函数形式和参数

通过针对不同类型的分子

进行优化和拟合得到

每个力场都有一定的适用范围

因此选择合适的力场

是分子动力学模拟中

非常重要的一个步骤

例如模拟蛋白质分子体系

我们可以选择Amber力场

Charmm力场

opls力场等

专门针对蛋白质和DNA

这些生物大分子开发的力场

模拟有机分子时

可以选择MM3MM4力场

或者Compass力场等等

选定力场后

通常分子动力学模拟过程

可以分成三步

首先通过对要模拟的体系

进行结构优化

来获得合理的初始结构

一般来说

当我们模拟一些大分子时

初始结构往往通过实验手段

比如X射线衍射

核磁共振等方法得到

由于实验技术的限制

其中的原子间可能还存在一些

能量较高的相互作用

通过能量最小化

可以将这些作用消除

从而得到相对合理的初始结构

第二步

经过优化后的结果

需要通过平衡相计算

来使模拟体系处于

指定系综下的平衡态

在平衡相计算过程中

我们可以监测体系的热力学参数

例如温度压力密度等等

以及分子的结构参数

例如RMSD均方根位移等等

来确保体系达到了平衡

第三步经过平衡后的体系

就可以真正用来生成分子动力学轨迹

并用于下一步的进一步分析了

分子动力学计算与分子结构优化计算一样

也存在收敛问题

分子动力学收敛的标准是

动力学轨迹已经收集到了

足够的相点可以近似

认为相空间已经得到了各态历经

但是这一标准

难以在实际操作中得到确认

因此在实际计算中

往往会选择一些更能反映

体系统计性质的指标作为判据

例如在对蛋白质分子进行模拟时

可以监测分子

在动力学过程中

相对于其晶体结构的均方根偏差

RMSD值

RMSD值一方面可以判断

模拟得到的结构与实验结构之间的偏差

同时根据RMSD值的波动范围

以及绝对值的大小

也可以判断动力学过程的收敛性

一般来说

在一段已经达到收敛的动力学轨迹中

RMSD值的波动范围应该非常小

同时为了避免

体系在某些介稳状态附近变化

而造成收敛的假相

在动力学计算过程中

还可以通过改变初始结构的方式

来获得多条相关性比较小的轨迹

来进行对比分析

在分子动力学模拟的过程中

有一个非常重要的问题就是

怎样仅对少数分子进行计算

来模拟包含有大量分子的宏观系统

显然仅仅针对数个水分子进行模拟

无法代表宏观状态下

一升水的真实情况

这是因为随着分子数量的减少

边界效应开始变得非常大

在模拟时应该对边界效应

做出适当的处理

周期性边界条件解决了这个问题

通过周期性边界条件的使用

人们能够基于对少量粒子进行模拟

来研究一个宏观体系的性质

周期性边界条件

将宏观体系

设想成一个个微观盒子的

周期性重复

而计算模拟

则仅需要对其中一个盒子

来完成就可以了

以图中的假想二维体系为例

每一个正方形的二维盒子

被八个相同的盒子包围

这样一个盒子里的粒子的行为

可以用来模拟整个宏观体系

而避免了因为使用较小的体系

而带来的边界效应

在模拟过程中

如果有分子从盒子的一侧离开

由于周期性设定

另外一侧会有同样的分子进入

从而也保证了体系粒子数的恒定

在三维空间中

通过同样的设定

每个盒子被26个相邻的盒子所包围

三维空间中采用立方体盒子

是最简单的边界条件

然而有时候也可以根据分子形状

从节省计算资源的角度

选择更合适的周期性边界条件

在三维空间中

还可以选择截角八面体

六棱柱，菱形12面体等

其他的周期性边界条件

例如对于接近于球形的分子

使用截角八面体边界条件时

所需要的溶剂分子就会更少

因此相比立方体边界条件

就减少了计算量

使用周期性边界条件

有一个需要注意的问题是

要避免不同盒子中

同一个原子之间产生作用

因为不同盒子互为映像

相同一个原子跟自己的映像

是不应该产生作用的

因此盒子大小需要跟势场的作用范围

尤其是非键作用的作用距离相互匹配

在实际计算中

可以通过对静电作用和范德华作用

设置截断距离来避免这一问题

非键作用的截断距离

一般应该小于周期性盒子的周期大小

针对非键作用设置截断距离

也是为了减少计算量

因为非键作用存在于体系中

任意一对原子之间

所以它的计算量

正比原子数N的平方

但是实际相距较远的原子间的作用

已经弱到可以忽略不计

因此在动力学模拟计算时

可以通过设定截断距离近似认为

超过截断距离的原子

对之间的相互作用等于0

从而减少计算量

通常为了避免

设定截断距离后

带来的势能函数的不连续

以及在截断距离附近

原子受力的波动

还可以通过开关函数

来对势能函数进行进一步的调整

原则上来说

由于静电作用的衰减

远小于范德华作用的衰减

因此针对静电作用

和范德华作用应该设置

不同的截断距离

通常我们一个我们

感兴趣的变化过程的时间长度

往往在微秒甚至毫秒级别

如果不经过任何处理

对这样的过程进行模拟计算

它的计算量将会过大

导致实际上无法完成

因此近些年来

发展了一些加速分子动力学计算的方法

它们包括通过周期性的

对体系进行加热退火

来避免模拟体系陷入局域极小

或者通过增加一定的偏移势

来引导动力学过程

确保针对感兴趣的结构区间进行采样

包括SMD TMD

Metadynamics等方法等等

另外副本交换

伞状取样等策略

也在现在的分子动力学计算中

经常采用

我们在这里就不再继续介绍了

今天的课我们就讲到这里

谢谢大家