gmx_hbond.c
1 
/* * mode: c; tabwidth: 4; indenttabsmode: nil; cbasicoffset: 4; cfilestyle: "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) 19912000, University of Groningen, The Netherlands.

13 
* Copyright (c) 20012008, 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 
{"nn","nn+1","nn+2","nn+3","nn+4","nn+5","nn>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_YESHB_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) = [1h(t)] H(t) where H(t) is one when the donor

143 
* acceptor distance is less than the userspecified 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 dapair.

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][framen0] */

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 hbinfo 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 
/* Selfinteraction */

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,framehb>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>len1] != p) { 
657 
pHist>len++; 
658 
srenew(pHist>frame, pHist>len); 
659 
srenew(pHist>p, pHist>len); 
660 
pHist>frame[pHist>len1] = frame;

661 
pHist>p[pHist>len1] = 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[i1]; 
681 
} 
682 

683 
/* It seems that frame is after the last periodic transition. Return the last periodicity. */

684 
return pHist.p[pHist.len1]; 
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 = framehb>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/atable!\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+2j];

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=nr11;

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 gridseach 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=DIM1; 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=DIM1; 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=DIM1; 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=DIM1; 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/yplane, 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 : (x1)) 
1360 
#define E(n,x,bTric,bEdge) ((n==1) ? x : bTric&&(bEdge) ? n1 : (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==0z==n[ZZ]1); \ 
1366 
yy<=E(n[YY],yo,bTric,z==0z==n[ZZ]1); yy++) { \ 
1367 
y=GRIDMOD(yy,n[YY]); \ 
1368 
for(xx=B(n[XX],xo,bTric,y==0y==n[YY]1z==0z==n[ZZ]1); \ 
1369 
xx<=E(n[XX],xo,bTric,y==0y==n[YY]1z==0z==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=DIM1; 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=DIM1; 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 cutoff.

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+n00nn0; 
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+n00nn0; 
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+n01nn0; 
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 hbmoutputloop 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[jj0]++; 
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 = nframes1;

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("debugac.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=1e16,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(k0.6)+0.7*sqr(kp1.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((N2)/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,1e3); 
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,1e3); 
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[n1],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<n2) && ((t[i0]t[0]) < fit_start); i0++) 
2132 
; 
2133 
if (i0 < n2) { 
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*sc2sqr(scn)); 
2144 
if ((tmp > 0) && (sn2 > 0)) { 
2145 
k = (sn2*sckscn*snk)/tmp; 
2146 
kp = (k*scnsnk)/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("Oneway %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[n1]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<n2); i++) { 
2198 
if ((ct[i] > e_1) && (ct[i+1] <= e_1)) { 
2199 
break;

2200 
} 
2201 
} 
2202 
if (i < n2) { 
2203 
/* Determine tau_relax from linear interpolation */

2204 
tau_rlx = t[i]t[0] + (e_1ct[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<nn1); j++) 
2220 
dydx[j] = (y[j+1]y[j1])/(x[j+1]x[j1]); 
2221 
/* Extrapolate endpoints */

2222 
dydx[0] = 2*dydx[1]  dydx[2]; 
2223 
dydx[nn1] = 2*dydx[nn2]  dydx[nn3]; 
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 = 1e3; 
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 NNloop. 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,bMergebContact,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 donorloop.\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 noone, 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 LevenbergMarquardt and nsimplex solvers that come with

2491 
* a BSDlicence 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 reallocating the arrays

2599 
* when taking on the next donoracceptor 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==(np1)) 
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][jpoff[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 LevenbergMarquardt and nsimplex solvers that come with

2722 
* a BSDlicence 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 hbmoutputloop 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 cutoff is provided, use it,

2813 
* otherwise use g(t) = 1h(t) */

2814 
if (!R2 && bContact)

2815 
gt[j] = 1ihb;

2816 
else

2817 
gt[j] = idist*(1ihb);

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(tail2tail*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)/(1tail);

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 ACFcalulation. 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],nframeshb>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,nhtotnbound);

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_framesnframes)); 
3063 
memset(&(p_hb>ndist[nframes]), 0, sizeof(int) * (p_hb>max_framesnframes)); 
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 resyncing. */

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 nonhydrogen atom.[PAR]",

3091 

3092 
"You need to specify two groups for analysis, which must be either",

3093 
"identical or nonoverlapping. 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 nn+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 nn+3, nn+4 and nn+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 DonorAcceptor (if TRUE) or HydrogenAcceptor (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 cutoff 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 
"Hbonds between the same donor and acceptor, but with different hydrogen are treated as a single Hbond. 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 nonpositive, 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 
/* NNloop? If so, what estimator to use ?*/

3269 
NN = 1;

3270 
/* Outcommented for now DvdS 20100713

3271 
while (NN < NN_NR && gmx_strcasecmp(NNtype[0], NNtype[NN])!=0)

3272 
NN++;

3273 
if (NN == NN_NR)

3274 
gmx_fatal(FARGS, "Invalid NNloop 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 donorH 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 DonorAcceptor 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 
"DonorHydrogenAcceptor",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 nonoverlapping 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 resyncing 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 = 1grp;

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 hbond */

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_hx1;

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 2931, 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 atompairs within %s distance\n",

3910 
hb>nrhb, bContact?"contacts":"hydrogen bonds", 
3911 
hb>nrdist,(r2cut>0)?"second cutoff":"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 alphahelix */

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 forstatement 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 hbmoutput

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 
} 