Amber與Gromacs讀入碳納米管的方法

這里說一下Amber和Gromacs讀入碳納米管的方法。只說操作方法，不涉及選用力場等問題。

Amber：

碳納米管結構比較大，而結構很有規律，一般僅有C原子。若不考慮邊界問題，可以將全部電荷當成0，這種近似并不會對結果真實性有多少影響。對于碳納米管這樣的特征，適合當成一般的小分子來處理，但不宜用antechamber而且也沒必要用。

先用Nanotube modeler生成一段碳納米管，含240個C原子，導出結構為tube.pdb。但是會發現，原子名都是相同的，而leap要求結構中相同殘基中每個原子都有獨立的名字，所以需要先改名。比較方便的方法是用ultraedit，打開后開啟列模式，選定從第1個原子到最后一個原子的15和16列(即原子名C后面的兩列)，選Column-Insert Number，OK，即把原子名改為了C1至C240，保存為tubename.pdb。

這里假定C原子使用gaff力場的c2原子類型對應的力場參數，依次運行
xleap -f leaprc.gaff
a=loadpdb /sob/tubename.pdb
bondbydistance a 默認是2埃內成鍵，故碳納米管的碳原子都會正常成鍵，可用edit a 檢查。
然后把碳納米管內全部240個原子都設成c2原子類型，可以用批處理腳本。由于leap只能一條一條運行批處理文件中的內容而不支持shell腳本的循環，需要將循環轉化成普通命令腳本再在leap運行。寫腳本：
for ((i=1;i<=240;i=i+1))
do
echo "set a.1.$i type c2" >> leapdo
done
在shell下運行，得到leapdo，然后在leap里運行source leapdo即可
check a 進行參數檢查，應該OK。如果想改力場參數或缺少某些參數，可以edit gaff往里添加和修改。
如果要加溶劑和普通情況無異，這里略過。
最后saveamberparm a tube.top tube.inpcrd

gromacs:

還是用上面的tube.pdb為例子。小分子在gromacs中常用prodrg來處理，但是理由和antechamber一樣并不適合，而且有更方便的辦法。gromacs提供了一個構建拓撲文件的工具x2top，專適用于構建結構規律性很強的體系。x2top中有一些bug和“規矩”，而且在不同gromacs版本中bug和“規矩”還不一樣，這里使用gromacs4.0.4的x2top。

首先需要寫n2t文件，比如使用ffgmx2力場，就在力場文件所在文件夾里寫一個ffgmx2.n2t，內容如下：
C    CX    0.0    12.011     3    C 0.14   C 0.14  C 0.14
C    CX    0.0    12.011     2    C 0.14   C 0.14
C    CX    0.0    12.011     1    C 0.14
第一行說明如果體系中任何一個C原子，與周圍3個C原子的距離都在0.14nm左右（判斷標準大概為正負10%），就把它當作CX原子類型（原子類型名字自定，可以不屬于力場包含的原子類型），電荷為0.0，相對原子質量為12.011。第二行和第三行與第一行意義類似，即代表與兩個C原子和與一個C原子成鍵的C原子被當成什么類型、多少電荷、多少原子質量，此例中將碳納米管所有原子都當成CX。

運行x2top -f tube.pdb -o tube.top -ff gmx2
這樣就得到tube.top，里面包含這個碳納米管全部鍵結項。
-ff代表用什么力場，此處即ffgmx2，如果-ff select則是出現列表自行選擇。這里如果使用ffgmx力場即便正確寫了ffgmx.n2t也可能無法正確執行，即不能正確根據距離判斷成鍵，此時應換用別的力場來執行x2top。
如果運行后程序卡住不動，可嘗試加-nopbc解決，這是一個bug。
默認情況下，x2top會自動在.top里相應鍵結項加上力場參數，這種方式自動加入的力場參數并不是根據力場中相應原子類型間的鍵結參數加入的。平衡距離/鍵角/二面角就是tube.pdb中的相應項的距離/鍵角/二面角，鍵的力常數/鍵角力常數/二面角力常數默認是400000/400/5，可以分別用-kb、-kt、-kp設定。如果不想讓其自動加入，可以加上-noparam。

如果n2t中把C原子轉換為ffgmx2力場中已有的類型，比如CB，那么此時這個.top搭配tube.pdb已經可以直接用于模擬了。
此例將C設為力場中沒有的原子類型CX，是因為假設要給它設定專門的力場參數。首先加入VDW參數，在ffgmx2nb.itp里的[atomtypes]和[nonbond_params]里面添加。
如果前面x2top時加了-noparam，還需要設定鍵結參數。分兩種情況：
1 直接在.top里面每一個鍵結項后面加入相應力場參數，類似于x2top自動在.top里添加力場參數的方式，用ultraedit等軟件的列模式批量加入效率最高。
2 不在.top里面加，而是在ff???bon.itp里面的[bondtypes] [angletypes] [dihedraltypes]里面加入相應的項，或者把這些內容單獨寫進一個.itp文件，之后include進.top里。這種方法比較省事、易于修改。
關于這兩種參數定義方式的討論，詳見《gromacs決定力場參數的方式》（http://www.shanxitv.org/4）

之后加盒子、填充溶劑等與一般體系無異。