分子动力学模拟及相关研究

1、分子动力学模拟及相关研究分子动力学模拟及相关研究MolecularMolecular Dynamics Dynamics Simulation and Related Research Simulation and Related Research 报告人：P.Li2021/10/131.粒子的运动服从Newton运动方程2.动力学机制决定运动方程的形式3.适用范围广泛1.分子运动方程及其数值求解 2.多体系统的基本概念与分子动力学方法3.分子力场1.模拟模型的设定2.给定初始条件3.趋于平衡 4.宏观物理量的计算 一一、分子动力学特点分子动力学特点二二、分子动力学基础知识、分子动力学基础知识

2、三、三、分子动力学模拟步骤及相关研究分子动力学模拟步骤及相关研究一一. . 特点特点1.1.在在分子动力学中，粒子的运动行为是通过经典的分子动力学中，粒子的运动行为是通过经典的NewtonNewton运动方程所描述。运动方程所描述。分子动力学是在原子、分子水平上求解多体问题的重要的计算机模拟方法，可以预测纳米尺度上的材料动力学特性。通过求解所有粒子的运动方程，分子动力学方法可以用于模拟与原子运动路径相关的基本过程。分子动力学模拟方法可以看作是体系在一段时间内的发展过程的模拟。在这样的处理过程中我们可以看出：分子动力学方法中不存在任何随机因素。在分子动力学方法处理过程中，方程组的建立是通过对物理

3、体系的微观数学描述给出的。在这个微观的物理体系中，每个分子都各自服从经典的牛顿力学。每个分子运动是用理论力学上的哈密顿量或者拉格朗日量来描述，也可以直接用牛顿运动方程来描述。2.系统系统的动力学机制决定运动方程的形式的动力学机制决定运动方程的形式这种方法可以处理与时间有关的过程，因而可以处理非平衡态问题。但是使用该方法的程序较复杂，计算量大，占内存也多。3.适用范围适用范围广泛广泛：原则上，分子动力学方法所适用的微观物理体系并无什么限制。这个方法适用的体系既可以是少体系统，也可以是多体系统；既可以是点粒子体系，也可以是具有内部结构的体系；处理的微观客体既可以是分子，也可以是其它的微观粒子。自五

4、十年代中期开始，分子动力学方法得到了广泛的应用。它与蒙特卡洛方法一起已经成为计算机模拟的重要方法。应用分子动力学方法取得了许多重要成果，例如气体或液体的状态方程、相变问题、吸附问题等，以及非平衡过程的研究。二. 分子动力学基础知识1、分子运动方程及其数值求解分子运动方程及其数值求解 采用分子动力学方法时，必须对一组分子运动微分方程做数值求解。从计算数学的角度来看，这是个求一个初值问题的微分方程的解。实际上计算数学为了求解这种问题已经发展了许多的算法，但是并不是所有的这些算法都可以用来解决物理问题。以一个一维谐振子为例，来看一下如何用计算机数值计算方法求解初值问题。一维谐振子的经典哈密顿量为：2

5、2212kxmpH这里的哈密顿量（即能量）为守恒量。假定初始条件为x（0）P（0），则它的哈密顿方程是对时间的一阶微分方程 mpPHdtdxkxxHdtdp计算在相空间中的运动轨迹(x(t),p(t)：采用有限差分法，将微分方程变为有限差分方程，以便在计算机上做数值求解，并得到空间坐标和动量随时间的演化关系。首先，取差分计算的时间步长为h，采用一阶微分形式的向前差商表示，直接运用展开到h的一阶泰勒展开公式： )()()(2hOdtdfhtfhtf即可得到：htfhtfdtdf)()(则微分方程改写为差分形式：mtphtxhtxdtdx)()()()()()(tkxhtphtpdtdp我们得到解

6、微分方程的欧拉(Euler)算法：)()()()()()(thkxtphtpmthptxhtx这是一组递推公式。有了初始条件，就可以一步一步地使用前一时刻的坐标、动量值确定下一时刻的坐标、动量值。这个方法是一步法的典型例子。 由于在实际数值计算时的大小是有限的，因而在上述算法中微分被离散化为差分形式来计算时总是有误差的。可以证明一步法的局部离散化误差与总体误差是相等的，都为的 量级。在实际应用中，适当地选择的大小是十分重要的。取得太大，得到的结果偏离也大，甚至于连能量都不守恒；取得太小，有可能结果仍然不够好。这就要求我们改进计算方法，进一步考虑二步法。)(2hO差分计算的二步法实际上泰勒展开式

7、的一般形式为：)()(!)()(1)(niniihOtfihtfhtf前面叙述的欧拉算法就是取n=1。取n=2 展开后可得到： )()(2)()()()(2)(2htxhtxhmtpmtFhhtxtxhtx这是二步法的一种, 称为Verlet方法。Verlet算法是分子动力学模拟中求解常微分方程最通用的方法2、多体系统的基本概念与分子动力学方法、多体系统的基本概念与分子动力学方法N体系统中，一个n体的密度函数一般可以写为：NnnnnrdrdrdRWnNNZrrr2121)()!(!1),(其中 是描写系统的几率函数。)(RW 通常为由系统中所有粒子的坐标、动量构成的相空间中的任意一点。R分子动

8、力学元胞分子动力学模拟方法往往用于研究大块物质在给定密度下的性质，而实际计算模拟不可能在几乎是无穷大的系统中进行。所以必须引进一个叫做分子动力学元胞的体积元, 以维持一个恒定的密度。对气体和液体，如果所占体积足够大，并且系统处于热平衡状态的情况下，那么这个体积的形状是无关紧要的。对于晶态的系统，元胞的形状是有影响的。为了计算简便，对于气体和液体，我们取一个立方形的体积为分子动力学元胞。设分子动力学元胞的线度大小为L，则其体积为L3。由于引进这样的立方体箱子，将产生六个我们不希望出现的表面。模拟中碰撞这些箱的表面的粒子应当被反射回到元胞内部，特别是对粒子数目很少的系统。然而这些表面的存在对系统的

9、任何一种性质都会有重大的影响。元胞的周期性边界条件为了将分子动力学元胞有限立方体内的模拟，扩展到真实大系统的模拟，我们通常采用周期性边界条件。采用这种边界条件, 我们就可以消除引入元胞后的表面效应，构造出一个准无穷大的体积来更精确地代表宏观系统。实际上，这里我们做了一个假定，即让这个小体积元胞镶嵌在一个无穷大的大块物质之中。周期性边界条件的数学表示形式为：),(),()(321nnnnLnxAxA其中 为任意的可观测量 ， 为任意整数。这个边界条件就是命令基本分子动力学元胞完全等同地重复无穷多次。A321,nnn该边界条件的具体实现是这样操作的：当有一个粒子穿过基本分子动力学元胞的六方体表面时

10、，就让这个粒子以相同的速度穿过此表面对面的表面重新进入分子动力学元胞内。不同分子动力学元胞盒子内粒子间的相互作用：对于不同分子动力学元胞盒子内粒子间的相互作用，如果相互作用是短程力，我们可以在长度rc处截断。这里V (rc)必须要足够小，以使截断不会显著地影响模拟结果。cijcijcijijSrrrrrr0-)()(VVV势函数直接截断势函数直接截断 典型的分子动力学元胞尺度L通常选得比rc大很多。我们往往选择元胞尺度满足不等式条件L/2rc ，使得距离大于L/2 的粒子的相互作用可以忽略，以避免有限尺寸效应。通常L的数值应当选得很大。在考虑粒子间的相互作用时，通常采用最小像力约定。最小像力约

11、定是在由无穷重复的分子动力学基本元胞中，每一个粒子只同它所在的基本元胞内的另外N -1个中（设在此元胞内有N 个粒子）的每个粒子或其最邻近的影像粒子发生相互作用。 分子力学，又叫力场方法（ force field method），是基于经典牛顿力学方程的一种计算分子的平衡结构和能量的方法。与量子力学不同，它求解的是Newton方程，而不是薛定谔方程。 基本假设：.The Born-Oppenheimer Approximation对势能面的经验性拟合。 量子力学中的薛定谔方程为:(非相对论和无时间依赖的情况下） )()(n);();(e),(),(RERRrERrrRErR体系的哈密顿算符体系

12、的哈密顿算符与原子核（与原子核（R）和电子（）和电子（r）位置相关的波函数）位置相关的波函数分子力学分子力学-简介简介3.3.分子力场分子力场基于Born-Oppenheimer近似，其物理模型可描述为：原子核的质量是电子质量的103105倍，电子速度远远大于原子核的运动速度，每当核的分布形式发生微小变更，电子立刻调整其运动状态以适应新的核场。这意味着，在任一确定的核分布形式下，电子都有相应的运动状态；同时核间的相对运动可视为所有电子运动的平均结果。所以电子的波函数只依赖于原子核的位置，而不是他们的动能。于是这个近似认为，电子的运动与原子核的运动可以分开处理，可以将上式分解为电子运动的波函数核

13、运动的波函数电子运动方程：电子运动方程：核运动方程：核运动方程：基本假设基本假设2()( , )( )( , )( )()( , )( , )()( )(0)0elNNeNeNelNNekeelekeelNNTHVR rRER rRHVR rER rTERERl方程中的能量方程中的能量Ee(Ee(势能面势能面) )仅仅是原子核坐标有关。相应的，方程仅仅是原子核坐标有关。相应的，方程所表示的为在核势能面所表示的为在核势能面E(R)E(R)上的核运动方程。上的核运动方程。l直接求解方程，采用的是从头算或者是半经验，这样的量化计算都直接求解方程，采用的是从头算或者是半经验，这样的量化计算都是把电子的

14、波函数和能量处理成原子核坐标的函数。由于量子化学求是把电子的波函数和能量处理成原子核坐标的函数。由于量子化学求解电子波函数和势能面耗时巨大，常常将势能面进行经验性的拟合，解电子波函数和势能面耗时巨大，常常将势能面进行经验性的拟合，成为力场，由此构成分子力学的基础。成为力场，由此构成分子力学的基础。l将方程用牛顿运动方程代替，势能面采用力场拟合，就构成了分子将方程用牛顿运动方程代替，势能面采用力场拟合，就构成了分子动力学的基础。动力学的基础。2()(, )()(, )()()(, )(, )()()(0)0elNNeNeNelNNekeelekeelNNTHVR rRER rRHVR rER r

15、TERER基本假设基本假设电子运动方程：电子运动方程：核运动方程：核运动方程： . .简单简单作用模型作用模型 对体系相互作用的贡献来自诸如键伸缩、键角的开合、单键的旋转等等。即使对体系相互作用的贡献来自诸如键伸缩、键角的开合、单键的旋转等等。即使使用类似使用类似HookeHooke定律这样的简单函数，也能令力场运转良好。定律这样的简单函数，也能令力场运转良好。 . .力场力场的可移植性（关键属性）的可移植性（关键属性） 仅在少数情况下通过测试的一套函数，可以用来解决更广范围内的问题。进一仅在少数情况下通过测试的一套函数，可以用来解决更广范围内的问题。进一步讲，从小分子得来的数据可以用来研究类

16、似高分子的大分子。步讲，从小分子得来的数据可以用来研究类似高分子的大分子。 简单分子力场简单分子力场 分子力场是分子力学的核心。分子力学的基本理论就是一个分子力场由分子内相分子力场是分子力学的核心。分子力学的基本理论就是一个分子力场由分子内相互作用和分子间相互作用两大部分构成，即力场的势能包括成键和非键相互作用，所互作用和分子间相互作用两大部分构成，即力场的势能包括成键和非键相互作用，所有的势能的总和即为分子的构象能。有的势能的总和即为分子的构象能。基本假设基本假设力场力场分类分类 力场分类规则 定义：把键长、键角、二面角分成不同类型，同种类型的立场参数相同 分类规则总是对原子进行。理论上说，

17、如果一个力场有N个原子类型，那么就有N(N-1)/2种键， N(N-1)(N-2)/2种键角， N(N-1)(N-2)(N-3)/2种二面角 举例：乙烷分子 1个C-C键，6个C-H键 6个C-C-H键角， 6个H-C-H键角 9个H-C-C-H二面角传统力场传统力场 AMBER力场 Kollman group, 1984 最初仅为蛋白质和核酸体系提供相应的原子类型和力场参数 1990，发展了适用于多糖模拟的力场参数（Homan 1990） 1995, 2000，加入了适用于有机小分子的原子类型和参数22001260()()1 cos()22244nbondsanglestorsionsiji

18、jijijijijijijVkkEllnq qrrr传统力场传统力场 CHARMm力场 （Chemistry at Harvard Macromolecular mechanics） Karplus group，1983 适用于各种分子性质的计算和模拟 对于从孤立的小分子到溶剂化的大生物体系的多种模拟体系都可以给出较好的结果 但不适合于有机金属配合物2200222201260()()cos()()(,)4bbijijijijonoffijijijijijEk rrkkknq qABksw rrrrrr传统力场传统力场 CVFF力场 （Consistent Valence Force Field

19、） Dauber-Osguthorpe group, 1988 适用范围包括有机小分子和蛋白质体系 扩展后可用于某些无机体系的模拟，如硅酸盐、铝硅酸盐、磷铝化合物 主要用于预测分子的结构和结合自由能0()220000000001261()1 cos()()()()()()()cos ()()2jb bbbbbbbbbijijEDeHHnHFbbbbFFbFFrrrr ijijijijq qr传统力场传统力场 MMX力场 Allinger group, 1989 包括MM2和MM3 主要针对有机小分子 函数形式比较复杂，包含交叉项 也可用于生物大分子体系，但是速度会比较慢sBTOOPSBVQEE

20、EEEEEEE传统力场传统力场第二代力场第二代力场 CFF力场 CFF91主要适用于模拟有机小分子、蛋白质以及小分子-蛋白质之间的相互作用 CFF95除了适用于蛋白质和有机小分子体系，还可用于有机高分子体系的模拟，如聚碳酸酯及多糖 PCFF在CFF91的基础上，还适用于聚碳酸酯、三聚氰胺甲醛树脂、多糖、核酸、分子筛等其他无机和有机材料体系的模拟。同时还提供20种金属离子的参数第二代力场第二代力场 COMPASS力场 （Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies用于原子水平模拟研究

21、的凝聚态优化的分子力场） Sun group, 1994 第一个能同时预测气态和凝聚态性质的从头力场 适用于常见的有机小分子、无机小分子和高分子 最新的力场也适合于金属、金属氧化物及金属卤化物的模拟 可较为准确的预测晶体的多种性质第二代力场第二代力场 MMFF94力场 （Merk Molecular Force Field） Hagler, 1996 目前最为准确的力场之一 定义了非常完备的原子类型 既适用于有机小分子，也适用于大分子体系，如蛋白质sBTOOPSBVQEEEEEEEEE第二代力场第二代力场通用力场通用力场 基于规则的力场 （Rule-based Force Fields） 较好

22、的可移植性，适用范围广 精度差 ESFF力场 （Extensive Systematic Force Field） 用于有机分子、无机分子及有机金属化合物的结构预测 UFF力场 （Universal Force Field） 覆盖了周期表中所有元素，应用最广泛的通用力场 Dreiding力场 适用于有机分子、生物大分子及位于主族的无机分子通用力场通用力场分子力场的选择分子力场的选择 蛋白质分子的模拟：首选AMBER力场和CHARMm力场，也可用CFF力场、CVFF力场和MMFF94力场 核酸分子的模拟：采用AMBER力场、CHARMm力场、MMFF94力场或用户自定义的力场 小分子-蛋白质复合

23、物体系的模拟：首选CHARMm力场和MMFF94力场，也可用CVFF力场和CFF力场 高分子的模拟：首选COMPASS力场，也可用PCFF力场和CFF95力场分子力场的选择分子力场的选择三.分子动力学模拟的基本步骤在计算机上对分子系统的分子动力学模拟的实际步骤可以划分为四步：首先是设定模拟所采用的模型设定模拟所采用的模型；第二，给定初始条件给定初始条件；第三，趋于平衡的计算趋于平衡的计算过程；最后是宏观物理量的计算宏观物理量的计算。1模拟模型的设定设定模型是分子动力学模拟的第一步工作。例如在一个分子系统中，假定两个分子间的相互作用势为硬球势，实际上，更常用的是Lennard-Jones型势。它

24、的势函数表示为:612)()(4)(rrrU其中是位势的最小值，r=时位势为零.第 i 个粒子与元胞内其它N-1个粒子或其最邻近的影像粒子的相互作用力在 x 方向的分量为:8142,)(21)()(48ijijNijjixirrxxFY 和Z分量的表达式与上式相同模型确定后,根据经典物理学的规律我们就可以知道在系综模拟中的守恒量。例如对在微正则系综的模拟中能量、动量和角动量均为守恒量。为了计算方便，取分子动力学元胞为立方形，元胞的线度大小为L。iiiiiiprMpPrVrmE)()(212例如，Verlet方法需要两组坐标来启动计算：一组是零时刻的坐标，另一组是前进一个时间步长时的坐标，或者是

25、一组零时刻的速度值。2给定初始条件分子动力学模拟的过程进入对系统微分方程组做数值求解时，需要知道粒子的初始位置和速度的数值。不同的算法要求不同的初始条件。但是，一般来说系统的初始条件都是不可能知道的。表面上看这是一个难题。实际上，精确选择待求系统的初始条件是没有什么意义的，因为模拟时间足够长时，系统就会忘掉初始条件。但是初始条件的合理选择将可以加快系统趋于平衡。 常用的初始条件可以选择为：（1）令初始位置在差分划分网格的格子上，初始速度则从玻尔兹曼分布随机抽样得到。（2）令初始位置随机地偏离差分划分网格的格子，初始速度为零。3趋于平衡 按照上面给出的运动方程、边界条件和初始条件，就可以进行分子

26、动力学模拟计算。但是，这样计算出的系统不会具有所要求的系统能量，并且这个状态本身也还不是一个平衡态。为了使系统达到平衡，模拟中需要一个趋衡过程。在这个过程中，我们增加或从系统中移出能量，直到系统具有所要求的能量。然后，再对运动方程中的时间向前积分若干步，使系统持续给出确定能量值。我们称这时系统已经达到平衡态。这段达到平衡所需的时间称为弛豫时间。在分子动力学模拟中，时间步长h的大小选择是十分重要的。为了减小误差，步长必须取得小一些；但是取得太小，系统模拟的弛豫时间就很长。这里需要积累一定的模拟经验，选择适当的时间步长h。4宏观物理量的计算 实际计算宏观物理量往往是在分子动力学模拟的最后阶段进行的

27、。它是沿着相空间轨迹求平均来计算得到的。例如对于一个宏观物理量A，它的测量值为平均值如果已知初始位置和动量为 和 (上标N表示系综N个粒子的对应坐标和动量参数），选择某种分子动力学算法求解具有初值问题的运动方程，便得到相空间轨迹:A)0()(Nr)0()(Np)(),()()(tptrNN对轨迹平均的宏观物理量A的表示为 :0)()()(),()(1limttNNtprAdttA如果宏观物理量为动能，它的平均为 :0)()()(1limttNktkpEdttE分子动力学相关软件CHARMm:由哈佛大学开发的分子力学计算工具。可以为用户提供处理各种小分子、大分子(包括蛋白质、核酸和糖)的经验化能量计算,模拟过程提供了有关分子结构、相互作用、能量等信息 。lAmber：用于分子（特别是生物大分子）的动力学模拟以及小