Run QM/MM calculations with CPMD

This tutorial explains how to run QM/MM calculations with CPMD starting from Amber system files (topology and coordinates).

Conversion from Amber to Gromos format

CPMD is making use of files in the Gromos format. We therefore have to convert the Amber files (system.top and system.rst) to the Gromos format.

The CPMD package comes with scripts that can be used to do that. They can be found in the following directory:

$CPMD_ROOT/modules/MM_Interface/Util_QMMM/

In particular we need the amber2gromos.x and estimate_gromos_size scripts. As first step:

ambertogromos.x system.top system.rst

This operation should produce three files:

  • gromos.crd
  • gromos.inp
  • gromos.top

Some modifications are required.

In gromos.inp, the sections SYSTEM, SUBMOLECULES and FORCES, need to be checked and maybe modified according to the solute/solvent and QM/MM separation. Here everything is not clear yet, especially when using non rigid solvent!

In the section ATOMTYPENAME of gromos.top, replace the names of the types of the atoms from the standard generic force field library gaff:

$AMBERHOME/dat/leap/parm/gaff.dat

to atom names from the Amber force field library:

$AMBERHOME/dat/leap/parm/parm99.dat

Create a CPMD input file for a QM/MM run

We now need to create an input file for CPMD. You can download a template here.

This input will perform a standard NVE Born-Oppenheimer molecular dynamics (BOMD). Places with XXX need to be updated. The rest of the file needs to be very carefully check and changed according to your system and the type of simulation you want to run.

First we need to update the array size at the end if the &QMMM section. for that run the following command and paste the results:

estimate_gromos_size gromos.top

We also need to determine the dimension of the cell. For that you can use the cell_size script from the comp_chem_py library:

cell_size XXX YYY

where XXX and YYY are the residue names of the QM part.

Alternatively, you can run the following commands:

grep XXX gromos.crd > QM.coor # XXX is the residue name of the QM part
sort -k 5 -n QM.coor # order x coordinates (-k 5)
echo "(D + 0.7)*10/0.529" | bc -l # where D is the max difference in x coordinates

The result should be taken as the dimension of the box in the x axis. repeat with -k 6 and -k 7 for the y and z axis. You can then take the maximum value of the three axis and paste it in place of the XXX after the CELL keyword in the &SYSTEM section of CPMD input.

Finally, in the &ATOMS section, add the number of QM atoms of a given type under the value for LMAX and on the next line all the indices of the QM atoms as given in the gromos.crd file.

Then run CPMD as usual. Just make sure CPMD has been compiled with the -qmmm option (see compile cpmd).