(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
A Coupled Potential QM/MM Simulation
By Ross Walker
Updated for AMBER
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:
Initial min of our structure QMMM
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
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:
|300K constant temp QMMMMD
cut=8.0, ntc=2, ntf=2,
$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.
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.
(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
Copyright Ross Walker 2015