src/gromacs/gmxana/gmx_hbond.c

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