#!/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 for GPU-accelerated molecular dynamics
module load GROMACS/2021.3-foss-2021a-CUDA-11.3.1

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



# ==============================================================================
# SECTION 1: Production MD run 
# ==============================================================================

# Prepare the production .tpr run file using:
#   - the NPT-equilibrated structure (*npt.gro)
#   - the NPT checkpoint to continue velocities (*npt.cpt)
#   - the production MD parameters (*md.mdp)
#   - the system topology and custom index file
gmx grompp -n *index.ndx -f *md.mdp \
           -c *npt.gro \
           -t *npt.cpt \
           -p "$TOP" \
           -o "${PREFIX}_md_0_10.tpr" 

# Launch the production MD simulation using the prepared .tpr file
# -v enables verbose output to monitor progress
gmx mdrun -v -deffnm "${PREFIX}_md_0_10" 


# ==============================================================================
# SECTION 2: Trajectory post-processing
# ==============================================================================

# Center the system and remove periodic boundary condition (PBC) artifacts
# from the raw trajectory before any analysis.
#   -pbc mol   : make molecules whole across periodic boundaries
#   -center    : re-center the system in the simulation box
#   -ur compact: use compact unitcell representation
# Input groups: 4 = Protein_LIG (centering reference), 0 = System (output)
#   -pbc nojump: ensures that atoms do not jump across the periodic boundary box between consecutive frames

gmx trjconv -s ${FOLDER}/*md_0_10.tpr \
            -f ${FOLDER}/*md_0_10.xtc \
            -o "${FOLDER}/${PREFIX}_nojump.xtc" \
            -pbc nojump <<< $'0'

gmx trjconv -s *md_0_10.tpr \
            -f *nojump.xtc \
            -o "${PREFIX}_traj_center_noPBC.xtc" \
            -pbc mol \
            -center \
            -n *index.ndx \
            -ur compact <<< $'4\n0' 

# Hand off to the analysis script to compute RMSD, RMSF, Rg, H-bonds, SASA and PCA
./MD_result.sh