(Note: These tutorials are meant to provide illustrative examples of how to use the AMBER software suite to carry out simulations that can be run on a simple workstation in a reasonable period of time. They do not necessarily provide the optimal choice of parameters or methods for the particular application area.)
Copyright Ross Walker 2015

A Simple Coupled Potential QM/MM Simulation

SECTION 3

A Coupled Potential QM/MM Simulation

By Ross Walker

Updated for AMBER 15

3) QM/MM MD

Now, let's try repeating the simulation using a coupled QM/MM potential. We will model the NMA molecule using the semi-empirical PM3 Hamiltonian while we model water classically. In this situation we have no bonds crossing the QM/MM boundary so do not need to worry about where sander will place hydrogen link atoms. There will be none. However, if we did have bonds that crossed the boundary sander would add the link atoms completely automatically. See the manual for more information.

For QM/MM MD the input files look very similar. There are a few notable differences, however, related to the specification of the atoms to be modeled quantum mechanically. Here is the minimization input file we will use:

The extra variables in this file are:

ifqnt = 1 - This is the flag that tells sander that we want a QM/MM run. It will then look for a &qmmm namelist.

Most of the options I have specified in the &qmmm namelist are the defaults but it is often good to include them in the input file to make it easier to work out what you did at the later date:

qmmask = ':1-2' - This specifies what to treat quantum mechanically using standard AMBER mask notation. Here it implies residue 1 and 2 which is the ACE and NME making up the NMA molecule. Alternatively you can manually specify the atom numbers using iqmatoms. E.g. iqmatoms=1,2,3,4,5,6,7,8,9,10,11,12,
qmcharge = 0   - The integer charge of the QM region (default = 0)
qm_theory = PM3   - Use the PM3 Hamiltonian (default = PM3)
qmshake = 1   - Shake QM hydrogen atoms (default = 1 if ntc=2)
qm_ewald = 1   - Use an Ewald type treatment for long range electrostatics (default = 1 if ntb>0)
qm_pme = 1   - Use an Particle Mesh Ewald method as Ewald type (default = 1 if qm_ewald=1 and use_pme=1)

For AMBER this is all we need:

$AMBERHOME/bin/sander -O -i min_qmmm.in -o min_qmmm.out -p NMA.prmtop -c NMA.inpcrd -r min_qmmm.rst

Input Files: min_qmmm.in, NMA.prmtop, NMA.inpcrd
Output Files: min_qmmm.out, min_qmmm.rst

As we did in the classical simulation above we will run 1,000fs of simulation. Again our input file is the same apart from the QM specifications. Here is the MD input file:

$AMBERHOME/bin/sander -O -i md_qmmm.in -o md_qmmm.out -p NMA.prmtop -c min_qmmm.rst -r md_qmmm.rst -x md_qmmm.mdcrd

You can track the output in another terminal window with "tail -f md_qmmm.out" if you wish. The first thing you should notice with this simulation is that each step takes slightly longer than the classical simulation. The reason for this is that we are doing a complete SCF at every step. However, the difference here should not be huge. In going to AMBER 9 the performance of the QM region of the code was increased more that 10 fold such that for systems up to about 50 QM atoms, the cost of doing QM/MM is minimal. This should take just over a minute to run.

Input Files: md_qmmm.in, NMA.prmtop, min_qmmm.rst
Output Files: md_qmmm.out, md_qmmm.rst, md_qmmm.mdcrd

Take a quick look at the md_qmmm.mdcrd file in VMD to check everything went okay. Remember to gunzip it first.

The final step is to compare the two simulations.

---

CLICK HERE TO GO TO SECTION 4

---

(Note: These tutorials are meant to provide illustrative examples of how to use the AMBER software suite to carry out simulations that can be run on a simple workstation in a reasonable period of time. They do not necessarily provide the optimal choice of parameters or methods for the particular application area.)
Copyright Ross Walker 2015