#!/bin/bash

name="MS-PBEl-rVV10_Ag"
name2="MS-PBEl"
name3="MS-PBEl-rVV10"
metal="Ag"
a3D_orig="4.090"
a3D_exp="4.062"
vdwZexp="1.98"
vdwEexp="-32.4"
yrange_a3D="[3.8:4.1]"
samples=2500
bvalues="3.0 4.0 5.0 6.0 6.3 7.0 8.0 9.0 10.0 11.0 12.0 13.0 14.0 16.0 17.0 18.0 19.0 20.0"

if [ -d PLOTS ];
then
  rm -r PLOTS
fi
mkdir PLOTS

#a3D
sort -n ${name}_a3d > PLOTS/${name}_a3d

#now sort and find vdw Zmin and Emin as function of b
for bparam in ${bvalues}
do
  echo "now at b = ${bparam}"
  #find gas ref value, which is on first line
  gasref=`awk 'NR == 1 {print $3}' ${name}_b${bparam}_VDW`
  #sort and shift by gasref
  sort -n ${name}_b${bparam}_VDW | awk '{print $2, ($3 - '${gasref}')*1000}' > PLOTS/${name}_b${bparam}_VDW

#use gnuplot to make dense table and find minimum
gnuplot<<EOF
set table "table_dummy"
set samples ${samples}
plot "PLOTS/${name}_b${bparam}_VDW" smooth cspline
EOF
  awk 'NR > 4 && NR <= '${samples}' {print $1, $2}' table_dummy > PLOTS/${name}_b${bparam}_VDW_cspline
  rm table_dummy

  Eold=1000
  Zold=2
  for (( i=1; i<${samples}; i++))
  do
    #awk 'NR == '${i}' { print $1, $2}' PLOTS/${name}_b${bparam}_VDW_cspline
    Zval=`awk 'NR == '${i}' { print $1}' PLOTS/${name}_b${bparam}_VDW_cspline`
    Eval=`awk 'NR == '${i}' { print $2}' PLOTS/${name}_b${bparam}_VDW_cspline`
    if (( $( echo "${Eval} > ${Eold}" | bc -l ) )); 
    then
      echo "${bparam} ${Zold} ${Eold}" >> PLOTS/${name}_VDWmin
      break
    else
    Eold=${Eval}
    Zold=${Zval}
    fi

  done #i

  #generate plot command all vdW wells
  if [ "$bparam" == "20.0" ];
  then
    plotcom_all_vdw+=" \"PLOTS/${name}_b${bparam}_VDW_cspline\" u 1:2 w l title\"b = ${bparam}\" "
  else
    plotcom_all_vdw+=" \"PLOTS/${name}_b${bparam}_VDW_cspline\" u 1:2 w l title\"b = ${bparam}\", "
  fi
done #bparam
#echo $plotcom_all_vdw

#plot all vdW wells
gnuplot<<EOF
set output "PLOTS/${name}_allVDW.pdf"
set terminal pdf enhanced 
set key font "Helvetica, 10"
set title "H_2 + ${metla}(111) $name3" font "Helvetica, 14"
set xlabel "Z [\305]" font "Helvetica, 14"
set ylabel "interaction energy [meV]" font "Helvetica, 14"
set yrange[-150:20]
set xrange[2.5:8.0]
set key bottom right
plot $plotcom_all_vdw
EOF

#plot a3D
gnuplot<<EOF
set output "PLOTS/${name}_a3D.pdf"
set terminal pdf enhanced
set key font "Helvetica, 10"
set key bottom right
set xtics nomirror
set ytics nomirror
set title"${metal}(111) ${name3} lattice constant" font "Helvetica, 14"
set ylabel "a3D [\305]" font "Helvetica, 14" 
set xlabel "b parameter" font "Helvetica, 14"
f(x) = ${a3D_exp}
g(x) = ${a3D_orig}
set xrange [2.0:20]
set yrange ${yrange_a3D}
plot "PLOTS/${name}_a3d" w lp lw 2 lt 1 lc 1 pt 7 pointsize 0.5 title"${name3}", g(x) lw 2 lt 4 lc 2 title"${name2}", f(x) lw 2 lt 1 lc 0 title"exp." 
EOF

gnuplot<<EOF
set output "PLOTS/${name}_multiplot.pdf"
set terminal pdf enhanced size 3in,6in
set multiplot layout 3,1 title"H_2 + ${metal}(111) ${name3}"
set title font "Helvetica, 10"

set title "(a) a3D"
set xtics nomirror
set ytics nomirror
set tics font "Helvetica, 10"
set ylabel "a3D [\305]" font "Helvetica, 10"
set xlabel "b parameter" font "Helvetica, 10"
f(x) = ${a3D_exp}
g(x) = ${a3D_orig}
set key bottom right
set key font "Helvetica, 10"
set xrange [2.0:20]
set yrange ${yrange_a3D}
plot "PLOTS/${name}_a3d" w lp lw 2 lt 1 lc 1 pt 7 pointsize 0.5 title"${name3}", g(x) lw 2 lt 4 lc 2 title"${name2}", f(x) lw 2 lt 1 lc 0 title"exp."

set title "(b) Depht of vdW well"
set ylabel "interaction energy [meV]" font "Helvetica, 10"
set xlabel "b parameter" font "Helvetica, 10"
set yrange [-150:20]
f(x) = ${vdwEexp}
plot "PLOTS/${name}_VDWmin" u 1:3 w lp lw 2 lt 1 lc 1 pt 7 pointsize 0.5 title"${name3}", f(x) lw 2 lc 0 title"exp."

set title "(c) Position of vdW well minimum"
set ylabel "Z [\305]"  font "Helvetica, 10"
set xlabel "b parameter" font "Helvetica, 10"
set yrange[1.8:5.0]
f(x) = ${vdwZexp}
plot "PLOTS/${name}_VDWmin" u 1:2 w lp lw 2 lt 1 lc 1 pt 7 pointsize 0.5 title"${name3}", f(x) lw 2 lc 0 title"exp."
EOF