?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Moved the real code of analysis programs to gmx_* and made g_* wrapper binaries, so
[gromacs.git] / src / tools / gmx_hbond.c
blob85bb37822b3cd143ba6c68b49b2b1ff3cc07580a
1 /*
2 * $Id$
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.2.0
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, 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 * Green Red Orange Magenta Azure Cyan Skyblue
35 */
36 #ifdef HAVE_CONFIG_H
37 #include <config.h>
38 #endif
39 #include <math.h>
40
41 #include "statutil.h"
42 #include "copyrite.h"
43 #include "sysstuff.h"
44 #include "txtdump.h"
45 #include "futil.h"
46 #include "tpxio.h"
47 #include "physics.h"
48 #include "macros.h"
49 #include "fatal.h"
50 #include "index.h"
51 #include "smalloc.h"
52 #include "assert.h"
53 #include "vec.h"
54 #include "xvgr.h"
55 #include "gstat.h"
56 #include "matio.h"
57 #include "string2.h"
58 #include "pbc.h"
59
60 #define max_hx 7
61 typedef int t_hx[max_hx];
62 #define NRHXTYPES max_hx
63 char *hxtypenames[NRHXTYPES]=
64 {"n-n","n-n+1","n-n+2","n-n+3","n-n+4","n-n+5","n-n>6"};
65
66 enum { gr0, gr1, grI, grNR };
67 #define OGRP (gr1-grp)
68
69 static char *grpnames[grNR] = {"0","1","I" };
70
71 #define HB_NO 0
72 #define HB_YES 1<<0
73 #define HB_INS 1<<1
74 #define HB_YESINS HB_YES|HB_INS
75 #define HB_NR (1<<2)
76 #define MAXHYDRO 4
77
78 typedef struct {
79 int nr;
80 int maxnr;
81 atom_id *atoms;
82 } t_ncell;
83
84 typedef struct {
85 t_ncell d[grNR];
86 t_ncell a[grNR];
87 } t_gridcell;
88
89 typedef int t_icell[grNR];
90 typedef atom_id h_id[MAXHYDRO];
91
92 typedef struct {
93 bool bExisted[MAXHYDRO]; /* has this hbond existed ever? */
94 bool bInsert; /* has this hbond been invaded? */
95 /* Bitmask array which tells whether a hbond is present
96 * at a given time. Either of these may be NULL
97 */
98 unsigned int **exist;
99 } t_hbond;
100
101 typedef struct {
102 int nra,max_nra;
103 atom_id *acc; /* Atom numbers of the acceptors */
104 int *aptr; /* Map atom number to acceptor index */
105 } t_acceptors;
106
107 typedef struct {
108 int nrd,max_nrd;
109 int *d; /* Atom numbers of the donors */
110 int *nhydro; /* Number of hydrogens for each donor */
111 h_id *hydro; /* The atom numbers of the hydrogens */
112 int *dptr; /* Map atom number to donor index */
113 } t_donors;
114
115 typedef struct {
116 bool bHBmap,bDAnr;
117 int wordlen;
118 /* The following arrays are nframes long */
119 int nframes,max_frames,maxhydro;
120 int *nhb;
121 real *time;
122 t_icell *danr;
123 t_hx *nhx;
124 /* These structures are initialized from the topology at start up */
125 t_donors *(d[grNR]);
126 t_acceptors *(a[grNR]);
127 /* This holds a matrix with all possible hydrogen bonds */
128 int nd,na,nrhb;
129 t_hbond **hbmap;
130 } t_hbdata;
131
132 static t_hbdata *mk_hbdata(bool bHBmap,bool bDAnr,bool bMerge)
133 {
134 t_hbdata *hb;
135
136 snew(hb,1);
137 hb->wordlen = 8*sizeof(unsigned int);
138 hb->bHBmap = bHBmap;
139 hb->bDAnr = bDAnr;
140 if (bMerge)
141 hb->maxhydro = 1;
142 else
143 hb->maxhydro = MAXHYDRO;
144
145
146 return hb;
147 }
148
149 static void mk_hbmap(t_hbdata *hb,bool bTwo,bool bInsert)
150 {
151 int i,j;
152
153 if (bInsert) {
154 hb->nd = hb->d[grI]->nrd;
155 hb->na = hb->a[grI]->nra;
156 }
157 else {
158 hb->nd = hb->d[gr0]->nrd;
159 hb->na = hb->a[gr0]->nra;
160 if (bTwo) {
161 hb->nd += hb->d[gr1]->nrd;
162 hb->na += hb->a[gr1]->nra;
163 }
164 }
165 fprintf(stderr,"Going to allocate %d kb of memory, and that's only the beginning\n",(int)((hb->nd*hb->na*sizeof(hb->hbmap[0][0]))/1024.0));
166
167 snew(hb->hbmap,hb->nd);
168 for(i=0; (i<hb->nd); i++) {
169 snew(hb->hbmap[i],hb->na);
170 for(j=0; (j<hb->na); j++)
171 snew(hb->hbmap[i][j].exist,hb->maxhydro);
172 }
173 }
174
175 static void add_frames(t_hbdata *hb,int nframes)
176 {
177 int i,j,k,l;
178
179 if (nframes >= hb->max_frames) {
180 hb->max_frames += 12000;
181 srenew(hb->time,hb->max_frames);
182 srenew(hb->nhb,hb->max_frames);
183 srenew(hb->nhx,hb->max_frames);
184 if (hb->bDAnr)
185 srenew(hb->danr,hb->max_frames);
186 if (hb->bHBmap) {
187 for(i=0; (i<hb->nd); i++)
188 for(j=0; (j<hb->na); j++)
189 for(k=0; (k<hb->maxhydro); k++)
190 if (hb->hbmap[i][j].exist[k]) {
191 srenew(hb->hbmap[i][j].exist[k],hb->max_frames/hb->wordlen);
192 for(l=hb->nframes/hb->wordlen; (l<hb->max_frames/hb->wordlen); l++)
193 hb->hbmap[i][j].exist[k][l] = 0;
194 }
195 }
196 }
197 hb->nframes=nframes;
198 }
199
200 #define OFFSET(frame) (frame / 32)
201 #define MASK(frame) (1 << (frame % 32))
202
203 static void set_hb(unsigned int hbexist[],unsigned int frame,bool bValue)
204 {
205 if (bValue)
206 hbexist[OFFSET(frame)] |= MASK(frame);
207 else
208 hbexist[OFFSET(frame)] &= ~MASK(frame);
209 }
210
211 static bool is_hb(unsigned int hbexist[],int frame)
212 {
213 return ((hbexist[OFFSET(frame)] & MASK(frame)) != 0) ? 1 : 0;
214 }
215
216 static void add_hbond(t_hbdata *hb,int d,int a,int h,int grpd,int grpa,
217 int frame,bool bInsert,bool bMerge)
218 {
219 int k,id,ia;
220
221 if ((id = hb->d[grpd]->dptr[d]) == NOTSET)
222 fatal_error(0,"No donor atom %d (%s %d)",d+1,__FILE__,__LINE__);
223 if ((ia = hb->a[grpa]->aptr[a]) == NOTSET)
224 fatal_error(0,"No acceptor atom %d (%s %d)",a+1,__FILE__,__LINE__);
225
226 /* Loop over hydrogens */
227 for(k=0; (k<hb->d[grpd]->nhydro[id]); k++)
228 if (hb->d[grpd]->hydro[id][k] == h)
229 break;
230 if (k == hb->d[grpd]->nhydro[id])
231 fatal_error(0,"Donor %d does not have hydrogen %d %s %d",d+1,h+1,
232 __FILE__,__LINE__);;
233
234 if (bMerge)
235 k = 0;
236 if (hb->bHBmap) {
237 if (!hb->hbmap[id][ia].exist[k])
238 snew(hb->hbmap[id][ia].exist[k],hb->max_frames/hb->wordlen);
239 if (frame >= 0)
240 set_hb(hb->hbmap[id][ia].exist[k],frame,TRUE);
241 }
242 hb->hbmap[id][ia].bInsert = bInsert || hb->hbmap[id][ia].bInsert;
243
244 if (frame >= 0) {
245 hb->nhb[frame]++;
246 if (!hb->hbmap[id][ia].bExisted[k]) {
247 hb->hbmap[id][ia].bExisted[k] = TRUE;
248 hb->nrhb++;
249 }
250 }
251 }
252
253 static bool in_list(atom_id selection,int isize,atom_id *index)
254 {
255 int i;
256 bool bFound;
257
258 bFound=FALSE;
259 for(i=0; (i<isize) && !bFound; i++)
260 if(selection == index[i])
261 bFound=TRUE;
262
263 return bFound;
264 }
265
266 static char *mkatomname(t_atoms *atoms,int i)
267 {
268 static char buf[32];
269 int rnr;
270
271 rnr = atoms->atom[i].resnr;
272 sprintf(buf,"%4s%d%-4s",*atoms->resname[rnr],rnr+1,*atoms->atomname[i]);
273
274 return buf;
275 }
276
277 static void add_acc(t_acceptors *a,int ia)
278 {
279 if (a->nra >= a->max_nra) {
280 a->max_nra += 16;
281 srenew(a->acc,a->max_nra);
282 }
283 a->acc[a->nra++] = ia;
284 }
285
286 static t_acceptors *search_acceptors(t_topology *top,int isize,
287 atom_id *index,bool bNitAcc,
288 bool bContact,bool bDoIt)
289 {
290 t_acceptors *a;
291 int i;
292
293 snew(a,1);
294 for (i=0; (i<top->atoms.nr); i++) {
295 if (bDoIt) {
296 if ((bContact ||
297 (((*top->atoms.atomname[i])[0] == 'O') ||
298 (bNitAcc && ((*top->atoms.atomname[i])[0] == 'N')))) &&
299 in_list(i,isize,index)) {
300 add_acc(a,i);
301 }
302 }
303 }
304 snew(a->aptr,top->atoms.nr);
305 for(i=0; (i<top->atoms.nr); i++)
306 a->aptr[i] = NOTSET;
307 for(i=0; (i<a->nra); i++)
308 a->aptr[a->acc[i]] = i;
309
310 return a;
311 }
312
313 static int acceptor_index(t_acceptors *a,atom_id i)
314 {
315 return a->aptr[i];
316 }
317
318 static int donor_index(t_donors *d,atom_id i)
319 {
320 return d->dptr[i];
321 }
322
323 static void add_h2d(int id,int ih,t_donors *ddd)
324 {
325 int i;
326
327 for(i=0; (i<ddd->nhydro[id]); i++)
328 if (ddd->hydro[id][i] == ih) {
329 fprintf(stderr,"Hm. This isn't first time I find this donor (%d,%d)\n",
330 ddd->d[id],ih);
331 break;
332 }
333 if (i == ddd->nhydro[id]) {
334 if (ddd->nhydro[id] >= MAXHYDRO)
335 fatal_error(0,"Donor %d has more than %d hydrogens! %s %d",
336 ddd->d[id],MAXHYDRO,__FILE__,__LINE__);
337 ddd->hydro[id][i] = ih;
338 ddd->nhydro[id]++;
339 }
340 }
341
342 static void add_dh(t_donors *ddd,int id,int ih)
343 {
344 int i;
345
346 for(i=0; (i<ddd->nrd); i++)
347 if (ddd->d[i] == id) {
348 add_h2d(i,ih,ddd);
349 break;
350 }
351 if (i == ddd->nrd) {
352 if (ddd->nrd >= ddd->max_nrd) {
353 ddd->max_nrd += 128;
354 srenew(ddd->d,ddd->max_nrd);
355 srenew(ddd->nhydro,ddd->max_nrd);
356 srenew(ddd->hydro,ddd->max_nrd);
357 }
358 ddd->d[ddd->nrd] = id;
359 ddd->nhydro[ddd->nrd] = 0;
360 add_h2d(ddd->nrd,ih,ddd);
361 ddd->nrd++;
362 }
363 }
364
365 static t_donors *search_donors(t_topology *top, int isize, atom_id *index,
366 bool bContact,bool bDoIt)
367 {
368 int i,j,nra;
369 t_functype func_type;
370 t_ilist *interaction;
371 atom_id nr1,nr2;
372 bool stop;
373 t_donors *ddd;
374
375 snew(ddd,1);
376
377 if (bContact) {
378 if (bDoIt)
379 for(i=0; (i<isize); i++)
380 add_dh(ddd,index[i],-1);
381 }
382 else {
383 for(func_type=0; (func_type < F_NRE); func_type++) {
384 interaction=&(top->idef.il[func_type]);
385 for(i=0; i < interaction->nr;
386 i+=interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1) {
387 /* next function */
388 if (func_type != top->idef.functype[interaction->iatoms[i]]) {
389 fatal_error(0,"Error in %s,%d func_type %s",__FILE__,__LINE__,
390 interaction_function[func_type].longname);
391 }
392
393 /* check out this functype */
394 if (func_type == F_SETTLE) {
395 nr1=interaction->iatoms[i+1];
396
397 if (in_list(nr1, isize,index)) {
398 if (in_list(nr1+1,isize,index))
399 add_dh(ddd,nr1,nr1+1);
400 if (in_list(nr1+2,isize,index))
401 add_dh(ddd,nr1,nr1+2);
402 }
403 }
404 else if (IS_CHEMBOND(func_type)) {
405 for (j=0; j<2; j++) {
406 nr1=interaction->iatoms[i+1+j];
407 nr2=interaction->iatoms[i+2-j];
408 if ((*top->atoms.atomname[nr1][0] == 'H') &&
409 ((*top->atoms.atomname[nr2][0] == 'O') ||
410 (*top->atoms.atomname[nr2][0] == 'N')) &&
411 in_list(nr1,isize,index) && in_list(nr2,isize,index))
412 add_dh(ddd,nr2,nr1);
413 }
414 }
415 }
416 }
417 for(func_type=0; func_type < F_NRE; func_type++) {
418 interaction=&top->idef.il[func_type];
419 for(i=0; i < interaction->nr;
420 i+=interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1) {
421 /* next function */
422 assert(func_type == top->idef.functype[interaction->iatoms[i]]);
423
424 if ( interaction_function[func_type].flags & IF_DUMMY ) {
425 nr1=interaction->iatoms[i+1];
426 if ( *top->atoms.atomname[nr1][0] == 'H') {
427 nr2=nr1-1;
428 stop=FALSE;
429 while (!stop && ( *top->atoms.atomname[nr2][0] == 'H'))
430 if (nr2)
431 nr2--;
432 else
433 stop=TRUE;
434 if ( !stop && ( ( *top->atoms.atomname[nr2][0] == 'O') ||
435 ( *top->atoms.atomname[nr2][0] == 'N') ) &&
436 in_list(nr1,isize,index) && in_list(nr2,isize,index) )
437 add_dh(ddd,nr2,nr1);
438 }
439 }
440 }
441 }
442 }
443 snew(ddd->dptr,top->atoms.nr);
444 for(i=0; (i<top->atoms.nr); i++)
445 ddd->dptr[i] = NOTSET;
446 for(i=0; (i<ddd->nrd); i++)
447 ddd->dptr[ddd->d[i]] = i;
448
449 return ddd;
450 }
451
452 static t_gridcell ***init_grid(bool bBox,rvec box[],real rcut,ivec ngrid)
453 {
454 t_gridcell ***grid;
455 int i,y,z;
456
457 if (bBox)
458 for(i=0; i<DIM; i++)
459 ngrid[i]=(box[i][i]/(1.2*rcut));
460
461 if ( !bBox || (ngrid[XX]<3) || (ngrid[YY]<3) || (ngrid[ZZ]<3) )
462 for(i=0; i<DIM; i++)
463 ngrid[i]=1;
464 printf("\nWill do grid-seach on %dx%dx%d grid, rcut=%g\n",
465 ngrid[XX],ngrid[YY],ngrid[ZZ],rcut);
466 snew(grid,ngrid[ZZ]);
467 for (z=0; z<ngrid[ZZ]; z++) {
468 snew((grid)[z],ngrid[YY]);
469 for (y=0; y<ngrid[YY]; y++)
470 snew((grid)[z][y],ngrid[XX]);
471 }
472 return grid;
473 }
474
475 static void build_grid(t_hbdata *hb,rvec x[], rvec xshell,
476 bool bBox, matrix box, rvec hbox,
477 real rcut, real rshell,
478 ivec ngrid, t_gridcell ***grid)
479 {
480 int i,m,gr,xi,yi,zi,nr;
481 atom_id *ad;
482 ivec grididx;
483 rvec invdelta,dshell;
484 t_ncell *newgrid;
485 bool bDoRshell,bInShell,bAcc;
486 real rshell2=0;
487
488 bDoRshell = (rshell > 0);
489 rshell2 = sqr(rshell);
490 bInShell = TRUE;
491
492 for(m=0; m<DIM; m++) {
493 hbox[m]=box[m][m]*0.5;
494 if (bBox) {
495 invdelta[m]=ngrid[m]/box[m][m];
496 if (1/invdelta[m] < rcut)
497 fatal_error(0,"box shrank too much to keep using this grid\n");
498 } else
499 invdelta[m]=0;
500 }
501 grididx[XX]=0;
502 grididx[YY]=0;
503 grididx[ZZ]=0;
504 for(gr=0; (gr<grNR); gr++) {
505 /* resetting atom counts */
506 for (zi=0; zi<ngrid[ZZ]; zi++)
507 for (yi=0; yi<ngrid[YY]; yi++)
508 for (xi=0; xi<ngrid[XX]; xi++) {
509 grid[zi][yi][xi].d[gr].nr=0;
510 grid[zi][yi][xi].a[gr].nr=0;
511 }
512 /*if (debug)
513 fprintf(debug,"Putting %d %s atoms into grid\n",nr[gr],grpnames[gr]);*/
514 /* put atoms in grid cells */
515 for(bAcc=FALSE; (bAcc<=TRUE); bAcc++) {
516 if (bAcc) {
517 nr = hb->a[gr]->nra;
518 ad = hb->a[gr]->acc;
519 }
520 else {
521 nr = hb->d[gr]->nrd;
522 ad = hb->d[gr]->d;
523 }
524 for(i=0; (i<nr); i++) {
525 /* check if we are inside the shell */
526 /* if bDoRshell=FALSE then bInShell=TRUE always */
527 if ( bDoRshell ) {
528 bInShell=TRUE;
529 rvec_sub(x[ad[i]],xshell,dshell);
530 if (bBox)
531 for(m=DIM-1; m>=0 && bInShell; m--) {
532 if ( dshell[m] < -hbox[m] )
533 rvec_inc(dshell,box[m]);
534 else if ( dshell[m] >= hbox[m] )
535 dshell[m] -= 2*hbox[m];
536 /* if we're outside the cube, we're outside the sphere also! */
537 if ( (dshell[m]>rshell) || (-dshell[m]>rshell) )
538 bInShell=FALSE;
539 }
540 /* if we're inside the cube, check if we're inside the sphere */
541 if (bInShell)
542 bInShell = norm2(dshell) < rshell2;
543 }
544 if ( bInShell ) {
545 if (bBox)
546 for(m=DIM-1; m>=0; m--) {
547 /* put atom in the box */
548 while( x[ad[i]][m] < 0 )
549 rvec_inc(x[ad[i]],box[m]);
550 while( x[ad[i]][m] >= box[m][m] )
551 rvec_dec(x[ad[i]],box[m]);
552 /* determine grid index of atom */
553 grididx[m]=x[ad[i]][m]*invdelta[m];
554 }
555 /* add atom to grid cell */
556 if (bAcc)
557 newgrid=&(grid[grididx[ZZ]][grididx[YY]][grididx[XX]].a[gr]);
558 else
559 newgrid=&(grid[grididx[ZZ]][grididx[YY]][grididx[XX]].d[gr]);
560 if (newgrid->nr >= newgrid->maxnr) {
561 newgrid->maxnr+=10;
562 srenew(newgrid->atoms, newgrid->maxnr);
563 }
564 newgrid->atoms[newgrid->nr]=ad[i];
565 newgrid->nr++;
566 }
567 }
568 }
569 }
570 }
571
572 static void count_da_grid(ivec ngrid, t_gridcell ***grid, t_icell danr)
573 {
574 int gr,xi,yi,zi;
575
576 for(gr=0; gr<grNR; gr++) {
577 danr[gr]=0;
578 for (zi=0; zi<ngrid[ZZ]; zi++)
579 for (yi=0; yi<ngrid[YY]; yi++)
580 for (xi=0; xi<ngrid[XX]; xi++) {
581 danr[gr] += grid[zi][yi][xi].d[gr].nr;
582 }
583 }
584 }
585
586 /* The grid loop.
587 * Without a box, the grid is 1x1x1, so all loops are 1 long.
588 * With a rectangular box (bTric==FALSE) all loops are 3 long.
589 * With a triclinic box all loops are 3 long, except when a cell is
590 * located next to one of the box edges which is not parallel to the
591 * x/y-plane, in that case all cells in a line or layer are searched.
592 * This could be implemented slightly more efficient, but the code
593 * would get much more complicated.
594 */
595 #define B(n,x,bTric,bEdge) ((n==1) ? x : bTric&&(bEdge) ? 0 : (x-1))
596 #define E(n,x,bTric,bEdge) ((n==1) ? x : bTric&&(bEdge) ? n-1 : (x+1))
597 #define GRIDMOD(j,n) (j+n)%(n)
598 #define LOOPGRIDINNER(x,y,z,xx,yy,zz,xo,yo,zo,n,bTric)\
599 for(zz=B(n[ZZ],zo,bTric,FALSE); zz<=E(n[ZZ],zo,bTric,FALSE); zz++) {\
600 z=GRIDMOD(zz,n[ZZ]);\
601 for(yy=B(n[YY],yo,bTric,z==0||z==n[ZZ]-1); \
602 yy<=E(n[YY],yo,bTric,z==0||z==n[ZZ]-1); yy++) {\
603 y=GRIDMOD(yy,n[YY]);\
604 for(xx=B(n[XX],xo,bTric,y==0||y==n[YY]-1||z==0||z==n[ZZ]-1); \
605 xx<=E(n[XX],xo,bTric,y==0||y==n[YY]-1||z==0||z==n[ZZ]-1); xx++) {\
606 x=GRIDMOD(xx,n[XX]);
607 #define ENDLOOPGRIDINNER\
608 }\
609 }\
610 }\
611
612 static void dump_grid(FILE *fp, ivec ngrid, t_gridcell ***grid)
613 {
614 int gr,x,y,z,sum[grNR];
615
616 fprintf(fp,"grid %dx%dx%d\n",ngrid[XX],ngrid[YY],ngrid[ZZ]);
617 for (gr=0; gr<grNR; gr++) {
618 sum[gr]=0;
619 fprintf(fp,"GROUP %d (%s)\n",gr,grpnames[gr]);
620 for (z=0; z<ngrid[ZZ]; z+=2) {
621 fprintf(fp,"Z=%d,%d\n",z,z+1);
622 for (y=0; y<ngrid[YY]; y++) {
623 for (x=0; x<ngrid[XX]; x++) {
624 fprintf(fp,"%3d",grid[x][y][z].d[gr].nr);
625 sum[gr]+=grid[z][y][x].d[gr].nr;
626 fprintf(fp,"%3d",grid[x][y][z].a[gr].nr);
627 sum[gr]+=grid[z][y][x].a[gr].nr;
628
629 }
630 fprintf(fp," | ");
631 if ( (z+1) < ngrid[ZZ] )
632 for (x=0; x<ngrid[XX]; x++) {
633 fprintf(fp,"%3d",grid[z+1][y][x].d[gr].nr);
634 sum[gr]+=grid[z+1][y][x].d[gr].nr;
635 fprintf(fp,"%3d",grid[z+1][y][x].a[gr].nr);
636 sum[gr]+=grid[z+1][y][x].a[gr].nr;
637 }
638 fprintf(fp,"\n");
639 }
640 }
641 }
642 fprintf(fp,"TOTALS:");
643 for (gr=0; gr<grNR; gr++)
644 fprintf(fp," %d=%d",gr,sum[gr]);
645 fprintf(fp,"\n");
646 }
647
648 /* New GMX record! 5 * in a row. Congratulations!
649 * Sorry, only four left.
650 */
651 static void free_grid(ivec ngrid, t_gridcell ****grid)
652 {
653 int y,z;
654 t_gridcell ***g = *grid;
655
656 for (z=0; z<ngrid[ZZ]; z++) {
657 for (y=0; y<ngrid[YY]; y++) {
658 sfree(g[z][y]);
659 }
660 sfree(g[z]);
661 }
662 sfree(g);
663 g=NULL;
664 }
665
666 static void pbc_correct(rvec dx,matrix box,rvec hbox)
667 {
668 int m;
669
670 for(m=DIM-1; m>=0; m--) {
671 if ( dx[m] < -hbox[m] )
672 rvec_inc(dx,box[m]);
673 else if ( dx[m] >= hbox[m] )
674 rvec_dec(dx,box[m]);
675 }
676 }
677
678 static bool low_is_hbond(int d,int a,int h,
679 rvec r_da,real rcut2, real ccut,
680 rvec x[], bool bBox, matrix box,rvec hbox,
681 real *d_ha, real *ang,real d2)
682 {
683 rvec r_dh,r_ha;
684 real ca;
685
686 if (d2 == 0) {
687 rvec_sub(x[h],x[a],r_ha);
688 if (bBox)
689 pbc_correct(r_ha,box,hbox);
690 d2 = iprod(r_ha,r_ha);
691 }
692
693 if ( d2 <= rcut2 ) {
694 rvec_sub(x[d],x[h],r_dh);
695 if (bBox)
696 pbc_correct(r_dh,box,hbox);
697
698 ca = cos_angle(r_dh,r_da);
699 /* if angle is smaller, cos is larger */
700 if (ca >= ccut) {
701 *d_ha = sqrt(d2);
702 *ang = acos(ca);
703 return TRUE;
704 }
705 }
706 return FALSE;
707 }
708
709 static bool is_hbond(t_hbdata *hb,int grpd,int grpa,int d,int a,
710 real rcut, real ccut,
711 rvec x[], bool bBox, matrix box,rvec hbox,
712 real *d_ha, real *ang,bool bDA,int *hhh,
713 bool bContact)
714 {
715 int h,hh,id,ja;
716 rvec r_da;
717 real rc2,d2;
718
719 if (d == a)
720 return FALSE;
721
722 id = donor_index(hb->d[grpd],d);
723 ja = acceptor_index(hb->a[grpa],a);
724
725 rvec_sub(x[d],x[a],r_da);
726 if (bBox)
727 pbc_correct(r_da,box,hbox);
728 d2 = iprod(r_da,r_da);
729 rc2 = rcut*rcut;
730
731 if (d2 < rc2) {
732 if (bContact) {
733 *hhh = -1;
734 return TRUE;
735 }
736 if (!bDA)
737 d2 = 0;
738
739 for(h=0; (h<hb->d[grpd]->nhydro[id]); h++) {
740 hh = hb->d[grpd]->hydro[id][h];
741 if (low_is_hbond(d,a,hh,r_da,rc2,ccut,x,bBox,box,hbox,d_ha,ang,d2)) {
742 *hhh = hh;
743 return TRUE;
744 }
745 }
746 }
747 return FALSE;
748 }
749
750 static void merge_hb(t_hbdata *hb,bool bTwo)
751 {
752 int i,inew,j,ii,jj,m,id,ia,ia0,id0,grp,ogrp,itest;
753 int *gNrD,*gNrA;
754
755 inew = hb->nrhb;
756
757 /* Check whether donors are also acceptors */
758 fprintf(stderr,"Merging hbonds with Acceptor and Donor swapped\n");
759 /* In the two arrays below we store the group number+1 (i.e. either
760 * gr0 or gr1). This is necessary in case we have two groups.
761 */
762 snew(gNrD,hb->na);
763 snew(gNrA,hb->nd);
764 ia0 = id0 = 0;
765 for(grp=gr0; (grp <= bTwo?gr1:gr0); grp++) {
766 for(i=0; (i<hb->d[grp]->nrd); i++) {
767 id = hb->d[grp]->d[i];
768 ia = hb->a[grp]->aptr[id];
769 gNrA[ia0++] = (ia != NOTSET) ? 1+grp : 0;
770 }
771 for(i=0; (i<hb->a[grp]->nra); i++) {
772 ia = hb->a[grp]->acc[i];
773 id = hb->d[grp]->dptr[ia];
774 gNrD[id0++] = (id != NOTSET) ? 1+grp : 0;
775 }
776 }
777 /* Consistency check */
778 if ((ia0 != hb->nd) || (id0 != hb->na))
779 fatal_error(0,"Unexepected inconsistency: ia0=%d, na=%d, id0=%d, nd=%d (%s,%d)",ia0,hb->na,id0,hb->nd,__FILE__,__LINE__);
780
781 itest = 0;
782 for(grp=gr0; (grp <= bTwo?gr1:gr0); grp++) {
783 ogrp = bTwo ? OGRP : grp;
784 for(i=0; (i<hb->nd); i++) {
785 if (gNrA[i] == (1+grp)) {
786 id = hb->d[grp]->d[i];
787 ii = hb->a[grp]->aptr[id];
788 for(j=0; (j<hb->na); j++) {
789 if (gNrD[j] == (1+ogrp)) {
790 ia = hb->a[ogrp]->acc[j];
791 jj = hb->d[ogrp]->dptr[ia];
792 if (id != ia) {
793 if (hb->hbmap[i][j].exist[0] && hb->hbmap[jj][ii].exist[0]) {
794 for(m=0; (m<hb->max_frames/hb->wordlen); m++)
795 hb->hbmap[i][j].exist[0][m] |= hb->hbmap[jj][ii].exist[0][m];
796 sfree(hb->hbmap[jj][ii].exist[0]);
797 hb->hbmap[jj][ii].exist[0] = NULL;
798 inew--;
799 }
800 itest++;
801 }
802 }
803 }
804 }
805 }
806 }
807 sfree(gNrD);
808 sfree(gNrA);
809 fprintf(stderr,"Reduced number of hbonds from %d to %d\n",hb->nrhb,inew);
810 fprintf(stderr,"tested %d pairs\n",itest);
811 hb->nrhb = inew;
812 }
813
814 static void do_hbac(char *fn,t_hbdata *hb,real aver_nhb,bool bDump,bool bMerge)
815 {
816 FILE *fp;
817 int i,j,k,m,ihb;
818 bool bNorm=FALSE;
819 double nhb = 0;
820 real *rhbex;
821 real *ct,tail,tail2,dtail,ct0;
822 const real tol = 1e-3;
823 int nframes = hb->nframes;
824 unsigned int **exist;
825 int nh,nhbonds,nhydro;
826
827 /* build hbexist matrix in reals for autocorr */
828 /* Allocate memory for computing ACF (rhbex) and aggregating the ACF (ct) */
829 snew(rhbex,nframes);
830 snew(ct,nframes);
831 snew(exist,hb->maxhydro);
832
833 /* Loop over hbonds */
834 if (bDump) {
835 fp = ffopen("debug-ac.xvg","w");
836 for(j=0; (j<nframes); j++) {
837 fprintf(fp,"%10g",hb->time[j]);
838 for(i=0; (i<hb->nd); i++) {
839 for(k=0; (k<hb->na); k++) {
840 if (bMerge) {
841 if (hb->hbmap[i][k].exist[0])
842 ihb = is_hb(hb->hbmap[i][k].exist[0],j);
843 else
844 ihb = 0;
845 }
846 else {
847 ihb = 0;
848 for(m=0; (m<hb->maxhydro) && !ihb ; m++)
849 ihb = ihb || ((hb->hbmap[i][k].exist[m]) &&
850 is_hb(hb->hbmap[i][k].exist[m],j));
851 }
852 fprintf(fp," %3d",ihb);
853 }
854 fprintf(fp,"\n");
855 }
856 ffclose(fp);
857 }
858 }
859 /* Total number of hbonds analyzed here */
860 nhbonds = 0;
861 for(i=0; (i<hb->nd); i++) {
862 for(k=0; (k<hb->na); k++) {
863 if (bMerge) {
864 if (hb->hbmap[i][k].exist[0]) {
865 exist[0] = hb->hbmap[i][k].exist[0];
866 nhydro = 1;
867 }
868 else
869 nhydro = 0;
870 }
871 else {
872 nhydro = 0;
873 for(m=0; (m<hb->maxhydro); m++)
874 if (hb->hbmap[i][k].exist[m]) {
875 exist[nhydro++] = hb->hbmap[i][k].exist[m];
876 }
877 }
878 for(nh=0; (nh<nhydro); nh++) {
879 if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == hb->nrhb))
880 fprintf(stderr,"\rACF %d/%d",nhbonds+1,hb->nrhb);
881 nhbonds++;
882 for(j=0; (j<nframes); j++) {
883 ihb = is_hb(exist[nh],j);
884 rhbex[j] = ihb-aver_nhb;
885 nhb += ihb;
886 }
887
888 /* The autocorrelation function is normalized after summation only */
889 low_do_autocorr(NULL,NULL,
890 nframes,1,-1,&rhbex,hb->time[1]-hb->time[0],eacNormal,1,
891 FALSE,TRUE,bNorm,FALSE,0,-1,0,1);
892 for(j=0; (j<nframes/2); j++) {
893 ct[j] += rhbex[j];
894 }
895 }
896 }
897 }
898 fprintf(stderr,"\n");
899 sfree(rhbex);
900
901 /* Normalize */
902 ct0 = ct[0];
903 for(j=0; (j<nframes/2); j++)
904 ct[j] /= ct0;
905
906 /* Determine tail value for statistics */
907 tail = 0;
908 tail2 = 0;
909 for(j=nframes/4; (j<nframes/2); j++) {
910 tail += ct[j];
911 tail2 += ct[j]*ct[j];
912 }
913 tail /= (nframes/2 - nframes/4);
914 tail2 /= (nframes/2 - nframes/4);
915 dtail = sqrt(tail2-tail*tail);
916
917 /* Check whether the ACF is long enough */
918 if (dtail > tol) {
919 printf("\nWARNING: Correlation function is probably not long enough\n"
920 "because the standard deviation in the tail of C(t) > %g\n"
921 "Total number of hbonds found: %g\n"
922 "Tail value (average C(t) over second half of acf): %g +/- %g\n",
923 tol,nhb,tail,dtail);
924 }
925 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation","Time (ps)","C(t)");
926 for(j=0; (j<nframes/2); j++)
927 fprintf(fp,"%10g %10g %10g\n",
928 hb->time[j]-hb->time[0],(ct[j]-tail)/(1-tail),ct[j]);
929
930 fclose(fp);
931 do_view(fn,NULL);
932 sfree(ct);
933 }
934
935 static void init_hbframe(t_hbdata *hb,int nframes,real t)
936 {
937 int i,j,m;
938
939 hb->time[nframes] = t;
940 hb->nhb[nframes] = 0;
941 for (i=0; (i<max_hx); i++)
942 hb->nhx[nframes][i]=0;
943 if (hb->bHBmap)
944 for (i=0; (i<hb->nd); i++)
945 for (j=0; (j<hb->na); j++)
946 for (m=0; (m<hb->maxhydro); m++)
947 if (hb->hbmap[i][j].exist[m])
948 set_hb(hb->hbmap[i][j].exist[m],nframes,HB_NO);
949 }
950
951
952 int gmx_hbond(int argc,char *argv[])
953 {
954 static char *desc[] = {
955 "g_hbond computes and analyzes hydrogen bonds. Hydrogen bonds are",
956 "determined based on cutoffs for the angle Donor - Hydrogen - Acceptor",
957 "(zero is extended) and the distance Hydrogen - Acceptor.",
958 "OH and NH groups are regarded as donors, O is an acceptor always,",
959 "N is an acceptor by default, but this can be switched using",
960 "[TT]-nitacc[tt]. Dummy hydrogen atoms are assumed to be connected",
961 "to the first preceding non-hydrogen atom.[PAR]",
962
963 "You need to specify two groups for analysis, which must be either",
964 "identical or non-overlapping. All hydrogen bonds between the two",
965 "groups are analyzed.[PAR]",
966
967 "If you set -shell, you will be asked for an additional index group",
968 "which should contain exactly one atom. In this case, only hydrogen",
969 "bonds between atoms within the shell distance from the one atom are",
970 "considered.[PAR]"
971
972 "It is also possible to analyse specific hydrogen bonds with",
973 "[TT]-sel[tt]. This index file must contain a group of atom triplets",
974 "Donor Hydrogen Acceptor, in the following way:[PAR]",
975
976 "[TT]",
977 "[ selected ][BR]",
978 " 20 21 24[BR]",
979 " 25 26 29[BR]",
980 " 1 3 6[BR]",
981 "[tt][BR]",
982 "Note that the triplets need not be on separate lines.",
983 "Each atom triplet specifies a hydrogen bond to be analyzed,",
984 "note also that no check is made for the types of atoms.[PAR]",
985
986 "[TT]-ins[tt] turns on computing solvent insertion into hydrogen bonds.",
987 "In this case an additional group must be selected, specifying the",
988 "solvent molecules.[PAR]",
989
990 "[BB]Output:[bb][BR]",
991 "[TT]-num[tt]: number of hydrogen bonds as a function of time.[BR]",
992 "[TT]-ac[tt]: average over all autocorrelations of the existence",
993 "functions (either 0 or 1) of all hydrogen bonds.[BR]",
994 "[TT]-dist[tt]: distance distribution of all hydrogen bonds.[BR]",
995 "[TT]-ang[tt]: angle distribution of all hydrogen bonds.[BR]",
996 "[TT]-hx[tt]: the number of n-n+i hydrogen bonds as a function of time",
997 "where n and n+i stand for residue numbers and i ranges from 0 to 6.",
998 "This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated",
999 "with helices in proteins.[BR]",
1000 "[TT]-hbn[tt]: all selected groups, donors, hydrogens and acceptors",
1001 "for selected groups, all hydrogen bonded atoms from all groups and",
1002 "all solvent atoms involved in insertion.[BR]",
1003 "[TT]-hbm[tt]: existence matrix for all hydrogen bonds over all",
1004 "frames, this also contains information on solvent insertion",
1005 "into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]",
1006 "index file.[BR]",
1007 "[TT]-dan[tt]: write out the number of donors and acceptors analyzed for",
1008 "each timeframe. This is especially usefull when using [TT]-shell[tt]."
1009 };
1010
1011 static real acut=30, abin=1, rcut=0.35, rbin=0.005, rshell=-1,maxnhb=0;
1012 static bool bNitAcc=TRUE,bInsert=FALSE,bDA=TRUE,bDump=FALSE,bMerge=TRUE;
1013 static bool bContact=FALSE;
1014 /* options */
1015 t_pargs pa [] = {
1016 { "-ins", FALSE, etBOOL, {&bInsert},
1017 "Analyze solvent insertion" },
1018 { "-a", FALSE, etREAL, {&acut},
1019 "Cutoff angle (degrees, Donor - Hydrogen - Acceptor)" },
1020 { "-r", FALSE, etREAL, {&rcut},
1021 "Cutoff radius (nm, X - Acceptor, see next option)" },
1022 { "-da", FALSE, etBOOL, {&bDA},
1023 "Use distance Donor-Acceptor (if TRUE) or Hydrogen-Acceptor (FALSE)" },
1024 { "-abin", FALSE, etREAL, {&abin},
1025 "Binwidth angle distribution (degrees)" },
1026 { "-rbin", FALSE, etREAL, {&rbin},
1027 "Binwidth distance distribution (nm)" },
1028 { "-nitacc",FALSE, etBOOL, {&bNitAcc},
1029 "Regard nitrogen atoms as acceptors" },
1030 { "-contact",FALSE,etBOOL, {&bContact},
1031 "Do not look for hydrogen bonds, but merely for contacts within the cut-off distance" },
1032 { "-shell", FALSE, etREAL, {&rshell},
1033 "when > 0, only calculate hydrogen bonds within # nm shell around "
1034 "one particle" },
1035 { "-dump", FALSE, etBOOL, {&bDump},
1036 "Dump all hydrogen bond ACFs (maximum 1000) in a single xvg file for debugging" },
1037 { "-max_hb",FALSE, etREAL, {&maxnhb},
1038 "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" },
1039 { "-merge", FALSE, etBOOL, {&bMerge},
1040 "H-bonds between the same donor and accepter, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." }
1041 };
1042
1043 t_filenm fnm[] = {
1044 { efTRX, "-f", NULL, ffREAD },
1045 { efTPX, NULL, NULL, ffREAD },
1046 { efNDX, NULL, NULL, ffOPTRD },
1047 { efLOG, "-g", "hbond", ffOPTWR },
1048 { efNDX, "-sel", "select", ffOPTRD },
1049 { efXVG, "-num", "hbnum", ffWRITE },
1050 { efXVG, "-ac", "hbac", ffOPTWR },
1051 { efXVG, "-dist","hbdist", ffOPTWR },
1052 { efXVG, "-ang", "hbang", ffOPTWR },
1053 { efXVG, "-hx", "hbhelix",ffOPTWR },
1054 { efNDX, "-hbn", "hbond", ffOPTWR },
1055 { efXPM, "-hbm", "hbmap", ffOPTWR },
1056 { efXVG, "-dan", "danum", ffOPTWR }
1057 };
1058 #define NFILE asize(fnm)
1059
1060 char hbmap [HB_NR]={ ' ', 'o', '-', '*' };
1061 char *hbdesc[HB_NR]={ "None", "Present", "Inserted", "Present & Inserted" };
1062 t_rgb hbrgb [HB_NR]={ {1,1,1},{1,0,0}, {0,0,1}, {1,0,1} };
1063
1064 int status;
1065 t_topology top;
1066 t_inputrec ir;
1067 int natoms,nframes=0,shatom;
1068 int *isize;
1069 char **grpnames;
1070 atom_id **index;
1071 rvec *x,hbox;
1072 matrix box;
1073 real t,ccut,dist,ang;
1074 double max_nhb,aver_nhb;
1075 int h,i,j,k,l,start,end,id,ja;
1076 int xi,yi,zi,ai;
1077 int xj,yj,zj,aj,xjj,yjj,zjj;
1078 int xk,yk,zk,ak,xkk,ykk,zkk;
1079 bool bSelected,bStop,bTwo,new,was,bBox,bTric;
1080 int *adist,*rdist;
1081 int grp,nabin,nrbin,bin,resdist;
1082 t_hbdata *hb;
1083 FILE *fp,*fpins=NULL,*fplog;
1084 t_gridcell ***grid;
1085 t_ncell *icell,*jcell,*kcell;
1086 ivec ngrid;
1087
1088 CopyRight(stderr,argv[0]);
1089 parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,NFILE,fnm,asize(pa),
1090 pa,asize(desc),desc,0,NULL);
1091
1092 /* process input */
1093 bSelected = opt2bSet("-sel",NFILE,fnm);
1094 ccut = cos(acut*DEG2RAD);
1095
1096 if (bContact) {
1097 if (bSelected)
1098 fatal_error(0,"Can not analyze selected contacts: turn off -sel");
1099 if (bInsert)
1100 fatal_error(0,"Can not analyze inserted contacts: turn off -ins");
1101 }
1102
1103 /* Initiate main data structure! */
1104 hb = mk_hbdata(opt2bSet("-hbm",NFILE,fnm) || opt2bSet("-ac",NFILE,fnm),
1105 opt2bSet("-da",NFILE,fnm),bMerge);
1106
1107 /* get topology */
1108 read_tpx(ftp2fn(efTPX,NFILE,fnm),&i,&t,&t,
1109 &ir,box,&natoms,NULL,NULL,NULL,&top);
1110
1111 snew(grpnames,grNR);
1112 snew(index,grNR);
1113 snew(isize,grNR);
1114 if (bSelected) {
1115 /* analyze selected hydrogen bonds */
1116 fprintf(stderr,"Select group with selected atoms:\n");
1117 get_index(&(top.atoms),opt2fn("-sel",NFILE,fnm),
1118 1,isize,index,grpnames);
1119 if (isize[0] % 3)
1120 fatal_error(0,"Number of atoms in group '%s' not a multiple of 3\n"
1121 "and therefore cannot contain triplets of "
1122 "Donor-Hydrogen-Acceptor",grpnames[0]);
1123 bTwo=FALSE;
1124
1125 snew(hb->d[gr0],1);
1126 snew(hb->a[gr0],1);
1127 for(i=0; (i<isize[0]); i+=3) {
1128 int dd = index[0][i];
1129 int hh = index[0][i+1];
1130 int aa = index[0][i+2];
1131 add_dh (hb->d[gr0],dd,hh);
1132 add_acc(hb->a[gr0],aa);
1133 /* Should this be here ? */
1134 add_hbond(hb,dd,hh,aa,gr0,gr0,-1,FALSE,bMerge);
1135 }
1136 fprintf(stderr,"Analyzing %d selected hydrogen bonds from '%s'\n",
1137 hb->nrhb,grpnames[0]);
1138 }
1139 else {
1140 /* analyze all hydrogen bonds: get group(s) */
1141 fprintf(stderr,"Specify 2 groups to analyze:\n");
1142 get_index(&(top.atoms),ftp2fn_null(efNDX,NFILE,fnm),
1143 2,isize,index,grpnames);
1144
1145 /* check if we have two identical or two non-overlapping groups */
1146 bTwo = isize[0] != isize[1];
1147 for (i=0; (i<isize[0]) && !bTwo; i++)
1148 bTwo = index[0][i] != index[1][i];
1149 if (bTwo) {
1150 fprintf(stderr,"Checking for overlap...\n");
1151 for (i=0; i<isize[0]; i++)
1152 for (j=0; j<isize[1]; j++)
1153 if (index[0][i] == index[1][j])
1154 fatal_error(0,"Partial overlap between groups '%s' and '%s'",
1155 grpnames[0],grpnames[1]);
1156 }
1157 if (bContact && !bTwo)
1158 fatal_error(0,"I need two different groups for a contact analysis");
1159 if (bTwo)
1160 fprintf(stderr,"Calculating %s "
1161 "between two groups of %d and %d atoms\n",
1162 bContact ? "contacts" : "hydrogen bonds",isize[0],isize[1]);
1163 else
1164 fprintf(stderr,"Calculating hydrogen bonds in one group of %d atoms\n",
1165 isize[0]);
1166 }
1167 if (bInsert) {
1168 fprintf(stderr,"Specify group for insertion analysis:\n");
1169 get_index(&(top.atoms),ftp2fn_null(efNDX,NFILE,fnm),
1170 1,&(isize[grI]),&(index[grI]),&(grpnames[grI]));
1171 fprintf(stderr,"Checking for overlap...\n");
1172 for (i=0; i<isize[grI]; i++)
1173 for (grp=0; grp<(bTwo?2:1); grp++)
1174 for (j=0; j<isize[grp]; j++)
1175 if (index[grI][i] == index[grp][j])
1176 fatal_error(0,"Partial overlap between groups '%s' and '%s'",
1177 grpnames[grp],grpnames[grI]);
1178 fpins=ffopen("insert.dat","w");
1179 fprintf(fpins,"%4s: %15s -> %15s (%7s) - %15s (%7s)\n",
1180 "time","insert","donor","distang","acceptor","distang");
1181 }
1182
1183 /* search donors and acceptors in groups */
1184 for (i=0; (i<grNR); i++)
1185 if ( ((i==gr0) && !bSelected ) ||
1186 ((i==gr1) && bTwo ) ||
1187 ((i==grI) && bInsert ) ) {
1188 if (bContact) {
1189 hb->a[i] = search_acceptors(&top,isize[i],index[i],bNitAcc,
1190 TRUE,(i==gr0));
1191 hb->d[i] = search_donors (&top,isize[i],index[i],
1192 TRUE,(i==gr1));
1193 }
1194 else {
1195 hb->a[i] = search_acceptors(&top,isize[i],index[i],bNitAcc,FALSE,TRUE);
1196 hb->d[i] = search_donors (&top,isize[i],index[i],FALSE,TRUE);
1197 }
1198 fprintf(stderr,"Found %d donors and %d acceptors in group '%s'\n",
1199 hb->d[i]->nrd,hb->a[i]->nra,grpnames[i]);
1200 }
1201 /*if (bSelected)
1202 snew(donors[gr0D], dons[gr0D].nrd);*/
1203 if (!bTwo) {
1204 hb->d[gr1] = hb->d[gr0];
1205 hb->a[gr1] = hb->a[gr0];
1206 /*donors[gr1D] = donors[gr0D];*/
1207 }
1208 if (!bInsert) {
1209 /* Initiate structures with zero entries */
1210 snew(hb->d[grI],1);
1211 snew(hb->a[grI],1);
1212 }
1213 /* Generate hbond data structure */
1214 mk_hbmap(hb,bTwo,bInsert);
1215
1216 /* check input */
1217 bStop=FALSE;
1218 for (grp=gr0; (grp <= (bTwo ? gr1 : gr0)); grp++)
1219 if (hb->d[grp]->nrd + hb->a[grp]->nra == 0) {
1220 fprintf(stderr,"No Donors or Acceptors found in group '%s'\n",
1221 grpnames[grp]);
1222 bStop=TRUE;
1223 }
1224 if (!bStop) {
1225 if (hb->d[gr0]->nrd + hb->d[gr1]->nrd == 0) {
1226 fprintf(stderr,"No Donors found\n");
1227 bStop=TRUE;
1228 }
1229 if (hb->a[gr0]->nra + hb->a[gr1]->nra == 0) {
1230 fprintf(stderr,"No Acceptors found\n");
1231 bStop=TRUE;
1232 }
1233 }
1234 if (bStop)
1235 fatal_error(0,"Nothing to be done");
1236 if ( bInsert && ((hb->d[grI]->nrd == 0) || (hb->a[grI]->nra == 0)) )
1237 fatal_error(0,"No %s%s%s found in insertion group '%s'\n",
1238 (hb->d[grI]->nrd == 0) ? "donors" : "",
1239 ((hb->d[grI]->nrd == 0) &&
1240 (hb->a[grI]->nra == 0)) ? " and " : "",
1241 (hb->a[grI]->nra == 0) ? "acceptors" : "",
1242 grpnames[grI]);
1243
1244 shatom=0;
1245 if (rshell > 0) {
1246 int shisz;
1247 atom_id *shidx;
1248 char *shgrpnm;
1249 /* get index group with atom for shell */
1250 do {
1251 fprintf(stderr,"Select atom for shell (1 atom):\n");
1252 get_index(&(top.atoms),ftp2fn_null(efNDX,NFILE,fnm),
1253 1,&shisz,&shidx,&shgrpnm);
1254 if (shisz != 1)
1255 fprintf(stderr,"group contains %d atoms, should be 1 (one)\n",shisz);
1256 } while(shisz != 1);
1257 shatom = shidx[0];
1258 fprintf(stderr,"Will calculate hydrogen bonds within a shell "
1259 "of %g nm around atom %i\n",rshell,shatom+1);
1260 }
1261
1262 /* Analyze trajectory */
1263 natoms=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
1264 if ( natoms > top.atoms.nr )
1265 fatal_error(0,"Topology (%d atoms) does not match trajectory (%d atoms)",
1266 top.atoms.nr,natoms);
1267
1268 bBox = ir.ePBC!=epbcNONE;
1269 grid = init_grid(bBox, box, rcut, ngrid);
1270 nabin = acut/abin;
1271 nrbin = rcut/rbin;
1272 snew(adist,nabin);
1273 snew(rdist,nrbin);
1274
1275 do {
1276 bTric = bBox && TRICLINIC(box);
1277 build_grid(hb,x,x[shatom], bBox,box,hbox, rcut, rshell, ngrid,grid);
1278 if (debug)
1279 dump_grid(debug, ngrid, grid);
1280
1281 add_frames(hb,nframes);
1282 init_hbframe(hb,nframes,t);
1283
1284 if (hb->bDAnr)
1285 count_da_grid(ngrid, grid, hb->danr[nframes]);
1286 /* loop over all gridcells (xi,yi,zi) */
1287 /* Removed confusing macro, DvdS 27/12/98 */
1288 for(xi=0; (xi<ngrid[XX]); xi++)
1289 for(yi=0; (yi<ngrid[YY]); yi++)
1290 for(zi=0; (zi<ngrid[ZZ]); zi++) {
1291
1292 /* loop over donor groups gr0 (always) and gr1 (if necessary) */
1293 for (grp=gr0; (grp <= (bTwo?gr1:gr0)); grp++) {
1294 icell=&(grid[zi][yi][xi].d[grp]);
1295
1296 /* loop over all hydrogen atoms from group (grp)
1297 * in this gridcell (icell)
1298 */
1299 for (ai=0; (ai<icell->nr); ai++) {
1300 i = icell->atoms[ai];
1301
1302 /* loop over all adjacent gridcells (xj,yj,zj) */
1303 /* This is a macro!!! */
1304 LOOPGRIDINNER(xj,yj,zj,xjj,yjj,zjj,xi,yi,zi,ngrid,bTric) {
1305 jcell=&(grid[zj][yj][xj].a[OGRP]);
1306 /* loop over acceptor atoms from other group (OGRP)
1307 * in this adjacent gridcell (jcell)
1308 */
1309 for (aj=0; (aj<jcell->nr); aj++) {
1310 j = jcell->atoms[aj];
1311
1312 if ((bSelected && (j == i)) || (!bSelected)) {
1313 /* check if this once was a h-bond */
1314 if (is_hbond(hb,grp,OGRP,i,j,rcut,ccut,x,
1315 bBox,box,hbox,&dist,&ang,bDA,&h,bContact)) {
1316 /* add to index if not already there */
1317 /* Add a hbond */
1318 add_hbond(hb,i,j,h,grp,OGRP,nframes,FALSE,bMerge);
1319
1320 /* make angle and distance distributions */
1321 ang*=RAD2DEG;
1322 adist[(int)( ang/abin)]++;
1323 rdist[(int)(dist/rbin)]++;
1324
1325 if (!bTwo) {
1326 int id = donor_index(hb->d[grp],i);
1327 int ia = acceptor_index(hb->a[OGRP],j);
1328 resdist=abs(top.atoms.atom[id].resnr-
1329 top.atoms.atom[ia].resnr);
1330 if (resdist >= max_hx)
1331 resdist = max_hx-1;
1332 hb->nhx[nframes][resdist]++;
1333 }
1334 }
1335 if (bInsert && bSelected) {
1336 /* this has been a h-bond, or we are analyzing
1337 selected bonds: check for inserted */
1338 bool ins_d, ins_a;
1339 real ins_d_dist, ins_d_ang, ins_a_dist, ins_a_ang;
1340 int ins_d_k=0,ins_a_k=0;
1341
1342 ins_d=ins_a=FALSE;
1343 ins_d_dist=ins_d_ang=ins_a_dist=ins_a_ang=1e6;
1344
1345 /* loop over gridcells adjacent to i (xk,yk,zk) */
1346 LOOPGRIDINNER(xk,yk,zk,xkk,ykk,zkk,xi,yi,zi,ngrid,bTric){
1347 kcell=&(grid[zk][yk][xk].a[grI]);
1348 /* loop over acceptor atoms from ins group
1349 in this adjacent gridcell (kcell) */
1350 for (ak=0; (ak<kcell->nr); ak++) {
1351 k=kcell->atoms[ak];
1352 if (is_hbond(hb,grp,grI,i,k,rcut,ccut,x,
1353 bBox,box,hbox,&dist,&ang,bDA,&h,
1354 bContact)) {
1355 if (dist < ins_d_dist) {
1356 ins_d=TRUE;
1357 ins_d_dist=dist;
1358 ins_d_ang =ang ;
1359 ins_d_k =k ;
1360 }
1361 }
1362 }
1363 }
1364 ENDLOOPGRIDINNER;
1365 /* loop over gridcells adjacent to j (xk,yk,zk) */
1366 LOOPGRIDINNER(xk,yk,zk,xkk,ykk,zkk,xj,yj,zj,ngrid,bTric){
1367 kcell=&grid[zk][yk][xk].d[grI];
1368 /* loop over hydrogen atoms from ins group
1369 in this adjacent gridcell (kcell) */
1370 for (ak=0; ak<kcell->nr; ak++) {
1371 k=kcell->atoms[ak];
1372 if (is_hbond(hb,grI,OGRP,k,j,rcut,ccut,x,
1373 bBox,box,hbox,&dist,&ang,bDA,&h,
1374 bContact) ) {
1375 if (dist<ins_a_dist) {
1376 ins_a=TRUE;
1377 ins_a_dist=dist;
1378 ins_a_ang =ang ;
1379 ins_a_k =k ;
1380 }
1381 }
1382 }
1383 }
1384 ENDLOOPGRIDINNER;
1385
1386 {
1387 if (ins_d && ins_a &&
1388 is_hbond(hb,grI,grI,ins_d_k,ins_a_k,rcut,ccut,x,
1389 bBox,box,hbox,&dist,&ang,bDA,&h,
1390 bContact)) {
1391 /* add to hbond index if not already there */
1392 add_hbond(hb,ins_d_k,ins_a_k,h,grI,OGRP,
1393 nframes,TRUE,bMerge);
1394
1395 /* print insertion info to file */
1396 /*fprintf(fpins,
1397 "%4g: %4u:%3.3s%4d%3.3s -> "
1398 "%4u:%3.3s%4d%3.3s (%4.2f,%2.0f) - "
1399 "%4u:%3.3s%4d%3.3s (%4.2f,%2.0f)\n",t,
1400 a[grIA][ins_d_k]+1,
1401 *top.atoms.resname[top.atoms.atom[a[grIA][ins_d_k]].resnr],
1402 top.atoms.atom[a[grIA][ins_d_k]].resnr+1,
1403 *top.atoms.atomname[a[grIA][ins_d_k]],
1404 a[grp+grD][i]+1,
1405 *top.atoms.resname[top.atoms.atom[a[grp+grD][i]].resnr],
1406 top.atoms.atom[a[grp+grD][i]].resnr+1,
1407 *top.atoms.atomname[a[grp+grD][i]],
1408 ins_d_dist,ins_d_ang*RAD2DEG,
1409 a[OGRP+grA][j]+1,
1410 *top.atoms.resname[top.atoms.atom[a[OGRP+grA][j]].resnr],
1411 top.atoms.atom[a[OGRP+grA][j]].resnr+1,
1412 *top.atoms.atomname[a[OGRP+grA][j]],
1413 ins_a_dist,ins_a_ang*RAD2DEG);*/
1414 }
1415 }
1416 }
1417 }
1418 } /* for aj */
1419 }
1420 ENDLOOPGRIDINNER;
1421 } /* for ai */
1422 } /* for grp */
1423 } /* for xi,yi,zi */
1424 nframes++;
1425 } while (read_next_x(status,&t,natoms,x,box));
1426
1427 free_grid(ngrid,&grid);
1428
1429 close_trj(status);
1430 if (bInsert)
1431 fclose(fpins);
1432
1433 if (hb->nrhb==0)
1434 fprintf(stderr,"No hydrogen bonds found!!\n");
1435 else {
1436 fprintf(stderr,"Found %d different hydrogen bonds in trajectory\n",
1437 hb->nrhb);
1438
1439 if (bMerge)
1440 merge_hb(hb,bTwo);
1441
1442 aver_nhb = 0;
1443 fp = xvgropen(opt2fn("-num",NFILE,fnm),bContact ? "Contacts" :
1444 "Hydrogen Bonds","Time","Number");
1445 for(i=0; (i<nframes); i++) {
1446 fprintf(fp,"%10g %10d\n",hb->time[i],hb->nhb[i]);
1447 aver_nhb += hb->nhb[i];
1448 }
1449 fclose(fp);
1450 aver_nhb /= nframes;
1451
1452 if (opt2bSet("-dist",NFILE,fnm)) {
1453 long sum;
1454
1455 sum=0;
1456 for(i=0; i<nrbin; i++)
1457 sum+=rdist[i];
1458
1459 fp = xvgropen(opt2fn("-dist",NFILE,fnm),
1460 "Hydrogen Bond Distribution",
1461 "Hydrogen - Acceptor Distance (nm)","");
1462 for(i=0; i<nrbin; i++)
1463 fprintf(fp,"%10g %10g\n",(i+0.5)*rbin,rdist[i]/(rbin*(real)sum));
1464 fclose(fp);
1465 }
1466
1467 if (opt2bSet("-ang",NFILE,fnm)) {
1468 long sum;
1469
1470 sum=0;
1471 for(i=0; i<nabin; i++)
1472 sum+=adist[i];
1473
1474 fp = xvgropen(opt2fn("-ang",NFILE,fnm),
1475 "Hydrogen Bond Distribution",
1476 "Donor - Hydrogen - Acceptor Angle (\\SO\\N)","");
1477 for(i=0; i<nabin; i++)
1478 fprintf(fp,"%10g %10g\n",(i+0.5)*abin,adist[i]/(abin*(real)sum));
1479 fclose(fp);
1480 }
1481
1482 if (opt2bSet("-hx",NFILE,fnm)) {
1483 fp = xvgropen(opt2fn("-hx",NFILE,fnm),
1484 "Hydrogen Bonds","Time(ps)","Count");
1485 xvgr_legend(fp,NRHXTYPES,hxtypenames);
1486 for(i=0; i<nframes; i++) {
1487 fprintf(fp,"%10g",hb->time[i]);
1488 for (j=0; j<max_hx; j++)
1489 fprintf(fp," %6d",hb->nhx[i][j]);
1490 fprintf(fp,"\n");
1491 }
1492 fclose(fp);
1493 }
1494
1495 /* Compute maximum possible number of different hbonds */
1496 if (maxnhb > 0)
1497 max_nhb = maxnhb;
1498 else {
1499 if (bTwo)
1500 max_nhb = 0.5*(hb->d[gr0]->nrd*hb->a[gr1]->nra +
1501 hb->d[gr1]->nrd*hb->a[gr0]->nra);
1502 else
1503 max_nhb = 0.5*(hb->d[gr0]->nrd*hb->a[gr0]->nra);
1504 }
1505 printf("Average number of hbonds per timeframe %.3f out of %g possible\n",
1506 aver_nhb,max_nhb);
1507 if (opt2bSet("-ac",NFILE,fnm))
1508 do_hbac(opt2fn("-ac",NFILE,fnm),hb,aver_nhb/max_nhb,bDump,bMerge);
1509
1510 if (opt2bSet("-hbm",NFILE,fnm)) {
1511 t_matrix mat;
1512 int id,ia,hh,x,y;
1513
1514 mat.nx=nframes;
1515 mat.ny=hb->nrhb;
1516 snew(mat.matrix,mat.nx);
1517 for(x=0; (x<mat.nx); x++){
1518 snew(mat.matrix[x],mat.ny);
1519 y=0;
1520 for(id=0; (id<hb->nd); id++)
1521 for(ia=0; (ia<hb->na); ia++)
1522 for(hh=0; (hh<hb->maxhydro); hh++)
1523 if (hb->hbmap[id][ia].exist[hh])
1524 for(j=0; (j<mat.ny); j++)
1525 mat.matrix[x][y++]=is_hb(hb->hbmap[id][ia].exist[hh],x);
1526
1527 mat.axis_x=hb->time;
1528 snew(mat.axis_y,mat.ny);
1529 for(j=0; j<mat.ny; j++)
1530 mat.axis_y[j]=j;
1531 sprintf(mat.title,"Hydrogen Bond Existence Map");
1532 sprintf(mat.legend,"Hydrogen Bonds");
1533 sprintf(mat.label_x,"time (ps)");
1534 sprintf(mat.label_y,"Hydrogen Bond Index");
1535 mat.bDiscrete=TRUE;
1536 if (bInsert)
1537 mat.nmap=HB_NR;
1538 else
1539 mat.nmap=2;
1540 snew(mat.map,mat.nmap);
1541 for(i=0; i<mat.nmap; i++) {
1542 mat.map[i].code.c1=hbmap[i];
1543 mat.map[i].desc=hbdesc[i];
1544 mat.map[i].rgb=hbrgb[i];
1545 }
1546 fp = opt2FILE("-hbm",NFILE,fnm,"w");
1547 write_xpm_m(fp, mat);
1548 fclose(fp);
1549 for(x=0; x<mat.nx; x++)
1550 sfree(mat.matrix[x]);
1551 sfree(mat.axis_y);
1552 sfree(mat.matrix);
1553 sfree(mat.map);
1554 }
1555 }
1556
1557 if (opt2bSet("-hbn",NFILE,fnm)) {
1558 fp = opt2FILE("-hbn",NFILE,fnm,"w");
1559 if (opt2bSet("-g",NFILE,fnm)) {
1560 fplog = ffopen(opt2fn("-g",NFILE,fnm),"w");
1561 fprintf(fplog,"# %10s %12s %12s\n","Donor","Hydrogen","Acceptor");
1562 }
1563 else
1564 fplog = NULL;
1565 for (grp=gr0; grp<=(bTwo?gr1:gr0); grp++) {
1566 fprintf(fp,"[ %s ]",grpnames[grp]);
1567 for (i=0; i<isize[grp]; i++) {
1568 fprintf(fp,(i%15)?" ":"\n");
1569 fprintf(fp,"%4u",index[grp][i]+1);
1570 }
1571 fprintf(fp,"\n");
1572 fprintf(fp,"[ donors_hydrogens_%s ]",grpnames[grp]);
1573 for (i=0; (i<hb->d[grp]->nrd); i++) {
1574 for(j=0; (j<hb->d[grp]->nhydro[i]); j++)
1575 fprintf(fp,"%4u %4u",hb->d[grp]->d[i]+1,
1576 hb->d[grp]->hydro[i][j]+1);
1577 fprintf(fp,"\n");
1578 }
1579 fprintf(fp,"[ acceptors_%s ]",grpnames[grp]);
1580 for (i=0; (i<hb->a[grp]->nra); i++) {
1581 fprintf(fp,(i%15)?" ":"\n");
1582 fprintf(fp,"%4u",hb->a[grp]->acc[i]+1);
1583 }
1584 fprintf(fp,"\n");
1585 }
1586 if (bTwo)
1587 fprintf(fp,"[ hbonds_%s-%s ]\n",grpnames[0],grpnames[1]);
1588 else
1589 fprintf(fp,"[ hbonds_%s ]\n",grpnames[0]);
1590 for(i=0; (i<hb->nrhb); i++) {
1591 int ddd,hhh,aaa;
1592 char ds[32],hs[32],as[32];
1593
1594 /* ddd = hb->hb[i].d;
1595 hhh = hb->hb[i].h;
1596 aaa = hb->hb[i].a;
1597 sprintf(ds,"%s",mkatomname(&top.atoms,ddd));
1598 sprintf(hs,"%s",mkatomname(&top.atoms,hhh));
1599 sprintf(as,"%s",mkatomname(&top.atoms,aaa));
1600 fprintf(fp,"%6u %6u %6u\n",ddd+1,hhh+1,aaa+1);
1601 if (fplog)
1602 fprintf(fplog,"%12s %12s %12s\n",ds,hs,as);
1603 */
1604 }
1605 if (bInsert) {
1606 if (bTwo)
1607 fprintf(fp,"[ insert_%s->%s-%s ]",
1608 grpnames[2],grpnames[0],grpnames[1]);
1609 else
1610 fprintf(fp,"[ insert_%s->%s ]",grpnames[2],grpnames[0]);
1611 j=0;
1612 /* for(i=0; (i<hb->nrhb); i++) {
1613 if (hb->hb[i].bInsert) {
1614 fprintf(fp,(j%15)?" ":"\n");
1615 fprintf(fp,"%4d",i+1);
1616 j++;
1617 }
1618 }*/
1619 fprintf(fp,"\n");
1620 }
1621 fclose(fp);
1622 if (fplog)
1623 fclose(fplog);
1624 }
1625
1626 if (hb->bDAnr) {
1627 int i,j,nleg;
1628 char **legnames;
1629 char buf[STRLEN];
1630
1631 #define USE_THIS_GROUP(j) ( (j == gr0) || (bTwo && (j == gr1)) || (bInsert && (j == grI)) )
1632
1633 fp = xvgropen(opt2fn("-da",NFILE,fnm),
1634 "Donors and Acceptors","Time(ps)","Count");
1635 nleg = (bTwo?2:1)*2 + (bInsert?2:0);
1636 snew(legnames,nleg);
1637 i=0;
1638 for(j=0; j<grNR; j++)
1639 if ( USE_THIS_GROUP(j) ) {
1640 sprintf(buf,"Donors %s",grpnames[j]);
1641 legnames[i++] = strdup(buf);
1642 sprintf(buf,"Acceptors %s",grpnames[j]);
1643 legnames[i++] = strdup(buf);
1644 }
1645 assert(i==nleg);
1646 xvgr_legend(fp,nleg,legnames);
1647 for(i=0; i<nframes; i++) {
1648 fprintf(fp,"%10g",hb->time[i]);
1649 for (j=0; (j<grNR); j++)
1650 if ( USE_THIS_GROUP(j) )
1651 fprintf(fp," %6d",hb->danr[i][j]);
1652 fprintf(fp,"\n");
1653 }
1654 fclose(fp);
1655 }
1656 }
1657
1658 thanx(stderr);
1659
1660 return 0;
1661 }
1662
The ultimate molecular dynamics simulation package