There is a good QM/MM tutorial which can be found here. However, the tutorial is a bit complex and it's difficult to reproduce the whole tutorial for some missing details. Here, we use a simple system to help beginners to setup a QM/MM simulation step by step.


    Before starting this tutorial, we assume that:
  • The GROMACS Software and the Gaussian package has been installed properly (the version of GROMACS used in this tutorial is 3.3.1).
  • The computer system used is Linux Operator System.
  • You have download the structural file of the molecular used here. Click peptide.pdb to download it.
  • You have the fundamental knowledge on GROMACS and Linux System.

Setup up the system in vacuo

    Step 1. Create the gro and top files

    Firstly, we create a new directory and save the peptide.pdb under this directory.

  • Type command: pdb2gmx -f peptide.pdb -p peptide -o peptide -ter
  • For force field, select 6 (OPLS force field)
    For N-terminus type, select 3 (None)
    For C-terminus type, select 3 (None)

    Then, we obtain the three files: the structural file peptide.gro, the topology file peptide.top and the position restriction file porse.itp.

    Step 2. Adding link atoms

    Here, we don't explain why we should introduce link atoms in QM/MM simulation, you should find the answer by yourself and it isn't a tough work.
    Before QM/MM simulation, we should modify the peptide.gro and peptide.top files, and specify which atoms should be included in QM region and specify the QMMM parameters in mdp file.

    The following figure show how we divide the system into QM part and MM part.

    Modify the peptide.gro and peptide.top files
    As illustrated in the figure above, we add two link atoms (LA) into the molecular, the atom number should be changed from 26 to 28. We should add the coordinates of these two atoms into the structural file peptide.gro. How to determine the coordinates of the link atoms? For the link atoms added between CA atom and C atom of GLY2 (The names of the atoms refer to the names in peptide.gro file), the link atoms should lie in the bond formated by CA atom and C atom of GLY2 and the bond length between LA atom and C atom should be proportional to the bond length between CA atom and C atom of GLY2. For the C-C bondlength is 0.153nm and the C-H bondlength is 0.108nm, we obtain the ratio is 0.108/0.153=0.706 (also see the "adding link atoms" part of this link). According to the coordinates of CA atom and C atom in peptide.gro, one can easily obtain the coordinates of the link atom. So we add the following content at the end of the structural file:

    5GLY LA 27 1.458 1.579 1.544
    6GLY LA 28 1.832 1.579 1.695

    We save the modified structural file as qmmm.gro. The group name of link atoms are not important, you can use other names instead of "5GLY" or "6GLY".

    We have modified the structure of the molecular, we should modify the topology file, too. Firstly, we add the following content at the end of the "atoms" part of the topology file:

    27 opls_997 5 GLY LA 9 0.00 0.000 ;
    28 opls_997 5 GLY LA 9 0.00 0.000 ;

    The "opls_997" is the atom type of link atom defined by ourselves in the ffoplsaanb.itp file (one of the OPLS force filed files). We added the following content into the ffoplsaanb.itp file:

    opls_997 LA 1 0.00000 0.000 A 0.00000e+00 0.00000e+00

    The modified ffoplsnb.itp is available here. We also add the following content into the ffoplsaa.atp file:

    opls_997 0.0 ; LA atoms in QMMM

    The modified ffoplsaa.atp is available here. Other related OPLS force field files are ffoplsaabon.itp and ffoplsaa.itp which are not modified at all. All the abovementioned OPLS force field files should be saved under the current directory.

    Other modifications are listed as following:

  • For the "bonds" part, we need to change the bond type of those bonds formed between two QM atoms into 5.
  • For the "angles" part, we should remove those angles within the QM part. We define an angle as an angle within the QM part if at least two of the three atoms which form this angle are QM atoms.
  • For the "dihedrals" part (improper or proper dihedrals), we should remove those dihedrals within the QM part. We define a dihedral as a dihedral within the QM part if at least three of the four atoms which form this dihedral are QM atoms.
  • Add the "dummies2" part and the "constraints" part into the topology file:
    [ dummies2 ]
    27 12 9 1 0.706
    28 16 19 1 0.706

    [ constraints ]
    9 12 2 0.153
    16 19 2 0.153

  • For details, also see this link. The modified topology file is saved as qmmm.top

    Specify which atoms should be included in QM region
  • Type command: make_ndx -f qmmm.gro
  • Tyep: q
    Using make_ndx, we generate a new index file (index.ndx) for the qmmm.gro, and add the following part:

    1 2 3 4 5 6 7 8 9 10 11 19 20 21
    22 23 24 25 26
    12 13 14 15 16 17 18 27 28

    to specify the atoms included in QM part as "QMatoms" and the atoms included in MM part as "MMatoms". The modified index.ndx is available here.

    Specify the QMMM parameters in mdp file
    To do QMMM simulation, we specify the QMMM parameters as following:
    QMMM = yes
    QMMM-grps = QMatoms
    QMmethod = B3LYP
    QMbasis = 6-31G*
    QMMMscheme = normal
    QMcharge = 0
    QMmult = 1
    The options of the QMMM parameters depends on the system. The mdp file we used is qmmm.mdp.

    Step 3. Generate tpr file

    Till now, we have prepared all the files for generating the tpr file for simulation.

  • Type command: grompp -f qmmm.mdp -p qmmm.top -n index.ndx -c qmmm.gro -o peptide.tpr
  • We obtain the peptide.tpr file. Using this tpr file, we can start QM/MM simulation. Here, we will find two warnings, because the incompleteness of the default force field. Since it's just tutorial to show how to setup a QMMM simulation, we ignored it. However, to setup a meaningful simulation, one need to learn more about molecular simulations, which is out of the scape of this tutorial.

Copyright©Wenjin Li. All Rights Reserved
Updated: 2011/4/11, By Wenjin Li