以下文章来源于网络,如有侵权,请联系管理员删除!
本教程内容适用于Gromacs,包括:
1. TIP4P/ICE水模型导入方法;
2. 冰水建模;转换TIP3P(SPC/E)水模型坐标文件为TIP4P的方法;
3. 两相共存状态下水结冰过程的模拟;
TIP4P是一种四点水模型,亦即除了一个O和2个H原子外还加入了一个虚原子(在氧附近),虚原子仅具有负电荷,改善了水分子周围的静电分布
其中,TIP4P/ICE是专门用来描述冰的一系列性质的水模型,对水的相图有较好的还原(1bar下冰-Ih熔点为272.2K)。而Gromacs的力场中没有自带该水模型,需要到gromacs官网获取(内容如下)。
;
; Note the strange order of atoms to make it faster in gromacs.
;
[ ]
; molnamenrexcl
SOL2
[ ]
; idat typeres nr residu nameat namecg nrcharge
1 O 1 SOL O 1 0.00
2 H 1 SOL H 1 0.5897
3 H 1 SOL H 1 0.5897
4 M 1 SOL M 1 -1.1794
[ ]
;i j funct doh dhh
1 2 1 0.09572
1 3 1 0.09572
2 3 1 0.15139
[ ]
1234
2134
3124
4123
; The position of the virtual site is computed as follows:
;
;O
;
; D
;
;HH
;
; const = distance (OD) / [ cos (angle(DOH)) * distance (OH) ]
; 0.01577 nm/ [ cos (52.26 deg)* 0.09572 nm]
; Vsite pos x4 = x1 + a*(x2-x1) + b*(x3-x1)
[ ]
; Vsite fromfunctab
412310.13458335 0.13458335
复 制以上内容另存为tip4p-ice.itp,保存到力场文件夹中(以OPLS-AA为例,gromacs根目录/share/gromacs/top/oplsaa.ff),然后就可以在拓扑文件中加入
来使用新的水模型。
注意你需要根据使用的力场来修改以上文件中的at type,如对于OPLS-AA,原子类型修改为
[ atoms ]
; idat typeres nr residu nameat namecg nrcharge
1 opls_113 1 SOL O 1 0.00
2 opls_114 1 SOL H 1 0.5897
3 opls_114 1 SOL H 1 0.5897
4 opls_115 1 SOL M 1 -1.1794
总之就是要修改为这个力场中有定义且定义正确的原子类型。你可以参照力场文件夹中其它水模型使用的原子类型来修改。
当然,你也可以直接简单粗暴地将以上水的拓扑文件直接放在当前模拟的拓扑文件中。
没有适当模板的结晶模拟是非常困难的,因此在本模拟中使用冰-Ih/水两相体系。
使用packmol建模,输入代码如下。上下填充两层冰,中间充满水
tolerance 2.0
output icebox.pdb
filetype pdb
structure ice.pdb
number 1
inside box 0. 0. 0. 28. 24. 23.
end structure
structure ice.pdb
number 1
inside box 0. 0. 46. 28. 23. 69.
end structure
structure water.pdb
number 600
inside box 0. 0. 23. 28. 23. 46.
end structure
得到的icebox.pdb若要直接用于模拟,还需要编辑模拟盒子信息,在最后一行REMARK后加入以下信息
CRYST1 28 23 69 90.00 90.00 90.00 P 1 1
该信息记载了模拟盒子的晶胞参数(三边长,角度等),可被gromacs读取
至此我们得到了三点水模型的坐标文件,先在spc/e模型下进行能量最小化
gmx grompp -f em.mdp -c icebox.pdb -p icebox.top -o em.tpr -maxwarn 2
gmx mdrun -v -deffnm em
这样就得到了经过能量最小化的em.gro。接着,我们要在不改变原坐标的情况下,给所有水引入虚原子,转化为四点水模型的坐标文件。运行
cat em.gro |awk '{print $0;if($2=="OW"){a=$1;b=$3;c=$4;d=$5;e=$6;}if($2=="HW2")printf("%8s MW %4d%8.3f%8.3f%8.3f\n",a,b,c,d,e)}' > newicebox.gro
这段代码意思是逐行读取em.gro文件,读到氧原子一行时,将氧原子的编号和坐标写入临时变量abcde,然后读到每个水分子的第二个氢原子时在下一行插入虚原子,编号和坐标与第一个氧原子相同。亦即在每个水的氧原子位置添加一个虚原子。
然后在修改过的使用tip4p/ice水模型的拓扑文件下重新进行能量最小化(由于新引入了虚原子,你需要手动修改gro文件的原子数量)
gmx grompp -f em.mdp -c newicebox.gro -p icebox4.top -o em.tpr -maxwarn 2
gmx mdrun -v -deffnm em
由于之前已经在三点水下进行过能量最小化,所以这一步能量最小化只需要调整虚原子位置即可,一步优化就完成了。
这样就得到了基于TIP4P的冰/水坐标文件。
模拟照常在能量最小化之后再进行NVT、NPT预平衡,然后进行产生相模拟。这里贴出产生相模拟所使用的mdp文件,其它mdp文件可在开头链接中的附件获取。
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 100000000 ; 1 * 100000000 = 10000 ps (10 ns)
dt = 0.001 ; 1 fs
; Output control
nstxout = 4000 ; save coordinates every 0.4 ps
nstvout = 4000 ; save velocities every 0.4 ps
nstenergy = 4000 ; save energies every 0.4 ps
nstlog = 4000 ; update log file every 0.4 ps
nstxtcout = 40000 ; xtc compressed trajectory output every 40 ps
; Bond parameters
continuation = yes ; Restarting after NVT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; only h-bonds constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
ns_type = grid ; search neighboring grid cels
nstlist = 5 ; 5 fs
rlist = 1.0 ; short-range neighborlist cutoff (in nm)
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = v-rescale ; More accurate thermostat
tc-grps = ICE SOL ; two coupling groups - more accurate
tau_t = 0.5 0.5 ; time constant, in ps
ref_t = 260 260 ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT
pcoupltype = semiisotropic ; uniform scaling of x-y box vectors, independent z
tau_p = 5.0 ; time constant, in ps
ref_p = 450.0 450.0 ; reference pressure, x-y, z (in bar)
compressibility = 4.5e-5 4.5e-5 ; isothermal compressibility, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation is off
; COM motion removal
; These options remove motion of the ice relative to the solvent/ions
nstcomm = 1
comm-mode = Linear
comm-grps = ICE SOL
; Scale COM of reference coordinates
refcoord_scaling = com
在这里,步长dt=1fs, 使用V-rescale控温,参考温度260K, Parrinello-Rahman 控压, 参考压力为 450bar 。
使用了半各向异性控压,并将冰和水分开来控温——不过这也许并不必要。
以上条件显然不算最优的,还需要按照自己需求进行调整参数。我实际模拟了300ns,能够观察到结冰过程
除了直观地观察冰晶的形成情况之外,还可以通过统计氢键平均寿命来衡量结晶的程度。