(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 2
A Coupled Potential QM/MM Simulation
By Ross Walker
Updated for AMBER 15
2) Classical MD
Now that we have our topology and coordinate files, we are ready to carry out classical molecular dynamics. First things first though we need to minimize our system to remove any bad contacts created by solvation and also to allow our artificially created NMA structure to relax. For this we will run 500 steps of minimization. Our input file is very similar to files we have used in the previous tutorials. Since AMBER supports periodic boundaries with PME for both classical and QM/MM MD we will use it (NTB=1). AMBER also supports shake for QM atoms and since we are not interested in any kind of reaction involving hydrogen in the QM/MM simulation we will use shake (NTC=2, NTF=2) for this simulation to allow us to use a 2 fs time step. Since this is a periodic simulation with PME we can safely use a cut off of 8 angstroms. Here is our input file:
Initial minimisation of our structure &cntrl imin=1, maxcyc=500, ncyc=200, cut=8.0, ntb=1, ntc=2, ntf=2 / |
|
$AMBERHOME/bin/sander -O -i min_classical.in -o min_classical.out -p NMA.prmtop -c NMA.inpcrd -r min_classical.rst
Input Files: min_classical.in,
NMA.prmtop, NMA.inpcrd
Output Files: min_classical.out,
min_classical.rst
Once we have minimized the system we will run 1,000 fs of molecular dynamics at a constant temperature of 300K. Normally we would slowly heat our system to 300K but since this is just a small NMA molecule in solution it should be fairly stable and should survive being started out at 300K. During the MD we will write to both our output file and mdcrd file every step. In this way we will be able to see the detailed motion of our system at the expense of slowing the simulation down while it writes a considerable amount of information to disk. Here is the input file for the MD run:
300K constant temp MD &cntrl imin=0, ntb=1, cut=8.0, ntc=2, ntf=2, tempi=300.0, temp0=300.0, ntt=3, gamma_ln=1.0, nstlim=1000, dt=0.002, ntpr=1, ntwx=1, / |
|
Now we run the md.
$AMBERHOME/bin/sander -O -i md_classical.in -o md_classical.out -p NMA.prmtop -c min_classical.rst -r md_classical.rst -x md_classical.mdcrd
This should take less than a minute to run.
Input Files: md_classical.in, NMA.prmtop,
min_classical.rst
Output Files: md_classical.out, md_classical.rst,
md_classical.mdcrd
Now, let's do two things. First of all we will visualize our MD trajectory in VMD and secondly we will use cpptraj to create an average structure for us.
Step 1 - VMD. Fire up VMD and load the NMA topology file (should be recognized as AMBER7 Parm) and then load the md_classical.mdcrd (remember that this coordinate file has box information in it so select "AMBER Coordinates with Periodic Box") file into this molecule. I gzipped it to save space so make sure you decompress it first.
While it is good to see our big box of water "floating" about we are more interested in the dynamics of the NMA residue so let's remove the water.
Graphics->Representations
In Selected Atoms type:
all not water
And hit apply
At the same time change the drawing method to CPK.
You should now be able to watch a movie of our NMA molecule "wiggling" about.
The thing to notice is how planar the amide C-N-H group remains during the simulation. This is largely artificial since the nitrogen atom has a lone pair which acts to force the hydrogen out of plane. This is a deficiency of the classical force field as it is currently implemented.
Lets repeat the simulation with QM/MM.
(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