src/gromacs/gmxana/gmx_hbond.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.3.3
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2008, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * Groningen Machine for Chemical Simulation
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39 #include <math.h>
  40
  41 /*#define HAVE_NN_LOOPS*/
  42
  43 #include "statutil.h"
  44 #include "copyrite.h"
  45 #include "sysstuff.h"
  46 #include "txtdump.h"
  47 #include "physics.h"
  48 #include "macros.h"
  49 #include "gmx_fatal.h"
  50 #include "index.h"
  51 #include "smalloc.h"
  52 #include "vec.h"
  53 #include "xvgr.h"
  54 #include "gstat.h"
  55 #include "string2.h"
  56 #include "pbc.h"
  57 #include "correl.h"
  58 #include "gmx_ana.h"
  59 #include "geminate.h"
  60
  61 #include "gromacs/fileio/futil.h"
  62 #include "gromacs/fileio/matio.h"
  63 #include "gromacs/fileio/tpxio.h"
  64 #include "gromacs/fileio/trxio.h"
  65 #include "gromacs/utility/gmxomp.h"
  66
  67 typedef short int t_E;
  68 typedef int t_EEst;
  69 #define max_hx 7
  70 typedef int t_hx[max_hx];
  71 #define NRHXTYPES max_hx
  72 const char *hxtypenames[NRHXTYPES] =
  73 {"n-n", "n-n+1", "n-n+2", "n-n+3", "n-n+4", "n-n+5", "n-n>6"};
  74 #define MAXHH 4
  75
  76 #ifdef GMX_OPENMP
  77 #define MASTER_THREAD_ONLY(threadNr) ((threadNr) == 0)
  78 #else
  79 #define MASTER_THREAD_ONLY(threadNr) ((threadNr) == (threadNr))
  80 #endif
  81
  82 /* -----------------------------------------*/
  83
  84 enum {
  85     gr0,  gr1,    grI,  grNR
  86 };
  87 enum {
  88     hbNo, hbDist, hbHB, hbNR, hbR2
  89 };
  90 enum {
  91     noDA, ACC, DON, DA, INGROUP
  92 };
  93 enum {
  94     NN_NULL, NN_NONE, NN_BINARY, NN_1_over_r3, NN_dipole, NN_NR
  95 };
  96
  97 static const char *grpnames[grNR] = {"0", "1", "I" };
  98
  99 static gmx_bool    bDebug = FALSE;
 100
 101 #define HB_NO 0
 102 #define HB_YES 1<<0
 103 #define HB_INS 1<<1
 104 #define HB_YESINS HB_YES|HB_INS
 105 #define HB_NR (1<<2)
 106 #define MAXHYDRO 4
 107
 108 #define ISHB(h)   (((h) & 2) == 2)
 109 #define ISDIST(h) (((h) & 1) == 1)
 110 #define ISDIST2(h) (((h) & 4) == 4)
 111 #define ISACC(h)  (((h) & 1) == 1)
 112 #define ISDON(h)  (((h) & 2) == 2)
 113 #define ISINGRP(h) (((h) & 4) == 4)
 114
 115 typedef struct {
 116     int      nr;
 117     int      maxnr;
 118     atom_id *atoms;
 119 } t_ncell;
 120
 121 typedef struct {
 122     t_ncell d[grNR];
 123     t_ncell a[grNR];
 124 } t_gridcell;
 125
 126 typedef int     t_icell[grNR];
 127 typedef atom_id h_id[MAXHYDRO];
 128
 129 typedef struct {
 130     int      history[MAXHYDRO];
 131     /* Has this hbond existed ever? If so as hbDist or hbHB or both.
 132      * Result is stored as a bitmap (1 = hbDist) || (2 = hbHB)
 133      */
 134     /* Bitmask array which tells whether a hbond is present
 135      * at a given time. Either of these may be NULL
 136      */
 137     int            n0;                 /* First frame a HB was found     */
 138     int            nframes, maxframes; /* Amount of frames in this hbond */
 139     unsigned int **h;
 140     unsigned int **g;
 141     /* See Xu and Berne, JPCB 105 (2001), p. 11929. We define the
 142      * function g(t) = [1-h(t)] H(t) where H(t) is one when the donor-
 143      * acceptor distance is less than the user-specified distance (typically
 144      * 0.35 nm).
 145      */
 146 } t_hbond;
 147
 148 typedef struct {
 149     int      nra, max_nra;
 150     atom_id *acc;             /* Atom numbers of the acceptors     */
 151     int     *grp;             /* Group index                       */
 152     int     *aptr;            /* Map atom number to acceptor index */
 153 } t_acceptors;
 154
 155 typedef struct {
 156     int       nrd, max_nrd;
 157     int      *don;               /* Atom numbers of the donors         */
 158     int      *grp;               /* Group index                        */
 159     int      *dptr;              /* Map atom number to donor index     */
 160     int      *nhydro;            /* Number of hydrogens for each donor */
 161     h_id     *hydro;             /* The atom numbers of the hydrogens  */
 162     h_id     *nhbonds;           /* The number of HBs per H at current */
 163 } t_donors;
 164
 165 /* Tune this to match memory requirements. It should be a signed integer type, e.g. signed char.*/
 166 #define PSTYPE int
 167
 168 typedef struct {
 169     int     len;   /* The length of frame and p. */
 170     int    *frame; /* The frames at which transitio*/
 171     PSTYPE *p;
 172 } t_pShift;
 173
 174 typedef struct {
 175     /* Periodicity history. Used for the reversible geminate recombination. */
 176     t_pShift **pHist; /* The periodicity of every hbond in t_hbdata->hbmap:
 177                        *   pHist[d][a]. We can safely assume that the same
 178                        *   periodic shift holds for all hydrogens of a da-pair.
 179                        *
 180                        * Nowadays it only stores TRANSITIONS, and not the shift at every frame.
 181                        *   That saves a LOT of memory, an hopefully kills a mysterious bug where
 182                        *   pHist gets contaminated. */
 183
 184     PSTYPE nper;      /* The length of p2i */
 185     ivec  *p2i;       /* Maps integer to periodic shift for a pair.*/
 186     matrix P;         /* Projection matrix to find the box shifts. */
 187     int    gemtype;   /* enumerated type */
 188 } t_gemPeriod;
 189
 190 typedef struct {
 191     int     nframes;
 192     int    *Etot; /* Total energy for each frame */
 193     t_E ****E;    /* Energy estimate for [d][a][h][frame-n0] */
 194 } t_hbEmap;
 195
 196 typedef struct {
 197     gmx_bool        bHBmap, bDAnr, bGem;
 198     int             wordlen;
 199     /* The following arrays are nframes long */
 200     int             nframes, max_frames, maxhydro;
 201     int            *nhb, *ndist;
 202     h_id           *n_bound;
 203     real           *time;
 204     t_icell        *danr;
 205     t_hx           *nhx;
 206     /* These structures are initialized from the topology at start up */
 207     t_donors        d;
 208     t_acceptors     a;
 209     /* This holds a matrix with all possible hydrogen bonds */
 210     int             nrhb, nrdist;
 211     t_hbond      ***hbmap;
 212 #ifdef HAVE_NN_LOOPS
 213     t_hbEmap        hbE;
 214 #endif
 215     /* For parallelization reasons this will have to be a pointer.
 216      * Otherwise discrepancies may arise between the periodicity data
 217      * seen by different threads. */
 218     t_gemPeriod *per;
 219 } t_hbdata;
 220
 221 static void clearPshift(t_pShift *pShift)
 222 {
 223     if (pShift->len > 0)
 224     {
 225         sfree(pShift->p);
 226         sfree(pShift->frame);
 227         pShift->len = 0;
 228     }
 229 }
 230
 231 static void calcBoxProjection(matrix B, matrix P)
 232 {
 233     const int vp[] = {XX, YY, ZZ};
 234     int       i, j;
 235     int       m, n;
 236     matrix    M, N, U;
 237
 238     for (i = 0; i < 3; i++)
 239     {
 240         m = vp[i];
 241         for (j = 0; j < 3; j++)
 242         {
 243             n       = vp[j];
 244             U[m][n] = i == j ? 1 : 0;
 245         }
 246     }
 247     m_inv(B, M);
 248     for (i = 0; i < 3; i++)
 249     {
 250         m = vp[i];
 251         mvmul(M, U[m], P[m]);
 252     }
 253     transpose(P, N);
 254 }
 255
 256 static void calcBoxDistance(matrix P, rvec d, ivec ibd)
 257 {
 258     /* returns integer distance in box coordinates.
 259      * P is the projection matrix from cartesian coordinates
 260      * obtained with calcBoxProjection(). */
 261     int  i;
 262     rvec bd;
 263     mvmul(P, d, bd);
 264     /* extend it by 0.5 in all directions since (int) rounds toward 0.*/
 265     for (i = 0; i < 3; i++)
 266     {
 267         bd[i] = bd[i] + (bd[i] < 0 ? -0.5 : 0.5);
 268     }
 269     ibd[XX] = (int)bd[XX];
 270     ibd[YY] = (int)bd[YY];
 271     ibd[ZZ] = (int)bd[ZZ];
 272 }
 273
 274 /* Changed argument 'bMerge' into 'oneHB' below,
 275  * since -contact should cause maxhydro to be 1,
 276  * not just -merge.
 277  * - Erik Marklund May 29, 2006
 278  */
 279
 280 static PSTYPE periodicIndex(ivec r, t_gemPeriod *per, gmx_bool daSwap)
 281 {
 282     /* Try to merge hbonds on the fly. That means that if the
 283      * acceptor and donor are mergable, then:
 284      * 1) store the hb-info so that acceptor id > donor id,
 285      * 2) add the periodic shift in pairs, so that [-x,-y,-z] is
 286      *    stored in per.p2i[] whenever acceptor id < donor id.
 287      * Note that [0,0,0] should already be the first element of per.p2i
 288      * by the time this function is called. */
 289
 290     /* daSwap is TRUE if the donor and acceptor were swapped.
 291      * If so, then the negative vector should be used. */
 292     PSTYPE i;
 293
 294     if (per->p2i == NULL || per->nper == 0)
 295     {
 296         gmx_fatal(FARGS, "'per' not initialized properly.");
 297     }
 298     for (i = 0; i < per->nper; i++)
 299     {
 300         if (r[XX] == per->p2i[i][XX] &&
 301             r[YY] == per->p2i[i][YY] &&
 302             r[ZZ] == per->p2i[i][ZZ])
 303         {
 304             return i;
 305         }
 306     }
 307     /* Not found apparently. Add it to the list! */
 308     /* printf("New shift found: %i,%i,%i\n",r[XX],r[YY],r[ZZ]); */
 309
 310 #pragma omp critical
 311     {
 312         if (!per->p2i)
 313         {
 314             fprintf(stderr, "p2i not initialized. This shouldn't happen!\n");
 315             snew(per->p2i, 1);
 316         }
 317         else
 318         {
 319             srenew(per->p2i, per->nper+2);
 320         }
 321         copy_ivec(r, per->p2i[per->nper]);
 322         (per->nper)++;
 323
 324         /* Add the mirror too. It's rather likely that it'll be needed. */
 325         per->p2i[per->nper][XX] = -r[XX];
 326         per->p2i[per->nper][YY] = -r[YY];
 327         per->p2i[per->nper][ZZ] = -r[ZZ];
 328         (per->nper)++;
 329     } /* omp critical */
 330     return per->nper - 1 - (daSwap ? 0 : 1);
 331 }
 332
 333 static t_hbdata *mk_hbdata(gmx_bool bHBmap, gmx_bool bDAnr, gmx_bool oneHB, gmx_bool bGem, int gemmode)
 334 {
 335     t_hbdata *hb;
 336
 337     snew(hb, 1);
 338     hb->wordlen = 8*sizeof(unsigned int);
 339     hb->bHBmap  = bHBmap;
 340     hb->bDAnr   = bDAnr;
 341     hb->bGem    = bGem;
 342     if (oneHB)
 343     {
 344         hb->maxhydro = 1;
 345     }
 346     else
 347     {
 348         hb->maxhydro = MAXHYDRO;
 349     }
 350     snew(hb->per, 1);
 351     hb->per->gemtype = bGem ? gemmode : 0;
 352
 353     return hb;
 354 }
 355
 356 static void mk_hbmap(t_hbdata *hb)
 357 {
 358     int  i, j;
 359
 360     snew(hb->hbmap, hb->d.nrd);
 361     for (i = 0; (i < hb->d.nrd); i++)
 362     {
 363         snew(hb->hbmap[i], hb->a.nra);
 364         if (hb->hbmap[i] == NULL)
 365         {
 366             gmx_fatal(FARGS, "Could not allocate enough memory for hbmap");
 367         }
 368         for (j = 0; (j > hb->a.nra); j++)
 369         {
 370             hb->hbmap[i][j] = NULL;
 371         }
 372     }
 373 }
 374
 375 /* Consider redoing pHist so that is only stores transitions between
 376  * periodicities and not the periodicity for all frames. This eats heaps of memory. */
 377 static void mk_per(t_hbdata *hb)
 378 {
 379     int i, j;
 380     if (hb->bGem)
 381     {
 382         snew(hb->per->pHist, hb->d.nrd);
 383         for (i = 0; i < hb->d.nrd; i++)
 384         {
 385             snew(hb->per->pHist[i], hb->a.nra);
 386             if (hb->per->pHist[i] == NULL)
 387             {
 388                 gmx_fatal(FARGS, "Could not allocate enough memory for per->pHist");
 389             }
 390             for (j = 0; j < hb->a.nra; j++)
 391             {
 392                 clearPshift(&(hb->per->pHist[i][j]));
 393             }
 394         }
 395         /* add the [0,0,0] shift to element 0 of p2i. */
 396         snew(hb->per->p2i, 1);
 397         clear_ivec(hb->per->p2i[0]);
 398         hb->per->nper = 1;
 399     }
 400 }
 401
 402 #ifdef HAVE_NN_LOOPS
 403 static void mk_hbEmap (t_hbdata *hb, int n0)
 404 {
 405     int i, j, k;
 406     hb->hbE.E       = NULL;
 407     hb->hbE.nframes = 0;
 408     snew(hb->hbE.E, hb->d.nrd);
 409     for (i = 0; i < hb->d.nrd; i++)
 410     {
 411         snew(hb->hbE.E[i], hb->a.nra);
 412         for (j = 0; j < hb->a.nra; j++)
 413         {
 414             snew(hb->hbE.E[i][j], MAXHYDRO);
 415             for (k = 0; k < MAXHYDRO; k++)
 416             {
 417                 hb->hbE.E[i][j][k] = NULL;
 418             }
 419         }
 420     }
 421     hb->hbE.Etot = NULL;
 422 }
 423
 424 static void free_hbEmap (t_hbdata *hb)
 425 {
 426     int i, j, k;
 427     for (i = 0; i < hb->d.nrd; i++)
 428     {
 429         for (j = 0; j < hb->a.nra; j++)
 430         {
 431             for (k = 0; k < MAXHYDRO; k++)
 432             {
 433                 sfree(hb->hbE.E[i][j][k]);
 434             }
 435             sfree(hb->hbE.E[i][j]);
 436         }
 437         sfree(hb->hbE.E[i]);
 438     }
 439     sfree(hb->hbE.E);
 440     sfree(hb->hbE.Etot);
 441 }
 442
 443 static void addFramesNN(t_hbdata *hb, int frame)
 444 {
 445
 446 #define DELTAFRAMES_HBE 10
 447
 448     int d, a, h, nframes;
 449
 450     if (frame >= hb->hbE.nframes)
 451     {
 452         nframes =  hb->hbE.nframes + DELTAFRAMES_HBE;
 453         srenew(hb->hbE.Etot, nframes);
 454
 455         for (d = 0; d < hb->d.nrd; d++)
 456         {
 457             for (a = 0; a < hb->a.nra; a++)
 458             {
 459                 for (h = 0; h < hb->d.nhydro[d]; h++)
 460                 {
 461                     srenew(hb->hbE.E[d][a][h], nframes);
 462                 }
 463             }
 464         }
 465
 466         hb->hbE.nframes += DELTAFRAMES_HBE;
 467     }
 468 }
 469
 470 static t_E calcHbEnergy(int d, int a, int h, rvec x[], t_EEst EEst,
 471                         matrix box, rvec hbox, t_donors *donors)
 472 {
 473     /* d     - donor atom
 474      * a     - acceptor atom
 475      * h     - hydrogen
 476      * alpha - angle between dipoles
 477      * x[]   - atomic positions
 478      * EEst  - the type of energy estimate (see enum in hbplugin.h)
 479      * box   - the box vectors   \
 480      * hbox  - half box lengths  _These two are only needed for the pbc correction
 481      */
 482
 483     t_E  E;
 484     rvec dist;
 485     rvec dipole[2], xmol[3], xmean[2];
 486     int  i;
 487     real r, realE;
 488
 489     if (d == a)
 490     {
 491         /* Self-interaction */
 492         return NONSENSE_E;
 493     }
 494
 495     switch (EEst)
 496     {
 497         case NN_BINARY:
 498             /* This is a simple binary existence function that sets E=1 whenever
 499              * the distance between the oxygens is equal too or less than 0.35 nm.
 500              */
 501             rvec_sub(x[d], x[a], dist);
 502             pbc_correct_gem(dist, box, hbox);
 503             if (norm(dist) <= 0.35)
 504             {
 505                 E = 1;
 506             }
 507             else
 508             {
 509                 E = 0;
 510             }
 511             break;
 512
 513         case NN_1_over_r3:
 514             /* Negative potential energy of a dipole.
 515              * E = -cos(alpha) * 1/r^3 */
 516
 517             copy_rvec(x[d], xmol[0]);                                 /* donor */
 518             copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */
 519             copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */
 520
 521             svmul(15.9994*(1/1.008), xmol[0], xmean[0]);
 522             rvec_inc(xmean[0], xmol[1]);
 523             rvec_inc(xmean[0], xmol[2]);
 524             for (i = 0; i < 3; i++)
 525             {
 526                 xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008;
 527             }
 528
 529             /* Assumes that all acceptors are also donors. */
 530             copy_rvec(x[a], xmol[0]);                                 /* acceptor */
 531             copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */
 532             copy_rvec(x[donors->hydro[donors->dptr[a]][1]], xmol[2]); /* hydrogen */
 533
 534
 535             svmul(15.9994*(1/1.008), xmol[0], xmean[1]);
 536             rvec_inc(xmean[1], xmol[1]);
 537             rvec_inc(xmean[1], xmol[2]);
 538             for (i = 0; i < 3; i++)
 539             {
 540                 xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008;
 541             }
 542
 543             rvec_sub(xmean[0], xmean[1], dist);
 544             pbc_correct_gem(dist, box, hbox);
 545             r = norm(dist);
 546
 547             realE = pow(r, -3.0);
 548             E     = (t_E)(SCALEFACTOR_E * realE);
 549             break;
 550
 551         case NN_dipole:
 552             /* Negative potential energy of a (unpolarizable) dipole.
 553              * E = -cos(alpha) * 1/r^3 */
 554             clear_rvec(dipole[1]);
 555             clear_rvec(dipole[0]);
 556
 557             copy_rvec(x[d], xmol[0]);                                 /* donor */
 558             copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */
 559             copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */
 560
 561             rvec_inc(dipole[0], xmol[1]);
 562             rvec_inc(dipole[0], xmol[2]);
 563             for (i = 0; i < 3; i++)
 564             {
 565                 dipole[0][i] *= 0.5;
 566             }
 567             rvec_dec(dipole[0], xmol[0]);
 568
 569             svmul(15.9994*(1/1.008), xmol[0], xmean[0]);
 570             rvec_inc(xmean[0], xmol[1]);
 571             rvec_inc(xmean[0], xmol[2]);
 572             for (i = 0; i < 3; i++)
 573             {
 574                 xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008;
 575             }
 576
 577             /* Assumes that all acceptors are also donors. */
 578             copy_rvec(x[a], xmol[0]);                                 /* acceptor */
 579             copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */
 580             copy_rvec(x[donors->hydro[donors->dptr[a]][2]], xmol[2]); /* hydrogen */
 581
 582
 583             rvec_inc(dipole[1], xmol[1]);
 584             rvec_inc(dipole[1], xmol[2]);
 585             for (i = 0; i < 3; i++)
 586             {
 587                 dipole[1][i] *= 0.5;
 588             }
 589             rvec_dec(dipole[1], xmol[0]);
 590
 591             svmul(15.9994*(1/1.008), xmol[0], xmean[1]);
 592             rvec_inc(xmean[1], xmol[1]);
 593             rvec_inc(xmean[1], xmol[2]);
 594             for (i = 0; i < 3; i++)
 595             {
 596                 xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008;
 597             }
 598
 599             rvec_sub(xmean[0], xmean[1], dist);
 600             pbc_correct_gem(dist, box, hbox);
 601             r = norm(dist);
 602
 603             double cosalpha = cos_angle(dipole[0], dipole[1]);
 604             realE = cosalpha * pow(r, -3.0);
 605             E     = (t_E)(SCALEFACTOR_E * realE);
 606             break;
 607
 608         default:
 609             printf("Can't do that type of energy estimate: %i\n.", EEst);
 610             E = NONSENSE_E;
 611     }
 612
 613     return E;
 614 }
 615
 616 static void storeHbEnergy(t_hbdata *hb, int d, int a, int h, t_E E, int frame)
 617 {
 618     /* hb - hbond data structure
 619        d  - donor
 620        a  - acceptor
 621        h  - hydrogen
 622        E  - estimate of the energy
 623        frame - the current frame.
 624      */
 625
 626     /* Store the estimated energy */
 627     if (E == NONSENSE_E)
 628     {
 629         E = 0;
 630     }
 631
 632     hb->hbE.E[d][a][h][frame] = E;
 633
 634 #pragma omp critical
 635     {
 636         hb->hbE.Etot[frame] += E;
 637     }
 638 }
 639 #endif /* HAVE_NN_LOOPS */
 640
 641
 642 /* Finds -v[] in the periodicity index */
 643 static int findMirror(PSTYPE p, ivec v[], PSTYPE nper)
 644 {
 645     PSTYPE i;
 646     ivec   u;
 647     for (i = 0; i < nper; i++)
 648     {
 649         if (v[i][XX] == -(v[p][XX]) &&
 650             v[i][YY] == -(v[p][YY]) &&
 651             v[i][ZZ] == -(v[p][ZZ]))
 652         {
 653             return (int)i;
 654         }
 655     }
 656     printf("Couldn't find mirror of [%i, %i, %i], index \n",
 657            v[p][XX],
 658            v[p][YY],
 659            v[p][ZZ]);
 660     return -1;
 661 }
 662
 663
 664 static void add_frames(t_hbdata *hb, int nframes)
 665 {
 666     int  i, j, k, l;
 667
 668     if (nframes >= hb->max_frames)
 669     {
 670         hb->max_frames += 4096;
 671         srenew(hb->time, hb->max_frames);
 672         srenew(hb->nhb, hb->max_frames);
 673         srenew(hb->ndist, hb->max_frames);
 674         srenew(hb->n_bound, hb->max_frames);
 675         srenew(hb->nhx, hb->max_frames);
 676         if (hb->bDAnr)
 677         {
 678             srenew(hb->danr, hb->max_frames);
 679         }
 680     }
 681     hb->nframes = nframes;
 682 }
 683
 684 #define OFFSET(frame) (frame / 32)
 685 #define MASK(frame)   (1 << (frame % 32))
 686
 687 static void _set_hb(unsigned int hbexist[], unsigned int frame, gmx_bool bValue)
 688 {
 689     if (bValue)
 690     {
 691         hbexist[OFFSET(frame)] |= MASK(frame);
 692     }
 693     else
 694     {
 695         hbexist[OFFSET(frame)] &= ~MASK(frame);
 696     }
 697 }
 698
 699 static gmx_bool is_hb(unsigned int hbexist[], int frame)
 700 {
 701     return ((hbexist[OFFSET(frame)] & MASK(frame)) != 0) ? 1 : 0;
 702 }
 703
 704 static void set_hb(t_hbdata *hb, int id, int ih, int ia, int frame, int ihb)
 705 {
 706     unsigned int *ghptr = NULL;
 707
 708     if (ihb == hbHB)
 709     {
 710         ghptr = hb->hbmap[id][ia]->h[ih];
 711     }
 712     else if (ihb == hbDist)
 713     {
 714         ghptr = hb->hbmap[id][ia]->g[ih];
 715     }
 716     else
 717     {
 718         gmx_fatal(FARGS, "Incomprehensible iValue %d in set_hb", ihb);
 719     }
 720
 721     _set_hb(ghptr, frame-hb->hbmap[id][ia]->n0, TRUE);
 722 }
 723
 724 static void addPshift(t_pShift *pHist, PSTYPE p, int frame)
 725 {
 726     if (pHist->len == 0)
 727     {
 728         snew(pHist->frame, 1);
 729         snew(pHist->p, 1);
 730         pHist->len      = 1;
 731         pHist->frame[0] = frame;
 732         pHist->p[0]     = p;
 733         return;
 734     }
 735     else
 736     if (pHist->p[pHist->len-1] != p)
 737     {
 738         pHist->len++;
 739         srenew(pHist->frame, pHist->len);
 740         srenew(pHist->p, pHist->len);
 741         pHist->frame[pHist->len-1] = frame;
 742         pHist->p[pHist->len-1]     = p;
 743     }     /* Otherwise, there is no transition. */
 744     return;
 745 }
 746
 747 static PSTYPE getPshift(t_pShift pHist, int frame)
 748 {
 749     int f, i;
 750
 751     if (pHist.len == 0
 752         || (pHist.len > 0 && pHist.frame[0] > frame))
 753     {
 754         return -1;
 755     }
 756
 757     for (i = 0; i < pHist.len; i++)
 758     {
 759         f = pHist.frame[i];
 760         if (f == frame)
 761         {
 762             return pHist.p[i];
 763         }
 764         if (f > frame)
 765         {
 766             return pHist.p[i-1];
 767         }
 768     }
 769
 770     /* It seems that frame is after the last periodic transition. Return the last periodicity. */
 771     return pHist.p[pHist.len-1];
 772 }
 773
 774 static void add_ff(t_hbdata *hbd, int id, int h, int ia, int frame, int ihb, PSTYPE p)
 775 {
 776     int         i, j, n;
 777     t_hbond    *hb       = hbd->hbmap[id][ia];
 778     int         maxhydro = min(hbd->maxhydro, hbd->d.nhydro[id]);
 779     int         wlen     = hbd->wordlen;
 780     int         delta    = 32*wlen;
 781     gmx_bool    bGem     = hbd->bGem;
 782
 783     if (!hb->h[0])
 784     {
 785         hb->n0        = frame;
 786         hb->maxframes = delta;
 787         for (i = 0; (i < maxhydro); i++)
 788         {
 789             snew(hb->h[i], hb->maxframes/wlen);
 790             snew(hb->g[i], hb->maxframes/wlen);
 791         }
 792     }
 793     else
 794     {
 795         hb->nframes = frame-hb->n0;
 796         /* We need a while loop here because hbonds may be returning
 797          * after a long time.
 798          */
 799         while (hb->nframes >= hb->maxframes)
 800         {
 801             n = hb->maxframes + delta;
 802             for (i = 0; (i < maxhydro); i++)
 803             {
 804                 srenew(hb->h[i], n/wlen);
 805                 srenew(hb->g[i], n/wlen);
 806                 for (j = hb->maxframes/wlen; (j < n/wlen); j++)
 807                 {
 808                     hb->h[i][j] = 0;
 809                     hb->g[i][j] = 0;
 810                 }
 811             }
 812
 813             hb->maxframes = n;
 814         }
 815     }
 816     if (frame >= 0)
 817     {
 818         set_hb(hbd, id, h, ia, frame, ihb);
 819         if (bGem)
 820         {
 821             if (p >= hbd->per->nper)
 822             {
 823                 gmx_fatal(FARGS, "invalid shift: p=%u, nper=%u", p, hbd->per->nper);
 824             }
 825             else
 826             {
 827                 addPshift(&(hbd->per->pHist[id][ia]), p, frame);
 828             }
 829
 830         }
 831     }
 832
 833 }
 834
 835 static void inc_nhbonds(t_donors *ddd, int d, int h)
 836 {
 837     int j;
 838     int dptr = ddd->dptr[d];
 839
 840     for (j = 0; (j < ddd->nhydro[dptr]); j++)
 841     {
 842         if (ddd->hydro[dptr][j] == h)
 843         {
 844             ddd->nhbonds[dptr][j]++;
 845             break;
 846         }
 847     }
 848     if (j == ddd->nhydro[dptr])
 849     {
 850         gmx_fatal(FARGS, "No such hydrogen %d on donor %d\n", h+1, d+1);
 851     }
 852 }
 853
 854 static int _acceptor_index(t_acceptors *a, int grp, atom_id i,
 855                            const char *file, int line)
 856 {
 857     int ai = a->aptr[i];
 858
 859     if (a->grp[ai] != grp)
 860     {
 861         if (debug && bDebug)
 862         {
 863             fprintf(debug, "Acc. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
 864                     ai, a->grp[ai], grp, file, line);
 865         }
 866         return NOTSET;
 867     }
 868     else
 869     {
 870         return ai;
 871     }
 872 }
 873 #define acceptor_index(a, grp, i) _acceptor_index(a, grp, i, __FILE__, __LINE__)
 874
 875 static int _donor_index(t_donors *d, int grp, atom_id i, const char *file, int line)
 876 {
 877     int di = d->dptr[i];
 878
 879     if (di == NOTSET)
 880     {
 881         return NOTSET;
 882     }
 883
 884     if (d->grp[di] != grp)
 885     {
 886         if (debug && bDebug)
 887         {
 888             fprintf(debug, "Don. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
 889                     di, d->grp[di], grp, file, line);
 890         }
 891         return NOTSET;
 892     }
 893     else
 894     {
 895         return di;
 896     }
 897 }
 898 #define donor_index(d, grp, i) _donor_index(d, grp, i, __FILE__, __LINE__)
 899
 900 static gmx_bool isInterchangable(t_hbdata *hb, int d, int a, int grpa, int grpd)
 901 {
 902     /* g_hbond doesn't allow overlapping groups */
 903     if (grpa != grpd)
 904     {
 905         return FALSE;
 906     }
 907     return
 908         donor_index(&hb->d, grpd, a) != NOTSET
 909         && acceptor_index(&hb->a, grpa, d) != NOTSET;
 910 }
 911
 912
 913 static void add_hbond(t_hbdata *hb, int d, int a, int h, int grpd, int grpa,
 914                       int frame, gmx_bool bMerge, int ihb, gmx_bool bContact, PSTYPE p)
 915 {
 916     int      k, id, ia, hh;
 917     gmx_bool daSwap = FALSE;
 918
 919     if ((id = hb->d.dptr[d]) == NOTSET)
 920     {
 921         gmx_fatal(FARGS, "No donor atom %d", d+1);
 922     }
 923     else if (grpd != hb->d.grp[id])
 924     {
 925         gmx_fatal(FARGS, "Inconsistent donor groups, %d iso %d, atom %d",
 926                   grpd, hb->d.grp[id], d+1);
 927     }
 928     if ((ia = hb->a.aptr[a]) == NOTSET)
 929     {
 930         gmx_fatal(FARGS, "No acceptor atom %d", a+1);
 931     }
 932     else if (grpa != hb->a.grp[ia])
 933     {
 934         gmx_fatal(FARGS, "Inconsistent acceptor groups, %d iso %d, atom %d",
 935                   grpa, hb->a.grp[ia], a+1);
 936     }
 937
 938     if (bMerge)
 939     {
 940
 941         if (isInterchangable(hb, d, a, grpd, grpa) && d > a)
 942         /* Then swap identity so that the id of d is lower then that of a.
 943          *
 944          * This should really be redundant by now, as is_hbond() now ought to return
 945          * hbNo in the cases where this conditional is TRUE. */
 946         {
 947             daSwap = TRUE;
 948             k      = d;
 949             d      = a;
 950             a      = k;
 951
 952             /* Now repeat donor/acc check. */
 953             if ((id = hb->d.dptr[d]) == NOTSET)
 954             {
 955                 gmx_fatal(FARGS, "No donor atom %d", d+1);
 956             }
 957             else if (grpd != hb->d.grp[id])
 958             {
 959                 gmx_fatal(FARGS, "Inconsistent donor groups, %d iso %d, atom %d",
 960                           grpd, hb->d.grp[id], d+1);
 961             }
 962             if ((ia = hb->a.aptr[a]) == NOTSET)
 963             {
 964                 gmx_fatal(FARGS, "No acceptor atom %d", a+1);
 965             }
 966             else if (grpa != hb->a.grp[ia])
 967             {
 968                 gmx_fatal(FARGS, "Inconsistent acceptor groups, %d iso %d, atom %d",
 969                           grpa, hb->a.grp[ia], a+1);
 970             }
 971         }
 972     }
 973
 974     if (hb->hbmap)
 975     {
 976         /* Loop over hydrogens to find which hydrogen is in this particular HB */
 977         if ((ihb == hbHB) && !bMerge && !bContact)
 978         {
 979             for (k = 0; (k < hb->d.nhydro[id]); k++)
 980             {
 981                 if (hb->d.hydro[id][k] == h)
 982                 {
 983                     break;
 984                 }
 985             }
 986             if (k == hb->d.nhydro[id])
 987             {
 988                 gmx_fatal(FARGS, "Donor %d does not have hydrogen %d (a = %d)",
 989                           d+1, h+1, a+1);
 990             }
 991         }
 992         else
 993         {
 994             k = 0;
 995         }
 996
 997         if (hb->bHBmap)
 998         {
 999
1000 #pragma omp critical
1001             {
1002                 if (hb->hbmap[id][ia] == NULL)
1003                 {
1004                     snew(hb->hbmap[id][ia], 1);
1005                     snew(hb->hbmap[id][ia]->h, hb->maxhydro);
1006                     snew(hb->hbmap[id][ia]->g, hb->maxhydro);
1007                 }
1008                 add_ff(hb, id, k, ia, frame, ihb, p);
1009             }
1010         }
1011
1012         /* Strange construction with frame >=0 is a relic from old code
1013          * for selected hbond analysis. It may be necessary again if that
1014          * is made to work again.
1015          */
1016         if (frame >= 0)
1017         {
1018             hh = hb->hbmap[id][ia]->history[k];
1019             if (ihb == hbHB)
1020             {
1021                 hb->nhb[frame]++;
1022                 if (!(ISHB(hh)))
1023                 {
1024                     hb->hbmap[id][ia]->history[k] = hh | 2;
1025                     hb->nrhb++;
1026                 }
1027             }
1028             else
1029             {
1030                 if (ihb == hbDist)
1031                 {
1032                     hb->ndist[frame]++;
1033                     if (!(ISDIST(hh)))
1034                     {
1035                         hb->hbmap[id][ia]->history[k] = hh | 1;
1036                         hb->nrdist++;
1037                     }
1038                 }
1039             }
1040         }
1041     }
1042     else
1043     {
1044         if (frame >= 0)
1045         {
1046             if (ihb == hbHB)
1047             {
1048                 hb->nhb[frame]++;
1049             }
1050             else
1051             {
1052                 if (ihb == hbDist)
1053                 {
1054                     hb->ndist[frame]++;
1055                 }
1056             }
1057         }
1058     }
1059     if (bMerge && daSwap)
1060     {
1061         h = hb->d.hydro[id][0];
1062     }
1063     /* Increment number if HBonds per H */
1064     if (ihb == hbHB && !bContact)
1065     {
1066         inc_nhbonds(&(hb->d), d, h);
1067     }
1068 }
1069
1070 static char *mkatomname(t_atoms *atoms, int i)
1071 {
1072     static char buf[32];
1073     int         rnr;
1074
1075     rnr = atoms->atom[i].resind;
1076     sprintf(buf, "%4s%d%-4s",
1077             *atoms->resinfo[rnr].name, atoms->resinfo[rnr].nr, *atoms->atomname[i]);
1078
1079     return buf;
1080 }
1081
1082 static void gen_datable(atom_id *index, int isize, unsigned char *datable, int natoms)
1083 {
1084     /* Generates table of all atoms and sets the ingroup bit for atoms in index[] */
1085     int i;
1086
1087     for (i = 0; i < isize; i++)
1088     {
1089         if (index[i] >= natoms)
1090         {
1091             gmx_fatal(FARGS, "Atom has index %d larger than number of atoms %d.", index[i], natoms);
1092         }
1093         datable[index[i]] |= INGROUP;
1094     }
1095 }
1096
1097 static void clear_datable_grp(unsigned char *datable, int size)
1098 {
1099     /* Clears group information from the table */
1100     int        i;
1101     const char mask = !(char)INGROUP;
1102     if (size > 0)
1103     {
1104         for (i = 0; i < size; i++)
1105         {
1106             datable[i] &= mask;
1107         }
1108     }
1109 }
1110
1111 static void add_acc(t_acceptors *a, int ia, int grp)
1112 {
1113     if (a->nra >= a->max_nra)
1114     {
1115         a->max_nra += 16;
1116         srenew(a->acc, a->max_nra);
1117         srenew(a->grp, a->max_nra);
1118     }
1119     a->grp[a->nra]   = grp;
1120     a->acc[a->nra++] = ia;
1121 }
1122
1123 static void search_acceptors(t_topology *top, int isize,
1124                              atom_id *index, t_acceptors *a, int grp,
1125                              gmx_bool bNitAcc,
1126                              gmx_bool bContact, gmx_bool bDoIt, unsigned char *datable)
1127 {
1128     int i, n;
1129
1130     if (bDoIt)
1131     {
1132         for (i = 0; (i < isize); i++)
1133         {
1134             n = index[i];
1135             if ((bContact ||
1136                  (((*top->atoms.atomname[n])[0] == 'O') ||
1137                   (bNitAcc && ((*top->atoms.atomname[n])[0] == 'N')))) &&
1138                 ISINGRP(datable[n]))
1139             {
1140                 datable[n] |= ACC; /* set the atom's acceptor flag in datable. */
1141                 add_acc(a, n, grp);
1142             }
1143         }
1144     }
1145     snew(a->aptr, top->atoms.nr);
1146     for (i = 0; (i < top->atoms.nr); i++)
1147     {
1148         a->aptr[i] = NOTSET;
1149     }
1150     for (i = 0; (i < a->nra); i++)
1151     {
1152         a->aptr[a->acc[i]] = i;
1153     }
1154 }
1155
1156 static void add_h2d(int id, int ih, t_donors *ddd)
1157 {
1158     int i;
1159
1160     for (i = 0; (i < ddd->nhydro[id]); i++)
1161     {
1162         if (ddd->hydro[id][i] == ih)
1163         {
1164             printf("Hm. This isn't the first time I found this donor (%d,%d)\n",
1165                    ddd->don[id], ih);
1166             break;
1167         }
1168     }
1169     if (i == ddd->nhydro[id])
1170     {
1171         if (ddd->nhydro[id] >= MAXHYDRO)
1172         {
1173             gmx_fatal(FARGS, "Donor %d has more than %d hydrogens!",
1174                       ddd->don[id], MAXHYDRO);
1175         }
1176         ddd->hydro[id][i] = ih;
1177         ddd->nhydro[id]++;
1178     }
1179 }
1180
1181 static void add_dh(t_donors *ddd, int id, int ih, int grp, unsigned char *datable)
1182 {
1183     int i;
1184
1185     if (ISDON(datable[id]) || !datable)
1186     {
1187         if (ddd->dptr[id] == NOTSET)   /* New donor */
1188         {
1189             i             = ddd->nrd;
1190             ddd->dptr[id] = i;
1191         }
1192         else
1193         {
1194             i = ddd->dptr[id];
1195         }
1196
1197         if (i == ddd->nrd)
1198         {
1199             if (ddd->nrd >= ddd->max_nrd)
1200             {
1201                 ddd->max_nrd += 128;
1202                 srenew(ddd->don, ddd->max_nrd);
1203                 srenew(ddd->nhydro, ddd->max_nrd);
1204                 srenew(ddd->hydro, ddd->max_nrd);
1205                 srenew(ddd->nhbonds, ddd->max_nrd);
1206                 srenew(ddd->grp, ddd->max_nrd);
1207             }
1208             ddd->don[ddd->nrd]    = id;
1209             ddd->nhydro[ddd->nrd] = 0;
1210             ddd->grp[ddd->nrd]    = grp;
1211             ddd->nrd++;
1212         }
1213         else
1214         {
1215             ddd->don[i] = id;
1216         }
1217         add_h2d(i, ih, ddd);
1218     }
1219     else
1220     if (datable)
1221     {
1222         printf("Warning: Atom %d is not in the d/a-table!\n", id);
1223     }
1224 }
1225
1226 static void search_donors(t_topology *top, int isize, atom_id *index,
1227                           t_donors *ddd, int grp, gmx_bool bContact, gmx_bool bDoIt,
1228                           unsigned char *datable)
1229 {
1230     int            i, j, nra, n;
1231     t_functype     func_type;
1232     t_ilist       *interaction;
1233     atom_id        nr1, nr2, nr3;
1234     gmx_bool       stop;
1235
1236     if (!ddd->dptr)
1237     {
1238         snew(ddd->dptr, top->atoms.nr);
1239         for (i = 0; (i < top->atoms.nr); i++)
1240         {
1241             ddd->dptr[i] = NOTSET;
1242         }
1243     }
1244
1245     if (bContact)
1246     {
1247         if (bDoIt)
1248         {
1249             for (i = 0; (i < isize); i++)
1250             {
1251                 datable[index[i]] |= DON;
1252                 add_dh(ddd, index[i], -1, grp, datable);
1253             }
1254         }
1255     }
1256     else
1257     {
1258         for (func_type = 0; (func_type < F_NRE); func_type++)
1259         {
1260             interaction = &(top->idef.il[func_type]);
1261             if (func_type == F_POSRES || func_type == F_FBPOSRES)
1262             {
1263                 /* The ilist looks strange for posre. Bug in grompp?
1264                  * We don't need posre interactions for hbonds anyway.*/
1265                 continue;
1266             }
1267             for (i = 0; i < interaction->nr;
1268                  i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1)
1269             {
1270                 /* next function */
1271                 if (func_type != top->idef.functype[interaction->iatoms[i]])
1272                 {
1273                     fprintf(stderr, "Error in func_type %s",
1274                             interaction_function[func_type].longname);
1275                     continue;
1276                 }
1277
1278                 /* check out this functype */
1279                 if (func_type == F_SETTLE)
1280                 {
1281                     nr1 = interaction->iatoms[i+1];
1282                     nr2 = interaction->iatoms[i+2];
1283                     nr3 = interaction->iatoms[i+3];
1284
1285                     if (ISINGRP(datable[nr1]))
1286                     {
1287                         if (ISINGRP(datable[nr2]))
1288                         {
1289                             datable[nr1] |= DON;
1290                             add_dh(ddd, nr1, nr1+1, grp, datable);
1291                         }
1292                         if (ISINGRP(datable[nr3]))
1293                         {
1294                             datable[nr1] |= DON;
1295                             add_dh(ddd, nr1, nr1+2, grp, datable);
1296                         }
1297                     }
1298                 }
1299                 else if (IS_CHEMBOND(func_type))
1300                 {
1301                     for (j = 0; j < 2; j++)
1302                     {
1303                         nr1 = interaction->iatoms[i+1+j];
1304                         nr2 = interaction->iatoms[i+2-j];
1305                         if ((*top->atoms.atomname[nr1][0] == 'H') &&
1306                             ((*top->atoms.atomname[nr2][0] == 'O') ||
1307                              (*top->atoms.atomname[nr2][0] == 'N')) &&
1308                             ISINGRP(datable[nr1]) && ISINGRP(datable[nr2]))
1309                         {
1310                             datable[nr2] |= DON;
1311                             add_dh(ddd, nr2, nr1, grp, datable);
1312                         }
1313                     }
1314                 }
1315             }
1316         }
1317 #ifdef SAFEVSITES
1318         for (func_type = 0; func_type < F_NRE; func_type++)
1319         {
1320             interaction = &top->idef.il[func_type];
1321             for (i = 0; i < interaction->nr;
1322                  i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1)
1323             {
1324                 /* next function */
1325                 if (func_type != top->idef.functype[interaction->iatoms[i]])
1326                 {
1327                     gmx_incons("function type in search_donors");
1328                 }
1329
1330                 if (interaction_function[func_type].flags & IF_VSITE)
1331                 {
1332                     nr1 = interaction->iatoms[i+1];
1333                     if (*top->atoms.atomname[nr1][0]  == 'H')
1334                     {
1335                         nr2  = nr1-1;
1336                         stop = FALSE;
1337                         while (!stop && ( *top->atoms.atomname[nr2][0] == 'H'))
1338                         {
1339                             if (nr2)
1340                             {
1341                                 nr2--;
1342                             }
1343                             else
1344                             {
1345                                 stop = TRUE;
1346                             }
1347                         }
1348                         if (!stop && ( ( *top->atoms.atomname[nr2][0] == 'O') ||
1349                                        ( *top->atoms.atomname[nr2][0] == 'N') ) &&
1350                             ISINGRP(datable[nr1]) && ISINGRP(datable[nr2]))
1351                         {
1352                             datable[nr2] |= DON;
1353                             add_dh(ddd, nr2, nr1, grp, datable);
1354                         }
1355                     }
1356                 }
1357             }
1358         }
1359 #endif
1360     }
1361 }
1362
1363 static t_gridcell ***init_grid(gmx_bool bBox, rvec box[], real rcut, ivec ngrid)
1364 {
1365     t_gridcell ***grid;
1366     int           i, y, z;
1367
1368     if (bBox)
1369     {
1370         for (i = 0; i < DIM; i++)
1371         {
1372             ngrid[i] = (box[i][i]/(1.2*rcut));
1373         }
1374     }
1375
1376     if (!bBox || (ngrid[XX] < 3) || (ngrid[YY] < 3) || (ngrid[ZZ] < 3) )
1377     {
1378         for (i = 0; i < DIM; i++)
1379         {
1380             ngrid[i] = 1;
1381         }
1382     }
1383     else
1384     {
1385         printf("\nWill do grid-seach on %dx%dx%d grid, rcut=%g\n",
1386                ngrid[XX], ngrid[YY], ngrid[ZZ], rcut);
1387     }
1388     snew(grid, ngrid[ZZ]);
1389     for (z = 0; z < ngrid[ZZ]; z++)
1390     {
1391         snew((grid)[z], ngrid[YY]);
1392         for (y = 0; y < ngrid[YY]; y++)
1393         {
1394             snew((grid)[z][y], ngrid[XX]);
1395         }
1396     }
1397     return grid;
1398 }
1399
1400 static void reset_nhbonds(t_donors *ddd)
1401 {
1402     int i, j;
1403
1404     for (i = 0; (i < ddd->nrd); i++)
1405     {
1406         for (j = 0; (j < MAXHH); j++)
1407         {
1408             ddd->nhbonds[i][j] = 0;
1409         }
1410     }
1411 }
1412
1413 void pbc_correct_gem(rvec dx, matrix box, rvec hbox);
1414
1415 static void build_grid(t_hbdata *hb, rvec x[], rvec xshell,
1416                        gmx_bool bBox, matrix box, rvec hbox,
1417                        real rcut, real rshell,
1418                        ivec ngrid, t_gridcell ***grid)
1419 {
1420     int         i, m, gr, xi, yi, zi, nr;
1421     atom_id    *ad;
1422     ivec        grididx;
1423     rvec        invdelta, dshell, xtemp = {0, 0, 0};
1424     t_ncell    *newgrid;
1425     gmx_bool    bDoRshell, bInShell, bAcc;
1426     real        rshell2 = 0;
1427     int         gx, gy, gz;
1428     int         dum = -1;
1429
1430     bDoRshell = (rshell > 0);
1431     rshell2   = sqr(rshell);
1432     bInShell  = TRUE;
1433
1434 #define DBB(x) if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n", __LINE__,#x, x)
1435     DBB(dum);
1436     for (m = 0; m < DIM; m++)
1437     {
1438         hbox[m] = box[m][m]*0.5;
1439         if (bBox)
1440         {
1441             invdelta[m] = ngrid[m]/box[m][m];
1442             if (1/invdelta[m] < rcut)
1443             {
1444                 gmx_fatal(FARGS, "Your computational box has shrunk too much.\n"
1445                           "%s can not handle this situation, sorry.\n",
1446                           ShortProgram());
1447             }
1448         }
1449         else
1450         {
1451             invdelta[m] = 0;
1452         }
1453     }
1454     grididx[XX] = 0;
1455     grididx[YY] = 0;
1456     grididx[ZZ] = 0;
1457     DBB(dum);
1458     /* resetting atom counts */
1459     for (gr = 0; (gr < grNR); gr++)
1460     {
1461         for (zi = 0; zi < ngrid[ZZ]; zi++)
1462         {
1463             for (yi = 0; yi < ngrid[YY]; yi++)
1464             {
1465                 for (xi = 0; xi < ngrid[XX]; xi++)
1466                 {
1467                     grid[zi][yi][xi].d[gr].nr = 0;
1468                     grid[zi][yi][xi].a[gr].nr = 0;
1469                 }
1470             }
1471         }
1472         DBB(dum);
1473
1474         /* put atoms in grid cells */
1475         for (bAcc = FALSE; (bAcc <= TRUE); bAcc++)
1476         {
1477             if (bAcc)
1478             {
1479                 nr = hb->a.nra;
1480                 ad = hb->a.acc;
1481             }
1482             else
1483             {
1484                 nr = hb->d.nrd;
1485                 ad = hb->d.don;
1486             }
1487             DBB(bAcc);
1488             for (i = 0; (i < nr); i++)
1489             {
1490                 /* check if we are inside the shell */
1491                 /* if bDoRshell=FALSE then bInShell=TRUE always */
1492                 DBB(i);
1493                 if (bDoRshell)
1494                 {
1495                     bInShell = TRUE;
1496                     rvec_sub(x[ad[i]], xshell, dshell);
1497                     if (bBox)
1498                     {
1499                         if (FALSE && !hb->bGem)
1500                         {
1501                             for (m = DIM-1; m >= 0 && bInShell; m--)
1502                             {
1503                                 if (dshell[m] < -hbox[m])
1504                                 {
1505                                     rvec_inc(dshell, box[m]);
1506                                 }
1507                                 else if (dshell[m] >= hbox[m])
1508                                 {
1509                                     dshell[m] -= 2*hbox[m];
1510                                 }
1511                                 /* if we're outside the cube, we're outside the sphere also! */
1512                                 if ( (dshell[m] > rshell) || (-dshell[m] > rshell) )
1513                                 {
1514                                     bInShell = FALSE;
1515                                 }
1516                             }
1517                         }
1518                         else
1519                         {
1520                             gmx_bool bDone = FALSE;
1521                             while (!bDone)
1522                             {
1523                                 bDone = TRUE;
1524                                 for (m = DIM-1; m >= 0 && bInShell; m--)
1525                                 {
1526                                     if (dshell[m] < -hbox[m])
1527                                     {
1528                                         bDone = FALSE;
1529                                         rvec_inc(dshell, box[m]);
1530                                     }
1531                                     if (dshell[m] >= hbox[m])
1532                                     {
1533                                         bDone      = FALSE;
1534                                         dshell[m] -= 2*hbox[m];
1535                                     }
1536                                 }
1537                             }
1538                             for (m = DIM-1; m >= 0 && bInShell; m--)
1539                             {
1540                                 /* if we're outside the cube, we're outside the sphere also! */
1541                                 if ( (dshell[m] > rshell) || (-dshell[m] > rshell) )
1542                                 {
1543                                     bInShell = FALSE;
1544                                 }
1545                             }
1546                         }
1547                     }
1548                     /* if we're inside the cube, check if we're inside the sphere */
1549                     if (bInShell)
1550                     {
1551                         bInShell = norm2(dshell) < rshell2;
1552                     }
1553                 }
1554                 DBB(i);
1555                 if (bInShell)
1556                 {
1557                     if (bBox)
1558                     {
1559                         if (hb->bGem)
1560                         {
1561                             copy_rvec(x[ad[i]], xtemp);
1562                         }
1563                         pbc_correct_gem(x[ad[i]], box, hbox);
1564                     }
1565                     for (m = DIM-1; m >= 0; m--)
1566                     {
1567                         if (TRUE || !hb->bGem)
1568                         {
1569                             /* put atom in the box */
1570                             while (x[ad[i]][m] < 0)
1571                             {
1572                                 rvec_inc(x[ad[i]], box[m]);
1573                             }
1574                             while (x[ad[i]][m] >= box[m][m])
1575                             {
1576                                 rvec_dec(x[ad[i]], box[m]);
1577                             }
1578                         }
1579                         /* determine grid index of atom */
1580                         grididx[m] = x[ad[i]][m]*invdelta[m];
1581                         grididx[m] = (grididx[m]+ngrid[m]) % ngrid[m];
1582                     }
1583                     if (hb->bGem)
1584                     {
1585                         copy_rvec(xtemp, x[ad[i]]); /* copy back */
1586                     }
1587                     gx = grididx[XX];
1588                     gy = grididx[YY];
1589                     gz = grididx[ZZ];
1590                     range_check(gx, 0, ngrid[XX]);
1591                     range_check(gy, 0, ngrid[YY]);
1592                     range_check(gz, 0, ngrid[ZZ]);
1593                     DBB(gx);
1594                     DBB(gy);
1595                     DBB(gz);
1596                     /* add atom to grid cell */
1597                     if (bAcc)
1598                     {
1599                         newgrid = &(grid[gz][gy][gx].a[gr]);
1600                     }
1601                     else
1602                     {
1603                         newgrid = &(grid[gz][gy][gx].d[gr]);
1604                     }
1605                     if (newgrid->nr >= newgrid->maxnr)
1606                     {
1607                         newgrid->maxnr += 10;
1608                         DBB(newgrid->maxnr);
1609                         srenew(newgrid->atoms, newgrid->maxnr);
1610                     }
1611                     DBB(newgrid->nr);
1612                     newgrid->atoms[newgrid->nr] = ad[i];
1613                     newgrid->nr++;
1614                 }
1615             }
1616         }
1617     }
1618 }
1619
1620 static void count_da_grid(ivec ngrid, t_gridcell ***grid, t_icell danr)
1621 {
1622     int gr, xi, yi, zi;
1623
1624     for (gr = 0; (gr < grNR); gr++)
1625     {
1626         danr[gr] = 0;
1627         for (zi = 0; zi < ngrid[ZZ]; zi++)
1628         {
1629             for (yi = 0; yi < ngrid[YY]; yi++)
1630             {
1631                 for (xi = 0; xi < ngrid[XX]; xi++)
1632                 {
1633                     danr[gr] += grid[zi][yi][xi].d[gr].nr;
1634                 }
1635             }
1636         }
1637     }
1638 }
1639
1640 /* The grid loop.
1641  * Without a box, the grid is 1x1x1, so all loops are 1 long.
1642  * With a rectangular box (bTric==FALSE) all loops are 3 long.
1643  * With a triclinic box all loops are 3 long, except when a cell is
1644  * located next to one of the box edges which is not parallel to the
1645  * x/y-plane, in that case all cells in a line or layer are searched.
1646  * This could be implemented slightly more efficient, but the code
1647  * would get much more complicated.
1648  */
1649 static inline gmx_bool grid_loop_begin(int n, int x, gmx_bool bTric, gmx_bool bEdge)
1650 {
1651     return ((n == 1) ? x : bTric && bEdge ? 0     : (x-1));
1652 }
1653 static inline gmx_bool grid_loop_end(int n, int x, gmx_bool bTric, gmx_bool bEdge)
1654 {
1655     return ((n == 1) ? x : bTric && bEdge ? (n-1) : (x+1));
1656 }
1657 static inline int grid_mod(int j, int n)
1658 {
1659     return (j+n) % (n);
1660 }
1661
1662 static void dump_grid(FILE *fp, ivec ngrid, t_gridcell ***grid)
1663 {
1664     int gr, x, y, z, sum[grNR];
1665
1666     fprintf(fp, "grid %dx%dx%d\n", ngrid[XX], ngrid[YY], ngrid[ZZ]);
1667     for (gr = 0; gr < grNR; gr++)
1668     {
1669         sum[gr] = 0;
1670         fprintf(fp, "GROUP %d (%s)\n", gr, grpnames[gr]);
1671         for (z = 0; z < ngrid[ZZ]; z += 2)
1672         {
1673             fprintf(fp, "Z=%d,%d\n", z, z+1);
1674             for (y = 0; y < ngrid[YY]; y++)
1675             {
1676                 for (x = 0; x < ngrid[XX]; x++)
1677                 {
1678                     fprintf(fp, "%3d", grid[x][y][z].d[gr].nr);
1679                     sum[gr] += grid[z][y][x].d[gr].nr;
1680                     fprintf(fp, "%3d", grid[x][y][z].a[gr].nr);
1681                     sum[gr] += grid[z][y][x].a[gr].nr;
1682
1683                 }
1684                 fprintf(fp, " | ");
1685                 if ( (z+1) < ngrid[ZZ])
1686                 {
1687                     for (x = 0; x < ngrid[XX]; x++)
1688                     {
1689                         fprintf(fp, "%3d", grid[z+1][y][x].d[gr].nr);
1690                         sum[gr] += grid[z+1][y][x].d[gr].nr;
1691                         fprintf(fp, "%3d", grid[z+1][y][x].a[gr].nr);
1692                         sum[gr] += grid[z+1][y][x].a[gr].nr;
1693                     }
1694                 }
1695                 fprintf(fp, "\n");
1696             }
1697         }
1698     }
1699     fprintf(fp, "TOTALS:");
1700     for (gr = 0; gr < grNR; gr++)
1701     {
1702         fprintf(fp, " %d=%d", gr, sum[gr]);
1703     }
1704     fprintf(fp, "\n");
1705 }
1706
1707 /* New GMX record! 5 * in a row. Congratulations!
1708  * Sorry, only four left.
1709  */
1710 static void free_grid(ivec ngrid, t_gridcell ****grid)
1711 {
1712     int           y, z;
1713     t_gridcell ***g = *grid;
1714
1715     for (z = 0; z < ngrid[ZZ]; z++)
1716     {
1717         for (y = 0; y < ngrid[YY]; y++)
1718         {
1719             sfree(g[z][y]);
1720         }
1721         sfree(g[z]);
1722     }
1723     sfree(g);
1724     g = NULL;
1725 }
1726
1727 void pbc_correct_gem(rvec dx, matrix box, rvec hbox)
1728 {
1729     int      m;
1730     gmx_bool bDone = FALSE;
1731     while (!bDone)
1732     {
1733         bDone = TRUE;
1734         for (m = DIM-1; m >= 0; m--)
1735         {
1736             if (dx[m] < -hbox[m])
1737             {
1738                 bDone = FALSE;
1739                 rvec_inc(dx, box[m]);
1740             }
1741             if (dx[m] >= hbox[m])
1742             {
1743                 bDone = FALSE;
1744                 rvec_dec(dx, box[m]);
1745             }
1746         }
1747     }
1748 }
1749
1750 /* Added argument r2cut, changed contact and implemented
1751  * use of second cut-off.
1752  * - Erik Marklund, June 29, 2006
1753  */
1754 static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a,
1755                     real rcut, real r2cut, real ccut,
1756                     rvec x[], gmx_bool bBox, matrix box, rvec hbox,
1757                     real *d_ha, real *ang, gmx_bool bDA, int *hhh,
1758                     gmx_bool bContact, gmx_bool bMerge, PSTYPE *p)
1759 {
1760     int      h, hh, id, ja, ihb;
1761     rvec     r_da, r_ha, r_dh, r = {0, 0, 0};
1762     ivec     ri;
1763     real     rc2, r2c2, rda2, rha2, ca;
1764     gmx_bool HAinrange = FALSE; /* If !bDA. Needed for returning hbDist in a correct way. */
1765     gmx_bool daSwap    = FALSE;
1766
1767     if (d == a)
1768     {
1769         return hbNo;
1770     }
1771
1772     if (((id = donor_index(&hb->d, grpd, d)) == NOTSET) ||
1773         ((ja = acceptor_index(&hb->a, grpa, a)) == NOTSET))
1774     {
1775         return hbNo;
1776     }
1777
1778     rc2  = rcut*rcut;
1779     r2c2 = r2cut*r2cut;
1780
1781     rvec_sub(x[d], x[a], r_da);
1782     /* Insert projection code here */
1783
1784     if (bMerge && d > a && isInterchangable(hb, d, a, grpd, grpa))
1785     {
1786         /* Then this hbond/contact will be found again, or it has already been found. */
1787         /*return hbNo;*/
1788     }
1789     if (bBox)
1790     {
1791         if (d > a && bMerge && isInterchangable(hb, d, a, grpd, grpa)) /* acceptor is also a donor and vice versa? */
1792         {                                                              /* return hbNo; */
1793             daSwap = TRUE;                                             /* If so, then their history should be filed with donor and acceptor swapped. */
1794         }
1795         if (hb->bGem)
1796         {
1797             copy_rvec(r_da, r); /* Save this for later */
1798             pbc_correct_gem(r_da, box, hbox);
1799         }
1800         else
1801         {
1802             pbc_correct_gem(r_da, box, hbox);
1803         }
1804     }
1805     rda2 = iprod(r_da, r_da);
1806
1807     if (bContact)
1808     {
1809         if (daSwap && grpa == grpd)
1810         {
1811             return hbNo;
1812         }
1813         if (rda2 <= rc2)
1814         {
1815             if (hb->bGem)
1816             {
1817                 calcBoxDistance(hb->per->P, r, ri);
1818                 *p = periodicIndex(ri, hb->per, daSwap);    /* find (or add) periodicity index. */
1819             }
1820             return hbHB;
1821         }
1822         else if (rda2 < r2c2)
1823         {
1824             return hbDist;
1825         }
1826         else
1827         {
1828             return hbNo;
1829         }
1830     }
1831     *hhh = NOTSET;
1832
1833     if (bDA && (rda2 > rc2))
1834     {
1835         return hbNo;
1836     }
1837
1838     for (h = 0; (h < hb->d.nhydro[id]); h++)
1839     {
1840         hh   = hb->d.hydro[id][h];
1841         rha2 = rc2+1;
1842         if (!bDA)
1843         {
1844             rvec_sub(x[hh], x[a], r_ha);
1845             if (bBox)
1846             {
1847                 pbc_correct_gem(r_ha, box, hbox);
1848             }
1849             rha2 = iprod(r_ha, r_ha);
1850         }
1851
1852         if (hb->bGem)
1853         {
1854             calcBoxDistance(hb->per->P, r, ri);
1855             *p = periodicIndex(ri, hb->per, daSwap);    /* find periodicity index. */
1856         }
1857
1858         if (bDA || (!bDA && (rha2 <= rc2)))
1859         {
1860             rvec_sub(x[d], x[hh], r_dh);
1861             if (bBox)
1862             {
1863                 pbc_correct_gem(r_dh, box, hbox);
1864             }
1865
1866             if (!bDA)
1867             {
1868                 HAinrange = TRUE;
1869             }
1870             ca = cos_angle(r_dh, r_da);
1871             /* if angle is smaller, cos is larger */
1872             if (ca >= ccut)
1873             {
1874                 *hhh  = hh;
1875                 *d_ha = sqrt(bDA ? rda2 : rha2);
1876                 *ang  = acos(ca);
1877                 return hbHB;
1878             }
1879         }
1880     }
1881     if (bDA || (!bDA && HAinrange))
1882     {
1883         return hbDist;
1884     }
1885     else
1886     {
1887         return hbNo;
1888     }
1889 }
1890
1891 /* Fixed previously undiscovered bug in the merge
1892    code, where the last frame of each hbond disappears.
1893    - Erik Marklund, June 1, 2006 */
1894 /* Added the following arguments:
1895  *   ptmp[] - temporary periodicity hisory
1896  *   a1     - identity of first acceptor/donor
1897  *   a2     - identity of second acceptor/donor
1898  * - Erik Marklund, FEB 20 2010 */
1899
1900 /* Merging is now done on the fly, so do_merge is most likely obsolete now.
1901  * Will do some more testing before removing the function entirely.
1902  * - Erik Marklund, MAY 10 2010 */
1903 static void do_merge(t_hbdata *hb, int ntmp,
1904                      unsigned int htmp[], unsigned int gtmp[], PSTYPE ptmp[],
1905                      t_hbond *hb0, t_hbond *hb1, int a1, int a2)
1906 {
1907     /* Here we need to make sure we're treating periodicity in
1908      * the right way for the geminate recombination kinetics. */
1909
1910     int       m, mm, n00, n01, nn0, nnframes;
1911     PSTYPE    pm;
1912     t_pShift *pShift;
1913
1914     /* Decide where to start from when merging */
1915     n00      = hb0->n0;
1916     n01      = hb1->n0;
1917     nn0      = min(n00, n01);
1918     nnframes = max(n00 + hb0->nframes, n01 + hb1->nframes) - nn0;
1919     /* Initiate tmp arrays */
1920     for (m = 0; (m < ntmp); m++)
1921     {
1922         htmp[m] = 0;
1923         gtmp[m] = 0;
1924         ptmp[m] = 0;
1925     }
1926     /* Fill tmp arrays with values due to first HB */
1927     /* Once again '<' had to be replaced with '<='
1928        to catch the last frame in which the hbond
1929        appears.
1930        - Erik Marklund, June 1, 2006 */
1931     for (m = 0; (m <= hb0->nframes); m++)
1932     {
1933         mm       = m+n00-nn0;
1934         htmp[mm] = is_hb(hb0->h[0], m);
1935         if (hb->bGem)
1936         {
1937             pm = getPshift(hb->per->pHist[a1][a2], m+hb0->n0);
1938             if (pm > hb->per->nper)
1939             {
1940                 gmx_fatal(FARGS, "Illegal shift!");
1941             }
1942             else
1943             {
1944                 ptmp[mm] = pm; /*hb->per->pHist[a1][a2][m];*/
1945             }
1946         }
1947     }
1948     /* If we're doing geminate recompbination we usually don't need the distances.
1949      * Let's save some memory and time. */
1950     if (TRUE || !hb->bGem || hb->per->gemtype == gemAD)
1951     {
1952         for (m = 0; (m <= hb0->nframes); m++)
1953         {
1954             mm       = m+n00-nn0;
1955             gtmp[mm] = is_hb(hb0->g[0], m);
1956         }
1957     }
1958     /* Next HB */
1959     for (m = 0; (m <= hb1->nframes); m++)
1960     {
1961         mm       = m+n01-nn0;
1962         htmp[mm] = htmp[mm] || is_hb(hb1->h[0], m);
1963         gtmp[mm] = gtmp[mm] || is_hb(hb1->g[0], m);
1964         if (hb->bGem /* && ptmp[mm] != 0 */)
1965         {
1966
1967             /* If this hbond has been seen before with donor and acceptor swapped,
1968              * then we need to find the mirrored (*-1) periodicity vector to truely
1969              * merge the hbond history. */
1970             pm = findMirror(getPshift(hb->per->pHist[a2][a1], m+hb1->n0), hb->per->p2i, hb->per->nper);
1971             /* Store index of mirror */
1972             if (pm > hb->per->nper)
1973             {
1974                 gmx_fatal(FARGS, "Illegal shift!");
1975             }
1976             ptmp[mm] = pm;
1977         }
1978     }
1979     /* Reallocate target array */
1980     if (nnframes > hb0->maxframes)
1981     {
1982         srenew(hb0->h[0], 4+nnframes/hb->wordlen);
1983         srenew(hb0->g[0], 4+nnframes/hb->wordlen);
1984     }
1985     if (NULL != hb->per->pHist)
1986     {
1987         clearPshift(&(hb->per->pHist[a1][a2]));
1988     }
1989
1990     /* Copy temp array to target array */
1991     for (m = 0; (m <= nnframes); m++)
1992     {
1993         _set_hb(hb0->h[0], m, htmp[m]);
1994         _set_hb(hb0->g[0], m, gtmp[m]);
1995         if (hb->bGem)
1996         {
1997             addPshift(&(hb->per->pHist[a1][a2]), ptmp[m], m+nn0);
1998         }
1999     }
2000
2001     /* Set scalar variables */
2002     hb0->n0        = nn0;
2003     hb0->maxframes = nnframes;
2004 }
2005
2006 /* Added argument bContact for nicer output.
2007  * Erik Marklund, June 29, 2006
2008  */
2009 static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact)
2010 {
2011     int           i, inrnew, indnew, j, ii, jj, m, id, ia, grp, ogrp, ntmp;
2012     unsigned int *htmp, *gtmp;
2013     PSTYPE       *ptmp;
2014     t_hbond      *hb0, *hb1;
2015
2016     inrnew = hb->nrhb;
2017     indnew = hb->nrdist;
2018
2019     /* Check whether donors are also acceptors */
2020     printf("Merging hbonds with Acceptor and Donor swapped\n");
2021
2022     ntmp = 2*hb->max_frames;
2023     snew(gtmp, ntmp);
2024     snew(htmp, ntmp);
2025     snew(ptmp, ntmp);
2026     for (i = 0; (i < hb->d.nrd); i++)
2027     {
2028         fprintf(stderr, "\r%d/%d", i+1, hb->d.nrd);
2029         id = hb->d.don[i];
2030         ii = hb->a.aptr[id];
2031         for (j = 0; (j < hb->a.nra); j++)
2032         {
2033             ia = hb->a.acc[j];
2034             jj = hb->d.dptr[ia];
2035             if ((id != ia) && (ii != NOTSET) && (jj != NOTSET) &&
2036                 (!bTwo || (bTwo && (hb->d.grp[i] != hb->a.grp[j]))))
2037             {
2038                 hb0 = hb->hbmap[i][j];
2039                 hb1 = hb->hbmap[jj][ii];
2040                 if (hb0 && hb1 && ISHB(hb0->history[0]) && ISHB(hb1->history[0]))
2041                 {
2042                     do_merge(hb, ntmp, htmp, gtmp, ptmp, hb0, hb1, i, j);
2043                     if (ISHB(hb1->history[0]))
2044                     {
2045                         inrnew--;
2046                     }
2047                     else if (ISDIST(hb1->history[0]))
2048                     {
2049                         indnew--;
2050                     }
2051                     else
2052                     if (bContact)
2053                     {
2054                         gmx_incons("No contact history");
2055                     }
2056                     else
2057                     {
2058                         gmx_incons("Neither hydrogen bond nor distance");
2059                     }
2060                     sfree(hb1->h[0]);
2061                     sfree(hb1->g[0]);
2062                     if (hb->bGem)
2063                     {
2064                         clearPshift(&(hb->per->pHist[jj][ii]));
2065                     }
2066                     hb1->h[0]       = NULL;
2067                     hb1->g[0]       = NULL;
2068                     hb1->history[0] = hbNo;
2069                 }
2070             }
2071         }
2072     }
2073     fprintf(stderr, "\n");
2074     printf("- Reduced number of hbonds from %d to %d\n", hb->nrhb, inrnew);
2075     printf("- Reduced number of distances from %d to %d\n", hb->nrdist, indnew);
2076     hb->nrhb   = inrnew;
2077     hb->nrdist = indnew;
2078     sfree(gtmp);
2079     sfree(htmp);
2080     sfree(ptmp);
2081 }
2082
2083 static void do_nhb_dist(FILE *fp, t_hbdata *hb, real t)
2084 {
2085     int  i, j, k, n_bound[MAXHH], nbtot;
2086     h_id nhb;
2087
2088
2089     /* Set array to 0 */
2090     for (k = 0; (k < MAXHH); k++)
2091     {
2092         n_bound[k] = 0;
2093     }
2094     /* Loop over possible donors */
2095     for (i = 0; (i < hb->d.nrd); i++)
2096     {
2097         for (j = 0; (j < hb->d.nhydro[i]); j++)
2098         {
2099             n_bound[hb->d.nhbonds[i][j]]++;
2100         }
2101     }
2102     fprintf(fp, "%12.5e", t);
2103     nbtot = 0;
2104     for (k = 0; (k < MAXHH); k++)
2105     {
2106         fprintf(fp, "  %8d", n_bound[k]);
2107         nbtot += n_bound[k]*k;
2108     }
2109     fprintf(fp, "  %8d\n", nbtot);
2110 }
2111
2112 /* Added argument bContact in do_hblife(...). Also
2113  * added support for -contact in function body.
2114  * - Erik Marklund, May 31, 2006 */
2115 /* Changed the contact code slightly.
2116  * - Erik Marklund, June 29, 2006
2117  */
2118 static void do_hblife(const char *fn, t_hbdata *hb, gmx_bool bMerge, gmx_bool bContact,
2119                       const output_env_t oenv)
2120 {
2121     FILE          *fp;
2122     const char    *leg[] = { "p(t)", "t p(t)" };
2123     int           *histo;
2124     int            i, j, j0, k, m, nh, ihb, ohb, nhydro, ndump = 0;
2125     int            nframes = hb->nframes;
2126     unsigned int **h;
2127     real           t, x1, dt;
2128     double         sum, integral;
2129     t_hbond       *hbh;
2130
2131     snew(h, hb->maxhydro);
2132     snew(histo, nframes+1);
2133     /* Total number of hbonds analyzed here */
2134     for (i = 0; (i < hb->d.nrd); i++)
2135     {
2136         for (k = 0; (k < hb->a.nra); k++)
2137         {
2138             hbh = hb->hbmap[i][k];
2139             if (hbh)
2140             {
2141                 if (bMerge)
2142                 {
2143                     if (hbh->h[0])
2144                     {
2145                         h[0]   = hbh->h[0];
2146                         nhydro = 1;
2147                     }
2148                     else
2149                     {
2150                         nhydro = 0;
2151                     }
2152                 }
2153                 else
2154                 {
2155                     nhydro = 0;
2156                     for (m = 0; (m < hb->maxhydro); m++)
2157                     {
2158                         if (hbh->h[m])
2159                         {
2160                             h[nhydro++] = bContact ? hbh->g[m] : hbh->h[m];
2161                         }
2162                     }
2163                 }
2164                 for (nh = 0; (nh < nhydro); nh++)
2165                 {
2166                     ohb = 0;
2167                     j0  = 0;
2168
2169                     /* Changed '<' into '<=' below, just like I
2170                        did in the hbm-output-loop in the main code.
2171                        - Erik Marklund, May 31, 2006
2172                      */
2173                     for (j = 0; (j <= hbh->nframes); j++)
2174                     {
2175                         ihb      = is_hb(h[nh], j);
2176                         if (debug && (ndump < 10))
2177                         {
2178                             fprintf(debug, "%5d  %5d\n", j, ihb);
2179                         }
2180                         if (ihb != ohb)
2181                         {
2182                             if (ihb)
2183                             {
2184                                 j0 = j;
2185                             }
2186                             else
2187                             {
2188                                 histo[j-j0]++;
2189                             }
2190                             ohb = ihb;
2191                         }
2192                     }
2193                     ndump++;
2194                 }
2195             }
2196         }
2197     }
2198     fprintf(stderr, "\n");
2199     if (bContact)
2200     {
2201         fp = xvgropen(fn, "Uninterrupted contact lifetime", output_env_get_xvgr_tlabel(oenv), "()", oenv);
2202     }
2203     else
2204     {
2205         fp = xvgropen(fn, "Uninterrupted hydrogen bond lifetime", output_env_get_xvgr_tlabel(oenv), "()",
2206                       oenv);
2207     }
2208
2209     xvgr_legend(fp, asize(leg), leg, oenv);
2210     j0 = nframes-1;
2211     while ((j0 > 0) && (histo[j0] == 0))
2212     {
2213         j0--;
2214     }
2215     sum = 0;
2216     for (i = 0; (i <= j0); i++)
2217     {
2218         sum += histo[i];
2219     }
2220     dt       = hb->time[1]-hb->time[0];
2221     sum      = dt*sum;
2222     integral = 0;
2223     for (i = 1; (i <= j0); i++)
2224     {
2225         t  = hb->time[i] - hb->time[0] - 0.5*dt;
2226         x1 = t*histo[i]/sum;
2227         fprintf(fp, "%8.3f  %10.3e  %10.3e\n", t, histo[i]/sum, x1);
2228         integral += x1;
2229     }
2230     integral *= dt;
2231     ffclose(fp);
2232     printf("%s lifetime = %.2f ps\n", bContact ? "Contact" : "HB", integral);
2233     printf("Note that the lifetime obtained in this manner is close to useless\n");
2234     printf("Use the -ac option instead and check the Forward lifetime\n");
2235     please_cite(stdout, "Spoel2006b");
2236     sfree(h);
2237     sfree(histo);
2238 }
2239
2240 /* Changed argument bMerge into oneHB to handle contacts properly.
2241  * - Erik Marklund, June 29, 2006
2242  */
2243 static void dump_ac(t_hbdata *hb, gmx_bool oneHB, int nDump)
2244 {
2245     FILE     *fp;
2246     int       i, j, k, m, nd, ihb, idist;
2247     int       nframes = hb->nframes;
2248     gmx_bool  bPrint;
2249     t_hbond  *hbh;
2250
2251     if (nDump <= 0)
2252     {
2253         return;
2254     }
2255     fp = ffopen("debug-ac.xvg", "w");
2256     for (j = 0; (j < nframes); j++)
2257     {
2258         fprintf(fp, "%10.3f", hb->time[j]);
2259         for (i = nd = 0; (i < hb->d.nrd) && (nd < nDump); i++)
2260         {
2261             for (k = 0; (k < hb->a.nra) && (nd < nDump); k++)
2262             {
2263                 bPrint = FALSE;
2264                 ihb    = idist = 0;
2265                 hbh    = hb->hbmap[i][k];
2266                 if (oneHB)
2267                 {
2268                     if (hbh->h[0])
2269                     {
2270                         ihb    = is_hb(hbh->h[0], j);
2271                         idist  = is_hb(hbh->g[0], j);
2272                         bPrint = TRUE;
2273                     }
2274                 }
2275                 else
2276                 {
2277                     for (m = 0; (m < hb->maxhydro) && !ihb; m++)
2278                     {
2279                         ihb   = ihb   || ((hbh->h[m]) && is_hb(hbh->h[m], j));
2280                         idist = idist || ((hbh->g[m]) && is_hb(hbh->g[m], j));
2281                     }
2282                     /* This is not correct! */
2283                     /* What isn't correct? -Erik M */
2284                     bPrint = TRUE;
2285                 }
2286                 if (bPrint)
2287                 {
2288                     fprintf(fp, "  %1d-%1d", ihb, idist);
2289                     nd++;
2290                 }
2291             }
2292         }
2293         fprintf(fp, "\n");
2294     }
2295     ffclose(fp);
2296 }
2297
2298 static real calc_dg(real tau, real temp)
2299 {
2300     real kbt;
2301
2302     kbt = BOLTZ*temp;
2303     if (tau <= 0)
2304     {
2305         return -666;
2306     }
2307     else
2308     {
2309         return kbt*log(kbt*tau/PLANCK);
2310     }
2311 }
2312
2313 typedef struct {
2314     int   n0, n1, nparams, ndelta;
2315     real  kkk[2];
2316     real *t, *ct, *nt, *kt, *sigma_ct, *sigma_nt, *sigma_kt;
2317 } t_luzar;
2318
2319 #ifdef HAVE_LIBGSL
2320 #include <gsl/gsl_multimin.h>
2321 #include <gsl/gsl_sf.h>
2322 #include <gsl/gsl_version.h>
2323
2324 static double my_f(const gsl_vector *v, void *params)
2325 {
2326     t_luzar *tl = (t_luzar *)params;
2327     int      i;
2328     double   tol = 1e-16, chi2 = 0;
2329     double   di;
2330     real     k, kp;
2331
2332     for (i = 0; (i < tl->nparams); i++)
2333     {
2334         tl->kkk[i] = gsl_vector_get(v, i);
2335     }
2336     k  = tl->kkk[0];
2337     kp = tl->kkk[1];
2338
2339     for (i = tl->n0; (i < tl->n1); i += tl->ndelta)
2340     {
2341         di = sqr(k*tl->sigma_ct[i]) + sqr(kp*tl->sigma_nt[i]) + sqr(tl->sigma_kt[i]);
2342         /*di = 1;*/
2343         if (di > tol)
2344         {
2345             chi2 += sqr(k*tl->ct[i]-kp*tl->nt[i]-tl->kt[i])/di;
2346         }
2347
2348         else
2349         {
2350             fprintf(stderr, "WARNING: sigma_ct = %g, sigma_nt = %g, sigma_kt = %g\n"
2351                     "di = %g k = %g kp = %g\n", tl->sigma_ct[i],
2352                     tl->sigma_nt[i], tl->sigma_kt[i], di, k, kp);
2353         }
2354     }
2355 #ifdef DEBUG
2356     chi2 = 0.3*sqr(k-0.6)+0.7*sqr(kp-1.3);
2357 #endif
2358     return chi2;
2359 }
2360
2361 static real optimize_luzar_parameters(FILE *fp, t_luzar *tl, int maxiter,
2362                                       real tol)
2363 {
2364     real   size, d2;
2365     int    iter   = 0;
2366     int    status = 0;
2367     int    i;
2368
2369     const gsl_multimin_fminimizer_type *T;
2370     gsl_multimin_fminimizer            *s;
2371
2372     gsl_vector                         *x, *dx;
2373     gsl_multimin_function               my_func;
2374
2375     my_func.f      = &my_f;
2376     my_func.n      = tl->nparams;
2377     my_func.params = (void *) tl;
2378
2379     /* Starting point */
2380     x = gsl_vector_alloc (my_func.n);
2381     for (i = 0; (i < my_func.n); i++)
2382     {
2383         gsl_vector_set (x, i, tl->kkk[i]);
2384     }
2385
2386     /* Step size, different for each of the parameters */
2387     dx = gsl_vector_alloc (my_func.n);
2388     for (i = 0; (i < my_func.n); i++)
2389     {
2390         gsl_vector_set (dx, i, 0.01*tl->kkk[i]);
2391     }
2392
2393     T = gsl_multimin_fminimizer_nmsimplex;
2394     s = gsl_multimin_fminimizer_alloc (T, my_func.n);
2395
2396     gsl_multimin_fminimizer_set (s, &my_func, x, dx);
2397     gsl_vector_free (x);
2398     gsl_vector_free (dx);
2399
2400     if (fp)
2401     {
2402         fprintf(fp, "%5s %12s %12s %12s %12s\n", "Iter", "k", "kp", "NM Size", "Chi2");
2403     }
2404
2405     do
2406     {
2407         iter++;
2408         status = gsl_multimin_fminimizer_iterate (s);
2409
2410         if (status != 0)
2411         {
2412             gmx_fatal(FARGS, "Something went wrong in the iteration in minimizer %s",
2413                       gsl_multimin_fminimizer_name(s));
2414         }
2415
2416         d2     = gsl_multimin_fminimizer_minimum(s);
2417         size   = gsl_multimin_fminimizer_size(s);
2418         status = gsl_multimin_test_size(size, tol);
2419
2420         if (status == GSL_SUCCESS)
2421         {
2422             if (fp)
2423             {
2424                 fprintf(fp, "Minimum found using %s at:\n",
2425                         gsl_multimin_fminimizer_name(s));
2426             }
2427         }
2428
2429         if (fp)
2430         {
2431             fprintf(fp, "%5d", iter);
2432             for (i = 0; (i < my_func.n); i++)
2433             {
2434                 fprintf(fp, " %12.4e", gsl_vector_get (s->x, i));
2435             }
2436             fprintf (fp, " %12.4e %12.4e\n", size, d2);
2437         }
2438     }
2439     while ((status == GSL_CONTINUE) && (iter < maxiter));
2440
2441     gsl_multimin_fminimizer_free (s);
2442
2443     return d2;
2444 }
2445
2446 static real quality_of_fit(real chi2, int N)
2447 {
2448     return gsl_sf_gamma_inc_Q((N-2)/2.0, chi2/2.0);
2449 }
2450
2451 #else
2452 static real optimize_luzar_parameters(FILE *fp, t_luzar *tl, int maxiter,
2453                                       real tol)
2454 {
2455     fprintf(stderr, "This program needs the GNU scientific library to work.\n");
2456
2457     return -1;
2458 }
2459 static real quality_of_fit(real chi2, int N)
2460 {
2461     fprintf(stderr, "This program needs the GNU scientific library to work.\n");
2462
2463     return -1;
2464 }
2465
2466 #endif
2467
2468 static real compute_weighted_rates(int n, real t[], real ct[], real nt[],
2469                                    real kt[], real sigma_ct[], real sigma_nt[],
2470                                    real sigma_kt[], real *k, real *kp,
2471                                    real *sigma_k, real *sigma_kp,
2472                                    real fit_start)
2473 {
2474 #define NK 10
2475     int      i, j;
2476     t_luzar  tl;
2477     real     kkk = 0, kkp = 0, kk2 = 0, kp2 = 0, chi2;
2478
2479     *sigma_k  = 0;
2480     *sigma_kp = 0;
2481
2482     for (i = 0; (i < n); i++)
2483     {
2484         if (t[i] >= fit_start)
2485         {
2486             break;
2487         }
2488     }
2489     tl.n0       = i;
2490     tl.n1       = n;
2491     tl.nparams  = 2;
2492     tl.ndelta   = 1;
2493     tl.t        = t;
2494     tl.ct       = ct;
2495     tl.nt       = nt;
2496     tl.kt       = kt;
2497     tl.sigma_ct = sigma_ct;
2498     tl.sigma_nt = sigma_nt;
2499     tl.sigma_kt = sigma_kt;
2500     tl.kkk[0]   = *k;
2501     tl.kkk[1]   = *kp;
2502
2503     chi2      = optimize_luzar_parameters(debug, &tl, 1000, 1e-3);
2504     *k        = tl.kkk[0];
2505     *kp       = tl.kkk[1] = *kp;
2506     tl.ndelta = NK;
2507     for (j = 0; (j < NK); j++)
2508     {
2509         (void) optimize_luzar_parameters(debug, &tl, 1000, 1e-3);
2510         kkk += tl.kkk[0];
2511         kkp += tl.kkk[1];
2512         kk2 += sqr(tl.kkk[0]);
2513         kp2 += sqr(tl.kkk[1]);
2514         tl.n0++;
2515     }
2516     *sigma_k  = sqrt(kk2/NK - sqr(kkk/NK));
2517     *sigma_kp = sqrt(kp2/NK - sqr(kkp/NK));
2518
2519     return chi2;
2520 }
2521
2522 static void smooth_tail(int n, real t[], real c[], real sigma_c[], real start,
2523                         const output_env_t oenv)
2524 {
2525     FILE *fp;
2526     real  e_1, fitparm[4];
2527     int   i;
2528
2529     e_1 = exp(-1);
2530     for (i = 0; (i < n); i++)
2531     {
2532         if (c[i] < e_1)
2533         {
2534             break;
2535         }
2536     }
2537     if (i < n)
2538     {
2539         fitparm[0] = t[i];
2540     }
2541     else
2542     {
2543         fitparm[0] = 10;
2544     }
2545     fitparm[1] = 0.95;
2546     do_lmfit(n, c, sigma_c, 0, t, start, t[n-1], oenv, bDebugMode(), effnEXP2, fitparm, 0);
2547 }
2548
2549 void analyse_corr(int n, real t[], real ct[], real nt[], real kt[],
2550                   real sigma_ct[], real sigma_nt[], real sigma_kt[],
2551                   real fit_start, real temp, real smooth_tail_start,
2552                   const output_env_t oenv)
2553 {
2554     int        i0, i;
2555     real       k = 1, kp = 1, kow = 1;
2556     real       Q = 0, chi22, chi2, dg, dgp, tau_hb, dtau, tau_rlx, e_1, dt, sigma_k, sigma_kp, ddg;
2557     double     tmp, sn2 = 0, sc2 = 0, sk2 = 0, scn = 0, sck = 0, snk = 0;
2558     gmx_bool   bError = (sigma_ct != NULL) && (sigma_nt != NULL) && (sigma_kt != NULL);
2559
2560     if (smooth_tail_start >= 0)
2561     {
2562         smooth_tail(n, t, ct, sigma_ct, smooth_tail_start, oenv);
2563         smooth_tail(n, t, nt, sigma_nt, smooth_tail_start, oenv);
2564         smooth_tail(n, t, kt, sigma_kt, smooth_tail_start, oenv);
2565     }
2566     for (i0 = 0; (i0 < n-2) && ((t[i0]-t[0]) < fit_start); i0++)
2567     {
2568         ;
2569     }
2570     if (i0 < n-2)
2571     {
2572         for (i = i0; (i < n); i++)
2573         {
2574             sc2 += sqr(ct[i]);
2575             sn2 += sqr(nt[i]);
2576             sk2 += sqr(kt[i]);
2577             sck += ct[i]*kt[i];
2578             snk += nt[i]*kt[i];
2579             scn += ct[i]*nt[i];
2580         }
2581         printf("Hydrogen bond thermodynamics at T = %g K\n", temp);
2582         tmp = (sn2*sc2-sqr(scn));
2583         if ((tmp > 0) && (sn2 > 0))
2584         {
2585             k    = (sn2*sck-scn*snk)/tmp;
2586             kp   = (k*scn-snk)/sn2;
2587             if (bError)
2588             {
2589                 chi2 = 0;
2590                 for (i = i0; (i < n); i++)
2591                 {
2592                     chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
2593                 }
2594                 chi22 = compute_weighted_rates(n, t, ct, nt, kt, sigma_ct, sigma_nt,
2595                                                sigma_kt, &k, &kp,
2596                                                &sigma_k, &sigma_kp, fit_start);
2597                 Q   = quality_of_fit(chi2, 2);
2598                 ddg = BOLTZ*temp*sigma_k/k;
2599                 printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n",
2600                        chi2, Q);
2601                 printf("The Rate and Delta G are followed by an error estimate\n");
2602                 printf("----------------------------------------------------------\n"
2603                        "Type      Rate (1/ps)  Sigma Time (ps)  DG (kJ/mol)  Sigma\n");
2604                 printf("Forward    %10.3f %6.2f   %8.3f  %10.3f %6.2f\n",
2605                        k, sigma_k, 1/k, calc_dg(1/k, temp), ddg);
2606                 ddg = BOLTZ*temp*sigma_kp/kp;
2607                 printf("Backward   %10.3f %6.2f   %8.3f  %10.3f %6.2f\n",
2608                        kp, sigma_kp, 1/kp, calc_dg(1/kp, temp), ddg);
2609             }
2610             else
2611             {
2612                 chi2 = 0;
2613                 for (i = i0; (i < n); i++)
2614                 {
2615                     chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
2616                 }
2617                 printf("Fitting parameters chi^2 = %10g\nQ = %10g\n",
2618                        chi2, Q);
2619                 printf("--------------------------------------------------\n"
2620                        "Type      Rate (1/ps) Time (ps)  DG (kJ/mol)  Chi^2\n");
2621                 printf("Forward    %10.3f   %8.3f  %10.3f  %10g\n",
2622                        k, 1/k, calc_dg(1/k, temp), chi2);
2623                 printf("Backward   %10.3f   %8.3f  %10.3f\n",
2624                        kp, 1/kp, calc_dg(1/kp, temp));
2625             }
2626         }
2627         if (sc2 > 0)
2628         {
2629             kow  = 2*sck/sc2;
2630             printf("One-way    %10.3f   %s%8.3f  %10.3f\n",
2631                    kow, bError ? "       " : "", 1/kow, calc_dg(1/kow, temp));
2632         }
2633         else
2634         {
2635             printf(" - Numerical problems computing HB thermodynamics:\n"
2636                    "sc2 = %g  sn2 = %g  sk2 = %g sck = %g snk = %g scn = %g\n",
2637                    sc2, sn2, sk2, sck, snk, scn);
2638         }
2639         /* Determine integral of the correlation function */
2640         tau_hb = evaluate_integral(n, t, ct, NULL, (t[n-1]-t[0])/2, &dtau);
2641         printf("Integral   %10.3f   %s%8.3f  %10.3f\n", 1/tau_hb,
2642                bError ? "       " : "", tau_hb, calc_dg(tau_hb, temp));
2643         e_1 = exp(-1);
2644         for (i = 0; (i < n-2); i++)
2645         {
2646             if ((ct[i] > e_1) && (ct[i+1] <= e_1))
2647             {
2648                 break;
2649             }
2650         }
2651         if (i < n-2)
2652         {
2653             /* Determine tau_relax from linear interpolation */
2654             tau_rlx = t[i]-t[0] + (e_1-ct[i])*(t[i+1]-t[i])/(ct[i+1]-ct[i]);
2655             printf("Relaxation %10.3f   %8.3f  %s%10.3f\n", 1/tau_rlx,
2656                    tau_rlx, bError ? "       " : "",
2657                    calc_dg(tau_rlx, temp));
2658         }
2659     }
2660     else
2661     {
2662         printf("Correlation functions too short to compute thermodynamics\n");
2663     }
2664 }
2665
2666 void compute_derivative(int nn, real x[], real y[], real dydx[])
2667 {
2668     int j;
2669
2670     /* Compute k(t) = dc(t)/dt */
2671     for (j = 1; (j < nn-1); j++)
2672     {
2673         dydx[j] = (y[j+1]-y[j-1])/(x[j+1]-x[j-1]);
2674     }
2675     /* Extrapolate endpoints */
2676     dydx[0]    = 2*dydx[1]   -  dydx[2];
2677     dydx[nn-1] = 2*dydx[nn-2] - dydx[nn-3];
2678 }
2679
2680 static void parallel_print(int *data, int nThreads)
2681 {
2682     /* This prints the donors on which each tread is currently working. */
2683     int i;
2684
2685     fprintf(stderr, "\r");
2686     for (i = 0; i < nThreads; i++)
2687     {
2688         fprintf(stderr, "%-7i", data[i]);
2689     }
2690 }
2691
2692 static void normalizeACF(real *ct, real *gt, int nhb, int len)
2693 {
2694     real ct_fac, gt_fac;
2695     int  i;
2696
2697     /* Xu and Berne use the same normalization constant */
2698
2699     ct_fac = 1.0/ct[0];
2700     gt_fac = (nhb == 0) ? 0 : 1.0/(real)nhb;
2701
2702     printf("Normalization for c(t) = %g for gh(t) = %g\n", ct_fac, gt_fac);
2703     for (i = 0; i < len; i++)
2704     {
2705         ct[i] *= ct_fac;
2706         if (gt != NULL)
2707         {
2708             gt[i] *= gt_fac;
2709         }
2710     }
2711 }
2712
2713 /* Added argument bContact in do_hbac(...). Also
2714  * added support for -contact in the actual code.
2715  * - Erik Marklund, May 31, 2006 */
2716 /* Changed contact code and added argument R2
2717  * - Erik Marklund, June 29, 2006
2718  */
2719 static void do_hbac(const char *fn, t_hbdata *hb,
2720                     int nDump, gmx_bool bMerge, gmx_bool bContact, real fit_start,
2721                     real temp, gmx_bool R2, real smooth_tail_start, const output_env_t oenv,
2722                     const char *gemType, int nThreads,
2723                     const int NN, const gmx_bool bBallistic, const gmx_bool bGemFit)
2724 {
2725     FILE          *fp;
2726     int            i, j, k, m, n, o, nd, ihb, idist, n2, nn, iter, nSets;
2727     const char    *legNN[]   = {
2728         "Ac(t)",
2729         "Ac'(t)"
2730     };
2731     static char  **legGem;
2732
2733     const char    *legLuzar[] = {
2734         "Ac\\sfin sys\\v{}\\z{}(t)",
2735         "Ac(t)",
2736         "Cc\\scontact,hb\\v{}\\z{}(t)",
2737         "-dAc\\sfs\\v{}\\z{}/dt"
2738     };
2739     gmx_bool       bNorm = FALSE, bOMP = FALSE;
2740     double         nhb   = 0;
2741     int            nhbi  = 0;
2742     real          *rhbex = NULL, *ht, *gt, *ght, *dght, *kt;
2743     real          *ct, *p_ct, tail, tail2, dtail, ct_fac, ght_fac, *cct;
2744     const real     tol     = 1e-3;
2745     int            nframes = hb->nframes, nf;
2746     unsigned int **h       = NULL, **g = NULL;
2747     int            nh, nhbonds, nhydro, ngh;
2748     t_hbond       *hbh;
2749     PSTYPE         p, *pfound = NULL, np;
2750     t_pShift      *pHist;
2751     int           *ptimes   = NULL, *poff = NULL, anhb, n0, mMax = INT_MIN;
2752     real         **rHbExGem = NULL;
2753     gmx_bool       c;
2754     int            acType;
2755     t_E           *E;
2756     double        *ctdouble, *timedouble, *fittedct;
2757     double         fittolerance = 0.1;
2758     int           *dondata      = NULL, thisThread;
2759
2760     enum {
2761         AC_NONE, AC_NN, AC_GEM, AC_LUZAR
2762     };
2763
2764 #ifdef GMX_OPENMP
2765     bOMP = TRUE;
2766 #else
2767     bOMP = FALSE;
2768 #endif
2769
2770     printf("Doing autocorrelation ");
2771
2772     /* Decide what kind of ACF calculations to do. */
2773     if (NN > NN_NONE && NN < NN_NR)
2774     {
2775 #ifdef HAVE_NN_LOOPS
2776         acType = AC_NN;
2777         printf("using the energy estimate.\n");
2778 #else
2779         acType = AC_NONE;
2780         printf("Can't do the NN-loop. Yet.\n");
2781 #endif
2782     }
2783     else if (hb->bGem)
2784     {
2785         acType = AC_GEM;
2786         printf("according to the reversible geminate recombination model by Omer Markowitch.\n");
2787
2788         nSets = 1 + (bBallistic ? 1 : 0) + (bGemFit ? 1 : 0);
2789         snew(legGem, nSets);
2790         for (i = 0; i < nSets; i++)
2791         {
2792             snew(legGem[i], 128);
2793         }
2794         sprintf(legGem[0], "Ac\\s%s\\v{}\\z{}(t)", gemType);
2795         if (bBallistic)
2796         {
2797             sprintf(legGem[1], "Ac'(t)");
2798         }
2799         if (bGemFit)
2800         {
2801             sprintf(legGem[(bBallistic ? 3 : 2)], "Ac\\s%s,fit\\v{}\\z{}(t)", gemType);
2802         }
2803
2804     }
2805     else
2806     {
2807         acType = AC_LUZAR;
2808         printf("according to the theory of Luzar and Chandler.\n");
2809     }
2810     fflush(stdout);
2811
2812     /* build hbexist matrix in reals for autocorr */
2813     /* Allocate memory for computing ACF (rhbex) and aggregating the ACF (ct) */
2814     n2 = 1;
2815     while (n2 < nframes)
2816     {
2817         n2 *= 2;
2818     }
2819
2820     nn = nframes/2;
2821
2822     if (acType != AC_NN || bOMP)
2823     {
2824         snew(h, hb->maxhydro);
2825         snew(g, hb->maxhydro);
2826     }
2827
2828     /* Dump hbonds for debugging */
2829     dump_ac(hb, bMerge || bContact, nDump);
2830
2831     /* Total number of hbonds analyzed here */
2832     nhbonds = 0;
2833     ngh     = 0;
2834     anhb    = 0;
2835
2836     if (acType != AC_LUZAR && bOMP)
2837     {
2838         nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
2839
2840         gmx_omp_set_num_threads(nThreads);
2841         snew(dondata, nThreads);
2842         for (i = 0; i < nThreads; i++)
2843         {
2844             dondata[i] = -1;
2845         }
2846         printf("ACF calculations parallelized with OpenMP using %i threads.\n"
2847                "Expect close to linear scaling over this donor-loop.\n", nThreads);
2848         fflush(stdout);
2849         fprintf(stderr, "Donors: [thread no]\n");
2850         {
2851             char tmpstr[7];
2852             for (i = 0; i < nThreads; i++)
2853             {
2854                 snprintf(tmpstr, 7, "[%i]", i);
2855                 fprintf(stderr, "%-7s", tmpstr);
2856             }
2857         }
2858         fprintf(stderr, "\n");
2859     }
2860
2861
2862     /* Build the ACF according to acType */
2863     switch (acType)
2864     {
2865
2866         case AC_NN:
2867 #ifdef HAVE_NN_LOOPS
2868             /* Here we're using the estimated energy for the hydrogen bonds. */
2869             snew(ct, nn);
2870
2871 #pragma omp parallel \
2872             private(i, j, k, nh, E, rhbex, thisThread) \
2873             default(shared)
2874             {
2875 #pragma omp barrier
2876                 thisThread = gmx_omp_get_thread_num();
2877                 rhbex      = NULL;
2878
2879                 snew(rhbex, n2);
2880                 memset(rhbex, 0, n2*sizeof(real)); /* Trust no-one, not even malloc()! */
2881
2882 #pragma omp barrier
2883 #pragma omp for schedule (dynamic)
2884                 for (i = 0; i < hb->d.nrd; i++) /* loop over donors */
2885                 {
2886                     if (bOMP)
2887                     {
2888 #pragma omp critical
2889                         {
2890                             dondata[thisThread] = i;
2891                             parallel_print(dondata, nThreads);
2892                         }
2893                     }
2894                     else
2895                     {
2896                         fprintf(stderr, "\r %i", i);
2897                     }
2898
2899                     for (j = 0; j < hb->a.nra; j++)              /* loop over acceptors */
2900                     {
2901                         for (nh = 0; nh < hb->d.nhydro[i]; nh++) /* loop over donors' hydrogens */
2902                         {
2903                             E = hb->hbE.E[i][j][nh];
2904                             if (E != NULL)
2905                             {
2906                                 for (k = 0; k < nframes; k++)
2907                                 {
2908                                     if (E[k] != NONSENSE_E)
2909                                     {
2910                                         rhbex[k] = (real)E[k];
2911                                     }
2912                                 }
2913
2914                                 low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &(rhbex), hb->time[1]-hb->time[0],
2915                                                 eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0, 1);
2916 #pragma omp critical
2917                                 {
2918                                     for (k = 0; (k < nn); k++)
2919                                     {
2920                                         ct[k] += rhbex[k];
2921                                     }
2922                                 }
2923                             }
2924                         } /* k loop */
2925                     }     /* j loop */
2926                 }         /* i loop */
2927                 sfree(rhbex);
2928 #pragma omp barrier
2929             }
2930
2931             if (bOMP)
2932             {
2933                 sfree(dondata);
2934             }
2935             normalizeACF(ct, NULL, 0, nn);
2936             snew(ctdouble, nn);
2937             snew(timedouble, nn);
2938             for (j = 0; j < nn; j++)
2939             {
2940                 timedouble[j] = (double)(hb->time[j]);
2941                 ctdouble[j]   = (double)(ct[j]);
2942             }
2943
2944             /* Remove ballistic term */
2945             /* Ballistic component removal and fitting to the reversible geminate recombination model
2946              * will be taken out for the time being. First of all, one can remove the ballistic
2947              * component with g_analyze afterwards. Secondly, and more importantly, there are still
2948              * problems with the robustness of the fitting to the model. More work is needed.
2949              * A third reason is that we're currently using gsl for this and wish to reduce dependence
2950              * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with
2951              * a BSD-licence that can do the job.
2952              *
2953              * - Erik Marklund, June 18 2010.
2954              */
2955 /*         if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */
2956 /*             takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */
2957 /*         else */
2958 /*             printf("\nNumber of data points is less than the number of parameters to fit\n." */
2959 /*                    "The system is underdetermined, hence no ballistic term can be found.\n\n"); */
2960
2961             fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)");
2962             xvgr_legend(fp, asize(legNN), legNN);
2963
2964             for (j = 0; (j < nn); j++)
2965             {
2966                 fprintf(fp, "%10g  %10g %10g\n",
2967                         hb->time[j]-hb->time[0],
2968                         ct[j],
2969                         ctdouble[j]);
2970             }
2971             xvgrclose(fp);
2972             sfree(ct);
2973             sfree(ctdouble);
2974             sfree(timedouble);
2975 #endif             /* HAVE_NN_LOOPS */
2976             break; /* case AC_NN */
2977
2978         case AC_GEM:
2979             snew(ct, 2*n2);
2980             memset(ct, 0, 2*n2*sizeof(real));
2981 #ifndef GMX_OPENMP
2982             fprintf(stderr, "Donor:\n");
2983 #define __ACDATA ct
2984 #else
2985 #define __ACDATA p_ct
2986 #endif
2987
2988 #pragma omp parallel \
2989             private(i, k, nh, hbh, pHist, h, g, n0, nf, np, j, m, \
2990             pfound, poff, rHbExGem, p, ihb, mMax, \
2991             thisThread, p_ct) \
2992             default(shared)
2993             { /* ##########  THE START OF THE ENORMOUS PARALLELIZED BLOCK!  ########## */
2994                 h          = NULL;
2995                 g          = NULL;
2996                 thisThread = gmx_omp_get_thread_num();
2997                 snew(h, hb->maxhydro);
2998                 snew(g, hb->maxhydro);
2999                 mMax     = INT_MIN;
3000                 rHbExGem = NULL;
3001                 poff     = NULL;
3002                 pfound   = NULL;
3003                 p_ct     = NULL;
3004                 snew(p_ct, 2*n2);
3005                 memset(p_ct, 0, 2*n2*sizeof(real));
3006
3007                 /* I'm using a chunk size of 1, since I expect      \
3008                  * the overhead to be really small compared         \
3009                  * to the actual calculations                       \ */
3010 #pragma omp for schedule(dynamic,1) nowait
3011                 for (i = 0; i < hb->d.nrd; i++)
3012                 {
3013
3014                     if (bOMP)
3015                     {
3016 #pragma omp critical
3017                         {
3018                             dondata[thisThread] = i;
3019                             parallel_print(dondata, nThreads);
3020                         }
3021                     }
3022                     else
3023                     {
3024                         fprintf(stderr, "\r %i", i);
3025                     }
3026                     for (k = 0; k < hb->a.nra; k++)
3027                     {
3028                         for (nh = 0; nh < ((bMerge || bContact) ? 1 : hb->d.nhydro[i]); nh++)
3029                         {
3030                             hbh = hb->hbmap[i][k];
3031                             if (hbh)
3032                             {
3033                                 /* Note that if hb->per->gemtype==gemDD, then distances will be stored in
3034                                  * hb->hbmap[d][a].h array anyway, because the contact flag will be set.
3035                                  * hence, it's only with the gemAD mode that hb->hbmap[d][a].g will be used. */
3036                                 pHist = &(hb->per->pHist[i][k]);
3037                                 if (ISHB(hbh->history[nh]) && pHist->len != 0)
3038                                 {
3039
3040                                     {
3041                                         h[nh] = hbh->h[nh];
3042                                         g[nh] = hb->per->gemtype == gemAD ? hbh->g[nh] : NULL;
3043                                     }
3044                                     n0 = hbh->n0;
3045                                     nf = hbh->nframes;
3046                                     /* count the number of periodic shifts encountered and store
3047                                      * them in separate arrays. */
3048                                     np = 0;
3049                                     for (j = 0; j < pHist->len; j++)
3050                                     {
3051                                         p = pHist->p[j];
3052                                         for (m = 0; m <= np; m++)
3053                                         {
3054                                             if (m == np) /* p not recognized in list. Add it and set up new array. */
3055                                             {
3056                                                 np++;
3057                                                 if (np > hb->per->nper)
3058                                                 {
3059                                                     gmx_fatal(FARGS, "Too many pshifts. Something's utterly wrong here.");
3060                                                 }
3061                                                 if (m >= mMax) /* Extend the arrays.
3062                                                                 * Doing it like this, using mMax to keep track of the sizes,
3063                                                                 * eleviates the need for freeing and re-allocating the arrays
3064                                                                 * when taking on the next donor-acceptor pair */
3065                                                 {
3066                                                     mMax = m;
3067                                                     srenew(pfound, np);   /* The list of found periodic shifts. */
3068                                                     srenew(rHbExGem, np); /* The hb existence functions (-aver_hb). */
3069                                                     snew(rHbExGem[m], 2*n2);
3070                                                     srenew(poff, np);
3071                                                 }
3072
3073                                                 {
3074                                                     if (rHbExGem != NULL && rHbExGem[m] != NULL)
3075                                                     {
3076                                                         /* This must be done, as this array was most likey
3077                                                          * used to store stuff in some previous iteration. */
3078                                                         memset(rHbExGem[m], 0, (sizeof(real)) * (2*n2));
3079                                                     }
3080                                                     else
3081                                                     {
3082                                                         fprintf(stderr, "rHbExGem not initialized! m = %i\n", m);
3083                                                     }
3084                                                 }
3085                                                 pfound[m] = p;
3086                                                 poff[m]   = -1;
3087
3088                                                 break;
3089                                             } /* m==np */
3090                                             if (p == pfound[m])
3091                                             {
3092                                                 break;
3093                                             }
3094                                         } /* m: Loop over found shifts */
3095                                     }     /* j: Loop over shifts */
3096
3097                                     /* Now unpack and disentangle the existence funtions. */
3098                                     for (j = 0; j < nf; j++)
3099                                     {
3100                                         /* i:       donor,
3101                                          * k:       acceptor
3102                                          * nh:      hydrogen
3103                                          * j:       time
3104                                          * p:       periodic shift
3105                                          * pfound:  list of periodic shifts found for this pair.
3106                                          * poff:    list of frame offsets; that is, the first
3107                                          *          frame a hbond has a particular periodic shift. */
3108                                         p = getPshift(*pHist, j+n0);
3109                                         if (p != -1)
3110                                         {
3111                                             for (m = 0; m < np; m++)
3112                                             {
3113                                                 if (pfound[m] == p)
3114                                                 {
3115                                                     break;
3116                                                 }
3117                                                 if (m == (np-1))
3118                                                 {
3119                                                     gmx_fatal(FARGS, "Shift not found, but must be there.");
3120                                                 }
3121                                             }
3122
3123                                             ihb = is_hb(h[nh], j) || ((hb->per->gemtype != gemAD || j == 0) ? FALSE : is_hb(g[nh], j));
3124                                             if (ihb)
3125                                             {
3126                                                 if (poff[m] == -1)
3127                                                 {
3128                                                     poff[m] = j; /* Here's where the first hbond with shift p is,
3129                                                                   * relative to the start of h[0].*/
3130                                                 }
3131                                                 if (j < poff[m])
3132                                                 {
3133                                                     gmx_fatal(FARGS, "j<poff[m]");
3134                                                 }
3135                                                 rHbExGem[m][j-poff[m]] += 1;
3136                                             }
3137                                         }
3138                                     }
3139
3140                                     /* Now, build ac. */
3141                                     for (m = 0; m < np; m++)
3142                                     {
3143                                         if (rHbExGem[m][0] > 0  && n0+poff[m] < nn /*  && m==0 */)
3144                                         {
3145                                             low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &(rHbExGem[m]), hb->time[1]-hb->time[0],
3146                                                             eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0);
3147                                             for (j = 0; (j < nn); j++)
3148                                             {
3149                                                 __ACDATA[j] += rHbExGem[m][j];
3150                                             }
3151                                         }
3152                                     } /* Building of ac. */
3153                                 }     /* if (ISHB(...*/
3154                             }         /* if (hbh) */
3155                         }             /* hydrogen loop */
3156                     }                 /* acceptor loop */
3157                 }                     /* donor loop */
3158
3159                 for (m = 0; m <= mMax; m++)
3160                 {
3161                     sfree(rHbExGem[m]);
3162                 }
3163                 sfree(pfound);
3164                 sfree(poff);
3165                 sfree(rHbExGem);
3166
3167                 sfree(h);
3168                 sfree(g);
3169
3170                 if (bOMP)
3171                 {
3172 #pragma omp critical
3173                     {
3174                         for (i = 0; i < nn; i++)
3175                         {
3176                             ct[i] += p_ct[i];
3177                         }
3178                     }
3179                     sfree(p_ct);
3180                 }
3181
3182             } /* ########## THE END OF THE ENORMOUS PARALLELIZED BLOCK ########## */
3183             if (bOMP)
3184             {
3185                 sfree(dondata);
3186             }
3187
3188             normalizeACF(ct, NULL, 0, nn);
3189
3190             fprintf(stderr, "\n\nACF successfully calculated.\n");
3191
3192             /* Use this part to fit to geminate recombination - JCP 129, 84505 (2008) */
3193
3194             snew(ctdouble, nn);
3195             snew(timedouble, nn);
3196             snew(fittedct, nn);
3197
3198             for (j = 0; j < nn; j++)
3199             {
3200                 timedouble[j] = (double)(hb->time[j]);
3201                 ctdouble[j]   = (double)(ct[j]);
3202             }
3203
3204             /* Remove ballistic term */
3205             /* Ballistic component removal and fitting to the reversible geminate recombination model
3206              * will be taken out for the time being. First of all, one can remove the ballistic
3207              * component with g_analyze afterwards. Secondly, and more importantly, there are still
3208              * problems with the robustness of the fitting to the model. More work is needed.
3209              * A third reason is that we're currently using gsl for this and wish to reduce dependence
3210              * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with
3211              * a BSD-licence that can do the job.
3212              *
3213              * - Erik Marklund, June 18 2010.
3214              */
3215 /*         if (bBallistic) { */
3216 /*             if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */
3217 /*                 takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */
3218 /*             else */
3219 /*                 printf("\nNumber of data points is less than the number of parameters to fit\n." */
3220 /*                        "The system is underdetermined, hence no ballistic term can be found.\n\n"); */
3221 /*         } */
3222 /*         if (bGemFit) */
3223 /*             fitGemRecomb(ctdouble, timedouble, &fittedct, nn, params); */
3224
3225
3226             if (bContact)
3227             {
3228                 fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3229             }
3230             else
3231             {
3232                 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3233             }
3234             xvgr_legend(fp, asize(legGem), (const char**)legGem, oenv);
3235
3236             for (j = 0; (j < nn); j++)
3237             {
3238                 fprintf(fp, "%10g  %10g", hb->time[j]-hb->time[0], ct[j]);
3239                 if (bBallistic)
3240                 {
3241                     fprintf(fp, "  %10g", ctdouble[j]);
3242                 }
3243                 if (bGemFit)
3244                 {
3245                     fprintf(fp, "  %10g", fittedct[j]);
3246                 }
3247                 fprintf(fp, "\n");
3248             }
3249             xvgrclose(fp);
3250
3251             sfree(ctdouble);
3252             sfree(timedouble);
3253             sfree(fittedct);
3254             sfree(ct);
3255
3256             break; /* case AC_GEM */
3257
3258         case AC_LUZAR:
3259             snew(rhbex, 2*n2);
3260             snew(ct, 2*n2);
3261             snew(gt, 2*n2);
3262             snew(ht, 2*n2);
3263             snew(ght, 2*n2);
3264             snew(dght, 2*n2);
3265
3266             snew(kt, nn);
3267             snew(cct, nn);
3268
3269             for (i = 0; (i < hb->d.nrd); i++)
3270             {
3271                 for (k = 0; (k < hb->a.nra); k++)
3272                 {
3273                     nhydro = 0;
3274                     hbh    = hb->hbmap[i][k];
3275
3276                     if (hbh)
3277                     {
3278                         if (bMerge || bContact)
3279                         {
3280                             if (ISHB(hbh->history[0]))
3281                             {
3282                                 h[0]   = hbh->h[0];
3283                                 g[0]   = hbh->g[0];
3284                                 nhydro = 1;
3285                             }
3286                         }
3287                         else
3288                         {
3289                             for (m = 0; (m < hb->maxhydro); m++)
3290                             {
3291                                 if (bContact ? ISDIST(hbh->history[m]) : ISHB(hbh->history[m]))
3292                                 {
3293                                     g[nhydro] = hbh->g[m];
3294                                     h[nhydro] = hbh->h[m];
3295                                     nhydro++;
3296                                 }
3297                             }
3298                         }
3299
3300                         nf = hbh->nframes;
3301                         for (nh = 0; (nh < nhydro); nh++)
3302                         {
3303                             int nrint = bContact ? hb->nrdist : hb->nrhb;
3304                             if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint))
3305                             {
3306                                 fprintf(stderr, "\rACF %d/%d", nhbonds+1, nrint);
3307                             }
3308                             nhbonds++;
3309                             for (j = 0; (j < nframes); j++)
3310                             {
3311                                 /* Changed '<' into '<=' below, just like I did in
3312                                    the hbm-output-loop in the gmx_hbond() block.
3313                                    - Erik Marklund, May 31, 2006 */
3314                                 if (j <= nf)
3315                                 {
3316                                     ihb   = is_hb(h[nh], j);
3317                                     idist = is_hb(g[nh], j);
3318                                 }
3319                                 else
3320                                 {
3321                                     ihb = idist = 0;
3322                                 }
3323                                 rhbex[j] = ihb;
3324                                 /* For contacts: if a second cut-off is provided, use it,
3325                                  * otherwise use g(t) = 1-h(t) */
3326                                 if (!R2 && bContact)
3327                                 {
3328                                     gt[j]  = 1-ihb;
3329                                 }
3330                                 else
3331                                 {
3332                                     gt[j]  = idist*(1-ihb);
3333                                 }
3334                                 ht[j]    = rhbex[j];
3335                                 nhb     += ihb;
3336                             }
3337
3338
3339                             /* The autocorrelation function is normalized after summation only */
3340                             low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &rhbex, hb->time[1]-hb->time[0],
3341                                             eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0);
3342
3343                             /* Cross correlation analysis for thermodynamics */
3344                             for (j = nframes; (j < n2); j++)
3345                             {
3346                                 ht[j] = 0;
3347                                 gt[j] = 0;
3348                             }
3349
3350                             cross_corr(n2, ht, gt, dght);
3351
3352                             for (j = 0; (j < nn); j++)
3353                             {
3354                                 ct[j]  += rhbex[j];
3355                                 ght[j] += dght[j];
3356                             }
3357                         }
3358                     }
3359                 }
3360             }
3361             fprintf(stderr, "\n");
3362             sfree(h);
3363             sfree(g);
3364             normalizeACF(ct, ght, nhb, nn);
3365
3366             /* Determine tail value for statistics */
3367             tail  = 0;
3368             tail2 = 0;
3369             for (j = nn/2; (j < nn); j++)
3370             {
3371                 tail  += ct[j];
3372                 tail2 += ct[j]*ct[j];
3373             }
3374             tail  /= (nn - nn/2);
3375             tail2 /= (nn - nn/2);
3376             dtail  = sqrt(tail2-tail*tail);
3377
3378             /* Check whether the ACF is long enough */
3379             if (dtail > tol)
3380             {
3381                 printf("\nWARNING: Correlation function is probably not long enough\n"
3382                        "because the standard deviation in the tail of C(t) > %g\n"
3383                        "Tail value (average C(t) over second half of acf): %g +/- %g\n",
3384                        tol, tail, dtail);
3385             }
3386             for (j = 0; (j < nn); j++)
3387             {
3388                 cct[j] = ct[j];
3389                 ct[j]  = (cct[j]-tail)/(1-tail);
3390             }
3391             /* Compute negative derivative k(t) = -dc(t)/dt */
3392             compute_derivative(nn, hb->time, ct, kt);
3393             for (j = 0; (j < nn); j++)
3394             {
3395                 kt[j] = -kt[j];
3396             }
3397
3398
3399             if (bContact)
3400             {
3401                 fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3402             }
3403             else
3404             {
3405                 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3406             }
3407             xvgr_legend(fp, asize(legLuzar), legLuzar, oenv);
3408
3409
3410             for (j = 0; (j < nn); j++)
3411             {
3412                 fprintf(fp, "%10g  %10g  %10g  %10g  %10g\n",
3413                         hb->time[j]-hb->time[0], ct[j], cct[j], ght[j], kt[j]);
3414             }
3415             ffclose(fp);
3416
3417             analyse_corr(nn, hb->time, ct, ght, kt, NULL, NULL, NULL,
3418                          fit_start, temp, smooth_tail_start, oenv);
3419
3420             do_view(oenv, fn, NULL);
3421             sfree(rhbex);
3422             sfree(ct);
3423             sfree(gt);
3424             sfree(ht);
3425             sfree(ght);
3426             sfree(dght);
3427             sfree(cct);
3428             sfree(kt);
3429             /* sfree(h); */
3430 /*         sfree(g); */
3431
3432             break; /* case AC_LUZAR */
3433
3434         default:
3435             gmx_fatal(FARGS, "Unrecognized type of ACF-calulation. acType = %i.", acType);
3436     } /* switch (acType) */
3437 }
3438
3439 static void init_hbframe(t_hbdata *hb, int nframes, real t)
3440 {
3441     int i, j, m;
3442
3443     hb->time[nframes]   = t;
3444     hb->nhb[nframes]    = 0;
3445     hb->ndist[nframes]  = 0;
3446     for (i = 0; (i < max_hx); i++)
3447     {
3448         hb->nhx[nframes][i] = 0;
3449     }
3450     /* Loop invalidated */
3451     if (hb->bHBmap && 0)
3452     {
3453         for (i = 0; (i < hb->d.nrd); i++)
3454         {
3455             for (j = 0; (j < hb->a.nra); j++)
3456             {
3457                 for (m = 0; (m < hb->maxhydro); m++)
3458                 {
3459                     if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[m])
3460                     {
3461                         set_hb(hb, i, m, j, nframes, HB_NO);
3462                     }
3463                 }
3464             }
3465         }
3466     }
3467     /*set_hb(hb->hbmap[i][j]->h[m],nframes-hb->hbmap[i][j]->n0,HB_NO);*/
3468 }
3469
3470 static void analyse_donor_props(const char *fn, t_hbdata *hb, int nframes, real t,
3471                                 const output_env_t oenv)
3472 {
3473     static FILE *fp    = NULL;
3474     const char  *leg[] = { "Nbound", "Nfree" };
3475     int          i, j, k, nbound, nb, nhtot;
3476
3477     if (!fn)
3478     {
3479         return;
3480     }
3481     if (!fp)
3482     {
3483         fp = xvgropen(fn, "Donor properties", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
3484         xvgr_legend(fp, asize(leg), leg, oenv);
3485     }
3486     nbound = 0;
3487     nhtot  = 0;
3488     for (i = 0; (i < hb->d.nrd); i++)
3489     {
3490         for (k = 0; (k < hb->d.nhydro[i]); k++)
3491         {
3492             nb = 0;
3493             nhtot++;
3494             for (j = 0; (j < hb->a.nra) && (nb == 0); j++)
3495             {
3496                 if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[k] &&
3497                     is_hb(hb->hbmap[i][j]->h[k], nframes))
3498                 {
3499                     nb = 1;
3500                 }
3501             }
3502             nbound += nb;
3503         }
3504     }
3505     fprintf(fp, "%10.3e  %6d  %6d\n", t, nbound, nhtot-nbound);
3506 }
3507
3508 static void dump_hbmap(t_hbdata *hb,
3509                        int nfile, t_filenm fnm[], gmx_bool bTwo,
3510                        gmx_bool bContact, int isize[], int *index[], char *grpnames[],
3511                        t_atoms *atoms)
3512 {
3513     FILE    *fp, *fplog;
3514     int      ddd, hhh, aaa, i, j, k, m, grp;
3515     char     ds[32], hs[32], as[32];
3516     gmx_bool first;
3517
3518     fp = opt2FILE("-hbn", nfile, fnm, "w");
3519     if (opt2bSet("-g", nfile, fnm))
3520     {
3521         fplog = ffopen(opt2fn("-g", nfile, fnm), "w");
3522         fprintf(fplog, "# %10s  %12s  %12s\n", "Donor", "Hydrogen", "Acceptor");
3523     }
3524     else
3525     {
3526         fplog = NULL;
3527     }
3528     for (grp = gr0; grp <= (bTwo ? gr1 : gr0); grp++)
3529     {
3530         fprintf(fp, "[ %s ]", grpnames[grp]);
3531         for (i = 0; i < isize[grp]; i++)
3532         {
3533             fprintf(fp, (i%15) ? " " : "\n");
3534             fprintf(fp, " %4u", index[grp][i]+1);
3535         }
3536         fprintf(fp, "\n");
3537         /*
3538            Added -contact support below.
3539            - Erik Marklund, May 29, 2006
3540          */
3541         if (!bContact)
3542         {
3543             fprintf(fp, "[ donors_hydrogens_%s ]\n", grpnames[grp]);
3544             for (i = 0; (i < hb->d.nrd); i++)
3545             {
3546                 if (hb->d.grp[i] == grp)
3547                 {
3548                     for (j = 0; (j < hb->d.nhydro[i]); j++)
3549                     {
3550                         fprintf(fp, " %4u %4u", hb->d.don[i]+1,
3551                                 hb->d.hydro[i][j]+1);
3552                     }
3553                     fprintf(fp, "\n");
3554                 }
3555             }
3556             first = TRUE;
3557             fprintf(fp, "[ acceptors_%s ]", grpnames[grp]);
3558             for (i = 0; (i < hb->a.nra); i++)
3559             {
3560                 if (hb->a.grp[i] == grp)
3561                 {
3562                     fprintf(fp, (i%15 && !first) ? " " : "\n");
3563                     fprintf(fp, " %4u", hb->a.acc[i]+1);
3564                     first = FALSE;
3565                 }
3566             }
3567             fprintf(fp, "\n");
3568         }
3569     }
3570     if (bTwo)
3571     {
3572         fprintf(fp, bContact ? "[ contacts_%s-%s ]\n" :
3573                 "[ hbonds_%s-%s ]\n", grpnames[0], grpnames[1]);
3574     }
3575     else
3576     {
3577         fprintf(fp, bContact ? "[ contacts_%s ]" : "[ hbonds_%s ]\n", grpnames[0]);
3578     }
3579
3580     for (i = 0; (i < hb->d.nrd); i++)
3581     {
3582         ddd = hb->d.don[i];
3583         for (k = 0; (k < hb->a.nra); k++)
3584         {
3585             aaa = hb->a.acc[k];
3586             for (m = 0; (m < hb->d.nhydro[i]); m++)
3587             {
3588                 if (hb->hbmap[i][k] && ISHB(hb->hbmap[i][k]->history[m]))
3589                 {
3590                     sprintf(ds, "%s", mkatomname(atoms, ddd));
3591                     sprintf(as, "%s", mkatomname(atoms, aaa));
3592                     if (bContact)
3593                     {
3594                         fprintf(fp, " %6u %6u\n", ddd+1, aaa+1);
3595                         if (fplog)
3596                         {
3597                             fprintf(fplog, "%12s  %12s\n", ds, as);
3598                         }
3599                     }
3600                     else
3601                     {
3602                         hhh = hb->d.hydro[i][m];
3603                         sprintf(hs, "%s", mkatomname(atoms, hhh));
3604                         fprintf(fp, " %6u %6u %6u\n", ddd+1, hhh+1, aaa+1);
3605                         if (fplog)
3606                         {
3607                             fprintf(fplog, "%12s  %12s  %12s\n", ds, hs, as);
3608                         }
3609                     }
3610                 }
3611             }
3612         }
3613     }
3614     ffclose(fp);
3615     if (fplog)
3616     {
3617         ffclose(fplog);
3618     }
3619 }
3620
3621 /* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template.
3622  * It mimics add_frames() and init_frame() to some extent. */
3623 static void sync_hbdata(t_hbdata *p_hb, int nframes)
3624 {
3625     int i;
3626     if (nframes >= p_hb->max_frames)
3627     {
3628         p_hb->max_frames += 4096;
3629         srenew(p_hb->nhb,   p_hb->max_frames);
3630         srenew(p_hb->ndist, p_hb->max_frames);
3631         srenew(p_hb->n_bound, p_hb->max_frames);
3632         srenew(p_hb->nhx, p_hb->max_frames);
3633         if (p_hb->bDAnr)
3634         {
3635             srenew(p_hb->danr, p_hb->max_frames);
3636         }
3637         memset(&(p_hb->nhb[nframes]),   0, sizeof(int) * (p_hb->max_frames-nframes));
3638         memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
3639         p_hb->nhb[nframes]   = 0;
3640         p_hb->ndist[nframes] = 0;
3641
3642     }
3643     p_hb->nframes = nframes;
3644 /*     for (i=0;) */
3645 /*     { */
3646 /*         p_hb->nhx[nframes][i] */
3647 /*     } */
3648     memset(&(p_hb->nhx[nframes]), 0, sizeof(int)*max_hx); /* zero the helix count for this frame */
3649
3650     /* hb->per will remain constant througout the frame loop,
3651      * even though the data its members point to will change,
3652      * hence no need for re-syncing. */
3653 }
3654
3655 int gmx_hbond(int argc, char *argv[])
3656 {
3657     const char        *desc[] = {
3658         "[THISMODULE] computes and analyzes hydrogen bonds. Hydrogen bonds are",
3659         "determined based on cutoffs for the angle Hydrogen - Donor - Acceptor",
3660         "(zero is extended) and the distance Donor - Acceptor",
3661         "(or Hydrogen - Acceptor using [TT]-noda[tt]).",
3662         "OH and NH groups are regarded as donors, O is an acceptor always,",
3663         "N is an acceptor by default, but this can be switched using",
3664         "[TT]-nitacc[tt]. Dummy hydrogen atoms are assumed to be connected",
3665         "to the first preceding non-hydrogen atom.[PAR]",
3666
3667         "You need to specify two groups for analysis, which must be either",
3668         "identical or non-overlapping. All hydrogen bonds between the two",
3669         "groups are analyzed.[PAR]",
3670
3671         "If you set [TT]-shell[tt], you will be asked for an additional index group",
3672         "which should contain exactly one atom. In this case, only hydrogen",
3673         "bonds between atoms within the shell distance from the one atom are",
3674         "considered.[PAR]",
3675
3676         "With option -ac, rate constants for hydrogen bonding can be derived with the model of Luzar and Chandler",
3677         "(Nature 394, 1996; J. Chem. Phys. 113:23, 2000) or that of Markovitz and Agmon (J. Chem. Phys 129, 2008).",
3678         "If contact kinetics are analyzed by using the -contact option, then",
3679         "n(t) can be defined as either all pairs that are not within contact distance r at time t",
3680         "(corresponding to leaving the -r2 option at the default value 0) or all pairs that",
3681         "are within distance r2 (corresponding to setting a second cut-off value with option -r2).",
3682         "See mentioned literature for more details and definitions."
3683         "[PAR]",
3684
3685         /*    "It is also possible to analyse specific hydrogen bonds with",
3686               "[TT]-sel[tt]. This index file must contain a group of atom triplets",
3687               "Donor Hydrogen Acceptor, in the following way:[PAR]",
3688          */
3689         "[TT]",
3690         "[ selected ][BR]",
3691         "     20    21    24[BR]",
3692         "     25    26    29[BR]",
3693         "      1     3     6[BR]",
3694         "[tt][BR]",
3695         "Note that the triplets need not be on separate lines.",
3696         "Each atom triplet specifies a hydrogen bond to be analyzed,",
3697         "note also that no check is made for the types of atoms.[PAR]",
3698
3699         "[BB]Output:[bb][BR]",
3700         "[TT]-num[tt]:  number of hydrogen bonds as a function of time.[BR]",
3701         "[TT]-ac[tt]:   average over all autocorrelations of the existence",
3702         "functions (either 0 or 1) of all hydrogen bonds.[BR]",
3703         "[TT]-dist[tt]: distance distribution of all hydrogen bonds.[BR]",
3704         "[TT]-ang[tt]:  angle distribution of all hydrogen bonds.[BR]",
3705         "[TT]-hx[tt]:   the number of n-n+i hydrogen bonds as a function of time",
3706         "where n and n+i stand for residue numbers and i ranges from 0 to 6.",
3707         "This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated",
3708         "with helices in proteins.[BR]",
3709         "[TT]-hbn[tt]:  all selected groups, donors, hydrogens and acceptors",
3710         "for selected groups, all hydrogen bonded atoms from all groups and",
3711         "all solvent atoms involved in insertion.[BR]",
3712         "[TT]-hbm[tt]:  existence matrix for all hydrogen bonds over all",
3713         "frames, this also contains information on solvent insertion",
3714         "into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]",
3715         "index file.[BR]",
3716         "[TT]-dan[tt]: write out the number of donors and acceptors analyzed for",
3717         "each timeframe. This is especially useful when using [TT]-shell[tt].[BR]",
3718         "[TT]-nhbdist[tt]: compute the number of HBonds per hydrogen in order to",
3719         "compare results to Raman Spectroscopy.",
3720         "[PAR]",
3721         "Note: options [TT]-ac[tt], [TT]-life[tt], [TT]-hbn[tt] and [TT]-hbm[tt]",
3722         "require an amount of memory proportional to the total numbers of donors",
3723         "times the total number of acceptors in the selected group(s)."
3724     };
3725
3726     static real        acut     = 30, abin = 1, rcut = 0.35, r2cut = 0, rbin = 0.005, rshell = -1;
3727     static real        maxnhb   = 0, fit_start = 1, fit_end = 60, temp = 298.15, smooth_tail_start = -1, D = -1;
3728     static gmx_bool    bNitAcc  = TRUE, bDA = TRUE, bMerge = TRUE;
3729     static int         nDump    = 0, nFitPoints = 100;
3730     static int         nThreads = 0, nBalExp = 4;
3731
3732     static gmx_bool    bContact     = FALSE, bBallistic = FALSE, bGemFit = FALSE;
3733     static real        logAfterTime = 10, gemBallistic = 0.2; /* ps */
3734     static const char *NNtype[]     = {NULL, "none", "binary", "oneOverR3", "dipole", NULL};
3735
3736     /* options */
3737     t_pargs     pa [] = {
3738         { "-a",    FALSE,  etREAL, {&acut},
3739           "Cutoff angle (degrees, Hydrogen - Donor - Acceptor)" },
3740         { "-r",    FALSE,  etREAL, {&rcut},
3741           "Cutoff radius (nm, X - Acceptor, see next option)" },
3742         { "-da",   FALSE,  etBOOL, {&bDA},
3743           "Use distance Donor-Acceptor (if TRUE) or Hydrogen-Acceptor (FALSE)" },
3744         { "-r2",   FALSE,  etREAL, {&r2cut},
3745           "Second cutoff radius. Mainly useful with [TT]-contact[tt] and [TT]-ac[tt]"},
3746         { "-abin", FALSE,  etREAL, {&abin},
3747           "Binwidth angle distribution (degrees)" },
3748         { "-rbin", FALSE,  etREAL, {&rbin},
3749           "Binwidth distance distribution (nm)" },
3750         { "-nitacc", FALSE, etBOOL, {&bNitAcc},
3751           "Regard nitrogen atoms as acceptors" },
3752         { "-contact", FALSE, etBOOL, {&bContact},
3753           "Do not look for hydrogen bonds, but merely for contacts within the cut-off distance" },
3754         { "-shell", FALSE, etREAL, {&rshell},
3755           "when > 0, only calculate hydrogen bonds within # nm shell around "
3756           "one particle" },
3757         { "-fitstart", FALSE, etREAL, {&fit_start},
3758           "Time (ps) from which to start fitting the correlation functions in order to obtain the forward and backward rate constants for HB breaking and formation. With [TT]-gemfit[tt] we suggest [TT]-fitstart 0[tt]" },
3759         { "-fitstart", FALSE, etREAL, {&fit_start},
3760           "Time (ps) to which to stop fitting the correlation functions in order to obtain the forward and backward rate constants for HB breaking and formation (only with [TT]-gemfit[tt])" },
3761         { "-temp",  FALSE, etREAL, {&temp},
3762           "Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" },
3763         { "-smooth", FALSE, etREAL, {&smooth_tail_start},
3764           "If >= 0, the tail of the ACF will be smoothed by fitting it to an exponential function: y = A exp(-x/[GRK]tau[grk])" },
3765         { "-dump",  FALSE, etINT, {&nDump},
3766           "Dump the first N hydrogen bond ACFs in a single [TT].xvg[tt] file for debugging" },
3767         { "-max_hb", FALSE, etREAL, {&maxnhb},
3768           "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" },
3769         { "-merge", FALSE, etBOOL, {&bMerge},
3770           "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." },
3771         { "-geminate", FALSE, etENUM, {gemType},
3772           "Use reversible geminate recombination for the kinetics/thermodynamics calclations. See Markovitch et al., J. Chem. Phys 129, 084505 (2008) for details."},
3773         { "-diff", FALSE, etREAL, {&D},
3774           "Dffusion coefficient to use in the reversible geminate recombination kinetic model. If negative, then it will be fitted to the ACF along with ka and kd."},
3775 #ifdef GMX_OPENMP
3776         { "-nthreads", FALSE, etINT, {&nThreads},
3777           "Number of threads used for the parallel loop over autocorrelations. nThreads <= 0 means maximum number of threads. Requires linking with OpenMP. The number of threads is limited by the number of processors (before OpenMP v.3 ) or environment variable OMP_THREAD_LIMIT (OpenMP v.3)"},
3778 #endif
3779     };
3780     const char *bugs[] = {
3781         "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and therefore not available for the time being."
3782     };
3783     t_filenm    fnm[] = {
3784         { efTRX, "-f",   NULL,     ffREAD  },
3785         { efTPX, NULL,   NULL,     ffREAD  },
3786         { efNDX, NULL,   NULL,     ffOPTRD },
3787         /*    { efNDX, "-sel", "select", ffOPTRD },*/
3788         { efXVG, "-num", "hbnum",  ffWRITE },
3789         { efLOG, "-g",   "hbond",  ffOPTWR },
3790         { efXVG, "-ac",  "hbac",   ffOPTWR },
3791         { efXVG, "-dist", "hbdist", ffOPTWR },
3792         { efXVG, "-ang", "hbang",  ffOPTWR },
3793         { efXVG, "-hx",  "hbhelix", ffOPTWR },
3794         { efNDX, "-hbn", "hbond",  ffOPTWR },
3795         { efXPM, "-hbm", "hbmap",  ffOPTWR },
3796         { efXVG, "-don", "donor",  ffOPTWR },
3797         { efXVG, "-dan", "danum",  ffOPTWR },
3798         { efXVG, "-life", "hblife", ffOPTWR },
3799         { efXVG, "-nhbdist", "nhbdist", ffOPTWR }
3800
3801     };
3802 #define NFILE asize(fnm)
3803
3804     char                  hbmap [HB_NR] = { ' ',    'o',      '-',       '*' };
3805     const char           *hbdesc[HB_NR] = { "None", "Present", "Inserted", "Present & Inserted" };
3806     t_rgb                 hbrgb [HB_NR] = { {1, 1, 1}, {1, 0, 0},   {0, 0, 1},    {1, 0, 1} };
3807
3808     t_trxstatus          *status;
3809     int                   trrStatus = 1;
3810     t_topology            top;
3811     t_inputrec            ir;
3812     t_pargs              *ppa;
3813     int                   npargs, natoms, nframes = 0, shatom;
3814     int                  *isize;
3815     char                **grpnames;
3816     atom_id             **index;
3817     rvec                 *x, hbox;
3818     matrix                box;
3819     real                  t, ccut, dist = 0.0, ang = 0.0;
3820     double                max_nhb, aver_nhb, aver_dist;
3821     int                   h = 0, i = 0, j, k = 0, l, start, end, id, ja, ogrp, nsel;
3822     int                   xi, yi, zi, ai;
3823     int                   xj, yj, zj, aj, xjj, yjj, zjj;
3824     int                   xk, yk, zk, ak, xkk, ykk, zkk;
3825     gmx_bool              bSelected, bHBmap, bStop, bTwo, was, bBox, bTric;
3826     int                  *adist, *rdist, *aptr, *rprt;
3827     int                   grp, nabin, nrbin, bin, resdist, ihb;
3828     char                **leg;
3829     t_hbdata             *hb, *hbptr;
3830     FILE                 *fp, *fpins = NULL, *fpnhb = NULL;
3831     t_gridcell         ***grid;
3832     t_ncell              *icell, *jcell, *kcell;
3833     ivec                  ngrid;
3834     unsigned char        *datable;
3835     output_env_t          oenv;
3836     int                   gemmode, NN;
3837     PSTYPE                peri = 0;
3838     t_E                   E;
3839     int                   ii, jj, hh, actual_nThreads;
3840     int                   threadNr = 0;
3841     gmx_bool              bGem, bNN, bParallel;
3842     t_gemParams          *params = NULL;
3843     gmx_bool              bEdge_yjj, bEdge_xjj, bOMP;
3844
3845     t_hbdata            **p_hb    = NULL;                   /* one per thread, then merge after the frame loop */
3846     int                 **p_adist = NULL, **p_rdist = NULL; /* a histogram for each thread. */
3847
3848 #ifdef GMX_OPENMP
3849     bOMP = TRUE;
3850 #else
3851     bOMP = FALSE;
3852 #endif
3853
3854     npargs = asize(pa);
3855     ppa    = add_acf_pargs(&npargs, pa);
3856
3857     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE, NFILE, fnm, npargs,
3858                            ppa, asize(desc), desc, asize(bugs), bugs, &oenv))
3859     {
3860         return 0;
3861     }
3862
3863     /* NN-loop? If so, what estimator to use ?*/
3864     NN = 1;
3865     /* Outcommented for now DvdS 2010-07-13
3866        while (NN < NN_NR && gmx_strcasecmp(NNtype[0], NNtype[NN])!=0)
3867         NN++;
3868        if (NN == NN_NR)
3869         gmx_fatal(FARGS, "Invalid NN-loop type.");
3870      */
3871     bNN = FALSE;
3872     for (i = 2; bNN == FALSE && i < NN_NR; i++)
3873     {
3874         bNN = bNN || NN == i;
3875     }
3876
3877     if (NN > NN_NONE && bMerge)
3878     {
3879         bMerge = FALSE;
3880     }
3881
3882     /* geminate recombination? If so, which flavor? */
3883     gemmode = 1;
3884     while (gemmode < gemNR && gmx_strcasecmp(gemType[0], gemType[gemmode]) != 0)
3885     {
3886         gemmode++;
3887     }
3888     if (gemmode == gemNR)
3889     {
3890         gmx_fatal(FARGS, "Invalid recombination type.");
3891     }
3892
3893     bGem = FALSE;
3894     for (i = 2; bGem == FALSE && i < gemNR; i++)
3895     {
3896         bGem = bGem || gemmode == i;
3897     }
3898
3899     if (bGem)
3900     {
3901         printf("Geminate recombination: %s\n", gemType[gemmode]);
3902 #ifndef HAVE_LIBGSL
3903         printf("Note that some aspects of reversible geminate recombination won't work without gsl.\n");
3904 #endif
3905         if (bContact)
3906         {
3907             if (gemmode != gemDD)
3908             {
3909                 printf("Turning off -contact option...\n");
3910                 bContact = FALSE;
3911             }
3912         }
3913         else
3914         {
3915             if (gemmode == gemDD)
3916             {
3917                 printf("Turning on -contact option...\n");
3918                 bContact = TRUE;
3919             }
3920         }
3921         if (bMerge)
3922         {
3923             if (gemmode == gemAA)
3924             {
3925                 printf("Turning off -merge option...\n");
3926                 bMerge = FALSE;
3927             }
3928         }
3929         else
3930         {
3931             if (gemmode != gemAA)
3932             {
3933                 printf("Turning on -merge option...\n");
3934                 bMerge = TRUE;
3935             }
3936         }
3937     }
3938
3939     /* process input */
3940     bSelected = FALSE;
3941     ccut      = cos(acut*DEG2RAD);
3942
3943     if (bContact)
3944     {
3945         if (bSelected)
3946         {
3947             gmx_fatal(FARGS, "Can not analyze selected contacts.");
3948         }
3949         if (!bDA)
3950         {
3951             gmx_fatal(FARGS, "Can not analyze contact between H and A: turn off -noda");
3952         }
3953     }
3954
3955     /* Initiate main data structure! */
3956     bHBmap = (opt2bSet("-ac", NFILE, fnm) ||
3957               opt2bSet("-life", NFILE, fnm) ||
3958               opt2bSet("-hbn", NFILE, fnm) ||
3959               opt2bSet("-hbm", NFILE, fnm) ||
3960               bGem);
3961
3962     if (opt2bSet("-nhbdist", NFILE, fnm))
3963     {
3964         const char *leg[MAXHH+1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" };
3965         fpnhb = xvgropen(opt2fn("-nhbdist", NFILE, fnm),
3966                          "Number of donor-H with N HBs", output_env_get_xvgr_tlabel(oenv), "N", oenv);
3967         xvgr_legend(fpnhb, asize(leg), leg, oenv);
3968     }
3969
3970     hb = mk_hbdata(bHBmap, opt2bSet("-dan", NFILE, fnm), bMerge || bContact, bGem, gemmode);
3971
3972     /* get topology */
3973     read_tpx_top(ftp2fn(efTPX, NFILE, fnm), &ir, box, &natoms, NULL, NULL, NULL, &top);
3974
3975     snew(grpnames, grNR);
3976     snew(index, grNR);
3977     snew(isize, grNR);
3978     /* Make Donor-Acceptor table */
3979     snew(datable, top.atoms.nr);
3980     gen_datable(index[0], isize[0], datable, top.atoms.nr);
3981
3982     if (bSelected)
3983     {
3984         /* analyze selected hydrogen bonds */
3985         printf("Select group with selected atoms:\n");
3986         get_index(&(top.atoms), opt2fn("-sel", NFILE, fnm),
3987                   1, &nsel, index, grpnames);
3988         if (nsel % 3)
3989         {
3990             gmx_fatal(FARGS, "Number of atoms in group '%s' not a multiple of 3\n"
3991                       "and therefore cannot contain triplets of "
3992                       "Donor-Hydrogen-Acceptor", grpnames[0]);
3993         }
3994         bTwo = FALSE;
3995
3996         for (i = 0; (i < nsel); i += 3)
3997         {
3998             int dd = index[0][i];
3999             int aa = index[0][i+2];
4000             /* int */ hh = index[0][i+1];
4001             add_dh (&hb->d, dd, hh, i, datable);
4002             add_acc(&hb->a, aa, i);
4003             /* Should this be here ? */
4004             snew(hb->d.dptr, top.atoms.nr);
4005             snew(hb->a.aptr, top.atoms.nr);
4006             add_hbond(hb, dd, aa, hh, gr0, gr0, 0, bMerge, 0, bContact, peri);
4007         }
4008         printf("Analyzing %d selected hydrogen bonds from '%s'\n",
4009                isize[0], grpnames[0]);
4010     }
4011     else
4012     {
4013         /* analyze all hydrogen bonds: get group(s) */
4014         printf("Specify 2 groups to analyze:\n");
4015         get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
4016                   2, isize, index, grpnames);
4017
4018         /* check if we have two identical or two non-overlapping groups */
4019         bTwo = isize[0] != isize[1];
4020         for (i = 0; (i < isize[0]) && !bTwo; i++)
4021         {
4022             bTwo = index[0][i] != index[1][i];
4023         }
4024         if (bTwo)
4025         {
4026             printf("Checking for overlap in atoms between %s and %s\n",
4027                    grpnames[0], grpnames[1]);
4028             for (i = 0; i < isize[1]; i++)
4029             {
4030                 if (ISINGRP(datable[index[1][i]]))
4031                 {
4032                     gmx_fatal(FARGS, "Partial overlap between groups '%s' and '%s'",
4033                               grpnames[0], grpnames[1]);
4034                 }
4035             }
4036             /*
4037                printf("Checking for overlap in atoms between %s and %s\n",
4038                grpnames[0],grpnames[1]);
4039                for (i=0; i<isize[0]; i++)
4040                for (j=0; j<isize[1]; j++)
4041                if (index[0][i] == index[1][j])
4042                gmx_fatal(FARGS,"Partial overlap between groups '%s' and '%s'",
4043                grpnames[0],grpnames[1]);
4044              */
4045         }
4046         if (bTwo)
4047         {
4048             printf("Calculating %s "
4049                    "between %s (%d atoms) and %s (%d atoms)\n",
4050                    bContact ? "contacts" : "hydrogen bonds",
4051                    grpnames[0], isize[0], grpnames[1], isize[1]);
4052         }
4053         else
4054         {
4055             fprintf(stderr, "Calculating %s in %s (%d atoms)\n",
4056                     bContact ? "contacts" : "hydrogen bonds", grpnames[0], isize[0]);
4057         }
4058     }
4059     sfree(datable);
4060
4061     /* search donors and acceptors in groups */
4062     snew(datable, top.atoms.nr);
4063     for (i = 0; (i < grNR); i++)
4064     {
4065         if ( ((i == gr0) && !bSelected ) ||
4066              ((i == gr1) && bTwo ))
4067         {
4068             gen_datable(index[i], isize[i], datable, top.atoms.nr);
4069             if (bContact)
4070             {
4071                 search_acceptors(&top, isize[i], index[i], &hb->a, i,
4072                                  bNitAcc, TRUE, (bTwo && (i == gr0)) || !bTwo, datable);
4073                 search_donors   (&top, isize[i], index[i], &hb->d, i,
4074                                  TRUE, (bTwo && (i == gr1)) || !bTwo, datable);
4075             }
4076             else
4077             {
4078                 search_acceptors(&top, isize[i], index[i], &hb->a, i, bNitAcc, FALSE, TRUE, datable);
4079                 search_donors   (&top, isize[i], index[i], &hb->d, i, FALSE, TRUE, datable);
4080             }
4081             if (bTwo)
4082             {
4083                 clear_datable_grp(datable, top.atoms.nr);
4084             }
4085         }
4086     }
4087     sfree(datable);
4088     printf("Found %d donors and %d acceptors\n", hb->d.nrd, hb->a.nra);
4089     /*if (bSelected)
4090        snew(donors[gr0D], dons[gr0D].nrd);*/
4091
4092     if (bHBmap)
4093     {
4094         printf("Making hbmap structure...");
4095         /* Generate hbond data structure */
4096         mk_hbmap(hb);
4097         printf("done.\n");
4098     }
4099
4100 #ifdef HAVE_NN_LOOPS
4101     if (bNN)
4102     {
4103         mk_hbEmap(hb, 0);
4104     }
4105 #endif
4106
4107     if (bGem)
4108     {
4109         printf("Making per structure...");
4110         /* Generate hbond data structure */
4111         mk_per(hb);
4112         printf("done.\n");
4113     }
4114
4115     /* check input */
4116     bStop = FALSE;
4117     if (hb->d.nrd + hb->a.nra == 0)
4118     {
4119         printf("No Donors or Acceptors found\n");
4120         bStop = TRUE;
4121     }
4122     if (!bStop)
4123     {
4124         if (hb->d.nrd == 0)
4125         {
4126             printf("No Donors found\n");
4127             bStop = TRUE;
4128         }
4129         if (hb->a.nra == 0)
4130         {
4131             printf("No Acceptors found\n");
4132             bStop = TRUE;
4133         }
4134     }
4135     if (bStop)
4136     {
4137         gmx_fatal(FARGS, "Nothing to be done");
4138     }
4139
4140     shatom = 0;
4141     if (rshell > 0)
4142     {
4143         int      shisz;
4144         atom_id *shidx;
4145         char    *shgrpnm;
4146         /* get index group with atom for shell */
4147         do
4148         {
4149             printf("Select atom for shell (1 atom):\n");
4150             get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
4151                       1, &shisz, &shidx, &shgrpnm);
4152             if (shisz != 1)
4153             {
4154                 printf("group contains %d atoms, should be 1 (one)\n", shisz);
4155             }
4156         }
4157         while (shisz != 1);
4158         shatom = shidx[0];
4159         printf("Will calculate hydrogen bonds within a shell "
4160                "of %g nm around atom %i\n", rshell, shatom+1);
4161     }
4162
4163     /* Analyze trajectory */
4164     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
4165     if (natoms > top.atoms.nr)
4166     {
4167         gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)",
4168                   top.atoms.nr, natoms);
4169     }
4170
4171     bBox  = ir.ePBC != epbcNONE;
4172     grid  = init_grid(bBox, box, (rcut > r2cut) ? rcut : r2cut, ngrid);
4173     nabin = acut/abin;
4174     nrbin = rcut/rbin;
4175     snew(adist, nabin+1);
4176     snew(rdist, nrbin+1);
4177
4178     if (bGem && !bBox)
4179     {
4180         gmx_fatal(FARGS, "Can't do geminate recombination without periodic box.");
4181     }
4182
4183     bParallel = FALSE;
4184
4185 #ifndef GMX_OPENMP
4186 #define __ADIST adist
4187 #define __RDIST rdist
4188 #define __HBDATA hb
4189 #else /* GMX_OPENMP ================================================== \
4190        * Set up the OpenMP stuff,                                       |
4191        * like the number of threads and such                            |
4192        * Also start the parallel loop.                                  |
4193        */
4194 #define __ADIST p_adist[threadNr]
4195 #define __RDIST p_rdist[threadNr]
4196 #define __HBDATA p_hb[threadNr]
4197 #endif
4198     if (bOMP)
4199     {
4200         bParallel = !bSelected;
4201
4202         if (bParallel)
4203         {
4204             actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
4205
4206             gmx_omp_set_num_threads(actual_nThreads);
4207             printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads);
4208             fflush(stdout);
4209         }
4210         else
4211         {
4212             actual_nThreads = 1;
4213         }
4214
4215         snew(p_hb,    actual_nThreads);
4216         snew(p_adist, actual_nThreads);
4217         snew(p_rdist, actual_nThreads);
4218         for (i = 0; i < actual_nThreads; i++)
4219         {
4220             snew(p_hb[i], 1);
4221             snew(p_adist[i], nabin+1);
4222             snew(p_rdist[i], nrbin+1);
4223
4224             p_hb[i]->max_frames = 0;
4225             p_hb[i]->nhb        = NULL;
4226             p_hb[i]->ndist      = NULL;
4227             p_hb[i]->n_bound    = NULL;
4228             p_hb[i]->time       = NULL;
4229             p_hb[i]->nhx        = NULL;
4230
4231             p_hb[i]->bHBmap     = hb->bHBmap;
4232             p_hb[i]->bDAnr      = hb->bDAnr;
4233             p_hb[i]->bGem       = hb->bGem;
4234             p_hb[i]->wordlen    = hb->wordlen;
4235             p_hb[i]->nframes    = hb->nframes;
4236             p_hb[i]->maxhydro   = hb->maxhydro;
4237             p_hb[i]->danr       = hb->danr;
4238             p_hb[i]->d          = hb->d;
4239             p_hb[i]->a          = hb->a;
4240             p_hb[i]->hbmap      = hb->hbmap;
4241             p_hb[i]->time       = hb->time; /* This may need re-syncing at every frame. */
4242             p_hb[i]->per        = hb->per;
4243
4244 #ifdef HAVE_NN_LOOPS
4245             p_hb[i]->hbE = hb->hbE;
4246 #endif
4247
4248             p_hb[i]->nrhb   = 0;
4249             p_hb[i]->nrdist = 0;
4250         }
4251     }
4252
4253     /* Make a thread pool here,
4254      * instead of forking anew at every frame. */
4255
4256 #pragma omp parallel \
4257     firstprivate(i) \
4258     private(j, h, ii, jj, hh, E, \
4259     xi, yi, zi, xj, yj, zj, threadNr, \
4260     dist, ang, peri, icell, jcell, \
4261     grp, ogrp, ai, aj, xjj, yjj, zjj, \
4262     xk, yk, zk, ihb, id,  resdist, \
4263     xkk, ykk, zkk, kcell, ak, k, bTric, \
4264     bEdge_xjj, bEdge_yjj) \
4265     default(shared)
4266     {    /* Start of parallel region */
4267         threadNr = gmx_omp_get_thread_num();
4268
4269         do
4270         {
4271
4272             bTric = bBox && TRICLINIC(box);
4273
4274             if (bOMP)
4275             {
4276                 sync_hbdata(p_hb[threadNr], nframes);
4277             }
4278 #pragma omp single
4279             {
4280                 build_grid(hb, x, x[shatom], bBox, box, hbox, (rcut > r2cut) ? rcut : r2cut,
4281                            rshell, ngrid, grid);
4282                 reset_nhbonds(&(hb->d));
4283
4284                 if (debug && bDebug)
4285                 {
4286                     dump_grid(debug, ngrid, grid);
4287                 }
4288
4289                 add_frames(hb, nframes);
4290                 init_hbframe(hb, nframes, output_env_conv_time(oenv, t));
4291
4292                 if (hb->bDAnr)
4293                 {
4294                     count_da_grid(ngrid, grid, hb->danr[nframes]);
4295                 }
4296             } /* omp single */
4297
4298             if (bOMP)
4299             {
4300                 p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */
4301             }
4302
4303             if (bNN)
4304             {
4305 #ifdef HAVE_NN_LOOPS /* Unlock this feature when testing */
4306                 /* Loop over all atom pairs and estimate interaction energy */
4307
4308 #pragma omp single
4309                 {
4310                     addFramesNN(hb, nframes);
4311                 }
4312
4313 #pragma omp barrier
4314 #pragma omp for schedule(dynamic)
4315                 for (i = 0; i < hb->d.nrd; i++)
4316                 {
4317                     for (j = 0; j < hb->a.nra; j++)
4318                     {
4319                         for (h = 0;
4320                              h < (bContact ? 1 : hb->d.nhydro[i]);
4321                              h++)
4322                         {
4323                             if (i == hb->d.nrd || j == hb->a.nra)
4324                             {
4325                                 gmx_fatal(FARGS, "out of bounds");
4326                             }
4327
4328                             /* Get the real atom ids */
4329                             ii = hb->d.don[i];
4330                             jj = hb->a.acc[j];
4331                             hh = hb->d.hydro[i][h];
4332
4333                             /* Estimate the energy from the geometry */
4334                             E = calcHbEnergy(ii, jj, hh, x, NN, box, hbox, &(hb->d));
4335                             /* Store the energy */
4336                             storeHbEnergy(hb, i, j, h, E, nframes);
4337                         }
4338                     }
4339                 }
4340 #endif        /* HAVE_NN_LOOPS */
4341             } /* if (bNN)*/
4342             else
4343             {
4344                 if (bSelected)
4345                 {
4346
4347 #pragma omp single
4348                     {
4349                         /* Do not parallelize this just yet. */
4350                         /* int ii; */
4351                         for (ii = 0; (ii < nsel); ii++)
4352                         {
4353                             int dd = index[0][i];
4354                             int aa = index[0][i+2];
4355                             /* int */ hh = index[0][i+1];
4356                             ihb          = is_hbond(hb, ii, ii, dd, aa, rcut, r2cut, ccut, x, bBox, box,
4357                                                     hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri);
4358
4359                             if (ihb)
4360                             {
4361                                 /* add to index if not already there */
4362                                 /* Add a hbond */
4363                                 add_hbond(hb, dd, aa, hh, ii, ii, nframes, bMerge, ihb, bContact, peri);
4364                             }
4365                         }
4366                     } /* omp single */
4367                 }     /* if (bSelected) */
4368                 else
4369                 {
4370
4371 #pragma omp single
4372                     {
4373                         if (bGem)
4374                         {
4375                             calcBoxProjection(box, hb->per->P);
4376                         }
4377
4378                         /* loop over all gridcells (xi,yi,zi)      */
4379                         /* Removed confusing macro, DvdS 27/12/98  */
4380
4381                     }
4382                     /* The outer grid loop will have to do for now. */
4383 #pragma omp for schedule(dynamic)
4384                     for (xi = 0; xi < ngrid[XX]; xi++)
4385                     {
4386                         for (yi = 0; (yi < ngrid[YY]); yi++)
4387                         {
4388                             for (zi = 0; (zi < ngrid[ZZ]); zi++)
4389                             {
4390
4391                                 /* loop over donor groups gr0 (always) and gr1 (if necessary) */
4392                                 for (grp = gr0; (grp <= (bTwo ? gr1 : gr0)); grp++)
4393                                 {
4394                                     icell = &(grid[zi][yi][xi].d[grp]);
4395
4396                                     if (bTwo)
4397                                     {
4398                                         ogrp = 1-grp;
4399                                     }
4400                                     else
4401                                     {
4402                                         ogrp = grp;
4403                                     }
4404
4405                                     /* loop over all hydrogen atoms from group (grp)
4406                                      * in this gridcell (icell)
4407                                      */
4408                                     for (ai = 0; (ai < icell->nr); ai++)
4409                                     {
4410                                         i  = icell->atoms[ai];
4411
4412                                         /* loop over all adjacent gridcells (xj,yj,zj) */
4413                                         for (zjj = grid_loop_begin(ngrid[ZZ], zi, bTric, FALSE);
4414                                              zjj <= grid_loop_end(ngrid[ZZ], zi, bTric, FALSE);
4415                                              zjj++)
4416                                         {
4417                                             zj        = grid_mod(zjj, ngrid[ZZ]);
4418                                             bEdge_yjj = (zj == 0) || (zj == ngrid[ZZ] - 1);
4419                                             for (yjj = grid_loop_begin(ngrid[YY], yi, bTric, bEdge_yjj);
4420                                                  yjj <= grid_loop_end(ngrid[YY], yi, bTric, bEdge_yjj);
4421                                                  yjj++)
4422                                             {
4423                                                 yj        = grid_mod(yjj, ngrid[YY]);
4424                                                 bEdge_xjj =
4425                                                     (yj == 0) || (yj == ngrid[YY] - 1) ||
4426                                                     (zj == 0) || (zj == ngrid[ZZ] - 1);
4427                                                 for (xjj = grid_loop_begin(ngrid[XX], xi, bTric, bEdge_xjj);
4428                                                      xjj <= grid_loop_end(ngrid[XX], xi, bTric, bEdge_xjj);
4429                                                      xjj++)
4430                                                 {
4431                                                     xj    = grid_mod(xjj, ngrid[XX]);
4432                                                     jcell = &(grid[zj][yj][xj].a[ogrp]);
4433                                                     /* loop over acceptor atoms from other group (ogrp)
4434                                                      * in this adjacent gridcell (jcell)
4435                                                      */
4436                                                     for (aj = 0; (aj < jcell->nr); aj++)
4437                                                     {
4438                                                         j = jcell->atoms[aj];
4439
4440                                                         /* check if this once was a h-bond */
4441                                                         peri = -1;
4442                                                         ihb  = is_hbond(__HBDATA, grp, ogrp, i, j, rcut, r2cut, ccut, x, bBox, box,
4443                                                                         hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri);
4444
4445                                                         if (ihb)
4446                                                         {
4447                                                             /* add to index if not already there */
4448                                                             /* Add a hbond */
4449                                                             add_hbond(__HBDATA, i, j, h, grp, ogrp, nframes, bMerge, ihb, bContact, peri);
4450
4451                                                             /* make angle and distance distributions */
4452                                                             if (ihb == hbHB && !bContact)
4453                                                             {
4454                                                                 if (dist > rcut)
4455                                                                 {
4456                                                                     gmx_fatal(FARGS, "distance is higher than what is allowed for an hbond: %f", dist);
4457                                                                 }
4458                                                                 ang *= RAD2DEG;
4459                                                                 __ADIST[(int)( ang/abin)]++;
4460                                                                 __RDIST[(int)(dist/rbin)]++;
4461                                                                 if (!bTwo)
4462                                                                 {
4463                                                                     int id, ia;
4464                                                                     if ((id = donor_index(&hb->d, grp, i)) == NOTSET)
4465                                                                     {
4466                                                                         gmx_fatal(FARGS, "Invalid donor %d", i);
4467                                                                     }
4468                                                                     if ((ia = acceptor_index(&hb->a, ogrp, j)) == NOTSET)
4469                                                                     {
4470                                                                         gmx_fatal(FARGS, "Invalid acceptor %d", j);
4471                                                                     }
4472                                                                     resdist = abs(top.atoms.atom[i].resind-
4473                                                                                   top.atoms.atom[j].resind);
4474                                                                     if (resdist >= max_hx)
4475                                                                     {
4476                                                                         resdist = max_hx-1;
4477                                                                     }
4478                                                                     __HBDATA->nhx[nframes][resdist]++;
4479                                                                 }
4480                                                             }
4481
4482                                                         }
4483                                                     } /* for aj  */
4484                                                 }     /* for xjj */
4485                                             }         /* for yjj */
4486                                         }             /* for zjj */
4487                                     }                 /* for ai  */
4488                                 }                     /* for grp */
4489                             }                         /* for xi,yi,zi */
4490                         }
4491                     }
4492                 } /* if (bSelected) {...} else */
4493
4494
4495                 /* Better wait for all threads to finnish using x[] before updating it. */
4496                 k = nframes;
4497 #pragma omp barrier
4498 #pragma omp critical
4499                 {
4500                     /* Sum up histograms and counts from p_hb[] into hb */
4501                     if (bOMP)
4502                     {
4503                         hb->nhb[k]   += p_hb[threadNr]->nhb[k];
4504                         hb->ndist[k] += p_hb[threadNr]->ndist[k];
4505                         for (j = 0; j < max_hx; j++)
4506                         {
4507                             hb->nhx[k][j]  += p_hb[threadNr]->nhx[k][j];
4508                         }
4509                     }
4510                 }
4511
4512                 /* Here are a handful of single constructs
4513                  * to share the workload a bit. The most
4514                  * important one is of course the last one,
4515                  * where there's a potential bottleneck in form
4516                  * of slow I/O.                    */
4517 #pragma omp barrier
4518 #pragma omp single
4519                 {
4520                     if (hb != NULL)
4521                     {
4522                         analyse_donor_props(opt2fn_null("-don", NFILE, fnm), hb, k, t, oenv);
4523                     }
4524                 }
4525
4526 #pragma omp single
4527                 {
4528                     if (fpnhb)
4529                     {
4530                         do_nhb_dist(fpnhb, hb, t);
4531                     }
4532                 }
4533             } /* if (bNN) {...} else  +   */
4534
4535 #pragma omp single
4536             {
4537                 trrStatus = (read_next_x(oenv, status, &t, x, box));
4538                 nframes++;
4539             }
4540
4541 #pragma omp barrier
4542         }
4543         while (trrStatus);
4544
4545         if (bOMP)
4546         {
4547 #pragma omp critical
4548             {
4549                 hb->nrhb   += p_hb[threadNr]->nrhb;
4550                 hb->nrdist += p_hb[threadNr]->nrdist;
4551             }
4552             /* Free parallel datastructures */
4553             sfree(p_hb[threadNr]->nhb);
4554             sfree(p_hb[threadNr]->ndist);
4555             sfree(p_hb[threadNr]->nhx);
4556
4557 #pragma omp for
4558             for (i = 0; i < nabin; i++)
4559             {
4560                 for (j = 0; j < actual_nThreads; j++)
4561                 {
4562
4563                     adist[i] += p_adist[j][i];
4564                 }
4565             }
4566 #pragma omp for
4567             for (i = 0; i <= nrbin; i++)
4568             {
4569                 for (j = 0; j < actual_nThreads; j++)
4570                 {
4571                     rdist[i] += p_rdist[j][i];
4572                 }
4573             }
4574
4575             sfree(p_adist[threadNr]);
4576             sfree(p_rdist[threadNr]);
4577         }
4578     } /* End of parallel region */
4579     if (bOMP)
4580     {
4581         sfree(p_adist);
4582         sfree(p_rdist);
4583     }
4584
4585     if (nframes < 2 && (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm)))
4586     {
4587         gmx_fatal(FARGS, "Cannot calculate autocorrelation of life times with less than two frames");
4588     }
4589
4590     free_grid(ngrid, &grid);
4591
4592     close_trj(status);
4593     if (fpnhb)
4594     {
4595         ffclose(fpnhb);
4596     }
4597
4598     /* Compute maximum possible number of different hbonds */
4599     if (maxnhb > 0)
4600     {
4601         max_nhb = maxnhb;
4602     }
4603     else
4604     {
4605         max_nhb = 0.5*(hb->d.nrd*hb->a.nra);
4606     }
4607     /* Added support for -contact below.
4608      * - Erik Marklund, May 29-31, 2006 */
4609     /* Changed contact code.
4610      * - Erik Marklund, June 29, 2006 */
4611     if (bHBmap && !bNN)
4612     {
4613         if (hb->nrhb == 0)
4614         {
4615             printf("No %s found!!\n", bContact ? "contacts" : "hydrogen bonds");
4616         }
4617         else
4618         {
4619             printf("Found %d different %s in trajectory\n"
4620                    "Found %d different atom-pairs within %s distance\n",
4621                    hb->nrhb, bContact ? "contacts" : "hydrogen bonds",
4622                    hb->nrdist, (r2cut > 0) ? "second cut-off" : "hydrogen bonding");
4623
4624             /*Control the pHist.*/
4625
4626             if (bMerge)
4627             {
4628                 merge_hb(hb, bTwo, bContact);
4629             }
4630
4631             if (opt2bSet("-hbn", NFILE, fnm))
4632             {
4633                 dump_hbmap(hb, NFILE, fnm, bTwo, bContact, isize, index, grpnames, &top.atoms);
4634             }
4635
4636             /* Moved the call to merge_hb() to a line BEFORE dump_hbmap
4637              * to make the -hbn and -hmb output match eachother.
4638              * - Erik Marklund, May 30, 2006 */
4639         }
4640     }
4641     /* Print out number of hbonds and distances */
4642     aver_nhb  = 0;
4643     aver_dist = 0;
4644     fp        = xvgropen(opt2fn("-num", NFILE, fnm), bContact ? "Contacts" :
4645                          "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
4646     snew(leg, 2);
4647     snew(leg[0], STRLEN);
4648     snew(leg[1], STRLEN);
4649     sprintf(leg[0], "%s", bContact ? "Contacts" : "Hydrogen bonds");
4650     sprintf(leg[1], "Pairs within %g nm", (r2cut > 0) ? r2cut : rcut);
4651     xvgr_legend(fp, 2, (const char**)leg, oenv);
4652     sfree(leg[1]);
4653     sfree(leg[0]);
4654     sfree(leg);
4655     for (i = 0; (i < nframes); i++)
4656     {
4657         fprintf(fp, "%10g  %10d  %10d\n", hb->time[i], hb->nhb[i], hb->ndist[i]);
4658         aver_nhb  += hb->nhb[i];
4659         aver_dist += hb->ndist[i];
4660     }
4661     ffclose(fp);
4662     aver_nhb  /= nframes;
4663     aver_dist /= nframes;
4664     /* Print HB distance distribution */
4665     if (opt2bSet("-dist", NFILE, fnm))
4666     {
4667         long sum;
4668
4669         sum = 0;
4670         for (i = 0; i < nrbin; i++)
4671         {
4672             sum += rdist[i];
4673         }
4674
4675         fp = xvgropen(opt2fn("-dist", NFILE, fnm),
4676                       "Hydrogen Bond Distribution",
4677                       bDA ?
4678                       "Donor - Acceptor Distance (nm)" :
4679                       "Hydrogen - Acceptor Distance (nm)", "", oenv);
4680         for (i = 0; i < nrbin; i++)
4681         {
4682             fprintf(fp, "%10g %10g\n", (i+0.5)*rbin, rdist[i]/(rbin*(real)sum));
4683         }
4684         ffclose(fp);
4685     }
4686
4687     /* Print HB angle distribution */
4688     if (opt2bSet("-ang", NFILE, fnm))
4689     {
4690         long sum;
4691
4692         sum = 0;
4693         for (i = 0; i < nabin; i++)
4694         {
4695             sum += adist[i];
4696         }
4697
4698         fp = xvgropen(opt2fn("-ang", NFILE, fnm),
4699                       "Hydrogen Bond Distribution",
4700                       "Hydrogen - Donor - Acceptor Angle (\\SO\\N)", "", oenv);
4701         for (i = 0; i < nabin; i++)
4702         {
4703             fprintf(fp, "%10g %10g\n", (i+0.5)*abin, adist[i]/(abin*(real)sum));
4704         }
4705         ffclose(fp);
4706     }
4707
4708     /* Print HB in alpha-helix */
4709     if (opt2bSet("-hx", NFILE, fnm))
4710     {
4711         fp = xvgropen(opt2fn("-hx", NFILE, fnm),
4712                       "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
4713         xvgr_legend(fp, NRHXTYPES, hxtypenames, oenv);
4714         for (i = 0; i < nframes; i++)
4715         {
4716             fprintf(fp, "%10g", hb->time[i]);
4717             for (j = 0; j < max_hx; j++)
4718             {
4719                 fprintf(fp, " %6d", hb->nhx[i][j]);
4720             }
4721             fprintf(fp, "\n");
4722         }
4723         ffclose(fp);
4724     }
4725     if (!bNN)
4726     {
4727         printf("Average number of %s per timeframe %.3f out of %g possible\n",
4728                bContact ? "contacts" : "hbonds",
4729                bContact ? aver_dist : aver_nhb, max_nhb);
4730     }
4731
4732     /* Do Autocorrelation etc. */
4733     if (hb->bHBmap)
4734     {
4735         /*
4736            Added support for -contact in ac and hbm calculations below.
4737            - Erik Marklund, May 29, 2006
4738          */
4739         ivec itmp;
4740         rvec rtmp;
4741         if (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm))
4742         {
4743             please_cite(stdout, "Spoel2006b");
4744         }
4745         if (opt2bSet("-ac", NFILE, fnm))
4746         {
4747             char *gemstring = NULL;
4748
4749             if (bGem || bNN)
4750             {
4751                 params = init_gemParams(rcut, D, hb->time, hb->nframes/2, nFitPoints, fit_start, fit_end,
4752                                         gemBallistic, nBalExp);
4753                 if (params == NULL)
4754                 {
4755                     gmx_fatal(FARGS, "Could not initiate t_gemParams params.");
4756                 }
4757             }
4758             gemstring = strdup(gemType[hb->per->gemtype]);
4759             do_hbac(opt2fn("-ac", NFILE, fnm), hb, nDump,
4760                     bMerge, bContact, fit_start, temp, r2cut > 0, smooth_tail_start, oenv,
4761                     gemstring, nThreads, NN, bBallistic, bGemFit);
4762         }
4763         if (opt2bSet("-life", NFILE, fnm))
4764         {
4765             do_hblife(opt2fn("-life", NFILE, fnm), hb, bMerge, bContact, oenv);
4766         }
4767         if (opt2bSet("-hbm", NFILE, fnm))
4768         {
4769             t_matrix mat;
4770             int      id, ia, hh, x, y;
4771
4772             if ((nframes > 0) && (hb->nrhb > 0))
4773             {
4774                 mat.nx = nframes;
4775                 mat.ny = hb->nrhb;
4776
4777                 snew(mat.matrix, mat.nx);
4778                 for (x = 0; (x < mat.nx); x++)
4779                 {
4780                     snew(mat.matrix[x], mat.ny);
4781                 }
4782                 y = 0;
4783                 for (id = 0; (id < hb->d.nrd); id++)
4784                 {
4785                     for (ia = 0; (ia < hb->a.nra); ia++)
4786                     {
4787                         for (hh = 0; (hh < hb->maxhydro); hh++)
4788                         {
4789                             if (hb->hbmap[id][ia])
4790                             {
4791                                 if (ISHB(hb->hbmap[id][ia]->history[hh]))
4792                                 {
4793                                     /* Changed '<' into '<=' in the for-statement below.
4794                                      * It fixed the previously undiscovered bug that caused
4795                                      * the last occurance of an hbond/contact to not be
4796                                      * set in mat.matrix. Have a look at any old -hbm-output
4797                                      * and you will notice that the last column is allways empty.
4798                                      * - Erik Marklund May 30, 2006
4799                                      */
4800                                     for (x = 0; (x <= hb->hbmap[id][ia]->nframes); x++)
4801                                     {
4802                                         int nn0 = hb->hbmap[id][ia]->n0;
4803                                         range_check(y, 0, mat.ny);
4804                                         mat.matrix[x+nn0][y] = is_hb(hb->hbmap[id][ia]->h[hh], x);
4805                                     }
4806                                     y++;
4807                                 }
4808                             }
4809                         }
4810                     }
4811                 }
4812                 mat.axis_x = hb->time;
4813                 snew(mat.axis_y, mat.ny);
4814                 for (j = 0; j < mat.ny; j++)
4815                 {
4816                     mat.axis_y[j] = j;
4817                 }
4818                 sprintf(mat.title, bContact ? "Contact Existence Map" :
4819                         "Hydrogen Bond Existence Map");
4820                 sprintf(mat.legend, bContact ? "Contacts" : "Hydrogen Bonds");
4821                 sprintf(mat.label_x, "%s", output_env_get_xvgr_tlabel(oenv));
4822                 sprintf(mat.label_y, bContact ? "Contact Index" : "Hydrogen Bond Index");
4823                 mat.bDiscrete = TRUE;
4824                 mat.nmap      = 2;
4825                 snew(mat.map, mat.nmap);
4826                 for (i = 0; i < mat.nmap; i++)
4827                 {
4828                     mat.map[i].code.c1 = hbmap[i];
4829                     mat.map[i].desc    = hbdesc[i];
4830                     mat.map[i].rgb     = hbrgb[i];
4831                 }
4832                 fp = opt2FILE("-hbm", NFILE, fnm, "w");
4833                 write_xpm_m(fp, mat);
4834                 ffclose(fp);
4835                 for (x = 0; x < mat.nx; x++)
4836                 {
4837                     sfree(mat.matrix[x]);
4838                 }
4839                 sfree(mat.axis_y);
4840                 sfree(mat.matrix);
4841                 sfree(mat.map);
4842             }
4843             else
4844             {
4845                 fprintf(stderr, "No hydrogen bonds/contacts found. No hydrogen bond map will be printed.\n");
4846             }
4847         }
4848     }
4849
4850     if (bGem)
4851     {
4852         fprintf(stderr, "There were %i periodic shifts\n", hb->per->nper);
4853         fprintf(stderr, "Freeing pHist for all donors...\n");
4854         for (i = 0; i < hb->d.nrd; i++)
4855         {
4856             fprintf(stderr, "\r%i", i);
4857             if (hb->per->pHist[i] != NULL)
4858             {
4859                 for (j = 0; j < hb->a.nra; j++)
4860                 {
4861                     clearPshift(&(hb->per->pHist[i][j]));
4862                 }
4863                 sfree(hb->per->pHist[i]);
4864             }
4865         }
4866         sfree(hb->per->pHist);
4867         sfree(hb->per->p2i);
4868         sfree(hb->per);
4869         fprintf(stderr, "...done.\n");
4870     }
4871
4872 #ifdef HAVE_NN_LOOPS
4873     if (bNN)
4874     {
4875         free_hbEmap(hb);
4876     }
4877 #endif
4878
4879     if (hb->bDAnr)
4880     {
4881         int    i, j, nleg;
4882         char **legnames;
4883         char   buf[STRLEN];
4884
4885 #define USE_THIS_GROUP(j) ( (j == gr0) || (bTwo && (j == gr1)) )
4886
4887         fp = xvgropen(opt2fn("-dan", NFILE, fnm),
4888                       "Donors and Acceptors", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
4889         nleg = (bTwo ? 2 : 1)*2;
4890         snew(legnames, nleg);
4891         i = 0;
4892         for (j = 0; j < grNR; j++)
4893         {
4894             if (USE_THIS_GROUP(j) )
4895             {
4896                 sprintf(buf, "Donors %s", grpnames[j]);
4897                 legnames[i++] = strdup(buf);
4898                 sprintf(buf, "Acceptors %s", grpnames[j]);
4899                 legnames[i++] = strdup(buf);
4900             }
4901         }
4902         if (i != nleg)
4903         {
4904             gmx_incons("number of legend entries");
4905         }
4906         xvgr_legend(fp, nleg, (const char**)legnames, oenv);
4907         for (i = 0; i < nframes; i++)
4908         {
4909             fprintf(fp, "%10g", hb->time[i]);
4910             for (j = 0; (j < grNR); j++)
4911             {
4912                 if (USE_THIS_GROUP(j) )
4913                 {
4914                     fprintf(fp, " %6d", hb->danr[i][j]);
4915                 }
4916             }
4917             fprintf(fp, "\n");
4918         }
4919         ffclose(fp);
4920     }
4921
4922     return 0;
4923 }