#!/bin/bash
#SBATCH --job-name=docking
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --output=output.log 
#SBATCH --partition=long30dq
#SBATCH --cpus-per-task=30

# Initialize Conda so that 'conda activate' works in non-interactive shell
eval "$(conda shell.bash hook)"

# Activate the conda environment that contains GNU parallel
conda activate parallel_env

# Load the AutoDock Vina module from the HPC module system
module load autodock-vina-1_1_2-gcc-10.2.0-2wyrl2q

# Create the top-level directory where all docking results will be stored
RESULTS_FOLDER="docking_resultats"
mkdir -p $RESULTS_FOLDER

# Define a function that runs Vina for a single ligand.
# This function will be exported so GNU parallel can call it in subshells.
dock_ligand() {
    ligand_file="$1"

    # Extract the ligand name by removing the path and the .pdbqt extension
    ligand_name=$(basename "$ligand_file" .pdbqt)

    echo "Processing ligand $ligand_name"

    # Create a dedicated output subdirectory for this ligand
    ligand_result_folder="$RESULTS_FOLDER/$ligand_name"
    mkdir -p "$ligand_result_folder"

    # Run AutoDock Vina using the shared config file (box coordinates, receptor, etc.)
    # Output poses are saved to out.pdbqt and the scoring log to log.txt
    vina --config config.txt \
         --ligand "$ligand_file" \
         --cpu 30 \
         --out "$ligand_result_folder/out.pdbqt" > "$ligand_result_folder/log.txt"
}

# Export the function so it is available to GNU parallel subprocesses
export -f dock_ligand

# Find all .pdbqt ligand files and run docking in parallel using 11 CPU cores
find Ligands_PDBQT/ -name "*.pdbqt" | parallel -j 11 dock_ligand {}

# Check whether any output files were generated before trying to parse them
if ls $RESULTS_FOLDER/*/out.pdbqt 1> /dev/null 2>&1; then

    # Write CSV header for the scores file
    echo "Ligand,Meilleur Score (kcal/mol)" > $RESULTS_FOLDER/docking_scores.csv

    # Loop over all output files and extract the best (lowest) binding affinity score
    for result in $RESULTS_FOLDER/*/out.pdbqt; do

        # Retrieve the ligand name from the parent directory name
        ligand_name=$(basename "$(dirname "$result")")

        # Parse the VINA RESULT line, extract the affinity column, sort numerically, take the best score
        best_score=$(grep "^REMARK VINA RESULT" "$result" | awk '{print $4}' | sort -n | head -1)

        # Keep only ligands with a binding affinity below -8.0 kcal/mol (strong binders)
        if (( $(echo "$best_score < -8.0" | bc -l) )); then
            echo "$ligand_name,$best_score" >> $RESULTS_FOLDER/docking_scores.csv
        fi
    done
fi