正确重复spce水的宏观介电常数计算

常数计算"/>

正确重复spce水的宏观介电常数计算

制作仿真盒子

 gmx solvate -box 2.5 2.5 5 -cs spc216.gro -o water_39_test.gro

再做top

vi water_39_test.top

#include "oplsaa.ff/forcefield.itp"
#include "oplsaa.ff/spce.itp"
#include "oplsaa.ff/ions.itp"[ system ]
; Name
water[ molecules ]
; Compound             #molsSOL         1040

做模拟

能量最小化

gmx grompp -f minim.mdp -c water_39_test.gro -p water_39_test.top -o water_39_minim.tpr
gmx mdrun -v -deffnm water_39_minim

nvt

gmx grompp -f nvt_38.mdp -c water_39_minim.gro -p water_39_test.top -o water_39_nvt.tpr
gmx mdrun -v -deffnm water_39_nvt

npt

gmx grompp -f npt_38.mdp -c water_39_nvt.gro -p water_39_test.top -o water_39_npt.tpr
gmx mdrun -v -deffnm water_39_npt

计算介电常数

gmx dipoles -corr total -c water_39_npt.xvg -f water_39_npt.trr -s water_39_npt.tpr

这里得出结论：截止半径越大越接近理论值，模拟时间越长越接近理论值
可以不做nvt，这里为了后续工作做了nvt

主要mdp

minim.mdp

; LINES STARTING WITH ';' ARE COMMENTS
title		= Minimization	; Title of run
;define      = -DPOSRES
; Parameters describing what to do, when to stop and what to save
integrator	= steep		; Algorithm (steep = steepest descent minimization)
emtol		= 500.0  	; Stop minimization when the maximum force < 10.0 kJ/mol
emstep      = 0.01      ; Energy step size
nsteps		= 50000	  	; Maximum number of (minimization) steps to perform
energygrps	= system	; Which energy group(s) to write to disk; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist		    = 1		    ; Frequency to update the neighbor list and long range forces
cutoff-scheme   = Verlet
ns_type		    = grid		; Method to determine neighbor list (simple, grid)
rlist		    = 0.8		; Cut-off for making neighbor list (short range forces)
coulombtype	    = PME		; Treatment of long range electrostatic interactions
rcoulomb	    = 0.8		; long range electrostatic cut-off
rvdw		    = 0.8		; long range Van der Waals cut-off
pbc             = xyz 		; Periodic Boundary Conditions

nvt_38.mdp

title       = Protein-ligand complex NVT equilibration  
;define      = -DPOSRES  ; position restrain the protein and ligand 
; Run parameters 
integrator  = md        ; leap-frog integrator 
nsteps      = 100000    ; 1 * 100 = 100 ps 
dt          = 0.001     ; 1 fs ; Output control 
nstxout     = 1000       ; save coordinates every 1.0 ps 
nstvout     = 1000       ; save velocities every 1.0 ps 
nstenergy   = 1000       ; save energies every 1.0 ps 
nstlog      = 1000       ; update log file every 1.0 ps 
energygrps  = system ; Bond parameters 
continuation    = no            ; first dynamics run 
constraint_algorithm = lincs    ; holonomic constraints  
constraints     = all-bonds     ; all bonds (even heavy atom-H bonds) constrained lincs_iter      = 1             ; accuracy of LINCS 
lincs_order     = 4             ; also related to accuracy 
; Neighborsearching 
cutoff-scheme   = Verlet 
ns_type         = grid      ; search neighboring grid cells 
nstlist         = 10        ; 20 fs, largely irrelevant with Verlet 
rcoulomb        = 1       ; short-range electrostatic cutoff (in nm) 
rvdw            = 1       ; 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 
tcoupl      = V-rescale                     ; modified Berendsen thermostat 
tc-grps     = system    ; two coupling groups - more accurate 
tau_t       = 0.1                        ; time constant, in ps 
ref_t       = 298                        ; reference temperature, one for each 
group, in K ; Pressure coupling 
pcoupl      = no        ; no pressure coupling in NVT 
; Periodic boundary conditions 
pbc         = xyz       ; 3-D PBC ; Dispersion correction 
DispCorr    = EnerPres  ; account for cut-off vdW scheme ; Velocity generation 
gen_vel     = yes       ; assign velocities from Maxwell distribution 
gen_temp    = 298       ; temperature for Maxwell distribution 
gen_seed    = -1        ; generate a random seed

npt_38.mdp

title       = Protein-ligand complex NPT equilibration  
;define      = -DPOSRES -DPOSRES_LIG  ; position restrain the protein and ligand 
; Run parameters 
integrator  = md        ; leap-frog integrator 
nsteps      = 2000000   ; 1.9ns 
dt          = 0.001     ; 1 fs ; Output control 
nstxout     = 1000       ; save coordinates every 1.0 ps 
nstvout     = 1000       ; save velocities every 1.0 ps 
nstenergy   = 1000       ; save energies every 1.0 ps 
nstlog      = 1000       ; update log file every 1.0 ps 
energygrps  = 
; Bond parameters continuation    = yes           ; first dynamics run constraint_algorithm = lincs    ; holonomic constraints  
constraints     = all-bonds     ; all bonds (even heavy atom-H bonds) constrained lincs_iter      = 1             ; accuracy of LINCS 
lincs_order     = 4             ; also related to accuracy 
; Neighborsearching 
cutoff-scheme   = Verlet 
ns_type         = grid      ; search neighboring grid cells 
nstlist         = 10        ; 20 fs, largely irrelevant with Verlet 
rcoulomb        = 1.5       ; short-range electrostatic cutoff (in nm) 
rvdw            = 1.5       ; 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 
tcoupl      = V-rescale                     ; modified Berendsen thermostat 
tc-grps     = system    ; two coupling groups - more accurate 
tau_t       = 0.1                        ; time constant, in ps 
ref_t       = 298                        ; reference temperature, one for each 
group, in K 
; Pressure coupling 
pcoupl      = Parrinello-Rahman             ; pressure coupling is on for NPT pcoupltype  = isotropic                     ; uniform scaling of box vectors 
tau_p       = 2.0                           ; time constant, in ps 
ref_p       = 1.0                           ; reference pressure, in bar compressibility = 4.5e-5                    ; isothermal compressibility of water, bar^-1 
refcoord_scaling    = com 
; 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 off after NVT