src/gromacs/gmxana/gmx_hbond.cpp

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