#!/bin/bash
module load amber/14.13
for j in ARA
do
	for i in 1 2 3 4 5 6 
	do      
		for k in 1 2 3
		do      
			cd ~/scratch
			cd Simulation4/Docked_Sims/State_${i}
			mkdir REPEATS 
			cd REPEATS 
			mkdir ${k}
			cd ${k} 
			mkdir 001.leap
			cp -t 001.leap/ ~/scratch/Simulation4/Docked_Sims/State_${i}/001.leap/*
			mkdir 002.min
			cd 002.min
			cat > min1.in << EOF
CYP2J2_${j}_${i}_${k}: Initial Minimisation with solute fixed
&cntrl
 imin=1, irest=0, ntx=1,	 !minimise, no restart, no vels
 ntmin=1, maxcyc=10000, ncyc=3000, !steepest followed by conjugated gradient minimisation
 ntb=1,				 !PBCs
 ntr=1, restraintmask=':1-461', restraint_wt=100.0, !apply restraint mask
 cut=10.0,                        !set nonbond cutoff at 10.0A
 ntxo=1			 	 !set output to ASCII
/
EOF
			cat > min2.in << EOF
CYP2J2_${j}_${i}_${k}: Minimisation 2 with heavy atoms fixed
&cntrl
 imin=1, irest=0, ntx=1,	 !minimise, no restart, no vels
 ntmin=1, maxcyc=10000, ncyc=3000,!steepest followed by conjugated gradient minimisation
 ntb=1,				 !PBCs
 ntr=1, restraintmask=':1-461 & !@H=', restraint_wt=10.0	 !restrain heavy atoms
 cut=10.0,                        !set nonbond cutoff at 10.0A
 ntxo=1			 	 !set output to ASCII
/
EOF
			cat > min3.in << EOF
CYP2J2_${j}_${i}_${k}: Minimisation 3 with backbone fixed
&cntrl
 imin=1, irest=0, ntx=1,	 !minimise, no restart, no vels
 ntmin=1, maxcyc=10000, ncyc=3000,!steepest followed by conjugated gradient minimisation
 ntb=1,				 !PBCs
 ntr=1, restraintmask=':1-459 & @CA,C,O,N,H', restraint_wt=10.0	 !restrain backbone
 cut=10.0,                        !set nonbond cutoff at 10.0A
 ntxo=1			 	 !set output to ASCII
/
EOF
			cat > min4.in << EOF
CYP2J2_${j}_${i}_${k}: Final long minimisation with no restraints
&cntrl
 imin=1, irest=0, ntx=1,	 !minimise, no restart, no vels
 ntmin=1, maxcyc=100000, ncyc=30000,!steepest followed by conjugated gradient minimisation
 ntb=1,				 !PBCs
 ntr=0,				 !no restraints
 cut=10.0,                        !set nonbond cutoff at 10.0A
 ntxo=1			 	 !set output to netCDF
/
EOF
			cd ..
			mkdir 003.heat
			cd 003.heat
			cat > heat.in << EOF
CYP2J2_${j}_${i}_${k}: Heat while restraining solute
&cntrl
 imin=0, irest=0, ntx=1,         !No minimise, no restart, no vels
 nstlim=100000, dt=0.002,	 !100,000 steps of 0.002ps
 ntf=2, ntc=2,			 !Set SHAKE to constain H
 tempi=0.0, temp0=310.0,	 !Set initial and final temps
 ntpr=500, ntwx=500,		 !Write energy and coordinates every 500 steps
 ntb=1, ntp=0, 			 !PBCs, no pressure scaling,
 ntt=3, gamma_ln=2.0, ig=-1,	 !Langevin dynamics, collision freq, 
 nmropt=1,			 !Use varying conditions
 cut=10.0,                        !Set nonbond cutoff to 10.0A
 ntr=1, restraintmask=':1-461', restraint_wt=10.0, !Apply restraint to solute
 ntxo=2                          !Output in netCDF format
/
&wt type='TEMP0', istep1=0, istep2=90000, value1=0.0, value2=310.0 /
&wt type='TEMP0', istep1=90001, istep2=100000, value1=310.0, value2=310.0 /
&wt type='END' /
EOF
			cd ..
			mkdir 004.equil
			cd 004.equil
			cat > equil1.in << EOF
CYP2J2_${j}_${i}_${k}: Equilibration phase 1 keeping solute restrained
&cntrl
 imin=0, irest=1, ntx=5,         !No minimise, use restart file
 nstlim=50000, dt=0.002,	 !50,000 steps of 0.002ps
 ntf=2, ntc=2,			 !Set SHAKE to constain H
 temp0=310.0,			 !Set temp
 ntpr=500, ntwx=500,		 !Write energy and coordinates every 500 steps
 ntwr=10000,			 !Write restart every 10000 steps
 ioutfm=1,			 !Output mdcrd and mdvel in netCDF format
 ntb=2, ntp=1, pres0=1.0, taup=2.0,	 !PBCs auto ntb, pressure scaling,
 ntt=3, gamma_ln=2.0, ig=-1,	 !Langevin dynamics, collision freq, 
 cut=10.0,                        !Set nonbond cutoff to 10.0A
 ntr=1, restraintmask=':1-461', restraint_wt=10.0, !Apply restraint 
 ntxo=2                          !Output in netCDF format
/
EOF
			cat > equil2.in << EOF
CYP2J2_${j}_${i}_${k}: Equil phase 2 keeping backbone fixed with very weak restraints
&cntrl
 imin=0, irest=1, ntx=5,         !No minimise, use restart file
 nstlim=20000, dt=0.002,	 !20,000 steps of 0.002ps
 ntf=2, ntc=2,			 !Set SHAKE to constain H
 temp0=310.0,			 !Set temp
 ntpr=500, ntwx=500,		 !Write energy and coordinates every 500 steps
 ntwr=10000,			 !Write restart every 10000 steps
 ioutfm=1,			 !Output mdcrd and mdvel in netCDF format
 ntb=2, ntp=1, pres0=1.0, taup=2.0,	 !PBCs auto ntb, pressure scaling,
 ntt=3, gamma_ln=2.0, ig=-1,	 !Langevin dynamics, collision freq, 
 cut=10.0,                        !Set nonbond cutoff to 10.0A
 ntr=1, restraintmask=':1-459 & @CA,C,O,N,H', restraint_wt=2.0, !Apply restraint 
 ntxo=2                          !Output in netCDF format
/
EOF
			cat > equil3.in << EOF
CYP2J2_${j}_${i}_${k}: Equil phase 3 without restraints
&cntrl
 imin=0, irest=1, ntx=5,         !No minimise, use restart file
 nstlim=50000, dt=0.002,	 !50,000 steps of 0.002ps
 ntf=2, ntc=2,			 !Set SHAKE to constain H
 temp0=310.0,			 !Set temp
 ntpr=500, ntwx=500,		 !Write energy and coordinates every 500 steps
 ntwr=10000,			 !Write restart every 10000 steps
 ioutfm=1,			 !Output mdcrd and mdvel in netCDF format
 ntb=2, ntp=1, pres0=1.0, taup=2.0,	 !PBCs auto ntb, pressure scaling,
 ntt=3, gamma_ln=2.0, ig=-1,	 !Langevin dynamics, collision freq, 
 cut=10.0,                        !Set nonbond cutoff to 10.0A
 ntxo=2                          !Output in netCDF format
/
EOF
			cd ..
			mkdir 005.md
			cd 005.md
			cat > md.in << EOF
CYP2J2_${j}_${i}_${k}: Production Phase
&cntrl
 imin=0, irest=1, ntx=5,         !No minimise, restart, vels
 nstlim=25000000, dt=0.002,	 !25,000,000 steps of 0.002ps
 ntf=2, ntc=2,			 !Set SHAKE to constain H
 temp0=310.0,			 !Set temp
 ntpr=10000, ntwx=10000,			 !Write energy and coordinates every 5000 stepsTHIS IS DIFFERENT HERE
 ntwr=100000,			 !Write restart every 100000 steps
 ioutfm=1,			 !Output mdcrd and mdvel in netCDF format
 ntb=2, ntp=1, pres0=1.0, taup=1.0,	 !PBCs auto ntb, pressure scaling
 ntt=3, gamma_ln=2.0, ig=-1,	 !Langevin dynamics, collision freq, 
 cut=10.0,                       !Set nonbond cutoff to 10.0A
 ntxo=2                          !Output in netCDF format
/
EOF
			cd ..
			mkdir 006.analysis
			mkdir 007.cpptraj
			mkdir 008.energy
			cd 008.energy
			cat > mmgbsa.in << EOF
Input file for running GB starting after 10ns (Frame 5000)
&general
	startframe=5000,
	endframe=25000, keep_files=2,
	interval=500, netcdf=1,
/
&gb
	igb=5, saltcon=0.100,
/
EOF
			cd ../007.cpptraj
			cat > strip.${i}.${k}.in << EOF
parm ../005.md/Sim4_min3_dock2_box.${i}.prmtop
trajin ../005.md/md.mdcrd 1 last 
autoimage
strip :WAT,Na+,Cl- nobox outprefix strip
rms first :1-460&!@H= mass
trajout strip.md.nc nobox
run
EOF
			cat > cluster_post_sm.${i}.${k}.in << EOF
parm strip.Sim4_min3_dock2_box.${i}.prmtop
trajin strip.md.nc
loadtraj name MYTRAJ
cluster crdset MYTRAJ :1-461&!@H= mass clusters 5 out cluster_${i}.${k}_out nofit averagelinkage summary summary_${i}.${k}_out info cluster_${i}.${k}_info repout State_${i}.${k}.rep repfmt pdb clusterout cluster.State_${i}.${k}.nc clusterfmt netcdf
run
EOF
			cat > rmsf_${j}.${i}.${k}.in <<EOF
parm strip.Sim4_min3_dock2_box.${i}.prmtop
trajin strip.md.nc 
atomicfluct out rmsf_${j}.${i}.${k}.agr :461&!@H=
run
EOF
			cat > rmsf.${i}.${k}.in << EOF
parm strip.Sim4_min3_dock2_box.${i}.prmtop
trajin strip.md.nc 
atomicfluct out rmsf.${i}.${k}.agr byres
run
EOF
			cat > rmsd.${j}.${i}.${k}.in << EOF
parm strip.Sim4_min3_dock2_box.${i}.prmtop
trajin strip.md.nc 
reference strip.md.nc 1 [frame1]
rms ${j}_nofit_${i}.${k} :461 reference out rmsd.${j}.nofit.${i}.${k}.agr nofit
run
EOF
			cat > rmsd.${i}.${k}.in << EOF
parm strip.Sim4_min3_dock2_box.${i}.prmtop
trajin strip.md.nc 
reference strip.md.nc 1 [frame1]
rms rmsd_${i}.${k} :1-460&!@H= reference out rmsd.protein.${i}.${k}.agr mass
run
EOF

			cd ..
			cat > amber_sub.sh << EOF
#!/bin/bash

# Job Name
#BSUB -J Sim4_PUFA_dock_${j}_${i}_${k}

# Hosts to use (3GPU or 8GPU); Leave blank for k80

# Queue Name; use emerald, emerald-short or emerald-k80
#BSUB -q emerald-k80

# Asks for exclusive host (not required unless doing MPI)

# Forces job to re-q and rerun from scratch if host fails
#BSUB -r

# Sends email when the job is done
#BSUB -N
#BSUB -u k.abelak@ucl.ac.uk

# Highest of the number of GPUs or CPUs needed
#BSUB -n 1

# Number of GPUs needed
#BSUB -R "span[ptile=1]"

# Output log and error files
#BSUB -o /home/ucl/eisuc224/job_logs/Sim4_PUFA_dock_${j}_${i}_${k}.%J.log
#BSUB -e /home/ucl/eisuc224/job_logs/Sim4_PUFA_dock_${j}_${i}_${k}.%J.err

# Wallclock time requested
#BSUB -W 168:00
#----------------------------------

# Load modules and change to appropriate scratch directory
cd /home/ucl/eisuc224/scratch/
module load amber/14.13
source /apps/amber/amber14u13/amber14/amber.sh

# Job commands start here; also, can use -i <input file> for commands
cd Simulation4/Docked_Sims/State_${i}/REPEATS/${k}
cd 001.leap
cp -t ../002.min Sim4_min3_dock2_box.${i}.prmtop Sim4_min3_dock2_box.${i}.rst7

cd ../002.min

$AMBERHOME/bin/pmemd.cuda -O -i min1.in -o min1.out -p Sim4_min3_dock2_box.${i}.prmtop -c Sim4_min3_dock2_box.${i}.rst7 -r min1.rst7 -ref Sim4_min3_dock2_box.${i}.rst7 -inf ../mdinfo
$AMBERHOME/bin/pmemd.cuda -O -i min2.in -o min2.out -p Sim4_min3_dock2_box.${i}.prmtop -c min1.rst7 -r min2.rst7 -ref min1.rst7 -inf ../mdinfo
$AMBERHOME/bin/pmemd.cuda -O -i min3.in -o min3.out -p Sim4_min3_dock2_box.${i}.prmtop -c min2.rst7 -r min3.rst7 -ref min2.rst7 -inf ../mdinfo
$AMBERHOME/bin/pmemd.cuda -O -i min4.in -o min4.out -p Sim4_min3_dock2_box.${i}.prmtop -c min3.rst7 -r min4.rst7 -inf ../mdinfo

cp -t ../003.heat Sim4_min3_dock2_box.${i}.prmtop min4.rst7
cd ../003.heat

$AMBERHOME/bin/pmemd.cuda -O -i heat.in -o heat.out -p Sim4_min3_dock2_box.${i}.prmtop -c min4.rst7 -r heat.rst7 -x heat.mdcrd -ref min4.rst7 -inf ../mdinfo

cp -t ../004.equil Sim4_min3_dock2_box.${i}.prmtop heat.rst7
cd ../004.equil

$AMBERHOME/bin/pmemd.cuda -O -i equil1.in -o equil1.out -p Sim4_min3_dock2_box.${i}.prmtop -c heat.rst7 -r equil1.rst7 -x equil1.mdcrd -ref heat.rst7 -inf ../mdinfo
$AMBERHOME/bin/pmemd.cuda -O -i equil2.in -o equil2.out -p Sim4_min3_dock2_box.${i}.prmtop -c equil1.rst7 -r equil2.rst7 -x equil2.mdcrd -ref equil1.rst7 -inf ../mdinfo
$AMBERHOME/bin/pmemd.cuda -O -i equil3.in -o equil3.out -p Sim4_min3_dock2_box.${i}.prmtop -c equil2.rst7 -r equil3.rst7 -x equil3.mdcrd -inf ../mdinfo

cp -t ../005.md Sim4_min3_dock2_box.${i}.prmtop equil3.rst7
cd ../005.md

$AMBERHOME/bin/pmemd.cuda -O -i md.in -o md.out -p Sim4_min3_dock2_box.${i}.prmtop -c equil3.rst7 -r md.rst7 -x md.mdcrd -inf ../mdinfo

cd ..
cd 006.analysis

$AMBERHOME/AmberTools/src/etc/process_mdout.perl ../003.heat/heat.out ../004.equil/equil1.out ../004.equil/equil2.out ../004.equil/equil3.out
sed -i -e 1,201d < summary.DENSITY
$AMBERHOME/bin/ambpdb -p ../002.min/Sim4_min3_dock2_box.${i}.prmtop -bres < ../002.min/min4.rst7 > minimised.Sim4_min3_dock2.${i}.${j}.${k}.pdb

cd ../007.cpptraj
$AMBERHOME/bin/cpptraj -i strip.${i}.${k}.in
$AMBERHOME/bin/cpptraj -i rmsf_${j}.${i}.${k}.in
$AMBERHOME/bin/cpptraj -i rmsf.${i}.${k}.in
$AMBERHOME/bin/cpptraj -i rmsd.${j}.${i}.${k}.in
$AMBERHOME/bin/cpptraj -i rmsd.${i}.${k}.in
#$AMBERHOME/bin/cpptraj -i cluster_post_sm.${i}.${k}.in

cd ../008.energy
#/apps/amber/amber14u13/amber14/bin/MMPBSA.py -O -i mmgbsa.in -o MMGBSA.${i}.${k}.dat -lp ../001.leap/${j}.${i}.${k}.prmtop -cp ../001.leap/C2J2_${j}_docked_gas_${i}.${k}.prmtop -rp ../001.leap/C2J2_${j}_receptor_${i}.${k}_gas.prmtop -sp ../001.leap/Sim4_min3_dock2_box.${i}.prmtop -y ../005.md/md.mdcrd > log.log 2>&1


EOF
			bsub < amber_sub.sh
		done
	done
done

