(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

Now that we have completed the equilibration, we will proceed to perform the ASMD simulations. The reaction coordinate is the end-to-end distance between the CA atoms of ALA1 and ALA10. We will pull the helical peptide using two different velocities, 10 Å/ns and 100 Å/ns. For simplicity, the reaction coordinate is partitioned into 5 equal segments (4 Å in length) and we use only 25 trajectories per stage. More trajectories would be needed to obtain a converged PMF but this tutorial is illustrated with this smaller number in order to allow the reader to walk through the steps quickly.

Create MD inputs and Distance restraint files

100 Å/ns velocity ASMD Simulation

At this pulling speed, each stage should be 40ps (nstlim=20000, dt=.002) in duration in order to match the length of the pull, etc. Each simulation is run with respect to the NVT ensemble (ntb=1). The temperature is maintained using Langevin dynamics (ntt=3, gamma_ln=5.0). Always remember to set ig=-1 in order to prevent synchronization! Setting jar=1 turns on steered molecular dynamics. Note that setting jar=1 internally sets nmropt=1 so there is no need to place nmropt in your MDIN. DISANG is the name of the file from which to read the distance restratints. DUMPAVE is the name of the output file which will contain the 4 columns: target distance, actual distance, force, and work This file will be written to every 1000 steps (istep1=1000).

ASMD Simulation Example MD Input
 &cntrl
   imin = 0, nstlim = 20000, dt = 0.002,
   ntx = 1, temp0 = 300.0,
   ntt = 3, gamma_ln=5.0
   ntc = 2, ntf = 2, ntb =1,
   ntwx =  1000, ntwr = 20000, ntpr = 1000,
   cut = 8.0, ig=-1, ioutfm=1,
   irest = 0, jar=1,
 /
 &wt type='DUMPFREQ', istep1=1000 /
 &wt type='END'   /
DISANG=dist.example.RST
DUMPAVE=dist.asmd.work.dat
LISTIN=POUT
LISTOUT=POUT
    

The distance restraint file resembles the file inputs used for NMR restraints. In the example, shown below, the parameters iat=9,99, r2=13, r2a=17 rk2=7.2 correspond to a request to change the distance between atoms 9 and 99 from 13 Å to 17 Å using a force constant of 7.2

 &rst
        iat=9,99,
        r2=13,
        r2a=17,
        rk2=7.2,

 &end
      

Create a directory for the asmd simulations that will use a velocity of a 100 Å/ns (e.g. mkdir 100ang_per_ns; cd 100ang_per_ns). The following bash script was used to create the directories and md inputs for the asmd tutorials:

#!/bin/bash

args=("$@")

#check the correct num of args
if [[ ${#args[@]} == 0 ]] || [[ ${#args[@]} > 2 ]]; then
  echo "create_ASMDinputs.sh {NUM of trajectories} {NUM of MDsteps}"
  exit
fi

num_asmd_sim=${args[0]}
mdsteps=${args[1]}

#create the ASMD directories
counter=1
for ((counter=1; counter<=$num_asmd_sim; counter++)); do
  if [ ! -d ASMD_$counter ]; then
    mkdir ASMD_$counter
  fi
done


#create the MDIN inputs & dist RST files for the ASMD simulation
for ((counter=1; counter<=$num_asmd_sim;counter++));do
  start_distance=13
  end_distance=17
  for stage in {1..5};do
    cat >>asmd_$counter.$stage.mdin<<EOF
ASMD simulation
 &cntrl
   imin = 0, nstlim = $mdsteps, dt = 0.002,
   ntx = 1, temp0 = 300.0,
   ntt = 3, gamma_ln=5.0
   ntc = 2, ntf = 2, ntb =1,
   ntwx =  1000, ntwr = $mdsteps, ntpr = 1000,
   cut = 8.0, ig=-1, ioutfm=1,
   irest = 0, jar=1, 
 /
 &wt type='DUMPFREQ', istep1=1000 /
 &wt type='END'   /
DISANG=dist.RST.dat.$stage
DUMPAVE=asmd_$counter.work.dat.$stage
LISTIN=POUT
LISTOUT=POUT
EOF
  mv asmd_$counter.$stage.mdin ASMD_$counter/
  done
done

#create the distance RST files
start_distance=13
end_distance=17
for stage in {1..5}; do
    cat >>dist.RST.dat.$stage<<EOF
 &rst
        iat=9,99,
        r2=$start_distance,
        r2a=$end_distance,
        rk2=7.2,

 &end
EOF
  start_distance=$((start_distance+4))
  end_distance=$((end_distance+4))
done
        

You may now type the following command, supplying create_ASMDinputs.sh with the additional choices of the number of trajectories as 25 and the number of MD steps as 20,000.

Running the ASMD simulations

Now it is time to run the first stage of ASMD simulations. Before you start, make sure that you have copied the prmtop and equilibrated rst7 file to the directory. Note that the bash script create_job.sh provided below creates the job submission scripts and run each SMD simulations sequentially. However, you could instead use multi-pmemd and create groupfiles to run all the SMD simulation simultaneously.

#!/bin/bash

args=("$@")

#check the correct num of args
if [[ ${#args[@]} == 0 ]] || [[ ${#args[@]} < 3 ]]; then
  echo "create_ASMDjobfile.sh {NUM of trajectories} {Coord/RST7} {Stage Num}"
  exit
fi

num_asmd_sim=${args[0]}

if [ ! -f ${args[1]} ]; then
  echo The file ${args[1]} does not exist
else
  coord=${args[1]}
  
fi

stage=${args[2]}

cat >>_job.sh<<EOF
#!/bin/sh
#PBS -m abe
#PBS -M username@gatech.edu
#PBS -N ASMD_stage$stage
#PBS -l walltime=4:00:00
#PBS -l nodes=2:ppn=8
#PBS -j oe

cd SSSPBS_O_WORKDIR

do_parallel="mpirun -np 16 pmemd.MPI"
prmtop="ala.asmd.prmtop"
coord="$coord"

EOF

for ((counter=1;counter<=$num_asmd_sim;counter+=1)); do
  echo SSSdo_parallel -O -i ASMD_$counter/asmd_$counter.$stage.mdin -p SSSprmtop -c SSScoord
-r ASMD_$counter/ASMD_$counter.$stage.rst7 -o ASMD_$counter/ASMD_$counter.$stage.mdout 
-x ASMD_$counter/ASMD_$counter.$stage.nc -inf ASMD_$counter/ASMD_$counter.$stage.info >> _job.sh
done
sed 's/SSS/$/g' _job.sh > job.$stage.sh; rm -f _job.sh
    

In order for create_job.sh to run, it needs the number of trajectories/simulations 25, coordinate file equil.rst7 and the stage number 1.

Stage 1

All of the work output files for the 10 Å/ns ASMD simulation can be downloaded here.

Submit the file job.1.sh to the queue

Note: you might have to edit the create_job.sh script to satisfy your own queueing system. The script is provided so as to illustrate how one might submit a job to a given queue.

On our computers (running on 16 cores of a 2.0 GHz AMD-Opteron-6128 processor, the average scaling was 23 ns/day), the job finished in approximately 1 hr.

Next, we used the program ASMD.py to determine, which of the ASMD simulations' work values is the closest to Jarzynski Average

Type the following command in order to obtain a list of available options and a description of them for ASMD.py.

In order to determine, which simulation has the closest work value execute the following command:

which should result in the following being displayed

Note: Do not be alarmed if after you run the simulation, the file chosen does not match the one listed above. (The result varies because of the use of Langevin dynamics with the random seeds and variations in the numerics of your computers random number generator).

Stage 2

Remember 25 for the number of simulations, asmd_1/ASMD_1.1.rst7 is the coordinate file and 2 for the stage number. Determine the simulation whose work value is closest to the Jarzynski Average

Stage 3

Using the following commands below: create the jobfile for the next stage, submit it to the queue and use ASMD.py to determine which trajectory's work is the closet to the Jarzynski average.

Stage 4

Using the following commands below: create the jobfile for the next stage, submit it to the queue and use ASMD.py to determine which trajectory's work is the closet to the Jarzynski average.

Stage 5

Using the following commands below: create the jobfile for the next stage, submit it to the queue and use ASMD.py to determine which trajectory's work is the closet to the Jarzynski average.

10 Å/ns velocity ASMD Simulation

Create a directory for the ASMD simulations using a velocity of 10 Å/ns (e.g. mkdir 10ang_per_ns; cd 10ang_per_ns).

Type the following command below supplying create_ASMDinputs.sh with the number of trajectories, 25 and the number of MD steps 200,000.

Follow the same procedure as outlined previously. Please remember that these simulations at a slower velocity of 10 Å/ns will require much more time complete.

All of the work output files for the 10 Å/ns ASMD simulation can be downloaded here.