基于分子動力學(xué)的常用力場、算法及結(jié)果分析

基于分子動力學(xué)的常用力場、算法及結(jié)果分析摘要：分子動力學(xué)模擬是隨著計算機技術(shù)發(fā)展而興起的一項新技術(shù)，結(jié)合了物理，數(shù)學(xué)和化學(xué)的模擬方法。分子動力學(xué)模擬具有溝通宏觀特性與微觀結(jié)構(gòu)的作用，對于許多在理論分析和實驗都難以了解和觀察到的現(xiàn)象能給出一定的微觀解釋，可以分析其中的物理現(xiàn)象和本質(zhì)。本文主要對分子動力學(xué)中常用力場，算法及結(jié)果分析方法進行綜述。最后基于計算機硬件的發(fā)展與模擬算法的優(yōu)化，對該領(lǐng)域未來發(fā)展進行展望。關(guān)鍵詞：分子動力學(xué)；力場；算法；結(jié)果分析；引言早在20世紀50年代末，人們就已經(jīng)開始進行計算模擬實驗。一開始，這些計算實驗的對象只是一些小的模型體系，比如幾個硬球的集合。但早在那時，Alder和Wainwright[1]就已經(jīng)預(yù)見到,這樣的計算機模擬實驗將會成為聯(lián)結(jié)宏觀實驗現(xiàn)象和其微觀本質(zhì)的重要橋梁。盡管解析理論進行模擬也可以部分解釋實驗，但是往往附帶了很多近似，而計算機模擬能夠提供體系的詳盡信息，具有更大的優(yōu)勢。由于宏觀和微觀在尺度上是截然不同的，為了能夠在兩者之間建立切實的關(guān)聯(lián)，即微觀模擬體系能夠反映相應(yīng)的宏觀實驗現(xiàn)象，需要通過周期性邊界條件對模擬對象體系進行周期性的復(fù)制，以避免在實際中并不存在的邊緣效應(yīng)。分子動力學(xué)模擬主要是依靠牛頓力學(xué)來模擬分子體系的運動，以在由分子體系的不同狀態(tài)構(gòu)成的系統(tǒng)中抽取樣本，從而計算體系的構(gòu)型積分，并以構(gòu)型積分的結(jié)果為基礎(chǔ)進一步計算體系的熱力學(xué)量和其他宏觀性質(zhì)。通過分子模擬技術(shù)，人們可以在分子、原子尺度上對所關(guān)心材料的物理、化學(xué)性質(zhì)以及與之相對應(yīng)的各種現(xiàn)象進行模擬，摒棄傳統(tǒng)材料科學(xué)研究中被稱之為“愛迪生方式”的炒菜式研究，利用計算機對材料的組成、性質(zhì)、制備工藝、加工工藝等環(huán)節(jié)進行虛擬實驗，獲得大量數(shù)據(jù)，再與少量驗證性實驗進行對照，分析物理現(xiàn)象和本質(zhì)，從而大大的降低材料的研發(fā)成本和時間。這種新材料的設(shè)計模式已經(jīng)很大程度上取代了傳統(tǒng)的新材料設(shè)計方式，成為新材料設(shè)計的主要方式之一。2011年，美國啟動一項價值超過5億美元的“先進制造業(yè)伙伴關(guān)系”計劃，而以材料計算與模擬為主的“材料基因組計劃”則是其中重要內(nèi)容之一[2]。目前，對分子動力學(xué)模擬的力場、算法及結(jié)果分析方法很多，但基本上是分散在英文文獻中，由于翻譯和理解上的差異，影響了國內(nèi)從事分子模擬科研工作者的理解交流。本文整理了目前分子動力學(xué)中主流的力場、算法及結(jié)果分析，并在最后對該領(lǐng)域進行了展望。1力場力場簡介分子動力學(xué)模擬是計算龐大復(fù)雜系統(tǒng)的有效方法，它以力場為依據(jù)，力場的完備與否決定計算的可靠程度。分子的總能量是動能與勢能之和，分子的勢能通常表示為簡單的幾何坐標(biāo)的函數(shù)。一般勢能中包括：（1）范德華力，與能量有關(guān)的非鍵相互作用交叉能量項，（2）構(gòu)成分子的各個化學(xué)鍵在鍵軸方向上的伸縮運動所引起的能量變化，（3）鍵角變化引起的分子能量變化，（4）單軸旋轉(zhuǎn)引起分子骨架扭曲所產(chǎn)生的能量變化，（5）離平面振動項，共平面原子的中心原子離平面小幅振動的勢能，（6）庫倫作用項，帶電荷粒子間存在的靜電吸引或排斥作用的勢能。力場可以看作是勢能面的經(jīng)驗表達式，它是分子動力學(xué)模擬的基礎(chǔ)、力場是通過原子位置計算體系能量的，與之前的量子力學(xué)方法相比，大大節(jié)約了計算時間，可用于計算包含上萬粒子數(shù)目的體系。勢能函數(shù)在大多數(shù)情況下將描述分子幾何形變最大程度地簡化為僅僅使用簡諧項和三角函數(shù)來實現(xiàn)，而非鍵原子之間的相互作用，則只采用庫倫相互作用和蘭納-瓊斯勢相結(jié)合來描述。勢能函數(shù)的可靠性主要取決于力場參數(shù)準確性，而力場參數(shù)通常通過擬合實驗觀測數(shù)據(jù)和量子力學(xué)從頭計算得到的數(shù)據(jù)。目前在生物大分子體系模擬中使用最為廣泛的分子力場是CHARMM力場[3]和AMBER力場[4]，也是早期研究生物大分子的分子力場，其現(xiàn)有的力場參數(shù)仍在不斷優(yōu)化，并且涵蓋的分子類型也在擴大。粗?；Ｐ驮谟嬎闵镂锢硌芯恐性絹碓揭鹑藗兊年P(guān)注[5,6]由于該模型中定義了粗粒化粒子，對應(yīng)于全原子模型中的若干原子或原子基團甚至分子，減少了體系中的粒子數(shù)和自由度，使得模擬的時間和空間尺度得以大幅度提高，雖然會丟失一些原子細節(jié)信息，但是這種模型是應(yīng)用于研究緩慢的生物現(xiàn)象或依賴于大組裝體的生物現(xiàn)象[7]，如生物膜的波動，對它的模擬需要巨大的膜片。1.2常見力場分子動力計算體系由最初的單原子分子系統(tǒng)延伸至多原子分子、聚合物分子、生化分子系統(tǒng)，力場也隨著系統(tǒng)復(fù)雜度的增加而增加其復(fù)雜性。針對特定的目的，力場分為許多不同的形式，具有不同的使用范圍與局限性，在執(zhí)行分子動力計算時，選擇適合的力場極為重要，往往決定了計算成果的優(yōu)劣與可靠性。下面就一些常用力場做簡單介紹。傳統(tǒng)力場：AMBER力場[4]：由Kollman課題組開發(fā)的力場，是目前使用比較廣泛的一種力場，適合處理較小的蛋白質(zhì)、核酸、多糖等生化分子，此力場參數(shù)全部來自計算結(jié)果與實驗值的比對。CHARMM力場[3]：由Karplus課題組開發(fā)，CHARM力場參數(shù)除了來自計算結(jié)果與實驗值的比對外，還引用了大量的量子計算結(jié)果為依據(jù)。此力場可應(yīng)用于研究許多分子系統(tǒng)，包括小的有機分子，溶液，聚合物，生化分子等。幾乎除了有機金屬分子外，通常皆可得到與實驗值相近的結(jié)果。MMX力場[8,9]：此力場為Allinger等人所發(fā)展，依其發(fā)展的先后順序分別為MM2，MM3，MM4，MM+等。MM力場將一些常見的原子細分，如將碳原子分為sp3，sp2，sp，酮基碳，環(huán)丙烷碳，碳自由基，碳陽離子等。這些不同形態(tài)的碳原子具有不同的力場參數(shù)。此力場適用于各種有機化合物，自由基，離子。在MM形式的力場中仔細考慮了許多交叉作用項，其結(jié)果往往優(yōu)于其他形式的力場。相對的，其力場形式較為復(fù)雜，比較不易程序化，計算也較費時。CVFF力場：為DauberOsguthope等所發(fā)展。此力場最初以生化分子為主，其力場參數(shù)適用于胺基酸，水，及各種官能基。其后，經(jīng)過不斷強化，CVFF力場可適用于各種多肽，蛋白質(zhì)和大量有機分子。第二代力場：第二代的勢能函數(shù)形式比傳統(tǒng)力場要更加復(fù)雜，涉及的力場參數(shù)更多，計算量也更大，當(dāng)然也相應(yīng)地更加準確。CFF力場：CFF力場是一個力場家族，包括了CFF91、PCFF、CFF95等很多力場，可以進行從有機小分子、生物大分子到分子篩等諸多體系的計算。COMPASS力場[10]：由MSI公司開發(fā)的力場，擅長進行高分子體系的計算。MMF94力場：Hagler開發(fā)的力場，是目前最準確的力場之一。SFF力場：SI公司開發(fā)的力場，可以進行有機、無機分子的計算。UFF力場：可以計算周期表上所有元素的參數(shù)。Dreiding力場[11]：用于有機小分子、大分子、主族元素的計算。2分子動力學(xué)關(guān)鍵算法分子動力學(xué)發(fā)展到現(xiàn)在成為一種穩(wěn)健的理論工具，并因此成為許多實驗方法必不可少的補充，離不開前人在熱力學(xué)平衡模擬和非平衡模擬方面所做的大量工作。有兩個算法非常值得我們回顧，那就是Nosé-Hoover恒溫算法[12,13]和Nosé-Andersen恒壓算法[14]。這兩個算法分別實現(xiàn)了正則體系的模擬以及恒壓恒焓體系的模擬。二者的結(jié)合，實現(xiàn)了目前常用的恒溫恒壓體系的模擬，而恒溫恒壓體系最接近大多數(shù)實驗條件。在這些算法中，格點加和算法[15]旨在處理長程靜電相互作用。而與其相關(guān)的，是一種使用不同時間步長[16,17]模擬運動演化的數(shù)值方法，即利用Trotter因子分解，將力場中短程力和長程力的貢獻分開處理，對于運動較快的和較慢的自由度也采用分別處理的方法。這樣一來，對于運動方程進行積分時就可以采用不同的時間步長。與多時間步長方法緊密相關(guān)的另一種方法是利用完整約束[18]，去除體系中一些較“硬”（具有較高振動頻率）的自由度，從而可增加時間步長。例如，可以通過Shake和Rattle算法[18]將含有氫原子的化學(xué)鍵保持在其平衡值不變，以消除其振動自由度，這樣就可以使用較大的積分步長，而又不致影響到系統(tǒng)總能量的守恒。這一方法的根位數(shù)是用來描述中心原子第一殼層內(nèi)原子的平均數(shù)目，反映的是中心原子與其它原子的結(jié)合能力和配位關(guān)系，描述的是體系中粒子排列的緊密程度，配位數(shù)越大，粒子排列越緊密。一般來說配位數(shù)是通過對徑向分布函數(shù)的第一峰進行積分來得到。對于Voronoi指數(shù)分析而言，由于泰森多邊形的算法就是計算中心原子與配位原子之間的連線平分面，因此，通過對Voronoi指數(shù)的累加，也能得到中心原子的配位數(shù)。對于晶體結(jié)構(gòu)，通過配位數(shù)可以判斷出晶體的結(jié)構(gòu)；對于液態(tài)的非晶態(tài)結(jié)構(gòu)，配位數(shù)可以作為一個發(fā)生結(jié)構(gòu)轉(zhuǎn)變的敏感參量，為結(jié)構(gòu)轉(zhuǎn)變的判斷提供依據(jù)。例如，對于液體中發(fā)生的液?液相變現(xiàn)象，眾多研究[39]都將配位數(shù)的變化作為判斷的依據(jù)之一。4展望如果分子體系足夠小，我們又有足夠的耐心，即便只使用為數(shù)不多的處理器，完成微秒甚至是毫秒級時間尺度的模擬已不成問題，幾年前就已經(jīng)實現(xiàn)。然而，只有大規(guī)模并行架構(gòu)，特別是千兆級運算速度的超級計算機的出現(xiàn),我們才得以處理數(shù)百萬個原子組成的分子體系。但現(xiàn)階段，毫秒級時間尺度或者百萬原子空間尺度的分子模擬仍然構(gòu)成技術(shù)上的巨大挑戰(zhàn)，這恰恰又昭示著一個新時代的來臨。隨著第一臺千兆級浮點運算速度的超級計算機投入使用，科學(xué)家已經(jīng)開始期待著億億級計算機的到來并開始考慮這種新型架構(gòu)將能夠處理的分子體系的種類。假使在未來的幾十年，分子動力學(xué)的發(fā)展仍然緊跟摩爾定律所預(yù)測的計算機發(fā)展趨勢，vanGunsteren[40]樂觀地預(yù)測，未來20年內(nèi)，我們將可以在納秒級時間尺度上模擬一個完整的細菌，如大腸桿菌，大約20年后，可以模擬一個完整的哺乳動物細胞。這些基于目前硬件性能和軟件發(fā)展水平所作出的推測固然很振奮人心，但我們應(yīng)該清楚地認識到橫亙在模擬時間尺度和空間尺度之間的鴻溝依然難以逾越。在未來的計算機模擬中，利用直接無偏分子動力學(xué)模擬研究生物大體系，有可能提供時間相關(guān)和動力學(xué)相關(guān)信息。而實現(xiàn)這一目標(biāo)，我們?nèi)匀贿€有很長的路要走。參考文獻：[1]AlderBJ,WainwrightTE.Phasetransitionforahardspheresystem[J].TheJournalofchemicalphysics,1957,27(5):1208.[2]趙繼成.材料基因組計劃簡介[J].自然雜志,2014,36(2):89-104.[3]BrooksBR,BruccoleriRE,OlafsonBD,etal.CHARMM:aprogramformacromolecularenergy,minimization,anddynamicscalculations[J].Journalofcomputationalchemistry,1983,4(2):187-217.[4]WeinerSJ,KollmanPA,CaseDA,etal.Anewforcefieldformolecularmechanicalsimulationofnucleicacidsandproteins[J].JournaloftheAmericanChemicalSociety,1984,106(3):765-784.[5]IzvekovS,VothGA.Amultiscalecoarse-grainingmethodforbiomolecularsystems[J].TheJournalofPhysicalChemistryB,2005,109(7):2469-2473.[6]MarrinkSJ,RisseladaHJ,YefimovS,etal.TheMARTINIforcefield:coarsegrainedmodelforbiomolecularsimulations[J].TheJournalofPhysicalChemistryB,2007,111(27):7812-7824.[7]KleinML,ShinodaW.Large-scalemoleculardynamicssimulationsofself-assemblingsystems[J].Science,2008,321(5890):798-800.[8]AllingerNL,YuhYH,LiiJH.JAmChemSoc111:8551;(b)LiiJH[J].AllingerNL(1989)JAmChemSoc,1989,111:8566.[9]AllingerNL,ChenK,LiiJH.Animprovedforcefield(MM4)forsaturatedhydrocarbons[J].Journalofcomputationalchemistry,1996,17(5‐6):642-668.[10]SunH.COMPASS:anabinitioforce-fieldoptimizedforcondensed-phaseapplicationsoverviewwithdetailsonalkaneandbenzenecompounds[J].TheJournalofPhysicalChemistryB,1998,102(38):7338-7364.[11]RappéAK,CasewitCJ,ColwellKS,etal.UFF,afullperiodictableforcefieldformolecularmechanicsandmoleculardynamicssimulations[J].JournaloftheAmericanchemicalsociety,1992,114(25):10024-10035.[12]NoséS.Amoleculardynamicsmethodforsimulationsinthecanonicalensemble[J].Molecularphysics,1984,52(2):255-268.[13]HooverWG.Canonicaldynamics:equilibriumphase-spacedistributions[J].PhysicalreviewA,1985,31(3):1695.[14]AndersenHC.Moleculardynamicssimulationsatconstantpressureand/ortemperature[J].TheJournalofchemicalphysics,1980,72(4):2384-2393.[15]DardenT,YorkD,PedersenL.ParticlemeshEwald:AnN?log(N)methodforEwaldsumsinlargesystems[J].TheJournalofchemicalphysics,1993,98(12):10089-10092.[16]TuckermanM,BerneBJ,MartynaGJ.Reversiblemultipletimescalemoleculardynamics[J].TheJournalofchemicalphysics,1992,97(3):1990-2001.[17]IzaguirreJA,ReichS,SkeelRD.Longertimestepsformoleculardynamics[J].TheJournalofchemicalphysics,1999,110(20):9853-9864.[18]MiyamotoS,KollmanPA.SETTLE:ananalyticalversionoftheSHAKEandRATTLEalgorithmforrigidwatermodels[J].Journalofcomputationalchemistry,1992,13(8):952-962.[19]BrownD,ClarkeJHR,OkudaM,etal.Adomaindecompositionparallelizationstrategyformoleculardynamicssimulationsondistributedmemorymachines[J].ComputerPhysicsCommunications,1993,74(1):67-80.[20]ChipotC,ángyánJG.Continuingchallengesintheparametrizationofintermolecularforcefields.Towardsanaccuratedescriptionofelectrostaticandinductionterms[J].Newjournalofchemistry,2005,29(3):411-420.[21]SaguiC,PedersenLG,DardenTA.Towardsanaccuraterepresentationofelectrostaticsinclassicalforcefields:Efficientimplementationofmultipolarinteractionsinbiomolecularsimulations[J].TheJournalofchemicalphysics,2004,120(1):73-87.[22]Nanoscience:Nanobiotechnologyandnanobiology[M].SpringerScience&BusinessMedia,2009.[23]PlimptonS.Fastparallelalgorithmsforshort-rangemoleculardynamics[J].Journalofcomputationalphysics,1995,117(1):1-19.[24]HessB,KutznerC,VanDerSpoelD,etal.GROMACS4:algorithmsforhighlyefficient,load-balanced,andscalablemolecularsimulation[J].Journalofchemicaltheoryandcomputation,2008,4(3):435-447.[25]StoneJE,PhillipsJC,FreddolinoPL,etal.Acceleratingmolecularmodelingapplicationswithgraphicsprocessors[J].Journalofcomputationalchemistry,2007,28(16):2618-2640.[26]AndersonJA,LorenzCD,TravessetA.Generalpurposemoleculardynamicssimulationsfullyimplementedongraphicsprocessingunits[J].JournalofComputationalPhysics,2008,227(10):5342-5359.[27]MarchNH,TosiMP,MarchNH.Introductiontoliquidstatephysics[M].Singapore:WorldScientific,2002.[28]AbrahamFF.Anisothermal–isobariccomputersimulationofthesupercooled‐liquid/glasstransitionregion:Istheshort‐rangeorderintheamorphoussolidfcc?[J].Thejournalofchemicalphysics,1980,72(1):359-365.[29]ChengYQ,MaE.Atomic-levelstructureandstructure–propertyrelationshipinmetallicglasses[J].ProgressinMaterialsScience,2011,56(4):379-473.[30]QiL,DongLF,ZhangSL,etal.GlassformationandlocalstructureevolutioninrapidlycooledPdNialloymeltunderhighpressure[J].PhysicsLettersA,2008,372(5):708-711.[31]KimTH,KeltonKF.Structuralstudyofsupercooledliquidtransitionmetals[J].TheJournalofchemicalphysics,2007,126(5):054513.[32]LaioA,VandeVondeleJ,RothlisbergerU.AHamiltonianelectrostaticcouplingschemeforhybridCar–Parrinellomoleculardynamicssimulations[J].TheJournalofchemicalphysics,2002,116(16):6941-6947.[33]孫民華.液態(tài)物理概論[M].科學(xué)出版社,2013.[34]MarchNH,TosiMP,MarchNH.Introductiontoliquidstatephysics[M].Singapore:WorldScientific,2002.[35]Qi