gmx_funcs.sh

#!/bin/bash

########### COLOURS

source $MWSHPATH/colours.sh
source $MWSHPATH/out.sh
colBold="\033[1m"
colHead="\033[1;4m"
colClear="\033[0m"

if [[ $GMXEXE != "" ]] ; then
  echo -e "$colBold""Using executable $colFunc$GMXEXE""$colClear"
else
  echo -e "$colBold""Using default executable $colFunc""gmx_mpi""$colClear"
  export GMXEXE=gmx_mpi
fi

########### NON-GROMACS

function finishUp {
  if [ $1 -eq 0 ] ; then 
    echo -e $colSuccess"Exiting $1"$colClear
  elif [ $1 -eq -1 ] ; then 
    echo -e $colError"Exiting"$colClear
  else 
    echo -e $colError"Exiting $1"$colClear 
  fi
  rm \#*
  headerOut "Issues:"
  grep -H --color=always -P "(?=There were).*" _*log
  grep -H --color=always -P "(?=There was).*" _*log
  grep -H --color=always -o "Fatal error" _*log
  if [[ -z ${SLURM_JOBID} ]] ; then
    if [[ $- == *i* ]] ; then
      return $1
    else
      exit $1
    fi
  else
    # mkdir ${SLURM_SUBMIT_DIR}/${SLURM_JOBID}_${SLURM_JOB_NAME}
    # cp -r * ${SLURM_SUBMIT_DIR}/${SLURM_JOBID}_${SLURM_JOB_NAME}/
    # mv ${SLURM_SUBMIT_DIR}/${SLURM_JOBID}_${SLURM_JOB_NAME}.? ${SLURM_SUBMIT_DIR}/${SLURM_JOBID}_${SLURM_JOB_NAME}/
    # rm -rf "/users/mw00368/parallel_scratch/"${SLURM_JOBID}_${SLURM_JOB_NAME}/
    exit $1
  fi
}

function cleanDir {
  rm -f \#*
  rm -f *.gro
  rm -f topol*.top
  rm -f topol*.itp
  rm -f posre*.itp
  rm -f ions.tpr
  rm -f *.log
  rm -f mdout.mdp
  rm -f *.xvg
  rm -f *.png
  rm -f *.gif
  rm -f *.cpt
  rm -f *.edr
  rm -f *.tpr
  rm -f *.trr
  rm -f *.xtc
  rm -f *.bak
  rm -f __temp__*
}

function gmxRet {
  GMX_RET=$1; echo -ne "Done, exit code: "
  if [ $GMX_RET -eq 0 ] ; then 
    echo -e $colSuccess"$GMX_RET "$colClear$(/usr/bin/date)
  else 
    echo -e $colError"$GMX_RET "$colClear$(/usr/bin/date)
  fi

  if [[ -z $THIS_SECTION ]] ; then
    THIS_SECTION=-1
  fi
  
  if [ $GMX_RET -ne 0 ] ; then finishUp $THIS_SECTION; fi
  return 0
}

function removeWater {
  echo -e $colBold"Removing water from PDB."$colClear
  grep -v HOH $1 > $2
}

# breakCheck $thisSection $START_AFTER $STOP_AFTER
function breakCheck {
  THIS_SECTION=$1
  # if this state is less than $START_AFTER
  if [ $1 -lt $2 ] ; then 
    # echo "Not started. Omitting $1."
    return 1
  fi
  # if this state is greater than $STOP_AFTER
  if [ $1 -gt $3 ] ; then 
    # echo "Program finished."
    finishUp 0 
  fi
  # echo "Running $1."
  return 0
}

# fancyOut <HEADER>
ECHO_STR=" >>> "
function fancyOut {
  echo -e $colBold$ECHO_STR$1" [ "$(/usr/bin/date)" ] "$colClear
}

########### GROMACS

function gromax {
  COMMAND=$1
  # echo $COMMAND
  # echo $@
  shift
  if [[ "$1" == "-l" ]] ; then
    shift
    LOGNUM=$1
    shift
  fi
  echo -ne $colBold"Running $COMMAND $LOGNUM... "
  case "$COMMAND" in
    # gyrate)
    #   VARIABLE=$1
    #   shift
    #   echo "$VARIABLE 0" | $GMXEXE $COMMAND $@ > _$COMMAND$LOGNUM.log 2>&1
    #   GMX_RET=$?
    #   ;;
    trjconv)
      if [[ "$1" == "-inter" ]] ; then 
        shift
        $GMXEXE $COMMAND $@
        GMX_RET=$?
      elif [[ "$1" == "-custom" ]] ; then
        shift
        STRING=$1
        shift
        echo -e "$STRING \n" | $GMXEXE $COMMAND $@ > _$COMMAND$LOGNUM.log 2>&1
        GMX_RET=$?
      elif [[ "$1" == "-2groups" ]] ; then
        shift
        VAR1=$1
        VAR2=$2
        shift
        shift
        echo "$VAR1 \n $VAR2 \n" | $GMXEXE $COMMAND $@ > _$COMMAND$LOGNUM.log 2>&1
        GMX_RET=$?
      else
        VARIABLE=$1
        shift
        echo -e "$VARIABLE \n $VARIABLE \n" | $GMXEXE $COMMAND $@ > _$COMMAND$LOGNUM.log 2>&1
        GMX_RET=$?
      fi
      ;;
    energy|gyrate|velacc)
      # echo $@
      if [[ "$1" == "-inter" ]] ; then 
        shift
        $GMXEXE $COMMAND $@
        GMX_RET=$?
      elif [[ "$1" == "-2groups" ]] ; then
        shift
        VAR1=$1
        VAR2=$2
        shift
        shift
        echo "$VAR1 $VAR2" | $GMXEXE $COMMAND $@ > _$COMMAND$LOGNUM.log 2>&1
        GMX_RET=$?
      else
        VARIABLE=$1
        shift
        echo -e "$VARIABLE \n" | $GMXEXE $COMMAND $@ > _$COMMAND$LOGNUM.log 2>&1
        GMX_RET=$?
      fi
      ;;
    genion)
      if [[ "$1" == "-inter" ]] ; then 
        shift
        $GMXEXE $COMMAND $@
        GMX_RET=$?
      else
        VARIABLE=$1
        shift
        echo "$VARIABLE" | $GMXEXE $COMMAND $@ > _$COMMAND$LOGNUM.log 2>&1
        GMX_RET=$?
      fi
      ;;
    rms)
      if [[ "$1" == "-inter" ]] ; then 
        shift
        $GMXEXE $COMMAND $@
        GMX_RET=$?
      else
        VARIABLE=$1
        shift
        echo -e "$VARIABLE $VARIABLE" | $GMXEXE $COMMAND $@ > _$COMMAND$LOGNUM.log 2>&1
        GMX_RET=$?
      fi
      ;;
    pdb2gmx)
      while test $# -gt 0; do
        case "$1" in
          -2ter)
            shift
            TER1=$1
            shift
            TER2=$1
            shift
            echo "$TER1 $TER2" | $GMXEXE pdb2gmx $@ -ter > _$COMMAND$LOGNUM.log 2>&1
            GMX_RET=$?
            break
            ;;
          -nter)
            shift
            TERS="$1"
            echo "<<<""$TERS"">>>"
            shift
            echo "$TERS" | $GMXEXE pdb2gmx $@ -ter > _$COMMAND$LOGNUM.log 2>&1
            GMX_RET=$?
            break
            ;;
          -inter)
            shift
            $GMXEXE pdb2gmx $@ -ter
            GMX_RET=$?
            break
            ;;
          *)
            $GMXEXE pdb2gmx $@ > _$COMMAND$LOGNUM.log 2>&1
            GMX_RET=$?
            break
            ;;
        esac
      done
      ;;
    make_ndx)
      if [[ "$1" == "-inter" ]] ; then 
        shift
        $GMXEXE $COMMAND $@
        GMX_RET=$?
      else
        # echo $COMMAND $@ "(_$COMMAND$LOGNUM.log)"
        echo "q" | $GMXEXE $COMMAND $@ > _$COMMAND$LOGNUM.log 2>&1
        GMX_RET=$?
      fi
      ;;
    *)
      if [[ "$1" == "-inter" ]] ; then 
        shift
        $GMXEXE $COMMAND $@
        GMX_RET=$?
      else
        # echo $COMMAND $@ "(_$COMMAND$LOGNUM.log)"
        $GMXEXE $COMMAND $@ > _$COMMAND$LOGNUM.log 2>&1
        GMX_RET=$?
      fi
      ;;
  esac
  gmxRet $GMX_RET
  return 0
}

########### ANALYSIS

function dnaBondStats {
  if [[ "$1" == "-l" ]] ; then
    shift
    LOGNUM=$1
    shift
  fi
  # echo amp_stats.py --dna-hydrogen-bonds $@ 
  echo -ne $colBold"Running amp_stats.py $LOGNUM... "
  amp_stats.py --dna-hydrogen-bonds $@ > _ampStats$LOGNUM.log 2>&1
  gmxRet $?
  return 0
}

function listGroups {
  echo -e $colBold"Atom Groups in $colFile"$1$colClear$colBold":"$colClear
  echo "q" | $GMXEXE make_ndx -f $1 > _make_ndx.log 2>&1
  echo -ne $colVarName
  grep -oP --color=never "(?<=\[ ).*(?= \])" index.ndx
  echo -ne $colClear
}

function totalCharge {
  QTOT=$(grep qtot $1 | tail -n1 | grep -oP '(?<=qtot).*')
  echo -e  $colFile$2$colClear$colVarName" Total charge $colClear="$colResult"$QTOT"$colClear
  return "$QTOT"
}

function numAtomsInResidue {
  nA=$(grep $1 $2 | wc -l)
  return "$nA"
}

function groupStats {
  nSys=$(grep "Group" $1 | grep " $2)" | grep -oP "(?<=has ).*(?= elements)")
  varOut $2 $nSys atoms
}

function genionSummary {
  nRep=$(grep "solute molecules in topology file" _genion*.log | sed "s/Replacing/Replaced/")
  echo $nRep
  # echo "Replaced "$nRep" with " $(grep -P "(?= by ).*(?=ions.)" _genion*.log)

#   Processing topology
# Replacing 46 solute molecules in topology file (topol.top)  by 46 NA and 0 CL ions.
}

function numAtomsInDNA {
  numAtomsInResidue "DG" $1; nG=$?
  numAtomsInResidue "DC" $1; nC=$?
  numAtomsInResidue "DT" $1; nT=$?
  numAtomsInResidue "DA" $1; nA=$?
  let "nTot = $nG + $nC + $nT + $nA"
  return $nTot
}

function numAtomsInProtein {
  numAtomsInResidue "ALA" $1; nTot=$?
  echo $nTot
  numAtomsInResidue "CYS" $1; let "nTot = $nTot + $?"
  echo $nTot
  numAtomsInResidue "ASP" $1; let "nTot = $nTot + $?"
  echo $nTot
  numAtomsInResidue "GLU" $1; let "nTot = $nTot + $?"
  echo $nTot
  numAtomsInResidue "PHE" $1; let "nTot = $nTot + $?"
  echo $nTot
  numAtomsInResidue "GLY" $1; let "nTot = $nTot + $?"
  echo $nTot
  numAtomsInResidue "HIS" $1; let "nTot = $nTot + $?"
  echo $nTot
  numAtomsInResidue "ILE" $1; let "nTot = $nTot + $?"
  echo $nTot
  numAtomsInResidue "LYS" $1; let "nTot = $nTot + $?"
  echo $nTot
  numAtomsInResidue "LEU" $1; let "nTot = $nTot + $?"
  echo $nTot
  numAtomsInResidue "MET" $1; let "nTot = $nTot + $?"
  echo $nTot
  numAtomsInResidue "ASN" $1; let "nTot = $nTot + $?"
  echo $nTot
  numAtomsInResidue "PRO" $1; let "nTot = $nTot + $?"
  echo $nTot
  numAtomsInResidue "GLN" $1; let "nTot = $nTot + $?"
  echo $nTot
  numAtomsInResidue "ARG" $1; let "nTot = $nTot + $?"
  echo $nTot
  numAtomsInResidue "SER" $1; let "nTot = $nTot + $?"
  echo $nTot
  numAtomsInResidue "THR" $1; let "nTot = $nTot + $?"
  echo $nTot
  numAtomsInResidue "VAL" $1; let "nTot = $nTot + $?"
  echo $nTot
  numAtomsInResidue "TRP" $1; let "nTot = $nTot + $?"
  echo $nTot
  numAtomsInResidue "TYR" $1; let "nTot = $nTot + $?"
  echo $nTot
  return $nTot
}

function minimStats {
  COMMAND=$1
  LOGNUM=$2
  tail -n7 _$COMMAND$LOGNUM.log | head -n4
}

function xvg2png {
  GPSTRING=""
  TITLE_SET=0
  while test $# -gt 0; do
    case "$1" in
      -h|--help)
        echo -e "Usage for "$colFunc"xvg2png"$colClear":"

        echo -e "\n""Basic:"
        flagOut "-f|--filename" "<STRING>" "Plot <STRING>.xvg"
        flagOut "-o|--output" "<STRING>" "Produce <STRING>.png"
        flagOut "-t|--title" "<STRING>" "Graph title"
        flagOut "-xl|--xlabel" "<STRING>" "Set the x-label"
        flagOut "-yl|--ylabel" "<STRING>" "Set the y-label"
        flagOut "-xt|--xtics" "<FLOAT>" "Set the frequency of the x ticks"
        flagOut "-yt|--ytics" "<FLOAT>" "Set the frequency of the y ticks"
        flagOut "-xs|--xscientific" "" "Show the x-axis values in scientific format"
        flagOut "-ys|--yscientific" "" "Show the y-axis values in scientific format"
        flagOut "-ylog|-logy" "" "Plot the y-axis in log scale"
        flagOut "-xmin" "<FLOAT>" "Set the lower bound of the x-axis range"
        flagOut "-xmax" "<FLOAT>" "Set the upper bound of the x-axis range"
        flagOut "-ymin" "<FLOAT>" "Set the lower bound of the y-axis range"
        flagOut "-ymax" "<FLOAT>" "Set the upper bound of the y-axis range"
        flagOut "-l|--lognum" "<INTEGER>" "Set the logfile to '_gp<INTEGER>.log'"
        
        echo -e "\n""Two xvg's on the same axes:"
        flagOut "-dp|--double-plot" "<STRING>" "Produce <STRING>.png"
        flagOut "-f2|--filename2" "<STRING>" "Also plot <STRING>.xvg"
        flagOut "-t1|--title1" "<STRING>" "Title for first datafile"
        flagOut "-t2|--title2" "<STRING>" "Title for second datafile"
        
        echo -e "\n""Analysis:"
        flagOut "-ra|--running-average" "" "Running average of last five data points"
        flagOut "-cf|--constfit" "" "Fit y=c to the data"
        flagOut "-fmin" "<FLOAT>" "Set the lower bound of the fitting range"
        flagOut "-fmax" "<FLOAT>" "Set the upper bound of the fitting range"

        # finishUp 0
        return
        exit
        ;;
      -f|-f1|--filename|--filename1)
        shift
        PLOTFILE=$1
        GPSTRING=$GPSTRING"filename='$1';"
        shift
        ;;
      -f2|--filename2)
        shift
        GPSTRING=$GPSTRING"filename2='$1';"
        shift
        ;;
      -o|--output)
        shift
        GPSTRING=$GPSTRING"output='$1';"
        shift
        ;;
      -t|--title)
        shift
        GPSTRING=$GPSTRING"title='$1';"
        TITLE_SET=1
        shift
        ;;
      -t1|--title1)
        shift
        GPSTRING=$GPSTRING"title1='$1';"
        shift
        ;;
      -t2|--title2)
        shift
        GPSTRING=$GPSTRING"title2='$1';"
        shift
        ;;
      -dp|--doubleplot)
        shift
        GPSTRING=$GPSTRING"doubleplot=1;"
        ;;
      -xl|--xlabel)
        shift
        GPSTRING=$GPSTRING"xlab='$1';"
        shift
        ;;
      -yl|--ylabel)
        shift
        GPSTRING=$GPSTRING"ylab='$1';"
        shift
        ;;
      -xt|--xtics)
        shift
        GPSTRING=$GPSTRING"xtic=$1;"
        shift
        ;;
      -yt|--ytics)
        shift
        GPSTRING=$GPSTRING"ytic=$1;"
        shift
        ;;
      -xs|--xscientific)
        shift
        GPSTRING=$GPSTRING"xsci=1;"
        ;;
      -ys|--yscientific)
        shift
        GPSTRING=$GPSTRING"ysci=1;"
        ;;
      -ylog|-logy)
        shift
        GPSTRING=$GPSTRING"ylog=1;"
        ;;
      -ra|--running-average)
        shift
        GPSTRING=$GPSTRING"runavg=1;"
        ;;
      -xmin)
        shift
        GPSTRING=$GPSTRING"xmin=$1;"
        shift
        ;;
      -xmax)
        shift
        GPSTRING=$GPSTRING"xmax=$1;"
        shift
        ;;
      -ymin)
        shift
        GPSTRING=$GPSTRING"ymin=$1;"
        shift
        ;;
      -ymax)
        shift
        GPSTRING=$GPSTRING"ymax=$1;"
        shift
        ;;
      -cf|--constfit)
        shift
        GPSTRING=$GPSTRING"constfit=1;"
        ;;
      -fmin|--fitmin)
        shift
        GPSTRING=$GPSTRING"fitmin=$1;"
        shift
        ;;
      -fmax|--fitmax)
        shift
        GPSTRING=$GPSTRING"fitmax=$1;"
        shift
        ;;
      -l|--lognum)
        shift
        LOGNUM=$1
        shift
        ;;
      *)
        echo -e $colError"Unrecognised flag: "$colArg$1$colClear
        shift
        finishUp 2
        ;;
    esac
  done

  NUM_MATCH=$(grep "@ s0 legend" $PLOTFILE".xvg" | grep -oP '(?<=").*(?=")' | wc -l)
  if [ $NUM_MATCH -eq 0 ] ; then
    SELECT=$(grep " title" $PLOTFILE".xvg" | grep -oP '(?<=").*(?=")')
  else
    SELECT=$(grep "@ s0 legend" $PLOTFILE".xvg" | grep -oP '(?<=").*(?=")')
  fi

  if [ $TITLE_SET -eq 0 ] ; then
    GPSTRING=$GPSTRING"title='$SELECT';"
  fi
  # echo $GPSTRING
  echo -ne $colBold"Plotting $PLOTFILE $colArg$SELECT$colClear$colBold... "
  gnuplot -e "$GPSTRING" $GROMAX/xvg2png.gp > _gp$LOGNUM.log 2>&1
  gmxRet $?
}