Gromacs制作能量收敛的水分子层

制作可用于计算的介电常数的水分子层

总的来说,需要用溶解的方式制作gro文件,然后自己制作top文件,如果直接用pdb2gmx和用packmol制作的由单体水分子构成的分子层会使得分子之间成H键。

制作gro文件

gmx solvate -box 11 11 2 -cs spc216.gro -o water_test.gro

盒子的大小是配合蛋白的最小外接和盒子尺寸的。

下面展示water_test.gro的前几行

Generated by gmx solvate
22188
    1SOL     OW    1   0.230   0.628   0.113
    1SOL    HW1    2   0.137   0.626   0.150
    1SOL    HW2    3   0.231   0.589   0.021
    2SOL     OW    4   0.225   0.275   0.996
    2SOL    HW1    5   0.260   0.258   1.088
    2SOL    HW2    6   0.137   0.230   0.984
    3SOL     OW    7   0.019   0.368   0.647
    3SOL    HW1    8  -0.063   0.411   0.686
    3SOL    HW2    9  -0.009   0.295   0.584
    4SOL     OW   10   0.569   1.275   1.165

这说明有22188个原子

制作top文件

22188/3=7396
所以这样写

#include "oplsaa.ff/forcefield.itp"
#include "oplsaa.ff/spce.itp"
#include "oplsaa.ff/ions.itp"

[ system ]
; Name
water

[ molecules ]
; Compound             #mols

SOL         7396

修改mdp文件

这里需要根据gro文件中仿真盒子尺寸的大小修改rlist,rcoulomb,rvdw这几个值,最好是以1A°为差一直尝试到gromacs不再提示调整rlist。还有约束键的类型要是H-bond,具体值如下:

em.mdp

; Run Control
integrator	= steep		; A steepest descent algorithm for energy minimization
emtol 		= 100		; Stop minimization when the maximum force < 100.0 kJ/mol/nm
emstep 		= 0.01		; [nm] initial step-size
nsteps 		= 50000		; Maximum number of (minimization) steps to perform

; Output Control
nstlog 		= 1		; # of steps that else between writing energies to the log file
nstenergy 	= 1		; # of steps that else between writing energies to energy file

; Neighbor searching and nonbonded interaction
cutoff-scheme	= verlet	; Cutoff-related parameters
nstlist		= 1		; Frequency to update the neighbor list and the long-range forces
ns_type		= grid		; Method to determine neighbor list (simple, grid)
pbc		= xyz		; Use periodic boundary conditions in all directions
rlist		= 0.9		; [nm] Cut-off for making neighbor list (short range forces)

; Electrostatics
coulombtype	= PME		; Treatment of long range electrostatic interactions
rcoulomb	= 0.9		; [nm] Short-range electrostatic cut-off

; Van der waals
vdwtype		= cutoff
rvdw		= 0.9		; [nm] Short-range Van der Waals cut-off

; Temperature and pressure coupling are off during EM
tcoupl          = no		; No temperature coupling
pcoupl          = no		; No pressure coupling (fixed box size)

; No velocities during EM 
gen_vel          = no		; Do not generate velocities

; options for bonds
constraints      = none		; No contraints except those defined in the topology

nvt.mdp

; Run Control
integrator	= md		; leap-frog integrator of Newton's equations of motion
tinit		= 0		; [ps] starting time for the run
dt		= 0.002		; 2 fs - [ps] time step for integration
nsteps		= 50000		; 100 ps

; Output Control
nstxout		= 500		; save coordinates every 1 ps
nstenergy	= 500		; save energies every 1 ps
nstlog		= 500		; update log file every 1 ps

; Bond parameters
continuation    = no		; apply constraints to the start configuration and reset shells
constraints     = h-bonds	; convert all bonds to constraints
constraint-algorithm = lincs	; holonomic constraints 

; Neighbor searching and nonbonded interactions
cutoff-scheme	= verlet	; Cutoff-related parameters
nstlist		= 20		; Frequency to update the neighbor list and the long-range forces
ns_type		= grid		; Method to determine neighbor list (simple, grid)
pbc		= xyz		; Use periodic boundary conditions in all directions
rlist		= 0.8		; [nm] Cut-off for making neighbor list (short range forces)

; Electrostatics
coulombtype	= PME		; Treatment of long range electrostatic interactions
rcoulomb	= 0.8		; [nm] Short-range electrostatic cut-off

; Van der waals
vdwtype		= cutoff
rvdw		= 0.8		; [nm] Short-range Van der Waals cut-off

;                    Do I need pme-order and fourier spacing?

;pme_order	= 4		; cubic interpolation
;fourierspacing	= 0.16		; grid spacing for FFT

; Temperature coupling
tcoupl		= v-rescale	; Temperature coupling, modified Berendsen thermostat‬
tc_grps         = system	; Group to couple to separate temperature baths
tau_t           = 0.1		; [ps] Time constant for coupling
ref_t           = 300		; [K] Reference temperature for coupling

; Pressure coupling is off for NVT
Pcoupl          = No		; No pressure coupling (fixed box size)

; Generate velocities to start
gen_vel		= yes		; assign velocities from Maxwell distribution
gen_temp	= 300		; temperature for Maxwell distribution
gen_seed	= -1		; generate a random seed

npt.mdp

; Run Control
integrator	= md		; leap-frog integrator of Newton's equations of motion
tinit		= 0		; [ps] starting time for the run
dt		= 0.002		; 2 fs - [ps] time step for integration
nsteps		= 500000	; 1000 ps

; Output Control
nstxout		= 5000		; save coordinates every 10 ps
nstenergy	= 500		; save energies every 1 ps
nstlog		= 50000		; update log file every 100 ps

; Bond parameters
constraints     = h-bonds	; convert all bonds to constraints
constraint-algorithm = lincs	; holonomic constraints 

; Neighbor searching and nonbonded interactions
cutoff-scheme	= verlet	; Cutoff-related parameters
nstlist		= 20		; Frequency to update the neighbor list and the long-range forces
ns_type		= grid		; Method to determine neighbor list (simple, grid)
pbc		= xyz		; Use periodic boundary conditions in all directions
rlist		= 0.9		; [nm] Cut-off for making neighbor list (short range forces)

; Electrostatics
coulombtype	= PME		; Treatment of long range electrostatic interactions
rcoulomb	= 0.9		; [nm] Short-range electrostatic cut-off

; Van der waals
vdwtype		= cutoff
rvdw		= 0.9		; [nm] Short-range Van der Waals cut-off

; Temperature coupling
tcoupl		= v-rescale	; Temperature coupling, modified Berendsen thermostat
tc_grps         = system	; Group to couple to separate temperature baths
tau_t           = 0.1		; [ps] Time constant for coupling
ref_t           = 300		; [K] Reference temperature for coupling

; Pressure coupling is off for NVT
Pcoupl          = berendsen	; Pressure coupling on in NPT
pcoupltype	= semiisotropic
tau_p           = 1.0		; [ps] time constant for coupling
ref_p           = 1.0 1.0	; [bar] reference pressure for coupling
compressibility = 0.0 4.5e-05 ; [bar^-1] isothermal compressibility of water

注意npt.mdp这个文件的这段和一般的模板不一样的。

; Pressure coupling is off for NVT
Pcoupl          = berendsen	; Pressure coupling on in NPT
pcoupltype	= semiisotropic
tau_p           = 1.0		; [ps] time constant for coupling
ref_p           = 1.0 1.0	; [bar] reference pressure for coupling
compressibility = 0.0 4.5e-05 ; [bar^-1] isothermal compressibility of water

能量最小化

gmx grompp -f em.mdp -c water_test.gro -p water_test.top -o water_test_em.tpr
gmx mdrun -v -deffnm water_test_em

NVT平衡

gmx grompp -f nvt.mdp -c water_test_em.gro -p water_test.top -o water_test_nvt.tpr
gmx mdrun -v -deffnm water_test_nvt

NPT平衡

gmx grompp -f npt.mdp -c water_test_nvt.gro -p water_test.top -o water_test_npt.tpr
gmx mdrun -v -deffnm water_test_npt

测试静止介电常数

gmx dipoles -corr total -c -f water_test_npt.trr -s  water_test_npt.tpr

这个输出的数应该是50,80这样的,和仿真盒子的大小有关系,但是一定不是1,建议做完每步看一下生成的gro文件,如果gromacs提示increase box size or decrease rlist, 不要增加盒子的大小,这样会使得一些分子跑到扩展盒子的界面上去,这样会严重影响分子的性质。在每步生成gro中如果出现分子分离了应该立刻停止更改结构或者调整mdp参数。

上一篇:amber模拟kcl水溶液


下一篇:UE4.26 水插件学习笔记