SMD of 1E0Q

In this page, we will run a SMD simulation on 1E0Q.

PATH: ~/1E0Q/amber/SMD/

copy 1E0Q.pdb, min.mdin, heat.mdin, md.mdin to ~/1E0Q/amber/SMD/

./pipeline.sh

the only two files we will be using are 1E0Q_md.rst and 1E0Q_box.prmtop

let's run 100A/ns SMD simulation as an example:

Input File Preparation

mkdir 100_ang_per_ns

copy create_ASMDinputs.sh and create_job.sh to ./100_ang_per_ns/

create_ASMDinputs.sh:

#!/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=7
  end_distance=57
  for stage in {1..1};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,iwrap=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=7
end_distance=57

for stage in {1..1}; do
    cat >>dist.RST.dat.$stage<<EOF
 &rst
        iat=9,275,
        r2=$start_distance,
        r2a=$end_distance,
        rk2=7.2,

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

./create_ASMDinputs.sh 5 250000

creat_job.sh

#!/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
  exit
else
  coord=${args[1]}
  
fi

stage=${args[2]}

cat >>_job.sh<<EOF
#!/bin/sh

do_parallel="pmemd.cuda"
prmtop="../1E0Q_box.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

./creat_job.sh 5 ../1E0Q_md.rst 1

Run the simulation

./job.1.sh

 

Repeat the same procedure for 10A/ns.

./create_ASMDinputs.sh 5 2500000

./creat_job.sh 5 ../1E0Q_md.rst 1

Data Analysis

./10ang_per_ns/:

../ASMD.py -i asmd*.work.dat.1 -o jar.stage1.dat

./100ang_per_ns/:

../ASMD.py -i asmd*.work.dat.1 -o jar.stage1.dat

plotting the PMF:

gnuplot:

set title "PMF of Steered Molecular Simulation of 1E0Q"

set xlabel "end-to-end distance(A)"

set ylabel "PMF (kcal/mol)"

plot "10ang_per_ns/jar.stage1.dat" w l ti "100A/ns","100ang_per_ns/jar.stage1.dat" w l ti "10A/ns"

 

As we can see, PMF for faster speed is significantly bigger than those of slower pulling speed, That is the reason when we were trying to explore the difference of gibbs free energy between two confomer using PMF, we need to emply as slow pulling speed as possible so as to make it close to reversible work.

Reference:

http://ambermd.org/tutorials/advanced/tutorial26/