#!/bin/bash

# This script:
# 1. Removes solvent from trajectory and outputs the new trajectory to "MD-trjconv.trr", also fixes periodic jumps.
# 2. Takes MD-trjconv.trr and removes all non-backbone atoms, and renames all backbone atoms to "CA"

# 2 is required because the coarse grain atom names are all "BB" (backbone) or SC1/SC2..etc (non-backbone).
# To visualize as tube representation or to do some analysis, VMD needs "CA" (alpha carbon). I.e. just change BB to CA.
# Don't really need HPC resource for this.

# INPUT
sysname=sym2_den
mdrun=MD313
# Trjconv stuff

module load gromacs



# Create an index file containing the indices for all protein atoms (and leave out all waters).
# Remember to change "MD.trr" to "MD.xtc" if the simulation has output .xtc file instead of .trr
gmx select -f "$mdrun.xtc" -s "$mdrun.tpr" -on protein-indices.ndx -select " (group protein) "
# Use the index file created above to select protein atoms only, fix periodic boundary jumps.
gmx trjconv -f "$mdrun.xtc" -o "$mdrun-noSolv.trr" -pbc nojump -n protein-indices.ndx

module load vmd
# Delete any pre-existing scripts
rm convert.tcl

# Write the script that loads the .pdb and trr, selects all backbone atoms and changes the backbone atom names to CA.
# We need to do this in separate instances of VMD otherwise the namechange in the middle will screw up the .trr
# Do the name change for the .pdb first.
echo "mol new" $sysname"-cg.pdb" >>convert.tcl
echo "set BB [atomselect top \"name BB\"] " >>convert.tcl
echo "\$BB set name CA" >>convert.tcl
# Write output pdb
echo "animate write pdb " $sysname"-CG-noSolv-backboneCAs.pdb sel \$BB" >>convert.tcl
echo "quit" >>convert.tcl
vmd -dispdev text -e convert.tcl

rm convert.tcl
# Do the name change for the trajectory.
echo "mol new" $sysname"-cg.pdb" >>convert.tcl
echo "mol addfile $mdrun-noSolv.trr waitfor all" >>convert.tcl
echo "set BB [atomselect top \"name BB\"] " >>convert.tcl
echo "\$BB set name CA" >>convert.tcl
# Write output trajectory
echo "animate write dcd $mdrun-CG-noSolv-backboneCAs.dcd sel \$BB" >>convert.tcl
echo "quit" >>convert.tcl
vmd -dispdev text -e convert.tcl

rm convert.tcl
# Align frames using align_frames.tcl
echo "mol new" $sysname"-CG-noSolv-backboneCAs.pdb" >>convert.tcl
echo "mol addfile $mdrun-CG-noSolv-backboneCAs.dcd waitfor all" >>convert.tcl
echo "source align_frames.tcl" >>convert.tcl
echo "align_frames top all" >>convert.tcl
echo "animate write dcd $mdrun-CG-noSolv-backboneCAs-ALIGNED.dcd" >>convert.tcl
echo "quit" >>convert.tcl
vmd -dispdev text -e convert.tcl


# get radius of gyration
printf '1' | gmx_mpi gyrate -f "$mdrun.xtc" -s "$mdrun.tpr" -o Rgyration.xvg

# Create a .tpr file with only alpha carbons and no solvent.
# Purpose: Allows input into chimera (which can allow us to cluster.). Chimera requires .tpr format. We only have .tpr of the fully solvated system so it wouldn't be compatible with our 'cleaned' MD file

printf '14\n3'| gmx_mpi pdb2gmx -f $sysname'-CG-noSolv-backboneCAs.pdb' -p $sysname'-CG-noSolv-backboneCAs.top' -i $sysname'-CG-noSolv-backboneCAs.itp' -missing
gmx_mpi editconf -f $sysname'-CG-noSolv-backboneCAs.pdb' -d 1.0 -bt cubic -o $sysname'-CG-noSolv-backboneCAs.gro'
gmx_mpi grompp -p $sysname'-CG-noSolv-backboneCAs.top' -f "$mdrun.mdp" -c $sysname'-CG-noSolv-backboneCAs.gro' -o $sysname'-CG-noSolv-backboneCAs.tpr' -maxwarn 1

rm *backboneCAs*.itp

#Extract potential energy (option no. 7) and total energy (no. 9) into energies.xvg file
printf '7 9 \n\n' | gmx_mpi energy -f "$mdrun.edr" -o energies.xvg
