#PBS -l nodes=8:ppn=12
#PBS -V
TOTAL_PROCS=96


        #************************************************************
        #                                                           *
	#  CALCULATION OF PHONON DISPERSION AND RAMAN SPECTROSCOPY  *
        #                                                           *
	#************************************************************

cd `echo $0 | sed 's/\(.*\)\/.*/\1/'` # extract pathname

EXAMPLE_DIR=/home2/soumya/qe-6.3/TEST/122/Ba22_ortho/raman

# check whether echo has the -e option
if test "`echo -e`" = "-e" ; then ECHO=echo ; else ECHO="echo -e" ; fi


# set the needed environment variables
. /home2/soumya/qe-6.3/environment_variables_$TOTAL_PROCS


# ----------- RESULT DIR ---------------------------
# result directory
for DIR in "$TMP_DIR" "$EXAMPLE_DIR/RAMAN" ; do
    if test ! -d $DIR ; then
        mkdir $DIR
    fi
done
cd $EXAMPLE_DIR/RAMAN
# ----------- RESULT DIR ---------------------------



# how to run executables
PW_COMMAND="$PARA_PREFIX $BIN_DIR/pw.x $PARA_POSTFIX"
$ECHO "  running pw.x as: $PW_COMMAND"
PH_COMMAND="$PARA_PREFIX $BIN_DIR/ph.x $PARA_POSTFIX"
$ECHO "  running ph.x as: $PH_COMMAND"
Q2R_COMMAND="$PARA_PREFIX $BIN_DIR/q2r.x $PARA_POSTFIX"
$ECHO "  running q2r.x as: $Q2R_COMMAND"
MATDYN_COMMAND="$PARA_PREFIX $BIN_DIR/matdyn.x $PARA_POSTFIX"
$ECHO "  running matdyn.x as: $MATDYN_COMMAND"
DM_COMMAND="$PARA_PREFIX $BIN_DIR/dynmat.x $PARA_POSTFIX"
$ECHO "  running dynmat.x as: $DM_COMMAND"
PLOTBAND_COMMAND="$BIN_DIR/plotband.x "
$ECHO "  running plotband.x as: $PLOTBAND_COMMAND"

$ECHO

# 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

## ===================================================
NAME='Ba22'
PREFIX='122'
Ecut=60
Erho=480
PSEUDO_DIR=/home2/soumya/qe-6.3/PSEUDO/QE_Special/Psilibrary-USPP-PBE
## ===================================================


# --------------------------------
# self-consistent calculation    |
# --------------------------------
cat > $NAME.scf.in << EOF
&CONTROL
  calculation = "scf",
  prefix = "$PREFIX",
  pseudo_dir = "$PSEUDO_DIR",
  outdir = "./",
  verbosity = "high"
  wf_collect = .true.
/
&system
    ibrav=0, 
    nat=5, 
    ntyp=3,
    occupations='smearing', 
    smearing='gaussian', 
    degauss=0.01,
    ecutwfc = $Ecut,
    ecutrho = $Erho,
                  ! nosym = .true.
 /
 &electrons
    electron_maxstep=1000
    mixing_mode='plain'
    mixing_beta = 0.5,
    startingwfc='random',
    conv_thr =  1.0d-8
    diagonalization = 'david'
               ! diago_full_acc = .true.
 /


ATOMIC_SPECIES
  Ba 137.327   Ba.pbe-spn-rrkjus_psl.1.0.0.UPF  
  Fe  55.845   Fe.pbe-spn-rrkjus_psl.1.0.0.UPF
  As  74.922   As.pbe-n-rrkjus_psl.1.0.0.UPF   


CELL_PARAMETERS (angstrom)
     0.0000000     2.7871000     6.4726500       
     2.8073000     0.0000000     6.4726500        
     2.8073000     2.7871000     0.0000000

ATOMIC_POSITIONS (crystal)
Fe	0.250000   0.250000   0.250000    
Fe     -0.250000  -0.250000  -0.250000    
As	0.353800   0.353800  -0.353800    
As     -0.353800  -0.353800   0.353800    
Ba	0.000000   0.000000   0.000000 

K_POINTS (automatic)
16 16 6 0 0 0

EOF

$ECHO "  running the SCF for $NAME...\c"
#$PW_COMMAND < $NAME.scf.in > $NAME.scf.out
check_failure $?
$ECHO " done"



		# ********************************
		#       RAMAN SPECTROSCOPY       #
		# ********************************


# -----------------------------------------------------------------
# self-consistent phonon calculation for RAMAN (Gamma only)       |
# -----------------------------------------------------------------
cat > $NAME.raman.in << EOF
Normal modes for $NAME
 &inputph
  tr2_ph=1.0d-14,     !Threshold for self-consistency
  prefix='$PREFIX',
  verbosity="high"
  outdir='./'
  amass(1)=137.327,   !Atomic mass [amu] of each atomic type
  amass(2)=55.845,
  amass(3)=74.921,
  !alpha_mix(1)=0.5,            !Mixing factor (for each iteration)
  !diagonalization='cg',
  fildvscf = 'dvscf',
  electron_phonon = 'simple'  !method of calculation
  el_ph_sigma = 0.02,               !spacing between double-delta smearing values
  el_ph_nsigma = 10                 !number of double-delta smearing values
  epsil=.false.,      !must be false for metallic system
  !lraman=.true.      !If .true. calculate non-resonant Raman coefficients
  trans=.true.,       !If .true. the phonons are computed.
  asr=.true.,          !Apply Acoustic Sum Rule to dynamical matrix
  fildyn='dmat-raman.$NAME' !File where the dynamical matrix is written
 /
 0.0 0.0 0.0          !list of q-points
EOF
$ECHO "  running phonon calculation for RAMAN of $NAME...\c"
#$PH_COMMAND < $NAME.raman.in > raman.out
check_failure $?
$ECHO " done"


# ----------------------------------------------
# extract phonon data with dynmat.x for RAMAN  |
# ----------------------------------------------
cat > $NAME.dm_raman.in << EOF
 &input 
   fildyn='dmat-raman.$NAME',   !File where the dynamical matrix is written
   asr='zero-dim'               !Indicates the type of Acoustic Sum Rule imposed
 /
EOF
$ECHO " Extracting phonon data with dynmat...\c"
$DM_COMMAND < $NAME.dm_raman.in > dm_raman.out
check_failure $?
$ECHO " done"



# --------------------------
# plotting using Gnuplot   |
# --------------------------
$ECHO "The data for spectrum is after '# mode   [cm-1]    [THz]      IR' on  $NAME.dm.out"
awk 'NR==1,/mode/{next}/DYNMAT/,NR==0{next}{print}' dm_raman.out > plotdata_$NAME.dat

$ECHO "Trying to plot it with gnuplot..."


cat > plot_command_$NAME.cmd << EOF
 set lmargin 8
 set rmargin 3

 set multiplot
 set key left top
 set origin 0.0,0.5
 set size 1,0.5
 set yrange [0:] 
 set format x ""
 set tmargin 1
 plot 'plotdata_$NAME.dat' u (\$2):(\$3) title ' $NAME-RAMAN' w i lw 2

 set key left bottom
 set origin 0.0,0.0
 set size 1,0.587
 set yrange [0:] reverse
 set format x
 set xlabel "Frequency [cm-1]"
 set bmargin 3
 set ylabel "Intensity" offset 0,5
 plot 'plotdata_$NAME.dat' u (\$2):(\$4) title '$NAME-IR' w i lw 2 lc 2 
 set nomultiplot
EOF


if [ "$GP_COMMAND" = "" ]; then
    $ECHO "No plot will be produced, because gnuplot was not found in the \$PATH "
    break
else
  $GP_COMMAND -persist plot_command_$NAME.cmd &
   
fi


$ECHO
$ECHO "$EXAMPLE_DIR: done"