src/gromacs/gmxana/gmx_hbond.c

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