首页/文章/ 详情

水结冰模拟及TIP4PICE水模型使用

刘十三613

2月前浏览1388

以下文章来源于网络，如有侵权,请联系管理员删除!

本教程内容适用于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.;[ moleculetype ]; molnamenrexclSOL2 [ atoms ]; idat typeres nr residu nameat namecg nrcharge1       O        1       SOL      O     1       0.002       H        1       SOL      H     1       0.58973       H        1       SOL      H     1       0.58974       M        1       SOL      M     1      -1.1794 [constraints];i j funct doh  dhh1       2       1       0.095721       3       1       0.095722       3       1       0.15139 [ exclusions ]1234213431244123 ; 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) [ virtual_sites3 ]; Vsite fromfunctab412310.13458335      0.13458335

复制以上内容另存为tip4p-ice.itp，保存到力场文件夹中（以OPLS-AA为例，gromacs根目录/share/gromacs/top/oplsaa.ff），然后就可以在拓扑文件中加入

#include "oplsaa.ff/forcefield.itp"#include "oplsaa.ff/tip4p-ice.itp"

来使用新的水模型。

注意你需要根据使用的力场来修改以上文件中的at type，如对于OPLS-AA，原子类型修改为

[ atoms ]; idat typeres nr residu nameat namecg nrcharge1       opls_113        1       SOL      O     1       0.002       opls_114        1       SOL      H     1       0.58973       opls_114        1       SOL      H     1       0.58974       opls_115        1       SOL      M     1      -1.1794

总之就是要修改为这个力场中有定义且定义正确的原子类型。你可以参照力场文件夹中其它水模型使用的原子类型来修改。

当然，你也可以直接简单粗暴地将以上水的拓扑文件直接放在当前模拟的拓扑文件中。

没有适当模板的结晶模拟是非常困难的，因此在本模拟中使用冰-Ih/水两相体系。

使用packmol建模，输入代码如下。上下填充两层冰，中间充满水

tolerance 2.0output icebox.pdbfiletype 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 2gmx 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 2gmx mdrun -v -deffnm em

由于之前已经在三点水下进行过能量最小化，所以这一步能量最小化只需要调整虚原子位置即可，一步优化就完成了。

这样就得到了基于TIP4P的冰/水坐标文件。

模拟照常在能量最小化之后再进行NVT、NPT预平衡，然后进行产生相模拟。这里贴出产生相模拟所使用的mdp文件，其它mdp文件可在开头链接中的附件获取。

; Run parametersintegrator  = md        ; leap-frog integratornsteps      = 100000000  ; 1 * 100000000 = 10000 ps (10 ns)dt          = 0.001     ; 1 fs; Output controlnstxout     = 4000       ; save coordinates every 0.4 psnstvout     = 4000       ; save velocities every 0.4 psnstenergy   = 4000       ; save energies every 0.4 psnstlog      = 4000       ; update log file every 0.4 psnstxtcout   = 40000      ; xtc compressed trajectory output every 40 ps; Bond parameterscontinuation            = yes       ; Restarting after NVTconstraint_algorithm    = lincs     ; holonomic constraintsconstraints             = h-bonds   ; only h-bonds constrainedlincs_iter              = 1         ; accuracy of LINCSlincs_order             = 4         ; also related to accuracy; Neighborsearchingns_type     = grid      ; search neighboring grid celsnstlist     = 5         ; 5 fsrlist       = 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); Electrostaticscoulombtype     = PME       ; Particle Mesh Ewald for long-range electrostaticspme_order       = 4         ; cubic interpolationfourierspacing  = 0.16      ; grid spacing for FFT; Temperature coupling is ontcoupl      = v-rescale                       ; More accurate thermostattc-grps     = ICE        SOL      ; two coupling groups - more accuratetau_t       = 0.5           0.5                 ; time constant, in psref_t       = 260           260                 ; reference temperature, one for each group, in K; Pressure coupling is onpcoupl      = Parrinello-Rahman     ; Pressure coupling on in NPTpcoupltype  = semiisotropic         ; uniform scaling of x-y box vectors, independent ztau_p       = 5.0                   ; time constant, in psref_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 conditionspbc         = xyz       ; 3-D PBC; Dispersion correctionDispCorr    = EnerPres  ; account for cut-off vdW scheme; Velocity generationgen_vel     = no        ; Velocity generation is off; COM motion removal; These options remove motion of the ice relative to the solvent/ionsnstcomm         = 1comm-mode       = Linearcomm-grps       = ICE SOL; Scale COM of reference coordinatesrefcoord_scaling = com