src/gromacs/gmxana/gmx_hbond.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2008, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "config.h"
  40
  41 #include <math.h>
  42
  43 #include "gromacs/commandline/pargs.h"
  44 #include "gromacs/correlationfunctions/autocorr.h"
  45 #include "gromacs/correlationfunctions/crosscorr.h"
  46 #include "gromacs/correlationfunctions/expfit.h"
  47 #include "gromacs/correlationfunctions/integrate.h"
  48 #include "gromacs/fileio/matio.h"
  49 #include "gromacs/fileio/tpxio.h"
  50 #include "gromacs/fileio/trxio.h"
  51 #include "gromacs/fileio/xvgr.h"
  52 #include "gromacs/gmxana/gmx_ana.h"
  53 #include "gromacs/legacyheaders/copyrite.h"
  54 #include "gromacs/legacyheaders/macros.h"
  55 #include "gromacs/legacyheaders/txtdump.h"
  56 #include "gromacs/legacyheaders/viewit.h"
  57 #include "gromacs/math/units.h"
  58 #include "gromacs/math/vec.h"
  59 #include "gromacs/pbcutil/pbc.h"
  60 #include "gromacs/topology/index.h"
  61 #include "gromacs/utility/cstringutil.h"
  62 #include "gromacs/utility/fatalerror.h"
  63 #include "gromacs/utility/futil.h"
  64 #include "gromacs/utility/gmxomp.h"
  65 #include "gromacs/utility/smalloc.h"
  66 #include "gromacs/utility/snprintf.h"
  67
  68 #define max_hx 7
  69 typedef int t_hx[max_hx];
  70 #define NRHXTYPES max_hx
  71 const char *hxtypenames[NRHXTYPES] =
  72 {"n-n", "n-n+1", "n-n+2", "n-n+3", "n-n+4", "n-n+5", "n-n>6"};
  73 #define MAXHH 4
  74
  75 #ifdef GMX_OPENMP
  76 #define MASTER_THREAD_ONLY(threadNr) ((threadNr) == 0)
  77 #else
  78 #define MASTER_THREAD_ONLY(threadNr) ((threadNr) == (threadNr))
  79 #endif
  80
  81 /* -----------------------------------------*/
  82
  83 enum {
  84     gr0,  gr1,    grI,  grNR
  85 };
  86 enum {
  87     hbNo, hbDist, hbHB, hbNR, hbR2
  88 };
  89 enum {
  90     noDA, ACC, DON, DA, INGROUP
  91 };
  92
  93 static const char *grpnames[grNR] = {"0", "1", "I" };
  94
  95 static gmx_bool    bDebug = FALSE;
  96
  97 #define HB_NO 0
  98 #define HB_YES 1<<0
  99 #define HB_INS 1<<1
 100 #define HB_YESINS HB_YES|HB_INS
 101 #define HB_NR (1<<2)
 102 #define MAXHYDRO 4
 103
 104 #define ISHB(h)   (((h) & 2) == 2)
 105 #define ISDIST(h) (((h) & 1) == 1)
 106 #define ISDIST2(h) (((h) & 4) == 4)
 107 #define ISACC(h)  (((h) & 1) == 1)
 108 #define ISDON(h)  (((h) & 2) == 2)
 109 #define ISINGRP(h) (((h) & 4) == 4)
 110
 111 typedef struct {
 112     int      nr;
 113     int      maxnr;
 114     atom_id *atoms;
 115 } t_ncell;
 116
 117 typedef struct {
 118     t_ncell d[grNR];
 119     t_ncell a[grNR];
 120 } t_gridcell;
 121
 122 typedef int     t_icell[grNR];
 123 typedef atom_id h_id[MAXHYDRO];
 124
 125 typedef struct {
 126     int      history[MAXHYDRO];
 127     /* Has this hbond existed ever? If so as hbDist or hbHB or both.
 128      * Result is stored as a bitmap (1 = hbDist) || (2 = hbHB)
 129      */
 130     /* Bitmask array which tells whether a hbond is present
 131      * at a given time. Either of these may be NULL
 132      */
 133     int            n0;                 /* First frame a HB was found     */
 134     int            nframes, maxframes; /* Amount of frames in this hbond */
 135     unsigned int **h;
 136     unsigned int **g;
 137     /* See Xu and Berne, JPCB 105 (2001), p. 11929. We define the
 138      * function g(t) = [1-h(t)] H(t) where H(t) is one when the donor-
 139      * acceptor distance is less than the user-specified distance (typically
 140      * 0.35 nm).
 141      */
 142 } t_hbond;
 143
 144 typedef struct {
 145     int      nra, max_nra;
 146     atom_id *acc;             /* Atom numbers of the acceptors     */
 147     int     *grp;             /* Group index                       */
 148     int     *aptr;            /* Map atom number to acceptor index */
 149 } t_acceptors;
 150
 151 typedef struct {
 152     int       nrd, max_nrd;
 153     int      *don;               /* Atom numbers of the donors         */
 154     int      *grp;               /* Group index                        */
 155     int      *dptr;              /* Map atom number to donor index     */
 156     int      *nhydro;            /* Number of hydrogens for each donor */
 157     h_id     *hydro;             /* The atom numbers of the hydrogens  */
 158     h_id     *nhbonds;           /* The number of HBs per H at current */
 159 } t_donors;
 160
 161 typedef struct {
 162     gmx_bool        bHBmap, bDAnr;
 163     int             wordlen;
 164     /* The following arrays are nframes long */
 165     int             nframes, max_frames, maxhydro;
 166     int            *nhb, *ndist;
 167     h_id           *n_bound;
 168     real           *time;
 169     t_icell        *danr;
 170     t_hx           *nhx;
 171     /* These structures are initialized from the topology at start up */
 172     t_donors        d;
 173     t_acceptors     a;
 174     /* This holds a matrix with all possible hydrogen bonds */
 175     int             nrhb, nrdist;
 176     t_hbond      ***hbmap;
 177 } t_hbdata;
 178
 179 /* Changed argument 'bMerge' into 'oneHB' below,
 180  * since -contact should cause maxhydro to be 1,
 181  * not just -merge.
 182  * - Erik Marklund May 29, 2006
 183  */
 184
 185 static t_hbdata *mk_hbdata(gmx_bool bHBmap, gmx_bool bDAnr, gmx_bool oneHB)
 186 {
 187     t_hbdata *hb;
 188
 189     snew(hb, 1);
 190     hb->wordlen = 8*sizeof(unsigned int);
 191     hb->bHBmap  = bHBmap;
 192     hb->bDAnr   = bDAnr;
 193     if (oneHB)
 194     {
 195         hb->maxhydro = 1;
 196     }
 197     else
 198     {
 199         hb->maxhydro = MAXHYDRO;
 200     }
 201     return hb;
 202 }
 203
 204 static void mk_hbmap(t_hbdata *hb)
 205 {
 206     int  i, j;
 207
 208     snew(hb->hbmap, hb->d.nrd);
 209     for (i = 0; (i < hb->d.nrd); i++)
 210     {
 211         snew(hb->hbmap[i], hb->a.nra);
 212         if (hb->hbmap[i] == NULL)
 213         {
 214             gmx_fatal(FARGS, "Could not allocate enough memory for hbmap");
 215         }
 216         for (j = 0; (j > hb->a.nra); j++)
 217         {
 218             hb->hbmap[i][j] = NULL;
 219         }
 220     }
 221 }
 222
 223 static void add_frames(t_hbdata *hb, int nframes)
 224 {
 225     int  i, j, k, l;
 226
 227     if (nframes >= hb->max_frames)
 228     {
 229         hb->max_frames += 4096;
 230         srenew(hb->time, hb->max_frames);
 231         srenew(hb->nhb, hb->max_frames);
 232         srenew(hb->ndist, hb->max_frames);
 233         srenew(hb->n_bound, hb->max_frames);
 234         srenew(hb->nhx, hb->max_frames);
 235         if (hb->bDAnr)
 236         {
 237             srenew(hb->danr, hb->max_frames);
 238         }
 239     }
 240     hb->nframes = nframes;
 241 }
 242
 243 #define OFFSET(frame) (frame / 32)
 244 #define MASK(frame)   (1 << (frame % 32))
 245
 246 static void _set_hb(unsigned int hbexist[], unsigned int frame, gmx_bool bValue)
 247 {
 248     if (bValue)
 249     {
 250         hbexist[OFFSET(frame)] |= MASK(frame);
 251     }
 252     else
 253     {
 254         hbexist[OFFSET(frame)] &= ~MASK(frame);
 255     }
 256 }
 257
 258 static gmx_bool is_hb(unsigned int hbexist[], int frame)
 259 {
 260     return ((hbexist[OFFSET(frame)] & MASK(frame)) != 0) ? 1 : 0;
 261 }
 262
 263 static void set_hb(t_hbdata *hb, int id, int ih, int ia, int frame, int ihb)
 264 {
 265     unsigned int *ghptr = NULL;
 266
 267     if (ihb == hbHB)
 268     {
 269         ghptr = hb->hbmap[id][ia]->h[ih];
 270     }
 271     else if (ihb == hbDist)
 272     {
 273         ghptr = hb->hbmap[id][ia]->g[ih];
 274     }
 275     else
 276     {
 277         gmx_fatal(FARGS, "Incomprehensible iValue %d in set_hb", ihb);
 278     }
 279
 280     _set_hb(ghptr, frame-hb->hbmap[id][ia]->n0, TRUE);
 281 }
 282
 283 static void add_ff(t_hbdata *hbd, int id, int h, int ia, int frame, int ihb)
 284 {
 285     int         i, j, n;
 286     t_hbond    *hb       = hbd->hbmap[id][ia];
 287     int         maxhydro = min(hbd->maxhydro, hbd->d.nhydro[id]);
 288     int         wlen     = hbd->wordlen;
 289     int         delta    = 32*wlen;
 290
 291     if (!hb->h[0])
 292     {
 293         hb->n0        = frame;
 294         hb->maxframes = delta;
 295         for (i = 0; (i < maxhydro); i++)
 296         {
 297             snew(hb->h[i], hb->maxframes/wlen);
 298             snew(hb->g[i], hb->maxframes/wlen);
 299         }
 300     }
 301     else
 302     {
 303         hb->nframes = frame-hb->n0;
 304         /* We need a while loop here because hbonds may be returning
 305          * after a long time.
 306          */
 307         while (hb->nframes >= hb->maxframes)
 308         {
 309             n = hb->maxframes + delta;
 310             for (i = 0; (i < maxhydro); i++)
 311             {
 312                 srenew(hb->h[i], n/wlen);
 313                 srenew(hb->g[i], n/wlen);
 314                 for (j = hb->maxframes/wlen; (j < n/wlen); j++)
 315                 {
 316                     hb->h[i][j] = 0;
 317                     hb->g[i][j] = 0;
 318                 }
 319             }
 320
 321             hb->maxframes = n;
 322         }
 323     }
 324     if (frame >= 0)
 325     {
 326         set_hb(hbd, id, h, ia, frame, ihb);
 327     }
 328
 329 }
 330
 331 static void inc_nhbonds(t_donors *ddd, int d, int h)
 332 {
 333     int j;
 334     int dptr = ddd->dptr[d];
 335
 336     for (j = 0; (j < ddd->nhydro[dptr]); j++)
 337     {
 338         if (ddd->hydro[dptr][j] == h)
 339         {
 340             ddd->nhbonds[dptr][j]++;
 341             break;
 342         }
 343     }
 344     if (j == ddd->nhydro[dptr])
 345     {
 346         gmx_fatal(FARGS, "No such hydrogen %d on donor %d\n", h+1, d+1);
 347     }
 348 }
 349
 350 static int _acceptor_index(t_acceptors *a, int grp, atom_id i,
 351                            const char *file, int line)
 352 {
 353     int ai = a->aptr[i];
 354
 355     if (a->grp[ai] != grp)
 356     {
 357         if (debug && bDebug)
 358         {
 359             fprintf(debug, "Acc. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
 360                     ai, a->grp[ai], grp, file, line);
 361         }
 362         return NOTSET;
 363     }
 364     else
 365     {
 366         return ai;
 367     }
 368 }
 369 #define acceptor_index(a, grp, i) _acceptor_index(a, grp, i, __FILE__, __LINE__)
 370
 371 static int _donor_index(t_donors *d, int grp, atom_id i, const char *file, int line)
 372 {
 373     int di = d->dptr[i];
 374
 375     if (di == NOTSET)
 376     {
 377         return NOTSET;
 378     }
 379
 380     if (d->grp[di] != grp)
 381     {
 382         if (debug && bDebug)
 383         {
 384             fprintf(debug, "Don. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
 385                     di, d->grp[di], grp, file, line);
 386         }
 387         return NOTSET;
 388     }
 389     else
 390     {
 391         return di;
 392     }
 393 }
 394 #define donor_index(d, grp, i) _donor_index(d, grp, i, __FILE__, __LINE__)
 395
 396 static gmx_bool isInterchangable(t_hbdata *hb, int d, int a, int grpa, int grpd)
 397 {
 398     /* g_hbond doesn't allow overlapping groups */
 399     if (grpa != grpd)
 400     {
 401         return FALSE;
 402     }
 403     return
 404         donor_index(&hb->d, grpd, a) != NOTSET
 405         && acceptor_index(&hb->a, grpa, d) != NOTSET;
 406 }
 407
 408
 409 static void add_hbond(t_hbdata *hb, int d, int a, int h, int grpd, int grpa,
 410                       int frame, gmx_bool bMerge, int ihb, gmx_bool bContact)
 411 {
 412     int      k, id, ia, hh;
 413     gmx_bool daSwap = FALSE;
 414
 415     if ((id = hb->d.dptr[d]) == NOTSET)
 416     {
 417         gmx_fatal(FARGS, "No donor atom %d", d+1);
 418     }
 419     else if (grpd != hb->d.grp[id])
 420     {
 421         gmx_fatal(FARGS, "Inconsistent donor groups, %d iso %d, atom %d",
 422                   grpd, hb->d.grp[id], d+1);
 423     }
 424     if ((ia = hb->a.aptr[a]) == NOTSET)
 425     {
 426         gmx_fatal(FARGS, "No acceptor atom %d", a+1);
 427     }
 428     else if (grpa != hb->a.grp[ia])
 429     {
 430         gmx_fatal(FARGS, "Inconsistent acceptor groups, %d iso %d, atom %d",
 431                   grpa, hb->a.grp[ia], a+1);
 432     }
 433
 434     if (bMerge)
 435     {
 436
 437         if (isInterchangable(hb, d, a, grpd, grpa) && d > a)
 438         /* Then swap identity so that the id of d is lower then that of a.
 439          *
 440          * This should really be redundant by now, as is_hbond() now ought to return
 441          * hbNo in the cases where this conditional is TRUE. */
 442         {
 443             daSwap = TRUE;
 444             k      = d;
 445             d      = a;
 446             a      = k;
 447
 448             /* Now repeat donor/acc check. */
 449             if ((id = hb->d.dptr[d]) == NOTSET)
 450             {
 451                 gmx_fatal(FARGS, "No donor atom %d", d+1);
 452             }
 453             else if (grpd != hb->d.grp[id])
 454             {
 455                 gmx_fatal(FARGS, "Inconsistent donor groups, %d iso %d, atom %d",
 456                           grpd, hb->d.grp[id], d+1);
 457             }
 458             if ((ia = hb->a.aptr[a]) == NOTSET)
 459             {
 460                 gmx_fatal(FARGS, "No acceptor atom %d", a+1);
 461             }
 462             else if (grpa != hb->a.grp[ia])
 463             {
 464                 gmx_fatal(FARGS, "Inconsistent acceptor groups, %d iso %d, atom %d",
 465                           grpa, hb->a.grp[ia], a+1);
 466             }
 467         }
 468     }
 469
 470     if (hb->hbmap)
 471     {
 472         /* Loop over hydrogens to find which hydrogen is in this particular HB */
 473         if ((ihb == hbHB) && !bMerge && !bContact)
 474         {
 475             for (k = 0; (k < hb->d.nhydro[id]); k++)
 476             {
 477                 if (hb->d.hydro[id][k] == h)
 478                 {
 479                     break;
 480                 }
 481             }
 482             if (k == hb->d.nhydro[id])
 483             {
 484                 gmx_fatal(FARGS, "Donor %d does not have hydrogen %d (a = %d)",
 485                           d+1, h+1, a+1);
 486             }
 487         }
 488         else
 489         {
 490             k = 0;
 491         }
 492
 493         if (hb->bHBmap)
 494         {
 495
 496 #pragma omp critical
 497             {
 498                 if (hb->hbmap[id][ia] == NULL)
 499                 {
 500                     snew(hb->hbmap[id][ia], 1);
 501                     snew(hb->hbmap[id][ia]->h, hb->maxhydro);
 502                     snew(hb->hbmap[id][ia]->g, hb->maxhydro);
 503                 }
 504                 add_ff(hb, id, k, ia, frame, ihb);
 505             }
 506         }
 507
 508         /* Strange construction with frame >=0 is a relic from old code
 509          * for selected hbond analysis. It may be necessary again if that
 510          * is made to work again.
 511          */
 512         if (frame >= 0)
 513         {
 514             hh = hb->hbmap[id][ia]->history[k];
 515             if (ihb == hbHB)
 516             {
 517                 hb->nhb[frame]++;
 518                 if (!(ISHB(hh)))
 519                 {
 520                     hb->hbmap[id][ia]->history[k] = hh | 2;
 521                     hb->nrhb++;
 522                 }
 523             }
 524             else
 525             {
 526                 if (ihb == hbDist)
 527                 {
 528                     hb->ndist[frame]++;
 529                     if (!(ISDIST(hh)))
 530                     {
 531                         hb->hbmap[id][ia]->history[k] = hh | 1;
 532                         hb->nrdist++;
 533                     }
 534                 }
 535             }
 536         }
 537     }
 538     else
 539     {
 540         if (frame >= 0)
 541         {
 542             if (ihb == hbHB)
 543             {
 544                 hb->nhb[frame]++;
 545             }
 546             else
 547             {
 548                 if (ihb == hbDist)
 549                 {
 550                     hb->ndist[frame]++;
 551                 }
 552             }
 553         }
 554     }
 555     if (bMerge && daSwap)
 556     {
 557         h = hb->d.hydro[id][0];
 558     }
 559     /* Increment number if HBonds per H */
 560     if (ihb == hbHB && !bContact)
 561     {
 562         inc_nhbonds(&(hb->d), d, h);
 563     }
 564 }
 565
 566 static char *mkatomname(t_atoms *atoms, int i)
 567 {
 568     static char buf[32];
 569     int         rnr;
 570
 571     rnr = atoms->atom[i].resind;
 572     sprintf(buf, "%4s%d%-4s",
 573             *atoms->resinfo[rnr].name, atoms->resinfo[rnr].nr, *atoms->atomname[i]);
 574
 575     return buf;
 576 }
 577
 578 static void gen_datable(atom_id *index, int isize, unsigned char *datable, int natoms)
 579 {
 580     /* Generates table of all atoms and sets the ingroup bit for atoms in index[] */
 581     int i;
 582
 583     for (i = 0; i < isize; i++)
 584     {
 585         if (index[i] >= natoms)
 586         {
 587             gmx_fatal(FARGS, "Atom has index %d larger than number of atoms %d.", index[i], natoms);
 588         }
 589         datable[index[i]] |= INGROUP;
 590     }
 591 }
 592
 593 static void clear_datable_grp(unsigned char *datable, int size)
 594 {
 595     /* Clears group information from the table */
 596     int        i;
 597     const char mask = !(char)INGROUP;
 598     if (size > 0)
 599     {
 600         for (i = 0; i < size; i++)
 601         {
 602             datable[i] &= mask;
 603         }
 604     }
 605 }
 606
 607 static void add_acc(t_acceptors *a, int ia, int grp)
 608 {
 609     if (a->nra >= a->max_nra)
 610     {
 611         a->max_nra += 16;
 612         srenew(a->acc, a->max_nra);
 613         srenew(a->grp, a->max_nra);
 614     }
 615     a->grp[a->nra]   = grp;
 616     a->acc[a->nra++] = ia;
 617 }
 618
 619 static void search_acceptors(t_topology *top, int isize,
 620                              atom_id *index, t_acceptors *a, int grp,
 621                              gmx_bool bNitAcc,
 622                              gmx_bool bContact, gmx_bool bDoIt, unsigned char *datable)
 623 {
 624     int i, n;
 625
 626     if (bDoIt)
 627     {
 628         for (i = 0; (i < isize); i++)
 629         {
 630             n = index[i];
 631             if ((bContact ||
 632                  (((*top->atoms.atomname[n])[0] == 'O') ||
 633                   (bNitAcc && ((*top->atoms.atomname[n])[0] == 'N')))) &&
 634                 ISINGRP(datable[n]))
 635             {
 636                 datable[n] |= ACC; /* set the atom's acceptor flag in datable. */
 637                 add_acc(a, n, grp);
 638             }
 639         }
 640     }
 641     snew(a->aptr, top->atoms.nr);
 642     for (i = 0; (i < top->atoms.nr); i++)
 643     {
 644         a->aptr[i] = NOTSET;
 645     }
 646     for (i = 0; (i < a->nra); i++)
 647     {
 648         a->aptr[a->acc[i]] = i;
 649     }
 650 }
 651
 652 static void add_h2d(int id, int ih, t_donors *ddd)
 653 {
 654     int i;
 655
 656     for (i = 0; (i < ddd->nhydro[id]); i++)
 657     {
 658         if (ddd->hydro[id][i] == ih)
 659         {
 660             printf("Hm. This isn't the first time I found this donor (%d,%d)\n",
 661                    ddd->don[id], ih);
 662             break;
 663         }
 664     }
 665     if (i == ddd->nhydro[id])
 666     {
 667         if (ddd->nhydro[id] >= MAXHYDRO)
 668         {
 669             gmx_fatal(FARGS, "Donor %d has more than %d hydrogens!",
 670                       ddd->don[id], MAXHYDRO);
 671         }
 672         ddd->hydro[id][i] = ih;
 673         ddd->nhydro[id]++;
 674     }
 675 }
 676
 677 static void add_dh(t_donors *ddd, int id, int ih, int grp, unsigned char *datable)
 678 {
 679     int i;
 680
 681     if (!datable || ISDON(datable[id]))
 682     {
 683         if (ddd->dptr[id] == NOTSET)   /* New donor */
 684         {
 685             i             = ddd->nrd;
 686             ddd->dptr[id] = i;
 687         }
 688         else
 689         {
 690             i = ddd->dptr[id];
 691         }
 692
 693         if (i == ddd->nrd)
 694         {
 695             if (ddd->nrd >= ddd->max_nrd)
 696             {
 697                 ddd->max_nrd += 128;
 698                 srenew(ddd->don, ddd->max_nrd);
 699                 srenew(ddd->nhydro, ddd->max_nrd);
 700                 srenew(ddd->hydro, ddd->max_nrd);
 701                 srenew(ddd->nhbonds, ddd->max_nrd);
 702                 srenew(ddd->grp, ddd->max_nrd);
 703             }
 704             ddd->don[ddd->nrd]    = id;
 705             ddd->nhydro[ddd->nrd] = 0;
 706             ddd->grp[ddd->nrd]    = grp;
 707             ddd->nrd++;
 708         }
 709         else
 710         {
 711             ddd->don[i] = id;
 712         }
 713         add_h2d(i, ih, ddd);
 714     }
 715     else
 716     if (datable)
 717     {
 718         printf("Warning: Atom %d is not in the d/a-table!\n", id);
 719     }
 720 }
 721
 722 static void search_donors(t_topology *top, int isize, atom_id *index,
 723                           t_donors *ddd, int grp, gmx_bool bContact, gmx_bool bDoIt,
 724                           unsigned char *datable)
 725 {
 726     int            i, j, nra, n;
 727     t_functype     func_type;
 728     t_ilist       *interaction;
 729     atom_id        nr1, nr2, nr3;
 730     gmx_bool       stop;
 731
 732     if (!ddd->dptr)
 733     {
 734         snew(ddd->dptr, top->atoms.nr);
 735         for (i = 0; (i < top->atoms.nr); i++)
 736         {
 737             ddd->dptr[i] = NOTSET;
 738         }
 739     }
 740
 741     if (bContact)
 742     {
 743         if (bDoIt)
 744         {
 745             for (i = 0; (i < isize); i++)
 746             {
 747                 datable[index[i]] |= DON;
 748                 add_dh(ddd, index[i], -1, grp, datable);
 749             }
 750         }
 751     }
 752     else
 753     {
 754         for (func_type = 0; (func_type < F_NRE); func_type++)
 755         {
 756             interaction = &(top->idef.il[func_type]);
 757             if (func_type == F_POSRES || func_type == F_FBPOSRES)
 758             {
 759                 /* The ilist looks strange for posre. Bug in grompp?
 760                  * We don't need posre interactions for hbonds anyway.*/
 761                 continue;
 762             }
 763             for (i = 0; i < interaction->nr;
 764                  i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1)
 765             {
 766                 /* next function */
 767                 if (func_type != top->idef.functype[interaction->iatoms[i]])
 768                 {
 769                     fprintf(stderr, "Error in func_type %s",
 770                             interaction_function[func_type].longname);
 771                     continue;
 772                 }
 773
 774                 /* check out this functype */
 775                 if (func_type == F_SETTLE)
 776                 {
 777                     nr1 = interaction->iatoms[i+1];
 778                     nr2 = interaction->iatoms[i+2];
 779                     nr3 = interaction->iatoms[i+3];
 780
 781                     if (ISINGRP(datable[nr1]))
 782                     {
 783                         if (ISINGRP(datable[nr2]))
 784                         {
 785                             datable[nr1] |= DON;
 786                             add_dh(ddd, nr1, nr1+1, grp, datable);
 787                         }
 788                         if (ISINGRP(datable[nr3]))
 789                         {
 790                             datable[nr1] |= DON;
 791                             add_dh(ddd, nr1, nr1+2, grp, datable);
 792                         }
 793                     }
 794                 }
 795                 else if (IS_CHEMBOND(func_type))
 796                 {
 797                     for (j = 0; j < 2; j++)
 798                     {
 799                         nr1 = interaction->iatoms[i+1+j];
 800                         nr2 = interaction->iatoms[i+2-j];
 801                         if ((*top->atoms.atomname[nr1][0] == 'H') &&
 802                             ((*top->atoms.atomname[nr2][0] == 'O') ||
 803                              (*top->atoms.atomname[nr2][0] == 'N')) &&
 804                             ISINGRP(datable[nr1]) && ISINGRP(datable[nr2]))
 805                         {
 806                             datable[nr2] |= DON;
 807                             add_dh(ddd, nr2, nr1, grp, datable);
 808                         }
 809                     }
 810                 }
 811             }
 812         }
 813 #ifdef SAFEVSITES
 814         for (func_type = 0; func_type < F_NRE; func_type++)
 815         {
 816             interaction = &top->idef.il[func_type];
 817             for (i = 0; i < interaction->nr;
 818                  i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1)
 819             {
 820                 /* next function */
 821                 if (func_type != top->idef.functype[interaction->iatoms[i]])
 822                 {
 823                     gmx_incons("function type in search_donors");
 824                 }
 825
 826                 if (interaction_function[func_type].flags & IF_VSITE)
 827                 {
 828                     nr1 = interaction->iatoms[i+1];
 829                     if (*top->atoms.atomname[nr1][0]  == 'H')
 830                     {
 831                         nr2  = nr1-1;
 832                         stop = FALSE;
 833                         while (!stop && ( *top->atoms.atomname[nr2][0] == 'H'))
 834                         {
 835                             if (nr2)
 836                             {
 837                                 nr2--;
 838                             }
 839                             else
 840                             {
 841                                 stop = TRUE;
 842                             }
 843                         }
 844                         if (!stop && ( ( *top->atoms.atomname[nr2][0] == 'O') ||
 845                                        ( *top->atoms.atomname[nr2][0] == 'N') ) &&
 846                             ISINGRP(datable[nr1]) && ISINGRP(datable[nr2]))
 847                         {
 848                             datable[nr2] |= DON;
 849                             add_dh(ddd, nr2, nr1, grp, datable);
 850                         }
 851                     }
 852                 }
 853             }
 854         }
 855 #endif
 856     }
 857 }
 858
 859 static t_gridcell ***init_grid(gmx_bool bBox, rvec box[], real rcut, ivec ngrid)
 860 {
 861     t_gridcell ***grid;
 862     int           i, y, z;
 863
 864     if (bBox)
 865     {
 866         for (i = 0; i < DIM; i++)
 867         {
 868             ngrid[i] = (box[i][i]/(1.2*rcut));
 869         }
 870     }
 871
 872     if (!bBox || (ngrid[XX] < 3) || (ngrid[YY] < 3) || (ngrid[ZZ] < 3) )
 873     {
 874         for (i = 0; i < DIM; i++)
 875         {
 876             ngrid[i] = 1;
 877         }
 878     }
 879     else
 880     {
 881         printf("\nWill do grid-seach on %dx%dx%d grid, rcut=%g\n",
 882                ngrid[XX], ngrid[YY], ngrid[ZZ], rcut);
 883     }
 884     snew(grid, ngrid[ZZ]);
 885     for (z = 0; z < ngrid[ZZ]; z++)
 886     {
 887         snew((grid)[z], ngrid[YY]);
 888         for (y = 0; y < ngrid[YY]; y++)
 889         {
 890             snew((grid)[z][y], ngrid[XX]);
 891         }
 892     }
 893     return grid;
 894 }
 895
 896 static void reset_nhbonds(t_donors *ddd)
 897 {
 898     int i, j;
 899
 900     for (i = 0; (i < ddd->nrd); i++)
 901     {
 902         for (j = 0; (j < MAXHH); j++)
 903         {
 904             ddd->nhbonds[i][j] = 0;
 905         }
 906     }
 907 }
 908
 909 void pbc_correct_gem(rvec dx, matrix box, rvec hbox);
 910 void pbc_in_gridbox(rvec dx, matrix box);
 911
 912 static void build_grid(t_hbdata *hb, rvec x[], rvec xshell,
 913                        gmx_bool bBox, matrix box, rvec hbox,
 914                        real rcut, real rshell,
 915                        ivec ngrid, t_gridcell ***grid)
 916 {
 917     int         i, m, gr, xi, yi, zi, nr;
 918     atom_id    *ad;
 919     ivec        grididx;
 920     rvec        invdelta, dshell, xtemp = {0, 0, 0};
 921     t_ncell    *newgrid;
 922     gmx_bool    bDoRshell, bInShell, bAcc;
 923     real        rshell2 = 0;
 924     int         gx, gy, gz;
 925     int         dum = -1;
 926
 927     bDoRshell = (rshell > 0);
 928     rshell2   = sqr(rshell);
 929     bInShell  = TRUE;
 930
 931 #define DBB(x) if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n", __LINE__,#x, x)
 932     DBB(dum);
 933     for (m = 0; m < DIM; m++)
 934     {
 935         hbox[m] = box[m][m]*0.5;
 936         if (bBox)
 937         {
 938             invdelta[m] = ngrid[m]/box[m][m];
 939             if (1/invdelta[m] < rcut)
 940             {
 941                 gmx_fatal(FARGS, "Your computational box has shrunk too much.\n"
 942                           "%s can not handle this situation, sorry.\n",
 943                           ShortProgram());
 944             }
 945         }
 946         else
 947         {
 948             invdelta[m] = 0;
 949         }
 950     }
 951     grididx[XX] = 0;
 952     grididx[YY] = 0;
 953     grididx[ZZ] = 0;
 954     DBB(dum);
 955     /* resetting atom counts */
 956     for (gr = 0; (gr < grNR); gr++)
 957     {
 958         for (zi = 0; zi < ngrid[ZZ]; zi++)
 959         {
 960             for (yi = 0; yi < ngrid[YY]; yi++)
 961             {
 962                 for (xi = 0; xi < ngrid[XX]; xi++)
 963                 {
 964                     grid[zi][yi][xi].d[gr].nr = 0;
 965                     grid[zi][yi][xi].a[gr].nr = 0;
 966                 }
 967             }
 968         }
 969         DBB(dum);
 970
 971         /* put atoms in grid cells */
 972         for (bAcc = FALSE; (bAcc <= TRUE); bAcc++)
 973         {
 974             if (bAcc)
 975             {
 976                 nr = hb->a.nra;
 977                 ad = hb->a.acc;
 978             }
 979             else
 980             {
 981                 nr = hb->d.nrd;
 982                 ad = hb->d.don;
 983             }
 984             DBB(bAcc);
 985             for (i = 0; (i < nr); i++)
 986             {
 987                 /* check if we are inside the shell */
 988                 /* if bDoRshell=FALSE then bInShell=TRUE always */
 989                 DBB(i);
 990                 if (bDoRshell)
 991                 {
 992                     bInShell = TRUE;
 993                     rvec_sub(x[ad[i]], xshell, dshell);
 994                     if (bBox)
 995                     {
 996                         gmx_bool bDone = FALSE;
 997                         while (!bDone)
 998                         {
 999                             bDone = TRUE;
1000                             for (m = DIM-1; m >= 0 && bInShell; m--)
1001                             {
1002                                 if (dshell[m] < -hbox[m])
1003                                 {
1004                                     bDone = FALSE;
1005                                     rvec_inc(dshell, box[m]);
1006                                 }
1007                                 if (dshell[m] >= hbox[m])
1008                                 {
1009                                     bDone      = FALSE;
1010                                     dshell[m] -= 2*hbox[m];
1011                                 }
1012                             }
1013                         }
1014                         for (m = DIM-1; m >= 0 && bInShell; m--)
1015                         {
1016                             /* if we're outside the cube, we're outside the sphere also! */
1017                             if ( (dshell[m] > rshell) || (-dshell[m] > rshell) )
1018                             {
1019                                 bInShell = FALSE;
1020                             }
1021                         }
1022                     }
1023                     /* if we're inside the cube, check if we're inside the sphere */
1024                     if (bInShell)
1025                     {
1026                         bInShell = norm2(dshell) < rshell2;
1027                     }
1028                 }
1029                 DBB(i);
1030                 if (bInShell)
1031                 {
1032                     if (bBox)
1033                     {
1034                         pbc_in_gridbox(x[ad[i]], box);
1035
1036                         for (m = DIM-1; m >= 0; m--)
1037                         {   /* determine grid index of atom */
1038                             grididx[m] = x[ad[i]][m]*invdelta[m];
1039                             grididx[m] = (grididx[m]+ngrid[m]) % ngrid[m];
1040                         }
1041                     }
1042
1043                     gx = grididx[XX];
1044                     gy = grididx[YY];
1045                     gz = grididx[ZZ];
1046                     range_check(gx, 0, ngrid[XX]);
1047                     range_check(gy, 0, ngrid[YY]);
1048                     range_check(gz, 0, ngrid[ZZ]);
1049                     DBB(gx);
1050                     DBB(gy);
1051                     DBB(gz);
1052                     /* add atom to grid cell */
1053                     if (bAcc)
1054                     {
1055                         newgrid = &(grid[gz][gy][gx].a[gr]);
1056                     }
1057                     else
1058                     {
1059                         newgrid = &(grid[gz][gy][gx].d[gr]);
1060                     }
1061                     if (newgrid->nr >= newgrid->maxnr)
1062                     {
1063                         newgrid->maxnr += 10;
1064                         DBB(newgrid->maxnr);
1065                         srenew(newgrid->atoms, newgrid->maxnr);
1066                     }
1067                     DBB(newgrid->nr);
1068                     newgrid->atoms[newgrid->nr] = ad[i];
1069                     newgrid->nr++;
1070                 }
1071             }
1072         }
1073     }
1074 }
1075
1076 static void count_da_grid(ivec ngrid, t_gridcell ***grid, t_icell danr)
1077 {
1078     int gr, xi, yi, zi;
1079
1080     for (gr = 0; (gr < grNR); gr++)
1081     {
1082         danr[gr] = 0;
1083         for (zi = 0; zi < ngrid[ZZ]; zi++)
1084         {
1085             for (yi = 0; yi < ngrid[YY]; yi++)
1086             {
1087                 for (xi = 0; xi < ngrid[XX]; xi++)
1088                 {
1089                     danr[gr] += grid[zi][yi][xi].d[gr].nr;
1090                 }
1091             }
1092         }
1093     }
1094 }
1095
1096 /* The grid loop.
1097  * Without a box, the grid is 1x1x1, so all loops are 1 long.
1098  * With a rectangular box (bTric==FALSE) all loops are 3 long.
1099  * With a triclinic box all loops are 3 long, except when a cell is
1100  * located next to one of the box edges which is not parallel to the
1101  * x/y-plane, in that case all cells in a line or layer are searched.
1102  * This could be implemented slightly more efficient, but the code
1103  * would get much more complicated.
1104  */
1105 static gmx_inline gmx_bool grid_loop_begin(int n, int x, gmx_bool bTric, gmx_bool bEdge)
1106 {
1107     return ((n == 1) ? x : bTric && bEdge ? 0     : (x-1));
1108 }
1109 static gmx_inline gmx_bool grid_loop_end(int n, int x, gmx_bool bTric, gmx_bool bEdge)
1110 {
1111     return ((n == 1) ? x : bTric && bEdge ? (n-1) : (x+1));
1112 }
1113 static gmx_inline int grid_mod(int j, int n)
1114 {
1115     return (j+n) % (n);
1116 }
1117
1118 static void dump_grid(FILE *fp, ivec ngrid, t_gridcell ***grid)
1119 {
1120     int gr, x, y, z, sum[grNR];
1121
1122     fprintf(fp, "grid %dx%dx%d\n", ngrid[XX], ngrid[YY], ngrid[ZZ]);
1123     for (gr = 0; gr < grNR; gr++)
1124     {
1125         sum[gr] = 0;
1126         fprintf(fp, "GROUP %d (%s)\n", gr, grpnames[gr]);
1127         for (z = 0; z < ngrid[ZZ]; z += 2)
1128         {
1129             fprintf(fp, "Z=%d,%d\n", z, z+1);
1130             for (y = 0; y < ngrid[YY]; y++)
1131             {
1132                 for (x = 0; x < ngrid[XX]; x++)
1133                 {
1134                     fprintf(fp, "%3d", grid[x][y][z].d[gr].nr);
1135                     sum[gr] += grid[z][y][x].d[gr].nr;
1136                     fprintf(fp, "%3d", grid[x][y][z].a[gr].nr);
1137                     sum[gr] += grid[z][y][x].a[gr].nr;
1138
1139                 }
1140                 fprintf(fp, " | ");
1141                 if ( (z+1) < ngrid[ZZ])
1142                 {
1143                     for (x = 0; x < ngrid[XX]; x++)
1144                     {
1145                         fprintf(fp, "%3d", grid[z+1][y][x].d[gr].nr);
1146                         sum[gr] += grid[z+1][y][x].d[gr].nr;
1147                         fprintf(fp, "%3d", grid[z+1][y][x].a[gr].nr);
1148                         sum[gr] += grid[z+1][y][x].a[gr].nr;
1149                     }
1150                 }
1151                 fprintf(fp, "\n");
1152             }
1153         }
1154     }
1155     fprintf(fp, "TOTALS:");
1156     for (gr = 0; gr < grNR; gr++)
1157     {
1158         fprintf(fp, " %d=%d", gr, sum[gr]);
1159     }
1160     fprintf(fp, "\n");
1161 }
1162
1163 /* New GMX record! 5 * in a row. Congratulations!
1164  * Sorry, only four left.
1165  */
1166 static void free_grid(ivec ngrid, t_gridcell ****grid)
1167 {
1168     int           y, z;
1169     t_gridcell ***g = *grid;
1170
1171     for (z = 0; z < ngrid[ZZ]; z++)
1172     {
1173         for (y = 0; y < ngrid[YY]; y++)
1174         {
1175             sfree(g[z][y]);
1176         }
1177         sfree(g[z]);
1178     }
1179     sfree(g);
1180     g = NULL;
1181 }
1182
1183 void pbc_correct_gem(rvec dx, matrix box, rvec hbox)
1184 {
1185     int      m;
1186     gmx_bool bDone = FALSE;
1187     while (!bDone)
1188     {
1189         bDone = TRUE;
1190         for (m = DIM-1; m >= 0; m--)
1191         {
1192             if (dx[m] < -hbox[m])
1193             {
1194                 bDone = FALSE;
1195                 rvec_inc(dx, box[m]);
1196             }
1197             if (dx[m] >= hbox[m])
1198             {
1199                 bDone = FALSE;
1200                 rvec_dec(dx, box[m]);
1201             }
1202         }
1203     }
1204 }
1205
1206 void pbc_in_gridbox(rvec dx, matrix box)
1207 {
1208     int      m;
1209     gmx_bool bDone = FALSE;
1210     while (!bDone)
1211     {
1212         bDone = TRUE;
1213         for (m = DIM-1; m >= 0; m--)
1214         {
1215             if (dx[m] < 0)
1216             {
1217                 bDone = FALSE;
1218                 rvec_inc(dx, box[m]);
1219             }
1220             if (dx[m] >= box[m][m])
1221             {
1222                 bDone = FALSE;
1223                 rvec_dec(dx, box[m]);
1224             }
1225         }
1226     }
1227 }
1228
1229 /* Added argument r2cut, changed contact and implemented
1230  * use of second cut-off.
1231  * - Erik Marklund, June 29, 2006
1232  */
1233 static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a,
1234                     real rcut, real r2cut, real ccut,
1235                     rvec x[], gmx_bool bBox, matrix box, rvec hbox,
1236                     real *d_ha, real *ang, gmx_bool bDA, int *hhh,
1237                     gmx_bool bContact, gmx_bool bMerge)
1238 {
1239     int      h, hh, id, ja, ihb;
1240     rvec     r_da, r_ha, r_dh, r = {0, 0, 0};
1241     ivec     ri;
1242     real     rc2, r2c2, rda2, rha2, ca;
1243     gmx_bool HAinrange = FALSE; /* If !bDA. Needed for returning hbDist in a correct way. */
1244     gmx_bool daSwap    = FALSE;
1245
1246     if (d == a)
1247     {
1248         return hbNo;
1249     }
1250
1251     if (((id = donor_index(&hb->d, grpd, d)) == NOTSET) ||
1252         ((ja = acceptor_index(&hb->a, grpa, a)) == NOTSET))
1253     {
1254         return hbNo;
1255     }
1256
1257     rc2  = rcut*rcut;
1258     r2c2 = r2cut*r2cut;
1259
1260     rvec_sub(x[d], x[a], r_da);
1261     /* Insert projection code here */
1262
1263     if (bMerge && d > a && isInterchangable(hb, d, a, grpd, grpa))
1264     {
1265         /* Then this hbond/contact will be found again, or it has already been found. */
1266         /*return hbNo;*/
1267     }
1268     if (bBox)
1269     {
1270         if (d > a && bMerge && isInterchangable(hb, d, a, grpd, grpa)) /* acceptor is also a donor and vice versa? */
1271         {                                                              /* return hbNo; */
1272             daSwap = TRUE;                                             /* If so, then their history should be filed with donor and acceptor swapped. */
1273         }
1274         pbc_correct_gem(r_da, box, hbox);
1275     }
1276     rda2 = iprod(r_da, r_da);
1277
1278     if (bContact)
1279     {
1280         if (daSwap && grpa == grpd)
1281         {
1282             return hbNo;
1283         }
1284         if (rda2 <= rc2)
1285         {
1286             return hbHB;
1287         }
1288         else if (rda2 < r2c2)
1289         {
1290             return hbDist;
1291         }
1292         else
1293         {
1294             return hbNo;
1295         }
1296     }
1297     *hhh = NOTSET;
1298
1299     if (bDA && (rda2 > rc2))
1300     {
1301         return hbNo;
1302     }
1303
1304     for (h = 0; (h < hb->d.nhydro[id]); h++)
1305     {
1306         hh   = hb->d.hydro[id][h];
1307         rha2 = rc2+1;
1308         if (!bDA)
1309         {
1310             rvec_sub(x[hh], x[a], r_ha);
1311             if (bBox)
1312             {
1313                 pbc_correct_gem(r_ha, box, hbox);
1314             }
1315             rha2 = iprod(r_ha, r_ha);
1316         }
1317
1318         if (bDA || (!bDA && (rha2 <= rc2)))
1319         {
1320             rvec_sub(x[d], x[hh], r_dh);
1321             if (bBox)
1322             {
1323                 pbc_correct_gem(r_dh, box, hbox);
1324             }
1325
1326             if (!bDA)
1327             {
1328                 HAinrange = TRUE;
1329             }
1330             ca = cos_angle(r_dh, r_da);
1331             /* if angle is smaller, cos is larger */
1332             if (ca >= ccut)
1333             {
1334                 *hhh  = hh;
1335                 *d_ha = sqrt(bDA ? rda2 : rha2);
1336                 *ang  = acos(ca);
1337                 return hbHB;
1338             }
1339         }
1340     }
1341     if (bDA || (!bDA && HAinrange))
1342     {
1343         return hbDist;
1344     }
1345     else
1346     {
1347         return hbNo;
1348     }
1349 }
1350
1351 /* Merging is now done on the fly, so do_merge is most likely obsolete now.
1352  * Will do some more testing before removing the function entirely.
1353  * - Erik Marklund, MAY 10 2010 */
1354 static void do_merge(t_hbdata *hb, int ntmp,
1355                      unsigned int htmp[], unsigned int gtmp[],
1356                      t_hbond *hb0, t_hbond *hb1)
1357 {
1358     /* Here we need to make sure we're treating periodicity in
1359      * the right way for the geminate recombination kinetics. */
1360
1361     int       m, mm, n00, n01, nn0, nnframes;
1362
1363     /* Decide where to start from when merging */
1364     n00      = hb0->n0;
1365     n01      = hb1->n0;
1366     nn0      = min(n00, n01);
1367     nnframes = max(n00 + hb0->nframes, n01 + hb1->nframes) - nn0;
1368     /* Initiate tmp arrays */
1369     for (m = 0; (m < ntmp); m++)
1370     {
1371         htmp[m] = 0;
1372         gtmp[m] = 0;
1373     }
1374     /* Fill tmp arrays with values due to first HB */
1375     /* Once again '<' had to be replaced with '<='
1376        to catch the last frame in which the hbond
1377        appears.
1378        - Erik Marklund, June 1, 2006 */
1379     for (m = 0; (m <= hb0->nframes); m++)
1380     {
1381         mm       = m+n00-nn0;
1382         htmp[mm] = is_hb(hb0->h[0], m);
1383     }
1384     for (m = 0; (m <= hb0->nframes); m++)
1385     {
1386         mm       = m+n00-nn0;
1387         gtmp[mm] = is_hb(hb0->g[0], m);
1388     }
1389     /* Next HB */
1390     for (m = 0; (m <= hb1->nframes); m++)
1391     {
1392         mm       = m+n01-nn0;
1393         htmp[mm] = htmp[mm] || is_hb(hb1->h[0], m);
1394         gtmp[mm] = gtmp[mm] || is_hb(hb1->g[0], m);
1395     }
1396     /* Reallocate target array */
1397     if (nnframes > hb0->maxframes)
1398     {
1399         srenew(hb0->h[0], 4+nnframes/hb->wordlen);
1400         srenew(hb0->g[0], 4+nnframes/hb->wordlen);
1401     }
1402
1403     /* Copy temp array to target array */
1404     for (m = 0; (m <= nnframes); m++)
1405     {
1406         _set_hb(hb0->h[0], m, htmp[m]);
1407         _set_hb(hb0->g[0], m, gtmp[m]);
1408     }
1409
1410     /* Set scalar variables */
1411     hb0->n0        = nn0;
1412     hb0->maxframes = nnframes;
1413 }
1414
1415 static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact)
1416 {
1417     int           i, inrnew, indnew, j, ii, jj, m, id, ia, grp, ogrp, ntmp;
1418     unsigned int *htmp, *gtmp;
1419     t_hbond      *hb0, *hb1;
1420
1421     inrnew = hb->nrhb;
1422     indnew = hb->nrdist;
1423
1424     /* Check whether donors are also acceptors */
1425     printf("Merging hbonds with Acceptor and Donor swapped\n");
1426
1427     ntmp = 2*hb->max_frames;
1428     snew(gtmp, ntmp);
1429     snew(htmp, ntmp);
1430     for (i = 0; (i < hb->d.nrd); i++)
1431     {
1432         fprintf(stderr, "\r%d/%d", i+1, hb->d.nrd);
1433         id = hb->d.don[i];
1434         ii = hb->a.aptr[id];
1435         for (j = 0; (j < hb->a.nra); j++)
1436         {
1437             ia = hb->a.acc[j];
1438             jj = hb->d.dptr[ia];
1439             if ((id != ia) && (ii != NOTSET) && (jj != NOTSET) &&
1440                 (!bTwo || (bTwo && (hb->d.grp[i] != hb->a.grp[j]))))
1441             {
1442                 hb0 = hb->hbmap[i][j];
1443                 hb1 = hb->hbmap[jj][ii];
1444                 if (hb0 && hb1 && ISHB(hb0->history[0]) && ISHB(hb1->history[0]))
1445                 {
1446                     do_merge(hb, ntmp, htmp, gtmp, hb0, hb1);
1447                     if (ISHB(hb1->history[0]))
1448                     {
1449                         inrnew--;
1450                     }
1451                     else if (ISDIST(hb1->history[0]))
1452                     {
1453                         indnew--;
1454                     }
1455                     else
1456                     if (bContact)
1457                     {
1458                         gmx_incons("No contact history");
1459                     }
1460                     else
1461                     {
1462                         gmx_incons("Neither hydrogen bond nor distance");
1463                     }
1464                     sfree(hb1->h[0]);
1465                     sfree(hb1->g[0]);
1466                     hb1->h[0]       = NULL;
1467                     hb1->g[0]       = NULL;
1468                     hb1->history[0] = hbNo;
1469                 }
1470             }
1471         }
1472     }
1473     fprintf(stderr, "\n");
1474     printf("- Reduced number of hbonds from %d to %d\n", hb->nrhb, inrnew);
1475     printf("- Reduced number of distances from %d to %d\n", hb->nrdist, indnew);
1476     hb->nrhb   = inrnew;
1477     hb->nrdist = indnew;
1478     sfree(gtmp);
1479     sfree(htmp);
1480 }
1481
1482 static void do_nhb_dist(FILE *fp, t_hbdata *hb, real t)
1483 {
1484     int  i, j, k, n_bound[MAXHH], nbtot;
1485     h_id nhb;
1486
1487
1488     /* Set array to 0 */
1489     for (k = 0; (k < MAXHH); k++)
1490     {
1491         n_bound[k] = 0;
1492     }
1493     /* Loop over possible donors */
1494     for (i = 0; (i < hb->d.nrd); i++)
1495     {
1496         for (j = 0; (j < hb->d.nhydro[i]); j++)
1497         {
1498             n_bound[hb->d.nhbonds[i][j]]++;
1499         }
1500     }
1501     fprintf(fp, "%12.5e", t);
1502     nbtot = 0;
1503     for (k = 0; (k < MAXHH); k++)
1504     {
1505         fprintf(fp, "  %8d", n_bound[k]);
1506         nbtot += n_bound[k]*k;
1507     }
1508     fprintf(fp, "  %8d\n", nbtot);
1509 }
1510
1511 static void do_hblife(const char *fn, t_hbdata *hb, gmx_bool bMerge, gmx_bool bContact,
1512                       const output_env_t oenv)
1513 {
1514     FILE          *fp;
1515     const char    *leg[] = { "p(t)", "t p(t)" };
1516     int           *histo;
1517     int            i, j, j0, k, m, nh, ihb, ohb, nhydro, ndump = 0;
1518     int            nframes = hb->nframes;
1519     unsigned int **h;
1520     real           t, x1, dt;
1521     double         sum, integral;
1522     t_hbond       *hbh;
1523
1524     snew(h, hb->maxhydro);
1525     snew(histo, nframes+1);
1526     /* Total number of hbonds analyzed here */
1527     for (i = 0; (i < hb->d.nrd); i++)
1528     {
1529         for (k = 0; (k < hb->a.nra); k++)
1530         {
1531             hbh = hb->hbmap[i][k];
1532             if (hbh)
1533             {
1534                 if (bMerge)
1535                 {
1536                     if (hbh->h[0])
1537                     {
1538                         h[0]   = hbh->h[0];
1539                         nhydro = 1;
1540                     }
1541                     else
1542                     {
1543                         nhydro = 0;
1544                     }
1545                 }
1546                 else
1547                 {
1548                     nhydro = 0;
1549                     for (m = 0; (m < hb->maxhydro); m++)
1550                     {
1551                         if (hbh->h[m])
1552                         {
1553                             h[nhydro++] = bContact ? hbh->g[m] : hbh->h[m];
1554                         }
1555                     }
1556                 }
1557                 for (nh = 0; (nh < nhydro); nh++)
1558                 {
1559                     ohb = 0;
1560                     j0  = 0;
1561
1562                     for (j = 0; (j <= hbh->nframes); j++)
1563                     {
1564                         ihb      = is_hb(h[nh], j);
1565                         if (debug && (ndump < 10))
1566                         {
1567                             fprintf(debug, "%5d  %5d\n", j, ihb);
1568                         }
1569                         if (ihb != ohb)
1570                         {
1571                             if (ihb)
1572                             {
1573                                 j0 = j;
1574                             }
1575                             else
1576                             {
1577                                 histo[j-j0]++;
1578                             }
1579                             ohb = ihb;
1580                         }
1581                     }
1582                     ndump++;
1583                 }
1584             }
1585         }
1586     }
1587     fprintf(stderr, "\n");
1588     if (bContact)
1589     {
1590         fp = xvgropen(fn, "Uninterrupted contact lifetime", output_env_get_xvgr_tlabel(oenv), "()", oenv);
1591     }
1592     else
1593     {
1594         fp = xvgropen(fn, "Uninterrupted hydrogen bond lifetime", output_env_get_xvgr_tlabel(oenv), "()",
1595                       oenv);
1596     }
1597
1598     xvgr_legend(fp, asize(leg), leg, oenv);
1599     j0 = nframes-1;
1600     while ((j0 > 0) && (histo[j0] == 0))
1601     {
1602         j0--;
1603     }
1604     sum = 0;
1605     for (i = 0; (i <= j0); i++)
1606     {
1607         sum += histo[i];
1608     }
1609     dt       = hb->time[1]-hb->time[0];
1610     sum      = dt*sum;
1611     integral = 0;
1612     for (i = 1; (i <= j0); i++)
1613     {
1614         t  = hb->time[i] - hb->time[0] - 0.5*dt;
1615         x1 = t*histo[i]/sum;
1616         fprintf(fp, "%8.3f  %10.3e  %10.3e\n", t, histo[i]/sum, x1);
1617         integral += x1;
1618     }
1619     integral *= dt;
1620     xvgrclose(fp);
1621     printf("%s lifetime = %.2f ps\n", bContact ? "Contact" : "HB", integral);
1622     printf("Note that the lifetime obtained in this manner is close to useless\n");
1623     printf("Use the -ac option instead and check the Forward lifetime\n");
1624     please_cite(stdout, "Spoel2006b");
1625     sfree(h);
1626     sfree(histo);
1627 }
1628
1629 static void dump_ac(t_hbdata *hb, gmx_bool oneHB, int nDump)
1630 {
1631     FILE     *fp;
1632     int       i, j, k, m, nd, ihb, idist;
1633     int       nframes = hb->nframes;
1634     gmx_bool  bPrint;
1635     t_hbond  *hbh;
1636
1637     if (nDump <= 0)
1638     {
1639         return;
1640     }
1641     fp = gmx_ffopen("debug-ac.xvg", "w");
1642     for (j = 0; (j < nframes); j++)
1643     {
1644         fprintf(fp, "%10.3f", hb->time[j]);
1645         for (i = nd = 0; (i < hb->d.nrd) && (nd < nDump); i++)
1646         {
1647             for (k = 0; (k < hb->a.nra) && (nd < nDump); k++)
1648             {
1649                 bPrint = FALSE;
1650                 ihb    = idist = 0;
1651                 hbh    = hb->hbmap[i][k];
1652                 if (oneHB)
1653                 {
1654                     if (hbh->h[0])
1655                     {
1656                         ihb    = is_hb(hbh->h[0], j);
1657                         idist  = is_hb(hbh->g[0], j);
1658                         bPrint = TRUE;
1659                     }
1660                 }
1661                 else
1662                 {
1663                     for (m = 0; (m < hb->maxhydro) && !ihb; m++)
1664                     {
1665                         ihb   = ihb   || ((hbh->h[m]) && is_hb(hbh->h[m], j));
1666                         idist = idist || ((hbh->g[m]) && is_hb(hbh->g[m], j));
1667                     }
1668                     /* This is not correct! */
1669                     /* What isn't correct? -Erik M */
1670                     bPrint = TRUE;
1671                 }
1672                 if (bPrint)
1673                 {
1674                     fprintf(fp, "  %1d-%1d", ihb, idist);
1675                     nd++;
1676                 }
1677             }
1678         }
1679         fprintf(fp, "\n");
1680     }
1681     gmx_ffclose(fp);
1682 }
1683
1684 static real calc_dg(real tau, real temp)
1685 {
1686     real kbt;
1687
1688     kbt = BOLTZ*temp;
1689     if (tau <= 0)
1690     {
1691         return -666;
1692     }
1693     else
1694     {
1695         return kbt*log(kbt*tau/PLANCK);
1696     }
1697 }
1698
1699 typedef struct {
1700     int   n0, n1, nparams, ndelta;
1701     real  kkk[2];
1702     real *t, *ct, *nt, *kt, *sigma_ct, *sigma_nt, *sigma_kt;
1703 } t_luzar;
1704
1705 static real compute_weighted_rates(int n, real t[], real ct[], real nt[],
1706                                    real kt[], real sigma_ct[], real sigma_nt[],
1707                                    real sigma_kt[], real *k, real *kp,
1708                                    real *sigma_k, real *sigma_kp,
1709                                    real fit_start)
1710 {
1711 #define NK 10
1712     int      i, j;
1713     t_luzar  tl;
1714     real     kkk = 0, kkp = 0, kk2 = 0, kp2 = 0, chi2;
1715
1716     *sigma_k  = 0;
1717     *sigma_kp = 0;
1718
1719     for (i = 0; (i < n); i++)
1720     {
1721         if (t[i] >= fit_start)
1722         {
1723             break;
1724         }
1725     }
1726     tl.n0       = i;
1727     tl.n1       = n;
1728     tl.nparams  = 2;
1729     tl.ndelta   = 1;
1730     tl.t        = t;
1731     tl.ct       = ct;
1732     tl.nt       = nt;
1733     tl.kt       = kt;
1734     tl.sigma_ct = sigma_ct;
1735     tl.sigma_nt = sigma_nt;
1736     tl.sigma_kt = sigma_kt;
1737     tl.kkk[0]   = *k;
1738     tl.kkk[1]   = *kp;
1739
1740     chi2      = 0; /*optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */
1741     *k        = tl.kkk[0];
1742     *kp       = tl.kkk[1] = *kp;
1743     tl.ndelta = NK;
1744     for (j = 0; (j < NK); j++)
1745     {
1746         /* (void) optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */
1747         kkk += tl.kkk[0];
1748         kkp += tl.kkk[1];
1749         kk2 += sqr(tl.kkk[0]);
1750         kp2 += sqr(tl.kkk[1]);
1751         tl.n0++;
1752     }
1753     *sigma_k  = sqrt(kk2/NK - sqr(kkk/NK));
1754     *sigma_kp = sqrt(kp2/NK - sqr(kkp/NK));
1755
1756     return chi2;
1757 }
1758
1759 void analyse_corr(int n, real t[], real ct[], real nt[], real kt[],
1760                   real sigma_ct[], real sigma_nt[], real sigma_kt[],
1761                   real fit_start, real temp)
1762 {
1763     int        i0, i;
1764     real       k = 1, kp = 1, kow = 1;
1765     real       Q = 0, chi22, chi2, dg, dgp, tau_hb, dtau, tau_rlx, e_1, dt, sigma_k, sigma_kp, ddg;
1766     double     tmp, sn2 = 0, sc2 = 0, sk2 = 0, scn = 0, sck = 0, snk = 0;
1767     gmx_bool   bError = (sigma_ct != NULL) && (sigma_nt != NULL) && (sigma_kt != NULL);
1768
1769     for (i0 = 0; (i0 < n-2) && ((t[i0]-t[0]) < fit_start); i0++)
1770     {
1771         ;
1772     }
1773     if (i0 < n-2)
1774     {
1775         for (i = i0; (i < n); i++)
1776         {
1777             sc2 += sqr(ct[i]);
1778             sn2 += sqr(nt[i]);
1779             sk2 += sqr(kt[i]);
1780             sck += ct[i]*kt[i];
1781             snk += nt[i]*kt[i];
1782             scn += ct[i]*nt[i];
1783         }
1784         printf("Hydrogen bond thermodynamics at T = %g K\n", temp);
1785         tmp = (sn2*sc2-sqr(scn));
1786         if ((tmp > 0) && (sn2 > 0))
1787         {
1788             k    = (sn2*sck-scn*snk)/tmp;
1789             kp   = (k*scn-snk)/sn2;
1790             if (bError)
1791             {
1792                 chi2 = 0;
1793                 for (i = i0; (i < n); i++)
1794                 {
1795                     chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
1796                 }
1797                 chi22 = compute_weighted_rates(n, t, ct, nt, kt, sigma_ct, sigma_nt,
1798                                                sigma_kt, &k, &kp,
1799                                                &sigma_k, &sigma_kp, fit_start);
1800                 Q   = 0; /* quality_of_fit(chi2, 2);*/
1801                 ddg = BOLTZ*temp*sigma_k/k;
1802                 printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n",
1803                        chi2, Q);
1804                 printf("The Rate and Delta G are followed by an error estimate\n");
1805                 printf("----------------------------------------------------------\n"
1806                        "Type      Rate (1/ps)  Sigma Time (ps)  DG (kJ/mol)  Sigma\n");
1807                 printf("Forward    %10.3f %6.2f   %8.3f  %10.3f %6.2f\n",
1808                        k, sigma_k, 1/k, calc_dg(1/k, temp), ddg);
1809                 ddg = BOLTZ*temp*sigma_kp/kp;
1810                 printf("Backward   %10.3f %6.2f   %8.3f  %10.3f %6.2f\n",
1811                        kp, sigma_kp, 1/kp, calc_dg(1/kp, temp), ddg);
1812             }
1813             else
1814             {
1815                 chi2 = 0;
1816                 for (i = i0; (i < n); i++)
1817                 {
1818                     chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
1819                 }
1820                 printf("Fitting parameters chi^2 = %10g\nQ = %10g\n",
1821                        chi2, Q);
1822                 printf("--------------------------------------------------\n"
1823                        "Type      Rate (1/ps) Time (ps)  DG (kJ/mol)  Chi^2\n");
1824                 printf("Forward    %10.3f   %8.3f  %10.3f  %10g\n",
1825                        k, 1/k, calc_dg(1/k, temp), chi2);
1826                 printf("Backward   %10.3f   %8.3f  %10.3f\n",
1827                        kp, 1/kp, calc_dg(1/kp, temp));
1828             }
1829         }
1830         if (sc2 > 0)
1831         {
1832             kow  = 2*sck/sc2;
1833             printf("One-way    %10.3f   %s%8.3f  %10.3f\n",
1834                    kow, bError ? "       " : "", 1/kow, calc_dg(1/kow, temp));
1835         }
1836         else
1837         {
1838             printf(" - Numerical problems computing HB thermodynamics:\n"
1839                    "sc2 = %g  sn2 = %g  sk2 = %g sck = %g snk = %g scn = %g\n",
1840                    sc2, sn2, sk2, sck, snk, scn);
1841         }
1842         /* Determine integral of the correlation function */
1843         tau_hb = evaluate_integral(n, t, ct, NULL, (t[n-1]-t[0])/2, &dtau);
1844         printf("Integral   %10.3f   %s%8.3f  %10.3f\n", 1/tau_hb,
1845                bError ? "       " : "", tau_hb, calc_dg(tau_hb, temp));
1846         e_1 = exp(-1);
1847         for (i = 0; (i < n-2); i++)
1848         {
1849             if ((ct[i] > e_1) && (ct[i+1] <= e_1))
1850             {
1851                 break;
1852             }
1853         }
1854         if (i < n-2)
1855         {
1856             /* Determine tau_relax from linear interpolation */
1857             tau_rlx = t[i]-t[0] + (e_1-ct[i])*(t[i+1]-t[i])/(ct[i+1]-ct[i]);
1858             printf("Relaxation %10.3f   %8.3f  %s%10.3f\n", 1/tau_rlx,
1859                    tau_rlx, bError ? "       " : "",
1860                    calc_dg(tau_rlx, temp));
1861         }
1862     }
1863     else
1864     {
1865         printf("Correlation functions too short to compute thermodynamics\n");
1866     }
1867 }
1868
1869 void compute_derivative(int nn, real x[], real y[], real dydx[])
1870 {
1871     int j;
1872
1873     /* Compute k(t) = dc(t)/dt */
1874     for (j = 1; (j < nn-1); j++)
1875     {
1876         dydx[j] = (y[j+1]-y[j-1])/(x[j+1]-x[j-1]);
1877     }
1878     /* Extrapolate endpoints */
1879     dydx[0]    = 2*dydx[1]   -  dydx[2];
1880     dydx[nn-1] = 2*dydx[nn-2] - dydx[nn-3];
1881 }
1882
1883 static void parallel_print(int *data, int nThreads)
1884 {
1885     /* This prints the donors on which each tread is currently working. */
1886     int i;
1887
1888     fprintf(stderr, "\r");
1889     for (i = 0; i < nThreads; i++)
1890     {
1891         fprintf(stderr, "%-7i", data[i]);
1892     }
1893 }
1894
1895 static void normalizeACF(real *ct, real *gt, int nhb, int len)
1896 {
1897     real ct_fac, gt_fac = 0;
1898     int  i;
1899
1900     /* Xu and Berne use the same normalization constant */
1901
1902     ct_fac = 1.0/ct[0];
1903     if (nhb != 0)
1904     {
1905         gt_fac = 1.0/(real)nhb;
1906     }
1907
1908     printf("Normalization for c(t) = %g for gh(t) = %g\n", ct_fac, gt_fac);
1909     for (i = 0; i < len; i++)
1910     {
1911         ct[i] *= ct_fac;
1912         if (gt != NULL)
1913         {
1914             gt[i] *= gt_fac;
1915         }
1916     }
1917 }
1918
1919 static void do_hbac(const char *fn, t_hbdata *hb,
1920                     int nDump, gmx_bool bMerge, gmx_bool bContact, real fit_start,
1921                     real temp, gmx_bool R2, const output_env_t oenv,
1922                     int nThreads)
1923 {
1924     FILE          *fp;
1925     int            i, j, k, m, n, o, nd, ihb, idist, n2, nn, iter, nSets;
1926     const char    *legNN[]   = {
1927         "Ac(t)",
1928         "Ac'(t)"
1929     };
1930     static char  **legGem;
1931
1932     const char    *legLuzar[] = {
1933         "Ac\\sfin sys\\v{}\\z{}(t)",
1934         "Ac(t)",
1935         "Cc\\scontact,hb\\v{}\\z{}(t)",
1936         "-dAc\\sfs\\v{}\\z{}/dt"
1937     };
1938     gmx_bool       bNorm = FALSE, bOMP = FALSE;
1939     double         nhb   = 0;
1940     int            nhbi  = 0;
1941     real          *rhbex = NULL, *ht, *gt, *ght, *dght, *kt;
1942     real          *ct, *p_ct, tail, tail2, dtail, ct_fac, ght_fac, *cct;
1943     const real     tol     = 1e-3;
1944     int            nframes = hb->nframes, nf;
1945     unsigned int **h       = NULL, **g = NULL;
1946     int            nh, nhbonds, nhydro, ngh;
1947     t_hbond       *hbh;
1948     int           *ptimes   = NULL, *poff = NULL, anhb, n0, mMax = INT_MIN;
1949     gmx_bool       c;
1950     int            acType;
1951     double        *ctdouble, *timedouble, *fittedct;
1952     double         fittolerance = 0.1;
1953     int           *dondata      = NULL, thisThread;
1954
1955     enum {
1956         AC_NONE, AC_NN, AC_GEM, AC_LUZAR
1957     };
1958
1959 #ifdef GMX_OPENMP
1960     bOMP = TRUE;
1961 #else
1962     bOMP = FALSE;
1963 #endif
1964
1965     printf("Doing autocorrelation ");
1966
1967     acType = AC_LUZAR;
1968     printf("according to the theory of Luzar and Chandler.\n");
1969     fflush(stdout);
1970
1971     /* build hbexist matrix in reals for autocorr */
1972     /* Allocate memory for computing ACF (rhbex) and aggregating the ACF (ct) */
1973     n2 = 1;
1974     while (n2 < nframes)
1975     {
1976         n2 *= 2;
1977     }
1978
1979     nn = nframes/2;
1980
1981     if (acType != AC_NN || bOMP)
1982     {
1983         snew(h, hb->maxhydro);
1984         snew(g, hb->maxhydro);
1985     }
1986
1987     /* Dump hbonds for debugging */
1988     dump_ac(hb, bMerge || bContact, nDump);
1989
1990     /* Total number of hbonds analyzed here */
1991     nhbonds = 0;
1992     ngh     = 0;
1993     anhb    = 0;
1994
1995     if (acType != AC_LUZAR && bOMP)
1996     {
1997         nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
1998
1999         gmx_omp_set_num_threads(nThreads);
2000         snew(dondata, nThreads);
2001         for (i = 0; i < nThreads; i++)
2002         {
2003             dondata[i] = -1;
2004         }
2005         printf("ACF calculations parallelized with OpenMP using %i threads.\n"
2006                "Expect close to linear scaling over this donor-loop.\n", nThreads);
2007         fflush(stdout);
2008         fprintf(stderr, "Donors: [thread no]\n");
2009         {
2010             char tmpstr[7];
2011             for (i = 0; i < nThreads; i++)
2012             {
2013                 snprintf(tmpstr, 7, "[%i]", i);
2014                 fprintf(stderr, "%-7s", tmpstr);
2015             }
2016         }
2017         fprintf(stderr, "\n");
2018     }
2019
2020
2021     /* Build the ACF */
2022     snew(rhbex, 2*n2);
2023     snew(ct, 2*n2);
2024     snew(gt, 2*n2);
2025     snew(ht, 2*n2);
2026     snew(ght, 2*n2);
2027     snew(dght, 2*n2);
2028
2029     snew(kt, nn);
2030     snew(cct, nn);
2031
2032     for (i = 0; (i < hb->d.nrd); i++)
2033     {
2034         for (k = 0; (k < hb->a.nra); k++)
2035         {
2036             nhydro = 0;
2037             hbh    = hb->hbmap[i][k];
2038
2039             if (hbh)
2040             {
2041                 if (bMerge || bContact)
2042                 {
2043                     if (ISHB(hbh->history[0]))
2044                     {
2045                         h[0]   = hbh->h[0];
2046                         g[0]   = hbh->g[0];
2047                         nhydro = 1;
2048                     }
2049                 }
2050                 else
2051                 {
2052                     for (m = 0; (m < hb->maxhydro); m++)
2053                     {
2054                         if (bContact ? ISDIST(hbh->history[m]) : ISHB(hbh->history[m]))
2055                         {
2056                             g[nhydro] = hbh->g[m];
2057                             h[nhydro] = hbh->h[m];
2058                             nhydro++;
2059                         }
2060                     }
2061                 }
2062
2063                 nf = hbh->nframes;
2064                 for (nh = 0; (nh < nhydro); nh++)
2065                 {
2066                     int nrint = bContact ? hb->nrdist : hb->nrhb;
2067                     if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint))
2068                     {
2069                         fprintf(stderr, "\rACF %d/%d", nhbonds+1, nrint);
2070                     }
2071                     nhbonds++;
2072                     for (j = 0; (j < nframes); j++)
2073                     {
2074                         if (j <= nf)
2075                         {
2076                             ihb   = is_hb(h[nh], j);
2077                             idist = is_hb(g[nh], j);
2078                         }
2079                         else
2080                         {
2081                             ihb = idist = 0;
2082                         }
2083                         rhbex[j] = ihb;
2084                         /* For contacts: if a second cut-off is provided, use it,
2085                          * otherwise use g(t) = 1-h(t) */
2086                         if (!R2 && bContact)
2087                         {
2088                             gt[j]  = 1-ihb;
2089                         }
2090                         else
2091                         {
2092                             gt[j]  = idist*(1-ihb);
2093                         }
2094                         ht[j]    = rhbex[j];
2095                         nhb     += ihb;
2096                     }
2097
2098                     /* The autocorrelation function is normalized after summation only */
2099                     low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &rhbex, hb->time[1]-hb->time[0],
2100                                     eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0);
2101
2102                     /* Cross correlation analysis for thermodynamics */
2103                     for (j = nframes; (j < n2); j++)
2104                     {
2105                         ht[j] = 0;
2106                         gt[j] = 0;
2107                     }
2108
2109                     cross_corr(n2, ht, gt, dght);
2110
2111                     for (j = 0; (j < nn); j++)
2112                     {
2113                         ct[j]  += rhbex[j];
2114                         ght[j] += dght[j];
2115                     }
2116                 }
2117             }
2118         }
2119     }
2120     fprintf(stderr, "\n");
2121     sfree(h);
2122     sfree(g);
2123     normalizeACF(ct, ght, nhb, nn);
2124
2125     /* Determine tail value for statistics */
2126     tail  = 0;
2127     tail2 = 0;
2128     for (j = nn/2; (j < nn); j++)
2129     {
2130         tail  += ct[j];
2131         tail2 += ct[j]*ct[j];
2132     }
2133     tail  /= (nn - nn/2);
2134     tail2 /= (nn - nn/2);
2135     dtail  = sqrt(tail2-tail*tail);
2136
2137     /* Check whether the ACF is long enough */
2138     if (dtail > tol)
2139     {
2140         printf("\nWARNING: Correlation function is probably not long enough\n"
2141                "because the standard deviation in the tail of C(t) > %g\n"
2142                "Tail value (average C(t) over second half of acf): %g +/- %g\n",
2143                tol, tail, dtail);
2144     }
2145     for (j = 0; (j < nn); j++)
2146     {
2147         cct[j] = ct[j];
2148         ct[j]  = (cct[j]-tail)/(1-tail);
2149     }
2150     /* Compute negative derivative k(t) = -dc(t)/dt */
2151     compute_derivative(nn, hb->time, ct, kt);
2152     for (j = 0; (j < nn); j++)
2153     {
2154         kt[j] = -kt[j];
2155     }
2156
2157
2158     if (bContact)
2159     {
2160         fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
2161     }
2162     else
2163     {
2164         fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
2165     }
2166     xvgr_legend(fp, asize(legLuzar), legLuzar, oenv);
2167
2168
2169     for (j = 0; (j < nn); j++)
2170     {
2171         fprintf(fp, "%10g  %10g  %10g  %10g  %10g\n",
2172                 hb->time[j]-hb->time[0], ct[j], cct[j], ght[j], kt[j]);
2173     }
2174     xvgrclose(fp);
2175
2176     analyse_corr(nn, hb->time, ct, ght, kt, NULL, NULL, NULL,
2177                  fit_start, temp);
2178
2179     do_view(oenv, fn, NULL);
2180     sfree(rhbex);
2181     sfree(ct);
2182     sfree(gt);
2183     sfree(ht);
2184     sfree(ght);
2185     sfree(dght);
2186     sfree(cct);
2187     sfree(kt);
2188 }
2189
2190 static void init_hbframe(t_hbdata *hb, int nframes, real t)
2191 {
2192     int i, j, m;
2193
2194     hb->time[nframes]   = t;
2195     hb->nhb[nframes]    = 0;
2196     hb->ndist[nframes]  = 0;
2197     for (i = 0; (i < max_hx); i++)
2198     {
2199         hb->nhx[nframes][i] = 0;
2200     }
2201 }
2202
2203 static FILE *open_donor_properties_file(const char        *fn,
2204                                         t_hbdata          *hb,
2205                                         const output_env_t oenv)
2206 {
2207     FILE       *fp    = NULL;
2208     const char *leg[] = { "Nbound", "Nfree" };
2209
2210     if (!fn || !hb)
2211     {
2212         return NULL;
2213     }
2214
2215     fp = xvgropen(fn, "Donor properties", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
2216     xvgr_legend(fp, asize(leg), leg, oenv);
2217
2218     return fp;
2219 }
2220
2221 static void analyse_donor_properties(FILE *fp, t_hbdata *hb, int nframes, real t)
2222 {
2223     int i, j, k, nbound, nb, nhtot;
2224
2225     if (!fp || !hb)
2226     {
2227         return;
2228     }
2229     nbound = 0;
2230     nhtot  = 0;
2231     for (i = 0; (i < hb->d.nrd); i++)
2232     {
2233         for (k = 0; (k < hb->d.nhydro[i]); k++)
2234         {
2235             nb = 0;
2236             nhtot++;
2237             for (j = 0; (j < hb->a.nra) && (nb == 0); j++)
2238             {
2239                 if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[k] &&
2240                     is_hb(hb->hbmap[i][j]->h[k], nframes))
2241                 {
2242                     nb = 1;
2243                 }
2244             }
2245             nbound += nb;
2246         }
2247     }
2248     fprintf(fp, "%10.3e  %6d  %6d\n", t, nbound, nhtot-nbound);
2249 }
2250
2251 static void dump_hbmap(t_hbdata *hb,
2252                        int nfile, t_filenm fnm[], gmx_bool bTwo,
2253                        gmx_bool bContact, int isize[], int *index[], char *grpnames[],
2254                        t_atoms *atoms)
2255 {
2256     FILE    *fp, *fplog;
2257     int      ddd, hhh, aaa, i, j, k, m, grp;
2258     char     ds[32], hs[32], as[32];
2259     gmx_bool first;
2260
2261     fp = opt2FILE("-hbn", nfile, fnm, "w");
2262     if (opt2bSet("-g", nfile, fnm))
2263     {
2264         fplog = gmx_ffopen(opt2fn("-g", nfile, fnm), "w");
2265         fprintf(fplog, "# %10s  %12s  %12s\n", "Donor", "Hydrogen", "Acceptor");
2266     }
2267     else
2268     {
2269         fplog = NULL;
2270     }
2271     for (grp = gr0; grp <= (bTwo ? gr1 : gr0); grp++)
2272     {
2273         fprintf(fp, "[ %s ]", grpnames[grp]);
2274         for (i = 0; i < isize[grp]; i++)
2275         {
2276             fprintf(fp, (i%15) ? " " : "\n");
2277             fprintf(fp, " %4d", index[grp][i]+1);
2278         }
2279         fprintf(fp, "\n");
2280
2281         if (!bContact)
2282         {
2283             fprintf(fp, "[ donors_hydrogens_%s ]\n", grpnames[grp]);
2284             for (i = 0; (i < hb->d.nrd); i++)
2285             {
2286                 if (hb->d.grp[i] == grp)
2287                 {
2288                     for (j = 0; (j < hb->d.nhydro[i]); j++)
2289                     {
2290                         fprintf(fp, " %4d %4d", hb->d.don[i]+1,
2291                                 hb->d.hydro[i][j]+1);
2292                     }
2293                     fprintf(fp, "\n");
2294                 }
2295             }
2296             first = TRUE;
2297             fprintf(fp, "[ acceptors_%s ]", grpnames[grp]);
2298             for (i = 0; (i < hb->a.nra); i++)
2299             {
2300                 if (hb->a.grp[i] == grp)
2301                 {
2302                     fprintf(fp, (i%15 && !first) ? " " : "\n");
2303                     fprintf(fp, " %4d", hb->a.acc[i]+1);
2304                     first = FALSE;
2305                 }
2306             }
2307             fprintf(fp, "\n");
2308         }
2309     }
2310     if (bTwo)
2311     {
2312         fprintf(fp, bContact ? "[ contacts_%s-%s ]\n" :
2313                 "[ hbonds_%s-%s ]\n", grpnames[0], grpnames[1]);
2314     }
2315     else
2316     {
2317         fprintf(fp, bContact ? "[ contacts_%s ]" : "[ hbonds_%s ]\n", grpnames[0]);
2318     }
2319
2320     for (i = 0; (i < hb->d.nrd); i++)
2321     {
2322         ddd = hb->d.don[i];
2323         for (k = 0; (k < hb->a.nra); k++)
2324         {
2325             aaa = hb->a.acc[k];
2326             for (m = 0; (m < hb->d.nhydro[i]); m++)
2327             {
2328                 if (hb->hbmap[i][k] && ISHB(hb->hbmap[i][k]->history[m]))
2329                 {
2330                     sprintf(ds, "%s", mkatomname(atoms, ddd));
2331                     sprintf(as, "%s", mkatomname(atoms, aaa));
2332                     if (bContact)
2333                     {
2334                         fprintf(fp, " %6d %6d\n", ddd+1, aaa+1);
2335                         if (fplog)
2336                         {
2337                             fprintf(fplog, "%12s  %12s\n", ds, as);
2338                         }
2339                     }
2340                     else
2341                     {
2342                         hhh = hb->d.hydro[i][m];
2343                         sprintf(hs, "%s", mkatomname(atoms, hhh));
2344                         fprintf(fp, " %6d %6d %6d\n", ddd+1, hhh+1, aaa+1);
2345                         if (fplog)
2346                         {
2347                             fprintf(fplog, "%12s  %12s  %12s\n", ds, hs, as);
2348                         }
2349                     }
2350                 }
2351             }
2352         }
2353     }
2354     gmx_ffclose(fp);
2355     if (fplog)
2356     {
2357         gmx_ffclose(fplog);
2358     }
2359 }
2360
2361 /* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template.
2362  * It mimics add_frames() and init_frame() to some extent. */
2363 static void sync_hbdata(t_hbdata *p_hb, int nframes)
2364 {
2365     int i;
2366     if (nframes >= p_hb->max_frames)
2367     {
2368         p_hb->max_frames += 4096;
2369         srenew(p_hb->nhb,   p_hb->max_frames);
2370         srenew(p_hb->ndist, p_hb->max_frames);
2371         srenew(p_hb->n_bound, p_hb->max_frames);
2372         srenew(p_hb->nhx, p_hb->max_frames);
2373         if (p_hb->bDAnr)
2374         {
2375             srenew(p_hb->danr, p_hb->max_frames);
2376         }
2377         memset(&(p_hb->nhb[nframes]),   0, sizeof(int) * (p_hb->max_frames-nframes));
2378         memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
2379         p_hb->nhb[nframes]   = 0;
2380         p_hb->ndist[nframes] = 0;
2381
2382     }
2383     p_hb->nframes = nframes;
2384 /*     for (i=0;) */
2385 /*     { */
2386 /*         p_hb->nhx[nframes][i] */
2387 /*     } */
2388     memset(&(p_hb->nhx[nframes]), 0, sizeof(int)*max_hx); /* zero the helix count for this frame */
2389
2390     /* hb->per will remain constant througout the frame loop,
2391      * even though the data its members point to will change,
2392      * hence no need for re-syncing. */
2393 }
2394
2395 int gmx_hbond(int argc, char *argv[])
2396 {
2397     const char        *desc[] = {
2398         "[THISMODULE] computes and analyzes hydrogen bonds. Hydrogen bonds are",
2399         "determined based on cutoffs for the angle Hydrogen - Donor - Acceptor",
2400         "(zero is extended) and the distance Donor - Acceptor",
2401         "(or Hydrogen - Acceptor using [TT]-noda[tt]).",
2402         "OH and NH groups are regarded as donors, O is an acceptor always,",
2403         "N is an acceptor by default, but this can be switched using",
2404         "[TT]-nitacc[tt]. Dummy hydrogen atoms are assumed to be connected",
2405         "to the first preceding non-hydrogen atom.[PAR]",
2406
2407         "You need to specify two groups for analysis, which must be either",
2408         "identical or non-overlapping. All hydrogen bonds between the two",
2409         "groups are analyzed.[PAR]",
2410
2411         "If you set [TT]-shell[tt], you will be asked for an additional index group",
2412         "which should contain exactly one atom. In this case, only hydrogen",
2413         "bonds between atoms within the shell distance from the one atom are",
2414         "considered.[PAR]",
2415
2416         "With option -ac, rate constants for hydrogen bonding can be derived with the model of Luzar and Chandler",
2417         "(Nature 394, 1996; J. Chem. Phys. 113:23, 2000) or that of Markovitz and Agmon (J. Chem. Phys 129, 2008).",
2418         "If contact kinetics are analyzed by using the -contact option, then",
2419         "n(t) can be defined as either all pairs that are not within contact distance r at time t",
2420         "(corresponding to leaving the -r2 option at the default value 0) or all pairs that",
2421         "are within distance r2 (corresponding to setting a second cut-off value with option -r2).",
2422         "See mentioned literature for more details and definitions."
2423         "[PAR]",
2424
2425         /*    "It is also possible to analyse specific hydrogen bonds with",
2426               "[TT]-sel[tt]. This index file must contain a group of atom triplets",
2427               "Donor Hydrogen Acceptor, in the following way::",
2428            "",
2429            "[ selected ]",
2430            "     20    21    24",
2431            "     25    26    29",
2432            "      1     3     6",
2433            "",
2434            "Note that the triplets need not be on separate lines.",
2435            "Each atom triplet specifies a hydrogen bond to be analyzed,",
2436            "note also that no check is made for the types of atoms.[PAR]",
2437          */
2438
2439         "[BB]Output:[bb]",
2440         "",
2441         " * [TT]-num[tt]:  number of hydrogen bonds as a function of time.",
2442         " * [TT]-ac[tt]:   average over all autocorrelations of the existence",
2443         "   functions (either 0 or 1) of all hydrogen bonds.",
2444         " * [TT]-dist[tt]: distance distribution of all hydrogen bonds.",
2445         " * [TT]-ang[tt]:  angle distribution of all hydrogen bonds.",
2446         " * [TT]-hx[tt]:   the number of n-n+i hydrogen bonds as a function of time",
2447         "   where n and n+i stand for residue numbers and i ranges from 0 to 6.",
2448         "   This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated",
2449         "   with helices in proteins.",
2450         " * [TT]-hbn[tt]:  all selected groups, donors, hydrogens and acceptors",
2451         "   for selected groups, all hydrogen bonded atoms from all groups and",
2452         "   all solvent atoms involved in insertion.",
2453         " * [TT]-hbm[tt]:  existence matrix for all hydrogen bonds over all",
2454         "   frames, this also contains information on solvent insertion",
2455         "   into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]",
2456         "   index file.",
2457         " * [TT]-dan[tt]: write out the number of donors and acceptors analyzed for",
2458         "   each timeframe. This is especially useful when using [TT]-shell[tt].",
2459         " * [TT]-nhbdist[tt]: compute the number of HBonds per hydrogen in order to",
2460         "   compare results to Raman Spectroscopy.",
2461         "",
2462         "Note: options [TT]-ac[tt], [TT]-life[tt], [TT]-hbn[tt] and [TT]-hbm[tt]",
2463         "require an amount of memory proportional to the total numbers of donors",
2464         "times the total number of acceptors in the selected group(s)."
2465     };
2466
2467     static real        acut     = 30, abin = 1, rcut = 0.35, r2cut = 0, rbin = 0.005, rshell = -1;
2468     static real        maxnhb   = 0, fit_start = 1, fit_end = 60, temp = 298.15, D = -1;
2469     static gmx_bool    bNitAcc  = TRUE, bDA = TRUE, bMerge = TRUE;
2470     static int         nDump    = 0, nFitPoints = 100;
2471     static int         nThreads = 0, nBalExp = 4;
2472
2473     static gmx_bool    bContact     = FALSE;
2474     static real        logAfterTime = 10, gemBallistic = 0.2; /* ps */
2475     static const char *NNtype[]     = {NULL, "none", "binary", "oneOverR3", "dipole", NULL};
2476
2477     /* options */
2478     t_pargs     pa [] = {
2479         { "-a",    FALSE,  etREAL, {&acut},
2480           "Cutoff angle (degrees, Hydrogen - Donor - Acceptor)" },
2481         { "-r",    FALSE,  etREAL, {&rcut},
2482           "Cutoff radius (nm, X - Acceptor, see next option)" },
2483         { "-da",   FALSE,  etBOOL, {&bDA},
2484           "Use distance Donor-Acceptor (if TRUE) or Hydrogen-Acceptor (FALSE)" },
2485         { "-r2",   FALSE,  etREAL, {&r2cut},
2486           "Second cutoff radius. Mainly useful with [TT]-contact[tt] and [TT]-ac[tt]"},
2487         { "-abin", FALSE,  etREAL, {&abin},
2488           "Binwidth angle distribution (degrees)" },
2489         { "-rbin", FALSE,  etREAL, {&rbin},
2490           "Binwidth distance distribution (nm)" },
2491         { "-nitacc", FALSE, etBOOL, {&bNitAcc},
2492           "Regard nitrogen atoms as acceptors" },
2493         { "-contact", FALSE, etBOOL, {&bContact},
2494           "Do not look for hydrogen bonds, but merely for contacts within the cut-off distance" },
2495         { "-shell", FALSE, etREAL, {&rshell},
2496           "when > 0, only calculate hydrogen bonds within # nm shell around "
2497           "one particle" },
2498         { "-fitstart", FALSE, etREAL, {&fit_start},
2499           "Time (ps) from which to start fitting the correlation functions in order to obtain the forward and backward rate constants for HB breaking and formation. With [TT]-gemfit[tt] we suggest [TT]-fitstart 0[tt]" },
2500         { "-fitend", FALSE, etREAL, {&fit_end},
2501           "Time (ps) to which to stop fitting the correlation functions in order to obtain the forward and backward rate constants for HB breaking and formation (only with [TT]-gemfit[tt])" },
2502         { "-temp",  FALSE, etREAL, {&temp},
2503           "Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" },
2504         { "-dump",  FALSE, etINT, {&nDump},
2505           "Dump the first N hydrogen bond ACFs in a single [REF].xvg[ref] file for debugging" },
2506         { "-max_hb", FALSE, etREAL, {&maxnhb},
2507           "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" },
2508         { "-merge", FALSE, etBOOL, {&bMerge},
2509           "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." },
2510 #ifdef GMX_OPENMP
2511         { "-nthreads", FALSE, etINT, {&nThreads},
2512           "Number of threads used for the parallel loop over autocorrelations. nThreads <= 0 means maximum number of threads. Requires linking with OpenMP. The number of threads is limited by the number of cores (before OpenMP v.3 ) or environment variable OMP_THREAD_LIMIT (OpenMP v.3)"},
2513 #endif
2514     };
2515     const char *bugs[] = {
2516         "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and therefore not available for the time being."
2517     };
2518     t_filenm    fnm[] = {
2519         { efTRX, "-f",   NULL,     ffREAD  },
2520         { efTPR, NULL,   NULL,     ffREAD  },
2521         { efNDX, NULL,   NULL,     ffOPTRD },
2522         /*    { efNDX, "-sel", "select", ffOPTRD },*/
2523         { efXVG, "-num", "hbnum",  ffWRITE },
2524         { efLOG, "-g",   "hbond",  ffOPTWR },
2525         { efXVG, "-ac",  "hbac",   ffOPTWR },
2526         { efXVG, "-dist", "hbdist", ffOPTWR },
2527         { efXVG, "-ang", "hbang",  ffOPTWR },
2528         { efXVG, "-hx",  "hbhelix", ffOPTWR },
2529         { efNDX, "-hbn", "hbond",  ffOPTWR },
2530         { efXPM, "-hbm", "hbmap",  ffOPTWR },
2531         { efXVG, "-don", "donor",  ffOPTWR },
2532         { efXVG, "-dan", "danum",  ffOPTWR },
2533         { efXVG, "-life", "hblife", ffOPTWR },
2534         { efXVG, "-nhbdist", "nhbdist", ffOPTWR }
2535
2536     };
2537 #define NFILE asize(fnm)
2538
2539     char                  hbmap [HB_NR] = { ' ',    'o',      '-',       '*' };
2540     const char           *hbdesc[HB_NR] = { "None", "Present", "Inserted", "Present & Inserted" };
2541     t_rgb                 hbrgb [HB_NR] = { {1, 1, 1}, {1, 0, 0},   {0, 0, 1},    {1, 0, 1} };
2542
2543     t_trxstatus          *status;
2544     int                   trrStatus = 1;
2545     t_topology            top;
2546     t_inputrec            ir;
2547     t_pargs              *ppa;
2548     int                   npargs, natoms, nframes = 0, shatom;
2549     int                  *isize;
2550     char                **grpnames;
2551     atom_id             **index;
2552     rvec                 *x, hbox;
2553     matrix                box;
2554     real                  t, ccut, dist = 0.0, ang = 0.0;
2555     double                max_nhb, aver_nhb, aver_dist;
2556     int                   h = 0, i = 0, j, k = 0, l, start, end, id, ja, ogrp, nsel;
2557     int                   xi, yi, zi, ai;
2558     int                   xj, yj, zj, aj, xjj, yjj, zjj;
2559     int                   xk, yk, zk, ak, xkk, ykk, zkk;
2560     gmx_bool              bSelected, bHBmap, bStop, bTwo, was, bBox, bTric;
2561     int                  *adist, *rdist, *aptr, *rprt;
2562     int                   grp, nabin, nrbin, bin, resdist, ihb;
2563     char                **leg;
2564     t_hbdata             *hb, *hbptr;
2565     FILE                 *fp, *fpins = NULL, *fpnhb = NULL, *donor_properties = NULL;
2566     t_gridcell         ***grid;
2567     t_ncell              *icell, *jcell, *kcell;
2568     ivec                  ngrid;
2569     unsigned char        *datable;
2570     output_env_t          oenv;
2571     int                   NN;
2572     int                   ii, jj, hh, actual_nThreads;
2573     int                   threadNr = 0;
2574     gmx_bool              bParallel;
2575     gmx_bool              bEdge_yjj, bEdge_xjj, bOMP;
2576
2577     t_hbdata            **p_hb    = NULL;                   /* one per thread, then merge after the frame loop */
2578     int                 **p_adist = NULL, **p_rdist = NULL; /* a histogram for each thread. */
2579
2580 #ifdef GMX_OPENMP
2581     bOMP = TRUE;
2582 #else
2583     bOMP = FALSE;
2584 #endif
2585
2586     npargs = asize(pa);
2587     ppa    = add_acf_pargs(&npargs, pa);
2588
2589     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT, NFILE, fnm, npargs,
2590                            ppa, asize(desc), desc, asize(bugs), bugs, &oenv))
2591     {
2592         return 0;
2593     }
2594
2595     /* process input */
2596     bSelected = FALSE;
2597     ccut      = cos(acut*DEG2RAD);
2598
2599     if (bContact)
2600     {
2601         if (bSelected)
2602         {
2603             gmx_fatal(FARGS, "Can not analyze selected contacts.");
2604         }
2605         if (!bDA)
2606         {
2607             gmx_fatal(FARGS, "Can not analyze contact between H and A: turn off -noda");
2608         }
2609     }
2610
2611     /* Initiate main data structure! */
2612     bHBmap = (opt2bSet("-ac", NFILE, fnm) ||
2613               opt2bSet("-life", NFILE, fnm) ||
2614               opt2bSet("-hbn", NFILE, fnm) ||
2615               opt2bSet("-hbm", NFILE, fnm));
2616
2617     if (opt2bSet("-nhbdist", NFILE, fnm))
2618     {
2619         const char *leg[MAXHH+1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" };
2620         fpnhb = xvgropen(opt2fn("-nhbdist", NFILE, fnm),
2621                          "Number of donor-H with N HBs", output_env_get_xvgr_tlabel(oenv), "N", oenv);
2622         xvgr_legend(fpnhb, asize(leg), leg, oenv);
2623     }
2624
2625     hb = mk_hbdata(bHBmap, opt2bSet("-dan", NFILE, fnm), bMerge || bContact);
2626
2627     /* get topology */
2628     read_tpx_top(ftp2fn(efTPR, NFILE, fnm), &ir, box, &natoms, NULL, NULL, NULL, &top);
2629
2630     snew(grpnames, grNR);
2631     snew(index, grNR);
2632     snew(isize, grNR);
2633     /* Make Donor-Acceptor table */
2634     snew(datable, top.atoms.nr);
2635     gen_datable(index[0], isize[0], datable, top.atoms.nr);
2636
2637     if (bSelected)
2638     {
2639         /* analyze selected hydrogen bonds */
2640         printf("Select group with selected atoms:\n");
2641         get_index(&(top.atoms), opt2fn("-sel", NFILE, fnm),
2642                   1, &nsel, index, grpnames);
2643         if (nsel % 3)
2644         {
2645             gmx_fatal(FARGS, "Number of atoms in group '%s' not a multiple of 3\n"
2646                       "and therefore cannot contain triplets of "
2647                       "Donor-Hydrogen-Acceptor", grpnames[0]);
2648         }
2649         bTwo = FALSE;
2650
2651         for (i = 0; (i < nsel); i += 3)
2652         {
2653             int dd = index[0][i];
2654             int aa = index[0][i+2];
2655             /* int */ hh = index[0][i+1];
2656             add_dh (&hb->d, dd, hh, i, datable);
2657             add_acc(&hb->a, aa, i);
2658             /* Should this be here ? */
2659             snew(hb->d.dptr, top.atoms.nr);
2660             snew(hb->a.aptr, top.atoms.nr);
2661             add_hbond(hb, dd, aa, hh, gr0, gr0, 0, bMerge, 0, bContact);
2662         }
2663         printf("Analyzing %d selected hydrogen bonds from '%s'\n",
2664                isize[0], grpnames[0]);
2665     }
2666     else
2667     {
2668         /* analyze all hydrogen bonds: get group(s) */
2669         printf("Specify 2 groups to analyze:\n");
2670         get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
2671                   2, isize, index, grpnames);
2672
2673         /* check if we have two identical or two non-overlapping groups */
2674         bTwo = isize[0] != isize[1];
2675         for (i = 0; (i < isize[0]) && !bTwo; i++)
2676         {
2677             bTwo = index[0][i] != index[1][i];
2678         }
2679         if (bTwo)
2680         {
2681             printf("Checking for overlap in atoms between %s and %s\n",
2682                    grpnames[0], grpnames[1]);
2683             for (i = 0; i < isize[1]; i++)
2684             {
2685                 if (ISINGRP(datable[index[1][i]]))
2686                 {
2687                     gmx_fatal(FARGS, "Partial overlap between groups '%s' and '%s'",
2688                               grpnames[0], grpnames[1]);
2689                 }
2690             }
2691             /*
2692                printf("Checking for overlap in atoms between %s and %s\n",
2693                grpnames[0],grpnames[1]);
2694                for (i=0; i<isize[0]; i++)
2695                for (j=0; j<isize[1]; j++)
2696                if (index[0][i] == index[1][j])
2697                gmx_fatal(FARGS,"Partial overlap between groups '%s' and '%s'",
2698                grpnames[0],grpnames[1]);
2699              */
2700         }
2701         if (bTwo)
2702         {
2703             printf("Calculating %s "
2704                    "between %s (%d atoms) and %s (%d atoms)\n",
2705                    bContact ? "contacts" : "hydrogen bonds",
2706                    grpnames[0], isize[0], grpnames[1], isize[1]);
2707         }
2708         else
2709         {
2710             fprintf(stderr, "Calculating %s in %s (%d atoms)\n",
2711                     bContact ? "contacts" : "hydrogen bonds", grpnames[0], isize[0]);
2712         }
2713     }
2714     sfree(datable);
2715
2716     /* search donors and acceptors in groups */
2717     snew(datable, top.atoms.nr);
2718     for (i = 0; (i < grNR); i++)
2719     {
2720         if ( ((i == gr0) && !bSelected ) ||
2721              ((i == gr1) && bTwo ))
2722         {
2723             gen_datable(index[i], isize[i], datable, top.atoms.nr);
2724             if (bContact)
2725             {
2726                 search_acceptors(&top, isize[i], index[i], &hb->a, i,
2727                                  bNitAcc, TRUE, (bTwo && (i == gr0)) || !bTwo, datable);
2728                 search_donors   (&top, isize[i], index[i], &hb->d, i,
2729                                  TRUE, (bTwo && (i == gr1)) || !bTwo, datable);
2730             }
2731             else
2732             {
2733                 search_acceptors(&top, isize[i], index[i], &hb->a, i, bNitAcc, FALSE, TRUE, datable);
2734                 search_donors   (&top, isize[i], index[i], &hb->d, i, FALSE, TRUE, datable);
2735             }
2736             if (bTwo)
2737             {
2738                 clear_datable_grp(datable, top.atoms.nr);
2739             }
2740         }
2741     }
2742     sfree(datable);
2743     printf("Found %d donors and %d acceptors\n", hb->d.nrd, hb->a.nra);
2744     /*if (bSelected)
2745        snew(donors[gr0D], dons[gr0D].nrd);*/
2746
2747     donor_properties = open_donor_properties_file(opt2fn_null("-don", NFILE, fnm), hb, oenv);
2748
2749     if (bHBmap)
2750     {
2751         printf("Making hbmap structure...");
2752         /* Generate hbond data structure */
2753         mk_hbmap(hb);
2754         printf("done.\n");
2755     }
2756
2757     /* check input */
2758     bStop = FALSE;
2759     if (hb->d.nrd + hb->a.nra == 0)
2760     {
2761         printf("No Donors or Acceptors found\n");
2762         bStop = TRUE;
2763     }
2764     if (!bStop)
2765     {
2766         if (hb->d.nrd == 0)
2767         {
2768             printf("No Donors found\n");
2769             bStop = TRUE;
2770         }
2771         if (hb->a.nra == 0)
2772         {
2773             printf("No Acceptors found\n");
2774             bStop = TRUE;
2775         }
2776     }
2777     if (bStop)
2778     {
2779         gmx_fatal(FARGS, "Nothing to be done");
2780     }
2781
2782     shatom = 0;
2783     if (rshell > 0)
2784     {
2785         int      shisz;
2786         atom_id *shidx;
2787         char    *shgrpnm;
2788         /* get index group with atom for shell */
2789         do
2790         {
2791             printf("Select atom for shell (1 atom):\n");
2792             get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
2793                       1, &shisz, &shidx, &shgrpnm);
2794             if (shisz != 1)
2795             {
2796                 printf("group contains %d atoms, should be 1 (one)\n", shisz);
2797             }
2798         }
2799         while (shisz != 1);
2800         shatom = shidx[0];
2801         printf("Will calculate hydrogen bonds within a shell "
2802                "of %g nm around atom %i\n", rshell, shatom+1);
2803     }
2804
2805     /* Analyze trajectory */
2806     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
2807     if (natoms > top.atoms.nr)
2808     {
2809         gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)",
2810                   top.atoms.nr, natoms);
2811     }
2812
2813     bBox  = ir.ePBC != epbcNONE;
2814     grid  = init_grid(bBox, box, (rcut > r2cut) ? rcut : r2cut, ngrid);
2815     nabin = acut/abin;
2816     nrbin = rcut/rbin;
2817     snew(adist, nabin+1);
2818     snew(rdist, nrbin+1);
2819
2820     bParallel = FALSE;
2821
2822 #ifndef GMX_OPENMP
2823 #define __ADIST adist
2824 #define __RDIST rdist
2825 #define __HBDATA hb
2826 #else /* GMX_OPENMP ================================================== \
2827        * Set up the OpenMP stuff,                                       |
2828        * like the number of threads and such                            |
2829        * Also start the parallel loop.                                  |
2830        */
2831 #define __ADIST p_adist[threadNr]
2832 #define __RDIST p_rdist[threadNr]
2833 #define __HBDATA p_hb[threadNr]
2834 #endif
2835     if (bOMP)
2836     {
2837         bParallel = !bSelected;
2838
2839         if (bParallel)
2840         {
2841             actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
2842
2843             gmx_omp_set_num_threads(actual_nThreads);
2844             printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads);
2845             fflush(stdout);
2846         }
2847         else
2848         {
2849             actual_nThreads = 1;
2850         }
2851
2852         snew(p_hb,    actual_nThreads);
2853         snew(p_adist, actual_nThreads);
2854         snew(p_rdist, actual_nThreads);
2855         for (i = 0; i < actual_nThreads; i++)
2856         {
2857             snew(p_hb[i], 1);
2858             snew(p_adist[i], nabin+1);
2859             snew(p_rdist[i], nrbin+1);
2860
2861             p_hb[i]->max_frames = 0;
2862             p_hb[i]->nhb        = NULL;
2863             p_hb[i]->ndist      = NULL;
2864             p_hb[i]->n_bound    = NULL;
2865             p_hb[i]->time       = NULL;
2866             p_hb[i]->nhx        = NULL;
2867
2868             p_hb[i]->bHBmap     = hb->bHBmap;
2869             p_hb[i]->bDAnr      = hb->bDAnr;
2870             p_hb[i]->wordlen    = hb->wordlen;
2871             p_hb[i]->nframes    = hb->nframes;
2872             p_hb[i]->maxhydro   = hb->maxhydro;
2873             p_hb[i]->danr       = hb->danr;
2874             p_hb[i]->d          = hb->d;
2875             p_hb[i]->a          = hb->a;
2876             p_hb[i]->hbmap      = hb->hbmap;
2877             p_hb[i]->time       = hb->time; /* This may need re-syncing at every frame. */
2878
2879             p_hb[i]->nrhb   = 0;
2880             p_hb[i]->nrdist = 0;
2881         }
2882     }
2883
2884     /* Make a thread pool here,
2885      * instead of forking anew at every frame. */
2886
2887 #pragma omp parallel \
2888     firstprivate(i) \
2889     private(j, h, ii, jj, hh, \
2890     xi, yi, zi, xj, yj, zj, threadNr, \
2891     dist, ang, icell, jcell, \
2892     grp, ogrp, ai, aj, xjj, yjj, zjj, \
2893     xk, yk, zk, ihb, id,  resdist, \
2894     xkk, ykk, zkk, kcell, ak, k, bTric, \
2895     bEdge_xjj, bEdge_yjj) \
2896     default(shared)
2897     {    /* Start of parallel region */
2898         threadNr = gmx_omp_get_thread_num();
2899
2900         do
2901         {
2902
2903             bTric = bBox && TRICLINIC(box);
2904
2905             if (bOMP)
2906             {
2907                 sync_hbdata(p_hb[threadNr], nframes);
2908             }
2909 #pragma omp single
2910             {
2911                 build_grid(hb, x, x[shatom], bBox, box, hbox, (rcut > r2cut) ? rcut : r2cut,
2912                            rshell, ngrid, grid);
2913                 reset_nhbonds(&(hb->d));
2914
2915                 if (debug && bDebug)
2916                 {
2917                     dump_grid(debug, ngrid, grid);
2918                 }
2919
2920                 add_frames(hb, nframes);
2921                 init_hbframe(hb, nframes, output_env_conv_time(oenv, t));
2922
2923                 if (hb->bDAnr)
2924                 {
2925                     count_da_grid(ngrid, grid, hb->danr[nframes]);
2926                 }
2927             } /* omp single */
2928
2929             if (bOMP)
2930             {
2931                 p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */
2932             }
2933
2934             if (bSelected)
2935             {
2936
2937 #pragma omp single
2938                 {
2939                     /* Do not parallelize this just yet. */
2940                     /* int ii; */
2941                     for (ii = 0; (ii < nsel); ii++)
2942                     {
2943                         int dd = index[0][i];
2944                         int aa = index[0][i+2];
2945                         /* int */ hh = index[0][i+1];
2946                         ihb          = is_hbond(hb, ii, ii, dd, aa, rcut, r2cut, ccut, x, bBox, box,
2947                                                 hbox, &dist, &ang, bDA, &h, bContact, bMerge);
2948
2949                         if (ihb)
2950                         {
2951                             /* add to index if not already there */
2952                             /* Add a hbond */
2953                             add_hbond(hb, dd, aa, hh, ii, ii, nframes, bMerge, ihb, bContact);
2954                         }
2955                     }
2956                 } /* omp single */
2957             }     /* if (bSelected) */
2958             else
2959             {
2960                 /* The outer grid loop will have to do for now. */
2961 #pragma omp for schedule(dynamic)
2962                 for (xi = 0; xi < ngrid[XX]; xi++)
2963                 {
2964                     for (yi = 0; (yi < ngrid[YY]); yi++)
2965                     {
2966                         for (zi = 0; (zi < ngrid[ZZ]); zi++)
2967                         {
2968
2969                             /* loop over donor groups gr0 (always) and gr1 (if necessary) */
2970                             for (grp = gr0; (grp <= (bTwo ? gr1 : gr0)); grp++)
2971                             {
2972                                 icell = &(grid[zi][yi][xi].d[grp]);
2973
2974                                 if (bTwo)
2975                                 {
2976                                     ogrp = 1-grp;
2977                                 }
2978                                 else
2979                                 {
2980                                     ogrp = grp;
2981                                 }
2982
2983                                 /* loop over all hydrogen atoms from group (grp)
2984                                  * in this gridcell (icell)
2985                                  */
2986                                 for (ai = 0; (ai < icell->nr); ai++)
2987                                 {
2988                                     i  = icell->atoms[ai];
2989
2990                                     /* loop over all adjacent gridcells (xj,yj,zj) */
2991                                     for (zjj = grid_loop_begin(ngrid[ZZ], zi, bTric, FALSE);
2992                                          zjj <= grid_loop_end(ngrid[ZZ], zi, bTric, FALSE);
2993                                          zjj++)
2994                                     {
2995                                         zj        = grid_mod(zjj, ngrid[ZZ]);
2996                                         bEdge_yjj = (zj == 0) || (zj == ngrid[ZZ] - 1);
2997                                         for (yjj = grid_loop_begin(ngrid[YY], yi, bTric, bEdge_yjj);
2998                                              yjj <= grid_loop_end(ngrid[YY], yi, bTric, bEdge_yjj);
2999                                              yjj++)
3000                                         {
3001                                             yj        = grid_mod(yjj, ngrid[YY]);
3002                                             bEdge_xjj =
3003                                                 (yj == 0) || (yj == ngrid[YY] - 1) ||
3004                                                 (zj == 0) || (zj == ngrid[ZZ] - 1);
3005                                             for (xjj = grid_loop_begin(ngrid[XX], xi, bTric, bEdge_xjj);
3006                                                  xjj <= grid_loop_end(ngrid[XX], xi, bTric, bEdge_xjj);
3007                                                  xjj++)
3008                                             {
3009                                                 xj    = grid_mod(xjj, ngrid[XX]);
3010                                                 jcell = &(grid[zj][yj][xj].a[ogrp]);
3011                                                 /* loop over acceptor atoms from other group (ogrp)
3012                                                  * in this adjacent gridcell (jcell)
3013                                                  */
3014                                                 for (aj = 0; (aj < jcell->nr); aj++)
3015                                                 {
3016                                                     j = jcell->atoms[aj];
3017
3018                                                     /* check if this once was a h-bond */
3019                                                     ihb  = is_hbond(__HBDATA, grp, ogrp, i, j, rcut, r2cut, ccut, x, bBox, box,
3020                                                                     hbox, &dist, &ang, bDA, &h, bContact, bMerge);
3021
3022                                                     if (ihb)
3023                                                     {
3024                                                         /* add to index if not already there */
3025                                                         /* Add a hbond */
3026                                                         add_hbond(__HBDATA, i, j, h, grp, ogrp, nframes, bMerge, ihb, bContact);
3027
3028                                                         /* make angle and distance distributions */
3029                                                         if (ihb == hbHB && !bContact)
3030                                                         {
3031                                                             if (dist > rcut)
3032                                                             {
3033                                                                 gmx_fatal(FARGS, "distance is higher than what is allowed for an hbond: %f", dist);
3034                                                             }
3035                                                             ang *= RAD2DEG;
3036                                                             __ADIST[(int)( ang/abin)]++;
3037                                                             __RDIST[(int)(dist/rbin)]++;
3038                                                             if (!bTwo)
3039                                                             {
3040                                                                 int id, ia;
3041                                                                 if ((id = donor_index(&hb->d, grp, i)) == NOTSET)
3042                                                                 {
3043                                                                     gmx_fatal(FARGS, "Invalid donor %d", i);
3044                                                                 }
3045                                                                 if ((ia = acceptor_index(&hb->a, ogrp, j)) == NOTSET)
3046                                                                 {
3047                                                                     gmx_fatal(FARGS, "Invalid acceptor %d", j);
3048                                                                 }
3049                                                                 resdist = abs(top.atoms.atom[i].resind-
3050                                                                               top.atoms.atom[j].resind);
3051                                                                 if (resdist >= max_hx)
3052                                                                 {
3053                                                                     resdist = max_hx-1;
3054                                                                 }
3055                                                                 __HBDATA->nhx[nframes][resdist]++;
3056                                                             }
3057                                                         }
3058
3059                                                     }
3060                                                 } /* for aj  */
3061                                             }     /* for xjj */
3062                                         }         /* for yjj */
3063                                     }             /* for zjj */
3064                                 }                 /* for ai  */
3065                             }                     /* for grp */
3066                         }                         /* for xi,yi,zi */
3067                     }
3068                 }
3069             } /* if (bSelected) {...} else */
3070
3071
3072             /* Better wait for all threads to finnish using x[] before updating it. */
3073             k = nframes;
3074 #pragma omp barrier
3075 #pragma omp critical
3076             {
3077                 /* Sum up histograms and counts from p_hb[] into hb */
3078                 if (bOMP)
3079                 {
3080                     hb->nhb[k]   += p_hb[threadNr]->nhb[k];
3081                     hb->ndist[k] += p_hb[threadNr]->ndist[k];
3082                     for (j = 0; j < max_hx; j++)
3083                     {
3084                         hb->nhx[k][j]  += p_hb[threadNr]->nhx[k][j];
3085                     }
3086                 }
3087             }
3088
3089             /* Here are a handful of single constructs
3090              * to share the workload a bit. The most
3091              * important one is of course the last one,
3092              * where there's a potential bottleneck in form
3093              * of slow I/O.                    */
3094 #pragma omp barrier
3095 #pragma omp single
3096             {
3097                 analyse_donor_properties(donor_properties, hb, k, t);
3098             }
3099
3100 #pragma omp single
3101             {
3102                 if (fpnhb)
3103                 {
3104                     do_nhb_dist(fpnhb, hb, t);
3105                 }
3106             }
3107
3108 #pragma omp single
3109             {
3110                 trrStatus = (read_next_x(oenv, status, &t, x, box));
3111                 nframes++;
3112             }
3113
3114 #pragma omp barrier
3115         }
3116         while (trrStatus);
3117
3118         if (bOMP)
3119         {
3120 #pragma omp critical
3121             {
3122                 hb->nrhb   += p_hb[threadNr]->nrhb;
3123                 hb->nrdist += p_hb[threadNr]->nrdist;
3124             }
3125             /* Free parallel datastructures */
3126             sfree(p_hb[threadNr]->nhb);
3127             sfree(p_hb[threadNr]->ndist);
3128             sfree(p_hb[threadNr]->nhx);
3129
3130 #pragma omp for
3131             for (i = 0; i < nabin; i++)
3132             {
3133                 for (j = 0; j < actual_nThreads; j++)
3134                 {
3135
3136                     adist[i] += p_adist[j][i];
3137                 }
3138             }
3139 #pragma omp for
3140             for (i = 0; i <= nrbin; i++)
3141             {
3142                 for (j = 0; j < actual_nThreads; j++)
3143                 {
3144                     rdist[i] += p_rdist[j][i];
3145                 }
3146             }
3147
3148             sfree(p_adist[threadNr]);
3149             sfree(p_rdist[threadNr]);
3150         }
3151     } /* End of parallel region */
3152     if (bOMP)
3153     {
3154         sfree(p_adist);
3155         sfree(p_rdist);
3156     }
3157
3158     if (nframes < 2 && (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm)))
3159     {
3160         gmx_fatal(FARGS, "Cannot calculate autocorrelation of life times with less than two frames");
3161     }
3162
3163     free_grid(ngrid, &grid);
3164
3165     close_trj(status);
3166
3167     if (donor_properties)
3168     {
3169         xvgrclose(donor_properties);
3170     }
3171
3172     if (fpnhb)
3173     {
3174         xvgrclose(fpnhb);
3175     }
3176
3177     /* Compute maximum possible number of different hbonds */
3178     if (maxnhb > 0)
3179     {
3180         max_nhb = maxnhb;
3181     }
3182     else
3183     {
3184         max_nhb = 0.5*(hb->d.nrd*hb->a.nra);
3185     }
3186
3187     if (bHBmap)
3188     {
3189         if (hb->nrhb == 0)
3190         {
3191             printf("No %s found!!\n", bContact ? "contacts" : "hydrogen bonds");
3192         }
3193         else
3194         {
3195             printf("Found %d different %s in trajectory\n"
3196                    "Found %d different atom-pairs within %s distance\n",
3197                    hb->nrhb, bContact ? "contacts" : "hydrogen bonds",
3198                    hb->nrdist, (r2cut > 0) ? "second cut-off" : "hydrogen bonding");
3199
3200             if (bMerge)
3201             {
3202                 merge_hb(hb, bTwo, bContact);
3203             }
3204
3205             if (opt2bSet("-hbn", NFILE, fnm))
3206             {
3207                 dump_hbmap(hb, NFILE, fnm, bTwo, bContact, isize, index, grpnames, &top.atoms);
3208             }
3209
3210             /* Moved the call to merge_hb() to a line BEFORE dump_hbmap
3211              * to make the -hbn and -hmb output match eachother.
3212              * - Erik Marklund, May 30, 2006 */
3213         }
3214     }
3215     /* Print out number of hbonds and distances */
3216     aver_nhb  = 0;
3217     aver_dist = 0;
3218     fp        = xvgropen(opt2fn("-num", NFILE, fnm), bContact ? "Contacts" :
3219                          "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
3220     snew(leg, 2);
3221     snew(leg[0], STRLEN);
3222     snew(leg[1], STRLEN);
3223     sprintf(leg[0], "%s", bContact ? "Contacts" : "Hydrogen bonds");
3224     sprintf(leg[1], "Pairs within %g nm", (r2cut > 0) ? r2cut : rcut);
3225     xvgr_legend(fp, 2, (const char**)leg, oenv);
3226     sfree(leg[1]);
3227     sfree(leg[0]);
3228     sfree(leg);
3229     for (i = 0; (i < nframes); i++)
3230     {
3231         fprintf(fp, "%10g  %10d  %10d\n", hb->time[i], hb->nhb[i], hb->ndist[i]);
3232         aver_nhb  += hb->nhb[i];
3233         aver_dist += hb->ndist[i];
3234     }
3235     xvgrclose(fp);
3236     aver_nhb  /= nframes;
3237     aver_dist /= nframes;
3238     /* Print HB distance distribution */
3239     if (opt2bSet("-dist", NFILE, fnm))
3240     {
3241         long sum;
3242
3243         sum = 0;
3244         for (i = 0; i < nrbin; i++)
3245         {
3246             sum += rdist[i];
3247         }
3248
3249         fp = xvgropen(opt2fn("-dist", NFILE, fnm),
3250                       "Hydrogen Bond Distribution",
3251                       bDA ?
3252                       "Donor - Acceptor Distance (nm)" :
3253                       "Hydrogen - Acceptor Distance (nm)", "", oenv);
3254         for (i = 0; i < nrbin; i++)
3255         {
3256             fprintf(fp, "%10g %10g\n", (i+0.5)*rbin, rdist[i]/(rbin*(real)sum));
3257         }
3258         xvgrclose(fp);
3259     }
3260
3261     /* Print HB angle distribution */
3262     if (opt2bSet("-ang", NFILE, fnm))
3263     {
3264         long sum;
3265
3266         sum = 0;
3267         for (i = 0; i < nabin; i++)
3268         {
3269             sum += adist[i];
3270         }
3271
3272         fp = xvgropen(opt2fn("-ang", NFILE, fnm),
3273                       "Hydrogen Bond Distribution",
3274                       "Hydrogen - Donor - Acceptor Angle (\\SO\\N)", "", oenv);
3275         for (i = 0; i < nabin; i++)
3276         {
3277             fprintf(fp, "%10g %10g\n", (i+0.5)*abin, adist[i]/(abin*(real)sum));
3278         }
3279         xvgrclose(fp);
3280     }
3281
3282     /* Print HB in alpha-helix */
3283     if (opt2bSet("-hx", NFILE, fnm))
3284     {
3285         fp = xvgropen(opt2fn("-hx", NFILE, fnm),
3286                       "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
3287         xvgr_legend(fp, NRHXTYPES, hxtypenames, oenv);
3288         for (i = 0; i < nframes; i++)
3289         {
3290             fprintf(fp, "%10g", hb->time[i]);
3291             for (j = 0; j < max_hx; j++)
3292             {
3293                 fprintf(fp, " %6d", hb->nhx[i][j]);
3294             }
3295             fprintf(fp, "\n");
3296         }
3297         xvgrclose(fp);
3298     }
3299
3300     printf("Average number of %s per timeframe %.3f out of %g possible\n",
3301            bContact ? "contacts" : "hbonds",
3302            bContact ? aver_dist : aver_nhb, max_nhb);
3303
3304     /* Do Autocorrelation etc. */
3305     if (hb->bHBmap)
3306     {
3307         /*
3308            Added support for -contact in ac and hbm calculations below.
3309            - Erik Marklund, May 29, 2006
3310          */
3311         ivec itmp;
3312         rvec rtmp;
3313         if (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm))
3314         {
3315             please_cite(stdout, "Spoel2006b");
3316         }
3317         if (opt2bSet("-ac", NFILE, fnm))
3318         {
3319             char *gemstring = NULL;
3320
3321             do_hbac(opt2fn("-ac", NFILE, fnm), hb, nDump,
3322                     bMerge, bContact, fit_start, temp, r2cut > 0, oenv,
3323                     nThreads);
3324         }
3325         if (opt2bSet("-life", NFILE, fnm))
3326         {
3327             do_hblife(opt2fn("-life", NFILE, fnm), hb, bMerge, bContact, oenv);
3328         }
3329         if (opt2bSet("-hbm", NFILE, fnm))
3330         {
3331             t_matrix mat;
3332             int      id, ia, hh, x, y;
3333             mat.flags = mat.y0 = 0;
3334
3335             if ((nframes > 0) && (hb->nrhb > 0))
3336             {
3337                 mat.nx = nframes;
3338                 mat.ny = hb->nrhb;
3339
3340                 snew(mat.matrix, mat.nx);
3341                 for (x = 0; (x < mat.nx); x++)
3342                 {
3343                     snew(mat.matrix[x], mat.ny);
3344                 }
3345                 y = 0;
3346                 for (id = 0; (id < hb->d.nrd); id++)
3347                 {
3348                     for (ia = 0; (ia < hb->a.nra); ia++)
3349                     {
3350                         for (hh = 0; (hh < hb->maxhydro); hh++)
3351                         {
3352                             if (hb->hbmap[id][ia])
3353                             {
3354                                 if (ISHB(hb->hbmap[id][ia]->history[hh]))
3355                                 {
3356                                     for (x = 0; (x <= hb->hbmap[id][ia]->nframes); x++)
3357                                     {
3358                                         int nn0 = hb->hbmap[id][ia]->n0;
3359                                         range_check(y, 0, mat.ny);
3360                                         mat.matrix[x+nn0][y] = is_hb(hb->hbmap[id][ia]->h[hh], x);
3361                                     }
3362                                     y++;
3363                                 }
3364                             }
3365                         }
3366                     }
3367                 }
3368                 mat.axis_x = hb->time;
3369                 snew(mat.axis_y, mat.ny);
3370                 for (j = 0; j < mat.ny; j++)
3371                 {
3372                     mat.axis_y[j] = j;
3373                 }
3374                 sprintf(mat.title, bContact ? "Contact Existence Map" :
3375                         "Hydrogen Bond Existence Map");
3376                 sprintf(mat.legend, bContact ? "Contacts" : "Hydrogen Bonds");
3377                 sprintf(mat.label_x, "%s", output_env_get_xvgr_tlabel(oenv));
3378                 sprintf(mat.label_y, bContact ? "Contact Index" : "Hydrogen Bond Index");
3379                 mat.bDiscrete = TRUE;
3380                 mat.nmap      = 2;
3381                 snew(mat.map, mat.nmap);
3382                 for (i = 0; i < mat.nmap; i++)
3383                 {
3384                     mat.map[i].code.c1 = hbmap[i];
3385                     mat.map[i].desc    = hbdesc[i];
3386                     mat.map[i].rgb     = hbrgb[i];
3387                 }
3388                 fp = opt2FILE("-hbm", NFILE, fnm, "w");
3389                 write_xpm_m(fp, mat);
3390                 gmx_ffclose(fp);
3391                 for (x = 0; x < mat.nx; x++)
3392                 {
3393                     sfree(mat.matrix[x]);
3394                 }
3395                 sfree(mat.axis_y);
3396                 sfree(mat.matrix);
3397                 sfree(mat.map);
3398             }
3399             else
3400             {
3401                 fprintf(stderr, "No hydrogen bonds/contacts found. No hydrogen bond map will be printed.\n");
3402             }
3403         }
3404     }
3405
3406     if (hb->bDAnr)
3407     {
3408         int    i, j, nleg;
3409         char **legnames;
3410         char   buf[STRLEN];
3411
3412 #define USE_THIS_GROUP(j) ( (j == gr0) || (bTwo && (j == gr1)) )
3413
3414         fp = xvgropen(opt2fn("-dan", NFILE, fnm),
3415                       "Donors and Acceptors", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
3416         nleg = (bTwo ? 2 : 1)*2;
3417         snew(legnames, nleg);
3418         i = 0;
3419         for (j = 0; j < grNR; j++)
3420         {
3421             if (USE_THIS_GROUP(j) )
3422             {
3423                 sprintf(buf, "Donors %s", grpnames[j]);
3424                 legnames[i++] = gmx_strdup(buf);
3425                 sprintf(buf, "Acceptors %s", grpnames[j]);
3426                 legnames[i++] = gmx_strdup(buf);
3427             }
3428         }
3429         if (i != nleg)
3430         {
3431             gmx_incons("number of legend entries");
3432         }
3433         xvgr_legend(fp, nleg, (const char**)legnames, oenv);
3434         for (i = 0; i < nframes; i++)
3435         {
3436             fprintf(fp, "%10g", hb->time[i]);
3437             for (j = 0; (j < grNR); j++)
3438             {
3439                 if (USE_THIS_GROUP(j) )
3440                 {
3441                     fprintf(fp, " %6d", hb->danr[i][j]);
3442                 }
3443             }
3444             fprintf(fp, "\n");
3445         }
3446         xvgrclose(fp);
3447     }
3448
3449     return 0;
3450 }