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