gmx_hbdat.bsh

echo -e "\
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   gmx_hbdat    <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>    Jicun Li    <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>>>>>>>>>>>>>>>>>>>>>>>>>>     2024-01-20 15:42:55     <<<<<<<<<<<<<<<<<<<<<<<<<\n
>>   Usage: gmx_hbdat -s     *.tpr  -n     *.ndx  -m     *.xpm
>> Default: gmx_hbdat -s topol.tpr  -n hbond.ndx  -m hbmap.xpm

--------------------------------------------------------------------------------
>> Log:
   2022-06-25: add #res for atoms
   2023-06-05: revise output
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"

tpr=topol.tpr
ndx=hbond.ndx
xpm=hbmap.xpm

gmx='gmx'        # /path/to/GMX/bin/gmx_mpi
dump="$gmx dump" # gmx dump

opt=($*); N=${#opt[@]}
for((i=0; i<N; i++)); do
	arg=${opt[$i]}; j=$((i+1)); val=${opt[$j]}
	[[ $arg =~ -s   ]] && { tpr=$val; }
	[[ $arg =~ -n   ]] && { ndx=$val; }
	[[ $arg =~ -m   ]] && { xpm=$val; }
done

$dump -s $tpr -quiet 2>>/dev/null | awk -v ndx=$ndx -v xpm=$xpm '
	BEGIN {
		isHB=0; nhb=0
		while(getline < ndx ) {
			if(index($0, "hbonds_")) isHB=1
			if(isHB && !index($1,"[")) {
				nhb++; don[nhb]=$1; hyd[nhb]=$2; acc[nhb]=$3
			}
		}
		close(ndx)

		isHB=0; nhb=0
		while(getline < xpm) {
			if(index($0, "y-axis")) isHB=1
			if(isHB && index($0, "\"")) {
				nhb++
				n=0; gsub(/[,"]/, "") #"
				for(i=1; i<=length($0); i++) if(substr($0, i, 1)=="o") n++
				occ[nhb]=n*100/length($0)
			}
		}
		close(xpm)
	}

	/#molblock/  { Ntyp=$3 }
	/moltype.+=/ { Imol=$3; getline; Nmol[Imol]=$3 }
	/moltype.+\(/ { Imol=$0; gsub(/[^0-9]/,"",Imol)
		getline txt; sub(/.*=/,"",txt); gsub(" ","_",txt)
		Name[Imol]=txt
		getline; getline txt;       gsub(/[^0-9]/,"",txt); Natm[Imol]=txt+0
		for(i=0; i<Natm[Imol]; i++) {
			getline; txt=$0; idx=$3; resID[Imol, i]=$(NF-2)+1+nres
		}
		getline
		for(i=0; i<Natm[Imol]; i++) {
			getline txt
			sub(/.+=./, "", txt); sub(/..$/, "", txt)
			Tatm[Imol, i]=txt
		}
	}

	/residue\[/ { nres++
		sub(/.*="/,"",$0); sub(/".*/,"",$0);
		resName[nres]=sprintf("%4d-%s", nres, $0)
	}

	END {

	n=split("ALA A CYS C ASP D GLU E PHE F GLY G HIS H ILE I LYS K LEU L MET M ASN N PRO P GLN Q ARG R SER S THR T VAL V TRP W TYR Y UNK X", arr)
	for(i=1; i<=n/2; i++) sname[arr[2*i-1]]=arr[2*i]

		Ntot=0; maxlen=0
		for(i=0; i<Ntyp; i++) {
			if(length(Name[i])>maxlen) maxlen=length(Name[i])
			for(n=0; n<Nmol[i]; n++) {
				for(j=0; j<Natm[i]; j++) {
					Ntot++
					Label[Ntot]=Ntot" "Name[i]" "resName[resID[i,j]]" "Tatm[i, j]
				}
			}
		}

		print "#HBond #atom                 mol.             #res Donor-H        mol.              #res Acceptor      Occupancy%   inter?  Label"
		fmt="%4d   %-20s %-"maxlen"s %8s %-9s  %-"maxlen"s %8s %-9s %9.3f    %4d    %-s\n"
		for(i=1; i<=nhb; i++) {
			split(Label[don[i]]" "Label[hyd[i]]" "Label[acc[i]], arr)
			
			n=arr[3]; sub(/[^0-9]+/, "", n)
			s=arr[3]; sub(/[0-9]+/,  "", s)
			tag=sname[s]""n"@"arr[4]"-"arr[8]
			n=arr[11]; sub(/[^0-9]+/, "", n)
			s=arr[11]; sub(/[0-9]+/,  "", s)
			tag=tag"..."sname[s]""n"@"arr[12]

			printf fmt, i,
				arr[1]"-"arr[5]"..."arr[9],
				arr[2],  arr[3], arr[4]"-"arr[8],
				arr[10], arr[11], arr[12], occ[nhb-i+1], 
				arr[3]==arr[11]? 0:1, tag
		}
	}
' > hbdat.dat