分子动力学方法

1、分子动力学方法分子动力学方法3.1 3.1 基本原理计算机模拟分类：(1)随机模拟方法。 优点: 随机模拟方法计算的程序简单，占内存少，但是该方法难于处理非平衡态的问题。(2)分子动力学方法（Molecular Dynamics或简称MD） 。 可以处理非平衡态问题。但是使用该方法的程序较复杂，计算量大，占内存也多。 分子动力学方法利用牛顿古典力学来计算许多分子在相空间中的轨迹。 3.1 3.1 基本原理l分子动力学(MD)方法的出发点是物理系统的确定的微观描述（哈密顿描述方程、拉格朗日方程或牛顿运动方程）。因此，分子动力学方法是用运动方程来计算多体或少体系统的性质，结果得到的既有系统的静态特

2、性，也有动态特性。蒙特卡罗方法只能得到系统的位形特性。lMD方法的具体做法是在计算机上求运动力程的数值解。为此，通过适当的格式对方程进行近似，使之适于在计算机上求数值解。其实质是计算一组分子的相空间轨道，其中每个分子各自都服从经典运动定律。这里的系统不仅是点粒子系统，也包括具有内部结构的粒子组成的系统。3.1 3.1 基本原理1.分子动力学是在原子、分子水平上求解多体问题的重要的计算机模拟方法。2.通过求解所有粒子的运动方程，分子动力学方法可以用于模拟与粒子运动路径相关的基本过程。3.在分子动力学中，粒子的运动行为是通过经典的Newton运动方程所描述。3.1.13.1.1粒子运动方程的数值求

3、解粒子体系的运动方程粒子体系的运动方程Lagrangian方程方程1.Lagrangian函数的定义为函数的定义为 则运动的则运动的Lagrangian方程为方程为0ddkkqqtLL3.1.13.1.1粒子运动方程的数值求解 粒子的运动方程粒子的运动方程-Lagrangian-Lagrangian方程方程kkkkkkkkkqHppHqLpqHqLp),(),(qqqp3.1.13.1.1粒子运动方程的数值求解l单原子的牛顿运动方程单原子的牛顿运动方程V.iiiiiiirmmrrfp3.1.13.1.1粒子运动方程的数值求解lVerlet 算法算法 r(t+t) = r(t) + v(t)t

4、+ (1/2)a(t)t 2 (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.13.1.1粒子运动方程的数值求解 其他求解算法：其他求解算法：lLeap-frog

5、 algorithm r(t+t) = r(t) + v(t+(1/2)t) t v(t+(1/2)t) = v(t-(1/2)t) + a(t) tlBeemans algorithm r(t+t) = r(t) + v(t)t + (2/3)a(t)t2 (1/6)a(t-t)t2 v(t+t) = v(t) + v(t)t + (1/3)a(t)t + (5/6)a(t)t (1/6)a(tt)t3.1.2 3.1.2 热力学量的计算热力学量的计算l在物理系统的计算机模拟中，系综平均必须用时间平均代替，在通常的模拟中，粒子数N和体积V是固定的。给定初始位置rN(0)和 初始动量pN(0)

6、后，一个MD算法将从运动方程生成轨道(rN(t), pN(t)，轨道平均的定义为l假定能量守恒，并且轨道在一切具有同一能量的相同休积内经历相同的时间，则轨道平均等于微正则系综平均100)();()()(limttNNttVtptrdtAttANVEAA 3.1.2 3.1.2 热力学量的计算热力学量的计算l孤立系统的总能量是一个守恒量，沿着分子动力学模拟生成的任何一条轨道，能量应保持不变，即 。l孤立系统的动能和势能不是守恒量，它们的大小沿着生成的轨道逐点变化 EE 101000,)()(lim,)()(limtttttktkdttrUttUdttEttEnnvvkkEnnE0013.1.2

7、3.1.2 热力学量的计算热力学量的计算l生成的动能路径是不连续的，必须在各个时间间断点上计算动能的值以求平均:l其中l从平均动能可以计算系统的温度，温度是一个重要的物理量，需要加以监测，特别是在模拟的起始阶段。nnvvkkEnnE001ivivkmE)(2123.1.2 3.1.2 能量均分定理能量均分定理l在热力学极限下，一切系综都是等同的，并且可以应用能量均分定理。l热力学极限是指粒子数（或体积）趋向无穷大时的极限。一般宏观物体包含了1023个粒子，可以认为是满足热力学极限的。相变只有在热力学极限下才会发生。l系统的哈密顿量为l则： ijiijiirumP221TkmBi212123.1

8、.2 3.1.2 能量均分定理能量均分定理l由于系统的每个粒子有三个自由度，因此 l假定位势在处被截断。系统内部的位形能量的轨道平均值为l其中TNkEBk23nnvvUnnU001jivijvruU)(3.1.2 3.1.2 能量均分定理能量均分定理l由于位势被截断，总能量和势能含有误差，为了估计必须作出改正。所有的内部位形能都加到截止距离为止，尾部改正可以取l其中是g(r)对关联函数，它是在原点r=0处有一个粒子时,在r周围的体积元dr内找到一个粒子的概率，令n(r)为离一个给定粒子的距离在r和r+dr之间的平均粒子数，于是crcdrrrgruU2)()(2rrrnNVrg24)()(3.1

9、.2 3.1.2 能量均分定理能量均分定理l其他的量也需要尾部改正，以压强的计算为例，这时维里(virial)状态方程成立。l至于势能的计算，可以把积分分成两项，一项是由相互作用力程之内的贡献引起的，一项是对位势截断的改正项：l长程改正项为 0324)(6drrrurgTPkBCjiijijkBPrurNTP6crCdrrrurgP324)(63.2 分子动力学模拟的基本步骤l分子动力学模拟的实际步骤可以划分为四步：l设定模拟所采用的模型；l给定初始条件；l趋于平衡的计算过程；l宏观物理量的计算。3.2 分子动力学模拟的基本步骤l1模拟模型的设定模拟模型的设定l硬球势硬球势lLennard-J

10、ones型势型势 .r,r如果如果 ,0,rU 6124rrrU0.60.81.01.21.41.61.82.02.22.42.6-3-2-101234 位势V(r)力F(r)吸引力排斥力v根据经典物理学的规律，可以知道在系综模拟中的守恒量。v微正则系综的模拟中能量、动量和角动量均为守恒量。在此系综中他们分别表示为： 元胞 周期性边界条件, iiirVrmE221iipPiiiprM3.2 分子动力学模拟的基本步骤分子动力学模拟的最小像力约定示意图 Lnrrminrjiij最小像力约定2.给定初始条件 给定粒子的初始位置和速度的数值： （1）令初始位置在差分网格格子上，初始速度从玻尔兹曼分布随

11、机抽样得到。（2）令初始位置随机地偏离差分网格格子，初始速度为零。（3）令初始位置随机地偏离差分网格格子，初始速度从玻尔兹曼分布随机抽样得到。 3趋于平衡 使系统达到平衡，模拟中需要一个趋衡过程。在这个过程中，我们增加或从系统中移出能量，直到系统具有所要求的能量。 3.2 分子动力学模拟的基本步骤3.3 平衡态分子动力学模拟1. 微正则系综的分子动力学模拟 粒子数恒定、体积恒定、能量恒定、整个系统的总动量恒等于零。 分子动力学模拟步骤如下： （1）给定初始空间位置。 （2）计算在第n步时粒子所受的力。 （3）计算在第n +1步时所有粒子所处的空间位置。 （4）计算第n步的速度。 （5）返回到步

12、骤（2），开始下一步的模拟计算。mhFrrrnininini/22)()1()()1(hrrvninini2/)1()1()(Verlet算法的速度形式： （1）给定初始空间位置 。 （2）给定初始速度 。 （3）利用公式： 计算在第n+1步时所有粒子所处的空间位置 。 （4） 计算在第n+1步时所有粒子的速度： （5） 返回到步骤（3），开始第n+2步的模拟计算。)1(iv)1(irmhFvhrrnininini2/2)()()()1()1( nirmFFhvvnininini2/)()1()()1(3.3 平衡态分子动力学模拟3.3 平衡态分子动力学模拟l一般对于能量确定的系统不可能给出精

13、确的初始条件。可以将初始位置固定在格子的格点上，而初始速度则由波尔兹曼分布得出。l以上述方法建立的系统一般不会具有所需要的能量，需要在模拟过程中逐渐调节系统能量达到给定值。步骤为：l(1)对运动方程求解若干步；l(2)计算动能和势能；l(3)若能量不等于所需的值，对速度进行标度；l(4)从第一步开始重复，直至系统达到平衡为止。速度标度因子: 2.正则系综的分子动力学模拟速度标度因子: 212161/iii*vm)N(kT2/12)43(iimvkTN3.3 3.3 平衡态分子动力学模拟3.3 3.3 平衡态分子动力学模拟 正则系综分子动力学的模拟具体步骤: (Verlet算法的速度形式) （1

14、）给定初始空间位置， （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(i2112)1(21inikvmE2/12)1()()43(inivmkTN11ninivv3.4 MD在材料科学中的应用在材料科学中的应用l3.4.1 高分子链动力学模拟l计明娟等用高温淬

15、火分子动力学模拟方法研究了甲硫氨酸一脑啡肽(Met-enkephalin)在真空中的构象性质，经聚类分析和能量优化得到了十三个低能构象，并与吗啡进行了空间拟合。3.4 MD在材料科学中的应用在材料科学中的应用l按图l的序列构造出Met-enkephalin分子并进行初始的优化，然后以每0.1 ps温度增加l0K进行动力学运算，加热到800K。在此温度下进行高温淬火分子动力学模拟计算。平衡时间为10ps，模拟时间为30 ps，时间步长为1.0fs，取样间隔为0.2ps，在正则(NTV)系综下运行对平衡后的构象作进一步分析，根据它们骨架( ， )构象分布的分析，以及所有结构中每一个残基Co-C骨架

16、片段之间均方根偏差不同性的计算产生五个等级，将均方根偏差最小间隔(0.1 nmrms0.2 nm)的结构分成十三个簇，从每个簇中选取最低能量构象作为此簇代表，最后将十三个所选结构中的每一个结构利用最陡下降法和共轭梯度法分别运行1000步进行完全的能量优化。直至rms能量偏差0.004 kJmol而收敛，然后将优化后的十三个构象与吗啡进行空间拟合。图3.2 经QMD模拟所得的十三个低能构象3.4 MD在材料科学中的应用在材料科学中的应用l模拟结果表明，Met-enkephalin是一种构象柔性大的多肽分子，可具有多种构象，主要为Gly2- Gly3和Gly3-Phe4之间的折叠型式。大多数的结构

17、具有由氢键连接两个末端残基所形成的假大环，侧链和芳香环定位在假大环的同侧。酪氨酸残基中的氨基N相当于吗啡中哌啶环上的季胺N，酪氨酸残基中的酚环相当于吗啡中的酚环，苯丙氨酸残基中的芳香环相当于吗啡中环己烯环，此肽与吗啡在空间结构(药效基团)上非常相似，都满足阿片受体的要求。所以，它们可作用于同一受体。3.4 MD在材料科学中的应用在材料科学中的应用l3.4.2 脆性断裂模拟l单德彬等为了研究在微观尺度上单晶铜的弯曲裂纹萌生和扩展机理，建立了单晶铜弯曲变形的分子动力学模型，应用速度标定法控制温度，采用Morse势进行了单晶铜的弯曲变形分子动力学模拟。3.4 MD在材料科学中的应用在材料科学中的应用

18、l二维分子动力学模型如图3.3。图3.3 弯曲变形仿真实验模型(无凹槽模型)图3.4弯曲变形仿真实验模型(无凹槽模型)3.4 MD在材料科学中的应用在材料科学中的应用l以上模拟晶面是(100)晶面，原子位置按理想点阵排列。原子运动按牛顿原子处理，牛顿原子的作用由Morse势函数计算。仿真实施后分别在模型的左右两端3列原子上施加向下的位移，同时施加弯矩(M)。计算采用速度-Verlet算法 。原子的初始位置设在fcc晶体的晶格上，原子的初始速度可利用Maxwell分布确定，也可直接置为零。温度控制采用速度标定法控制在绝对零度，以避免原子热激活的复杂影响。图3.5单晶Cu(100)晶面无凹槽模型弯

19、曲变形仿真实验图3.6单晶Cu(100)晶面有凹槽模型弯曲变形仿真实验3.4 MD在材料科学中的应用在材料科学中的应用l通过模拟试验结果图3.5发现，在位移和弯矩的作用下，个别原子间间隙的不断增大和应变能的不断积累最终使晶体内部出现空位，在128 ps已经很明显的看到空位，这些空位位于模型中部，空位的不断合并萌生裂纹，微观裂纹的顶端又形成空位，随着时间步的推移，塑性变形加剧，裂纹逐渐萌生和扩展，发展成为贯穿模型的裂纹，形成裂纹的扩展，微观裂纹在裂纹顶端形核而造成主裂纹的连续扩展，裂尖向下萌生，后续微观裂纹的扩展类似于宏观裂纹扩展，与无凹槽模型相比，图3.6有凹槽模型其裂纹萌生和扩展的时间明显缩

20、短，表明裂纹缺陷对断裂过程起促进作用。3.4 MD在材料科学中的应用在材料科学中的应用l研究结果表明：应变能的不断积累使晶体内部产生空位，材料的裂纹萌生于空位，空位的合并形成纳米级裂纹，后续微观裂纹的扩展类似于宏观裂纹；裂纹缺陷促进了裂纹的萌生和扩展。而且对无凹槽模型和有凹槽模型裂纹萌生和扩展的时间步长的对比研究表明，裂纹缺陷对断裂过程起促进作用。3.4 MD在材料科学中的应用在材料科学中的应用l铁中裂纹扩展的结构演化l吴映飞用分子动力学模拟的方法研究了bccFe中I型裂纹在应力及温度场下裂尖区的结构演化问题。图3.7裂纹的几何结构示意图 3.4 MD在材料科学中的应用在材料科学中的应用l裂纹

21、前沿沿101方向，裂纹面垂直于0l0方向，裂纹扩展沿101。加载通过增加应力强度因子(K)实现，K1由临界应力强度因子KIC度量。l基于各向同性线弹性力学的平面应变条件得到初始裂纹，进行了两组模拟。(1)第一组模拟选取初始温度T=5 K、应力强度因子K =1.0 KIC，设在所对应的外载下存在初始裂纹，随后温度增加到l00，300，和500 K，外载增加到K =2.8 KIC。在这组模拟中，显示出堆垛层错，裂纹形状尖锐，且扩展速度较快，(2)第二组模拟选取初始温度T=100 K，应力场条件同于第一组，温度加到300和500 K，外载加到K =2.8 KIC 。在这组模拟中，显示出了堆垛层错、位

22、错发射、裂尖钝化和分枝，以及孪晶带；该组裂纹的扩展速度慢于第一组。 图 3.9 KI=1.8KIC后，裂尖区出现了孪晶带(T=100K,KI=2.5KIC)图3.8 低温下裂尖区原子的ABAB型堆垛层错(T=5K,KI=1.8KIC)3.4 MD在材料科学中的应用在材料科学中的应用l研究结果表明：在裂纹扩展过程中，裂尖区结构演化与初始裂纹的工作条件相关，并且依赖于温度与外载的协同作用。在相同条件下裂纹周围局域原子的能量也受初始裂纹的工作条件影响。 3.4 MD在材料科学中的应用在材料科学中的应用l3.4.3表面效应对于晶界附近原子迁移的影响lK.Saitoh等采用分子动力学模拟方法（对微电子器件中铝芯丝材内晶界附近原子的迁移进行了分析。该模拟是基于EMT理论建立的模型，这一理论结合大量电子浓度的体效应可以验证出更加准确的晶界或表面处原子结构的研究结果。l主要工作分为两部分，即晶界附近扩散性能的模拟和晶界附近原子的重新排列。它的算法是对MD时间段（2.0fs）和运动方程的积分，在该模型中，将=5（100）晶界和自由表面分别设定为xz面和yz面。图3.10 扩散性能研究中的原子模型 3.4 MD