#fraction of Fock exchange and the screening length for the system
AEXX=
HFSCREEN=

#VASP executable
VASP="srun vasp_std"

# model dielectric function 
msX="AEXX=$AEXX ; HFSCREEN=$HFSCREEN ; LMODELHF = .TRUE. "

# number of k-points
NKPTS=6





#
# shifts
#
data[1]="0.000000  0.000000  0.000000"     
data[2]="0.250000 -0.000000  0.000000"     
data[3]="0.500000 -0.000000  0.000000"     
data[4]="0.250000  0.250000  0.000000"     
data[5]="0.500000  0.250000  0.00000"    
data[6]="-0.250000  0.250000  0.00000" 
data[7]="0.500000  0.500000  0.000000" 
data[8]="-0.250000  0.500000  0.250000"
#



#
# TDHF related flags
#
TDHF="ALGO = TDHF ; PRECFOCK = Fast ; NEDOS = 2001 ; MAXMEM = 6000 ; OMEGAMAX = 20 ;  NMAXFOCKAE = 1 ; LFXC = .TRUE."

#------------------------------------------------------------------------------------
#  loop
#------------------------------------------------------------------------------------
cd $1
for i in {1..8}
do
  rm WAVECAR
  cat >KPOINTS <<!
Automatically generated mesh
       0
Gamma
 $NKPTS $NKPTS $NKPTS
   ${data[$i]}
!
  cat >INCAR.bak <<!
ISMEAR = 0 ; SIGMA = 0.02
ISYM=-1
!

#------------------------------------------------------------------------------------
#  DFT calculation
#------------------------------------------------------------------------------------
  
  echo "EDIFF = 1E-8 "  >INCAR
  echo "LOPTICS = .TRUE." >>INCAR
  echo "LPEAD = .TRUE." >>INCAR
  echo "LORBIT  = 11 "  >>INCAR
  cat INCAR.bak >>INCAR
  $VASP
  cat INCAR  OUTCAR > OUTCAR.$NKPTS.$i.PBE

  #
  # Number of unoccupied bands
  # 
  NBANDS=$(awk -v xx=$NCPU '/maximum number of plane-waves:/ { print int(($5*(240/400)^(3/2))/xx)*xx }' <OUTCAR)

#------------------------------------------------------------------------------------
#  CAM calculation
#------------------------------------------------------------------------------------
  if [ -f WAVECAR.$NKPTS.$i.CAM -a -f WAVEDER.$NKPTS.$i.CAM ]
  then 
    cp WAVECAR.$NKPTS.$i.CAM WAVECAR
    cp WAVEDER.$NKPTS.$i.CAM WAVEDER
  elif [  -f WAVECAR.$NKPTS.$i.CAM ]
  then
    cp WAVECAR.$NKPTS.$i.CAM WAVECAR
    echo "NBANDS = $NBANDS" >INCAR
    cat INCAR.bak >>INCAR
    echo "LOPTICS = .TRUE." >>INCAR
    echo "LPEAD = .TRUE." >>INCAR
    echo "LHFCALC = .TRUE. ; ALGO = Exact ; PRECFOCK = Normal ; ISIF = 0" >>INCAR
    echo "$msX" >> INCAR
    $VASP
    cat INCAR  OUTCAR > OUTCAR.$NKPTS.$i.CAM
    cp WAVECAR WAVECAR.$NKPTS.$i.CAM
    cp WAVEDER WAVEDER.$NKPTS.$i.CAM
  else
    cat INCAR.bak >INCAR
    echo "LHFCALC = .TRUE. ; ALGO = A ; TIME = 0.3 ; PRECFOCK = Normal ; ISIF = 0 ; NKRED = 2" >>INCAR
    echo "$msX" >> INCAR
    $VASP
    cat INCAR  OUTCAR > OUTCAR.$NKPTS.$i.CAM
 
    cat INCAR.bak >INCAR
    echo "LHFCALC = .TRUE. ; ALGO = A ; TIME = 0.3 ; PRECFOCK = Normal ; ISIF = 0 ; EDIFF = 1E-7" >>INCAR
    echo "$msX" >> INCAR
    $VASP
    cat INCAR  OUTCAR > OUTCAR.$NKPTS.$i.CAM

    echo "NBANDS = $NBANDS" >INCAR
    cat INCAR.bak >>INCAR
    echo "LOPTICS = .TRUE." >>INCAR
    echo "LPEAD = .TRUE." >>INCAR
    echo "LHFCALC = .TRUE. ; ALGO = Exact ; PRECFOCK = Normal ; ISIF = 0" >>INCAR
    echo "$msX" >> INCAR
    $VASP
    cat INCAR  OUTCAR > OUTCAR.$NKPTS.$i.CAM
    cp WAVECAR WAVECAR.$NKPTS.$i.CAM
    cp WAVEDER WAVEDER.$NKPTS.$i.CAM
  fi
#------------------------------------------------------------------------------------
#  CAM based TDHF calculations
#------------------------------------------------------------------------------------
  echo "NBANDS = $NBANDS" >INCAR
  cat INCAR.bak >>INCAR
  echo " $TDHF " >> INCAR
  echo "NBANDSO = 4 ; NBANDSV = 8 ;"  >>INCAR
  echo "NCSHMEM = 8 " >>INCAR
  echo "$msX" >> INCAR
  $VASP
  cat INCAR  OUTCAR > OUTCAR.$NKPTS.$NBANDS.$i.TDHF.CAM
  mv vasprun.xml  vasprun.xml.$NKPTS.$NBANDS.$i.TDHF.CAM


#------------------------------------------------------------------------------------
#  finish
#------------------------------------------------------------------------------------
  rm TWOEINT*
 rm WAVECAR WAVEDER *.tmp

  mv INCAR.bak INCAR
done

