gromacs讀入碳納米管的方法

1、與gromacs讀入碳納米管的方法讀入碳納米管的方法有人問(wèn)我leap怎么讀入碳納米管，這里說(shuō)一下，順便也說(shuō)一下gromacs讀入的方法。只說(shuō)操作方法，不涉及選用力場(chǎng)等問(wèn)題。碳納米管結(jié)構(gòu)比較大，而結(jié)構(gòu)很有規(guī)律，一般僅有C原子。若不考慮邊界問(wèn)題，可以將全部電荷當(dāng)成0，這種近似并不會(huì)對(duì)結(jié)果真實(shí)性有多少影響。對(duì)于碳納米管這樣的特征，適合當(dāng)成一般的小分子來(lái)處理，但不宜用antechamber而且也沒(méi)必要用。先用Nanotubemodeler生成一

2、段碳納米管，含240個(gè)C原子，導(dǎo)出結(jié)構(gòu)為tube.pdb。但是會(huì)發(fā)現(xiàn)，原子名都是相同的，而leap要求結(jié)構(gòu)中相同殘基中每個(gè)原子都有獨(dú)立的名字，所以需要先改名。比較方便的方法是用ultraedit，打開(kāi)后開(kāi)啟列模式，選定從第1個(gè)原子到最后一個(gè)原子的15和16列(即原子名C后面的兩列)，選ColumnNumber，OK，即把原子名改為了C1至C240，保存為tubename.pdb。這里假定C原子使用gaff力場(chǎng)的c2原子類型對(duì)應(yīng)的力場(chǎng)參數(shù)

3、，依次運(yùn)行xleapfleaprc.gaffa=loadpdbsobtubename.pdbbondbydistancea默認(rèn)是2埃內(nèi)成鍵，故碳納米管的碳原子都會(huì)正常成鍵，可用edita檢查。然后把碳納米管內(nèi)全部240個(gè)原子都設(shè)成c2原子類型，可以用批處理腳本。由于leap只能一條一條運(yùn)行批處理文件中的內(nèi)容而不支持shell腳本的循環(huán)，需要將循環(huán)轉(zhuǎn)化成普通命令腳本再在leap運(yùn)行。寫(xiě)腳本：f((i=1ileapdodone在shell下

4、運(yùn)行，得到leapdo，然后在leap里運(yùn)行sourceleapdo即可checka進(jìn)行參數(shù)檢查，應(yīng)該OK。如果想改力場(chǎng)參數(shù)或缺少某些參數(shù)，可以editgaff往里添加和修改。如果要加溶劑和普通情況無(wú)異，這里略過(guò)。最后saveamberparmatube.toptube.inpcrd=============gromacs讀入碳納米管還是用上面的tube.pdb為例子。小分子在gromacs中常用prodrg來(lái)處理，但是理由和antec

5、hamber一樣并不適合，而且有更方便的辦法。gromacs提供了一個(gè)構(gòu)建拓?fù)湮募墓ぞ選2top，專適用于構(gòu)建結(jié)構(gòu)規(guī)律性很強(qiáng)的體系。x2top中有一些bug和“規(guī)矩”，而且在不同gromacs版本中bug和“規(guī)矩”還不一樣，這里使用gromacs4.0.4的x2top。首先需要寫(xiě)n2t文件，比如使用ffgmx2力場(chǎng)，就在力場(chǎng)文件所在文件夾里寫(xiě)一個(gè)ffgmx2.n2t，內(nèi)容如下：CCX0.012.0113C0.14C0.14C0.14C

6、CX0.012.0112C0.14C0.14CCX0.012.0111C0.14第一行說(shuō)明如果體系中任何一個(gè)C原子，與周圍3個(gè)C原子的距離都在0.14nm左右（判斷標(biāo)準(zhǔn)大概為正負(fù)10%），就把它當(dāng)作CX原子類型（原子類型名字自定，可以不屬于力場(chǎng)包含的原子類型），電荷為0.0，相對(duì)原子質(zhì)量為12.011。第二行和第三行與第一行意義類似，即代表與兩個(gè)C原子和與一個(gè)C原子成鍵的C原子被當(dāng)成什么類型、多少電荷、多少原子質(zhì)量，此例中將碳納米管所有

7、原子都當(dāng)成CX。運(yùn)行x2topftube.pdbotube.topffgmx2這樣就得到tube.top，里面包含這個(gè)碳納米管全部鍵結(jié)項(xiàng)。ff代表用什么力場(chǎng)，此處即ffgmx2，如果ff則是出現(xiàn)列表自行選擇。這里如果使用ffgmx力場(chǎng)即便正確寫(xiě)了ffgmx.n2t也可能無(wú)法正確執(zhí)行，即不能正確根據(jù)距離判斷成鍵，此時(shí)應(yīng)換用別的力場(chǎng)來(lái)執(zhí)行x2top。如果運(yùn)行后程序卡住不動(dòng)，可嘗試加nopbc解決，這是一個(gè)bug。默認(rèn)情況下，x2top會(huì)自動(dòng)

8、在.top里相應(yīng)鍵結(jié)項(xiàng)加上力場(chǎng)參數(shù)，這種方式自動(dòng)加入的力場(chǎng)參數(shù)并不是根據(jù)力場(chǎng)中相應(yīng)原子類型間的鍵結(jié)參數(shù)加入的。平衡距離鍵角二面角就是tube.pdb中的相應(yīng)項(xiàng)的距離鍵角二面角，鍵的力常數(shù)鍵角力常數(shù)二面角力常數(shù)默認(rèn)是4000004005，可以分別用kb、kt、kp設(shè)定。如果不想讓其自動(dòng)加入，可以加上noparam。如果n2t中把C原子轉(zhuǎn)換為ffgmx2力場(chǎng)中已有的類型，比如CB，那么此時(shí)這個(gè).top搭配tube.pdb已經(jīng)可以直接用于模擬

9、了。此例將C設(shè)為力場(chǎng)中沒(méi)有的原子類型CX是為了假設(shè)需要加入專門的力場(chǎng)參數(shù)的情況，此時(shí)需要定義新的原子類型。首先加入VDW和14VDW參數(shù)，分別在ffgmx2nb.itp里的[atomtypes]和[nonbond_params]里面添加。如果前面x2top時(shí)加了noparam，還需要設(shè)定鍵結(jié)參數(shù)。分兩種情況：1直接在.top里面每一個(gè)鍵結(jié)項(xiàng)后面加入相應(yīng)力場(chǎng)參數(shù)，類似于x2top自動(dòng)在.top里添加力場(chǎng)參數(shù)的方式，用ultraedit等軟