毛英臣，張德鵬，劉佳慧，孫 甍，王詩(shī)佳

（遼寧師范大學(xué) 物理與電子技術(shù)學(xué)院，遼寧 大連 116029）

粗?；Ｐ椭械腉ay－Berne相互作用勢(shì)

毛英臣，張德鵬，劉佳慧，孫 甍，王詩(shī)佳

（遼寧師范大學(xué) 物理與電子技術(shù)學(xué)院，遼寧 大連 116029）

合理描述非球形粗粒化粒子間的相互作用是提高粗?；肿觿?dòng)力學(xué)模擬速度的關(guān)鍵．為此本文介紹了 Gay-Berne勢(shì)．將之應(yīng)用于兩種有機(jī)小分子體系，在合理選擇構(gòu)象集后，由遺傳算法得到了粗?；Ｗ拥?GB參數(shù)，并通過(guò)對(duì)粗?；Ｐ秃腿幽Ｐ偷玫降姆兜峦叨瓜嗷プ饔玫膶?duì)比檢驗(yàn)了GB力場(chǎng)參數(shù)．最后，指出如何處理作用位點(diǎn)是粗?；Ｐ桶l(fā)展的一個(gè)關(guān)鍵問(wèn)題．

粗?；肿觿?dòng)力學(xué)模擬；Gay-Berne勢(shì)；范德瓦耳斯相互作用；作用位點(diǎn)

對(duì)生物大體系和軟物質(zhì)體系的分子模擬而言，隨著模擬體系粒子數(shù)目的增加，基于微觀層次的全原子模擬的運(yùn)算速度會(huì)極速下降，而若采用介觀模擬，某些微觀特征將會(huì)被“模糊”掉，故此介于兩種方法間的粗粒化分子模擬常被采用［1］．所謂“粗?；本褪茄芯矿w系中的小分子或者小原子集團(tuán)（如官能團(tuán)、殘基等）視為一個(gè) “粗?；Ｗ印保瑢?duì)于“粒子”相互作用則類似于全原子分子力場(chǎng)那樣進(jìn)行構(gòu)建，運(yùn)動(dòng)方程則一般常采用（質(zhì)點(diǎn)或剛體）牛頓方程．應(yīng)當(dāng)指出的是在實(shí)際模擬中，對(duì)分子間相互作用勢(shì)的合理描述是準(zhǔn)確而高效模擬體系性質(zhì)的關(guān)鍵所在．本文將介紹目前在粗?；肿幽M中采用的 Gay-Berne（GB）勢(shì)［2］，并結(jié)合兩種有機(jī)小分子（C5H5N、DMF：C3H7NO）體系進(jìn)行詳細(xì)的討論．

1 Gay－Berne相互作用勢(shì)

任何一個(gè)好的分子模型必須具備兩個(gè)主要特征，一方面數(shù)學(xué)上要簡(jiǎn)單并易于處理，另一方面能夠?qū)δM體系作較好的物理描述［3］．基于這兩方面的考慮，Berne和 Pechukas近似地將分子視為繞其主軸旋轉(zhuǎn)的橢球體，得到了由高斯函數(shù)描述的經(jīng)驗(yàn)勢(shì)．在對(duì)非球形分子的分子模擬中，短距離吸引和排斥相互作用可通過(guò)如下方法來(lái)完成，即在分子內(nèi)定義多個(gè)位點(diǎn)，每?jī)蓚€(gè)位點(diǎn)間的范德瓦爾斯相互作用勢(shì)由勒納德-瓊斯勢(shì)描述．需指出的是對(duì)較大分子，計(jì)算勢(shì)的時(shí)間會(huì)隨位點(diǎn)數(shù)目的平方而遞增，從而使得計(jì)算效率大為降低．為了解決這一問(wèn)題，Gay和 Berne對(duì)描述非球形分子間的短距離排斥和吸引相互作用的高斯重疊勢(shì)進(jìn)行了修正，把分子處理成軟的、單軸的橢球體，把 GB作用位點(diǎn)放置在橢球體的質(zhì)心上，用單位點(diǎn)勢(shì)替代多位點(diǎn)勢(shì)，這樣就得到了 GB勢(shì)［2］．在慣性系中，橢球體在位形空間中由其質(zhì)心坐標(biāo)和3個(gè)歐拉角描述，對(duì)應(yīng)于由一套 GB參數(shù)來(lái)描述．目前，GB勢(shì)在描述液晶體系中得到了廣泛應(yīng)用［4-6］．為了更好地控制 GB勢(shì)的軟硬度，Brene和Pechukas建議增加參數(shù)來(lái)控制經(jīng)驗(yàn)勢(shì)的柔性，由此引入了 dw參數(shù)．

兩個(gè)粗?；Ｗ娱g的GB勢(shì)可表示為

其中

而l和d分別描述分子的長(zhǎng)和寬，這樣分子的形狀就可以用任意的棒狀、盤狀或者球狀等來(lái)表示.

式（1）中的阱深參數(shù)取如下形式：

其中μ和ν均為可調(diào)指數(shù)，借鑒文獻(xiàn)［7］工作，本文分別將其取值為2.0和1.0.而ε1和ε2可分別表示為：

式（5）中的χ′和α′分別取如下形式：

其中 ε0描述了交叉型結(jié)構(gòu)的阱深，εE和 εS分別表示尾對(duì)尾和邊對(duì)邊結(jié)構(gòu)的阱深（形象描述參見(jiàn)表1）.當(dāng)兩個(gè)分子相同或者其中一個(gè)分子是球形分子時(shí)，勢(shì)函數(shù)用 Berne和 Pechukas的原始形式表示．如果兩個(gè)分子都是球形的，勢(shì)函數(shù)將進(jìn)一步簡(jiǎn)化為勒納德-瓊斯勢(shì)形式．

表1　C5H5N和DMF二聚物的4種基本構(gòu)象

2 粗?；Ｐ土?chǎng)參數(shù)的確定

為了獲得GB勢(shì)參數(shù)，我們需要對(duì)研究體系進(jìn)行合理取樣，然后利用描述有機(jī)分子的AMBER 力場(chǎng)（GAFF）［8］，通過(guò)與全原子的范德瓦爾斯相互作用勢(shì)進(jìn)行比較后擬合得到GB參數(shù)．選取GAFF力場(chǎng)是由于它的函數(shù)形式簡(jiǎn)單，原子類型有限．

在研究中，我們發(fā)現(xiàn)構(gòu)象的選取會(huì)直接影響GB參數(shù)的擬合，進(jìn)而決定了GB勢(shì)是否能夠準(zhǔn)確地模擬全原子的范德瓦爾斯相互作用勢(shì)，因此選擇可以反映真實(shí)分子液體環(huán)境的構(gòu)象是十分重要的．由于在實(shí)際的溶液環(huán)境中，二聚物的面對(duì)面（或尾對(duì)尾）、邊對(duì)邊、T型和交叉結(jié)構(gòu)這 4種構(gòu)象出現(xiàn)的概率相對(duì)較大［5］，另外考慮到在真實(shí)的液體環(huán)境中分子的構(gòu)象趨于取能量最低的構(gòu)象，所以本文在篩選用于擬合GB參數(shù)的構(gòu)象時(shí)，對(duì)上述4種構(gòu)象分別基于玻爾茲曼分布采用了蒙特卡洛方法獲得分子的構(gòu)象集．這樣才可以確保在一定的能量范圍內(nèi)，構(gòu)象數(shù)目隨能量的升高而減少．表1給出了C5H5N和 DMF二聚物的4種基本構(gòu)象．

在參數(shù)擬合的具體過(guò)程中，5個(gè)自由參數(shù)（l，d，ε0，εE/εS和 dw）由遺傳算法擬合得到，兩種有機(jī)小分子的GB參數(shù)見(jiàn)表2.

表2　C5H5N和DMF的GB參數(shù)

3 在有機(jī)小分子體系的應(yīng)用

為了檢驗(yàn)本文對(duì)分子粗?；Ｐ虶B參數(shù)的擬合，我們對(duì)C5H5N、DMF兩種有機(jī)分子溶劑分別進(jìn)行了粗?；暮腿拥姆肿觿?dòng)力學(xué)模擬．為了在計(jì)算的過(guò)程中可任意地在粗粒化模型和全原子模型間進(jìn)行轉(zhuǎn)換，方便比較粗粒化模型和全原子模型的結(jié)果，故本文在計(jì)算中采用了TINKER軟件．此外還基于如下考慮：TINKER還包含了對(duì)粗?；Ｗ拥腉B相互作用勢(shì)、力和力矩部分的計(jì)算，并可在慣性系中非常容易地整合GB相互作用勢(shì)的計(jì)算結(jié)果．兩分子間的靜電相互作用由電多極展開勢(shì)（EMP）計(jì)算．

C5H5N和 DMF有機(jī)溶劑的立方體盒子由Materials Studio（MS）5．0軟件生成，其邊長(zhǎng)均約為30 ?．對(duì)GB截?cái)喟霃蕉既?2．0 ?（這考慮到在通常情況下，要求截?cái)喟霃街挡淮笥诤凶舆呴L(zhǎng)的一半）．分子動(dòng)力學(xué)模擬在NVT系綜下進(jìn)行，溫度取為298 K．由于自由度的減小，以及分子內(nèi)高頻運(yùn)動(dòng)的 “屏蔽”，在粗?；肿觿?dòng)力學(xué)模擬中時(shí)間步長(zhǎng)取為5 fs．

在圖1和圖2中，我們比較了利用兩種模型計(jì)算得到的 C5H5N和 DMF的范德瓦爾斯相互作用勢(shì)，可以看到粗?；Ｐ徒Y(jié)果整體上較好地符合了全原子模型結(jié)果．但對(duì)于T型構(gòu)象，可以看出相比于全原子模型，粗?；Ｐ偷玫搅溯^小的阱深和距離參數(shù)，同時(shí)阱的位置也有一定的偏移．一方面主要是由于粗?；Ｐ陀脙蓚€(gè)橢球體（長(zhǎng)橢球或者扁橢球）并不能很好地表示這兩種分子二聚物的T型構(gòu)象．另一方面，更主要的是我們?cè)谟?jì)算時(shí)，將所用的作用位點(diǎn)都放在其質(zhì)心（中心）上，而實(shí)際上由于原子的電負(fù)性（C5H5N中的N原子，DMF中的N原子和O原子）往往會(huì)使靜電相互作用位點(diǎn)偏離其中心，且由于電多極矩效應(yīng)，需要用到多個(gè)作用位點(diǎn)．對(duì)比圖1與圖2，我們還發(fā)現(xiàn)DMF的這種偏移要較C5H5N更大．產(chǎn)生這樣結(jié)果的原因是，對(duì)于DMF分子而言，靜電作用位點(diǎn)與GB作用位點(diǎn)的距離要相較C5H5N分子更大，所以造成了勢(shì)阱位置的更大偏移．由圖1和圖2我們還可看到，對(duì)于二聚物的邊靠邊型和交叉型兩種構(gòu)象，利用兩種模型得到的范德瓦爾斯相互作用非常接近，這同樣是由于在目前的計(jì)算中僅考慮了單作用位點(diǎn)的原因．對(duì)于上述問(wèn)題，文獻(xiàn)［9］提出了多作用位點(diǎn)的GB-EMP模型，很好地改進(jìn)了計(jì)算結(jié)果．

圖1　基于粗粒化模型和全原子模型對(duì) C5H5N分子的范德瓦爾斯相互作用曲線的比較．二聚物的4種構(gòu)象：面對(duì)面（1、2），邊對(duì)邊（7、8），交叉型（5、6）和 T形（3、4）．實(shí)線代表了粗?；Ｐ徒Y(jié)果，虛線代表了全原子模型結(jié)果．

圖2　基于粗?；Ｐ秃腿幽Ｐ蛯?duì) DMF分子范德瓦爾斯相互作用曲線的比較（同圖1）

4 結(jié)論

加快分子動(dòng)力學(xué)模擬速度的一種有效的辦法是采用粗?；Ｐ停瑸榇诵枰侠砬腋咝У孛鑼懘至；Ｗ娱g的相互作用．一般地粗粒化粒子都不是球形的，為了精確描述非球形分子間的范德瓦耳斯相互作用，Gay-Berne勢(shì)將原來(lái)的球形分子改進(jìn)為對(duì)稱橢球形分子．本文詳細(xì)地介紹了Gay-Berne勢(shì)，首先對(duì)C5H5N和DMF兩種分子的二聚體進(jìn)行了合理取樣，然后由遺傳算法擬合得到了兩種分子的粗粒化粒子的GB參數(shù)．為了檢驗(yàn)所獲得的GB參數(shù)，我們對(duì)兩種有機(jī)小分子溶劑進(jìn)行了分子動(dòng)力學(xué)模擬，并同全原子模擬結(jié)果進(jìn)行了比較．粗粒化模型結(jié)果較好地符合了全原子模型結(jié)果．但是對(duì)于T型構(gòu)象，由于本工作將所有作用位點(diǎn)放在粗粒化粒子的中心上，所以勢(shì)阱參數(shù)和位置有了較大的偏差．更準(zhǔn)確地描述粒子間的相互作用，需考慮多作用位點(diǎn)，但應(yīng)當(dāng)指出的是作用位點(diǎn)的增加意味著計(jì)算速度的降低，因此如何平衡計(jì)算速度和精度是粗?；Ｐ桶l(fā)展必須考慮的一個(gè)問(wèn)題．

［1］ Voth G A．Coarse-Graining of Condensed Phase and Biomolecular Systems［M］．Boca Raton：CRC Press，2009．

［2］ Gay J G，Berne B J．Modification of the overlap potential to mimic a linear site-site potential［J］．Journal of Chemical Physics，1981，74：3316-3319．

［3］ Berne B J，Pechukas P．Gaussian model potentials for molecular interactions［J］．Journal of Chemical Physics，1972，56：4213-4216．

［4］ Cleaver D J，Care C M，Allen M P，et al．Extension and generalization of the Gay-Berne potential［J］．Physical Review E，1996，54：559-567．

［5］ Care C M，Cleaver D J．Computer simulation of liquid crystals［J］．Reports on Progress in Physics，2005，68：2665-2700．

［6］ 陳正隆，徐為人，湯立達(dá)．分子模擬的理論與實(shí)踐 ［M］．北京：化學(xué)工業(yè)出版社，2007：18-23．

［7］ Golubkov P A，Ren P Y．Generalized coarse-grained model based on point multipole and Gay-Berne potentials［J］．Journal of Chemical Physics，2006，125：064103．

［8］ Wang J M，Wolf R M，Caldwell W J，et al．Development and testing of a general AMBER force field［J］．Journal of Computational Chemistry，2004，25：1157-1174．

［9］ Golubkov P A，Wu J C，Ren P Y．A transferable coarsegrained model for hydrogen-bonding liquids［J］．Physical Chemistry Chemical Physics，2008，10：2050-2057．

Gay-Berne interaction potential used in coarse-graining molecular simulation

MAO Ying-chen，ZHANG De-peng，LIU Jia-hui，SUN Meng，WANG Shi-jia
（School of Physics and Electronic Technology，Liaoning Normal University，Dalian，Liaoning 116029，China）

How to compute effectively the interactions of non-spherical coarse-graining particles is critical for the improvement of computational speed of the coarse-grained molecular simulation，so we introduce Gay-Berne potential in this paper，and carry out applications for two small organic molecular systems．We derive the GB parameters of coarse-graining particles with genetic algorithm based on the reasonable conformational set．Moreover，we make some tests about the GB parameters through making comparison with the results from coarse-grained model and all-atom model，respectively．Finally，we point out how to deal with the interaction site is an essential problem for the development of coarse-graining model．

coarse-grained molecular dynamics simulation；Gay-Berne interaction；van der Waals potential；interaction site

教學(xué)討論

O 469

A

1000-0712（2016）10-0020-03

2015-12-05；

2016-02-25

國(guó)家自然科學(xué)基金項(xiàng)目（11447019）、教育部高等學(xué)校力學(xué)課程教學(xué)研究項(xiàng)目（JZW-15-LX-18）資助

毛英臣（1977—），男，山東萊州人，博士，遼寧師范大學(xué)物理與電子技術(shù)學(xué)院副教授，主要從事理論物理的研究和教學(xué)工作．