分 子 模 拟
牛继南
njn0516@cumt.edu.cn
2011.3
第玖章 限制性分子动力学
虽然分子动力模拟的计算能力随着计算
机速度日益剧增，但通常能计算只有数纳
秒内上万个分子体系的变化，与实际的体
系相距还很远，因此，改进计算方法为分
子动力模拟中最重要的问题之一。
9.1 刚性线性分子的计算
• 执行分子动力学计算时，若系统的热力学
性质或宏观性质为主要的研究对象，则可
以忽略分子内部的相对运动。
• 通常采用刚性分子计算法，即固定分子中
原子的相对位置，移动时分子作为一个整
体。
• 通过这种方法可以增加计算的效率，并可
采用较大的步长。
平
动
转
动
• 分子质心的坐标为：
其中n为原子个数，ma和 分别为a原子的质
量和其位置矢量，M为分子的总质量。
• 对于线性分子，如HCN和CO2，其角速度
和力矩 必定垂直其分子轴，设 为分子轴
向的单位矢量，则：
• 其中 取决于分子间作用力。
具有以下特性：
可以依此来求分子
力矩平方的平均值。
• 分子中原子相对于质心的位置为：
其中 为a原子相对于质心的距离。
• 则：
• 为原子a所受的力。式子中的 仅有垂直
于分子轴的分量 起作用，即：
• 对于 ，可由下式求出：
• 线性分子的转动可写成：
其中 为角速度， 为分子的转动惯量。
为拉格朗日未定乘子（Lagrange
undetermined multiplier）。由于
和 ，以LF法解此微分式可得:
• 分子的移动即为质心的移动，而质心所受
的力为分子中原子所受力的总和，即：
• 由LF法解得的质心运动方程为：
• 式中M为分子的总质量。
9.2 刚性非线性分子的计算
• 对于非线性刚性分子，造成分子转动的力
矩为：
• 分子的位向定义了分子固定轴和空间固定
轴之间的关系。一般分子固定轴的选取是
使分子的转动惯量矩阵成对角化状态。
• 设空间中一个单位向量在分子轴坐标系中
的位置是 ，在空间坐标系中的位置是
，两者之间的关系为：
• 式中，转置矩阵和欧拉角欧拉角欧拉角欧拉角(Ruler angle)
相关。
两坐标系
O-ξ η ζ
O-x y z
交线
ON
称为节线
ON-Oξ φ 进动角
ON-Ox ψ 自转角
Oξ -Oz θ 章动角
• 由欧拉角确定的矩阵 为：
• 得到分子坐标轴的运动方程为：
• 式中 为分子固定轴的3个转动惯
量。
• 分子固定轴和空间固定轴中力矩和角速度
的关系为：
• 式子中，
• 可得欧拉角的运动方程：
θ趋近于0时，分
母有很大的误差
• Evans和Murad发展出了新的计算法，其定
义 为含有4个分量的向量，即：
• 四个分量满足：
• 定义它们与欧拉角的关系为：
• 将转置矩阵 用q来表示，则：
• 每一个分子的 向量满足下面的运动方程：
• 以LF求出角动量的运动方程：
• 可以利用
• 求得分子轴的角速度：
• 带入前面的式子还可求得：
这就改进了直接用欧拉
角求解的缺点。
• 利用刚性分子动力计算方法可以节省很多
时间，如水分子系统中，由于羟基伸缩振
动频率很大约为3660cm-1，因此仅能采用约
为5×10-16s的积分步长，但若采用刚性分子
可增大至6×10-15s，即增加约10倍的计算效
率。
• 因为在大多数化学体系中都是以水溶液为
基础，所以采用刚性水分子模型能大大节
省实验时间。
9.3 键长限制法
• 在分子的各种内部运动中，以键伸缩和键
角弯曲速度最快，但多数分子的键长和键
角仅在平衡位置附近呈小幅度振动，键长
的变化约为3%，键角的变化约为5~10°。
• 实际当中分子动力计算的一些性质与键长
或键角的变化并没有太大关系，因此在计
算中固定键长和键角，可以采用较长的积
分步长，并获得更多的动力学资料。
• 将键长和键角固定的动力计算也是限制计
算法中的一种。
• 最通常的键长限制法由Ryckaert等
设计
领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计
，这
种方法是在原来的势能中加入一项限制。
设欲限制键长在dij，则势能为：
• 其中U0为分子的一般势能项，第二项为键
长限制项，λij为拉格朗日未定乘子。
• 原子受到的力为：
• 式中F
0
为无限制时原子所受的力。
• 按Verlet法得到的运动方程为：
• 根据键长限制，无论何时
• 因此
• 又根据梯度关系式：
-
得到：
• 带入
得到:
其中 表示δt的四次方，可以忽略，因此：
根据此式可以解出λ，可以进而求出t+δt时
原子的位置。
9.4 限制键长和键角的计算
• 以非线性三原子分子为例，将连接1、2、3
原子的键长r12、r23和键角θ123同时限制住，
限制项为：
θ
按照前面的算法可导出下
面的关系式：
• 类似根据
得方程组：
• 将其写成矩阵的形式：
• 可以通过逆矩阵求得各个λ值，带入前面的
式子即可求得t+δt时原子的位置。
• 上面的矩阵除了用通常的方式求解外，还
可以用迭代计算法，如SHAKE方法。这种
方法中将矩阵中所有非对角元素忽略，然
后求出λ12、λ23和λ13，进而求出r(t+δt)后，检
查新原子位置处的键长是否满足限制条件
， 如果不满足，将r(t+δ
t)视为新的起点计算新的λ和新的键长，知
道满足为止。通常要经过7~10次左右才能
找到原子的正确位置，但不需求逆矩阵，
因而适用于限制较多的情况。
• 除了SHAKE法以外，还有RATTLE算法。
9.5 聚合物的限制动力学计算
• Toxvard将有机正葵烷分子的甲基和亚甲基
均视为粒子（称为联合原子模型），每个
分子含有10个这样的粒子，计算时分为全无
限制情况和限制情况，限制时分为两类。
一是限制所有键长，二是限制所有键长和
键角。
• 系统的密度定为0.63g/ml,温度定为481K，全
无限制的计算时间为20ns，加入限制时为
6ns，动力计算采用LF方法，积分步骤δt为
2fs。
• 计算时，对聚合物比较重要的性质有：
1. 末端距离(end-to-end distance)平方的平均值
，R为链分子头端基团与末端基团的距
离;
2. 回旋半径(gyration radius)，回旋半径的
计算
公式
小学单位换算公式大全免费下载公式下载行测公式大全下载excel公式下载逻辑回归公式下载
为：
mi为每个粒子的质量，di为各粒
子距质心的距离。
3.此外，还需计算聚合物的转动惯量：
式中x，y，z为粒子的笛卡尔坐标，通常还
可以选取另外的坐标系使得矩阵对角化。
• 即
Iaa，Ibb和Icc相当于绕着a，b，c轴旋转的转
动惯量。
• 由分子模拟得到的末端距离平方、回转半
径和转动惯量可以反映聚合物的形状。
• 对于含N根聚合物的系统，任何物理量A的
平均值为：
T为计算的总时间，Ai,t为
第j根链在t时刻的A值。
• 从图上可以看出，若
仅限制键长，则计算
的扭转角分布与未限
制系统的结果十分相
近，若加上键角的限
制，则误差较大，但
仍处在相似的分布。
• 从图中可以看出增
加键角限制也会导
致较大的末端距离
误差。
• 聚合物中每四个依次相邻的C链中（C-C-C-
C）扭转角均有邻位交叉式（gauche）和反
式（trans）构象。
• 邻位交叉式有两种，g+和g-，反式则为最稳
定的构型。
• 如果用ng和nt表示分子中邻位和反式构象的
含量，其平均值的比为/，此比值
的大小反映了聚合物的稳定性，值越大表
明结构越不稳定，越小反式构型越多，聚
合物也越稳定。
• 从表6-2中可以看出 均小于
• ，这表示聚合物的平均形状略显扁盘
形。
• 从上面的分析可知：
1.限制系统和无限制系统的变化趋势是一致
的；仅限制键长系统的计算结果无论是从
分子形状、角度分布、物理性质等都非常
接近无限制系统的结果；但加入键长后可
允许分子动力计算中引入更长的积分步骤
，节省计算时间和计算更长的运动轨迹。
2.因此在长时间、大体系中宜采用限制计算
，尤其是键长限制法。因为对复杂体系，
同时限制键长和键角并不容易实现，需要
更好的数学方法才能解决。
浙公网安备 33021202002483