#!/bin/bash
#
# Code by  venkentom@gmail.com
#
# Licenced by the GPL General Public License
# Have a look at:
# http://www.gnu.org/copyleft/gpl.html
# For more information about the GPL License
#
# I am not responsible for what this code does to your hw/sw/mind
# Use at your own risk
# 

#Put gromacs commands after this line

DIR=`pwd`

#Config
NPROCS=2

mkdir hbonds

#make index file
make_ndx -f system.pdb -o index_aa.ndx <<EOL
splitres 1
q
EOL

#loop over all residues
for i in {01..19}
do
k=`echo $i|sed 's/^0*//'`
j=$(($k + 17))

echo $j

#calculate number of contacts between POPC and protein residues
g_hbond -s system.tpr -f prod_time.xtc -dt 0.001 -num hbnum_r${i}-POPC.xvg -tu us -xvg none -n index_aa.ndx <<EOL 
$j
13
EOL

#calculate number of contacts between CHOL and protein residues
g_hbond -s system.tpr -f prod_time.xtc -dt 0.001 -num hbnum_r${i}-CHOL.xvg -tu us -xvg none -n index_aa.ndx <<EOL
$j
14
EOL

#calculate number of contacts between PW and protein residues
g_hbond -s system.tpr -f prod_time.xtc -dt 0.001 -num hbnum_r${i}-PW.xvg -tu us -xvg none -n index_aa.ndx <<EOL
$j
15
EOL

#rm mindist.xvg

done

#putting everything together
pr -m -t -s\  hbnum*-POPC.xvg | gawk '{print $2,$4,$6,$8,$10,$12,$14,$16,$18,$20,$22,$24,$26,$28,$30,$32,$34,$36,$38,$40}' > hbnum_POPC_all.xvg

pr -m -t -s\  hbnum*-CHOL.xvg | gawk '{print $2,$4,$6,$8,$10,$12,$14,$16,$18,$20,$22,$24,$26,$28,$30,$32,$34,$36,$38,$40}' > hbnum_CHOL_all.xvg

pr -m -t -s\  hbnum*-PW.xvg | gawk '{print $2,$4,$6,$8,$10,$12,$14,$16,$18,$20,$22,$24,$26,$28,$30,$32,$34,$36,$38,$40}' > hbnum_PW_all.xvg

#tranpose file for plotting
gawk '
{ 
    for (i=1; i<=NF; i++)  {
        a[NR,i] = $i
    }
}
NF>p { p = NF }
END {    
    for(j=1; j<=p; j++) {
        str=a[1,j]
        for(i=2; i<=NR; i++){
            str=str" "a[i,j];
        }
        print str
    }
}' hbnum_POPC_all.xvg > hbnum_POPC_all_transpose.xvg

gawk '
{
    for (i=1; i<=NF; i++)  {
        a[NR,i] = $i
    }
}
NF>p { p = NF }
END {
    for(j=1; j<=p; j++) {
        str=a[1,j]
        for(i=2; i<=NR; i++){
            str=str" "a[i,j];
        }
        print str
    }
}' hbnum_CHOL_all.xvg > hbnum_CHOL_all_transpose.xvg

gawk '
{ 
    for (i=1; i<=NF; i++)  {
        a[NR,i] = $i
    }
}
NF>p { p = NF }
END {    
    for(j=1; j<=p; j++) {
        str=a[1,j]
        for(i=2; i<=NR; i++){
            str=str" "a[i,j];
        }
        print str
    }
}' hbnum_PW_all.xvg > hbnum_PW_all_transpose.xvg

#plot
#gnuplot -e "datafile='mindistnum_POPC_all_transpose.xvg'; outputname='mindistnum_POPC.png'" analysis_contacts.plg
#gnuplot -e "datafile='mindistnum_CHOL_all_transpose.xvg'; outputname='mindistnum_CHOL.png'" analysis_contacts.plg
#gnuplot -e "datafile='mindistnum_PW_all_transpose.xvg'; outputname='mindistnum_PW.png'" analysis_contacts.plg


mv hbnum_* hbonds
rm \#*