#!/bin/bash

name="MS-PBEl-rVV10"
name2="MS-PBEl"
a3D_orig="3.580"
a3D_exp="3.596"
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 + Cu(111) $name" font "Helvetica, 14"
set xlabel "Z [\305]" font "Helvetica, 14"
set ylabel "interaction energy [meV]" font "Helvetica, 14"
set yrange[-120: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"Cu(111) ${name} 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 [3.35:3.65]
plot "PLOTS/${name}_a3d" w lp lw 2 lt 1 lc 1 pt 7 pointsize 0.5 title"${name}", 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 + Cu(111) MS-PBEl-rVV10"
set title font "Helvetica, 10"

set lmargin at screen 0.3
set rmargin at screen 0.95


set title "(a)" offset -12,-2.5
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 [3.35:3.65]
plot "PLOTS/${name}_a3d" w lp lw 2 lt 1 lc 1 pt 7 pointsize 0.5 title"${name}", g(x) lw 2 lt 4 lc 2 title"${name2}", f(x) lw 2 lt 1 lc 0 title"exp."

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

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