#!/bin/bash
name="MS-PBEl-rVV10_Ag"
directory="/home/guido/work/VASP/MS.PBE-rVV10/test_Ag"

for bparam in 3.0 4.0 5.0 6.0 6.3 7.0 8.0 9.0 10.0 11.0 12.0 13.0 14.0 15.0 16.0 17.0 18.0 19.0 20.0 25.0 30.0 35.0
do

mkdir b${bparam}
cd ${directory}/b${bparam}

cat >Job<<!
#!/bin/bash
#PBS -lwalltime=168:00:00  
#PBS -S /bin/bash
#PBS -lnodes=1:ppn=16
#PBS -N vasp544-${name}

#go to directory
cd ${directory}/b${bparam}

#first the bulk
mkdir BULK
cd BULK

cp /home/guido/VASP/PseudopotentialsVASP/PAWPBE/Ag/POTCAR .
cp /home/guido/VASP/POSCARfiles/POSCAR.fcc.bulk.1x1x1.vasp_Ag POSCAR

cat >KPOINTS<<EOF
Automatic Mesh
0
Gamma
28 28 28
0. 0. 0.
EOF

cat >INCAR<<BLA
SYSTEM = Cu Bulk
PREC=Accurate
ALGO=Normal
LREAL=Auto


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


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

IBRION=2                #most accurate algorithm for minimization
ISIF=3                  #relax ions -change shape -change volume

EDIFF=1.E-7             #global break conidition SC-loop
EDIFFG=-0.0001          #relaxation will stop if all forces are smaller then ediffg
ENCUT=600             #plain wave cutoff [eV]

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

NSW = 100               #maximum number of ionic steps


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


ISMEAR=1                #MP smearing
SIGMA=0.2               

LUSE_VDW=.TRUE.
BPARAM=$bparam  

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

BLA

#run vasp
mpirun /home/guido/VASP/test/bin/vasp_std

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

echo "${bparam} \${a3D}" >> ${directory}/${name}_a3d


#-----------------------------------------------
#-----------------------------------------------
#layer
cd ${directory}/b${bparam}
mkdir LAYER_RELAX
cd LAYER_RELAX

cp /home/guido/VASP/PseudopotentialsVASP/PAWPBE/Ag/POTCAR .
cp /home/guido/VASP/scripts/poscar .
uc=1
layer=6
Dvac=16

totatom=\`echo " (\$uc * \$uc * \$layer ) + 9 " | bc -l \`
./poscar \$a3D \$uc \$layer \$Dvac nofix | head -\${totatom} > POSCAR
echo " " >> POSCAR

cat >KPOINTS<<BLI
Automatic Mesh
0
Gamma
28 28 1
0. 0. 0.
BLI

cat >INCAR<<BLO
SYSTEM = LA_REL

PREC=Accurate
ALGO=Normal               #this selects a mixture of algorithms, inlcuding ialgo38
LREAL=Auto       #projection operators

NELM=120                 #number of 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=2                #most accurate algorithm for minimization
ISIF=2                  #relax ions only 

EDIFF=1.E-7             #global break conidition SC-loop
ENCUT=600             #plain wave cutoff [eV]

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

NSW = 100               #maximum number of ionic steps

NPAR = 2                #recommended for parallelization
LMIXTAU=.TRUE.
LASPH=.TRUE.

LUSE_VDW=.TRUE.
BPARAM=${bparam}  

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

BLO

mpirun /home/guido/VASP/test/bin/vasp_std

#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 * \$layer " | 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
#Zatoms="0.0 1.0 2.0 3.0 4.0 5.0"
echo "${bparam} \${a3D} \${a3Dcalc} \$Zatoms" >> ${directory}/${name}_LAYER_RELAX

#-----------------------
#-----------------------
#VDW
cd ${directory}/b${bparam}
mkdir VDW
cd VDW

for Zvdw in 8.00 2.00 2.25 2.50 2.75 3.00 3.25 3.50 3.75 4.00 4.25 4.50 4.75 5.00 6.00 7.00
do
cd ${directory}/b${bparam}/VDW
mkdir Z\${Zvdw}
cd Z\${Zvdw}

cat /home/guido/VASP/PseudopotentialsVASP/PAWPBE/Ag/POTCAR /home/guido/VASP/PseudopotentialsVASP/PAWPBE/H/POTCAR > POTCAR
cp /home/guido/VASP/scripts/pes .

uc=2
layer=6
Dvac=16

totatom=\`echo " (\$uc * \$uc * \$layer ) + 9 + 2" | bc -l \`
./pes \$a3D \$uc \$layer \$Dvac fix \${Zatoms} Molecule 0 0 90 90 \$Zvdw 0.75 > POSCAR
echo " " >> POSCAR

cat >KPOINTS<<BLU
Automatic Mesh
0
Gamma
11 11 1
0. 0. 0.
BLU

cat >INCAR<<BLE
SYSTEM = Cu MOL
PREC=Accurate
ALGO=Normal
LREAL=Auto


NPAR=2
LMIXTAU=.TRUE.
LASPH=.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               

LUSE_VDW=.TRUE.
BPARAM=${bparam}   

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

BLE

mpirun /home/guido/VASP/test/bin/vasp_std

ES0=\`grep sigma OUTCAR | tail -1 | awk '{print \$7}'\`

echo "${bparam} \${Zvdw} \${ES0}" >> ${directory}/${name}_b${bparam}_VDW
done #Zvdw
!

qsub Job
cd ${directory}

done #bparam