#!/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 and postprocessing codes to make a"
$ECHO "contour plot in the [110] plane of the charge density for Si, and to"
$ECHO "plot the band structure of Si."

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

# required executables and pseudopotentials
BIN_LIST="pw.x pp.x plotrho.x bands.x plotband.x"
PSEUDO_LIST="Si.pz-vbc.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 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 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"
PP_COMMAND="$PARA_PREFIX $BIN_DIR/pp.x $PARA_POSTFIX"
PLOTRHO_COMMAND="$BIN_DIR/plotrho.x"
BANDS_COMMAND="$PARA_PREFIX $BIN_DIR/bands.x $PARA_POSTFIX"
PLOTBAND_COMMAND="$BIN_DIR/plotband.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 bands.x as:    $BANDS_COMMAND"
$ECHO "  running plotband.x as: $PLOTBAND_COMMAND"
$ECHO

# self-consistent calculation
cat > si.scf.in << EOF
 &control
    calculation='scf'
    restart_mode='from_scratch',
    prefix='si'
    pseudo_dir = '$PSEUDO_DIR/',
    outdir='$TMP_DIR/'
 /
 &system
    ibrav= 2, celldm(1)= 10.2, nat= 2, ntyp= 1,
    ecutwfc =18.0
 /
 &electrons
    conv_thr =  1.0d-8
    mixing_beta = 0.7
 /
ATOMIC_SPECIES
 Si  28.086  Si.pz-vbc.UPF
ATOMIC_POSITIONS alat
 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...\c"
$PW_COMMAND < si.scf.in > si.scf.out
check_failure $?
$ECHO " done"

# post-processing for charge density
cat > si.pp_rho.in << EOF
 &inputpp
    prefix  = 'si'
    outdir = '$TMP_DIR/'
    filplot = 'sicharge'
    plot_num= 0
 /
 &plot
    nfile = 1
    filepp(1) = 'sicharge'
    weight(1) = 1.0
    iflag = 2
    output_format = 2
    fileout = 'si.rho.dat'
    e1(1) =1.0, e1(2)=1.0, e1(3) = 0.0,
    e2(1) =0.0, e2(2)=0.0, e2(3) = 1.0,
    nx=56, ny=40
 /
EOF
$ECHO "  running pp.x to do a 2-d plot of the charge density...\c"
$PP_COMMAND < si.pp_rho.in > si.pp_rho.out
check_failure $?
$ECHO " done"

# plotrho
cat > si.plotrho.in << EOF
si.rho.dat
si.rho.ps
n
0 0.09 6
EOF
$ECHO "  running plotrho.x to generate rho.ps...\c"
$PLOTRHO_COMMAND < si.plotrho.in > si.plotrho.out
$ECHO " done"

# post-processing for charge density
cat > si.pp_rho_new.in << EOF
 &inputpp
    prefix  = 'si'
    outdir = '$TMP_DIR/'
    filplot = 'sicharge'
    plot_num= 0
 /
 &plot
    nfile = 1
    filepp(1) = 'sicharge'
    weight(1) = 1.0
    iflag = 2
    output_format = 7
    fileout = 'si.rho_new.dat'
    e1(1) =1.0, e1(2)=1.0, e1(3) = 0.0,
    e2(1) =0.0, e2(2)=0.0, e2(3) = 1.0,
    nx=141, ny=100
 /
EOF
$ECHO
$ECHO "  running pp.x to do another 2-d plot of the charge density...\c"
$PP_COMMAND < si.pp_rho_new.in > si.pp_rho_new.out
check_failure $?
$ECHO " done"


if [ "$GP_COMMAND" = "" ]; then
    break
else
cat > gnuplot.tmp <<EOF
#!$GP_COMMAND
#
set term png font ",18" enh size 1000,707 
set pm3d
set palette model HSV functions gray*0.75, 1, 0.9
set view 0,0
#
alat=10.2
set xra[0:1.4142136*alat]
set yra [0.:alat]
set xtics out nomirror
set ytics axis in offset -4.0,0 nomirror 
set label "r (a.u)" at 6.8,-2.2 center
set label "r (a.u)" at -1.7,5.0 rotate by 90 center
unset ztics
unset key
set colorbox 
#
set out 'si.charge.png'
set title "Si charge"
splot 'si.rho_new.dat' u 1:2:3 w pm3d
EOF
$ECHO "  generating si.charge.png...\c"
$GP_COMMAND < gnuplot.tmp
$ECHO " done"
#rm gnuplot.tmp

cat > gnuplot1.tmp <<EOF
#!$GP_COMMAND
set view map
set size square
unset surface
unset clabel
set contour
set dgrid3d 100,141
set cntrparam cubicspline
set table
#
#  Define here the countour values. Each set of countours will have the same
#  color and is written in a different file
#
set cntrparam levels discrete 0.01,0.02,0.03
set output "table.dat"
splot 'si.rho_new.dat' using 1:2:3 w l 

set cntrparam levels discrete 0.04,0.05,0.06
set output "table1.dat"
splot 'si.rho_new.dat' using 1:2:3 w l 

set cntrparam levels discrete 0.07,0.08
set output "table2.dat"
splot 'si.rho_new.dat' using 1:2:3 w l 
#
unset table
#
#  Now define a postcript terminal 
#
set encoding iso_8859_15
set terminal postscript enhanced solid color "Helvetica" 20
set output "si.contour.ps"
#
#  prepare the countour graph
#
set size ratio 1./1.4142
set key off
alat=10.2
set border lw 3
set label "Si" at 10.6,8.7 
set xrange [0:1.4142136*alat]
set yrange [0:alat]
set xlabel "r (a.u.)"
set ylabel "r (a.u.)"
#
#  Set contour labels
#
dato="0.08"
set obj 9 rect at 6.,1.3 size char strlen(dato)*0.6, char 0.6 
set obj 9 fillstyle solid noborder front
set label at 6.,1.3 dato front center font "Helvetica,12" tc rgb "blue"

dato="0.07"
set obj 10 rect at 7.9,1.3 size char strlen(dato)*0.6, char 0.6 
set obj 10 fillstyle solid noborder front
set label at 7.9,1.3 dato front center font "Helvetica,12" tc rgb "blue"

dato="0.06"
set obj 11 rect at 3.3,1.0 size char strlen(dato)*0.6, char 0.6 
set obj 11 fillstyle solid noborder front
set label at 3.3,1.0 dato front center font "Helvetica,12" tc rgb "green"

dato="0.01"
set obj 12 rect at 3.6,6.0 size char strlen(dato)*0.6, char 0.6 
set obj 12 fillstyle solid noborder front
set label at 3.6,6.0 dato front center font "Helvetica,12" tc rgb "red"

dato="0.02"
set obj 13 rect at 3.6,5.4 size char strlen(dato)*0.6, char 0.6 
set obj 13 fillstyle solid noborder front
set label at 3.6,5.4 dato front center font "Helvetica,12" tc rgb "red"

dato="0.03"
set obj 14 rect at 3.6,4.9 size char strlen(dato)*0.6, char 0.6 
set obj 14 fillstyle solid noborder front
set label at 3.6,4.9 dato front center font "Helvetica,12" tc rgb "red"

dato="0.04"
set obj 15 rect at 2.3,3.7 size char strlen(dato)*0.6, char 0.6 
set obj 15 fillstyle solid noborder front
set label at 2.3,3.7 dato front center font "Helvetica,12" tc rgb "green"

dato="0.05"
set obj 16 rect at 3.6,3.7 size char strlen(dato)*0.6, char 0.6 
set obj 16 fillstyle solid noborder front
set label at 3.6,3.7 dato front center font "Helvetica,12" tc rgb "green"

dato="0.05"
set obj 17 rect at 7.2,1.9 size char strlen(dato)*0.6, char 0.6 
set obj 17 fillstyle solid noborder front
set label at 7.2,1.9 dato front center font "Helvetica,12" tc rgb "green"

dato="0.01"
set obj 18 rect at 10.8,2.3 size char strlen(dato)*0.6, char 0.6 
set obj 18 fillstyle solid noborder front
set label at 10.8,2.3 dato front center font "Helvetica,12" tc rgb "red"

#
# Print the countour
#
plot "table.dat" u 1:2 w l lw 3 lc rgb "red", "table1.dat" u 1:2 w l lw 3 lc rgb "green", "table2.dat" u 1:2 w l lw 3 lc rgb "blue"
EOF
$ECHO "  generating contour plot of the charge si.contour.ps...\c"
$GP_COMMAND < gnuplot1.tmp
$ECHO " done"
rm gnuplot1.tmp table.dat table1.dat table2.dat
$ECHO 
fi


# band structure calculation along high-symmetry lines
cat > si.band.in << EOF
 &control
    calculation='bands'
    pseudo_dir = '$PSEUDO_DIR/',
    outdir='$TMP_DIR/',
    prefix='si'
 /
 &system
    ibrav=  2, celldm(1) =10.20, nat=  2, ntyp= 1,
    ecutwfc =18.0, nbnd = 8,
 /
 &electrons
 /
ATOMIC_SPECIES
 Si  28.086  Si.pz-vbc.UPF
ATOMIC_POSITIONS alat
 Si 0.00 0.00 0.00
 Si 0.25 0.25 0.25
K_POINTS tpiba_b
  5
   L 20
   gG 20
X 0
   1.0 1.0 0.0 30
   gG  1
EOF
$ECHO "  running the band-structure calculation for Si...\c"
$PW_COMMAND < si.band.in > si.band.out
check_failure $?
$ECHO " done"

# post-processing for band structure
cat > si.bands.in << EOF
 &bands
    prefix  = 'si'
    outdir = '$TMP_DIR/'
    filband = 'sibands.dat'
    lsym=.true.,
 /
EOF
$ECHO "  running the post-processing for band structure...\c"
$BANDS_COMMAND < si.bands.in > si.bands.out
check_failure $?
$ECHO " done"

# plotband.x
cat > si.plotband.in << EOF
sibands.dat
-6.0 10
sibands.xmgr
sibands.ps
6.255
1.0 6.255
EOF
$ECHO "  running plotband.x to generate sibands.ps...\c"
$PLOTBAND_COMMAND < si.plotband.in > si.plotband.out
check_failure $?
$ECHO " done"


if [ "$GP_COMMAND" = "" ]; then
    break
else
cat > gnuplot1.tmp <<EOF
set encoding iso_8859_15
set terminal postscript enhanced solid color "Helvetica" 20
set output "si.bands.ps"
#
set key off

dim1=-12.5
dim2=6
sqrt32=sqrt(3)*0.5
sqrt2=sqrt(2)
set xrange [0:sqrt2+sqrt32+1]
set yrange [dim1:dim2]
set arrow from sqrt32,dim1 to sqrt32,dim2 nohead lw 2 front
set arrow from sqrt32+1,dim1 to sqrt32+1,dim2 nohead lw 2 front
set arrow from 0,0 to sqrt2+sqrt32+1,0 nohead lw 1 front
set ylabel "Energy (eV)"
set label "Si" at 0.3,5 center
unset xtics
set border lw 2
lpos=dim1-0.45
set label "L" at 0.,lpos center
set label "{/Symbol G}" at sqrt32,lpos center
set label "X" at sqrt32+1,lpos center
set label "{/Symbol G}" at sqrt2+sqrt32+1,lpos center

set label "L_1" at 0.3,-10. center tc lt 1
set label "L_3" at 0.3,-2. center  tc lt 2

set label "{/Symbol G}_1" at sqrt32+0.02,-11.5 left
set label "{/Symbol G}_{25}'" at sqrt32+0.02,0.5 left
set label "{/Symbol G}_{15}" at sqrt32+0.02,2.0 left
set label "{/Symbol G}_{2}'" at sqrt32+0.02,4.6 left

set label "{/Symbol D}_1" at sqrt32+0.6,-9.5 center tc lt 1
set label "{/Symbol D}_2'" at sqrt32+0.6,-5.5 center tc lt 2
set label "{/Symbol D}_5" at sqrt32+0.6,-1.7 center tc lt 3

set label "{/Symbol S}_1" at sqrt32+1.7,-11 center tc lt 1
set label "{/Symbol S}_4" at sqrt32+1.7,-6.5 center tc lt 3
set label "{/Symbol S}_2" at sqrt32+1.7,-2.2 center tc lt 2
set label "{/Symbol S}_3" at sqrt32+1.7,4. center tc lt 5

hpos=0.6
set label "C_{3v}" at sqrt32*0.5,dim2+hpos center
set label "O_h" at sqrt32,dim2+hpos center
set label "C_{4v}" at sqrt32+0.5,dim2+hpos center
set label "C_{2v}" at sqrt32+1.0+sqrt2/2,dim2+hpos center
set label "O_h" at sqrt32+1.0+sqrt2,dim2+hpos center

vb_max=6.255

plot "sibands.xmgr.1.1"   u 1:(\$2-vb_max) w l lw 3  lt 1 ,\
"sibands.xmgr.1.3" u 1:(\$2-vb_max) w l lw 3  lt 2 ,\
"sibands.xmgr.2.1" u 1:(\$2-vb_max) w l lw 3  lt 1 ,\
"sibands.xmgr.2.4" u 1:(\$2-vb_max) w l lw 3  lt 2 ,\
"sibands.xmgr.2.5" u 1:(\$2-vb_max) w l lw 3  lt 3 ,\
"sibands.xmgr.3"   u 1:(\$2-vb_max) w l lw 3  lt 1 ,\
"sibands.xmgr.4.1" u 1:(\$2-vb_max) w l lw 3  lt 1 ,\
"sibands.xmgr.4.2" u 1:(\$2-vb_max) w l lw 3  lt 2 ,\
"sibands.xmgr.4.3" u 1:(\$2-vb_max) w l lw 3  lt 5 ,\
"sibands.xmgr.4.4" u 1:(\$2-vb_max) w l lw 3  lt 3 
EOF
$ECHO "  creating the postscript file si.bands.ps...\c"
$GP_COMMAND gnuplot1.tmp
check_failure $?
$ECHO " done"
fi
$ECHO

$ECHO "  cleaning $TMP_DIR...\c"
rm -rf $TMP_DIR/si.*

$ECHO
$ECHO "$EXAMPLE_DIR: done"