A freeze group can now be allowed to move rigidly in some dimension by using "freezed...

[gromacs/rigid-bodies.git] / src / tools / gmx_hbond.c

/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
 * 
 * 
 *                This source code is part of
 * 
 *                 G   R   O   M   A   C   S
 * 
 *          GROningen MAchine for Chemical Simulations
 * 
 *                        VERSION 3.3.3
 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
 * Copyright (c) 2001-2008, The GROMACS development team,
 * check out http://www.gromacs.org for more information.

 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 * 
 * If you want to redistribute modifications, please consider that
 * scientific software is very special. Version control is crucial -
 * bugs must be traceable. We will be happy to consider code for
 * inclusion in the official distribution, but derived work must not
 * be called official GROMACS. Details are found in the README & COPYING
 * files - if they are missing, get the official version at www.gromacs.org.
 * 
 * To help us fund GROMACS development, we humbly ask that you cite
 * the papers on the package - you can find them in the top README file.
 * 
 * For more info, check our website at http://www.gromacs.org
 * 
 * And Hey:
 * Groningen Machine for Chemical Simulation
 */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <math.h>

/*#define HAVE_NN_LOOPS*/
/* Set environment variable CFLAGS = "-fopenmp" when running
 * configure and define DOUSEOPENMP to make use of parallelized
 * calculation of autocorrelation function.
 * It also adds a new option -nthreads which sets the number of threads.
 * */
/*#define DOUSEOPENMP*/

#ifdef DOUSEOPENMP
#define HAVE_OPENMP
#include "omp.h"
#endif

#include "statutil.h"
#include "copyrite.h"
#include "sysstuff.h"
#include "txtdump.h"
#include "futil.h"
#include "tpxio.h"
#include "physics.h"
#include "macros.h"
#include "gmx_fatal.h"
#include "index.h"
#include "smalloc.h"
#include "vec.h"
#include "xvgr.h"
#include "gstat.h"
#include "matio.h"
#include "string2.h"
#include "pbc.h"
#include "correl.h"
#include "gmx_ana.h"
#include "geminate.h"

typedef short int t_E;
typedef int t_EEst;
#define max_hx 7
typedef int t_hx[max_hx];
#define NRHXTYPES max_hx
const char *hxtypenames[NRHXTYPES]=
{"n-n","n-n+1","n-n+2","n-n+3","n-n+4","n-n+5","n-n>6"};
#define MAXHH 4

#ifdef HAVE_OPENMP
#define MASTER_THREAD_ONLY(threadNr) ((threadNr)==0)
#else
#define MASTER_THREAD_ONLY(threadNr) ((threadNr)==(threadNr))
#endif

/* -----------------------------------------*/

enum { gr0,  gr1,    grI,  grNR };
enum { hbNo, hbDist, hbHB, hbNR, hbR2}; 
enum { noDA, ACC, DON, DA, INGROUP};
enum {NN_NULL, NN_NONE, NN_BINARY, NN_1_over_r3, NN_dipole, NN_NR};
  
static const char *grpnames[grNR] = {"0","1","I" };

static gmx_bool bDebug = FALSE;

#define HB_NO 0
#define HB_YES 1<<0
#define HB_INS 1<<1
#define HB_YESINS HB_YES|HB_INS
#define HB_NR (1<<2)
#define MAXHYDRO 4

#define ISHB(h)   (((h) & 2) == 2)
#define ISDIST(h) (((h) & 1) == 1)
#define ISDIST2(h) (((h) & 4) == 4)
#define ISACC(h)  (((h) & 1) == 1)
#define ISDON(h)  (((h) & 2) == 2)
#define ISINGRP(h) (((h) & 4) == 4)

typedef struct {
    int nr;
    int maxnr;
    atom_id *atoms;
} t_ncell;

typedef struct {
    t_ncell d[grNR];
    t_ncell a[grNR];
} t_gridcell;

typedef int     t_icell[grNR];
typedef atom_id h_id[MAXHYDRO];
 
typedef struct {
    int      history[MAXHYDRO]; 
    /* Has this hbond existed ever? If so as hbDist or hbHB or both.
     * Result is stored as a bitmap (1 = hbDist) || (2 = hbHB)
     */
    /* Bitmask array which tells whether a hbond is present
     * at a given time. Either of these may be NULL 
     */
    int      n0;       /* First frame a HB was found     */ 
    int      nframes,maxframes;  /* Amount of frames in this hbond */
    unsigned int **h; 
    unsigned int **g; 
    /* See Xu and Berne, JPCB 105 (2001), p. 11929. We define the
     * function g(t) = [1-h(t)] H(t) where H(t) is one when the donor-
     * acceptor distance is less than the user-specified distance (typically
     * 0.35 nm).
     */
} t_hbond;

typedef struct {
    int     nra,max_nra;
    atom_id *acc;             /* Atom numbers of the acceptors     */
    int     *grp;             /* Group index                       */
    int     *aptr;            /* Map atom number to acceptor index */
} t_acceptors;

typedef struct {
    int      nrd,max_nrd;
    int      *don;               /* Atom numbers of the donors         */
    int      *grp;               /* Group index                        */
    int      *dptr;              /* Map atom number to donor index     */
    int      *nhydro;            /* Number of hydrogens for each donor */
    h_id     *hydro;             /* The atom numbers of the hydrogens  */
    h_id     *nhbonds;           /* The number of HBs per H at current */
} t_donors;

/* Tune this to match memory requirements. It should be a signed integer type, e.g. signed char.*/
#define PSTYPE int

typedef struct {
    int len;     /* The length of frame and p. */
    int *frame;  /* The frames at which transitio*/
    PSTYPE *p;
} t_pShift;

typedef struct {
    /* Periodicity history. Used for the reversible geminate recombination. */
    t_pShift **pHist;            /* The periodicity of every hbond in t_hbdata->hbmap:
                                  *   pHist[d][a]. We can safely assume that the same
                                  *   periodic shift holds for all hydrogens of a da-pair.
                                  *
                                  * Nowadays it only stores TRANSITIONS, and not the shift at every frame.
                                  *   That saves a LOT of memory, an hopefully kills a mysterious bug where
                                  *   pHist gets contaminated. */
  
    PSTYPE nper;    /* The length of p2i */
    ivec *p2i;       /* Maps integer to periodic shift for a pair.*/
    matrix P;        /* Projection matrix to find the box shifts. */
    int gemtype;     /* enumerated type */
} t_gemPeriod;

typedef struct {
    int nframes;
    int *Etot;  /* Total energy for each frame */
    t_E ****E;  /* Energy estimate for [d][a][h][frame-n0] */
} t_hbEmap;

typedef struct {
    gmx_bool        bHBmap,bDAnr,bGem;
    int         wordlen;
    /* The following arrays are nframes long */
    int         nframes,max_frames,maxhydro;
    int         *nhb,*ndist;
    h_id        *n_bound;
    real        *time;
    t_icell     *danr;
    t_hx        *nhx;
    /* These structures are initialized from the topology at start up */
    t_donors    d;
    t_acceptors a;
    /* This holds a matrix with all possible hydrogen bonds */
    int         nrhb,nrdist;
    t_hbond     ***hbmap;
#ifdef HAVE_NN_LOOPS
    t_hbEmap    hbE;
#endif
    /* For parallelization reasons this will have to be a pointer.
     * Otherwise discrepancies may arise between the periodicity data
     * seen by different threads. */
    t_gemPeriod *per;
} t_hbdata;

static void clearPshift(t_pShift *pShift)
{
    if (pShift->len > 0) {
        sfree(pShift->p);
        sfree(pShift->frame);
        pShift->len = 0;
    }
}

static void calcBoxProjection(matrix B, matrix P)
{
    const int vp[] = {XX,YY,ZZ};
    int i,j;
    int m,n;
    matrix M, N, U;

    for (i=0; i<3; i++){ m = vp[i];
        for (j=0; j<3; j++){ n = vp[j];
            U[m][n] = i==j ? 1:0;
        }
    }
    m_inv(B,M);
    for (i=0; i<3; i++){ m = vp[i];
        mvmul(M, U[m], P[m]);
    }
    transpose(P,N);
}

static void calcBoxDistance(matrix P, rvec d, ivec ibd){
    /* returns integer distance in box coordinates.
     * P is the projection matrix from cartesian coordinates
     * obtained with calcBoxProjection(). */
    int i;
    rvec bd;
    mvmul(P, d, bd);
    /* extend it by 0.5 in all directions since (int) rounds toward 0.*/
    for (i=0;i<3;i++)
        bd[i]=bd[i] + (bd[i]<0 ? -0.5 : 0.5);
    ibd[XX] = (int)bd[XX];
    ibd[YY] = (int)bd[YY];
    ibd[ZZ] = (int)bd[ZZ];
}

/* Changed argument 'bMerge' into 'oneHB' below,
 * since -contact should cause maxhydro to be 1,
 * not just -merge.
 * - Erik Marklund May 29, 2006
 */

static PSTYPE periodicIndex(ivec r, t_gemPeriod *per, gmx_bool daSwap) {
    /* Try to merge hbonds on the fly. That means that if the
     * acceptor and donor are mergable, then:
     * 1) store the hb-info so that acceptor id > donor id,
     * 2) add the periodic shift in pairs, so that [-x,-y,-z] is
     *    stored in per.p2i[] whenever acceptor id < donor id.
     * Note that [0,0,0] should already be the first element of per.p2i
     * by the time this function is called. */

    /* daSwap is TRUE if the donor and acceptor were swapped.
     * If so, then the negative vector should be used. */
    PSTYPE i;

    if (per->p2i == NULL || per->nper == 0)
        gmx_fatal(FARGS, "'per' not initialized properly.");
    for (i=0; i<per->nper; i++) {
        if (r[XX] == per->p2i[i][XX] &&
            r[YY] == per->p2i[i][YY] &&
            r[ZZ] == per->p2i[i][ZZ])
            return i;
    }
    /* Not found apparently. Add it to the list! */
    /* printf("New shift found: %i,%i,%i\n",r[XX],r[YY],r[ZZ]); */

/* Unfortunately this needs to be critical it seems. */
#ifdef HAVE_OPENMP
#pragma omp critical
#endif
    {
        if (!per->p2i) {
            fprintf(stderr, "p2i not initialized. This shouldn't happen!\n");
            snew(per->p2i, 1);
        }
        else
            srenew(per->p2i, per->nper+2);
        copy_ivec(r, per->p2i[per->nper]);
        (per->nper)++;

        /* Add the mirror too. It's rather likely that it'll be needed. */
        per->p2i[per->nper][XX] = -r[XX];
        per->p2i[per->nper][YY] = -r[YY];
        per->p2i[per->nper][ZZ] = -r[ZZ];
        (per->nper)++;
    } /* omp critical */
    return per->nper - 1 - (daSwap ? 0:1);
}

static t_hbdata *mk_hbdata(gmx_bool bHBmap,gmx_bool bDAnr,gmx_bool oneHB, gmx_bool bGem, int gemmode)
{
    t_hbdata *hb;
  
    snew(hb,1);
    hb->wordlen = 8*sizeof(unsigned int);
    hb->bHBmap  = bHBmap;
    hb->bDAnr   = bDAnr;
    hb->bGem    = bGem;
    if (oneHB)
        hb->maxhydro = 1;
    else
        hb->maxhydro = MAXHYDRO;
    snew(hb->per, 1);
    hb->per->gemtype = bGem ? gemmode : 0;
  
    return hb;
}

static void mk_hbmap(t_hbdata *hb,gmx_bool bTwo)
{
    int  i,j;

    snew(hb->hbmap,hb->d.nrd);
    for(i=0; (i<hb->d.nrd); i++) {
        snew(hb->hbmap[i],hb->a.nra);
        if (hb->hbmap[i] == NULL)
            gmx_fatal(FARGS,"Could not allocate enough memory for hbmap");
        for (j=0; (j>hb->a.nra); j++)
            hb->hbmap[i][j] = NULL;
    }
}

/* Consider redoing pHist so that is only stores transitions between
 * periodicities and not the periodicity for all frames. This eats heaps of memory. */
static void mk_per(t_hbdata *hb)
{
    int i,j;
    if (hb->bGem) {
        snew(hb->per->pHist, hb->d.nrd);
        for (i=0; i<hb->d.nrd; i++) {
            snew(hb->per->pHist[i], hb->a.nra);
            if (hb->per->pHist[i]==NULL)
                gmx_fatal(FARGS,"Could not allocate enough memory for per->pHist");
            for (j=0; j<hb->a.nra; j++) {
                clearPshift(&(hb->per->pHist[i][j]));
            }
        }
        /* add the [0,0,0] shift to element 0 of p2i. */
        snew(hb->per->p2i, 1);
        clear_ivec(hb->per->p2i[0]);
        hb->per->nper = 1;
    }
}

#ifdef HAVE_NN_LOOPS
static void mk_hbEmap (t_hbdata *hb, int n0)
{
    int i, j, k;
    hb->hbE.E = NULL;
    hb->hbE.nframes = 0;
    snew(hb->hbE.E, hb->d.nrd);
    for (i=0; i<hb->d.nrd; i++)
    {
        snew(hb->hbE.E[i], hb->a.nra);
        for (j=0; j<hb->a.nra; j++)
        {
            snew(hb->hbE.E[i][j], MAXHYDRO);
            for (k=0; k<MAXHYDRO; k++)
                hb->hbE.E[i][j][k] = NULL;
        }
    }
    hb->hbE.Etot = NULL;
}

static void free_hbEmap (t_hbdata *hb)
{
    int i, j, k;
    for (i=0; i<hb->d.nrd; i++)
    {
        for (j=0; j<hb->a.nra; j++)
        {
            for (k=0; k<MAXHYDRO; k++)
                sfree(hb->hbE.E[i][j][k]);
            sfree(hb->hbE.E[i][j]);
        }
        sfree(hb->hbE.E[i]);
    }
    sfree(hb->hbE.E);
    sfree(hb->hbE.Etot);
}

static void addFramesNN(t_hbdata *hb, int frame)
{

#define DELTAFRAMES_HBE 10

    int d,a,h,nframes;

    if (frame >= hb->hbE.nframes) {
        nframes =  hb->hbE.nframes + DELTAFRAMES_HBE;
        srenew(hb->hbE.Etot, nframes);

        for (d=0; d<hb->d.nrd; d++)
            for (a=0; a<hb->a.nra; a++)
                for (h=0; h<hb->d.nhydro[d]; h++)
                    srenew(hb->hbE.E[d][a][h], nframes);
        
        hb->hbE.nframes += DELTAFRAMES_HBE;
    }
}

static t_E calcHbEnergy(int d, int a, int h, rvec x[], t_EEst EEst,
                        matrix box, rvec hbox, t_donors *donors){
    /* d     - donor atom
     * a     - acceptor atom
     * h     - hydrogen
     * alpha - angle between dipoles
     * x[]   - atomic positions
     * EEst  - the type of energy estimate (see enum in hbplugin.h)
     * box   - the box vectors   \
     * hbox  - half box lengths  _These two are only needed for the pbc correction
     */

    t_E E;
    rvec dist;
    rvec dipole[2], xmol[3], xmean[2]; 
    int i;
    real r, realE;

    if (d == a)
        /* Self-interaction */
        return NONSENSE_E;

    switch (EEst)
    {
    case NN_BINARY:
        /* This is a simple binary existence function that sets E=1 whenever
         * the distance between the oxygens is equal too or less than 0.35 nm.
         */
        rvec_sub(x[d], x[a], dist);
        pbc_correct_gem(dist, box, hbox);
        if (norm(dist) <= 0.35)
            E = 1;
        else
            E = 0;
        break;

    case NN_1_over_r3:
        /* Negative potential energy of a dipole.
         * E = -cos(alpha) * 1/r^3 */     
     
        copy_rvec(x[d], xmol[0]); /* donor */
        copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */
        copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */

        svmul(15.9994*(1/1.008), xmol[0], xmean[0]);
        rvec_inc(xmean[0], xmol[1]);
        rvec_inc(xmean[0], xmol[2]);
        for(i=0; i<3; i++)
            xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008;

        /* Assumes that all acceptors are also donors. */
        copy_rvec(x[a], xmol[0]); /* acceptor */
        copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */
        copy_rvec(x[donors->hydro[donors->dptr[a]][1]], xmol[2]); /* hydrogen */


        svmul(15.9994*(1/1.008), xmol[0], xmean[1]);
        rvec_inc(xmean[1], xmol[1]);
        rvec_inc(xmean[1], xmol[2]);
        for(i=0; i<3; i++)
            xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008;

        rvec_sub(xmean[0], xmean[1], dist);
        pbc_correct_gem(dist, box, hbox);
        r = norm(dist);

        realE = pow(r, -3.0);
        E = (t_E)(SCALEFACTOR_E * realE);
        break;

    case NN_dipole:
        /* Negative potential energy of a (unpolarizable) dipole.
         * E = -cos(alpha) * 1/r^3 */
        clear_rvec(dipole[1]);
        clear_rvec(dipole[0]);
     
        copy_rvec(x[d], xmol[0]); /* donor */
        copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */
        copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */

        rvec_inc(dipole[0], xmol[1]);
        rvec_inc(dipole[0], xmol[2]);
        for (i=0; i<3; i++)
            dipole[0][i] *= 0.5;
        rvec_dec(dipole[0], xmol[0]);

        svmul(15.9994*(1/1.008), xmol[0], xmean[0]);
        rvec_inc(xmean[0], xmol[1]);
        rvec_inc(xmean[0], xmol[2]);
        for(i=0; i<3; i++)
            xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008;

        /* Assumes that all acceptors are also donors. */
        copy_rvec(x[a], xmol[0]); /* acceptor */
        copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */
        copy_rvec(x[donors->hydro[donors->dptr[a]][2]], xmol[2]); /* hydrogen */


        rvec_inc(dipole[1], xmol[1]);
        rvec_inc(dipole[1], xmol[2]);
        for (i=0; i<3; i++)
            dipole[1][i] *= 0.5;
        rvec_dec(dipole[1], xmol[0]);

        svmul(15.9994*(1/1.008), xmol[0], xmean[1]);
        rvec_inc(xmean[1], xmol[1]);
        rvec_inc(xmean[1], xmol[2]);
        for(i=0; i<3; i++)
            xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008;

        rvec_sub(xmean[0], xmean[1], dist);
        pbc_correct_gem(dist, box, hbox);
        r = norm(dist);

        double cosalpha = cos_angle(dipole[0],dipole[1]);
        realE = cosalpha * pow(r, -3.0);
        E = (t_E)(SCALEFACTOR_E * realE);
        break;
      
    default:
        printf("Can't do that type of energy estimate: %i\n.", EEst);
        E = NONSENSE_E;
    }

    return E;
}

static void storeHbEnergy(t_hbdata *hb, int d, int a, int h, t_E E, int frame){
    /* hb - hbond data structure
       d  - donor
       a  - acceptor
       h  - hydrogen
       E  - estimate of the energy
       frame - the current frame.
    */

    /* Store the estimated energy */
    if (E == NONSENSE_E)
        E = 0;

    hb->hbE.E[d][a][h][frame] = E;
#ifdef HAVE_OPENMP
#pragma omp critical
#endif
    {
        hb->hbE.Etot[frame] += E;
    }
}
#endif /* HAVE_NN_LOOPS */


/* Finds -v[] in the periodicity index */
static int findMirror(PSTYPE p, ivec v[], PSTYPE nper)
{
    PSTYPE i;
    ivec u;
    for (i=0; i<nper; i++){
        if (v[i][XX] == -(v[p][XX]) &&
            v[i][YY] == -(v[p][YY]) &&
            v[i][ZZ] == -(v[p][ZZ]))
            return (int)i;
    }
    printf("Couldn't find mirror of [%i, %i, %i], index \n",
           v[p][XX],
           v[p][YY],
           v[p][ZZ]);
    return -1;
}
  

static void add_frames(t_hbdata *hb,int nframes)
{
    int  i,j,k,l;
  
    if (nframes >= hb->max_frames) {
        hb->max_frames += 4096;
        srenew(hb->time,hb->max_frames);
        srenew(hb->nhb,hb->max_frames);
        srenew(hb->ndist,hb->max_frames);
        srenew(hb->n_bound,hb->max_frames);
        srenew(hb->nhx,hb->max_frames);
        if (hb->bDAnr)
            srenew(hb->danr,hb->max_frames);
    }
    hb->nframes=nframes;
}

#define OFFSET(frame) (frame / 32)
#define MASK(frame)   (1 << (frame % 32))

static void _set_hb(unsigned int hbexist[],unsigned int frame,gmx_bool bValue)
{
    if (bValue)
        hbexist[OFFSET(frame)] |= MASK(frame);
    else
        hbexist[OFFSET(frame)] &= ~MASK(frame);
}

static gmx_bool is_hb(unsigned int hbexist[],int frame)
{
    return ((hbexist[OFFSET(frame)] & MASK(frame)) != 0) ? 1 : 0;
}

static void set_hb(t_hbdata *hb,int id,int ih, int ia,int frame,int ihb)
{
    unsigned int *ghptr=NULL;
  
    if (ihb == hbHB)
        ghptr = hb->hbmap[id][ia]->h[ih];
    else if (ihb == hbDist)
        ghptr = hb->hbmap[id][ia]->g[ih];
    else
        gmx_fatal(FARGS,"Incomprehensible iValue %d in set_hb",ihb);

    _set_hb(ghptr,frame-hb->hbmap[id][ia]->n0,TRUE);
}

static void addPshift(t_pShift *pHist, PSTYPE p, int frame)
{
    if (pHist->len == 0) {
        snew(pHist->frame, 1);
        snew(pHist->p, 1);
        pHist->len      = 1;
        pHist->frame[0] = frame;
        pHist->p[0]     = p;
        return;
    } else
        if (pHist->p[pHist->len-1] != p) {
            pHist->len++;
            srenew(pHist->frame, pHist->len);
            srenew(pHist->p, pHist->len);
            pHist->frame[pHist->len-1] = frame;
            pHist->p[pHist->len-1]     = p;
        } /* Otherwise, there is no transition. */
    return;
}

static PSTYPE getPshift(t_pShift pHist, int frame)
{
    int f, i;

    if (pHist.len == 0
        || (pHist.len > 0 && pHist.frame[0]>frame))
        return -1;
  
    for (i=0; i<pHist.len; i++)
    {
        f = pHist.frame[i];
        if (f==frame)
            return pHist.p[i];
        if (f>frame)
            return pHist.p[i-1];
    }
  
    /* It seems that frame is after the last periodic transition. Return the last periodicity. */
    return pHist.p[pHist.len-1];
}

static void add_ff(t_hbdata *hbd,int id,int h,int ia,int frame,int ihb, PSTYPE p)
{
    int     i,j,n;
    t_hbond *hb      = hbd->hbmap[id][ia];
    int     maxhydro = min(hbd->maxhydro,hbd->d.nhydro[id]);
    int     wlen     = hbd->wordlen;
    int     delta    = 32*wlen;
    gmx_bool    bGem     = hbd->bGem;

    if (!hb->h[0]) {
        hb->n0        = frame;
        hb->maxframes = delta;
        for(i=0; (i<maxhydro); i++) {
            snew(hb->h[i],hb->maxframes/wlen);
            snew(hb->g[i],hb->maxframes/wlen);
        }
    } else {
        hb->nframes = frame-hb->n0;
        /* We need a while loop here because hbonds may be returning
         * after a long time.
         */
        while (hb->nframes >= hb->maxframes) {
            n = hb->maxframes + delta;
            for(i=0; (i<maxhydro); i++) {
                srenew(hb->h[i],n/wlen);
                srenew(hb->g[i],n/wlen);
                for(j=hb->maxframes/wlen; (j<n/wlen); j++) {
                    hb->h[i][j] = 0;
                    hb->g[i][j] = 0;
                }
            }

            hb->maxframes = n;
        }
    }
    if (frame >= 0) {
        set_hb(hbd,id,h,ia,frame,ihb);
        if (bGem) {
            if (p>=hbd->per->nper)
                gmx_fatal(FARGS, "invalid shift: p=%u, nper=%u", p, hbd->per->nper);
            else
                addPshift(&(hbd->per->pHist[id][ia]), p, frame);
      
        }
    }
}

static void inc_nhbonds(t_donors *ddd,int d, int h)
{
    int j;
    int dptr = ddd->dptr[d];
  
    for(j=0; (j<ddd->nhydro[dptr]); j++)
        if (ddd->hydro[dptr][j] == h) {
            ddd->nhbonds[dptr][j]++;
            break;
        }
    if (j == ddd->nhydro[dptr])
        gmx_fatal(FARGS,"No such hydrogen %d on donor %d\n",h+1,d+1);
}

static int _acceptor_index(t_acceptors *a,int grp,atom_id i,
                           const char *file,int line)
{
    int ai = a->aptr[i];

    if (a->grp[ai] != grp) {
        if (debug && bDebug) 
            fprintf(debug,"Acc. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
                    ai,a->grp[ai],grp,file,line);
        return NOTSET;
    }
    else
        return ai;
}
#define acceptor_index(a,grp,i) _acceptor_index(a,grp,i,__FILE__,__LINE__)

static int _donor_index(t_donors *d,int grp,atom_id i,const char *file,int line)
{
    int di = d->dptr[i];
  
    if (di == NOTSET)
        return NOTSET;

    if (d->grp[di] != grp) {
        if (debug && bDebug)
            fprintf(debug,"Don. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
                    di,d->grp[di],grp,file,line);
        return NOTSET;
    }
    else
        return di;
}
#define donor_index(d,grp,i) _donor_index(d,grp,i,__FILE__,__LINE__)

static gmx_bool isInterchangable(t_hbdata *hb, int d, int a, int grpa, int grpd)
{
    /* g_hbond doesn't allow overlapping groups */
    if (grpa!=grpd)
        return FALSE;
    return
        donor_index(&hb->d,grpd,a) != NOTSET
        && acceptor_index(&hb->a,grpa,d) != NOTSET;
}


static void add_hbond(t_hbdata *hb,int d,int a,int h,int grpd,int grpa,
                      int frame,gmx_bool bMerge,int ihb,gmx_bool bContact, PSTYPE p)
{ 
    int k,id,ia,hh;
    gmx_bool daSwap = FALSE;

    if ((id = hb->d.dptr[d]) == NOTSET)
        gmx_fatal(FARGS,"No donor atom %d",d+1);
    else if (grpd != hb->d.grp[id])
        gmx_fatal(FARGS,"Inconsistent donor groups, %d iso %d, atom %d",
                  grpd,hb->d.grp[id],d+1);
    if ((ia = hb->a.aptr[a]) == NOTSET)
        gmx_fatal(FARGS,"No acceptor atom %d",a+1);
    else if (grpa != hb->a.grp[ia])
        gmx_fatal(FARGS,"Inconsistent acceptor groups, %d iso %d, atom %d",
                  grpa,hb->a.grp[ia],a+1);

    if (bMerge)
        if ((daSwap = isInterchangable(hb, d, a, grpd, grpa) || bContact) && d>a)
            /* Then swap identity so that the id of d is lower then that of a.
             *
             * This should really be redundant by now, as is_hbond() now ought to return
             * hbNo in the cases where this conditional is TRUE. */
        {
            k = d;
            d = a;
            a = k;
        
            /* Now repeat donor/acc check. */
            if ((id = hb->d.dptr[d]) == NOTSET)
                gmx_fatal(FARGS,"No donor atom %d",d+1);
            else if (grpd != hb->d.grp[id])
                gmx_fatal(FARGS,"Inconsistent donor groups, %d iso %d, atom %d",
                          grpd,hb->d.grp[id],d+1);
            if ((ia = hb->a.aptr[a]) == NOTSET)
                gmx_fatal(FARGS,"No acceptor atom %d",a+1);
            else if (grpa != hb->a.grp[ia])
                gmx_fatal(FARGS,"Inconsistent acceptor groups, %d iso %d, atom %d",
                          grpa,hb->a.grp[ia],a+1);
        }

    if (hb->hbmap) {
        /* Loop over hydrogens to find which hydrogen is in this particular HB */
        if ((ihb == hbHB) && !bMerge && !bContact) {
            for(k=0; (k<hb->d.nhydro[id]); k++) 
                if (hb->d.hydro[id][k] == h)
                    break;
            if (k == hb->d.nhydro[id])
                gmx_fatal(FARGS,"Donor %d does not have hydrogen %d (a = %d)",
                          d+1,h+1,a+1);
        }
        else
            k = 0;
    
        if (hb->bHBmap) {
            if (hb->hbmap[id][ia] == NULL) {
                snew(hb->hbmap[id][ia],1);
                snew(hb->hbmap[id][ia]->h,hb->maxhydro);
                snew(hb->hbmap[id][ia]->g,hb->maxhydro);
            }
            add_ff(hb,id,k,ia,frame,ihb,p);
        }
    
        /* Strange construction with frame >=0 is a relic from old code
         * for selected hbond analysis. It may be necessary again if that
         * is made to work again.
         */
        if (frame >= 0) {
            hh = hb->hbmap[id][ia]->history[k];
            if (ihb == hbHB) {
                hb->nhb[frame]++;
                if (!(ISHB(hh))) {
                    hb->hbmap[id][ia]->history[k] = hh | 2;
                    hb->nrhb++;
                }
            }
            else
            {
                if (ihb == hbDist) {
                    hb->ndist[frame]++;
                    if (!(ISDIST(hh))) {
                        hb->hbmap[id][ia]->history[k] = hh | 1;
                        hb->nrdist++;
                    }
                }
            }
        }
    } else {
        if (frame >= 0) {
            if (ihb == hbHB) {
                hb->nhb[frame]++;
            } else {
                if (ihb == hbDist) {
                    hb->ndist[frame]++;
                }
            }
        }
    }
    if (bMerge && daSwap)
        h = hb->d.hydro[id][0];
    /* Increment number if HBonds per H */
    if (ihb == hbHB && !bContact)
        inc_nhbonds(&(hb->d),d,h);
}

/* Now a redundant function. It might find use at some point though. */
static gmx_bool in_list(atom_id selection,int isize,atom_id *index)
{
    int i;
    gmx_bool bFound;
  
    bFound=FALSE;
    for(i=0; (i<isize) && !bFound; i++)
        if(selection == index[i])
            bFound=TRUE;
  
    return bFound;
}

static char *mkatomname(t_atoms *atoms,int i)
{
    static char buf[32];
    int rnr;
  
    rnr = atoms->atom[i].resind;
    sprintf(buf,"%4s%d%-4s",
            *atoms->resinfo[rnr].name,atoms->resinfo[rnr].nr,*atoms->atomname[i]);
  
    return buf;
}

static void gen_datable(atom_id *index, int isize, unsigned char *datable, int natoms){
    /* Generates table of all atoms and sets the ingroup bit for atoms in index[] */
    int i;

    for (i=0; i<isize; i++){
        if (index[i] >= natoms)
            gmx_fatal(FARGS,"Atom has index %d larger than number of atoms %d.",index[i],natoms);
        datable[index[i]] |= INGROUP;
    }
}

static void clear_datable_grp(unsigned char *datable, int size){
    /* Clears group information from the table */
    int i;
    const char mask = !(char)INGROUP;
    if (size > 0)
        for (i=0;i<size;i++)
            datable[i] &= mask;
}

static void add_acc(t_acceptors *a,int ia,int grp)
{
    if (a->nra >= a->max_nra) {
        a->max_nra += 16;
        srenew(a->acc,a->max_nra);
        srenew(a->grp,a->max_nra);
    }
    a->grp[a->nra]   = grp;
    a->acc[a->nra++] = ia;
}

static void search_acceptors(t_topology *top,int isize, 
                             atom_id *index,t_acceptors *a,int grp,
                             gmx_bool bNitAcc,
                             gmx_bool bContact,gmx_bool bDoIt, unsigned char *datable)
{
    int i,n;
  
    if (bDoIt) {
        for (i=0; (i<isize); i++) {
            n = index[i];
            if ((bContact ||
                 (((*top->atoms.atomname[n])[0] == 'O') || 
                  (bNitAcc && ((*top->atoms.atomname[n])[0] == 'N')))) &&
                ISINGRP(datable[n])) {
                datable[n] |= ACC; /* set the atom's acceptor flag in datable. */
                add_acc(a,n,grp);
            }
        }
    }
    snew(a->aptr,top->atoms.nr);
    for(i=0; (i<top->atoms.nr); i++)
        a->aptr[i] = NOTSET;
    for(i=0; (i<a->nra); i++)
        a->aptr[a->acc[i]] = i;
}

static void add_h2d(int id,int ih,t_donors *ddd)
{
    int i;
  
    for(i=0; (i<ddd->nhydro[id]); i++) 
        if (ddd->hydro[id][i] == ih) {
            printf("Hm. This isn't the first time I found this donor (%d,%d)\n",
                   ddd->don[id],ih);
            break;
        }
    if (i == ddd->nhydro[id]) {
        if (ddd->nhydro[id] >= MAXHYDRO)
            gmx_fatal(FARGS,"Donor %d has more than %d hydrogens!",
                      ddd->don[id],MAXHYDRO);
        ddd->hydro[id][i] = ih;
        ddd->nhydro[id]++;
    }
}
  
static void add_dh(t_donors *ddd,int id,int ih,int grp, unsigned char *datable)
{
    int i;

    if (ISDON(datable[id]) || !datable) {
        if (ddd->dptr[id] == NOTSET) { /* New donor */
            i = ddd->nrd;
            ddd->dptr[id] = i;
        } else 
            i = ddd->dptr[id];
  
        if (i == ddd->nrd) {
            if (ddd->nrd >= ddd->max_nrd) {
                ddd->max_nrd += 128;
                srenew(ddd->don,ddd->max_nrd);
                srenew(ddd->nhydro,ddd->max_nrd);
                srenew(ddd->hydro,ddd->max_nrd);
                srenew(ddd->nhbonds,ddd->max_nrd);
                srenew(ddd->grp,ddd->max_nrd);
            }
            ddd->don[ddd->nrd] = id;
            ddd->nhydro[ddd->nrd] = 0;
            ddd->grp[ddd->nrd] = grp;
            ddd->nrd++;
        } else
            ddd->don[i] = id;
        add_h2d(i,ih,ddd);
    } else
        if (datable)
            printf("Warning: Atom %d is not in the d/a-table!\n", id);
}

static void search_donors(t_topology *top, int isize, atom_id *index,
                          t_donors *ddd,int grp,gmx_bool bContact,gmx_bool bDoIt,
                          unsigned char *datable)
{
    int        i,j,nra,n;
    t_functype func_type;
    t_ilist    *interaction;
    atom_id    nr1,nr2;
    gmx_bool       stop;

    if (!ddd->dptr) {
        snew(ddd->dptr,top->atoms.nr);
        for(i=0; (i<top->atoms.nr); i++)
            ddd->dptr[i] = NOTSET;
    }

    if (bContact) {
        if (bDoIt)
            for(i=0; (i<isize); i++) {
                datable[index[i]] |= DON;
                add_dh(ddd,index[i],-1,grp,datable);
            }
    }
    else {
        for(func_type=0; (func_type < F_NRE); func_type++) {
            interaction=&(top->idef.il[func_type]);
            if (func_type == F_POSRES)
            { /* The ilist looks strange for posre. Bug in grompp?
               * We don't need posre interactions for hbonds anyway.*/
                continue;
            }
            for(i=0; i < interaction->nr; 
                i+=interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1) {
                /* next function */
                if (func_type != top->idef.functype[interaction->iatoms[i]]) {
                    fprintf(stderr,"Error in func_type %s",
                            interaction_function[func_type].longname);
                    continue;
                }
        
                /* check out this functype */
                if (func_type == F_SETTLE) {
                    nr1=interaction->iatoms[i+1];
          
                    if (ISINGRP(datable[nr1])) {
                        if (ISINGRP(datable[nr1+1])) {
                            datable[nr1] |= DON;
                            add_dh(ddd,nr1,nr1+1,grp,datable);
                        }
                        if (ISINGRP(datable[nr1+2])) {
                            datable[nr1] |= DON;
                            add_dh(ddd,nr1,nr1+2,grp,datable);
                        }
                    }
                } 
                else if (IS_CHEMBOND(func_type)) {
                    for (j=0; j<2; j++) {
                        nr1=interaction->iatoms[i+1+j];
                        nr2=interaction->iatoms[i+2-j];
                        if ((*top->atoms.atomname[nr1][0] == 'H') && 
                            ((*top->atoms.atomname[nr2][0] == 'O') ||
                             (*top->atoms.atomname[nr2][0] == 'N')) &&
                            ISINGRP(datable[nr1]) && ISINGRP(datable[nr2])) {
                            datable[nr2] |= DON;
                            add_dh(ddd,nr2,nr1,grp,datable);
                        }
                    }
                }
            }
        }
#ifdef SAFEVSITES
        for(func_type=0; func_type < F_NRE; func_type++) {
            interaction=&top->idef.il[func_type];
            for(i=0; i < interaction->nr; 
                i+=interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1) {
                /* next function */
                if (func_type != top->idef.functype[interaction->iatoms[i]])
                    gmx_incons("function type in search_donors");
        
                if ( interaction_function[func_type].flags & IF_VSITE ) {
                    nr1=interaction->iatoms[i+1];
                    if ( *top->atoms.atomname[nr1][0]  == 'H') {
                        nr2=nr1-1;
                        stop=FALSE;
                        while (!stop && ( *top->atoms.atomname[nr2][0] == 'H'))
                            if (nr2)
                                nr2--;
                            else
                                stop=TRUE;
                        if ( !stop && ( ( *top->atoms.atomname[nr2][0] == 'O') ||
                                        ( *top->atoms.atomname[nr2][0] == 'N') ) &&
                             ISINGRP(datable[nr1]) && ISINGRP(datable[nr2])) {
                            datable[nr2] |= DON;
                            add_dh(ddd,nr2,nr1,grp,datable);
                        }
                    }
                }
            }
        }
#endif
    }
}

static t_gridcell ***init_grid(gmx_bool bBox,rvec box[],real rcut,ivec ngrid)
{
    t_gridcell ***grid;
    int i,y,z;
  
    if (bBox)
        for(i=0; i<DIM; i++)
            ngrid[i]=(box[i][i]/(1.2*rcut));
  
    if ( !bBox || (ngrid[XX]<3) || (ngrid[YY]<3) || (ngrid[ZZ]<3) )
        for(i=0; i<DIM; i++)
            ngrid[i]=1;
    else 
        printf("\nWill do grid-seach on %dx%dx%d grid, rcut=%g\n",
               ngrid[XX],ngrid[YY],ngrid[ZZ],rcut);
    snew(grid,ngrid[ZZ]);
    for (z=0; z<ngrid[ZZ]; z++) {
        snew((grid)[z],ngrid[YY]);
        for (y=0; y<ngrid[YY]; y++)
            snew((grid)[z][y],ngrid[XX]);
    }
    return grid;
}

static void control_pHist(t_hbdata *hb, int nframes)
{
    int i,j,k;
    PSTYPE p;
    for (i=0;i<hb->d.nrd;i++)
        for (j=0;j<hb->a.nra;j++)
            if (hb->per->pHist[i][j].len != 0)
                for (k=hb->hbmap[i][j][0].n0; k<nframes; k++) {
                    p = getPshift(hb->per->pHist[i][j], k);
                    if (p>hb->per->nper)
                        fprintf(stderr, "Weird stuff in pHist[%i][%i].p at frame %i: p=%i\n",
                                i,j,k,p);
                }
}

static void reset_nhbonds(t_donors *ddd)
{
    int i,j;
  
    for(i=0; (i<ddd->nrd); i++) 
        for(j=0; (j<MAXHH); j++)
            ddd->nhbonds[i][j] = 0;
}

void pbc_correct_gem(rvec dx,matrix box,rvec hbox);

static void build_grid(t_hbdata *hb,rvec x[], rvec xshell,
                       gmx_bool bBox, matrix box, rvec hbox,
                       real rcut, real rshell,
                       ivec ngrid, t_gridcell ***grid)
{
    int     i,m,gr,xi,yi,zi,nr;
    atom_id *ad;
    ivec    grididx;
    rvec    invdelta,dshell,xtemp={0,0,0};
    t_ncell *newgrid;
    gmx_bool    bDoRshell,bInShell,bAcc;
    real    rshell2=0;
    int     gx,gy,gz;
    int     dum = -1;
  
    bDoRshell = (rshell > 0);
    rshell2   = sqr(rshell);
    bInShell  = TRUE;
  
#define DBB(x) if (debug && bDebug) fprintf(debug,"build_grid, line %d, %s = %d\n",__LINE__,#x,x)
    DBB(dum);
    for(m=0; m<DIM; m++) {
        hbox[m]=box[m][m]*0.5;
        if (bBox) {
            invdelta[m]=ngrid[m]/box[m][m];
            if (1/invdelta[m] < rcut)
                gmx_fatal(FARGS,"Your computational box has shrunk too much.\n"
                          "%s can not handle this situation, sorry.\n",
                          ShortProgram());
        } else
            invdelta[m]=0;
    }
    grididx[XX]=0;
    grididx[YY]=0;
    grididx[ZZ]=0;
    DBB(dum);
    /* resetting atom counts */
    for(gr=0; (gr<grNR); gr++) {
        for (zi=0; zi<ngrid[ZZ]; zi++)
            for (yi=0; yi<ngrid[YY]; yi++)
                for (xi=0; xi<ngrid[XX]; xi++) {
                    grid[zi][yi][xi].d[gr].nr=0;
                    grid[zi][yi][xi].a[gr].nr=0;
                }
        DBB(dum);
    
        /* put atoms in grid cells */
        for(bAcc=FALSE; (bAcc<=TRUE); bAcc++) {
            if (bAcc) {
                nr = hb->a.nra;
                ad = hb->a.acc;
            }
            else {
                nr = hb->d.nrd;
                ad = hb->d.don;
            }
            DBB(bAcc);
            for(i=0; (i<nr); i++) {
                /* check if we are inside the shell */
                /* if bDoRshell=FALSE then bInShell=TRUE always */
                DBB(i);
                if ( bDoRshell ) {
                    bInShell=TRUE;
                    rvec_sub(x[ad[i]],xshell,dshell);
                    if (bBox) {
                        if (FALSE && !hb->bGem) {
                            for(m=DIM-1; m>=0 && bInShell; m--) {
                                if ( dshell[m] < -hbox[m] )
                                    rvec_inc(dshell,box[m]);
                                else if ( dshell[m] >= hbox[m] ) 
                                    dshell[m] -= 2*hbox[m];
                                /* if we're outside the cube, we're outside the sphere also! */
                                if ( (dshell[m]>rshell) || (-dshell[m]>rshell) )
                                    bInShell=FALSE;
                            }
                        } else {
                            gmx_bool bDone = FALSE;
                            while (!bDone)
                            {
                                bDone = TRUE;
                                for(m=DIM-1; m>=0 && bInShell; m--) {
                                    if ( dshell[m] < -hbox[m] ) {
                                        bDone = FALSE;
                                        rvec_inc(dshell,box[m]);
                                    }
                                    if ( dshell[m] >= hbox[m] ) {
                                        bDone = FALSE;
                                        dshell[m] -= 2*hbox[m];
                                    }
                                }
                            }
                            for(m=DIM-1; m>=0 && bInShell; m--) {
                                /* if we're outside the cube, we're outside the sphere also! */
                                if ( (dshell[m]>rshell) || (-dshell[m]>rshell) )
                                    bInShell=FALSE;
                            }
                        }
                    }
                    /* if we're inside the cube, check if we're inside the sphere */
                    if (bInShell)
                        bInShell = norm2(dshell) < rshell2;
                }
                DBB(i);
                if ( bInShell ) {
                    if (bBox) {
                        if (hb->bGem)
                            copy_rvec(x[ad[i]], xtemp);
                        pbc_correct_gem(x[ad[i]], box, hbox);
                    }
                    for(m=DIM-1; m>=0; m--) {
                        if (TRUE || !hb->bGem){
                            /* put atom in the box */
                            while( x[ad[i]][m] < 0 )
                                rvec_inc(x[ad[i]],box[m]);
                            while( x[ad[i]][m] >= box[m][m] )
                                rvec_dec(x[ad[i]],box[m]);
                        }
                        /* determine grid index of atom */
                        grididx[m]=x[ad[i]][m]*invdelta[m];
                        grididx[m] = (grididx[m]+ngrid[m]) % ngrid[m];
                    }
                    if (hb->bGem)
                        copy_rvec(xtemp, x[ad[i]]); /* copy back */
                    gx = grididx[XX];
                    gy = grididx[YY];
                    gz = grididx[ZZ];
                    range_check(gx,0,ngrid[XX]);
                    range_check(gy,0,ngrid[YY]);
                    range_check(gz,0,ngrid[ZZ]);
                    DBB(gx);
                    DBB(gy);
                    DBB(gz);
                    /* add atom to grid cell */
                    if (bAcc)
                        newgrid=&(grid[gz][gy][gx].a[gr]);
                    else
                        newgrid=&(grid[gz][gy][gx].d[gr]);
                    if (newgrid->nr >= newgrid->maxnr) {
                        newgrid->maxnr+=10;
                        DBB(newgrid->maxnr);
                        srenew(newgrid->atoms, newgrid->maxnr);
                    }
                    DBB(newgrid->nr);
                    newgrid->atoms[newgrid->nr]=ad[i];
                    newgrid->nr++;
                }
            }
        }
    }
}

static void count_da_grid(ivec ngrid, t_gridcell ***grid, t_icell danr)
{
    int gr,xi,yi,zi;
  
    for(gr=0; (gr<grNR); gr++) {
        danr[gr]=0;
        for (zi=0; zi<ngrid[ZZ]; zi++)
            for (yi=0; yi<ngrid[YY]; yi++)
                for (xi=0; xi<ngrid[XX]; xi++) {
                    danr[gr] += grid[zi][yi][xi].d[gr].nr;
                }
    }
}

/* The grid loop.
 * Without a box, the grid is 1x1x1, so all loops are 1 long.
 * With a rectangular box (bTric==FALSE) all loops are 3 long.
 * With a triclinic box all loops are 3 long, except when a cell is
 * located next to one of the box edges which is not parallel to the
 * x/y-plane, in that case all cells in a line or layer are searched.
 * This could be implemented slightly more efficient, but the code
 * would get much more complicated.
 */
#define B(n,x,bTric,bEdge) ((n==1) ? x : bTric&&(bEdge) ? 0   : (x-1))
#define E(n,x,bTric,bEdge) ((n==1) ? x : bTric&&(bEdge) ? n-1 : (x+1))
#define GRIDMOD(j,n) (j+n)%(n)
#define LOOPGRIDINNER(x,y,z,xx,yy,zz,xo,yo,zo,n,bTric)                  \
    for(zz=B(n[ZZ],zo,bTric,FALSE); zz<=E(n[ZZ],zo,bTric,FALSE); zz++) { \
    z=GRIDMOD(zz,n[ZZ]);                                                \
    for(yy=B(n[YY],yo,bTric,z==0||z==n[ZZ]-1);                          \
        yy<=E(n[YY],yo,bTric,z==0||z==n[ZZ]-1); yy++) {                 \
    y=GRIDMOD(yy,n[YY]);                                                \
    for(xx=B(n[XX],xo,bTric,y==0||y==n[YY]-1||z==0||z==n[ZZ]-1);                \
        xx<=E(n[XX],xo,bTric,y==0||y==n[YY]-1||z==0||z==n[ZZ]-1); xx++) { \
    x=GRIDMOD(xx,n[XX]);
#define ENDLOOPGRIDINNER                                                \
    }                                                                   \
        }                                                               \
                                        }                                                               \
                                                                        \

static void dump_grid(FILE *fp, ivec ngrid, t_gridcell ***grid)
{
    int gr,x,y,z,sum[grNR];
  
    fprintf(fp,"grid %dx%dx%d\n",ngrid[XX],ngrid[YY],ngrid[ZZ]);
    for (gr=0; gr<grNR; gr++) {
        sum[gr]=0;
        fprintf(fp,"GROUP %d (%s)\n",gr,grpnames[gr]);
        for (z=0; z<ngrid[ZZ]; z+=2) {
            fprintf(fp,"Z=%d,%d\n",z,z+1);
            for (y=0; y<ngrid[YY]; y++) {
                for (x=0; x<ngrid[XX]; x++) {
                    fprintf(fp,"%3d",grid[x][y][z].d[gr].nr);
                    sum[gr]+=grid[z][y][x].d[gr].nr;
                    fprintf(fp,"%3d",grid[x][y][z].a[gr].nr);
                    sum[gr]+=grid[z][y][x].a[gr].nr;
          
                }
                fprintf(fp," | ");
                if ( (z+1) < ngrid[ZZ] )
                    for (x=0; x<ngrid[XX]; x++) {
                        fprintf(fp,"%3d",grid[z+1][y][x].d[gr].nr);
                        sum[gr]+=grid[z+1][y][x].d[gr].nr;
                        fprintf(fp,"%3d",grid[z+1][y][x].a[gr].nr);
                        sum[gr]+=grid[z+1][y][x].a[gr].nr;
                    }
                fprintf(fp,"\n");
            }
        }
    }
    fprintf(fp,"TOTALS:");
    for (gr=0; gr<grNR; gr++)
        fprintf(fp," %d=%d",gr,sum[gr]);
    fprintf(fp,"\n");
}

/* New GMX record! 5 * in a row. Congratulations! 
 * Sorry, only four left.
 */
static void free_grid(ivec ngrid, t_gridcell ****grid)
{
    int y,z;
    t_gridcell ***g = *grid;
  
    for (z=0; z<ngrid[ZZ]; z++) {
        for (y=0; y<ngrid[YY]; y++) {
            sfree(g[z][y]);
        }
        sfree(g[z]);
    }
    sfree(g);
    g=NULL;
}

static void pbc_correct(rvec dx,matrix box,rvec hbox)
{
    int m;
    for(m=DIM-1; m>=0; m--) {
        if ( dx[m] < -hbox[m] )
            rvec_inc(dx,box[m]);
        else if ( dx[m] >= hbox[m] )
            rvec_dec(dx,box[m]);
    }
}

void pbc_correct_gem(rvec dx,matrix box,rvec hbox)
{
    int m;
    gmx_bool bDone = FALSE;
    while (!bDone) {
        bDone = TRUE;
        for(m=DIM-1; m>=0; m--) {
            if ( dx[m] < -hbox[m] ) {
                bDone = FALSE;
                rvec_inc(dx,box[m]);
            }
            if ( dx[m] >= hbox[m] ) {
                bDone = FALSE;
                rvec_dec(dx,box[m]);
            }
        }
    }
}

/* Added argument r2cut, changed contact and implemented 
 * use of second cut-off.
 * - Erik Marklund, June 29, 2006
 */
static int is_hbond(t_hbdata *hb,int grpd,int grpa,int d,int a,
                    real rcut, real r2cut, real ccut, 
                    rvec x[], gmx_bool bBox, matrix box,rvec hbox,
                    real *d_ha, real *ang,gmx_bool bDA,int *hhh,
                    gmx_bool bContact, gmx_bool bMerge, PSTYPE *p)
{
    int  h,hh,id,ja,ihb;
    rvec r_da,r_ha,r_dh, r={0, 0, 0};
    ivec ri;
    real rc2,r2c2,rda2,rha2,ca;
    gmx_bool HAinrange = FALSE; /* If !bDA. Needed for returning hbDist in a correct way. */
    gmx_bool daSwap = FALSE;

    if (d == a)
        return hbNo;

    if (((id = donor_index(&hb->d,grpd,d)) == NOTSET) ||
        ((ja = acceptor_index(&hb->a,grpa,a)) == NOTSET))
        return hbNo;
  
    rc2  = rcut*rcut;
    r2c2 = r2cut*r2cut;
  
    rvec_sub(x[d],x[a],r_da);
    /* Insert projection code here */

    /* if (d>a && ((isInterchangable(hb, d, a, grpd, grpa) && bMerge) || bContact)) */
/*         /\* Then this hbond will be found again, or it has already been found. *\/ */
/*         return hbNo; */

    if (bBox){
        if (d>a && bMerge && (bContact || isInterchangable(hb, d, a, grpd, grpa))) { /* acceptor is also a donor and vice versa? */
            return hbNo;
            daSwap = TRUE; /* If so, then their history should be filed with donor and acceptor swapped. */
        }
        if (hb->bGem) {
            copy_rvec(r_da, r); /* Save this for later */
            pbc_correct_gem(r_da,box,hbox);
        } else {
            pbc_correct_gem(r_da,box,hbox);    
        }
    }
    rda2 = iprod(r_da,r_da);
  
    if (bContact) {
        if (daSwap)
            return hbNo;
        if (rda2 <= rc2){
            if (hb->bGem){
                calcBoxDistance(hb->per->P, r, ri);
                *p = periodicIndex(ri, hb->per, daSwap);        /* find (or add) periodicity index. */
            }
            return hbHB;
        }
        else if (rda2 < r2c2)
            return hbDist;
        else
            return hbNo;
    }
    *hhh = NOTSET;
  
    if (bDA && (rda2 > rc2))
        return hbNo;
  
    for(h=0; (h < hb->d.nhydro[id]); h++) {
        hh = hb->d.hydro[id][h];
        rha2 = rc2+1;
        if (!bDA) {
            rvec_sub(x[hh],x[a],r_ha);
            if (bBox) {
                pbc_correct_gem(r_ha,box,hbox);
            }
            rha2 = iprod(r_ha,r_ha);
        }

        if (hb->bGem) {
            calcBoxDistance(hb->per->P, r, ri);
            *p = periodicIndex(ri, hb->per, daSwap);    /* find periodicity index. */
        }

        if (bDA || (!bDA && (rha2 <= rc2))) {
            rvec_sub(x[d],x[hh],r_dh);
            if (bBox) {
                if (hb->bGem)
                    pbc_correct_gem(r_dh,box,hbox);
                else
                    pbc_correct_gem(r_dh,box,hbox);
            }
        
            if (!bDA)
                HAinrange = TRUE;
            ca = cos_angle(r_dh,r_da);
            /* if angle is smaller, cos is larger */
            if (ca >= ccut) {
                *hhh  = hh;
                *d_ha = sqrt(bDA?rda2:rha2);
                *ang  = acos(ca);
                return hbHB;
            }
        }
    }
    if (bDA || (!bDA && HAinrange))
        return hbDist;
    else
        return hbNo;
}

/* Fixed previously undiscovered bug in the merge
   code, where the last frame of each hbond disappears.
   - Erik Marklund, June 1, 2006 */
/* Added the following arguments:
 *   ptmp[] - temporary periodicity hisory
 *   a1     - identity of first acceptor/donor
 *   a2     - identity of second acceptor/donor
 * - Erik Marklund, FEB 20 2010 */

/* Merging is now done on the fly, so do_merge is most likely obsolete now.
 * Will do some more testing before removing the function entirely.
 * - Erik Marklund, MAY 10 2010 */
static void do_merge(t_hbdata *hb,int ntmp,
                     unsigned int htmp[],unsigned int gtmp[],PSTYPE ptmp[],
                     t_hbond *hb0,t_hbond *hb1, int a1, int a2)
{
    /* Here we need to make sure we're treating periodicity in
     * the right way for the geminate recombination kinetics. */

    int m,mm,n00,n01,nn0,nnframes;
    PSTYPE pm;
    t_pShift *pShift;

    /* Decide where to start from when merging */
    n00      = hb0->n0;
    n01      = hb1->n0;
    nn0      = min(n00,n01);
    nnframes = max(n00 + hb0->nframes,n01 + hb1->nframes) - nn0;
    /* Initiate tmp arrays */
    for(m=0; (m<ntmp); m++) {
        htmp[m] = 0;
        gtmp[m] = 0;
        ptmp[m] = 0;
    }
    /* Fill tmp arrays with values due to first HB */
    /* Once again '<' had to be replaced with '<='
       to catch the last frame in which the hbond
       appears.
       - Erik Marklund, June 1, 2006 */  
    for(m=0; (m<=hb0->nframes); m++) {
        mm = m+n00-nn0;
        htmp[mm] = is_hb(hb0->h[0],m);
        if (hb->bGem) {
            pm = getPshift(hb->per->pHist[a1][a2], m+hb0->n0);
            if (pm > hb->per->nper)
                gmx_fatal(FARGS, "Illegal shift!");
            else
                ptmp[mm] = pm; /*hb->per->pHist[a1][a2][m];*/
        }
    }
    /* If we're doing geminate recompbination we usually don't need the distances.
     * Let's save some memory and time. */
    if (TRUE || !hb->bGem || hb->per->gemtype == gemAD){
        for(m=0; (m<=hb0->nframes); m++) {
            mm = m+n00-nn0;
            gtmp[mm] = is_hb(hb0->g[0],m);
        }
    }
    /* Next HB */
    for(m=0; (m<=hb1->nframes); m++) {
        mm = m+n01-nn0;
        htmp[mm] = htmp[mm] || is_hb(hb1->h[0],m);
        gtmp[mm] = gtmp[mm] || is_hb(hb1->g[0],m);
        if (hb->bGem /* && ptmp[mm] != 0 */) {

            /* If this hbond has been seen before with donor and acceptor swapped,
             * then we need to find the mirrored (*-1) periodicity vector to truely
             * merge the hbond history. */
            pm = findMirror(getPshift(hb->per->pHist[a2][a1],m+hb1->n0), hb->per->p2i, hb->per->nper);
            /* Store index of mirror */
            if (pm > hb->per->nper)
                gmx_fatal(FARGS, "Illegal shift!");
            ptmp[mm] = pm;
        }
    }
    /* Reallocate target array */
    if (nnframes > hb0->maxframes) {
        srenew(hb0->h[0],4+nnframes/hb->wordlen);
        srenew(hb0->g[0],4+nnframes/hb->wordlen);  
    }
    if (NULL != hb->per->pHist)
    {
        clearPshift(&(hb->per->pHist[a1][a2]));
    }

    /* Copy temp array to target array */
    for(m=0; (m<=nnframes); m++) {
        _set_hb(hb0->h[0],m,htmp[m]);
        _set_hb(hb0->g[0],m,gtmp[m]);
        if (hb->bGem)
            addPshift(&(hb->per->pHist[a1][a2]), ptmp[m], m+nn0);
    }
  
    /* Set scalar variables */
    hb0->n0        = nn0;
    hb0->maxframes = nnframes;
}

/* Added argument bContact for nicer output.
 * Erik Marklund, June 29, 2006
 */
static void merge_hb(t_hbdata *hb,gmx_bool bTwo, gmx_bool bContact){
    int  i,inrnew,indnew,j,ii,jj,m,id,ia,grp,ogrp,ntmp;
    unsigned int *htmp,*gtmp;
    PSTYPE *ptmp;
    t_hbond *hb0,*hb1;

    inrnew = hb->nrhb;
    indnew = hb->nrdist;
    
    /* Check whether donors are also acceptors */
    printf("Merging hbonds with Acceptor and Donor swapped\n");

    ntmp = 2*hb->max_frames;
    snew(gtmp,ntmp);
    snew(htmp,ntmp);
    snew(ptmp,ntmp);
    for(i=0; (i<hb->d.nrd); i++) {
        fprintf(stderr,"\r%d/%d",i+1,hb->d.nrd);
        id = hb->d.don[i];
        ii = hb->a.aptr[id];
        for(j=0; (j<hb->a.nra); j++) {
            ia = hb->a.acc[j];
            jj = hb->d.dptr[ia];
            if ((id != ia) && (ii != NOTSET) && (jj != NOTSET) &&
                (!bTwo || (bTwo && (hb->d.grp[i] != hb->a.grp[j])))) {
                hb0 = hb->hbmap[i][j];
                hb1 = hb->hbmap[jj][ii];
                if (hb0 && hb1 && ISHB(hb0->history[0]) && ISHB(hb1->history[0])) {
                    do_merge(hb,ntmp,htmp,gtmp,ptmp,hb0,hb1,i,j);
                    if (ISHB(hb1->history[0])) 
                        inrnew--;
                    else if (ISDIST(hb1->history[0])) 
                        indnew--;
                    else
                        if (bContact) 
                            gmx_incons("No contact history");
                        else
                            gmx_incons("Neither hydrogen bond nor distance");
                    sfree(hb1->h[0]);
                    sfree(hb1->g[0]);
                    if (hb->bGem) {
                        clearPshift(&(hb->per->pHist[jj][ii]));
                    }
                    hb1->h[0] = NULL;
                    hb1->g[0] = NULL;
                    hb1->history[0] = hbNo;
                }
            }
        }
    }
    fprintf(stderr,"\n");
    printf("- Reduced number of hbonds from %d to %d\n",hb->nrhb,inrnew);
    printf("- Reduced number of distances from %d to %d\n",hb->nrdist,indnew);
    hb->nrhb   = inrnew;
    hb->nrdist = indnew;
    sfree(gtmp);
    sfree(htmp);
    sfree(ptmp);
}

static void do_nhb_dist(FILE *fp,t_hbdata *hb,real t) 
{
    int  i,j,k,n_bound[MAXHH],nbtot;
    h_id nhb;

  
    /* Set array to 0 */
    for(k=0; (k<MAXHH); k++)
        n_bound[k] = 0;
    /* Loop over possible donors */
    for(i=0; (i<hb->d.nrd); i++) {
        for(j=0; (j<hb->d.nhydro[i]); j++)
            n_bound[hb->d.nhbonds[i][j]]++;
    }      
    fprintf(fp,"%12.5e",t);
    nbtot = 0;
    for(k=0; (k<MAXHH); k++) {
        fprintf(fp,"  %8d",n_bound[k]);
        nbtot += n_bound[k]*k;
    }
    fprintf(fp,"  %8d\n",nbtot);
}

/* Added argument bContact in do_hblife(...). Also
 * added support for -contact in function body.
 * - Erik Marklund, May 31, 2006 */
/* Changed the contact code slightly.
 * - Erik Marklund, June 29, 2006
 */
static void do_hblife(const char *fn,t_hbdata *hb,gmx_bool bMerge,gmx_bool bContact,
                      const output_env_t oenv)
{
    FILE *fp;
    const char *leg[] = { "p(t)", "t p(t)" };
    int  *histo;
    int  i,j,j0,k,m,nh,ihb,ohb,nhydro,ndump=0;
    int   nframes = hb->nframes;
    unsigned int **h;
    real   t,x1,dt;
    double sum,integral;
    t_hbond *hbh;
  
    snew(h,hb->maxhydro);
    snew(histo,nframes+1);
    /* Total number of hbonds analyzed here */
    for(i=0; (i<hb->d.nrd); i++) {
        for(k=0; (k<hb->a.nra); k++) {
            hbh = hb->hbmap[i][k];
            if (hbh) {
                if (bMerge) {
                    if (hbh->h[0]) {
                        h[0] = hbh->h[0];
                        nhydro = 1;
                    }
                    else
                        nhydro = 0;
                }
                else {
                    nhydro = 0;
                    for(m=0; (m<hb->maxhydro); m++)
                        if (hbh->h[m]) {
                            h[nhydro++] = bContact ? hbh->g[m] : hbh->h[m];
                        }
                }
                for(nh=0; (nh<nhydro); nh++) {
                    ohb = 0;
                    j0  = 0;

                    /* Changed '<' into '<=' below, just like I
                       did in the hbm-output-loop in the main code.
                       - Erik Marklund, May 31, 2006
                    */
                    for(j=0; (j<=hbh->nframes); j++) {
                        ihb      = is_hb(h[nh],j);
                        if (debug && (ndump < 10))
                            fprintf(debug,"%5d  %5d\n",j,ihb);
                        if (ihb != ohb) {
                            if (ihb) {
                                j0 = j;
                            }
                            else {
                                histo[j-j0]++;
                            }
                            ohb = ihb;
                        }
                    }
                    ndump++;
                }
            }
        }
    }
    fprintf(stderr,"\n");
    if (bContact)
        fp = xvgropen(fn,"Uninterrupted contact lifetime",output_env_get_xvgr_tlabel(oenv),"()",oenv);
    else
        fp = xvgropen(fn,"Uninterrupted hydrogen bond lifetime",output_env_get_xvgr_tlabel(oenv),"()",
                      oenv);

    xvgr_legend(fp,asize(leg),leg,oenv);
    j0 = nframes-1;
    while ((j0 > 0) && (histo[j0] == 0))
        j0--;
    sum = 0;
    for(i=0; (i<=j0); i++)
        sum+=histo[i];
    dt       = hb->time[1]-hb->time[0];
    sum      = dt*sum;
    integral = 0;
    for(i=1; (i<=j0); i++) {
        t  = hb->time[i] - hb->time[0] - 0.5*dt;
        x1 = t*histo[i]/sum;
        fprintf(fp,"%8.3f  %10.3e  %10.3e\n",t,histo[i]/sum,x1);
        integral += x1;
    }
    integral *= dt;
    ffclose(fp);
    printf("%s lifetime = %.2f ps\n", bContact?"Contact":"HB", integral);
    printf("Note that the lifetime obtained in this manner is close to useless\n");
    printf("Use the -ac option instead and check the Forward lifetime\n");
    please_cite(stdout,"Spoel2006b");
    sfree(h);
    sfree(histo);
}

/* Changed argument bMerge into oneHB to handle contacts properly.
 * - Erik Marklund, June 29, 2006
 */
static void dump_ac(t_hbdata *hb,gmx_bool oneHB,int nDump)
{
    FILE  *fp;
    int   i,j,k,m,nd,ihb,idist;
    int   nframes = hb->nframes;
    gmx_bool  bPrint;
    t_hbond *hbh;

    if (nDump <= 0)
        return;
    fp = ffopen("debug-ac.xvg","w");
    for(j=0; (j<nframes); j++) {
        fprintf(fp,"%10.3f",hb->time[j]);
        for(i=nd=0; (i<hb->d.nrd) && (nd < nDump); i++) {
            for(k=0; (k<hb->a.nra) && (nd < nDump); k++) {
                bPrint = FALSE;
                ihb = idist = 0;
                hbh = hb->hbmap[i][k];
                if (oneHB) {
                    if (hbh->h[0]) {
                        ihb   = is_hb(hbh->h[0],j);
                        idist = is_hb(hbh->g[0],j);
                        bPrint = TRUE;
                    }
                } 
                else {
                    for(m=0; (m<hb->maxhydro) && !ihb ; m++) {
                        ihb   = ihb   || ((hbh->h[m]) && is_hb(hbh->h[m],j));
                        idist = idist || ((hbh->g[m]) && is_hb(hbh->g[m],j));
                    }
                    /* This is not correct! */
                    /* What isn't correct? -Erik M */
                    bPrint = TRUE;
                }
                if (bPrint) {
                    fprintf(fp,"  %1d-%1d",ihb,idist);
                    nd++;
                }
            }
        }
        fprintf(fp,"\n");
    }
    ffclose(fp);
}

static real calc_dg(real tau,real temp)
{
    real kbt;
  
    kbt = BOLTZ*temp;
    if (tau <= 0)
        return -666;
    else
        return kbt*log(kbt*tau/PLANCK);  
}

typedef struct {
    int  n0,n1,nparams,ndelta;
    real kkk[2];
    real *t,*ct,*nt,*kt,*sigma_ct,*sigma_nt,*sigma_kt;
} t_luzar;

#ifdef HAVE_LIBGSL
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_version.h>

static double my_f(const gsl_vector *v,void *params)
{
    t_luzar *tl = (t_luzar *)params;
    int    i;
    double tol=1e-16,chi2=0;
    double di;
    real   k,kp;
  
    for(i=0; (i<tl->nparams); i++) {
        tl->kkk[i] = gsl_vector_get(v, i);
    }
    k  = tl->kkk[0];
    kp = tl->kkk[1];
  
    for(i=tl->n0; (i<tl->n1); i+=tl->ndelta) {
        di=sqr(k*tl->sigma_ct[i]) + sqr(kp*tl->sigma_nt[i]) + sqr(tl->sigma_kt[i]);
        /*di = 1;*/
        if (di > tol)
            chi2 += sqr(k*tl->ct[i]-kp*tl->nt[i]-tl->kt[i])/di;
      
        else
            fprintf(stderr,"WARNING: sigma_ct = %g, sigma_nt = %g, sigma_kt = %g\n"
                    "di = %g k = %g kp = %g\n",tl->sigma_ct[i],
                    tl->sigma_nt[i],tl->sigma_kt[i],di,k,kp);
    }
#ifdef DEBUG
    chi2 = 0.3*sqr(k-0.6)+0.7*sqr(kp-1.3);
#endif
    return chi2;
}

static real optimize_luzar_parameters(FILE *fp,t_luzar *tl,int maxiter,
                                      real tol)
{
    real   size,d2;
    int    iter   = 0;
    int    status = 0;
    int    i;

    const gsl_multimin_fminimizer_type *T;
    gsl_multimin_fminimizer *s;

    gsl_vector *x,*dx;
    gsl_multimin_function my_func;

    my_func.f      = &my_f;
    my_func.n      = tl->nparams;
    my_func.params = (void *) tl;

    /* Starting point */
    x = gsl_vector_alloc (my_func.n);
    for(i=0; (i<my_func.n); i++)
        gsl_vector_set (x, i, tl->kkk[i]);
  
    /* Step size, different for each of the parameters */
    dx = gsl_vector_alloc (my_func.n);
    for(i=0; (i<my_func.n); i++)
        gsl_vector_set (dx, i, 0.01*tl->kkk[i]);

    T = gsl_multimin_fminimizer_nmsimplex;
    s = gsl_multimin_fminimizer_alloc (T, my_func.n);

    gsl_multimin_fminimizer_set (s, &my_func, x, dx);
    gsl_vector_free (x);
    gsl_vector_free (dx);

    if (fp)
        fprintf(fp,"%5s %12s %12s %12s %12s\n","Iter","k","kp","NM Size","Chi2");
  
    do  {
        iter++;
        status = gsl_multimin_fminimizer_iterate (s);
    
        if (status != 0)
            gmx_fatal(FARGS,"Something went wrong in the iteration in minimizer %s",
                      gsl_multimin_fminimizer_name(s));
    
        d2     = gsl_multimin_fminimizer_minimum(s);
        size   = gsl_multimin_fminimizer_size(s);
        status = gsl_multimin_test_size(size,tol);
    
        if (status == GSL_SUCCESS)
            if (fp) 
                fprintf(fp,"Minimum found using %s at:\n",
                        gsl_multimin_fminimizer_name(s));
  
        if (fp) {
            fprintf(fp,"%5d", iter);
            for(i=0; (i<my_func.n); i++) 
                fprintf(fp," %12.4e",gsl_vector_get (s->x,i));
            fprintf (fp," %12.4e %12.4e\n",size,d2);
        }
    }
    while ((status == GSL_CONTINUE) && (iter < maxiter));
  
    gsl_multimin_fminimizer_free (s);
  
    return d2;
}

static real quality_of_fit(real chi2,int N)
{
    return gsl_sf_gamma_inc_Q((N-2)/2.0,chi2/2.0);
}

#else
static real optimize_luzar_parameters(FILE *fp,t_luzar *tl,int maxiter,
                                      real tol)
{
    fprintf(stderr,"This program needs the GNU scientific library to work.\n");
  
    return -1;
}
static real quality_of_fit(real chi2,int N)
{
    fprintf(stderr,"This program needs the GNU scientific library to work.\n");
  
    return -1;
}

#endif

static real compute_weighted_rates(int n,real t[],real ct[],real nt[],
                                   real kt[],real sigma_ct[],real sigma_nt[],
                                   real sigma_kt[],real *k,real *kp,
                                   real *sigma_k,real *sigma_kp,
                                   real fit_start)
{
#define NK 10
    int      i,j;
    t_luzar  tl;
    real     kkk=0,kkp=0,kk2=0,kp2=0,chi2;
  
    *sigma_k  = 0;
    *sigma_kp = 0;
  
    for(i=0; (i<n); i++) {
        if (t[i] >= fit_start)
            break;
    }
    tl.n0      = i;
    tl.n1      = n;
    tl.nparams = 2;
    tl.ndelta  = 1;
    tl.t  = t;
    tl.ct = ct;
    tl.nt = nt;
    tl.kt = kt;
    tl.sigma_ct = sigma_ct;
    tl.sigma_nt = sigma_nt;
    tl.sigma_kt = sigma_kt;
    tl.kkk[0] = *k;
    tl.kkk[1] = *kp;
  
    chi2 = optimize_luzar_parameters(debug,&tl,1000,1e-3);
    *k  = tl.kkk[0];
    *kp = tl.kkk[1] = *kp;
    tl.ndelta = NK;
    for(j=0; (j<NK); j++) {
        (void) optimize_luzar_parameters(debug,&tl,1000,1e-3);
        kkk += tl.kkk[0];
        kkp += tl.kkk[1];
        kk2 += sqr(tl.kkk[0]);
        kp2 += sqr(tl.kkk[1]);
        tl.n0++;
    }
    *sigma_k  = sqrt(kk2/NK - sqr(kkk/NK));
    *sigma_kp = sqrt(kp2/NK - sqr(kkp/NK));
  
    return chi2;
}

static void smooth_tail(int n,real t[],real c[],real sigma_c[],real start,
                        const output_env_t oenv)
{
    FILE *fp;
    real e_1,fitparm[4];
    int  i;
  
    e_1 = exp(-1);
    for(i=0; (i<n); i++)
        if (c[i] < e_1)
            break;
    if (i < n)
        fitparm[0] = t[i];
    else
        fitparm[0] = 10;
    fitparm[1] = 0.95;
    do_lmfit(n,c,sigma_c,0,t,start,t[n-1],oenv,bDebugMode(),effnEXP2,fitparm,0);
}

void analyse_corr(int n,real t[],real ct[],real nt[],real kt[],
                  real sigma_ct[],real sigma_nt[],real sigma_kt[],
                  real fit_start,real temp,real smooth_tail_start,
                  const output_env_t oenv)
{
    int    i0,i;
    real   k=1,kp=1,kow=1;
    real   Q=0,chi22,chi2,dg,dgp,tau_hb,dtau,tau_rlx,e_1,dt,sigma_k,sigma_kp,ddg;
    double tmp,sn2=0,sc2=0,sk2=0,scn=0,sck=0,snk=0;
    gmx_bool   bError = (sigma_ct != NULL) && (sigma_nt != NULL) && (sigma_kt != NULL);
  
    if (smooth_tail_start >= 0) {
        smooth_tail(n,t,ct,sigma_ct,smooth_tail_start,oenv);
        smooth_tail(n,t,nt,sigma_nt,smooth_tail_start,oenv);
        smooth_tail(n,t,kt,sigma_kt,smooth_tail_start,oenv);
    }
    for(i0=0; (i0<n-2) && ((t[i0]-t[0]) < fit_start); i0++)
        ;
    if (i0 < n-2) { 
        for(i=i0; (i<n); i++) {
            sc2 += sqr(ct[i]);
            sn2 += sqr(nt[i]);
            sk2 += sqr(kt[i]);
            sck += ct[i]*kt[i];
            snk += nt[i]*kt[i];
            scn += ct[i]*nt[i];
        }
        printf("Hydrogen bond thermodynamics at T = %g K\n",temp);
        tmp = (sn2*sc2-sqr(scn));
        if ((tmp > 0) && (sn2 > 0)) {
            k    = (sn2*sck-scn*snk)/tmp;
            kp   = (k*scn-snk)/sn2;
            if (bError) {
                chi2 = 0;
                for(i=i0; (i<n); i++) {
                    chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
                }
                chi22 = compute_weighted_rates(n,t,ct,nt,kt,sigma_ct,sigma_nt,
                                               sigma_kt,&k,&kp,
                                               &sigma_k,&sigma_kp,fit_start);
                Q = quality_of_fit(chi2,2);
                ddg = BOLTZ*temp*sigma_k/k;
                printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n",
                       chi2,Q);
                printf("The Rate and Delta G are followed by an error estimate\n");
                printf("----------------------------------------------------------\n"
                       "Type      Rate (1/ps)  Sigma Time (ps)  DG (kJ/mol)  Sigma\n");
                printf("Forward    %10.3f %6.2f   %8.3f  %10.3f %6.2f\n",
                       k,sigma_k,1/k,calc_dg(1/k,temp),ddg);
                ddg = BOLTZ*temp*sigma_kp/kp;
                printf("Backward   %10.3f %6.2f   %8.3f  %10.3f %6.2f\n",
                       kp,sigma_kp,1/kp,calc_dg(1/kp,temp),ddg);
            }
            else {
                chi2 = 0;
                for(i=i0; (i<n); i++) {
                    chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
                }
                printf("Fitting parameters chi^2 = %10g\nQ = %10g\n",
                       chi2,Q);
                printf("--------------------------------------------------\n"
                       "Type      Rate (1/ps) Time (ps)  DG (kJ/mol)  Chi^2\n");
                printf("Forward    %10.3f   %8.3f  %10.3f  %10g\n",
                       k,1/k,calc_dg(1/k,temp),chi2);
                printf("Backward   %10.3f   %8.3f  %10.3f\n",
                       kp,1/kp,calc_dg(1/kp,temp));
            }
        }
        if (sc2 > 0) {
            kow  = 2*sck/sc2;
            printf("One-way    %10.3f   %s%8.3f  %10.3f\n",
                   kow,bError ? "       " : "",1/kow,calc_dg(1/kow,temp));
        }
        else 
            printf(" - Numerical problems computing HB thermodynamics:\n"
                   "sc2 = %g  sn2 = %g  sk2 = %g sck = %g snk = %g scn = %g\n",
                   sc2,sn2,sk2,sck,snk,scn);
        /* Determine integral of the correlation function */
        tau_hb = evaluate_integral(n,t,ct,NULL,(t[n-1]-t[0])/2,&dtau);
        printf("Integral   %10.3f   %s%8.3f  %10.3f\n",1/tau_hb,
               bError ? "       " : "",tau_hb,calc_dg(tau_hb,temp));
        e_1 = exp(-1);
        for(i=0; (i<n-2); i++) {
            if ((ct[i] > e_1) && (ct[i+1] <= e_1)) {
                break;
            }
        }
        if (i < n-2) {
            /* Determine tau_relax from linear interpolation */
            tau_rlx = t[i]-t[0] + (e_1-ct[i])*(t[i+1]-t[i])/(ct[i+1]-ct[i]);
            printf("Relaxation %10.3f   %8.3f  %s%10.3f\n",1/tau_rlx,
                   tau_rlx,bError ? "       " : "",
                   calc_dg(tau_rlx,temp));
        }
    }
    else 
        printf("Correlation functions too short to compute thermodynamics\n");
}

void compute_derivative(int nn,real x[],real y[],real dydx[])
{
    int j;
  
    /* Compute k(t) = dc(t)/dt */
    for(j=1; (j<nn-1); j++)
        dydx[j] = (y[j+1]-y[j-1])/(x[j+1]-x[j-1]);
    /* Extrapolate endpoints */
    dydx[0]    = 2*dydx[1]   -  dydx[2];
    dydx[nn-1] = 2*dydx[nn-2] - dydx[nn-3];
}

static void parallel_print(int *data, int nThreads)
{
    /* This prints the donors on which each tread is currently working. */
    int i;

    fprintf(stderr, "\r");
    for (i=0; i<nThreads; i++)
        fprintf(stderr, "%-7i",data[i]);
}

static void normalizeACF(real *ct, real *gt, int len)
{
    real ct_fac, gt_fac;
    int i;

    /* Xu and Berne use the same normalization constant */

    ct_fac = 1.0/ct[0];
    gt_fac = (gt!=NULL && gt[0]!=0) ? 1.0/gt[0] : 0;
    printf("Normalization for c(t) = %g for gh(t) = %g\n",ct_fac,gt_fac);
    for (i=0; i<len; i++)
    {
        ct[i] *= ct_fac;
        if (gt != NULL)
            gt[i] *= gt_fac;
    }
}

/* Added argument bContact in do_hbac(...). Also
 * added support for -contact in the actual code.
 * - Erik Marklund, May 31, 2006 */
/* Changed contact code and added argument R2 
 * - Erik Marklund, June 29, 2006
 */
static void do_hbac(const char *fn,t_hbdata *hb,
                    int nDump,gmx_bool bMerge,gmx_bool bContact, real fit_start,
                    real temp,gmx_bool R2,real smooth_tail_start, const output_env_t oenv,
                    t_gemParams *params, const char *gemType, int nThreads,
                    const int NN, const gmx_bool bBallistic, const gmx_bool bGemFit)
{
    FILE *fp;
    int  i,j,k,m,n,o,nd,ihb,idist,n2,nn,iter,nSets;
    const char *legNN[]   = { "Ac(t)",
                               "Ac'(t)"};
    static char **legGem;
                              
    const char *legLuzar[] = { "Ac\\sfin sys\\v{}\\z{}(t)",
                                "Ac(t)",
                                "Cc\\scontact,hb\\v{}\\z{}(t)",
                                "-dAc\\sfs\\v{}\\z{}/dt" };
    gmx_bool bNorm=FALSE;
    double nhb = 0;
    int nhbi=0;
    real *rhbex=NULL,*ht,*gt,*ght,*dght,*kt;
    real *ct,*p_ct,tail,tail2,dtail,ct_fac,ght_fac,*cct;
    const real tol = 1e-3;
    int   nframes = hb->nframes,nf;
    unsigned int **h,**g;
    int   nh,nhbonds,nhydro,ngh;
    t_hbond *hbh;
    PSTYPE p, *pfound = NULL, np;
    t_pShift *pHist;
    int *ptimes=NULL, *poff=NULL, anhb, n0, mMax=INT_MIN;
    real **rHbExGem = NULL;
    gmx_bool c;
    int acType;
    t_E *E;
    double *ctdouble, *timedouble, *fittedct;
    double fittolerance=0.1;

    enum {AC_NONE, AC_NN, AC_GEM, AC_LUZAR};


#ifdef HAVE_OPENMP
    int *dondata=NULL, thisThread;
#endif


    printf("Doing autocorrelation ");

    /* Decide what kind of ACF calculations to do. */
    if (NN > NN_NONE && NN < NN_NR) {
#ifdef HAVE_NN_LOOPS
        acType = AC_NN;
        printf("using the energy estimate.\n");
#else
        acType = AC_NONE;
        printf("Can't do the NN-loop. Yet.\n");
#endif
    } else if (hb->bGem) {
        acType = AC_GEM;
        printf("according to the reversible geminate recombination model by Omer Markowitch.\n");

        nSets = 1 + (bBallistic ? 1:0) + (bGemFit ? 1:0);
        snew(legGem, nSets);
        for (i=0;i<nSets;i++)
            snew(legGem[i], 128);
        sprintf(legGem[0], "Ac\\s%s\\v{}\\z{}(t)", gemType);
        if (bBallistic)
            sprintf(legGem[1], "Ac'(t)");
        if (bGemFit)
            sprintf(legGem[(bBallistic ? 3:2)], "Ac\\s%s,fit\\v{}\\z{}(t)", gemType);

    } else {
        acType = AC_LUZAR;
        printf("according to the theory of Luzar and Chandler.\n");
    }
    fflush(stdout);

    /* build hbexist matrix in reals for autocorr */
    /* Allocate memory for computing ACF (rhbex) and aggregating the ACF (ct) */
    n2=1;
    while (n2 < nframes)
        n2*=2;
  
    nn = nframes/2;
  
    if (acType != AC_NN ||
#ifndef HAVE_OPENMP
        TRUE
#else
        FALSE
#endif
        ) {
        snew(h,hb->maxhydro);
        snew(g,hb->maxhydro);
    }

    /* Dump hbonds for debugging */
    dump_ac(hb,bMerge||bContact,nDump);
  
    /* Total number of hbonds analyzed here */
    nhbonds = 0;
    ngh     = 0;
    anhb    = 0;

    /* ------------------------------------------------
     * I got tired of waiting for the acf calculations
     * and parallelized it with openMP
     * set environment variable CFLAGS = "-fopenmp" when running
     * configure and define DOUSEOPENMP to make use of it.
     */

#ifdef HAVE_OPENMP  /* ================================================= \
                     * Set up the OpenMP stuff,                           |
                     * like the number of threads and such                |
                     */
    if (acType != AC_LUZAR)
    {
#if (_OPENMP >= 200805) /* =====================\ */
        nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_thread_limit());
#else
        nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_num_procs());
#endif /* _OPENMP >= 200805 ====================/ */

        omp_set_num_threads(nThreads);
        snew(dondata, nThreads);
        for (i=0; i<nThreads; i++)
            dondata[i] = -1;
        printf("ACF calculations parallelized with OpenMP using %i threads.\n"
               "Expect close to linear scaling over this donor-loop.\n", nThreads);
        fflush(stdout);
        fprintf(stderr, "Donors: [thread no]\n");
        {
            char tmpstr[7];
            for (i=0;i<nThreads;i++) {
                snprintf(tmpstr, 7, "[%i]", i);
                fprintf(stderr, "%-7s", tmpstr);
            }
        }
        fprintf(stderr, "\n"); /*                                         | */
    }  /*                                                                 | */
#endif /* HAVE_OPENMP ===================================================/  */


    /* Build the ACF according to acType */
    switch (acType)
    {
      
    case AC_NN:
#ifdef HAVE_NN_LOOPS
        /* Here we're using the estimated energy for the hydrogen bonds. */
        snew(ct,nn);
#ifdef HAVE_OPENMP /* ==================================\ */      
#pragma omp parallel                            \
    private(i, j, k, nh, E, rhbex, thisThread), \
    default(shared)
        {
#pragma omp barrier
            thisThread = omp_get_thread_num();
            rhbex = NULL;
#endif /* ==============================================/ */

            snew(rhbex, n2);
            memset(rhbex, 0, n2*sizeof(real)); /* Trust no-one, not even malloc()! */

#ifdef HAVE_OPENMP /* ################################################## \
                    *                                                    #
                    *                                                    #
                    */
#pragma omp barrier
#pragma omp for schedule (dynamic)
#endif
            for (i=0; i<hb->d.nrd; i++) /* loop over donors */
            {
#ifdef HAVE_OPENMP /* ====== Write some output ======\ */
#pragma omp critical
                {
                    dondata[thisThread] = i;
                    parallel_print(dondata, nThreads);
                }
#else
                fprintf(stderr, "\r %i", i);
#endif /* ===========================================/ */

                for (j=0; j<hb->a.nra; j++) /* loop over acceptors */
                {
                    for (nh=0; nh<hb->d.nhydro[i]; nh++) /* loop over donors' hydrogens */
                    {
                        E = hb->hbE.E[i][j][nh];
                        if (E != NULL)
                        {
                            for (k=0; k<nframes; k++)
                            {
                                if (E[k] != NONSENSE_E)
                                    rhbex[k] = (real)E[k];
                            }
                      
                            low_do_autocorr(NULL,oenv,NULL,nframes,1,-1,&(rhbex),hb->time[1]-hb->time[0],
                                            eacNormal,1,FALSE,bNorm,FALSE,0,-1,0,1);
#ifdef HAVE_OPENMP
#pragma omp critical
#endif
                            {
                                for(k=0; (k<nn); k++)
                                    ct[k] += rhbex[k];
                            }
                        }
                    }   /* k loop */
                }       /* j loop */
            }           /* i loop */
            sfree(rhbex);
#pragma omp barrier
#ifdef HAVE_OPENMP 
            /*                                                           # */
        } /* End of parallel block                                       # */
        /* ##############################################################/ */
        sfree(dondata);
#endif
        normalizeACF(ct, NULL, nn);
        snew(ctdouble, nn);
        snew(timedouble, nn);
        for (j=0; j<nn; j++)
        {
            timedouble[j]=(double)(hb->time[j]);
            ctdouble[j]=(double)(ct[j]);
        }

        /* Remove ballistic term */
        /* Ballistic component removal and fitting to the reversible geminate recombination model
         * will be taken out for the time being. First of all, one can remove the ballistic
         * component with g_analyze afterwards. Secondly, and more importantly, there are still
         * problems with the robustness of the fitting to the model. More work is needed.
         * A third reason is that we're currently using gsl for this and wish to reduce dependence
         * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with
         * a BSD-licence that can do the job.
         *
         * - Erik Marklund, June 18 2010. 
         */
/*         if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */
/*             takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */
/*         else */
/*             printf("\nNumber of data points is less than the number of parameters to fit\n." */
/*                    "The system is underdetermined, hence no ballistic term can be found.\n\n"); */

        fp = xvgropen(fn, "Hydrogen Bond Autocorrelation",output_env_get_xvgr_tlabel(oenv),"C(t)");
        xvgr_legend(fp,asize(legNN),legNN);
      
        for(j=0; (j<nn); j++)
            fprintf(fp,"%10g  %10g %10g\n",
                    hb->time[j]-hb->time[0],
                    ct[j],
                    ctdouble[j]);
        fclose(fp);
        sfree(ct);
        sfree(ctdouble);
        sfree(timedouble);
#endif /* HAVE_NN_LOOPS */
        break; /* case AC_NN */

    case AC_GEM:
        snew(ct,2*n2);
        memset(ct,0,2*n2*sizeof(real));
#ifndef HAVE_OPENMP
        fprintf(stderr, "Donor:\n");
#define __ACDATA ct
#else
#define __ACDATA p_ct
#endif

#ifdef HAVE_OPENMP /*  =========================================\
                    *                                          */
#pragma omp parallel default(none)                              \
    private(i, k, nh, hbh, pHist, h, g, n0, nf, np, j, m,               \
            pfound, poff, rHbExGem, p, ihb, mMax,               \
            thisThread, p_ct)                                   \
    shared(hb, dondata, ct, nn, nThreads, n2, stderr, bNorm,    \
           nframes, bMerge, bContact)
        { /* ##########  THE START OF THE ENORMOUS PARALLELIZED BLOCK!  ########## */
            h = NULL;
            g = NULL;
            thisThread = omp_get_thread_num();
            snew(h,hb->maxhydro);
            snew(g,hb->maxhydro);
            mMax = INT_MIN;
            rHbExGem = NULL;
            poff = NULL;
            pfound = NULL;
            p_ct = NULL;
            snew(p_ct,2*n2);
            memset(p_ct,0,2*n2*sizeof(real));

            /* I'm using a chunk size of 1, since I expect      \
             * the overhead to be really small compared         \
             * to the actual calculations                       \ */
#pragma omp for schedule(dynamic,1) nowait /*                   \ */
#endif /* HAVE_OPENMP  =========================================/ */
      
            for (i=0; i<hb->d.nrd; i++) {
#ifdef HAVE_OPENMP
#pragma omp critical
                {
                    dondata[thisThread] = i;
                    parallel_print(dondata, nThreads);
                }
#else
                fprintf(stderr, "\r %i", i);
#endif
        
                for (k=0; k<hb->a.nra; k++) {
                    for (nh=0; nh < ((bMerge || bContact) ? 1 : hb->d.nhydro[i]); nh++) {
                        hbh = hb->hbmap[i][k];
                        if (hbh) {
                            /* Note that if hb->per->gemtype==gemDD, then distances will be stored in
                             * hb->hbmap[d][a].h array anyway, because the contact flag will be set.
                             * hence, it's only with the gemAD mode that hb->hbmap[d][a].g will be used. */
                            pHist = &(hb->per->pHist[i][k]);
                            if (ISHB(hbh->history[nh]) && pHist->len != 0) {

/* No need for a critical section */
/* #ifdef HAVE_OPENMP */
/* #pragma omp critical */
/* #endif */
                                {
                                    h[nh] = hbh->h[nh];
                                    g[nh] = hb->per->gemtype==gemAD ? hbh->g[nh] : NULL;
                                }
                                n0 = hbh->n0;
                                nf = hbh->nframes;
                                /* count the number of periodic shifts encountered and store
                                 * them in separate arrays. */
                                np = 0;
                                for (j=0; j<pHist->len; j++)
                                {
                                    p = pHist->p[j];
                                    for (m=0; m<=np; m++) {
                                        if (m == np) { /* p not recognized in list. Add it and set up new array. */ 
                                            np++;
                                            if (np>hb->per->nper)
                                                gmx_fatal(FARGS, "Too many pshifts. Something's utterly wrong here.");
                                            if (m>=mMax) { /* Extend the arrays.
                                                            * Doing it like this, using mMax to keep track of the sizes,
                                                            * eleviates the need for freeing and re-allocating the arrays
                                                            * when taking on the next donor-acceptor pair */
                                                mMax = m;
                                                srenew(pfound,np);   /* The list of found periodic shifts. */
                                                srenew(rHbExGem,np); /* The hb existence functions (-aver_hb). */
                                                snew(rHbExGem[m],2*n2);
                                                srenew(poff,np);
                                            }

/* This shouldn't have to be critical, right? */
/* #ifdef HAVE_OPENMP */
/* #pragma omp critical */
/* #endif */
                                            {
                                                if (rHbExGem != NULL && rHbExGem[m] != NULL) {
                                                    /* This must be done, as this array was most likey
                                                     * used to store stuff in some previous iteration. */                           
                                                    memset(rHbExGem[m], 0, (sizeof(real)) * (2*n2));
                                                }
                                                else
                                                    fprintf(stderr, "rHbExGem not initialized! m = %i\n", m);
                                            }
                                            pfound[m] = p;
                                            poff[m] = -1;
                      
                                            break;
                                        } /* m==np */
                                        if (p == pfound[m])
                                            break;
                                    } /* m: Loop over found shifts */
                                }   /* j: Loop over shifts */

                                /* Now unpack and disentangle the existence funtions. */
                                for (j=0; j<nf; j++) {
                                    /* i:       donor,
                                     * k:       acceptor
                                     * nh:      hydrogen
                                     * j:       time
                                     * p:       periodic shift
                                     * pfound:  list of periodic shifts found for this pair.
                                     * poff:    list of frame offsets; that is, the first
                                     *          frame a hbond has a particular periodic shift. */
                                    p = getPshift(*pHist, j+n0);
                                    if (p != -1)
                                    {
                                        for (m=0; m<np; m++)
                                        {
                                            if (pfound[m] == p)
                                                break;
                                            if (m==(np-1))
                                                gmx_fatal(FARGS,"Shift not found, but must be there.");
                                        }

                                        ihb = is_hb(h[nh],j) || ((hb->per->gemtype!=gemAD || j==0) ? FALSE : is_hb(g[nh],j));
                                        if (ihb)
                                        {
                                            if (poff[m] == -1)
                                                poff[m] = j; /* Here's where the first hbond with shift p is,
                                                              * relative to the start of h[0].*/
                                            if (j<poff[m])
                                                gmx_fatal(FARGS, "j<poff[m]");
                                            rHbExGem[m][j-poff[m]] += 1;
                                        }
                                    }
                                }

                                /* Now, build ac. */
                                for (m=0; m<np; m++) {
                                    if (rHbExGem[m][0]>0  && n0+poff[m]<nn/*  && m==0 */) {
                                        low_do_autocorr(NULL,oenv,NULL,nframes,1,-1,&(rHbExGem[m]),hb->time[1]-hb->time[0],
                                                        eacNormal,1,FALSE,bNorm,FALSE,0,-1,0,1);
                                        for(j=0; (j<nn); j++)
                                            __ACDATA[j] += rHbExGem[m][j];
                                    }
                                } /* Building of ac. */
                            } /* if (ISHB(...*/
                        } /* if (hbh) */
                    } /* hydrogen loop */
                } /* acceptor loop */
            } /* donor loop */

            for (m=0; m<=mMax; m++) {
                sfree(rHbExGem[m]);
            }
            sfree(pfound);
            sfree(poff);
            sfree(rHbExGem);

            sfree(h);
            sfree(g);
#ifdef HAVE_OPENMP /* =======================================\ */
#pragma omp critical
            {
                for (i=0; i<nn; i++)
                    ct[i] += p_ct[i];
            }
            sfree(p_ct);

        } /* ########## THE END OF THE ENORMOUS PARALLELIZED BLOCK ########## */
        sfree(dondata);
#endif /* HAVE_OPENMP =======================================/ */

        normalizeACF(ct, NULL, nn);

        fprintf(stderr, "\n\nACF successfully calculated.\n");

        /* Use this part to fit to geminate recombination - JCP 129, 84505 (2008) */
      
        snew(ctdouble, nn);
        snew(timedouble, nn);
        snew(fittedct, nn);
    
        for (j=0;j<nn;j++){
            timedouble[j]=(double)(hb->time[j]);
            ctdouble[j]=(double)(ct[j]);
        }

        /* Remove ballistic term */
        /* Ballistic component removal and fitting to the reversible geminate recombination model
         * will be taken out for the time being. First of all, one can remove the ballistic
         * component with g_analyze afterwards. Secondly, and more importantly, there are still
         * problems with the robustness of the fitting to the model. More work is needed.
         * A third reason is that we're currently using gsl for this and wish to reduce dependence
         * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with
         * a BSD-licence that can do the job.
         *
         * - Erik Marklund, June 18 2010. 
         */
/*         if (bBallistic) { */
/*             if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */
/*                 takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */
/*             else */
/*                 printf("\nNumber of data points is less than the number of parameters to fit\n." */
/*                        "The system is underdetermined, hence no ballistic term can be found.\n\n"); */
/*         } */
/*         if (bGemFit) */
/*             fitGemRecomb(ctdouble, timedouble, &fittedct, nn, params); */
    

        if (bContact)
            fp = xvgropen(fn, "Contact Autocorrelation",output_env_get_xvgr_tlabel(oenv),"C(t)",oenv);
        else
            fp = xvgropen(fn, "Hydrogen Bond Autocorrelation",output_env_get_xvgr_tlabel(oenv),"C(t)",oenv);
        xvgr_legend(fp,asize(legGem),(const char**)legGem,oenv);

        for(j=0; (j<nn); j++)
        {
            fprintf(fp, "%10g  %10g", hb->time[j]-hb->time[0],ct[j]);
            if (bBallistic)
                fprintf(fp,"  %10g", ctdouble[j]);
            if (bGemFit)
                fprintf(fp,"  %10g", fittedct[j]);
            fprintf(fp,"\n");
        }
        fclose(fp);

        sfree(ctdouble);
        sfree(timedouble);
        sfree(fittedct);
        sfree(ct);

        break; /* case AC_GEM */

    case AC_LUZAR:
        snew(rhbex,2*n2);
        snew(ct,2*n2);
        snew(gt,2*n2);
        snew(ht,2*n2);
        snew(ght,2*n2);
        snew(dght,2*n2);

        snew(kt,nn);
        snew(cct,nn);
    
        for(i=0; (i<hb->d.nrd); i++) {
            for(k=0; (k<hb->a.nra); k++) {
                nhydro = 0;
                hbh = hb->hbmap[i][k];
   
                if (hbh) {
                    if (bMerge || bContact) {
                        if (ISHB(hbh->history[0])) {
                            h[0] = hbh->h[0];
                            g[0] = hbh->g[0];
                            nhydro = 1;
                        }
                    }
                    else {
                        for(m=0; (m<hb->maxhydro); m++)
                            if (bContact ? ISDIST(hbh->history[m]) : ISHB(hbh->history[m])) {
                                g[nhydro] = hbh->g[m];
                                h[nhydro] = hbh->h[m];
                                nhydro++;
                            }
                    }
        
                    nf = hbh->nframes;
                    for(nh=0; (nh<nhydro); nh++) {
                        int nrint = bContact ? hb->nrdist : hb->nrhb;
                        if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint))
                            fprintf(stderr,"\rACF %d/%d",nhbonds+1,nrint);
                        nhbonds++;
                        for(j=0; (j<nframes); j++) {
                            /* Changed '<' into '<=' below, just like I did in
                               the hbm-output-loop in the gmx_hbond() block.
                               - Erik Marklund, May 31, 2006 */
                            if (j <= nf) {
                                ihb   = is_hb(h[nh],j);
                                idist = is_hb(g[nh],j);
                            }
                            else {
                                ihb = idist = 0;
                            }
                            rhbex[j] = ihb;
                            /* For contacts: if a second cut-off is provided, use it,
                             * otherwise use g(t) = 1-h(t) */
                            if (!R2 && bContact)
                                gt[j]  = 1-ihb;
                            else
                                gt[j]  = idist*(1-ihb); 
                            ht[j]    = rhbex[j];
                            nhb     += ihb;
                        }
          

                        /* The autocorrelation function is normalized after summation only */
                        low_do_autocorr(NULL,oenv,NULL,nframes,1,-1,&rhbex,hb->time[1]-hb->time[0],
                                        eacNormal,1,FALSE,bNorm,FALSE,0,-1,0,1);
          
                        /* Cross correlation analysis for thermodynamics */
                        for(j=nframes; (j<n2); j++) {
                            ht[j] = 0;
                            gt[j] = 0;
                        }

                        cross_corr(n2,ht,gt,dght);
          
                        for(j=0; (j<nn); j++) {
                            ct[j]  += rhbex[j];
                            ght[j] += dght[j];
                        }
                    }
                }
            }
        }
        fprintf(stderr,"\n");
        sfree(h);
        sfree(g);
        normalizeACF(ct, gt, nn);

        /* Determine tail value for statistics */
        tail  = 0;
        tail2 = 0;
        for(j=nn/2; (j<nn); j++) {
            tail  += ct[j];
            tail2 += ct[j]*ct[j];
        }
        tail  /= (nn - nn/2);
        tail2 /= (nn - nn/2);
        dtail  = sqrt(tail2-tail*tail);
  
        /* Check whether the ACF is long enough */
        if (dtail > tol) {
            printf("\nWARNING: Correlation function is probably not long enough\n"
                   "because the standard deviation in the tail of C(t) > %g\n"
                   "Tail value (average C(t) over second half of acf): %g +/- %g\n",
                   tol,tail,dtail);
        }
        for(j=0; (j<nn); j++) {
            cct[j] = ct[j];
            ct[j]  = (cct[j]-tail)/(1-tail); 
        }
        /* Compute negative derivative k(t) = -dc(t)/dt */
        compute_derivative(nn,hb->time,ct,kt);
        for(j=0; (j<nn); j++)
            kt[j] = -kt[j];


        if (bContact)
            fp = xvgropen(fn, "Contact Autocorrelation",output_env_get_xvgr_tlabel(oenv),"C(t)",oenv);
        else
            fp = xvgropen(fn, "Hydrogen Bond Autocorrelation",output_env_get_xvgr_tlabel(oenv),"C(t)",oenv);
        xvgr_legend(fp,asize(legLuzar),legLuzar, oenv);

      
        for(j=0; (j<nn); j++)
            fprintf(fp,"%10g  %10g  %10g  %10g  %10g\n",
                    hb->time[j]-hb->time[0],ct[j],cct[j],ght[j],kt[j]);
        ffclose(fp);

        analyse_corr(nn,hb->time,ct,ght,kt,NULL,NULL,NULL,
                     fit_start,temp,smooth_tail_start,oenv);
  
        do_view(oenv,fn,NULL);
        sfree(rhbex);
        sfree(ct);
        sfree(gt);
        sfree(ht);
        sfree(ght);
        sfree(dght);
        sfree(cct);
        sfree(kt);
        /* sfree(h); */
/*         sfree(g); */

        break; /* case AC_LUZAR */

    default:
        gmx_fatal(FARGS, "Unrecognized type of ACF-calulation. acType = %i.", acType);
    } /* switch (acType) */
}

static void init_hbframe(t_hbdata *hb,int nframes,real t)
{
    int i,j,m;
  
    hb->time[nframes]   = t;
    hb->nhb[nframes]    = 0;
    hb->ndist[nframes]  = 0;
    for (i=0; (i<max_hx); i++)
        hb->nhx[nframes][i]=0;
    /* Loop invalidated */
    if (hb->bHBmap && 0)
        for (i=0; (i<hb->d.nrd); i++)
            for (j=0; (j<hb->a.nra); j++)
                for (m=0; (m<hb->maxhydro); m++)
                    if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[m])
                        set_hb(hb,i,m,j,nframes,HB_NO);
    /*set_hb(hb->hbmap[i][j]->h[m],nframes-hb->hbmap[i][j]->n0,HB_NO);*/
}

static void analyse_donor_props(const char *fn,t_hbdata *hb,int nframes,real t,
                                const output_env_t oenv)
{
    static FILE *fp = NULL;
    const char *leg[] = { "Nbound", "Nfree" };
    int i,j,k,nbound,nb,nhtot;
  
    if (!fn)
        return;
    if (!fp) {
        fp = xvgropen(fn,"Donor properties",output_env_get_xvgr_tlabel(oenv),"Number",oenv);
        xvgr_legend(fp,asize(leg),leg,oenv);
    }
    nbound = 0;
    nhtot  = 0;
    for(i=0; (i<hb->d.nrd); i++) {
        for(k=0; (k<hb->d.nhydro[i]); k++) {
            nb = 0;
            nhtot ++;
            for(j=0; (j<hb->a.nra) && (nb == 0); j++) {
                if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[k] && 
                    is_hb(hb->hbmap[i][j]->h[k],nframes)) 
                    nb = 1;
            }
            nbound += nb;
        }
    }
    fprintf(fp,"%10.3e  %6d  %6d\n",t,nbound,nhtot-nbound);
}

static void dump_hbmap(t_hbdata *hb,
                       int nfile,t_filenm fnm[],gmx_bool bTwo,
                       gmx_bool bContact, int isize[],int *index[],char *grpnames[],
                       t_atoms *atoms)
{
    FILE *fp,*fplog;
    int  ddd,hhh,aaa,i,j,k,m,grp;
    char ds[32],hs[32],as[32];
    gmx_bool first;
  
    fp = opt2FILE("-hbn",nfile,fnm,"w");
    if (opt2bSet("-g",nfile,fnm)) {
        fplog = ffopen(opt2fn("-g",nfile,fnm),"w");
        if (bContact)
            fprintf(fplog,"# %10s  %12s  %12s\n","Donor","Hydrogen","Acceptor");
        else
            fprintf(fplog,"# %10s  %12s  %12s\n","Donor","Hydrogen","Acceptor");
    }
    else
        fplog = NULL;
    for (grp=gr0; grp<=(bTwo?gr1:gr0); grp++) {
        fprintf(fp,"[ %s ]",grpnames[grp]);
        for (i=0; i<isize[grp]; i++) {
            fprintf(fp,(i%15)?" ":"\n");
            fprintf(fp," %4u",index[grp][i]+1);
        }
        fprintf(fp,"\n");
        /*
          Added -contact support below.
          - Erik Marklund, May 29, 2006
        */
        if (!bContact) {
            fprintf(fp,"[ donors_hydrogens_%s ]\n",grpnames[grp]);
            for (i=0; (i<hb->d.nrd); i++) {
                if (hb->d.grp[i] == grp) { 
                    for(j=0; (j<hb->d.nhydro[i]); j++)
                        fprintf(fp," %4u %4u",hb->d.don[i]+1,
                                hb->d.hydro[i][j]+1);
                    fprintf(fp,"\n");
                }
            }
            first = TRUE;
            fprintf(fp,"[ acceptors_%s ]",grpnames[grp]);
            for (i=0; (i<hb->a.nra); i++) {
                if (hb->a.grp[i] == grp) { 
                    fprintf(fp,(i%15 && !first)?" ":"\n");
                    fprintf(fp," %4u",hb->a.acc[i]+1);
                    first = FALSE;
                }
            }
            fprintf(fp,"\n");
        }
    }
    if (bTwo)
        fprintf(fp,bContact ? "[ contacts_%s-%s ]\n" :
                "[ hbonds_%s-%s ]\n",grpnames[0],grpnames[1]);
    else
        fprintf(fp,bContact ? "[ contacts_%s ]" : "[ hbonds_%s ]\n",grpnames[0]);
  
    for(i=0; (i<hb->d.nrd); i++) {
        ddd = hb->d.don[i];
        for(k=0; (k<hb->a.nra); k++) {
            aaa = hb->a.acc[k];
            for(m=0; (m<hb->d.nhydro[i]); m++) {
                if (hb->hbmap[i][k] && ISHB(hb->hbmap[i][k]->history[m])) {
                    sprintf(ds,"%s",mkatomname(atoms,ddd));
                    sprintf(as,"%s",mkatomname(atoms,aaa));
                    if (bContact) {
                        fprintf(fp," %6u %6u\n",ddd+1,aaa+1);
                        if (fplog) 
                            fprintf(fplog,"%12s  %12s\n",ds,as);
                    } else {
                        hhh = hb->d.hydro[i][m];
                        sprintf(hs,"%s",mkatomname(atoms,hhh));
                        fprintf(fp," %6u %6u %6u\n",ddd+1,hhh+1,aaa+1);
                        if (fplog) 
                            fprintf(fplog,"%12s  %12s  %12s\n",ds,hs,as);
                    }
                }
            }
        }
    }
    ffclose(fp);
    if (fplog)
        ffclose(fplog);
}

#ifdef HAVE_OPENMP
/* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template.
 * It mimics add_frames() and init_frame() to some extent. */
static void sync_hbdata(t_hbdata *hb, t_hbdata *p_hb,
                        int nframes, real t)
{
    int i;
    if (nframes >= p_hb->max_frames)
    {
        p_hb->max_frames += 4096;
        srenew(p_hb->nhb,   p_hb->max_frames);
        srenew(p_hb->ndist, p_hb->max_frames);
        srenew(p_hb->n_bound, p_hb->max_frames);
        srenew(p_hb->nhx, p_hb->max_frames);
        if (p_hb->bDAnr)
            srenew(p_hb->danr, p_hb->max_frames);
        memset(&(p_hb->nhb[nframes]),   0, sizeof(int) * (p_hb->max_frames-nframes));
        memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
        p_hb->nhb[nframes] = 0;
        p_hb->ndist[nframes] = 0;

    }
    p_hb->nframes = nframes;
/*     for (i=0;) */
/*     { */
/*         p_hb->nhx[nframes][i] */
/*     } */
    memset(&(p_hb->nhx[nframes]), 0 ,sizeof(int)*max_hx); /* zero the helix count for this frame */

    /* hb->per will remain constant througout the frame loop,
     * even though the data its members point to will change,
     * hence no need for re-syncing. */
}
#endif

int gmx_hbond(int argc,char *argv[])
{
    const char *desc[] = {
        "[TT]g_hbond[tt] computes and analyzes hydrogen bonds. Hydrogen bonds are",
        "determined based on cutoffs for the angle Acceptor - Donor - Hydrogen",
        "(zero is extended) and the distance Hydrogen - Acceptor.",
        "OH and NH groups are regarded as donors, O is an acceptor always,",
        "N is an acceptor by default, but this can be switched using",
        "[TT]-nitacc[tt]. Dummy hydrogen atoms are assumed to be connected",
        "to the first preceding non-hydrogen atom.[PAR]",
    
        "You need to specify two groups for analysis, which must be either",
        "identical or non-overlapping. All hydrogen bonds between the two",
        "groups are analyzed.[PAR]",
    
        "If you set [TT]-shell[tt], you will be asked for an additional index group",
        "which should contain exactly one atom. In this case, only hydrogen",
        "bonds between atoms within the shell distance from the one atom are",
        "considered.[PAR]",
    
        /*    "It is also possible to analyse specific hydrogen bonds with",
              "[TT]-sel[tt]. This index file must contain a group of atom triplets",
              "Donor Hydrogen Acceptor, in the following way:[PAR]",
        */
        "[TT]",
        "[ selected ][BR]",
        "     20    21    24[BR]",
        "     25    26    29[BR]",
        "      1     3     6[BR]",
        "[tt][BR]",
        "Note that the triplets need not be on separate lines.",
        "Each atom triplet specifies a hydrogen bond to be analyzed,",
        "note also that no check is made for the types of atoms.[PAR]",
     
        "[BB]Output:[bb][BR]",
        "[TT]-num[tt]:  number of hydrogen bonds as a function of time.[BR]",
        "[TT]-ac[tt]:   average over all autocorrelations of the existence",
        "functions (either 0 or 1) of all hydrogen bonds.[BR]",
        "[TT]-dist[tt]: distance distribution of all hydrogen bonds.[BR]",
        "[TT]-ang[tt]:  angle distribution of all hydrogen bonds.[BR]",
        "[TT]-hx[tt]:   the number of n-n+i hydrogen bonds as a function of time",
        "where n and n+i stand for residue numbers and i ranges from 0 to 6.",
        "This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated",
        "with helices in proteins.[BR]",
        "[TT]-hbn[tt]:  all selected groups, donors, hydrogens and acceptors",
        "for selected groups, all hydrogen bonded atoms from all groups and",
        "all solvent atoms involved in insertion.[BR]",
        "[TT]-hbm[tt]:  existence matrix for all hydrogen bonds over all",
        "frames, this also contains information on solvent insertion",
        "into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]",
        "index file.[BR]",
        "[TT]-dan[tt]: write out the number of donors and acceptors analyzed for",
        "each timeframe. This is especially useful when using [TT]-shell[tt].[BR]",
        "[TT]-nhbdist[tt]: compute the number of HBonds per hydrogen in order to",
        "compare results to Raman Spectroscopy.",
        "[PAR]",
        "Note: options [TT]-ac[tt], [TT]-life[tt], [TT]-hbn[tt] and [TT]-hbm[tt]",
        "require an amount of memory proportional to the total numbers of donors",
        "times the total number of acceptors in the selected group(s)."
    };
  
    static real acut=30, abin=1, rcut=0.35, r2cut=0, rbin=0.005, rshell=-1;
    static real maxnhb=0,fit_start=1,fit_end=60,temp=298.15,smooth_tail_start=-1, D=-1;
    static gmx_bool bNitAcc=TRUE,bDA=TRUE,bMerge=TRUE;
    static int  nDump=0, nFitPoints=100;
    static int nThreads = 0, nBalExp=4;

    static gmx_bool bContact=FALSE, bBallistic=FALSE, bBallisticDt=FALSE, bGemFit=FALSE;
    static real logAfterTime = 10, gemBallistic = 0.2; /* ps */
    static const char *NNtype[] = {NULL, "none", "binary", "oneOverR3", "dipole", NULL};

    /* options */
    t_pargs pa [] = {
        { "-a",    FALSE,  etREAL, {&acut},
          "Cutoff angle (degrees, Acceptor - Donor - Hydrogen)" },
        { "-r",    FALSE,  etREAL, {&rcut},
          "Cutoff radius (nm, X - Acceptor, see next option)" },
        { "-da",   FALSE,  etBOOL, {&bDA},
          "Use distance Donor-Acceptor (if TRUE) or Hydrogen-Acceptor (FALSE)" },
        { "-r2",   FALSE,  etREAL, {&r2cut},
          "Second cutoff radius. Mainly useful with [TT]-contact[tt] and [TT]-ac[tt]"},
        { "-abin", FALSE,  etREAL, {&abin},
          "Binwidth angle distribution (degrees)" },
        { "-rbin", FALSE,  etREAL, {&rbin},
          "Binwidth distance distribution (nm)" },
        { "-nitacc",FALSE, etBOOL, {&bNitAcc},
          "Regard nitrogen atoms as acceptors" },
        { "-contact",FALSE,etBOOL, {&bContact},
          "Do not look for hydrogen bonds, but merely for contacts within the cut-off distance" },
        { "-shell", FALSE, etREAL, {&rshell},
          "when > 0, only calculate hydrogen bonds within # nm shell around "
          "one particle" },
        { "-fitstart", FALSE, etREAL, {&fit_start},
          "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]" },
        { "-fitstart", FALSE, etREAL, {&fit_start},
          "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])" },
        { "-temp",  FALSE, etREAL, {&temp},
          "Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" },
        { "-smooth",FALSE, etREAL, {&smooth_tail_start},
          "If >= 0, the tail of the ACF will be smoothed by fitting it to an exponential function: y = A exp(-x/[GRK]tau[grk])" },
        { "-dump",  FALSE, etINT, {&nDump},
          "Dump the first N hydrogen bond ACFs in a single [TT].xvg[tt] file for debugging" },
        { "-max_hb",FALSE, etREAL, {&maxnhb},
          "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" },
        { "-merge", FALSE, etBOOL, {&bMerge},
          "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." },
        { "-geminate", FALSE, etENUM, {gemType},
          "Use reversible geminate recombination for the kinetics/thermodynamics calclations. See Markovitch et al., J. Chem. Phys 129, 084505 (2008) for details."},
        { "-diff", FALSE, etREAL, {&D},
          "Dffusion coefficient to use in the reversible geminate recombination kinetic model. If negative, then it will be fitted to the ACF along with ka and kd."},
#ifdef HAVE_OPENMP
        { "-nthreads", FALSE, etINT, {&nThreads},
          "Number of threads used for the parallel loop over autocorrelations. nThreads <= 0 means maximum number of threads. Requires linking with OpenMP. The number of threads is limited by the number of processors (before OpenMP v.3 ) or environment variable OMP_THREAD_LIMIT (OpenMP v.3)"},
#endif
    };
    const char *bugs[] = {
        "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and therefore not available for the time being."
    };
    t_filenm fnm[] = {
        { efTRX, "-f",   NULL,     ffREAD  },
        { efTPX, NULL,   NULL,     ffREAD  },
        { efNDX, NULL,   NULL,     ffOPTRD },
        /*    { efNDX, "-sel", "select", ffOPTRD },*/
        { efXVG, "-num", "hbnum",  ffWRITE },
        { efLOG, "-g",   "hbond",  ffOPTWR },
        { efXVG, "-ac",  "hbac",   ffOPTWR },
        { efXVG, "-dist","hbdist", ffOPTWR },
        { efXVG, "-ang", "hbang",  ffOPTWR },
        { efXVG, "-hx",  "hbhelix",ffOPTWR },
        { efNDX, "-hbn", "hbond",  ffOPTWR },
        { efXPM, "-hbm", "hbmap",  ffOPTWR },
        { efXVG, "-don", "donor",  ffOPTWR },
        { efXVG, "-dan", "danum",  ffOPTWR },
        { efXVG, "-life","hblife", ffOPTWR },
        { efXVG, "-nhbdist", "nhbdist",ffOPTWR }
    
    };
#define NFILE asize(fnm)
  
    char  hbmap [HB_NR]={ ' ',    'o',      '-',       '*' };
    const char *hbdesc[HB_NR]={ "None", "Present", "Inserted", "Present & Inserted" };
    t_rgb hbrgb [HB_NR]={ {1,1,1},{1,0,0},   {0,0,1},    {1,0,1} };

    t_trxstatus *status;
    int trrStatus=1;
    t_topology top;
    t_inputrec ir;
    t_pargs *ppa;
    int     npargs,natoms,nframes=0,shatom;
    int     *isize;
    char    **grpnames;
    atom_id **index;
    rvec    *x,hbox;
    matrix  box;
    real    t,ccut,dist=0.0,ang=0.0;
    double  max_nhb,aver_nhb,aver_dist;
    int     h=0,i,j,k=0,l,start,end,id,ja,ogrp,nsel;
    int     xi,yi,zi,ai;
    int     xj,yj,zj,aj,xjj,yjj,zjj;
    int     xk,yk,zk,ak,xkk,ykk,zkk;
    gmx_bool    bSelected,bHBmap,bStop,bTwo,was,bBox,bTric;
    int     *adist,*rdist;
    int        grp,nabin,nrbin,bin,resdist,ihb;
    char       **leg;
    t_hbdata   *hb;
    FILE       *fp,*fpins=NULL,*fpnhb=NULL;
    t_gridcell ***grid;
    t_ncell    *icell,*jcell,*kcell;
    ivec       ngrid;
    unsigned char        *datable;
    output_env_t oenv;
    int     gemmode, NN;
    PSTYPE  peri=0;
    t_E     E;
    int     ii, jj, hh, actual_nThreads;
    int     threadNr=0;
    gmx_bool    bGem, bNN, bParallel;
    t_gemParams *params=NULL;
    
    CopyRight(stdout,argv[0]);

    npargs = asize(pa);  
    ppa    = add_acf_pargs(&npargs,pa);
  
    parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,NFILE,fnm,npargs,
                      ppa,asize(desc),desc,asize(bugs),bugs,&oenv);

    /* NN-loop? If so, what estimator to use ?*/
    NN = 1;
    /* Outcommented for now DvdS 2010-07-13
    while (NN < NN_NR && gmx_strcasecmp(NNtype[0], NNtype[NN])!=0)
        NN++;
    if (NN == NN_NR)
        gmx_fatal(FARGS, "Invalid NN-loop type.");
    */
    bNN = FALSE;
    for (i=2; bNN==FALSE && i<NN_NR; i++)
        bNN = bNN || NN == i;

    if (NN > NN_NONE && bMerge)
        bMerge = FALSE;

    /* geminate recombination? If so, which flavor? */
    gemmode = 1;
    while (gemmode < gemNR && gmx_strcasecmp(gemType[0], gemType[gemmode])!=0)
        gemmode++;
    if (gemmode == gemNR)
        gmx_fatal(FARGS, "Invalid recombination type.");
  
    bGem = FALSE;
    for (i=2; bGem==FALSE && i<gemNR; i++)
        bGem = bGem || gemmode == i;
  
    if (bGem) {
        printf("Geminate recombination: %s\n" ,gemType[gemmode]);
#ifndef HAVE_LIBGSL
        printf("Note that some aspects of reversible geminate recombination won't work without gsl.\n");
#endif
        if (bContact) {
            if (gemmode != gemDD) {
                printf("Turning off -contact option...\n");
                bContact = FALSE;
            }
        } else {
            if (gemmode == gemDD) {
                printf("Turning on -contact option...\n");
                bContact = TRUE;
            }
        }
        if (bMerge) {
            if (gemmode == gemAA) {
                printf("Turning off -merge option...\n");
                bMerge = FALSE;
            }
        } else {
            if (gemmode != gemAA) {
                printf("Turning on -merge option...\n");
                bMerge = TRUE;
            }
        }
    } 
    
    /* process input */
    bSelected = FALSE;
    ccut = cos(acut*DEG2RAD);
  
    if (bContact) {
        if (bSelected)
            gmx_fatal(FARGS,"Can not analyze selected contacts.");
        if (!bDA) {
            gmx_fatal(FARGS,"Can not analyze contact between H and A: turn off -noda");
        }
    }

#ifndef HAVE_LIBGSL
    /* Don't pollute stdout with information about external libraries.
     *
     * printf("NO GSL! Can't find and take away ballistic term in ACF without GSL\n.");
     */
#endif
  
    /* Initiate main data structure! */
    bHBmap = (opt2bSet("-ac",NFILE,fnm) ||
              opt2bSet("-life",NFILE,fnm) ||
              opt2bSet("-hbn",NFILE,fnm) || 
              opt2bSet("-hbm",NFILE,fnm) ||
              bGem);
  
#ifdef HAVE_OPENMP
    /* Same thing here. There is no reason whatsoever to write the specific version of
     * OpenMP used for compilation to stdout for normal usage.
     *
     * printf("Compiled with OpenMP (%i)\n", _OPENMP);
     */
#endif

    /*   if (bContact && bGem) */
    /*     gmx_fatal(FARGS, "Can't do reversible geminate recombination with -contact yet."); */

    if (opt2bSet("-nhbdist",NFILE,fnm)) {
        const char *leg[MAXHH+1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" };
        fpnhb = xvgropen(opt2fn("-nhbdist",NFILE,fnm),
                         "Number of donor-H with N HBs",output_env_get_xvgr_tlabel(oenv),"N",oenv);
        xvgr_legend(fpnhb,asize(leg),leg,oenv);
    }
  
    hb = mk_hbdata(bHBmap,opt2bSet("-dan",NFILE,fnm),bMerge || bContact, bGem, gemmode);

    /* get topology */
    read_tpx_top(ftp2fn(efTPX,NFILE,fnm),&ir,box,&natoms,NULL,NULL,NULL,&top);
  
    snew(grpnames,grNR);
    snew(index,grNR);
    snew(isize,grNR);
    /* Make Donor-Acceptor table */
    snew(datable, top.atoms.nr);
    gen_datable(index[0],isize[0],datable,top.atoms.nr);
  
    if (bSelected) {
        /* analyze selected hydrogen bonds */
        printf("Select group with selected atoms:\n");
        get_index(&(top.atoms),opt2fn("-sel",NFILE,fnm),
                  1,&nsel,index,grpnames);
        if (nsel % 3)
            gmx_fatal(FARGS,"Number of atoms in group '%s' not a multiple of 3\n"
                      "and therefore cannot contain triplets of "
                      "Donor-Hydrogen-Acceptor",grpnames[0]);
        bTwo=FALSE;
    
        for(i=0; (i<nsel); i+=3) {
            int dd = index[0][i];
            int aa = index[0][i+2];
            /* int */ hh = index[0][i+1];
            add_dh (&hb->d,dd,hh,i,datable);
            add_acc(&hb->a,aa,i);
            /* Should this be here ? */
            snew(hb->d.dptr,top.atoms.nr);
            snew(hb->a.aptr,top.atoms.nr);
            add_hbond(hb,dd,aa,hh,gr0,gr0,0,bMerge,0,bContact,peri);
        }
        printf("Analyzing %d selected hydrogen bonds from '%s'\n",
               isize[0],grpnames[0]);
    } 
    else {
        /* analyze all hydrogen bonds: get group(s) */
        printf("Specify 2 groups to analyze:\n");
        get_index(&(top.atoms),ftp2fn_null(efNDX,NFILE,fnm),
                  2,isize,index,grpnames);
    
        /* check if we have two identical or two non-overlapping groups */
        bTwo = isize[0] != isize[1];
        for (i=0; (i<isize[0]) && !bTwo; i++)
            bTwo = index[0][i] != index[1][i];
        if (bTwo) {
            printf("Checking for overlap in atoms between %s and %s\n",
                   grpnames[0],grpnames[1]);
            for (i=0; i<isize[1];i++)
                if (ISINGRP(datable[index[1][i]]))
                    gmx_fatal(FARGS,"Partial overlap between groups '%s' and '%s'",
                              grpnames[0],grpnames[1]);
            /*
              printf("Checking for overlap in atoms between %s and %s\n",
              grpnames[0],grpnames[1]);
              for (i=0; i<isize[0]; i++)
              for (j=0; j<isize[1]; j++)
              if (index[0][i] == index[1][j]) 
              gmx_fatal(FARGS,"Partial overlap between groups '%s' and '%s'",
              grpnames[0],grpnames[1]);
            */
        }
        if (bTwo)
            printf("Calculating %s "
                   "between %s (%d atoms) and %s (%d atoms)\n",
                   bContact ? "contacts" : "hydrogen bonds",
                   grpnames[0],isize[0],grpnames[1],isize[1]);
        else
            fprintf(stderr,"Calculating %s in %s (%d atoms)\n",
                    bContact?"contacts":"hydrogen bonds",grpnames[0],isize[0]);
    }
    sfree(datable);
  
    /* search donors and acceptors in groups */
    snew(datable, top.atoms.nr);
    for (i=0; (i<grNR); i++)
        if ( ((i==gr0) && !bSelected ) ||
             ((i==gr1) && bTwo )) {
            gen_datable(index[i],isize[i],datable,top.atoms.nr);
            if (bContact) {
                search_acceptors(&top,isize[i],index[i],&hb->a,i,
                                 bNitAcc,TRUE,(bTwo && (i==gr0)) || !bTwo, datable);
                search_donors   (&top,isize[i],index[i],&hb->d,i,
                                 TRUE,(bTwo && (i==gr1)) || !bTwo, datable);
            }
            else {
                search_acceptors(&top,isize[i],index[i],&hb->a,i,bNitAcc,FALSE,TRUE, datable);
                search_donors   (&top,isize[i],index[i],&hb->d,i,FALSE,TRUE, datable);
            }
            if (bTwo)
                clear_datable_grp(datable,top.atoms.nr);
        }
    sfree(datable);
    printf("Found %d donors and %d acceptors\n",hb->d.nrd,hb->a.nra);
    /*if (bSelected)
      snew(donors[gr0D], dons[gr0D].nrd);*/

    if (bHBmap) {
        printf("Making hbmap structure...");
        /* Generate hbond data structure */
        mk_hbmap(hb,bTwo);
        printf("done.\n");
    }

#ifdef HAVE_NN_LOOPS
    if (bNN)
        mk_hbEmap(hb, 0);
#endif

    if (bGem) {
        printf("Making per structure...");
        /* Generate hbond data structure */
        mk_per(hb);
        printf("done.\n");
    }
  
    /* check input */
    bStop=FALSE;
    if (hb->d.nrd + hb->a.nra == 0) {
        printf("No Donors or Acceptors found\n");
        bStop=TRUE;
    }
    if (!bStop) {
        if (hb->d.nrd == 0) {
            printf("No Donors found\n");
            bStop=TRUE;
        }
        if (hb->a.nra == 0) {
            printf("No Acceptors found\n");
            bStop=TRUE;
        }
    }
    if (bStop)
        gmx_fatal(FARGS,"Nothing to be done");

    shatom=0;
    if (rshell > 0) {
        int shisz;
        atom_id *shidx;
        char *shgrpnm;
        /* get index group with atom for shell */
        do {
            printf("Select atom for shell (1 atom):\n");
            get_index(&(top.atoms),ftp2fn_null(efNDX,NFILE,fnm),
                      1,&shisz,&shidx,&shgrpnm);
            if (shisz != 1)
                printf("group contains %d atoms, should be 1 (one)\n",shisz);
        } while(shisz != 1);
        shatom = shidx[0];
        printf("Will calculate hydrogen bonds within a shell "
               "of %g nm around atom %i\n",rshell,shatom+1);
    }

    /* Analyze trajectory */
    natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
    if ( natoms > top.atoms.nr )
        gmx_fatal(FARGS,"Topology (%d atoms) does not match trajectory (%d atoms)",
                  top.atoms.nr,natoms);
                
    bBox  = ir.ePBC!=epbcNONE;
    grid  = init_grid(bBox, box, (rcut>r2cut)?rcut:r2cut, ngrid);
    nabin = acut/abin;
    nrbin = rcut/rbin;
    snew(adist,nabin+1);
    snew(rdist,nrbin+1);

    if (bGem && !bBox)
        gmx_fatal(FARGS, "Can't do geminate recombination without periodic box.");

    bParallel = FALSE;

#ifndef HAVE_OPENMP
#define __ADIST adist
#define __RDIST rdist
#define __HBDATA hb
#else /* HAVE_OPENMP ================================================== \
       * Set up the OpenMP stuff,                                       |
       * like the number of threads and such                            |
       * Also start the parallel loop.                                  |
       */
#define __ADIST p_adist[threadNr]
#define __RDIST p_rdist[threadNr]
#define __HBDATA p_hb[threadNr]

    bParallel = !bSelected;

    if (bParallel)
    {
#if (_OPENMP > 200805)
        actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_thread_limit());
#else
        actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_num_procs());
#endif
        omp_set_num_threads(actual_nThreads);
        printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads);
        fflush(stdout);
    }
    else
    {
        actual_nThreads = 1;
    }

    t_hbdata **p_hb;          /* one per thread, then merge after the frame loop */
    int **p_adist, **p_rdist; /* a histogram for each thread. */
    snew(p_hb,    actual_nThreads);
    snew(p_adist, actual_nThreads);
    snew(p_rdist, actual_nThreads);
    for (i=0; i<actual_nThreads; i++)
    {
        snew(p_hb[i], 1);
        snew(p_adist[i], nabin+1);
        snew(p_rdist[i], nrbin+1);

        p_hb[i]->max_frames = 0;
        p_hb[i]->nhb = NULL;
        p_hb[i]->ndist = NULL;
        p_hb[i]->n_bound = NULL;
        p_hb[i]->time = NULL;
        p_hb[i]->nhx = NULL;

        p_hb[i]->bHBmap     = hb->bHBmap;
        p_hb[i]->bDAnr      = hb->bDAnr;
        p_hb[i]->bGem       = hb->bGem;
        p_hb[i]->wordlen    = hb->wordlen;
        p_hb[i]->nframes    = hb->nframes;
        p_hb[i]->maxhydro   = hb->maxhydro;
        p_hb[i]->danr       = hb->danr;
        p_hb[i]->d          = hb->d;
        p_hb[i]->a          = hb->a;
        p_hb[i]->hbmap      = hb->hbmap;
        p_hb[i]->time       = hb->time; /* This may need re-syncing at every frame. */
        p_hb[i]->per        = hb->per;

#ifdef HAVE_NN_LOOPS
        p_hb[i]->hbE = hb->hbE;
#endif

        p_hb[i]->nrhb   = 0;
        p_hb[i]->nrdist = 0;
    }
  
    /* Make a thread pool here,
     * instead of forking anew at every frame. */
  
#pragma omp parallel                                    \
    private(i, j, h, ii, jj, hh, E,                     \
            xi, yi, zi, xj, yj, zj, threadNr,           \
            dist, ang, peri, icell, jcell,              \
            grp, ogrp, ai, aj, xjj, yjj, zjj,           \
            xk, yk, zk, ihb, id,  resdist,              \
            xkk, ykk, zkk, kcell, ak, k, bTric)         \
    default(none)                                       \
    shared(hb, p_hb, p_adist, p_rdist, actual_nThreads, \
           x, bBox, box, hbox, rcut, r2cut, rshell,     \
           shatom, ngrid, grid, nframes, t,             \
           bParallel, bNN, index, bMerge, bContact,     \
           bTwo, bDA,ccut, abin, rbin, top,             \
           bSelected, bDebug, stderr, nsel,             \
           bGem, oenv, fnm, fpnhb, trrStatus, natoms,   \
           status, nabin, nrbin, adist, rdist, debug)
    {    /* Start of parallel region */
        threadNr = omp_get_thread_num();
#endif /* HAVE_OPENMP ================================================= */
        do
        {
            bTric = bBox && TRICLINIC(box);

#ifdef HAVE_OPENMP
            sync_hbdata(hb, p_hb[threadNr], nframes, t);
#pragma omp single
#endif
            {
                build_grid(hb,x,x[shatom], bBox,box,hbox, (rcut>r2cut)?rcut:r2cut, 
                           rshell, ngrid,grid);
                reset_nhbonds(&(hb->d));
        
                if (debug && bDebug)
                    dump_grid(debug, ngrid, grid);
                
                add_frames(hb,nframes);
                init_hbframe(hb,nframes,output_env_conv_time(oenv,t));
        
                if (hb->bDAnr)
                    count_da_grid(ngrid, grid, hb->danr[nframes]);
            } /* omp single */

#ifdef HAVE_OPENMP
            p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */
#endif
            if (bNN)
            {
#ifdef HAVE_NN_LOOPS /* Unlock this feature when testing */
                /* Loop over all atom pairs and estimate interaction energy */
#ifdef HAVE_OPENMP /* ------- */
#pragma omp single
#endif /* HAVE_OPENMP ------- */
                {
                    addFramesNN(hb, nframes);
                }
#ifdef HAVE_OPENMP /* ---------------- */
#pragma omp barrier
#pragma omp for schedule(dynamic)
#endif /* HAVE_OPENMP ---------------- */
                for (i=0; i<hb->d.nrd; i++)
                {
                    for(j=0;j<hb->a.nra; j++)
                    {
                        for (h=0;
                             h < (bContact ? 1 : hb->d.nhydro[i]);
                             h++)
                        {
                            if (i==hb->d.nrd || j==hb->a.nra)
                                gmx_fatal(FARGS, "out of bounds");

                            /* Get the real atom ids */
                            ii = hb->d.don[i];
                            jj = hb->a.acc[j];
                            hh = hb->d.hydro[i][h];
                    
                            /* Estimate the energy from the geometry */
                            E = calcHbEnergy(ii, jj, hh, x, NN, box, hbox, &(hb->d));
                            /* Store the energy */
                            storeHbEnergy(hb, i, j, h, E, nframes);
                        }
                    }
                }
#endif /* HAVE_NN_LOOPS */
            } /* if (bNN)*/
            else
            {
                if (bSelected)
                {
#ifdef HAVE_OPENMP
#pragma omp single
#endif
                    {
                        /* Do not parallelize this just yet. */
                        /* int ii; */
                        for(ii=0; (ii<nsel); ii++) {
                            int dd = index[0][i];
                            int aa = index[0][i+2];
                            /* int */ hh = index[0][i+1];
                            ihb = is_hbond(hb,ii,ii,dd,aa,rcut,r2cut,ccut,x,bBox,box,
                                           hbox,&dist,&ang,bDA,&h,bContact,bMerge,&peri);
          
                            if (ihb) {
                                /* add to index if not already there */
                                /* Add a hbond */
                                add_hbond(hb,dd,aa,hh,ii,ii,nframes,bMerge,ihb,bContact,peri);
                            }
                        }
                    } /* omp single */
                } /* if (bSelected) */
                else
                {
#ifdef HAVE_OPENMP
#pragma omp single
                    {
#endif
                    if (bGem)
                        calcBoxProjection(box, hb->per->P);

                    /* loop over all gridcells (xi,yi,zi)      */
                    /* Removed confusing macro, DvdS 27/12/98  */
#ifdef HAVE_OPENMP
                    }
                    /* The outer grid loop will have to do for now. */
#pragma omp for schedule(dynamic)
#endif
                    for(xi=0; (xi<ngrid[XX]); xi++)
                        for(yi=0; (yi<ngrid[YY]); yi++)
                            for(zi=0; (zi<ngrid[ZZ]); zi++) {
              
                                /* loop over donor groups gr0 (always) and gr1 (if necessary) */
                                for (grp=gr0; (grp <= (bTwo?gr1:gr0)); grp++) {
                                    icell=&(grid[zi][yi][xi].d[grp]);
                
                                    if (bTwo)
                                        ogrp = 1-grp;
                                    else
                                        ogrp = grp;
                
                                    /* loop over all hydrogen atoms from group (grp) 
                                     * in this gridcell (icell) 
                                     */
                                    for (ai=0; (ai<icell->nr); ai++) {
                                        i  = icell->atoms[ai];
                
                                        /* loop over all adjacent gridcells (xj,yj,zj) */
                                        /* This is a macro!!! */
                                        LOOPGRIDINNER(xj,yj,zj,xjj,yjj,zjj,xi,yi,zi,ngrid,bTric) {
                                            jcell=&(grid[zj][yj][xj].a[ogrp]);
                                            /* loop over acceptor atoms from other group (ogrp) 
                                             * in this adjacent gridcell (jcell) 
                                             */
                                            for (aj=0; (aj<jcell->nr); aj++) {
                                                j = jcell->atoms[aj];
                  
                                                /* check if this once was a h-bond */
                                                peri = -1;
                                                ihb = is_hbond(__HBDATA,grp,ogrp,i,j,rcut,r2cut,ccut,x,bBox,box,
                                                               hbox,&dist,&ang,bDA,&h,bContact,bMerge,&peri);
                    
                                                if (ihb) {
                                                    /* add to index if not already there */
                                                    /* Add a hbond */
                                                    add_hbond(__HBDATA,i,j,h,grp,ogrp,nframes,bMerge,ihb,bContact,peri);
                      
                                                    /* make angle and distance distributions */
                                                    if (ihb == hbHB && !bContact) {
                                                        if (dist>rcut)
                                                            gmx_fatal(FARGS,"distance is higher than what is allowed for an hbond: %f",dist);
                                                        ang*=RAD2DEG;
                                                        __ADIST[(int)( ang/abin)]++;
                                                        __RDIST[(int)(dist/rbin)]++;
                                                        if (!bTwo) {
                                                            int id,ia;
                                                            if ((id = donor_index(&hb->d,grp,i)) == NOTSET)
                                                                gmx_fatal(FARGS,"Invalid donor %d",i);
                                                            if ((ia = acceptor_index(&hb->a,ogrp,j)) == NOTSET)
                                                                gmx_fatal(FARGS,"Invalid acceptor %d",j);
                                                            resdist=abs(top.atoms.atom[i].resind-
                                                                        top.atoms.atom[j].resind);
                                                            if (resdist >= max_hx)
                                                                resdist = max_hx-1;
                                                            __HBDATA->nhx[nframes][resdist]++;
                                                        }
                                                    }

                                                }
                                            } /* for aj  */
                                        }
                                        ENDLOOPGRIDINNER;
                                    } /* for ai  */
                                } /* for grp */
                            } /* for xi,yi,zi */
                } /* if (bSelected) {...} else */ 

#ifdef HAVE_OPENMP /* ---------------------------- */
                /* Better wait for all threads to finnish using x[] before updating it. */
                k = nframes;            /*         */
#pragma omp barrier                     /*         */
#pragma omp critical                    /*         */
                {                       /*         */
                    /* Sum up histograms and counts from p_hb[] into hb */
                    {                   /*         */
                        hb->nhb[k]   += p_hb[threadNr]->nhb[k];
                        hb->ndist[k] += p_hb[threadNr]->ndist[k];
                        for (j=0; j<max_hx; j++) /**/
                            hb->nhx[k][j]  += p_hb[threadNr]->nhx[k][j];
                    }                   /*         */
                }                       /*         */
                /*                                 */
                /* Here are a handful of single constructs
                 * to share the workload a bit. The most
                 * important one is of course the last one,
                 * where there's a potential bottleneck in form
                 * of slow I/O.                    */
#pragma omp single /* ++++++++++++++++,            */
#endif /* HAVE_OPENMP ----------------+------------*/
                { /*                  +   */
                    if (hb != NULL)  /*   */
                    { /*              +   */
                        analyse_donor_props(opt2fn_null("-don",NFILE,fnm),hb,k,t,oenv);
                    } /*              +   */
                } /*                  +   */
#ifdef HAVE_OPENMP /*                 +   */
#pragma omp single /* +++           +++   */
#endif       /*                       +   */
                {  /*                 +   */
                    if (fpnhb)  /*    +   */
                        do_nhb_dist(fpnhb,hb,t);
                }  /*                 +   */
            } /* if (bNN) {...} else  +   */
#ifdef HAVE_OPENMP /*                 +   */
#pragma omp single /* +++           +++   */
#endif       /*                       +   */
            {      /*                 +   */
                trrStatus = (read_next_x(oenv,status,&t,natoms,x,box));
                nframes++;      /*    +   */
            }      /*                 +   */
#ifdef HAVE_OPENMP /* +++++++++++++++++   */
#pragma omp barrier
#endif
        } while (trrStatus);

#ifdef HAVE_OPENMP
#pragma omp critical
        {
            hb->nrhb += p_hb[threadNr]->nrhb;
            hb->nrdist += p_hb[threadNr]->nrdist;
        }
        /* Free parallel datastructures */
        sfree(p_hb[threadNr]->nhb);
        sfree(p_hb[threadNr]->ndist);
        sfree(p_hb[threadNr]->nhx);

#pragma omp for
        for (i=0; i<nabin; i++)
            for (j=0; j<actual_nThreads; j++)

                adist[i] += p_adist[j][i];
#pragma omp for
        for (i=0; i<=nrbin; i++)
            for (j=0; j<actual_nThreads; j++)
                rdist[i] += p_rdist[j][i];
    
        sfree(p_adist[threadNr]);
        sfree(p_rdist[threadNr]);
    } /* End of parallel region */
    sfree(p_adist);
    sfree(p_rdist);
#endif
  
    if(nframes <2 && (opt2bSet("-ac",NFILE,fnm) || opt2bSet("-life",NFILE,fnm)))
    {
        gmx_fatal(FARGS,"Cannot calculate autocorrelation of life times with less than two frames");
    }
  
    free_grid(ngrid,&grid);
  
    close_trj(status);
    if (fpnhb)
        ffclose(fpnhb);
    
    /* Compute maximum possible number of different hbonds */
    if (maxnhb > 0)
        max_nhb = maxnhb;
    else {
        max_nhb = 0.5*(hb->d.nrd*hb->a.nra);
    }
    /* Added support for -contact below.
     * - Erik Marklund, May 29-31, 2006 */
    /* Changed contact code.
     * - Erik Marklund, June 29, 2006 */
    if (bHBmap && !bNN) {
        if (hb->nrhb==0) {
            printf("No %s found!!\n", bContact ? "contacts" : "hydrogen bonds");
        } else {
            printf("Found %d different %s in trajectory\n"
                   "Found %d different atom-pairs within %s distance\n",
                   hb->nrhb, bContact?"contacts":"hydrogen bonds",
                   hb->nrdist,(r2cut>0)?"second cut-off":"hydrogen bonding");

            /*Control the pHist.*/

            if (bMerge)
                merge_hb(hb,bTwo,bContact);

            if (opt2bSet("-hbn",NFILE,fnm)) 
                dump_hbmap(hb,NFILE,fnm,bTwo,bContact,isize,index,grpnames,&top.atoms);

            /* Moved the call to merge_hb() to a line BEFORE dump_hbmap
             * to make the -hbn and -hmb output match eachother. 
             * - Erik Marklund, May 30, 2006 */
        }
    }
    /* Print out number of hbonds and distances */
    aver_nhb  = 0;    
    aver_dist = 0;
    fp = xvgropen(opt2fn("-num",NFILE,fnm),bContact ? "Contacts" :
                  "Hydrogen Bonds",output_env_get_xvgr_tlabel(oenv),"Number",oenv);
    snew(leg,2);
    snew(leg[0],STRLEN);
    snew(leg[1],STRLEN);
    sprintf(leg[0],"%s",bContact?"Contacts":"Hydrogen bonds");
    sprintf(leg[1],"Pairs within %g nm",(r2cut>0)?r2cut:rcut);
    xvgr_legend(fp,2,(const char**)leg,oenv);
    sfree(leg[1]);
    sfree(leg[0]);
    sfree(leg);
    for(i=0; (i<nframes); i++) {
        fprintf(fp,"%10g  %10d  %10d\n",hb->time[i],hb->nhb[i],hb->ndist[i]);
        aver_nhb  += hb->nhb[i];
        aver_dist += hb->ndist[i];
    }
    ffclose(fp);
    aver_nhb  /= nframes;
    aver_dist /= nframes;
    /* Print HB distance distribution */
    if (opt2bSet("-dist",NFILE,fnm)) {
        long sum;
    
        sum=0;
        for(i=0; i<nrbin; i++)
            sum+=rdist[i];
    
        fp = xvgropen(opt2fn("-dist",NFILE,fnm),
                      "Hydrogen Bond Distribution",
                      "Hydrogen - Acceptor Distance (nm)","",oenv);
        for(i=0; i<nrbin; i++)
            fprintf(fp,"%10g %10g\n",(i+0.5)*rbin,rdist[i]/(rbin*(real)sum));
        ffclose(fp);
    }
  
    /* Print HB angle distribution */
    if (opt2bSet("-ang",NFILE,fnm)) {
        long sum;
    
        sum=0;
        for(i=0; i<nabin; i++)
            sum+=adist[i];
    
        fp = xvgropen(opt2fn("-ang",NFILE,fnm),
                      "Hydrogen Bond Distribution",
                      "Donor - Hydrogen - Acceptor Angle (\\SO\\N)","",oenv);
        for(i=0; i<nabin; i++)
            fprintf(fp,"%10g %10g\n",(i+0.5)*abin,adist[i]/(abin*(real)sum));
        ffclose(fp);
    }
    
    /* Print HB in alpha-helix */
    if (opt2bSet("-hx",NFILE,fnm)) {
        fp = xvgropen(opt2fn("-hx",NFILE,fnm),
                      "Hydrogen Bonds",output_env_get_xvgr_tlabel(oenv),"Count",oenv);
        xvgr_legend(fp,NRHXTYPES,hxtypenames,oenv);
        for(i=0; i<nframes; i++) {
            fprintf(fp,"%10g",hb->time[i]);
            for (j=0; j<max_hx; j++)
                fprintf(fp," %6d",hb->nhx[i][j]);
            fprintf(fp,"\n");
        }
        ffclose(fp);
    }
    if (!bNN)
        printf("Average number of %s per timeframe %.3f out of %g possible\n",
               bContact ? "contacts" : "hbonds",
               bContact ? aver_dist : aver_nhb, max_nhb);
         
    /* Do Autocorrelation etc. */
    if (hb->bHBmap) {
        /* 
           Added support for -contact in ac and hbm calculations below.
           - Erik Marklund, May 29, 2006
        */
        ivec itmp;
        rvec rtmp;
        if (opt2bSet("-ac",NFILE,fnm) || opt2bSet("-life",NFILE,fnm))
            please_cite(stdout,"Spoel2006b");
        if (opt2bSet("-ac",NFILE,fnm)) {
            char *gemstring=NULL;

            if (bGem || bNN) {
                params = init_gemParams(rcut, D, hb->time, hb->nframes/2, nFitPoints, fit_start, fit_end,
                                        gemBallistic, nBalExp, bBallisticDt);
                if (params == NULL)
                    gmx_fatal(FARGS, "Could not initiate t_gemParams params.");
            }
            gemstring = strdup(gemType[hb->per->gemtype]);
            do_hbac(opt2fn("-ac",NFILE,fnm),hb,nDump,
                    bMerge,bContact,fit_start,temp,r2cut>0,smooth_tail_start,oenv,
                    params, gemstring, nThreads, NN, bBallistic, bGemFit);
        }
        if (opt2bSet("-life",NFILE,fnm))
            do_hblife(opt2fn("-life",NFILE,fnm),hb,bMerge,bContact,oenv);
        if (opt2bSet("-hbm",NFILE,fnm)) {
            t_matrix mat;
            int id,ia,hh,x,y;
      
            mat.nx=nframes;
            mat.ny=(bContact ? hb->nrdist : hb->nrhb);

            snew(mat.matrix,mat.nx);
            for(x=0; (x<mat.nx); x++) 
                snew(mat.matrix[x],mat.ny);
            y=0;
            for(id=0; (id<hb->d.nrd); id++) 
                for(ia=0; (ia<hb->a.nra); ia++) {
                    for(hh=0; (hh<hb->maxhydro); hh++) {
                        if (hb->hbmap[id][ia]) {
                            if (ISHB(hb->hbmap[id][ia]->history[hh])) {
                                /* Changed '<' into '<=' in the for-statement below.
                                 * It fixed the previously undiscovered bug that caused
                                 * the last occurance of an hbond/contact to not be
                                 * set in mat.matrix. Have a look at any old -hbm-output
                                 * and you will notice that the last column is allways empty.
                                 * - Erik Marklund May 30, 2006
                                 */
                                for(x=0; (x<=hb->hbmap[id][ia]->nframes); x++) {
                                    int nn0 = hb->hbmap[id][ia]->n0;
                                    range_check(y,0,mat.ny);
                                    mat.matrix[x+nn0][y] = is_hb(hb->hbmap[id][ia]->h[hh],x);
                                }
                                y++;
                            }
                        }
                    }
                }
            mat.axis_x=hb->time;
            snew(mat.axis_y,mat.ny);
            for(j=0; j<mat.ny; j++)
                mat.axis_y[j]=j;
            sprintf(mat.title,bContact ? "Contact Existence Map":
                    "Hydrogen Bond Existence Map");
            sprintf(mat.legend,bContact ? "Contacts" : "Hydrogen Bonds");
            sprintf(mat.label_x,"%s",output_env_get_xvgr_tlabel(oenv));
            sprintf(mat.label_y, bContact ? "Contact Index" : "Hydrogen Bond Index");
            mat.bDiscrete=TRUE;
            mat.nmap=2;
            snew(mat.map,mat.nmap);
            for(i=0; i<mat.nmap; i++) {
                mat.map[i].code.c1=hbmap[i];
                mat.map[i].desc=hbdesc[i];
                mat.map[i].rgb=hbrgb[i];
            }
            fp = opt2FILE("-hbm",NFILE,fnm,"w");
            write_xpm_m(fp, mat);
            ffclose(fp);
            for(x=0; x<mat.nx; x++)
                sfree(mat.matrix[x]);
            sfree(mat.axis_y);
            sfree(mat.matrix);
            sfree(mat.map);
        }
    }

    if (bGem) {
        fprintf(stderr, "There were %i periodic shifts\n", hb->per->nper);
        fprintf(stderr, "Freeing pHist for all donors...\n");
        for (i=0; i<hb->d.nrd; i++) {
            fprintf(stderr, "\r%i",i);
            if (hb->per->pHist[i] != NULL) {
                for (j=0; j<hb->a.nra; j++)
                    clearPshift(&(hb->per->pHist[i][j]));
                sfree(hb->per->pHist[i]);
            }
        }
        sfree(hb->per->pHist);
        sfree(hb->per->p2i);
        sfree(hb->per);
        fprintf(stderr, "...done.\n");
    }

#ifdef HAVE_NN_LOOPS
    if (bNN)
        free_hbEmap(hb);
#endif
    
    if (hb->bDAnr) {
        int  i,j,nleg;
        char **legnames;
        char buf[STRLEN];
    
#define USE_THIS_GROUP(j) ( (j == gr0) || (bTwo && (j == gr1)) )
    
        fp = xvgropen(opt2fn("-dan",NFILE,fnm),
                      "Donors and Acceptors",output_env_get_xvgr_tlabel(oenv),"Count",oenv);
        nleg = (bTwo?2:1)*2;
        snew(legnames,nleg);
        i=0;
        for(j=0; j<grNR; j++)
            if ( USE_THIS_GROUP(j) ) {
                sprintf(buf,"Donors %s",grpnames[j]);
                legnames[i++] = strdup(buf);
                sprintf(buf,"Acceptors %s",grpnames[j]);
                legnames[i++] = strdup(buf);
            }
        if (i != nleg)
            gmx_incons("number of legend entries");
        xvgr_legend(fp,nleg,(const char**)legnames,oenv);
        for(i=0; i<nframes; i++) {
            fprintf(fp,"%10g",hb->time[i]);
            for (j=0; (j<grNR); j++)
                if ( USE_THIS_GROUP(j) )
                    fprintf(fp," %6d",hb->danr[i][j]);
            fprintf(fp,"\n");
        }
        ffclose(fp);
    }
  
    thanx(stdout);
  
    return 0;
}