#!/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 a calculation of STM maps."

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

# required executables and pseudopotentials
BIN_LIST="pw.x pp.x plotrho.x projwfc.x sumpdos.x"
PSEUDO_LIST="Al.pz-vbc.UPF As.pz-bhs.UPF"

$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 gnuplot
GP_COMMAND=`which gnuplot 2>/dev/null`
if [ "$GP_COMMAND" = "" ]; then
        $ECHO
        $ECHO "gnuplot not in PATH"
        $ECHO "Results will not be plotted"
fi

# 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"
PP_COMMAND="$PARA_PREFIX $BIN_DIR/pp.x $PARA_POSTFIX"
PLOTRHO_COMMAND="$BIN_DIR/plotrho.x"
PROJWFC_COMMAND="$PARA_PREFIX $BIN_DIR/projwfc.x $PARA_POSTFIX"
SUMPDOS_COMMAND="$BIN_DIR/sumpdos.x"
$ECHO
$ECHO "  running pw.x as:      $PW_COMMAND"
$ECHO "  running pp.x as:      $PP_COMMAND"
$ECHO "  running plotrho.x as: $PLOTRHO_COMMAND"
$ECHO "  running projwfc.x as: $PROJWFC_COMMAND"
$ECHO "  running sumpdos.x as: $SUMPDOS_COMMAND"
$ECHO "  running gnuplot as:   $GP_COMMAND"
$ECHO

# self-consistent calculation
cat > AlAs110re.scf.in << EOF
 &control
    calculation = 'scf'
    restart_mode='from_scratch',
    title='AlAs 110 surface slab, relaxed (central plane fixed)'
    pseudo_dir = '$PSEUDO_DIR/',
    outdir='$TMP_DIR/',
    prefix='AlAs110'
 /
 &system
    ibrav=  8, celldm(1) =7.424621202, celldm(2)=1.414213576,
    celldm(3)= 6.00000,
    nat= 14, ntyp= 2,
    ecutwfc =14.0,
 /
 &electrons
    mixing_mode = 'plain'
    mixing_beta = 0.7
    conv_thr =  1.0d-6
 /
 ATOMIC_SPECIES
 Al  26.98  Al.pz-vbc.UPF
 As  74.92  As.pz-bhs.UPF
ATOMIC_POSITIONS alat
 As  0.000000000  -0.044777195  -0.058158722
 Al  0.500000000   0.251460333   0.113525468
 As  0.500000000   0.712279009   0.504183885
 Al  0.000000000   1.067633546   0.480460620
 As  0.000000000  -0.003937059   0.995826731
 Al  0.500000000   0.351400965   1.004220212
 As  0.000000000  -0.003937059   2.004173269
 Al  0.500000000   0.351400965   1.995779788
 As  0.500000000   0.712279009   2.495816115
 Al  0.000000000   1.067633546   2.519539380
 As  0.000000000  -0.044777195   3.058158722
 Al  0.500000000   0.251460333   2.886474532
 As  0.500000000   0.707106800   1.500000000
 Al  0.000000000   1.060660200   1.500000000
K_POINTS {automatic}
 6 2 1 0 0 0
EOF
$ECHO "  running the scf calculation...\c"
$PW_COMMAND < AlAs110re.scf.in > AlAs110re.scf.out
check_failure $?
$ECHO " done"

cat > AlAs110re.nonscf.in << EOF
 &control
    calculation = 'nscf'
    restart_mode='from_scratch',
    title='AlAs 110 surface slab, relaxed (central plane fixed)'
    pseudo_dir = '$PSEUDO_DIR/',
    outdir='$TMP_DIR/',
    prefix='AlAs110'
 /
 &system
    ibrav=  8, celldm(1) =7.424621202, celldm(2)=1.414213576,
    celldm(3)= 6.00000,
    nat= 14, ntyp= 2,
    occupations='smearing', smearing='gaussian', degauss=0.01, nbnd=34,
    ecutwfc =14.0,
 /
 &electrons
    mixing_mode = 'plain'
    mixing_beta = 0.7
    conv_thr =  1.0d-6
 /
 ATOMIC_SPECIES
 Al  26.98  Al.pz-vbc.UPF
 As  74.92  As.pz-bhs.UPF
ATOMIC_POSITIONS alat
 As  0.000000000  -0.044777195  -0.058158722
 Al  0.500000000   0.251460333   0.113525468
 As  0.500000000   0.712279009   0.504183885
 Al  0.000000000   1.067633546   0.480460620
 As  0.000000000  -0.003937059   0.995826731
 Al  0.500000000   0.351400965   1.004220212
 As  0.000000000  -0.003937059   2.004173269
 Al  0.500000000   0.351400965   1.995779788
 As  0.500000000   0.712279009   2.495816115
 Al  0.000000000   1.067633546   2.519539380
 As  0.000000000  -0.044777195   3.058158722
 Al  0.500000000   0.251460333   2.886474532
 As  0.500000000   0.707106800   1.500000000
 Al  0.000000000   1.060660200   1.500000000
K_POINTS {automatic}
 12 4 1 0 0 0
EOF
$ECHO "  running the non-scf calculation...\c"
$PW_COMMAND < AlAs110re.nonscf.in > AlAs110re.nonscf.out
check_failure $?
$ECHO " done"

# post-processing for stm images (sample bias given in Ry!)
cat > AlAs110.pp_stm-.in << EOF
 &inputpp
    prefix  = 'AlAs110'
    outdir='$TMP_DIR/',
    filplot = 'AlAsresm-1.0'
    sample_bias=-0.0735d0,
    plot_num= 5
 /
 &plot
   nfile=1
   filepp(1)='AlAsresm-1.0'
   weight(1)=1.0
   iflag=2
   output_format=2
   e1(1)=7.0, e1(2)=0.0,     e1(3)=0.0
   e2(1)=0.0, e2(2)=7.07107, e2(3)=0.0
   x0(1)=0.0, x0(2)=-0.18,   x0(3)=3.25
   nx=36 ,ny=56
   fileout='AlAs110-1.0'
 /
EOF
$ECHO
$ECHO "  running the post-processing phase, negative bias...\c"
$PP_COMMAND < AlAs110.pp_stm-.in > AlAs110.pp_stm-.out
check_failure $?
$ECHO " done"

# run plotrho to do the figure
cat > AlAs110.plotrho-.in << EOF
AlAs110-1.0
AlAs110-1.0eV.ps
n
0.00005 0.0078  8
EOF
$ECHO "  running plotrho on negative bias data...\c"
$PLOTRHO_COMMAND < AlAs110.plotrho-.in > AlAs110.plotrho-.out
check_failure $?
$ECHO " done"

# post-processing for stm images (negative bias, constant current)
cat > AlAs110.pp_isostm-.in << EOF
 &inputpp
 /
 &plot
    nfile=1
    filepp(1)='AlAsresm-1.0'
    weight(1)=1.0
    iflag=2
    output_format=7    
    fileout='AlAs110.pp_isostm-.dat'
    e1(1)=7.0, e1(2)=0.0,     e1(3)=0.0
    e2(1)=0.0, e2(2)=7.07107, e2(3)=0.0
    nx=150, ny=150
    isostm_flag=.true.
    isovalue=0.00005
    heightmin=0.4
    heightmax=0.75
    direction=1
 /
EOF
$ECHO
$ECHO "  STM image, negative bias and constant current...\c"
$PP_COMMAND < AlAs110.pp_isostm-.in > AlAs110.pp_isostm-.out
check_failure $?
$ECHO " done"

# run gnuplot to do the figure
if [ "$GP_COMMAND" = "" ]; then
    break
else
cat > gnuplot.tmp <<EOF
set term postscript enhanced color solid lw 3 24
set output 'AlAs110-1.0eV.isoplot.ps'
set xlabel "x (bohr)"
set ylabel "y (bohr)"
set pm3d map
set size ratio -1
set palette rgb 21,22,23
set tics out
unset key
splot [0:51.972][0:52.500] 'AlAs110.pp_isostm-.dat'
EOF
$ECHO
$ECHO "  plotting results ...\c"
$GP_COMMAND < gnuplot.tmp
$ECHO " done"
rm gnuplot.tmp
fi

# post-processing for stm images (as before, but positive bias)
cat > AlAs110.pp_stm+.in << EOF
 &inputpp
    prefix  = 'AlAs110'
    outdir='$TMP_DIR/',
    filplot = 'AlAsresm+1.0'
    sample_bias=0.0735d0,
    plot_num= 5
 /
 &plot
   nfile=1
   filepp(1)='AlAsresm+1.0'
   weight(1)=1.0
   iflag=2
   output_format=2
   e1(1)=7.0, e1(2)=0.0,     e1(3)=0.0
   e2(1)=0.0, e2(2)=7.07107, e2(3)=0.0
   x0(1)=0.0, x0(2)=-0.18,   x0(3)=3.25
   nx=36 ,ny=56
   fileout='AlAs110+1.0'
 /
EOF
$ECHO "  running the post-processing phase, positive bias...\c"
$PP_COMMAND < AlAs110.pp_stm+.in > AlAs110.pp_stm+.out
check_failure $?
$ECHO " done"

# plotrho
cat > AlAs110.plotrho+.in << EOF
AlAs110+1.0
AlAs110+1.0eV.ps
n
0.00002 0.0021  8
EOF
$ECHO "  running plotrho on positive bias data...\c"
$PLOTRHO_COMMAND < AlAs110.plotrho+.in > AlAs110.plotrho+.out
check_failure $?
$ECHO " done"

# post-processing for stm images (positive bias, constant current)
cat > AlAs110.pp_isostm+.in << EOF
 &inputpp
 /
 &plot
    nfile=1
    filepp(1)='AlAsresm+1.0'
    weight(1)=1.0
    iflag=2
    output_format=7    
    fileout='AlAs110.pp_isostm+.dat'
    e1(1)=7.0, e1(2)=0.0,     e1(3)=0.0
    e2(1)=0.0, e2(2)=7.07107, e2(3)=0.0
    nx=150, ny=150
    isostm_flag=.true.
    isovalue=0.00005
    heightmin=0.4
    heightmax=0.75
    direction=1
 /
EOF
$ECHO
$ECHO "  STM image, positive bias and constant current...\c"
$PP_COMMAND < AlAs110.pp_isostm+.in > AlAs110.pp_isostm+.out
check_failure $?
$ECHO " done"

# run gnuplot to do the figure
if [ "$GP_COMMAND" = "" ]; then
    break
else
cat > gnuplot.tmp <<EOF
set term postscript enhanced color solid lw 3 24
set output 'AlAs110+1.0eV.isoplot.ps'
set xlabel "x (bohr)"
set ylabel "y (bohr)"
set pm3d map
set size ratio -1
set palette rgb 21,22,23
set tics out
unset key
splot [0:51.972][0:52.500] 'AlAs110.pp_isostm+.dat'
EOF
$ECHO
$ECHO "  plotting results ...\c"
$GP_COMMAND < gnuplot.tmp
$ECHO " done"
rm gnuplot.tmp
fi

# Projection of the DOS on volumes (boxes)
cat > AlAs110.box.projwfc.in << EOF
 &projwfc
    prefix  = 'AlAs110'
    outdir='$TMP_DIR/',
    ngauss=0
    degauss=0.01
    DeltaE=0.02
    tdosinboxes=.true.
    plotboxes=.true.
    n_proj_boxes=8

!! Boxes centered on the first vacuum layer:
  !! 1) above the surface Al
    irmin(1,1)= 0, irmax(1,1)= 2, irmin(2,1)= 0, irmax(2,1)= 2, irmin(3,1)=63, irmax(3,1)=65,
  !! 2) above the surface As
    irmin(1,2)= 9, irmax(1,2)=11, irmin(2,2)= 5, irmax(2,2)= 7, irmin(3,2)=63, irmax(3,2)=65,
  !! 3) above the 2nd layer Al
    irmin(1,3)= 9, irmax(1,3)=11, irmin(2,3)=14, irmax(2,3)=16, irmin(3,3)=63, irmax(3,3)=65,
  !! 4) as large as the surface unit cell
    irmin(1,4)= 1, irmax(1,4)=18, irmin(2,4)= 1, irmax(2,4)=27, irmin(3,4)=63, irmax(3,4)=65,

!! Same as above, centered on the second vacuum layer:
    irmin(1,5)= 0, irmax(1,5)= 2, irmin(2,5)= 0, irmax(2,5)= 2, irmin(3,5)=72, irmax(3,5)=74,
    irmin(1,6)= 9, irmax(1,6)=11, irmin(2,6)= 5, irmax(2,6)= 7, irmin(3,6)=72, irmax(3,6)=74,
    irmin(1,7)= 9, irmax(1,7)=11, irmin(2,7)=14, irmax(2,7)=16, irmin(3,7)=72, irmax(3,7)=74,
    irmin(1,8)= 1, irmax(1,8)=18, irmin(2,8)= 1, irmax(2,8)=27, irmin(3,8)=72, irmax(3,8)=74,
 /
EOF
$ECHO
$ECHO "  running local DOS calculation...\c"
$PROJWFC_COMMAND < AlAs110.box.projwfc.in > AlAs110.box.projwfc.out
check_failure $?
$ECHO " done"

# Projection of the DOS on atomic wavefunctions
cat > AlAs110.projwfc.in << EOF
 &projwfc
    prefix  = 'AlAs110'
    outdir='$TMP_DIR/',
    ngauss=0
    degauss=0.01
    DeltaE=0.02
    tdosinboxes=.false.
 /
EOF
$ECHO
$ECHO "  running projected DOS calculation...\c"
$PROJWFC_COMMAND < AlAs110.projwfc.in > AlAs110.projwfc.out
check_failure $?
$ECHO " done"

$ECHO
$ECHO "  summing the atomic PDOS...\c"
$SUMPDOS_COMMAND AlAs110.pdos_atm\#10\(Al\)_wfc* > "AlAs110.pdos_atm#10(Al)" 2> /dev/null
$SUMPDOS_COMMAND AlAs110.pdos_atm\#11\(As\)_wfc* > "AlAs110.pdos_atm#11(As)" 2> /dev/null
$ECHO " done"

#
#  if gnuplot was found, the results are plotted
#
if [ "$GP_COMMAND" = "" ]; then
    break
else
eFermi=`grep "Fermi" AlAs110re.nonscf.out | cut -d  \  -f 14`
cat > gnuplot.tmp <<EOF
set term postscript enhanced color solid lw 3 24
set output 'AlAs110.box.projwfc.ps'
ef=$eFermi
set xlabel "Energy - E_F (eV)"
set ylabel "Local DOS (states/eV)"
set style data lines
set key top left Left reverse
set border 31 lw 0.2
set title "Projected DOS"
plot \\
 './AlAs110.pdos_atm#10(Al)' u (\$1-ef):2 t "Surface Al" ,\\
 './AlAs110.pdos_atm#11(As)' u (\$1-ef):2 t "Surface As"
set title "Local DOS centered in the first vacuum layer"
plot \\
 './AlAs110.ldos_boxes' u (\$1-ef):4 t "Above Al" ,\\
 './AlAs110.ldos_boxes' u (\$1-ef):5 t "Above As" ,\\
 './AlAs110.ldos_boxes' u (\$1-ef):(\$7/54) t "Surface average"
set title "Local DOS centered in the second vacuum layer"
plot \\
 './AlAs110.ldos_boxes' u (\$1-ef):8 t "Above Al" ,\\
 './AlAs110.ldos_boxes' u (\$1-ef):9 t "Above As" ,\\
 './AlAs110.ldos_boxes' u (\$1-ef):(\$11/54) t "Surface average"
EOF
$ECHO
$ECHO "  plotting DOS results ...\c"
$GP_COMMAND < gnuplot.tmp
$ECHO " done"
rm gnuplot.tmp
fi

$ECHO
$ECHO "  To visualize a volume in which the DOS is integrated, execute:"
$ECHO "    xcrysden --xsf results/AlAs110.box#1.xsf"
$ECHO "  and plot the isosurface corresponding to isovalue 0.5"

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