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

基于分子動(dòng)力學(xué)的常用力場(chǎng)、算法及結(jié)果分析摘要：分子動(dòng)力學(xué)模擬是隨著計(jì)算機(jī)技術(shù)發(fā)展而興起的一項(xiàng)新技術(shù)，結(jié)合了物理，數(shù)學(xué)和化學(xué)的模擬方法。分子動(dòng)力學(xué)模擬具有溝通宏觀特性與微觀結(jié)構(gòu)的作用，對(duì)于許多在理論分析和實(shí)驗(yàn)都難以了解和觀察到的現(xiàn)象能給出一定的微觀解釋，可以分析其中的物理現(xiàn)象和本質(zhì)。本文主要對(duì)分子動(dòng)力學(xué)中常用力場(chǎng)，算法及結(jié)果分析方法進(jìn)行綜述。最后基于計(jì)算機(jī)硬件的發(fā)展與模擬算法的優(yōu)化，對(duì)該領(lǐng)域未來發(fā)展進(jìn)行展望。關(guān)鍵詞：分子動(dòng)力學(xué)；力場(chǎng)；算法；結(jié)果分析；引言早在20世紀(jì)50年代末，人們就已經(jīng)開始進(jìn)行計(jì)算模擬實(shí)驗(yàn)。一開始，這些計(jì)算實(shí)驗(yàn)的對(duì)象只是一些小的模型體系，比如幾個(gè)硬球的集合。但早在那時(shí)，Alder和Wainwright[1]就已經(jīng)預(yù)見到,這樣的計(jì)算機(jī)模擬實(shí)驗(yàn)將會(huì)成為聯(lián)結(jié)宏觀實(shí)驗(yàn)現(xiàn)象和其微觀本質(zhì)的重要橋梁。盡管解析理論進(jìn)行模擬也可以部分解釋實(shí)驗(yàn)，但是往往附帶了很多近似，而計(jì)算機(jī)模擬能夠提供體系的詳盡信息，具有更大的優(yōu)勢(shì)。由于宏觀和微觀在尺度上是截然不同的，為了能夠在兩者之間建立切實(shí)的關(guān)聯(lián)，即微觀模擬體系能夠反映相應(yīng)的宏觀實(shí)驗(yàn)現(xiàn)象，需要通過周期性邊界條件對(duì)模擬對(duì)象體系進(jìn)行周期性的復(fù)制，以避免在實(shí)際中并不存在的邊緣效應(yīng)。分子動(dòng)力學(xué)模擬主要是依靠牛頓力學(xué)來模擬分子體系的運(yùn)動(dòng)，以在由分子體系的不同狀態(tài)構(gòu)成的系統(tǒng)中抽取樣本，從而計(jì)算體系的構(gòu)型積分，并以構(gòu)型積分的結(jié)果為基礎(chǔ)進(jìn)一步計(jì)算體系的熱力學(xué)量和其他宏觀性質(zhì)。通過分子模擬技術(shù)，人們可以在分子、原子尺度上對(duì)所關(guān)心材料的物理、化學(xué)性質(zhì)以及與之相對(duì)應(yīng)的各種現(xiàn)象進(jìn)行模擬，摒棄傳統(tǒng)材料科學(xué)研究中被稱之為“愛迪生方式”的炒菜式研究，利用計(jì)算機(jī)對(duì)材料的組成、性質(zhì)、制備工藝、加工工藝等環(huán)節(jié)進(jìn)行虛擬實(shí)驗(yàn)，獲得大量數(shù)據(jù)，再與少量驗(yàn)證性實(shí)驗(yàn)進(jìn)行對(duì)照，分析物理現(xiàn)象和本質(zhì)，從而大大的降低材料的研發(fā)成本和時(shí)間。這種新材料的設(shè)計(jì)模式已經(jīng)很大程度上取代了傳統(tǒng)的新材料設(shè)計(jì)方式，成為新材料設(shè)計(jì)的主要方式之一。2011年，美國(guó)啟動(dòng)一項(xiàng)價(jià)值超過5億美元的“先進(jìn)制造業(yè)伙伴關(guān)系”計(jì)劃，而以材料計(jì)算與模擬為主的“材料基因組計(jì)劃”則是其中重要內(nèi)容之一[2]。目前，對(duì)分子動(dòng)力學(xué)模擬的力場(chǎng)、算法及結(jié)果分析方法很多，但基本上是分散在英文文獻(xiàn)中，由于翻譯和理解上的差異，影響了國(guó)內(nèi)從事分子模擬科研工作者的理解交流。本文整理了目前分子動(dòng)力學(xué)中主流的力場(chǎng)、算法及結(jié)果分析，并在最后對(duì)該領(lǐng)域進(jìn)行了展望。1力場(chǎng)力場(chǎng)簡(jiǎn)介分子動(dòng)力學(xué)模擬是計(jì)算龐大復(fù)雜系統(tǒng)的有效方法，它以力場(chǎng)為依據(jù)，力場(chǎng)的完備與否決定計(jì)算的可靠程度。分子的總能量是動(dòng)能與勢(shì)能之和，分子的勢(shì)能通常表示為簡(jiǎn)單的幾何坐標(biāo)的函數(shù)。一般勢(shì)能中包括：（1）范德華力，與能量有關(guān)的非鍵相互作用交叉能量項(xiàng)，（2）構(gòu)成分子的各個(gè)化學(xué)鍵在鍵軸方向上的伸縮運(yùn)動(dòng)所引起的能量變化，（3）鍵角變化引起的分子能量變化，（4）單軸旋轉(zhuǎn)引起分子骨架扭曲所產(chǎn)生的能量變化，（5）離平面振動(dòng)項(xiàng)，共平面原子的中心原子離平面小幅振動(dòng)的勢(shì)能，（6）庫(kù)倫作用項(xiàng)，帶電荷粒子間存在的靜電吸引或排斥作用的勢(shì)能。力場(chǎng)可以看作是勢(shì)能面的經(jīng)驗(yàn)表達(dá)式，它是分子動(dòng)力學(xué)模擬的基礎(chǔ)、力場(chǎng)是通過原子位置計(jì)算體系能量的，與之前的量子力學(xué)方法相比，大大節(jié)約了計(jì)算時(shí)間，可用于計(jì)算包含上萬(wàn)粒子數(shù)目的體系。勢(shì)能函數(shù)在大多數(shù)情況下將描述分子幾何形變最大程度地簡(jiǎn)化為僅僅使用簡(jiǎn)諧項(xiàng)和三角函數(shù)來實(shí)現(xiàn)，而非鍵原子之間的相互作用，則只采用庫(kù)倫相互作用和蘭納-瓊斯勢(shì)相結(jié)合來描述。勢(shì)能函數(shù)的可靠性主要取決于力場(chǎng)參數(shù)準(zhǔn)確性，而力場(chǎng)參數(shù)通常通過擬合實(shí)驗(yàn)觀測(cè)數(shù)據(jù)和量子力學(xué)從頭計(jì)算得到的數(shù)據(jù)。目前在生物大分子體系模擬中使用最為廣泛的分子力場(chǎng)是CHARMM力場(chǎng)[3]和AMBER力場(chǎng)[4]，也是早期研究生物大分子的分子力場(chǎng)，其現(xiàn)有的力場(chǎng)參數(shù)仍在不斷優(yōu)化，并且涵蓋的分子類型也在擴(kuò)大。粗粒化模型在計(jì)算生物物理研究中越來越引起人們的關(guān)注[5,6]由于該模型中定義了粗粒化粒子，對(duì)應(yīng)于全原子模型中的若干原子或原子基團(tuán)甚至分子，減少了體系中的粒子數(shù)和自由度，使得模擬的時(shí)間和空間尺度得以大幅度提高，雖然會(huì)丟失一些原子細(xì)節(jié)信息，但是這種模型是應(yīng)用于研究緩慢的生物現(xiàn)象或依賴于大組裝體的生物現(xiàn)象[7]，如生物膜的波動(dòng)，對(duì)它的模擬需要巨大的膜片。1.2常見力場(chǎng)分子動(dòng)力計(jì)算體系由最初的單原子分子系統(tǒng)延伸至多原子分子、聚合物分子、生化分子系統(tǒng)，力場(chǎng)也隨著系統(tǒng)復(fù)雜度的增加而增加其復(fù)雜性。針對(duì)特定的目的，力場(chǎng)分為許多不同的形式，具有不同的使用范圍與局限性，在執(zhí)行分子動(dòng)力計(jì)算時(shí)，選擇適合的力場(chǎng)極為重要，往往決定了計(jì)算成果的優(yōu)劣與可靠性。下面就一些常用力場(chǎng)做簡(jiǎn)單介紹。傳統(tǒng)力場(chǎng)：AMBER力場(chǎng)[4]：由Kollman課題組開發(fā)的力場(chǎng)，是目前使用比較廣泛的一種力場(chǎng)，適合處理較小的蛋白質(zhì)、核酸、多糖等生化分子，此力場(chǎng)參數(shù)全部來自計(jì)算結(jié)果與實(shí)驗(yàn)值的比對(duì)。CHARMM力場(chǎng)[3]：由Karplus課題組開發(fā)，CHARM力場(chǎng)參數(shù)除了來自計(jì)算結(jié)果與實(shí)驗(yàn)值的比對(duì)外，還引用了大量的量子計(jì)算結(jié)果為依據(jù)。此力場(chǎng)可應(yīng)用于研究許多分子系統(tǒng)，包括小的有機(jī)分子，溶液，聚合物，生化分子等。幾乎除了有機(jī)金屬分子外，通常皆可得到與實(shí)驗(yàn)值相近的結(jié)果。MMX力場(chǎng)[8,9]：此力場(chǎng)為Allinger等人所發(fā)展，依其發(fā)展的先后順序分別為MM2，MM3，MM4，MM+等。MM力場(chǎng)將一些常見的原子細(xì)分，如將碳原子分為sp3，sp2，sp，酮基碳，環(huán)丙烷碳，碳自由基，碳陽(yáng)離子等。這些不同形態(tài)的碳原子具有不同的力場(chǎng)參數(shù)。此力場(chǎng)適用于各種有機(jī)化合物，自由基，離子。在MM形式的力場(chǎng)中仔細(xì)考慮了許多交叉作用項(xiàng)，其結(jié)果往往優(yōu)于其他形式的力場(chǎng)。相對(duì)的，其力場(chǎng)形式較為復(fù)雜，比較不易程序化，計(jì)算也較費(fèi)時(shí)。CVFF力場(chǎng)：為DauberOsguthope等所發(fā)展。此力場(chǎng)最初以生化分子為主，其力場(chǎng)參數(shù)適用于胺基酸，水，及各種官能基。其后，經(jīng)過不斷強(qiáng)化，CVFF力場(chǎng)可適用于各種多肽，蛋白質(zhì)和大量有機(jī)分子。第二代力場(chǎng)：第二代的勢(shì)能函數(shù)形式比傳統(tǒng)力場(chǎng)要更加復(fù)雜，涉及的力場(chǎng)參數(shù)更多，計(jì)算量也更大，當(dāng)然也相應(yīng)地更加準(zhǔn)確。CFF力場(chǎng)：CFF力場(chǎng)是一個(gè)力場(chǎng)家族，包括了CFF91、PCFF、CFF95等很多力場(chǎng)，可以進(jìn)行從有機(jī)小分子、生物大分子到分子篩等諸多體系的計(jì)算。COMPASS力場(chǎng)[10]：由MSI公司開發(fā)的力場(chǎng)，擅長(zhǎng)進(jìn)行高分子體系的計(jì)算。MMF94力場(chǎng)：Hagler開發(fā)的力場(chǎng)，是目前最準(zhǔn)確的力場(chǎng)之一。SFF力場(chǎng)：SI公司開發(fā)的力場(chǎng)，可以進(jìn)行有機(jī)、無(wú)機(jī)分子的計(jì)算。UFF力場(chǎng)：可以計(jì)算周期表上所有元素的參數(shù)。Dreiding力場(chǎng)[11]：用于有機(jī)小分子、大分子、主族元素的計(jì)算。2分子動(dòng)力學(xué)關(guān)鍵算法分子動(dòng)力學(xué)發(fā)展到現(xiàn)在成為一種穩(wěn)健的理論工具，并因此成為許多實(shí)驗(yàn)方法必不可少的補(bǔ)充，離不開前人在熱力學(xué)平衡模擬和非平衡模擬方面所做的大量工作。有兩個(gè)算法非常值得我們回顧，那就是Nosé-Hoover恒溫算法[12,13]和Nosé-Andersen恒壓算法[14]。這兩個(gè)算法分別實(shí)現(xiàn)了正則體系的模擬以及恒壓恒焓體系的模擬。二者的結(jié)合，實(shí)現(xiàn)了目前常用的恒溫恒壓體系的模擬，而恒溫恒壓體系最接近大多數(shù)實(shí)驗(yàn)條件。在這些算法中，格點(diǎn)加和算法[15]旨在處理長(zhǎng)程靜電相互作用。而與其相關(guān)的，是一種使用不同時(shí)間步長(zhǎng)[16,17]模擬運(yùn)動(dòng)演化的數(shù)值方法，即利用Trotter因子分解，將力場(chǎng)中短程力和長(zhǎng)程力的貢獻(xiàn)分開處理，對(duì)于運(yùn)動(dòng)較快的和較慢的自由度也采用分別處理的方法。這樣一來，對(duì)于運(yùn)動(dòng)方程進(jìn)行積分時(shí)就可以采用不同的時(shí)間步長(zhǎng)。與多時(shí)間步長(zhǎng)方法緊密相關(guān)的另一種方法是利用完整約束[18]，去除體系中一些較“硬”（具有較高振動(dòng)頻率）的自由度，從而可增加時(shí)間步長(zhǎng)。例如，可以通過Shake和Rattle算法[18]將含有氫原子的化學(xué)鍵保持在其平衡值不變，以消除其振動(dòng)自由度，這樣就可以使用較大的積分步長(zhǎng)，而又不致影響到系統(tǒng)總能量的守恒。這一方法的根位數(shù)是用來描述中心原子第一殼層內(nèi)原子的平均數(shù)目，反映的是中心原子與其它原子的結(jié)合能力和配位關(guān)系，描述的是體系中粒子排列的緊密程度，配位數(shù)越大，粒子排列越緊密。一般來說配位數(shù)是通過對(duì)徑向分布函數(shù)的第一峰進(jìn)行積分來得到。對(duì)于Voronoi指數(shù)分析而言，由于泰森多邊形的算法就是計(jì)算中心原子與配位原子之間的連線平分面，因此，通過對(duì)Voronoi指數(shù)的累加，也能得到中心原子的配位數(shù)。對(duì)于晶體結(jié)構(gòu)，通過配位數(shù)可以判斷出晶體的結(jié)構(gòu)；對(duì)于液態(tài)的非晶態(tài)結(jié)構(gòu)，配位數(shù)可以作為一個(gè)發(fā)生結(jié)構(gòu)轉(zhuǎn)變的敏感參量，為結(jié)構(gòu)轉(zhuǎn)變的判斷提供依據(jù)。例如，對(duì)于液體中發(fā)生的液?液相變現(xiàn)象，眾多研究[39]都將配位數(shù)的變化作為判斷的依據(jù)之一。4展望如果分子體系足夠小，我們又有足夠的耐心，即便只使用為數(shù)不多的處理器，完成微秒甚至是毫秒級(jí)時(shí)間尺度的模擬已不成問題，幾年前就已經(jīng)實(shí)現(xiàn)。然而，只有大規(guī)模并行架構(gòu)，特別是千兆級(jí)運(yùn)算速度的超級(jí)計(jì)算機(jī)的出現(xiàn),我們才得以處理數(shù)百萬(wàn)個(gè)原子組成的分子體系。但現(xiàn)階段，毫秒級(jí)時(shí)間尺度或者百萬(wàn)原子空間尺度的分子模擬仍然構(gòu)成技術(shù)上的巨大挑戰(zhàn)，這恰恰又昭示著一個(gè)新時(shí)代的來臨。隨著第一臺(tái)千兆級(jí)浮點(diǎn)運(yùn)算速度的超級(jí)計(jì)算機(jī)投入使用，科學(xué)家已經(jīng)開始期待著億億級(jí)計(jì)算機(jī)的到來并開始考慮這種新型架構(gòu)將能夠處理的分子體系的種類。假使在未來的幾十年，分子動(dòng)力學(xué)的發(fā)展仍然緊跟摩爾定律所預(yù)測(cè)的計(jì)算機(jī)發(fā)展趨勢(shì)，vanGunsteren[40]樂觀地預(yù)測(cè)，未來20年內(nèi)，我們將可以在納秒級(jí)時(shí)間尺度上模擬一個(gè)完整的細(xì)菌，如大腸桿菌，大約20年后，可以模擬一個(gè)完整的哺乳動(dòng)物細(xì)胞。這些基于目前硬件性能和軟件發(fā)展水平所作出的推測(cè)固然很振奮人心，但我們應(yīng)該清楚地認(rèn)識(shí)到橫亙?cè)谀M時(shí)間尺度和空間尺度之間的鴻溝依然難以逾越。在未來的計(jì)算機(jī)模擬中，利用直接無(wú)偏分子動(dòng)力學(xué)模擬研究生物大體系，有可能提供時(shí)間相關(guān)和動(dòng)力學(xué)相關(guān)信息。而實(shí)現(xiàn)這一目標(biāo)，我們?nèi)匀贿€有很長(zhǎng)的路要走。參考文獻(xiàn)：[1]AlderBJ,WainwrightTE.Phasetransitionforahardspheresystem[J].TheJournalofchemicalphysics,1957,27(5):1208.[2]趙繼成.材料基因組計(jì)劃簡(jiǎn)介[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