#!/bin/bash
#relevant paths and variable
runname="WELL"            #run name for sub folder naming
smearing="MP"           #eg FE=FERMI MP=Methfessel Paxton TH=tetrahedron
sigma="0.2"
functional="DF2"        #functional to use: DF1=vdW-DF1 or DF2=vdW-DF2
uc=3                    #unit cell type eg. 1 2 or 3 for 1x1x1 2x2x2 or 3x3x3 bulk cell
kptype="Gamma"          #type of kpoint mesh Gamma, Monckhorst --> GAMMA for HEXAGONAL
Dvac=20                 #vacuum distance
walltime="48"           #walltime (h) 
LREALtag="Auto"         #projection scheme, can be: .TRUE. .FALSE. or Auto 
GGAtag="S2"
PARAM1="0.6265"
a3D="3.64427"		#SRPsol62.65
layerpos="0.00000 2.10402 4.20374 6.31019 8.41056 10.50654"	#SRPsol62.65
#LAYER_RELAX SUMMARY:
#summary="SRPsol75-DF2.2x2.MP.e450-s0.20.summary"

#select which input to use. WELL use Req and Zrange. ELBOW use Rrange and Zrange. BAR use Zbar and Rbar only.
WellElbow="ELBOW"        #make WELL ELBOW or BAR calculation

#Polar angles of the H2 molecule --> used for WELL or ELBOW calculations
theta=90		#90=parallel   0=perpendicular
phi=0

#sites that can be used --> The site refers to the COM coordinate of the H2 molecule 
sites="TOP"		#position in unitcell -- use HCP FCC BtH TOP TtF TtH 

#number of layers that can be used for the calculation
#WARNING 
layers="6"              #at the moment only 5 6 7 layes can be retrieved from summary file 

#Kpoint range for which to perform calculations
kpoints="13"  				#can be multiple numbers

#barrier calculations
#barriers based on Diaz et al. 10.1126/science.1178722
ZbarBtH="1.164"		; 	phiBtH=90			#BtH Barrier and phi
RbarBtH="1.032"		; 	thetaBtH=90			#dissociation above the BRIDGE site, with H pointing towards hollows

ZbarFCC="1.270"		; 	phiFCC=0			#FCC 
RbarFCC="1.5875"	; 	thetaFCC=90			#dissociation above the FCC site

ZbarHCP="1.270"		;	phiHCP=0			#HCP
RbarHCP="1.5875"	;	thetaHCP=90			#dissociation above the HCP site

ZbarTtB="1.386"		;	phiTtB=0			#Top to Bridge
RbarTtB="1.397"		;	thetaTtB=90			#Dissociation above the TOP site, with H pointing towards the Bridge.

ZbarTtH="1.27"		;	phiTtH=0			#Top to HCP
RbarTtH="1.27"		;       thetaTtH=90			#dissociation above TtH site

ZbarTtF="1.27"		;	phiTtF=120			#Top to FCC
RbarTtF="1.27"		;	thetaTtF=90			#dissociation above TtF site

#WELL  -->When WELL is selected the Rrange is Req from gasphase calculation
Zrange="10.0"
Rrelax="NO"

#when ELBOW is selected --> The grid is specified by Zrange and Rrange
#small grid ELBOW 
#Zrange="1.152 1.156 1.160 1.164 1.168 1.172 1.176"
#Rrange="1.026 1.028 1.030 1.032 1.034 1.036 1.038"
#LARge grid ELBOW
#Zrange="10.0"
Rrange="0.4 0.5 0.6 0.65 0.7 0.75 0.8 0.85 1.0 1.25 1.5 1.75 2.0 2.3"

#check if Req.dat is available, then store equilibrium bond distance
#if [ ! -f ../../GAS/Req.dat ]; then
#  echo "Req not available from gasphase calculation!"
#  exit 1
#fi
#Req=`awk '{print}' ../../GAS/Req.dat`
Req=".7458825761"	#SRPsol62.65-DF2
if [ $WellElbow == "WELL" ]; then
  Rrange=$Req				#when computing well --> use gasphase Req
elif [ $WellElbow == "BAR" ]; then
  Rrange=$Rbar
  Zrange=$Zbar
elif [ $WellElbow != "ELBOW" ]; then
  echo " ERROR"
  exit 1
fi

#BackupPath="/home/guido/VASP/BACKUP/LAYER_RELAX/"
#if [ ! -f ${BackupPath}${summary} ]; then
#  echo "ERROR summary file not present!"
#  exit 1
#fi
#a3D=`awk '{print $5}' ${BackupPath}${summary} | tail -1` 
#ecut=`awk '{print $2}' ${BackupPath}${summary} | tail -1`
ecut="450"
DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )"

#provide some checks before submitting
echo
echo "-------------------------------------------------------------------------------------"
echo "We are in folder: $DIR"
echo "path to LAYER_RELAX summary: ${BackupPath}${summary}"
echo "Calculation typ: $WellElbow "
echo "on sites: $sites"
echo
echo "Req=$Req"
echo "DVAC=$Dvac"
echo "a3D=$a3D"
echo "ecut=$ecut"
echo "KPoints=$kpoints"
echo "walltime=${walltime}"
echo "Unit Cell=${uc}"
echo "GGA=${GGAtag}"
echo "PARAM1=${PARAM1}"
echo
echo "Number of layers: $layers"
echo "LAYER_RELAX summary"
echo $layerpos
#cat ${BackupPath}${summary}
echo "--------------------------------------------------------------------------------------"
echo
echo "CORRECT?  YES/NO"
WhereYouAnIdiot=oops
read WhereYouAnIdiot
#echo " WhereYouAnIdiot = $WhereYouAnIdiot "
if [ "$WhereYouAnIdiot" != "YES" ]
then 
  echo "exiting script"
  exit 1
fi
echo

#structure of summary file
echo "#site | Z | Rin[A] | Rcalc[A] | DFx | GGA | param1 | uc | KP | ENCUT[eV] | smear | SIGMA | a3Din[A] | a3Dcalc[A] | layers | Dvac[A] | TotalE[eV] | E(s->0)[eV]| entropy[eV] | pressure[KB] | h | theta | phi | Z(A) col=lays " >> ${runname}-summary.dat
kpmax=0
counter=0
for site in $sites
do
#set Z and R according to site here when WellElbow==BAR
  if [ $WellElbow == "BAR" ]; then
    if [ $site == "BtH" ]; then
      Zrange=$ZbarBtH;		Rrange=$RbarBtH;	
      phi=$phiBtH;		theta=$thetaBtH;
    elif [ $site == "FCC" ]; then
      Zrange=$ZbarFCC;		Rrange=$RbarFCC;
      phi=$phiFCC;		theta=$thetaFCC;
    elif [ $site == "HCP" ]; then
      Zrange=$ZbarHCP;		Rrange=$RbarHCP;
      phi=$phiHCP;		theta=$thetaHCP;	
    elif [ $site == "TOP" ]; then
      Zrange=$ZbarTtB;		Rrange=$RbarTtB;
      phi=$phiTtB;		theta=$thetaTtB;
    elif [ $site == "TtH" ]; then
      Zrange=$ZbarTtH;		Rrange=$RbarTtH;
      phi=$phiTtH;		theta=$thetaTtH;
    elif [ $site == "TtF" ]; then
      Zrange=$ZbarTtF;		Rrange=$RbarTtF;
      phi=$phiTtF;		theta=$thetaTtF;
    else
      echo "error: site not supported for barrier calculation!"
      exit 1	
    fi
  fi
  for Z in $Zrange  
  do
    for R in $Rrange
    do
      for kp in $kpoints
      do
        for lay in $layers
        do
#retrieve layer spacings from summary file
#if [ $lay == 5 ]; then
#  layerpos=`awk '$7 == "'$lay'" {print}' ${BackupPath}${summary} | sort -n -k 1 | tail -1 | awk '{print $13, $14, $15, $16, $17}'`
#elif [ $lay == 6 ]; then
#  layerpos=`awk '$7 == "'$lay'" {print}' ${BackupPath}${summary} | sort -n -k 1 | tail -1 | awk '{print $13, $14, $15, $16, $17, $18}'`
#elif [ $lay == 7 ]; then
#  layerpos=`awk '$7 == "'$lay'" {print}' ${BackupPath}${summary} | sort -n -k 1 | tail -1 | awk '{print $13, $14, $15, $16, $17, $18, $19}'`
#else
#  echo "ERROR layers not available"
#  exit 1
#fi

#systematic name of the sub folders
naming="${runname}-${site}-Z${Z}-R${R}-t${theta}-p${phi}-l${lay}-k${kp}-e${ecut}-s${sigma}" #CHECK GASPHASE NAMING
mkdir ${naming}
cd ${naming}
echo ${naming}

if [ $GGAtag == "XC" ]; then
  cat >XCFUNC<<!
2
1.0 101  
1.0 130 
!
fi

if [ "$walltime" -le "8" ]
then
  que="small"
else
  que="long"
fi

#create Job file
cat >Job<<!
#!/bin/bash
#PBS -lwalltime=${walltime}:00:00  
#PBS -S /bin/bash
#PBS -lnodes=1:ppn=16
#PBS -N ${naming}
#PBS -q ${que}

# workdirectory
WORKDIR=\$PBS_O_WORKDIR
#make all nodes go to workdir
cd \${WORKDIR}

# make scratch directory
SCRATCH=/scratch/\${USER}/\${PBS_JOBID}
mkdir \${SCRATCH}
# copy files from one to other
scp -r \${WORKDIR}/* \${SCRATCH}/.

#go to scratch
cd \${SCRATCH}

#execute vasp with the generated input fules
mpiexec vasp-LIBXC

rm vasp-LIBXC

# copy files back
scp -r \${SCRATCH}/* \${WORKDIR}/.
#clean scratch
rm -r \${SCRATCH}/*

#go to workdir
cd \${WORKDIR}
rm vasp-LIBXC

#save relevant data to file summary.dat
#DO NOT FORGET escape characters!

#total energy 
E=\`tail -1 OSZICAR | awk '{print \$3}'\`
#entropy
entropy=\`grep EENTRO OUTCAR | tail -1 | awk '{print \$5}'\`
#E(s->0)
ES0=\`grep sigma OUTCAR | tail -1 | awk '{print \$7}'\`

#calculate a3D
a3D1=\`awk "NR==2" CONTCAR\`
a3D2=\`awk "NR==3" CONTCAR | awk '{print \$1}'\`
a3Dcalc=\`echo "scale=5; ( \$a3D1 * \$a3D2 * sqrt(2) / $uc ) / 1" | bc -l \`

#calculate z of atoms--------------------------------
#first the scaling factor
z3D1=\`awk "NR==2" CONTCAR\`
z3D2=\`awk "NR==5" CONTCAR | awk '{print \$3}'\`
Ztot=\`echo "( \$z3D1 * \$z3D2 ) " | bc -l\`

#then the other atoms
nr=10     #line with first atom position in CONTCAR
i=1
#each layer = 1 atom in primitive (111) cell -> totat = uc *uc *layer
totatoms=\`echo " $uc * $uc * $lay " | bc -l \`
Zatoms=\`awk "NR==10" CONTCAR | awk '{print \$3}'\`
Zold=\$Zatoms
while [ "\$i" -lt "\$totatoms" ]
do
  nr=\$((\$nr + 1))
  i=\$((\$i + 1))
  Znew=\`awk "NR==\$nr" CONTCAR | awk '{print \$3}'\`

  if [ "\$Znew" != "\$Zold" ]
  then
    Zold=\$Znew
    Znew=\`echo "scale=5; ( \$Znew * \$Ztot ) / 1" | bc -l \`
    Zatoms="\$Zatoms \$Znew"
  fi
done #while
#----------------------------------------------------

#pressure
P=\`grep -m 1 pressure OUTCAR | awk '{print \$4}'\`
#----------------------------------------------------

#calculate distance H1 - H2 --CONTCAR will be in direct coordinates
a3D=\`awk "NR==2" CONTCAR\`
z3D=\`awk "NR==5" CONTCAR | awk '{print \$3}'\`
Ztot=\`echo "( \$a3D * \$z3D ) " | bc -l\`

#the script for generating poscars places H1 and H2 at boundaries unit cell
#need to calculate R by using periodic conditions if H1y-H2y is to big 
line=\`echo " $uc * $uc * $lay + 9 + 1 " | bc -l\`
H1x=\`head -\$line CONTCAR | tail -1 | awk '{print \$1}'\`
H1y=\`head -\$line CONTCAR | tail -1 | awk '{print \$2}'\`
H1z=\`head -\$line CONTCAR | tail -1 | awk '{print \$3}'\`
line=\`echo " $uc * $uc * $lay + 9 + 2 " | bc -l\`
H2x=\`head -\$line CONTCAR | tail -1 | awk '{print \$1}'\`
H2y=\`head -\$line CONTCAR | tail -1 | awk '{print \$2}'\`
H2z=\`head -\$line CONTCAR | tail -1 | awk '{print \$3}'\`

#check if we need periodic H image to calculate R
#CONTCAR will not have negative coordinates
if (( \$(echo "sqrt((\$H1y - \$H2y)*(\$H1y - \$H2y)) > 0.5 " | bc -l ) )); then
  H1y=\`echo " \$H1y + 1.0" | bc -l\`
fi

#transform to cartesian
H1x=\`echo " (\$a3D / sqrt(2)) * \$H1x + (\$a3D / sqrt(2))*c(120*3.14159265/180)*\$H1y " | bc -l \`
H1y=\`echo " (\$a3D / sqrt(2)) * s(120*3.14159265/180)*\$H1y " | bc -l \`
H1z=\`echo " \$Ztot * \$H1z " | bc -l \`
H2x=\`echo " (\$a3D / sqrt(2)) * \$H2x + (\$a3D / sqrt(2))*c(120*3.14159265/180)*\$H2y " | bc -l \`
H2y=\`echo " (\$a3D / sqrt(2)) * s(120*3.14159265/180)*\$H2y " | bc -l \`
H2z=\`echo " \$Ztot * \$H2z " | bc -l \`
Rcalc=\`echo "scale=5; sqrt(  (\$H1x - \$H2x)*(\$H1x - \$H2x) + (\$H1y - \$H2y)*(\$H1y - \$H2y) + (\$H1z - \$H2z)*(\$H1z - \$H2z) )/1" | bc -l \`
#----------------------------------------------------

calctime=\`grep "System time (sec):" OUTCAR | awk '{print \$4}'\`
calctime=\`echo "scale=1; ((\$calctime * 1000) / 60 ) / 60 " | bc -l \`


#print results to summary
echo "$site $Z $R \$Rcalc $functional $GGAtag $PARAM1 $uc $kp $ecut $smearing $sigma $a3D \$a3Dcalc $lay $Dvac \$E \$ES0 \$entropy \$P \$calctime $theta $phi \$Zatoms " >> ../${runname}-summary.dat
!

#create kpoints file [MP for bulk / GAMMA for hexagonal]
cat >KPOINTS<<!
Automatic Mesh
0
$kptype
$kp $kp 1
0. 0. 0.
!

if [ $Rrelax == "YES" ]; then
  IBRION="1 "
else
  IBRION="-1 "
fi
#create INCAR file - DO NOT USE TABS IN INCAR FILE!!!!! 
cat >INCAR<<!
SYSTEM = Cu-LA_REL SPBXX-vdW-${functional} ${uc}x${uc} Dvac=${Dvac}

PREC=Accurate
ALGO=Fast               #this selects a mixture of algorithms, inlcuding ialgo38
LREAL=${LREALtag}       #projection operators

NELM=60                 #number of SCF steps 
NELMIN=3                #minimum SCF steps
NELMDL=-10              #number of steps H is not updated on start

ISTART=0                #start from scratch
INIWAV=1                #fill wave function with random numbers


IBRION=$IBRION
ISIF=1                 

EDIFF=1.E-5             #global break conidition SC-loop
ENCUT=$ecut             #plain wave cutoff [eV]

LCHARGE = .FALSE.       #do not write charge density dile
LWAVE = .FALSE.         #do write wave function file

NPAR = 2                #recommended for parallelization
LPLANE=.TRUE.          #recommended for parallelization
NSIM=1                  #recommended for parallelization


!
#add smearing type to INCAR
if [ "$smearing" == "FE" ]
then
  echo "#Fermi smearing " >> INCAR
  echo "ISMEAR=-1 " >> INCAR 
  echo "SIGMA=${sigma} " >> INCAR
elif [ "$smearing" == "MP" ]
then
  echo "#Methfessel Paxton smearing " >> INCAR
  echo "ISMEAR=1 " >> INCAR
  echo "SIGMA=${sigma} " >> INCAR
elif [ "$smearing" == "TH" ] 
then
  echo "#Tetrahedron smearing " >> INCAR
  echo "ISMEAR=-5 " >> INCAR
else
  echo "NO VALID SMEARING SELECTED! "
  exit 1
fi

#add DF type to INCAR
if [ "$functional" == "DF1" ]
then
  echo "#vdW-DF1 functional - Dion et al. " >> INCAR
  echo "GGA=${GGAtag}  " >> INCAR
  echo "LUSE_VDW=.TRUE. " >> INCAR
  echo "AGGAC=0.0000 " >> INCAR
  echo "PARAM1=${PARAM1}  " >> INCAR
elif [ "$functional" == "DF2" ]
then
  echo "#vdW-DF2 functional - Lee et al. " >> INCAR
  echo "GGA=${GGAtag}  " >> INCAR
  echo "LUSE_VDW=.TRUE. " >> INCAR
  echo "Zab_vdW=-1.8867 " >> INCAR
  echo "AGGAC=0.0000 " >> INCAR
  echo "PARAM1=${PARAM1}  " >> INCAR
elif [ "$functional" == "DF2.2" ]
then 
  echo "#good luck with your DF experiment! ">> INCAR
  echo "GGA=${GGAtag}  " >> INCAR
  echo "LUSE_VDW=.TRUE.  " >> INCAR
  echo "Zab_vdW=-2.2 " >> INCAR
  echo "AGGAC=0.0000 " >> INCAR
  echo "PARAM1=${PARAM1}  " >> INCAR
elif [ "$functional" == "DF1.4" ]
then
  echo "#good luck with your DF experiment! ">> INCAR
  echo "GGA=${GGAtag}  " >> INCAR
  echo "Zab_vdW=-1.4 " >> INCAR
  echo "LUSE_VDW=.TRUE.  " >> INCAR
  echo "AGGAC=0.0000 " >> INCAR
  echo "PARAM1=${PARAM1}  " >> INCAR
elif [ "$functional" == "DFoptPBE" ] 
then
  echo "#optPBE vdW " >> INCAR
  echo "GGA=BO " >> INCAR
  echo "PARAM1 = 0.1833333333 " >> INCAR
  echo "PARAM2 = 0.2200000000 " >> INCAR
  echo "LUSE_VDW = .TRUE. " >> INCAR
  echo "AGGAC=0.0000 " >> INCAR
elif [ "$functional" == "NO" ]
then
  echo "#no vdW " >> INCAR
  echo "GGA=${GGAtag} " >>INCAR
  echo "PARAM1=${PARAM1}  " >> INCAR
else
  echo "NO VALID vdW-DF selected! "
  exit 1
fi

#PAW potentials needed Copy correct POTCAR file
cat /home/guido/VASP/PseudopotentialsVASP/PAWPBE/Cu/POTCAR /home/guido/VASP/PseudopotentialsVASP/PAWPBE/H/POTCAR > POTCAR
#copy vdW-kernel file
cp /home/guido/VASP/PseudopotentialsVASP/vdw_kernel.bindat .
#copy vasp executable
cp /home/guido/VASP/VASP-LIBXC/vasp.5.3.5/vasp-LIBXC .

#create POSCAR using POSCAR script
cp /home/guido/VASP/scripts/vdW .
totatom=`echo " ($uc * $uc * $lay ) + 9 + 2" | bc -l `
./vdW $a3D $uc $lay $Dvac fix $layerpos $R $site $Z $theta $phi | head -${totatom} > POSCAR
echo " " >> POSCAR
#rm vdW

if [ $kp -gt $kpmax ]; then
  kpmax=$kp
fi

qsub Job
((counter++))
cd .. 
        done #lay
      done #kp
    done #R
  done #Z
done #site

exit 1
#retrieve highest kp
#loop over layers -> gas reference for each amount of layers
#theta and phi should not matter
#site should not matter
#cp everything except kp file -> create new KP file according to kpmax
#gasphase reference calculation
Zgas=`echo "scale=2; $Dvac / 2" | bc -l `
for lay in $layers
do 
  ((counter++))
  mkdir ${runname}-GAS-Z${Zgas}-R${Req}-t${theta}-p${phi}-l${lay}-k${kp}-e${ecut}-s${sigma}
  cd ${runname}-GAS-Z${Zgas}-R${Req}-t${theta}-p${phi}-l${lay}-k${kp}-e${ecut}-s${sigma}
  echo "${runname}-GAS-Z${Zgas}-R${Req}-t${theta}-p${phi}-l${lay}-k${kp}-e${ecut}-s${sigma}"
  cp ../${naming}/INCAR .
  cp ../${naming}/POTCAR .
  #cp ../${naming}/KPOINTS .
  cp ../${naming}/vdw_kernel.bindat .
  cp ../${naming}/vasp-LIBXC .
cat >KPOINTS<<!
Automatic Mesh
0
$kptype
$kpmax $kpmax 1
0. 0. 0.
!
  if [ $GGAtag == "XC" ]; then
    cp ../${naming}/XCFUNC .
  fi
  less ../${naming}/Job | head -4 > Job
  echo "#PBS -N GAS" >> Job
  Jobl=`wc -l ../${naming}/Job | awk '{print $1}'`
  Jobl=`echo "$Jobl -5" | bc -l `
  less ../${naming}/Job | tail -$Jobl | head -$(($Jobl-1)) >> Job
  echo "echo \"#GAS $Zgas $Req \$Rcalc $functional $GGAtag $PARAM1 $uc $kp $ecut $smearing $sigma $a3D \$a3Dcalc $lay $Dvac \$E \$ES0 \$entropy \$P \$Zatoms \" >> ../${runname}-summary.dat" >> Job
  #bottom two lines of poscar script have fixed position for H
  totatom=`echo " ($uc * $uc * $lay ) + 9 " | bc -l `
  cp /home/guido/VASP/scripts/vdW .
  ./vdW $a3D $uc $lay $Dvac fix $layerpos $Req BtH $Zgas $theta $phi | head -${totatom} > POSCAR
  ./vdW $a3D $uc $lay $Dvac fix $layerpos $Req BtH $Zgas $theta $phi | tail -2 >> POSCAR
  rm vdW
  #qsub Job
  cd ..
done
echo "total jobs = $counter "
