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