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