(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 Dwight McGee, Hailey Bureau, Caley Allen, Rigoberto Hernandez
Section 1: Preparing the system
In this tutorial, we focus on a helical peptide deca-alanine, whose structure is similiar to that used by Ozer et al. except with different end caps( ACE and NH2).
Download a pdb of the peptide: ALA.pdb
Generating the topology/coordinate files
Next, we load the peptide into tleap in order to solvate the structure and create the necessary files to run our ASMD simulations. The following script (below) generates the topology and coordinate files
| tleap.in |
#Load the Force Field
source oldff/leaprc.ff14SB
#Load the PDB
alapdb=loadpdb ALA.pdb
#solvate the structure
solvateoct alapdb TIP3PBOX 15.0
#Check & LOADOFF the topology and coordinate files
check alapdb
saveamberparm alapdb ala.asmd.prmtop ala.asmd.rst7
quit
|
Run the following command or type each of the commands above line by line in tleap/xleap.
| >tleap -f tleap.in |
Note: There is no need to add counterions because our system has total net charge of 0
You should now have generated the files, ala.asmd.prmtop and ala.asmd.rst7
With the topology and coordinate files in hand, we need to equilibrate our system. The system can be equilibrated through the following series of steps: minimization, heating, density, slowly releasing restraints, and finally an unconstrained molecular dynamics simulation. We assume that the reader is familar with performing minimizations and running MD simulations so we will not describe each input variable in the scripts outlined below. (Information about each input variable is described in the AMBER manual.)
Minimization
The minization is performed in two parts: First minimize the waters and then minimize the entire system.
| min1.in | min2.in |
Minimize solvent
&cntrl
imin=1, maxcyc=10000, ntmin=1,
ntb=1, ntc=1, ntf=1,
ntr=1, restraintmask="!(:WAT)", restraint_wt=20.0
ntwx=100, ntpr=100, ntwr=100, iwrap=1,
ioutfm=1,
/
|
Minimize solute and solvent
&cntrl
imin=1, maxcyc=10000, ntmin=1,
ntb=1, ntc=1, ntf=1,
ntwx=100, ntpr=100, ntwr=100,
ioutfm=1, iwrap=1,
/
|
Heating
The system is heated from 100 K to 300 K incrementally.
| heat.mdin |
Heat system (constant volume)
&cntrl
imin=0, nstlim=250000, dt=0.002,
ntb=1, ntc=2,ntf=2,
cut=8.0, iwrap=1,
ntpr=10000, ntwx=10000, ntwr=10000,
ntt=3, gamma_ln=5.0, ig=-1,
tempi=100.0, temp0=300.0,
ntr=1, restraintmask=':1-12', restraint_wt=20.0,
nmropt=1, ioutfm=1,
/
&wt TYPE='TEMP0', istep1=0, istep2=250000,
value1=100.0, value2=300.0, /
&wt TYPE='END' /
|
Density Equilibration
Equilibrate the density of our system.
| density.mdin |
Density Equilibration
&cntrl
imin=0, nstlim=250000, dt=0.002
irest=1, ntx=5,
ntc=2, ntf=2, ntb=2, ntp=1, taup=1.0,
cut=8.0, iwrap=1, ig=-1
ntpr=10000, ntwx=10000, ntwr=10000,
ntt=3, gamma_ln=5.0, temp0=300.0,
ntr=1, restraintmask=":1-12", restraint_wt=20.0
ioutfm=1, barostat=2,
/
|
System Equilibration
The restraints of the backbone atoms "C,CA,N,O" are reduced from 10 to 0 through a series of molecular dynamics simulations.
| equil1.mdin | equil2.mdin |
MD1 simulation backbone restrained &cntrl imin=0, nstlim=100000, dt=0.002, irest=1, ntx=5, ntc=2, ntf=2, ntb=1, cut=8.0, iwrap=1, ig=-1, ntpr=5000, ntwx=5000, ntwr=5000, ntt=3, gamma_ln=5.0, temp0=300.0, ntr=1, restraintmask='@CA,C,N,O', restraint_wt=10.0, ioutfm=1, / |
MD2 simulation1 backbone restrained &cntrl imin=0, nstlim=100000, dt=0.002, irest=1, ntx=5, ntc=2, ntf=2, ntb=1, cut=8.0, iwrap=1, ig=-1, ntpr=5000, ntwx=5000, ntwr=5000, ntt=3, gamma_ln=5.0, temp0=300.0, ntr=1, restraintmask='@CA,C,N,O', restraint_wt=5.0, ioutfm=1, / |
| equil3.mdin | equil4.mdin |
MD3 simulation backbone restrained &cntrl imin=0, nstlim=100000, dt=0.002, irest=1, ntx=5, ntc=2, ntf=2, ntb=1, cut=8.0, iwrap=1, ig=-1, ntpr=5000, ntwx=5000, ntwr=5000, ntt=3, gamma_ln=5.0, temp0=300.0, ntr=1, restraintmask='@CA,C,N,O', restraint_wt=1.0, ioutfm=1, / |
MD4 simulation Unconstrained &cntrl imin=0, nstlim=1000000, dt=0.002 irest=1, ntx=5, ntc=2, ntf=2, ntb=1, cut=8.0, iwrap=1, ig=-1 ntpr=10000, ntwx=10000, ntwr=100000, ntt=3, gamma_ln=5.0, temp0=300.0, ioutfm=1, nscm=10000, / |
All of the commands needed to perform the equilibration were combined into one script and run on a local cluster using 16 processors. The equilibration took just over 3 hours to complete.
The inputs/outputs from the equilibration runs can be download here: equilibration.tar.gz
| equilibration.sh |
#!/bin/sh
#PBS -m abe
#PBS -M username@gatech.edu
#PBS -N ASMD_equilibration
#PBS -l walltime=4:00:00
#PBS -l nodes=2:ppn=8
#PBS -j oe
cd $PBS_O_WORKDIR
do_parallel="mpirun -np 16 pmemd.MPI "
prmtop="ala.asmd.prmtop"
#Run the ASMD equilibration
$do_parallel -O -i min1.in -p $prmtop -c ala.asmd.rst7 -x min1.nc -o min1.out -inf min1.info -r min1.rst7 -ref ala.asmd.rst7
$do_parallel -O -i min2.in -p $prmtop -c min1.rst7 -x min2.nc -o min2.out -inf min2.info -r min2.rst7
$do_parallel -O -i heat.mdin -p $prmtop -c min2.rst7 -x heat.nc -o heat.mdout -inf heat.info -r heat.rst7 -ref min2.rst7
$do_parallel -O -i density.mdin -p $prmtop -c heat.rst7 -x density.nc -o density.mdout -inf density.info -r density.rst7 -ref heat.rst7
$do_parallel -O -i equil1.mdin -p $prmtop -c density.rst7 -x equil1.nc -o equil1.mdout -inf equil1.info -r equil1.rst7 -ref density.rst7
$do_parallel -O -i equil2.mdin -p $prmtop -c equil1.rst7 -x equil2.nc -o equil2.mdout -inf equil2.info -r equil2.rst7 -ref equil1.rst7
$do_parallel -O -i equil3.mdin -p $prmtop -c equil2.rst7 -x equil3.nc -o equil3.mdout -inf equil3.info -r equil3.rst7 -ref equil2.rst7
$do_parallel -O -i equil4.mdin -p $prmtop -c equil3.rst7 -x equil4.nc -o equil4.mdout -inf equil4.info -r equil4.rst7
|