#!/bin/bash
#SBATCH -J xxxxxx
#SBATCH --partition=gpu-prodq
#SBATCH --account=gpu_users
#SBATCH --gres=gpu:1
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=1
#SBATCH --cpus-per-task=20
#SBATCH --qos=gpu
#SBATCH -o %x-%j.out
#SBATCH -e %x-%j.err

# Load GROMACS with CUDA support
module load GROMACS/2021.3-foss-2021a-CUDA-11.3.1

# Define global variables
PREFIX=""
TOP="topol.top"



# ==============================================================================
# SECTION 1: RMSD — Root Mean Square Deviation
# Measures structural drift of the protein backbone and ligand over time
# relative to a reference structure (here the initial .tpr structure)
# ==============================================================================

# RMSD of the protein backbone (group 7 = backbone, used as both fit and analysis group)
# Output in nanoseconds (-tu ns)
echo -e "7\n7" | gmx rms -s *md_0_10.tpr \
                    -f *traj_center_noPBC.xtc \
                    -o "${PREFIX}_rmsd_prot.xvg" \
                    -tu ns 

# RMSD of the ligand (group 2 = ligand, used as both fit and analysis group)
echo -e "2\n2" | gmx rms -s *md_0_10.tpr \
                    -f *traj_center_noPBC.xtc \
                    -o "${PREFIX}_rmsd_lig.xvg" \
                    -tu ns 


# ==============================================================================
# SECTION 2: RMSF — Root Mean Square Fluctuation
# Measures per-residue atomic flexibility averaged over the trajectory.
# Group 4 = protein Cα atoms; -res averages fluctuations per residue
# ==============================================================================

echo "4" | gmx rmsf -s *md_0_10.tpr \
                    -f *traj_center_noPBC.xtc \
                    -o "${PREFIX}_rmsf.xvg" \
                    -res 


# ==============================================================================
# SECTION 3: Radius of Gyration (Rg)
# Measures protein compactness over time.
# A stable Rg indicates the protein maintains its folded conformation.
# Group 4 = protein
# ==============================================================================

echo "4" | gmx gyrate -s *md_0_10.tpr \
                      -f *traj_center_noPBC.xtc \
                      -o "${PREFIX}_rg.xvg"


# ==============================================================================
# SECTION 4: Hydrogen Bond Analysis (protein–ligand interface)
# Tracks the number, distance, and angle of H-bonds between protein and ligand.
# Group 4 = protein, Group 2 = ligand
# ==============================================================================

echo -e "4\n2" | gmx hbond -s *md_0_10.tpr \
                            -f *traj_center_noPBC.xtc \
                            -dist "${PREFIX}_hbdist.xvg" \
                            -num "${PREFIX}_hbonds.xvg" \
                            -ang "${PREFIX}_hbang.xvg" 


# ==============================================================================
# SECTION 5: Solvent Accessible Surface Area (SASA)
# Computes the surface area of the protein exposed to solvent over time.
# Group 22 = ligand (as defined in the custom index file)
# ==============================================================================

echo -e "22" | \
gmx sasa \
    -s *md_0_10.tpr \
    -f *traj_center_noPBC.xtc \
    -o "${PREFIX}_sasa.xvg" \
    -n *index.ndx \
    -tu ns


# ==============================================================================
# SECTION 6: Principal Component Analysis (PCA)
# Captures the dominant collective motions of the protein during the simulation
# ==============================================================================

# Step 6.1: Covariance matrix diagonalization
# Group 6 = backbone (fit reference), Group 4 = protein (analysis group)
# Outputs eigenvalues (-o) and eigenvectors (-v)
echo -e "6\n4" | gmx covar -s *md_0_10.tpr \
                     -f *traj_center_noPBC.xtc \
                     -o "${PREFIX}_eigenval.xvg" \
                     -v "${PREFIX}_eigenvec.trr"

# Step 6.2: Project trajectory onto the first two principal components (PC1 and PC2)
# -comp : projection of each frame onto eigenvectors
# -2d   : 2D scatter plot (PC1 vs PC2) to visualize conformational sampling
# -b 50 : skip the first 50 ps (equilibration period)
# -first 1 -last 2 : use eigenvectors 1 and 2 only
echo -e "6\n4" | gmx anaeig -s *md_0_10.tpr \
                      -f *traj_center_noPBC.xtc \
                      -v *eigenvec.trr \
                      -comp "${PREFIX}_comp.xvg" \
                      -2d "${PREFIX}_pca.xvg" \
                      -b 50 \
                      -tu ps \
                      -first 1 -last 2
# ==============================================================================
# SECTION 7: Potential Energy
# Extracts the potential energy of the system over the production trajectory.
# The .edr file is the binary energy file written by gmx mdrun.
# Output is an .xvg file that can be plotted to verify system stability.
# ==============================================================================

echo "Potential" | gmx energy \
    -f *md_0_10.edr \
    -o "${PREFIX}_potential_energy.xvg"


# ==============================================================================
# SECTION 8: Free Energy Landscape (FEL)
# Computes the Gibbs free energy landscape from the PC1/PC2 projection
# generated in Section 6 (PCA). The landscape reveals the conformational
# states sampled during the simulation and their relative stabilities.
#   -f      : 2D projection file (PC1 vs PC2) from gmx anaeig
#   -ls     : output free energy landscape in .xpm format
#   -notime : ignore the time column in the input .xvg file
# ==============================================================================

gmx sham \
    -f "${PREFIX}_pca.xvg" \
    -ls "${PREFIX}_FEL.xpm" \
    -notime