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

1、基于分子動(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í)，Alde

2、r和Wainwright1就已經(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é)量和

3、其他宏觀性質(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)

4、、算法及結(jié)果分析方法很多，但基本上是分散在英文文獻(xiàn)中，由于翻譯和理解上的差異，影響了國(guó)內(nèi)從事分子模擬科研工作者的理解交流。本文整理了目前分子動(dòng)力學(xué)中主流的力場(chǎng)、算法及結(jié)果分析，并在最后對(duì)該領(lǐng)域進(jìn)行了展望。1 力場(chǎng)1.1 力場(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、能量變化，（5）離平面振動(dòng)項(xiàng)，共平面原子的中心原子離平面小幅振動(dòng)的勢(shì)能，（6）庫倫作用項(xiàng)，帶電荷粒子間存在的靜電吸引或排斥作用的勢(shì)能。力場(chǎng)可以看作是勢(shì)能面的經(jīng)驗(yàn)表達(dá)式，它是分子動(dòng)力學(xué)模擬的基礎(chǔ)、力場(chǎng)是通過原子位置計(jì)算體系能量的，與之前的量子力學(xué)方法相比，大大節(jié)約了計(jì)算時(shí)間，可用于計(jì)算包含上萬粒子數(shù)目的體系。勢(shì)能函數(shù)在大多數(shù)情況下將描述分子幾何形變最大程度地簡(jiǎn)化為僅僅使用簡(jiǎn)諧項(xiàng)和三角函數(shù)來實(shí)現(xiàn)，而非鍵原子之間的相互作用，則只采用庫倫相互作用和蘭納-瓊斯勢(shì)相結(jié)合來描述。勢(shì)能函數(shù)的可靠性主要取決于力場(chǎng)參數(shù)準(zhǔn)確性，而力場(chǎng)參數(shù)通常通過擬合實(shí)驗(yàn)觀測(cè)數(shù)據(jù)和量子力學(xué)從頭計(jì)算得到的數(shù)據(jù)。目前在生物大分子體系模擬

6、中使用最為廣泛的分子力場(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)

7、也隨著系統(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ī)分子，溶液，聚合物，生化分子

8、等。幾乎除了有機(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)丙烷碳，碳自由基，碳陽離子等。這些不同形態(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)：為Dauber Osguthope等所發(fā)展。此力場(chǎng)最初以生化分子為主，其力場(chǎng)參數(shù)適用于胺基

9、酸，水，及各種官能基。其后，經(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ī)、無機(jī)分子的計(jì)算。UFF力場(chǎng)：可以計(jì)算周期表上所有元素的參數(shù)。Dre

10、iding力場(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ù)值

11、方法，即利用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)總能量的守恒。這一方法的根本目的是在不增加計(jì)算消耗的情況下能夠獲得更多的采樣。但是，由于增加了額外的計(jì)算求解約束方程，總體上所提高的計(jì)算

12、效率并不理想。但是，在人們孜孜以求不斷增大時(shí)間和空間尺度的過程中，真正革命性的成就莫過于并行架構(gòu)和空間分解算法19的出現(xiàn)，使得計(jì)算時(shí)間可以隨著處理器數(shù)目的增加而線性地減少。 然而，這一偉大的成就并沒有讓人們追求最長(zhǎng)模擬時(shí)間或最大模擬體系的欲望得到徹底的滿足，反而變得更加強(qiáng)烈，同時(shí)，在一定程度上也遮住了在改進(jìn)原子間相互作用力的表達(dá)方面所取得的輝煌成就的光芒，特別是極化率20和分布式多極矩21的引入，當(dāng)然這也相應(yīng)地增加了計(jì)算的開銷，但這是為了提高計(jì)算精度而必須付出的代價(jià)，另一方面也多少忽視了在表征運(yùn)動(dòng)方程數(shù)值積分的誤差方面所付出的巨大努力。分子動(dòng)力學(xué)模擬要求按一定的時(shí)間步長(zhǎng)對(duì)經(jīng)典運(yùn)動(dòng)方程進(jìn)行不連續(xù)

13、的積分，而時(shí)間步長(zhǎng)受制于分子體系中運(yùn)動(dòng)最快的自由度22，在對(duì)生物體系的全原子描述中，由于要描述包括氫鍵在內(nèi)的化學(xué)鍵的振動(dòng)，為了確保體系能量守恒，時(shí)間步長(zhǎng)只能在飛秒數(shù)量級(jí)。因此，需要進(jìn)行數(shù)百萬甚至數(shù)十億次積分才能在與生物體系相應(yīng)的時(shí)間尺度上描述體系的行為。如前所述，盡可能減少計(jì)算能量時(shí)的開銷是非常重要的。另一方面，在對(duì)勢(shì)加和的方法中，理論上的計(jì)算開銷正比于N2，即體系中粒子數(shù)的平方。在實(shí)際計(jì)算中15，對(duì)靜電相互作用中的短程力部分可采用球形截?cái)嗖呗?，而求解Poisson方程時(shí)則可采用離散化的處理方案，這樣，就可以在倒易空間里處理長(zhǎng)程相互作用，從而使實(shí)際的計(jì)算開銷減少到N logN。在利用大型中央處

14、理器CPU陣列運(yùn)行大規(guī)模并行分子動(dòng)力學(xué)程序的做法出現(xiàn)之后大約二十年，數(shù)值模擬領(lǐng)域又成為了另一場(chǎng)革新的舞臺(tái)，這一革新是由通用圖形處理器GPU的出現(xiàn)而引發(fā)的。如今，通用圖形處理器已經(jīng)變成一種并不昂貴的多核通用處理器，它能夠以并行的方式處理浮點(diǎn)數(shù)據(jù)運(yùn)算，而制造成本的降低得益于其在消費(fèi)類電子產(chǎn)品，尤其是在視頻游戲控制平臺(tái)方面的廣泛應(yīng)用。曾經(jīng)風(fēng)靡一時(shí)的程序，如NAMD23和GROMACS24等都很快作出了調(diào)整以適應(yīng)這種新型的異構(gòu)化并行架構(gòu)25，同時(shí)，新的專門針對(duì)圖形處理器而設(shè)計(jì)的分子動(dòng)力學(xué)代碼如HOOMD-blue26等也應(yīng)運(yùn)而生。3 結(jié)果分析方法3.1 常用宏觀統(tǒng)計(jì)分析方法對(duì)分布函數(shù)、徑向分布函數(shù)和靜

15、態(tài)結(jié)構(gòu)因子都是表征體系結(jié)構(gòu)的宏觀統(tǒng)計(jì)參量，三者之間可以進(jìn)行數(shù)值的相互轉(zhuǎn)化。由于靜態(tài)結(jié)構(gòu)因子也可以通過X衍射或中子衍射在實(shí)驗(yàn)中得到，因此對(duì)分布函數(shù)、徑向分布函數(shù)和靜態(tài)結(jié)構(gòu)因子這3種方法就可以作為連接計(jì)算機(jī)模擬與實(shí)驗(yàn)分析的橋梁，成為分子動(dòng)力學(xué)模擬中最常使用的結(jié)構(gòu)分析方法。對(duì)分布函數(shù)（PDF）反映的是從一個(gè)任意指定的“中心”粒子出發(fā)，到半徑為的位置上出現(xiàn)其它粒子的幾率（單位體積內(nèi)的粒子數(shù)目）27。如圖1所示，取一原子為觀測(cè)中心，以1 為單位向外畫出一層層的同心球殼，然后統(tǒng)計(jì)每一殼層里的原子數(shù)密度（粒子數(shù)n(r)/體積V）與平均數(shù)密度（0總粒子數(shù)/體積V）比值，即為對(duì)分布函數(shù)，常用g(r)表示。圖1

16、對(duì)分布函數(shù)與原子結(jié)構(gòu)關(guān)系示意圖對(duì)于g(r)的分析，主要是基于以下幾點(diǎn)：峰的分布規(guī)律。根據(jù)對(duì)分布函數(shù)的物理含義，對(duì)于晶體，由于原子的周期性排列，每層原子會(huì)對(duì)應(yīng)出現(xiàn)分布峰值，而在層與層的原子之間出現(xiàn)原子的幾率為0。對(duì)于液體結(jié)構(gòu)，由于短程有序，在第一近鄰處會(huì)出現(xiàn)較為尖銳的峰值，在中遠(yuǎn)程的范圍里，液體結(jié)構(gòu)的對(duì)分布函數(shù)會(huì)出現(xiàn)漫高峰，并且趨近于1。這說明在液態(tài)結(jié)構(gòu)中，相對(duì)于一個(gè)中心原子，在無窮遠(yuǎn)處會(huì)總會(huì)有一個(gè)原子的存在。第一峰的高度和形狀。第一峰代表第一近鄰原子之間的結(jié)合強(qiáng)度，如果它的外形比較尖銳，表明在此半徑范圍內(nèi)的原子數(shù)密度比平均密度要高很多，中心原子與最近鄰原子的相互結(jié)合強(qiáng)度也比較大。第一峰與第一谷

17、的關(guān)系。Abraham28等人提出一種經(jīng)驗(yàn)性方法，根據(jù)對(duì)分布函數(shù)中第一谷值gmin與第一峰值gmax的比例，定義參 量R=gmin/gmax，用來確定玻璃轉(zhuǎn)變點(diǎn)Tg，如圖2(a)所示。第二峰的形狀與分布。第二峰代表著中程序的連接強(qiáng)度，如果第一峰和第二峰比較明顯，同時(shí)第一峰與第二峰之間有一個(gè)很深的波谷，基本上表明原子是以一種共價(jià)鍵的形式存在。如果第二峰出現(xiàn)了劈裂，則這是出現(xiàn)非晶結(jié)構(gòu)的典型特征29, 30。圖2 PDF參量、原子體積與溫度關(guān)系徑向分布函數(shù)（RDF）表示在半徑為r處厚度r的球殼內(nèi)的平均原子數(shù)，其數(shù)值等于該位置出現(xiàn)粒子的幾率g(r)與球殼面積之積。常用符號(hào)G(r)表示。由于語言轉(zhuǎn)換和理

18、解上的差異，徑向分布函數(shù)與對(duì)分布函數(shù)在中外文的文獻(xiàn)中經(jīng)常被互用，有時(shí)使用的是徑向分布函數(shù)G(r)概念，表達(dá)出來的卻是對(duì)分布函數(shù)g(r)含義。其實(shí)二者無論是在圖像上還是物理意義上都有明確的區(qū)別31, 32。如圖所示，對(duì)分布函數(shù)g(r)表示的是概率，對(duì)于液體或非晶結(jié)構(gòu)，分布函數(shù)曲線隨著半徑r的增加應(yīng)該是歸一的，如圖3(a)；而徑向分布函數(shù)G(r)則表示在球面上粒子的統(tǒng)計(jì)值，曲線隨著半徑的增加整體上升，如圖3(b)。圖3 對(duì)分布函數(shù)與徑向分布函數(shù)的關(guān)系靜態(tài)結(jié)構(gòu)因子（SSF）是一個(gè)連接實(shí)驗(yàn)分析與模擬分析的重要參量，對(duì)于衍射分析來說，它表征的是材料對(duì)射線的散射能力，反映結(jié)構(gòu)的平均信息。靜態(tài)結(jié)構(gòu)因子(S(

19、k)-1)與對(duì)分布函數(shù)(g(r)-1)互為傅里葉變換，在粒子的空間結(jié)構(gòu)信息中，靜態(tài)結(jié)構(gòu)因子為倒易空間，對(duì)分布函數(shù)為實(shí)空間。峰位置的相對(duì)關(guān)系。由于結(jié)構(gòu)因子是倒易空間，變量k與距中心原子的半徑r成反比關(guān)系。根據(jù)第二峰位與第一峰位k2/k1的比值，可以將液體分為3種類型。第一種為大部分的液態(tài)金屬，比值為1.86；第二種為Ga, Sn, Bi等液體，比值為1.96，其特征為在第一峰右邊存在一個(gè)肩峰；第三種為In和Tl，比值為1.8833；預(yù)峰的特征。在第一峰的前面，如果出現(xiàn)小的預(yù)峰，反映的是體系中存在著中程序結(jié)構(gòu)特征33；熔點(diǎn)溫度峰值特征。通過對(duì)多種液體的統(tǒng)計(jì)發(fā)現(xiàn)，液態(tài)結(jié)構(gòu)因子的第一峰高度在熔點(diǎn)附近都

20、會(huì)達(dá)到一個(gè)同樣的高度，即S(k) 2.8，這個(gè)規(guī)律叫做Hansen-Verlet凝固判據(jù)，是可以用來作為熔化和凝固的判據(jù)之一34。3.2 常用共近鄰分析方法共近鄰分析（CNA）是用來描述一個(gè)原子周圍近鄰原子的排列特征的方法，根據(jù)計(jì)算方法的不同，可包括鍵對(duì)分析、鍵角分析、泰森多邊形指數(shù)分析、配位數(shù)等。最近鄰的距離一般是采用對(duì)分布函數(shù)的第一低谷處的原子位置值rmin。鍵對(duì)分析是研究液體和固體中團(tuán)簇的一種常用的有效方法35。這種方法用4個(gè)指數(shù)i, j, l, m對(duì)原子進(jìn)行分類標(biāo)定。i=1表示研究的兩個(gè)原子對(duì)能夠成鍵，i=2表示兩個(gè)原子對(duì)不能成鍵，一般我們都關(guān)注于成鍵原子之間的關(guān)系。j表示兩個(gè)被研究原

21、子對(duì)中相同的近鄰原子數(shù)。1則表示j個(gè)相同的近鄰原子中成鍵的原子對(duì)數(shù)。m表示l中成鍵近鄰原子的排布情況，如果成鍵相連，則m為1，否則m為236。鍵角分布就是體系中所有原子與其它最近鄰原子成鍵所構(gòu)成的鍵角的分布，所獲取的是原子在三維空間上的排布信息，因而可以彌補(bǔ)對(duì)分布函數(shù)g(r)僅僅依據(jù)原子距離結(jié)構(gòu)獲取信息的不足。然而由于溫度的作用，原子之間的振動(dòng)會(huì)影響實(shí)際的角度分布情況，因此從角度的分布變化上，可以看出體系中的二十面體或四面體的扭曲和破壞情況37。泰森多邊形是一種幾何算法，在氣象、物流、商業(yè)等多個(gè)行業(yè)和領(lǐng)域中已被廣泛利用，也常被稱為Voronoi網(wǎng)格。它是由兩個(gè)鄰點(diǎn)連接直線的垂直平分線組成的連續(xù)

22、多邊形組成，二維泰森多邊形原理示意圖如圖所示。這種方法可以描述在平面上的N個(gè)點(diǎn)，按照最鄰近原則劃分與它的最近鄰區(qū)域的相互關(guān)系。利用Voronoi算法，將離散的原子進(jìn)行三維空間上的劃分，利用所得到的空間結(jié)構(gòu)來表征材料的微觀結(jié)構(gòu)，目前也已成為結(jié)構(gòu)表征的一種常用手段。泰森多邊形指數(shù)常用的形式表示，其中ni代表著Voronoi 多面體具有的i邊形數(shù)38。第一近鄰的配位數(shù)也是一種傳統(tǒng)的分析微觀結(jié)構(gòu)的方式，此概念首先由阿爾弗雷德維爾納提出，原來更多地用于化學(xué)結(jié)構(gòu)的表征，是配位化學(xué)的基礎(chǔ)。在微觀結(jié)構(gòu)的表征中，配位數(shù)是用來描述中心原子第一殼層內(nèi)原子的平均數(shù)目，反映的是中心原子與其它原子的結(jié)合能力和配位關(guān)系，描

23、述的是體系中粒子排列的緊密程度，配位數(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í)間尺

24、度的模擬已不成問題，幾年前就已經(jīng)實(shí)現(xiàn)。然而，只有大規(guī)模并行架構(gòu)，特別是千兆級(jí)運(yùn)算速度的超級(jí)計(jì)算機(jī)的出現(xiàn), 我們才得以處理數(shù)百萬個(gè)原子組成的分子體系。但現(xiàn)階段，毫秒級(jí)時(shí)間尺度或者百萬原子空間尺度的分子模擬仍然構(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ì)，van Gunsteren40樂觀地預(yù)測(cè)，未來20年內(nèi)，我們將可以在納秒級(jí)時(shí)間尺度上模擬一個(gè)完整的細(xì)菌，如大腸桿菌，大約20年

25、后，可以模擬一個(gè)完整的哺乳動(dòng)物細(xì)胞。這些基于目前硬件性能和軟件發(fā)展水平所作出的推測(cè)固然很振奮人心，但我們應(yīng)該清楚地認(rèn)識(shí)到橫亙?cè)谀M時(shí)間尺度和空間尺度之間的鴻溝依然難以逾越。在未來的計(jì)算機(jī)模擬中，利用直接無偏分子動(dòng)力學(xué)模擬研究生物大體系，有可能提供時(shí)間相關(guān)和動(dòng)力學(xué)相關(guān)信息。而實(shí)現(xiàn)這一目標(biāo)，我們?nèi)匀贿€有很長(zhǎng)的路要走。參考文獻(xiàn)：1 Alder B J, Wainwright T E. Phase transition for a hard sphere systemJ. The Journal of chemical physics, 1957, 27(5): 1208.2 趙繼成. 材料基因組計(jì)劃

26、簡(jiǎn)介J. 自然雜志, 2014, 36(2): 89-104.3 Brooks B R, Bruccoleri R E, Olafson B D, et al. CHARMM: a program for macromolecular energy, minimization, and dynamics calculationsJ. Journal of computational chemistry, 1983, 4(2): 187-217.4 Weiner S J, Kollman P A, Case D A, et al. A new force field for molecular

27、mechanical simulation of nucleic acids and proteinsJ. Journal of the American Chemical Society, 1984, 106(3): 765-784.5 Izvekov S, Voth G A. A multiscale coarse-graining method for biomolecular systemsJ. The Journal of Physical Chemistry B, 2005, 109(7): 2469-2473.6 Marrink S J, Risselada H J, Yefim

28、ov S, et al. The MARTINI force field: coarse grained model for biomolecular simulationsJ. The Journal of Physical Chemistry B, 2007, 111(27): 7812-7824.7 Klein M L, Shinoda W. Large-scale molecular dynamics simulations of self-assembling systemsJ. Science, 2008, 321(5890): 798-800.8 Allinger N L, Yu

29、h Y H, Lii J H. J Am Chem Soc 111: 8551;(b) Lii JHJ. Allinger NL (1989) J Am Chem Soc, 1989, 111: 8566.9 Allinger N L, Chen K, Lii J H. An improved force field (MM4) for saturated hydrocarbonsJ. Journal of computational chemistry, 1996, 17(56): 642-668.10 Sun H. COMPASS: an ab initio force-field opt

30、imized for condensed-phase applications overview with details on alkane and benzene compoundsJ. The Journal of Physical Chemistry B, 1998, 102(38): 7338-7364.11 Rapp A K, Casewit C J, Colwell K S, et al. UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulation

31、sJ. Journal of the American chemical society, 1992, 114(25): 10024-10035.12 Nos S. A molecular dynamics method for simulations in the canonical ensembleJ. Molecular physics, 1984, 52(2): 255-268.13 Hoover W G. Canonical dynamics: equilibrium phase-space distributionsJ. Physical review A, 1985, 31(3)

32、: 1695.14 Andersen H C. Molecular dynamics simulations at constant pressure and/or temperatureJ. The Journal of chemical physics, 1980, 72(4): 2384-2393.15 Darden T, York D, Pedersen L. Particle mesh Ewald: An N log (N) method for Ewald sums in large systemsJ. The Journal of chemical physics, 1993,

33、98(12): 10089-10092.16 Tuckerman M, Berne B J, Martyna G J. Reversible multiple time scale molecular dynamicsJ. The Journal of chemical physics, 1992, 97(3): 1990-2001.17 Izaguirre J A, Reich S, Skeel R D. Longer time steps for molecular dynamicsJ. The Journal of chemical physics, 1999, 110(20): 985

34、3-9864.18 Miyamoto S, Kollman P A. SETTLE: an analytical version of the SHAKE and RATTLE algorithm for rigid water modelsJ. Journal of computational chemistry, 1992, 13(8): 952-962.19 Brown D, Clarke J H R, Okuda M, et al. A domain decomposition parallelization strategy for molecular dynamics simula

35、tions on distributed memory machinesJ. Computer Physics Communications, 1993, 74(1): 67-80.20 Chipot C, ngyn J G. Continuing challenges in the parametrization of intermolecular force fields. Towards an accurate description of electrostatic and induction termsJ. New journal of chemistry, 2005, 29(3):

36、 411-420.21 Sagui C, Pedersen L G, Darden T A. Towards an accurate representation of electrostatics in classical force fields: Efficient implementation of multipolar interactions in biomolecular simulationsJ. The Journal of chemical physics, 2004, 120(1): 73-87.22 Nanoscience: Nanobiotechnology and

37、nanobiologyM. Springer Science & Business Media, 2009.23 Plimpton S. Fast parallel algorithms for short-range molecular dynamicsJ. Journal of computational physics, 1995, 117(1): 1-19.24 Hess B, Kutzner C, Van Der Spoel D, et al. GROMACS 4: algorithms for highly efficient, load-balanced, and scalabl

38、e molecular simulationJ. Journal of chemical theory and computation, 2008, 4(3): 435-447.25 Stone J E, Phillips J C, Freddolino P L, et al. Accelerating molecular modeling applications with graphics processorsJ. Journal of computational chemistry, 2007, 28(16): 2618-2640.26 Anderson J A, Lorenz C D,

39、 Travesset A. General purpose molecular dynamics simulations fully implemented on graphics processing unitsJ. Journal of Computational Physics, 2008, 227(10): 5342-5359.27 March N H, Tosi M P, March N H. Introduction to liquid state physicsM. Singapore: World Scientific, 2002.28 Abraham F F. An isot

40、hermalisobaric computer simulation of the supercooledliquid/glass transition region: Is the shortrange order in the amorphous solid fcc?J. The journal of chemical physics, 1980, 72(1): 359-365.29 Cheng Y Q, Ma E. Atomic-level structure and structureproperty relationship in metallic glassesJ. Progres

41、s in Materials Science, 2011, 56(4): 379-473.30 Qi L, Dong L F, Zhang S L, et al. Glass formation and local structure evolution in rapidly cooled Pd Ni alloy melt under high pressureJ. Physics Letters A, 2008, 372(5): 708-711.31 Kim T H, Kelton K F. Structural study of supercooled liquid transition metalsJ. The Journal of chemical physics, 2007, 126(5): 054513.32 Laio A, Va