#!/bin/sh

###############################################################################
##
##  HIGH VERBOSITY EXAMPLE
##
###############################################################################

# 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

# function to test the exit status of a job
. ../check_failure.sh

$ECHO
$ECHO "$EXAMPLE_DIR : starting"
$ECHO
$ECHO "This example shows how to calculate electron-phonon interaction"
$ECHO "coefficients for fcc Al."

# set the needed environment variables
. ../environment_variables
# required executables and pseudopotentials
BIN_LIST="pw.x ph.x q2r.x matdyn.x lambda.x"
PSEUDO_LIST="Al.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 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
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 "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"
PH_COMMAND="$PARA_PREFIX $BIN_DIR/ph.x $PARA_POSTFIX"
Q2R_COMMAND="$PARA_PREFIX $BIN_DIR/q2r.x $PARA_POSTFIX"
MATDYN_COMMAND="$PARA_PREFIX $BIN_DIR/matdyn.x $PARA_POSTFIX"
LAMBDA_COMMAND="$BIN_DIR/lambda.x "
$ECHO
$ECHO "  running pw.x as: $PW_COMMAND"
$ECHO "  running ph.x as: $PH_COMMAND"
$ECHO "  running q2r.x as: $Q2R_COMMAND"
$ECHO "  running matdyn.x as: $MATDYN_COMMAND"
$ECHO "  running lambda.x as: $LAMBDA_COMMAND"
$ECHO

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

#
# SCF at dense k-mesh, good enough for electronic DOS
#
cat > al.scf.fit.in << EOF
 &control
    calculation='scf'
    restart_mode='from_scratch',
    prefix='al',
    pseudo_dir = '$PSEUDO_DIR/',
    outdir='$TMP_DIR/'
 /
 &system
    ibrav=  2, celldm(1) =7.5, nat= 1, ntyp= 1,
    ecutwfc =15.0,
    occupations='smearing', smearing='methfessel-paxton', degauss=0.05,
    la2F = .true.,
 /
 &electrons
    conv_thr =  1.0d-8
    mixing_beta = 0.7
 /
ATOMIC_SPECIES
 Al  26.98 Al.pz-vbc.UPF
ATOMIC_POSITIONS
 Al 0.00 0.00 0.00
K_POINTS {automatic}
 16 16 16  0 0 0
EOF
$ECHO "  running the scf calculation with dense k-point grid...\c"
$PW_COMMAND  < al.scf.fit.in > al.scf.fit.out
check_failure $?
$ECHO "  done"
#
#  SCF at k-mesh good enough for phonons
#
cat > al.scf.in << EOF
 &control
    calculation='scf'
    restart_mode='from_scratch',
    prefix='al',
    pseudo_dir = '$PSEUDO_DIR/',
    outdir='$TMP_DIR/'
 /
 &system
    ibrav=  2, celldm(1) =7.5, nat= 1, ntyp= 1,
    ecutwfc =15.0,
    occupations='smearing', smearing='methfessel-paxton', degauss=0.05
 /
 &electrons
    conv_thr =  1.0d-8
    mixing_beta = 0.7
 /
ATOMIC_SPECIES
 Al  26.98 Al.pz-vbc.UPF
ATOMIC_POSITIONS
 Al 0.00 0.00 0.00
K_POINTS {automatic}
 8 8 8  0 0 0
EOF
$ECHO "  running the scf calculation...\c"
$PW_COMMAND < al.scf.in > al.scf.out
check_failure $?
$ECHO "  done"
#
cat > al.ph.in << EOF
Electron-phonon coefficients for Al
 &inputph
  tr2_ph=1.0d-10,
  prefix='al',
  fildvscf='aldv',
  amass(1)=26.98,
  outdir='$TMP_DIR/',
  fildyn='al.dyn',
  elph=.false.,
  trans=.true.,
  ldisp=.true.
  nq1=4, nq2=4, nq3=4
 /
EOF
$ECHO "  running the ph calculation...\c"
$PH_COMMAND < al.ph.in > al.elph.out
check_failure $?
$ECHO "  done"
#
cat > al.elph.in << EOF
Electron-phonon coefficients for Al
 &inputph
  tr2_ph=1.0d-10,
  prefix='al',
  fildvscf='aldv',
  amass(1)=26.98,
  outdir='$TMP_DIR/',
  fildyn='al.dyn',
  elph=.true.,
  trans=.false.,
  ldisp=.true.
  nq1=4, nq2=4, nq3=4
 /
EOF
$ECHO "  running the el-ph calculation...\c"
$PH_COMMAND < al.elph.in > al.elph.out
check_failure $?
$ECHO "  done"
#
#   q2r and matdyn
#
cat > q2r.in << EOF
 &input
  zasr='simple',  fildyn='al.dyn', flfrc='Al444.fc', la2F=.true.
 /
EOF
$ECHO "  running q2r...\c"
$Q2R_COMMAND < q2r.in > q2r.out
check_failure $?
$ECHO "  done"
#
#
#
cat > matdyn.in.freq << EOF
 &input
    asr='simple',  amass(1)=26.98,
    flfrc='Al444.fc', flfrq='Al444.freq', la2F=.true., dos=.false.
 /
  19
  0.000 0.0 0.0     0.0
  0.125 0.0 0.0     0.0
  0.250 0.0 0.0     0.0
  0.375 0.0 0.0     0.0
  0.500 0.0 0.0     0.0
  0.750 0.0 0.0     0.0
  1.000 0.0 0.0     0.0
  0.825 0.125 0.125 0.0
  0.750 0.250 0.250 0.0
  0.625 0.375 0.375 0.0
  0.500 0.500 0.500 0.0
  0.325 0.325 0.325 0.0
  0.250 0.250 0.250 0.0
  0.125 0.125 0.125 0.0
  0.000 0.000 0.000 0.0
  0.125 0.125 0.000 0.0
  0.250 0.250 0.000 0.0
  0.325 0.325 0.000 0.0
  0.500 0.500 0.000 0.0
EOF
$ECHO "  running matdyn for frequency calculation...\c"
$MATDYN_COMMAND < matdyn.in.freq > matdyn.out.freq
check_failure $?
$ECHO "  done"
#
#
#
cat > matdyn.in.dos << EOF
 &input
    asr='simple',  amass(1)=26.98,
    flfrc='Al444.fc', flfrq='Al444.freq', la2F=.true., dos=.true.
    fldos='phonon.dos', nk1=10, nk2=10, nk3=10, ndos=50
 /
EOF
$ECHO "  running matdyn for a2F(omega) calculation...\c"
$MATDYN_COMMAND < matdyn.in.dos > matdyn.out.dos
check_failure $?
$ECHO "  done"
#
# 
#
cat > lambda.in << EOF
10  0.12  1    ! emax (something more than highest phonon mode in THz), degauss, smearing method 
    8          ! Number of q-points for which EPC is calculated, 
    0.0000000  0.0000000  0.0000000   1.00  ! the first q-point, use kpoints.x program to calculate 
    -0.2500000 -0.2500000  0.2500000   8.00  ! q-points and their weight
    -0.5000000 -0.5000000  0.5000000   4.00  !
     0.0000000  0.0000000  0.5000000   6.00  ! 4th q-point, qx,qy,qz
    -0.2500000 -0.2500000  0.7500000  24.00  !
    -0.5000000 -0.5000000  1.0000000  12.00  !
     0.0000000  0.0000000  1.0000000   3.00  !
    -0.5000000  0.0000000  1.0000000   6.00  ! the last q-point 
elph. 0.000000. 0.000000. 0.000000 ! elph output file names,
elph.-0.250000. 0.250000.-0.250000 ! in the same order as the q-points before  
elph. 0.500000.-0.500000. 0.500000
elph. 0.000000. 0.500000. 0.000000
elph. 0.750000.-0.250000. 0.750000
elph. 0.500000. 0.000000. 0.500000
elph. 0.000000.-1.000000. 0.000000
elph.-0.500000.-1.000000. 0.000000
0.10                     ! \mu the Coloumb coefficient in the modified
                         ! Allen-Dynes formula for T_c (via \omega_log)
EOF
$ECHO "  running lambda.x for lambda calculation...\c"
$LAMBDA_COMMAND < lambda.in > lambda.out
check_failure $?
$ECHO "  done"