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