在Gromacs中模擬金納米線拉伸過程（含視頻）

1 前言

前一陣子在多尺度材料模擬會議上看見有人演示了一段拉伸納米線的動畫，令我有些愉悅。雖然寡人并不是搞材料模擬的，不過還是想模擬一下玩玩，就用熟悉的gromacs做了一下拉伸金納米線的過程，結果看上去還挺合理，這里把流程介紹下。當然了，納米線拉伸這類問題用lammps去做最適合，gromacs做這個就是寓教于樂罷了。

gromacs模擬金屬最大的困難就是不支持模擬金屬最常用的EAM勢及其變種，gromacs只能直接支持對兒勢，然而一般普遍認為對兒勢對于模擬金屬不好。不過，這個共識也并不完全正確，對兒勢的形式用于金屬，至少是特定的金屬，不代表結果一定不好，只要能把參數搞好的話。在JPCC,112,17281就給出了一套很好的LJ 12-6或9-6對兒勢形式的適合fcc金屬（包括Ag, Al, Au, Cu, Ni, Pb, Pd, Pt）的參數，結果不遜于EAM勢，而且計算量還大大降低，還可以直接結合AMBER, CHARMM, COMPASS, CVFF, OPLS-AA, GROMOS, PCFF這些用LJ勢的分子力場在amber、gromacs、NAMD、Forcite等程序中研究與金屬混合的分子體系，比如研究水在金屬表面的行為。

本文用的gromacs為4.5.5版。本文涉及的文件都可以從這里下載：http://pan.baidu.com/s/1dFeEGEX

獲得的拉伸過程的軌跡如下所示：

也順便做了一下擠壓納米線的軌跡

2 建模

我們首先建立一個金的納米線的結構，具體來說，就是個25nm*6nm*6nm的金的圓棒。

最方便的方法是用Material studio+VMD來建。首先在MS中選File-import，在...Materials Studio 6.0\share\Structures\metals\pure-metals目錄下選中Au.msi，可知此晶胞的邊長為4.078埃。我們要將它平移復制成一個金屬條，然后從中挖出一個金屬棒。我們選Build-symmetry-Supercell，為了能達到25nm*6nm*6nm的尺寸，在三個框分別填65,15,15。選File-export，將得到的超胞導出到Au.pdb文件中。

用VMD打開Au.pdb，我們先找出這個超胞的中心。在控制臺輸入
atomselect top all
measure center atomselect0
得到中心的坐標為（埃）：
131.52516174316406 29.5676326751709 29.5676326751709

將Au.pdb保存為新的pdb文件Au_rod.pdb，在保存時的selected atoms框里面輸入：
abs(x-131.525)<125 and (y-29.567)^2+(z-29.567)^2<30^2
這樣就以超胞中心作為納米線的中心，挖出來了25nm*6nm*6nm的圓柱形。

將剛保存出來的Au_rod.pdb載入VMD。此納米線包含41480個原子。在gromacs中進行拉伸時，我們要將下圖中黃色的區域原子同時向兩側拉，每個黃色區域厚12埃。（因為我們按照原文用12埃的cut-off方式計算VDW勢，因此這樣可以讓黃色區域外的原子都不會感到邊緣是被人為截斷的）
[1.png]
為了能在gromacs中拉它們，必須將左右黃色區域的原子的序號分別寫進index文件中，作為名為left和right的group。為此，我們先找出這個棒的最左端和最右端的坐標，即x最小值和x最大值。運行：
measure minmax atomselect0
返回結果
{8.156999588012695 0.0 0.0} {254.8939971923828 59.1349983215332 59.1349983215332}

下面這個tcl小腳本是將指定區域的原子的序號按照index文件允許的格式輸出到c:\a.txt中。將這段腳本拷進控制臺運行
proc rangeindex {range {fps 0}} {
set sel [atomselect top $range frame $fps]
set result [open c:\\a.txt w]
set k 0
foreach i [$sel list] {
incr k
puts -nonewline $result [expr $i+1]\
if {$k%6==0} {puts $result " "}
}
close $result
$sel delete
puts "total $k atoms"
}

接下來運行rangeindex {x<12+8.156999588012695}。將生成的a.txt改名為left.txt
再運行rangeindex {x>254.8939971923828-12}。將a.txt改名為right.txt

把Au_rod.pdb里"AU1 MOL"全都替換為"AU  AU "。這代表每個金原子作為一個分子對待，金原子的原子名和分子名都叫AU。

3 在gromacs中設定參數文件并模擬

將下面這行
   Au  196.9665
添加到gromacs/top/gromos53a6.ff/atomtypes.atp當中

編寫一個拓撲文件，名為Au_rod.top，內容如下。其中金的參數來自JPCC,112,17281的表1的Au的B、A值。文中這兩個參數的能量和長度單位是kcal和埃，這里轉換為了gromacs用的KJ和nm。
#include "gromos53a6.ff/forcefield.itp"

[ atomtypes ]
   AU   79      0.000      0.000     A  0.029247582 9.657102e-06

#include "gromos53a6.ff/spce.itp"

[ moleculetype ]
; molname   nrexcl
AU         1

[ atoms ]
; id    at type res nr  residu name at name  cg nr  charge   mass
1       AU    1       AU         AU       1      0.0        196.9665

[ system ]
; Name
Au_rod

[ molecules ]
; Compound        #mols
AU  41480
實際上我們沒必要把整個gromos53a6的力場文件以及spce.itp給include進去，但是這樣這個拓撲文件更為普適，以后可以直接往當前體系里加入水、小分子之類的東西一起模擬。如果不把gromos53a6力場include進去的話，至少得把gromos53a6.ff/forcefield.itp中的[ defaults ]段落替換到#include "gromos53a6.ff/forcefield.itp"所在的行。

這里之所以用gromos力場，而不是OPLS、amber、charmm27等，是因為gromos力場在gromacs中用的是C12/r^12-C6/r^6形式，其中金的C12和C6參數在JPCC,112,17281里面給出了。雖然這篇文章中也給了epsilon+r0形式的LJ勢的參數，但是函數形式不標準，和gromacs中的形式并不一樣，而OPLS、amber、charmm27等用的都是epsilon+r0參數表達LJ勢，因此不兼容而沒法直接用。（注：決定用C12+C6還是epsilon+r0參數表達LJ勢取決于力場的forcefield.itp文件里的comb-rule項）

運行Make_ndx -f Au_rod.pdb
把left.txt和right.txt里的原子序號都寫到index.ndx最后，開頭設的group名分別為[left]和[right]。

構建md1.mdp，用于平衡模擬，內容如下
define =
integrator = md
dt = 0.002
nsteps = 50000 !100ps
comm_mode = ANGULAR
nstcomm = 10
comm-grps = system
nstxout = 1000
nstvout = 1000
nstfout = 0
nstlog = 500
nstenergy = 500
nstxtcout = 1000
xtc_grps = system
freezegrps = left right
freezedim = Y Y Y Y Y Y
;
nstlist = 10
ns_type = grid
pbc = no
rlist = 1.2
coulombtype = cut-off
rcoulomb = 1.2
vdwtype = cut-off
rvdw = 1.2
DispCorr = no
;
Tcoupl = Berendsen
tau_t = 5
tc_grps = system
ref_t = 300
本例是在真空下模擬。我們先做100ps的平衡模擬，讓體系升溫到300K，在此過程中納米棒兩端1.2nm厚度的原子都凍結不動。

運行
grompp -f md1.mdp -c Au_rod.pdb -p Au_rod.top -o Au_md1.tpr -maxwarn 5 -n index.ndx
mdrun -v -deffnm Au_md1
通過觀看軌跡，會看到此體系很結實，原子排布沒有變化，只是由于熱運動原子有微微晃動。溫度也早已平衡到300K了。

接下來編輯拉伸過程的mdp文件。將md1.mdp復制為SMD.mdp，把nsteps設為250000，即模擬500ps。freezedim設為N Y Y N Y Y，即拉伸中只讓兩端1.2nm厚度內的原子在拉伸方向運動（否則可能斜著散架）。在末尾添加拉伸的設定：
pull=constraint
pull_dim=Y N N
pull_geometry=distance
pull_group0=left
pull_group1=right
pull_ngroups=1
pull_init1=23.8
pull_rate1=0.025
這代表勻速x方向拉伸，讓名為left和right的group的質心距離從23.8nm（這是開始時它們的質心在當前結構下實際的距離）開始以0.025nm/ps的速度逐漸加大（折合每秒25m...不要介意，本文目的只是帶來愉悅）。這樣在500ps內可以拉長12.5nm，即當前體系長度的一半。

運行
grompp -f SMD.mdp -c Au_md1.gro -p Au_rod.top -o Au_SMD.tpr -maxwarn 5 -n index.ndx
mdrun -v -deffnm Au_SMD
最終的軌跡應當是類似于前面的視頻那樣。如預期的，拉開后斷面是尖的。有興趣的話可以用JPCC,112,17281里的其它金屬的參數也拉伸其它金屬，或者將納米線的直徑改大點再跑跑，以及看看不同速度、溫度下拉伸結果的差異。

我們再做個擠壓過程，也就是將pull_rate1=0.025改成pull_rate1=-0.025，模擬900ps，這樣left和right的質心距離就剩23.8-900*0.025=1.3nm了，差0.1nm就相接了。