# script_harm.sh for computing the band gap at 300 K, 500 K and 700 K using harmonic phonons
cd $PWD
QE=path_to_qe_bin

# copy inputs
cp ../inputs/ZG_harm.in .
cp ../inputs/Kpoints.txt .
ibrun -np 1 $QE/ZG.x < ZG_harm.in > ZG_harm.out

# run scf for equilibrium structure (no el-ph coupling)
ibrun -np 48 $QE/pw.x -nk 3 < equil-scf_222.in > equil-scf_222.out
#
# prepare file to run nscf by adding empty bands 
cp equil-scf_222.in equil-nscf_222.in
sed -i 's/scf/nscf/g' equil-nscf_222.in
sed -i 's/-nscf/-scf/g' equil-nscf_222.in
sed -i 's/ecutwfc/nosym = .true., nbnd = 140, ecutwfc/g' equil-nscf_222.in
sed -i '/K_P/d' equil-nscf_222.in
sed -i '/3   3   3/d' equil-nscf_222.in
# K point as specified in Kpoints.txt. We use Gamma point which captures the direct gap 
# of PbTe in a 2x2x2 supercell. Note that for, e.g., a 3x3x3 supercell this is not the 
# case and one needs to specify a different K point. 
cat equil-nscf_222.in Kpoints.txt > tmp
mv tmp equil-nscf_222.in
# run nscf
ibrun -n 48 $QE/pw.x < equil-nscf_222.in > equil-nscf_222.out
rm *wfc* *.save/*wfc*
#
grep highest equil-nscf_222.out > tmp
awk '{print $8-$7}' tmp > equil_gap.dat

# repeat the step above but now for T-dependent ZG files
# Those are generated relying on the harmonic phonons --> see ZG_harm.in
declare -i j
array=(300 500 700)
for j in ${array[@]}
do
  file=ZG-nscf_222_"$j".00K
  # run scf
  ibrun -n 48 $QE/pw.x -nk 8 < ZG-scf_222_"$j".00K.in > ZG-scf_222_"$j".00K.out
  cp ZG-scf_222_"$j".00K.in "$file".in
  # prepare nscf files
  sed -i 's/scf/nscf/g' "$file".in
  sed -i 's/-nscf/-scf/g' "$file".in
  sed -i 's/ecutwfc/nosym = .true., nbnd = 140, ecutwfc/g' "$file".in
  sed -i '/K_P/d' "$file".in
  sed -i '/3   3   3/d' "$file".in
  cat "$file".in Kpoints.txt > tmp
  mv tmp "$file".in
  # run nscf
  ibrun -n 48 $QE/pw.x < "$file".in > "$file".out
  rm *wfc* *.save/*wfc*
done

exit