分子动力学模拟方法的基本原理与应用.docx
- 文档编号:4684431
- 上传时间:2022-12-07
- 格式:DOCX
- 页数:8
- 大小:26.24KB
分子动力学模拟方法的基本原理与应用.docx
《分子动力学模拟方法的基本原理与应用.docx》由会员分享,可在线阅读,更多相关《分子动力学模拟方法的基本原理与应用.docx(8页珍藏版)》请在冰豆网上搜索。
分子动力学模拟方法的基本原理与应用
分子动力学模拟方法的基本原理与应用
摘要:
介绍了分子动力学模拟的基本原理及常用的原子间相互作用势,如Leonard-Jones势:
论述了几种常用的有限差分算法,如VeCet算法;说明了分子动力学模拟的几种系综及感兴趣的宏观统计量的提取。
关键词:
分子动力学模拟;原子间相互作用势;有限差分算法;
分子动力学是一门结合物理,数学和化学的综合技术。
分子动力学是一套分子模拟方法,该方法主要是依靠牛顿力学来模拟分子体系的运动,以在由分子体系的不同状态构成的系统中抽取样本,从而计算体系的构型积分,并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。
从统计物理学中衍生出来的分子动力学模拟方法(MolecularDynamicsSimulation,MDS),实践证明是一种描述纳米科技研究对象的有效方法,得到越来越广泛的重视。
所谓分子动力学模拟,是指对于原子核和电子所构成的多体系统,用计算机模拟原子核的运动过程,从而计算系统的结构和性质,其中每一个原子核被视为在全部其他原子核和电子所提供的经验势场作用下按牛顿定律运动。
它被认为是本世纪以来除理论分析和实验观察之外的第三种科学研究手段,称之为“计算机实验”手段,在物理学、化学、生物学和材料科学等许多领域中得到广泛地应用。
科学工作者在长期的科学研究实践中发现,当实验研究方法不能满足研究工作的需求时,用计算机模拟却可以提供实验上尚无法获得或很难获得的重要信息;尽管计算机模拟不能完全取代实验,但可以用来指导实验,并验证某些理论假设,从而促进理论和实验的发展。
特别是在材料形成过程中许多与原子有关的微观细节,在实验中基本上是无法获得的,而在计算机模拟中即可以方便地得到。
这种优点使分子动力学模拟在材料研究中显得非常有吸引力。
分子动力学模拟就是用计算机方法来表示统计力学,作为实验的一个辅助手段。
分子模拟就是对于原子核和电子所构成的多体系统,求解运动方程(如牛顿方程、哈密顿方程或拉格朗日方程),其中每一个原子核被视为在全部其它原子核和电子作用下运动,通过分析系统中各粒子的受力情况,用经典或量子的方法求解系统中各粒子在某时刻的位置和速度,以确定粒子的运动状态,进而计算系统的结构和性质。
该模拟技术主要涉及粒子运动的动力学问题,与蒙特卡罗模拟方法(简称MC)相比,分子动力学是一种“确定性方法”,它所计算的是时间平均,而MC进行的是系综平均。
然而按照统计力学各态历经假设,时间平均等价于系综平均。
因此,两种方法严格的比较计算能给岀几乎相同的结果。
经典的分子动力学方法是Alder等于1957年提岀并首先在“硬球”液体模型下应用,发现了由Kirkwood在1939年根据统计力学预言的“刚性球组成的集合系统会发生有液相到结晶相的转变”。
后来人们称这种相变为Alder相变oRahman于1963年采用连续势模型研究了液体的分子动力学模拟。
1972年Less等发展了该方法并扩展了存在速度梯度的非平衡系统。
1980年Andersen等创造了恒压分子动力学方法。
1983年GilIan等将该方法推广到具有温度梯度的非平衡系统,从而形成了非平衡系统分子动力学方法体系。
1984年Nose等完成了恒温分子动力学方法的创建。
1985年针对势函数模型化比较困难的半导体和金属等,Car等提出了将电子论与分子动力学方法有机统一起来的第一性原理分子动力学方法。
1991年Cagin等进一步提岀了应用于处理吸附问题的巨正则系综分子动力学方法。
20世纪80年代后期,计算机技术飞速发展,加上多体势函数的提岀与发展,使分子动力学模拟技术有了进一步的发展。
1、分子动力学的运动方程:
分子动力学模拟的岀发点是假定粒子的运动可以用经典动力学来处理,对一个由N个粒子构成的孤立体系,粒子的运动由牛顿运动方程决定,也就是:
mid2ri/dt2=-VjV(rbr2,rM),式中,G分别为第i个原子的质量和位置。
▽i二-?
/?
「,V(rbr2,rN)为体系所处的势。
2、运动方程的数值积分:
计算机模拟方法的基点是利用现代计算机高速和精确的优点,对几百个以至上千个分子的运动方程进行数值积分有许多不同的积分方法,它们的效率和方便程度各异问题基本上就是用有限差分法来对二阶常微分方程进行积分常用的有以下几种方法:
、Verlet算法⑴:
Verlet算法是在60年代后期岀现的,对扩散分子的质心运动的积分是最稳定的也是最常用的数值方法。
它运用t时刻的位置和加速度以及t时刻的位置来预测t+6t位置,其积分方案,以三阶Taylor展开为基础,由以下方程给岀:
r(t+6t)=2r(t)-r(t-6t)+6t2a(t)
这里,为简单计,省略了i。
速度可按微分的基本法则得岀:
V(t)=[r(t+6t)-r(t-6t)]/26to这种算法的优点是占有计算机的内存小,并且很容易编程。
但它的缺点是位置r(t+6t)要通过小项(6t2)与非常大的两项2r(t)和r(t-6t)的差相加得到,这容易造成精度损失。
并且从式中可以看岀,这种算法不是一个自启动算法,新位置必须由t和t-6t时刻的位置得到。
、Gear的预侧-校正算法:
这种算法分为三步来完成:
首先,根据Taylor展开,预测新的位置、速度和加速度。
然后,根据新的计算的力计算加速度。
这个加速度再由与Taylor级数展开式中的加速度进行比较,两者之差在校正步里用来校正位置和速度项。
这种方法的缺点就是占有计算机的内存大。
、“蛙跳”(Leap-frog)算法⑵:
Hockey提岀的Leap-frog算法是Verlet算法的变化,这种方法设计半时间间隔的速度,即:
r(t+6t)=r(t)+6tv(t+6t/2),v(t+6t/2)=v(t-6t/2)+6ta(t)«t时刻的速度由下式给岀:
V(t)=[v(t+6t/2)+v(t-6t/2)]/2o这种算法与Verlet算法相比有两个优点:
(1)包括显速度项;
(2)收敛速度快,计算量小。
这种算法明显的缺陷是位置和速度不同步。
除了上述提及的几种方法外,还有Beeman算法、Rahman等。
3、周期性边界条件和长程力:
即使是使用现代的巨型计算机,MD方法还是只能用于粒子数大约是几百到几千的系统。
这就引起一个问题:
用这样少量的粒子,如何来模拟宏观体系为了解决这个问题,引入了周期性边界条件⑶。
采用这种方法,模拟体系实际上是由基本单元(也称为模拟计算元胞)在各个方向上重复叠合而成。
但在模拟中只需保留基本单元,所有其它单元与基本单元由平移对称性关联。
在处理粒子之间的相互作用时,通常采用“最小影像”约定。
这个约定是在由无穷重复的MD基本模拟计算元胞中,一个粒子只与它所在的基本元胞内的另外N-1个(设在此元胞内有N个粒子)中的每个粒子或其最邻近影像粒子发生相互作用。
实际上,这个约定就是通过满足不等式条件rc 4、势函数: MD模拟结果准确与否的关键在于对系统内的原子之间相互作用势函数的选取,总的来说,原子(或分子)之间的相互作用势的研究进展一直很缓慢,在一定程度上制约了MD方法在实际研究中的应用,原子间的势函数的发展经历了从对势到多体势的过程,对势认为原子之间的相互作用是两两之间的作用,与其它原子的位置无关,而实际上,在多原子体系中,一个原子的位置不同,将影响其它原子间的有效相互作用。 所以,多体势能更准确地表示多原子体系势函数。 、对势: 在分子动力学模拟的初期,经常采用的就是对势。 对势可以分为间断对势⑴和连续对势,而连续对势主要有以下几种: Lennard-Jones(L-J)势、BornTande(B-L)势、Morse势和Johnson势等,其中,L7势是为描述惰性气体分子之间的相互作用而建立的,因此它表达的力比较弱,描述材料的行为也比较柔韧,也可以用来描述过渡金属原子之间的相互作用°B-L势是用来描述离子晶体离子之间的相互作用的。 Morse势和Johnson势多用于描述金属原子之间的相互作用。 对势虽然简单,得到的结果往往也符合某些宏观物理规律,但其缺点是必然导致Cauchy关系,所以,对势实际上不能准确地描述晶体的弹性性质。 、多体势: 多体势是在20世纪80年代初期开始岀现的,1984年Daw和Baskes首次提岀了原子嵌入(EAM)势⑷。 此势的基本思想是: 把晶体的总势能分成位于晶格点阵上的原子核之间的相互作用对势和原子核镶嵌在电子云背景中的嵌入能(多体相互作用)两部分,其中,对势和多体作用势的函数形式往往根据经验选取。 基于EAM势的势函数还有很多种,这些多体势大多用于金属的微观模拟,此外,还有许多形式的多体势函数形式,1987年Jacobsen等人在等效介质原理的基础上提岀的另一种多体势函数形式,由于其简单、有效,也得到了广泛地应用。 5、分子动力学模拟的系综: 平衡态分子动力学模拟是在一定的系综瞎进行的,经常用的平衡系综是NVT或NPT系综。 在这两种系综中,牵涉到控制温度和压力的几种技术,分别如下: 、控温方法: (1)速度标定: 系统的温度与动能存在如下关系: Ek=为呼=勢輕,式中,N是原子数,汕是约束数,焉是Boltzmann常数,v,是原子i的速度。 由于系统的温度和动能存在这样的关系,所以一种最简单和最直观的方法是直接对速度进行标定。 这种方法的基本思想是: 如果t时刻的温度是T(t),速度乘以因子入后,温度的变化为△T=(X2-1)T(t),其中,X二、/T/T(t),T为所控制体系的温度。 (2)Berendsen恒温槽方法: 这种控温方法假设所模拟的体系与一个恒温槽连载一起,则两者之间就可以通过热交换而使模拟体系达到恒温的目的,方法如下,定义一个参数入: 入二J+2(T二TQ/t,式中T表征系统与恒温槽之间的热交换速率,At为MD的时间步长,那么通过v"FZw校正即可保持体系的温度在T。 附近震动,而参数T(通常取为~)则可用于控制这个震动幅度。 (3)Nose-Hoover方法: 这种方法是通过改变模拟体系的HamiItonian来实现控温的,因而有更强的物理意义,其基本思想就是在Hamiltonian加入一个假想的项来代表一个恒温源,具体做法如下: H二丽疋/2+V(r)+Qg2/2+gKTInS,式中,$和§分别是假想项的坐标和动量。 这样体系的微分方程就变为: v-dr/dr,ai=-(dV/dr+mivi©)/mi,d芽=(Emjv|-gKBT)/Qo式中,g为体系的自由度,Q为一个可调参量,表征着假想项的质量,T为温度。 在这三种控温方法中,速度标定是一种非真实的物理效应。 但这种方法可以使系统很快达到平衡,在经典分子动力学方法中,这是一种比较常用的控温方法。 Berendsen控温方法是通过系统和恒温槽进行热交换来控制温度,此方法的优点是非常简单,应用起来非常方便。 Nose-Hoover控温方法是基于统计力学而提岀来的一种控温方法,当系统和恒温槽进行热交换时,在系统中粒子出现的几率遵从统计力学规律,所以这是一种真实的物理效应。 、控压方法: 用MD方法研究压力的诱导相变和结构重构,在等压下模拟比在等体积下更容易实现。 在等压模拟下,可以通过改变模拟元胞的三个方向或一个方向的尺寸来实现体积的变化,类似于温度控制的方法,也有许多方法用于压力控制,在实际应用中,有以下两种较为常用的压力控制方法。 (1)Berendsen方法: 这种控压方法的基本思路与基本控制方程均与他的控温方法类似,如下: H=[1+At(P-Po)Y/t]1/3,式中Y是一个可调参数,Po为所控制的压力。 然后,在每一次的MD迭代中,粒子坐标x,y,z均用口相乘,得到新的坐标,即可实现对压力P的控制。 该方法的缺点是可能导致在很长的模拟时间步长内,保持着起伏。 (2)Parrinello-Rahaman方法: 这种方法是通过改变体系的Lagrangian来实现的,新的Lagrangian定义如下: 1-二》miSiGh/2-&=>jGij)+WTrh'h/2-pQ,其中,h二{a,b,c}是模拟单胞的基矢,G=h'h,原子坐标I为n=hs,V为原子之间的相互作用势,0为单胞体积,Tr表示矩阵的迹,W是一个可调参量。 这种方法的有点是模拟单胞不仅大小可变,形状也可变。 6、热力学性质的计算: 经典力学的一个基本前提是,只要知道了物理体系在相空间的运动,就可以导岀其所有的宏观性质。 在分子动力学中,对模拟得到的位形进行结构分析是很重要的。 比如在研究晶体的熔化时,需要及时分析它的结构变化。 就结构分析而言,目前常用的方法有: 径向分布函数、静态结构因子和配位数等方法。 、径向分布函数: 径向分布函数可以由下式来定义: g(r)=^n(r)/4nr2Ar,式中,n(r)表示距 离原点r到r+Ar之间的平均粒子数,V是模拟的体积。 RDF表示的物理含义是: 再空间位置r点周围的体积元中单位体积内发现另一个粒子的几率。 由此看出,RDF表征着结构的无序化程度。 、静态结构因子: 静态结构因子表示式为: S(k)=|^=1exp(iKrj)|/N,其中,N是总原子数,K是倒格矢,rj为原子位置矢量。 对理想晶体而言,SSF为1,而对理想流体,则为0。 、配位数: 所谓配位数,就是原子的第一最近邻原子的个数配位数的多少,表示了该原子周围的原子分布的密度大小,它经常在一些分析中作为结构分析的一种辅助手段。 7、MD方法的进一步研究方向: 经典分子动力学方法存在两个缺陷: (1)元胞体积和形状保持不变,限制MD方法的应用; (2)不适合含有自由电子的系统,对金属等系统的计算的结果不理想。 为克服经典分子动力学方法的局限性,20世纪80年代初AndersorivParrinelIo等人发展了可变元胞分子动力学方法。 由于能带论中的密度泛函(DF)方法是研究凝聚态物质电子结构的有效方法,1985年Car等人提岀了第一原理的MD方法。 但是,进一步拓宽分子动力学方法的应用面显然是一个热门的发展方向。 对MD方法本身而言,最重要的两个要素是初始位形的给定和相互作用势函数的选取。 MD模拟时间足够长,初始条件的选择不会影响计算的结果,但是选择合理的初始条件可以加快系统趋于平衡,可以节省机时。 影响MD计算结果精确程度的最主要的因素是相互作用势的精确性,人们对提高势函数的精确性进行了大量的研究提岀了许多势函数的形式,但是,对大多数原子而言,努力找到即准确而形式又不太复杂的势函数,仍然是一个很有挑战性的课题。 我们认为,把第一原理、量子化学分析和参数拟合相结合是势函数研究的最好的方法,NanxianChen等人就用此方法研究了氯化钠晶体离子之间的相互作用势,并研究了微结构的变化,这是一个值得进一步研究的方向。 ⑸ 8、分子动力学模拟方法的应用: 、氨的自扩散系数的分子动力学模拟研究: 扩散系数是描述传递现象的重要的流体热物理性质之一,较宽温度和压力范围的扩散系数在化工生产中具有重要应用价值。 扩散系数的实验测量是一项很困难的工作,而分子动力学(MD)模拟撒获得分子扩散系数的有效方法。 氨,作为大气中最重要的碱性气体之一,它在大气化学、物质传输与沉降过程中发挥着重要作用。 采用MD模拟方法计算氨在较宽温度和压力范围的自扩散系数,通过将两种不同的控温方式得到的模拟值与大量的实验值进行比较,验证了所采用的模拟方法是合理的。 从而可有效地预测极端条件下实验难以测量的扩散系数。 MD模拟采用TINKER程序包,OPLS-AA全原子力场。 模拟体系包含300个分子。 分子间的作用势采用Lennard-Jones势,体系中的长程作用力采用Ewald加和的形式。 在模拟过程中均采用周期性边界条件、Beeman算法求解运动方程。 温控方式分别采用Andersen和Berendsen。 体系的截断半径为1nm,时间步长为1fs。 MD模拟在NVT系综进行。 通过模拟发现: 从常温到高温,采用Berendsen控温方式得到的氨的自扩散系数的模拟值与实验值吻合的很好。 于是,可以采用分子动力学模拟来代替实验,获得高温高压条件下实验难以测量的分子自扩散系数。 另一方面,氨的自扩散系数受温度的影响比受压力的影响更明显。 ⑹ 、半晶态聚合物的分子动力学模拟: 聚合物材料区别于其它材料的一个典型特征是它具有半结晶形态.当从熔点以上温度开始缓慢冷却到室温时,很多聚合物将发生部分区域结晶、部分保留为非晶态的转变,而形成半晶态的结构。 实验观测到折叠链结构为聚合物结晶的基本形式,由折叠链片晶可形成多层堆叠片晶、球晶、串晶和枝晶等。 针对半晶态聚合物的总体形态,人们提岀了多种多样的聚合物结构模型,如缨状微束模型、折叠链模型和折叠链缨状胶粒模型等。 采用Meyer和Muler-PIathe提岀的粗粒化聚乙烯醇模型,应用分子动力学方法,模拟熔融态聚合物经过缓慢冷却、局部结晶形成半晶态聚合物的过程,研究冷却过程中聚合物的结晶和凝固行为,揭示半晶态聚合物的结构特征及其形成机制。 为半晶态聚合物物理、力学性能的分子模拟研究提供基础。 应用分子动力学方法,模拟熔融态聚合物经过缓慢冷却,局部结晶形成半晶态聚合物的过程。 静态结构因子的演变显示岀,在结晶行为的初期小角散射强度的增大先于布拉格峰的岀现,这个现象与小角/大角x射线散射实验结果相一致。 模拟得到了半晶态的聚合物结构,其中的微小晶区由折叠链构成,它们随机取向分布在非晶态区中,总体形态可用缨状微束结构模型描述。 模拟得到的半晶态聚合物的结晶度为45%,分子链持久长度(即晶区平均厚度)约为。 通过比较结晶度和分子链持久长度随温度的变化曲线,发现在不同的冷却阶段具有不同的有序结构形成机制。 从随机缠绕的分子链到折叠链晶区,需要两种形式的结构转变: 分子链的伸展和伸直分子链之间的平行排列。 在从结晶温度到玻璃化温度的凝固转变过程中,存在上述两种形式的结构转变;而在玻璃化温度之后,材料的活性只能允许调整伸直分子链之间的相对排列位置,但由此也导致了聚合物结晶度明显的增大。 ⑺ 、多肽抑制剂抑制淀粉质多肽42构象转换的分子动力学模拟和结合自由能计算: 应用分子动力学模拟和结合自由能计算方法研究多肽抑制剂KLVFF.VVIA和LPFFD抑制淀粉质多肽42(A(342)构象转换的分子机理。 结果表明,三种多肽抑制剂均能够有效抑制AB42的二级结构由a-螺旋向B-折叠的构象转换。 另外,多肽抑制剂降低了A042分子内的疏水相互作用,减少了多肽分子内远距离的接触,有效抑制了AB42的疏水塌缩,从而起到稳定其初始构象的作用。 这些抑制剂与Ap42之间的疏水和静电相互作用(包括氢键)均有利于它们抑制A042的构象转换。 此外,抑制剂中的带电氨基酸残基可以增强其和AP42之间的静电相互作用(包括氢键)并降低抑制剂之间的聚集,从而大大增强对AM2构象转换的抑制能力。 但脯氨酸的引入会破坏多肽的线性结构,从而大大降低其与AM2之间的作用力。 上述分子模拟的结果揭示了多肽抑制剂紅VFFvVVIA和LPFFD抑制Ap42构象转换的分子机理,对进一步合理设计AB高效短肽抑制剂具有非常重要的理论指导意义。 ⑻ 、超临界水中醋酸锌水解反应的分子动力学模拟: 超临界水具有特殊的溶解度、可变的密度、较低的粘度和表面张力以及较高的扩散特性,利用超临界水热合成法制备纳米ZnO逐渐成为常用的方法之一。 在这种方法中,超临界水(SCW,温度T二647K,压力p二22MPa,密度p=g.cm_3)同时作为反应介质和反应物,含Zn离子的盐如Zn(N03)2vZnS04vZn(CH3COO)2的水溶液在超临界状态下与碱(如KOH)或纯水先发生水解反应生成Zn(OH)2然后Zn(0H)2脱水得到纳米级的ZnO颗粒。 利用超临界水热合成还成功制备岀不同种类和形貌的纳米颗粒,如MgO、TiO、Zr02vSnO2XZnS和FesO.等。 利用分子动力学模拟Zn(CH3COO)2与SCW反应的能量变化,以及产物CH3COOH在SCW的分布特性,为实验上通过控制温度或压力生成纳米ZnO颗粒提供参考。 利用分子动力学模拟研究Zn(CH3COO)2在SCW中的平衡特征以及水解反应过程的能量变化,定性分析了水分子在Zn(CH3C00)2团簇内部由于静电场的影响而发生解离的可能性,得岀以下结论: (1)Zn(CH3COO)2在SCW中易聚集成团,1个Zn周围平均结合5个CHXXXy和1个H20满足Zn形成6配位八面体的空间构型的需要。 Zn"在Zn(CH3C00)2团簇内部和表明的配位情况有所不同,更多的H20在Zn(CH3C00)2团簇表面参与Zn"配位。 (2)在标准状? S,Zn(CH3COO)冰解生成Zn(0H)2和CH3COOH的反应焙变为・moL,从热力学的角度看难以发生,但在超临界条件下,模拟发现上述水解反应使得体系势能降低,表明反应容易发生,并且水解反应发生的同时伴随Zn(CH3COO)2团簇结构的改变。 水解产物0H「容易进入Zn(CH3COO)2团簇内部,富集Zn,而产物CH3COOH扩散进入水相,并倾向分布在水相一气相界面。 (3)定性分析处于Zn(CH3C00)2团簇内部的水分子周围静电势的分布发现,水分子周围分布有很强的负的静电势,容易使水分子极化变形进而解离,促进Zn(CH3COO)2在SCW中的水解反应。 ⑼ 参考文献: [1]fIuids.I.ThermodynamicaIpropertiesofLennard-JonesmolecuIes[J]. Phys.Rev.r1967,159: 98-103 [2]ThepotentialcaIcuIationandsomeappIications[J].Methodsin ComputationaIPhysics,1970,9: 136-211 [3].ofstatecaIculationsbyfastcomputingmachines[J]・atommethodderivationandappIicationtoimpurities,surfaces,andotherdefectsinmetaIs[J].崔守鑫,胡海泉,肖效光,黄海军,分子动力学模拟基本原理和主要技术[0].聊城大学学报,2005,18 (1): 30-34 [6]周昌林,黎多来,李东凯,等,氨的自扩散系数的分子动力学模拟究[J].广东化工,2012,39(13): 183-184 [刀段芳莉,颜世铛,半晶态聚合物的分子动力学模拟[A].计算物理,2012,29(5): 759-765 [8]董晓燕,都文婕,刘夫锋,多肽抑制剂抑制淀粉质多肽42构象转换的分子动力学模拟和结合自由能计算[0],物理化学学报,2012,28(11),2734-2744 [9]王晓娟,李志义,刘志军,超临界水中醋酸锌水解反应的分子动力学模拟[0],物理化学学报,2013,29 (1)
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 分子 动力学 模拟 方法 基本原理 应用