从头计算分子动力学方法及其应用

收稿日期:2004-12-21
基金项目:山东省自然科学基金资助项目(Y2003A01)和石油科技中青年创新基金(04E7038)作者简介:蓝建慧(1979-),女(汉族),山东即墨人,硕士研究生,专业方向为计算物理。

文章编号:100025870(2005)0420143204
综述
从头计算分子动力学方法及其应用
蓝建慧,卢贵武,黄乔松,李英峰,朱　阁
(中国石油大学物理科学与技术学院,山东东营257061)
摘要:从头计算分子动力学方法把密度泛函理论和分子动力学方法有机地结合起来,使电子的极化效应及化学键的本质均可用计算机分子模拟方法进行研究,是目前计算机模拟实验中最先进、最重要的方法之一。

文章简述了从头计算分子动力学方法的基本原理,介绍了该方法在水、水溶液及其他氢键液体的结构与动力学研究中的应用。

关键词:从头计算;密度泛函理论;分子动力学;计算机分子模拟中图分类号:O 35 文献标识码:A
Method of ab initio molecular dynamics and its applications
LAN Jian 2hui ,L U Gui 2wu ,HUAN G Qiao 2song ,L I Y ing 2feng ,ZHU G e
(College of Physics Science and Technology in China U niversity of Pet roleum ,Dongying 257061,China )
Abstract :The ab initio molecular dynamics method ,which combines the density functional theory with the molecular dy 2namics methodology ,made it convenient to study the electronic polarization effects and the nature of the chemical bonds in term of the computer molecular simulation.The method is one of the most im portant and advanced com puter simulation ex 2periment methods.The basic principle of the ab initio molecular dynamics method and its applications in structure and dy 2namics research of liquid water ,aqueous solutions and other hydrogen 2bond liquids were introduced.K ey w ords :ab initio ;density functional theory ;molecular dynamics ;computer molecular simulation
现代凝聚态理论研究应用最普遍的方法之一是
分子动力学(MD )方法。

MD 方法把原子的运动与特定的轨道联系在一起,通过求解原子的牛顿运动方程得到体系的热力学性质。

MD 计算的核心是选择合适的力场,即在系统中引入简单数学模型来描述原子间的结合、弯曲和二面角势及原子间的范德华力和静电作用,并预设实验数据或模型参数进行计算。

利用该方法既可得到原子的运动轨迹,也可在计算过程中对平衡或非平衡系统的微观细节进行研究,因此被广泛应用于研究各种物理、化学问题,在处理凝聚态体系时获得了极大的成功[1]。

但是MD 方法是基于力场的,因此存在严重的缺陷。

首先它忽略了电子极化效应。

后来有人提出了校正该缺陷的极化模型,但极化模型只适用于特定的问题,目前还没有得到广泛的应用。

此外,MD 方法无法描述化学键形成或断裂的本质问题,若采用经验价键方法或其他的半经验方法作近似处理,即使对不
同反应设置不同的参数,模拟结果仍可能偏离反应路径[2]。

为了解决上述问题,从头计算分子动力学方法(ab initio molecular dynamics ,A IMD )被提出。

本文中重点介绍从头计算分子动力学方法的基本原理及其在水、水溶液和其他氢键液体研究中的应用。

1　从头计算分子动力学方法
从头计算分子动力学(A IMD )方法主要基于以下3个假设:(1)忽略系统的核量子效应;(2)认为系统满足轨道近似(即单电子近似);(3)认为系统满足绝热近似。

其中电子基态本征函数和本征值的计算是A IMD 的核心内容。

电子基态计算属于复杂的量子多体问题,需引入密度泛函理论(densisty functional theory ,DF T )以简化计算量,把复杂的多体问题转化为一组自洽的单电子轨道方程,即K ohn 2Sham 方程,并根据电子和原子核的相互作用对电子密度的影响程度,对交
2005年　第29卷 石油大学学报(自然科学版) Vol.29　No.4　第4期 Journal of the University of Petroleum ,China Aug.2005
换势采用局域密度近似[1](local density approxima 2
tion ,LDA )或广义梯度近似(general gradient ap 2proximation ,GG A ),从而方程可解。

这也是目前凝聚态物理中计算电子结构普遍采用的方法。

电子基态计算不仅得到原子间相互作用势,也为分子动力学研究提供了精确的力场,
这种将密度函数理论和分子动力学结合起来的方法,
是目前从头计算分子动力学中普遍采用的算法之一(即Car 2Parrinello ,CP 算法)[2]。

该算法的基本步骤见图1。

由图可见,A IMD 模拟可归结为电子系统动力学计算、电子和离子耦合系统动力学计算两个过程。

图1　从头计算分子动力学模拟流程
第一过程主要有以下几个步骤:①根据电子系统动力学计算原子初始构型的基态单电子轨道波函数。

通常采用尝试法,先用少量的平面波基函数将分子轨道展开,按矩阵对角化技术得到分子轨道本征波函数,再以按照扩展平面波基函数展开的分子轨道作尝试波函数进行虚拟动力学模拟。

②根据密度泛函理论计算电子体系总能量E 和离子间的相互作用势V 。

体系的基态能量对应于E 的极小值,是电子基态密度的泛函。

离子间的相互作用势V
(即电子和离子耦合系统的势能)以及V 对离子坐
标的微分根据密度泛函理论和Hellman 2Feynman 定理计算。

③计算电子波函数正交化产生的“虚拟力”,求解运动方程,得到新的电子基态波函数。

如果运动方程中未包含由波函数正交化产生的“虚拟力”,需要对新的波函数正交化。

④将新的波函数作
为t =t +Δt 时刻的输入,进行下一步的电子动力学模拟,直至t >t max (t max 即A IMD 模拟的最大设定时间)。

第二过程算法步骤如下:得到初始本征波函数
后,设定虚拟系统广义拉格朗日欧拉方程中的折合
质量等参量,使电子的运动标度远小于离子的运动标度,从而保证电子随着离子位形的改变尽量趋于基态(即电子在等能面E 附近作简谐振动),这样无论高温还是低温条件下,整个耦合系统就处在BO (Born 2Oppenheimer )面上,然后开始电子和离子耦合系统的动力学模拟。

根据第一过程的计算结果即可计算离子运动方程中的“虚拟力”,求解离子运动方程并得到t =t +Δt 时刻的离子坐标。

2　从头计算分子动力学方法的应用
从头计算分子动力学方法在凝聚态体系的微观几何结构和电子结构以及反应动力学的研究中应用广泛,本文中主要介绍该方法在液态水及其他氢键液体的结构与动力学研究中的应用。

2.1　液态水的从头计算分子动力学模拟
由于300K 水的特殊性,对该温度下水的结构的A IMD 研究[3-7]比较多。

文献[8]中研究了由64个水分子构成的小系统,模拟布里渊区中心Γ点,采用Troullier 2Martins 型原子赝势,电子相关效应用BL YP 方法处理。

运行周期为11ps ,时间步长为0.17fs ,得到的径向分布函数与实验结果一致。

文献[7]采用Perdew 交换关联函数,当运行周期为715ps 时,得到的O -O 径向分布函数与实验吻合。

A IMD 模拟还被用于计算水的红外光谱和平均偶极矩。

文献[4]和[9]中分别对含32个水分子和64个水分子的体系进行了研究,采用BL YP 函数且运行周期为10～12ps 时,水的红外光谱和平均偶极矩的计算结果与实验结果吻合较好。

另外,将路径取样法[10,11]和A IMD 技术结合起来,可以阐述水的自身离解以及H 3O -/OH -离子对的形成问题。

应当指出,由于A IMD 模拟中系统和基组大小的影响[3-7]以及直接从中子和X 2射线衍射数据中提取这些量尚存在困难,通过A IMD 方法得到300K 水的准确径向分布函数仍是一个具有挑战性的课题。

2.2　液态氨的从头计算分子动力学模拟
液氨化学类似于水化学,水中的反应类型在液氨中基本上都可以见到。

液氨分子中氮原子的电负性比水中氧原子的小,相应的非键电子更容易被利用,故被广泛地应用于大量普通的有机反应和与金属反应的溶液中。

其结构因子和径向分布函数已于最近通过中子衍射技术在实验中得到确认,通过A IMD 研究结果和纯液态氨实验数据的比较,可以
・441・石油大学学报(自然科学版) 2005年8月
对A IMD方法做出评价。

目前,基于CP运动方程的A IMD模拟被应用于含32个氨分子、箱长为11.27!且满足周期性边界条件的样品,实验温度分别为260K和273 K[12,13]。

交换关联函数用BL YP GG A方法[14,15]处理,核心电子用Troullier2Martins型赝势处理。

系统平衡态可保持的最长时间为2.2ps,运行周期为6.0ps,时间步长为0.12fs。

260K液氨的结构因子和径向分布函数的计算结果与实验结果[16]吻合极好。

另外根据计算[17]确定的自扩散常数为1.1×10-4cm2・s-1,比较接近实验值1.0×10-4cm2・s-1。

2.3　300K液态甲醇的从头计算分子动力学模拟
另一种重要的氢键液体是甲醇(CH3OH)。

和液氨一样,甲醇作为溶剂被广泛地应用于许多普通的有机反应中,尤其在新兴的燃料电池技术中发挥着重要作用,是一种重要的工业用液体。

利用中子衍射技术,液体甲醇的结构最近得到确认[18,19]。

与从头计算相比,利用该技术更容易得到结构因子和径向分布函数。

Morrone等[20]对液体甲醇中质子的输运特性进行了A IMD研究,模拟采用含32个甲醇分子、周期性箱的尺寸为12.93!的样品,交换与关联用B IYP GG A方法处理。

分别采用两种不同的模拟方案进行研究。

第一方案:选用25个平面波并采用超软赝势法[21],运行周期为20ps;第二方案:采用量子力学/分子力学(QM/MM)方法[22], OH基在QM级别上用BL YP函数处理[14,15],CH3基在MM级别上用AMB ER力场处理。

研究结果发现,对于结构因子,QM/MM模型、超软赝势法与实验三者吻合较好,而对于OO,OH和CC的径向分布函数,QM/MM模型比超软赝势法好得多。

2.4　K OD溶液的从头计算分子动力学模拟
了解水环境中由于H-或OH-的增加而产生的电荷转移行为在酸碱化学、生物工程、燃料电池、皂化和工业催化剂等领域至关重要。

含水合氢离子H3O-和氢氧根离子OH-的水溶液是良好的质子导体,其输运过程的微观细节目前仍存在争议。

部分争议在于难以找到能证明其特殊机制的明确的实验证据。

高水平A IMD模拟在揭开这些机制的细节方面可以发挥特殊作用。

在文献[23]的研究中,1.5 mol/L的溶液含32个水分子和一个孤立的KOD分子,置于尺寸为10.25!的周期性立方箱中;13 mol/L的溶液由27个水分子和8个孤立的KOD分子构成,置于尺寸为10.15!的周期性立方箱中。

为了减小轻氢原子核量子效应的影响,用KOD代替KOH。

交换关联采用BL YP方法处理,核心电子用Troullier2Martins赝势和一个钾的半核赝势处理。

模拟周期为10ps,采用600au虚拟电子群,时间步长为5ps。

此模拟长度可容纳2ps的平衡态。

模拟计算了偶极矩并给出红外光谱。

由于模拟长度相对较短,采用最大熵技术[24]提取红外光谱。

比较计算得到的KOD水溶液的红外光谱和14mol/L KOH溶液的实验测量结果[2]发现,前者再现了14 mol/L KOH溶液的实验所得光谱的所有特征。

3　结束语
介绍了结合有限温度分子动力学(MD)与虚拟电子结构计算的从头计算分子动力学模拟方法,以及液态水、液态氨、液态甲醇和KOD水溶液中的A IMD模拟应用。

此外,A IMD方法在液体结构、酸化学、工业和生物催化剂以及材料学等研究领域也得到广泛应用,揭示了许多经验模型无法解释的现象,同时启发人们对一些实验数据做出新的解释,设计新的实验研究方案。

应当指出,从头计算分子动力学方法使基于密度泛函理论的第一性原理计算直接应用于统计力学模拟成为可能,极大地扩展了计算机分子模拟的广度和深度,成为计算机分子模拟最先进和最重要的方法之一。

但与以力场为基础的传统计算机模拟方法相比,A IMD方法的计算量显著增加。

传统计算机分子模拟可应用于时间标度在几十个纳秒数量级、原子数约为104～106的大分子体系,A IMD计算则主要应用于时间标度在几十个皮秒数量级而原子数在几十到几百的小分子体系。

这些特点在实际应用中需格外注意。

参考文献:
[1]　赵宇军,姜明,曹培林.从头计算分子动力学[J].物理
学进展,1998,18(1):47-75.
ZHAO Yu2jun,J IAN G Ming,CAO Pei2lin.Ab initio
molecular dynamics[J].Progress in Physics,1998,18
(1):47-75.
[2]　MAR K E Tuckerman.Ab initio molecular dynamics:ba2
sic concepts,current trends and novel applications[J].J
Phys,2002,14(50):1297-1355.
[3]　LAASON EN K,SPRIK M,PARRIN ELLO M,et al.
Ab initio liquid water[J].J Chem Phys,1993,99(11):
9080-9089.
[4]　SILV ESTRELL I P L,BERNASCON I M and PAR2
・
5
4
1
・
第29卷　第4期 蓝建慧等:从头计算分子动力学方法及其应用
RIN ELLO M.Ab initio infrared s pectrum of liquid wa 2ter [J ].Chem Phys Lett ,1997,277(5-6):478-482.[5]　TROU T B L and PARRIN ELLO M.X 2ray spectroscopic
and quantum 2chemical study of carbon tubes produced in arc 2discharge [J ].Chem Phys Lett ,1998,288(3-4):343-347.
[6]　SPRIK putation of the p K of liquid water using
coordination constraints [J ].Chem Phys ,2000,258(2-3):139-150.
[7]　SCHWEG L ER E ,G ALL I G ,GYGI F ,et al .Dissocia 2
tion of water under pressure [J ].Phys Rev Lett ,2001,87(26):265-501.
[8]　IZV EKOV S and VO TH G.A Car 2Parrinello molecular
dynamics simulation of liquid water :new results [J ].J Chem Phys ,2002,116(23):10372-10376.
[9]　SILV ESTRELL I P L and PARRINDLLO M.Water
molecule dipole in the gas and in the liquid phase [J ].Phys Rev Lett ,1999,82(16):3308-3311.
[10]　DELLA G O C ,BOL HU IS P G ,CSAJ K A F S ,et al .
Transition path sampling and the calculation of rate con 2stants[J ].J Chem Phys ,1998,108(5):1964-1977.[11]　DELLA GE C ,BOL HU IS P G and CHANDL ER D.On
the calculation of reaction rate constants in the transition path ensemble [J ].J Chem Phys ,1999,110(14):6617-6625.
[12]　DIRAISON M ,MART YNA G J and TUCKERMAN
M E.Simulation studies of liquid ammonia by classical ab initio ,classical ,and path 2integral molecular dynam 2ics[J ].J Chem Phys ,1999,111(3):1096-1103.
[13]　L IU Y and TUCKERMAN M E.Protonic defects in
hydrogen bonded liquids :structure and dynamics in am 2monia and comparison with water[J ].J Phys Chem B ,2001,105(28):6598-6610.
[14]　BECKE A D.Density 2functional exchange 2energy ap 2
proximation with correct asymptotic behavior[J ].Phys Rev A ,1988,38(6):3098-3100.
[15]　L EE W Y C and PARR R C.Development of the Colle 2
Salvetti correlation 2energy formula into a functional of the electron density[J ].Phys Rev B ,1988,37(2):785-789.
[16]　RICCI M A ,NARDON E M ,RICCI F P ,et al .Mi 2
croscopic structure of low tem perature liquid ammonia :a neutron diffraction experiment [J ].J Chem Phys 1995,102(19):7650-7655.
[17]　O ’REILL Y D E ,PETERSON E M and SCHEIE C E.
Self 2diffusion in liquid ammonia and deuteroammonia [J ].J Chem Phys ,1973,58(10):4072-4075.
[18]　Y ANA GUCHI T ,HIDA K A K and SOPER A K.The
structure of liquid methanol revisited :a neutron diffrac 2tion experiment at -80℃and +25℃[J ].Mol Phys ,1999,96(8):1159-1168.
[19]　Y ANA GUCHI T ,HIDA K A K and SOPER A K.Erra 2
tum :the structure of liquid methanol revisited :a neu 2tron diffraction experiment at -80℃and +25℃[J ].Mol Phys ,1999,97(4):603-605.
[20]　MORRONE J A and TUCKERMAN M E.Ab initio m olec 2
ular dynamics study of proton m obility in liquid methanol[J ].J Chem Phys ,2002,117(9):4403-4413.
[21]　LAASON EN K ,PASQUARELLO A ,CAR R ,et al.
Car 2Parrinello molecular dynamics with Vanderbilt ul 2trasoft pseudopotentials [J ].Phys Rev B ,1993,47(16):10142-10153.
[22]　TSUCHIDA E ,K ANADA Y and TSU K ADA M.Den 2
sity 2functional study of liquid methanol [J ].Chem Phys Lett ,1999,311(3):236-240.
[23]　ZHU Z and TUCKERMAN M E.Ab initio molecular
dynamics investigation of the concentration dependence of charged defect transport in basic solutions via calcula 2tion of the infrared s pectrum [J ].J Phys Chem B ,2002,106(33):8009-8018.
[24]　NAILON G W.The Maximum Entropy Method [M ].
New Y ork :S pringer ,1997.
(编辑　陈淑娴)
“油田开发过程中油层保护及解堵技术研究”达国际领先水平
由中国石油大学石油工程学院蒋官澄博士为首的课题组与中国石化胜利油田有限公司孤东采油厂共同
完成的“油田开发过程中油层保护及解堵技术研究”,近日通过了由山东省科技厅组织的专家鉴定。

鉴定委员会一致认为,该项研究达到了同类研究的国际领先水平。

目前,该项研究成果已在孤东油田进行了1877井次的现场试验,提高了施工成功率和有效期,保护了油层,共创直接经济效益29664.67万元,总的投入产出比达到1∶2.7。

(摘自中国石油大学校园网)
・641・石油大学学报(自然科学版) 2005年8月。