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 void pbc_in_gridbox(rvec dx, matrix box);
1421
1422 static void build_grid(t_hbdata *hb, rvec x[], rvec xshell,
1423                        gmx_bool bBox, matrix box, rvec hbox,
1424                        real rcut, real rshell,
1425                        ivec ngrid, t_gridcell ***grid)
1426 {
1427     int         i, m, gr, xi, yi, zi, nr;
1428     atom_id    *ad;
1429     ivec        grididx;
1430     rvec        invdelta, dshell, xtemp = {0, 0, 0};
1431     t_ncell    *newgrid;
1432     gmx_bool    bDoRshell, bInShell, bAcc;
1433     real        rshell2 = 0;
1434     int         gx, gy, gz;
1435     int         dum = -1;
1436
1437     bDoRshell = (rshell > 0);
1438     rshell2   = sqr(rshell);
1439     bInShell  = TRUE;
1440
1441 #define DBB(x) if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n", __LINE__,#x, x)
1442     DBB(dum);
1443     for (m = 0; m < DIM; m++)
1444     {
1445         hbox[m] = box[m][m]*0.5;
1446         if (bBox)
1447         {
1448             invdelta[m] = ngrid[m]/box[m][m];
1449             if (1/invdelta[m] < rcut)
1450             {
1451                 gmx_fatal(FARGS, "Your computational box has shrunk too much.\n"
1452                           "%s can not handle this situation, sorry.\n",
1453                           ShortProgram());
1454             }
1455         }
1456         else
1457         {
1458             invdelta[m] = 0;
1459         }
1460     }
1461     grididx[XX] = 0;
1462     grididx[YY] = 0;
1463     grididx[ZZ] = 0;
1464     DBB(dum);
1465     /* resetting atom counts */
1466     for (gr = 0; (gr < grNR); gr++)
1467     {
1468         for (zi = 0; zi < ngrid[ZZ]; zi++)
1469         {
1470             for (yi = 0; yi < ngrid[YY]; yi++)
1471             {
1472                 for (xi = 0; xi < ngrid[XX]; xi++)
1473                 {
1474                     grid[zi][yi][xi].d[gr].nr = 0;
1475                     grid[zi][yi][xi].a[gr].nr = 0;
1476                 }
1477             }
1478         }
1479         DBB(dum);
1480
1481         /* put atoms in grid cells */
1482         for (bAcc = FALSE; (bAcc <= TRUE); bAcc++)
1483         {
1484             if (bAcc)
1485             {
1486                 nr = hb->a.nra;
1487                 ad = hb->a.acc;
1488             }
1489             else
1490             {
1491                 nr = hb->d.nrd;
1492                 ad = hb->d.don;
1493             }
1494             DBB(bAcc);
1495             for (i = 0; (i < nr); i++)
1496             {
1497                 /* check if we are inside the shell */
1498                 /* if bDoRshell=FALSE then bInShell=TRUE always */
1499                 DBB(i);
1500                 if (bDoRshell)
1501                 {
1502                     bInShell = TRUE;
1503                     rvec_sub(x[ad[i]], xshell, dshell);
1504                     if (bBox)
1505                     {
1506                         if (FALSE && !hb->bGem)
1507                         {
1508                             for (m = DIM-1; m >= 0 && bInShell; m--)
1509                             {
1510                                 if (dshell[m] < -hbox[m])
1511                                 {
1512                                     rvec_inc(dshell, box[m]);
1513                                 }
1514                                 else if (dshell[m] >= hbox[m])
1515                                 {
1516                                     dshell[m] -= 2*hbox[m];
1517                                 }
1518                                 /* if we're outside the cube, we're outside the sphere also! */
1519                                 if ( (dshell[m] > rshell) || (-dshell[m] > rshell) )
1520                                 {
1521                                     bInShell = FALSE;
1522                                 }
1523                             }
1524                         }
1525                         else
1526                         {
1527                             gmx_bool bDone = FALSE;
1528                             while (!bDone)
1529                             {
1530                                 bDone = TRUE;
1531                                 for (m = DIM-1; m >= 0 && bInShell; m--)
1532                                 {
1533                                     if (dshell[m] < -hbox[m])
1534                                     {
1535                                         bDone = FALSE;
1536                                         rvec_inc(dshell, box[m]);
1537                                     }
1538                                     if (dshell[m] >= hbox[m])
1539                                     {
1540                                         bDone      = FALSE;
1541                                         dshell[m] -= 2*hbox[m];
1542                                     }
1543                                 }
1544                             }
1545                             for (m = DIM-1; m >= 0 && bInShell; m--)
1546                             {
1547                                 /* if we're outside the cube, we're outside the sphere also! */
1548                                 if ( (dshell[m] > rshell) || (-dshell[m] > rshell) )
1549                                 {
1550                                     bInShell = FALSE;
1551                                 }
1552                             }
1553                         }
1554                     }
1555                     /* if we're inside the cube, check if we're inside the sphere */
1556                     if (bInShell)
1557                     {
1558                         bInShell = norm2(dshell) < rshell2;
1559                     }
1560                 }
1561                 DBB(i);
1562                 if (bInShell)
1563                 {
1564                     if (bBox)
1565                     {
1566                         if (hb->bGem)
1567                         {
1568                             copy_rvec(x[ad[i]], xtemp);
1569                         }
1570                         pbc_in_gridbox(x[ad[i]], box);
1571
1572                         for (m = DIM-1; m >= 0; m--)
1573                         {   /* determine grid index of atom */
1574                             grididx[m] = x[ad[i]][m]*invdelta[m];
1575                             grididx[m] = (grididx[m]+ngrid[m]) % ngrid[m];
1576                         }
1577                         if (hb->bGem)
1578                         {
1579                             copy_rvec(xtemp, x[ad[i]]); /* copy back */
1580                         }
1581                     }
1582
1583                     gx = grididx[XX];
1584                     gy = grididx[YY];
1585                     gz = grididx[ZZ];
1586                     range_check(gx, 0, ngrid[XX]);
1587                     range_check(gy, 0, ngrid[YY]);
1588                     range_check(gz, 0, ngrid[ZZ]);
1589                     DBB(gx);
1590                     DBB(gy);
1591                     DBB(gz);
1592                     /* add atom to grid cell */
1593                     if (bAcc)
1594                     {
1595                         newgrid = &(grid[gz][gy][gx].a[gr]);
1596                     }
1597                     else
1598                     {
1599                         newgrid = &(grid[gz][gy][gx].d[gr]);
1600                     }
1601                     if (newgrid->nr >= newgrid->maxnr)
1602                     {
1603                         newgrid->maxnr += 10;
1604                         DBB(newgrid->maxnr);
1605                         srenew(newgrid->atoms, newgrid->maxnr);
1606                     }
1607                     DBB(newgrid->nr);
1608                     newgrid->atoms[newgrid->nr] = ad[i];
1609                     newgrid->nr++;
1610                 }
1611             }
1612         }
1613     }
1614 }
1615
1616 static void count_da_grid(ivec ngrid, t_gridcell ***grid, t_icell danr)
1617 {
1618     int gr, xi, yi, zi;
1619
1620     for (gr = 0; (gr < grNR); gr++)
1621     {
1622         danr[gr] = 0;
1623         for (zi = 0; zi < ngrid[ZZ]; zi++)
1624         {
1625             for (yi = 0; yi < ngrid[YY]; yi++)
1626             {
1627                 for (xi = 0; xi < ngrid[XX]; xi++)
1628                 {
1629                     danr[gr] += grid[zi][yi][xi].d[gr].nr;
1630                 }
1631             }
1632         }
1633     }
1634 }
1635
1636 /* The grid loop.
1637  * Without a box, the grid is 1x1x1, so all loops are 1 long.
1638  * With a rectangular box (bTric==FALSE) all loops are 3 long.
1639  * With a triclinic box all loops are 3 long, except when a cell is
1640  * located next to one of the box edges which is not parallel to the
1641  * x/y-plane, in that case all cells in a line or layer are searched.
1642  * This could be implemented slightly more efficient, but the code
1643  * would get much more complicated.
1644  */
1645 static gmx_inline gmx_bool grid_loop_begin(int n, int x, gmx_bool bTric, gmx_bool bEdge)
1646 {
1647     return ((n == 1) ? x : bTric && bEdge ? 0     : (x-1));
1648 }
1649 static gmx_inline gmx_bool grid_loop_end(int n, int x, gmx_bool bTric, gmx_bool bEdge)
1650 {
1651     return ((n == 1) ? x : bTric && bEdge ? (n-1) : (x+1));
1652 }
1653 static gmx_inline int grid_mod(int j, int n)
1654 {
1655     return (j+n) % (n);
1656 }
1657
1658 static void dump_grid(FILE *fp, ivec ngrid, t_gridcell ***grid)
1659 {
1660     int gr, x, y, z, sum[grNR];
1661
1662     fprintf(fp, "grid %dx%dx%d\n", ngrid[XX], ngrid[YY], ngrid[ZZ]);
1663     for (gr = 0; gr < grNR; gr++)
1664     {
1665         sum[gr] = 0;
1666         fprintf(fp, "GROUP %d (%s)\n", gr, grpnames[gr]);
1667         for (z = 0; z < ngrid[ZZ]; z += 2)
1668         {
1669             fprintf(fp, "Z=%d,%d\n", z, z+1);
1670             for (y = 0; y < ngrid[YY]; y++)
1671             {
1672                 for (x = 0; x < ngrid[XX]; x++)
1673                 {
1674                     fprintf(fp, "%3d", grid[x][y][z].d[gr].nr);
1675                     sum[gr] += grid[z][y][x].d[gr].nr;
1676                     fprintf(fp, "%3d", grid[x][y][z].a[gr].nr);
1677                     sum[gr] += grid[z][y][x].a[gr].nr;
1678
1679                 }
1680                 fprintf(fp, " | ");
1681                 if ( (z+1) < ngrid[ZZ])
1682                 {
1683                     for (x = 0; x < ngrid[XX]; x++)
1684                     {
1685                         fprintf(fp, "%3d", grid[z+1][y][x].d[gr].nr);
1686                         sum[gr] += grid[z+1][y][x].d[gr].nr;
1687                         fprintf(fp, "%3d", grid[z+1][y][x].a[gr].nr);
1688                         sum[gr] += grid[z+1][y][x].a[gr].nr;
1689                     }
1690                 }
1691                 fprintf(fp, "\n");
1692             }
1693         }
1694     }
1695     fprintf(fp, "TOTALS:");
1696     for (gr = 0; gr < grNR; gr++)
1697     {
1698         fprintf(fp, " %d=%d", gr, sum[gr]);
1699     }
1700     fprintf(fp, "\n");
1701 }
1702
1703 /* New GMX record! 5 * in a row. Congratulations!
1704  * Sorry, only four left.
1705  */
1706 static void free_grid(ivec ngrid, t_gridcell ****grid)
1707 {
1708     int           y, z;
1709     t_gridcell ***g = *grid;
1710
1711     for (z = 0; z < ngrid[ZZ]; z++)
1712     {
1713         for (y = 0; y < ngrid[YY]; y++)
1714         {
1715             sfree(g[z][y]);
1716         }
1717         sfree(g[z]);
1718     }
1719     sfree(g);
1720     g = NULL;
1721 }
1722
1723 void pbc_correct_gem(rvec dx, matrix box, rvec hbox)
1724 {
1725     int      m;
1726     gmx_bool bDone = FALSE;
1727     while (!bDone)
1728     {
1729         bDone = TRUE;
1730         for (m = DIM-1; m >= 0; m--)
1731         {
1732             if (dx[m] < -hbox[m])
1733             {
1734                 bDone = FALSE;
1735                 rvec_inc(dx, box[m]);
1736             }
1737             if (dx[m] >= hbox[m])
1738             {
1739                 bDone = FALSE;
1740                 rvec_dec(dx, box[m]);
1741             }
1742         }
1743     }
1744 }
1745
1746 void pbc_in_gridbox(rvec dx, matrix box)
1747 {
1748     int      m;
1749     gmx_bool bDone = FALSE;
1750     while (!bDone)
1751     {
1752         bDone = TRUE;
1753         for (m = DIM-1; m >= 0; m--)
1754         {
1755             if (dx[m] < 0)
1756             {
1757                 bDone = FALSE;
1758                 rvec_inc(dx, box[m]);
1759             }
1760             if (dx[m] >= box[m][m])
1761             {
1762                 bDone = FALSE;
1763                 rvec_dec(dx, box[m]);
1764             }
1765         }
1766     }
1767 }
1768
1769 /* Added argument r2cut, changed contact and implemented
1770  * use of second cut-off.
1771  * - Erik Marklund, June 29, 2006
1772  */
1773 static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a,
1774                     real rcut, real r2cut, real ccut,
1775                     rvec x[], gmx_bool bBox, matrix box, rvec hbox,
1776                     real *d_ha, real *ang, gmx_bool bDA, int *hhh,
1777                     gmx_bool bContact, gmx_bool bMerge, PSTYPE *p)
1778 {
1779     int      h, hh, id, ja, ihb;
1780     rvec     r_da, r_ha, r_dh, r = {0, 0, 0};
1781     ivec     ri;
1782     real     rc2, r2c2, rda2, rha2, ca;
1783     gmx_bool HAinrange = FALSE; /* If !bDA. Needed for returning hbDist in a correct way. */
1784     gmx_bool daSwap    = FALSE;
1785
1786     if (d == a)
1787     {
1788         return hbNo;
1789     }
1790
1791     if (((id = donor_index(&hb->d, grpd, d)) == NOTSET) ||
1792         ((ja = acceptor_index(&hb->a, grpa, a)) == NOTSET))
1793     {
1794         return hbNo;
1795     }
1796
1797     rc2  = rcut*rcut;
1798     r2c2 = r2cut*r2cut;
1799
1800     rvec_sub(x[d], x[a], r_da);
1801     /* Insert projection code here */
1802
1803     if (bMerge && d > a && isInterchangable(hb, d, a, grpd, grpa))
1804     {
1805         /* Then this hbond/contact will be found again, or it has already been found. */
1806         /*return hbNo;*/
1807     }
1808     if (bBox)
1809     {
1810         if (d > a && bMerge && isInterchangable(hb, d, a, grpd, grpa)) /* acceptor is also a donor and vice versa? */
1811         {                                                              /* return hbNo; */
1812             daSwap = TRUE;                                             /* If so, then their history should be filed with donor and acceptor swapped. */
1813         }
1814         if (hb->bGem)
1815         {
1816             copy_rvec(r_da, r); /* Save this for later */
1817             pbc_correct_gem(r_da, box, hbox);
1818         }
1819         else
1820         {
1821             pbc_correct_gem(r_da, box, hbox);
1822         }
1823     }
1824     rda2 = iprod(r_da, r_da);
1825
1826     if (bContact)
1827     {
1828         if (daSwap && grpa == grpd)
1829         {
1830             return hbNo;
1831         }
1832         if (rda2 <= rc2)
1833         {
1834             if (hb->bGem)
1835             {
1836                 calcBoxDistance(hb->per->P, r, ri);
1837                 *p = periodicIndex(ri, hb->per, daSwap);    /* find (or add) periodicity index. */
1838             }
1839             return hbHB;
1840         }
1841         else if (rda2 < r2c2)
1842         {
1843             return hbDist;
1844         }
1845         else
1846         {
1847             return hbNo;
1848         }
1849     }
1850     *hhh = NOTSET;
1851
1852     if (bDA && (rda2 > rc2))
1853     {
1854         return hbNo;
1855     }
1856
1857     for (h = 0; (h < hb->d.nhydro[id]); h++)
1858     {
1859         hh   = hb->d.hydro[id][h];
1860         rha2 = rc2+1;
1861         if (!bDA)
1862         {
1863             rvec_sub(x[hh], x[a], r_ha);
1864             if (bBox)
1865             {
1866                 pbc_correct_gem(r_ha, box, hbox);
1867             }
1868             rha2 = iprod(r_ha, r_ha);
1869         }
1870
1871         if (hb->bGem)
1872         {
1873             calcBoxDistance(hb->per->P, r, ri);
1874             *p = periodicIndex(ri, hb->per, daSwap);    /* find periodicity index. */
1875         }
1876
1877         if (bDA || (!bDA && (rha2 <= rc2)))
1878         {
1879             rvec_sub(x[d], x[hh], r_dh);
1880             if (bBox)
1881             {
1882                 pbc_correct_gem(r_dh, box, hbox);
1883             }
1884
1885             if (!bDA)
1886             {
1887                 HAinrange = TRUE;
1888             }
1889             ca = cos_angle(r_dh, r_da);
1890             /* if angle is smaller, cos is larger */
1891             if (ca >= ccut)
1892             {
1893                 *hhh  = hh;
1894                 *d_ha = sqrt(bDA ? rda2 : rha2);
1895                 *ang  = acos(ca);
1896                 return hbHB;
1897             }
1898         }
1899     }
1900     if (bDA || (!bDA && HAinrange))
1901     {
1902         return hbDist;
1903     }
1904     else
1905     {
1906         return hbNo;
1907     }
1908 }
1909
1910 /* Fixed previously undiscovered bug in the merge
1911    code, where the last frame of each hbond disappears.
1912    - Erik Marklund, June 1, 2006 */
1913 /* Added the following arguments:
1914  *   ptmp[] - temporary periodicity hisory
1915  *   a1     - identity of first acceptor/donor
1916  *   a2     - identity of second acceptor/donor
1917  * - Erik Marklund, FEB 20 2010 */
1918
1919 /* Merging is now done on the fly, so do_merge is most likely obsolete now.
1920  * Will do some more testing before removing the function entirely.
1921  * - Erik Marklund, MAY 10 2010 */
1922 static void do_merge(t_hbdata *hb, int ntmp,
1923                      unsigned int htmp[], unsigned int gtmp[], PSTYPE ptmp[],
1924                      t_hbond *hb0, t_hbond *hb1, int a1, int a2)
1925 {
1926     /* Here we need to make sure we're treating periodicity in
1927      * the right way for the geminate recombination kinetics. */
1928
1929     int       m, mm, n00, n01, nn0, nnframes;
1930     PSTYPE    pm;
1931     t_pShift *pShift;
1932
1933     /* Decide where to start from when merging */
1934     n00      = hb0->n0;
1935     n01      = hb1->n0;
1936     nn0      = min(n00, n01);
1937     nnframes = max(n00 + hb0->nframes, n01 + hb1->nframes) - nn0;
1938     /* Initiate tmp arrays */
1939     for (m = 0; (m < ntmp); m++)
1940     {
1941         htmp[m] = 0;
1942         gtmp[m] = 0;
1943         ptmp[m] = 0;
1944     }
1945     /* Fill tmp arrays with values due to first HB */
1946     /* Once again '<' had to be replaced with '<='
1947        to catch the last frame in which the hbond
1948        appears.
1949        - Erik Marklund, June 1, 2006 */
1950     for (m = 0; (m <= hb0->nframes); m++)
1951     {
1952         mm       = m+n00-nn0;
1953         htmp[mm] = is_hb(hb0->h[0], m);
1954         if (hb->bGem)
1955         {
1956             pm = getPshift(hb->per->pHist[a1][a2], m+hb0->n0);
1957             if (pm > hb->per->nper)
1958             {
1959                 gmx_fatal(FARGS, "Illegal shift!");
1960             }
1961             else
1962             {
1963                 ptmp[mm] = pm; /*hb->per->pHist[a1][a2][m];*/
1964             }
1965         }
1966     }
1967     /* If we're doing geminate recompbination we usually don't need the distances.
1968      * Let's save some memory and time. */
1969     if (TRUE || !hb->bGem || hb->per->gemtype == gemAD)
1970     {
1971         for (m = 0; (m <= hb0->nframes); m++)
1972         {
1973             mm       = m+n00-nn0;
1974             gtmp[mm] = is_hb(hb0->g[0], m);
1975         }
1976     }
1977     /* Next HB */
1978     for (m = 0; (m <= hb1->nframes); m++)
1979     {
1980         mm       = m+n01-nn0;
1981         htmp[mm] = htmp[mm] || is_hb(hb1->h[0], m);
1982         gtmp[mm] = gtmp[mm] || is_hb(hb1->g[0], m);
1983         if (hb->bGem /* && ptmp[mm] != 0 */)
1984         {
1985
1986             /* If this hbond has been seen before with donor and acceptor swapped,
1987              * then we need to find the mirrored (*-1) periodicity vector to truely
1988              * merge the hbond history. */
1989             pm = findMirror(getPshift(hb->per->pHist[a2][a1], m+hb1->n0), hb->per->p2i, hb->per->nper);
1990             /* Store index of mirror */
1991             if (pm > hb->per->nper)
1992             {
1993                 gmx_fatal(FARGS, "Illegal shift!");
1994             }
1995             ptmp[mm] = pm;
1996         }
1997     }
1998     /* Reallocate target array */
1999     if (nnframes > hb0->maxframes)
2000     {
2001         srenew(hb0->h[0], 4+nnframes/hb->wordlen);
2002         srenew(hb0->g[0], 4+nnframes/hb->wordlen);
2003     }
2004     if (NULL != hb->per->pHist)
2005     {
2006         clearPshift(&(hb->per->pHist[a1][a2]));
2007     }
2008
2009     /* Copy temp array to target array */
2010     for (m = 0; (m <= nnframes); m++)
2011     {
2012         _set_hb(hb0->h[0], m, htmp[m]);
2013         _set_hb(hb0->g[0], m, gtmp[m]);
2014         if (hb->bGem)
2015         {
2016             addPshift(&(hb->per->pHist[a1][a2]), ptmp[m], m+nn0);
2017         }
2018     }
2019
2020     /* Set scalar variables */
2021     hb0->n0        = nn0;
2022     hb0->maxframes = nnframes;
2023 }
2024
2025 /* Added argument bContact for nicer output.
2026  * Erik Marklund, June 29, 2006
2027  */
2028 static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact)
2029 {
2030     int           i, inrnew, indnew, j, ii, jj, m, id, ia, grp, ogrp, ntmp;
2031     unsigned int *htmp, *gtmp;
2032     PSTYPE       *ptmp;
2033     t_hbond      *hb0, *hb1;
2034
2035     inrnew = hb->nrhb;
2036     indnew = hb->nrdist;
2037
2038     /* Check whether donors are also acceptors */
2039     printf("Merging hbonds with Acceptor and Donor swapped\n");
2040
2041     ntmp = 2*hb->max_frames;
2042     snew(gtmp, ntmp);
2043     snew(htmp, ntmp);
2044     snew(ptmp, ntmp);
2045     for (i = 0; (i < hb->d.nrd); i++)
2046     {
2047         fprintf(stderr, "\r%d/%d", i+1, hb->d.nrd);
2048         id = hb->d.don[i];
2049         ii = hb->a.aptr[id];
2050         for (j = 0; (j < hb->a.nra); j++)
2051         {
2052             ia = hb->a.acc[j];
2053             jj = hb->d.dptr[ia];
2054             if ((id != ia) && (ii != NOTSET) && (jj != NOTSET) &&
2055                 (!bTwo || (bTwo && (hb->d.grp[i] != hb->a.grp[j]))))
2056             {
2057                 hb0 = hb->hbmap[i][j];
2058                 hb1 = hb->hbmap[jj][ii];
2059                 if (hb0 && hb1 && ISHB(hb0->history[0]) && ISHB(hb1->history[0]))
2060                 {
2061                     do_merge(hb, ntmp, htmp, gtmp, ptmp, hb0, hb1, i, j);
2062                     if (ISHB(hb1->history[0]))
2063                     {
2064                         inrnew--;
2065                     }
2066                     else if (ISDIST(hb1->history[0]))
2067                     {
2068                         indnew--;
2069                     }
2070                     else
2071                     if (bContact)
2072                     {
2073                         gmx_incons("No contact history");
2074                     }
2075                     else
2076                     {
2077                         gmx_incons("Neither hydrogen bond nor distance");
2078                     }
2079                     sfree(hb1->h[0]);
2080                     sfree(hb1->g[0]);
2081                     if (hb->bGem)
2082                     {
2083                         clearPshift(&(hb->per->pHist[jj][ii]));
2084                     }
2085                     hb1->h[0]       = NULL;
2086                     hb1->g[0]       = NULL;
2087                     hb1->history[0] = hbNo;
2088                 }
2089             }
2090         }
2091     }
2092     fprintf(stderr, "\n");
2093     printf("- Reduced number of hbonds from %d to %d\n", hb->nrhb, inrnew);
2094     printf("- Reduced number of distances from %d to %d\n", hb->nrdist, indnew);
2095     hb->nrhb   = inrnew;
2096     hb->nrdist = indnew;
2097     sfree(gtmp);
2098     sfree(htmp);
2099     sfree(ptmp);
2100 }
2101
2102 static void do_nhb_dist(FILE *fp, t_hbdata *hb, real t)
2103 {
2104     int  i, j, k, n_bound[MAXHH], nbtot;
2105     h_id nhb;
2106
2107
2108     /* Set array to 0 */
2109     for (k = 0; (k < MAXHH); k++)
2110     {
2111         n_bound[k] = 0;
2112     }
2113     /* Loop over possible donors */
2114     for (i = 0; (i < hb->d.nrd); i++)
2115     {
2116         for (j = 0; (j < hb->d.nhydro[i]); j++)
2117         {
2118             n_bound[hb->d.nhbonds[i][j]]++;
2119         }
2120     }
2121     fprintf(fp, "%12.5e", t);
2122     nbtot = 0;
2123     for (k = 0; (k < MAXHH); k++)
2124     {
2125         fprintf(fp, "  %8d", n_bound[k]);
2126         nbtot += n_bound[k]*k;
2127     }
2128     fprintf(fp, "  %8d\n", nbtot);
2129 }
2130
2131 /* Added argument bContact in do_hblife(...). Also
2132  * added support for -contact in function body.
2133  * - Erik Marklund, May 31, 2006 */
2134 /* Changed the contact code slightly.
2135  * - Erik Marklund, June 29, 2006
2136  */
2137 static void do_hblife(const char *fn, t_hbdata *hb, gmx_bool bMerge, gmx_bool bContact,
2138                       const output_env_t oenv)
2139 {
2140     FILE          *fp;
2141     const char    *leg[] = { "p(t)", "t p(t)" };
2142     int           *histo;
2143     int            i, j, j0, k, m, nh, ihb, ohb, nhydro, ndump = 0;
2144     int            nframes = hb->nframes;
2145     unsigned int **h;
2146     real           t, x1, dt;
2147     double         sum, integral;
2148     t_hbond       *hbh;
2149
2150     snew(h, hb->maxhydro);
2151     snew(histo, nframes+1);
2152     /* Total number of hbonds analyzed here */
2153     for (i = 0; (i < hb->d.nrd); i++)
2154     {
2155         for (k = 0; (k < hb->a.nra); k++)
2156         {
2157             hbh = hb->hbmap[i][k];
2158             if (hbh)
2159             {
2160                 if (bMerge)
2161                 {
2162                     if (hbh->h[0])
2163                     {
2164                         h[0]   = hbh->h[0];
2165                         nhydro = 1;
2166                     }
2167                     else
2168                     {
2169                         nhydro = 0;
2170                     }
2171                 }
2172                 else
2173                 {
2174                     nhydro = 0;
2175                     for (m = 0; (m < hb->maxhydro); m++)
2176                     {
2177                         if (hbh->h[m])
2178                         {
2179                             h[nhydro++] = bContact ? hbh->g[m] : hbh->h[m];
2180                         }
2181                     }
2182                 }
2183                 for (nh = 0; (nh < nhydro); nh++)
2184                 {
2185                     ohb = 0;
2186                     j0  = 0;
2187
2188                     /* Changed '<' into '<=' below, just like I
2189                        did in the hbm-output-loop in the main code.
2190                        - Erik Marklund, May 31, 2006
2191                      */
2192                     for (j = 0; (j <= hbh->nframes); j++)
2193                     {
2194                         ihb      = is_hb(h[nh], j);
2195                         if (debug && (ndump < 10))
2196                         {
2197                             fprintf(debug, "%5d  %5d\n", j, ihb);
2198                         }
2199                         if (ihb != ohb)
2200                         {
2201                             if (ihb)
2202                             {
2203                                 j0 = j;
2204                             }
2205                             else
2206                             {
2207                                 histo[j-j0]++;
2208                             }
2209                             ohb = ihb;
2210                         }
2211                     }
2212                     ndump++;
2213                 }
2214             }
2215         }
2216     }
2217     fprintf(stderr, "\n");
2218     if (bContact)
2219     {
2220         fp = xvgropen(fn, "Uninterrupted contact lifetime", output_env_get_xvgr_tlabel(oenv), "()", oenv);
2221     }
2222     else
2223     {
2224         fp = xvgropen(fn, "Uninterrupted hydrogen bond lifetime", output_env_get_xvgr_tlabel(oenv), "()",
2225                       oenv);
2226     }
2227
2228     xvgr_legend(fp, asize(leg), leg, oenv);
2229     j0 = nframes-1;
2230     while ((j0 > 0) && (histo[j0] == 0))
2231     {
2232         j0--;
2233     }
2234     sum = 0;
2235     for (i = 0; (i <= j0); i++)
2236     {
2237         sum += histo[i];
2238     }
2239     dt       = hb->time[1]-hb->time[0];
2240     sum      = dt*sum;
2241     integral = 0;
2242     for (i = 1; (i <= j0); i++)
2243     {
2244         t  = hb->time[i] - hb->time[0] - 0.5*dt;
2245         x1 = t*histo[i]/sum;
2246         fprintf(fp, "%8.3f  %10.3e  %10.3e\n", t, histo[i]/sum, x1);
2247         integral += x1;
2248     }
2249     integral *= dt;
2250     xvgrclose(fp);
2251     printf("%s lifetime = %.2f ps\n", bContact ? "Contact" : "HB", integral);
2252     printf("Note that the lifetime obtained in this manner is close to useless\n");
2253     printf("Use the -ac option instead and check the Forward lifetime\n");
2254     please_cite(stdout, "Spoel2006b");
2255     sfree(h);
2256     sfree(histo);
2257 }
2258
2259 /* Changed argument bMerge into oneHB to handle contacts properly.
2260  * - Erik Marklund, June 29, 2006
2261  */
2262 static void dump_ac(t_hbdata *hb, gmx_bool oneHB, int nDump)
2263 {
2264     FILE     *fp;
2265     int       i, j, k, m, nd, ihb, idist;
2266     int       nframes = hb->nframes;
2267     gmx_bool  bPrint;
2268     t_hbond  *hbh;
2269
2270     if (nDump <= 0)
2271     {
2272         return;
2273     }
2274     fp = gmx_ffopen("debug-ac.xvg", "w");
2275     for (j = 0; (j < nframes); j++)
2276     {
2277         fprintf(fp, "%10.3f", hb->time[j]);
2278         for (i = nd = 0; (i < hb->d.nrd) && (nd < nDump); i++)
2279         {
2280             for (k = 0; (k < hb->a.nra) && (nd < nDump); k++)
2281             {
2282                 bPrint = FALSE;
2283                 ihb    = idist = 0;
2284                 hbh    = hb->hbmap[i][k];
2285                 if (oneHB)
2286                 {
2287                     if (hbh->h[0])
2288                     {
2289                         ihb    = is_hb(hbh->h[0], j);
2290                         idist  = is_hb(hbh->g[0], j);
2291                         bPrint = TRUE;
2292                     }
2293                 }
2294                 else
2295                 {
2296                     for (m = 0; (m < hb->maxhydro) && !ihb; m++)
2297                     {
2298                         ihb   = ihb   || ((hbh->h[m]) && is_hb(hbh->h[m], j));
2299                         idist = idist || ((hbh->g[m]) && is_hb(hbh->g[m], j));
2300                     }
2301                     /* This is not correct! */
2302                     /* What isn't correct? -Erik M */
2303                     bPrint = TRUE;
2304                 }
2305                 if (bPrint)
2306                 {
2307                     fprintf(fp, "  %1d-%1d", ihb, idist);
2308                     nd++;
2309                 }
2310             }
2311         }
2312         fprintf(fp, "\n");
2313     }
2314     gmx_ffclose(fp);
2315 }
2316
2317 static real calc_dg(real tau, real temp)
2318 {
2319     real kbt;
2320
2321     kbt = BOLTZ*temp;
2322     if (tau <= 0)
2323     {
2324         return -666;
2325     }
2326     else
2327     {
2328         return kbt*log(kbt*tau/PLANCK);
2329     }
2330 }
2331
2332 typedef struct {
2333     int   n0, n1, nparams, ndelta;
2334     real  kkk[2];
2335     real *t, *ct, *nt, *kt, *sigma_ct, *sigma_nt, *sigma_kt;
2336 } t_luzar;
2337
2338 static real compute_weighted_rates(int n, real t[], real ct[], real nt[],
2339                                    real kt[], real sigma_ct[], real sigma_nt[],
2340                                    real sigma_kt[], real *k, real *kp,
2341                                    real *sigma_k, real *sigma_kp,
2342                                    real fit_start)
2343 {
2344 #define NK 10
2345     int      i, j;
2346     t_luzar  tl;
2347     real     kkk = 0, kkp = 0, kk2 = 0, kp2 = 0, chi2;
2348
2349     *sigma_k  = 0;
2350     *sigma_kp = 0;
2351
2352     for (i = 0; (i < n); i++)
2353     {
2354         if (t[i] >= fit_start)
2355         {
2356             break;
2357         }
2358     }
2359     tl.n0       = i;
2360     tl.n1       = n;
2361     tl.nparams  = 2;
2362     tl.ndelta   = 1;
2363     tl.t        = t;
2364     tl.ct       = ct;
2365     tl.nt       = nt;
2366     tl.kt       = kt;
2367     tl.sigma_ct = sigma_ct;
2368     tl.sigma_nt = sigma_nt;
2369     tl.sigma_kt = sigma_kt;
2370     tl.kkk[0]   = *k;
2371     tl.kkk[1]   = *kp;
2372
2373     chi2      = 0; /*optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */
2374     *k        = tl.kkk[0];
2375     *kp       = tl.kkk[1] = *kp;
2376     tl.ndelta = NK;
2377     for (j = 0; (j < NK); j++)
2378     {
2379         /* (void) optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */
2380         kkk += tl.kkk[0];
2381         kkp += tl.kkk[1];
2382         kk2 += sqr(tl.kkk[0]);
2383         kp2 += sqr(tl.kkk[1]);
2384         tl.n0++;
2385     }
2386     *sigma_k  = sqrt(kk2/NK - sqr(kkk/NK));
2387     *sigma_kp = sqrt(kp2/NK - sqr(kkp/NK));
2388
2389     return chi2;
2390 }
2391
2392 static void smooth_tail(int n, real t[], real c[], real sigma_c[], real start,
2393                         const output_env_t oenv)
2394 {
2395     FILE  *fp;
2396     real   e_1;
2397     double fitparm[4];
2398     int    i;
2399
2400     e_1 = exp(-1);
2401     for (i = 0; (i < n); i++)
2402     {
2403         if (c[i] < e_1)
2404         {
2405             break;
2406         }
2407     }
2408     if (i < n)
2409     {
2410         fitparm[0] = t[i];
2411     }
2412     else
2413     {
2414         fitparm[0] = 10;
2415     }
2416     fitparm[1] = 0.95;
2417     do_lmfit(n, c, sigma_c, 0, t, start, t[n-1], oenv, bDebugMode(),
2418              effnEXP2, fitparm, 0);
2419 }
2420
2421 void analyse_corr(int n, real t[], real ct[], real nt[], real kt[],
2422                   real sigma_ct[], real sigma_nt[], real sigma_kt[],
2423                   real fit_start, real temp, real smooth_tail_start,
2424                   const output_env_t oenv)
2425 {
2426     int        i0, i;
2427     real       k = 1, kp = 1, kow = 1;
2428     real       Q = 0, chi22, chi2, dg, dgp, tau_hb, dtau, tau_rlx, e_1, dt, sigma_k, sigma_kp, ddg;
2429     double     tmp, sn2 = 0, sc2 = 0, sk2 = 0, scn = 0, sck = 0, snk = 0;
2430     gmx_bool   bError = (sigma_ct != NULL) && (sigma_nt != NULL) && (sigma_kt != NULL);
2431
2432     if (smooth_tail_start >= 0)
2433     {
2434         smooth_tail(n, t, ct, sigma_ct, smooth_tail_start, oenv);
2435         smooth_tail(n, t, nt, sigma_nt, smooth_tail_start, oenv);
2436         smooth_tail(n, t, kt, sigma_kt, smooth_tail_start, oenv);
2437     }
2438     for (i0 = 0; (i0 < n-2) && ((t[i0]-t[0]) < fit_start); i0++)
2439     {
2440         ;
2441     }
2442     if (i0 < n-2)
2443     {
2444         for (i = i0; (i < n); i++)
2445         {
2446             sc2 += sqr(ct[i]);
2447             sn2 += sqr(nt[i]);
2448             sk2 += sqr(kt[i]);
2449             sck += ct[i]*kt[i];
2450             snk += nt[i]*kt[i];
2451             scn += ct[i]*nt[i];
2452         }
2453         printf("Hydrogen bond thermodynamics at T = %g K\n", temp);
2454         tmp = (sn2*sc2-sqr(scn));
2455         if ((tmp > 0) && (sn2 > 0))
2456         {
2457             k    = (sn2*sck-scn*snk)/tmp;
2458             kp   = (k*scn-snk)/sn2;
2459             if (bError)
2460             {
2461                 chi2 = 0;
2462                 for (i = i0; (i < n); i++)
2463                 {
2464                     chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
2465                 }
2466                 chi22 = compute_weighted_rates(n, t, ct, nt, kt, sigma_ct, sigma_nt,
2467                                                sigma_kt, &k, &kp,
2468                                                &sigma_k, &sigma_kp, fit_start);
2469                 Q   = 0; /* quality_of_fit(chi2, 2);*/
2470                 ddg = BOLTZ*temp*sigma_k/k;
2471                 printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n",
2472                        chi2, Q);
2473                 printf("The Rate and Delta G are followed by an error estimate\n");
2474                 printf("----------------------------------------------------------\n"
2475                        "Type      Rate (1/ps)  Sigma Time (ps)  DG (kJ/mol)  Sigma\n");
2476                 printf("Forward    %10.3f %6.2f   %8.3f  %10.3f %6.2f\n",
2477                        k, sigma_k, 1/k, calc_dg(1/k, temp), ddg);
2478                 ddg = BOLTZ*temp*sigma_kp/kp;
2479                 printf("Backward   %10.3f %6.2f   %8.3f  %10.3f %6.2f\n",
2480                        kp, sigma_kp, 1/kp, calc_dg(1/kp, temp), ddg);
2481             }
2482             else
2483             {
2484                 chi2 = 0;
2485                 for (i = i0; (i < n); i++)
2486                 {
2487                     chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
2488                 }
2489                 printf("Fitting parameters chi^2 = %10g\nQ = %10g\n",
2490                        chi2, Q);
2491                 printf("--------------------------------------------------\n"
2492                        "Type      Rate (1/ps) Time (ps)  DG (kJ/mol)  Chi^2\n");
2493                 printf("Forward    %10.3f   %8.3f  %10.3f  %10g\n",
2494                        k, 1/k, calc_dg(1/k, temp), chi2);
2495                 printf("Backward   %10.3f   %8.3f  %10.3f\n",
2496                        kp, 1/kp, calc_dg(1/kp, temp));
2497             }
2498         }
2499         if (sc2 > 0)
2500         {
2501             kow  = 2*sck/sc2;
2502             printf("One-way    %10.3f   %s%8.3f  %10.3f\n",
2503                    kow, bError ? "       " : "", 1/kow, calc_dg(1/kow, temp));
2504         }
2505         else
2506         {
2507             printf(" - Numerical problems computing HB thermodynamics:\n"
2508                    "sc2 = %g  sn2 = %g  sk2 = %g sck = %g snk = %g scn = %g\n",
2509                    sc2, sn2, sk2, sck, snk, scn);
2510         }
2511         /* Determine integral of the correlation function */
2512         tau_hb = evaluate_integral(n, t, ct, NULL, (t[n-1]-t[0])/2, &dtau);
2513         printf("Integral   %10.3f   %s%8.3f  %10.3f\n", 1/tau_hb,
2514                bError ? "       " : "", tau_hb, calc_dg(tau_hb, temp));
2515         e_1 = exp(-1);
2516         for (i = 0; (i < n-2); i++)
2517         {
2518             if ((ct[i] > e_1) && (ct[i+1] <= e_1))
2519             {
2520                 break;
2521             }
2522         }
2523         if (i < n-2)
2524         {
2525             /* Determine tau_relax from linear interpolation */
2526             tau_rlx = t[i]-t[0] + (e_1-ct[i])*(t[i+1]-t[i])/(ct[i+1]-ct[i]);
2527             printf("Relaxation %10.3f   %8.3f  %s%10.3f\n", 1/tau_rlx,
2528                    tau_rlx, bError ? "       " : "",
2529                    calc_dg(tau_rlx, temp));
2530         }
2531     }
2532     else
2533     {
2534         printf("Correlation functions too short to compute thermodynamics\n");
2535     }
2536 }
2537
2538 void compute_derivative(int nn, real x[], real y[], real dydx[])
2539 {
2540     int j;
2541
2542     /* Compute k(t) = dc(t)/dt */
2543     for (j = 1; (j < nn-1); j++)
2544     {
2545         dydx[j] = (y[j+1]-y[j-1])/(x[j+1]-x[j-1]);
2546     }
2547     /* Extrapolate endpoints */
2548     dydx[0]    = 2*dydx[1]   -  dydx[2];
2549     dydx[nn-1] = 2*dydx[nn-2] - dydx[nn-3];
2550 }
2551
2552 static void parallel_print(int *data, int nThreads)
2553 {
2554     /* This prints the donors on which each tread is currently working. */
2555     int i;
2556
2557     fprintf(stderr, "\r");
2558     for (i = 0; i < nThreads; i++)
2559     {
2560         fprintf(stderr, "%-7i", data[i]);
2561     }
2562 }
2563
2564 static void normalizeACF(real *ct, real *gt, int nhb, int len)
2565 {
2566     real ct_fac, gt_fac = 0;
2567     int  i;
2568
2569     /* Xu and Berne use the same normalization constant */
2570
2571     ct_fac = 1.0/ct[0];
2572     if (nhb != 0)
2573     {
2574         gt_fac = 1.0/(real)nhb;
2575     }
2576
2577     printf("Normalization for c(t) = %g for gh(t) = %g\n", ct_fac, gt_fac);
2578     for (i = 0; i < len; i++)
2579     {
2580         ct[i] *= ct_fac;
2581         if (gt != NULL)
2582         {
2583             gt[i] *= gt_fac;
2584         }
2585     }
2586 }
2587
2588 /* Added argument bContact in do_hbac(...). Also
2589  * added support for -contact in the actual code.
2590  * - Erik Marklund, May 31, 2006 */
2591 /* Changed contact code and added argument R2
2592  * - Erik Marklund, June 29, 2006
2593  */
2594 static void do_hbac(const char *fn, t_hbdata *hb,
2595                     int nDump, gmx_bool bMerge, gmx_bool bContact, real fit_start,
2596                     real temp, gmx_bool R2, real smooth_tail_start, const output_env_t oenv,
2597                     const char *gemType, int nThreads,
2598                     const int NN, const gmx_bool bBallistic, const gmx_bool bGemFit)
2599 {
2600     FILE          *fp;
2601     int            i, j, k, m, n, o, nd, ihb, idist, n2, nn, iter, nSets;
2602     const char    *legNN[]   = {
2603         "Ac(t)",
2604         "Ac'(t)"
2605     };
2606     static char  **legGem;
2607
2608     const char    *legLuzar[] = {
2609         "Ac\\sfin sys\\v{}\\z{}(t)",
2610         "Ac(t)",
2611         "Cc\\scontact,hb\\v{}\\z{}(t)",
2612         "-dAc\\sfs\\v{}\\z{}/dt"
2613     };
2614     gmx_bool       bNorm = FALSE, bOMP = FALSE;
2615     double         nhb   = 0;
2616     int            nhbi  = 0;
2617     real          *rhbex = NULL, *ht, *gt, *ght, *dght, *kt;
2618     real          *ct, *p_ct, tail, tail2, dtail, ct_fac, ght_fac, *cct;
2619     const real     tol     = 1e-3;
2620     int            nframes = hb->nframes, nf;
2621     unsigned int **h       = NULL, **g = NULL;
2622     int            nh, nhbonds, nhydro, ngh;
2623     t_hbond       *hbh;
2624     PSTYPE         p, *pfound = NULL, np;
2625     t_pShift      *pHist;
2626     int           *ptimes   = NULL, *poff = NULL, anhb, n0, mMax = INT_MIN;
2627     real         **rHbExGem = NULL;
2628     gmx_bool       c;
2629     int            acType;
2630     t_E           *E;
2631     double        *ctdouble, *timedouble, *fittedct;
2632     double         fittolerance = 0.1;
2633     int           *dondata      = NULL, thisThread;
2634
2635     enum {
2636         AC_NONE, AC_NN, AC_GEM, AC_LUZAR
2637     };
2638
2639 #ifdef GMX_OPENMP
2640     bOMP = TRUE;
2641 #else
2642     bOMP = FALSE;
2643 #endif
2644
2645     printf("Doing autocorrelation ");
2646
2647     /* Decide what kind of ACF calculations to do. */
2648     if (NN > NN_NONE && NN < NN_NR)
2649     {
2650 #ifdef HAVE_NN_LOOPS
2651         acType = AC_NN;
2652         printf("using the energy estimate.\n");
2653 #else
2654         acType = AC_NONE;
2655         printf("Can't do the NN-loop. Yet.\n");
2656 #endif
2657     }
2658     else if (hb->bGem)
2659     {
2660         acType = AC_GEM;
2661         printf("according to the reversible geminate recombination model by Omer Markowitch.\n");
2662
2663         nSets = 1 + (bBallistic ? 1 : 0) + (bGemFit ? 1 : 0);
2664         snew(legGem, nSets);
2665         for (i = 0; i < nSets; i++)
2666         {
2667             snew(legGem[i], 128);
2668         }
2669         sprintf(legGem[0], "Ac\\s%s\\v{}\\z{}(t)", gemType);
2670         if (bBallistic)
2671         {
2672             sprintf(legGem[1], "Ac'(t)");
2673         }
2674         if (bGemFit)
2675         {
2676             sprintf(legGem[(bBallistic ? 3 : 2)], "Ac\\s%s,fit\\v{}\\z{}(t)", gemType);
2677         }
2678
2679     }
2680     else
2681     {
2682         acType = AC_LUZAR;
2683         printf("according to the theory of Luzar and Chandler.\n");
2684     }
2685     fflush(stdout);
2686
2687     /* build hbexist matrix in reals for autocorr */
2688     /* Allocate memory for computing ACF (rhbex) and aggregating the ACF (ct) */
2689     n2 = 1;
2690     while (n2 < nframes)
2691     {
2692         n2 *= 2;
2693     }
2694
2695     nn = nframes/2;
2696
2697     if (acType != AC_NN || bOMP)
2698     {
2699         snew(h, hb->maxhydro);
2700         snew(g, hb->maxhydro);
2701     }
2702
2703     /* Dump hbonds for debugging */
2704     dump_ac(hb, bMerge || bContact, nDump);
2705
2706     /* Total number of hbonds analyzed here */
2707     nhbonds = 0;
2708     ngh     = 0;
2709     anhb    = 0;
2710
2711     if (acType != AC_LUZAR && bOMP)
2712     {
2713         nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
2714
2715         gmx_omp_set_num_threads(nThreads);
2716         snew(dondata, nThreads);
2717         for (i = 0; i < nThreads; i++)
2718         {
2719             dondata[i] = -1;
2720         }
2721         printf("ACF calculations parallelized with OpenMP using %i threads.\n"
2722                "Expect close to linear scaling over this donor-loop.\n", nThreads);
2723         fflush(stdout);
2724         fprintf(stderr, "Donors: [thread no]\n");
2725         {
2726             char tmpstr[7];
2727             for (i = 0; i < nThreads; i++)
2728             {
2729                 snprintf(tmpstr, 7, "[%i]", i);
2730                 fprintf(stderr, "%-7s", tmpstr);
2731             }
2732         }
2733         fprintf(stderr, "\n");
2734     }
2735
2736
2737     /* Build the ACF according to acType */
2738     switch (acType)
2739     {
2740
2741         case AC_NN:
2742 #ifdef HAVE_NN_LOOPS
2743             /* Here we're using the estimated energy for the hydrogen bonds. */
2744             snew(ct, nn);
2745
2746 #pragma omp parallel \
2747             private(i, j, k, nh, E, rhbex, thisThread) \
2748             default(shared)
2749             {
2750 #pragma omp barrier
2751                 thisThread = gmx_omp_get_thread_num();
2752                 rhbex      = NULL;
2753
2754                 snew(rhbex, n2);
2755                 memset(rhbex, 0, n2*sizeof(real)); /* Trust no-one, not even malloc()! */
2756
2757 #pragma omp barrier
2758 #pragma omp for schedule (dynamic)
2759                 for (i = 0; i < hb->d.nrd; i++) /* loop over donors */
2760                 {
2761                     if (bOMP)
2762                     {
2763 #pragma omp critical
2764                         {
2765                             dondata[thisThread] = i;
2766                             parallel_print(dondata, nThreads);
2767                         }
2768                     }
2769                     else
2770                     {
2771                         fprintf(stderr, "\r %i", i);
2772                     }
2773
2774                     for (j = 0; j < hb->a.nra; j++)              /* loop over acceptors */
2775                     {
2776                         for (nh = 0; nh < hb->d.nhydro[i]; nh++) /* loop over donors' hydrogens */
2777                         {
2778                             E = hb->hbE.E[i][j][nh];
2779                             if (E != NULL)
2780                             {
2781                                 for (k = 0; k < nframes; k++)
2782                                 {
2783                                     if (E[k] != NONSENSE_E)
2784                                     {
2785                                         rhbex[k] = (real)E[k];
2786                                     }
2787                                 }
2788
2789                                 low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &(rhbex), hb->time[1]-hb->time[0],
2790                                                 eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0, 1);
2791 #pragma omp critical
2792                                 {
2793                                     for (k = 0; (k < nn); k++)
2794                                     {
2795                                         ct[k] += rhbex[k];
2796                                     }
2797                                 }
2798                             }
2799                         } /* k loop */
2800                     }     /* j loop */
2801                 }         /* i loop */
2802                 sfree(rhbex);
2803 #pragma omp barrier
2804             }
2805
2806             if (bOMP)
2807             {
2808                 sfree(dondata);
2809             }
2810             normalizeACF(ct, NULL, 0, nn);
2811             snew(ctdouble, nn);
2812             snew(timedouble, nn);
2813             for (j = 0; j < nn; j++)
2814             {
2815                 timedouble[j] = (double)(hb->time[j]);
2816                 ctdouble[j]   = (double)(ct[j]);
2817             }
2818
2819             /* Remove ballistic term */
2820             /* Ballistic component removal and fitting to the reversible geminate recombination model
2821              * will be taken out for the time being. First of all, one can remove the ballistic
2822              * component with g_analyze afterwards. Secondly, and more importantly, there are still
2823              * problems with the robustness of the fitting to the model. More work is needed.
2824              * A third reason is that we're currently using gsl for this and wish to reduce dependence
2825              * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with
2826              * a BSD-licence that can do the job.
2827              *
2828              * - Erik Marklund, June 18 2010.
2829              */
2830 /*         if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */
2831 /*             takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */
2832 /*         else */
2833 /*             printf("\nNumber of data points is less than the number of parameters to fit\n." */
2834 /*                    "The system is underdetermined, hence no ballistic term can be found.\n\n"); */
2835
2836             fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)");
2837             xvgr_legend(fp, asize(legNN), legNN);
2838
2839             for (j = 0; (j < nn); j++)
2840             {
2841                 fprintf(fp, "%10g  %10g %10g\n",
2842                         hb->time[j]-hb->time[0],
2843                         ct[j],
2844                         ctdouble[j]);
2845             }
2846             xvgrclose(fp);
2847             sfree(ct);
2848             sfree(ctdouble);
2849             sfree(timedouble);
2850 #endif             /* HAVE_NN_LOOPS */
2851             break; /* case AC_NN */
2852
2853         case AC_GEM:
2854             snew(ct, 2*n2);
2855             memset(ct, 0, 2*n2*sizeof(real));
2856 #ifndef GMX_OPENMP
2857             fprintf(stderr, "Donor:\n");
2858 #define __ACDATA ct
2859 #else
2860 #define __ACDATA p_ct
2861 #endif
2862
2863 #pragma omp parallel \
2864             private(i, k, nh, hbh, pHist, h, g, n0, nf, np, j, m, \
2865             pfound, poff, rHbExGem, p, ihb, mMax, \
2866             thisThread, p_ct) \
2867             default(shared)
2868             { /* ##########  THE START OF THE ENORMOUS PARALLELIZED BLOCK!  ########## */
2869                 h          = NULL;
2870                 g          = NULL;
2871                 thisThread = gmx_omp_get_thread_num();
2872                 snew(h, hb->maxhydro);
2873                 snew(g, hb->maxhydro);
2874                 mMax     = INT_MIN;
2875                 rHbExGem = NULL;
2876                 poff     = NULL;
2877                 pfound   = NULL;
2878                 p_ct     = NULL;
2879                 snew(p_ct, 2*n2);
2880                 memset(p_ct, 0, 2*n2*sizeof(real));
2881
2882                 /* I'm using a chunk size of 1, since I expect      \
2883                  * the overhead to be really small compared         \
2884                  * to the actual calculations                       \ */
2885 #pragma omp for schedule(dynamic,1) nowait
2886                 for (i = 0; i < hb->d.nrd; i++)
2887                 {
2888
2889                     if (bOMP)
2890                     {
2891 #pragma omp critical
2892                         {
2893                             dondata[thisThread] = i;
2894                             parallel_print(dondata, nThreads);
2895                         }
2896                     }
2897                     else
2898                     {
2899                         fprintf(stderr, "\r %i", i);
2900                     }
2901                     for (k = 0; k < hb->a.nra; k++)
2902                     {
2903                         for (nh = 0; nh < ((bMerge || bContact) ? 1 : hb->d.nhydro[i]); nh++)
2904                         {
2905                             hbh = hb->hbmap[i][k];
2906                             if (hbh)
2907                             {
2908                                 /* Note that if hb->per->gemtype==gemDD, then distances will be stored in
2909                                  * hb->hbmap[d][a].h array anyway, because the contact flag will be set.
2910                                  * hence, it's only with the gemAD mode that hb->hbmap[d][a].g will be used. */
2911                                 pHist = &(hb->per->pHist[i][k]);
2912                                 if (ISHB(hbh->history[nh]) && pHist->len != 0)
2913                                 {
2914
2915                                     {
2916                                         h[nh] = hbh->h[nh];
2917                                         g[nh] = hb->per->gemtype == gemAD ? hbh->g[nh] : NULL;
2918                                     }
2919                                     n0 = hbh->n0;
2920                                     nf = hbh->nframes;
2921                                     /* count the number of periodic shifts encountered and store
2922                                      * them in separate arrays. */
2923                                     np = 0;
2924                                     for (j = 0; j < pHist->len; j++)
2925                                     {
2926                                         p = pHist->p[j];
2927                                         for (m = 0; m <= np; m++)
2928                                         {
2929                                             if (m == np) /* p not recognized in list. Add it and set up new array. */
2930                                             {
2931                                                 np++;
2932                                                 if (np > hb->per->nper)
2933                                                 {
2934                                                     gmx_fatal(FARGS, "Too many pshifts. Something's utterly wrong here.");
2935                                                 }
2936                                                 if (m >= mMax) /* Extend the arrays.
2937                                                                 * Doing it like this, using mMax to keep track of the sizes,
2938                                                                 * eleviates the need for freeing and re-allocating the arrays
2939                                                                 * when taking on the next donor-acceptor pair */
2940                                                 {
2941                                                     mMax = m;
2942                                                     srenew(pfound, np);   /* The list of found periodic shifts. */
2943                                                     srenew(rHbExGem, np); /* The hb existence functions (-aver_hb). */
2944                                                     snew(rHbExGem[m], 2*n2);
2945                                                     srenew(poff, np);
2946                                                 }
2947
2948                                                 {
2949                                                     if (rHbExGem != NULL && rHbExGem[m] != NULL)
2950                                                     {
2951                                                         /* This must be done, as this array was most likey
2952                                                          * used to store stuff in some previous iteration. */
2953                                                         memset(rHbExGem[m], 0, (sizeof(real)) * (2*n2));
2954                                                     }
2955                                                     else
2956                                                     {
2957                                                         fprintf(stderr, "rHbExGem not initialized! m = %i\n", m);
2958                                                     }
2959                                                 }
2960                                                 pfound[m] = p;
2961                                                 poff[m]   = -1;
2962
2963                                                 break;
2964                                             } /* m==np */
2965                                             if (p == pfound[m])
2966                                             {
2967                                                 break;
2968                                             }
2969                                         } /* m: Loop over found shifts */
2970                                     }     /* j: Loop over shifts */
2971
2972                                     /* Now unpack and disentangle the existence funtions. */
2973                                     for (j = 0; j < nf; j++)
2974                                     {
2975                                         /* i:       donor,
2976                                          * k:       acceptor
2977                                          * nh:      hydrogen
2978                                          * j:       time
2979                                          * p:       periodic shift
2980                                          * pfound:  list of periodic shifts found for this pair.
2981                                          * poff:    list of frame offsets; that is, the first
2982                                          *          frame a hbond has a particular periodic shift. */
2983                                         p = getPshift(*pHist, j+n0);
2984                                         if (p != -1)
2985                                         {
2986                                             for (m = 0; m < np; m++)
2987                                             {
2988                                                 if (pfound[m] == p)
2989                                                 {
2990                                                     break;
2991                                                 }
2992                                                 if (m == (np-1))
2993                                                 {
2994                                                     gmx_fatal(FARGS, "Shift not found, but must be there.");
2995                                                 }
2996                                             }
2997
2998                                             ihb = is_hb(h[nh], j) || ((hb->per->gemtype != gemAD || j == 0) ? FALSE : is_hb(g[nh], j));
2999                                             if (ihb)
3000                                             {
3001                                                 if (poff[m] == -1)
3002                                                 {
3003                                                     poff[m] = j; /* Here's where the first hbond with shift p is,
3004                                                                   * relative to the start of h[0].*/
3005                                                 }
3006                                                 if (j < poff[m])
3007                                                 {
3008                                                     gmx_fatal(FARGS, "j<poff[m]");
3009                                                 }
3010                                                 rHbExGem[m][j-poff[m]] += 1;
3011                                             }
3012                                         }
3013                                     }
3014
3015                                     /* Now, build ac. */
3016                                     for (m = 0; m < np; m++)
3017                                     {
3018                                         if (rHbExGem[m][0] > 0  && n0+poff[m] < nn /*  && m==0 */)
3019                                         {
3020                                             low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &(rHbExGem[m]), hb->time[1]-hb->time[0],
3021                                                             eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0);
3022                                             for (j = 0; (j < nn); j++)
3023                                             {
3024                                                 __ACDATA[j] += rHbExGem[m][j];
3025                                             }
3026                                         }
3027                                     } /* Building of ac. */
3028                                 }     /* if (ISHB(...*/
3029                             }         /* if (hbh) */
3030                         }             /* hydrogen loop */
3031                     }                 /* acceptor loop */
3032                 }                     /* donor loop */
3033
3034                 for (m = 0; m <= mMax; m++)
3035                 {
3036                     sfree(rHbExGem[m]);
3037                 }
3038                 sfree(pfound);
3039                 sfree(poff);
3040                 sfree(rHbExGem);
3041
3042                 sfree(h);
3043                 sfree(g);
3044
3045                 if (bOMP)
3046                 {
3047 #pragma omp critical
3048                     {
3049                         for (i = 0; i < nn; i++)
3050                         {
3051                             ct[i] += p_ct[i];
3052                         }
3053                     }
3054                     sfree(p_ct);
3055                 }
3056
3057             } /* ########## THE END OF THE ENORMOUS PARALLELIZED BLOCK ########## */
3058             if (bOMP)
3059             {
3060                 sfree(dondata);
3061             }
3062
3063             normalizeACF(ct, NULL, 0, nn);
3064
3065             fprintf(stderr, "\n\nACF successfully calculated.\n");
3066
3067             /* Use this part to fit to geminate recombination - JCP 129, 84505 (2008) */
3068
3069             snew(ctdouble, nn);
3070             snew(timedouble, nn);
3071             snew(fittedct, nn);
3072
3073             for (j = 0; j < nn; j++)
3074             {
3075                 timedouble[j] = (double)(hb->time[j]);
3076                 ctdouble[j]   = (double)(ct[j]);
3077             }
3078
3079             /* Remove ballistic term */
3080             /* Ballistic component removal and fitting to the reversible geminate recombination model
3081              * will be taken out for the time being. First of all, one can remove the ballistic
3082              * component with g_analyze afterwards. Secondly, and more importantly, there are still
3083              * problems with the robustness of the fitting to the model. More work is needed.
3084              * A third reason is that we're currently using gsl for this and wish to reduce dependence
3085              * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with
3086              * a BSD-licence that can do the job.
3087              *
3088              * - Erik Marklund, June 18 2010.
3089              */
3090 /*         if (bBallistic) { */
3091 /*             if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */
3092 /*                 takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */
3093 /*             else */
3094 /*                 printf("\nNumber of data points is less than the number of parameters to fit\n." */
3095 /*                        "The system is underdetermined, hence no ballistic term can be found.\n\n"); */
3096 /*         } */
3097 /*         if (bGemFit) */
3098 /*             fitGemRecomb(ctdouble, timedouble, &fittedct, nn, params); */
3099
3100
3101             if (bContact)
3102             {
3103                 fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3104             }
3105             else
3106             {
3107                 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3108             }
3109             xvgr_legend(fp, asize(legGem), (const char**)legGem, oenv);
3110
3111             for (j = 0; (j < nn); j++)
3112             {
3113                 fprintf(fp, "%10g  %10g", hb->time[j]-hb->time[0], ct[j]);
3114                 if (bBallistic)
3115                 {
3116                     fprintf(fp, "  %10g", ctdouble[j]);
3117                 }
3118                 if (bGemFit)
3119                 {
3120                     fprintf(fp, "  %10g", fittedct[j]);
3121                 }
3122                 fprintf(fp, "\n");
3123             }
3124             xvgrclose(fp);
3125
3126             sfree(ctdouble);
3127             sfree(timedouble);
3128             sfree(fittedct);
3129             sfree(ct);
3130
3131             break; /* case AC_GEM */
3132
3133         case AC_LUZAR:
3134             snew(rhbex, 2*n2);
3135             snew(ct, 2*n2);
3136             snew(gt, 2*n2);
3137             snew(ht, 2*n2);
3138             snew(ght, 2*n2);
3139             snew(dght, 2*n2);
3140
3141             snew(kt, nn);
3142             snew(cct, nn);
3143
3144             for (i = 0; (i < hb->d.nrd); i++)
3145             {
3146                 for (k = 0; (k < hb->a.nra); k++)
3147                 {
3148                     nhydro = 0;
3149                     hbh    = hb->hbmap[i][k];
3150
3151                     if (hbh)
3152                     {
3153                         if (bMerge || bContact)
3154                         {
3155                             if (ISHB(hbh->history[0]))
3156                             {
3157                                 h[0]   = hbh->h[0];
3158                                 g[0]   = hbh->g[0];
3159                                 nhydro = 1;
3160                             }
3161                         }
3162                         else
3163                         {
3164                             for (m = 0; (m < hb->maxhydro); m++)
3165                             {
3166                                 if (bContact ? ISDIST(hbh->history[m]) : ISHB(hbh->history[m]))
3167                                 {
3168                                     g[nhydro] = hbh->g[m];
3169                                     h[nhydro] = hbh->h[m];
3170                                     nhydro++;
3171                                 }
3172                             }
3173                         }
3174
3175                         nf = hbh->nframes;
3176                         for (nh = 0; (nh < nhydro); nh++)
3177                         {
3178                             int nrint = bContact ? hb->nrdist : hb->nrhb;
3179                             if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint))
3180                             {
3181                                 fprintf(stderr, "\rACF %d/%d", nhbonds+1, nrint);
3182                             }
3183                             nhbonds++;
3184                             for (j = 0; (j < nframes); j++)
3185                             {
3186                                 /* Changed '<' into '<=' below, just like I did in
3187                                    the hbm-output-loop in the gmx_hbond() block.
3188                                    - Erik Marklund, May 31, 2006 */
3189                                 if (j <= nf)
3190                                 {
3191                                     ihb   = is_hb(h[nh], j);
3192                                     idist = is_hb(g[nh], j);
3193                                 }
3194                                 else
3195                                 {
3196                                     ihb = idist = 0;
3197                                 }
3198                                 rhbex[j] = ihb;
3199                                 /* For contacts: if a second cut-off is provided, use it,
3200                                  * otherwise use g(t) = 1-h(t) */
3201                                 if (!R2 && bContact)
3202                                 {
3203                                     gt[j]  = 1-ihb;
3204                                 }
3205                                 else
3206                                 {
3207                                     gt[j]  = idist*(1-ihb);
3208                                 }
3209                                 ht[j]    = rhbex[j];
3210                                 nhb     += ihb;
3211                             }
3212
3213
3214                             /* The autocorrelation function is normalized after summation only */
3215                             low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &rhbex, hb->time[1]-hb->time[0],
3216                                             eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0);
3217
3218                             /* Cross correlation analysis for thermodynamics */
3219                             for (j = nframes; (j < n2); j++)
3220                             {
3221                                 ht[j] = 0;
3222                                 gt[j] = 0;
3223                             }
3224
3225                             cross_corr(n2, ht, gt, dght);
3226
3227                             for (j = 0; (j < nn); j++)
3228                             {
3229                                 ct[j]  += rhbex[j];
3230                                 ght[j] += dght[j];
3231                             }
3232                         }
3233                     }
3234                 }
3235             }
3236             fprintf(stderr, "\n");
3237             sfree(h);
3238             sfree(g);
3239             normalizeACF(ct, ght, nhb, nn);
3240
3241             /* Determine tail value for statistics */
3242             tail  = 0;
3243             tail2 = 0;
3244             for (j = nn/2; (j < nn); j++)
3245             {
3246                 tail  += ct[j];
3247                 tail2 += ct[j]*ct[j];
3248             }
3249             tail  /= (nn - nn/2);
3250             tail2 /= (nn - nn/2);
3251             dtail  = sqrt(tail2-tail*tail);
3252
3253             /* Check whether the ACF is long enough */
3254             if (dtail > tol)
3255             {
3256                 printf("\nWARNING: Correlation function is probably not long enough\n"
3257                        "because the standard deviation in the tail of C(t) > %g\n"
3258                        "Tail value (average C(t) over second half of acf): %g +/- %g\n",
3259                        tol, tail, dtail);
3260             }
3261             for (j = 0; (j < nn); j++)
3262             {
3263                 cct[j] = ct[j];
3264                 ct[j]  = (cct[j]-tail)/(1-tail);
3265             }
3266             /* Compute negative derivative k(t) = -dc(t)/dt */
3267             compute_derivative(nn, hb->time, ct, kt);
3268             for (j = 0; (j < nn); j++)
3269             {
3270                 kt[j] = -kt[j];
3271             }
3272
3273
3274             if (bContact)
3275             {
3276                 fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3277             }
3278             else
3279             {
3280                 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3281             }
3282             xvgr_legend(fp, asize(legLuzar), legLuzar, oenv);
3283
3284
3285             for (j = 0; (j < nn); j++)
3286             {
3287                 fprintf(fp, "%10g  %10g  %10g  %10g  %10g\n",
3288                         hb->time[j]-hb->time[0], ct[j], cct[j], ght[j], kt[j]);
3289             }
3290             xvgrclose(fp);
3291
3292             analyse_corr(nn, hb->time, ct, ght, kt, NULL, NULL, NULL,
3293                          fit_start, temp, smooth_tail_start, oenv);
3294
3295             do_view(oenv, fn, NULL);
3296             sfree(rhbex);
3297             sfree(ct);
3298             sfree(gt);
3299             sfree(ht);
3300             sfree(ght);
3301             sfree(dght);
3302             sfree(cct);
3303             sfree(kt);
3304             /* sfree(h); */
3305 /*         sfree(g); */
3306
3307             break; /* case AC_LUZAR */
3308
3309         default:
3310             gmx_fatal(FARGS, "Unrecognized type of ACF-calulation. acType = %i.", acType);
3311     } /* switch (acType) */
3312 }
3313
3314 static void init_hbframe(t_hbdata *hb, int nframes, real t)
3315 {
3316     int i, j, m;
3317
3318     hb->time[nframes]   = t;
3319     hb->nhb[nframes]    = 0;
3320     hb->ndist[nframes]  = 0;
3321     for (i = 0; (i < max_hx); i++)
3322     {
3323         hb->nhx[nframes][i] = 0;
3324     }
3325     /* Loop invalidated */
3326     if (hb->bHBmap && 0)
3327     {
3328         for (i = 0; (i < hb->d.nrd); i++)
3329         {
3330             for (j = 0; (j < hb->a.nra); j++)
3331             {
3332                 for (m = 0; (m < hb->maxhydro); m++)
3333                 {
3334                     if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[m])
3335                     {
3336                         set_hb(hb, i, m, j, nframes, HB_NO);
3337                     }
3338                 }
3339             }
3340         }
3341     }
3342     /*set_hb(hb->hbmap[i][j]->h[m],nframes-hb->hbmap[i][j]->n0,HB_NO);*/
3343 }
3344
3345 static FILE *open_donor_properties_file(const char        *fn,
3346                                         t_hbdata          *hb,
3347                                         const output_env_t oenv)
3348 {
3349     FILE       *fp    = NULL;
3350     const char *leg[] = { "Nbound", "Nfree" };
3351
3352     if (!fn || !hb)
3353     {
3354         return NULL;
3355     }
3356
3357     fp = xvgropen(fn, "Donor properties", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
3358     xvgr_legend(fp, asize(leg), leg, oenv);
3359
3360     return fp;
3361 }
3362
3363 static void analyse_donor_properties(FILE *fp, t_hbdata *hb, int nframes, real t)
3364 {
3365     int i, j, k, nbound, nb, nhtot;
3366
3367     if (!fp || !hb)
3368     {
3369         return;
3370     }
3371     nbound = 0;
3372     nhtot  = 0;
3373     for (i = 0; (i < hb->d.nrd); i++)
3374     {
3375         for (k = 0; (k < hb->d.nhydro[i]); k++)
3376         {
3377             nb = 0;
3378             nhtot++;
3379             for (j = 0; (j < hb->a.nra) && (nb == 0); j++)
3380             {
3381                 if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[k] &&
3382                     is_hb(hb->hbmap[i][j]->h[k], nframes))
3383                 {
3384                     nb = 1;
3385                 }
3386             }
3387             nbound += nb;
3388         }
3389     }
3390     fprintf(fp, "%10.3e  %6d  %6d\n", t, nbound, nhtot-nbound);
3391 }
3392
3393 static void dump_hbmap(t_hbdata *hb,
3394                        int nfile, t_filenm fnm[], gmx_bool bTwo,
3395                        gmx_bool bContact, int isize[], int *index[], char *grpnames[],
3396                        t_atoms *atoms)
3397 {
3398     FILE    *fp, *fplog;
3399     int      ddd, hhh, aaa, i, j, k, m, grp;
3400     char     ds[32], hs[32], as[32];
3401     gmx_bool first;
3402
3403     fp = opt2FILE("-hbn", nfile, fnm, "w");
3404     if (opt2bSet("-g", nfile, fnm))
3405     {
3406         fplog = gmx_ffopen(opt2fn("-g", nfile, fnm), "w");
3407         fprintf(fplog, "# %10s  %12s  %12s\n", "Donor", "Hydrogen", "Acceptor");
3408     }
3409     else
3410     {
3411         fplog = NULL;
3412     }
3413     for (grp = gr0; grp <= (bTwo ? gr1 : gr0); grp++)
3414     {
3415         fprintf(fp, "[ %s ]", grpnames[grp]);
3416         for (i = 0; i < isize[grp]; i++)
3417         {
3418             fprintf(fp, (i%15) ? " " : "\n");
3419             fprintf(fp, " %4d", index[grp][i]+1);
3420         }
3421         fprintf(fp, "\n");
3422         /*
3423            Added -contact support below.
3424            - Erik Marklund, May 29, 2006
3425          */
3426         if (!bContact)
3427         {
3428             fprintf(fp, "[ donors_hydrogens_%s ]\n", grpnames[grp]);
3429             for (i = 0; (i < hb->d.nrd); i++)
3430             {
3431                 if (hb->d.grp[i] == grp)
3432                 {
3433                     for (j = 0; (j < hb->d.nhydro[i]); j++)
3434                     {
3435                         fprintf(fp, " %4d %4d", hb->d.don[i]+1,
3436                                 hb->d.hydro[i][j]+1);
3437                     }
3438                     fprintf(fp, "\n");
3439                 }
3440             }
3441             first = TRUE;
3442             fprintf(fp, "[ acceptors_%s ]", grpnames[grp]);
3443             for (i = 0; (i < hb->a.nra); i++)
3444             {
3445                 if (hb->a.grp[i] == grp)
3446                 {
3447                     fprintf(fp, (i%15 && !first) ? " " : "\n");
3448                     fprintf(fp, " %4d", hb->a.acc[i]+1);
3449                     first = FALSE;
3450                 }
3451             }
3452             fprintf(fp, "\n");
3453         }
3454     }
3455     if (bTwo)
3456     {
3457         fprintf(fp, bContact ? "[ contacts_%s-%s ]\n" :
3458                 "[ hbonds_%s-%s ]\n", grpnames[0], grpnames[1]);
3459     }
3460     else
3461     {
3462         fprintf(fp, bContact ? "[ contacts_%s ]" : "[ hbonds_%s ]\n", grpnames[0]);
3463     }
3464
3465     for (i = 0; (i < hb->d.nrd); i++)
3466     {
3467         ddd = hb->d.don[i];
3468         for (k = 0; (k < hb->a.nra); k++)
3469         {
3470             aaa = hb->a.acc[k];
3471             for (m = 0; (m < hb->d.nhydro[i]); m++)
3472             {
3473                 if (hb->hbmap[i][k] && ISHB(hb->hbmap[i][k]->history[m]))
3474                 {
3475                     sprintf(ds, "%s", mkatomname(atoms, ddd));
3476                     sprintf(as, "%s", mkatomname(atoms, aaa));
3477                     if (bContact)
3478                     {
3479                         fprintf(fp, " %6d %6d\n", ddd+1, aaa+1);
3480                         if (fplog)
3481                         {
3482                             fprintf(fplog, "%12s  %12s\n", ds, as);
3483                         }
3484                     }
3485                     else
3486                     {
3487                         hhh = hb->d.hydro[i][m];
3488                         sprintf(hs, "%s", mkatomname(atoms, hhh));
3489                         fprintf(fp, " %6d %6d %6d\n", ddd+1, hhh+1, aaa+1);
3490                         if (fplog)
3491                         {
3492                             fprintf(fplog, "%12s  %12s  %12s\n", ds, hs, as);
3493                         }
3494                     }
3495                 }
3496             }
3497         }
3498     }
3499     gmx_ffclose(fp);
3500     if (fplog)
3501     {
3502         gmx_ffclose(fplog);
3503     }
3504 }
3505
3506 /* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template.
3507  * It mimics add_frames() and init_frame() to some extent. */
3508 static void sync_hbdata(t_hbdata *p_hb, int nframes)
3509 {
3510     int i;
3511     if (nframes >= p_hb->max_frames)
3512     {
3513         p_hb->max_frames += 4096;
3514         srenew(p_hb->nhb,   p_hb->max_frames);
3515         srenew(p_hb->ndist, p_hb->max_frames);
3516         srenew(p_hb->n_bound, p_hb->max_frames);
3517         srenew(p_hb->nhx, p_hb->max_frames);
3518         if (p_hb->bDAnr)
3519         {
3520             srenew(p_hb->danr, p_hb->max_frames);
3521         }
3522         memset(&(p_hb->nhb[nframes]),   0, sizeof(int) * (p_hb->max_frames-nframes));
3523         memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
3524         p_hb->nhb[nframes]   = 0;
3525         p_hb->ndist[nframes] = 0;
3526
3527     }
3528     p_hb->nframes = nframes;
3529 /*     for (i=0;) */
3530 /*     { */
3531 /*         p_hb->nhx[nframes][i] */
3532 /*     } */
3533     memset(&(p_hb->nhx[nframes]), 0, sizeof(int)*max_hx); /* zero the helix count for this frame */
3534
3535     /* hb->per will remain constant througout the frame loop,
3536      * even though the data its members point to will change,
3537      * hence no need for re-syncing. */
3538 }
3539
3540 int gmx_hbond(int argc, char *argv[])
3541 {
3542     const char        *desc[] = {
3543         "[THISMODULE] computes and analyzes hydrogen bonds. Hydrogen bonds are",
3544         "determined based on cutoffs for the angle Hydrogen - Donor - Acceptor",
3545         "(zero is extended) and the distance Donor - Acceptor",
3546         "(or Hydrogen - Acceptor using [TT]-noda[tt]).",
3547         "OH and NH groups are regarded as donors, O is an acceptor always,",
3548         "N is an acceptor by default, but this can be switched using",
3549         "[TT]-nitacc[tt]. Dummy hydrogen atoms are assumed to be connected",
3550         "to the first preceding non-hydrogen atom.[PAR]",
3551
3552         "You need to specify two groups for analysis, which must be either",
3553         "identical or non-overlapping. All hydrogen bonds between the two",
3554         "groups are analyzed.[PAR]",
3555
3556         "If you set [TT]-shell[tt], you will be asked for an additional index group",
3557         "which should contain exactly one atom. In this case, only hydrogen",
3558         "bonds between atoms within the shell distance from the one atom are",
3559         "considered.[PAR]",
3560
3561         "With option -ac, rate constants for hydrogen bonding can be derived with the model of Luzar and Chandler",
3562         "(Nature 394, 1996; J. Chem. Phys. 113:23, 2000) or that of Markovitz and Agmon (J. Chem. Phys 129, 2008).",
3563         "If contact kinetics are analyzed by using the -contact option, then",
3564         "n(t) can be defined as either all pairs that are not within contact distance r at time t",
3565         "(corresponding to leaving the -r2 option at the default value 0) or all pairs that",
3566         "are within distance r2 (corresponding to setting a second cut-off value with option -r2).",
3567         "See mentioned literature for more details and definitions."
3568         "[PAR]",
3569
3570         /*    "It is also possible to analyse specific hydrogen bonds with",
3571               "[TT]-sel[tt]. This index file must contain a group of atom triplets",
3572               "Donor Hydrogen Acceptor, in the following way::",
3573            "",
3574            "[ selected ]",
3575            "     20    21    24",
3576            "     25    26    29",
3577            "      1     3     6",
3578            "",
3579            "Note that the triplets need not be on separate lines.",
3580            "Each atom triplet specifies a hydrogen bond to be analyzed,",
3581            "note also that no check is made for the types of atoms.[PAR]",
3582          */
3583
3584         "[BB]Output:[bb]",
3585         "",
3586         " * [TT]-num[tt]:  number of hydrogen bonds as a function of time.",
3587         " * [TT]-ac[tt]:   average over all autocorrelations of the existence",
3588         "   functions (either 0 or 1) of all hydrogen bonds.",
3589         " * [TT]-dist[tt]: distance distribution of all hydrogen bonds.",
3590         " * [TT]-ang[tt]:  angle distribution of all hydrogen bonds.",
3591         " * [TT]-hx[tt]:   the number of n-n+i hydrogen bonds as a function of time",
3592         "   where n and n+i stand for residue numbers and i ranges from 0 to 6.",
3593         "   This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated",
3594         "   with helices in proteins.",
3595         " * [TT]-hbn[tt]:  all selected groups, donors, hydrogens and acceptors",
3596         "   for selected groups, all hydrogen bonded atoms from all groups and",
3597         "   all solvent atoms involved in insertion.",
3598         " * [TT]-hbm[tt]:  existence matrix for all hydrogen bonds over all",
3599         "   frames, this also contains information on solvent insertion",
3600         "   into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]",
3601         "   index file.",
3602         " * [TT]-dan[tt]: write out the number of donors and acceptors analyzed for",
3603         "   each timeframe. This is especially useful when using [TT]-shell[tt].",
3604         " * [TT]-nhbdist[tt]: compute the number of HBonds per hydrogen in order to",
3605         "   compare results to Raman Spectroscopy.",
3606         "",
3607         "Note: options [TT]-ac[tt], [TT]-life[tt], [TT]-hbn[tt] and [TT]-hbm[tt]",
3608         "require an amount of memory proportional to the total numbers of donors",
3609         "times the total number of acceptors in the selected group(s)."
3610     };
3611
3612     static real        acut     = 30, abin = 1, rcut = 0.35, r2cut = 0, rbin = 0.005, rshell = -1;
3613     static real        maxnhb   = 0, fit_start = 1, fit_end = 60, temp = 298.15, smooth_tail_start = -1, D = -1;
3614     static gmx_bool    bNitAcc  = TRUE, bDA = TRUE, bMerge = TRUE;
3615     static int         nDump    = 0, nFitPoints = 100;
3616     static int         nThreads = 0, nBalExp = 4;
3617
3618     static gmx_bool    bContact     = FALSE, bBallistic = FALSE, bGemFit = FALSE;
3619     static real        logAfterTime = 10, gemBallistic = 0.2; /* ps */
3620     static const char *NNtype[]     = {NULL, "none", "binary", "oneOverR3", "dipole", NULL};
3621
3622     /* options */
3623     t_pargs     pa [] = {
3624         { "-a",    FALSE,  etREAL, {&acut},
3625           "Cutoff angle (degrees, Hydrogen - Donor - Acceptor)" },
3626         { "-r",    FALSE,  etREAL, {&rcut},
3627           "Cutoff radius (nm, X - Acceptor, see next option)" },
3628         { "-da",   FALSE,  etBOOL, {&bDA},
3629           "Use distance Donor-Acceptor (if TRUE) or Hydrogen-Acceptor (FALSE)" },
3630         { "-r2",   FALSE,  etREAL, {&r2cut},
3631           "Second cutoff radius. Mainly useful with [TT]-contact[tt] and [TT]-ac[tt]"},
3632         { "-abin", FALSE,  etREAL, {&abin},
3633           "Binwidth angle distribution (degrees)" },
3634         { "-rbin", FALSE,  etREAL, {&rbin},
3635           "Binwidth distance distribution (nm)" },
3636         { "-nitacc", FALSE, etBOOL, {&bNitAcc},
3637           "Regard nitrogen atoms as acceptors" },
3638         { "-contact", FALSE, etBOOL, {&bContact},
3639           "Do not look for hydrogen bonds, but merely for contacts within the cut-off distance" },
3640         { "-shell", FALSE, etREAL, {&rshell},
3641           "when > 0, only calculate hydrogen bonds within # nm shell around "
3642           "one particle" },
3643         { "-fitstart", FALSE, etREAL, {&fit_start},
3644           "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]" },
3645         { "-fitend", FALSE, etREAL, {&fit_end},
3646           "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])" },
3647         { "-temp",  FALSE, etREAL, {&temp},
3648           "Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" },
3649         { "-smooth", FALSE, etREAL, {&smooth_tail_start},
3650           "If >= 0, the tail of the ACF will be smoothed by fitting it to an exponential function: y = A exp(-x/[GRK]tau[grk])" },
3651         { "-dump",  FALSE, etINT, {&nDump},
3652           "Dump the first N hydrogen bond ACFs in a single [REF].xvg[ref] file for debugging" },
3653         { "-max_hb", FALSE, etREAL, {&maxnhb},
3654           "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" },
3655         { "-merge", FALSE, etBOOL, {&bMerge},
3656           "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." },
3657         { "-geminate", FALSE, etENUM, {gemType},
3658           "HIDDENUse reversible geminate recombination for the kinetics/thermodynamics calclations. See Markovitch et al., J. Chem. Phys 129, 084505 (2008) for details."},
3659         { "-diff", FALSE, etREAL, {&D},
3660           "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."},
3661 #ifdef GMX_OPENMP
3662         { "-nthreads", FALSE, etINT, {&nThreads},
3663           "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)"},
3664 #endif
3665     };
3666     const char *bugs[] = {
3667         "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and therefore not available for the time being."
3668     };
3669     t_filenm    fnm[] = {
3670         { efTRX, "-f",   NULL,     ffREAD  },
3671         { efTPR, NULL,   NULL,     ffREAD  },
3672         { efNDX, NULL,   NULL,     ffOPTRD },
3673         /*    { efNDX, "-sel", "select", ffOPTRD },*/
3674         { efXVG, "-num", "hbnum",  ffWRITE },
3675         { efLOG, "-g",   "hbond",  ffOPTWR },
3676         { efXVG, "-ac",  "hbac",   ffOPTWR },
3677         { efXVG, "-dist", "hbdist", ffOPTWR },
3678         { efXVG, "-ang", "hbang",  ffOPTWR },
3679         { efXVG, "-hx",  "hbhelix", ffOPTWR },
3680         { efNDX, "-hbn", "hbond",  ffOPTWR },
3681         { efXPM, "-hbm", "hbmap",  ffOPTWR },
3682         { efXVG, "-don", "donor",  ffOPTWR },
3683         { efXVG, "-dan", "danum",  ffOPTWR },
3684         { efXVG, "-life", "hblife", ffOPTWR },
3685         { efXVG, "-nhbdist", "nhbdist", ffOPTWR }
3686
3687     };
3688 #define NFILE asize(fnm)
3689
3690     char                  hbmap [HB_NR] = { ' ',    'o',      '-',       '*' };
3691     const char           *hbdesc[HB_NR] = { "None", "Present", "Inserted", "Present & Inserted" };
3692     t_rgb                 hbrgb [HB_NR] = { {1, 1, 1}, {1, 0, 0},   {0, 0, 1},    {1, 0, 1} };
3693
3694     t_trxstatus          *status;
3695     int                   trrStatus = 1;
3696     t_topology            top;
3697     t_inputrec            ir;
3698     t_pargs              *ppa;
3699     int                   npargs, natoms, nframes = 0, shatom;
3700     int                  *isize;
3701     char                **grpnames;
3702     atom_id             **index;
3703     rvec                 *x, hbox;
3704     matrix                box;
3705     real                  t, ccut, dist = 0.0, ang = 0.0;
3706     double                max_nhb, aver_nhb, aver_dist;
3707     int                   h = 0, i = 0, j, k = 0, l, start, end, id, ja, ogrp, nsel;
3708     int                   xi, yi, zi, ai;
3709     int                   xj, yj, zj, aj, xjj, yjj, zjj;
3710     int                   xk, yk, zk, ak, xkk, ykk, zkk;
3711     gmx_bool              bSelected, bHBmap, bStop, bTwo, was, bBox, bTric;
3712     int                  *adist, *rdist, *aptr, *rprt;
3713     int                   grp, nabin, nrbin, bin, resdist, ihb;
3714     char                **leg;
3715     t_hbdata             *hb, *hbptr;
3716     FILE                 *fp, *fpins = NULL, *fpnhb = NULL, *donor_properties = NULL;
3717     t_gridcell         ***grid;
3718     t_ncell              *icell, *jcell, *kcell;
3719     ivec                  ngrid;
3720     unsigned char        *datable;
3721     output_env_t          oenv;
3722     int                   gemmode, NN;
3723     PSTYPE                peri = 0;
3724     t_E                   E;
3725     int                   ii, jj, hh, actual_nThreads;
3726     int                   threadNr = 0;
3727     gmx_bool              bGem, bNN, bParallel;
3728     t_gemParams          *params = NULL;
3729     gmx_bool              bEdge_yjj, bEdge_xjj, bOMP;
3730
3731     t_hbdata            **p_hb    = NULL;                   /* one per thread, then merge after the frame loop */
3732     int                 **p_adist = NULL, **p_rdist = NULL; /* a histogram for each thread. */
3733
3734 #ifdef GMX_OPENMP
3735     bOMP = TRUE;
3736 #else
3737     bOMP = FALSE;
3738 #endif
3739
3740     npargs = asize(pa);
3741     ppa    = add_acf_pargs(&npargs, pa);
3742
3743     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT, NFILE, fnm, npargs,
3744                            ppa, asize(desc), desc, asize(bugs), bugs, &oenv))
3745     {
3746         return 0;
3747     }
3748
3749     /* NN-loop? If so, what estimator to use ?*/
3750     NN = 1;
3751     /* Outcommented for now DvdS 2010-07-13
3752        while (NN < NN_NR && gmx_strcasecmp(NNtype[0], NNtype[NN])!=0)
3753         NN++;
3754        if (NN == NN_NR)
3755         gmx_fatal(FARGS, "Invalid NN-loop type.");
3756      */
3757     bNN = FALSE;
3758     for (i = 2; bNN == FALSE && i < NN_NR; i++)
3759     {
3760         bNN = bNN || NN == i;
3761     }
3762
3763     if (NN > NN_NONE && bMerge)
3764     {
3765         bMerge = FALSE;
3766     }
3767
3768     /* geminate recombination? If so, which flavor? */
3769     gemmode = 1;
3770     while (gemmode < gemNR && gmx_strcasecmp(gemType[0], gemType[gemmode]) != 0)
3771     {
3772         gemmode++;
3773     }
3774     if (gemmode == gemNR)
3775     {
3776         gmx_fatal(FARGS, "Invalid recombination type.");
3777     }
3778
3779     bGem = FALSE;
3780     for (i = 2; bGem == FALSE && i < gemNR; i++)
3781     {
3782         bGem = bGem || gemmode == i;
3783     }
3784
3785     if (bGem)
3786     {
3787         printf("Geminate recombination: %s\n", gemType[gemmode]);
3788         if (bContact)
3789         {
3790             if (gemmode != gemDD)
3791             {
3792                 printf("Turning off -contact option...\n");
3793                 bContact = FALSE;
3794             }
3795         }
3796         else
3797         {
3798             if (gemmode == gemDD)
3799             {
3800                 printf("Turning on -contact option...\n");
3801                 bContact = TRUE;
3802             }
3803         }
3804         if (bMerge)
3805         {
3806             if (gemmode == gemAA)
3807             {
3808                 printf("Turning off -merge option...\n");
3809                 bMerge = FALSE;
3810             }
3811         }
3812         else
3813         {
3814             if (gemmode != gemAA)
3815             {
3816                 printf("Turning on -merge option...\n");
3817                 bMerge = TRUE;
3818             }
3819         }
3820     }
3821
3822     /* process input */
3823     bSelected = FALSE;
3824     ccut      = cos(acut*DEG2RAD);
3825
3826     if (bContact)
3827     {
3828         if (bSelected)
3829         {
3830             gmx_fatal(FARGS, "Can not analyze selected contacts.");
3831         }
3832         if (!bDA)
3833         {
3834             gmx_fatal(FARGS, "Can not analyze contact between H and A: turn off -noda");
3835         }
3836     }
3837
3838     /* Initiate main data structure! */
3839     bHBmap = (opt2bSet("-ac", NFILE, fnm) ||
3840               opt2bSet("-life", NFILE, fnm) ||
3841               opt2bSet("-hbn", NFILE, fnm) ||
3842               opt2bSet("-hbm", NFILE, fnm) ||
3843               bGem);
3844
3845     if (opt2bSet("-nhbdist", NFILE, fnm))
3846     {
3847         const char *leg[MAXHH+1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" };
3848         fpnhb = xvgropen(opt2fn("-nhbdist", NFILE, fnm),
3849                          "Number of donor-H with N HBs", output_env_get_xvgr_tlabel(oenv), "N", oenv);
3850         xvgr_legend(fpnhb, asize(leg), leg, oenv);
3851     }
3852
3853     hb = mk_hbdata(bHBmap, opt2bSet("-dan", NFILE, fnm), bMerge || bContact, bGem, gemmode);
3854
3855     /* get topology */
3856     read_tpx_top(ftp2fn(efTPR, NFILE, fnm), &ir, box, &natoms, NULL, NULL, NULL, &top);
3857
3858     snew(grpnames, grNR);
3859     snew(index, grNR);
3860     snew(isize, grNR);
3861     /* Make Donor-Acceptor table */
3862     snew(datable, top.atoms.nr);
3863     gen_datable(index[0], isize[0], datable, top.atoms.nr);
3864
3865     if (bSelected)
3866     {
3867         /* analyze selected hydrogen bonds */
3868         printf("Select group with selected atoms:\n");
3869         get_index(&(top.atoms), opt2fn("-sel", NFILE, fnm),
3870                   1, &nsel, index, grpnames);
3871         if (nsel % 3)
3872         {
3873             gmx_fatal(FARGS, "Number of atoms in group '%s' not a multiple of 3\n"
3874                       "and therefore cannot contain triplets of "
3875                       "Donor-Hydrogen-Acceptor", grpnames[0]);
3876         }
3877         bTwo = FALSE;
3878
3879         for (i = 0; (i < nsel); i += 3)
3880         {
3881             int dd = index[0][i];
3882             int aa = index[0][i+2];
3883             /* int */ hh = index[0][i+1];
3884             add_dh (&hb->d, dd, hh, i, datable);
3885             add_acc(&hb->a, aa, i);
3886             /* Should this be here ? */
3887             snew(hb->d.dptr, top.atoms.nr);
3888             snew(hb->a.aptr, top.atoms.nr);
3889             add_hbond(hb, dd, aa, hh, gr0, gr0, 0, bMerge, 0, bContact, peri);
3890         }
3891         printf("Analyzing %d selected hydrogen bonds from '%s'\n",
3892                isize[0], grpnames[0]);
3893     }
3894     else
3895     {
3896         /* analyze all hydrogen bonds: get group(s) */
3897         printf("Specify 2 groups to analyze:\n");
3898         get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
3899                   2, isize, index, grpnames);
3900
3901         /* check if we have two identical or two non-overlapping groups */
3902         bTwo = isize[0] != isize[1];
3903         for (i = 0; (i < isize[0]) && !bTwo; i++)
3904         {
3905             bTwo = index[0][i] != index[1][i];
3906         }
3907         if (bTwo)
3908         {
3909             printf("Checking for overlap in atoms between %s and %s\n",
3910                    grpnames[0], grpnames[1]);
3911             for (i = 0; i < isize[1]; i++)
3912             {
3913                 if (ISINGRP(datable[index[1][i]]))
3914                 {
3915                     gmx_fatal(FARGS, "Partial overlap between groups '%s' and '%s'",
3916                               grpnames[0], grpnames[1]);
3917                 }
3918             }
3919             /*
3920                printf("Checking for overlap in atoms between %s and %s\n",
3921                grpnames[0],grpnames[1]);
3922                for (i=0; i<isize[0]; i++)
3923                for (j=0; j<isize[1]; j++)
3924                if (index[0][i] == index[1][j])
3925                gmx_fatal(FARGS,"Partial overlap between groups '%s' and '%s'",
3926                grpnames[0],grpnames[1]);
3927              */
3928         }
3929         if (bTwo)
3930         {
3931             printf("Calculating %s "
3932                    "between %s (%d atoms) and %s (%d atoms)\n",
3933                    bContact ? "contacts" : "hydrogen bonds",
3934                    grpnames[0], isize[0], grpnames[1], isize[1]);
3935         }
3936         else
3937         {
3938             fprintf(stderr, "Calculating %s in %s (%d atoms)\n",
3939                     bContact ? "contacts" : "hydrogen bonds", grpnames[0], isize[0]);
3940         }
3941     }
3942     sfree(datable);
3943
3944     /* search donors and acceptors in groups */
3945     snew(datable, top.atoms.nr);
3946     for (i = 0; (i < grNR); i++)
3947     {
3948         if ( ((i == gr0) && !bSelected ) ||
3949              ((i == gr1) && bTwo ))
3950         {
3951             gen_datable(index[i], isize[i], datable, top.atoms.nr);
3952             if (bContact)
3953             {
3954                 search_acceptors(&top, isize[i], index[i], &hb->a, i,
3955                                  bNitAcc, TRUE, (bTwo && (i == gr0)) || !bTwo, datable);
3956                 search_donors   (&top, isize[i], index[i], &hb->d, i,
3957                                  TRUE, (bTwo && (i == gr1)) || !bTwo, datable);
3958             }
3959             else
3960             {
3961                 search_acceptors(&top, isize[i], index[i], &hb->a, i, bNitAcc, FALSE, TRUE, datable);
3962                 search_donors   (&top, isize[i], index[i], &hb->d, i, FALSE, TRUE, datable);
3963             }
3964             if (bTwo)
3965             {
3966                 clear_datable_grp(datable, top.atoms.nr);
3967             }
3968         }
3969     }
3970     sfree(datable);
3971     printf("Found %d donors and %d acceptors\n", hb->d.nrd, hb->a.nra);
3972     /*if (bSelected)
3973        snew(donors[gr0D], dons[gr0D].nrd);*/
3974
3975     donor_properties = open_donor_properties_file(opt2fn_null("-don", NFILE, fnm), hb, oenv);
3976
3977     if (bHBmap)
3978     {
3979         printf("Making hbmap structure...");
3980         /* Generate hbond data structure */
3981         mk_hbmap(hb);
3982         printf("done.\n");
3983     }
3984
3985 #ifdef HAVE_NN_LOOPS
3986     if (bNN)
3987     {
3988         mk_hbEmap(hb, 0);
3989     }
3990 #endif
3991
3992     if (bGem)
3993     {
3994         printf("Making per structure...");
3995         /* Generate hbond data structure */
3996         mk_per(hb);
3997         printf("done.\n");
3998     }
3999
4000     /* check input */
4001     bStop = FALSE;
4002     if (hb->d.nrd + hb->a.nra == 0)
4003     {
4004         printf("No Donors or Acceptors found\n");
4005         bStop = TRUE;
4006     }
4007     if (!bStop)
4008     {
4009         if (hb->d.nrd == 0)
4010         {
4011             printf("No Donors found\n");
4012             bStop = TRUE;
4013         }
4014         if (hb->a.nra == 0)
4015         {
4016             printf("No Acceptors found\n");
4017             bStop = TRUE;
4018         }
4019     }
4020     if (bStop)
4021     {
4022         gmx_fatal(FARGS, "Nothing to be done");
4023     }
4024
4025     shatom = 0;
4026     if (rshell > 0)
4027     {
4028         int      shisz;
4029         atom_id *shidx;
4030         char    *shgrpnm;
4031         /* get index group with atom for shell */
4032         do
4033         {
4034             printf("Select atom for shell (1 atom):\n");
4035             get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
4036                       1, &shisz, &shidx, &shgrpnm);
4037             if (shisz != 1)
4038             {
4039                 printf("group contains %d atoms, should be 1 (one)\n", shisz);
4040             }
4041         }
4042         while (shisz != 1);
4043         shatom = shidx[0];
4044         printf("Will calculate hydrogen bonds within a shell "
4045                "of %g nm around atom %i\n", rshell, shatom+1);
4046     }
4047
4048     /* Analyze trajectory */
4049     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
4050     if (natoms > top.atoms.nr)
4051     {
4052         gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)",
4053                   top.atoms.nr, natoms);
4054     }
4055
4056     bBox  = ir.ePBC != epbcNONE;
4057     grid  = init_grid(bBox, box, (rcut > r2cut) ? rcut : r2cut, ngrid);
4058     nabin = acut/abin;
4059     nrbin = rcut/rbin;
4060     snew(adist, nabin+1);
4061     snew(rdist, nrbin+1);
4062
4063     if (bGem && !bBox)
4064     {
4065         gmx_fatal(FARGS, "Can't do geminate recombination without periodic box.");
4066     }
4067
4068     bParallel = FALSE;
4069
4070 #ifndef GMX_OPENMP
4071 #define __ADIST adist
4072 #define __RDIST rdist
4073 #define __HBDATA hb
4074 #else /* GMX_OPENMP ================================================== \
4075        * Set up the OpenMP stuff,                                       |
4076        * like the number of threads and such                            |
4077        * Also start the parallel loop.                                  |
4078        */
4079 #define __ADIST p_adist[threadNr]
4080 #define __RDIST p_rdist[threadNr]
4081 #define __HBDATA p_hb[threadNr]
4082 #endif
4083     if (bOMP)
4084     {
4085         bParallel = !bSelected;
4086
4087         if (bParallel)
4088         {
4089             actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
4090
4091             gmx_omp_set_num_threads(actual_nThreads);
4092             printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads);
4093             fflush(stdout);
4094         }
4095         else
4096         {
4097             actual_nThreads = 1;
4098         }
4099
4100         snew(p_hb,    actual_nThreads);
4101         snew(p_adist, actual_nThreads);
4102         snew(p_rdist, actual_nThreads);
4103         for (i = 0; i < actual_nThreads; i++)
4104         {
4105             snew(p_hb[i], 1);
4106             snew(p_adist[i], nabin+1);
4107             snew(p_rdist[i], nrbin+1);
4108
4109             p_hb[i]->max_frames = 0;
4110             p_hb[i]->nhb        = NULL;
4111             p_hb[i]->ndist      = NULL;
4112             p_hb[i]->n_bound    = NULL;
4113             p_hb[i]->time       = NULL;
4114             p_hb[i]->nhx        = NULL;
4115
4116             p_hb[i]->bHBmap     = hb->bHBmap;
4117             p_hb[i]->bDAnr      = hb->bDAnr;
4118             p_hb[i]->bGem       = hb->bGem;
4119             p_hb[i]->wordlen    = hb->wordlen;
4120             p_hb[i]->nframes    = hb->nframes;
4121             p_hb[i]->maxhydro   = hb->maxhydro;
4122             p_hb[i]->danr       = hb->danr;
4123             p_hb[i]->d          = hb->d;
4124             p_hb[i]->a          = hb->a;
4125             p_hb[i]->hbmap      = hb->hbmap;
4126             p_hb[i]->time       = hb->time; /* This may need re-syncing at every frame. */
4127             p_hb[i]->per        = hb->per;
4128
4129 #ifdef HAVE_NN_LOOPS
4130             p_hb[i]->hbE = hb->hbE;
4131 #endif
4132
4133             p_hb[i]->nrhb   = 0;
4134             p_hb[i]->nrdist = 0;
4135         }
4136     }
4137
4138     /* Make a thread pool here,
4139      * instead of forking anew at every frame. */
4140
4141 #pragma omp parallel \
4142     firstprivate(i) \
4143     private(j, h, ii, jj, hh, E, \
4144     xi, yi, zi, xj, yj, zj, threadNr, \
4145     dist, ang, peri, icell, jcell, \
4146     grp, ogrp, ai, aj, xjj, yjj, zjj, \
4147     xk, yk, zk, ihb, id,  resdist, \
4148     xkk, ykk, zkk, kcell, ak, k, bTric, \
4149     bEdge_xjj, bEdge_yjj) \
4150     default(shared)
4151     {    /* Start of parallel region */
4152         threadNr = gmx_omp_get_thread_num();
4153
4154         do
4155         {
4156
4157             bTric = bBox && TRICLINIC(box);
4158
4159             if (bOMP)
4160             {
4161                 sync_hbdata(p_hb[threadNr], nframes);
4162             }
4163 #pragma omp single
4164             {
4165                 build_grid(hb, x, x[shatom], bBox, box, hbox, (rcut > r2cut) ? rcut : r2cut,
4166                            rshell, ngrid, grid);
4167                 reset_nhbonds(&(hb->d));
4168
4169                 if (debug && bDebug)
4170                 {
4171                     dump_grid(debug, ngrid, grid);
4172                 }
4173
4174                 add_frames(hb, nframes);
4175                 init_hbframe(hb, nframes, output_env_conv_time(oenv, t));
4176
4177                 if (hb->bDAnr)
4178                 {
4179                     count_da_grid(ngrid, grid, hb->danr[nframes]);
4180                 }
4181             } /* omp single */
4182
4183             if (bOMP)
4184             {
4185                 p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */
4186             }
4187
4188             if (bNN)
4189             {
4190 #ifdef HAVE_NN_LOOPS /* Unlock this feature when testing */
4191                 /* Loop over all atom pairs and estimate interaction energy */
4192
4193 #pragma omp single
4194                 {
4195                     addFramesNN(hb, nframes);
4196                 }
4197
4198 #pragma omp barrier
4199 #pragma omp for schedule(dynamic)
4200                 for (i = 0; i < hb->d.nrd; i++)
4201                 {
4202                     for (j = 0; j < hb->a.nra; j++)
4203                     {
4204                         for (h = 0;
4205                              h < (bContact ? 1 : hb->d.nhydro[i]);
4206                              h++)
4207                         {
4208                             if (i == hb->d.nrd || j == hb->a.nra)
4209                             {
4210                                 gmx_fatal(FARGS, "out of bounds");
4211                             }
4212
4213                             /* Get the real atom ids */
4214                             ii = hb->d.don[i];
4215                             jj = hb->a.acc[j];
4216                             hh = hb->d.hydro[i][h];
4217
4218                             /* Estimate the energy from the geometry */
4219                             E = calcHbEnergy(ii, jj, hh, x, NN, box, hbox, &(hb->d));
4220                             /* Store the energy */
4221                             storeHbEnergy(hb, i, j, h, E, nframes);
4222                         }
4223                     }
4224                 }
4225 #endif        /* HAVE_NN_LOOPS */
4226             } /* if (bNN)*/
4227             else
4228             {
4229                 if (bSelected)
4230                 {
4231
4232 #pragma omp single
4233                     {
4234                         /* Do not parallelize this just yet. */
4235                         /* int ii; */
4236                         for (ii = 0; (ii < nsel); ii++)
4237                         {
4238                             int dd = index[0][i];
4239                             int aa = index[0][i+2];
4240                             /* int */ hh = index[0][i+1];
4241                             ihb          = is_hbond(hb, ii, ii, dd, aa, rcut, r2cut, ccut, x, bBox, box,
4242                                                     hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri);
4243
4244                             if (ihb)
4245                             {
4246                                 /* add to index if not already there */
4247                                 /* Add a hbond */
4248                                 add_hbond(hb, dd, aa, hh, ii, ii, nframes, bMerge, ihb, bContact, peri);
4249                             }
4250                         }
4251                     } /* omp single */
4252                 }     /* if (bSelected) */
4253                 else
4254                 {
4255
4256 #pragma omp single
4257                     {
4258                         if (bGem)
4259                         {
4260                             calcBoxProjection(box, hb->per->P);
4261                         }
4262
4263                         /* loop over all gridcells (xi,yi,zi)      */
4264                         /* Removed confusing macro, DvdS 27/12/98  */
4265
4266                     }
4267                     /* The outer grid loop will have to do for now. */
4268 #pragma omp for schedule(dynamic)
4269                     for (xi = 0; xi < ngrid[XX]; xi++)
4270                     {
4271                         for (yi = 0; (yi < ngrid[YY]); yi++)
4272                         {
4273                             for (zi = 0; (zi < ngrid[ZZ]); zi++)
4274                             {
4275
4276                                 /* loop over donor groups gr0 (always) and gr1 (if necessary) */
4277                                 for (grp = gr0; (grp <= (bTwo ? gr1 : gr0)); grp++)
4278                                 {
4279                                     icell = &(grid[zi][yi][xi].d[grp]);
4280
4281                                     if (bTwo)
4282                                     {
4283                                         ogrp = 1-grp;
4284                                     }
4285                                     else
4286                                     {
4287                                         ogrp = grp;
4288                                     }
4289
4290                                     /* loop over all hydrogen atoms from group (grp)
4291                                      * in this gridcell (icell)
4292                                      */
4293                                     for (ai = 0; (ai < icell->nr); ai++)
4294                                     {
4295                                         i  = icell->atoms[ai];
4296
4297                                         /* loop over all adjacent gridcells (xj,yj,zj) */
4298                                         for (zjj = grid_loop_begin(ngrid[ZZ], zi, bTric, FALSE);
4299                                              zjj <= grid_loop_end(ngrid[ZZ], zi, bTric, FALSE);
4300                                              zjj++)
4301                                         {
4302                                             zj        = grid_mod(zjj, ngrid[ZZ]);
4303                                             bEdge_yjj = (zj == 0) || (zj == ngrid[ZZ] - 1);
4304                                             for (yjj = grid_loop_begin(ngrid[YY], yi, bTric, bEdge_yjj);
4305                                                  yjj <= grid_loop_end(ngrid[YY], yi, bTric, bEdge_yjj);
4306                                                  yjj++)
4307                                             {
4308                                                 yj        = grid_mod(yjj, ngrid[YY]);
4309                                                 bEdge_xjj =
4310                                                     (yj == 0) || (yj == ngrid[YY] - 1) ||
4311                                                     (zj == 0) || (zj == ngrid[ZZ] - 1);
4312                                                 for (xjj = grid_loop_begin(ngrid[XX], xi, bTric, bEdge_xjj);
4313                                                      xjj <= grid_loop_end(ngrid[XX], xi, bTric, bEdge_xjj);
4314                                                      xjj++)
4315                                                 {
4316                                                     xj    = grid_mod(xjj, ngrid[XX]);
4317                                                     jcell = &(grid[zj][yj][xj].a[ogrp]);
4318                                                     /* loop over acceptor atoms from other group (ogrp)
4319                                                      * in this adjacent gridcell (jcell)
4320                                                      */
4321                                                     for (aj = 0; (aj < jcell->nr); aj++)
4322                                                     {
4323                                                         j = jcell->atoms[aj];
4324
4325                                                         /* check if this once was a h-bond */
4326                                                         peri = -1;
4327                                                         ihb  = is_hbond(__HBDATA, grp, ogrp, i, j, rcut, r2cut, ccut, x, bBox, box,
4328                                                                         hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri);
4329
4330                                                         if (ihb)
4331                                                         {
4332                                                             /* add to index if not already there */
4333                                                             /* Add a hbond */
4334                                                             add_hbond(__HBDATA, i, j, h, grp, ogrp, nframes, bMerge, ihb, bContact, peri);
4335
4336                                                             /* make angle and distance distributions */
4337                                                             if (ihb == hbHB && !bContact)
4338                                                             {
4339                                                                 if (dist > rcut)
4340                                                                 {
4341                                                                     gmx_fatal(FARGS, "distance is higher than what is allowed for an hbond: %f", dist);
4342                                                                 }
4343                                                                 ang *= RAD2DEG;
4344                                                                 __ADIST[(int)( ang/abin)]++;
4345                                                                 __RDIST[(int)(dist/rbin)]++;
4346                                                                 if (!bTwo)
4347                                                                 {
4348                                                                     int id, ia;
4349                                                                     if ((id = donor_index(&hb->d, grp, i)) == NOTSET)
4350                                                                     {
4351                                                                         gmx_fatal(FARGS, "Invalid donor %d", i);
4352                                                                     }
4353                                                                     if ((ia = acceptor_index(&hb->a, ogrp, j)) == NOTSET)
4354                                                                     {
4355                                                                         gmx_fatal(FARGS, "Invalid acceptor %d", j);
4356                                                                     }
4357                                                                     resdist = abs(top.atoms.atom[i].resind-
4358                                                                                   top.atoms.atom[j].resind);
4359                                                                     if (resdist >= max_hx)
4360                                                                     {
4361                                                                         resdist = max_hx-1;
4362                                                                     }
4363                                                                     __HBDATA->nhx[nframes][resdist]++;
4364                                                                 }
4365                                                             }
4366
4367                                                         }
4368                                                     } /* for aj  */
4369                                                 }     /* for xjj */
4370                                             }         /* for yjj */
4371                                         }             /* for zjj */
4372                                     }                 /* for ai  */
4373                                 }                     /* for grp */
4374                             }                         /* for xi,yi,zi */
4375                         }
4376                     }
4377                 } /* if (bSelected) {...} else */
4378
4379
4380                 /* Better wait for all threads to finnish using x[] before updating it. */
4381                 k = nframes;
4382 #pragma omp barrier
4383 #pragma omp critical
4384                 {
4385                     /* Sum up histograms and counts from p_hb[] into hb */
4386                     if (bOMP)
4387                     {
4388                         hb->nhb[k]   += p_hb[threadNr]->nhb[k];
4389                         hb->ndist[k] += p_hb[threadNr]->ndist[k];
4390                         for (j = 0; j < max_hx; j++)
4391                         {
4392                             hb->nhx[k][j]  += p_hb[threadNr]->nhx[k][j];
4393                         }
4394                     }
4395                 }
4396
4397                 /* Here are a handful of single constructs
4398                  * to share the workload a bit. The most
4399                  * important one is of course the last one,
4400                  * where there's a potential bottleneck in form
4401                  * of slow I/O.                    */
4402 #pragma omp barrier
4403 #pragma omp single
4404                 {
4405                     analyse_donor_properties(donor_properties, hb, k, t);
4406                 }
4407
4408 #pragma omp single
4409                 {
4410                     if (fpnhb)
4411                     {
4412                         do_nhb_dist(fpnhb, hb, t);
4413                     }
4414                 }
4415             } /* if (bNN) {...} else  +   */
4416
4417 #pragma omp single
4418             {
4419                 trrStatus = (read_next_x(oenv, status, &t, x, box));
4420                 nframes++;
4421             }
4422
4423 #pragma omp barrier
4424         }
4425         while (trrStatus);
4426
4427         if (bOMP)
4428         {
4429 #pragma omp critical
4430             {
4431                 hb->nrhb   += p_hb[threadNr]->nrhb;
4432                 hb->nrdist += p_hb[threadNr]->nrdist;
4433             }
4434             /* Free parallel datastructures */
4435             sfree(p_hb[threadNr]->nhb);
4436             sfree(p_hb[threadNr]->ndist);
4437             sfree(p_hb[threadNr]->nhx);
4438
4439 #pragma omp for
4440             for (i = 0; i < nabin; i++)
4441             {
4442                 for (j = 0; j < actual_nThreads; j++)
4443                 {
4444
4445                     adist[i] += p_adist[j][i];
4446                 }
4447             }
4448 #pragma omp for
4449             for (i = 0; i <= nrbin; i++)
4450             {
4451                 for (j = 0; j < actual_nThreads; j++)
4452                 {
4453                     rdist[i] += p_rdist[j][i];
4454                 }
4455             }
4456
4457             sfree(p_adist[threadNr]);
4458             sfree(p_rdist[threadNr]);
4459         }
4460     } /* End of parallel region */
4461     if (bOMP)
4462     {
4463         sfree(p_adist);
4464         sfree(p_rdist);
4465     }
4466
4467     if (nframes < 2 && (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm)))
4468     {
4469         gmx_fatal(FARGS, "Cannot calculate autocorrelation of life times with less than two frames");
4470     }
4471
4472     free_grid(ngrid, &grid);
4473
4474     close_trj(status);
4475
4476     if (donor_properties)
4477     {
4478         xvgrclose(donor_properties);
4479     }
4480
4481     if (fpnhb)
4482     {
4483         xvgrclose(fpnhb);
4484     }
4485
4486     /* Compute maximum possible number of different hbonds */
4487     if (maxnhb > 0)
4488     {
4489         max_nhb = maxnhb;
4490     }
4491     else
4492     {
4493         max_nhb = 0.5*(hb->d.nrd*hb->a.nra);
4494     }
4495     /* Added support for -contact below.
4496      * - Erik Marklund, May 29-31, 2006 */
4497     /* Changed contact code.
4498      * - Erik Marklund, June 29, 2006 */
4499     if (bHBmap && !bNN)
4500     {
4501         if (hb->nrhb == 0)
4502         {
4503             printf("No %s found!!\n", bContact ? "contacts" : "hydrogen bonds");
4504         }
4505         else
4506         {
4507             printf("Found %d different %s in trajectory\n"
4508                    "Found %d different atom-pairs within %s distance\n",
4509                    hb->nrhb, bContact ? "contacts" : "hydrogen bonds",
4510                    hb->nrdist, (r2cut > 0) ? "second cut-off" : "hydrogen bonding");
4511
4512             /*Control the pHist.*/
4513
4514             if (bMerge)
4515             {
4516                 merge_hb(hb, bTwo, bContact);
4517             }
4518
4519             if (opt2bSet("-hbn", NFILE, fnm))
4520             {
4521                 dump_hbmap(hb, NFILE, fnm, bTwo, bContact, isize, index, grpnames, &top.atoms);
4522             }
4523
4524             /* Moved the call to merge_hb() to a line BEFORE dump_hbmap
4525              * to make the -hbn and -hmb output match eachother.
4526              * - Erik Marklund, May 30, 2006 */
4527         }
4528     }
4529     /* Print out number of hbonds and distances */
4530     aver_nhb  = 0;
4531     aver_dist = 0;
4532     fp        = xvgropen(opt2fn("-num", NFILE, fnm), bContact ? "Contacts" :
4533                          "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
4534     snew(leg, 2);
4535     snew(leg[0], STRLEN);
4536     snew(leg[1], STRLEN);
4537     sprintf(leg[0], "%s", bContact ? "Contacts" : "Hydrogen bonds");
4538     sprintf(leg[1], "Pairs within %g nm", (r2cut > 0) ? r2cut : rcut);
4539     xvgr_legend(fp, 2, (const char**)leg, oenv);
4540     sfree(leg[1]);
4541     sfree(leg[0]);
4542     sfree(leg);
4543     for (i = 0; (i < nframes); i++)
4544     {
4545         fprintf(fp, "%10g  %10d  %10d\n", hb->time[i], hb->nhb[i], hb->ndist[i]);
4546         aver_nhb  += hb->nhb[i];
4547         aver_dist += hb->ndist[i];
4548     }
4549     xvgrclose(fp);
4550     aver_nhb  /= nframes;
4551     aver_dist /= nframes;
4552     /* Print HB distance distribution */
4553     if (opt2bSet("-dist", NFILE, fnm))
4554     {
4555         long sum;
4556
4557         sum = 0;
4558         for (i = 0; i < nrbin; i++)
4559         {
4560             sum += rdist[i];
4561         }
4562
4563         fp = xvgropen(opt2fn("-dist", NFILE, fnm),
4564                       "Hydrogen Bond Distribution",
4565                       bDA ?
4566                       "Donor - Acceptor Distance (nm)" :
4567                       "Hydrogen - Acceptor Distance (nm)", "", oenv);
4568         for (i = 0; i < nrbin; i++)
4569         {
4570             fprintf(fp, "%10g %10g\n", (i+0.5)*rbin, rdist[i]/(rbin*(real)sum));
4571         }
4572         xvgrclose(fp);
4573     }
4574
4575     /* Print HB angle distribution */
4576     if (opt2bSet("-ang", NFILE, fnm))
4577     {
4578         long sum;
4579
4580         sum = 0;
4581         for (i = 0; i < nabin; i++)
4582         {
4583             sum += adist[i];
4584         }
4585
4586         fp = xvgropen(opt2fn("-ang", NFILE, fnm),
4587                       "Hydrogen Bond Distribution",
4588                       "Hydrogen - Donor - Acceptor Angle (\\SO\\N)", "", oenv);
4589         for (i = 0; i < nabin; i++)
4590         {
4591             fprintf(fp, "%10g %10g\n", (i+0.5)*abin, adist[i]/(abin*(real)sum));
4592         }
4593         xvgrclose(fp);
4594     }
4595
4596     /* Print HB in alpha-helix */
4597     if (opt2bSet("-hx", NFILE, fnm))
4598     {
4599         fp = xvgropen(opt2fn("-hx", NFILE, fnm),
4600                       "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
4601         xvgr_legend(fp, NRHXTYPES, hxtypenames, oenv);
4602         for (i = 0; i < nframes; i++)
4603         {
4604             fprintf(fp, "%10g", hb->time[i]);
4605             for (j = 0; j < max_hx; j++)
4606             {
4607                 fprintf(fp, " %6d", hb->nhx[i][j]);
4608             }
4609             fprintf(fp, "\n");
4610         }
4611         xvgrclose(fp);
4612     }
4613     if (!bNN)
4614     {
4615         printf("Average number of %s per timeframe %.3f out of %g possible\n",
4616                bContact ? "contacts" : "hbonds",
4617                bContact ? aver_dist : aver_nhb, max_nhb);
4618     }
4619
4620     /* Do Autocorrelation etc. */
4621     if (hb->bHBmap)
4622     {
4623         /*
4624            Added support for -contact in ac and hbm calculations below.
4625            - Erik Marklund, May 29, 2006
4626          */
4627         ivec itmp;
4628         rvec rtmp;
4629         if (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm))
4630         {
4631             please_cite(stdout, "Spoel2006b");
4632         }
4633         if (opt2bSet("-ac", NFILE, fnm))
4634         {
4635             char *gemstring = NULL;
4636
4637             if (bGem || bNN)
4638             {
4639                 params = init_gemParams(rcut, D, hb->time, hb->nframes/2, nFitPoints, fit_start, fit_end,
4640                                         gemBallistic, nBalExp);
4641                 if (params == NULL)
4642                 {
4643                     gmx_fatal(FARGS, "Could not initiate t_gemParams params.");
4644                 }
4645             }
4646             gemstring = gmx_strdup(gemType[hb->per->gemtype]);
4647             do_hbac(opt2fn("-ac", NFILE, fnm), hb, nDump,
4648                     bMerge, bContact, fit_start, temp, r2cut > 0, smooth_tail_start, oenv,
4649                     gemstring, nThreads, NN, bBallistic, bGemFit);
4650         }
4651         if (opt2bSet("-life", NFILE, fnm))
4652         {
4653             do_hblife(opt2fn("-life", NFILE, fnm), hb, bMerge, bContact, oenv);
4654         }
4655         if (opt2bSet("-hbm", NFILE, fnm))
4656         {
4657             t_matrix mat;
4658             int      id, ia, hh, x, y;
4659             mat.flags = mat.y0 = 0;
4660
4661             if ((nframes > 0) && (hb->nrhb > 0))
4662             {
4663                 mat.nx = nframes;
4664                 mat.ny = hb->nrhb;
4665
4666                 snew(mat.matrix, mat.nx);
4667                 for (x = 0; (x < mat.nx); x++)
4668                 {
4669                     snew(mat.matrix[x], mat.ny);
4670                 }
4671                 y = 0;
4672                 for (id = 0; (id < hb->d.nrd); id++)
4673                 {
4674                     for (ia = 0; (ia < hb->a.nra); ia++)
4675                     {
4676                         for (hh = 0; (hh < hb->maxhydro); hh++)
4677                         {
4678                             if (hb->hbmap[id][ia])
4679                             {
4680                                 if (ISHB(hb->hbmap[id][ia]->history[hh]))
4681                                 {
4682                                     /* Changed '<' into '<=' in the for-statement below.
4683                                      * It fixed the previously undiscovered bug that caused
4684                                      * the last occurance of an hbond/contact to not be
4685                                      * set in mat.matrix. Have a look at any old -hbm-output
4686                                      * and you will notice that the last column is allways empty.
4687                                      * - Erik Marklund May 30, 2006
4688                                      */
4689                                     for (x = 0; (x <= hb->hbmap[id][ia]->nframes); x++)
4690                                     {
4691                                         int nn0 = hb->hbmap[id][ia]->n0;
4692                                         range_check(y, 0, mat.ny);
4693                                         mat.matrix[x+nn0][y] = is_hb(hb->hbmap[id][ia]->h[hh], x);
4694                                     }
4695                                     y++;
4696                                 }
4697                             }
4698                         }
4699                     }
4700                 }
4701                 mat.axis_x = hb->time;
4702                 snew(mat.axis_y, mat.ny);
4703                 for (j = 0; j < mat.ny; j++)
4704                 {
4705                     mat.axis_y[j] = j;
4706                 }
4707                 sprintf(mat.title, bContact ? "Contact Existence Map" :
4708                         "Hydrogen Bond Existence Map");
4709                 sprintf(mat.legend, bContact ? "Contacts" : "Hydrogen Bonds");
4710                 sprintf(mat.label_x, "%s", output_env_get_xvgr_tlabel(oenv));
4711                 sprintf(mat.label_y, bContact ? "Contact Index" : "Hydrogen Bond Index");
4712                 mat.bDiscrete = TRUE;
4713                 mat.nmap      = 2;
4714                 snew(mat.map, mat.nmap);
4715                 for (i = 0; i < mat.nmap; i++)
4716                 {
4717                     mat.map[i].code.c1 = hbmap[i];
4718                     mat.map[i].desc    = hbdesc[i];
4719                     mat.map[i].rgb     = hbrgb[i];
4720                 }
4721                 fp = opt2FILE("-hbm", NFILE, fnm, "w");
4722                 write_xpm_m(fp, mat);
4723                 gmx_ffclose(fp);
4724                 for (x = 0; x < mat.nx; x++)
4725                 {
4726                     sfree(mat.matrix[x]);
4727                 }
4728                 sfree(mat.axis_y);
4729                 sfree(mat.matrix);
4730                 sfree(mat.map);
4731             }
4732             else
4733             {
4734                 fprintf(stderr, "No hydrogen bonds/contacts found. No hydrogen bond map will be printed.\n");
4735             }
4736         }
4737     }
4738
4739     if (bGem)
4740     {
4741         fprintf(stderr, "There were %i periodic shifts\n", hb->per->nper);
4742         fprintf(stderr, "Freeing pHist for all donors...\n");
4743         for (i = 0; i < hb->d.nrd; i++)
4744         {
4745             fprintf(stderr, "\r%i", i);
4746             if (hb->per->pHist[i] != NULL)
4747             {
4748                 for (j = 0; j < hb->a.nra; j++)
4749                 {
4750                     clearPshift(&(hb->per->pHist[i][j]));
4751                 }
4752                 sfree(hb->per->pHist[i]);
4753             }
4754         }
4755         sfree(hb->per->pHist);
4756         sfree(hb->per->p2i);
4757         sfree(hb->per);
4758         fprintf(stderr, "...done.\n");
4759     }
4760
4761 #ifdef HAVE_NN_LOOPS
4762     if (bNN)
4763     {
4764         free_hbEmap(hb);
4765     }
4766 #endif
4767
4768     if (hb->bDAnr)
4769     {
4770         int    i, j, nleg;
4771         char **legnames;
4772         char   buf[STRLEN];
4773
4774 #define USE_THIS_GROUP(j) ( (j == gr0) || (bTwo && (j == gr1)) )
4775
4776         fp = xvgropen(opt2fn("-dan", NFILE, fnm),
4777                       "Donors and Acceptors", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
4778         nleg = (bTwo ? 2 : 1)*2;
4779         snew(legnames, nleg);
4780         i = 0;
4781         for (j = 0; j < grNR; j++)
4782         {
4783             if (USE_THIS_GROUP(j) )
4784             {
4785                 sprintf(buf, "Donors %s", grpnames[j]);
4786                 legnames[i++] = gmx_strdup(buf);
4787                 sprintf(buf, "Acceptors %s", grpnames[j]);
4788                 legnames[i++] = gmx_strdup(buf);
4789             }
4790         }
4791         if (i != nleg)
4792         {
4793             gmx_incons("number of legend entries");
4794         }
4795         xvgr_legend(fp, nleg, (const char**)legnames, oenv);
4796         for (i = 0; i < nframes; i++)
4797         {
4798             fprintf(fp, "%10g", hb->time[i]);
4799             for (j = 0; (j < grNR); j++)
4800             {
4801                 if (USE_THIS_GROUP(j) )
4802                 {
4803                     fprintf(fp, " %6d", hb->danr[i][j]);
4804                 }
4805             }
4806             fprintf(fp, "\n");
4807         }
4808         xvgrclose(fp);
4809     }
4810
4811     return 0;
4812 }