#!/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
location+=("Molecule 0 0 90 0")                         #31 R   for gasphase
location+=("Molecule 0 0 90 0")                         #32 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"
uc=2                    #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=15                 #vacuum distance
walltime="24"           #walltime (h) 
a3D="3.58005" #MS6.PBE
#loop over layers
layers="6"              
 
layerpos="0.000000 2.06694 4.10673 6.16382 8.21068 10.25783" #MS6.PBE
		
#Kpoint range for which to perform calculations
kpoints="11"  				

Zgasphase="8.0"
Rgasphase="0.4 0.5 0.6 0.65 0.7 0.725 0.75 0.775 0.8 0.85 1.0 1.25 1.5 1.75 2.0 2.3"

Zextrapol="8.0 7.5 7.0 6.5 6.0 5.5 5.0 4.75 4.5 4.25 4.0 3.75 3.5 3.25 3.0 2.75 2.5 2.0"
Rextrapol="0.742924" #MS6.RPBE



#ZR Cut
Rrange="0.4 0.5 0.6 0.65 0.7 0.75 0.8 0.85 0.95 1.05 1.15 1.25 1.5 1.75 2.0 2.3"
Zrange="0.25 0.5 0.75 1.0 1.125 1.25 1.375 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 5.0 5.75"

#Rrange="0.4 0.5 0.6 0.65 0.7 0.75 0.8 0.85 1.0 1.1 1.2 1.25 1.3 1.5 1.75 2.0 2.3"
#Zrange="0.25 0.5 0.75 0.95 1.0 1.05 1.1 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 5.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

if [ ! -d "${naming}" ]; then
mkdir ${naming}
cd ${naming}
echo ${naming}


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}
#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}


#execute vasp with the generated input fules

mpirun /home/guido/VASP/test/bin/vasp_std
rm CHG* EIGENVAL IBZKPT PCDAT REPORT WAVECAR XDATCAR

# 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.
!

#create INCAR file - DO NOT USE TABS IN INCAR FILE!!!!! 
cat >INCAR<<!
SYSTEM = Cu BAR SRMGGA MS6 RPBE
PREC=Accurate
ALGO=VeryFast
LREAL=Auto


NPAR=2
LMIXTAU=.TRUE.
LASPH=.TRUE.
ADDGRID=.TRUE.

NELM=120                 #number of SCF steps
NELMIN=4                #minimum SCF steps
NELMDL=-10              #number of steps H is not updated on start
NSW = 1               #maximum number of ionic steps

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

IBRION=-1                #most accurate algorithm for minimization
ISIF=0                  #relax ions -change shape -change volume

EDIFF=1.E-6             #global break conidition SC-loop
ENCUT=600

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

ISMEAR=1                #MP smearing
"SIGMA=0.2               

METAGGA=SRMGGA
SRVAL=1.0
FINTB=1.0 
SRAC=1.0

SRTA0=0
SRA0MU=0.12345679012345678
SRA0KP=0.804
SRA0C=0.1036199

SRTA1=0
SRA1MU=0.12345679012345678
SRA1KP=0.804

SRTB0=0
SRB0MU=0.12345679012345678
SRB0KP=0.804
SRB0C=0.1036199

SRTB1=0
SRB1MU=0.12345679012345678
SRB1KP=0.804

!
#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 .. 
fi #skip if folder naming exists
cd ..
        done #lay
      done #kp
    done #R
  done #Z
done #site

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