許佩軍 唐媛媛,2 張 靜 張知博 王 昆

邵 穎3 沈虎峻2,* 毛英臣1,*

(1遼寧師范大學(xué)物理與電子技術(shù)學(xué)院,遼寧大連116029;2中國(guó)科學(xué)院大連化學(xué)物理研究所,遼寧大連116023; 3大連海事大學(xué)物理系,遼寧大連116026)

基于粗?；Ｐ蛯?duì)有機(jī)溶劑的分子動(dòng)力學(xué)模擬

許佩軍1唐媛媛1,2張 靜1張知博1王 昆1

邵 穎3沈虎峻2,*毛英臣1,*

(1遼寧師范大學(xué)物理與電子技術(shù)學(xué)院,遼寧大連116029;2中國(guó)科學(xué)院大連化學(xué)物理研究所,遼寧大連116023;3大連海事大學(xué)物理系,遼寧大連116026)

在基于Boltzmann分布對(duì)四種基本構(gòu)象進(jìn)行Monte Carlo取樣后,通過與全原子模型的范德華勢(shì)比較得到了Gay-Berne(GB)參數(shù).又在對(duì)用量化計(jì)算得到的分子體系的電勢(shì)進(jìn)行電荷、偶極矩和四極矩的擬合后,得到了電多極展開勢(shì)(EMP)參數(shù).利用得到的粗粒化參數(shù),基于粗?；Ｐ?對(duì)CHCl3及四氫呋喃(THF)兩種有機(jī)溶劑進(jìn)行了分子動(dòng)力學(xué)模擬(MDS),并將結(jié)果同全原子模擬進(jìn)行了比較.計(jì)算結(jié)果表明用粗粒化模型從整體上能重復(fù)全原子模型的模擬結(jié)果,但在某些細(xì)節(jié)的計(jì)算與全原子模型有偏差,其原因可能是目前工作僅考慮了單位點(diǎn)情況,為此今后在對(duì)具有復(fù)雜結(jié)構(gòu)的分子進(jìn)行粗?；M時(shí)還應(yīng)考慮合理放置及增加相互作用位點(diǎn).

粗粒化模型;Gay-Berne勢(shì);電多極展開勢(shì);徑向分布函數(shù);分子動(dòng)力學(xué)模擬

1 引言

分子動(dòng)力學(xué)模擬(MDS)已經(jīng)成為了研究復(fù)雜分子體系的重要工具,盡管一些計(jì)算方法和技術(shù)手段已經(jīng)取得了突破,1-3但當(dāng)前基于全原子模型對(duì)某些較大體系的動(dòng)力學(xué)模擬時(shí)間尺度仍然非常短,并不能滿足與實(shí)驗(yàn)比較的迫切需要.4-6已知在利用經(jīng)驗(yàn)力場(chǎng)進(jìn)行的全原子分子動(dòng)力學(xué)模擬主要包括兩部分:一是對(duì)分子的構(gòu)象進(jìn)行統(tǒng)計(jì)分析;二是結(jié)合相應(yīng)的構(gòu)象分析體系的能量及受力情況.7-9而對(duì)某一分子體系,由于非鍵結(jié)的范德華相互作用項(xiàng)和靜電相互作用項(xiàng)的數(shù)目遠(yuǎn)遠(yuǎn)大于鍵長(zhǎng)、鍵角和二面角相互作用項(xiàng)的數(shù)目,因此對(duì)范德華相互作用和靜電相互作用的計(jì)算將要耗費(fèi)大量的計(jì)算時(shí)間.為了簡(jiǎn)化對(duì)這兩種相互作用的計(jì)算,研究人員已提出了多種粗?；Ｐ?CG模型).10-14CG模型的主要特點(diǎn)在于它將研究體系中的分子或分子的某一部分近似看成一個(gè)原子團(tuán)簇,該原子團(tuán)簇內(nèi)部信息(鍵長(zhǎng)、鍵角、二面角等)被屏蔽,即該團(tuán)簇被處理為沒有內(nèi)部結(jié)構(gòu)的“粒子”,該粒子對(duì)應(yīng)于全原子模型中的單個(gè)原子.而在利用CG模型進(jìn)行分子動(dòng)力學(xué)模擬之前還必須獲得描述這些粒子間范德華相互作用和靜電相互作用的力場(chǎng)參數(shù),應(yīng)當(dāng)指出這一部分工作也是構(gòu)建CG模型的關(guān)鍵部分.CG模型結(jié)合了日益發(fā)展的計(jì)算機(jī)水平,使得對(duì)更大更復(fù)雜體系的長(zhǎng)時(shí)模擬成為可能.

該方法的主要缺陷就在于它的計(jì)算精度要低于基于全原子模型的分子動(dòng)力學(xué)模擬.10此外還需指出的是大多數(shù)的CG模型僅對(duì)某些特殊體系適用,而不能應(yīng)用于其他體系,即模型的可移植性有待提高.產(chǎn)生這一問題的主要原因在于多數(shù)CG模型將粗?；Ｗ涌闯墒乔蛐?這樣導(dǎo)致了這些模型沒能更好地描述分子間的范德華相互作用,以及未能對(duì)由此而導(dǎo)致的電勢(shì)空間分布改變的靜電相互作用作出更好的描述.為了克服這些缺陷,改進(jìn)CG模型,Golubkov和Ren14提出了在計(jì)算粗?；Ｗ娱g相互作用時(shí),必須考慮由粗粒化粒子的形狀而引起的Gay-Berne(GB)效應(yīng)及電多極矩效應(yīng),這些改進(jìn)在他們對(duì)水、苯及甲醇三種構(gòu)象相對(duì)簡(jiǎn)單的分子體系的動(dòng)力學(xué)模擬中得到了較理想體現(xiàn).此外也應(yīng)認(rèn)識(shí)到的是他們工作的基石就在于較好地?cái)M合得到了描述上述三種分子的粗粒化力場(chǎng)參數(shù).為了檢驗(yàn)如何利用CG模型對(duì)復(fù)雜分子體系進(jìn)行高效計(jì)算,本文采用類似方法對(duì)CHCl3及四氫呋喃(THF)分子的兩種構(gòu)象在相對(duì)復(fù)雜的有機(jī)溶劑中進(jìn)行了分子動(dòng)力學(xué)模擬.

2 計(jì)算與模擬方法

2.1 描述非球形分子的Gay-Berne勢(shì)和電多極展開勢(shì)(EMP)

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

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

其中

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

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

其中μ和ν均為可調(diào)指數(shù),借鑒Golubkov和Ren14的工作,本文分別將其取值為2.0和1.0.而ε1和ε2可分別表示為

(9)式中的χ?和α?分別取如下形式

其中,ε0描述了交叉結(jié)構(gòu)的阱深,εE和εS分別表示尾對(duì)尾和邊對(duì)邊結(jié)構(gòu)的阱深.當(dāng)兩個(gè)分子相同或者其中一個(gè)分子是球形分子時(shí),勢(shì)函數(shù)用Berne和Pechukas的原始形式表示.如果兩個(gè)分子都是球形的,勢(shì)函數(shù)將進(jìn)一步簡(jiǎn)化為L(zhǎng)ennard-Jones勢(shì).

長(zhǎng)距離靜電相互作用可以由電多極展開勢(shì)來表示.在大于分子尺寸的距離處,分子的電勢(shì)可準(zhǔn)確地通過對(duì)其質(zhì)心的電多極展開勢(shì)來表示23

其中q、μ、Θ分別表示電荷、偶極矩和四極矩.電多極勢(shì)可表示為電荷-電荷相互作用勢(shì)、電荷-偶極相互作用勢(shì)、電荷-四極相互作用勢(shì)等項(xiàng)之和.本文在分子的質(zhì)心上放置了一個(gè)電多極展開位點(diǎn),同時(shí)選取了和GB位點(diǎn)相同的慣性系.

兩個(gè)多極位點(diǎn)間的相互作用可通過電多極展開勢(shì)表示為

其矩陣形式可參見文獻(xiàn).14,24本文對(duì)靜電相互作用截?cái)嘀档奶幚矸椒?采用了TINKER軟件的標(biāo)準(zhǔn)方法.25-27

此外還需指出的是,考慮到由于粗粒化粒子具有各向異性的形狀而使得在進(jìn)入排斥的范德華相互作用區(qū)域前電多極展開方法將給出不正確的相互作用勢(shì),類似于文獻(xiàn),14,28本文同樣加入了一個(gè)阻尼函數(shù)以確保對(duì)短距離多極相互作用勢(shì)的精確描述.

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

2.2.1 Gay-Berne(GB)參數(shù)的擬合

本文利用粗粒化模型進(jìn)行分子模擬的前提是要獲得GB-EMP模型的參數(shù),為此我們首先對(duì)研究體系進(jìn)行了合理取樣選取,然后利用描述有機(jī)分子的AMBER力場(chǎng)(GAFF),29,30通過與全原子的范德華勢(shì)進(jìn)行比較后擬合得到了GB參數(shù).選取GAFF力場(chǎng)是由于它的函數(shù)形式簡(jiǎn)單,原子類型有限,還可利用實(shí)驗(yàn)和理論模型計(jì)算獲得力常數(shù)和部分原子電荷.在研究中,我們發(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)這四種構(gòu)象出現(xiàn)的幾率相對(duì)較大,17,19另外考慮到在真實(shí)的液體環(huán)境中分子的構(gòu)象趨于取能量最低的構(gòu)象,所以本文在篩選用于擬合GB參數(shù)的構(gòu)象時(shí),對(duì)上述四種構(gòu)象分別基于Boltzmann分布采用Monte Carlo方法獲得到分子的構(gòu)象集.這樣才可確保在一定的能量范圍內(nèi),構(gòu)象數(shù)目隨能量的升高而減少.表1給出了CHCl3和THF這兩種分子的四種基本構(gòu)象.

表1 氯仿和THF分子的四種基本構(gòu)象Table 1 Four reference configurations of CHCl3and THF molecules

圖1及圖2給出了氯仿和THF分子的四種基本構(gòu)象在各能量區(qū)間內(nèi)的構(gòu)象數(shù)目.從圖中我們可以清楚地看到所選的分子構(gòu)象在較低的能量范圍內(nèi)較多,且構(gòu)象數(shù)隨著能量的升高而逐漸減小.

表2 氯仿和THF分子的GB參數(shù)Table 2 GB parameters of CHCl3and THF molecules

在參數(shù)擬合的具體過程中,5個(gè)自由參數(shù)(l,d, ε0,εE/εS,dw)由遺傳算法擬合得到.表2給出用上述方案擬合得到的氯仿和THF分子的GB參數(shù).

2.2.2 EMP的參數(shù)擬合

在擬合粗?；Ｐ偷腅MP參數(shù)時(shí),本文首先應(yīng)用Gaussian 03軟件包31的MP2方法采用6-311G** (2d,2p)機(jī)組計(jì)算得到了兩種分子的電勢(shì),然后利用GDMA程序32對(duì)電勢(shì)進(jìn)行電荷、偶極和四極的擬合,最后得到了兩種分子粗?；Ｐ偷腅MP參數(shù).表3中給出了氯仿和THF的EMP參數(shù).

2.3 模擬體系與動(dòng)力學(xué)模擬方法

為了檢驗(yàn)本文對(duì)氯仿和THF分子粗粒化模型GB-EMP參數(shù)的擬合,我們利用粗?；Ｐ蛯?duì)兩種有機(jī)分子溶劑進(jìn)行了粗?；姆肿觿?dòng)力學(xué)模擬.為了在計(jì)算的過程中可任意地在粗粒化模型和全原子模型間進(jìn)行轉(zhuǎn)換,方便容易地比較CG模型和全原子模型的結(jié)果,故本文在計(jì)算中采用了TINKER軟件包.此外還基于如下考慮,即TINKER也包含了對(duì)粗?；Ｗ覩B-EMP能量、力和力矩部分的計(jì)算,并可在慣性系中非常容易地整合GB和EMP兩部分的計(jì)算結(jié)果.

氯仿和THF有機(jī)溶劑的立方體盒子由Materials Studio(MS)5.0軟件生成,其邊長(zhǎng)均為3 nm.對(duì)GB和EMP兩部分的截?cái)喟霃蕉既?.2 nm(這考慮到在通常情況下,要求截?cái)喟霃街挡淮笥诤凶舆呴L(zhǎng)的一半).CG模型和全原子模型的分子動(dòng)力學(xué)模擬都是采用NVT系綜在室溫(298 K)下進(jìn)行的.基于全原子模型的分子動(dòng)力學(xué)模擬的計(jì)算步長(zhǎng)取為1 fs.文獻(xiàn)14對(duì)考慮了GB-EMP效應(yīng)的CG模型的時(shí)間步長(zhǎng)作了實(shí)驗(yàn),推薦CG模型的計(jì)算步長(zhǎng)可取5-20 fs.考慮到由于自由度的減小,以及分子內(nèi)高頻運(yùn)動(dòng)的“屏蔽”,在本文我們?nèi)r(shí)間步長(zhǎng)為5 fs.在對(duì)運(yùn)動(dòng)方程進(jìn)行數(shù)值化處理時(shí),兩種模型的計(jì)算都采用了TINKER推薦的BERMAN積分算法.25

3 結(jié)果與討論

3.1 對(duì)范德華相互作用的分析

在圖3中,我們比較了分別利用CG模型和全原子模型計(jì)算得到的氯仿分子的范德華相互作用勢(shì),可以看到CG模型結(jié)果與全原子模型結(jié)果總體符合得較好.對(duì)面對(duì)面及交叉結(jié)構(gòu),兩種方法得到的作用曲線幾乎完全相同;而對(duì)于邊對(duì)邊及T型構(gòu)象,可以看出相比于全原子模型,粗?；Ｐ偷玫搅溯^小的阱深和距離參數(shù),這一現(xiàn)象對(duì)T型構(gòu)象尤為明顯.此外,基于粗?；Ｐ蛯?duì)四種構(gòu)象的計(jì)算結(jié)果整體趨勢(shì)較好地符合了Gay及Berne16的結(jié)論.

通過將分子構(gòu)象的取樣圖(圖1)與范德華相互作用勢(shì)能曲線圖(圖3)進(jìn)行比較,可以發(fā)現(xiàn)分子構(gòu)象數(shù)最多的能量區(qū)域正好對(duì)應(yīng)于范德華相互作用勢(shì)能曲線圖的最低點(diǎn)區(qū)域,這也印證了分子構(gòu)象都趨于位于能量最低的構(gòu)象狀態(tài)的結(jié)論.19

圖4給出了基于兩種模型對(duì)THF分子范德華相互作用的比較.從圖中可以看出,相比于結(jié)構(gòu)簡(jiǎn)單的氯仿分子,對(duì)THF分子進(jìn)行粗?；Ｐ偷挠?jì)算結(jié)果與全原子結(jié)果有較大的偏差.這表現(xiàn)在對(duì)應(yīng)于四種基本構(gòu)象的阱深參數(shù)和距離參數(shù)都一致性地小于全原子模型的計(jì)算結(jié)果.從圖3和圖4中,我們也觀察到兩種分子的交叉結(jié)構(gòu)和邊對(duì)邊結(jié)構(gòu)的范德華作用曲線非常接近,幾乎重合,這一現(xiàn)象很好地支持了文獻(xiàn)16的結(jié)論.在對(duì)粗粒化參數(shù)做多次擬合后,計(jì)算結(jié)果仍然如此,由此我們認(rèn)為出現(xiàn)這一現(xiàn)象的主要原因是我們僅考慮了單位點(diǎn)情況.這也建議在用粗?；Ｐ蛯?duì)復(fù)雜結(jié)構(gòu)的分子來計(jì)算分子間范德華相互作用時(shí)應(yīng)考慮增加相互作用位點(diǎn).

圖3 基于粗粒化模型和全原子模型對(duì)氯仿分子的范德華相互作用曲線的比較Fig.3 Comparison of van der Waals potential of CHCl3 based on the CG model and the all-atom model The curves of four configurations,namely face to face(lines 1 and2),side by side(lines 5 and 6),T-shape(lines 3 and 4),and cross shape(lines 7 and 8),are also displayed,respectively.

圖4 基于粗?；Ｐ秃腿幽Ｐ蛯?duì)THF分子的范德華相互作用曲線的比較Fig.4 Comparison of van der Waals potential of THF based on the CG model and the all-atom modelThe curves of four configurations,namely face to face(lines 1 and 2),side by side(lines 5 and 6),T-shape(lines 3 and 4)and cross shape(lines 7 and 8),are also displayed,respectively.The results of the CG model are represented with solid lines,and the dash lines are the results of the all-atom model.

考慮到粗粒化模型的普遍特征,10在目前階段的工作中,要進(jìn)一步檢驗(yàn)粗?；Ｐ偷挠?jì)算效率,還需結(jié)合其他物理量,如溶劑的徑向分布曲線等來做進(jìn)一步研究.

3.2 徑向分布函數(shù)的比較分析

需要指出,為了比較兩種模型得到的徑向分布函數(shù)曲線,首先需要把用粗?；Ｐ偷玫降腃G軌跡轉(zhuǎn)換為全原子軌跡.圖5給出了氯仿分子的徑向分布函數(shù)曲線.從圖中可以看出,盡管時(shí)間步長(zhǎng)增大為5 fs,粗?；姆肿觿?dòng)力學(xué)模擬仍能給出一條穩(wěn)定的軌跡,這說明CG模型能夠較好地重復(fù)全原子的模擬結(jié)果.與全原子模型相比,CG模型的徑向分布函數(shù)曲線起伏較少,這符合粗?；Ｐ偷幕咎卣?從C－H徑向分布曲線可看出,在0.65 nm處CG模型得到一波峰,此外從Cl－Cl及H－Cl徑向分布曲線還可看出CG模型不能重現(xiàn)0.30 nm處的波峰,但當(dāng)距離較大時(shí),CG模型和全原子模型的結(jié)果符合得較好.在距離較小處,兩者的結(jié)果符合得較差,其原因可能與EMP相互作用位點(diǎn)的放置有關(guān).在本文中,與GB相互作用位點(diǎn)的放置相同, EMP作用位點(diǎn)放置在粗粒化粒子的質(zhì)心上,這樣就造成了對(duì)電多極展開勢(shì)空間分布描述的偏差,從而影響了對(duì)電負(fù)性較強(qiáng)元素(Cl元素)的相關(guān)徑向分布函數(shù)曲線的模擬.而我們注意到Golubkov和Ren等33在最近實(shí)驗(yàn)工作中把EMP作用位點(diǎn)放到了其他位置上,相比于他們之前工作,14這樣更好地模擬了所研究體系的相關(guān)性質(zhì).這也提示在以后的工作中,應(yīng)針對(duì)分子的具體結(jié)構(gòu),考慮把EMP位點(diǎn)放在其他位置(如C原子上)而不是分子的質(zhì)心,或者考慮增加額外的相互作用位點(diǎn).

圖6給出了THF溶劑的徑向分布函數(shù)曲線.可以看出CG模型能較好地重復(fù)全原子的模擬結(jié)果,然而必須指出的是,與兩種方法得到的范德華相互作用曲線的較大差異不同,兩種模型的模擬結(jié)果符合得較好,沒有出現(xiàn)類似于氯仿的徑向分布曲線中的較大差異.鑒于分子結(jié)構(gòu)的差異,相比于將EMP相互位點(diǎn)放置在氯仿粗?；Ｗ拥馁|(zhì)心上,將該位點(diǎn)放置在THF粗粒化粒子的質(zhì)心上,目前我們認(rèn)為是一種可以接受的處理方式.對(duì)于CG模型仍然不能較好地重復(fù)全原子模型在0.30 nm處的較小起伏,我們認(rèn)為考慮增加相互作用位點(diǎn)或者適當(dāng)?shù)卣{(diào)整EMP作用位點(diǎn)的位置會(huì)提高粗?；Ｐ偷挠?jì)算精度.

圖5 對(duì)基于兩種模型得到的氯仿溶劑的徑向分布函數(shù)(RDF)曲線的比較Fig.5 Comparison of radial distribution functions(RDF)of CHCl3solvent based on the two models(a)C－H,(b)Cl－Cl,(c)H－Cl;Theresultsof theCG modelarerepresentedwithsolidlines,andthedash linesaretheresultsof theall-atommodel.

圖6 對(duì)基于兩種模型得到的THF溶劑的徑向分布函數(shù)(RDF)曲線的比較Fig.6 Comparison of radial distribution functions(RDF)of THF solvent based on the two models(a)C－H,(b)H－H,(c)O－H;The results of the CG model are represented with solid lines,and the dash lines are the results of the all-atom model.

通過CG模型進(jìn)行分子動(dòng)力學(xué)模擬的主要優(yōu)點(diǎn)是可以在保證一定精度的條件下較大幅度地提高計(jì)算的效率,這主要是由于用粗粒化的GB-EMP相互作用代替了多對(duì)原子間的相互作用,使得自由度數(shù)得到了較大地減少.此外一重要因素便是由于忽略了鍵長(zhǎng)、鍵角和二面角的高頻運(yùn)動(dòng),而使得模擬的時(shí)間步長(zhǎng)增大了.

還需指出的是,由于在處理Gay-Berne效應(yīng)中引入的dw參數(shù)較好地控制了GB勢(shì)的柔性,所以用本文中的CG模型還可用來模擬更大的復(fù)雜分子體系.此外在考慮了能更精確地描述靜電相互作用的電多極勢(shì)效應(yīng)后,使得基于CG模型模擬那些帶高電量的生物大分子體系成為可能.

4 結(jié)論

為了精確描述非球形分子間的相互作用,考慮了Gay-Berne效應(yīng)和電多極勢(shì)效應(yīng),基于粗?；Ｐ蛯?duì)氯仿和THF兩種有機(jī)溶劑進(jìn)行了分子動(dòng)力學(xué)模擬,并同全原子模擬結(jié)果進(jìn)行了比較.為此本文首先對(duì)研究體系的分子構(gòu)象進(jìn)行了合理取樣,然后利用描述有機(jī)分子的GAFF力場(chǎng),通過與全原子的范德華勢(shì)進(jìn)行比較后,用遺傳算法擬合得到了GB參數(shù).在擬合GB參數(shù)的過程中,我們也注意到構(gòu)象選取對(duì)計(jì)算結(jié)果的影響,針對(duì)CHCl3和THF兩種有機(jī)溶劑中面對(duì)面、邊對(duì)邊、T型及交叉結(jié)構(gòu)這四種構(gòu)象出現(xiàn)幾率相對(duì)較大,及分子構(gòu)象都趨于位于能量最低的構(gòu)象狀態(tài)的特征,利用Monte Carlo方法基于Boltzmann分布得到了擬合參數(shù)的構(gòu)象集,最后利用遺傳算法得到了描述粗?；Ｗ娱g相互作用的GB參數(shù).結(jié)合利用量化計(jì)算得到的EMP參數(shù),我們就獲得到了精確描述兩種有機(jī)分子的GB-EMP參數(shù).利用得到的粗?；瘏?shù),基于粗?；Ｐ蛯?duì)范德華相互作用曲線和徑向分布函數(shù)曲線進(jìn)行的計(jì)算表明粗粒化模型能以較高的效率從整體上重復(fù)全原子的計(jì)算結(jié)果,但是在某些細(xì)節(jié)方面,粗?；Ｐ筒荒芎芎玫胤先拥姆肿觿?dòng)力學(xué)模擬結(jié)果,這主要與目前工作描述相互作用的位點(diǎn)較少及相互作用位點(diǎn)放置簡(jiǎn)單有關(guān).為此我們應(yīng)在下一步工作中增加相互作用位點(diǎn),并應(yīng)適當(dāng)?shù)卣{(diào)整GB位點(diǎn)和EMP位點(diǎn)的具體位置.

(1) Dror,R.O.;Jensen,M.?.;Borhani,D.W.;Shaw,D.E.J.Gen. Physiol.2010,135,555.

(2) Karplus,M.;McCammon,J.A.Nature Struct.Biol.2002,9, 646.

(3) Shaw,D.E.;Maragakis,P.;Lindorff-Larsen,K.;Piana,S.;Dror, R.O.;Eastwood,M.P.;Bank,J.A.;Jumper,J.M.;Salmon,J. K.;Shan,Y.B.;Wriggers,W.Science 2010,330,341.

(4) van Gunsteren,W.F.;Bakowies,D.;Baron,R.;Chandrasekhar, I.;Christen,M.;Daura,X.;Gee,P.;Geerke,D.P.;Gl?ttli,A.; Hünenberger,P.H.;Kastenholz,M.A.;Oostenbrink,C.; Schenk,M.;Trzesniak,D.;van der Vegt,N.F.A.;Yu,H.B.B. Angew.Chem.Int.Edit.2006,45,4064.

(5) Klepeis,J.L.;Lindorff-Larsen,K.;Dror,R.O.;Shaw,D.E. Curr.Opin.Struct.Biol.2009,19,120.

(6) Freddolino,P.L.;Harrison,C.B.;Liu,Y.X.;Schulten,K. Nature Phys.2010,6,751.

(7)Xu,X.J.;Hou,T.J.;Qiao,X.B.;Zhang,W.Computer Aided Drug Design;Chemical Industry Press:Beijing,2004;pp 169-172.[徐筱杰,侯廷軍,喬學(xué)斌,章 威.計(jì)算機(jī)輔助藥物分子設(shè)計(jì).北京:化學(xué)工業(yè)出版社,2004:169-172.]

(8) Leach,A.P.Molecular Modeling,Principles and Application; Person Education Limited:England,2001;pp 165-245.

(9) Schlick,T.Molecular Modeling and Simulation:An Interdisciplinary Guide,2nd ed.;Springer:New York,2010; pp 265-343.

(10) Voth,G.A.Coarse-Graining of Condensed Phase and Biomolecular Systems;CRC Press:England,2009.

(11) Chu,J.W.;Izvekov,S.;Voth,G.A.Mol.Sim.2006,32,211.

(12) Marrink,S.J.;de Vries,A.H.;Mark,A.E.J.Phys.Chem.B 2004,108,750.

(13) Tozzini,V.Curr.Opin.Struc.Biol.2005,15,114.

(14) Golubkov,P.A.;Ren,P.Y.J.Chem.Phys.2006,125,064103.

(15) Berne,B.J.;Pechukas,P.J.Chem.Phys.1972,56,4213.

(16) Gay,J.G.;Berne,B.J.J.Chem.Phys.1981,74,3316.

(17) Cleaver,D.J.;Care,C.M.;Allen,M.P.;Neal,M.P.Phys.Rev. E 1996,54,559.

(18) Wilson,M.R.J.Chem.Phys.1997,107,8654.

(19) Care,C.M.;Cleaver,D.J.Rep.Prog.Phys.2005,68,2665.

(20) Paramonov,L.;Yaliraki,S.N.J.Chem.Phys.2005,123, 194111.

(21) Kabadi,V.N.;Steele,W.A.Ber.Bunsenges.Phys.Chem.1985, 89,2.

(22) Kabadi,V.N.Ber.Bunsenges.Phys.Chem.1986,90,327.

(23) Jackson,J.D.Classical Electrodynamics,3rd ed;John Wiley& Sons Inc.:New York,1999;pp 145-150.

(24) Applequist,J.J.Phys.A 1989,22,4303.

(25) Ponder,J.W.TINKER Molecular Modeling,Package 5.1; Washington University Medical School.

(26) Darden,T.;York,D.;Pedersen,L.G.J.Chem.Phys.1993,98, 10089.

(27) Sagui,C.;Pedersen,L.G.;Darden,T.A.J.Chem.Phys.2004, 120,73.

(28) Ren,P.Y.;Ponder,J.W.J.Phys.Chem.B 2003,107,5933.

(29)Wang,J.M.;Wolf,R.M.;Caldwell,W.J.;Kollman,P.A.;Case, D.A.J.Comput.Chem.2004,25,1157.

(30)Yang,L.J.;Tan,C.H.;Hsieh,M.J.;Wang,J.M.;Duan,Y.; Cieplak,P.;Caldwell,W.J.;Kollman,P.A.;Luo,R.J.Phys. Chem.B 2006,110,13166.

(31) Frisch,M.J.;Trucks,G.W.;Schlegel,H.B.;et al.Gaussian 03, RevisionA.01;Gaussian Inc.:Pittsburgh,PA,2003.

(32) Stone,A.J.J.Chem.Theory Comput.2005,1,1128.

(33)Golubkov,P.A.;Wu,J.C.;Ren,P.Y.Phys.Chem.Chem.Phys. 2008,10,2050.

Febraury 28,2011;Revised:May 9,2011;Published on Web:June 14,2011.

Molecular Dynamics Simulation of Organic Solvents Based on the Coarse-Grained Model

XU Pei-Jun1TANG Yuan-Yuan1,2ZHANG Jing1ZHANG Zhi-Bo1WANG Kun1SHAO Ying3SHEN Hu-Jun2,*MAO Ying-Chen1,*
(1School of Physics and Electronic Technology,Liaoning Normal University,Dalian 116029,Liaoning Province,P.R.China;
2Dalian Institute of Chemical Physics,Chinese Academy of Sciences,Dalian 116023,Liaoning Province,P.R.China;
3Department of Physics,Dalian Maritime University,Dalian 116026,Liaoning Province,P.R.China)

To obtain Gay-Berne(GB)parameters,we carried out Monte Carlo sampling of four reference configurations based on the Boltzmann distribution.After comparing with the van der Waals potential within the all-atom model we obtained the GB parameters.Also by fitting the charge,dipole,and quadrupole with the electric potential obtained from quantum chemical computations with Gaussian 03 we obtained the electric multipole potential(EMP)parameters.With the GB-EMP parameters we then carried out molecular dynamics simulations(MDS)for CHCl3and tetrahydrofuran(THF)based on the coarsegrained(CG)model.Compared with the all-atom model,the CG model can reproduce the simulation results on the whole,but there are some deviations in the simulations in some details.The reason is that we only take one interaction site into account in this work.Therefore,for more complicated molecules it is necessary to take the placement of the interaction sites into account.Additionally,the multi-sites situation is also considered in the MDS within the frame of the coarse-grained model.

Coarse-grained model;Gay-Berne potential;Electric multipole potential;Radial distribution function;Molecular dynamics simulation

O645;O641

*Corresponding authors.MAO Ying-Chen,Email:myc@lnnu.edu.cn;Tel:+86-411-82158367.SHEN Hu-Jun,Email:hshen@dicp.ac.cn; Tel:+86-411-84379875.

The project was supported by the Fundamental Research Funds for the Central Universities,China(2009QN069).

中央高校專項(xiàng)資金優(yōu)秀青年教師基金(2009QN069)資助項(xiàng)目