分子动力学方法

分子动力学方法3.1基本原理计算机模拟分类：(1)随机模拟方法。优点:随机模拟方法计算的程序简单，占内存少，但是该方法难于处理非平衡态的问题。(2)分子动力学方法（MolecularDynamics或简称MD）。可以处理非平衡态问题。但是使用该方法的程序较复杂，计算量大，占内存也多。分子动力学方法利用牛顿古典力学来计算许多分子在相空间中的轨迹。3.1基本原理分子动力学(MD)方法的出发点是物理系统的确定的微观描述（哈密顿描述方程、拉格朗日方程或牛顿运动方程）。因此，分子动力学方法是用运动方程来计算多体或少体系统的性质，结果得到的既有系统的静态特性，也有动态特性。蒙特卡罗方法只能得到系统的位形特性。MD方法的具体做法是在计算机上求运动力程的数值解。为此，通过适当的格式对方程进行近似，使之适于在计算机上求数值解。其实质是计算一组分子的相空间轨道，其中每个分子各自都服从经典运动定律。这里的系统不仅是点粒子系统，也包括具有内部结构的粒子组成的系统。3.1基本原理1.分子动力学是在原子、分子水平上求解多体问题的重要的计算机模拟方法。2.通过求解所有粒子的运动方程，分子动力学方法可以用于模拟与粒子运动路径相关的基本过程。3.在分子动力学中，粒子的运动行为是通过经典的Newton运动方程所描述。3.1.1粒子运动方程的数值求解粒子体系的运动方程Lagrangian方程1.Lagrangian函数的定义为则运动的Lagrangian方程为0ddkkqqtLL3.1.1粒子运动方程的数值求解粒子的运动方程-Lagrangian方程kkkkkkkkkqHppHqLpqHqLp),(),(qqqp3.1.1粒子运动方程的数值求解单原子的牛顿运动方程V...iiiiiiirmmrrfp3.1.1粒子运动方程的数值求解Verlet算法r(t+∆t)=r(t)+v(t)∆t+(1/2)a(t)∆t2(1)r(t-∆t)=r(t)–v(t)∆t+(1/2)a(t)∆t2(2)将上面两式相加，得到：r(t+∆t)=2r(t)–r(t-∆t)+a(t)∆t2(3)v(t+∆t)=v(t)+a(t)∆t+(1/2)b(t)∆t2(4)a(t+∆t)=a(t)+b(t)∆t(5)将(5)式的b(t)代入(4)，得到：v(t+∆t)=v(t)+(1/2)[a(t)+a(t+∆t)]∆t(6)3.1.1粒子运动方程的数值求解其他求解算法：Leap-frogalgorithmr(t+∆t)=r(t)+v(t+(1/2)∆t)∆tv(t+(1/2)∆t)=v(t-(1/2)∆t)+a(t)∆tBeeman’salgorithmr(t+∆t)=r(t)+v(t)∆t+(2/3)a(t)∆t2–(1/6)a(t-∆t)∆t2v(t+∆t)=v(t)+v(t)∆t+(1/3)a(t)∆t+(5/6)a(t)∆t–(1/6)a(t∆t)∆t3.1.2热力学量的计算在物理系统的计算机模拟中，系综平均必须用时间平均代替，在通常的模拟中，粒子数N和体积V是固定的。给定初始位置rN(0)和初始动量pN(0)后，一个MD算法将从运动方程生成轨道(rN(t),pN(t))，轨道平均的定义为假定能量守恒，并且轨道在一切具有同一能量的相同休积内经历相同的时间，则轨道平均等于微正则系综平均'10'0))();()(()'(limttNNttVtptrdtAttANVEAA3.1.2热力学量的计算孤立系统的总能量是一个守恒量，沿着分子动力学模拟生成的任何一条轨道，能量应保持不变，即。孤立系统的动能和势能不是守恒量，它们的大小沿着生成的轨道逐点变化EE'10''10'00,))(()'(lim,))(()'(limtttttktkdttrUttUdttEttEnnvvkkEnnE0013.1.2热力学量的计算生成的动能路径是不连续的，必须在各个时间间断点上计算动能的值以求平均:其中从平均动能可以计算系统的温度，温度是一个重要的物理量，需要加以监测，特别是在模拟的起始阶段。nnvvkkEnnE001ivivkmE)(2123.1.2能量均分定理在热力学极限下，一切系综都是等同的，并且可以应用能量均分定理。热力学极限是指粒子数（或体积）趋向无穷大时的极限。一般宏观物体包含了1023个粒子，可以认为是满足热力学极限的。相变只有在热力学极限下才会发生。系统的哈密顿量为则：ijiijiirumP221TkmBi212123.1.2能量均分定理由于系统的每个粒子有三个自由度，因此假定位势在处被截断。系统内部的位形能量的轨道平均值为其中TNkEBk23nnvvUnnU001jivijvruU)(3.1.2能量均分定理由于位势被截断，总能量和势能含有误差，为了估计必须作出改正。所有的内部位形能都加到截止距离为止，尾部改正可以取其中是g(r)对关联函数，它是在原点r=0处有一个粒子时,在r周围的体积元dr内找到一个粒子的概率，令n(r)为离一个给定粒子的距离在r和r+dr之间的平均粒子数，于是crcdrrrgruU2)()(2rrrnNVrg24)()(3.1.2能量均分定理其他的量也需要尾部改正，以压强的计算为例，这时维里(virial)状态方程成立。至于势能的计算，可以把积分分成两项，一项是由相互作用力程之内的贡献引起的，一项是对位势截断的改正项：长程改正项为0324)(6drrrurgTPkBCjiijijkBPrurNTP6crCdrrrurgP324)(63.2分子动力学模拟的基本步骤分子动力学模拟的实际步骤可以划分为四步：设定模拟所采用的模型；给定初始条件；趋于平衡的计算过程；宏观物理量的计算。3.2分子动力学模拟的基本步骤1．模拟模型的设定硬球势Lennard-Jones型势.r,r如果如果,0,rU6124rrrU0.60.81.01.21.41.61.82.02.22.42.6-3-2-101234位势V(r)力F(r)吸引力排斥力根据经典物理学的规律，可以知道在系综模拟中的守恒量。微正则系综的模拟中能量、动量和角动量均为守恒量。在此系综中他们分别表示为：元胞周期性边界条件,iiirVrmE221iipPiiiprM3.2分子动力学模拟的基本步骤分子动力学模拟的最小像力约定示意图Lnrrminrjiij最小像力约定2.给定初始条件给定粒子的初始位置和速度的数值：（1）令初始位置在差分网格格子上，初始速度从玻尔兹曼分布随机抽样得到。（2）令初始位置随机地偏离差分网格格子，初始速度为零。（3）令初始位置随机地偏离差分网格格子，初始速度从玻尔兹曼分布随机抽样得到。3．趋于平衡使系统达到平衡，模拟中需要一个趋衡过程。在这个过程中，我们增加或从系统中移出能量，直到系统具有所要求的能量。3.2分子动力学模拟的基本步骤3.3平衡态分子动力学模拟1.微正则系综的分子动力学模拟粒子数恒定、体积恒定、能量恒定、整个系统的总动量恒等于零。分子动力学模拟步骤如下（Verlet算法）：（1）给定初始空间位置。（2）计算在第n步时粒子所受的力。（3）计算在第n+1步时所有粒子所处的空间位置。（4）计算第n步的速度。（5）返回到步骤（2），开始下一步的模拟计算。mhFrrrnininini/22)()1()()1(hrrvninini2/)1()1()(Verlet算法的速度形式：（1）给定初始空间位置。（2）给定初始速度。（3）利用公式：计算在第n+1步时所有粒子所处的空间位置。（4）计算在第n+1步时所有粒子的速度：（5）返回到步骤（3），开始第n+2步的模拟计算。)1(iv)1(irmhFvhrrnininini2/2)()()()1()1(nirmFFhvvnininini2/)()1()()1(3.3平衡态分子动力学模拟3.3平衡态分子动力学模拟一般对于能量确定的系统不可能给出精确的初始条件。可以将初始位置固定在格子的格点上，而初始速度则由波尔兹曼分布得出。以上述方法建立的系统一般不会具有所需要的能量，需要在模拟过程中逐渐调节系统能量达到给定值。步骤为：(1)对运动方程求解若干步；(2)计算动能和势能；(3)若能量不等于所需的值，对速度进行标度；(4)从第一步开始重复，直至系统达到平衡为止。速度标度因子:2.正则系综的分子动力学模拟速度标度因子:212161/iii*vm)N(kT2/12)43(iimvkTN3.3平衡态分子动力学模拟3.3平衡态分子动力学模拟正则系综分子动力学的模拟具体步骤:(Verlet算法的速度形式)（1）给定初始空间位置，（2）给定初始速度，（3）利用公式：计算在第n+1步时所有粒子所处的空间位置，（4）计算在第n+1步时所有粒子的速度：动能和速度标度因子：（5）计算将速度乘以标度因子的值，并让该值作为下一次计算时，第n+1步粒子的速度：（6）返回到步骤（3），开始第n+2步的模拟计算。mhFvhrrnininini2/2)()()()1(m/FFhvv)n(i)n(i)n(i)n(i2112)1(21inikvmE2/12)1()()43(inivmkTN11ninivv3.4MD在材料科学中的应用3.4.1高分子链动力学模拟计明娟等用高温淬火分子动力学模拟方法研究了甲硫氨酸一脑啡肽(Met-enkephalin)在真空中的构象性质，经聚类分析和能量优化得到了十三个低能构象，并与吗啡进行了空间拟合。3.4MD在材料科学中的应用按图l的序列构造出Met-enkephalin分子并进行初始的优化，然后以每0.1ps温度增加l0K进行动力学运算，加热到800K。在此温度下进行高温淬火分子动力学模拟计算。平衡时间为10ps，模拟时间为30ps，时间步长为1.0fs，取样间隔为0.2ps，在正则(NTV)系综下运行．对平衡后的构象作进一步分析，根据它们骨架(φ，Ψ)构象分布的分析，以及所有结构中每一个残基Co-C骨架片段之间均方根偏差不同性的计算产生五个等级，将均方根偏差最小间隔(0.1nm≤rms≤0.2nm)的结构分成十三个簇，从每个簇中选取最低能量构象作为此簇代表，最后将十三个所选结构中的每一个结构利用最陡下降法和共轭梯度法分别运行1000步进行完全的能量优化。直至rms能量偏差0.004kJ•mol而收敛，然后将优化后的十三个构象与吗啡进行空间拟合。图3.2经QMD模拟所得的十三个低能构象3.4MD在材料科学中的应用模拟结果表明，Met-enkephalin是一种构象柔性大的多肽分子，可具有多种构象，主要为Gly2-Gly3和Gly3-Phe4之间的β折叠型式。大多数的结构具有由氢键连接两个末端残基所形成的假大环，侧链和芳香环定位在假大环的同侧。酪氨酸残基中的氨基N相当于吗啡中哌啶环上的季胺N，酪氨酸残基中的酚环相当于吗啡中的酚环，苯丙氨酸残基中的芳香环相当于吗啡中环己烯环，此肽与吗啡在空间结构(药效基团)上非常相似，都满足阿片受体的要求。所以，它们可作用于同一受体。3.4MD在材料科学中的应用3.4.2脆性断裂模拟单德彬等为了研究在微观尺度上单晶铜的弯曲裂纹萌生和扩展机理，建立了单晶铜弯曲变形的分子动力学模型，应用速度标定法控制温度，采用Morse势进行了单晶