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 "gromacs/commandline/pargs.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 gmx_unused    *fp,
2453                                       t_luzar gmx_unused *tl,
2454                                       int gmx_unused      maxiter,
2455                                       real gmx_unused     tol)
2456 {
2457     fprintf(stderr, "This program needs the GNU scientific library to work.\n");
2458
2459     return -1;
2460 }
2461 static real quality_of_fit(real gmx_unused chi2, int gmx_unused N)
2462 {
2463     fprintf(stderr, "This program needs the GNU scientific library to work.\n");
2464
2465     return -1;
2466 }
2467
2468 #endif
2469
2470 static real compute_weighted_rates(int n, real t[], real ct[], real nt[],
2471                                    real kt[], real sigma_ct[], real sigma_nt[],
2472                                    real sigma_kt[], real *k, real *kp,
2473                                    real *sigma_k, real *sigma_kp,
2474                                    real fit_start)
2475 {
2476 #define NK 10
2477     int      i, j;
2478     t_luzar  tl;
2479     real     kkk = 0, kkp = 0, kk2 = 0, kp2 = 0, chi2;
2480
2481     *sigma_k  = 0;
2482     *sigma_kp = 0;
2483
2484     for (i = 0; (i < n); i++)
2485     {
2486         if (t[i] >= fit_start)
2487         {
2488             break;
2489         }
2490     }
2491     tl.n0       = i;
2492     tl.n1       = n;
2493     tl.nparams  = 2;
2494     tl.ndelta   = 1;
2495     tl.t        = t;
2496     tl.ct       = ct;
2497     tl.nt       = nt;
2498     tl.kt       = kt;
2499     tl.sigma_ct = sigma_ct;
2500     tl.sigma_nt = sigma_nt;
2501     tl.sigma_kt = sigma_kt;
2502     tl.kkk[0]   = *k;
2503     tl.kkk[1]   = *kp;
2504
2505     chi2      = optimize_luzar_parameters(debug, &tl, 1000, 1e-3);
2506     *k        = tl.kkk[0];
2507     *kp       = tl.kkk[1] = *kp;
2508     tl.ndelta = NK;
2509     for (j = 0; (j < NK); j++)
2510     {
2511         (void) optimize_luzar_parameters(debug, &tl, 1000, 1e-3);
2512         kkk += tl.kkk[0];
2513         kkp += tl.kkk[1];
2514         kk2 += sqr(tl.kkk[0]);
2515         kp2 += sqr(tl.kkk[1]);
2516         tl.n0++;
2517     }
2518     *sigma_k  = sqrt(kk2/NK - sqr(kkk/NK));
2519     *sigma_kp = sqrt(kp2/NK - sqr(kkp/NK));
2520
2521     return chi2;
2522 }
2523
2524 static void smooth_tail(int n, real t[], real c[], real sigma_c[], real start,
2525                         const output_env_t oenv)
2526 {
2527     FILE *fp;
2528     real  e_1, fitparm[4];
2529     int   i;
2530
2531     e_1 = exp(-1);
2532     for (i = 0; (i < n); i++)
2533     {
2534         if (c[i] < e_1)
2535         {
2536             break;
2537         }
2538     }
2539     if (i < n)
2540     {
2541         fitparm[0] = t[i];
2542     }
2543     else
2544     {
2545         fitparm[0] = 10;
2546     }
2547     fitparm[1] = 0.95;
2548     do_lmfit(n, c, sigma_c, 0, t, start, t[n-1], oenv, bDebugMode(), effnEXP2, fitparm, 0);
2549 }
2550
2551 void analyse_corr(int n, real t[], real ct[], real nt[], real kt[],
2552                   real sigma_ct[], real sigma_nt[], real sigma_kt[],
2553                   real fit_start, real temp, real smooth_tail_start,
2554                   const output_env_t oenv)
2555 {
2556     int        i0, i;
2557     real       k = 1, kp = 1, kow = 1;
2558     real       Q = 0, chi22, chi2, dg, dgp, tau_hb, dtau, tau_rlx, e_1, dt, sigma_k, sigma_kp, ddg;
2559     double     tmp, sn2 = 0, sc2 = 0, sk2 = 0, scn = 0, sck = 0, snk = 0;
2560     gmx_bool   bError = (sigma_ct != NULL) && (sigma_nt != NULL) && (sigma_kt != NULL);
2561
2562     if (smooth_tail_start >= 0)
2563     {
2564         smooth_tail(n, t, ct, sigma_ct, smooth_tail_start, oenv);
2565         smooth_tail(n, t, nt, sigma_nt, smooth_tail_start, oenv);
2566         smooth_tail(n, t, kt, sigma_kt, smooth_tail_start, oenv);
2567     }
2568     for (i0 = 0; (i0 < n-2) && ((t[i0]-t[0]) < fit_start); i0++)
2569     {
2570         ;
2571     }
2572     if (i0 < n-2)
2573     {
2574         for (i = i0; (i < n); i++)
2575         {
2576             sc2 += sqr(ct[i]);
2577             sn2 += sqr(nt[i]);
2578             sk2 += sqr(kt[i]);
2579             sck += ct[i]*kt[i];
2580             snk += nt[i]*kt[i];
2581             scn += ct[i]*nt[i];
2582         }
2583         printf("Hydrogen bond thermodynamics at T = %g K\n", temp);
2584         tmp = (sn2*sc2-sqr(scn));
2585         if ((tmp > 0) && (sn2 > 0))
2586         {
2587             k    = (sn2*sck-scn*snk)/tmp;
2588             kp   = (k*scn-snk)/sn2;
2589             if (bError)
2590             {
2591                 chi2 = 0;
2592                 for (i = i0; (i < n); i++)
2593                 {
2594                     chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
2595                 }
2596                 chi22 = compute_weighted_rates(n, t, ct, nt, kt, sigma_ct, sigma_nt,
2597                                                sigma_kt, &k, &kp,
2598                                                &sigma_k, &sigma_kp, fit_start);
2599                 Q   = quality_of_fit(chi2, 2);
2600                 ddg = BOLTZ*temp*sigma_k/k;
2601                 printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n",
2602                        chi2, Q);
2603                 printf("The Rate and Delta G are followed by an error estimate\n");
2604                 printf("----------------------------------------------------------\n"
2605                        "Type      Rate (1/ps)  Sigma Time (ps)  DG (kJ/mol)  Sigma\n");
2606                 printf("Forward    %10.3f %6.2f   %8.3f  %10.3f %6.2f\n",
2607                        k, sigma_k, 1/k, calc_dg(1/k, temp), ddg);
2608                 ddg = BOLTZ*temp*sigma_kp/kp;
2609                 printf("Backward   %10.3f %6.2f   %8.3f  %10.3f %6.2f\n",
2610                        kp, sigma_kp, 1/kp, calc_dg(1/kp, temp), ddg);
2611             }
2612             else
2613             {
2614                 chi2 = 0;
2615                 for (i = i0; (i < n); i++)
2616                 {
2617                     chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
2618                 }
2619                 printf("Fitting parameters chi^2 = %10g\nQ = %10g\n",
2620                        chi2, Q);
2621                 printf("--------------------------------------------------\n"
2622                        "Type      Rate (1/ps) Time (ps)  DG (kJ/mol)  Chi^2\n");
2623                 printf("Forward    %10.3f   %8.3f  %10.3f  %10g\n",
2624                        k, 1/k, calc_dg(1/k, temp), chi2);
2625                 printf("Backward   %10.3f   %8.3f  %10.3f\n",
2626                        kp, 1/kp, calc_dg(1/kp, temp));
2627             }
2628         }
2629         if (sc2 > 0)
2630         {
2631             kow  = 2*sck/sc2;
2632             printf("One-way    %10.3f   %s%8.3f  %10.3f\n",
2633                    kow, bError ? "       " : "", 1/kow, calc_dg(1/kow, temp));
2634         }
2635         else
2636         {
2637             printf(" - Numerical problems computing HB thermodynamics:\n"
2638                    "sc2 = %g  sn2 = %g  sk2 = %g sck = %g snk = %g scn = %g\n",
2639                    sc2, sn2, sk2, sck, snk, scn);
2640         }
2641         /* Determine integral of the correlation function */
2642         tau_hb = evaluate_integral(n, t, ct, NULL, (t[n-1]-t[0])/2, &dtau);
2643         printf("Integral   %10.3f   %s%8.3f  %10.3f\n", 1/tau_hb,
2644                bError ? "       " : "", tau_hb, calc_dg(tau_hb, temp));
2645         e_1 = exp(-1);
2646         for (i = 0; (i < n-2); i++)
2647         {
2648             if ((ct[i] > e_1) && (ct[i+1] <= e_1))
2649             {
2650                 break;
2651             }
2652         }
2653         if (i < n-2)
2654         {
2655             /* Determine tau_relax from linear interpolation */
2656             tau_rlx = t[i]-t[0] + (e_1-ct[i])*(t[i+1]-t[i])/(ct[i+1]-ct[i]);
2657             printf("Relaxation %10.3f   %8.3f  %s%10.3f\n", 1/tau_rlx,
2658                    tau_rlx, bError ? "       " : "",
2659                    calc_dg(tau_rlx, temp));
2660         }
2661     }
2662     else
2663     {
2664         printf("Correlation functions too short to compute thermodynamics\n");
2665     }
2666 }
2667
2668 void compute_derivative(int nn, real x[], real y[], real dydx[])
2669 {
2670     int j;
2671
2672     /* Compute k(t) = dc(t)/dt */
2673     for (j = 1; (j < nn-1); j++)
2674     {
2675         dydx[j] = (y[j+1]-y[j-1])/(x[j+1]-x[j-1]);
2676     }
2677     /* Extrapolate endpoints */
2678     dydx[0]    = 2*dydx[1]   -  dydx[2];
2679     dydx[nn-1] = 2*dydx[nn-2] - dydx[nn-3];
2680 }
2681
2682 static void parallel_print(int *data, int nThreads)
2683 {
2684     /* This prints the donors on which each tread is currently working. */
2685     int i;
2686
2687     fprintf(stderr, "\r");
2688     for (i = 0; i < nThreads; i++)
2689     {
2690         fprintf(stderr, "%-7i", data[i]);
2691     }
2692 }
2693
2694 static void normalizeACF(real *ct, real *gt, int nhb, int len)
2695 {
2696     real ct_fac, gt_fac;
2697     int  i;
2698
2699     /* Xu and Berne use the same normalization constant */
2700
2701     ct_fac = 1.0/ct[0];
2702     gt_fac = (nhb == 0) ? 0 : 1.0/(real)nhb;
2703
2704     printf("Normalization for c(t) = %g for gh(t) = %g\n", ct_fac, gt_fac);
2705     for (i = 0; i < len; i++)
2706     {
2707         ct[i] *= ct_fac;
2708         if (gt != NULL)
2709         {
2710             gt[i] *= gt_fac;
2711         }
2712     }
2713 }
2714
2715 /* Added argument bContact in do_hbac(...). Also
2716  * added support for -contact in the actual code.
2717  * - Erik Marklund, May 31, 2006 */
2718 /* Changed contact code and added argument R2
2719  * - Erik Marklund, June 29, 2006
2720  */
2721 static void do_hbac(const char *fn, t_hbdata *hb,
2722                     int nDump, gmx_bool bMerge, gmx_bool bContact, real fit_start,
2723                     real temp, gmx_bool R2, real smooth_tail_start, const output_env_t oenv,
2724                     const char *gemType, int nThreads,
2725                     const int NN, const gmx_bool bBallistic, const gmx_bool bGemFit)
2726 {
2727     FILE          *fp;
2728     int            i, j, k, m, n, o, nd, ihb, idist, n2, nn, iter, nSets;
2729     const char    *legNN[]   = {
2730         "Ac(t)",
2731         "Ac'(t)"
2732     };
2733     static char  **legGem;
2734
2735     const char    *legLuzar[] = {
2736         "Ac\\sfin sys\\v{}\\z{}(t)",
2737         "Ac(t)",
2738         "Cc\\scontact,hb\\v{}\\z{}(t)",
2739         "-dAc\\sfs\\v{}\\z{}/dt"
2740     };
2741     gmx_bool       bNorm = FALSE, bOMP = FALSE;
2742     double         nhb   = 0;
2743     int            nhbi  = 0;
2744     real          *rhbex = NULL, *ht, *gt, *ght, *dght, *kt;
2745     real          *ct, *p_ct, tail, tail2, dtail, ct_fac, ght_fac, *cct;
2746     const real     tol     = 1e-3;
2747     int            nframes = hb->nframes, nf;
2748     unsigned int **h       = NULL, **g = NULL;
2749     int            nh, nhbonds, nhydro, ngh;
2750     t_hbond       *hbh;
2751     PSTYPE         p, *pfound = NULL, np;
2752     t_pShift      *pHist;
2753     int           *ptimes   = NULL, *poff = NULL, anhb, n0, mMax = INT_MIN;
2754     real         **rHbExGem = NULL;
2755     gmx_bool       c;
2756     int            acType;
2757     t_E           *E;
2758     double        *ctdouble, *timedouble, *fittedct;
2759     double         fittolerance = 0.1;
2760     int           *dondata      = NULL, thisThread;
2761
2762     enum {
2763         AC_NONE, AC_NN, AC_GEM, AC_LUZAR
2764     };
2765
2766 #ifdef GMX_OPENMP
2767     bOMP = TRUE;
2768 #else
2769     bOMP = FALSE;
2770 #endif
2771
2772     printf("Doing autocorrelation ");
2773
2774     /* Decide what kind of ACF calculations to do. */
2775     if (NN > NN_NONE && NN < NN_NR)
2776     {
2777 #ifdef HAVE_NN_LOOPS
2778         acType = AC_NN;
2779         printf("using the energy estimate.\n");
2780 #else
2781         acType = AC_NONE;
2782         printf("Can't do the NN-loop. Yet.\n");
2783 #endif
2784     }
2785     else if (hb->bGem)
2786     {
2787         acType = AC_GEM;
2788         printf("according to the reversible geminate recombination model by Omer Markowitch.\n");
2789
2790         nSets = 1 + (bBallistic ? 1 : 0) + (bGemFit ? 1 : 0);
2791         snew(legGem, nSets);
2792         for (i = 0; i < nSets; i++)
2793         {
2794             snew(legGem[i], 128);
2795         }
2796         sprintf(legGem[0], "Ac\\s%s\\v{}\\z{}(t)", gemType);
2797         if (bBallistic)
2798         {
2799             sprintf(legGem[1], "Ac'(t)");
2800         }
2801         if (bGemFit)
2802         {
2803             sprintf(legGem[(bBallistic ? 3 : 2)], "Ac\\s%s,fit\\v{}\\z{}(t)", gemType);
2804         }
2805
2806     }
2807     else
2808     {
2809         acType = AC_LUZAR;
2810         printf("according to the theory of Luzar and Chandler.\n");
2811     }
2812     fflush(stdout);
2813
2814     /* build hbexist matrix in reals for autocorr */
2815     /* Allocate memory for computing ACF (rhbex) and aggregating the ACF (ct) */
2816     n2 = 1;
2817     while (n2 < nframes)
2818     {
2819         n2 *= 2;
2820     }
2821
2822     nn = nframes/2;
2823
2824     if (acType != AC_NN || bOMP)
2825     {
2826         snew(h, hb->maxhydro);
2827         snew(g, hb->maxhydro);
2828     }
2829
2830     /* Dump hbonds for debugging */
2831     dump_ac(hb, bMerge || bContact, nDump);
2832
2833     /* Total number of hbonds analyzed here */
2834     nhbonds = 0;
2835     ngh     = 0;
2836     anhb    = 0;
2837
2838     if (acType != AC_LUZAR && bOMP)
2839     {
2840         nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
2841
2842         gmx_omp_set_num_threads(nThreads);
2843         snew(dondata, nThreads);
2844         for (i = 0; i < nThreads; i++)
2845         {
2846             dondata[i] = -1;
2847         }
2848         printf("ACF calculations parallelized with OpenMP using %i threads.\n"
2849                "Expect close to linear scaling over this donor-loop.\n", nThreads);
2850         fflush(stdout);
2851         fprintf(stderr, "Donors: [thread no]\n");
2852         {
2853             char tmpstr[7];
2854             for (i = 0; i < nThreads; i++)
2855             {
2856                 snprintf(tmpstr, 7, "[%i]", i);
2857                 fprintf(stderr, "%-7s", tmpstr);
2858             }
2859         }
2860         fprintf(stderr, "\n");
2861     }
2862
2863
2864     /* Build the ACF according to acType */
2865     switch (acType)
2866     {
2867
2868         case AC_NN:
2869 #ifdef HAVE_NN_LOOPS
2870             /* Here we're using the estimated energy for the hydrogen bonds. */
2871             snew(ct, nn);
2872
2873 #pragma omp parallel \
2874             private(i, j, k, nh, E, rhbex, thisThread) \
2875             default(shared)
2876             {
2877 #pragma omp barrier
2878                 thisThread = gmx_omp_get_thread_num();
2879                 rhbex      = NULL;
2880
2881                 snew(rhbex, n2);
2882                 memset(rhbex, 0, n2*sizeof(real)); /* Trust no-one, not even malloc()! */
2883
2884 #pragma omp barrier
2885 #pragma omp for schedule (dynamic)
2886                 for (i = 0; i < hb->d.nrd; i++) /* loop over donors */
2887                 {
2888                     if (bOMP)
2889                     {
2890 #pragma omp critical
2891                         {
2892                             dondata[thisThread] = i;
2893                             parallel_print(dondata, nThreads);
2894                         }
2895                     }
2896                     else
2897                     {
2898                         fprintf(stderr, "\r %i", i);
2899                     }
2900
2901                     for (j = 0; j < hb->a.nra; j++)              /* loop over acceptors */
2902                     {
2903                         for (nh = 0; nh < hb->d.nhydro[i]; nh++) /* loop over donors' hydrogens */
2904                         {
2905                             E = hb->hbE.E[i][j][nh];
2906                             if (E != NULL)
2907                             {
2908                                 for (k = 0; k < nframes; k++)
2909                                 {
2910                                     if (E[k] != NONSENSE_E)
2911                                     {
2912                                         rhbex[k] = (real)E[k];
2913                                     }
2914                                 }
2915
2916                                 low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &(rhbex), hb->time[1]-hb->time[0],
2917                                                 eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0, 1);
2918 #pragma omp critical
2919                                 {
2920                                     for (k = 0; (k < nn); k++)
2921                                     {
2922                                         ct[k] += rhbex[k];
2923                                     }
2924                                 }
2925                             }
2926                         } /* k loop */
2927                     }     /* j loop */
2928                 }         /* i loop */
2929                 sfree(rhbex);
2930 #pragma omp barrier
2931             }
2932
2933             if (bOMP)
2934             {
2935                 sfree(dondata);
2936             }
2937             normalizeACF(ct, NULL, 0, nn);
2938             snew(ctdouble, nn);
2939             snew(timedouble, nn);
2940             for (j = 0; j < nn; j++)
2941             {
2942                 timedouble[j] = (double)(hb->time[j]);
2943                 ctdouble[j]   = (double)(ct[j]);
2944             }
2945
2946             /* Remove ballistic term */
2947             /* Ballistic component removal and fitting to the reversible geminate recombination model
2948              * will be taken out for the time being. First of all, one can remove the ballistic
2949              * component with g_analyze afterwards. Secondly, and more importantly, there are still
2950              * problems with the robustness of the fitting to the model. More work is needed.
2951              * A third reason is that we're currently using gsl for this and wish to reduce dependence
2952              * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with
2953              * a BSD-licence that can do the job.
2954              *
2955              * - Erik Marklund, June 18 2010.
2956              */
2957 /*         if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */
2958 /*             takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */
2959 /*         else */
2960 /*             printf("\nNumber of data points is less than the number of parameters to fit\n." */
2961 /*                    "The system is underdetermined, hence no ballistic term can be found.\n\n"); */
2962
2963             fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)");
2964             xvgr_legend(fp, asize(legNN), legNN);
2965
2966             for (j = 0; (j < nn); j++)
2967             {
2968                 fprintf(fp, "%10g  %10g %10g\n",
2969                         hb->time[j]-hb->time[0],
2970                         ct[j],
2971                         ctdouble[j]);
2972             }
2973             xvgrclose(fp);
2974             sfree(ct);
2975             sfree(ctdouble);
2976             sfree(timedouble);
2977 #endif             /* HAVE_NN_LOOPS */
2978             break; /* case AC_NN */
2979
2980         case AC_GEM:
2981             snew(ct, 2*n2);
2982             memset(ct, 0, 2*n2*sizeof(real));
2983 #ifndef GMX_OPENMP
2984             fprintf(stderr, "Donor:\n");
2985 #define __ACDATA ct
2986 #else
2987 #define __ACDATA p_ct
2988 #endif
2989
2990 #pragma omp parallel \
2991             private(i, k, nh, hbh, pHist, h, g, n0, nf, np, j, m, \
2992             pfound, poff, rHbExGem, p, ihb, mMax, \
2993             thisThread, p_ct) \
2994             default(shared)
2995             { /* ##########  THE START OF THE ENORMOUS PARALLELIZED BLOCK!  ########## */
2996                 h          = NULL;
2997                 g          = NULL;
2998                 thisThread = gmx_omp_get_thread_num();
2999                 snew(h, hb->maxhydro);
3000                 snew(g, hb->maxhydro);
3001                 mMax     = INT_MIN;
3002                 rHbExGem = NULL;
3003                 poff     = NULL;
3004                 pfound   = NULL;
3005                 p_ct     = NULL;
3006                 snew(p_ct, 2*n2);
3007                 memset(p_ct, 0, 2*n2*sizeof(real));
3008
3009                 /* I'm using a chunk size of 1, since I expect      \
3010                  * the overhead to be really small compared         \
3011                  * to the actual calculations                       \ */
3012 #pragma omp for schedule(dynamic,1) nowait
3013                 for (i = 0; i < hb->d.nrd; i++)
3014                 {
3015
3016                     if (bOMP)
3017                     {
3018 #pragma omp critical
3019                         {
3020                             dondata[thisThread] = i;
3021                             parallel_print(dondata, nThreads);
3022                         }
3023                     }
3024                     else
3025                     {
3026                         fprintf(stderr, "\r %i", i);
3027                     }
3028                     for (k = 0; k < hb->a.nra; k++)
3029                     {
3030                         for (nh = 0; nh < ((bMerge || bContact) ? 1 : hb->d.nhydro[i]); nh++)
3031                         {
3032                             hbh = hb->hbmap[i][k];
3033                             if (hbh)
3034                             {
3035                                 /* Note that if hb->per->gemtype==gemDD, then distances will be stored in
3036                                  * hb->hbmap[d][a].h array anyway, because the contact flag will be set.
3037                                  * hence, it's only with the gemAD mode that hb->hbmap[d][a].g will be used. */
3038                                 pHist = &(hb->per->pHist[i][k]);
3039                                 if (ISHB(hbh->history[nh]) && pHist->len != 0)
3040                                 {
3041
3042                                     {
3043                                         h[nh] = hbh->h[nh];
3044                                         g[nh] = hb->per->gemtype == gemAD ? hbh->g[nh] : NULL;
3045                                     }
3046                                     n0 = hbh->n0;
3047                                     nf = hbh->nframes;
3048                                     /* count the number of periodic shifts encountered and store
3049                                      * them in separate arrays. */
3050                                     np = 0;
3051                                     for (j = 0; j < pHist->len; j++)
3052                                     {
3053                                         p = pHist->p[j];
3054                                         for (m = 0; m <= np; m++)
3055                                         {
3056                                             if (m == np) /* p not recognized in list. Add it and set up new array. */
3057                                             {
3058                                                 np++;
3059                                                 if (np > hb->per->nper)
3060                                                 {
3061                                                     gmx_fatal(FARGS, "Too many pshifts. Something's utterly wrong here.");
3062                                                 }
3063                                                 if (m >= mMax) /* Extend the arrays.
3064                                                                 * Doing it like this, using mMax to keep track of the sizes,
3065                                                                 * eleviates the need for freeing and re-allocating the arrays
3066                                                                 * when taking on the next donor-acceptor pair */
3067                                                 {
3068                                                     mMax = m;
3069                                                     srenew(pfound, np);   /* The list of found periodic shifts. */
3070                                                     srenew(rHbExGem, np); /* The hb existence functions (-aver_hb). */
3071                                                     snew(rHbExGem[m], 2*n2);
3072                                                     srenew(poff, np);
3073                                                 }
3074
3075                                                 {
3076                                                     if (rHbExGem != NULL && rHbExGem[m] != NULL)
3077                                                     {
3078                                                         /* This must be done, as this array was most likey
3079                                                          * used to store stuff in some previous iteration. */
3080                                                         memset(rHbExGem[m], 0, (sizeof(real)) * (2*n2));
3081                                                     }
3082                                                     else
3083                                                     {
3084                                                         fprintf(stderr, "rHbExGem not initialized! m = %i\n", m);
3085                                                     }
3086                                                 }
3087                                                 pfound[m] = p;
3088                                                 poff[m]   = -1;
3089
3090                                                 break;
3091                                             } /* m==np */
3092                                             if (p == pfound[m])
3093                                             {
3094                                                 break;
3095                                             }
3096                                         } /* m: Loop over found shifts */
3097                                     }     /* j: Loop over shifts */
3098
3099                                     /* Now unpack and disentangle the existence funtions. */
3100                                     for (j = 0; j < nf; j++)
3101                                     {
3102                                         /* i:       donor,
3103                                          * k:       acceptor
3104                                          * nh:      hydrogen
3105                                          * j:       time
3106                                          * p:       periodic shift
3107                                          * pfound:  list of periodic shifts found for this pair.
3108                                          * poff:    list of frame offsets; that is, the first
3109                                          *          frame a hbond has a particular periodic shift. */
3110                                         p = getPshift(*pHist, j+n0);
3111                                         if (p != -1)
3112                                         {
3113                                             for (m = 0; m < np; m++)
3114                                             {
3115                                                 if (pfound[m] == p)
3116                                                 {
3117                                                     break;
3118                                                 }
3119                                                 if (m == (np-1))
3120                                                 {
3121                                                     gmx_fatal(FARGS, "Shift not found, but must be there.");
3122                                                 }
3123                                             }
3124
3125                                             ihb = is_hb(h[nh], j) || ((hb->per->gemtype != gemAD || j == 0) ? FALSE : is_hb(g[nh], j));
3126                                             if (ihb)
3127                                             {
3128                                                 if (poff[m] == -1)
3129                                                 {
3130                                                     poff[m] = j; /* Here's where the first hbond with shift p is,
3131                                                                   * relative to the start of h[0].*/
3132                                                 }
3133                                                 if (j < poff[m])
3134                                                 {
3135                                                     gmx_fatal(FARGS, "j<poff[m]");
3136                                                 }
3137                                                 rHbExGem[m][j-poff[m]] += 1;
3138                                             }
3139                                         }
3140                                     }
3141
3142                                     /* Now, build ac. */
3143                                     for (m = 0; m < np; m++)
3144                                     {
3145                                         if (rHbExGem[m][0] > 0  && n0+poff[m] < nn /*  && m==0 */)
3146                                         {
3147                                             low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &(rHbExGem[m]), hb->time[1]-hb->time[0],
3148                                                             eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0);
3149                                             for (j = 0; (j < nn); j++)
3150                                             {
3151                                                 __ACDATA[j] += rHbExGem[m][j];
3152                                             }
3153                                         }
3154                                     } /* Building of ac. */
3155                                 }     /* if (ISHB(...*/
3156                             }         /* if (hbh) */
3157                         }             /* hydrogen loop */
3158                     }                 /* acceptor loop */
3159                 }                     /* donor loop */
3160
3161                 for (m = 0; m <= mMax; m++)
3162                 {
3163                     sfree(rHbExGem[m]);
3164                 }
3165                 sfree(pfound);
3166                 sfree(poff);
3167                 sfree(rHbExGem);
3168
3169                 sfree(h);
3170                 sfree(g);
3171
3172                 if (bOMP)
3173                 {
3174 #pragma omp critical
3175                     {
3176                         for (i = 0; i < nn; i++)
3177                         {
3178                             ct[i] += p_ct[i];
3179                         }
3180                     }
3181                     sfree(p_ct);
3182                 }
3183
3184             } /* ########## THE END OF THE ENORMOUS PARALLELIZED BLOCK ########## */
3185             if (bOMP)
3186             {
3187                 sfree(dondata);
3188             }
3189
3190             normalizeACF(ct, NULL, 0, nn);
3191
3192             fprintf(stderr, "\n\nACF successfully calculated.\n");
3193
3194             /* Use this part to fit to geminate recombination - JCP 129, 84505 (2008) */
3195
3196             snew(ctdouble, nn);
3197             snew(timedouble, nn);
3198             snew(fittedct, nn);
3199
3200             for (j = 0; j < nn; j++)
3201             {
3202                 timedouble[j] = (double)(hb->time[j]);
3203                 ctdouble[j]   = (double)(ct[j]);
3204             }
3205
3206             /* Remove ballistic term */
3207             /* Ballistic component removal and fitting to the reversible geminate recombination model
3208              * will be taken out for the time being. First of all, one can remove the ballistic
3209              * component with g_analyze afterwards. Secondly, and more importantly, there are still
3210              * problems with the robustness of the fitting to the model. More work is needed.
3211              * A third reason is that we're currently using gsl for this and wish to reduce dependence
3212              * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with
3213              * a BSD-licence that can do the job.
3214              *
3215              * - Erik Marklund, June 18 2010.
3216              */
3217 /*         if (bBallistic) { */
3218 /*             if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */
3219 /*                 takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */
3220 /*             else */
3221 /*                 printf("\nNumber of data points is less than the number of parameters to fit\n." */
3222 /*                        "The system is underdetermined, hence no ballistic term can be found.\n\n"); */
3223 /*         } */
3224 /*         if (bGemFit) */
3225 /*             fitGemRecomb(ctdouble, timedouble, &fittedct, nn, params); */
3226
3227
3228             if (bContact)
3229             {
3230                 fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3231             }
3232             else
3233             {
3234                 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3235             }
3236             xvgr_legend(fp, asize(legGem), (const char**)legGem, oenv);
3237
3238             for (j = 0; (j < nn); j++)
3239             {
3240                 fprintf(fp, "%10g  %10g", hb->time[j]-hb->time[0], ct[j]);
3241                 if (bBallistic)
3242                 {
3243                     fprintf(fp, "  %10g", ctdouble[j]);
3244                 }
3245                 if (bGemFit)
3246                 {
3247                     fprintf(fp, "  %10g", fittedct[j]);
3248                 }
3249                 fprintf(fp, "\n");
3250             }
3251             xvgrclose(fp);
3252
3253             sfree(ctdouble);
3254             sfree(timedouble);
3255             sfree(fittedct);
3256             sfree(ct);
3257
3258             break; /* case AC_GEM */
3259
3260         case AC_LUZAR:
3261             snew(rhbex, 2*n2);
3262             snew(ct, 2*n2);
3263             snew(gt, 2*n2);
3264             snew(ht, 2*n2);
3265             snew(ght, 2*n2);
3266             snew(dght, 2*n2);
3267
3268             snew(kt, nn);
3269             snew(cct, nn);
3270
3271             for (i = 0; (i < hb->d.nrd); i++)
3272             {
3273                 for (k = 0; (k < hb->a.nra); k++)
3274                 {
3275                     nhydro = 0;
3276                     hbh    = hb->hbmap[i][k];
3277
3278                     if (hbh)
3279                     {
3280                         if (bMerge || bContact)
3281                         {
3282                             if (ISHB(hbh->history[0]))
3283                             {
3284                                 h[0]   = hbh->h[0];
3285                                 g[0]   = hbh->g[0];
3286                                 nhydro = 1;
3287                             }
3288                         }
3289                         else
3290                         {
3291                             for (m = 0; (m < hb->maxhydro); m++)
3292                             {
3293                                 if (bContact ? ISDIST(hbh->history[m]) : ISHB(hbh->history[m]))
3294                                 {
3295                                     g[nhydro] = hbh->g[m];
3296                                     h[nhydro] = hbh->h[m];
3297                                     nhydro++;
3298                                 }
3299                             }
3300                         }
3301
3302                         nf = hbh->nframes;
3303                         for (nh = 0; (nh < nhydro); nh++)
3304                         {
3305                             int nrint = bContact ? hb->nrdist : hb->nrhb;
3306                             if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint))
3307                             {
3308                                 fprintf(stderr, "\rACF %d/%d", nhbonds+1, nrint);
3309                             }
3310                             nhbonds++;
3311                             for (j = 0; (j < nframes); j++)
3312                             {
3313                                 /* Changed '<' into '<=' below, just like I did in
3314                                    the hbm-output-loop in the gmx_hbond() block.
3315                                    - Erik Marklund, May 31, 2006 */
3316                                 if (j <= nf)
3317                                 {
3318                                     ihb   = is_hb(h[nh], j);
3319                                     idist = is_hb(g[nh], j);
3320                                 }
3321                                 else
3322                                 {
3323                                     ihb = idist = 0;
3324                                 }
3325                                 rhbex[j] = ihb;
3326                                 /* For contacts: if a second cut-off is provided, use it,
3327                                  * otherwise use g(t) = 1-h(t) */
3328                                 if (!R2 && bContact)
3329                                 {
3330                                     gt[j]  = 1-ihb;
3331                                 }
3332                                 else
3333                                 {
3334                                     gt[j]  = idist*(1-ihb);
3335                                 }
3336                                 ht[j]    = rhbex[j];
3337                                 nhb     += ihb;
3338                             }
3339
3340
3341                             /* The autocorrelation function is normalized after summation only */
3342                             low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &rhbex, hb->time[1]-hb->time[0],
3343                                             eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0);
3344
3345                             /* Cross correlation analysis for thermodynamics */
3346                             for (j = nframes; (j < n2); j++)
3347                             {
3348                                 ht[j] = 0;
3349                                 gt[j] = 0;
3350                             }
3351
3352                             cross_corr(n2, ht, gt, dght);
3353
3354                             for (j = 0; (j < nn); j++)
3355                             {
3356                                 ct[j]  += rhbex[j];
3357                                 ght[j] += dght[j];
3358                             }
3359                         }
3360                     }
3361                 }
3362             }
3363             fprintf(stderr, "\n");
3364             sfree(h);
3365             sfree(g);
3366             normalizeACF(ct, ght, nhb, nn);
3367
3368             /* Determine tail value for statistics */
3369             tail  = 0;
3370             tail2 = 0;
3371             for (j = nn/2; (j < nn); j++)
3372             {
3373                 tail  += ct[j];
3374                 tail2 += ct[j]*ct[j];
3375             }
3376             tail  /= (nn - nn/2);
3377             tail2 /= (nn - nn/2);
3378             dtail  = sqrt(tail2-tail*tail);
3379
3380             /* Check whether the ACF is long enough */
3381             if (dtail > tol)
3382             {
3383                 printf("\nWARNING: Correlation function is probably not long enough\n"
3384                        "because the standard deviation in the tail of C(t) > %g\n"
3385                        "Tail value (average C(t) over second half of acf): %g +/- %g\n",
3386                        tol, tail, dtail);
3387             }
3388             for (j = 0; (j < nn); j++)
3389             {
3390                 cct[j] = ct[j];
3391                 ct[j]  = (cct[j]-tail)/(1-tail);
3392             }
3393             /* Compute negative derivative k(t) = -dc(t)/dt */
3394             compute_derivative(nn, hb->time, ct, kt);
3395             for (j = 0; (j < nn); j++)
3396             {
3397                 kt[j] = -kt[j];
3398             }
3399
3400
3401             if (bContact)
3402             {
3403                 fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3404             }
3405             else
3406             {
3407                 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3408             }
3409             xvgr_legend(fp, asize(legLuzar), legLuzar, oenv);
3410
3411
3412             for (j = 0; (j < nn); j++)
3413             {
3414                 fprintf(fp, "%10g  %10g  %10g  %10g  %10g\n",
3415                         hb->time[j]-hb->time[0], ct[j], cct[j], ght[j], kt[j]);
3416             }
3417             ffclose(fp);
3418
3419             analyse_corr(nn, hb->time, ct, ght, kt, NULL, NULL, NULL,
3420                          fit_start, temp, smooth_tail_start, oenv);
3421
3422             do_view(oenv, fn, NULL);
3423             sfree(rhbex);
3424             sfree(ct);
3425             sfree(gt);
3426             sfree(ht);
3427             sfree(ght);
3428             sfree(dght);
3429             sfree(cct);
3430             sfree(kt);
3431             /* sfree(h); */
3432 /*         sfree(g); */
3433
3434             break; /* case AC_LUZAR */
3435
3436         default:
3437             gmx_fatal(FARGS, "Unrecognized type of ACF-calulation. acType = %i.", acType);
3438     } /* switch (acType) */
3439 }
3440
3441 static void init_hbframe(t_hbdata *hb, int nframes, real t)
3442 {
3443     int i, j, m;
3444
3445     hb->time[nframes]   = t;
3446     hb->nhb[nframes]    = 0;
3447     hb->ndist[nframes]  = 0;
3448     for (i = 0; (i < max_hx); i++)
3449     {
3450         hb->nhx[nframes][i] = 0;
3451     }
3452     /* Loop invalidated */
3453     if (hb->bHBmap && 0)
3454     {
3455         for (i = 0; (i < hb->d.nrd); i++)
3456         {
3457             for (j = 0; (j < hb->a.nra); j++)
3458             {
3459                 for (m = 0; (m < hb->maxhydro); m++)
3460                 {
3461                     if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[m])
3462                     {
3463                         set_hb(hb, i, m, j, nframes, HB_NO);
3464                     }
3465                 }
3466             }
3467         }
3468     }
3469     /*set_hb(hb->hbmap[i][j]->h[m],nframes-hb->hbmap[i][j]->n0,HB_NO);*/
3470 }
3471
3472 static void analyse_donor_props(const char *fn, t_hbdata *hb, int nframes, real t,
3473                                 const output_env_t oenv)
3474 {
3475     static FILE *fp    = NULL;
3476     const char  *leg[] = { "Nbound", "Nfree" };
3477     int          i, j, k, nbound, nb, nhtot;
3478
3479     if (!fn)
3480     {
3481         return;
3482     }
3483     if (!fp)
3484     {
3485         fp = xvgropen(fn, "Donor properties", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
3486         xvgr_legend(fp, asize(leg), leg, oenv);
3487     }
3488     nbound = 0;
3489     nhtot  = 0;
3490     for (i = 0; (i < hb->d.nrd); i++)
3491     {
3492         for (k = 0; (k < hb->d.nhydro[i]); k++)
3493         {
3494             nb = 0;
3495             nhtot++;
3496             for (j = 0; (j < hb->a.nra) && (nb == 0); j++)
3497             {
3498                 if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[k] &&
3499                     is_hb(hb->hbmap[i][j]->h[k], nframes))
3500                 {
3501                     nb = 1;
3502                 }
3503             }
3504             nbound += nb;
3505         }
3506     }
3507     fprintf(fp, "%10.3e  %6d  %6d\n", t, nbound, nhtot-nbound);
3508 }
3509
3510 static void dump_hbmap(t_hbdata *hb,
3511                        int nfile, t_filenm fnm[], gmx_bool bTwo,
3512                        gmx_bool bContact, int isize[], int *index[], char *grpnames[],
3513                        t_atoms *atoms)
3514 {
3515     FILE    *fp, *fplog;
3516     int      ddd, hhh, aaa, i, j, k, m, grp;
3517     char     ds[32], hs[32], as[32];
3518     gmx_bool first;
3519
3520     fp = opt2FILE("-hbn", nfile, fnm, "w");
3521     if (opt2bSet("-g", nfile, fnm))
3522     {
3523         fplog = ffopen(opt2fn("-g", nfile, fnm), "w");
3524         fprintf(fplog, "# %10s  %12s  %12s\n", "Donor", "Hydrogen", "Acceptor");
3525     }
3526     else
3527     {
3528         fplog = NULL;
3529     }
3530     for (grp = gr0; grp <= (bTwo ? gr1 : gr0); grp++)
3531     {
3532         fprintf(fp, "[ %s ]", grpnames[grp]);
3533         for (i = 0; i < isize[grp]; i++)
3534         {
3535             fprintf(fp, (i%15) ? " " : "\n");
3536             fprintf(fp, " %4u", index[grp][i]+1);
3537         }
3538         fprintf(fp, "\n");
3539         /*
3540            Added -contact support below.
3541            - Erik Marklund, May 29, 2006
3542          */
3543         if (!bContact)
3544         {
3545             fprintf(fp, "[ donors_hydrogens_%s ]\n", grpnames[grp]);
3546             for (i = 0; (i < hb->d.nrd); i++)
3547             {
3548                 if (hb->d.grp[i] == grp)
3549                 {
3550                     for (j = 0; (j < hb->d.nhydro[i]); j++)
3551                     {
3552                         fprintf(fp, " %4u %4u", hb->d.don[i]+1,
3553                                 hb->d.hydro[i][j]+1);
3554                     }
3555                     fprintf(fp, "\n");
3556                 }
3557             }
3558             first = TRUE;
3559             fprintf(fp, "[ acceptors_%s ]", grpnames[grp]);
3560             for (i = 0; (i < hb->a.nra); i++)
3561             {
3562                 if (hb->a.grp[i] == grp)
3563                 {
3564                     fprintf(fp, (i%15 && !first) ? " " : "\n");
3565                     fprintf(fp, " %4u", hb->a.acc[i]+1);
3566                     first = FALSE;
3567                 }
3568             }
3569             fprintf(fp, "\n");
3570         }
3571     }
3572     if (bTwo)
3573     {
3574         fprintf(fp, bContact ? "[ contacts_%s-%s ]\n" :
3575                 "[ hbonds_%s-%s ]\n", grpnames[0], grpnames[1]);
3576     }
3577     else
3578     {
3579         fprintf(fp, bContact ? "[ contacts_%s ]" : "[ hbonds_%s ]\n", grpnames[0]);
3580     }
3581
3582     for (i = 0; (i < hb->d.nrd); i++)
3583     {
3584         ddd = hb->d.don[i];
3585         for (k = 0; (k < hb->a.nra); k++)
3586         {
3587             aaa = hb->a.acc[k];
3588             for (m = 0; (m < hb->d.nhydro[i]); m++)
3589             {
3590                 if (hb->hbmap[i][k] && ISHB(hb->hbmap[i][k]->history[m]))
3591                 {
3592                     sprintf(ds, "%s", mkatomname(atoms, ddd));
3593                     sprintf(as, "%s", mkatomname(atoms, aaa));
3594                     if (bContact)
3595                     {
3596                         fprintf(fp, " %6u %6u\n", ddd+1, aaa+1);
3597                         if (fplog)
3598                         {
3599                             fprintf(fplog, "%12s  %12s\n", ds, as);
3600                         }
3601                     }
3602                     else
3603                     {
3604                         hhh = hb->d.hydro[i][m];
3605                         sprintf(hs, "%s", mkatomname(atoms, hhh));
3606                         fprintf(fp, " %6u %6u %6u\n", ddd+1, hhh+1, aaa+1);
3607                         if (fplog)
3608                         {
3609                             fprintf(fplog, "%12s  %12s  %12s\n", ds, hs, as);
3610                         }
3611                     }
3612                 }
3613             }
3614         }
3615     }
3616     ffclose(fp);
3617     if (fplog)
3618     {
3619         ffclose(fplog);
3620     }
3621 }
3622
3623 /* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template.
3624  * It mimics add_frames() and init_frame() to some extent. */
3625 static void sync_hbdata(t_hbdata *p_hb, int nframes)
3626 {
3627     int i;
3628     if (nframes >= p_hb->max_frames)
3629     {
3630         p_hb->max_frames += 4096;
3631         srenew(p_hb->nhb,   p_hb->max_frames);
3632         srenew(p_hb->ndist, p_hb->max_frames);
3633         srenew(p_hb->n_bound, p_hb->max_frames);
3634         srenew(p_hb->nhx, p_hb->max_frames);
3635         if (p_hb->bDAnr)
3636         {
3637             srenew(p_hb->danr, p_hb->max_frames);
3638         }
3639         memset(&(p_hb->nhb[nframes]),   0, sizeof(int) * (p_hb->max_frames-nframes));
3640         memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
3641         p_hb->nhb[nframes]   = 0;
3642         p_hb->ndist[nframes] = 0;
3643
3644     }
3645     p_hb->nframes = nframes;
3646 /*     for (i=0;) */
3647 /*     { */
3648 /*         p_hb->nhx[nframes][i] */
3649 /*     } */
3650     memset(&(p_hb->nhx[nframes]), 0, sizeof(int)*max_hx); /* zero the helix count for this frame */
3651
3652     /* hb->per will remain constant througout the frame loop,
3653      * even though the data its members point to will change,
3654      * hence no need for re-syncing. */
3655 }
3656
3657 int gmx_hbond(int argc, char *argv[])
3658 {
3659     const char        *desc[] = {
3660         "[THISMODULE] computes and analyzes hydrogen bonds. Hydrogen bonds are",
3661         "determined based on cutoffs for the angle Hydrogen - Donor - Acceptor",
3662         "(zero is extended) and the distance Donor - Acceptor",
3663         "(or Hydrogen - Acceptor using [TT]-noda[tt]).",
3664         "OH and NH groups are regarded as donors, O is an acceptor always,",
3665         "N is an acceptor by default, but this can be switched using",
3666         "[TT]-nitacc[tt]. Dummy hydrogen atoms are assumed to be connected",
3667         "to the first preceding non-hydrogen atom.[PAR]",
3668
3669         "You need to specify two groups for analysis, which must be either",
3670         "identical or non-overlapping. All hydrogen bonds between the two",
3671         "groups are analyzed.[PAR]",
3672
3673         "If you set [TT]-shell[tt], you will be asked for an additional index group",
3674         "which should contain exactly one atom. In this case, only hydrogen",
3675         "bonds between atoms within the shell distance from the one atom are",
3676         "considered.[PAR]",
3677
3678         "With option -ac, rate constants for hydrogen bonding can be derived with the model of Luzar and Chandler",
3679         "(Nature 394, 1996; J. Chem. Phys. 113:23, 2000) or that of Markovitz and Agmon (J. Chem. Phys 129, 2008).",
3680         "If contact kinetics are analyzed by using the -contact option, then",
3681         "n(t) can be defined as either all pairs that are not within contact distance r at time t",
3682         "(corresponding to leaving the -r2 option at the default value 0) or all pairs that",
3683         "are within distance r2 (corresponding to setting a second cut-off value with option -r2).",
3684         "See mentioned literature for more details and definitions."
3685         "[PAR]",
3686
3687         /*    "It is also possible to analyse specific hydrogen bonds with",
3688               "[TT]-sel[tt]. This index file must contain a group of atom triplets",
3689               "Donor Hydrogen Acceptor, in the following way:[PAR]",
3690          */
3691         "[TT]",
3692         "[ selected ][BR]",
3693         "     20    21    24[BR]",
3694         "     25    26    29[BR]",
3695         "      1     3     6[BR]",
3696         "[tt][BR]",
3697         "Note that the triplets need not be on separate lines.",
3698         "Each atom triplet specifies a hydrogen bond to be analyzed,",
3699         "note also that no check is made for the types of atoms.[PAR]",
3700
3701         "[BB]Output:[bb][BR]",
3702         "[TT]-num[tt]:  number of hydrogen bonds as a function of time.[BR]",
3703         "[TT]-ac[tt]:   average over all autocorrelations of the existence",
3704         "functions (either 0 or 1) of all hydrogen bonds.[BR]",
3705         "[TT]-dist[tt]: distance distribution of all hydrogen bonds.[BR]",
3706         "[TT]-ang[tt]:  angle distribution of all hydrogen bonds.[BR]",
3707         "[TT]-hx[tt]:   the number of n-n+i hydrogen bonds as a function of time",
3708         "where n and n+i stand for residue numbers and i ranges from 0 to 6.",
3709         "This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated",
3710         "with helices in proteins.[BR]",
3711         "[TT]-hbn[tt]:  all selected groups, donors, hydrogens and acceptors",
3712         "for selected groups, all hydrogen bonded atoms from all groups and",
3713         "all solvent atoms involved in insertion.[BR]",
3714         "[TT]-hbm[tt]:  existence matrix for all hydrogen bonds over all",
3715         "frames, this also contains information on solvent insertion",
3716         "into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]",
3717         "index file.[BR]",
3718         "[TT]-dan[tt]: write out the number of donors and acceptors analyzed for",
3719         "each timeframe. This is especially useful when using [TT]-shell[tt].[BR]",
3720         "[TT]-nhbdist[tt]: compute the number of HBonds per hydrogen in order to",
3721         "compare results to Raman Spectroscopy.",
3722         "[PAR]",
3723         "Note: options [TT]-ac[tt], [TT]-life[tt], [TT]-hbn[tt] and [TT]-hbm[tt]",
3724         "require an amount of memory proportional to the total numbers of donors",
3725         "times the total number of acceptors in the selected group(s)."
3726     };
3727
3728     static real        acut     = 30, abin = 1, rcut = 0.35, r2cut = 0, rbin = 0.005, rshell = -1;
3729     static real        maxnhb   = 0, fit_start = 1, fit_end = 60, temp = 298.15, smooth_tail_start = -1, D = -1;
3730     static gmx_bool    bNitAcc  = TRUE, bDA = TRUE, bMerge = TRUE;
3731     static int         nDump    = 0, nFitPoints = 100;
3732     static int         nThreads = 0, nBalExp = 4;
3733
3734     static gmx_bool    bContact     = FALSE, bBallistic = FALSE, bGemFit = FALSE;
3735     static real        logAfterTime = 10, gemBallistic = 0.2; /* ps */
3736     static const char *NNtype[]     = {NULL, "none", "binary", "oneOverR3", "dipole", NULL};
3737
3738     /* options */
3739     t_pargs     pa [] = {
3740         { "-a",    FALSE,  etREAL, {&acut},
3741           "Cutoff angle (degrees, Hydrogen - Donor - Acceptor)" },
3742         { "-r",    FALSE,  etREAL, {&rcut},
3743           "Cutoff radius (nm, X - Acceptor, see next option)" },
3744         { "-da",   FALSE,  etBOOL, {&bDA},
3745           "Use distance Donor-Acceptor (if TRUE) or Hydrogen-Acceptor (FALSE)" },
3746         { "-r2",   FALSE,  etREAL, {&r2cut},
3747           "Second cutoff radius. Mainly useful with [TT]-contact[tt] and [TT]-ac[tt]"},
3748         { "-abin", FALSE,  etREAL, {&abin},
3749           "Binwidth angle distribution (degrees)" },
3750         { "-rbin", FALSE,  etREAL, {&rbin},
3751           "Binwidth distance distribution (nm)" },
3752         { "-nitacc", FALSE, etBOOL, {&bNitAcc},
3753           "Regard nitrogen atoms as acceptors" },
3754         { "-contact", FALSE, etBOOL, {&bContact},
3755           "Do not look for hydrogen bonds, but merely for contacts within the cut-off distance" },
3756         { "-shell", FALSE, etREAL, {&rshell},
3757           "when > 0, only calculate hydrogen bonds within # nm shell around "
3758           "one particle" },
3759         { "-fitstart", FALSE, etREAL, {&fit_start},
3760           "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]" },
3761         { "-fitstart", FALSE, etREAL, {&fit_start},
3762           "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])" },
3763         { "-temp",  FALSE, etREAL, {&temp},
3764           "Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" },
3765         { "-smooth", FALSE, etREAL, {&smooth_tail_start},
3766           "If >= 0, the tail of the ACF will be smoothed by fitting it to an exponential function: y = A exp(-x/[GRK]tau[grk])" },
3767         { "-dump",  FALSE, etINT, {&nDump},
3768           "Dump the first N hydrogen bond ACFs in a single [TT].xvg[tt] file for debugging" },
3769         { "-max_hb", FALSE, etREAL, {&maxnhb},
3770           "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" },
3771         { "-merge", FALSE, etBOOL, {&bMerge},
3772           "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." },
3773         { "-geminate", FALSE, etENUM, {gemType},
3774           "Use reversible geminate recombination for the kinetics/thermodynamics calclations. See Markovitch et al., J. Chem. Phys 129, 084505 (2008) for details."},
3775         { "-diff", FALSE, etREAL, {&D},
3776           "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."},
3777 #ifdef GMX_OPENMP
3778         { "-nthreads", FALSE, etINT, {&nThreads},
3779           "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)"},
3780 #endif
3781     };
3782     const char *bugs[] = {
3783         "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and therefore not available for the time being."
3784     };
3785     t_filenm    fnm[] = {
3786         { efTRX, "-f",   NULL,     ffREAD  },
3787         { efTPX, NULL,   NULL,     ffREAD  },
3788         { efNDX, NULL,   NULL,     ffOPTRD },
3789         /*    { efNDX, "-sel", "select", ffOPTRD },*/
3790         { efXVG, "-num", "hbnum",  ffWRITE },
3791         { efLOG, "-g",   "hbond",  ffOPTWR },
3792         { efXVG, "-ac",  "hbac",   ffOPTWR },
3793         { efXVG, "-dist", "hbdist", ffOPTWR },
3794         { efXVG, "-ang", "hbang",  ffOPTWR },
3795         { efXVG, "-hx",  "hbhelix", ffOPTWR },
3796         { efNDX, "-hbn", "hbond",  ffOPTWR },
3797         { efXPM, "-hbm", "hbmap",  ffOPTWR },
3798         { efXVG, "-don", "donor",  ffOPTWR },
3799         { efXVG, "-dan", "danum",  ffOPTWR },
3800         { efXVG, "-life", "hblife", ffOPTWR },
3801         { efXVG, "-nhbdist", "nhbdist", ffOPTWR }
3802
3803     };
3804 #define NFILE asize(fnm)
3805
3806     char                  hbmap [HB_NR] = { ' ',    'o',      '-',       '*' };
3807     const char           *hbdesc[HB_NR] = { "None", "Present", "Inserted", "Present & Inserted" };
3808     t_rgb                 hbrgb [HB_NR] = { {1, 1, 1}, {1, 0, 0},   {0, 0, 1},    {1, 0, 1} };
3809
3810     t_trxstatus          *status;
3811     int                   trrStatus = 1;
3812     t_topology            top;
3813     t_inputrec            ir;
3814     t_pargs              *ppa;
3815     int                   npargs, natoms, nframes = 0, shatom;
3816     int                  *isize;
3817     char                **grpnames;
3818     atom_id             **index;
3819     rvec                 *x, hbox;
3820     matrix                box;
3821     real                  t, ccut, dist = 0.0, ang = 0.0;
3822     double                max_nhb, aver_nhb, aver_dist;
3823     int                   h = 0, i = 0, j, k = 0, l, start, end, id, ja, ogrp, nsel;
3824     int                   xi, yi, zi, ai;
3825     int                   xj, yj, zj, aj, xjj, yjj, zjj;
3826     int                   xk, yk, zk, ak, xkk, ykk, zkk;
3827     gmx_bool              bSelected, bHBmap, bStop, bTwo, was, bBox, bTric;
3828     int                  *adist, *rdist, *aptr, *rprt;
3829     int                   grp, nabin, nrbin, bin, resdist, ihb;
3830     char                **leg;
3831     t_hbdata             *hb, *hbptr;
3832     FILE                 *fp, *fpins = NULL, *fpnhb = NULL;
3833     t_gridcell         ***grid;
3834     t_ncell              *icell, *jcell, *kcell;
3835     ivec                  ngrid;
3836     unsigned char        *datable;
3837     output_env_t          oenv;
3838     int                   gemmode, NN;
3839     PSTYPE                peri = 0;
3840     t_E                   E;
3841     int                   ii, jj, hh, actual_nThreads;
3842     int                   threadNr = 0;
3843     gmx_bool              bGem, bNN, bParallel;
3844     t_gemParams          *params = NULL;
3845     gmx_bool              bEdge_yjj, bEdge_xjj, bOMP;
3846
3847     t_hbdata            **p_hb    = NULL;                   /* one per thread, then merge after the frame loop */
3848     int                 **p_adist = NULL, **p_rdist = NULL; /* a histogram for each thread. */
3849
3850 #ifdef GMX_OPENMP
3851     bOMP = TRUE;
3852 #else
3853     bOMP = FALSE;
3854 #endif
3855
3856     npargs = asize(pa);
3857     ppa    = add_acf_pargs(&npargs, pa);
3858
3859     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE, NFILE, fnm, npargs,
3860                            ppa, asize(desc), desc, asize(bugs), bugs, &oenv))
3861     {
3862         return 0;
3863     }
3864
3865     /* NN-loop? If so, what estimator to use ?*/
3866     NN = 1;
3867     /* Outcommented for now DvdS 2010-07-13
3868        while (NN < NN_NR && gmx_strcasecmp(NNtype[0], NNtype[NN])!=0)
3869         NN++;
3870        if (NN == NN_NR)
3871         gmx_fatal(FARGS, "Invalid NN-loop type.");
3872      */
3873     bNN = FALSE;
3874     for (i = 2; bNN == FALSE && i < NN_NR; i++)
3875     {
3876         bNN = bNN || NN == i;
3877     }
3878
3879     if (NN > NN_NONE && bMerge)
3880     {
3881         bMerge = FALSE;
3882     }
3883
3884     /* geminate recombination? If so, which flavor? */
3885     gemmode = 1;
3886     while (gemmode < gemNR && gmx_strcasecmp(gemType[0], gemType[gemmode]) != 0)
3887     {
3888         gemmode++;
3889     }
3890     if (gemmode == gemNR)
3891     {
3892         gmx_fatal(FARGS, "Invalid recombination type.");
3893     }
3894
3895     bGem = FALSE;
3896     for (i = 2; bGem == FALSE && i < gemNR; i++)
3897     {
3898         bGem = bGem || gemmode == i;
3899     }
3900
3901     if (bGem)
3902     {
3903         printf("Geminate recombination: %s\n", gemType[gemmode]);
3904 #ifndef HAVE_LIBGSL
3905         printf("Note that some aspects of reversible geminate recombination won't work without gsl.\n");
3906 #endif
3907         if (bContact)
3908         {
3909             if (gemmode != gemDD)
3910             {
3911                 printf("Turning off -contact option...\n");
3912                 bContact = FALSE;
3913             }
3914         }
3915         else
3916         {
3917             if (gemmode == gemDD)
3918             {
3919                 printf("Turning on -contact option...\n");
3920                 bContact = TRUE;
3921             }
3922         }
3923         if (bMerge)
3924         {
3925             if (gemmode == gemAA)
3926             {
3927                 printf("Turning off -merge option...\n");
3928                 bMerge = FALSE;
3929             }
3930         }
3931         else
3932         {
3933             if (gemmode != gemAA)
3934             {
3935                 printf("Turning on -merge option...\n");
3936                 bMerge = TRUE;
3937             }
3938         }
3939     }
3940
3941     /* process input */
3942     bSelected = FALSE;
3943     ccut      = cos(acut*DEG2RAD);
3944
3945     if (bContact)
3946     {
3947         if (bSelected)
3948         {
3949             gmx_fatal(FARGS, "Can not analyze selected contacts.");
3950         }
3951         if (!bDA)
3952         {
3953             gmx_fatal(FARGS, "Can not analyze contact between H and A: turn off -noda");
3954         }
3955     }
3956
3957     /* Initiate main data structure! */
3958     bHBmap = (opt2bSet("-ac", NFILE, fnm) ||
3959               opt2bSet("-life", NFILE, fnm) ||
3960               opt2bSet("-hbn", NFILE, fnm) ||
3961               opt2bSet("-hbm", NFILE, fnm) ||
3962               bGem);
3963
3964     if (opt2bSet("-nhbdist", NFILE, fnm))
3965     {
3966         const char *leg[MAXHH+1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" };
3967         fpnhb = xvgropen(opt2fn("-nhbdist", NFILE, fnm),
3968                          "Number of donor-H with N HBs", output_env_get_xvgr_tlabel(oenv), "N", oenv);
3969         xvgr_legend(fpnhb, asize(leg), leg, oenv);
3970     }
3971
3972     hb = mk_hbdata(bHBmap, opt2bSet("-dan", NFILE, fnm), bMerge || bContact, bGem, gemmode);
3973
3974     /* get topology */
3975     read_tpx_top(ftp2fn(efTPX, NFILE, fnm), &ir, box, &natoms, NULL, NULL, NULL, &top);
3976
3977     snew(grpnames, grNR);
3978     snew(index, grNR);
3979     snew(isize, grNR);
3980     /* Make Donor-Acceptor table */
3981     snew(datable, top.atoms.nr);
3982     gen_datable(index[0], isize[0], datable, top.atoms.nr);
3983
3984     if (bSelected)
3985     {
3986         /* analyze selected hydrogen bonds */
3987         printf("Select group with selected atoms:\n");
3988         get_index(&(top.atoms), opt2fn("-sel", NFILE, fnm),
3989                   1, &nsel, index, grpnames);
3990         if (nsel % 3)
3991         {
3992             gmx_fatal(FARGS, "Number of atoms in group '%s' not a multiple of 3\n"
3993                       "and therefore cannot contain triplets of "
3994                       "Donor-Hydrogen-Acceptor", grpnames[0]);
3995         }
3996         bTwo = FALSE;
3997
3998         for (i = 0; (i < nsel); i += 3)
3999         {
4000             int dd = index[0][i];
4001             int aa = index[0][i+2];
4002             /* int */ hh = index[0][i+1];
4003             add_dh (&hb->d, dd, hh, i, datable);
4004             add_acc(&hb->a, aa, i);
4005             /* Should this be here ? */
4006             snew(hb->d.dptr, top.atoms.nr);
4007             snew(hb->a.aptr, top.atoms.nr);
4008             add_hbond(hb, dd, aa, hh, gr0, gr0, 0, bMerge, 0, bContact, peri);
4009         }
4010         printf("Analyzing %d selected hydrogen bonds from '%s'\n",
4011                isize[0], grpnames[0]);
4012     }
4013     else
4014     {
4015         /* analyze all hydrogen bonds: get group(s) */
4016         printf("Specify 2 groups to analyze:\n");
4017         get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
4018                   2, isize, index, grpnames);
4019
4020         /* check if we have two identical or two non-overlapping groups */
4021         bTwo = isize[0] != isize[1];
4022         for (i = 0; (i < isize[0]) && !bTwo; i++)
4023         {
4024             bTwo = index[0][i] != index[1][i];
4025         }
4026         if (bTwo)
4027         {
4028             printf("Checking for overlap in atoms between %s and %s\n",
4029                    grpnames[0], grpnames[1]);
4030             for (i = 0; i < isize[1]; i++)
4031             {
4032                 if (ISINGRP(datable[index[1][i]]))
4033                 {
4034                     gmx_fatal(FARGS, "Partial overlap between groups '%s' and '%s'",
4035                               grpnames[0], grpnames[1]);
4036                 }
4037             }
4038             /*
4039                printf("Checking for overlap in atoms between %s and %s\n",
4040                grpnames[0],grpnames[1]);
4041                for (i=0; i<isize[0]; i++)
4042                for (j=0; j<isize[1]; j++)
4043                if (index[0][i] == index[1][j])
4044                gmx_fatal(FARGS,"Partial overlap between groups '%s' and '%s'",
4045                grpnames[0],grpnames[1]);
4046              */
4047         }
4048         if (bTwo)
4049         {
4050             printf("Calculating %s "
4051                    "between %s (%d atoms) and %s (%d atoms)\n",
4052                    bContact ? "contacts" : "hydrogen bonds",
4053                    grpnames[0], isize[0], grpnames[1], isize[1]);
4054         }
4055         else
4056         {
4057             fprintf(stderr, "Calculating %s in %s (%d atoms)\n",
4058                     bContact ? "contacts" : "hydrogen bonds", grpnames[0], isize[0]);
4059         }
4060     }
4061     sfree(datable);
4062
4063     /* search donors and acceptors in groups */
4064     snew(datable, top.atoms.nr);
4065     for (i = 0; (i < grNR); i++)
4066     {
4067         if ( ((i == gr0) && !bSelected ) ||
4068              ((i == gr1) && bTwo ))
4069         {
4070             gen_datable(index[i], isize[i], datable, top.atoms.nr);
4071             if (bContact)
4072             {
4073                 search_acceptors(&top, isize[i], index[i], &hb->a, i,
4074                                  bNitAcc, TRUE, (bTwo && (i == gr0)) || !bTwo, datable);
4075                 search_donors   (&top, isize[i], index[i], &hb->d, i,
4076                                  TRUE, (bTwo && (i == gr1)) || !bTwo, datable);
4077             }
4078             else
4079             {
4080                 search_acceptors(&top, isize[i], index[i], &hb->a, i, bNitAcc, FALSE, TRUE, datable);
4081                 search_donors   (&top, isize[i], index[i], &hb->d, i, FALSE, TRUE, datable);
4082             }
4083             if (bTwo)
4084             {
4085                 clear_datable_grp(datable, top.atoms.nr);
4086             }
4087         }
4088     }
4089     sfree(datable);
4090     printf("Found %d donors and %d acceptors\n", hb->d.nrd, hb->a.nra);
4091     /*if (bSelected)
4092        snew(donors[gr0D], dons[gr0D].nrd);*/
4093
4094     if (bHBmap)
4095     {
4096         printf("Making hbmap structure...");
4097         /* Generate hbond data structure */
4098         mk_hbmap(hb);
4099         printf("done.\n");
4100     }
4101
4102 #ifdef HAVE_NN_LOOPS
4103     if (bNN)
4104     {
4105         mk_hbEmap(hb, 0);
4106     }
4107 #endif
4108
4109     if (bGem)
4110     {
4111         printf("Making per structure...");
4112         /* Generate hbond data structure */
4113         mk_per(hb);
4114         printf("done.\n");
4115     }
4116
4117     /* check input */
4118     bStop = FALSE;
4119     if (hb->d.nrd + hb->a.nra == 0)
4120     {
4121         printf("No Donors or Acceptors found\n");
4122         bStop = TRUE;
4123     }
4124     if (!bStop)
4125     {
4126         if (hb->d.nrd == 0)
4127         {
4128             printf("No Donors found\n");
4129             bStop = TRUE;
4130         }
4131         if (hb->a.nra == 0)
4132         {
4133             printf("No Acceptors found\n");
4134             bStop = TRUE;
4135         }
4136     }
4137     if (bStop)
4138     {
4139         gmx_fatal(FARGS, "Nothing to be done");
4140     }
4141
4142     shatom = 0;
4143     if (rshell > 0)
4144     {
4145         int      shisz;
4146         atom_id *shidx;
4147         char    *shgrpnm;
4148         /* get index group with atom for shell */
4149         do
4150         {
4151             printf("Select atom for shell (1 atom):\n");
4152             get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
4153                       1, &shisz, &shidx, &shgrpnm);
4154             if (shisz != 1)
4155             {
4156                 printf("group contains %d atoms, should be 1 (one)\n", shisz);
4157             }
4158         }
4159         while (shisz != 1);
4160         shatom = shidx[0];
4161         printf("Will calculate hydrogen bonds within a shell "
4162                "of %g nm around atom %i\n", rshell, shatom+1);
4163     }
4164
4165     /* Analyze trajectory */
4166     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
4167     if (natoms > top.atoms.nr)
4168     {
4169         gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)",
4170                   top.atoms.nr, natoms);
4171     }
4172
4173     bBox  = ir.ePBC != epbcNONE;
4174     grid  = init_grid(bBox, box, (rcut > r2cut) ? rcut : r2cut, ngrid);
4175     nabin = acut/abin;
4176     nrbin = rcut/rbin;
4177     snew(adist, nabin+1);
4178     snew(rdist, nrbin+1);
4179
4180     if (bGem && !bBox)
4181     {
4182         gmx_fatal(FARGS, "Can't do geminate recombination without periodic box.");
4183     }
4184
4185     bParallel = FALSE;
4186
4187 #ifndef GMX_OPENMP
4188 #define __ADIST adist
4189 #define __RDIST rdist
4190 #define __HBDATA hb
4191 #else /* GMX_OPENMP ================================================== \
4192        * Set up the OpenMP stuff,                                       |
4193        * like the number of threads and such                            |
4194        * Also start the parallel loop.                                  |
4195        */
4196 #define __ADIST p_adist[threadNr]
4197 #define __RDIST p_rdist[threadNr]
4198 #define __HBDATA p_hb[threadNr]
4199 #endif
4200     if (bOMP)
4201     {
4202         bParallel = !bSelected;
4203
4204         if (bParallel)
4205         {
4206             actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
4207
4208             gmx_omp_set_num_threads(actual_nThreads);
4209             printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads);
4210             fflush(stdout);
4211         }
4212         else
4213         {
4214             actual_nThreads = 1;
4215         }
4216
4217         snew(p_hb,    actual_nThreads);
4218         snew(p_adist, actual_nThreads);
4219         snew(p_rdist, actual_nThreads);
4220         for (i = 0; i < actual_nThreads; i++)
4221         {
4222             snew(p_hb[i], 1);
4223             snew(p_adist[i], nabin+1);
4224             snew(p_rdist[i], nrbin+1);
4225
4226             p_hb[i]->max_frames = 0;
4227             p_hb[i]->nhb        = NULL;
4228             p_hb[i]->ndist      = NULL;
4229             p_hb[i]->n_bound    = NULL;
4230             p_hb[i]->time       = NULL;
4231             p_hb[i]->nhx        = NULL;
4232
4233             p_hb[i]->bHBmap     = hb->bHBmap;
4234             p_hb[i]->bDAnr      = hb->bDAnr;
4235             p_hb[i]->bGem       = hb->bGem;
4236             p_hb[i]->wordlen    = hb->wordlen;
4237             p_hb[i]->nframes    = hb->nframes;
4238             p_hb[i]->maxhydro   = hb->maxhydro;
4239             p_hb[i]->danr       = hb->danr;
4240             p_hb[i]->d          = hb->d;
4241             p_hb[i]->a          = hb->a;
4242             p_hb[i]->hbmap      = hb->hbmap;
4243             p_hb[i]->time       = hb->time; /* This may need re-syncing at every frame. */
4244             p_hb[i]->per        = hb->per;
4245
4246 #ifdef HAVE_NN_LOOPS
4247             p_hb[i]->hbE = hb->hbE;
4248 #endif
4249
4250             p_hb[i]->nrhb   = 0;
4251             p_hb[i]->nrdist = 0;
4252         }
4253     }
4254
4255     /* Make a thread pool here,
4256      * instead of forking anew at every frame. */
4257
4258 #pragma omp parallel \
4259     firstprivate(i) \
4260     private(j, h, ii, jj, hh, E, \
4261     xi, yi, zi, xj, yj, zj, threadNr, \
4262     dist, ang, peri, icell, jcell, \
4263     grp, ogrp, ai, aj, xjj, yjj, zjj, \
4264     xk, yk, zk, ihb, id,  resdist, \
4265     xkk, ykk, zkk, kcell, ak, k, bTric, \
4266     bEdge_xjj, bEdge_yjj) \
4267     default(shared)
4268     {    /* Start of parallel region */
4269         threadNr = gmx_omp_get_thread_num();
4270
4271         do
4272         {
4273
4274             bTric = bBox && TRICLINIC(box);
4275
4276             if (bOMP)
4277             {
4278                 sync_hbdata(p_hb[threadNr], nframes);
4279             }
4280 #pragma omp single
4281             {
4282                 build_grid(hb, x, x[shatom], bBox, box, hbox, (rcut > r2cut) ? rcut : r2cut,
4283                            rshell, ngrid, grid);
4284                 reset_nhbonds(&(hb->d));
4285
4286                 if (debug && bDebug)
4287                 {
4288                     dump_grid(debug, ngrid, grid);
4289                 }
4290
4291                 add_frames(hb, nframes);
4292                 init_hbframe(hb, nframes, output_env_conv_time(oenv, t));
4293
4294                 if (hb->bDAnr)
4295                 {
4296                     count_da_grid(ngrid, grid, hb->danr[nframes]);
4297                 }
4298             } /* omp single */
4299
4300             if (bOMP)
4301             {
4302                 p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */
4303             }
4304
4305             if (bNN)
4306             {
4307 #ifdef HAVE_NN_LOOPS /* Unlock this feature when testing */
4308                 /* Loop over all atom pairs and estimate interaction energy */
4309
4310 #pragma omp single
4311                 {
4312                     addFramesNN(hb, nframes);
4313                 }
4314
4315 #pragma omp barrier
4316 #pragma omp for schedule(dynamic)
4317                 for (i = 0; i < hb->d.nrd; i++)
4318                 {
4319                     for (j = 0; j < hb->a.nra; j++)
4320                     {
4321                         for (h = 0;
4322                              h < (bContact ? 1 : hb->d.nhydro[i]);
4323                              h++)
4324                         {
4325                             if (i == hb->d.nrd || j == hb->a.nra)
4326                             {
4327                                 gmx_fatal(FARGS, "out of bounds");
4328                             }
4329
4330                             /* Get the real atom ids */
4331                             ii = hb->d.don[i];
4332                             jj = hb->a.acc[j];
4333                             hh = hb->d.hydro[i][h];
4334
4335                             /* Estimate the energy from the geometry */
4336                             E = calcHbEnergy(ii, jj, hh, x, NN, box, hbox, &(hb->d));
4337                             /* Store the energy */
4338                             storeHbEnergy(hb, i, j, h, E, nframes);
4339                         }
4340                     }
4341                 }
4342 #endif        /* HAVE_NN_LOOPS */
4343             } /* if (bNN)*/
4344             else
4345             {
4346                 if (bSelected)
4347                 {
4348
4349 #pragma omp single
4350                     {
4351                         /* Do not parallelize this just yet. */
4352                         /* int ii; */
4353                         for (ii = 0; (ii < nsel); ii++)
4354                         {
4355                             int dd = index[0][i];
4356                             int aa = index[0][i+2];
4357                             /* int */ hh = index[0][i+1];
4358                             ihb          = is_hbond(hb, ii, ii, dd, aa, rcut, r2cut, ccut, x, bBox, box,
4359                                                     hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri);
4360
4361                             if (ihb)
4362                             {
4363                                 /* add to index if not already there */
4364                                 /* Add a hbond */
4365                                 add_hbond(hb, dd, aa, hh, ii, ii, nframes, bMerge, ihb, bContact, peri);
4366                             }
4367                         }
4368                     } /* omp single */
4369                 }     /* if (bSelected) */
4370                 else
4371                 {
4372
4373 #pragma omp single
4374                     {
4375                         if (bGem)
4376                         {
4377                             calcBoxProjection(box, hb->per->P);
4378                         }
4379
4380                         /* loop over all gridcells (xi,yi,zi)      */
4381                         /* Removed confusing macro, DvdS 27/12/98  */
4382
4383                     }
4384                     /* The outer grid loop will have to do for now. */
4385 #pragma omp for schedule(dynamic)
4386                     for (xi = 0; xi < ngrid[XX]; xi++)
4387                     {
4388                         for (yi = 0; (yi < ngrid[YY]); yi++)
4389                         {
4390                             for (zi = 0; (zi < ngrid[ZZ]); zi++)
4391                             {
4392
4393                                 /* loop over donor groups gr0 (always) and gr1 (if necessary) */
4394                                 for (grp = gr0; (grp <= (bTwo ? gr1 : gr0)); grp++)
4395                                 {
4396                                     icell = &(grid[zi][yi][xi].d[grp]);
4397
4398                                     if (bTwo)
4399                                     {
4400                                         ogrp = 1-grp;
4401                                     }
4402                                     else
4403                                     {
4404                                         ogrp = grp;
4405                                     }
4406
4407                                     /* loop over all hydrogen atoms from group (grp)
4408                                      * in this gridcell (icell)
4409                                      */
4410                                     for (ai = 0; (ai < icell->nr); ai++)
4411                                     {
4412                                         i  = icell->atoms[ai];
4413
4414                                         /* loop over all adjacent gridcells (xj,yj,zj) */
4415                                         for (zjj = grid_loop_begin(ngrid[ZZ], zi, bTric, FALSE);
4416                                              zjj <= grid_loop_end(ngrid[ZZ], zi, bTric, FALSE);
4417                                              zjj++)
4418                                         {
4419                                             zj        = grid_mod(zjj, ngrid[ZZ]);
4420                                             bEdge_yjj = (zj == 0) || (zj == ngrid[ZZ] - 1);
4421                                             for (yjj = grid_loop_begin(ngrid[YY], yi, bTric, bEdge_yjj);
4422                                                  yjj <= grid_loop_end(ngrid[YY], yi, bTric, bEdge_yjj);
4423                                                  yjj++)
4424                                             {
4425                                                 yj        = grid_mod(yjj, ngrid[YY]);
4426                                                 bEdge_xjj =
4427                                                     (yj == 0) || (yj == ngrid[YY] - 1) ||
4428                                                     (zj == 0) || (zj == ngrid[ZZ] - 1);
4429                                                 for (xjj = grid_loop_begin(ngrid[XX], xi, bTric, bEdge_xjj);
4430                                                      xjj <= grid_loop_end(ngrid[XX], xi, bTric, bEdge_xjj);
4431                                                      xjj++)
4432                                                 {
4433                                                     xj    = grid_mod(xjj, ngrid[XX]);
4434                                                     jcell = &(grid[zj][yj][xj].a[ogrp]);
4435                                                     /* loop over acceptor atoms from other group (ogrp)
4436                                                      * in this adjacent gridcell (jcell)
4437                                                      */
4438                                                     for (aj = 0; (aj < jcell->nr); aj++)
4439                                                     {
4440                                                         j = jcell->atoms[aj];
4441
4442                                                         /* check if this once was a h-bond */
4443                                                         peri = -1;
4444                                                         ihb  = is_hbond(__HBDATA, grp, ogrp, i, j, rcut, r2cut, ccut, x, bBox, box,
4445                                                                         hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri);
4446
4447                                                         if (ihb)
4448                                                         {
4449                                                             /* add to index if not already there */
4450                                                             /* Add a hbond */
4451                                                             add_hbond(__HBDATA, i, j, h, grp, ogrp, nframes, bMerge, ihb, bContact, peri);
4452
4453                                                             /* make angle and distance distributions */
4454                                                             if (ihb == hbHB && !bContact)
4455                                                             {
4456                                                                 if (dist > rcut)
4457                                                                 {
4458                                                                     gmx_fatal(FARGS, "distance is higher than what is allowed for an hbond: %f", dist);
4459                                                                 }
4460                                                                 ang *= RAD2DEG;
4461                                                                 __ADIST[(int)( ang/abin)]++;
4462                                                                 __RDIST[(int)(dist/rbin)]++;
4463                                                                 if (!bTwo)
4464                                                                 {
4465                                                                     int id, ia;
4466                                                                     if ((id = donor_index(&hb->d, grp, i)) == NOTSET)
4467                                                                     {
4468                                                                         gmx_fatal(FARGS, "Invalid donor %d", i);
4469                                                                     }
4470                                                                     if ((ia = acceptor_index(&hb->a, ogrp, j)) == NOTSET)
4471                                                                     {
4472                                                                         gmx_fatal(FARGS, "Invalid acceptor %d", j);
4473                                                                     }
4474                                                                     resdist = abs(top.atoms.atom[i].resind-
4475                                                                                   top.atoms.atom[j].resind);
4476                                                                     if (resdist >= max_hx)
4477                                                                     {
4478                                                                         resdist = max_hx-1;
4479                                                                     }
4480                                                                     __HBDATA->nhx[nframes][resdist]++;
4481                                                                 }
4482                                                             }
4483
4484                                                         }
4485                                                     } /* for aj  */
4486                                                 }     /* for xjj */
4487                                             }         /* for yjj */
4488                                         }             /* for zjj */
4489                                     }                 /* for ai  */
4490                                 }                     /* for grp */
4491                             }                         /* for xi,yi,zi */
4492                         }
4493                     }
4494                 } /* if (bSelected) {...} else */
4495
4496
4497                 /* Better wait for all threads to finnish using x[] before updating it. */
4498                 k = nframes;
4499 #pragma omp barrier
4500 #pragma omp critical
4501                 {
4502                     /* Sum up histograms and counts from p_hb[] into hb */
4503                     if (bOMP)
4504                     {
4505                         hb->nhb[k]   += p_hb[threadNr]->nhb[k];
4506                         hb->ndist[k] += p_hb[threadNr]->ndist[k];
4507                         for (j = 0; j < max_hx; j++)
4508                         {
4509                             hb->nhx[k][j]  += p_hb[threadNr]->nhx[k][j];
4510                         }
4511                     }
4512                 }
4513
4514                 /* Here are a handful of single constructs
4515                  * to share the workload a bit. The most
4516                  * important one is of course the last one,
4517                  * where there's a potential bottleneck in form
4518                  * of slow I/O.                    */
4519 #pragma omp barrier
4520 #pragma omp single
4521                 {
4522                     if (hb != NULL)
4523                     {
4524                         analyse_donor_props(opt2fn_null("-don", NFILE, fnm), hb, k, t, oenv);
4525                     }
4526                 }
4527
4528 #pragma omp single
4529                 {
4530                     if (fpnhb)
4531                     {
4532                         do_nhb_dist(fpnhb, hb, t);
4533                     }
4534                 }
4535             } /* if (bNN) {...} else  +   */
4536
4537 #pragma omp single
4538             {
4539                 trrStatus = (read_next_x(oenv, status, &t, x, box));
4540                 nframes++;
4541             }
4542
4543 #pragma omp barrier
4544         }
4545         while (trrStatus);
4546
4547         if (bOMP)
4548         {
4549 #pragma omp critical
4550             {
4551                 hb->nrhb   += p_hb[threadNr]->nrhb;
4552                 hb->nrdist += p_hb[threadNr]->nrdist;
4553             }
4554             /* Free parallel datastructures */
4555             sfree(p_hb[threadNr]->nhb);
4556             sfree(p_hb[threadNr]->ndist);
4557             sfree(p_hb[threadNr]->nhx);
4558
4559 #pragma omp for
4560             for (i = 0; i < nabin; i++)
4561             {
4562                 for (j = 0; j < actual_nThreads; j++)
4563                 {
4564
4565                     adist[i] += p_adist[j][i];
4566                 }
4567             }
4568 #pragma omp for
4569             for (i = 0; i <= nrbin; i++)
4570             {
4571                 for (j = 0; j < actual_nThreads; j++)
4572                 {
4573                     rdist[i] += p_rdist[j][i];
4574                 }
4575             }
4576
4577             sfree(p_adist[threadNr]);
4578             sfree(p_rdist[threadNr]);
4579         }
4580     } /* End of parallel region */
4581     if (bOMP)
4582     {
4583         sfree(p_adist);
4584         sfree(p_rdist);
4585     }
4586
4587     if (nframes < 2 && (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm)))
4588     {
4589         gmx_fatal(FARGS, "Cannot calculate autocorrelation of life times with less than two frames");
4590     }
4591
4592     free_grid(ngrid, &grid);
4593
4594     close_trj(status);
4595     if (fpnhb)
4596     {
4597         ffclose(fpnhb);
4598     }
4599
4600     /* Compute maximum possible number of different hbonds */
4601     if (maxnhb > 0)
4602     {
4603         max_nhb = maxnhb;
4604     }
4605     else
4606     {
4607         max_nhb = 0.5*(hb->d.nrd*hb->a.nra);
4608     }
4609     /* Added support for -contact below.
4610      * - Erik Marklund, May 29-31, 2006 */
4611     /* Changed contact code.
4612      * - Erik Marklund, June 29, 2006 */
4613     if (bHBmap && !bNN)
4614     {
4615         if (hb->nrhb == 0)
4616         {
4617             printf("No %s found!!\n", bContact ? "contacts" : "hydrogen bonds");
4618         }
4619         else
4620         {
4621             printf("Found %d different %s in trajectory\n"
4622                    "Found %d different atom-pairs within %s distance\n",
4623                    hb->nrhb, bContact ? "contacts" : "hydrogen bonds",
4624                    hb->nrdist, (r2cut > 0) ? "second cut-off" : "hydrogen bonding");
4625
4626             /*Control the pHist.*/
4627
4628             if (bMerge)
4629             {
4630                 merge_hb(hb, bTwo, bContact);
4631             }
4632
4633             if (opt2bSet("-hbn", NFILE, fnm))
4634             {
4635                 dump_hbmap(hb, NFILE, fnm, bTwo, bContact, isize, index, grpnames, &top.atoms);
4636             }
4637
4638             /* Moved the call to merge_hb() to a line BEFORE dump_hbmap
4639              * to make the -hbn and -hmb output match eachother.
4640              * - Erik Marklund, May 30, 2006 */
4641         }
4642     }
4643     /* Print out number of hbonds and distances */
4644     aver_nhb  = 0;
4645     aver_dist = 0;
4646     fp        = xvgropen(opt2fn("-num", NFILE, fnm), bContact ? "Contacts" :
4647                          "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
4648     snew(leg, 2);
4649     snew(leg[0], STRLEN);
4650     snew(leg[1], STRLEN);
4651     sprintf(leg[0], "%s", bContact ? "Contacts" : "Hydrogen bonds");
4652     sprintf(leg[1], "Pairs within %g nm", (r2cut > 0) ? r2cut : rcut);
4653     xvgr_legend(fp, 2, (const char**)leg, oenv);
4654     sfree(leg[1]);
4655     sfree(leg[0]);
4656     sfree(leg);
4657     for (i = 0; (i < nframes); i++)
4658     {
4659         fprintf(fp, "%10g  %10d  %10d\n", hb->time[i], hb->nhb[i], hb->ndist[i]);
4660         aver_nhb  += hb->nhb[i];
4661         aver_dist += hb->ndist[i];
4662     }
4663     ffclose(fp);
4664     aver_nhb  /= nframes;
4665     aver_dist /= nframes;
4666     /* Print HB distance distribution */
4667     if (opt2bSet("-dist", NFILE, fnm))
4668     {
4669         long sum;
4670
4671         sum = 0;
4672         for (i = 0; i < nrbin; i++)
4673         {
4674             sum += rdist[i];
4675         }
4676
4677         fp = xvgropen(opt2fn("-dist", NFILE, fnm),
4678                       "Hydrogen Bond Distribution",
4679                       bDA ?
4680                       "Donor - Acceptor Distance (nm)" :
4681                       "Hydrogen - Acceptor Distance (nm)", "", oenv);
4682         for (i = 0; i < nrbin; i++)
4683         {
4684             fprintf(fp, "%10g %10g\n", (i+0.5)*rbin, rdist[i]/(rbin*(real)sum));
4685         }
4686         ffclose(fp);
4687     }
4688
4689     /* Print HB angle distribution */
4690     if (opt2bSet("-ang", NFILE, fnm))
4691     {
4692         long sum;
4693
4694         sum = 0;
4695         for (i = 0; i < nabin; i++)
4696         {
4697             sum += adist[i];
4698         }
4699
4700         fp = xvgropen(opt2fn("-ang", NFILE, fnm),
4701                       "Hydrogen Bond Distribution",
4702                       "Hydrogen - Donor - Acceptor Angle (\\SO\\N)", "", oenv);
4703         for (i = 0; i < nabin; i++)
4704         {
4705             fprintf(fp, "%10g %10g\n", (i+0.5)*abin, adist[i]/(abin*(real)sum));
4706         }
4707         ffclose(fp);
4708     }
4709
4710     /* Print HB in alpha-helix */
4711     if (opt2bSet("-hx", NFILE, fnm))
4712     {
4713         fp = xvgropen(opt2fn("-hx", NFILE, fnm),
4714                       "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
4715         xvgr_legend(fp, NRHXTYPES, hxtypenames, oenv);
4716         for (i = 0; i < nframes; i++)
4717         {
4718             fprintf(fp, "%10g", hb->time[i]);
4719             for (j = 0; j < max_hx; j++)
4720             {
4721                 fprintf(fp, " %6d", hb->nhx[i][j]);
4722             }
4723             fprintf(fp, "\n");
4724         }
4725         ffclose(fp);
4726     }
4727     if (!bNN)
4728     {
4729         printf("Average number of %s per timeframe %.3f out of %g possible\n",
4730                bContact ? "contacts" : "hbonds",
4731                bContact ? aver_dist : aver_nhb, max_nhb);
4732     }
4733
4734     /* Do Autocorrelation etc. */
4735     if (hb->bHBmap)
4736     {
4737         /*
4738            Added support for -contact in ac and hbm calculations below.
4739            - Erik Marklund, May 29, 2006
4740          */
4741         ivec itmp;
4742         rvec rtmp;
4743         if (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm))
4744         {
4745             please_cite(stdout, "Spoel2006b");
4746         }
4747         if (opt2bSet("-ac", NFILE, fnm))
4748         {
4749             char *gemstring = NULL;
4750
4751             if (bGem || bNN)
4752             {
4753                 params = init_gemParams(rcut, D, hb->time, hb->nframes/2, nFitPoints, fit_start, fit_end,
4754                                         gemBallistic, nBalExp);
4755                 if (params == NULL)
4756                 {
4757                     gmx_fatal(FARGS, "Could not initiate t_gemParams params.");
4758                 }
4759             }
4760             gemstring = strdup(gemType[hb->per->gemtype]);
4761             do_hbac(opt2fn("-ac", NFILE, fnm), hb, nDump,
4762                     bMerge, bContact, fit_start, temp, r2cut > 0, smooth_tail_start, oenv,
4763                     gemstring, nThreads, NN, bBallistic, bGemFit);
4764         }
4765         if (opt2bSet("-life", NFILE, fnm))
4766         {
4767             do_hblife(opt2fn("-life", NFILE, fnm), hb, bMerge, bContact, oenv);
4768         }
4769         if (opt2bSet("-hbm", NFILE, fnm))
4770         {
4771             t_matrix mat;
4772             int      id, ia, hh, x, y;
4773
4774             if ((nframes > 0) && (hb->nrhb > 0))
4775             {
4776                 mat.nx = nframes;
4777                 mat.ny = hb->nrhb;
4778
4779                 snew(mat.matrix, mat.nx);
4780                 for (x = 0; (x < mat.nx); x++)
4781                 {
4782                     snew(mat.matrix[x], mat.ny);
4783                 }
4784                 y = 0;
4785                 for (id = 0; (id < hb->d.nrd); id++)
4786                 {
4787                     for (ia = 0; (ia < hb->a.nra); ia++)
4788                     {
4789                         for (hh = 0; (hh < hb->maxhydro); hh++)
4790                         {
4791                             if (hb->hbmap[id][ia])
4792                             {
4793                                 if (ISHB(hb->hbmap[id][ia]->history[hh]))
4794                                 {
4795                                     /* Changed '<' into '<=' in the for-statement below.
4796                                      * It fixed the previously undiscovered bug that caused
4797                                      * the last occurance of an hbond/contact to not be
4798                                      * set in mat.matrix. Have a look at any old -hbm-output
4799                                      * and you will notice that the last column is allways empty.
4800                                      * - Erik Marklund May 30, 2006
4801                                      */
4802                                     for (x = 0; (x <= hb->hbmap[id][ia]->nframes); x++)
4803                                     {
4804                                         int nn0 = hb->hbmap[id][ia]->n0;
4805                                         range_check(y, 0, mat.ny);
4806                                         mat.matrix[x+nn0][y] = is_hb(hb->hbmap[id][ia]->h[hh], x);
4807                                     }
4808                                     y++;
4809                                 }
4810                             }
4811                         }
4812                     }
4813                 }
4814                 mat.axis_x = hb->time;
4815                 snew(mat.axis_y, mat.ny);
4816                 for (j = 0; j < mat.ny; j++)
4817                 {
4818                     mat.axis_y[j] = j;
4819                 }
4820                 sprintf(mat.title, bContact ? "Contact Existence Map" :
4821                         "Hydrogen Bond Existence Map");
4822                 sprintf(mat.legend, bContact ? "Contacts" : "Hydrogen Bonds");
4823                 sprintf(mat.label_x, "%s", output_env_get_xvgr_tlabel(oenv));
4824                 sprintf(mat.label_y, bContact ? "Contact Index" : "Hydrogen Bond Index");
4825                 mat.bDiscrete = TRUE;
4826                 mat.nmap      = 2;
4827                 snew(mat.map, mat.nmap);
4828                 for (i = 0; i < mat.nmap; i++)
4829                 {
4830                     mat.map[i].code.c1 = hbmap[i];
4831                     mat.map[i].desc    = hbdesc[i];
4832                     mat.map[i].rgb     = hbrgb[i];
4833                 }
4834                 fp = opt2FILE("-hbm", NFILE, fnm, "w");
4835                 write_xpm_m(fp, mat);
4836                 ffclose(fp);
4837                 for (x = 0; x < mat.nx; x++)
4838                 {
4839                     sfree(mat.matrix[x]);
4840                 }
4841                 sfree(mat.axis_y);
4842                 sfree(mat.matrix);
4843                 sfree(mat.map);
4844             }
4845             else
4846             {
4847                 fprintf(stderr, "No hydrogen bonds/contacts found. No hydrogen bond map will be printed.\n");
4848             }
4849         }
4850     }
4851
4852     if (bGem)
4853     {
4854         fprintf(stderr, "There were %i periodic shifts\n", hb->per->nper);
4855         fprintf(stderr, "Freeing pHist for all donors...\n");
4856         for (i = 0; i < hb->d.nrd; i++)
4857         {
4858             fprintf(stderr, "\r%i", i);
4859             if (hb->per->pHist[i] != NULL)
4860             {
4861                 for (j = 0; j < hb->a.nra; j++)
4862                 {
4863                     clearPshift(&(hb->per->pHist[i][j]));
4864                 }
4865                 sfree(hb->per->pHist[i]);
4866             }
4867         }
4868         sfree(hb->per->pHist);
4869         sfree(hb->per->p2i);
4870         sfree(hb->per);
4871         fprintf(stderr, "...done.\n");
4872     }
4873
4874 #ifdef HAVE_NN_LOOPS
4875     if (bNN)
4876     {
4877         free_hbEmap(hb);
4878     }
4879 #endif
4880
4881     if (hb->bDAnr)
4882     {
4883         int    i, j, nleg;
4884         char **legnames;
4885         char   buf[STRLEN];
4886
4887 #define USE_THIS_GROUP(j) ( (j == gr0) || (bTwo && (j == gr1)) )
4888
4889         fp = xvgropen(opt2fn("-dan", NFILE, fnm),
4890                       "Donors and Acceptors", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
4891         nleg = (bTwo ? 2 : 1)*2;
4892         snew(legnames, nleg);
4893         i = 0;
4894         for (j = 0; j < grNR; j++)
4895         {
4896             if (USE_THIS_GROUP(j) )
4897             {
4898                 sprintf(buf, "Donors %s", grpnames[j]);
4899                 legnames[i++] = strdup(buf);
4900                 sprintf(buf, "Acceptors %s", grpnames[j]);
4901                 legnames[i++] = strdup(buf);
4902             }
4903         }
4904         if (i != nleg)
4905         {
4906             gmx_incons("number of legend entries");
4907         }
4908         xvgr_legend(fp, nleg, (const char**)legnames, oenv);
4909         for (i = 0; i < nframes; i++)
4910         {
4911             fprintf(fp, "%10g", hb->time[i]);
4912             for (j = 0; (j < grNR); j++)
4913             {
4914                 if (USE_THIS_GROUP(j) )
4915                 {
4916                     fprintf(fp, " %6d", hb->danr[i][j]);
4917                 }
4918             }
4919             fprintf(fp, "\n");
4920         }
4921         ffclose(fp);
4922     }
4923
4924     return 0;
4925 }