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