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