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