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).