[gmx-developers] Need to export hydrogen bond indices per frame

[YOLANDA SMALL]

Hello Gmx developers,

Does anyone know which lines need to be changed in the gmx_hbond.c code to
determine which donor-hydrogen-acceptor hbond indices exist at each frame?

In other words, I would like to export lines 1483 to 1505 (see below) for every
frame where a corresponding hydrogen bond exists (lines 2069 to 2073):

------begin line 1483----------
 if (bTwo)
    fprintf(fp,"[ hbonds_%s-%s ]\n",grpnames[0],grpnames[1]);
  else
    fprintf(fp,"[ hbonds_%s ]\n",grpnames[0]);

  for(i=0; (i<hb->d.nrd); i++) {
    ddd = hb->d.don[i];
    for(k=0; (k<hb->a.nra); k++) {
      aaa = hb->a.acc[k];
      for(m=0; (m<hb->d.nhydro[i]); m++) {
        if (hb->hbmap[i][k] && ISHB(hb->hbmap[i][k]->history[m])) {
          hhh = hb->d.hydro[i][m];

          sprintf(ds,"%s",mkatomname(atoms,ddd));
          sprintf(hs,"%s",mkatomname(atoms,hhh));
          sprintf(as,"%s",mkatomname(atoms,aaa));
          fprintf(fp,"%6u %6u %6u\n",ddd+1,hhh+1,aaa+1);
          if (fplog)
            fprintf(fplog,"%12s  %12s  %12s\n",ds,hs,as);
        }
      }
    }
  }
--------end line 1505------------


--------begin line 2069----------
  for(i=0; (i<nframes); i++) {
    fprintf(fp,"%10g  %10d  %10d\n",hb->time[i],hb->nhb[i],hb->ndist[i]);
    aver_nhb  += hb->nhb[i];
    aver_dist += hb->ndist[i];
  }
--------end line 2073------------

A little more background on why I need this....I ran a simulation for 400ns on a
large protein and I need to determine hydrogen bonding pathways between certain
residues.  I wrote an external code in Fortran 90 to calculate all possible
paths between residue A and B via other protein atoms or solvent. The input for
my external code relies on a file which I hope to make g_hbond export.  The
file is a four column format, for example:

Time donor-index hydrogen-index acceptor-index
0     1           2               25
0     3           4               25
0     25          26               5
0     16          17               33
0.025 1           2               25
0.025  3           4               25
0.025 25          26               5
0.050 16          17               33
0.050 1           2               25
0.075

This represents the identity of the hydrogen bonding triad at each frame and
will be blank if no hydrogen bonds exist for a frame.  Does anyone see a way
that I can change the gmx_hbond.c code to give me the output in this format in
a single file? 

Thanks in advance for your help,
Yolanda