(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 2004

TUTORIAL 6 - AMBER 9 VERSION - SECTION 3

A Coupled Potential QM/MM Simulation

By Ross Walker

PLEASE NOTE THIS TUTORIAL IS FOR AMBER 9 ONLY.

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 semiempirical PM3 Hamiltonian while all the water we will model classically. In this situation we have no bonds crossing the QM/MM boundary and 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 info.

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 QMMM 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)
qmtheory = 1   - Use the PM3 Hamiltonian (default = 1)
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)

Unlike AMBER 8 where we needed extra input files and a special executable for AMBER 9 this is all we need.

$AMBERHOME/exe/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/exe/sander -O -i md_qmmm.in -o md_qmmm.out -p NMA.prmtop -c NMA_min_qmmm.crd -r NMA_md_qmmm.rst -x NMA_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 has been 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 around 3 minutes or so to run.