#!/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 will calculate the water dipole and calculate the work"
$ECHO "function on a Ni slab with a CO molecule adsorbed using the dipole"
$ECHO "correction"

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

# required executables and pseudopotentials
BIN_LIST="pw.x pp.x average.x"
PSEUDO_LIST="C.pbe-rrkjus.UPF O.pbe-rrkjus.UPF Ni.pbe-nd-rrkjus.UPF HUSPBE.RRKJ3"

$ECHO
$ECHO "  executables directory: $BIN_DIR"
$ECHO "  pseudo directory:      $PSEUDO_DIR"
$ECHO "  temporary directory:   $TMP_DIR"
$ECHO
$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"
AVERAGE_COMMAND="$BIN_DIR/average.x"
$ECHO
$ECHO "  running pw.x as: $PW_COMMAND"
$ECHO "  running pp.x as: $PP_COMMAND"
$ECHO "  running average.x as: $AVERAGE_COMMAND"
$ECHO "  running gnuplot as:    $GP_COMMAND"
$ECHO

#
# Run the calculation for the Ni+CO slab with the dipole correction
# emaxpos, the starting of the potential inversion, has to be carefully
# placed in a position where there's almost no charge.
#
cat > ni+co.scf.in << EOF
 &control
    calculation='scf',
    restart_mode='from_scratch',
    prefix='ni+co',
    pseudo_dir = '$PSEUDO_DIR',
    outdir='$TMP_DIR'
    
    tefield=.true.,
    dipfield=.true.,
 /
&system

	nat=5, ntyp=3,
	ibrav=0, celldm(1)=4.70366666,

    ecutwfc = 30.0
	occupations='smearing', smearing='m-v', degauss=0.03

    edir=3,
    emaxpos=0.55,
    eopreg=0.06,
    eamp=0,
	
/
&electrons
	mixing_beta = 0.3
	conv_thr =  1.0d-6
/

ATOMIC_SPECIES
C	1.0	C.pbe-rrkjus.UPF
O	1.0	O.pbe-rrkjus.UPF
Ni	1.0	Ni.pbe-nd-rrkjus.UPF  

CELL_PARAMETERS alat
1.00000000 0.00000000 0.00000000
0.00000000 1.41421356 0.00000000
0.00000000 0.00000000 9.10000001

ATOMIC_POSITIONS (alat)

C       -0.00364039     0.02119538      1.54673745
O       -0.00634860     0.04192428      2.02021975
Ni      0.48527378      0.00197332      0.97713547
Ni      -0.00049546     0.70236680      0.45417840
Ni      0.50000000      0.00000000      0.00000000

K_POINTS {gamma}
EOF
$ECHO
$ECHO "  running pw.x for Ni+CO slab...\c"
$PW_COMMAND < ni+co.scf.in > ni+co.scf.out
check_failure $?
$ECHO " done"

#
# Extract the potential from the SCF calculation 
#
cat > ni+co.pp.in << EOF
&inputpp
    prefix='ni+co',
    outdir='$TMP_DIR',
    filplot = 'ni+co.vpot'
    plot_num= 11
/
EOF
$ECHO
$ECHO "  running pp.x to extract ni+co potential...\c"
$PP_COMMAND < ni+co.pp.in > ni+co.pp.out
check_failure $?
$ECHO " done"

#
# Average the potential to obtain the planar average along edir
# 
cat > ni+co.avg.in << EOF
1
ni+co.vpot
1.D0
150
3
3.000000
EOF
$ECHO
$ECHO "  running average.x to obtain the potential along Z...\c"
$AVERAGE_COMMAND < ni+co.avg.in > ni+co.vpot-z
check_failure $?
$ECHO " done"
mv avg.dat ni+co.vpot-z 

#
# If gnuplot is present plot the potential. In the graph it's clear where
# the potential reaches a constant value thanks to the dipole correction.
#

if [ "$GP_COMMAND" = "" ]; then
    break
else    
cat > ni+co.gnuplot <<EOF
set terminal postscript enhanced color
set output 'ni+co.eps'

set multiplot

set size 1,1   
set origin 0,0

set style line 1 lt 1
set style line 2 lt 2
set style line 3 lt 3 

set arrow 1 from 8,-2.2976 to 8,4.75 ls 3
set arrow 2 from 8,4.75 to 8,-2.2976 ls 3
#set arrow 3 from 18,4.488 to 16.5,-4  nohead ls 3

#set object 1 rect from 16.5,4.488 to 19.05048,4.896 fc ls 3 fs empty

set xlabel "Length (Angstroms)"
set ylabel "Energy (eV)"

set label 'W_{funct} = 7.04 eV' at 8.2,0

plot '< tail -150 ni+co.vpot-z' u (\$1*0.52918):(\$2*13.6) w l ls 1 title 'Pot.' , '< tail -150 ni+co.vpot-z' u (\$1*0.52918):(-2.2976) w l ls 2 title 'E_{f}'



set size 0.4,0.4 
set origin 0.45,0.1
set xrange [8:18]
set yrange[4.74:4.76]

set xtics 8.2,2.8,18
set ytics 4.742, 0.005, 4.758

set style line 1 lt 1
set style line 2 lt 2 
set style line 3 lt 3 

set arrow 1 from 10,4.7566 to 10,4.7529 ls 3
set arrow 2 from 10,4.7529 to 10,4.7566 ls 3
#unset arrow 3

set label "Dipole corr." at 10.8,4.748

unset xlabel 
unset ylabel 

plot [8:18][4.74:4.76] "< tail -150 ni+co.vpot-z" u (\$1*0.52918):(\$3*13.6) w l ls 1 notitle, "< tail -150 ni+co.vpot-z" u (\$1*0.52918):(4.7566) w l ls 2 notitle, "< tail -150 ni+co.vpot-z" u (\$1*0.52918):(4.7529) w l ls 2 notitle

unset multiplot
set output
reset
EOF
$ECHO "  Plotting averages and work function results ...\c"
$GP_COMMAND < ni+co.gnuplot
$ECHO " done"
rm ni+co.gnuplot
fi


###############################
#        WATER
###############################


#
# Run the calculation for the water molecule with the dipole correction
#
cat > water.scf.in << EOF
 &control
    calculation='scf',
    restart_mode='from_scratch',
    prefix='water',
    pseudo_dir = '$PSEUDO_DIR',
    outdir='$TMP_DIR'
    
   tefield=.true., 
    dipfield=.true.
/
&SYSTEM
    ibrav=1, celldm(1) = 15 
    nat=3, ntyp=2

    ecutwfc=30.0
    occupations='smearing', degauss=0.01
    
    edir=3
    eamp=0.D0
    eopreg=0.1
    emaxpos=0.6
/
&ELECTRONS
    mixing_beta = 0.7
    conv_thr =  1.0d-8
/
ATOMIC_SPECIES
    O    15.9994    O.pbe-rrkjus.UPF
    H    1.00794    HUSPBE.RRKJ3
ATOMIC_POSITIONS { Angstrom }
O    0.0  0.0  0.0
H    0.77 0.0  0.62
H   -0.77 0.0  0.62
K_POINTS { gamma }
EOF
$ECHO
$ECHO "  running pw.x for water molecule ...\c"
$PW_COMMAND < water.scf.in > water.scf.out
check_failure $?
$ECHO " done"

#
# Extract the potential from the SCF calculation  
#
cat > water.pp.in << EOF
&inputpp
    prefix='water',
    outdir='$TMP_DIR',
    filplot = 'water.vpot'
    plot_num= 11
/
EOF
$ECHO
$ECHO "  running pp.x to extract water potential...\c"
$PP_COMMAND < water.pp.in > water.pp.out
check_failure $?
$ECHO " done"

#
# If gnuplot is present plot the potential. In the graph it's clear how the
# dipole influences the potential slope
#
cat > water.avg.in << EOF
1
water.vpot
1.D0
90
3
3.000000
EOF
$ECHO
$ECHO "  running average.x to obtain the potential along Z...\c"
$AVERAGE_COMMAND < water.avg.in > water.vpot-z
check_failure $?
$ECHO " done"
mv avg.dat water.vpot-z


if [ "$GP_COMMAND" = "" ]; then
    break
else    
cat > water.gnuplot <<EOF
set terminal postscript enhanced color
set output 'water.eps'

set style line 1 lt 1
set style line 2 lt 2
set style line 3 lt 3

set arrow 1 from 2.5,-0.15 to 2.5,0.93 ls 3
set arrow 2 from 2.5,0.93 to 2.5,-0.15 ls 3

set label "Dipole correction" at 3,0.5

set yrange [-1.2:1.2]

set xlabel "Length (Angstroms)" 
set ylabel "Energy (eV)" 

plot "< tail -90 water.vpot-z" u (\$1*0.52918):(\$2*13.6) w l ls 1 title "Pot.", "< tail -90 water.vpot-z" u (\$1*0.52918):(-0.15) w l ls 2 notitle, "< tail -90 water.vpot-z" u (\$1*0.52918):(0.93) w l ls 2 notitle 

set output
reset

EOF
$ECHO "  Plotting averages and work function results ...\c"
$GP_COMMAND < water.gnuplot
$ECHO " done"
rm water.gnuplot
fi

$ECHO "  cleaning $TMP_DIR...\c"
rm -rf $TMP_DIR/ni+co.* $TMP_DIR/water.*

$ECHO "Calculations done."
$ECHO ""
$ECHO "----------------------------------------------------"
$ECHO " results/ni+co.eps -> Graph for the slab calculation"
$ECHO " results/water.eps -> Graph for the water molecule"
$ECHO "----------------------------------------------------"