制作可用于计算的介电常数的水分子层
总的来说，需要用溶解的方式制作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参数。