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