src/tools/gmx_hbond.c

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