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