Biomodeling: water box

Hello, water box! 1_water_box in the tutorial files

Input files for MD simulation topology file

grompp.mdp: MD Option file

The first MD simulation

topol.tpr: Input file for MDRUN

The first MDRUN

Minimize the system

Try MD again

MD for 10 ns

check box size and temperature

Structural properties of water: Radial distribution function, g(r)

Calculation of g(r) using Gromacs

set xlabel "r (nm)"

set ylabel "g(r)"

set grid

pl [0:1] 'rdf.xvg' u 1:2 w l

The relationship between g(r) and the free energy G(r)

set xlabel "r (nm)"

set ylabel "Energy (kT)"

set grid

pl [0:1][-1:2] 'rdf.xvg' u 1:(-log($2)) w l

Calculation of g(r) using VMD

Potential of mean force by Umbrella sampling

Umbrella sampling

gmx make_ndx -f topol.tpr -o pull.ndx

ri 1

ri 2


Perform simulation in each window

pull-ngroups             = 2

pull-ncoords             = 1

pull-group1-name         = r_1

pull-group2-name         = r_2

pull-coord1-type         = umbrella

pull-coord1-geometry     = distance

pull-coord1-groups       = 1 2

pull-coord1-dim          = Y Y Y

pull-coord1-init         = 0.2              ; the equilibrium length of spring is 0.2 nm

pull-coord1-k            = 1000             ; spring constant in kJ/mol nm

for i in 03 04 05 06 07 08 09 10


gmx grompp -c 2.pdb -f pull.$i.mdp -n pull.ndx -o pull.$i.tpr

gmx mdrun -s pull.$i.tpr -nt 4 -gpu_id 0 -cpi -noappend -deffnm pull.$i -nsteps 500000


Weighted Histogram Analysis Method (WHAM)

gmx wham

set xlabel "r (nm)"

set ylabel "Probability"

set grid

pl 'histo.xvg' u 1:2 w l title "Window r = 0.2 nm" ,\

   'histo.xvg' u 1:3 w l title "Window r = 0.3 nm" ,\

   'histo.xvg' u 1:4 w l title "Window r = 0.4 nm" ,\

   'histo.xvg' u 1:5 w l title "Window r = 0.5 nm" ,\

   'histo.xvg' u 1:6 w l title "Window r = 0.6 nm" ,\

   'histo.xvg' u 1:7 w l title "Window r = 0.7 nm" ,\

   'histo.xvg' u 1:8 w l title "Window r = 0.8 nm" ,\

   'histo.xvg' u 1:9 w l title "Window r = 0.9 nm" ,\

   'histo.xvg' u 1:10 w l title "Window r = 1.0 nm"

set xlabel "Distance r (nm)"

set ylabel "Energy (kJ/mol nm)"

set grid

pl [0:1] 'rdf.xvg' u 1:(-2.5 * log($2)) w l title "-kT ln g(r)", 'profile.xvg' u 1:($2+9) w l title "Delta G from WHAM", 'profile.xvg' u 1:($2 + 2.5*2*log($1) +8.5) w l title "Delta G corrected by 2r ln(r)"


Dynamic properties of water: Diffusion


What is the pairwise interaction free energy between a pair of Na and Cl ions?

[ system]

SOL                  NNNN        ; the number of SOL must be decreased by 2

Na                     1

Cl                     1