#!/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 grafene."

# 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="C.tpss-mt.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_grafene" ; do
    if test ! -d $DIR ; then
        mkdir $DIR
    fi
done
cd $EXAMPLE_DIR/results_grafene

# 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
PARA_PREFIX="mpirun -np 8"
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 > grafene.scf.fit.in << EOF
 &control
    calculation='scf'
    restart_mode='from_scratch',
    prefix='grafene',
    tstress = .true.
    tprnfor = .true.
    etot_conv_thr = 1.0D-6
    forc_conv_thr = 1.0D-6
    pseudo_dir = '$PSEUDO_DIR/',
    outdir='$TMP_DIR/'
 /
 &system
    ibrav=  4, celldm(1) =4.65815, celldm(3)=4.05679d0, nat=  2, ntyp= 1,
    ecutwfc =80.0, ecutrho = 400.d0, occupations = 'smearing', smearing = 'cs', degauss= 0.001D0,
    input_dft='pbe',la2F=.true.,
 /
 &electrons
    conv_thr =  1.0d-8
    mixing_beta = 0.7
 /
ATOMIC_SPECIES
 C  12.00 C.tpss-mt.UPF
ATOMIC_POSITIONS
C        0.000000000   0.000000000   0.000318168
C        0.000000000   0.577350270   0.000681832
K_POINTS {automatic}
 16 16 4  0 0 0
EOF
$ECHO "  running the scf calculation with dense k-point grid...\c"
$PW_COMMAND  < grafene.scf.fit.in > grafene.scf.fit.out
check_failure $?
$ECHO "  done"
#
#  SCF at k-mesh good enough for phonons
#
cat > grafene.scf.in << EOF
 &control
    calculation='scf'
    restart_mode='from_scratch',
    prefix='grafene',
    tstress = .true.
    tprnfor = .true.
    etot_conv_thr = 1.0D-6
    forc_conv_thr = 1.0D-6
    pseudo_dir = '$PSEUDO_DIR/',
    outdir='$TMP_DIR/'
 /
 &system
    ibrav=  4, celldm(1) =4.65815, celldm(3)=4.05679d0, nat=  2, ntyp= 1,
    ecutwfc =80.0, ecutrho = 400.d0, occupations = 'smearing', smearing = 'cs', degauss= 0.001D0,
    input_dft='pbe',la2F = .true.,
 /
 &electrons
    conv_thr =  1.0d-8
    mixing_beta = 0.7
 /
ATOMIC_SPECIES
 C  12.00 C.tpss-mt.UPF
ATOMIC_POSITIONS
C        0.000000000   0.000000000   0.000318168
C        0.000000000   0.577350270   0.000681832
K_POINTS {automatic}
8 8 2 0 0 0 
EOF
$ECHO "  running the scf calculation...\c"
$PW_COMMAND < grafene.scf.in > grafene.scf.out
check_failure $?
$ECHO "  done"
#
cat > grafene.elph.in << EOF
Electron-phonon coefficients for graphene 
 &inputph
  tr2_ph=1.0d-10,
  prefix='grafene',
  fildvscf='grafenedv',
  amass(1)=12.00,
  outdir='$TMP_DIR/',
  fildyn='grafene.dyn',
  elph=.true.,
  trans=.true.,
  ldisp=.true.
  nq1=4, nq2=4, nq3=1
 /
EOF
$ECHO "  running the el-ph calculation...\c"
$PH_COMMAND < grafene.elph.in > grafene.elph.out
check_failure $?
$ECHO "  done"

# copy the dvscf info from TMP:
$ECHO " "
$ECHO " copy filedvscf into the current directory ."
cp $TMP_DIR/
#
#   q2r and matdyn
#
cat > q2r.in << EOF
 &input
  zasr='simple',  fildyn='grafene.dyn', flfrc='grafene881.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)=12.00,
    flfrc='grafene881.fc', flfrq='grafene881.freq', la2F=.true., dos=.false.
 /
91
 0.016667  0.009622  0.000000
 0.033333  0.019245  0.000000
 0.050000  0.028867  0.000000
 0.066667  0.038489  0.000000
 0.083333  0.048112  0.000000
 0.100000  0.057734  0.000000
 0.116667  0.067356  0.000000
 0.133333  0.076979  0.000000
 0.150000  0.086601  0.000000
 0.166667  0.096223  0.000000
 0.183333  0.105846  0.000000
 0.200000  0.115468  0.000000
 0.216667  0.125090  0.000000
 0.233333  0.134713  0.000000
 0.250000  0.144335  0.000000
 0.266667  0.153957  0.000000
 0.283333  0.163580  0.000000
 0.300000  0.173202  0.000000
 0.316667  0.182824  0.000000
 0.333333  0.192447  0.000000
 0.350000  0.202069  0.000000
 0.366667  0.211691  0.000000
 0.383333  0.221314  0.000000
 0.400000  0.230936  0.000000
 0.416667  0.240558  0.000000
 0.433333  0.250181  0.000000
 0.450000  0.259803  0.000000
 0.466667  0.269425  0.000000
 0.483333  0.279048  0.000000
 0.500000  0.288670  0.000000
 0.494444  0.298293  0.000000
 0.488889  0.307915  0.000000
 0.483333  0.317538  0.000000
 0.477777  0.327161  0.000000
 0.472222  0.336783  0.000000
 0.466666  0.346406  0.000000
 0.461110  0.356029  0.000000
 0.455555  0.365651  0.000000
 0.449999  0.375274  0.000000
 0.444443  0.384897  0.000000
 0.438888  0.394519  0.000000
 0.433332  0.404142  0.000000
 0.427776  0.413765  0.000000
 0.422221  0.423387  0.000000
 0.416665  0.433010  0.000000
 0.411109  0.442633  0.000000
 0.405554  0.452255  0.000000
 0.399998  0.461878  0.000000
 0.394442  0.471501  0.000000
 0.388887  0.481123  0.000000
 0.383331  0.490746  0.000000
 0.377775  0.500369  0.000000
 0.372220  0.509991  0.000000
 0.366664  0.519614  0.000000
 0.361108  0.529237  0.000000
 0.355553  0.538859  0.000000
 0.349997  0.548482  0.000000
 0.344441  0.558105  0.000000
 0.338886  0.567727  0.000000
 0.333330  0.577350  0.000000
 0.322219  0.558105  0.000000
 0.311108  0.538860  0.000000
 0.299997  0.519615  0.000000
 0.288886  0.500370  0.000000
 0.277775  0.481125  0.000000
 0.266664  0.461880  0.000000
 0.255553  0.442635  0.000000
 0.244442  0.423390  0.000000
 0.233331  0.404145  0.000000
 0.222220  0.384900  0.000000
 0.211109  0.365655  0.000000
 0.199998  0.346410  0.000000
 0.188887  0.327165  0.000000
 0.177776  0.307920  0.000000
 0.166665  0.288675  0.000000
 0.155554  0.269430  0.000000
 0.144443  0.250185  0.000000
 0.133332  0.230940  0.000000
 0.122221  0.211695  0.000000
 0.111110  0.192450  0.000000
 0.099999  0.173205  0.000000
 0.088888  0.153960  0.000000
 0.077777  0.134715  0.000000
 0.066666  0.115470  0.000000
 0.055555  0.096225  0.000000
 0.044444  0.076980  0.000000
 0.033333  0.057735  0.000000
 0.022222  0.038490  0.000000
 0.011111  0.019245  0.000000
 0.000000  0.000000  0.000000
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)=12.00,
    flfrc='grafene881.fc', flfrq='grafene881.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 
   10
0.0000000  0.0000000  0.0000000   1.00
0.1250000  0.0721688  0.0000000   6.00
0.2500000  0.1443376  0.0000000   6.00
0.3750000  0.2165064  0.0000000   6.00
0.5000000  0.2886751  0.0000000   3.00
0.1250000  0.2165064  0.0000000   6.00
0.2500000  0.2886751  0.0000000  12.00
0.3750000  0.3608439  0.0000000  12.00
0.2500000  0.4330127  0.0000000   6.00
0.3750000  0.5051815  0.0000000   6.00
elph. 0.000000. 0.000000. 0.000000
elph. 0.000000. 0.144338. 0.000000
elph. 0.000000. 0.288675. 0.000000
elph. 0.000000. 0.433013. 0.000000
elph. 0.000000.-0.577350. 0.000000
elph. 0.125000. 0.216506. 0.000000
elph. 0.125000. 0.360844. 0.000000
elph. 0.125000. 0.505181. 0.000000
elph. 0.250000. 0.433013. 0.000000
elph. 0.250000. 0.577350. 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"