#!/bin/bash
#modified WELL-elbow barr script

#number is commandline argument
number=$1


LOG="/home/guido/LOG-PES-B86SRP68-DF2/"
if [ $number -ge 31 ]; then
  exit 1
fi
#changed site loop -->loop over the 30 configuration
#used pes program instead of vdW

location=()
#location should contain u v theta phi
location+=("Molecule 0 0 0 0")				#0 R
location+=("Molecule 0 0 90 0")				#1 R
location+=("Molecule 0 0 90 30")			#2 R
location+=("Molecule 0.5 0 0 0")			#3 R 
location+=("Molecule 0.5 0 90 0")			#4 R
location+=("Molecule 0.5 0 90 60")			#5 R
location+=("Molecule 0.5 0 90 90")			#6 R
location+=("Molecule 0.66666666 0.33333333 0 0")	#7 R
location+=("Molecule 0.66666666 0.33333333 45 30")	#8 R
location+=("Molecule 0.66666666 0.33333333 45 210")	#9 R
location+=("Molecule 0.66666666 0.33333333 90 0")	#10 R
location+=("Molecule 0.66666666 0.33333333 90 30")	#11 R
location+=("Molecule 0.33333333 0.16666666 0 0")	#12 R
location+=("Molecule 0.33333333 0.16666666 45 30")	#13 R
location+=("Molecule 0.33333333 0.16666666 45 120")	#14 R
location+=("Molecule 0.33333333 0.16666666 45 210")	#15 R
location+=("Molecule 0.33333333 0.16666666 90 30")	#16 R
location+=("Molecule 0.33333333 0.16666666 90 120")	#17 R octo
location+=("Molecule 0.33333333 -0.33333333 0 0")	#18 R octo
location+=("Molecule 0.33333333 -0.33333333 45 150")	#19 R octo 
location+=("Molecule 0.33333333 -0.33333333 45 330")	#20 R
location+=("Molecule 0.33333333 -0.33333333 90 0")	#21 R
location+=("Molecule 0.33333333 -0.33333333 90 330")	#22 R
location+=("Molecule 0.16666666 -0.16666666 0 0")	#23 R
location+=("Molecule 0.16666666 -0.16666666 45 150")	#24 R
location+=("Molecule 0.16666666 -0.16666666 45 240")	#25 R
location+=("Molecule 0.16666666 -0.16666666 45 330")	#26 R
location+=("Molecule 0.16666666 -0.16666666 90 240")	#27 R octo
location+=("Molecule 0.16666666 -0.16666666 90 330")	#28 R octo
location+=("Molecule 0 0 90 0")				#29 R	for gasphase 
location+=("Molecule 0 0 90 0")				#30 R	for extrapol

#for k in {0..28}
#do
#  echo ${location[$k]} | awk '{print $2, $3, $4, $5}'
#done 
#exit 1

#Checked positions using VESTA -->correct
#cp /home/guido/VASP/scripts/pes .
#for i in {0..9}
#do 
#  ./pes 3.6287 3 6 20 fix 0.00000 2.09503 4.18699 6.28615 8.37945 10.46748 ${location[$i]} 12 0 | tail -1
#done
#exit 1

#relevant paths and variable
runname="MOL"
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="12"           #walltime (h) 
LREALtag="Auto"         #projection scheme, can be: .TRUE. .FALSE. or Auto 
GGAtag="XC"
PARAM1="0.0"
ecut=450
a3D="3.63999" #B86SRP68-DF2
Dvac=20
#loop over layers
layers="6"              
 
layerpos="0.000000 2.10154 4.19923 6.30387 8.40225 10.49567" #B86SRP68-DF2
		
#Kpoint range for which to perform calculations
kpoints="11"  				

Zgasphase="10.0"
Rgasphase="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"

Zextrapol="10.0 9.0 8.5 8.0 7.5 7.0 6.5 6.0 5.5 5.0 4.5"
Rextrapol="0.7465449432007" #B86SRP-DF2



#ZR Cut
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"
Zrange="0.25 0.5 1.0 1.25 1.50 1.75 2.0 2.25 2.50 2.75 3.0 3.25 3.50 3.75 4.0 4.25 4.5 4.75 5.0 5.5 6.0 7.0"


Rrelax="NO"
count=0
#echo "sites= $sites"
echo "Zrange= $Zrange "
for site in $number
  do
  #mkdir site${site}
  echo ${location[$site]}
  for Z in $Zrange
  do
    for R in $Rrange
    do  
      ((count++))
    done
  done
done
echo "total amount of jobs = $count "
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 {0..30}
for site in $number  
  do
  if [ $site == "29" ]; then
    Zrange=$Zgasphase
    Rrange=$Rgasphase
  fi
  if [ $site == "30" ]; then
    Zrange=$Zextrapol
    Rrange=$Rextrapol
  fi
  for Z in $Zrange
  do
    for R in $Rrange
    do
      for kp in $kpoints
      do
        for lay in $layers
        do

#systematic name of the sub folders
naming="${runname}-site${site}-Z${Z}-R${R}-l${lay}-k${kp}-e${ecut}-s${sigma}" #CHECK GASPHASE NAMING
if [ -d "site${site}" ]; then
  cd site${site}
else  
  mkdir site${site}
  cd site${site}
fi


mkdir ${naming}
cd ${naming}
echo ${naming}

if [ $GGAtag == "XC" ]; then
  cat >XCFUNC<<!
3
0.68 41
0.32 117  
1.0 12 
!
fi

if [ $Z == "6.0" ]; then
  REALwalltime="8"
else
  REALwalltime=$walltime
fi

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


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

# 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
#cp -r \${WORKDIR}/* \${SCRATCH}/.

#go to scratch
cd \${SCRATCH}

cp ${LOG}standaard/* .
cp ${LOG}input/POSCAR.site${site}.R${R}.Z${Z} POSCAR
cp /home/guido/VASP/VASP-LIBXC/vasp.5.3.5/vasp-LIBXC .

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

rm vasp-LIBXC

cp OSZICAR ${LOG}out/OSZICAR.site${site}.R${R}.Z${Z}

# copy files back
cp -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}'\`
#----------------------------------------------------

#----------------------------------------------------

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 $functional $GGAtag $PARAM1 $uc $kp $ecut $smearing $sigma $a3D \$a3Dcalc $lay $Dvac \$E \$ES0 \$entropy \$P \$calctime \$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=-8              #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/pes .
totatom=`echo " ($uc * $uc * $lay ) + 9 + 2" | bc -l `
#  ./pes 3.6287 3 6 20 fix 0.00000 2.09503 4.18699 6.28615 8.37945 10.46748 ${location[$i]} 12 0 | tail -1
./pes $a3D $uc $lay $Dvac fix ${layerpos} ${location[$site]} $Z $R > POSCAR
echo " " >> POSCAR
rm pes

cp POSCAR ${LOG}input/POSCAR.site${site}.R${R}.Z${Z}
#cp XCFUNC ${LOG}standaard/.
#cp KPOINTS ${LOG}standaard/.
#cp POTCAR ${LOG}standaard/.
#cp INCAR ${LOG}standaard/.
#cp vdw_kernel.bindat ${LOG}standaard/.

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

echo "total amount of jobs= $counter "
echo "exiting ....."
echo
