分子动力学模拟方法的基本原理与应用Word下载.docx
《分子动力学模拟方法的基本原理与应用Word下载.docx》由会员分享,可在线阅读,更多相关《分子动力学模拟方法的基本原理与应用Word下载.docx(8页珍藏版)》请在冰豆网上搜索。
Rahman于1963年采用连续势模型研究了液体的分子动力学模拟。
1972年Less等发展了该方法并扩展了存在速度梯度的非平衡系统。
1980年Andersen等创造了恒压分子动力学方法。
1983年Gillan等将该方法推广到具有温度梯度的非平衡系统,从而形成了非平衡系统分子动力学方法体系。
1984年Nose等完成了恒温分子动力学方法的创建。
1985年针对势函数模型化比较困难的半导体和金屈等,Car等提岀了将电子论与分子动力学方法有机统一起來的第一性原理分子动力学方法。
1991年Cagin等进一步提出了应用于处理吸附问题的巨正则系综分子动力学方法。
20世纪80年代后期,计算机技术飞速发展,加上多体势函数的提出与发展,使分子动力学模拟技术有了进一步的发展。
1、分子动力学的运动方程:
分子动力学模拟的出发点是假定粒子的运动可以用经典动力学来处理,对一个由N个粒子构成的孤立体系,粒子的运动由牛顿运动方程决定,也就是:
mid2ri/dt2=-r2,rf),式中,nu,“分别为第i个原子的质量和位置。
▽产-0/3“,V(mr2,rf)为体系所处的势。
2、运动方程的数值积分:
计算机模拟方法的基点是利用现代计算机高速和精确的优点,对几百个以至上千个分子的运动方程进行数值积分有许务不同的积分方法,它们的效率和方便程度各界问题基本上就是用河限差分法来对二阶常微分方程进行积分常用的有以卜•几种方法:
21、Verlot算法叫
Verlet算法是在60年代后期岀现的,对扩散分子的质心运动的积分是最稳定的也是最常用的数值方法。
它运用t时刻的位置和加速度以及t时刻的位置來预测t+8t位置,其积分方案,以三阶Taylor•展开为基础,由以下方程给岀:
r(t+8t)=2r(t)-r(t-8t)+8t2a(t)
这里,为简单计,省略Ti•速度可按微分的基本法则得出:
V(t)=[r(t+5t)-r(t-8t)]/25t0这种算法的优点是占有计算机的内存小,并且很容易编程。
但它的缺点是位Br(t+8t)要通过小项©
2)与非常人的两项2r(t)和心&
)的差相加得到,这容易造成精度损失。
并且从式中可以看出,这种算法不是一个自启动算法,新位置必须由t和tdt时刻的位置得到。
2.2、的预侧-校正算法:
这种算法分为三步来完成:
首先,根据Taylor展开,预测新的位置、速度和加速度。
然后,根据新的计算的力计算加速度。
这个加速度再由与Taylor级数展开式中的加速度进行比较,两者之差在校正步里用來校正位置和速度项。
这种方法的缺点就是占有计篦机的内存大。
2.3、“蛙跳”(Leap-frog)算法⑵:
Hockey提出的Leap-frog算法是Verlet算法的变化,这种方法设计半时间间隔的速度,即:
r(t+8t)=r(t)+8tv(t+5t/2),v(t+8t/2)=v(t-8t/2)+5ta(t)«
t时刻的速度由下式给Hl:
v(t)=[v(t+8t/2)+v(t-5t/2)]/2。
这种算法与Verier算法相比冇两个优点:
(1)包括显速度项;
(2)收敛速度快,计算最小。
这种算法明显的缺陷是位置和速度不同步。
除了上述提及的儿种方法外,还右°
Beeman隽法、Rahman等。
3、周期性边界条件和长程力:
即使是便用现代的巨型计SZ机,MD方法还星只能用于粒子数大约是几百到几千的系统。
这就引起一个问题:
用这样少量的粒子,如何来模拟宏观体系?
为了解决这个问题,引入了周期性边界条件卩】。
采用这种方法,模拟体系实际上是由基本单元(也称为模拟计算元胞)在各个方向上重复叠合而成。
但在模拟中只需保留基本单元,所有其它单元与基本单元由平移对称性关联。
在处理粒子之间的相互作用时,通常采用“最小影像”约定。
这个约定是在由无穷重复的MD基本模拟计算元胞中,一个粒子只与它所在的基本元胞内的另外N-1个(设在此元胞内有N个粒子)中的每个粒子或其垠邻近影像粒子发生相互作用。
实际上,这个约定就是通过满足不等式条件rc<
L/2來截断位势(门为截断半径,L是元胞的边长)通常L的数值应当选得很人,以避免有限尺寸效应,但这样会增人计算量,同城采用对相互作用势的修正来近似处理。
4、势函数:
MD模拟结果准确与否的关键在于对系统内的原子之间相互作用势函数的选取,总的来说,原子(或分子)之间的相互作用势的研究进展一直很缓慢,在一定程度上制约了MD方法在实际研究中的应用,原子间的势函数的发展经历了从对坍到多体势的过程,对势认为原子之间的相互作用是两两之间的作用,与其它原子的位置无关,而实际上,在多原子体系中,一个原子的位置不同,将影响其它原子间的有效相互作用。
所以,多体势能更准确地表示多原子体系势函数。
4.1、对势:
在分子动力学模拟的初期,经常采用的就是对势。
对势可以分为间断对势[1】和连续对势,而连续对势主要有以卜几种:
LennardJones(L-J)势、Boni-lande(B-L)势、Morse势和Johnson孙等,其中,L-J势是为描述惰性气体分子之间的相互作用而建立的,因此它表达的力比较弱,描述材料的行为也比较柔韧,也町以用来描述过渡金属原子之间的相互作用。
B丄势是用来描込离子晶体离子之间的相互作用的。
MorsP势和Johnson坍篦川丁•描述金属原子之间的相互作用°
对孙虽然简单•・得到的结果往往也符合某些宏观物理规律,但其缺点是必然导致Cauchy关系,所以,对势实际上不能准确地描述晶体的弹性性质。
4.2、多体势:
多体势是在20世纪80年代初期开始出现的,1984年Daw和Baskes首次提出了原子嵌入(EAM)势囹。
此势的基本思想是:
把晶体的总势能分成位于晶格点阵上的原子核之间的相互作用对坍和原子核镶嵌在电子云背景中的嵌入能(多体相互作用)两部分,其中,对势和多体作用势的函数形式往往根据经验选取。
基于EAM势的势函数还有很多种,这些多体势人多用丁•金属的微观模拟,此外,还有许多形式的多体势函数形式,1987年Jacobsen等人在等效介质原理的基础上提出的另一种多体势函数形式,由于其简单、有效,也得到了广泛地应用。
5、分子动力学模拟的系综:
平衡态分子动力学模拟是在一定的系综瞎进行的,经常用的平衡系综是NVT或NPT系综。
在这两种系综中,牵涉到控制温度和压力的几种技术,分别如2
5.1、控温方法:
(1)速度标定:
系统的温度与动能存在如下关系:
比=丫?
零=竺泸口,式中,N是原子数,2是约束数,Kb是Boltzmann常数,w是原子i的速度。
由于系统的温皮和动能存在这样的关系,所以一种最简单和最直观的方法是直接对速度进行标定。
这种方法的基本思想是:
如果t时刻的温度是T(t),速度乘以因子入后,温度的变化为AT二(入F)T(t),其中,入=07丽,T为所控制体系的温度。
(2)Berendsen恒温槽方法:
这种控温方法假设所模拟的体系与一个恒温槽连载一起,则两者之间就可以通过热交换而便模拟体系达到恒温的目的,方法如下,定义一个参数入:
入=+加(T一乙)/「式中T表征系统与恒温槽之间的热交换速率,At为MD的时间步长,那么通过Vnew=AVold校正即可保持体系的温度在T。
附近震动,而参数丫(通常取为0.l~0.4ps)则可用于控制这个震动幅度。
(3)Nose-Hoover方法:
这种方法是通过改变模拟体系的Hamiltonian来实现控温的,因而有更强的物理意义,其基本思想就是在Haniiltoiiian加入一个假想的项来代表一个恒温源,具体做法如下:
H=》0评/2+V(r)+Q^/2+gKTlnS,式中,S和«
分别是假想项的坐标和动量。
这样体系的微分方程就变为:
vi=dn/dr,ain(dV/dr+mivi0/mi,d#=山2—9KbT)/Q。
式中,g为体系的自由度,Q为
一个可调参量,表征着假想项的质量,T为温度。
在这三种控温方法中,速度标定是一种年真实的物理效应。
但这种方法可以使系统很快达到平衡,在经典分子动力学方法中,这是一种比较常用的控温方法。
Berendsen控温方法是通过系统和恒温槽进行热交换来控制温度,此方法的优点是非常简单,应用起来非常方便。
Nose-Hoover控温方法是基于统计力学而提出来的一种控温方法,当系统和恒温槽进行热交换时,在系统中粒子出现的儿率遵从统计力学规律,所以这是一种真实的物理效应。
5.2、控压方法:
用MD方法研究压力的诱导和变和结构重构,在等压下模拟比在等体积下更容易实现。
在等压模拟下,可以通过改变模拟元胞的三个方向或一个方向的尺寸來实现体积的变化,类似于温度控制的方法,也有许多方法用于压力控制,在实际应用中,右以下两种较为常用的压力控制方法。
(1)Berendsen方法:
这种控压方法的基本思路与皐本控制方程均与他的控温方法类似,如下:
卩=[1+M(P-Po)Y/t]1/3,式中丫是一个可调参数,Po为所控制的压力。
然后,在每一次的MD迭代中,粒子坐标x,y,z均用甘H乘,得到新的坐标,即可实现对压力P的控制。
该方法的缺点是可能导致在很长的模拟时间步长内,保持着起伏。
(2)Parrinello-Rahaman方法:
这种方法是通过改变体系的Lagrangian来实现的,
新的Lagrangian定义如下:
L=S—S;
>
1"
(中)+WTrh'
h/2—pO.
其中,/2={a,b,c}是模拟单胞的基矢,G=h'
h,原子坐标耳为r,=hs,V为原子之间的相互作用势,Q为单胞体积,Tr表示矩阵的迹,W是一个可调参量。
这种方法的有点是模拟单胞不仅大小可变,形状也可变。
6、热力学性质的计算:
经典力学的一个基本前提是,只要知道了物理体系在相空间的运动,就可以导出其所行的宏观性质。
在分子动力学中,对模拟得到的位形进行结构分析是很巫要的。
比如在研究晶体的熔化时,需要及时分析它的结构变化。
就结构分析而言,目前常用的方法育:
径向分布函数、静态结构因子和配位数等方法。
6.1、径向分布函数:
径向分布函数可以由下式来定义:
g(r)=^n(r)/4irr2Ar,式中,n(r)表示距离原点N
「到"
△「之间的平均粒子数,▼是模拟的体积。
RDF表示的物理含义是:
再空间位置「点周围的体积元中单位体积内发现另一个粒子的几率。
由此看出,RDF表征着结构的无序化程度。
6.2、静态结构因子:
静态结构因子表示式为:
S(k)=|M〔]exp(iK7j)|/N,其中,N是总原子数,K是倒格矢,z;
为原子位置矢量。
对理想晶体而言,SSF为1,而对理想流体,则为0。
6.3、配位数:
所谓配位数,就是原子的第一最近邻原子的个数配位数的多少,表示了该原子
周围的原子分布的密度大小,它经常在一些分析中作为结构分析的一种辅助手段。
7.MD方法的进一步研究方向:
经典分子动力学方法存在两个缺陷:
(1)元胞体积和形状保持不变,限制MD方法的应用:
(2)不适合含有自由电子的系统,对金属等系统的计算的结果不理想。
为克服经典分子动力学方法的局限性,20世纪80年代初Anderson、Parrinello等人发展了可变元胞分子动力学方法°
由于能带论中的密度泛函(DF)方決是研究凝聚态物质电子结构的有效方法,1985年Cai•等人提出了第一原理的MD方法。
但是,进一步拓宽分子动力学方法的应用而显然是一个热门的发展方向。
对MD方法本身而言,最重要的两个要素是初始位形的给定和相互作用势函数的选取。
MD模拟时间足够长,初始条件的选择不会影响计算的结果,但是选择合理的初始条件口J以加快系统趋于平衡,町以节省机时。
影响MD计算结果精确程度的垠主要的因索是相互作用坍的精确性,人们对提高势函数的精确性进行了人量的研究提出了许多势函数的形式,但是,对人多数原子而言,努力找到即准确而形式又不太复杂的势函数,仍然是一个很有挑战性的课题。
我们认为,把第一原理、量子化学分析和参数拟合相结合是势函数研究的最好的方法,NaiixianChen等人就用此方法研究了氯化钠晶体离子之间的相互作用势,并研究了微结构的变化,这是一个值得进一步研究的方向。
P]
8、分子动力学模拟方法的应用:
8.1、氨的自扩散系数的分子动力学模拟研究:
扩散系数是描述传递现象的巫要的流体热物理性质么一,较宽温度和压力范圉的扩散系数在化工生产中具右重要应用价值。
扩散系数的实验测量是一项很困难的工作,而分子动力学(MD)模拟撒获得分子扩故系数的有效方法。
氨,作为大气中最重要的碱性气体之一,它在大气化学、物质传输与沉降过程中发挥着重要作用。
采用MD模拟方法计算氨在较宽温度和压力范圉的自扩散系数,通过将两种不同的控温方式得到的模拟值与大最的实验值进行比较,验证了所采用的模拟方法是合理的。
从而可有效地预测极端条件卜•实验难以测量的扩散系数。
MD模拟采用TINKERv5.0程序包,OPLS・AA全原子力场。
模拟体系包含300个分子。
分子间的作用势采用LennaidJones势,体系中的长程作用力采用Ewald加和的形式°
祚模拟过程中均采用周期性边界条件、Beeman^法求解运动方程。
温控方式分別采用Andersen和Berendsen。
体系的截断半径为limi,时间步长为lfs。
MD模拟在NVT系综进行。
通过模拟发现:
从常温到髙温,采用Berendsen控温方式得到的氨的自扩散系数的模拟值与实验值吻合的很好。
于是,可以釆用分子动力学模拟来代替实验,获得高温高压条件卜•实验难以测量的分子自扩散系数。
另一方面,氨的自扩散系数受温度的影响比受压力的影响更明显。
⑹
8.2、半晶态聚合物的分子动力学模拟:
聚合物材料区别于其它材料的一个典型特征是它具有半结晶形态•当从熔点以上温度开始缓慢冷却到室温时,很多聚合物将发生部分区域结晶、部分保留为非晶态的转变,而形成半晶态的结构。
实验观测到折叠链结构为聚合物结晶的基本形式,由折穂链片晶可形成多层堆穂片晶、球晶、串晶和枝晶等。
针对半晶态聚合物的总体形态,人们提出了多种多样的聚合物结构模型,如缨状微束模型、折亮链模型利折叠链缨状胶粒模型等。
采用Meyer和Mule「P4the提出的粗粒化聚乙烯醇模型,应用分子动力学方法,模拟熔融态聚合物经过缓慢冷却、局部结晶形成半晶态聚合物的过程,研究冷却过程中聚合物的结晶和凝固行为,揭示半晶态聚合物的结构特征及其形成机制。
为半晶态聚合物物理、力学性能的分子模拟研究提供基础。
应用分子动力学方法,模拟熔融态聚合物经过缓慢冷却,局部结晶形成半晶态聚合物的过程。
静态结构因子的演变显示出,在结晶行为的初期小角散射强度的增大先于布拉格峰的出现,这个现象与小角/大角X射线散射实验结果相一致。
模拟得到了半品态的聚合物结构,其中的微小晶区由折叠链构成,它们随机取向分布在非晶态区中,总体形态可用缨状微束结构模型描述。
模拟得到的半晶态聚合物的结晶度为45%,分子链持久长度(即晶区平均厚度)约为6.5nm。
通过比较结晶度和分子链持久长度随温度的变化曲线•发现在不同的冷却阶段具育不同的仃序结构形成机制C从随机缠绕的分子链到折叠链晶区,需要两种形式的结构转变:
分子链的伸展和伸直分子链之间的平行排列。
在从结晶温度到玻璃化温度的凝固转变过程中,存在上述两种形式的结构转变;
而在玻璃化温度之后,材料的活性只能允许调整伸直分子链之间的相对排列位置,但由此也导致了聚合物结晶度明显的增大。
卩】
8.3、多肽抑制剂抑制淀粉质多肽42构象转换的分子动力学模拟和结合自由能计算:
应用分子动力学模拟和结合自由能计算方法研究务肽抑制剂KLVFF、VVIA利LPFFD抑制淀粉质务肽42(Ag42)构象转换的分子机理。
结果表明,三种多肽抑制剂均能够有效抑制Ap42的二级结构由8螺旋向P■折叠的构象转换。
另外,多肽抑制剂降低了Ap42分子内的疏水相互作用,减少了多肽分子内远距离的接触,有效抑制了Ap42的疏水塌缩,从而起到稳定其初始构象的作用。
这些抑制剂与Ap42之间的疏水和静电相互作用(包括氢键)均有利于它们抑制A卩42的构象转换。
此外,抑制剂中的带电氨基酸残基町以增强其和AR42之间的静电相互作用(包括氢键)并降低抑制剂之间的聚集,从而大大增强对Ap42构象转换的抑制能力。
但脯氨酸的引入会破坏多肽的线性结构,从而人人降低其与A042之间的作用力。
上述分子模拟的结果揭示了多肽抑制剂KLVFF.VVIA和LPFFD抑制A|342构象转换的分子机理,对进一步合理设计AB高效短肽抑制剂具有非常重要的理论指导意义。
同
8.4、超临界水中醋酸锌水解反应的分子动力学模拟:
超临界水具冇特殊的溶解度、可变的密度、较低的粘度和表面张力以及较高的扩散特性,利用超临界水热合成法制备纳米ZnO逐渐成为常用的方法之一。
在这种方法中,超临界水(SCW,温度T二647K,压力p二22MPa,密度p二T?
)同时作为反应介质和反应物,含Zn离子的盐如Zn(NO3)2、ZnS04>
Zn(CH3COO)2的水溶液在超临界状态下与碱(如KOH)或纯水先发生水解反应生成Zn(OH)2然后Zn(OH)2脱水得到纳米级的颗料°
利用超临界水热合成还成功制备出不同种类和形貌的纳米颗粒,如MgO、TiO、ZrO2>
SnO2>
ZnS和FesO斗等。
利用分子动力学模拟Zn(CH3COO)2与SCW反应的能量变化,以及产物CHsCOOH在SCW的分布特性,为实验上通过控制温度或压力生成纳米ZnO颗粒提供参考。
利用分子动力学模拟研究Zn(CH3COO)2在SCW中的平衡特征以及水解反应过程的能量变化,定性分析了水分子在Zn(CH3COO)2团簇内部由于静电场的影响而发生解离的可能性,得出以卜•结论:
(1)Zn(CH3COO)2在SCW中易聚集成团,1个Zn周围平均结合5个CH3COO和1个H20满足Zn形成6配位八面体的空间构型的碍要。
Zi0在Zn(CH3COO)2团簇内部和农明的配位情况有所不同,更多的氏0在Zn(CH3COO)2团簇衷而参巧ZM♦配位。
(2)在标准状态,Zn(CH3COO)2水解生成Zn(OH)2和CH3COOH的反应焙变为60.4kJ.moli,从热力学的角度看难以发生,但在超临界条件下,模拟发现上述水解反应使得体系势能降低,表明反应容易发生,并且水解反应发生的同时伴随Zn(CH3COO)2团簇结构的改变。
水解产物0H-容易进入Zn(CH3COO)2团簇内部,富集Zn,而产物CH3COOH扩散进入水相,并倾向分布在水和一气相界面。
(3)定性分析处于Zn(CH3COO)2团簇内部的水分子周国静电势的分布发现,水分子周田分布有很强的负的静电势,容易使水分子极化变形进而解离,促进Zn(CH3COO)2在SCW中的水解反应。
叨
参考文献:
[1]LVerlet.Coniputer"
experiments'
classicalfluids.I・Therniodynamicalpropertiesof
Lennard-Jonesmolecules[J].Phys.Rev;
1967,159:
98-103
[2]RW.Honeycutt.Thepotentialcalculationandsomeapplications[J].MethodsinComputationalPhysics,1970,9:
136-211
[3]N.MetropolisAWRosenbluth,M.N.RosenbluthfA.H.Tellei;
E.Teller.EquationofstatecalculationsbyfastcomputingniachinesQJ.J.ChenLPliys.,1953,21:
1087-1092
[4]M.S.Daw;
M丄Baskes.Embeddedatommethodderivationandapplicationtoimpuiities^suifaces.andotherdefectsinmetals[J].Phys.RevB,1984,29:
8485-8495
[5]崔守鑫,胡海泉,肖效光,黄海军,分子动力学模拟基本原理和主要技术[0].聊
城大学学报,2005.18
(1):
30-34
[6]周昌林,黎多来,李东凯,等,氨的口扩散系数的分子动力学模拟究[J].广东化工,
2012,39(13):
183-184
[7]段芳莉,颜世铛,半晶态聚合物的分子动力学模拟[A].计算物理,2012,29(5):
759