molecular_dynamics/trajectory_analysis/11-rmsd-equilibration.sh

#!/usr/bin/env bash
export VMDNOCUDA=1
#export VMDNOOPTIX=1
export VMDNOOSPRAY=1

. functions.sh

#######################################################################################
## ADJUSTABLE PARAMETERS                                                             ##
#######################################################################################

# set atom selection for alignment
align_sel=${1?"chain A and name CA and resid 81 to 103 118 to 138 195 to 215 221 to 241 308 to 328 331 to 350 860 to 88 912 to 932 991 to 1011 1014 to 1034 1103 to 1123 1129 to 1149"}
# set atom selection for rmsd calculation
rmsd_sel=${2?"chain A and name CA and resid 81 to 103 118 to 138 195 to 215 221 to 241 308 to 328 331 to 350 860 to 88 912 to 932 991 to 1011 1014 to 1034 1103 to 1123 1129 to 1149"}

#######################################################################################
## SET VARS                                                                          ##
#######################################################################################

# set folder name for output files
set_and_make_out_dir
# Get array of concatenated, indvidual jobs
get_jobs_individual

#######################################################################################
## CALCULATE RMSD                                                                    ##
#######################################################################################

function calc_rmsd() {

  for _job in ${_jobs_individual}; do
  
    get_system_and_date ${_job}
    psf=$(ls ${_job}/*.psf)
    dcd=$(ls ${_job}/*.dcd)
  
    # load trajectory
    echo "mol new ${psf}" > temp.tcl
    echo "mol addfile ${dcd} waitfor all" >> temp.tcl
  
    # ALIGN
    echo "set align_ref [atomselect top \""${align_sel}"\" frame 0]" >> temp.tcl
    echo "set align_sel [atomselect top \""${align_sel}"\"]" >> temp.tcl
    echo "set rmsd_ref [atomselect top \""${rmsd_sel}"\" frame 0]" >> temp.tcl
    echo "set rmsd_sel [atomselect top \""${rmsd_sel}"\"]" >> temp.tcl
    echo 'set all [atomselect top "all"]' >> temp.tcl
    echo 'set n [molinfo top get numframes]' >> temp.tcl
    echo "set fp [open "${out_dir}/${_system}-${_date}.csv" w]" >> temp.tcl
    echo 'for { set i 1 } { $i < $n } { incr i } {' >> temp.tcl
    echo '  $align_sel frame $i' >> temp.tcl
    echo '  $rmsd_sel frame $i' >> temp.tcl
    echo '  $all frame $i' >> temp.tcl
    echo '  $all move [measure fit $align_sel $align_ref]' >> temp.tcl
    # CALC RMSD
    echo '  puts $fp "$i [measure rmsd $rmsd_sel $rmsd_ref]"' >> temp.tcl
    echo '}' >> temp.tcl
    echo 'close $fp'  >> temp.tcl
    echo 'exit' >> temp.tcl
  
    vmd -dispdev text -eofexit < "temp.tcl"
    sync
    rm "temp.tcl"
    
    plot_rmsd ${_system} "${out_dir}/${_system}-${_date}"
    
  done
  
}


function plot_rmsd() {
gnuplot -persist <<-GNUPLOT_INPUT
  set terminal pngcairo enhanced size 1200,1200 font "arial,24" linewidth 2
  set output "${2}.png"
  set autoscale
  #set title "Backbone RMSD" font "arial,40"
  #set xlabel "Time (ns)" font "arial,32"
  #set ylabel "RMSD (Å)" font "arial,32"
  set border 1+2
  #set xtics nomirror
  set ytics nomirror 1
  set grid
  set grid noxtics
  set key box off
  #set yrange [0:4]
  #set xrange [0:120]
  #set format y "%+-2.f"
  set style data lines
  set xzeroaxis ls 8
  
  # linetypes
  #set style line 1 pi 0 lw 1 ps 1 pt 0 lc rgb '#000000'

  plot "${2}.csv" using (\$1 / 50):(\$2 + 1) with linespoints ls 8 pi 0 pt 0
  
  exit

GNUPLOT_INPUT
}
  


calc_rmsd