#!/bin/bash
#SBATCH --mail-user=icamps@gmail.com
#SBATCH --mail-type=ALL

#SBATCH --time=7-0:0
#SBATCH --job-name="BNMobPb3D"

#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=10
#SBATCH --mem=10G


module load StdEnv/2020  
source ~/bin/xtb/share/xtb/config_env.bash

export MKL_NUM_THREADS=${SLURM_CPUS_PER_TASK}
export OMP_NUM_THREADS=${SLURM_CPUS_PER_TASK},1
export OMP_STACKSIZE=4G
ulimit -s unlimited

filename1="MBNNB"
filename2="PbM4_3D"
T="300"

#level Econv/Eh Gconv/Eh·α⁻¹ Accuracy
#crude 5 × 10⁻⁴ 1 × 10⁻² 3.00
#sloppy 1 × 10⁻⁴ 6 × 10⁻³ 3.00
#loose 5 × 10⁻⁵ 4 × 10⁻³ 2.00
#lax 2 × 10⁻⁵ 2 × 10⁻³ 2.00
#normal 5 × 10⁻⁶ 1 × 10⁻³ 1.00
#tight 1 × 10⁻⁶ 8 × 10⁻⁴ 0.20
#vtight 1 × 10⁻⁷ 2 × 10⁻⁴ 0.05
#extreme 5 × 10⁻⁸ 5 × 10⁻⁵ 0.01
nOpt="extreme"

# --gfn 0
# --gfn 1
# --gfn 2
# --gfnff
hamilton="--gfn 2"

nIter="500" #default 250


# Single molecule geometry optimization
mkdir out_opt
cd out_opt

echo "Optimizing molecule 1..."
mkdir opt_mol_1
cd  opt_mol_1

cat > in_geo-opt.inp <<!
\$scc
    temp=${T}
\$write
    output file=_.out
    esp=true
    density=true
    spin population=true
    spin density=true
    mos=true
    wiberg=true
    charges=true
    mulliken=false
\$opt
    engine=rf
!

xtb ${hamilton} ./../../in_XYZ/${filename1}.xyz --input in_geo-opt.inp --molden --iterations ${nIter} --opt ${nOpt} -P ${nCPU} --namespace ${filename1} > opt_${filename1}.log
cd ..

echo "Optimizing molecule 2..."
mkdir opt_mol_2
cd  opt_mol_2

cat > in_geo-opt.inp <<!
\$scc
    temp=${T}
\$write
    output file=_.out
    esp=true
    density=true
    spin population=true
    spin density=true
    mos=true
    wiberg=true
    charges=true
    mulliken=false
\$opt
    engine=rf
!

xtb ${hamilton} ./../../in_XYZ/${filename2}.xyz --input in_geo-opt.inp --molden --iterations ${nIter} --opt ${nOpt}  -P ${nCPU} --namespace ${filename2} > opt_${filename2}.log
cd ../../

# DOCKING
mkdir out_dock
cd out_dock

cat > in_dock.inp <<!
\$dock
   pocket
   stack
   maxparent = 100
   nfinal = 10
   atm
\$end
!

echo "Docking..."

Complex=${filename1}+${filename2}
echo ${Complex}
xtb dock ./../in_XYZ/${filename1}.xyz ./../in_XYZ/${filename2}.xyz --input in_dock.inp --opt ${nOpt} --etemp ${T}> dock.log
dock_filename=dock_$Complex
mv best.xyz ${dock_filename}.xyz
cd ..

# Complex geometry optimization
mkdir out_opt-complx
cd out_opt-complx


echo "Optimizing Complex..."
opt_filename=opt_$Complex
cp ../out_dock/${dock_filename}.xyz .

cat > in_geo-opt.inp <<!
\$scc
    temp=${T}
\$write
    output file=_.out
    esp=true
    density=true
    spin population=true
    spin density=true
    mos=true
    wiberg=true
    charges=true
    mulliken=false
\$opt
    engine=rf
!

xtb ${hamilton} ${dock_filename}.xyz --input in_geo-opt.inp --molden --iterations ${nIter} --opt ${nOpt} -P ${nCPU} --namespace ${opt_filename} > opt_Complex.log
cd ..



# Molecular Dynamics
mkdir out_md
cd out_md
md_filename=md_$Complex
cp ../out_dock/${dock_filename}.xyz .

cat > in_md.inp <<!
\$md
   temp=298.15 # in K
   time= 100.0  # in ps
   dump= 50.0  # in fs
   step=  2.0  # in fs
   velo= false
   nve = true
   hmass=4
   shake=2 # constrain all bonds
   sccacc=1.0
\$end
!

echo "Molecular Dynamics..."
echo ${md_filename}
xtb $hamilton ${dock_filename}.xyz --input in_md.inp --md --namespace ${md_filename} -P ${nCPU} --iterations ${nIter} > md.log
mkdir scoord
mv *.scoord* scoord