使用Gromacs模擬碳納米管的一個簡單例子
注：本文的做法如今已經完全過時，如今強烈建議用Sobtop（http://www.shanxitv.org/soft/Sobtop）建立納米球/管/板的拓撲文件，比x2top靈活得多，參考網頁里的例子。
使用Gromacs模擬碳納米管的一個簡單例子
文/Sobereva @北京科音 2014-Dec-15
很早以前寫過一個帖子《amber與gromacs讀入碳納米管的方法》（http://www.shanxitv.org/39）介紹了怎么用gromacs處理碳納米管，不過貌似還是有很多人不會。這里就提供一個在真空中和在水中模擬碳納米管的簡單且完整的例子。這個小教程本來是給一個國際友人寫的，所以都是英文。本文用的是Gromacs 4.6.5，文中的CNT10指的是一個很普通的碳納米管，cnt10.pdb是其結構文件，這可以用Nanotube Modeler很容易地產生。
文中涉及到的文件在此：/usr/uploads/file/20150610/20150610055407_71312.rar。
Create a new file named atomname2type.n2t in "share/gromacs/top/gromos54a7.ff" folder in gromacs directory, and then copy below content into it
C CR1 0.0 12.011 3 C 0.142 C 0.142 C 0.142
C CR1 0.0 12.011 2 C 0.142 C 0.142
C CR1 0.0 12.011 1 C 0.142
That means when we use x2top to generate topology file under G54A7 forcefield, if a carbon atom has one or two or three neighbours with a distance about 0.142nm (normal C-C bond length in CNT), then this carbon will be recognized as CR1 atom, which is the atomtype corresponding to the carbon in aromatic CH-group of G54A7. This atomtype is suitable for representing vdW interaction between CNT carbons and environmental atoms.
Run below command to generate the CNT10.top:
g_x2top -f cnt10.pdb -o CNT10.top -ff select -nopbc -name CNT10 -kb 400000 -kt 600 -kp 150
Select G54A7 from the forcefield list, then CNT10.top will be yielded in current folder, and the molecular name is CNT10 (specified by "-name"). -nopbc have to be specified, otherwise x2top can't work normally. -kb, -kt and -kp are used to define the force constant of bond, angle and dihedral terms.
In the CNT10.top, the [ bonds ] field looks like below
1 3 1 1.420000e-01 4.000000e+05 1.420000e-01 4.000000e+05
2 3 1 1.420000e-01 4.000000e+05 1.420000e-01 4.000000e+05
2 4 1 1.420000e-01 4.000000e+05 1.420000e-01 4.000000e+05
2 27 1 1.420000e-01 4.000000e+05 1.420000e-01 4.000000e+05
3 218 1 1.420000e-01 4.000000e+05 1.420000e-01 4.000000e+05
...
The 4th column is the equilibrium length determined based on the input geometry (cnt10.pdb), the 5th column is the force constant set by -kb. The last two columns are redundant, you can simply ignore or even directly delete them. The content of [ angles ] and [ dihedrals ] fields are similar to [ bonds ].
The bond, angle and dihedral forcefield parameters we set above are not important, they are mainly used to keep the nearly rigid structure of CNT during simulation, so the force constant can be simply set to a large value; however, too large values will cause high-frequency vibrations and make the simulation unstable. If you want to make the CNT more flexible (rigid), you can decrease (increase) the dihedral force contant. Currently, AFAIK, no forcefield contains bond, angle and dihedral parameters specifically optimized for CNT modelling.
Then there are two cases, you can follow either one
1 MD simulation in vaccum
grompp -f md_vacuum.mdp -c cnt10.pdb -p CNT10.top -o md_vacuum.tpr
mdrun -v -deffnm md_vacuum
2 MD simulation in solvated box
editconf -bt triclinic -f cnt10.pdb -o cnt10_box.gro -d 2
genbox -cp cnt10_box.gro -cs spc216.gro -o cnt10_wat.gro -p CNT10.top
Add below sentence at the head of CNT10.top
#include "gromos54a7.ff/spce.itp"
Start simulation:
grompp -f md.mdp -c cnt10_wat.gro -p CNT10.top -o md.tpr
mdrun -v -deffnm md
模擬和石墨烯、碳球的過程和此文沒有任何區別，把文中的cnt10.pdb替換為相應的結構文件即可。石墨烯結構也可以由Nanotube Modeler產生。Nanotube Modeler自帶的富勒烯庫里有大量碳球結構，各種原子數的各種富勒烯異構體也可以用CaGe產生，CaGe的基本使用方法見《生成富勒烯的螺旋算法簡介以及使CaGe中的編號與Fowler-Manolopoulos編號相符的方法》（http://www.shanxitv.org/104）