molecular_dynamics/auto_namd/functions.sh

#!/usr/bin/env bash

function get_next_step() {

  # get system name from path
  system=$(basename "${job}")

  # get prev_step and prev_step_num from parsing coor filenames
  last_step_coorname=$(find ${job} -name "*.coor" -type f | sort -V | tail -n 1)
  if [[ "${last_step_coorname}" == "" ]]; then
    prev_step=""
  else
    filename=$(basename ${last_step_coorname})
    filename=${filename%%.*}
    prev_step=${filename##${system}_}
    prev_step=${prev_step%%_*}
  fi
  case "${prev_step}" in
    "4-sim")
      this_step="4-sim"
      prev_step_num=${filename##*[^0-9]}
      final_step_num=$((${prev_step_num}+${step_size})) 
      ;;
    "3-heat")
      this_step="4-sim"
      prev_step_num=0
      final_step_num=${step_size} 
      ;;
    "2-min")
      this_step="3-heat"
      prev_step_num=0
      final_step_num=0 
      ;;
    "1-min")
      this_step="2-min"
      prev_step_num=0
      final_step_num=0 
      ;;
    "")
      if [[ -f "${job}/${system}.pdb" && -f "${job}/${system}.psf" ]]; then
        this_step="1-min"
        prev_step_num=0
        final_step_num=0
      else
        echo "No requisite pdb or psf files, exiting..."
        exit 1
      fi
      ;;
  esac
}


function create_conf() {

  step_size=$((${final_step_num} - ${prev_step_num}))
  
  echo "Creating conf file using job:"${job}" system:"${system}" this_step:"$this_step" prev_step:$prev_step_num final_step:$final_step_num"

  #######################################################################################
  ## CREATE NAMD CONFIGURATION FILES                                                   ##
  #######################################################################################

  case "${this_step}" in
  
    ##########################
    ## MINIMIZATION 1       ##
    ##########################
    "1-min")
    
      conf="${job}/${system}_${this_step}.conf"
      out="${job}/${system}_${this_step}.out"
      
      echo "structure ${system}.psf" > ${conf}
      echo "coordinates ${system}.pdb" >> ${conf}
      echo "paraTypeCharmm on" >> ${conf}
      for ff in "${ffs[@]}"; do
        echo "parameters "${ff}"" >> ${conf}
      done
      echo "temperature 0" >> ${conf}
      echo "" >> ${conf}
      echo "# Force-Field Parameters" >> ${conf}
      echo "exclude scaled1-4" >> ${conf}
      echo "1-4scaling 1.0" >> ${conf}
      echo "cutoff 12." >> ${conf}
      echo "switching on" >> ${conf}
      echo "switchdist 10." >> ${conf}
      echo "pairlistdist 14" >> ${conf}
      echo "" >> ${conf}
      echo "# Integrator Parameters" >> ${conf}
      echo "timestep 1.0  ;# 1fs/step" >> ${conf}
      echo "nonbondedFreq 1" >> ${conf}
      echo "fullElectFrequency 2" >> ${conf}
      echo "stepspercycle 10" >> ${conf}
      echo "" >> ${conf}
      echo "# Output" >> ${conf}
      echo "outputName ${system}_${this_step}" >> ${conf}
      echo "outputEnergies 100" >> ${conf}
      echo "outputPressure 100" >> ${conf}
      # Minimize
      echo "minimize 15000" >> ${conf}
      ;;
    ##########################
    ## MINIMIZATION 2       ##
    ##########################
    "2-min")
     
      conf="${job}/${system}_${this_step}.conf"
      out="${job}/${system}_${this_step}.out"
     
      # Solvate and ionize
      solv_ion
      
      # create pcell.txt
      calc_pcell_min2

      # Source pcell.txt for PBC values
      source "${job}/${system}_${this_step}_pcell.txt"

      # Begin writing NAMD conf file
      echo "structure ${system}_solv_ion.psf" > ${conf}
      echo "coordinates ${system}_solv_ion.pdb" >> ${conf}
      echo "paraTypeCharmm on" >> ${conf}
      for ff in "${ffs[@]}"; do
        echo "parameters "${ff}"" >> ${conf}
      done
      echo "temperature 0" >> ${conf}
      echo "" >> ${conf}
      echo "# Periodic boundary conditions" >> ${conf}
      echo "wrapWater on" >> ${conf}
      echo "wrapAll on" >> ${conf}
      echo "cellOrigin	$CORX	$CORY	$CORZ" >> ${conf}
      echo "cellBasisVector1	$CBVX	0.0	0.0" >> ${conf}
      echo "cellBasisVector2	0.0	$CBVY	0.0" >> ${conf}
      echo "cellBasisVector3	0.0	0.0	$CBVZ" >> ${conf}
      echo "# Force-Field Parameters" >> ${conf}
      echo "exclude scaled1-4" >> ${conf}
      echo "1-4scaling 1.0" >> ${conf}
      echo "cutoff 12." >> ${conf}
      echo "switching on" >> ${conf}
      echo "switchdist 10." >> ${conf}
      echo "pairlistdist 14" >> ${conf}
      echo "" >> ${conf}
      echo "# Integrator Parameters" >> ${conf}
      echo "timestep 2.0  ;# 2fs/step" >> ${conf}
      echo "rigidBonds all  ;# needed for 2fs steps" >> ${conf}
      echo "nonbondedFreq 1" >> ${conf}
      echo "fullElectFrequency 2" >> ${conf}
      echo "stepspercycle 10" >> ${conf}
      echo "" >> ${conf}
      echo "# Output" >> ${conf}
      echo "outputName ${system}_${this_step}" >> ${conf}
      echo "outputEnergies 10000" >> ${conf}
      echo "outputPressure 10000" >> ${conf}
      # Minimize
      echo "minimize 15000" >> ${conf}
      ;;
    ##########################
    ## EQUILIBRATION        ##
    ##########################
    "3-heat")
    
      conf="${job}/${system}_${this_step}.conf"
      out="${job}/${system}_${this_step}.out"
      
      # create pcell.txt
      calc_pcell_heat

      # Source pcell.txt for PBC values
      source "${job}/${system}_${this_step}_pcell.txt"

      echo "structure ${system}_solv_ion.psf" > ${conf}
      echo "coordinates ${system}_solv_ion.pdb" >> ${conf}
      echo "bincoordinates ${system}_2-min.coor" >> ${conf}
      echo "#extendedSystem ${system}_2-min.xsc" >> ${conf}
      echo "cellOrigin	$CORX	$CORY	$CORZ" >> ${conf}
      echo "cellBasisVector1	$CBVX	0.0	0.0" >> ${conf}
      echo "cellBasisVector2	0.0	$CBVY	0.0" >> ${conf}
      echo "cellBasisVector3	0.0	0.0	$CBVZ" >> ${conf}
      echo "" >> ${conf}
      echo "firsttimestep ${prev_step_num}" >> ${conf}
      echo "" >> ${conf}
      echo "paraTypeCharmm on" >> ${conf}
      for ff in "${ffs[@]}"; do
        echo "parameters "${ff}"" >> ${conf}
      done
      echo "temperature 0" >> ${conf}
      echo "#margin 2" >> ${conf}
      echo "" >> ${conf}
      echo "# Periodic boundary conditions" >> ${conf}
      echo "wrapWater on" >> ${conf}
      echo "wrapAll on" >> ${conf}
      echo "" >> ${conf}
      echo "# Full Electrostatics" >> ${conf}
      echo "PME on" >> ${conf}
      echo "PMEGridSpacing 1.0" >> ${conf}
      echo "# Force-Field Parameters" >> ${conf}
      echo "exclude scaled1-4" >> ${conf}
      echo "1-4scaling 1.0" >> ${conf}
      echo "cutoff 12." >> ${conf}
      echo "switching on" >> ${conf}
      echo "switchdist 10." >> ${conf}
      echo "pairlistdist 14" >> ${conf}
      echo "" >> ${conf}
      echo "# Integrator Parameters" >> ${conf}
      echo "timestep 2.0  ;# 2fs/step" >> ${conf}
      echo "rigidBonds all  ;# needed for 2fs steps" >> ${conf}
      echo "nonbondedFreq 1" >> ${conf}
      echo "fullElectFrequency 2  " >> ${conf}
      echo "stepspercycle 10" >> ${conf}
      echo "" >> ${conf}
      echo "# Output" >> ${conf}
      echo "outputName ${system}_${this_step}" >> ${conf}
      echo "outputEnergies 1000" >> ${conf}
      echo "outputPressure 1000" >> ${conf}
      echo "dcdfreq 1000" >> ${conf}
      echo "" >> ${conf}
      echo "# Constant Temperature Control" >> ${conf}
      echo "langevin on    ;# do langevin dynamics" >> ${conf}
      echo "langevinDamping 0.5   ;# damping coefficient (gamma) of 0.5/ps" >> ${conf}
      echo "langevinTemp 310" >> ${conf}
      echo "langevinHydrogen no    ;# don't couple langevin bath to hydrogens" >> ${conf}
      echo "" >> ${conf}
      echo "# Constant Pressure Control (variable volume)" >> ${conf}
      echo "useGroupPressure yes ;# needed for 2fs steps" >> ${conf}
      echo "useFlexibleCell yes ;# no for water box, yes for membrane" >> ${conf}
      echo "useConstantRatio yes  ;# no for water box, yes for membrane" >> ${conf}
      echo "langevinPiston on" >> ${conf}
      echo "langevinPistonTarget 1.01325 ;#  in bar -> 1 atm" >> ${conf}
      echo "langevinPistonPeriod 100." >> ${conf}
      echo "langevinPistonDecay 50." >> ${conf}
      echo "langevinPistonTemp 310" >> ${conf}
      echo "" >> ${conf}
      # Equilibration loop
      echo "set freq 100" >> ${conf}
      echo 'for {set i 10} {$i <= 310} {incr i 10} {' >> ${conf}
      echo 'reinitvels $i' >> ${conf}
      echo 'langevinTemp $i' >> ${conf}
      echo 'run $freq' >> ${conf}
      echo '}' >> ${conf}
      # Stabilize
      echo 'run 10000' >> ${conf}
      ;;
    ##########################
    ## SIMULATION           ##
    ##########################
    "4-sim")
    
      conf="${job}/${system}_${this_step}_${final_step_num}.conf"
      out="${job}/${system}_${this_step}_${final_step_num}.out"
      
      # Begin writing NAMD conf file
      echo "structure ${system}_solv_ion.psf" > ${conf}
      echo "coordinates ${system}_solv_ion.pdb" >> ${conf}

      # if this is the first timestep, begin with heat
      if [[ ${prev_step_num} == 0 ]]; then
        echo "temperature 310" >> ${conf}
        echo "bincoordinates ${system}_3-heat.coor" >> ${conf}
        # create pcell.txt
        calc_pcell_heat
        source "${job}/${system}_${this_step}_pcell.txt"
        echo "cellOrigin	$CORX	$CORY	$CORZ" >> ${conf}
        echo "cellBasisVector1	$CBVX	0.0	0.0" >> ${conf}
        echo "cellBasisVector2	0.0	$CBVY	0.0" >> ${conf}
        echo "cellBasisVector3	0.0	0.0	$CBVZ" >> ${conf}
      # else resume
      else
        echo "binvelocities ${system}_${this_step}_${prev_step_num}.vel" >> ${conf}
        echo "bincoordinates ${system}_${this_step}_${prev_step_num}.coor" >> ${conf}
        echo "extendedSystem ${system}_${this_step}_${prev_step_num}.xsc" >> ${conf}
      fi

      echo "" >> ${conf}
      echo "firsttimestep ${prev_step_num}" >> ${conf}
      echo "" >> ${conf}
      echo "#margin 2" >> ${conf}
      echo "paraTypeCharmm on" >> ${conf}
      for ff in "${ffs[@]}"; do
        echo "parameters "${ff}"" >> ${conf}
      done
      echo "" >> ${conf}
      echo "# Periodic boundary conditions" >> ${conf}
      echo "wrapWater on" >> ${conf}
      echo "wrapAll on" >> ${conf}
      echo "" >> ${conf}
      echo "# Full Electrostatics" >> ${conf}
      echo "PME on" >> ${conf}
      echo "PMEGridSpacing 1.0" >> ${conf}
      echo "" >> ${conf}
      echo "# Force-Field Parameters" >> ${conf}
      echo "exclude scaled1-4" >> ${conf}
      echo "1-4scaling 1.0" >> ${conf}
      echo "cutoff 12." >> ${conf}
      echo "switching on" >> ${conf}
      echo "switchdist 10." >> ${conf}
      echo "pairlistdist 14" >> ${conf}
      echo "" >> ${conf}
      echo "# Integrator Parameters" >> ${conf}
      echo "timestep 2.0  ;# 2fs/step" >> ${conf}
      echo "rigidBonds all  ;# needed for 2fs steps" >> ${conf}
      echo "nonbondedFreq 1" >> ${conf}
      echo "fullElectFrequency 2  " >> ${conf}
      echo "stepspercycle 10" >> ${conf}
      echo "" >> ${conf}
      echo "# Output" >> ${conf}
      echo "outputName ${system}_${this_step}_${final_step_num}" >> ${conf}
      echo "outputEnergies 10000" >> ${conf}
      echo "outputPressure 10000" >> ${conf}
      echo "dcdfreq 10000" >> ${conf}
      echo "" >> ${conf}
      echo "# Constant Temperature Control" >> ${conf}
      echo "langevin on    ;# do langevin dynamics" >> ${conf}
      echo "langevinDamping 0.5   ;# damping coefficient (gamma) of 0.5/ps" >> ${conf}
      echo "langevinTemp 310" >> ${conf}
      echo "langevinHydrogen no    ;# don't couple langevin bath to hydrogens" >> ${conf}
      echo "" >> ${conf}
      echo "# Constant Pressure Control (variable volume)" >> ${conf}
      echo "useGroupPressure yes ;# needed for 2fs steps" >> ${conf}
      echo "useFlexibleCell yes ;# no for water box, yes for membrane" >> ${conf}
      echo "useConstantRatio yes  ;# no for water box, yes for membrane" >> ${conf}
      echo "langevinPiston on" >> ${conf}
      echo "langevinPistonTarget 1.01325 ;#  in bar -> 1 atm" >> ${conf}
      echo "langevinPistonPeriod 100." >> ${conf}
      echo "langevinPistonDecay 50." >> ${conf}
      echo "langevinPistonTemp 310" >> ${conf}
      echo "" >> ${conf}
      echo "run ${step_size}" >> ${conf}
      ;;
  esac
}

function solv_ion() {
  # Create temporary VMD script to solvate, ionize, and generate PBC measurements
  # Load molecule
  echo "mol load psf [glob ${job}/${system}.psf] namdbin [glob ${job}/${system}_${prev_step}.coor]" > ${job}/temp_solv_ion.tcl
  echo "set all [atomselect top all]" >> ${job}/temp_solv_ion.tcl
  echo "\$all writepdb ${job}/${system}_${prev_step}.pdb" >> ${job}/temp_solv_ion.tcl
  # Solvate
  echo "package require solvate" >> ${job}/temp_solv_ion.tcl
  echo "solvate ${job}/${system}.psf ${job}/${system}_${prev_step}.pdb -o ${job}/${system}_solv -s WT -x 13 -y 13 -z 13 +x 13 +y 13 +z 13 -b 2.4" >> ${job}/temp_solv_ion.tcl
  # Ionize
  echo "package require autoionize" >> ${job}/temp_solv_ion.tcl
  echo "autoionize -psf ${job}/${system}_solv.psf -pdb ${job}/${system}_solv.pdb -o ${job}/${system}_solv_ion -sc 0.15" >> ${job}/temp_solv_ion.tcl

  # Run VMD
  vmd -dispdev text -eofexit < "${job}/temp_solv_ion.tcl"
  rm ${job}/temp_solv_ion.tcl
}

function calc_pcell_min2() {
  # Load molecule
  echo "mol load psf [glob ${job}/${system}_solv_ion.psf] pdb [glob ${job}/${system}_solv_ion.pdb]" > ${job}/temp_build_box.tcl

  # Get PBC box measurements from VMD and output to pcell.txt
  echo "set fileId [open ${job}/${system}_${this_step}_pcell.txt "w"]" >> ${job}/temp_build_box.tcl
  echo "set all [atomselect top "not lipid"]" >> ${job}/temp_build_box.tcl
  echo 'set minmax [measure minmax $all]' >> ${job}/temp_build_box.tcl
  echo 'set vec [vecsub [lindex $minmax 1] [lindex $minmax 0]]' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "#!/usr/bin/env bash"' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "CBVX=[lindex $vec 0]"' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "CBVY=[lindex $vec 1]"' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "CBVZ=[lindex $vec 2]"' >> ${job}/temp_build_box.tcl
  echo 'set center [measure center $all]' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "CORX=[lindex $center 0]"' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "CORY=[lindex $center 1]"' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "CORZ=[lindex $center 2]"' >> ${job}/temp_build_box.tcl
  echo 'close $fileId' >> ${job}/temp_build_box.tcl

  # Run VMD
  vmd -dispdev text -eofexit < "${job}/temp_build_box.tcl"
  rm "${job}/temp_build_box.tcl"
}

function calc_pcell_heat() {
  # Load molecule
  echo "mol load psf [glob ${job}/${system}_solv_ion.psf] namdbin [glob ${job}/${system}_${prev_step}.coor]" > ${job}/temp_build_box.tcl

  # Get PBC box measurements from VMD and output to pcell.txt
  echo "set fileId [open ${job}/${system}_${this_step}_pcell.txt "w"]" >> ${job}/temp_build_box.tcl
  echo "set all [atomselect top all]" >> ${job}/temp_build_box.tcl
  echo 'set minmax [measure minmax $all]' >> ${job}/temp_build_box.tcl
  echo 'set vec [vecsub [lindex $minmax 1] [lindex $minmax 0]]' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "#!/usr/bin/env bash"' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "CBVX=[lindex $vec 0]"' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "CBVY=[lindex $vec 1]"' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "CBVZ=[lindex $vec 2]"' >> ${job}/temp_build_box.tcl
  echo 'set center [measure center $all]' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "CORX=[lindex $center 0]"' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "CORY=[lindex $center 1]"' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "CORZ=[lindex $center 2]"' >> ${job}/temp_build_box.tcl
  echo 'close $fileId' >> ${job}/temp_build_box.tcl

  # Run VMD
  vmd -dispdev text -eofexit < "${job}/temp_build_box.tcl"
  rm "${job}/temp_build_box.tcl"
}

function calc_pcell_sim() {
  # Load molecule
  echo "mol load psf [glob ${job}/${system}_solv_ion.psf] namdbin [glob ${job}/${system}_${prev_step}_${prev_step_num}.coor]" > ${job}/temp_build_box.tcl

  # Get PBC box measurements from VMD and output to pcell.txt
  echo "set fileId [open ${job}/${system}_${this_step}_${prev_step_num}_pcell.txt "w"]" >> ${job}/temp_build_box.tcl
  echo "set all [atomselect top all]" >> ${job}/temp_build_box.tcl
  echo 'set minmax [measure minmax $all]' >> ${job}/temp_build_box.tcl
  echo 'set vec [vecsub [lindex $minmax 1] [lindex $minmax 0]]' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "#!/usr/bin/env bash"' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "CBVX=[lindex $vec 0]"' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "CBVY=[lindex $vec 1]"' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "CBVZ=[lindex $vec 2]"' >> ${job}/temp_build_box.tcl
  echo 'set center [measure center $all]' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "CORX=[lindex $center 0]"' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "CORY=[lindex $center 1]"' >> ${job}/temp_build_box.tcl
  echo 'puts $fileId "CORZ=[lindex $center 2]"' >> ${job}/temp_build_box.tcl
  echo 'close $fileId' >> ${job}/temp_build_box.tcl

  # Run VMD
  vmd -dispdev text -eofexit < "${job}/temp_build_box.tcl"
  rm "${job}/temp_build_box.tcl"
}

function copy_ffs() {
  ffs=()
  for ff in ../0-forcefields/*; do
    if [[ -f "${ff}" ]]; then
      # strip path
      ff_name=${ff##*/}
      # copy ff and append to ffs array
      cp -f "${ff}" "${job}/${ff_name}"
      ffs+=(${ff##*/})
    fi
  done
}

function backup() {
  if [[ -f "${out}" ]]; then
   mv "${out}" "${out}.bak"
  fi
}