src/gromacs/gmxana/gmx_hbond.c

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