#!/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 mindistnum_hydro_withoutca

#make index file
make_ndx -f system.pdb -o index_aa.ndx <<EOL
"POP" & !a C21 C22 O* P N C11 C12 C13 H11* H12* H13* H14* H15* C1 C2 C3 HA HB HX HY HS
name 18 lipidtails
splitres 1
19 & !a N HN C O CA HA HT* NZ HZ*
20 & !a N HN C O CA HA SD
21 & !a N HN C O CA HA NZ HZ*
22 & !a N HN C O CA HA
23 & !a N HN C O CA HA
24 & !a N HN C O CA HA SD
25 & !a N HN C O CA HA
26 & !a N HN C O CA HA1 HA2
27 & !a N HN C O CA HA
28 & !a N HN C O CA HA SG HG1
29 & !a N HN C O CA HA
30 & !a N HN C O CA HA
31 & !a N HN C O CA HA
32 & !a N HN C O CA HA CG OD* HD2
33 & !a N HN C O CA HA
34 & !a N HN C O CA HA SD
35 & !a N HN C O CA HA
36 & !a N HN C O CA HA
37 & !a N HN C O CA HA CD OE* OT*
q
EOL

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

echo $j

#calculate number of contacts between POPC and protein residues
g_mindist -s system.pdb -f prod_time.xtc -n index_aa.ndx -on hydro_r${i}-POPC.xvg -tu us -xvg none -d 0.3 <<EOL
$j
18
EOL

rm mindist.xvg

done

#putting everything together
pr -m -t -s\  hydro*-POPC.xvg | gawk '{print $2,$4,$6,$8,$10,$12,$14,$16,$18,$20,$22,$24,$26,$28,$30,$32,$34,$36}' > hydro_POPC_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
    }
}' hydro_POPC_all.xvg > hydro_POPC_all_transpose.xvg

#plot
gnuplot -e "datafile='hydro_POPC_all_transpose.xvg'; outputname='hydro_POPC.png'" analysis_contacts.plg


mv hydro_* mindistnum_hydro_withoutca
rm \#*