src/gromacs/gmxana/gmx_hbond.cpp

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