#!/bin/sh

# run from directory where this script is
cd `echo $0 | sed 's/\(.*\)\/.*/\1/'` # extract pathname
EXAMPLE_DIR=`pwd`

# check whether echo has the -e option
if test "`echo -e`" = "-e" ; then ECHO=echo ; else ECHO="echo -e" ; fi

$ECHO
$ECHO "$EXAMPLE_DIR : starting"
$ECHO
$ECHO "This example shows how to use pw.x to calculate the total energy"
$ECHO "of silicon and of a few small molecules using hybrid functionals."

# set the needed environment variables
. ../../../environment_variables

# required executables and pseudopotentials
BIN_LIST="pw.x"
PSEUDO_LIST="Si.pz-vbc.UPF"

x_gamma_extrapolation=".TRUE."
exxdiv_treatment="gygi-baldereschi"
if [ ! -z "$1" ] ; then exxdiv_treatment="$1" ; fi
if [ "$exxdiv_treatment" = "vcut_ws" ] ; then x_gamma_extrapolation=.FALSE. ; fi
if [ "$exxdiv_treatment" = "vcut_spherical" ] ; then x_gamma_extrapolation=.FALSE. ; fi
ecutvcut=0.7

$ECHO
$ECHO "  executables directory: $BIN_DIR"
$ECHO "  pseudo directory:      $PSEUDO_DIR"
$ECHO "  temporary directory:   $TMP_DIR"
$ECHO "  checking that needed directories and files exist...\c"

# check for directories
for DIR in "$BIN_DIR" "$PSEUDO_DIR" ; do
    if test ! -d $DIR ; then 
        $ECHO
        $ECHO "ERROR: $DIR not existent or not a directory"
        $ECHO "Aborting"
        exit 1
    fi
done
for DIR in "$TMP_DIR" "$EXAMPLE_DIR/results" ; do
    if test ! -d $DIR ; then
        mkdir $DIR
    fi
done
cd $EXAMPLE_DIR/results

# check for executables
for FILE in $BIN_LIST ; do
    if test ! -x $BIN_DIR/$FILE ; then
        $ECHO
        $ECHO "ERROR: $BIN_DIR/$FILE not existent or not executable"
        $ECHO "Aborting"
        exit 1
    fi
done

# check for pseudopotentials
for FILE in $PSEUDO_LIST ; do
    if test ! -r $PSEUDO_DIR/$FILE ; then
       $ECHO
       $ECHO "Downloading $FILE to $PSEUDO_DIR...\c"
            $WGET $PSEUDO_DIR/$FILE $NETWORK_PSEUDO/$FILE 2> /dev/null
    fi
    if test $? != 0; then
        $ECHO
        $ECHO "ERROR: $PSEUDO_DIR/$FILE not existent or not readable"
        $ECHO "Aborting"
        exit 1
    fi
done
$ECHO " done"

# how to run executables
PW_COMMAND="$PARA_PREFIX $BIN_DIR/pw.x $PARA_POSTFIX"
$ECHO
$ECHO "  running pw.x as: $PW_COMMAND"
$ECHO

$ECHO
$ECHO "  running PBE0 calculation for Si with nq=1,2,4 \c"
$ECHO

for nq in 1 2 4 ; do

# self-consistent calculation
cat > si.in << EOF
 &control
    calculation = 'scf'
    restart_mode='from_scratch',
    prefix='silicon',
    pseudo_dir = '$PSEUDO_DIR/',
    outdir='$TMP_DIR/'
    tstress = .TRUE.
 /
 &system    
    ibrav=  2, celldm(1) =10.20, nat=  2, ntyp= 1,
    ecutwfc =12.0,  nbnd = 8,
    input_dft='pbe0', nqx1 = $nq, nqx2 = $nq, nqx3 = $nq, 
    exxdiv_treatment='$exxdiv_treatment'
    ecutvcut=$ecutvcut
    x_gamma_extrapolation = $x_gamma_extrapolation
 /
 &electrons
    mixing_beta = 0.7 
 /
ATOMIC_SPECIES
 Si  28.086  Si.pz-vbc.UPF
ATOMIC_POSITIONS
 Si 0.00 0.00 0.00 
 Si 0.25 0.25 0.25 
K_POINTS
  10
   0.1250000  0.1250000  0.1250000   1.00
   0.1250000  0.1250000  0.3750000   3.00
   0.1250000  0.1250000  0.6250000   3.00
   0.1250000  0.1250000  0.8750000   3.00
   0.1250000  0.3750000  0.3750000   3.00
   0.1250000  0.3750000  0.6250000   6.00
   0.1250000  0.3750000  0.8750000   6.00
   0.1250000  0.6250000  0.6250000   3.00
   0.3750000  0.3750000  0.3750000   1.00
   0.3750000  0.3750000  0.6250000   3.00
EOF
$ECHO "  running the scf calculation for Si with nq = $nq ...\c"
$PW_COMMAND < si.in > si.PBE0_nq=${nq}.out
$ECHO " done"
grep -e ! si.PBE0_nq=${nq}.out | tail -1

# clean TMP_DIR
$ECHO "  cleaning $TMP_DIR...\c"
rm -rf $TMP_DIR/silicon*
$ECHO " done"

done

$ECHO
$ECHO "  running HSE calculation for Si with nq=1,2,4 \c"
$ECHO

for nq in 1 2 4 ; do

# self-consistent calculation
cat > si.in << EOF
 &control
    calculation = 'scf'
    restart_mode='from_scratch',
    prefix='silicon',
    pseudo_dir = '$PSEUDO_DIR/',
    outdir='$TMP_DIR/'
 /
 &system    
    ibrav=  2, celldm(1) =10.20, nat=  2, ntyp= 1,
    ecutwfc =12.0,  nbnd = 8,
    input_dft='hse', nqx1 = $nq, nqx2 = $nq, nqx3 = $nq, 
    x_gamma_extrapolation = $x_gamma_extrapolation
    exxdiv_treatment = '$exxdiv_treatment'
 /
 &electrons
    mixing_beta = 0.7 
 /
ATOMIC_SPECIES
 Si  28.086  Si.pz-vbc.UPF
ATOMIC_POSITIONS
 Si 0.00 0.00 0.00 
 Si 0.25 0.25 0.25 
K_POINTS
  10
   0.1250000  0.1250000  0.1250000   1.00
   0.1250000  0.1250000  0.3750000   3.00
   0.1250000  0.1250000  0.6250000   3.00
   0.1250000  0.1250000  0.8750000   3.00
   0.1250000  0.3750000  0.3750000   3.00
   0.1250000  0.3750000  0.6250000   6.00
   0.1250000  0.3750000  0.8750000   6.00
   0.1250000  0.6250000  0.6250000   3.00
   0.3750000  0.3750000  0.3750000   1.00
   0.3750000  0.3750000  0.6250000   3.00
EOF
$ECHO "  running the scf calculation for Si with nq = $nq ...\c"
$PW_COMMAND < si.in > si.hse_nq=${nq}.out
$ECHO " done"
grep -e ! si.hse_nq=${nq}.out | tail -1

# clean TMP_DIR
$ECHO "  cleaning $TMP_DIR...\c"
rm -rf $TMP_DIR/silicon*
$ECHO " done"

done


$ECHO
$ECHO "  running Gau-PBE calculation for Si with nq=1,2,4 \c"
$ECHO

for nq in 1 2 4 ; do

# self-consistent calculation
cat > si.in << EOF
 &control
    calculation = 'scf'
    restart_mode='from_scratch',
    prefix='silicon',
    pseudo_dir = '$PSEUDO_DIR/',
    outdir='$TMP_DIR/'
 /
 &system    
    ibrav=  2, celldm(1) =10.20, nat=  2, ntyp= 1,
    ecutwfc =12.0,  nbnd = 8,
    input_dft='gaupbe', nqx1 = $nq, nqx2 = $nq, nqx3 = $nq, 
    x_gamma_extrapolation = .false.
    exxdiv_treatment = 'none'
 /
 &electrons
    mixing_beta = 0.7 
 /
ATOMIC_SPECIES
 Si  28.086  Si.pz-vbc.UPF
ATOMIC_POSITIONS
 Si 0.00 0.00 0.00 
 Si 0.25 0.25 0.25 
K_POINTS
  10
   0.1250000  0.1250000  0.1250000   1.00
   0.1250000  0.1250000  0.3750000   3.00
   0.1250000  0.1250000  0.6250000   3.00
   0.1250000  0.1250000  0.8750000   3.00
   0.1250000  0.3750000  0.3750000   3.00
   0.1250000  0.3750000  0.6250000   6.00
   0.1250000  0.3750000  0.8750000   6.00
   0.1250000  0.6250000  0.6250000   3.00
   0.3750000  0.3750000  0.3750000   1.00
   0.3750000  0.3750000  0.6250000   3.00
EOF
$ECHO "  running the scf calculation for Si with nq = $nq ...\c"
$PW_COMMAND < si.in > si.gaupbe_nq=${nq}.out
$ECHO " done"
grep -e ! si.gaupbe_nq=${nq}.out | tail -1

# clean TMP_DIR
$ECHO "  cleaning $TMP_DIR...\c"
rm -rf $TMP_DIR/silicon*
$ECHO " done"

done

$ECHO
$ECHO "  running now a few molecules with Gamma sampling ...\c"
$ECHO
PSEUDO_DIR=$EXAMPLE_DIR/Pseudo
$ECHO "  pseudo directory changed to:      $PSEUDO_DIR"
$ECHO


for xc in pbe0 hse gaupbe; do

$ECHO "  Exchange and correlation is: " $xc "...\c"
$ECHO

if [ "$xc" = "gaupbe" ] ; then x_gamma_extrapolation=.FALSE. ; fi
if [ "$xc" = "gaupbe" ] ; then exxdiv_treatment=none ; fi
if [ "$xc" = "gaupbe" ] ; then ecutvcut=0.0 ; fi

ps=1nlcc
ecut=80

cat > o.inp << EOF
&CONTROL
                 calculation = 'scf' ,
                restart_mode = 'from_scratch' ,
                      outdir = '$TMP_DIR/' ,
                  pseudo_dir = '$PSEUDO_DIR/',
                      prefix = 'o',
                     disk_io = 'minimal' ,
                      iprint = 1
                     tprnfor = .true.
 /
 &SYSTEM
                       ibrav = 1,
                   celldm(1) = 12.0,
                         nat = 1,
                        ntyp = 1,
                     ecutwfc = $ecut ,
                   input_dft = '$xc'
                       nspin = 2
                        nbnd = 4
           tot_magnetization = 2.0
            exxdiv_treatment = '$exxdiv_treatment'
                    ecutvcut = $ecutvcut
       x_gamma_extrapolation = $x_gamma_extrapolation
 /
 &ELECTRONS
                    conv_thr = 0.5d-3
 /
ATOMIC_SPECIES
  O 16.0  OPBE$ps.RRKJ3
ATOMIC_POSITIONS angstrom 
  O  0.1 0.2 0.3 
K_POINTS gamma
#automatic 
#1 1 1  0 0 0 
EOF
$ECHO "  running oxygen atom..\c"
$PW_COMMAND < o.inp > o.$xc.$ps.out-$ecut
$ECHO " done"

# clean TMP_DIR
$ECHO "  cleaning $TMP_DIR...\c"
rm -rf $TMP_DIR/o.*
$ECHO " done"

cat > c.inp << EOF
&CONTROL
                 calculation = 'scf' ,
                restart_mode = 'from_scratch' ,
                      outdir = '$TMP_DIR/' ,
                  pseudo_dir = '$PSEUDO_DIR/',
                      prefix = 'c',
                     disk_io = 'minimal' ,
                      iprint = 1
                     tprnfor = .true.
 /
 &SYSTEM
                       ibrav = 1,
                   celldm(1) = 12.0,
                         nat = 1,
                        ntyp = 1,
                     ecutwfc = $ecut ,
                        input_dft='$xc'
                       nspin = 2
                        nbnd = 4
           tot_magnetization = 2.0
            exxdiv_treatment = '$exxdiv_treatment'
                    ecutvcut = $ecutvcut
       x_gamma_extrapolation = $x_gamma_extrapolation
 /
 &ELECTRONS
                    conv_thr = 0.5d-3
 /
ATOMIC_SPECIES
  C 16.0  CPBE$ps.RRKJ3
ATOMIC_POSITIONS angstrom 
  C  0.1 0.2 0.3 
K_POINTS gamma
#automatic 
#1 1 1  0 0 0 
EOF
$ECHO "  running carbon atom..\c"
$PW_COMMAND < c.inp > c.$xc.$ps.out-$ecut
$ECHO " done"

# clean TMP_DIR
$ECHO "  cleaning $TMP_DIR...\c"
rm -rf $TMP_DIR/c.*
$ECHO " done"

cat > n.inp << EOF
&CONTROL
                 calculation = 'scf' ,
                restart_mode = 'from_scratch' ,
                      outdir = '$TMP_DIR/' ,
                  pseudo_dir = '$PSEUDO_DIR/',
                      prefix = 'n',
                     disk_io = 'minimal' ,
                      iprint = 1
                     tprnfor = .true.
 /
 &SYSTEM
                       ibrav = 1,
                   celldm(1) = 12.0,
                         nat = 1,
                        ntyp = 1,
                     ecutwfc = $ecut ,
                        input_dft='$xc'
                       nspin = 2
                        nbnd = 4
           tot_magnetization = 3.0
            exxdiv_treatment = '$exxdiv_treatment'
                    ecutvcut = $ecutvcut
       x_gamma_extrapolation = $x_gamma_extrapolation
 /
 &ELECTRONS
                    conv_thr = 0.5d-4
 /
ATOMIC_SPECIES
  N 16.0  NPBE$ps.RRKJ3
ATOMIC_POSITIONS angstrom 
  N  0.1 0.2 0.3 
K_POINTS gamma
#automatic 
#1 1 1  0 0 0 
EOF
$ECHO "  running nitrogen atom..\c"
$PW_COMMAND < n.inp > n.$xc.$ps.out-$ecut
$ECHO " done"

# clean TMP_DIR
$ECHO "  cleaning $TMP_DIR...\c"
rm -rf $TMP_DIR/n.*
$ECHO " done"

b=0.3169 
cat > n2.inp << EOF
&CONTROL
                 calculation = 'scf' ,
                restart_mode = 'from_scratch' ,
                      outdir = '$TMP_DIR/' ,
                  pseudo_dir = '$PSEUDO_DIR/',
                      prefix = 'n2',
                     disk_io = 'minimal' ,
                      iprint = 1
                     tprnfor = .true.
 /
 &SYSTEM
                       ibrav = 1,
                   celldm(1) = 12.0,
                         nat = 2,
                        ntyp = 1,
                     ecutwfc = $ecut ,
                        input_dft='$xc'
                        nbnd = 8
            exxdiv_treatment = '$exxdiv_treatment'
                    ecutvcut = $ecutvcut
       x_gamma_extrapolation = $x_gamma_extrapolation
 /
 &ELECTRONS
                    conv_thr = 1.d-4
 /
 &IONS
 /
ATOMIC_SPECIES
  N 16.0  NPBE$ps.RRKJ3
ATOMIC_POSITIONS angstrom 
  N  $b  $b  $b
  N -$b -$b -$b
K_POINTS gamma
#automatic 
#1 1 1  0 0 0 
EOF
$ECHO "  running n2 molecule..\c"
$PW_COMMAND < n2.inp > n2.$xc.$ps.out-$ecut
$ECHO " done"

# clean TMP_DIR
$ECHO "  cleaning $TMP_DIR...\c"
rm -rf $TMP_DIR/n2.*
$ECHO " done"

b=0.3256
cat > co.inp << EOF
&CONTROL
                 calculation = 'scf' ,
                restart_mode = 'from_scratch' ,
                      outdir = '$TMP_DIR/' ,
                  pseudo_dir = '$PSEUDO_DIR/',
                      prefix = 'co',
                     disk_io = 'minimal' ,
                      iprint = 1
                     tprnfor = .true.
 /
 &SYSTEM
                       ibrav = 1,
                   celldm(1) = 12.0,
                         nat = 2,
                        ntyp = 2,
                     ecutwfc = $ecut ,
                        input_dft='$xc'
                        nbnd = 8
            exxdiv_treatment = '$exxdiv_treatment'
                    ecutvcut = $ecutvcut
       x_gamma_extrapolation = $x_gamma_extrapolation
 /
 &ELECTRONS
                    conv_thr = 0.5d-3
 /
 &IONS
 /
ATOMIC_SPECIES
  C 16.0  CPBE$ps.RRKJ3
  O 16.0  OPBE$ps.RRKJ3
ATOMIC_POSITIONS angstrom 
  C  $b  $b  $b
  O -$b -$b -$b
K_POINTS gamma
#automatic 
#1 1 1  0 0 0 
EOF
$ECHO "  running co molecule..\c"
$PW_COMMAND < co.inp > co.$xc.$ps.out-$ecut
$ECHO " done" 

# clean TMP_DIR
$ECHO "  cleaning $TMP_DIR...\c"
rm -rf $TMP_DIR/n2.*
$ECHO " done"

b=0.3478
cat > o2.inp << EOF
&CONTROL
                 calculation = 'scf' ,
                restart_mode = 'from_scratch' ,
                      outdir = '$TMP_DIR/' ,
                  pseudo_dir = '$PSEUDO_DIR/',
                      prefix = 'o2',
                     disk_io = 'minimal' ,
                      iprint = 1
                     tprnfor = .true.
 /
 &SYSTEM
                       ibrav = 1,
                   celldm(1) = 12.0,
                         nat = 2,
                        ntyp = 1,
                     ecutwfc = $ecut ,
                        input_dft='$xc'
                       nspin = 2
                        nbnd = 8
           tot_magnetization = 2.0
            exxdiv_treatment = '$exxdiv_treatment'
                    ecutvcut = $ecutvcut
       x_gamma_extrapolation = $x_gamma_extrapolation
 /
 &ELECTRONS
                    conv_thr = 0.5d-3
 /
 &IONS
 /
ATOMIC_SPECIES
  O 16.0  OPBE$ps.RRKJ3
ATOMIC_POSITIONS angstrom 
  O  $b  $b  $b
  O -$b -$b -$b
K_POINTS gamma
#automatic 
#1 1 1  0 0 0 
EOF
$ECHO "  running o2 molecule..\c"
$PW_COMMAND < o2.inp > o2.$xc.$ps.out-$ecut
$ECHO " done"
# clean TMP_DIR
$ECHO "  cleaning $TMP_DIR...\c"
rm -rf $TMP_DIR/o2.*
$ECHO " done"

$ECHO
cat > summarize << EOF
grep -e ! n.$xc.$ps.out-$ecut  | tail -1 | awk '{print \$5}' > N
grep -e ! n2.$xc.$ps.out-$ecut | tail -1 | awk '{print \$5}' > N2
paste N2 N  | awk '{be= (\$1-\$2*2.0) * 13.6058 * 23.06; print "N2 : ",be}'
grep -e ! o.$xc.$ps.out-$ecut  | tail -1 | awk '{print \$5}' > O
grep -e ! o2.$xc.$ps.out-$ecut | tail -1 | awk '{print \$5}' > O2
paste O2 O  | awk '{be= (\$1-\$2*2.0) * 13.6058 * 23.06 ; print "O2 : ",be}'
grep -e ! c.$xc.$ps.out-$ecut  | tail -1 | awk '{print \$5}' > C
grep -e ! co.$xc.$ps.out-$ecut | tail -1 | awk '{print \$5}' > CO
paste CO O C | awk '{be= (\$1-\$2-\$3) * 13.6058 * 23.06; print "CO : ",be}'
rm C N O CO O2 N2
EOF
sh summarize
$ECHO 
done

$ECHO "$EXAMPLE_DIR : done"