src/gromacs/gmxana/gmx_hbond.c

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