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