src/tools/gmx_hbond.c

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