Project

General

Profile

gmx_hbond.c

Modified gmx_hbond.c - Justin Lemkul, 11/10/2010 04:42 PM

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

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

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

    
49
#ifdef DOUSEOPENMP
50
#define HAVE_OPENMP
51
#include "omp.h"
52
#endif
53

    
54
#include "statutil.h"
55
#include "copyrite.h"
56
#include "sysstuff.h"
57
#include "txtdump.h"
58
#include "futil.h"
59
#include "tpxio.h"
60
#include "physics.h"
61
#include "macros.h"
62
#include "gmx_fatal.h"
63
#include "index.h"
64
#include "smalloc.h"
65
#include "vec.h"
66
#include "xvgr.h"
67
#include "gstat.h"
68
#include "matio.h"
69
#include "string2.h"
70
#include "pbc.h"
71
#include "correl.h"
72
#include "gmx_ana.h"
73
#include "geminate.h"
74

    
75
typedef short int t_E;
76
typedef int t_EEst;
77
#define max_hx 7
78
typedef int t_hx[max_hx];
79
#define NRHXTYPES max_hx
80
const char *hxtypenames[NRHXTYPES]=
81
{"n-n","n-n+1","n-n+2","n-n+3","n-n+4","n-n+5","n-n>6"};
82
#define MAXHH 4
83

    
84
#ifdef HAVE_OPENMP
85
#define MASTER_THREAD_ONLY(threadNr) ((threadNr)==0)
86
#else
87
#define MASTER_THREAD_ONLY(threadNr) ((threadNr)==(threadNr))
88
#endif
89

    
90
/* -----------------------------------------*/
91

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

    
99
static gmx_bool bDebug = FALSE;
100

    
101
#define HB_NO 0
102
#define HB_YES 1<<0
103
#define HB_INS 1<<1
104
#define HB_YESINS HB_YES|HB_INS
105
#define HB_NR (1<<2)
106
#define MAXHYDRO 4
107

    
108
#define ISHB(h)   (((h) & 2) == 2)
109
#define ISDIST(h) (((h) & 1) == 1)
110
#define ISDIST2(h) (((h) & 4) == 4)
111
#define ISACC(h)  (((h) & 1) == 1)
112
#define ISDON(h)  (((h) & 2) == 2)
113
#define ISINGRP(h) (((h) & 4) == 4)
114

    
115
typedef struct {
116
    int nr;
117
    int maxnr;
118
    atom_id *atoms;
119
} t_ncell;
120

    
121
typedef struct {
122
    t_ncell d[grNR];
123
    t_ncell a[grNR];
124
} t_gridcell;
125

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

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

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

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

    
168
typedef struct {
169
    int len;     /* The length of frame and p. */
170
    int *frame;  /* The frames at which transitio*/
171
    PSTYPE *p;
172
} t_pShift;
173

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

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

    
196
typedef struct {
197
    gmx_bool        bHBmap,bDAnr,bGem;
198
    int         wordlen;
199
    /* The following arrays are nframes long */
200
    int         nframes,max_frames,maxhydro;
201
    int         *nhb,*ndist;
202
    h_id        *n_bound;
203
    real        *time;
204
    t_icell     *danr;
205
    t_hx        *nhx;
206
    /* These structures are initialized from the topology at start up */
207
    t_donors    d;
208
    t_acceptors a;
209
    /* This holds a matrix with all possible hydrogen bonds */
210
    int         nrhb,nrdist;
211
    t_hbond     ***hbmap;
212
#ifdef HAVE_NN_LOOPS
213
    t_hbEmap    hbE;
214
#endif
215
    /* For parallelization reasons this will have to be a pointer.
216
     * Otherwise discrepancies may arise between the periodicity data
217
     * seen by different threads. */
218
    t_gemPeriod *per;
219
} t_hbdata;
220

    
221
static void clearPshift(t_pShift *pShift)
222
{
223
    if (pShift->len > 0) {
224
        sfree(pShift->p);
225
        sfree(pShift->frame);
226
        pShift->len = 0;
227
    }
228
}
229

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

    
237
    for (i=0; i<3; i++){ m = vp[i];
238
        for (j=0; j<3; j++){ n = vp[j];
239
            U[m][n] = i==j ? 1:0;
240
        }
241
    }
242
    m_inv(B,M);
243
    for (i=0; i<3; i++){ m = vp[i];
244
        mvmul(M, U[m], P[m]);
245
    }
246
    transpose(P,N);
247
}
248

    
249
static void calcBoxDistance(matrix P, rvec d, ivec ibd){
250
    /* returns integer distance in box coordinates.
251
     * P is the projection matrix from cartesian coordinates
252
     * obtained with calcBoxProjection(). */
253
    int i;
254
    rvec bd;
255
    mvmul(P, d, bd);
256
    /* extend it by 0.5 in all directions since (int) rounds toward 0.*/
257
    for (i=0;i<3;i++)
258
        bd[i]=bd[i] + (bd[i]<0 ? -0.5 : 0.5);
259
    ibd[XX] = (int)bd[XX];
260
    ibd[YY] = (int)bd[YY];
261
    ibd[ZZ] = (int)bd[ZZ];
262
}
263

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

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

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

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

    
294
/* Unfortunately this needs to be critical it seems. */
295
#ifdef HAVE_OPENMP
296
#pragma omp critical
297
#endif
298
    {
299
        if (!per->p2i) {
300
            fprintf(stderr, "p2i not initialized. This shouldn't happen!\n");
301
            snew(per->p2i, 1);
302
        }
303
        else
304
            srenew(per->p2i, per->nper+2);
305
        copy_ivec(r, per->p2i[per->nper]);
306
        (per->nper)++;
307

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

    
317
static t_hbdata *mk_hbdata(gmx_bool bHBmap,gmx_bool bDAnr,gmx_bool oneHB, gmx_bool bGem, int gemmode)
318
{
319
    t_hbdata *hb;
320
  
321
    snew(hb,1);
322
    hb->wordlen = 8*sizeof(unsigned int);
323
    hb->bHBmap  = bHBmap;
324
    hb->bDAnr   = bDAnr;
325
    hb->bGem    = bGem;
326
    if (oneHB)
327
        hb->maxhydro = 1;
328
    else
329
        hb->maxhydro = MAXHYDRO;
330
    snew(hb->per, 1);
331
    hb->per->gemtype = bGem ? gemmode : 0;
332
  
333
    return hb;
334
}
335

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

    
340
    snew(hb->hbmap,hb->d.nrd);
341
    for(i=0; (i<hb->d.nrd); i++) {
342
        snew(hb->hbmap[i],hb->a.nra);
343
        if (hb->hbmap[i] == NULL)
344
            gmx_fatal(FARGS,"Could not allocate enough memory for hbmap");
345
        for (j=0; (j>hb->a.nra); j++)
346
            hb->hbmap[i][j] = NULL;
347
    }
348
}
349

    
350
/* Consider redoing pHist so that is only stores transitions between
351
 * periodicities and not the periodicity for all frames. This eats heaps of memory. */
352
static void mk_per(t_hbdata *hb)
353
{
354
    int i,j;
355
    if (hb->bGem) {
356
        snew(hb->per->pHist, hb->d.nrd);
357
        for (i=0; i<hb->d.nrd; i++) {
358
            snew(hb->per->pHist[i], hb->a.nra);
359
            if (hb->per->pHist[i]==NULL)
360
                gmx_fatal(FARGS,"Could not allocate enough memory for per->pHist");
361
            for (j=0; j<hb->a.nra; j++) {
362
                clearPshift(&(hb->per->pHist[i][j]));
363
            }
364
        }
365
        /* add the [0,0,0] shift to element 0 of p2i. */
366
        snew(hb->per->p2i, 1);
367
        clear_ivec(hb->per->p2i[0]);
368
        hb->per->nper = 1;
369
    }
370
}
371

    
372
#ifdef HAVE_NN_LOOPS
373
static void mk_hbEmap (t_hbdata *hb, int n0)
374
{
375
    int i, j, k;
376
    hb->hbE.E = NULL;
377
    hb->hbE.nframes = 0;
378
    snew(hb->hbE.E, hb->d.nrd);
379
    for (i=0; i<hb->d.nrd; i++)
380
    {
381
        snew(hb->hbE.E[i], hb->a.nra);
382
        for (j=0; j<hb->a.nra; j++)
383
        {
384
            snew(hb->hbE.E[i][j], MAXHYDRO);
385
            for (k=0; k<MAXHYDRO; k++)
386
                hb->hbE.E[i][j][k] = NULL;
387
        }
388
    }
389
    hb->hbE.Etot = NULL;
390
}
391

    
392
static void free_hbEmap (t_hbdata *hb)
393
{
394
    int i, j, k;
395
    for (i=0; i<hb->d.nrd; i++)
396
    {
397
        for (j=0; j<hb->a.nra; j++)
398
        {
399
            for (k=0; k<MAXHYDRO; k++)
400
                sfree(hb->hbE.E[i][j][k]);
401
            sfree(hb->hbE.E[i][j]);
402
        }
403
        sfree(hb->hbE.E[i]);
404
    }
405
    sfree(hb->hbE.E);
406
    sfree(hb->hbE.Etot);
407
}
408

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

    
412
#define DELTAFRAMES_HBE 10
413

    
414
    int d,a,h,nframes;
415

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

    
420
        for (d=0; d<hb->d.nrd; d++)
421
            for (a=0; a<hb->a.nra; a++)
422
                for (h=0; h<hb->d.nhydro[d]; h++)
423
                    srenew(hb->hbE.E[d][a][h], nframes);
424
    
425
        hb->hbE.nframes += DELTAFRAMES_HBE;
426
    }
427
}
428

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

    
441
    t_E E;
442
    rvec dist;
443
    rvec dipole[2], xmol[3], xmean[2]; 
444
    int i;
445
    real r, realE;
446

    
447
    if (d == a)
448
        /* Self-interaction */
449
        return NONSENSE_E;
450

    
451
    switch (EEst)
452
    {
453
    case NN_BINARY:
454
        /* This is a simple binary existence function that sets E=1 whenever
455
         * the distance between the oxygens is equal too or less than 0.35 nm.
456
         */
457
        rvec_sub(x[d], x[a], dist);
458
        pbc_correct_gem(dist, box, hbox);
459
        if (norm(dist) <= 0.35)
460
            E = 1;
461
        else
462
            E = 0;
463
        break;
464

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

    
473
        svmul(15.9994*(1/1.008), xmol[0], xmean[0]);
474
        rvec_inc(xmean[0], xmol[1]);
475
        rvec_inc(xmean[0], xmol[2]);
476
        for(i=0; i<3; i++)
477
            xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008;
478

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

    
484

    
485
        svmul(15.9994*(1/1.008), xmol[0], xmean[1]);
486
        rvec_inc(xmean[1], xmol[1]);
487
        rvec_inc(xmean[1], xmol[2]);
488
        for(i=0; i<3; i++)
489
            xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008;
490

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

    
495
        realE = pow(r, -3.0);
496
        E = (t_E)(SCALEFACTOR_E * realE);
497
        break;
498

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

    
509
        rvec_inc(dipole[0], xmol[1]);
510
        rvec_inc(dipole[0], xmol[2]);
511
        for (i=0; i<3; i++)
512
            dipole[0][i] *= 0.5;
513
        rvec_dec(dipole[0], xmol[0]);
514

    
515
        svmul(15.9994*(1/1.008), xmol[0], xmean[0]);
516
        rvec_inc(xmean[0], xmol[1]);
517
        rvec_inc(xmean[0], xmol[2]);
518
        for(i=0; i<3; i++)
519
            xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008;
520

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

    
526

    
527
        rvec_inc(dipole[1], xmol[1]);
528
        rvec_inc(dipole[1], xmol[2]);
529
        for (i=0; i<3; i++)
530
            dipole[1][i] *= 0.5;
531
        rvec_dec(dipole[1], xmol[0]);
532

    
533
        svmul(15.9994*(1/1.008), xmol[0], xmean[1]);
534
        rvec_inc(xmean[1], xmol[1]);
535
        rvec_inc(xmean[1], xmol[2]);
536
        for(i=0; i<3; i++)
537
            xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008;
538

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

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

    
553
    return E;
554
}
555

    
556
static void storeHbEnergy(t_hbdata *hb, int d, int a, int h, t_E E, int frame){
557
    /* hb - hbond data structure
558
       d  - donor
559
       a  - acceptor
560
       h  - hydrogen
561
       E  - estimate of the energy
562
       frame - the current frame.
563
    */
564

    
565
    /* Store the estimated energy */
566
    if (E == NONSENSE_E)
567
        E = 0;
568

    
569
    hb->hbE.E[d][a][h][frame] = E;
570
#ifdef HAVE_OPENMP
571
#pragma omp critical
572
#endif
573
    {
574
        hb->hbE.Etot[frame] += E;
575
    }
576
}
577
#endif /* HAVE_NN_LOOPS */
578

    
579

    
580
/* Finds -v[] in the periodicity index */
581
static int findMirror(PSTYPE p, ivec v[], PSTYPE nper)
582
{
583
    PSTYPE i;
584
    ivec u;
585
    for (i=0; i<nper; i++){
586
        if (v[i][XX] == -(v[p][XX]) &&
587
            v[i][YY] == -(v[p][YY]) &&
588
            v[i][ZZ] == -(v[p][ZZ]))
589
            return (int)i;
590
    }
591
    printf("Couldn't find mirror of [%i, %i, %i], index \n",
592
           v[p][XX],
593
           v[p][YY],
594
           v[p][ZZ]);
595
    return -1;
596
}
597
  
598

    
599
static void add_frames(t_hbdata *hb,int nframes)
600
{
601
    int  i,j,k,l;
602
  
603
    if (nframes >= hb->max_frames) {
604
        hb->max_frames += 4096;
605
        srenew(hb->time,hb->max_frames);
606
        srenew(hb->nhb,hb->max_frames);
607
        srenew(hb->ndist,hb->max_frames);
608
        srenew(hb->n_bound,hb->max_frames);
609
        srenew(hb->nhx,hb->max_frames);
610
        if (hb->bDAnr)
611
            srenew(hb->danr,hb->max_frames);
612
    }
613
    hb->nframes=nframes;
614
}
615

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

    
619
static void _set_hb(unsigned int hbexist[],unsigned int frame,gmx_bool bValue)
620
{
621
    if (bValue)
622
        hbexist[OFFSET(frame)] |= MASK(frame);
623
    else
624
        hbexist[OFFSET(frame)] &= ~MASK(frame);
625
}
626

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

    
632
static void set_hb(t_hbdata *hb,int id,int ih, int ia,int frame,int ihb)
633
{
634
    unsigned int *ghptr=NULL;
635
  
636
    if (ihb == hbHB)
637
        ghptr = hb->hbmap[id][ia]->h[ih];
638
    else if (ihb == hbDist)
639
        ghptr = hb->hbmap[id][ia]->g[ih];
640
    else
641
        gmx_fatal(FARGS,"Incomprehensible iValue %d in set_hb",ihb);
642

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

    
646
static void addPshift(t_pShift *pHist, PSTYPE p, int frame)
647
{
648
    if (pHist->len == 0) {
649
        snew(pHist->frame, 1);
650
        snew(pHist->p, 1);
651
        pHist->len      = 1;
652
        pHist->frame[0] = frame;
653
        pHist->p[0]     = p;
654
        return;
655
    } else
656
        if (pHist->p[pHist->len-1] != p) {
657
            pHist->len++;
658
            srenew(pHist->frame, pHist->len);
659
            srenew(pHist->p, pHist->len);
660
            pHist->frame[pHist->len-1] = frame;
661
            pHist->p[pHist->len-1]     = p;
662
        } /* Otherwise, there is no transition. */
663
    return;
664
}
665

    
666
static PSTYPE getPshift(t_pShift pHist, int frame)
667
{
668
    int f, i;
669

    
670
    if (pHist.len == 0
671
        || (pHist.len > 0 && pHist.frame[0]>frame))
672
        return -1;
673
  
674
    for (i=0; i<pHist.len; i++)
675
    {
676
        f = pHist.frame[i];
677
        if (f==frame)
678
            return pHist.p[i];
679
        if (f>frame)
680
            return pHist.p[i-1];
681
    }
682
  
683
    /* It seems that frame is after the last periodic transition. Return the last periodicity. */
684
    return pHist.p[pHist.len-1];
685
}
686

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

    
696
    if (!hb->h[0]) {
697
        hb->n0        = frame;
698
        hb->maxframes = delta;
699
        for(i=0; (i<maxhydro); i++) {
700
            snew(hb->h[i],hb->maxframes/wlen);
701
            snew(hb->g[i],hb->maxframes/wlen);
702
        }
703
    } else {
704
        hb->nframes = frame-hb->n0;
705
        /* We need a while loop here because hbonds may be returning
706
         * after a long time.
707
         */
708
        while (hb->nframes >= hb->maxframes) {
709
            n = hb->maxframes + delta;
710
            for(i=0; (i<maxhydro); i++) {
711
                srenew(hb->h[i],n/wlen);
712
                srenew(hb->g[i],n/wlen);
713
                for(j=hb->maxframes/wlen; (j<n/wlen); j++) {
714
                    hb->h[i][j] = 0;
715
                    hb->g[i][j] = 0;
716
                }
717
            }
718

    
719
            hb->maxframes = n;
720
        }
721
    }
722
    if (frame >= 0) {
723
        set_hb(hbd,id,h,ia,frame,ihb);
724
        if (bGem) {
725
            if (p>=hbd->per->nper)
726
                gmx_fatal(FARGS, "invalid shift: p=%u, nper=%u", p, hbd->per->nper);
727
            else
728
                addPshift(&(hbd->per->pHist[id][ia]), p, frame);
729
      
730
        }
731
    }
732
}
733

    
734
static void inc_nhbonds(t_donors *ddd,int d, int h)
735
{
736
    int j;
737
    int dptr = ddd->dptr[d];
738
  
739
    for(j=0; (j<ddd->nhydro[dptr]); j++)
740
        if (ddd->hydro[dptr][j] == h) {
741
            ddd->nhbonds[dptr][j]++;
742
            break;
743
        }
744
    if (j == ddd->nhydro[dptr])
745
        gmx_fatal(FARGS,"No such hydrogen %d on donor %d\n",h+1,d+1);
746
}
747

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

    
753
    if (a->grp[ai] != grp) {
754
        if (debug && bDebug) 
755
            fprintf(debug,"Acc. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
756
                    ai,a->grp[ai],grp,file,line);
757
        return NOTSET;
758
    }
759
    else
760
        return ai;
761
}
762
#define acceptor_index(a,grp,i) _acceptor_index(a,grp,i,__FILE__,__LINE__)
763

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

    
771
    if (d->grp[di] != grp) {
772
        if (debug && bDebug)
773
            fprintf(debug,"Don. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
774
                    di,d->grp[di],grp,file,line);
775
        return NOTSET;
776
    }
777
    else
778
        return di;
779
}
780
#define donor_index(d,grp,i) _donor_index(d,grp,i,__FILE__,__LINE__)
781

    
782
static gmx_bool isInterchangable(t_hbdata *hb, int d, int a, int grpa, int grpd)
783
{
784
    /* g_hbond doesn't allow overlapping groups */
785
    if (grpa!=grpd)
786
        return FALSE;
787
    return
788
        donor_index(&hb->d,grpd,a) != NOTSET
789
        && acceptor_index(&hb->a,grpa,d) != NOTSET;
790
}
791

    
792

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

    
799
    if ((id = hb->d.dptr[d]) == NOTSET)
800
        gmx_fatal(FARGS,"No donor atom %d",d+1);
801
    else if (grpd != hb->d.grp[id])
802
        gmx_fatal(FARGS,"Inconsistent donor groups, %d iso %d, atom %d",
803
                  grpd,hb->d.grp[id],d+1);
804
    if ((ia = hb->a.aptr[a]) == NOTSET)
805
        gmx_fatal(FARGS,"No acceptor atom %d",a+1);
806
    else if (grpa != hb->a.grp[ia])
807
        gmx_fatal(FARGS,"Inconsistent acceptor groups, %d iso %d, atom %d",
808
                  grpa,hb->a.grp[ia],a+1);
809

    
810
    if (bMerge)
811
        if ((daSwap = isInterchangable(hb, d, a, grpd, grpa) || bContact) && d>a)
812
            /* Then swap identity so that the id of d is lower then that of a.
813
             *
814
             * This should really be redundant by now, as is_hbond() now ought to return
815
             * hbNo in the cases where this conditional is TRUE. */
816
        {
817
            k = d;
818
            d = a;
819
            a = k;
820
    
821
            /* Now repeat donor/acc check. */
822
            if ((id = hb->d.dptr[d]) == NOTSET)
823
                gmx_fatal(FARGS,"No donor atom %d",d+1);
824
            else if (grpd != hb->d.grp[id])
825
                gmx_fatal(FARGS,"Inconsistent donor groups, %d iso %d, atom %d",
826
                          grpd,hb->d.grp[id],d+1);
827
            if ((ia = hb->a.aptr[a]) == NOTSET)
828
                gmx_fatal(FARGS,"No acceptor atom %d",a+1);
829
            else if (grpa != hb->a.grp[ia])
830
                gmx_fatal(FARGS,"Inconsistent acceptor groups, %d iso %d, atom %d",
831
                          grpa,hb->a.grp[ia],a+1);
832
        }
833

    
834
    if (hb->hbmap) {
835
        /* Loop over hydrogens to find which hydrogen is in this particular HB */
836
        if ((ihb == hbHB) && !bMerge && !bContact) {
837
            for(k=0; (k<hb->d.nhydro[id]); k++) 
838
                if (hb->d.hydro[id][k] == h)
839
                    break;
840
            if (k == hb->d.nhydro[id])
841
                gmx_fatal(FARGS,"Donor %d does not have hydrogen %d (a = %d)",
842
                          d+1,h+1,a+1);
843
        }
844
        else
845
            k = 0;
846
    
847
        if (hb->bHBmap) {
848
            if (hb->hbmap[id][ia] == NULL) {
849
                snew(hb->hbmap[id][ia],1);
850
                snew(hb->hbmap[id][ia]->h,hb->maxhydro);
851
                snew(hb->hbmap[id][ia]->g,hb->maxhydro);
852
            }
853
            add_ff(hb,id,k,ia,frame,ihb,p);
854
        }
855
    
856
        /* Strange construction with frame >=0 is a relic from old code
857
         * for selected hbond analysis. It may be necessary again if that
858
         * is made to work again.
859
         */
860
        if (frame >= 0) {
861
            hh = hb->hbmap[id][ia]->history[k];
862
            if (ihb == hbHB) {
863
                hb->nhb[frame]++;
864
                if (!(ISHB(hh))) {
865
                    hb->hbmap[id][ia]->history[k] = hh | 2;
866
                    hb->nrhb++;
867
                }
868
            }
869
            else
870
            {
871
                if (ihb == hbDist) {
872
                    hb->ndist[frame]++;
873
                    if (!(ISDIST(hh))) {
874
                        hb->hbmap[id][ia]->history[k] = hh | 1;
875
                        hb->nrdist++;
876
                    }
877
                }
878
            }
879
        }
880
    } else {
881
        if (frame >= 0) {
882
            if (ihb == hbHB) {
883
                hb->nhb[frame]++;
884
            } else {
885
                if (ihb == hbDist) {
886
                    hb->ndist[frame]++;
887
                }
888
            }
889
        }
890
    }
891
    if (bMerge && daSwap)
892
        h = hb->d.hydro[id][0];
893
    /* Increment number if HBonds per H */
894
    if (ihb == hbHB && !bContact)
895
        inc_nhbonds(&(hb->d),d,h);
896
}
897

    
898
/* Now a redundant function. It might find use at some point though. */
899
static gmx_bool in_list(atom_id selection,int isize,atom_id *index)
900
{
901
    int i;
902
    gmx_bool bFound;
903
  
904
    bFound=FALSE;
905
    for(i=0; (i<isize) && !bFound; i++)
906
        if(selection == index[i])
907
            bFound=TRUE;
908
  
909
    return bFound;
910
}
911

    
912
static char *mkatomname(t_atoms *atoms,int i)
913
{
914
    static char buf[32];
915
    int rnr;
916
  
917
    rnr = atoms->atom[i].resind;
918
    sprintf(buf,"%4s%d%-4s",
919
            *atoms->resinfo[rnr].name,atoms->resinfo[rnr].nr,*atoms->atomname[i]);
920
  
921
    return buf;
922
}
923

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

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

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

    
944
static void add_acc(t_acceptors *a,int ia,int grp)
945
{
946
    if (a->nra >= a->max_nra) {
947
        a->max_nra += 16;
948
        srenew(a->acc,a->max_nra);
949
        srenew(a->grp,a->max_nra);
950
    }
951
    a->grp[a->nra]   = grp;
952
    a->acc[a->nra++] = ia;
953
}
954

    
955
static void search_acceptors(t_topology *top,int isize, 
956
                             atom_id *index,t_acceptors *a,int grp,
957
                             gmx_bool bNitAcc,
958
                             gmx_bool bContact,gmx_bool bDoIt, unsigned char *datable)
959
{
960
    int i,n;
961
  
962
    if (bDoIt) {
963
        for (i=0; (i<isize); i++) {
964
            n = index[i];
965
            if ((bContact ||
966
                 (((*top->atoms.atomname[n])[0] == 'O') || 
967
                  (bNitAcc && ((*top->atoms.atomname[n])[0] == 'N')))) &&
968
                ISINGRP(datable[n])) {
969
                datable[n] |= ACC; /* set the atom's acceptor flag in datable. */
970
                add_acc(a,n,grp);
971
            }
972
        }
973
    }
974
    snew(a->aptr,top->atoms.nr);
975
    for(i=0; (i<top->atoms.nr); i++)
976
        a->aptr[i] = NOTSET;
977
    for(i=0; (i<a->nra); i++)
978
        a->aptr[a->acc[i]] = i;
979
}
980

    
981
static void add_h2d(int id,int ih,t_donors *ddd)
982
{
983
    int i;
984
  
985
    for(i=0; (i<ddd->nhydro[id]); i++) 
986
        if (ddd->hydro[id][i] == ih) {
987
            printf("Hm. This isn't the first time I found this donor (%d,%d)\n",
988
                   ddd->don[id],ih);
989
            break;
990
        }
991
    if (i == ddd->nhydro[id]) {
992
        if (ddd->nhydro[id] >= MAXHYDRO)
993
            gmx_fatal(FARGS,"Donor %d has more than %d hydrogens!",
994
                      ddd->don[id],MAXHYDRO);
995
        ddd->hydro[id][i] = ih;
996
        ddd->nhydro[id]++;
997
    }
998
}
999
  
1000
static void add_dh(t_donors *ddd,int id,int ih,int grp, unsigned char *datable)
1001
{
1002
    int i;
1003

    
1004
    if (ISDON(datable[id]) || !datable) {
1005
        if (ddd->dptr[id] == NOTSET) { /* New donor */
1006
            i = ddd->nrd;
1007
            ddd->dptr[id] = i;
1008
        } else 
1009
            i = ddd->dptr[id];
1010
  
1011
        if (i == ddd->nrd) {
1012
            if (ddd->nrd >= ddd->max_nrd) {
1013
                ddd->max_nrd += 128;
1014
                srenew(ddd->don,ddd->max_nrd);
1015
                srenew(ddd->nhydro,ddd->max_nrd);
1016
                srenew(ddd->hydro,ddd->max_nrd);
1017
                srenew(ddd->nhbonds,ddd->max_nrd);
1018
                srenew(ddd->grp,ddd->max_nrd);
1019
            }
1020
            ddd->don[ddd->nrd] = id;
1021
            ddd->nhydro[ddd->nrd] = 0;
1022
            ddd->grp[ddd->nrd] = grp;
1023
            ddd->nrd++;
1024
        } else
1025
            ddd->don[i] = id;
1026
        add_h2d(i,ih,ddd);
1027
    } else
1028
        if (datable)
1029
            printf("Warning: Atom %d is not in the d/a-table!\n", id);
1030
}
1031

    
1032
static void search_donors(t_topology *top, int isize, atom_id *index,
1033
                          t_donors *ddd,int grp,gmx_bool bContact,gmx_bool bDoIt,
1034
                          unsigned char *datable)
1035
{
1036
    int        i,j,nra,n;
1037
    t_functype func_type;
1038
    t_ilist    *interaction;
1039
    atom_id    nr1,nr2;
1040
    gmx_bool       stop;
1041

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

    
1048
    if (bContact) {
1049
        if (bDoIt)
1050
            for(i=0; (i<isize); i++) {
1051
                datable[index[i]] |= DON;
1052
                add_dh(ddd,index[i],-1,grp,datable);
1053
            }
1054
    }
1055
    else {
1056
        for(func_type=0; (func_type < F_NRE); func_type++) {
1057
            interaction=&(top->idef.il[func_type]);
1058
            if (func_type == F_POSRES)
1059
            { /* The ilist looks strange for posre. Bug in grompp?
1060
               * We don't need posre interactions for hbonds anyway.*/
1061
                continue;
1062
            }
1063
            for(i=0; i < interaction->nr; 
1064
                i+=interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1) {
1065
                /* next function */
1066
                if (func_type != top->idef.functype[interaction->iatoms[i]]) {
1067
                    fprintf(stderr,"Error in func_type %s",
1068
                            interaction_function[func_type].longname);
1069
                    continue;
1070
                }
1071
    
1072
                /* check out this functype */
1073
                if (func_type == F_SETTLE) {
1074
                    nr1=interaction->iatoms[i+1];
1075
      
1076
                    if (ISINGRP(datable[nr1])) {
1077
                        if (ISINGRP(datable[nr1+1])) {
1078
                            datable[nr1] |= DON;
1079
                            add_dh(ddd,nr1,nr1+1,grp,datable);
1080
                        }
1081
                        if (ISINGRP(datable[nr1+2])) {
1082
                            datable[nr1] |= DON;
1083
                            add_dh(ddd,nr1,nr1+2,grp,datable);
1084
                        }
1085
                    }
1086
                } 
1087
                else if (IS_CHEMBOND(func_type)) {
1088
                    for (j=0; j<2; j++) {
1089
                        nr1=interaction->iatoms[i+1+j];
1090
                        nr2=interaction->iatoms[i+2-j];
1091
                        if ((*top->atoms.atomname[nr1][0] == 'H') && 
1092
                            ((*top->atoms.atomname[nr2][0] == 'O') ||
1093
                             (*top->atoms.atomname[nr2][0] == 'N')) &&
1094
                            ISINGRP(datable[nr1]) && ISINGRP(datable[nr2])) {
1095
                            datable[nr2] |= DON;
1096
                            add_dh(ddd,nr2,nr1,grp,datable);
1097
                        }
1098
                    }
1099
                }
1100
            }
1101
        }
1102
#ifdef SAFEVSITES
1103
        for(func_type=0; func_type < F_NRE; func_type++) {
1104
            interaction=&top->idef.il[func_type];
1105
            for(i=0; i < interaction->nr; 
1106
                i+=interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1) {
1107
                /* next function */
1108
                if (func_type != top->idef.functype[interaction->iatoms[i]])
1109
                    gmx_incons("function type in search_donors");
1110
    
1111
                if ( interaction_function[func_type].flags & IF_VSITE ) {
1112
                    nr1=interaction->iatoms[i+1];
1113
                    if ( *top->atoms.atomname[nr1][0]  == 'H') {
1114
                        nr2=nr1-1;
1115
                        stop=FALSE;
1116
                        while (!stop && ( *top->atoms.atomname[nr2][0] == 'H'))
1117
                            if (nr2)
1118
                                nr2--;
1119
                            else
1120
                                stop=TRUE;
1121
                        if ( !stop && ( ( *top->atoms.atomname[nr2][0] == 'O') ||
1122
                                        ( *top->atoms.atomname[nr2][0] == 'N') ) &&
1123
                             ISINGRP(datable[nr1]) && ISINGRP(datable[nr2])) {
1124
                            datable[nr2] |= DON;
1125
                            add_dh(ddd,nr2,nr1,grp,datable);
1126
                        }
1127
                    }
1128
                }
1129
            }
1130
        }
1131
#endif
1132
    }
1133
}
1134

    
1135
static t_gridcell ***init_grid(gmx_bool bBox,rvec box[],real rcut,ivec ngrid)
1136
{
1137
    t_gridcell ***grid;
1138
    int i,y,z;
1139
  
1140
    if (bBox)
1141
        for(i=0; i<DIM; i++)
1142
            ngrid[i]=(box[i][i]/(1.2*rcut));
1143
  
1144
    if ( !bBox || (ngrid[XX]<3) || (ngrid[YY]<3) || (ngrid[ZZ]<3) )
1145
        for(i=0; i<DIM; i++)
1146
            ngrid[i]=1;
1147
    else 
1148
        printf("\nWill do grid-seach on %dx%dx%d grid, rcut=%g\n",
1149
               ngrid[XX],ngrid[YY],ngrid[ZZ],rcut);
1150
    snew(grid,ngrid[ZZ]);
1151
    for (z=0; z<ngrid[ZZ]; z++) {
1152
        snew((grid)[z],ngrid[YY]);
1153
        for (y=0; y<ngrid[YY]; y++)
1154
            snew((grid)[z][y],ngrid[XX]);
1155
    }
1156
    return grid;
1157
}
1158

    
1159
static void control_pHist(t_hbdata *hb, int nframes)
1160
{
1161
    int i,j,k;
1162
    PSTYPE p;
1163
    for (i=0;i<hb->d.nrd;i++)
1164
        for (j=0;j<hb->a.nra;j++)
1165
            if (hb->per->pHist[i][j].len != 0)
1166
                for (k=hb->hbmap[i][j][0].n0; k<nframes; k++) {
1167
                    p = getPshift(hb->per->pHist[i][j], k);
1168
                    if (p>hb->per->nper)
1169
                        fprintf(stderr, "Weird stuff in pHist[%i][%i].p at frame %i: p=%i\n",
1170
                                i,j,k,p);
1171
                }
1172
}
1173

    
1174
static void reset_nhbonds(t_donors *ddd)
1175
{
1176
    int i,j;
1177
  
1178
    for(i=0; (i<ddd->nrd); i++) 
1179
        for(j=0; (j<MAXHH); j++)
1180
            ddd->nhbonds[i][j] = 0;
1181
}
1182

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

    
1185
static void build_grid(t_hbdata *hb,rvec x[], rvec xshell,
1186
                       gmx_bool bBox, matrix box, rvec hbox,
1187
                       real rcut, real rshell,
1188
                       ivec ngrid, t_gridcell ***grid)
1189
{
1190
    int     i,m,gr,xi,yi,zi,nr;
1191
    atom_id *ad;
1192
    ivec    grididx;
1193
    rvec    invdelta,dshell,xtemp={0,0,0};
1194
    t_ncell *newgrid;
1195
    gmx_bool    bDoRshell,bInShell,bAcc;
1196
    real    rshell2=0;
1197
    int     gx,gy,gz;
1198
    int     dum = -1;
1199
  
1200
    bDoRshell = (rshell > 0);
1201
    rshell2   = sqr(rshell);
1202
    bInShell  = TRUE;
1203
  
1204
#define DBB(x) if (debug && bDebug) fprintf(debug,"build_grid, line %d, %s = %d\n",__LINE__,#x,x)
1205
    DBB(dum);
1206
    for(m=0; m<DIM; m++) {
1207
        hbox[m]=box[m][m]*0.5;
1208
        if (bBox) {
1209
            invdelta[m]=ngrid[m]/box[m][m];
1210
            if (1/invdelta[m] < rcut)
1211
                gmx_fatal(FARGS,"Your computational box has shrunk too much.\n"
1212
                          "%s can not handle this situation, sorry.\n",
1213
                          ShortProgram());
1214
        } else
1215
            invdelta[m]=0;
1216
    }
1217
    grididx[XX]=0;
1218
    grididx[YY]=0;
1219
    grididx[ZZ]=0;
1220
    DBB(dum);
1221
    /* resetting atom counts */
1222
    for(gr=0; (gr<grNR); gr++) {
1223
        for (zi=0; zi<ngrid[ZZ]; zi++)
1224
            for (yi=0; yi<ngrid[YY]; yi++)
1225
                for (xi=0; xi<ngrid[XX]; xi++) {
1226
                    grid[zi][yi][xi].d[gr].nr=0;
1227
                    grid[zi][yi][xi].a[gr].nr=0;
1228
                }
1229
        DBB(dum);
1230
    
1231
        /* put atoms in grid cells */
1232
        for(bAcc=FALSE; (bAcc<=TRUE); bAcc++) {
1233
            if (bAcc) {
1234
                nr = hb->a.nra;
1235
                ad = hb->a.acc;
1236
            }
1237
            else {
1238
                nr = hb->d.nrd;
1239
                ad = hb->d.don;
1240
            }
1241
            DBB(bAcc);
1242
            for(i=0; (i<nr); i++) {
1243
                /* check if we are inside the shell */
1244
                /* if bDoRshell=FALSE then bInShell=TRUE always */
1245
                DBB(i);
1246
                if ( bDoRshell ) {
1247
                    bInShell=TRUE;
1248
                    rvec_sub(x[ad[i]],xshell,dshell);
1249
                    if (bBox) {
1250
                        if (FALSE && !hb->bGem) {
1251
                            for(m=DIM-1; m>=0 && bInShell; m--) {
1252
                                if ( dshell[m] < -hbox[m] )
1253
                                    rvec_inc(dshell,box[m]);
1254
                                else if ( dshell[m] >= hbox[m] ) 
1255
                                    dshell[m] -= 2*hbox[m];
1256
                                /* if we're outside the cube, we're outside the sphere also! */
1257
                                if ( (dshell[m]>rshell) || (-dshell[m]>rshell) )
1258
                                    bInShell=FALSE;
1259
                            }
1260
                        } else {
1261
                            gmx_bool bDone = FALSE;
1262
                            while (!bDone)
1263
                            {
1264
                                bDone = TRUE;
1265
                                for(m=DIM-1; m>=0 && bInShell; m--) {
1266
                                    if ( dshell[m] < -hbox[m] ) {
1267
                                        bDone = FALSE;
1268
                                        rvec_inc(dshell,box[m]);
1269
                                    }
1270
                                    if ( dshell[m] >= hbox[m] ) {
1271
                                        bDone = FALSE;
1272
                                        dshell[m] -= 2*hbox[m];
1273
                                    }
1274
                                }
1275
                            }
1276
                            for(m=DIM-1; m>=0 && bInShell; m--) {
1277
                                /* if we're outside the cube, we're outside the sphere also! */
1278
                                if ( (dshell[m]>rshell) || (-dshell[m]>rshell) )
1279
                                    bInShell=FALSE;
1280
                            }
1281
                        }
1282
                    }
1283
                    /* if we're inside the cube, check if we're inside the sphere */
1284
                    if (bInShell)
1285
                        bInShell = norm2(dshell) < rshell2;
1286
                }
1287
                DBB(i);
1288
                if ( bInShell ) {
1289
                    if (bBox) {
1290
                        if (hb->bGem)
1291
                            copy_rvec(x[ad[i]], xtemp);
1292
                        pbc_correct_gem(x[ad[i]], box, hbox);
1293
                    }
1294
                    for(m=DIM-1; m>=0; m--) {
1295
                        if (TRUE || !hb->bGem){
1296
                            /* put atom in the box */
1297
                            while( x[ad[i]][m] < 0 )
1298
                                rvec_inc(x[ad[i]],box[m]);
1299
                            while( x[ad[i]][m] >= box[m][m] )
1300
                                rvec_dec(x[ad[i]],box[m]);
1301
                        }
1302
                        /* determine grid index of atom */
1303
                        grididx[m]=x[ad[i]][m]*invdelta[m];
1304
                        grididx[m] = (grididx[m]+ngrid[m]) % ngrid[m];
1305
                    }
1306
                    if (hb->bGem)
1307
                        copy_rvec(xtemp, x[ad[i]]); /* copy back */
1308
                    gx = grididx[XX];
1309
                    gy = grididx[YY];
1310
                    gz = grididx[ZZ];
1311
                    range_check(gx,0,ngrid[XX]);
1312
                    range_check(gy,0,ngrid[YY]);
1313
                    range_check(gz,0,ngrid[ZZ]);
1314
                    DBB(gx);
1315
                    DBB(gy);
1316
                    DBB(gz);
1317
                    /* add atom to grid cell */
1318
                    if (bAcc)
1319
                        newgrid=&(grid[gz][gy][gx].a[gr]);
1320
                    else
1321
                        newgrid=&(grid[gz][gy][gx].d[gr]);
1322
                    if (newgrid->nr >= newgrid->maxnr) {
1323
                        newgrid->maxnr+=10;
1324
                        DBB(newgrid->maxnr);
1325
                        srenew(newgrid->atoms, newgrid->maxnr);
1326
                    }
1327
                    DBB(newgrid->nr);
1328
                    newgrid->atoms[newgrid->nr]=ad[i];
1329
                    newgrid->nr++;
1330
                }
1331
            }
1332
        }
1333
    }
1334
}
1335

    
1336
static void count_da_grid(ivec ngrid, t_gridcell ***grid, t_icell danr)
1337
{
1338
    int gr,xi,yi,zi;
1339
  
1340
    for(gr=0; (gr<grNR); gr++) {
1341
        danr[gr]=0;
1342
        for (zi=0; zi<ngrid[ZZ]; zi++)
1343
            for (yi=0; yi<ngrid[YY]; yi++)
1344
                for (xi=0; xi<ngrid[XX]; xi++) {
1345
                    danr[gr] += grid[zi][yi][xi].d[gr].nr;
1346
                }
1347
    }
1348
}
1349

    
1350
/* The grid loop.
1351
 * Without a box, the grid is 1x1x1, so all loops are 1 long.
1352
 * With a rectangular box (bTric==FALSE) all loops are 3 long.
1353
 * With a triclinic box all loops are 3 long, except when a cell is
1354
 * located next to one of the box edges which is not parallel to the
1355
 * x/y-plane, in that case all cells in a line or layer are searched.
1356
 * This could be implemented slightly more efficient, but the code
1357
 * would get much more complicated.
1358
 */
1359
#define B(n,x,bTric,bEdge) ((n==1) ? x : bTric&&(bEdge) ? 0   : (x-1))
1360
#define E(n,x,bTric,bEdge) ((n==1) ? x : bTric&&(bEdge) ? n-1 : (x+1))
1361
#define GRIDMOD(j,n) (j+n)%(n)
1362
#define LOOPGRIDINNER(x,y,z,xx,yy,zz,xo,yo,zo,n,bTric)                  \
1363
    for(zz=B(n[ZZ],zo,bTric,FALSE); zz<=E(n[ZZ],zo,bTric,FALSE); zz++) { \
1364
    z=GRIDMOD(zz,n[ZZ]);                                                \
1365
    for(yy=B(n[YY],yo,bTric,z==0||z==n[ZZ]-1);                          \
1366
        yy<=E(n[YY],yo,bTric,z==0||z==n[ZZ]-1); yy++) {                 \
1367
    y=GRIDMOD(yy,n[YY]);                                                \
1368
    for(xx=B(n[XX],xo,bTric,y==0||y==n[YY]-1||z==0||z==n[ZZ]-1);        \
1369
        xx<=E(n[XX],xo,bTric,y==0||y==n[YY]-1||z==0||z==n[ZZ]-1); xx++) { \
1370
    x=GRIDMOD(xx,n[XX]);
1371
#define ENDLOOPGRIDINNER                                                \
1372
    }                                                                   \
1373
        }                                                               \
1374
                                        }                               \
1375
                                                                        \
1376

    
1377
static void dump_grid(FILE *fp, ivec ngrid, t_gridcell ***grid)
1378
{
1379
    int gr,x,y,z,sum[grNR];
1380
  
1381
    fprintf(fp,"grid %dx%dx%d\n",ngrid[XX],ngrid[YY],ngrid[ZZ]);
1382
    for (gr=0; gr<grNR; gr++) {
1383
        sum[gr]=0;
1384
        fprintf(fp,"GROUP %d (%s)\n",gr,grpnames[gr]);
1385
        for (z=0; z<ngrid[ZZ]; z+=2) {
1386
            fprintf(fp,"Z=%d,%d\n",z,z+1);
1387
            for (y=0; y<ngrid[YY]; y++) {
1388
                for (x=0; x<ngrid[XX]; x++) {
1389
                    fprintf(fp,"%3d",grid[x][y][z].d[gr].nr);
1390
                    sum[gr]+=grid[z][y][x].d[gr].nr;
1391
                    fprintf(fp,"%3d",grid[x][y][z].a[gr].nr);
1392
                    sum[gr]+=grid[z][y][x].a[gr].nr;
1393
      
1394
                }
1395
                fprintf(fp," | ");
1396
                if ( (z+1) < ngrid[ZZ] )
1397
                    for (x=0; x<ngrid[XX]; x++) {
1398
                        fprintf(fp,"%3d",grid[z+1][y][x].d[gr].nr);
1399
                        sum[gr]+=grid[z+1][y][x].d[gr].nr;
1400
                        fprintf(fp,"%3d",grid[z+1][y][x].a[gr].nr);
1401
                        sum[gr]+=grid[z+1][y][x].a[gr].nr;
1402
                    }
1403
                fprintf(fp,"\n");
1404
            }
1405
        }
1406
    }
1407
    fprintf(fp,"TOTALS:");
1408
    for (gr=0; gr<grNR; gr++)
1409
        fprintf(fp," %d=%d",gr,sum[gr]);
1410
    fprintf(fp,"\n");
1411
}
1412

    
1413
/* New GMX record! 5 * in a row. Congratulations! 
1414
 * Sorry, only four left.
1415
 */
1416
static void free_grid(ivec ngrid, t_gridcell ****grid)
1417
{
1418
    int y,z;
1419
    t_gridcell ***g = *grid;
1420
  
1421
    for (z=0; z<ngrid[ZZ]; z++) {
1422
        for (y=0; y<ngrid[YY]; y++) {
1423
            sfree(g[z][y]);
1424
        }
1425
        sfree(g[z]);
1426
    }
1427
    sfree(g);
1428
    g=NULL;
1429
}
1430

    
1431
static void pbc_correct(rvec dx,matrix box,rvec hbox)
1432
{
1433
    int m;
1434
    for(m=DIM-1; m>=0; m--) {
1435
        if ( dx[m] < -hbox[m] )
1436
            rvec_inc(dx,box[m]);
1437
        else if ( dx[m] >= hbox[m] )
1438
            rvec_dec(dx,box[m]);
1439
    }
1440
}
1441

    
1442
void pbc_correct_gem(rvec dx,matrix box,rvec hbox)
1443
{
1444
    int m;
1445
    gmx_bool bDone = FALSE;
1446
    while (!bDone) {
1447
        bDone = TRUE;
1448
        for(m=DIM-1; m>=0; m--) {
1449
            if ( dx[m] < -hbox[m] ) {
1450
                bDone = FALSE;
1451
                rvec_inc(dx,box[m]);
1452
            }
1453
            if ( dx[m] >= hbox[m] ) {
1454
                bDone = FALSE;
1455
                rvec_dec(dx,box[m]);
1456
            }
1457
        }
1458
    }
1459
}
1460

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

    
1478
    if (d == a)
1479
        return hbNo;
1480

    
1481
    if (((id = donor_index(&hb->d,grpd,d)) == NOTSET) ||
1482
        ((ja = acceptor_index(&hb->a,grpa,a)) == NOTSET))
1483
        return hbNo;
1484
  
1485
    rc2  = rcut*rcut;
1486
    r2c2 = r2cut*r2cut;
1487
  
1488
    rvec_sub(x[d],x[a],r_da);
1489
    /* Insert projection code here */
1490

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

    
1495
    if (bBox){
1496
        if (d>a && bMerge && (bContact || isInterchangable(hb, d, a, grpd, grpa))) { /* acceptor is also a donor and vice versa? */
1497
            return hbNo;
1498
            daSwap = TRUE; /* If so, then their history should be filed with donor and acceptor swapped. */
1499
        }
1500
        if (hb->bGem) {
1501
            copy_rvec(r_da, r); /* Save this for later */
1502
            pbc_correct_gem(r_da,box,hbox);
1503
        } else {
1504
            pbc_correct_gem(r_da,box,hbox);    
1505
        }
1506
    }
1507
    rda2 = iprod(r_da,r_da);
1508
  
1509
    if (bContact) {
1510
        if (daSwap)
1511
            return hbNo;
1512
        if (rda2 <= rc2){
1513
            if (hb->bGem){
1514
                calcBoxDistance(hb->per->P, r, ri);
1515
                *p = periodicIndex(ri, hb->per, daSwap);    /* find (or add) periodicity index. */
1516
            }
1517
            return hbHB;
1518
        }
1519
        else if (rda2 < r2c2)
1520
            return hbDist;
1521
        else
1522
            return hbNo;
1523
    }
1524
    *hhh = NOTSET;
1525
  
1526
    if (bDA && (rda2 > rc2))
1527
        return hbNo;
1528
  
1529
    for(h=0; (h < hb->d.nhydro[id]); h++) {
1530
        hh = hb->d.hydro[id][h];
1531
        rha2 = rc2+1;
1532
        if (!bDA) {
1533
            rvec_sub(x[hh],x[a],r_ha);
1534
            if (bBox) {
1535
                pbc_correct_gem(r_ha,box,hbox);
1536
            }
1537
            rha2 = iprod(r_ha,r_ha);
1538
        }
1539

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

    
1545
        if (bDA || (!bDA && (rha2 <= rc2))) {
1546
            rvec_sub(x[d],x[hh],r_dh);
1547
            if (bBox) {
1548
                if (hb->bGem)
1549
                    pbc_correct_gem(r_dh,box,hbox);
1550
                else
1551
                    pbc_correct_gem(r_dh,box,hbox);
1552
            }
1553
    
1554
            if (!bDA)
1555
                HAinrange = TRUE;
1556
            ca = cos_angle(r_dh,r_da);
1557
            /* if angle is smaller, cos is larger */
1558
            if (ca >= ccut) {
1559
                *hhh  = hh;
1560
                *d_ha = sqrt(bDA?rda2:rha2);
1561
                *ang  = acos(ca);
1562
                return hbHB;
1563
            }
1564
        }
1565
    }
1566
    if (bDA || (!bDA && HAinrange))
1567
        return hbDist;
1568
    else
1569
        return hbNo;
1570
}
1571

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

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

    
1591
    int m,mm,n00,n01,nn0,nnframes;
1592
    PSTYPE pm;
1593
    t_pShift *pShift;
1594

    
1595
    /* Decide where to start from when merging */
1596
    n00      = hb0->n0;
1597
    n01      = hb1->n0;
1598
    nn0      = min(n00,n01);
1599
    nnframes = max(n00 + hb0->nframes,n01 + hb1->nframes) - nn0;
1600
    /* Initiate tmp arrays */
1601
    for(m=0; (m<ntmp); m++) {
1602
        htmp[m] = 0;
1603
        gtmp[m] = 0;
1604
        ptmp[m] = 0;
1605
    }
1606
    /* Fill tmp arrays with values due to first HB */
1607
    /* Once again '<' had to be replaced with '<='
1608
       to catch the last frame in which the hbond
1609
       appears.
1610
       - Erik Marklund, June 1, 2006 */  
1611
    for(m=0; (m<=hb0->nframes); m++) {
1612
        mm = m+n00-nn0;
1613
        htmp[mm] = is_hb(hb0->h[0],m);
1614
        if (hb->bGem) {
1615
            pm = getPshift(hb->per->pHist[a1][a2], m+hb0->n0);
1616
            if (pm > hb->per->nper)
1617
                gmx_fatal(FARGS, "Illegal shift!");
1618
            else
1619
                ptmp[mm] = pm; /*hb->per->pHist[a1][a2][m];*/
1620
        }
1621
    }
1622
    /* If we're doing geminate recompbination we usually don't need the distances.
1623
     * Let's save some memory and time. */
1624
    if (TRUE || !hb->bGem || hb->per->gemtype == gemAD){
1625
        for(m=0; (m<=hb0->nframes); m++) {
1626
            mm = m+n00-nn0;
1627
            gtmp[mm] = is_hb(hb0->g[0],m);
1628
        }
1629
    }
1630
    /* Next HB */
1631
    for(m=0; (m<=hb1->nframes); m++) {
1632
        mm = m+n01-nn0;
1633
        htmp[mm] = htmp[mm] || is_hb(hb1->h[0],m);
1634
        gtmp[mm] = gtmp[mm] || is_hb(hb1->g[0],m);
1635
        if (hb->bGem /* && ptmp[mm] != 0 */) {
1636

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

    
1654
    /* Copy temp array to target array */
1655
    for(m=0; (m<=nnframes); m++) {
1656
        _set_hb(hb0->h[0],m,htmp[m]);
1657
        _set_hb(hb0->g[0],m,gtmp[m]);
1658
        if (hb->bGem)
1659
            addPshift(&(hb->per->pHist[a1][a2]), ptmp[m], m+nn0);
1660
    }
1661
  
1662
    /* Set scalar variables */
1663
    hb0->n0        = nn0;
1664
    hb0->maxframes = nnframes;
1665
}
1666

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

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

    
1682
    ntmp = 2*hb->max_frames;
1683
    snew(gtmp,ntmp);
1684
    snew(htmp,ntmp);
1685
    snew(ptmp,ntmp);
1686
    for(i=0; (i<hb->d.nrd); i++) {
1687
        fprintf(stderr,"\r%d/%d",i+1,hb->d.nrd);
1688
        id = hb->d.don[i];
1689
        ii = hb->a.aptr[id];
1690
        for(j=0; (j<hb->a.nra); j++) {
1691
            ia = hb->a.acc[j];
1692
            jj = hb->d.dptr[ia];
1693
            if ((id != ia) && (ii != NOTSET) && (jj != NOTSET) &&
1694
                (!bTwo || (bTwo && (hb->d.grp[i] != hb->a.grp[j])))) {
1695
                hb0 = hb->hbmap[i][j];
1696
                hb1 = hb->hbmap[jj][ii];
1697
                if (hb0 && hb1 && ISHB(hb0->history[0]) && ISHB(hb1->history[0])) {
1698
                    do_merge(hb,ntmp,htmp,gtmp,ptmp,hb0,hb1,i,j);
1699
                    if (ISHB(hb1->history[0])) 
1700
                        inrnew--;
1701
                    else if (ISDIST(hb1->history[0])) 
1702
                        indnew--;
1703
                    else
1704
                        if (bContact) 
1705
                            gmx_incons("No contact history");
1706
                        else
1707
                            gmx_incons("Neither hydrogen bond nor distance");
1708
                    sfree(hb1->h[0]);
1709
                    sfree(hb1->g[0]);
1710
                    if (hb->bGem) {
1711
                        clearPshift(&(hb->per->pHist[jj][ii]));
1712
                    }
1713
                    hb1->h[0] = NULL;
1714
                    hb1->g[0] = NULL;
1715
                    hb1->history[0] = hbNo;
1716
                }
1717
            }
1718
        }
1719
    }
1720
    fprintf(stderr,"\n");
1721
    printf("- Reduced number of hbonds from %d to %d\n",hb->nrhb,inrnew);
1722
    printf("- Reduced number of distances from %d to %d\n",hb->nrdist,indnew);
1723
    hb->nrhb   = inrnew;
1724
    hb->nrdist = indnew;
1725
    sfree(gtmp);
1726
    sfree(htmp);
1727
    sfree(ptmp);
1728
}
1729

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

    
1735
  
1736
    /* Set array to 0 */
1737
    for(k=0; (k<MAXHH); k++)
1738
        n_bound[k] = 0;
1739
    /* Loop over possible donors */
1740
    for(i=0; (i<hb->d.nrd); i++) {
1741
        for(j=0; (j<hb->d.nhydro[i]); j++)
1742
            n_bound[hb->d.nhbonds[i][j]]++;
1743
    }      
1744
    fprintf(fp,"%12.5e",t);
1745
    nbtot = 0;
1746
    for(k=0; (k<MAXHH); k++) {
1747
        fprintf(fp,"  %8d",n_bound[k]);
1748
        nbtot += n_bound[k]*k;
1749
    }
1750
    fprintf(fp,"  %8d\n",nbtot);
1751
}
1752

    
1753
/* Added argument bContact in do_hblife(...). Also
1754
 * added support for -contact in function body.
1755
 * - Erik Marklund, May 31, 2006 */
1756
/* Changed the contact code slightly.
1757
 * - Erik Marklund, June 29, 2006
1758
 */
1759
static void do_hblife(const char *fn,t_hbdata *hb,gmx_bool bMerge,gmx_bool bContact,
1760
                      const output_env_t oenv)
1761
{
1762
    FILE *fp;
1763
    const char *leg[] = { "p(t)", "t p(t)" };
1764
    int  *histo;
1765
    int  i,j,j0,k,m,nh,ihb,ohb,nhydro,ndump=0;
1766
    int   nframes = hb->nframes;
1767
    unsigned int **h;
1768
    real   t,x1,dt;
1769
    double sum,integral;
1770
    t_hbond *hbh;
1771
  
1772
    snew(h,hb->maxhydro);
1773
    snew(histo,nframes+1);
1774
    /* Total number of hbonds analyzed here */
1775
    for(i=0; (i<hb->d.nrd); i++) {
1776
        for(k=0; (k<hb->a.nra); k++) {
1777
            hbh = hb->hbmap[i][k];
1778
            if (hbh) {
1779
                if (bMerge) {
1780
                    if (hbh->h[0]) {
1781
                        h[0] = hbh->h[0];
1782
                        nhydro = 1;
1783
                    }
1784
                    else
1785
                        nhydro = 0;
1786
                }
1787
                else {
1788
                    nhydro = 0;
1789
                    for(m=0; (m<hb->maxhydro); m++)
1790
                        if (hbh->h[m]) {
1791
                            h[nhydro++] = bContact ? hbh->g[m] : hbh->h[m];
1792
                        }
1793
                }
1794
                for(nh=0; (nh<nhydro); nh++) {
1795
                    ohb = 0;
1796
                    j0  = 0;
1797

    
1798
                    /* Changed '<' into '<=' below, just like I
1799
                       did in the hbm-output-loop in the main code.
1800
                       - Erik Marklund, May 31, 2006
1801
                    */
1802
                    for(j=0; (j<=hbh->nframes); j++) {
1803
                        ihb      = is_hb(h[nh],j);
1804
                        if (debug && (ndump < 10))
1805
                            fprintf(debug,"%5d  %5d\n",j,ihb);
1806
                        if (ihb != ohb) {
1807
                            if (ihb) {
1808
                                j0 = j;
1809
                            }
1810
                            else {
1811
                                histo[j-j0]++;
1812
                            }
1813
                            ohb = ihb;
1814
                        }
1815
                    }
1816
                    ndump++;
1817
                }
1818
            }
1819
        }
1820
    }
1821
    fprintf(stderr,"\n");
1822
    if (bContact)
1823
        fp = xvgropen(fn,"Uninterrupted contact lifetime",output_env_get_xvgr_tlabel(oenv),"()",oenv);
1824
    else
1825
        fp = xvgropen(fn,"Uninterrupted hydrogen bond lifetime",output_env_get_xvgr_tlabel(oenv),"()",
1826
                      oenv);
1827

    
1828
    xvgr_legend(fp,asize(leg),leg,oenv);
1829
    j0 = nframes-1;
1830
    while ((j0 > 0) && (histo[j0] == 0))
1831
        j0--;
1832
    sum = 0;
1833
    for(i=0; (i<=j0); i++)
1834
        sum+=histo[i];
1835
    dt       = hb->time[1]-hb->time[0];
1836
    sum      = dt*sum;
1837
    integral = 0;
1838
    for(i=1; (i<=j0); i++) {
1839
        t  = hb->time[i] - hb->time[0] - 0.5*dt;
1840
        x1 = t*histo[i]/sum;
1841
        fprintf(fp,"%8.3f  %10.3e  %10.3e\n",t,histo[i]/sum,x1);
1842
        integral += x1;
1843
    }
1844
    integral *= dt;
1845
    ffclose(fp);
1846
    printf("%s lifetime = %.2f ps\n", bContact?"Contact":"HB", integral);
1847
    printf("Note that the lifetime obtained in this manner is close to useless\n");
1848
    printf("Use the -ac option instead and check the Forward lifetime\n");
1849
    please_cite(stdout,"Spoel2006b");
1850
    sfree(h);
1851
    sfree(histo);
1852
}
1853

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

    
1865
    if (nDump <= 0)
1866
        return;
1867
    fp = ffopen("debug-ac.xvg","w");
1868
    for(j=0; (j<nframes); j++) {
1869
        fprintf(fp,"%10.3f",hb->time[j]);
1870
        for(i=nd=0; (i<hb->d.nrd) && (nd < nDump); i++) {
1871
            for(k=0; (k<hb->a.nra) && (nd < nDump); k++) {
1872
                bPrint = FALSE;
1873
                ihb = idist = 0;
1874
                hbh = hb->hbmap[i][k];
1875
                if (oneHB) {
1876
                    if (hbh->h[0]) {
1877
                        ihb   = is_hb(hbh->h[0],j);
1878
                        idist = is_hb(hbh->g[0],j);
1879
                        bPrint = TRUE;
1880
                    }
1881
                } 
1882
                else {
1883
                    for(m=0; (m<hb->maxhydro) && !ihb ; m++) {
1884
                        ihb   = ihb   || ((hbh->h[m]) && is_hb(hbh->h[m],j));
1885
                        idist = idist || ((hbh->g[m]) && is_hb(hbh->g[m],j));
1886
                    }
1887
                    /* This is not correct! */
1888
                    /* What isn't correct? -Erik M */
1889
                    bPrint = TRUE;
1890
                }
1891
                if (bPrint) {
1892
                    fprintf(fp,"  %1d-%1d",ihb,idist);
1893
                    nd++;
1894
                }
1895
            }
1896
        }
1897
        fprintf(fp,"\n");
1898
    }
1899
    ffclose(fp);
1900
}
1901

    
1902
static real calc_dg(real tau,real temp)
1903
{
1904
    real kbt;
1905
  
1906
    kbt = BOLTZ*temp;
1907
    if (tau <= 0)
1908
        return -666;
1909
    else
1910
        return kbt*log(kbt*tau/PLANCK);  
1911
}
1912

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

    
1919
#ifdef HAVE_LIBGSL
1920
#include <gsl/gsl_multimin.h>
1921
#include <gsl/gsl_sf.h>
1922
#include <gsl/gsl_version.h>
1923

    
1924
static double my_f(const gsl_vector *v,void *params)
1925
{
1926
    t_luzar *tl = (t_luzar *)params;
1927
    int    i;
1928
    double tol=1e-16,chi2=0;
1929
    double di;
1930
    real   k,kp;
1931
  
1932
    for(i=0; (i<tl->nparams); i++) {
1933
        tl->kkk[i] = gsl_vector_get(v, i);
1934
    }
1935
    k  = tl->kkk[0];
1936
    kp = tl->kkk[1];
1937
  
1938
    for(i=tl->n0; (i<tl->n1); i+=tl->ndelta) {
1939
        di=sqr(k*tl->sigma_ct[i]) + sqr(kp*tl->sigma_nt[i]) + sqr(tl->sigma_kt[i]);
1940
        /*di = 1;*/
1941
        if (di > tol)
1942
            chi2 += sqr(k*tl->ct[i]-kp*tl->nt[i]-tl->kt[i])/di;
1943
      
1944
        else
1945
            fprintf(stderr,"WARNING: sigma_ct = %g, sigma_nt = %g, sigma_kt = %g\n"
1946
                    "di = %g k = %g kp = %g\n",tl->sigma_ct[i],
1947
                    tl->sigma_nt[i],tl->sigma_kt[i],di,k,kp);
1948
    }
1949
#ifdef DEBUG
1950
    chi2 = 0.3*sqr(k-0.6)+0.7*sqr(kp-1.3);
1951
#endif
1952
    return chi2;
1953
}
1954

    
1955
static real optimize_luzar_parameters(FILE *fp,t_luzar *tl,int maxiter,
1956
                                      real tol)
1957
{
1958
    real   size,d2;
1959
    int    iter   = 0;
1960
    int    status = 0;
1961
    int    i;
1962

    
1963
    const gsl_multimin_fminimizer_type *T;
1964
    gsl_multimin_fminimizer *s;
1965

    
1966
    gsl_vector *x,*dx;
1967
    gsl_multimin_function my_func;
1968

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

    
1973
    /* Starting point */
1974
    x = gsl_vector_alloc (my_func.n);
1975
    for(i=0; (i<my_func.n); i++)
1976
        gsl_vector_set (x, i, tl->kkk[i]);
1977
  
1978
    /* Step size, different for each of the parameters */
1979
    dx = gsl_vector_alloc (my_func.n);
1980
    for(i=0; (i<my_func.n); i++)
1981
        gsl_vector_set (dx, i, 0.01*tl->kkk[i]);
1982

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

    
1986
    gsl_multimin_fminimizer_set (s, &my_func, x, dx);
1987
    gsl_vector_free (x);
1988
    gsl_vector_free (dx);
1989

    
1990
    if (fp)
1991
        fprintf(fp,"%5s %12s %12s %12s %12s\n","Iter","k","kp","NM Size","Chi2");
1992
  
1993
    do  {
1994
        iter++;
1995
        status = gsl_multimin_fminimizer_iterate (s);
1996
    
1997
        if (status != 0)
1998
            gmx_fatal(FARGS,"Something went wrong in the iteration in minimizer %s",
1999
                      gsl_multimin_fminimizer_name(s));
2000
    
2001
        d2     = gsl_multimin_fminimizer_minimum(s);
2002
        size   = gsl_multimin_fminimizer_size(s);
2003
        status = gsl_multimin_test_size(size,tol);
2004
    
2005
        if (status == GSL_SUCCESS)
2006
            if (fp) 
2007
                fprintf(fp,"Minimum found using %s at:\n",
2008
                        gsl_multimin_fminimizer_name(s));
2009
  
2010
        if (fp) {
2011
            fprintf(fp,"%5d", iter);
2012
            for(i=0; (i<my_func.n); i++) 
2013
                fprintf(fp," %12.4e",gsl_vector_get (s->x,i));
2014
            fprintf (fp," %12.4e %12.4e\n",size,d2);
2015
        }
2016
    }
2017
    while ((status == GSL_CONTINUE) && (iter < maxiter));
2018
  
2019
    gsl_multimin_fminimizer_free (s);
2020
  
2021
    return d2;
2022
}
2023

    
2024
static real quality_of_fit(real chi2,int N)
2025
{
2026
    return gsl_sf_gamma_inc_Q((N-2)/2.0,chi2/2.0);
2027
}
2028

    
2029
#else
2030
static real optimize_luzar_parameters(FILE *fp,t_luzar *tl,int maxiter,
2031
                                      real tol)
2032
{
2033
    fprintf(stderr,"This program needs the GNU scientific library to work.\n");
2034
  
2035
    return -1;
2036
}
2037
static real quality_of_fit(real chi2,int N)
2038
{
2039
    fprintf(stderr,"This program needs the GNU scientific library to work.\n");
2040
  
2041
    return -1;
2042
}
2043

    
2044
#endif
2045

    
2046
static real compute_weighted_rates(int n,real t[],real ct[],real nt[],
2047
                                   real kt[],real sigma_ct[],real sigma_nt[],
2048
                                   real sigma_kt[],real *k,real *kp,
2049
                                   real *sigma_k,real *sigma_kp,
2050
                                   real fit_start)
2051
{
2052
#define NK 10
2053
    int      i,j;
2054
    t_luzar  tl;
2055
    real     kkk=0,kkp=0,kk2=0,kp2=0,chi2;
2056
  
2057
    *sigma_k  = 0;
2058
    *sigma_kp = 0;
2059
  
2060
    for(i=0; (i<n); i++) {
2061
        if (t[i] >= fit_start)
2062
            break;
2063
    }
2064
    tl.n0      = i;
2065
    tl.n1      = n;
2066
    tl.nparams = 2;
2067
    tl.ndelta  = 1;
2068
    tl.t  = t;
2069
    tl.ct = ct;
2070
    tl.nt = nt;
2071
    tl.kt = kt;
2072
    tl.sigma_ct = sigma_ct;
2073
    tl.sigma_nt = sigma_nt;
2074
    tl.sigma_kt = sigma_kt;
2075
    tl.kkk[0] = *k;
2076
    tl.kkk[1] = *kp;
2077
  
2078
    chi2 = optimize_luzar_parameters(debug,&tl,1000,1e-3);
2079
    *k  = tl.kkk[0];
2080
    *kp = tl.kkk[1] = *kp;
2081
    tl.ndelta = NK;
2082
    for(j=0; (j<NK); j++) {
2083
        (void) optimize_luzar_parameters(debug,&tl,1000,1e-3);
2084
        kkk += tl.kkk[0];
2085
        kkp += tl.kkk[1];
2086
        kk2 += sqr(tl.kkk[0]);
2087
        kp2 += sqr(tl.kkk[1]);
2088
        tl.n0++;
2089
    }
2090
    *sigma_k  = sqrt(kk2/NK - sqr(kkk/NK));
2091
    *sigma_kp = sqrt(kp2/NK - sqr(kkp/NK));
2092
  
2093
    return chi2;
2094
}
2095

    
2096
static void smooth_tail(int n,real t[],real c[],real sigma_c[],real start,
2097
                        const output_env_t oenv)
2098
{
2099
    FILE *fp;
2100
    real e_1,fitparm[4];
2101
    int  i;
2102
  
2103
    e_1 = exp(-1);
2104
    for(i=0; (i<n); i++)
2105
        if (c[i] < e_1)
2106
            break;
2107
    if (i < n)
2108
        fitparm[0] = t[i];
2109
    else
2110
        fitparm[0] = 10;
2111
    fitparm[1] = 0.95;
2112
    do_lmfit(n,c,sigma_c,0,t,start,t[n-1],oenv,bDebugMode(),effnEXP2,fitparm,0);
2113
}
2114

    
2115
void analyse_corr(int n,real t[],real ct[],real nt[],real kt[],
2116
                  real sigma_ct[],real sigma_nt[],real sigma_kt[],
2117
                  real fit_start,real temp,real smooth_tail_start,
2118
                  const output_env_t oenv)
2119
{
2120
    int    i0,i;
2121
    real   k=1,kp=1,kow=1;
2122
    real   Q=0,chi22,chi2,dg,dgp,tau_hb,dtau,tau_rlx,e_1,dt,sigma_k,sigma_kp,ddg;
2123
    double tmp,sn2=0,sc2=0,sk2=0,scn=0,sck=0,snk=0;
2124
    gmx_bool   bError = (sigma_ct != NULL) && (sigma_nt != NULL) && (sigma_kt != NULL);
2125
  
2126
    if (smooth_tail_start >= 0) {
2127
        smooth_tail(n,t,ct,sigma_ct,smooth_tail_start,oenv);
2128
        smooth_tail(n,t,nt,sigma_nt,smooth_tail_start,oenv);
2129
        smooth_tail(n,t,kt,sigma_kt,smooth_tail_start,oenv);
2130
    }
2131
    for(i0=0; (i0<n-2) && ((t[i0]-t[0]) < fit_start); i0++)
2132
        ;
2133
    if (i0 < n-2) { 
2134
        for(i=i0; (i<n); i++) {
2135
            sc2 += sqr(ct[i]);
2136
            sn2 += sqr(nt[i]);
2137
            sk2 += sqr(kt[i]);
2138
            sck += ct[i]*kt[i];
2139
            snk += nt[i]*kt[i];
2140
            scn += ct[i]*nt[i];
2141
        }
2142
        printf("Hydrogen bond thermodynamics at T = %g K\n",temp);
2143
        tmp = (sn2*sc2-sqr(scn));
2144
        if ((tmp > 0) && (sn2 > 0)) {
2145
            k    = (sn2*sck-scn*snk)/tmp;
2146
            kp   = (k*scn-snk)/sn2;
2147
            if (bError) {
2148
                chi2 = 0;
2149
                for(i=i0; (i<n); i++) {
2150
                    chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
2151
                }
2152
                chi22 = compute_weighted_rates(n,t,ct,nt,kt,sigma_ct,sigma_nt,
2153
                                               sigma_kt,&k,&kp,
2154
                                               &sigma_k,&sigma_kp,fit_start);
2155
                Q = quality_of_fit(chi2,2);
2156
                ddg = BOLTZ*temp*sigma_k/k;
2157
                printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n",
2158
                       chi2,Q);
2159
                printf("The Rate and Delta G are followed by an error estimate\n");
2160
                printf("----------------------------------------------------------\n"
2161
                       "Type      Rate (1/ps)  Sigma Time (ps)  DG (kJ/mol)  Sigma\n");
2162
                printf("Forward    %10.3f %6.2f   %8.3f  %10.3f %6.2f\n",
2163
                       k,sigma_k,1/k,calc_dg(1/k,temp),ddg);
2164
                ddg = BOLTZ*temp*sigma_kp/kp;
2165
                printf("Backward   %10.3f %6.2f   %8.3f  %10.3f %6.2f\n",
2166
                       kp,sigma_kp,1/kp,calc_dg(1/kp,temp),ddg);
2167
            }
2168
            else {
2169
                chi2 = 0;
2170
                for(i=i0; (i<n); i++) {
2171
                    chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
2172
                }
2173
                printf("Fitting parameters chi^2 = %10g\nQ = %10g\n",
2174
                       chi2,Q);
2175
                printf("--------------------------------------------------\n"
2176
                       "Type      Rate (1/ps) Time (ps)  DG (kJ/mol)  Chi^2\n");
2177
                printf("Forward    %10.3f   %8.3f  %10.3f  %10g\n",
2178
                       k,1/k,calc_dg(1/k,temp),chi2);
2179
                printf("Backward   %10.3f   %8.3f  %10.3f\n",
2180
                       kp,1/kp,calc_dg(1/kp,temp));
2181
            }
2182
        }
2183
        if (sc2 > 0) {
2184
            kow  = 2*sck/sc2;
2185
            printf("One-way    %10.3f   %s%8.3f  %10.3f\n",
2186
                   kow,bError ? "       " : "",1/kow,calc_dg(1/kow,temp));
2187
        }
2188
        else 
2189
            printf(" - Numerical problems computing HB thermodynamics:\n"
2190
                   "sc2 = %g  sn2 = %g  sk2 = %g sck = %g snk = %g scn = %g\n",
2191
                   sc2,sn2,sk2,sck,snk,scn);
2192
        /* Determine integral of the correlation function */
2193
        tau_hb = evaluate_integral(n,t,ct,NULL,(t[n-1]-t[0])/2,&dtau);
2194
        printf("Integral   %10.3f   %s%8.3f  %10.3f\n",1/tau_hb,
2195
               bError ? "       " : "",tau_hb,calc_dg(tau_hb,temp));
2196
        e_1 = exp(-1);
2197
        for(i=0; (i<n-2); i++) {
2198
            if ((ct[i] > e_1) && (ct[i+1] <= e_1)) {
2199
                break;
2200
            }
2201
        }
2202
        if (i < n-2) {
2203
            /* Determine tau_relax from linear interpolation */
2204
            tau_rlx = t[i]-t[0] + (e_1-ct[i])*(t[i+1]-t[i])/(ct[i+1]-ct[i]);
2205
            printf("Relaxation %10.3f   %8.3f  %s%10.3f\n",1/tau_rlx,
2206
                   tau_rlx,bError ? "       " : "",
2207
                   calc_dg(tau_rlx,temp));
2208
        }
2209
    }
2210
    else 
2211
        printf("Correlation functions too short to compute thermodynamics\n");
2212
}
2213

    
2214
void compute_derivative(int nn,real x[],real y[],real dydx[])
2215
{
2216
    int j;
2217
  
2218
    /* Compute k(t) = dc(t)/dt */
2219
    for(j=1; (j<nn-1); j++)
2220
        dydx[j] = (y[j+1]-y[j-1])/(x[j+1]-x[j-1]);
2221
    /* Extrapolate endpoints */
2222
    dydx[0]    = 2*dydx[1]   -  dydx[2];
2223
    dydx[nn-1] = 2*dydx[nn-2] - dydx[nn-3];
2224
}
2225

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

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

    
2236
static void normalizeACF(real *ct, real *gt, int len)
2237
{
2238
    real ct_fac, gt_fac;
2239
    int i;
2240

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

    
2243
    ct_fac = 1.0/ct[0];
2244
    gt_fac = (gt!=NULL && gt[0]!=0) ? 1.0/gt[0] : 0;
2245
    printf("Normalization for c(t) = %g for gh(t) = %g\n",ct_fac,gt_fac);
2246
    for (i=0; i<len; i++)
2247
    {
2248
        ct[i] *= ct_fac;
2249
        if (gt != NULL)
2250
            gt[i] *= gt_fac;
2251
    }
2252
}
2253

    
2254
/* Added argument bContact in do_hbac(...). Also
2255
 * added support for -contact in the actual code.
2256
 * - Erik Marklund, May 31, 2006 */
2257
/* Changed contact code and added argument R2 
2258
 * - Erik Marklund, June 29, 2006
2259
 */
2260
static void do_hbac(const char *fn,t_hbdata *hb,
2261
                    int nDump,gmx_bool bMerge,gmx_bool bContact, real fit_start,
2262
                    real temp,gmx_bool R2,real smooth_tail_start, const output_env_t oenv,
2263
                    t_gemParams *params, const char *gemType, int nThreads,
2264
                    const int NN, const gmx_bool bBallistic, const gmx_bool bGemFit)
2265
{
2266
    FILE *fp;
2267
    int  i,j,k,m,n,o,nd,ihb,idist,n2,nn,iter,nSets;
2268
    const char *legNN[]   = { "Ac(t)",
2269
                               "Ac'(t)"};
2270
    static char **legGem;
2271
                  
2272
    const char *legLuzar[] = { "Ac\\sfin sys\\v{}\\z{}(t)",
2273
                                "Ac(t)",
2274
                                "Cc\\scontact,hb\\v{}\\z{}(t)",
2275
                                "-dAc\\sfs\\v{}\\z{}/dt" };
2276
    gmx_bool bNorm=FALSE;
2277
    double nhb = 0;
2278
    int nhbi=0;
2279
    real *rhbex=NULL,*ht,*gt,*ght,*dght,*kt;
2280
    real *ct,*p_ct,tail,tail2,dtail,ct_fac,ght_fac,*cct;
2281
    const real tol = 1e-3;
2282
    int   nframes = hb->nframes,nf;
2283
    unsigned int **h,**g;
2284
    int   nh,nhbonds,nhydro,ngh;
2285
    t_hbond *hbh;
2286
    PSTYPE p, *pfound = NULL, np;
2287
    t_pShift *pHist;
2288
    int *ptimes=NULL, *poff=NULL, anhb, n0, mMax=INT_MIN;
2289
    real **rHbExGem = NULL;
2290
    gmx_bool c;
2291
    int acType;
2292
    t_E *E;
2293
    double *ctdouble, *timedouble, *fittedct;
2294
    double fittolerance=0.1;
2295

    
2296
    enum {AC_NONE, AC_NN, AC_GEM, AC_LUZAR};
2297

    
2298

    
2299
#ifdef HAVE_OPENMP
2300
    int *dondata=NULL, thisThread;
2301
#endif
2302

    
2303

    
2304
    printf("Doing autocorrelation ");
2305

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

    
2319
        nSets = 1 + (bBallistic ? 1:0) + (bGemFit ? 1:0);
2320
        snew(legGem, nSets);
2321
        for (i=0;i<nSets;i++)
2322
            snew(legGem[i], 128);
2323
        sprintf(legGem[0], "Ac\\s%s\\v{}\\z{}(t)", gemType);
2324
        if (bBallistic)
2325
            sprintf(legGem[1], "Ac'(t)");
2326
        if (bGemFit)
2327
            sprintf(legGem[(bBallistic ? 3:2)], "Ac\\s%s,fit\\v{}\\z{}(t)", gemType);
2328

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

    
2335
    /* build hbexist matrix in reals for autocorr */
2336
    /* Allocate memory for computing ACF (rhbex) and aggregating the ACF (ct) */
2337
    n2=1;
2338
    while (n2 < nframes)
2339
        n2*=2;
2340
  
2341
    nn = nframes/2;
2342
  
2343
    if (acType != AC_NN ||
2344
#ifndef HAVE_OPENMP
2345
        TRUE
2346
#else
2347
        FALSE
2348
#endif
2349
        ) {
2350
        snew(h,hb->maxhydro);
2351
        snew(g,hb->maxhydro);
2352
    }
2353

    
2354
    /* Dump hbonds for debugging */
2355
    dump_ac(hb,bMerge||bContact,nDump);
2356
  
2357
    /* Total number of hbonds analyzed here */
2358
    nhbonds = 0;
2359
    ngh     = 0;
2360
    anhb    = 0;
2361

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

    
2369
#ifdef HAVE_OPENMP  /* ================================================= \
2370
                     * Set up the OpenMP stuff,                           |
2371
                     * like the number of threads and such                |
2372
                     */
2373
    if (acType != AC_LUZAR)
2374
    {
2375
#if (_OPENMP >= 200805) /* =====================\ */
2376
        nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_thread_limit());
2377
#else
2378
        nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_num_procs());
2379
#endif /* _OPENMP >= 200805 ====================/ */
2380

    
2381
        omp_set_num_threads(nThreads);
2382
        snew(dondata, nThreads);
2383
        for (i=0; i<nThreads; i++)
2384
            dondata[i] = -1;
2385
        printf("ACF calculations parallelized with OpenMP using %i threads.\n"
2386
               "Expect close to linear scaling over this donor-loop.\n", nThreads);
2387
        fflush(stdout);
2388
        fprintf(stderr, "Donors: [thread no]\n");
2389
        {
2390
            char tmpstr[7];
2391
            for (i=0;i<nThreads;i++) {
2392
                snprintf(tmpstr, 7, "[%i]", i);
2393
                fprintf(stderr, "%-7s", tmpstr);
2394
            }
2395
        }
2396
        fprintf(stderr, "\n"); /*                                         | */
2397
    }  /*                                                                 | */
2398
#endif /* HAVE_OPENMP ===================================================/  */
2399

    
2400

    
2401
    /* Build the ACF according to acType */
2402
    switch (acType)
2403
    {
2404
      
2405
    case AC_NN:
2406
#ifdef HAVE_NN_LOOPS
2407
        /* Here we're using the estimated energy for the hydrogen bonds. */
2408
        snew(ct,nn);
2409
#ifdef HAVE_OPENMP /* ==================================\ */      
2410
#pragma omp parallel                            \
2411
    private(i, j, k, nh, E, rhbex, thisThread), \
2412
    default(shared)
2413
        {
2414
#pragma omp barrier
2415
            thisThread = omp_get_thread_num();
2416
            rhbex = NULL;
2417
#endif /* ==============================================/ */
2418

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

    
2422
#ifdef HAVE_OPENMP /* ################################################## \
2423
                    *                                                    #
2424
                    *                                                    #
2425
                    */
2426
#pragma omp barrier
2427
#pragma omp for schedule (dynamic)
2428
#endif
2429
            for (i=0; i<hb->d.nrd; i++) /* loop over donors */
2430
            {
2431
#ifdef HAVE_OPENMP /* ====== Write some output ======\ */
2432
#pragma omp critical
2433
                {
2434
                    dondata[thisThread] = i;
2435
                    parallel_print(dondata, nThreads);
2436
                }
2437
#else
2438
                fprintf(stderr, "\r %i", i);
2439
#endif /* ===========================================/ */
2440

    
2441
                for (j=0; j<hb->a.nra; j++) /* loop over acceptors */
2442
                {
2443
                    for (nh=0; nh<hb->d.nhydro[i]; nh++) /* loop over donors' hydrogens */
2444
                    {
2445
                        E = hb->hbE.E[i][j][nh];
2446
                        if (E != NULL)
2447
                        {
2448
                            for (k=0; k<nframes; k++)
2449
                            {
2450
                                if (E[k] != NONSENSE_E)
2451
                                    rhbex[k] = (real)E[k];
2452
                            }
2453
              
2454
                            low_do_autocorr(NULL,oenv,NULL,nframes,1,-1,&(rhbex),hb->time[1]-hb->time[0],
2455
                                            eacNormal,1,FALSE,bNorm,FALSE,0,-1,0,1);
2456
#ifdef HAVE_OPENMP
2457
#pragma omp critical
2458
#endif
2459
                            {
2460
                                for(k=0; (k<nn); k++)
2461
                                    ct[k] += rhbex[k];
2462
                            }
2463
                        }
2464
                    }   /* k loop */
2465
                }       /* j loop */
2466
            }           /* i loop */
2467
            sfree(rhbex);
2468
#pragma omp barrier
2469
#ifdef HAVE_OPENMP 
2470
            /*                                                           # */
2471
        } /* End of parallel block                                       # */
2472
        /* ##############################################################/ */
2473
        sfree(dondata);
2474
#endif
2475
        normalizeACF(ct, NULL, nn);
2476
        snew(ctdouble, nn);
2477
        snew(timedouble, nn);
2478
        for (j=0; j<nn; j++)
2479
        {
2480
            timedouble[j]=(double)(hb->time[j]);
2481
            ctdouble[j]=(double)(ct[j]);
2482
        }
2483

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

    
2501
        fp = xvgropen(fn, "Hydrogen Bond Autocorrelation",output_env_get_xvgr_tlabel(oenv),"C(t)");
2502
        xvgr_legend(fp,asize(legNN),legNN);
2503
      
2504
        for(j=0; (j<nn); j++)
2505
            fprintf(fp,"%10g  %10g %10g\n",
2506
                    hb->time[j]-hb->time[0],
2507
                    ct[j],
2508
                    ctdouble[j]);
2509
        fclose(fp);
2510
        sfree(ct);
2511
        sfree(ctdouble);
2512
        sfree(timedouble);
2513
#endif /* HAVE_NN_LOOPS */
2514
        break; /* case AC_NN */
2515

    
2516
    case AC_GEM:
2517
        snew(ct,2*n2);
2518
        memset(ct,0,2*n2*sizeof(real));
2519
#ifndef HAVE_OPENMP
2520
        fprintf(stderr, "Donor:\n");
2521
#define __ACDATA ct
2522
#else
2523
#define __ACDATA p_ct
2524
#endif
2525

    
2526
#ifdef HAVE_OPENMP /*  =========================================\
2527
                    *                                          */
2528
#pragma omp parallel default(none)                              \
2529
    private(i, k, nh, hbh, pHist, h, g, n0, nf, np, j, m,       \
2530
            pfound, poff, rHbExGem, p, ihb, mMax,               \
2531
            thisThread, p_ct)                                   \
2532
    shared(hb, dondata, ct, nn, nThreads, n2, stderr, bNorm,    \
2533
           nframes, bMerge, bContact)
2534
        { /* ##########  THE START OF THE ENORMOUS PARALLELIZED BLOCK!  ########## */
2535
            h = NULL;
2536
            g = NULL;
2537
            thisThread = omp_get_thread_num();
2538
            snew(h,hb->maxhydro);
2539
            snew(g,hb->maxhydro);
2540
            mMax = INT_MIN;
2541
            rHbExGem = NULL;
2542
            poff = NULL;
2543
            pfound = NULL;
2544
            p_ct = NULL;
2545
            snew(p_ct,2*n2);
2546
            memset(p_ct,0,2*n2*sizeof(real));
2547

    
2548
            /* I'm using a chunk size of 1, since I expect      \
2549
             * the overhead to be really small compared         \
2550
             * to the actual calculations                       \ */
2551
#pragma omp for schedule(dynamic,1) nowait /*                   \ */
2552
#endif /* HAVE_OPENMP  =========================================/ */
2553
      
2554
            for (i=0; i<hb->d.nrd; i++) {
2555
#ifdef HAVE_OPENMP
2556
#pragma omp critical
2557
                {
2558
                    dondata[thisThread] = i;
2559
                    parallel_print(dondata, nThreads);
2560
                }
2561
#else
2562
                fprintf(stderr, "\r %i", i);
2563
#endif
2564
    
2565
                for (k=0; k<hb->a.nra; k++) {
2566
                    for (nh=0; nh < ((bMerge || bContact) ? 1 : hb->d.nhydro[i]); nh++) {
2567
                        hbh = hb->hbmap[i][k];
2568
                        if (hbh) {
2569
                            /* Note that if hb->per->gemtype==gemDD, then distances will be stored in
2570
                             * hb->hbmap[d][a].h array anyway, because the contact flag will be set.
2571
                             * hence, it's only with the gemAD mode that hb->hbmap[d][a].g will be used. */
2572
                            pHist = &(hb->per->pHist[i][k]);
2573
                            if (ISHB(hbh->history[nh]) && pHist->len != 0) {
2574

    
2575
/* No need for a critical section */
2576
/* #ifdef HAVE_OPENMP */
2577
/* #pragma omp critical */
2578
/* #endif */
2579
                                {
2580
                                    h[nh] = hbh->h[nh];
2581
                                    g[nh] = hb->per->gemtype==gemAD ? hbh->g[nh] : NULL;
2582
                                }
2583
                                n0 = hbh->n0;
2584
                                nf = hbh->nframes;
2585
                                /* count the number of periodic shifts encountered and store
2586
                                 * them in separate arrays. */
2587
                                np = 0;
2588
                                for (j=0; j<pHist->len; j++)
2589
                                {
2590
                                    p = pHist->p[j];
2591
                                    for (m=0; m<=np; m++) {
2592
                                        if (m == np) { /* p not recognized in list. Add it and set up new array. */ 
2593
                                            np++;
2594
                                            if (np>hb->per->nper)
2595
                                                gmx_fatal(FARGS, "Too many pshifts. Something's utterly wrong here.");
2596
                                            if (m>=mMax) { /* Extend the arrays.
2597
                                                            * Doing it like this, using mMax to keep track of the sizes,
2598
                                                            * eleviates the need for freeing and re-allocating the arrays
2599
                                                            * when taking on the next donor-acceptor pair */
2600
                                                mMax = m;
2601
                                                srenew(pfound,np);   /* The list of found periodic shifts. */
2602
                                                srenew(rHbExGem,np); /* The hb existence functions (-aver_hb). */
2603
                                                snew(rHbExGem[m],2*n2);
2604
                                                srenew(poff,np);
2605
                                            }
2606

    
2607
/* This shouldn't have to be critical, right? */
2608
/* #ifdef HAVE_OPENMP */
2609
/* #pragma omp critical */
2610
/* #endif */
2611
                                            {
2612
                                                if (rHbExGem != NULL && rHbExGem[m] != NULL) {
2613
                                                    /* This must be done, as this array was most likey
2614
                                                     * used to store stuff in some previous iteration. */               
2615
                                                    memset(rHbExGem[m], 0, (sizeof(real)) * (2*n2));
2616
                                                }
2617
                                                else
2618
                                                    fprintf(stderr, "rHbExGem not initialized! m = %i\n", m);
2619
                                            }
2620
                                            pfound[m] = p;
2621
                                            poff[m] = -1;
2622
              
2623
                                            break;
2624
                                        } /* m==np */
2625
                                        if (p == pfound[m])
2626
                                            break;
2627
                                    } /* m: Loop over found shifts */
2628
                                }   /* j: Loop over shifts */
2629

    
2630
                                /* Now unpack and disentangle the existence funtions. */
2631
                                for (j=0; j<nf; j++) {
2632
                                    /* i:       donor,
2633
                                     * k:       acceptor
2634
                                     * nh:      hydrogen
2635
                                     * j:       time
2636
                                     * p:       periodic shift
2637
                                     * pfound:  list of periodic shifts found for this pair.
2638
                                     * poff:    list of frame offsets; that is, the first
2639
                                     *          frame a hbond has a particular periodic shift. */
2640
                                    p = getPshift(*pHist, j+n0);
2641
                                    if (p != -1)
2642
                                    {
2643
                                        for (m=0; m<np; m++)
2644
                                        {
2645
                                            if (pfound[m] == p)
2646
                                                break;
2647
                                            if (m==(np-1))
2648
                                                gmx_fatal(FARGS,"Shift not found, but must be there.");
2649
                                        }
2650

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

    
2664
                                /* Now, build ac. */
2665
                                for (m=0; m<np; m++) {
2666
                                    if (rHbExGem[m][0]>0  && n0+poff[m]<nn/*  && m==0 */) {
2667
                                        low_do_autocorr(NULL,oenv,NULL,nframes,1,-1,&(rHbExGem[m]),hb->time[1]-hb->time[0],
2668
                                                        eacNormal,1,FALSE,bNorm,FALSE,0,-1,0,1);
2669
                                        for(j=0; (j<nn); j++)
2670
                                            __ACDATA[j] += rHbExGem[m][j];
2671
                                    }
2672
                                } /* Building of ac. */
2673
                            } /* if (ISHB(...*/
2674
                        } /* if (hbh) */
2675
                    } /* hydrogen loop */
2676
                } /* acceptor loop */
2677
            } /* donor loop */
2678

    
2679
            for (m=0; m<=mMax; m++) {
2680
                sfree(rHbExGem[m]);
2681
            }
2682
            sfree(pfound);
2683
            sfree(poff);
2684
            sfree(rHbExGem);
2685

    
2686
            sfree(h);
2687
            sfree(g);
2688
#ifdef HAVE_OPENMP /* =======================================\ */
2689
#pragma omp critical
2690
            {
2691
                for (i=0; i<nn; i++)
2692
                    ct[i] += p_ct[i];
2693
            }
2694
            sfree(p_ct);
2695

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

    
2700
        normalizeACF(ct, NULL, nn);
2701

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

    
2704
        /* Use this part to fit to geminate recombination - JCP 129, 84505 (2008) */
2705
      
2706
        snew(ctdouble, nn);
2707
        snew(timedouble, nn);
2708
        snew(fittedct, nn);
2709
    
2710
        for (j=0;j<nn;j++){
2711
            timedouble[j]=(double)(hb->time[j]);
2712
            ctdouble[j]=(double)(ct[j]);
2713
        }
2714

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

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

    
2743
        for(j=0; (j<nn); j++)
2744
        {
2745
            fprintf(fp, "%10g  %10g", hb->time[j]-hb->time[0],ct[j]);
2746
            if (bBallistic)
2747
                fprintf(fp,"  %10g", ctdouble[j]);
2748
            if (bGemFit)
2749
                fprintf(fp,"  %10g", fittedct[j]);
2750
            fprintf(fp,"\n");
2751
        }
2752
        fclose(fp);
2753

    
2754
        sfree(ctdouble);
2755
        sfree(timedouble);
2756
        sfree(fittedct);
2757
        sfree(ct);
2758

    
2759
        break; /* case AC_GEM */
2760

    
2761
    case AC_LUZAR:
2762
        snew(rhbex,2*n2);
2763
        snew(ct,2*n2);
2764
        snew(gt,2*n2);
2765
        snew(ht,2*n2);
2766
        snew(ght,2*n2);
2767
        snew(dght,2*n2);
2768

    
2769
        snew(kt,nn);
2770
        snew(cct,nn);
2771
    
2772
        for(i=0; (i<hb->d.nrd); i++) {
2773
            for(k=0; (k<hb->a.nra); k++) {
2774
                nhydro = 0;
2775
                hbh = hb->hbmap[i][k];
2776
   
2777
                if (hbh) {
2778
                    if (bMerge || bContact) {
2779
                        if (ISHB(hbh->history[0])) {
2780
                            h[0] = hbh->h[0];
2781
                            g[0] = hbh->g[0];
2782
                            nhydro = 1;
2783
                        }
2784
                    }
2785
                    else {
2786
                        for(m=0; (m<hb->maxhydro); m++)
2787
                            if (bContact ? ISDIST(hbh->history[m]) : ISHB(hbh->history[m])) {
2788
                                g[nhydro] = hbh->g[m];
2789
                                h[nhydro] = hbh->h[m];
2790
                                nhydro++;
2791
                            }
2792
                    }
2793
    
2794
                    nf = hbh->nframes;
2795
                    for(nh=0; (nh<nhydro); nh++) {
2796
                        int nrint = bContact ? hb->nrdist : hb->nrhb;
2797
                        if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint))
2798
                            fprintf(stderr,"\rACF %d/%d",nhbonds+1,nrint);
2799
                        nhbonds++;
2800
                        for(j=0; (j<nframes); j++) {
2801
                            /* Changed '<' into '<=' below, just like I did in
2802
                               the hbm-output-loop in the gmx_hbond() block.
2803
                               - Erik Marklund, May 31, 2006 */
2804
                            if (j <= nf) {
2805
                                ihb   = is_hb(h[nh],j);
2806
                                idist = is_hb(g[nh],j);
2807
                            }
2808
                            else {
2809
                                ihb = idist = 0;
2810
                            }
2811
                            rhbex[j] = ihb;
2812
                            /* For contacts: if a second cut-off is provided, use it,
2813
                             * otherwise use g(t) = 1-h(t) */
2814
                            if (!R2 && bContact)
2815
                                gt[j]  = 1-ihb;
2816
                            else
2817
                                gt[j]  = idist*(1-ihb); 
2818
                            ht[j]    = rhbex[j];
2819
                            nhb     += ihb;
2820
                        }
2821
      
2822

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

    
2833
                        cross_corr(n2,ht,gt,dght);
2834
      
2835
                        for(j=0; (j<nn); j++) {
2836
                            ct[j]  += rhbex[j];
2837
                            ght[j] += dght[j];
2838
                        }
2839
                    }
2840
                }
2841
            }
2842
        }
2843
        fprintf(stderr,"\n");
2844
        sfree(h);
2845
        sfree(g);
2846
        normalizeACF(ct, gt, nn);
2847

    
2848
        /* Determine tail value for statistics */
2849
        tail  = 0;
2850
        tail2 = 0;
2851
        for(j=nn/2; (j<nn); j++) {
2852
            tail  += ct[j];
2853
            tail2 += ct[j]*ct[j];
2854
        }
2855
        tail  /= (nn - nn/2);
2856
        tail2 /= (nn - nn/2);
2857
        dtail  = sqrt(tail2-tail*tail);
2858
  
2859
        /* Check whether the ACF is long enough */
2860
        if (dtail > tol) {
2861
            printf("\nWARNING: Correlation function is probably not long enough\n"
2862
                   "because the standard deviation in the tail of C(t) > %g\n"
2863
                   "Tail value (average C(t) over second half of acf): %g +/- %g\n",
2864
                   tol,tail,dtail);
2865
        }
2866
        for(j=0; (j<nn); j++) {
2867
            cct[j] = ct[j];
2868
            ct[j]  = (cct[j]-tail)/(1-tail); 
2869
        }
2870
        /* Compute negative derivative k(t) = -dc(t)/dt */
2871
        compute_derivative(nn,hb->time,ct,kt);
2872
        for(j=0; (j<nn); j++)
2873
            kt[j] = -kt[j];
2874

    
2875

    
2876
        if (bContact)
2877
            fp = xvgropen(fn, "Contact Autocorrelation",output_env_get_xvgr_tlabel(oenv),"C(t)",oenv);
2878
        else
2879
            fp = xvgropen(fn, "Hydrogen Bond Autocorrelation",output_env_get_xvgr_tlabel(oenv),"C(t)",oenv);
2880
        xvgr_legend(fp,asize(legLuzar),legLuzar, oenv);
2881

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

    
2888
        analyse_corr(nn,hb->time,ct,ght,kt,NULL,NULL,NULL,
2889
                     fit_start,temp,smooth_tail_start,oenv);
2890
  
2891
        do_view(oenv,fn,NULL);
2892
        sfree(rhbex);
2893
        sfree(ct);
2894
        sfree(gt);
2895
        sfree(ht);
2896
        sfree(ght);
2897
        sfree(dght);
2898
        sfree(cct);
2899
        sfree(kt);
2900
        /* sfree(h); */
2901
/*         sfree(g); */
2902

    
2903
        break; /* case AC_LUZAR */
2904

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

    
2910
static void init_hbframe(t_hbdata *hb,int nframes,real t)
2911
{
2912
    int i,j,m;
2913
  
2914
    hb->time[nframes]   = t;
2915
    hb->nhb[nframes]    = 0;
2916
    hb->ndist[nframes]  = 0;
2917
    for (i=0; (i<max_hx); i++)
2918
        hb->nhx[nframes][i]=0;
2919
    /* Loop invalidated */
2920
    if (hb->bHBmap && 0)
2921
        for (i=0; (i<hb->d.nrd); i++)
2922
            for (j=0; (j<hb->a.nra); j++)
2923
                for (m=0; (m<hb->maxhydro); m++)
2924
                    if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[m])
2925
                        set_hb(hb,i,m,j,nframes,HB_NO);
2926
    /*set_hb(hb->hbmap[i][j]->h[m],nframes-hb->hbmap[i][j]->n0,HB_NO);*/
2927
}
2928

    
2929
static void analyse_donor_props(const char *fn,t_hbdata *hb,int nframes,real t,
2930
                                const output_env_t oenv)
2931
{
2932
    static FILE *fp = NULL;
2933
    const char *leg[] = { "Nbound", "Nfree" };
2934
    int i,j,k,nbound,nb,nhtot;
2935
  
2936
    if (!fn)
2937
        return;
2938
    if (!fp) {
2939
        fp = xvgropen(fn,"Donor properties",output_env_get_xvgr_tlabel(oenv),"Number",oenv);
2940
        xvgr_legend(fp,asize(leg),leg,oenv);
2941
    }
2942
    nbound = 0;
2943
    nhtot  = 0;
2944
    for(i=0; (i<hb->d.nrd); i++) {
2945
        for(k=0; (k<hb->d.nhydro[i]); k++) {
2946
            nb = 0;
2947
            nhtot ++;
2948
            for(j=0; (j<hb->a.nra) && (nb == 0); j++) {
2949
                if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[k] && 
2950
                    is_hb(hb->hbmap[i][j]->h[k],nframes)) 
2951
                    nb = 1;
2952
            }
2953
            nbound += nb;
2954
        }
2955
    }
2956
    fprintf(fp,"%10.3e  %6d  %6d\n",t,nbound,nhtot-nbound);
2957
}
2958

    
2959
static void dump_hbmap(t_hbdata *hb,
2960
                       int nfile,t_filenm fnm[],gmx_bool bTwo,
2961
                       gmx_bool bContact, int isize[],int *index[],char *grpnames[],
2962
                       t_atoms *atoms)
2963
{
2964
    FILE *fp,*fplog;
2965
    int  ddd,hhh,aaa,i,j,k,m,grp;
2966
    char ds[32],hs[32],as[32];
2967
    gmx_bool first;
2968
  
2969
    fp = opt2FILE("-hbn",nfile,fnm,"w");
2970
    if (opt2bSet("-g",nfile,fnm)) {
2971
        fplog = ffopen(opt2fn("-g",nfile,fnm),"w");
2972
        if (bContact)
2973
            fprintf(fplog,"# %10s  %12s  %12s\n","Donor","Hydrogen","Acceptor");
2974
        else
2975
            fprintf(fplog,"# %10s  %12s  %12s\n","Donor","Hydrogen","Acceptor");
2976
    }
2977
    else
2978
        fplog = NULL;
2979
    for (grp=gr0; grp<=(bTwo?gr1:gr0); grp++) {
2980
        fprintf(fp,"[ %s ]",grpnames[grp]);
2981
        for (i=0; i<isize[grp]; i++) {
2982
            fprintf(fp,(i%15)?" ":"\n");
2983
            fprintf(fp," %4u",index[grp][i]+1);
2984
        }
2985
        fprintf(fp,"\n");
2986
        /*
2987
          Added -contact support below.
2988
          - Erik Marklund, May 29, 2006
2989
        */
2990
        if (!bContact) {
2991
            fprintf(fp,"[ donors_hydrogens_%s ]\n",grpnames[grp]);
2992
            for (i=0; (i<hb->d.nrd); i++) {
2993
                if (hb->d.grp[i] == grp) { 
2994
                    for(j=0; (j<hb->d.nhydro[i]); j++)
2995
                        fprintf(fp," %4u %4u",hb->d.don[i]+1,
2996
                                hb->d.hydro[i][j]+1);
2997
                    fprintf(fp,"\n");
2998
                }
2999
            }
3000
            first = TRUE;
3001
            fprintf(fp,"[ acceptors_%s ]",grpnames[grp]);
3002
            for (i=0; (i<hb->a.nra); i++) {
3003
                if (hb->a.grp[i] == grp) { 
3004
                    fprintf(fp,(i%15 && !first)?" ":"\n");
3005
                    fprintf(fp," %4u",hb->a.acc[i]+1);
3006
                    first = FALSE;
3007
                }
3008
            }
3009
            fprintf(fp,"\n");
3010
        }
3011
    }
3012
    if (bTwo)
3013
        fprintf(fp,bContact ? "[ contacts_%s-%s ]\n" :
3014
                "[ hbonds_%s-%s ]\n",grpnames[0],grpnames[1]);
3015
    else
3016
        fprintf(fp,bContact ? "[ contacts_%s ]" : "[ hbonds_%s ]\n",grpnames[0]);
3017
  
3018
    for(i=0; (i<hb->d.nrd); i++) {
3019
        ddd = hb->d.don[i];
3020
        for(k=0; (k<hb->a.nra); k++) {
3021
            aaa = hb->a.acc[k];
3022
            for(m=0; (m<hb->d.nhydro[i]); m++) {
3023
                if (hb->hbmap[i][k] && ISHB(hb->hbmap[i][k]->history[m])) {
3024
                    sprintf(ds,"%s",mkatomname(atoms,ddd));
3025
                    sprintf(as,"%s",mkatomname(atoms,aaa));
3026
                    if (bContact) {
3027
                        fprintf(fp," %6u %6u\n",ddd+1,aaa+1);
3028
                        if (fplog) 
3029
                            fprintf(fplog,"%12s  %12s\n",ds,as);
3030
                    } else {
3031
                        hhh = hb->d.hydro[i][m];
3032
                        sprintf(hs,"%s",mkatomname(atoms,hhh));
3033
                        fprintf(fp," %6u %6u %6u\n",ddd+1,hhh+1,aaa+1);
3034
                        if (fplog) 
3035
                            fprintf(fplog,"%12s  %12s  %12s\n",ds,hs,as);
3036
                    }
3037
                }
3038
            }
3039
        }
3040
    }
3041
    ffclose(fp);
3042
    if (fplog)
3043
        ffclose(fplog);
3044
}
3045

    
3046
#ifdef HAVE_OPENMP
3047
/* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template.
3048
 * It mimics add_frames() and init_frame() to some extent. */
3049
static void sync_hbdata(t_hbdata *hb, t_hbdata *p_hb,
3050
                        int nframes, real t)
3051
{
3052
    int i;
3053
    if (nframes >= p_hb->max_frames)
3054
    {
3055
        p_hb->max_frames += 4096;
3056
        srenew(p_hb->nhb,   p_hb->max_frames);
3057
        srenew(p_hb->ndist, p_hb->max_frames);
3058
        srenew(p_hb->n_bound, p_hb->max_frames);
3059
        srenew(p_hb->nhx, p_hb->max_frames);
3060
        if (p_hb->bDAnr)
3061
            srenew(p_hb->danr, p_hb->max_frames);
3062
        memset(&(p_hb->nhb[nframes]),   0, sizeof(int) * (p_hb->max_frames-nframes));
3063
        memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
3064
        p_hb->nhb[nframes] = 0;
3065
        p_hb->ndist[nframes] = 0;
3066

    
3067
    }
3068
    p_hb->nframes = nframes;
3069
/*     for (i=0;) */
3070
/*     { */
3071
/*         p_hb->nhx[nframes][i] */
3072
/*     } */
3073
    memset(&(p_hb->nhx[nframes]), 0 ,sizeof(int)*max_hx); /* zero the helix count for this frame */
3074

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

    
3081
int gmx_hbond(int argc,char *argv[])
3082
{
3083
    const char *desc[] = {
3084
        "g_hbond computes and analyzes hydrogen bonds. Hydrogen bonds are",
3085
        "determined based on cutoffs for the angle Acceptor - Donor - Hydrogen",
3086
        "(zero is extended) and the distance Hydrogen - Acceptor.",
3087
        "OH and NH groups are regarded as donors, O is an acceptor always,",
3088
        "N is an acceptor by default, but this can be switched using",
3089
        "[TT]-nitacc[tt]. Dummy hydrogen atoms are assumed to be connected",
3090
        "to the first preceding non-hydrogen atom.[PAR]",
3091
    
3092
        "You need to specify two groups for analysis, which must be either",
3093
        "identical or non-overlapping. All hydrogen bonds between the two",
3094
        "groups are analyzed.[PAR]",
3095
    
3096
        "If you set -shell, you will be asked for an additional index group",
3097
        "which should contain exactly one atom. In this case, only hydrogen",
3098
        "bonds between atoms within the shell distance from the one atom are",
3099
        "considered.[PAR]",
3100
    
3101
        /*    "It is also possible to analyse specific hydrogen bonds with",
3102
              "[TT]-sel[tt]. This index file must contain a group of atom triplets",
3103
              "Donor Hydrogen Acceptor, in the following way:[PAR]",
3104
        */
3105
        "[TT]",
3106
        "[ selected ][BR]",
3107
        "     20    21    24[BR]",
3108
        "     25    26    29[BR]",
3109
        "      1     3     6[BR]",
3110
        "[tt][BR]",
3111
        "Note that the triplets need not be on separate lines.",
3112
        "Each atom triplet specifies a hydrogen bond to be analyzed,",
3113
        "note also that no check is made for the types of atoms.[PAR]",
3114
     
3115
        "[BB]Output:[bb][BR]",
3116
        "[TT]-num[tt]:  number of hydrogen bonds as a function of time.[BR]",
3117
        "[TT]-ac[tt]:   average over all autocorrelations of the existence",
3118
        "functions (either 0 or 1) of all hydrogen bonds.[BR]",
3119
        "[TT]-dist[tt]: distance distribution of all hydrogen bonds.[BR]",
3120
        "[TT]-ang[tt]:  angle distribution of all hydrogen bonds.[BR]",
3121
        "[TT]-hx[tt]:   the number of n-n+i hydrogen bonds as a function of time",
3122
        "where n and n+i stand for residue numbers and i ranges from 0 to 6.",
3123
        "This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated",
3124
        "with helices in proteins.[BR]",
3125
        "[TT]-hbn[tt]:  all selected groups, donors, hydrogens and acceptors",
3126
        "for selected groups, all hydrogen bonded atoms from all groups and",
3127
        "all solvent atoms involved in insertion.[BR]",
3128
        "[TT]-hbm[tt]:  existence matrix for all hydrogen bonds over all",
3129
        "frames, this also contains information on solvent insertion",
3130
        "into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]",
3131
        "index file.[BR]",
3132
        "[TT]-dan[tt]: write out the number of donors and acceptors analyzed for",
3133
        "each timeframe. This is especially useful when using [TT]-shell[tt].[BR]",
3134
        "[TT]-nhbdist[tt]: compute the number of HBonds per hydrogen in order to",
3135
        "compare results to Raman Spectroscopy.",
3136
        "[PAR]",
3137
        "Note: options [TT]-ac[tt], [TT]-life[tt], [TT]-hbn[tt] and [TT]-hbm[tt]",
3138
        "require an amount of memory proportional to the total numbers of donors",
3139
        "times the total number of acceptors in the selected group(s)."
3140
    };
3141
  
3142
    static real acut=30, abin=1, rcut=0.35, r2cut=0, rbin=0.005, rshell=-1;
3143
    static real maxnhb=0,fit_start=1,fit_end=60,temp=298.15,smooth_tail_start=-1, D=-1;
3144
    static gmx_bool bNitAcc=TRUE,bDA=TRUE,bMerge=TRUE;
3145
    static int  nDump=0, nFitPoints=100;
3146
    static int nThreads = 0, nBalExp=4;
3147

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

    
3152
    /* options */
3153
    t_pargs pa [] = {
3154
        { "-a",    FALSE,  etREAL, {&acut},
3155
          "Cutoff angle (degrees, Acceptor - Donor - Hydrogen)" },
3156
        { "-r",    FALSE,  etREAL, {&rcut},
3157
          "Cutoff radius (nm, X - Acceptor, see next option)" },
3158
        { "-da",   FALSE,  etBOOL, {&bDA},
3159
          "Use distance Donor-Acceptor (if TRUE) or Hydrogen-Acceptor (FALSE)" },
3160
        { "-r2",   FALSE,  etREAL, {&r2cut},
3161
          "Second cutoff radius. Mainly useful with -contact and -ac"},
3162
        { "-abin", FALSE,  etREAL, {&abin},
3163
          "Binwidth angle distribution (degrees)" },
3164
        { "-rbin", FALSE,  etREAL, {&rbin},
3165
          "Binwidth distance distribution (nm)" },
3166
        { "-nitacc",FALSE, etBOOL, {&bNitAcc},
3167
          "Regard nitrogen atoms as acceptors" },
3168
        { "-contact",FALSE,etBOOL, {&bContact},
3169
          "Do not look for hydrogen bonds, but merely for contacts within the cut-off distance" },
3170
        { "-shell", FALSE, etREAL, {&rshell},
3171
          "when > 0, only calculate hydrogen bonds within # nm shell around "
3172
          "one particle" },
3173
        { "-fitstart", FALSE, etREAL, {&fit_start},
3174
          "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 -gemfit we suggest -fitstart 0" },
3175
        { "-fitstart", FALSE, etREAL, {&fit_start},
3176
          "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 -gemfit)" },
3177
        { "-temp",  FALSE, etREAL, {&temp},
3178
          "Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" },
3179
        { "-smooth",FALSE, etREAL, {&smooth_tail_start},
3180
          "If >= 0, the tail of the ACF will be smoothed by fitting it to an exponential function: y = A exp(-x/tau)" },
3181
        { "-dump",  FALSE, etINT, {&nDump},
3182
          "Dump the first N hydrogen bond ACFs in a single xvg file for debugging" },
3183
        { "-max_hb",FALSE, etREAL, {&maxnhb},
3184
          "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" },
3185
        { "-merge", FALSE, etBOOL, {&bMerge},
3186
          "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." },
3187
        { "-geminate", FALSE, etENUM, {gemType},
3188
          "Use reversible geminate recombination for the kinetics/thermodynamics calclations. See Markovitch et al., J. Chem. Phys 129, 084505 (2008) for details."},
3189
        { "-diff", FALSE, etREAL, {&D},
3190
          "Dffusion coefficient to use in the rev. gem. recomb. kinetic model. If non-positive, then it will be fitted to the ACF along with ka and kd."},
3191
#ifdef HAVE_OPENMP
3192
        { "-nthreads", FALSE, etINT, {&nThreads},
3193
          "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)"},
3194
#endif
3195
    };
3196
    const char *bugs[] = {
3197
        "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and therefore not available for the time being."
3198
    };
3199
    t_filenm fnm[] = {
3200
        { efTRX, "-f",   NULL,     ffREAD  },
3201
        { efTPX, NULL,   NULL,     ffREAD  },
3202
        { efNDX, NULL,   NULL,     ffOPTRD },
3203
        /*    { efNDX, "-sel", "select", ffOPTRD },*/
3204
        { efXVG, "-num", "hbnum",  ffWRITE },
3205
        { efLOG, "-g",   "hbond",  ffOPTWR },
3206
        { efXVG, "-ac",  "hbac",   ffOPTWR },
3207
        { efXVG, "-dist","hbdist", ffOPTWR },
3208
        { efXVG, "-ang", "hbang",  ffOPTWR },
3209
        { efXVG, "-hx",  "hbhelix",ffOPTWR },
3210
        { efNDX, "-hbn", "hbond",  ffOPTWR },
3211
        { efXPM, "-hbm", "hbmap",  ffOPTWR },
3212
        { efXVG, "-don", "donor",  ffOPTWR },
3213
        { efXVG, "-dan", "danum",  ffOPTWR },
3214
        { efXVG, "-life","hblife", ffOPTWR },
3215
        { efXVG, "-nhbdist", "nhbdist",ffOPTWR }
3216
    
3217
    };
3218
#define NFILE asize(fnm)
3219
  
3220
    char  hbmap [HB_NR]={ ' ',    'o',      '-',       '*' };
3221
    const char *hbdesc[HB_NR]={ "None", "Present", "Inserted", "Present & Inserted" };
3222
    t_rgb hbrgb [HB_NR]={ {1,1,1},{1,0,0},   {0,0,1},    {1,0,1} };
3223

    
3224
    t_trxstatus *status;
3225
    int trrStatus=1;
3226
    t_topology top;
3227
    t_inputrec ir;
3228
    t_pargs *ppa;
3229
    int     npargs,natoms,nframes=0,shatom;
3230
    int     *isize;
3231
    char    **grpnames;
3232
    atom_id **index;
3233
    rvec    *x,hbox;
3234
    matrix  box;
3235
    real    t,ccut,dist=0.0,ang=0.0;
3236
    double  max_nhb,aver_nhb,aver_dist;
3237
    int     h=0,i,j,k=0,l,start,end,id,ja,ogrp,nsel;
3238
    int     xi,yi,zi,ai;
3239
    int     xj,yj,zj,aj,xjj,yjj,zjj;
3240
    int     xk,yk,zk,ak,xkk,ykk,zkk;
3241
    gmx_bool    bSelected,bHBmap,bStop,bTwo,was,bBox,bTric;
3242
    int     *adist,*rdist;
3243
    int        grp,nabin,nrbin,bin,resdist,ihb;
3244
    char       **leg;
3245
    t_hbdata   *hb;
3246
    FILE       *fp,*fpins=NULL,*fpnhb=NULL;
3247
    t_gridcell ***grid;
3248
    t_ncell    *icell,*jcell,*kcell;
3249
    ivec       ngrid;
3250
    unsigned char        *datable;
3251
    output_env_t oenv;
3252
    int     gemmode, NN;
3253
    PSTYPE  peri=0;
3254
    t_E     E;
3255
    int     ii, jj, hh, actual_nThreads;
3256
    int     threadNr=0;
3257
    gmx_bool    bGem, bNN, bParallel;
3258
    t_gemParams *params=NULL;
3259
    
3260
    CopyRight(stdout,argv[0]);
3261

    
3262
    npargs = asize(pa);  
3263
    ppa    = add_acf_pargs(&npargs,pa);
3264
  
3265
    parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,NFILE,fnm,npargs,
3266
                      ppa,asize(desc),desc,asize(bugs),bugs,&oenv);
3267

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

    
3280
    if (NN > NN_NONE && bMerge)
3281
        bMerge = FALSE;
3282

    
3283
    /* geminate recombination? If so, which flavor? */
3284
    gemmode = 1;
3285
    while (gemmode < gemNR && gmx_strcasecmp(gemType[0], gemType[gemmode])!=0)
3286
        gemmode++;
3287
    if (gemmode == gemNR)
3288
        gmx_fatal(FARGS, "Invalid recombination type.");
3289
  
3290
    bGem = FALSE;
3291
    for (i=2; bGem==FALSE && i<gemNR; i++)
3292
        bGem = bGem || gemmode == i;
3293
  
3294
    if (bGem) {
3295
        printf("Geminate recombination: %s\n" ,gemType[gemmode]);
3296
#ifndef HAVE_LIBGSL
3297
        printf("Note that some aspects of reversible geminate recombination won't work without gsl.\n");
3298
#endif
3299
        if (bContact) {
3300
            if (gemmode != gemDD) {
3301
                printf("Turning off -contact option...\n");
3302
                bContact = FALSE;
3303
            }
3304
        } else {
3305
            if (gemmode == gemDD) {
3306
                printf("Turning on -contact option...\n");
3307
                bContact = TRUE;
3308
            }
3309
        }
3310
        if (bMerge) {
3311
            if (gemmode == gemAA) {
3312
                printf("Turning off -merge option...\n");
3313
                bMerge = FALSE;
3314
            }
3315
        } else {
3316
            if (gemmode != gemAA) {
3317
                printf("Turning on -merge option...\n");
3318
                bMerge = TRUE;
3319
            }
3320
        }
3321
    } 
3322
    
3323
    /* process input */
3324
    bSelected = FALSE;
3325
    ccut = cos(acut*DEG2RAD);
3326
  
3327
    if (bContact) {
3328
        if (bSelected)
3329
            gmx_fatal(FARGS,"Can not analyze selected contacts.");
3330
        if (!bDA) {
3331
            gmx_fatal(FARGS,"Can not analyze contact between H and A: turn off -noda");
3332
        }
3333
    }
3334

    
3335
#ifndef HAVE_LIBGSL
3336
    /* Don't pollute stdout with information about external libraries.
3337
     *
3338
     * printf("NO GSL! Can't find and take away ballistic term in ACF without GSL\n.");
3339
     */
3340
#endif
3341
  
3342
    /* Initiate main data structure! */
3343
    bHBmap = (opt2bSet("-ac",NFILE,fnm) ||
3344
              opt2bSet("-life",NFILE,fnm) ||
3345
              opt2bSet("-hbn",NFILE,fnm) || 
3346
              opt2bSet("-hbm",NFILE,fnm) ||
3347
              bGem);
3348
  
3349
#ifdef HAVE_OPENMP
3350
    /* Same thing here. There is no reason whatsoever to write the specific version of
3351
     * OpenMP used for compilation to stdout for normal usage.
3352
     *
3353
     * printf("Compiled with OpenMP (%i)\n", _OPENMP);
3354
     */
3355
#endif
3356

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

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

    
3369
    /* get topology */
3370
    read_tpx_top(ftp2fn(efTPX,NFILE,fnm),&ir,box,&natoms,NULL,NULL,NULL,&top);
3371
  
3372
    snew(grpnames,grNR);
3373
    snew(index,grNR);
3374
    snew(isize,grNR);
3375
    /* Make Donor-Acceptor table */
3376
    snew(datable, top.atoms.nr);
3377
    gen_datable(index[0],isize[0],datable,top.atoms.nr);
3378
  
3379
    if (bSelected) {
3380
        /* analyze selected hydrogen bonds */
3381
        printf("Select group with selected atoms:\n");
3382
        get_index(&(top.atoms),opt2fn("-sel",NFILE,fnm),
3383
                  1,&nsel,index,grpnames);
3384
        if (nsel % 3)
3385
            gmx_fatal(FARGS,"Number of atoms in group '%s' not a multiple of 3\n"
3386
                      "and therefore cannot contain triplets of "
3387
                      "Donor-Hydrogen-Acceptor",grpnames[0]);
3388
        bTwo=FALSE;
3389
    
3390
        for(i=0; (i<nsel); i+=3) {
3391
            int dd = index[0][i];
3392
            int aa = index[0][i+2];
3393
            /* int */ hh = index[0][i+1];
3394
            add_dh (&hb->d,dd,hh,i,datable);
3395
            add_acc(&hb->a,aa,i);
3396
            /* Should this be here ? */
3397
            snew(hb->d.dptr,top.atoms.nr);
3398
            snew(hb->a.aptr,top.atoms.nr);
3399
            add_hbond(hb,dd,aa,hh,gr0,gr0,0,bMerge,0,bContact,peri);
3400
        }
3401
        printf("Analyzing %d selected hydrogen bonds from '%s'\n",
3402
               isize[0],grpnames[0]);
3403
    } 
3404
    else {
3405
        /* analyze all hydrogen bonds: get group(s) */
3406
        printf("Specify 2 groups to analyze:\n");
3407
        get_index(&(top.atoms),ftp2fn_null(efNDX,NFILE,fnm),
3408
                  2,isize,index,grpnames);
3409
    
3410
        /* check if we have two identical or two non-overlapping groups */
3411
        bTwo = isize[0] != isize[1];
3412
        for (i=0; (i<isize[0]) && !bTwo; i++)
3413
            bTwo = index[0][i] != index[1][i];
3414
        if (bTwo) {
3415
            printf("Checking for overlap in atoms between %s and %s\n",
3416
                   grpnames[0],grpnames[1]);
3417
            for (i=0; i<isize[1];i++)
3418
                if (ISINGRP(datable[index[1][i]]))
3419
                    gmx_fatal(FARGS,"Partial overlap between groups '%s' and '%s'",
3420
                              grpnames[0],grpnames[1]);
3421
            /*
3422
              printf("Checking for overlap in atoms between %s and %s\n",
3423
              grpnames[0],grpnames[1]);
3424
              for (i=0; i<isize[0]; i++)
3425
              for (j=0; j<isize[1]; j++)
3426
              if (index[0][i] == index[1][j]) 
3427
              gmx_fatal(FARGS,"Partial overlap between groups '%s' and '%s'",
3428
              grpnames[0],grpnames[1]);
3429
            */
3430
        }
3431
        if (bTwo)
3432
            printf("Calculating %s "
3433
                   "between %s (%d atoms) and %s (%d atoms)\n",
3434
                   bContact ? "contacts" : "hydrogen bonds",
3435
                   grpnames[0],isize[0],grpnames[1],isize[1]);
3436
        else
3437
            fprintf(stderr,"Calculating %s in %s (%d atoms)\n",
3438
                    bContact?"contacts":"hydrogen bonds",grpnames[0],isize[0]);
3439
    }
3440
    sfree(datable);
3441
  
3442
    /* search donors and acceptors in groups */
3443
    snew(datable, top.atoms.nr);
3444
    for (i=0; (i<grNR); i++)
3445
        if ( ((i==gr0) && !bSelected ) ||
3446
             ((i==gr1) && bTwo )) {
3447
            gen_datable(index[i],isize[i],datable,top.atoms.nr);
3448
            if (bContact) {
3449
                search_acceptors(&top,isize[i],index[i],&hb->a,i,
3450
                                 bNitAcc,TRUE,(bTwo && (i==gr0)) || !bTwo, datable);
3451
                search_donors   (&top,isize[i],index[i],&hb->d,i,
3452
                                 TRUE,(bTwo && (i==gr1)) || !bTwo, datable);
3453
            }
3454
            else {
3455
                search_acceptors(&top,isize[i],index[i],&hb->a,i,bNitAcc,FALSE,TRUE, datable);
3456
                search_donors   (&top,isize[i],index[i],&hb->d,i,FALSE,TRUE, datable);
3457
            }
3458
            if (bTwo)
3459
                clear_datable_grp(datable,top.atoms.nr);
3460
        }
3461
    sfree(datable);
3462
    printf("Found %d donors and %d acceptors\n",hb->d.nrd,hb->a.nra);
3463
    /*if (bSelected)
3464
      snew(donors[gr0D], dons[gr0D].nrd);*/
3465

    
3466
    if (bHBmap) {
3467
        printf("Making hbmap structure...");
3468
        /* Generate hbond data structure */
3469
        mk_hbmap(hb,bTwo);
3470
        printf("done.\n");
3471
    }
3472

    
3473
#ifdef HAVE_NN_LOOPS
3474
    if (bNN)
3475
        mk_hbEmap(hb, 0);
3476
#endif
3477

    
3478
    if (bGem) {
3479
        printf("Making per structure...");
3480
        /* Generate hbond data structure */
3481
        mk_per(hb);
3482
        printf("done.\n");
3483
    }
3484
  
3485
    /* check input */
3486
    bStop=FALSE;
3487
    if (hb->d.nrd + hb->a.nra == 0) {
3488
        printf("No Donors or Acceptors found\n");
3489
        bStop=TRUE;
3490
    }
3491
    if (!bStop) {
3492
        if (hb->d.nrd == 0) {
3493
            printf("No Donors found\n");
3494
            bStop=TRUE;
3495
        }
3496
        if (hb->a.nra == 0) {
3497
            printf("No Acceptors found\n");
3498
            bStop=TRUE;
3499
        }
3500
    }
3501
    if (bStop)
3502
        gmx_fatal(FARGS,"Nothing to be done");
3503

    
3504
    shatom=0;
3505
    if (rshell > 0) {
3506
        int shisz;
3507
        atom_id *shidx;
3508
        char *shgrpnm;
3509
        /* get index group with atom for shell */
3510
        do {
3511
            printf("Select atom for shell (1 atom):\n");
3512
            get_index(&(top.atoms),ftp2fn_null(efNDX,NFILE,fnm),
3513
                      1,&shisz,&shidx,&shgrpnm);
3514
            if (shisz != 1)
3515
                printf("group contains %d atoms, should be 1 (one)\n",shisz);
3516
        } while(shisz != 1);
3517
        shatom = shidx[0];
3518
        printf("Will calculate hydrogen bonds within a shell "
3519
               "of %g nm around atom %i\n",rshell,shatom+1);
3520
    }
3521

    
3522
    /* Analyze trajectory */
3523
    natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
3524
    if ( natoms > top.atoms.nr )
3525
        gmx_fatal(FARGS,"Topology (%d atoms) does not match trajectory (%d atoms)",
3526
                  top.atoms.nr,natoms);
3527

    
3528
    bBox  = ir.ePBC!=epbcNONE;
3529
    grid  = init_grid(bBox, box, (rcut>r2cut)?rcut:r2cut, ngrid);
3530
    nabin = acut/abin;
3531
    nrbin = rcut/rbin;
3532
    snew(adist,nabin+1);
3533
    snew(rdist,nrbin+1);
3534

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

    
3538
    bParallel = FALSE;
3539

    
3540
#ifndef HAVE_OPENMP
3541
#define __ADIST adist
3542
#define __RDIST rdist
3543
#define __HBDATA hb
3544
#else /* HAVE_OPENMP ================================================== \
3545
       * Set up the OpenMP stuff,                                       |
3546
       * like the number of threads and such                            |
3547
       * Also start the parallel loop.                                  |
3548
       */
3549
#define __ADIST p_adist[threadNr]
3550
#define __RDIST p_rdist[threadNr]
3551
#define __HBDATA p_hb[threadNr]
3552

    
3553
    bParallel = !bSelected;
3554

    
3555
    if (bParallel)
3556
    {
3557
#if (_OPENMP > 200805)
3558
        actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_thread_limit());
3559
#else
3560
        actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_num_procs());
3561
#endif
3562
        omp_set_num_threads(actual_nThreads);
3563
        printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads);
3564
        fflush(stdout);
3565
    }
3566
    else
3567
    {
3568
        actual_nThreads = 1;
3569
    }
3570

    
3571
    t_hbdata **p_hb;          /* one per thread, then merge after the frame loop */
3572
    int **p_adist, **p_rdist; /* a histogram for each thread. */
3573
    snew(p_hb,    actual_nThreads);
3574
    snew(p_adist, actual_nThreads);
3575
    snew(p_rdist, actual_nThreads);
3576
    for (i=0; i<actual_nThreads; i++)
3577
    {
3578
        snew(p_hb[i], 1);
3579
        snew(p_adist[i], nabin+1);
3580
        snew(p_rdist[i], nrbin+1);
3581

    
3582
        p_hb[i]->max_frames = 0;
3583
        p_hb[i]->nhb = NULL;
3584
        p_hb[i]->ndist = NULL;
3585
        p_hb[i]->n_bound = NULL;
3586
        p_hb[i]->time = NULL;
3587
        p_hb[i]->nhx = NULL;
3588

    
3589
        p_hb[i]->bHBmap     = hb->bHBmap;
3590
        p_hb[i]->bDAnr      = hb->bDAnr;
3591
        p_hb[i]->bGem       = hb->bGem;
3592
        p_hb[i]->wordlen    = hb->wordlen;
3593
        p_hb[i]->nframes    = hb->nframes;
3594
        p_hb[i]->maxhydro   = hb->maxhydro;
3595
        p_hb[i]->danr       = hb->danr;
3596
        p_hb[i]->d          = hb->d;
3597
        p_hb[i]->a          = hb->a;
3598
        p_hb[i]->hbmap      = hb->hbmap;
3599
        p_hb[i]->time       = hb->time; /* This may need re-syncing at every frame. */
3600
        p_hb[i]->per        = hb->per;
3601

    
3602
#ifdef HAVE_NN_LOOPS
3603
        p_hb[i]->hbE = hb->hbE;
3604
#endif
3605

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

    
3636
#ifdef HAVE_OPENMP
3637
            sync_hbdata(hb, p_hb[threadNr], nframes, t);
3638
#pragma omp single
3639
#endif
3640
            {
3641
                build_grid(hb,x,x[shatom], bBox,box,hbox, (rcut>r2cut)?rcut:r2cut, 
3642
                           rshell, ngrid,grid);
3643
                reset_nhbonds(&(hb->d));
3644
    
3645
                if (debug && bDebug)
3646
                    dump_grid(debug, ngrid, grid);
3647
                
3648
                add_frames(hb,nframes);
3649
                init_hbframe(hb,nframes,output_env_conv_time(oenv,t));    
3650

    
3651
                if (hb->bDAnr)
3652
                    count_da_grid(ngrid, grid, hb->danr[nframes]);
3653
            } /* omp single */
3654

    
3655
#ifdef HAVE_OPENMP
3656
            p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */
3657
#endif
3658
            if (bNN)
3659
            {
3660
#ifdef HAVE_NN_LOOPS /* Unlock this feature when testing */
3661
                /* Loop over all atom pairs and estimate interaction energy */
3662
#ifdef HAVE_OPENMP /* ------- */
3663
#pragma omp single
3664
#endif /* HAVE_OPENMP ------- */
3665
                {
3666
                    addFramesNN(hb, nframes);
3667
                }
3668
#ifdef HAVE_OPENMP /* ---------------- */
3669
#pragma omp barrier
3670
#pragma omp for schedule(dynamic)
3671
#endif /* HAVE_OPENMP ---------------- */
3672
                for (i=0; i<hb->d.nrd; i++)
3673
                {
3674
                    for(j=0;j<hb->a.nra; j++)
3675
                    {
3676
                        for (h=0;
3677
                             h < (bContact ? 1 : hb->d.nhydro[i]);
3678
                             h++)
3679
                        {
3680
                            if (i==hb->d.nrd || j==hb->a.nra)
3681
                                gmx_fatal(FARGS, "out of bounds");
3682

    
3683
                            /* Get the real atom ids */
3684
                            ii = hb->d.don[i];
3685
                            jj = hb->a.acc[j];
3686
                            hh = hb->d.hydro[i][h];
3687
            
3688
                            /* Estimate the energy from the geometry */
3689
                            E = calcHbEnergy(ii, jj, hh, x, NN, box, hbox, &(hb->d));
3690
                            /* Store the energy */
3691
                            storeHbEnergy(hb, i, j, h, E, nframes);
3692
                        }
3693
                    }
3694
                }
3695
#endif /* HAVE_NN_LOOPS */
3696
            } /* if (bNN)*/
3697
            else
3698
            {
3699
                if (bSelected)
3700
                {
3701
#ifdef HAVE_OPENMP
3702
#pragma omp single
3703
#endif
3704
                    {
3705
                        /* Do not parallelize this just yet. */
3706
                        /* int ii; */
3707
                        for(ii=0; (ii<nsel); ii++) {
3708
                            int dd = index[0][i];
3709
                            int aa = index[0][i+2];
3710
                            /* int */ hh = index[0][i+1];
3711
                            ihb = is_hbond(hb,ii,ii,dd,aa,rcut,r2cut,ccut,x,bBox,box,
3712
                                           hbox,&dist,&ang,bDA,&h,bContact,bMerge,&peri);
3713
      
3714
                            if (ihb) {
3715
                                /* add to index if not already there */
3716
                                /* Add a hbond */
3717
                                add_hbond(hb,dd,aa,hh,ii,ii,nframes,bMerge,ihb,bContact,peri);
3718
                            }
3719
                        }
3720
                    } /* omp single */
3721
                } /* if (bSelected) */
3722
                else
3723
                {
3724
#ifdef HAVE_OPENMP
3725
#pragma omp single
3726
                    {
3727
#endif
3728
                    if (bGem)
3729
                        calcBoxProjection(box, hb->per->P);
3730

    
3731
                    /* loop over all gridcells (xi,yi,zi)      */
3732
                    /* Removed confusing macro, DvdS 27/12/98  */
3733
#ifdef HAVE_OPENMP
3734
                    }
3735
                    /* The outer grid loop will have to do for now. */
3736
#pragma omp for schedule(dynamic)
3737
#endif
3738
                    for(xi=0; (xi<ngrid[XX]); xi++)
3739
                        for(yi=0; (yi<ngrid[YY]); yi++)
3740
                            for(zi=0; (zi<ngrid[ZZ]); zi++) {
3741
          
3742
                                /* loop over donor groups gr0 (always) and gr1 (if necessary) */
3743
                                for (grp=gr0; (grp <= (bTwo?gr1:gr0)); grp++) {
3744
                                    icell=&(grid[zi][yi][xi].d[grp]);
3745
        
3746
                                    if (bTwo)
3747
                                        ogrp = 1-grp;
3748
                                    else
3749
                                        ogrp = grp;
3750
        
3751
                                    /* loop over all hydrogen atoms from group (grp) 
3752
                                     * in this gridcell (icell) 
3753
                                     */
3754
                                    for (ai=0; (ai<icell->nr); ai++) {
3755
                                        i  = icell->atoms[ai];
3756
        
3757
                                        /* loop over all adjacent gridcells (xj,yj,zj) */
3758
                                        /* This is a macro!!! */
3759
                                        LOOPGRIDINNER(xj,yj,zj,xjj,yjj,zjj,xi,yi,zi,ngrid,bTric) {
3760
                                            jcell=&(grid[zj][yj][xj].a[ogrp]);
3761
                                            /* loop over acceptor atoms from other group (ogrp) 
3762
                                             * in this adjacent gridcell (jcell) 
3763
                                             */
3764
                                            for (aj=0; (aj<jcell->nr); aj++) {
3765
                                                j = jcell->atoms[aj];
3766
          
3767
                                                /* check if this once was a h-bond */
3768
                                                peri = -1;
3769
                                                ihb = is_hbond(__HBDATA,grp,ogrp,i,j,rcut,r2cut,ccut,x,bBox,box,
3770
                                                               hbox,&dist,&ang,bDA,&h,bContact,bMerge,&peri);
3771
            
3772
                                                if (ihb) {
3773
                                                    /* add to index if not already there */
3774
                                                    /* Add a hbond */
3775
                                                    add_hbond(__HBDATA,i,j,h,grp,ogrp,nframes,bMerge,ihb,bContact,peri);
3776
              
3777
                                                    /* make angle and distance distributions */
3778
                                                    if (ihb == hbHB && !bContact) {
3779
                                                        if (dist>rcut)
3780
                                                            gmx_fatal(FARGS,"distance is higher than what is allowed for an hbond: %f",dist);
3781
                                                        ang*=RAD2DEG;
3782
                                                        __ADIST[(int)( ang/abin)]++;
3783
                                                        __RDIST[(int)(dist/rbin)]++;
3784
                                                        if (!bTwo) {
3785
                                                            int id,ia;
3786
                                                            if ((id = donor_index(&hb->d,grp,i)) == NOTSET)
3787
                                                                gmx_fatal(FARGS,"Invalid donor %d",i);
3788
                                                            if ((ia = acceptor_index(&hb->a,ogrp,j)) == NOTSET)
3789
                                                                gmx_fatal(FARGS,"Invalid acceptor %d",j);
3790
                                                            resdist=abs(top.atoms.atom[i].resind-
3791
                                                                        top.atoms.atom[j].resind);
3792
                                                            if (resdist >= max_hx)
3793
                                                                resdist = max_hx-1;
3794
                                                            __HBDATA->nhx[nframes][resdist]++;
3795
                                                        }
3796
                                                    }
3797

    
3798
                                                }
3799
                                            } /* for aj  */
3800
                                        }
3801
                                        ENDLOOPGRIDINNER;
3802
                                    } /* for ai  */
3803
                                } /* for grp */
3804
                            } /* for xi,yi,zi */
3805
                } /* if (bSelected) {...} else */ 
3806

    
3807
#ifdef HAVE_OPENMP /* ---------------------------- */
3808
                /* Better wait for all threads to finnish using x[] before updating it. */
3809
                k = nframes;            /*         */
3810
#pragma omp barrier                     /*         */
3811
#pragma omp critical                    /*         */
3812
                {                       /*         */
3813
                    /* Sum up histograms and counts from p_hb[] into hb */
3814
                    {                   /*         */
3815
                        hb->nhb[k]   += p_hb[threadNr]->nhb[k];
3816
                        hb->ndist[k] += p_hb[threadNr]->ndist[k];
3817
                        for (j=0; j<max_hx; j++) /**/
3818
                            hb->nhx[k][j]  += p_hb[threadNr]->nhx[k][j];
3819
                    }                   /*         */
3820
                }                       /*         */
3821
                /*                                 */
3822
                /* Here are a handful of single constructs
3823
                 * to share the workload a bit. The most
3824
                 * important one is of course the last one,
3825
                 * where there's a potential bottleneck in form
3826
                 * of slow I/O.                    */
3827
#pragma omp single /* ++++++++++++++++,            */
3828
#endif /* HAVE_OPENMP ----------------+------------*/
3829
                { /*                  +   */
3830
                    if (hb != NULL)  /*   */
3831
                    { /*              +   */
3832
                        analyse_donor_props(opt2fn_null("-don",NFILE,fnm),hb,k,t,oenv);
3833
                    } /*              +   */
3834
                } /*                  +   */
3835
#ifdef HAVE_OPENMP /*                 +   */
3836
#pragma omp single /* +++           +++   */
3837
#endif       /*                       +   */
3838
                {  /*                 +   */
3839
                    if (fpnhb)  /*    +   */
3840
                        do_nhb_dist(fpnhb,hb,t);
3841
                }  /*                 +   */
3842
            } /* if (bNN) {...} else  +   */
3843
#ifdef HAVE_OPENMP /*                 +   */
3844
#pragma omp single /* +++           +++   */
3845
#endif       /*                       +   */
3846
            {      /*                 +   */
3847
                trrStatus = (read_next_x(oenv,status,&t,natoms,x,box));
3848
                nframes++;      /*    +   */
3849
            }      /*                 +   */
3850
#ifdef HAVE_OPENMP /* ++++++++++++++++?   */
3851
#pragma omp barrier
3852
#endif
3853
        } while (trrStatus);
3854

    
3855
#ifdef HAVE_OPENMP
3856
#pragma omp critical
3857
        {
3858
            hb->nrhb += p_hb[threadNr]->nrhb;
3859
            hb->nrdist += p_hb[threadNr]->nrdist;
3860
        }
3861
        /* Free parallel datastructures */
3862
        sfree(p_hb[threadNr]->nhb);
3863
        sfree(p_hb[threadNr]->ndist);
3864
        sfree(p_hb[threadNr]->nhx);
3865

    
3866
#pragma omp for
3867
        for (i=0; i<nabin; i++)
3868
            for (j=0; j<actual_nThreads; j++)
3869

    
3870
                adist[i] += p_adist[j][i];
3871
#pragma omp for
3872
        for (i=0; i<=nrbin; i++)
3873
            for (j=0; j<actual_nThreads; j++)
3874
                rdist[i] += p_rdist[j][i];
3875
    
3876
        sfree(p_adist[threadNr]);
3877
        sfree(p_rdist[threadNr]);
3878
    } /* End of parallel region */
3879
    sfree(p_adist);
3880
    sfree(p_rdist);
3881
#endif
3882
  
3883
    if(nframes <2 && (opt2bSet("-ac",NFILE,fnm) || opt2bSet("-life",NFILE,fnm)))
3884
    {
3885
        gmx_fatal(FARGS,"Cannot calculate autocorrelation of life times with less than two frames");
3886
    }
3887
  
3888
    free_grid(ngrid,&grid);
3889
  
3890
    close_trj(status);
3891
    if (fpnhb)
3892
        ffclose(fpnhb);
3893
    
3894
    /* Compute maximum possible number of different hbonds */
3895
    if (maxnhb > 0)
3896
        max_nhb = maxnhb;
3897
    else {
3898
        max_nhb = 0.5*(hb->d.nrd*hb->a.nra);
3899
    }
3900
    /* Added support for -contact below.
3901
     * - Erik Marklund, May 29-31, 2006 */
3902
    /* Changed contact code.
3903
     * - Erik Marklund, June 29, 2006 */
3904
    if (bHBmap && !bNN) {
3905
        if (hb->nrhb==0) {
3906
            printf("No %s found!!\n", bContact ? "contacts" : "hydrogen bonds");
3907
        } else {
3908
            printf("Found %d different %s in trajectory\n"
3909
                   "Found %d different atom-pairs within %s distance\n",
3910
                   hb->nrhb, bContact?"contacts":"hydrogen bonds",
3911
                   hb->nrdist,(r2cut>0)?"second cut-off":"hydrogen bonding");
3912

    
3913
            /*Control the pHist.*/
3914

    
3915
            if (bMerge)
3916
                merge_hb(hb,bTwo,bContact);
3917

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

    
3921
            /* Moved the call to merge_hb() to a line BEFORE dump_hbmap
3922
             * to make the -hbn and -hmb output match eachother. 
3923
             * - Erik Marklund, May 30, 2006 */
3924
        }
3925
    }
3926
    /* Print out number of hbonds and distances */
3927
    aver_nhb  = 0;    
3928
    aver_dist = 0;
3929
    fp = xvgropen(opt2fn("-num",NFILE,fnm),bContact ? "Contacts" :
3930
                  "Hydrogen Bonds",output_env_get_xvgr_tlabel(oenv),"Number",oenv);
3931
    snew(leg,2);
3932
    snew(leg[0],STRLEN);
3933
    snew(leg[1],STRLEN);
3934
    sprintf(leg[0],"%s",bContact?"Contacts":"Hydrogen bonds");
3935
    sprintf(leg[1],"Pairs within %g nm",(r2cut>0)?r2cut:rcut);
3936
    xvgr_legend(fp,2,(const char**)leg,oenv);
3937
    sfree(leg[1]);
3938
    sfree(leg[0]);
3939
    sfree(leg);
3940
    for(i=0; (i<nframes); i++) {
3941
        fprintf(fp,"%10g  %10d  %10d\n",hb->time[i],hb->nhb[i],hb->ndist[i]);
3942
        aver_nhb  += hb->nhb[i];
3943
        aver_dist += hb->ndist[i];
3944
    }
3945
    ffclose(fp);
3946
    aver_nhb  /= nframes;
3947
    aver_dist /= nframes;
3948
    /* Print HB distance distribution */
3949
    if (opt2bSet("-dist",NFILE,fnm)) {
3950
        long sum;
3951
    
3952
        sum=0;
3953
        for(i=0; i<nrbin; i++)
3954
            sum+=rdist[i];
3955
    
3956
        fp = xvgropen(opt2fn("-dist",NFILE,fnm),
3957
                      "Hydrogen Bond Distribution",
3958
                      "Hydrogen - Acceptor Distance (nm)","",oenv);
3959
        for(i=0; i<nrbin; i++)
3960
            fprintf(fp,"%10g %10g\n",(i+0.5)*rbin,rdist[i]/(rbin*(real)sum));
3961
        ffclose(fp);
3962
    }
3963
  
3964
    /* Print HB angle distribution */
3965
    if (opt2bSet("-ang",NFILE,fnm)) {
3966
        long sum;
3967
    
3968
        sum=0;
3969
        for(i=0; i<nabin; i++)
3970
            sum+=adist[i];
3971
    
3972
        fp = xvgropen(opt2fn("-ang",NFILE,fnm),
3973
                      "Hydrogen Bond Distribution",
3974
                      "Donor - Hydrogen - Acceptor Angle (\\SO\\N)","",oenv);
3975
        for(i=0; i<nabin; i++)
3976
            fprintf(fp,"%10g %10g\n",(i+0.5)*abin,adist[i]/(abin*(real)sum));
3977
        ffclose(fp);
3978
    }
3979
    
3980
    /* Print HB in alpha-helix */
3981
    if (opt2bSet("-hx",NFILE,fnm)) {
3982
        fp = xvgropen(opt2fn("-hx",NFILE,fnm),
3983
                      "Hydrogen Bonds",output_env_get_xvgr_tlabel(oenv),"Count",oenv);
3984
        xvgr_legend(fp,NRHXTYPES,hxtypenames,oenv);
3985
        for(i=0; i<nframes; i++) {
3986
            fprintf(fp,"%10g",hb->time[i]);
3987
            for (j=0; j<max_hx; j++)
3988
                fprintf(fp," %6d",hb->nhx[i][j]);
3989
            fprintf(fp,"\n");
3990
        }
3991
        ffclose(fp);
3992
    }
3993
    if (!bNN)
3994
        printf("Average number of %s per timeframe %.3f out of %g possible\n",
3995
               bContact ? "contacts" : "hbonds",
3996
               bContact ? aver_dist : aver_nhb, max_nhb);
3997
     
3998
    /* Do Autocorrelation etc. */
3999
    if (hb->bHBmap) {
4000
        /* 
4001
           Added support for -contact in ac and hbm calculations below.
4002
           - Erik Marklund, May 29, 2006
4003
        */
4004
        ivec itmp;
4005
        rvec rtmp;
4006
        if (opt2bSet("-ac",NFILE,fnm) || opt2bSet("-life",NFILE,fnm))
4007
            please_cite(stdout,"Spoel2006b");
4008
        if (opt2bSet("-ac",NFILE,fnm)) {
4009
            char *gemstring=NULL;
4010

    
4011
            if (bGem || bNN) {
4012
                params = init_gemParams(rcut, D, hb->time, hb->nframes/2, nFitPoints, fit_start, fit_end,
4013
                                        gemBallistic, nBalExp, bBallisticDt);
4014
                if (params == NULL)
4015
                    gmx_fatal(FARGS, "Could not initiate t_gemParams params.");
4016
            }
4017
            gemstring = strdup(gemType[hb->per->gemtype]);
4018
            do_hbac(opt2fn("-ac",NFILE,fnm),hb,nDump,
4019
                    bMerge,bContact,fit_start,temp,r2cut>0,smooth_tail_start,oenv,
4020
                    params, gemstring, nThreads, NN, bBallistic, bGemFit);
4021
        }
4022
        if (opt2bSet("-life",NFILE,fnm))
4023
            do_hblife(opt2fn("-life",NFILE,fnm),hb,bMerge,bContact,oenv);
4024
        if (opt2bSet("-hbm",NFILE,fnm)) {
4025
            t_matrix mat;
4026
            int id,ia,hh,x,y;
4027
      
4028
            mat.nx=nframes;
4029
            mat.ny=(bContact ? hb->nrdist : hb->nrhb);
4030

    
4031
            snew(mat.matrix,mat.nx);
4032
            for(x=0; (x<mat.nx); x++) 
4033
                snew(mat.matrix[x],mat.ny);
4034
            y=0;
4035
            for(id=0; (id<hb->d.nrd); id++) 
4036
                for(ia=0; (ia<hb->a.nra); ia++) {
4037
                    for(hh=0; (hh<hb->maxhydro); hh++) {
4038
                        if (hb->hbmap[id][ia]) {
4039
                            if (ISHB(hb->hbmap[id][ia]->history[hh])) {
4040
                                /* Changed '<' into '<=' in the for-statement below.
4041
                                 * It fixed the previously undiscovered bug that caused
4042
                                 * the last occurance of an hbond/contact to not be
4043
                                 * set in mat.matrix. Have a look at any old -hbm-output
4044
                                 * and you will notice that the last column is allways empty.
4045
                                 * - Erik Marklund May 30, 2006
4046
                                 */
4047
                                for(x=0; (x<=hb->hbmap[id][ia]->nframes); x++) {
4048
                                    int nn0 = hb->hbmap[id][ia]->n0;
4049
                                    range_check(y,0,mat.ny);
4050
                                    mat.matrix[x+nn0][y] = is_hb(hb->hbmap[id][ia]->h[hh],x);
4051
                                }
4052
                                y++;
4053
                            }
4054
                        }
4055
                    }
4056
                }
4057
            mat.axis_x=hb->time;
4058
            snew(mat.axis_y,mat.ny);
4059
            for(j=0; j<mat.ny; j++)
4060
                mat.axis_y[j]=j;
4061
            sprintf(mat.title,bContact ? "Contact Existence Map":
4062
                    "Hydrogen Bond Existence Map");
4063
            sprintf(mat.legend,bContact ? "Contacts" : "Hydrogen Bonds");
4064
            sprintf(mat.label_x,output_env_get_xvgr_tlabel(oenv));
4065
            sprintf(mat.label_y, bContact ? "Contact Index" : "Hydrogen Bond Index");
4066
            mat.bDiscrete=TRUE;
4067
            mat.nmap=2;
4068
            snew(mat.map,mat.nmap);
4069
            for(i=0; i<mat.nmap; i++) {
4070
                mat.map[i].code.c1=hbmap[i];
4071
                mat.map[i].desc=hbdesc[i];
4072
                mat.map[i].rgb=hbrgb[i];
4073
            }
4074
            fp = opt2FILE("-hbm",NFILE,fnm,"w");
4075
            write_xpm_m(fp, mat);
4076
            ffclose(fp);
4077
            for(x=0; x<mat.nx; x++)
4078
                sfree(mat.matrix[x]);
4079
            sfree(mat.axis_y);
4080
            sfree(mat.matrix);
4081
            sfree(mat.map);
4082
        }
4083
    }
4084

    
4085
    if (bGem) {
4086
        fprintf(stderr, "There were %i periodic shifts\n", hb->per->nper);
4087
        fprintf(stderr, "Freeing pHist for all donors...\n");
4088
        for (i=0; i<hb->d.nrd; i++) {
4089
            fprintf(stderr, "\r%i",i);
4090
            if (hb->per->pHist[i] != NULL) {
4091
                for (j=0; j<hb->a.nra; j++)
4092
                    clearPshift(&(hb->per->pHist[i][j]));
4093
                sfree(hb->per->pHist[i]);
4094
            }
4095
        }
4096
        sfree(hb->per->pHist);
4097
        sfree(hb->per->p2i);
4098
        sfree(hb->per);
4099
        fprintf(stderr, "...done.\n");
4100
    }
4101

    
4102
#ifdef HAVE_NN_LOOPS
4103
    if (bNN)
4104
        free_hbEmap(hb);
4105
#endif
4106
    
4107
    if (hb->bDAnr) {
4108
        int  i,j,nleg;
4109
        char **legnames;
4110
        char buf[STRLEN];
4111
    
4112
#define USE_THIS_GROUP(j) ( (j == gr0) || (bTwo && (j == gr1)) )
4113
    
4114
        fp = xvgropen(opt2fn("-dan",NFILE,fnm),
4115
                      "Donors and Acceptors",output_env_get_xvgr_tlabel(oenv),"Count",oenv);
4116
        nleg = (bTwo?2:1)*2;
4117
        snew(legnames,nleg);
4118
        i=0;
4119
        for(j=0; j<grNR; j++)
4120
            if ( USE_THIS_GROUP(j) ) {
4121
                sprintf(buf,"Donors %s",grpnames[j]);
4122
                legnames[i++] = strdup(buf);
4123
                sprintf(buf,"Acceptors %s",grpnames[j]);
4124
                legnames[i++] = strdup(buf);
4125
            }
4126
        if (i != nleg)
4127
            gmx_incons("number of legend entries");
4128
        xvgr_legend(fp,nleg,(const char**)legnames,oenv);
4129
        for(i=0; i<nframes; i++) {
4130
            fprintf(fp,"%10g",hb->time[i]);
4131
            for (j=0; (j<grNR); j++)
4132
                if ( USE_THIS_GROUP(j) )
4133
                    fprintf(fp," %6d",hb->danr[i][j]);
4134
            fprintf(fp,"\n");
4135
        }
4136
        ffclose(fp);
4137
    }
4138
  
4139
    thanx(stdout);
4140
  
4141
    return 0;
4142
}