Biomodeling: water box

• Generate a water box in a 5-nm cube

  • gmx solvate -box 5 5 5 -cs -o 1.pdb

    • -box: box dimension

    • -cs: type of water. By default, a 3-point water type.

    • -o: output file name

  • Open the output file using a text editor and see the format of PDB file

  • How many water molecules do we have in this 5-nm cubic box?

    • Log message will tell the number.

    • Otherwise, we can count by using wc command

      • grep OW 1.pdb | wc -l

  • Use VMD to see the water box.

topol.top: topology file

•  topol. top file contains information of force field and the simulation system.

• copy files/topol.top to PWD (present working directory)

  • cp files/topol.top .

    • "." indicates the current directory.

• A reusable part can be saved in a separate file, which can included in topol.top 

  • #include REUSABLE_FILE

• To see the installed force fields, see $GMXDATA/top directory

  • ls $GMXDATA/top

• In this example, we will use amber99sb.ff.

  • amber99sb.ff/forcefield.itp contains files with all force field parameters.

    • See the files that forcefield.itp includes.

  • amber99sb.ff/tip3p.itp contains "moleculetype" definition for the TIP3P water model.

grompp.mdp: MD Option file

• MDP file contains MD simulation Options.

  • See https://manual.gromacs.org/current/user-guide/mdp-options.html for details.

• Copy a sample MDP file to PWD (present working directory)

  • cp files/*.mdp .

• Open grompp.mdp using a text editor

  • See each section of the file.

topol.tpr: Input file for MDRUN

• Now, we have all input files for simulation:

  • Structure file: 1.pdb

  • topol.top

  • grompp.mdp

• gmx grompp program assembles information from all three files into a topol.tpr file

  • gmx grompp -f grompp.mdp -c 1.pdb -p topol.top -pp processed.top  -o topol.tpr

    • -pp: write the topology file interpreted by grompp. It writes a long file with all the parameters in a single file without using #include.

    • 
      

  • topol.tpr is a binary file. 

    • To see the content of topol.tpr use gmx dump command.

      • gmx dump -s topol.tpr

The first MDRUN

• Now, we perform MD by gmx mdrun program

  • gmx mdrun -nt 4 -gpu_id 0

    • -nt: the number of CPU cores for this run.

    • -gpu_id: GPU ID for this run if you have GPU in your system, 0, 1, 2, etc. If you have only one GPU in the system, use 0 or omit this option.

• The simulation will fail because LJ interaction forces are too high (i.e., too much acceleration).

  • see md.log file

    • It is important to check the energy values

    • In normal MD simulations in a water box, the total energy of the system should be a huge negative value.

• Prepare a new topol.tpr file for minimization:

  • gmx grompp -f mini.mdp -c 1.pdb

    • Compare mini.mdp with grompp.mdp:

      • diff mini.mdp grompp.mdp

  • gmx mdrun -nt 4 -gpu_id 0 -c 2.pdb

    • -c: The minimized system will be written t0 2.pdb

  • see md.log and ener.edr files to see change in energy.

    • gmx energy -f ener.edr

• Open grompp.mdp and change the settings:

  • nsteps                   = 2500  ; Do MD for 2500 steps

  • nstxout-compressed       = 1     ; write trajectory every step.

• Generate topol.tpr

  • gmx grompp -c 2.pdb -f grompp.mdp

    • -c: Now, we are using the minimized structure.

• Perform MD.

  • gmx mdrun -nt 4 -gpu_id 0

    •  it will generate some output files:

      • md.log : text log file

      • ener.edr : binary energy file

      • traj_comp.xtc : trajectory file

      • state.cpt : check point file for continuation

• See traj_comp.xtc using VMD

  • traj_comp.xtc will look weird because water molecules at the boundary are broken.

  • This happens because gmx mdrun wants to have all atoms in a unit cell.

• We can fix this by using gmx trjconv

  • gmx trjconv -f traj_comp.xtc -pbc mol -ur compact -o tmp1.xtc

  • Load tmp1.xtc to VMD, and see the difference.

• We can fix PDB file using gmx trjconv

  • 2.pdb file has broken water molecules because it is an output from gmx mdrun.

  • gmx trjconv -f 2.pdb -pbc mol -ur compact -o 2.pdb

• To prevent molecules from box crossing, use -pbc nojump option.

  • gmx trjconv -s 1.pdb -f traj_comp.xtc -pbc nojump -o tmp2.xtc

    • the starting conformation should be whole. 1.pdb file in this case.

• Change grompp.mdp:

  • nsteps                   = 5000000  ; 10 ns

  • nstxout-compressed       = 10000    ; write trajectory every 20 ps

• gmx grompp -c 2.pdb -f grompp.mdp

• gmx mdrun -nt 4 -gpu_id 0 -cpi -noappend

check box size and temperature

• gmx traj -ob -f traj_comp.part0002.xtc

  • see box.xvg

• gmx energy -f ener.part0002.edr

  • see energy.xvg

Calculation of g(r) using VMD

• Try

• g(r) tells us the ensemble-averaged interaction for water pairs.

  • This is possible we can sample an ensemble of water configurations sufficiently.

    • In other words, this is now a rare event.

  • But, to calculate interaction ∆G for the hard-to-sample pairs, we need a technique that makes sampling rare events possible.

    • Umbrella sampling.

Umbrella sampling

• First, we need to select two water molecules, for which we want to calculate intermolecular interaction free energy.

  • The chance that these two specific water molecules meet is RARE!

    • Thus, it will take a very long time to sample the rare events sufficiently.

  • In umbrella sampling, we constraint the water-water pair by using a harmonic spring.

    • See PULL section in MDP options.

  • Generate index file for two water molecules

    • We will just select the first and the second water molecules. The selection will not change the results. All water molecules are equal!

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

ri 1

ri 2

q

• • • Caution: The initial distance between chosen water molecules must be smaller than half of the box. Otherwise, it violated PBC.

Perform simulation in each window

• 0.2 nm window

  • See PULL section in pull.02.mdp

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

• • Generate pull.02.tpr for gmx mdrun

    • gmx grompp -c 2.pdb -f pull.02.mdp -n pull.ndx -o pull.02.tpr -t state.cpt

      • -t state.cpt: use information from the check point file.

  • Perform MDRUN for 1 ns

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

      • pullx files contain the length of spring

      • pullf files contain the force on the spring (F = -kx)

  • Repeat same things for other windows

for i in 03 04 05 06 07 08 09 10

do

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

done

Discussions

• How to estimate error?

  • gmx wham has bootstrap options.

  • More straightforward way is getting multiple results (or splitting simulation into multiple blocks) and calculate stardard errors.

• How can we tell if the Umbrella sampling is long enough?

  • For many real problems, very long simulations are necessary.

    • Keep continuing each window.

    • When you have multiple pullx files for a window, you need to combine them into a pullx file. For example, awk command.

    • Calculate gmx wham frequently.

      • If PMF converges (i.e., PMF does not change any more), quit.

      • Otherwise, keep simulating.

• Mean squared displacement (MSD) = 6D∆t

• gmx msd -f traj_comp.part0002.xtc -n rdf.xvg