#!/bin/bash
# computes  work from dhdl file
# syntax:
# $ works.bash project_name 
#

if [ $# -eq 0 ] ; then 
    echo " Extract work data from bound and ubound directories generated with gmx/multidir" 
    echo " syntax: $BIN_PATH/works.bash DIRE"
    echo "          where DIRE is a user-defined name for the drug-receptor pair"
    exit
fi 

if [ ! -d Results ]; then
    mkdir Results 
fi

# process bound state dhdl files 

histog=bin/histog
bound_dir=bound
# vdw work
rm blj.wrk >& /dev/null
for i in  $bound_dir/fsdam/b*/dhdlvdw.xvg; do
    # find the pass
    dl=`grep -v -e "[#@]" $i | wc -l | awk '{print dl=1/($1-1)}'`
    # do the integral
    grep -v -e "[#@]" $i | awk -v dl=$dl '{a+=$2*dl}END{print -a}' >> blj.wrk
done
echo "lj bound done"

# qq work 
rm bqq.wrk >& /dev/null
for i in  $bound_dir/fsdam/b*/dhdlQ.xvg; do
    # find the pass
    dl=`grep -v -e "[#@]" $i | wc -l | awk '{print dl=1/($1-1)}'`
    # do the integral
    grep -v -e "[#@]" $i | awk -v dl=$dl '{a+=$2*dl}END{print -a}' >> bqq.wrk
done
echo "qq bound done"

unbound_dir=unbound
# vdw work
rm ulj.wrk >& /dev/null
for i in $unbound_dir/fsdam/u*/dhdlvdw.xvg; do
    # find the pass
    dl=`grep -v -e "[#@]" $i | wc -l | awk '{print dl=1/($1-1)}'`
    # do the integral
    grep -v -e "[#@]" $i | awk -v dl=$dl '{a+=$2*dl}END{print  a}' >> ulj.wrk
done
echo "lj unbound done"

# qq work 
rm uqq.wrk >& /dev/null
for i in  $unbound_dir/fsdam/u*/dhdlQ.xvg; do
    # find the pass
    dl=`grep -v -e "[#@]" $i | wc -l | awk '{print dl=1/($1-1)}'`
    # do the integral
    grep -v -e "[#@]" $i | awk -v dl=$dl '{a+=$2*dl}END{print  a}' >> uqq.wrk
done
echo "qq unbound done"

paste ulj.wrk uqq.wrk   | awk '{print $1+$2}' >  temp-u.wrk
paste blj.wrk bqq.wrk   | awk '{print $1+$2}' >  temp-b.wrk

bin/zscore.awk temp-u.wrk | grep -v "#" > temp1-u.wrk
bin/zscore.awk temp1-u.wrk | grep -v "#" > temp2-u.wrk
bin/zscore.awk temp2-u.wrk | grep -v "#" > u.wrk

rm temp-u.wrk temp1-u.wrk temp2-u.wrk


bin/zscore.awk temp-b.wrk | grep -v "#" > temp1-b.wrk
bin/zscore.awk temp1-b.wrk | grep -v "#" > temp2-b.wrk
bin/zscore.awk temp2-b.wrk | grep -v "#" > b.wrk

rm temp-b.wrk temp1-b.wrk temp2-b.wrk


dire=$1

# now make the convolution of b and u wrk values as
# P(z) int P_b(z-x)P_u(x) dx 
# bootstrap (50 samples would do) loop to get the error
# we will do bootstrap on u and b work *indipendently* prior to compute
# the convoilution. This gives a better estimate of the error

rm -fr em.tmp jar.tmp
for i in {1..20} ; do  # bootstrap loop  
    fileu=u.wrk
    fileb=b.wrk
    bin/bs_res.awk -v seed=$i $fileu > u.tmp
    bin/bs_res.awk -v seed=$i $fileb > b.tmp
    nu=`wc -l $fileu | awk '{print $1}'`
    nb=`wc -l $fileb | awk '{print $1}'`
    cat u.tmp b.tmp  > filebu.tmp  # combine bootstrapped files
    awk -v nu=$nu -v nb=$nb '\
  {\
   if(NR<=nu)  {wu[NR]=$1} else {wb[NR-nu]=$1}\
       }\
  END{\
      for(i=1; i<=nu; i++) {\
          for(j=1; j<=nb; j++){\
             print  wu[i]+wb[j]\
          }\
      }\
}' filebu.tmp > $dire.BU.wrk
    awk '{print $1/4.184}' $dire.BU.wrk | $histog +i0.1 > $dire.BU.hst   
    #computes Jarzynski free energy on the bound-undbound convolution
    bin/free.awk $dire.BU.wrk | awk '{print $6}' >>  jar.tmp
    #computes EexpecMaxim estimate  free energy on the convolution
    bin/em/em 0.1 $i $dire.BU.wrk  >> em.tmp
done




grep "Components=  1" em.tmp |  awk '{a+=$4; a2+=$4^2}END{print a/NR,1.96*sqrt(a2/NR -(a/NR)^2)}' > DG.1 
grep "Components=  2" em.tmp |  awk '{a+=$4; a2+=$4^2}END{print a/NR,1.96*sqrt(a2/NR -(a/NR)^2)}' > DG.2
grep "Components=  3" em.tmp |  awk '{a+=$4; a2+=$4^2}END{print a/NR,1.96*sqrt(a2/NR -(a/NR)^2)}' > DG.3

# print free energy estimates 
DG1=`awk -v vol=$volcor -v qc=$corr '{print $1+vol+qc}' DG.1`
DG1err=`awk '{print $2}' DG.1`

printf " DG1= %4.1f" $DG1
printf " %4.1f" $DG1err


DG2=`awk -v vol=$volcor -v qc=$corr '{print $1+vol+qc}' DG.2`
DG2err=`awk '{print $2}' DG.2`

printf " DG2= %4.1f" $DG2
printf " %4.1f" $DG2err

DG3=`awk -v vol=$volcor -v qc=$corr '{print $1+vol+qc}' DG.3`
DG3err=`awk '{print $2}' DG.3`

printf " DG3= %4.1f" $DG3
printf " %4.1f" $DG3err

DGj=`awk '{a+=$1}END{print a/NR}' jar.tmp`
DGjerr=`awk '{a+=$1; a2+=$1^2}END{print 1.96*sqrt((a2/NR) - (a/NR)^2)}' jar.tmp`

printf " DGj= %4.1f" $DGj
printf " %4.1f\n" $DGjerr


nu=`wc -l $fileu | awk '{print $1}'`
nb=`wc -l $fileb | awk '{print $1}'`
cat u.wrk b.wrk  > filebu.tmp  # combine bootstrapped files
awk -v nu=$nu -v nb=$nb '\
  {\
   if(NR<=nu)  {wu[NR]=$1} else {wb[NR-nu]=$1}\
       }\
  END{\
      for(i=1; i<=nu; i++) {\
          for(j=1; j<=nb; j++){\
             print  wu[i]+wb[j]\
          }\
      }\
}' filebu.tmp > $dire.BU.wrk
awk '{print $1/4.184}' $dire.BU.wrk | $histog +i0.1 > $dire.BU.hst   

mv $dire.BU.hst Results
mv u.wrk Results/${dire}_u.wrk
mv b.wrk Results/${dire}_b.wrk
mv em.tmp Results/$dire.em
mv jar.tmp Results/$dire.jrz