water molecular dynamics

em1.mdp

define                   = -DFLEXIBLE
integrator               = steep
nsteps                   = 1000

nstenergy                = 500
nstlog                   = 500
nstxout-compressed       = 1000

constraint-algorithm     = lincs
constraints              = h-bonds

cutoff-scheme            = Verlet

coulombtype              = PME
rcoulomb                 = 0.9

vdwtype                  = Cut-off
rvdw                     = 0.9
DispCorr                 = EnerPres

em2.mdp

integrator               = steep
nsteps                   = 50000

nstenergy                = 500
nstlog                   = 500
nstxout-compressed       = 1000

cutoff-scheme            = Verlet

coulombtype              = PME
rcoulomb                 = 0.9

vdwtype                  = Cut-off
rvdw                     = 0.9
DispCorr                 = EnerPres

npt.mdp

integrator               = md        
dt                       = 0.002     ; 2 fs
nsteps                   = 1000000    ; 2.0 ns
pbc                      = xyz

nstenergy                = 100
nstxout                  = 100
nstvout                  = 100
nstfout                  = 100
nstlog                   = 1000
nstxout-compressed       = 10000

constraint-algorithm     = lincs
constraints              = h-bonds

cutoff-scheme            = Verlet
verlet-buffer-tolerance  = 1.0e-4

coulombtype              = PME
rcoulomb                 = 0.9

vdwtype                  = Cut-off
rvdw                     = 0.9
DispCorr                 = EnerPres

tcoupl                   = Nose-Hoover
tc-grps                  = System
tau-t                    = 2.0
ref-t                    = 298.15
nhchainlength            = 1

pcoupl                   = Parrinello-Rahman 
tau_p                    = 2.0
compressibility          = 4.46e-5
ref_p                    = 1.0 

;continuation             = yes
gen-vel                  = yes
gen-temp                 = 298.15

nve.mdp

integrator		= md 
dt			= 0.002
nsteps			= 200000    ;400ps
pbc			= xyz

;
; Removing CM Translation and Rotation
;comm_mode		=  Linear 
;nstcomm			=  1000
;
; Output Control 
nstcalcenergy           = 10
nstlog			= 10
nstenergy		= 10
nstxout			= 10
nstvout			= 10
nstfout			= 10
;
; Neighbour Searching
;nstlist			= 10
;ns_type			= grid
;rlist			= 1.4 
;
cutoff-scheme            = Verlet
verlet-buffer-tolerance  = 1.0e-4

coulombtype		= pme
rcoulomb		= 0.9

vdwtype            	= Cut-off
rvdw                	= 0.9
DispCorr                 = EnerPres


; Constraints
constraints		= h-bonds 
constraint-algorithm	= lincs

tcoupl			= no 
pcoupl			= no

; Initial Velocities
gen_vel			= yes
gen_temp		= 298.15
gen_seed		= -1

nvt.mdp

integrator               = md        
dt                       = 0.002     ; 2 fs
nsteps                   = 1000000     ; 2.0 ns
pbc                      = xyz

nstenergy                = 100
nstlog                   = 100
nstxout                  = 100
nstvout                  = 100
nstfout                  = 100
nstxout-compressed       = 10000

constraint-algorithm     = lincs
constraints              = h-bonds

cutoff-scheme            = Verlet
verlet-buffer-tolerance  = 1.0e-4

coulombtype              = PME
rcoulomb                 = 0.9

vdwtype                  = Cut-off
rvdw                     = 0.9
DispCorr                 = EnerPres

tcoupl                   = Nose-Hoover
tc-grps                  = System
tau-t                    = 2.0
ref-t                    = 298.15
nhchainlength            = 1

continuation             = yes
;gen-vel                  = yes
;gen-temp                 = 298.15

topol.top

#include "oplsaa.ff/forcefield.itp"
#include "oplsaa.ff/tip4pew.itp"

[ System ]

[ Molecules ]

water_molecular_dynamics.md

Determining the Bulk Viscosity of Rigid Water Models

http://pubs.acs.org/doi/pdf/10.1021/jp211952y

の追試用ファイルたち(for Gromacs 5.1.2)。

#-- We will get 250 molecules
gmx solvate -cs tip4p -o water.gro -box 2.1 2.0 2.0 -p conf/topol.top

#-- Energy minimization
gmx grompp -f conf/em1.mdp -o em1 -c water.gro -p conf/topol.top
gmx mdrun -s em1.tpr -deffnm em1 -v 
gmx grompp -f conf/em2.mdp -o em2 -pp em2 -po em2 -c em1 -t em1 -p conf/topol.top
gmx mdrun -s em2.tpr -deffnm em2

##gmx energy -f em1.edr -o min-energy.xvg
##gmx energy -f em2.edr -o min2-energy.xvg 

#-- NPT ensemble (to determine the density of liquid water)
gmx grompp -f conf/npt.mdp -o nptin -c em2.gro -p conf/topol.top -maxwarn 1
gmx mdrun -s nptin.tpr -deffnm nptout

#-- check density ,pressure, energy, temparature
echo -e "7\n8\n10\n14\n15\n" | gmx energy -f nptout.edr -o nptout.xvg

#-- NVT ensemble (to determine the average energy)
gmx grompp -f conf/nvt.mdp -o nvtin -c nptout -t nptout -p conf/topol.top
gmx mdrun -s nvtin.tpr -deffnm nvtout

#-- check energy
echo -e "7\n9\n11\n" | gmx energy -f nvtout.edr -o nvtout.xvg

#-- NVE ensemble
for i in `seq 1 100`;do
  gmx grompp -f conf/nve.mdp -o nvein -c nvtout -t nvtout -p conf/topol.top
  gmx mdrun -s nvein.tpr -deffnm nveout_$i
  echo -e "7\n8\n10\n" | gmx energy -f nveout_$i.edr -o nveout_$i.xvg    #--energy,temperature,pressure
done