QM/MM tutorial

QM/MM calculations on thymine dimer repair.

I. Building a model of a solvated double stranded DNA molecule with a dimer lesion

DNA building

B (formely know as Biomer) is a Java-based, on-line biomolecular modeling package. It is useful for generating initial structures of biopolymers and small organic molecules. We will use B to create double-stranded DNA oligomer. You are free to choose any sequence, as long as there is at least one TT repeat on one of the strands. The example files that we provide in this tutorial are based on a A20T20 oligomer.

Please consult the B manual and faq page for help on the model building. Save your DNA molecule in pdb format. The pdb file of the A20T20 oligomer can be found here.

Topology Building

To model the interactions we make use of the Amber99 forcefield. The Amber port for gromacs can be downloaded from Eric Sorin's webpage.

To allow the pdb2gmx tool to recognize the thymine dimer, we need to modify the ffamber99.rtp residue database. We introduce two new thymine residues dtA and dtB, and define additional interbase bonds between the C5 atoms and the C6 atoms. Here is what your modified ffamber99.rtp should look like.

With the pdb2gmx command we create the topology (topol.top) and the configuration (conf.gro):

localhost:~>pdb2gmx -ignh -f Bmodel.pdb

The -ignh flag ignores hydrogenatoms in the input file, and forces pdb2gmx to construct the hydrogens based on the ffamber99.hdb database.

Structure Preparation

We first perform a minimization in vacuum, using the steepest descent algorithm. Check steep.mdp to familiarize yourself with the minimization protocol that we follow.

We perform the minimization by executing grompp and then mdrun:

localhost:~>grompp -f steep.mdp

localhost:~>mdrun -v -c minimized.pdb

After minimization the DNA structure should looks like this.

The next step is to add water and neutralize the system. We use the gromacs tools editconf, genbox, and genion.

We first place the minimized structure in the center of a rectangular box:

localhost:~>editconf -c -d 1.6 -f minimized.pdb -o boxed.pdb

Then, we add water, for which we use the tip4p model

localhost:~>genbox -cs tip4p.gro -cp boxed.pdb -p topol -o solved.pdb

Finally, we neutralize the overall negative charge of -38 by adding 38 sodium ions. We first need to create a tpr file, for which we use neutralize.mdp:

localhost:~>grompp -f neutralize.mdp -c solved.pdb -p topol

With the topol.tpr file, we can now use genion to place the 38 counter-ions at random positions in the box:

localhost:~>genion -np 38 -pname Na -o neutral.pdb -random

We now have to modify the topology by adding 38 sodium ions (NA) and removing 38 tip4p water molecules manually. The final topology file can be found here. In what follows we also need an index file that contain the atomnumbers of the DNA, water and ions.


Because of possible overlap between water molecules and DNA atoms, we have to minimize the model system. We use again the steepest descent algortihm to remove possible strain in the initial structure. The mdp file can be found here.

localhost:~>grompp -f mini_sol.mdp -p neutral.top -c neutral.pdb -n neutral.ndx

localhost:~>mdrun -v -c neutral_mini.pdb

Then, we equilibrate the water and ions for 100 ps, using the parameter in equi_sol.mdp. The heavy atoms of the DNA are restrained during the solvent equilibration.

localhost:~>grompp -f equi_sol.mdp -p neutral.top -c neutral_mini.pdb -n neutral.ndx

localhost:~>mdrun -v -c neutral_sol_equi.pdb

The result of the solvent equilibration on the A20T20 oligomer is available here.

Next:II. Equilibration of the DNA with gromacs
Previous: Introduction

back to top

updated 28/10/08