src/gromacs/gmxana/gmx_hbond.cpp

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