#!/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 \
                http://www.quantum-espresso.org/pseudo/1.3/UPF/$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

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

# self-consistent calculation
cat > AlAs110re.scf.xml << EOF
<?xml version="1.0" encoding="UTF-8"?>


<input calculation="scf" prefix="AlAs110" title="AlAs 110 surface slab, relaxed (central plane fixed)">

	<cell type="qecell">
		<qecell ibrav="8" alat="7.424621202">
			<real rank="1" n1="5">
				1.414213576 6.00000 0.0 0.0 0.0
			</real>
		</qecell>
	</cell>

	<atomic_species ntyp="2">
		<specie name="Al">
			<property name="mass">
				<real>26.98</real>
			</property>
			<property name="pseudofile">
				<string>Al.pz-vbc.UPF</string>
			</property>
		</specie>
		<specie name="As">
			<property name="mass">
				<real>74.92</real>
			</property>
			<property name="pseudofile">
				<string>As.pz-bhs.UPF</string>
			</property>
		</specie>
	</atomic_species>

	<atomic_list units="alat" nat="14" >
		<atom name="As">
			<position>	
				<real rank="1" n1="3">
					0.000000000  -0.044777195  -0.058158722
				</real>
			</position>
		</atom>	
		<atom name="Al">
			<position>
				<real rank="1" n1="3">
					0.500000000   0.251460333   0.113525468
				</real>
			</position>
		</atom>
		<atom name="As">
			<position>
				<real rank="1" n1="3">
					0.500000000   0.712279009   0.504183885
				</real>
			</position>
		</atom>	
		<atom name="Al">
			<position>
				<real rank="1" n1="3">
					0.000000000   1.067633546   0.480460620
				</real>
			</position>
		</atom>	
		<atom name="As">
			<position>
				<real rank="1" n1="3">
					0.000000000  -0.003937059   0.995826731
				</real>
			</position>
		</atom>	
		<atom name="Al">
			<position>
				<real rank="1" n1="3">
					0.500000000   0.351400965   1.004220212
				</real>
			</position>
		</atom>	
		<atom name="As">
			<position>
				<real rank="1" n1="3">
					0.000000000  -0.003937059   2.004173269
				</real>
			</position>
		</atom>	
		<atom name="Al">
			<position>
				<real rank="1" n1="3">
					0.500000000   0.351400965   1.995779788
				</real>
			</position>
		</atom>	
		<atom name="As">
			<position>
				<real rank="1" n1="3">
					0.500000000   0.712279009   2.495816115
				</real>
			</position>
		</atom>	
		<atom name="Al">
			<position>
				<real rank="1" n1="3">
					0.000000000   1.067633546   2.519539380
				</real>
			</position>
		</atom>	
		<atom name="As">
			<position>
				<real rank="1" n1="3">
					0.000000000  -0.044777195   3.058158722
				</real>
			</position>
		</atom>	
		<atom name="Al">
			<position>
				<real rank="1" n1="3">
					0.500000000   0.251460333   2.886474532
				</real>
			</position>
		</atom>	
		<atom name="As">
			<position>
				<real rank="1" n1="3">
					0.500000000   0.707106800   1.500000000
				</real>
			</position>
		</atom>	
		<atom name="Al">
			<position>
				<real rank="1" n1="3">
					0.000000000   1.060660200   1.500000000
				</real>
			</position>
		</atom>				
	</atomic_list>		
	
	
	<field name="InputOutput">

		<parameter name="restart_mode">
			<string>
				from_scratch
			</string>
		</parameter>

		<parameter name="pseudo_dir">
			<string>
				$PSEUDO_DIR/
			</string>
		</parameter>
		
		<parameter name="outdir">
			<string>
				$TMP_DIR/
			</string>
		</parameter>
					
	</field>
	
	<field name="Numerics">

		<parameter name="ecutwfc">
			<real>
				14.0
			</real>
		</parameter>
		
		<parameter name="mixing_mode">
			<string>
				plain
			</string>
		</parameter>
		
		<parameter name="mixing_beta">
			<real>
				0.7
			</real>
		</parameter>
		
		<parameter name="conv_thr">
			<real>
				1.0d-6
			</real>
		</parameter>
	
	</field>
	
	<k_points type="automatic">
		<mesh>
			<integer rank="1" n1="6">
				6 2 1 0 0 0
			</integer>
		</mesh> 
	</k_points>
</input>
EOF
$ECHO "  running the scf calculation...\c"
$PW_COMMAND < AlAs110re.scf.xml > AlAs110re.scf.out
check_failure $?
$ECHO " done"

cat > AlAs110re.nonscf.xml << EOF
<?xml version="1.0" encoding="UTF-8"?>


<input calculation="nscf" prefix="AlAs110" title="AlAs 110 surface slab, relaxed (central plane fixed)">

	<cell type="qecell">
		<qecell ibrav="8" alat="7.424621202">
			<real rank="1" n1="5">
				1.414213576 6.00000 0.0 0.0 0.0
			</real>
		</qecell>
	</cell>

	<atomic_species ntyp="2">
		<specie name="Al">
			<property name="mass">
				<real>26.98</real>
			</property>
			<property name="pseudofile">
				<string>Al.pz-vbc.UPF</string>
			</property>
		</specie>
		<specie name="As">
			<property name="mass">
				<real>74.92</real>
			</property>
			<property name="pseudofile">
				<string>As.pz-bhs.UPF</string>
			</property>
		</specie>
	</atomic_species>

	<atomic_list units="alat" nat="14" >
		<atom name="As">
			<position>
				<real rank="1" n1="3">
					0.000000000  -0.044777195  -0.058158722
				</real>
			</position>
		</atom>	
		<atom name="Al">
			<position>
				<real rank="1" n1="3">
					0.500000000   0.251460333   0.113525468
				</real>
			</position>
		</atom>
		<atom name="As">
			<position>
				<real rank="1" n1="3">
					0.500000000   0.712279009   0.504183885
				</real>
			</position>
		</atom>	
		<atom name="Al">
			<position>
				<real rank="1" n1="3">
					0.000000000   1.067633546   0.480460620
				</real>
			</position>
		</atom>	
		<atom name="As">
			<position>
				<real rank="1" n1="3">
					0.000000000  -0.003937059   0.995826731
				</real>
			</position>
		</atom>	
		<atom name="Al">
			<position>
				<real rank="1" n1="3">
					0.500000000   0.351400965   1.004220212
				</real>
			</position>
		</atom>	
		<atom name="As">
			<position>
				<real rank="1" n1="3">
					0.000000000  -0.003937059   2.004173269
				</real>
			</position>
		</atom>	
		<atom name="Al">
			<position>
				<real rank="1" n1="3">
					0.500000000   0.351400965   1.995779788
				</real>
			</position>
		</atom>	
		<atom name="As">
			<position>
				<real rank="1" n1="3">
					0.500000000   0.712279009   2.495816115
				</real>
			</position>
		</atom>	
		<atom name="Al">
			<position>
				<real rank="1" n1="3">
					0.000000000   1.067633546   2.519539380
				</real>
			</position>
		</atom>	
		<atom name="As">
			<position>
				<real rank="1" n1="3">
					0.000000000  -0.044777195   3.058158722
				</real>
			</position>
		</atom>	
		<atom name="Al">
			<position>
				<real rank="1" n1="3">
					0.500000000   0.251460333   2.886474532
				</real>
			</position>
		</atom>	
		<atom name="As">
			<position>
				<real rank="1" n1="3">
					0.500000000   0.707106800   1.500000000
				</real>
			</position>
		</atom>	
		<atom name="Al">
			<position>
				<real rank="1" n1="3">
					0.000000000   1.060660200   1.500000000
				</real>
			</position>
		</atom>				
	</atomic_list>	
	
	
	<field name="InputOutput">

		<parameter name="restart_mode">
			<string>
				from_scratch
			</string>
		</parameter>

		<parameter name="pseudo_dir">
			<string>
				$PSEUDO_DIR/
			</string>
		</parameter>
		
		<parameter name="outdir">
			<string>
				$TMP_DIR/
			</string>
		</parameter>
					
	</field>
	
	<field name="Numerics">

		<parameter name="ecutwfc">
			<real>
				14.0
			</real>
		</parameter>
		
		<parameter name="mixing_mode">
			<string>
				plain
			</string>
		</parameter>
		
		<parameter name="mixing_beta">
			<real>
				0.7
			</real>
		</parameter>
		
		<parameter name="conv_thr">
			<real>
				1.0d-6
			</real>
		</parameter>
	
	</field>
	
	<field name="Options">

		<parameter name="occupations">
			<string>
				smearing
			</string>
		</parameter>
		
		<parameter name="smearing">
			<string>
				gaussian
			</string>
		</parameter>
		
		<parameter name="degauss">
			<real>
				0.01
			</real>
		</parameter>
		
		<parameter name="nbnd">
			<integer>
				34
			</integer>
		</parameter>
	
	</field>	
	
	<k_points type="automatic">
		<mesh>
			<integer rank="1" n1="6">
				12 4 1 0 0 0
			</integer>
		</mesh> 
	</k_points>
</input>
EOF
$ECHO "  running the non-scf calculation...\c"
$PW_COMMAND < AlAs110re.nonscf.xml > 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 (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"

# Projection of the DOS on volumes (boxes)
cat > AlAs110.box.projwfc.in << EOF
 &inputpp
    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
 &inputpp
    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
$ECHO " $EXAMPLE_DIR: done"