#! /bin/bash

# This script writes a Gromacs topology template from the final pdb file generated by MEMGEN
#
# http://memgen.uni-goettingen.de/
#
# Please cite:
#
# Knight and Hub, MemGen: A general web server for the setup of lipid membrane simulation systems, submitted
#

if [ $# -ne 1 ]; then
    echo -e "\nUsage: finalpdb2GromacsTop.sh final.pdb\n"
    exit 1
fi

pdb=$1

if [ ! -e "$pdb" ]; then
    echo "Error, file $pdb not found"
    exit 1
fi

nSOL=$(grep REMARK "$pdb" | grep -w Water | tr -s ' ' | cut -d" " -f3)
nNA=$(grep REMARK "$pdb" | grep -w "Na+" | tr -s ' ' | cut -d" " -f3)
nCL=$(grep REMARK "$pdb" | grep -w "Cl-" | tr -s ' ' | cut -d" " -f3)

grep REMARK "$pdb" | sed -n '/THIS PDB FILE CONTAINS/,$ p' | grep -v "REMARK$" | \
    egrep -v 'Na\+|Cl-|Water' | sed 's/THIS PDB FILE CONTAINS://' > lipids.tmp

lipids=()
N=()
while read x lipid n; do
    lipids+=$lipid
    N+=$n
done < lipids.tmp
nTypes=${#lipids[@]}

echo Found $nTypes different lipids in $pdb. 

{
    cat <<EOF
;
; Topology template for GROMACS
;
; Generated by MEMGEN, http://memgen.uni-goettingen.de
;

; Please fill in your force field
#include "YOUR-FORCE-FIELD.ff/forcefield.itp"

; Include lipid itp files. Please replace by your itp file names:
EOF

    for ((i=0; i<nTypes; i++)); do
	printf '#include "%s.itp"\n' ${lipids[i]}
    done
    echo
    
    cat <<EOF
; Include water topology. 
; Replace force field and water model (such as tip3p, spc, ...):
#include "YOUR-FORCE-FIELD.ff/WATERMODEL.itp"
#include "YOUR-FORCE-FIELD.ff/ions.itp"

[ system ]
; Name
EOF
    printf "Membrane of "
    for ((i=0; i<nTypes; i++)); do
	printf '%d %s' ${N[i]} ${lipids[i]}
	[ $i -lt $[nTypes-1] ] && printf ', '
    done
    echo 
    echo
    echo '[ molecules ]'
    echo '; Compound        #mols'
    echo '; Replace compound name by moleculetype(s) in your itp file(s)'
    for ((i=0; i<nTypes; i++)); do
        printf '%-18s %d\n' ${lipids[i]} ${N[i]}
    done
    printf '%-18s %d\n' SOL $nSOL
    [ $nNA -gt 0 ] && printf '%-18s %d\n' NA $nNA
    [ $nCL -gt 0 ] && printf '%-18s %d\n' CL $nCL
} > topol.top