[gmx-users] Hydrogen bonding energy calculation (Ben)

[leoneverlasting]

Hi Dr. Spoel, 

Thanks for your reply! I set my script as follows:

gmx pdb2gmx -f water.pdb -o water.gro -water tip3p << EOF #water.pdb is a pdb file contains a single water molecule
5		#OPLS
2		#TIP3P
EOF
gmx editconf -f water.gro -o water_newbox.gro -c -d 1.5 -bt cubic
gmx solvate -cp water_newbox.gro -cs spc216.gro -o water_solv.gro -p topol.top
gmx grompp -f em.mdp -c water_solv.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em
gmx grompp -f nvt.mdp -c em.gro -p topol.top -o nvt.tpr
gmx mdrun -deffnm nvt
gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
gmx mdrun -deffnm npt
gmx grompp -f md1.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_1.tpr  #a short 200 ps run
gmx mdrun -deffnm md_0_1
gmx grompp -f md2.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_2.tpr  #a longer 2 ns run
gmx mdrun -deffnm md_0_2
gmx trjconv -s md_0_2.tpr -f md_0_2.trr -o md_0_2_noPBC.trr -pbc mol -ur compact << EOF
2           #SOL (water)
EOF
gmx hbond -f md_0_2_noPBC.trr -s md_0_2.tpr -num hbnum.xvg -temp 298.15 -ac hbac.xvg << EOF
2	    #SOL (water)
2           #SOL (water)
EOF

Here is the specification of my md*.mdp file

; Run parameter
integrator          =  md               ; leap-frog integrator
dt                  =  0.002    ; ps !
nsteps              =  1000000  ; total 2 ns
; Output control
nstxout             =  50       ; save coordinates every 100 fs
nstvout             =  0
nstenergy           =  50       ; save energy every 100 fs
nstlog              =  50
; Bond parameters
continuation        = yes
constraint_algorithm   = SHAKE
shake-tol           =  0.0001
constraints         =  all-bonds
; Neighborsearching
cutoff-scheme   = Verlet
ns_type             =  grid
nstlist             =  5
rcoulomb            =  1.0
rvdw                =  1.0
; Electrostatics
coulombtype         =  PME
pme_order           = 4
fourierspacing      = 0.16
; Temperature coupling is on
Tcoupl              =  berendsen
tc-grps             =  SOL
tau_t               =  0.1
ref_t               =  298.15
; Isotropic pressure coupling is now on
Pcoupl              =  berendsen
Pcoupltype          =  isotropic
compressibility     =  4.5e-5
tau_p               =  1.0
ref_p               =  1.0
refcoord_scaling    = com
; Periodic boundary conditions
pbc                 = xyz
; Dispersion correction
DispCorr            = EnerPres
gen_vel             =  no

Could you take a look at my parameters and script and let me know if there is anything that might cause the deviation? Thanks!!

Ben