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

2005年第29卷石油大学学报(自然科Vol. 29 No 4第4期Journal of the University of Petroleum, ChinaAug.2005文章编号:1000-5870(2005)04014304综述t从头计算分子动力学方法及其应用形1蓝建慧,卢贵武,黄乔松,李英峰,朱阁(中国石油大学物理科学与技术学院,山东东营257061)摘要:从头计算分子动力学方法把密度泛函理论和分子动力学方法有机地结合起来,使电子的极化效应及化学键的本质均可用计算机分子模拟方法进行研究,是目前计算机模拟实验中最先进、最重要的方法之一。文章简述了从头计算分子动力学方法的基本原理,介绍了该方法在水、水溶液及其他氢键液体的结构与动力学研究中的应用关键词:从头计算;密度泛函理论;分子动力学;计算机分子模拟中图分类号:O35文献标识码:AMethod of ab initio molecular dynamics and its applicationsLAN Jian-hui, LU Gui-wu, HUANG Qiaosong, LI Ying-feng, ZHU GeCollege of Physics Science and Technology in China University of Petroleum, Dongying 257061, China)Abstract: The ab initio molecular dynamics method, which combines the density functional theory with the molecular dynamics methodology, made it convenient to study the electronic polarization effects and the nature of the chemical bonds interm of the computer molecular simulation. The method is one of the most important and advanced computer simulation experiment methods. The basic principle of the ab initio molecular dynamics method and its applications in structure and dynamics research of liquid water, aqueous solutions and other hydrogen-bond liquids were introducedKey words: ab initio density functional theory: molecular dynamics; computer molecular simulation现代凝聚态理论研究应用最普遍的方法之一是同反应设置不同的参数,模拟结果仍可能偏离反应分子动力学(MD)方法。MD方法把原子的运动与路径{2。为了解决上述问题,从头计算分子动力学特定的轨道联系在一起,通过求解原子的牛顿运动方法( ab initio molecular dynamics,AIMD)被提出。方程得到体系的热力学性质。MD计算的核心是选本文中重点介绍从头计算分子动力学方法的基本原择合适的力场,即在系统中引入简单数学模型来描理及其在水、水溶液和其他氢键液体研究中的应用。述原子间的结合、弯曲和二面角势及原子间的范德华力和静电作用并预设实验数据或模型参数进行1从头计算分子动力学方法计算。利用该方法既可得到原子的运动轨迹,也可从头计算分子动力学(AIMD)方法主要基于以在计算过程中对平衡或非平衡系统的微观细节进行下3个假设:(1)忽略系统的核量子效应;(2)认为系研究,因此被广泛应用于研究各种物理、化学问题,统满足轨道近似(即单电子近似);(3)认为系统满足在处理凝聚态体系时获得了极大的成功。但是绝热近似。其中电子基态本征函数和本征值的计算MD方法是基于力场的,因此存在严重的缺陷。首是AIMD的核心内容。先它忽略了电子极化效应。后来有人提出了校正该电子基态计算属于复杂的量子多体问题,需引缺陷的极化模型,但极化模型只适用于特定的问题,人密度泛函理论( density functional theory,DFT)以目前还没有得到广泛的应用。此外,MD方法无法简化计算量,把复杂的多体问题转化为一组自洽的描述化学键形成或断裂的本质问题,若采用经验价单键方法或其他的半经验方法作近似处理即使对不和原中国煤化T1方程并根据电子的影响程度,对交CNMHG收稿日期:2004-12-21基金项目:山东省自然科学基金资助项目(Y2003A01)和石油科技中青年创新基金(04E7038)作者简介:蓝建慧(1979-),女(汉族),山东即墨人,硕士研究生,专业方向为计算物理。144石油大学学报(自然科学版)2005年8月换势采用局域密度近似( local density approxima-后设定虚拟系统广义拉格朗日-欧拉方程中的折合tion,LDA)或广义梯度近似( general gradient ap-质量等参量,使电子的运动标度远小于离子的运动proximation,GGA),从而方程可解。这也是目前凝标度,从而保证电子随着离子位形的改变尽量趋于聚态物理中计算电子结构普遍采用的方法。电子基基态(即电子在等能面E附近作简谐振动),这样无态计算不仅得到原子间相互作用势,也为分子动力论高温还是低温条件下,整个耦合系统就处在BO学研究提供了精确的力场,这种将密度函数理论和( Born-Oppenheimer)面上,然后开始电子和离子耦分子动力学结合起来的方法,是目前从头计算分子合系统的动力学模拟。根据第一过程的计算结果即动力学中普遍采用的算法之一(即 Car-Parrinello,可计算离子运动方程中的“虚拟力”,求解离子运动CP算法)2。该算法的基本步骤见图1。由图可方程并得到t=t+△t时刻的离子坐标。见,AIM模拟可归结为电子系统动力学计算、电子和离子耦合系统动力学计算两个过程。2从头计算分子动力学方法的应用始坐标从头计算分子动力学方法在凝聚态体系的微观YES几何结构和电子结构以及反应动力学的研究中应用计算基态单电子波函数广泛,本文中主要介绍该方法在液态水及其他氢键液体的结构与动力学研究中的应用。计算离子间相互作用势V及电子系统总能量E2.1液态水的从头计算分子动力学模拟计算离子虚拟力计算电子的虚拟力由于300K水的特殊性,对该温度下水的结构的AMD研究13-7比较多。文献[8]中研究了由64求解离子运动方程求解电子运动方程个水分子构成的小系统,模拟布里渊区中心点采用 Troullier- Martins型原子赝势,电子相关效应新的离子坐标新的体系波函数用BLYP方法处理。运行周期为11ps,时间步长为0.17fs,得到的径向分布函数与实验结果一致。文t=t+△t献[7]采用 Perdew交换关联函数,当运行周期75ps时,得到的O-O径向分布函数与实验吻合。图1从头计算分子动力学模拟流程AMD模拟还被用于计算水的红外光谱和平均第一过程主要有以下几个步骤:①根据电子系偶极矩。文献[4]和[9]中分别对含32个水分子和统动力学计算原子初始构型的基态单电子轨道波函64个水分子的体系进行了研究采用BLYP函数且数。通常采用尝试法,先用少量的平面波基函数将运行周期为10~12p时,水的红外光谱和平均偶分子轨道展开,按矩阵对角化技术得到分子轨道本极矩的计算结果与实验结果吻合较好。另外,将路征波函数,再以按照扩展平面波基函数展开的分子径取样法10和AMD技术结合起来,可以阐述水轨道作尝试波函数进行虚拟动力学模拟。②根据密的自身离解以及H3OOH离子对的形成问题。度泛函理论计算电子体系总能量E和离子间的相应当指出,由于AIMD模拟中系统和基组大小互作用势v。体系的基态能量对应于E的极小值,的影响3-以及直接从中子和x射线衍射数据中是电子基态密度的泛函。离子间的相互作用势v提取这些量尚存在困难通过AIMm方法得到300(即电子和离子耦合系统的势能)以及V对离子坐K水的准确径向分布函数仍是一个具有挑战性的课标的微分根据密度泛函理论和 Hellman- Feynman定题。理计算。③计算电子波函数正交化产生的“虚拟2.2液态氨的从头计算分子动力学模拟力”,求解运动方程,得到新的电子基态波函数。如液氨化学类似于水化学,水中的反应类型在液果运动方程中未包含由波函数正交化产生的“虚拟氨中基本上都可以见到。液氨分子中氮原子的电负力”需要对新的波函数正交化。④将新的波函数作性比才V山中国煤化工建电子更容易被利为t=t+△时刻的输人,进行下一步的电子动力用故CNMH(有机反应和与金学模拟直至t>tm(tm即AIMD模拟的最大设属反应的溶液中。其结构因子和径向分布函数已于定时间)。最近通过中子衍射技术在实验中得到确认,通过第二过程算法步骤如下:得到初始本征波函数AIMD研究结果和纯液态氨实验数据的比较,可以第29卷第4期蓝建慧等:从头计算分子动力学方法及其应用145对AIMD方法做出评价。为了减小轻氢原子核量子效应的影响,用KOD代目前,基于CP运动方程的AMD模拟被应用替KOH。交换关联采用BLYP方法处理,核心电子于含32个氨分子、箱长为1127A且满足周期性边用 Troullier-Martins腰势和一个钾的半核赝势处界条件的样品,实验温度分别为260K和273理。模拟周期为10ps,采用600au虚拟电子群时K12.13]。交换关联函数用 BLYP GGA方法4.1处间步长为5ps此模拟长度可容纳2ps的平衡态。理,核心电子用 Troullier Martins型赝势处理。系模拟计算了偶极矩并给出红外光谱。由于模拟长度统平衡态可保持的最长时间为22p,运行周期为相对较短,采用最大熵技术4提取红外光谱。比较60ps,时间步长为0.12fs。260K液氨的结构因计算得到的KOD水溶液的红外光谱和14moL子和径向分布函数的计算结果与实验结果1)吻合KOH溶液的实验测量结果2发现,前者再现了14极好。另外根据计算1确定的自扩散常数为1.1 mol/L Koh溶液的实验所得光谱的所有特征×10-4cm2·s-1,比较接近实验值1.0×10-4cm23结束语2.3300K液态甲醇的从头计算分子动力学模拟介绍了结合有限温度分子动力学(MD)与虚拟另一种重要的氢键液体是甲醇(CH3OH)。和电子结构计算的从头计算分子动力学模拟方法,以液氨一样,甲醇作为溶剂被广泛地应用于许多普通及液态水、液态氨、液态甲醇和KOD水溶液中的的有机反应中,尤其在新兴的燃料电池技术中发挥AM模拟应用。此外,AMD方法在液体结构、酸着重要作用,是一种重要的工业用液体。利用中子化学、工业和生物催化剂以及材料学等研究领域也衍射技术,液体甲醇的结构最近得到确认18.9。与得到广泛应用揭示了许多经验模型无法解释的现从头计算相比,利用该技术更容易得到结构因子和象,同时启发人们对一些实验数据做出新的解释,设径向分布函数。 Morrone等201对液体甲醇中质子计新的实验研究方案的输运特性进行了AIMD研究,模拟采用含32个应当指出,从头计算分子动力学方法使基于密甲醇分子、周期性箱的尺寸为12.93A的样品,交换度泛函理论的第一性原理计算直接应用于统计力学与关联用 BIYP GGA方法处理。分别采用两种不模拟成为可能极大地扩展了计算机分子模拟的广同的模拟方案进行研究。第一方案:选用25个平面度和深度,成为计算机分子模拟最先进和最重要的波并采用超软赝势法(21,运行周期为20ps;第二方方法之一。但与以力场为基础的传统计算机模拟方案:采用量子力学/分子力学(QMMM)方法2),法相比,AIMD方法的计算量显著增加。传统计算OH基在QM级别上用BLYP函数处理14.1,CH3机分子模拟可应用于时间标度在几十个纳秒数量基在MM级别上用 AMBER力场处理。研究结果级原子数约为104~106的大分子体系,AMD计发现,对于结构因子,QMMM模型超软赝势法与算则主要应用于时间标度在几十个皮秒数量级而原实验三者吻合较好,而对于OO,OH和CC的径向子数在几十到几百的小分子体系。这些特点在实际分布函数,QM^MM模型比超软赝势法好得多应用中需格外注意。2.4KOD溶液的从头计算分子动力学模拟了解水环境中由于H或OH的增加而产生的参考文献电荷转移行为在酸碱化学、生物工程燃料电池、皂[1]赵宇军姜明,曹培林从头计算分子动力学J物理化和工业催化剂等领域至关重要。含水合氢离子学进展,1998,18(1):47H3O和氢氧根离子OH的水溶液是良好的质子导ZHAO Yu-jun, JIANG Ming, CAO Pei-lin. Ab initio体,其输运过程的微观细节目前仍存在争议。部分molecular dynamics[J]. Progress in Physics, 1998, 18争议在于难以找到能证明其特殊机制的明确的实验(1):47-75证据。高水平AMD模拟在揭开这些机制的细节2] MARK E Tuckerman. Ab initio molecular dynamics:tend novel applications [J].J方面可以发挥特殊作用。在文献[23]的研究中,1.5中国煤化工5mL的溶液含32个水分子和一个孤立的KOD分[3C N M H GRRINELLO M, et al子,置于尺寸为10.25A的周期性立方箱中;13Ab initio liquid water [J]. J Chem Phys, 1993, 99(11)mo/L的溶液由27个水分子和8个孤立的KOD分9080-9089子构成,置于尺寸为10.15A的周期性立方箱中。[4] SILVESTRELLI P L, BERNASCONI M and PAR石油大学学报(自然科学版)2005年8月RINELLO M. Ab initio infrared spectrum of liquidSalvetti correlation-energy formula into a functional ofter[J. Chem Phys Lett,1997,277(5-6):478-482the electron density[j]. Phys Rev B, 1988, 37(2): 785[5] TROUT B L and PARRINELLO M. X-ray spectroscopicnd quantum-chemical study of carbon tubes produced in [16] RICCI M A, NARDONE M, RICCI F P, et al.Miarc-discharge [J]. Chem Phys Lett, 1998, 288 (3-4)croscopic structure of low temperature liquid ammonia: aneutron diffraction experiment [J]. J Chem Phys 1995[6 SPRIK M. Computation of the pK of liquid water using102(19):7650-7655coordination constraints [J]. Chem Phys, 2000, 258(2-[17]O'REILLY D E, PETERSON E M and SCHEIE CE.Self-diffusion in liquid ammonia and deuteroammonia[7] SCHWEGLER E, GALLI G, GYGI F, et al. Dissocia[J]. J Chem Phys,1973,58(10):4072-4075tion of water under pressure [J]. Phys Rev Lett, 2001, [18] YANAGUCHI T, HIDAKA K and SOPER AK.The87(26):265-501structure of liquid methanol revisited: a neutron diffrac[ 8] IZVEKOV S and VOTH G. A Car- Parrinello moleculartion expenment at-80℃and+25℃[J].Moldynamics simulation of liquid water: new results [J].JPhys,1999,96(8):1159-1168Chem Phys,202,116(23):10372-10376[19] YANAGUCHI T, HIDAKA K and SOPER AK. Erra9] SILVESTRELLI P L and PARRINDLLO M. Watum: the structure of liquid methanol revisited a neu-[J]ron diffraction experiment at-80℃and+25℃[J]hys Rev Lett,99,82(16):3308-3311Mol Phys,1999,97(4):603-605[10] DELLAGO C, BOLHUIS P G, CSAJKA F S, et al. [20] MORRONE J A and TUKERMANME. Ab initio mlecTransition path sampling and the calculation of rate con-ular dynamics study of proton bility in liquid methanol[J]stants[J]. J Chem Phys,1998,108(5):1964-1977J Chem Phys,2002,117(9:4403-4413[11] DELLAGE C, BOLHUIS P G and CHANDLER D On [21] LAASONEN K, PASQUARELLO A, CAR R, et althe calculation of reaction rate constants in the transitionCar-Parrinello molecular dynamics with Vanderbilt ulpath ensemble [J]. J Chem Phys, 1999, 110(14): 6617rasoft pseudopotentials [J]. Phys Rev B, 1993, 47(16):10142-1015[12] DIRAISON M, MARTYNA G J and TUCKERMAN [22] TSUCHIDA E, KANADA Y and TSUKADA M.Den-ME. Simulation studies of liquid ammoniasity-functional study of liquid methanol[J]. Chem Phab initio, classical, and path-integral molecular dynam-tt,1999,311(3):236-240ics[J]. j Chem Phys,199,01(3):1096-1103[23] ZHU Z and TUCKERMAN M E. Ab initio molecular[13] LIU Y and TUCKERMAN M E. Protonic defectsdynamics investigation of the concentration dependencehydrogen bonded liquids: structure and dynamics in am-ofmonia and comparison with water[J]. J Phys Chem B,tion of the infrared spectrum [J]. J Phys Chen B2001,105(28):6598-66102002,106(33):8009-8018[14] BECKE A D. Density-functional exchange-energy ap- [24] NAILONG W. The Maximum Entropy Method [M]proximation with correct asymptotic behavior[J]. PhysNew York: Springer, 1997RevA,1988,38(6):3098-3100(编辑陈淑娴」[15] LEE W Y C and PARR R C. Development of the Colle-油田开发过程中油层保护及解堵技术研究”达国际领先水平由中国石油大学石油工程学院蒋官澄博士为首的课题组与中国石化胜利油田有限公司孤东采油厂共同完成的油田开发过程中油层保护及解堵技术研究”,近日通过V凵中国煤化工家鉴定。鉴定委员会一致认为,该项研究达到了同类研究的国际领先水平。CNMHG目前,该项研究成果已在孤东油田进行了1877井次的现场试验,提高了施工成功率和有效期,保护了油层,共创直接经济效益29664.67万元,总的投入产出比达到1:2.7。(摘自中国石油大学校园网)