123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441 |
- #!/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
- }
|