?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
changed reading hint
[gromacs/adressmacs.git] / src / tools / g_hbond.c
blobd1f5ecae7f070ad9d0699b4bc156127d9cbcc847
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 2.0
11 *
12 * Copyright (c) 1991-1999
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
15 *
16 * Please refer to:
17 * GROMACS: A message-passing parallel molecular dynamics implementation
18 * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
19 * Comp. Phys. Comm. 91, 43-56 (1995)
20 *
21 * Also check out our WWW page:
22 * http://md.chem.rug.nl/~gmx
23 * or e-mail to:
24 * gromacs@chem.rug.nl
25 *
26 * And Hey:
27 * GRowing Old MAkes el Chrono Sweat
28 */
29 static char *SRCID_g_hbond_c = "$Id$";
30
31 #ifdef HAVE_IDENT
32 #ident "@(#) g_hbond.cc 1.29 9/30/97"
33 #endif /* HAVE_IDENT */
34
35 #include "statutil.h"
36 #include "copyrite.h"
37 #include "sysstuff.h"
38 #include "futil.h"
39 #include "tpxio.h"
40 #include "physics.h"
41 #include "macros.h"
42 #include "fatal.h"
43 #include "rdgroup.h"
44 #include "smalloc.h"
45 #include "assert.h"
46 #include "vec.h"
47 #include "xvgr.h"
48 #include "gstat.h"
49 #include "matio.h"
50 #include "lutab.h"
51 #include "string2.h"
52 #include <math.h>
53
54 #define max_hx 7
55 typedef int t_hx[max_hx];
56 #define NRHXTYPES max_hx
57 char *hxtypenames[NRHXTYPES]=
58 {"n-n","n-n+1","n-n+2","n-n+3","n-n+4","n-n+5","n-n>6"};
59
60 enum { gr0D, gr0H, gr0A, gr1D, gr1H, gr1A, grID, grIH, grIA, grNR };
61 #define grD gr0D
62 #define grH gr0H
63 #define grA gr0A
64 #define gr0 gr0D
65 #define gr1 gr1D
66 #define grI grID
67 #define grINC (gr1-gr0)
68
69 typedef struct {
70 int nr;
71 int maxnr;
72 atom_id *atoms;
73 } t_ncell;
74
75 typedef t_ncell t_gridcell[grNR];
76 typedef int t_icell[grNR];
77
78 typedef struct {
79 int a;
80 int nr;
81 } t_acceptor;
82
83 typedef struct {
84 int nrhb;
85 int maxnr;
86 t_acceptor *hb;
87 } t_donor;
88
89 bool in_list(atom_id selection,int isize,atom_id *index)
90 {
91 int i;
92 bool bFound;
93
94 bFound=FALSE;
95 for(i=0; (i<isize) && !bFound; i++)
96 if(selection == index[i])
97 bFound=TRUE;
98
99 return bFound;
100 }
101
102 void add_acc(int i, int *max_nr,int *nr,atom_id **a)
103 {
104 (*nr)++;
105 if ( *nr > *max_nr ) {
106 (*max_nr)+=10;
107 srenew(*a,*max_nr);
108 }
109 (*a)[*nr-1]=i;
110 }
111
112 void search_acceptors(t_topology *top, int isize, atom_id *index,
113 int *nr_a, atom_id **a, bool bNitAcc)
114 {
115 int i,max_nr_a;
116
117 max_nr_a=*nr_a;
118 for (i=0; i<top->atoms.nr; i++)
119 if ( ( *top->atoms.atomname[i][0] == 'O' ||
120 ( bNitAcc && ( *top->atoms.atomname[i][0] == 'N' ) ) ) &&
121 in_list(i,isize,index) )
122 add_acc(i,&max_nr_a,nr_a,a);
123 srenew(*a,*nr_a);
124 }
125
126 void add_dh(int id, int ih, int *max_nr,int *nr,atom_id **d,atom_id **h)
127 {
128 (*nr)++;
129 if ( *nr > *max_nr ) {
130 (*max_nr)+=10;
131 srenew(*d,*max_nr);
132 srenew(*h,*max_nr);
133 }
134 (*d)[*nr-1]=id;
135 (*h)[*nr-1]=ih;
136 }
137
138 void search_donors(t_topology *top, int isize, atom_id *index,
139 int *nr_d, atom_id **d, atom_id **h)
140 {
141 int i,j,max_nr_d;
142 t_functype func_type;
143 t_ilist *interaction;
144 atom_id nr1,nr2;
145 bool stop;
146
147
148 max_nr_d=*nr_d;
149 for(func_type=0; func_type < F_NRE; func_type++) {
150 interaction=&top->idef.il[func_type];
151 for(i=0; i < interaction->nr; i+=interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1 /* next function */) {
152 assert(func_type == top->idef.functype[interaction->iatoms[i]]);
153
154 /* check out this functype */
155 if (func_type == F_SETTLE) {
156 nr1=interaction->iatoms[i+1];
157
158 if (in_list(nr1, isize,index) &&
159 in_list(nr1+1,isize,index) &&
160 in_list(nr1+2,isize,index) ) {
161 add_dh(nr1,nr1+1,&max_nr_d,nr_d,d,h);
162 add_dh(nr1,nr1+2,&max_nr_d,nr_d,d,h);
163 }
164 } else if ( interaction_function[func_type].flags & IF_CONNECT ) {
165 for (j=0; j<2; j++) {
166 nr1=interaction->iatoms[i+1+j];
167 nr2=interaction->iatoms[i+2-j];
168 if ( ( *top->atoms.atomname[nr1][0] == 'H' ) &&
169 ( ( *top->atoms.atomname[nr2][0] == 'O' ) ||
170 ( *top->atoms.atomname[nr2][0] == 'N' )) &&
171 in_list(nr1,isize,index) && in_list(nr2,isize,index))
172 add_dh(nr2,nr1,&max_nr_d,nr_d,d,h);
173 }
174 } else if ( interaction_function[func_type].flags & IF_DUMMY ) {
175 nr1=interaction->iatoms[i+1];
176 if ( *top->atoms.atomname[nr1][0] == 'H') {
177 nr2=nr1-1;
178 stop=FALSE;
179 while (!stop && ( *top->atoms.atomname[nr2][0] == 'H'))
180 if (nr2)
181 nr2--;
182 else
183 stop=TRUE;
184 if ( !stop && ( ( *top->atoms.atomname[nr2][0] == 'O') ||
185 ( *top->atoms.atomname[nr2][0] == 'N') ) &&
186 in_list(nr1,isize,index) && in_list(nr2,isize,index) )
187 add_dh(nr2,nr1,&max_nr_d,nr_d,d,h);
188 }
189 }
190 }
191 }
192 srenew(*d,*nr_d);
193 srenew(*h,*nr_d);
194 }
195
196 void init_grid(bool bBox, matrix box, real rcut,
197 ivec ngrid, t_gridcell ****grid)
198 {
199 int i,x,y,z;
200
201 if (bBox)
202 for(i=0; i<DIM; i++)
203 ngrid[i]=box[i][i]/(1.2*rcut);
204
205 if ( !bBox || (ngrid[XX]<3) || (ngrid[YY]<3) || (ngrid[ZZ]<3) )
206 for(i=0; i<DIM; i++)
207 ngrid[i]=1;
208 if (debug)
209 fprintf(debug,"Will do grid-seach on %dx%dx%d grid\n",
210 ngrid[XX],ngrid[YY],ngrid[ZZ]);
211 snew(*grid,ngrid[XX]);
212 for (x=0; x<ngrid[XX]; x++) {
213 snew((*grid)[x],ngrid[YY]);
214 for (y=0; y<ngrid[YY]; y++)
215 snew((*grid)[x][y],ngrid[ZZ]);
216 }
217 }
218
219 char *grpnames[grNR] = {"0D","0H","0A","1D","1H","1A","iD","iH","iA"};
220
221 void build_grid(int *nr, atom_id **a, rvec x[], rvec xshell,
222 bool bBox, matrix box, rvec hbox, real rcut, real rshell,
223 ivec ngrid, t_gridcell ***grid)
224 {
225 int i,m,gr,xi,yi,zi;
226 ivec grididx;
227 rvec invdelta,dshell;
228 t_ncell *newgrid;
229 bool bDoRshell,bInShell;
230 real rshell2;
231
232 bDoRshell = rshell > 0;
233 if (bDoRshell)
234 rshell2 = sqr(rshell);
235 bInShell = TRUE;
236
237 for(m=0; m<DIM; m++) {
238 hbox[m]=box[m][m]*0.5;
239 if (bBox) {
240 invdelta[m]=ngrid[m]/box[m][m];
241 if (1/invdelta[m] < rcut)
242 fatal_error(0,"box shrank too much to keep using this grid\n");
243 } else
244 invdelta[m]=0;
245 }
246 for(gr=0; gr<grNR; gr++) {
247 /* resetting atom counts */
248 for (xi=0; xi<ngrid[XX]; xi++)
249 for (yi=0; yi<ngrid[YY]; yi++)
250 for (zi=0; zi<ngrid[ZZ]; zi++)
251 grid[xi][yi][zi][gr].nr=0;
252 if (debug)
253 fprintf(debug,"Putting %d %s atoms into grid\n",nr[gr],grpnames[gr]);
254 /* put atoms in grid cells */
255 for(i=0; i<nr[gr]; i++) {
256 /* check if we are inside the shell */
257 /* if bDoRshell=FALSE then bInShell=TRUE always */
258 if ( bDoRshell ) {
259 bInShell=TRUE;
260 for(m=0; (m<DIM) && bInShell; m++) {
261 dshell[m] = x[a[gr][i]][m] - xshell[m];
262 if (bBox)
263 if ( dshell[m] < -hbox[m] )
264 dshell[m] += 2*hbox[m];
265 else if ( dshell[m] >= hbox[m] )
266 dshell[m] -= 2*hbox[m];
267 /* if we're outside the cube, we're outside the sphere also! */
268 if ( (dshell[m]>rshell) || (-dshell[m]>rshell) )
269 bInShell=FALSE;
270 }
271 /* if we're inside the cube, check if we're inside the sphere */
272 if (bInShell)
273 bInShell = norm2(dshell) < rshell2;
274 }
275 if ( bInShell ) {
276 for(m=0; m<DIM; m++) {
277 /* put atom in the box */
278 if (bBox) {
279 while( x[a[gr][i]][m] < 0 )
280 x[a[gr][i]][m] += box[m][m];
281 while( x[a[gr][i]][m] >= box[m][m] )
282 x[a[gr][i]][m] -= box[m][m];
283 }
284 /* get grid index of atom */
285 grididx[m]=x[a[gr][i]][m]*invdelta[m];
286 }
287 /* add atom to grid cell */
288 newgrid=&(grid[grididx[XX]][grididx[YY]][grididx[ZZ]][gr]);
289 if (newgrid->nr >= newgrid->maxnr) {
290 newgrid->maxnr+=10;
291 srenew(newgrid->atoms, newgrid->maxnr);
292 }
293 newgrid->atoms[newgrid->nr]=i;
294 newgrid->nr++;
295 }
296 }
297 }
298 }
299
300 void count_da_grid(ivec ngrid, t_gridcell ***grid, t_icell danr)
301 {
302 int gr,xi,yi,zi;
303
304 for(gr=0; gr<grNR; gr++) {
305 danr[gr]=0;
306 for (xi=0; xi<ngrid[XX]; xi++)
307 for (yi=0; yi<ngrid[YY]; yi++)
308 for (zi=0; zi<ngrid[ZZ]; zi++)
309 danr[gr] += grid[xi][yi][zi][gr].nr;
310 }
311 }
312
313 #define B(n,x)((n!=1)?(x-1):x)
314 #define E(n,x)((n!=1)?(x+1):x)
315 #define GRIDMOD(j,n) (j+n)%(n)
316 #define LOOPGRIDINNER(x,y,z,xx,yy,zz,xo,yo,zo,n)\
317 for(xx=B(n[XX],xo); xx<=E(n[XX],xo); xx++) {\
318 x=GRIDMOD(xx,n[XX]);\
319 for(yy=B(n[YY],yo); yy<=E(n[YY],yo); yy++) {\
320 y=GRIDMOD(yy,n[YY]);\
321 for(zz=B(n[ZZ],zo); zz<=E(n[ZZ],zo); zz++) {\
322 z=GRIDMOD(zz,n[ZZ]);
323 #define ENDLOOPGRIDINNER\
324 }\
325 }\
326 }\
327
328 void dump_grid(FILE *fp, ivec ngrid, t_gridcell ***grid)
329 {
330 int gr,x,y,z,sum[grNR];
331
332 fprintf(fp,"grid %dx%dx%d\n",ngrid[XX],ngrid[YY],ngrid[ZZ]);
333 for (gr=0; gr<grNR; gr++) {
334 sum[gr]=0;
335 fprintf(fp,"GROUP %d (%s)\n",gr,grpnames[gr]);
336 for (z=0; z<ngrid[ZZ]; z+=2) {
337 fprintf(fp,"Z=%d,%d\n",z,z+1);
338 for (y=0; y<ngrid[YY]; y++) {
339 for (x=0; x<ngrid[XX]; x++) {
340 fprintf(fp,"%3d",grid[x][y][z][gr].nr);
341 sum[gr]+=grid[x][y][z][gr].nr;
342 }
343 fprintf(fp," | ");
344 if ( (z+1) < ngrid[ZZ] )
345 for (x=0; x<ngrid[XX]; x++) {
346 fprintf(fp,"%3d",grid[x][y][z+1][gr].nr);
347 sum[gr]+=grid[x][y][z+1][gr].nr;
348 }
349 fprintf(fp,"\n");
350 }
351 }
352 }
353 fprintf(fp,"TOTALS:");
354 for (gr=0; gr<grNR; gr++)
355 fprintf(fp," %d=%d",gr,sum[gr]);
356 fprintf(fp,"\n");
357 }
358
359 /* New GMX record! 5 * in a row. Congratulations!
360 * Sorry, only four left.
361 */
362 void free_grid(ivec ngrid, t_gridcell ****grid)
363 {
364 int x,y,z,i;
365 t_gridcell ***g = *grid;
366
367 for (x=0; x<ngrid[XX]; x++) {
368 for (y=0; y<ngrid[YY]; y++) {
369 sfree(g[x][y]);
370 }
371 sfree(g[x]);
372 }
373 sfree(g);
374 g=NULL;
375 }
376
377 bool is_hbond(atom_id d, atom_id h, atom_id a,
378 real rcut, real ccut,
379 rvec x[], bool bBox, rvec hbox, real *d_ha, real *ang)
380 {
381 rvec r_dh,r_ha;
382 real d_ha2,ca;
383 int m;
384
385 for(m=0; m<DIM; m++) {
386 r_ha[m] = x[a][m] - x[h][m];
387 if (bBox) {
388 if ( r_ha[m] < -hbox[m] )
389 r_ha[m] += 2*hbox[m];
390 else if ( r_ha[m] >= hbox[m] )
391 r_ha[m] -= 2*hbox[m];
392 }
393 if ( (r_ha[m]>rcut) || (-r_ha[m]>rcut) )
394 return FALSE;
395 }
396 d_ha2 = iprod(r_ha,r_ha);
397 if ( d_ha2 <= rcut*rcut ) {
398 rvec_sub(x[h],x[d],r_dh);
399 if (bBox)
400 for(m=0; m<DIM; m++) {
401 if ( r_dh[m] < -hbox[m] )
402 r_dh[m] += 2*hbox[m];
403 else if ( r_dh[m] >= hbox[m] )
404 r_dh[m] -= 2*hbox[m];
405 }
406 ca = cos_angle(r_dh,r_ha);
407 /* if angle is smaller, cos is larger */
408 if ( ca>=ccut ) {
409 *d_ha = sqrt(d_ha2);
410 *ang = acos(ca);
411 return TRUE;
412 }
413 }
414 return FALSE;
415 }
416
417 void sort_dha(int nr_d, atom_id *d, atom_id *h, int nr_a, atom_id *a)
418 {
419 int i,j;
420 atom_id temp;
421
422 for(i=0; i<nr_d; i++)
423 for(j=i+1; j<nr_d; j++)
424 if ( (d[i] > d[j]) ||
425 ( (d[i] == d[j]) && (h[i] > h[j]) ) ) {
426 temp=d[i];
427 d[i]=d[j];
428 d[j]=temp;
429 temp=h[i];
430 h[i]=h[j];
431 h[j]=temp;
432 }
433 for(i=0; i<nr_a; i++)
434 for(j=i+1; j<nr_a; j++)
435 if (a[i] > a[j]) {
436 temp=a[i];
437 a[i]=a[j];
438 a[j]=temp;
439 }
440 }
441
442 void sort_hb(int *nr_a, t_donor **donors)
443 {
444 int gr,i,j,k;
445 t_acceptor ta;
446
447 fprintf(stderr,"Sorting hydrogen bonds\n");
448 for (gr=0; gr<grNR; gr+=grINC)
449 for (i=0; i<nr_a[gr+grD]; i++)
450 for (j=0; j<donors[gr+grD][i].nrhb; j++)
451 for (k=j+1; k<donors[gr+grD][i].nrhb; k++)
452 if (donors[gr+grD][i].hb[j].a > donors[gr+grD][i].hb[k].a) {
453 ta = donors[gr+grD][i].hb[k];
454 donors[gr+grD][i].hb[k] = donors[gr+grD][i].hb[j];
455 donors[gr+grD][i].hb[j] = ta;
456 }
457 }
458
459 int main(int argc,char *argv[])
460 {
461 static char *desc[] = {
462 "g_hbond computes and analyzes hydrogen bonds. Hydrogen bonds are",
463 "determined based on cutoffs for the angle Donor - Hydrogen - Acceptor",
464 "(zero is extended) and the distance Hydrogen - Acceptor.",
465 "OH and NH groups are regarded as donors, O is an acceptor always,",
466 "N is an acceptor by default, but this can be switched using",
467 "[TT]-nitacc[tt]. Dummy hydrogen atoms are assumed to be connected",
468 "to the first preceding non-hydrogen atom.[PAR]",
469
470 "You need to specify two groups for analysis, which must be either",
471 "identical or non-overlapping. All hydrogen bonds between the two",
472 "groups are analyzed.[PAR]",
473
474 "If you set -shell, you will be asked for an additional index group",
475 "which should contain exactly one atom. In this case, only hydrogen",
476 "bonds between atoms within the shell distance from the one atom are",
477 "considered.[PAR]"
478
479 "It is also possible to analyse specific hydrogen bonds with",
480 "[TT]-sel[tt]. This index file must contain a group of atom triplets",
481 "Donor Hydrogen Acceptor, in the following way:[PAR]",
482
483 "[TT]",
484 "[ selected ][BR]",
485 " 20 21 24[BR]",
486 " 25 26 29[BR]",
487 " 1 3 6[BR]",
488 "[tt][BR]",
489 "Note that the triplets need not be on separate lines.",
490 "Each atom triplet specifies a hydrogen bond to be analyzed,",
491 "note also that no check is made for the types of atoms.[PAR]",
492
493 "[TT]-ins[tt] turns on computing solvent insertion into hydrogen bonds.",
494 "In this case an additional group must be selected, specifying the",
495 "solvent molecules.[PAR]",
496
497 "[BB]Output:[bb][BR]",
498 "[TT]-num[tt]: number of hydrogen bonds as a function of time.[BR]",
499 "[TT]-ac[tt]: average over all autocorrelations of the existence",
500 "functions (either 0 or 1) of all hydrogen bonds.[BR]",
501 "[TT]-dist[tt]: distance distribution of all hydrogen bonds.[BR]",
502 "[TT]-ang[tt]: angle distribution of all hydrogen bonds.[BR]",
503 "[TT]-hx[tt]: the number of n-n+i hydrogen bonds as a function of time",
504 "where n and n+i stand for residue numbers and i ranges from 0 to 6.",
505 "This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated",
506 "with helices in proteins.[BR]",
507 "[TT]-hbn[tt]: all selected groups, donors, hydrogens and acceptors",
508 "for selected groups, all hydrogen bonded atoms from all groups and",
509 "all solvent atoms involved in insertion.[BR]",
510 "[TT]-hbm[tt]: existence matrix for all hydrogen bonds over all",
511 "frames, this also contains information on solvent insertion",
512 "into hydrogen bonds.[BR]",
513 "[TT}-da[tt]: write out the number of donors and acceptors analyzed for",
514 "each timeframe. This is especially usefull when using [TT]-shell[tt]."
515 };
516
517 static real acut=60, abin=1, rcut=0.25, rbin=0.005, rshell=-1;
518 static bool bNitAcc=TRUE,bInsert=FALSE;
519 /* options */
520 t_pargs pa [] = {
521 { "-ins", FALSE, etBOOL, {&bInsert},
522 "analyze solvent insertion" },
523 { "-a", FALSE, etREAL, {&acut},
524 "cutoff angle (degrees, Donor - Hydrogen - Acceptor)" },
525 { "-r", FALSE, etREAL, {&rcut},
526 "cutoff radius (nm, Hydrogen - Acceptor)" },
527 { "-abin", FALSE, etREAL, {&abin},
528 "binwidth angle distribution (degrees)" },
529 { "-rbin", FALSE, etREAL, {&rbin},
530 "binwidth distance distribution (nm)" },
531 { "-nitacc",FALSE,etBOOL, {&bNitAcc},
532 "regard nitrogen atoms as acceptors" },
533 { "-shell", FALSE,etREAL, {&rshell},
534 "only calculate hydrogen bonds within shell around one particle" }
535 };
536
537 t_filenm fnm[] = {
538 { efTRX, "-f", NULL, ffREAD },
539 { efTPX, NULL, NULL, ffREAD },
540 { efNDX, NULL, NULL, ffOPTRD },
541 { efNDX, "-sel", "select", ffOPTRD },
542 { efXVG, "-num", "hbnum", ffWRITE },
543 { efXVG, "-ac", "hbac", ffOPTWR },
544 { efXVG, "-dist","hbdist", ffOPTWR },
545 { efXVG, "-ang", "hbang", ffOPTWR },
546 { efXVG, "-hx", "hbhelix",ffOPTWR },
547 { efNDX, "-hbn", "hbond", ffOPTWR },
548 { efXPM, "-hbm", "hbmap", ffOPTWR },
549 { efXVG, "-da", "danum", ffOPTWR }
550 };
551 #define NFILE asize(fnm)
552
553 #define FRINC 100
554 #define HBINC 100
555
556 #define NRGRPS 3
557 #define OGRP (gr1-grp)
558 #define INSGRP 2
559
560 #define HB_NO 0
561 #define HB_YES 1<<0
562 #define HB_INS 1<<1
563 #define HB_YESINS HB_YES|HB_INS
564 #define HB_NR (1<<2)
565 char hbmap [HB_NR]={ ' ', 'o', '-', '*' };
566 char *hbdesc[HB_NR]={ "None", "Present", "Inserted", "Present & Inserted" };
567 t_rgb hbrgb [HB_NR]={ {1,1,1},{1,0,0}, {0,0,1}, {1,0,1} };
568
569 int status;
570 t_topology top;
571 t_inputrec ir;
572 int natoms,nframes,max_nframes,shatom;
573 int *isize;
574 char **grpnames;
575 atom_id **index;
576 rvec *x,hbox;
577 matrix box;
578 real t,ccut,dist,ang;
579 real *time;
580 int i,j,k,l,start,end;
581 int xi,yi,zi,ai;
582 int xj,yj,zj,aj,xjj,yjj,zjj;
583 int xk,yk,zk,ak,xkk,ykk,zkk;
584 bool bSelected,bStop,bTwo,bHBMap,new,was,bBox,bDAnr;
585 bool *insert=NULL;
586 int nr_a[grNR];
587 atom_id *a[grNR];
588 int grp;
589 int *nhb,*adist,*rdist;
590 t_hx *nhx;
591 int max_nrhb,nrhb,nabin,nrbin,bin,resdist,idx;
592 unsigned char **hbexist;
593 FILE *fp,*fpins=NULL;
594 t_gridcell ***grid;
595 t_ncell *icell,*jcell,*kcell;
596 t_icell *danr;
597 ivec ngrid;
598 t_donor *donors[grNR];
599 real **rhbex;
600
601 CopyRight(stderr,argv[0]);
602 parse_common_args(&argc,argv,PCA_CAN_TIME,TRUE,NFILE,fnm,asize(pa),pa,
603 asize(desc),desc,0,NULL);
604 /* Initiate lookup table for sqrt calculations */
605 init_lookup_table(stdout);
606
607 /* process input */
608 bHBMap = opt2bSet("-hbm",NFILE,fnm) || opt2bSet("-ac",NFILE,fnm);
609 bDAnr = opt2bSet("-da",NFILE,fnm);
610 bSelected = opt2bSet("-sel",NFILE,fnm);
611 ccut = cos(acut*DEG2RAD);
612
613 /* get topology */
614 read_tpx(ftp2fn(efTPX,NFILE,fnm),&i,&t,&t,
615 &ir,box,&natoms,NULL,NULL,NULL,&top);
616
617 /* initialize h-bond atom groups: */
618 for (i=gr0; i<grNR; i+=grINC) {
619 nr_a[i+grD] = nr_a[i+grH] = nr_a[i+grA] = 0;
620 a[i+grD] = a[i+grH] = a[i+grA] = NULL;
621 }
622
623 snew(grpnames,NRGRPS);
624 snew(index,NRGRPS);
625 snew(isize,NRGRPS);
626 if (bSelected) {
627 /* analyze selected hydrogen bonds */
628 fprintf(stderr,"Select group with selected atoms:\n");
629 get_index(&(top.atoms),opt2fn("-sel",NFILE,fnm),
630 1,isize,index,grpnames);
631 if (isize[0] % 3)
632 fatal_error(0,"Number of atoms in group '%s' not a multiple of 3\n"
633 "and therefore cannot contain triplets of "
634 "Donor-Hydrogen-Acceptor",grpnames[0]);
635 bTwo=FALSE;
636 for(grp=gr0; grp<gr0+grINC; grp++) {
637 nr_a[grp]=isize[0]/3;
638 snew(a[grp],nr_a[grp]);
639 for(i=0; i<nr_a[grp]; i++)
640 a[grp][i]=index[0][i*3+grp];
641 }
642 fprintf(stderr,"Analyzing %d selected hydrogen bonds from '%s'\n",
643 nr_a[gr0D],grpnames[0]);
644 } else {
645 /* analyze all hydrogen bonds: get group(s) */
646 fprintf(stderr,"Specify 2 groups to analyze:\n");
647 get_index(&(top.atoms),ftp2fn_null(efNDX,NFILE,fnm),
648 2,isize,index,grpnames);
649
650 /* check if we have two identical or two non-overlapping groups */
651 bTwo = isize[0] != isize[1];
652 for (i=0; (i<isize[0]) && !bTwo; i++)
653 bTwo = index[0][i] != index[1][i];
654 if (bTwo) {
655 fprintf(stderr,"Checking for overlap...\n");
656 for (i=0; i<isize[0]; i++)
657 for (j=0; j<isize[1]; j++)
658 if (index[0][i] == index[1][j])
659 fatal_error(0,"Partial overlap between groups '%s' and '%s'",
660 grpnames[0],grpnames[1]);
661 }
662 if (bTwo)
663 fprintf(stderr,"Calculating hydrogen bonds "
664 "between two groups of %d and %d atoms\n",
665 isize[0],isize[1]);
666 else
667 fprintf(stderr,"Calculating hydrogen bonds in one group of %d atoms\n",
668 isize[0]);
669 }
670 if (bInsert) {
671 fprintf(stderr,"Specify group for insertion analysis:\n");
672 get_index(&(top.atoms),ftp2fn_null(efNDX,NFILE,fnm),
673 1,&(isize[INSGRP]),&(index[INSGRP]),&(grpnames[INSGRP]));
674 fprintf(stderr,"Checking for overlap...\n");
675 for (i=0; i<isize[INSGRP]; i++)
676 for (grp=0; grp<(bTwo?2:1); grp++)
677 for (j=0; j<isize[grp]; j++)
678 if (index[INSGRP][i] == index[grp][j])
679 fatal_error(0,"Partial overlap between groups '%s' and '%s'",
680 grpnames[grp],grpnames[INSGRP]);
681 fpins=ffopen("insert.dat","w");
682 fprintf(fpins,"%4s: %15s -> %15s (%7s) - %15s (%7s)\n",
683 "time","insert","donor","distang","acceptor","distang");
684 }
685
686 /* search donors and acceptors in groups */
687 for (i=gr0; i<grNR; i+=grINC)
688 if ( ((i==gr0) && !bSelected ) ||
689 ((i==gr1) && bTwo ) ||
690 ((i==grI) && bInsert ) ) {
691 search_acceptors(&top, isize[i/grINC], index[i/grINC],
692 &nr_a[i+grA], &a[i+grA], bNitAcc);
693 search_donors (&top, isize[i/grINC], index[i/grINC],
694 &nr_a[i+grH], &a[i+grD], &a[i+grH]);
695 nr_a[i+grD] = nr_a[i+grH];
696 fprintf(stderr,"Found %d donors and %d acceptors in group '%s'\n",
697 nr_a[i+grD], nr_a[i+grA], grpnames[i/grINC]);
698 snew(donors[i+grD], nr_a[i+grD] );
699 sort_dha(nr_a[i+grD], a[i+grD], a[i+grH], nr_a[i+grA], a[i+grA]);
700 }
701 if (bSelected)
702 snew(donors[gr0D], nr_a[gr0D]);
703 if (!bTwo) {
704 for(i=0; i<grINC; i++) {
705 nr_a[gr1+i] = nr_a[gr0+i];
706 a[gr1+i] = a[gr0+i];
707 }
708 donors[gr1D] = donors[gr0D];
709 }
710
711 /* check input */
712 bStop=FALSE;
713 for (grp=gr0; grp<=(bTwo?gr1:gr0); grp+=grINC)
714 if (nr_a[grp+grD]+nr_a[grp+grA]==0) {
715 fprintf(stderr,"No Donors or Acceptors found in group '%s'\n",
716 grpnames[grp/grINC]);
717 bStop=TRUE;
718 }
719 if (!bStop) {
720 if (nr_a[gr0D]+nr_a[gr1D]==0) {
721 fprintf(stderr,"No Donors found\n");
722 bStop=TRUE;
723 }
724 if (nr_a[gr0A]+nr_a[gr1A]==0) {
725 fprintf(stderr,"No Acceptors found\n");
726 bStop=TRUE;
727 }
728 }
729 if (bStop)
730 fatal_error(0,"Nothing to be done");
731 if ( bInsert && ((nr_a[grID]==0) || (nr_a[grIA]==0)) )
732 fatal_error(0,"No %s%s%s found in insertion group '%s'\n",
733 (nr_a[grID]==0)?"donors":"",
734 ((nr_a[grID]==0)&&(nr_a[grIA]==0))?" and ":"",
735 (nr_a[grIA]==0)?"acceptors":"",
736 grpnames[grI/grINC]);
737
738 shatom=0;
739 if (rshell > 0) {
740 int shisz;
741 atom_id *shidx;
742 char *shgrpnm;
743 /* get index group with atom for shell */
744 do {
745 fprintf(stderr,"Select atom for shell (1 atom):\n");
746 get_index(&(top.atoms),ftp2fn_null(efNDX,NFILE,fnm),
747 1,&shisz,&shidx,&shgrpnm);
748 if (shisz != 1)
749 fprintf(stderr,"group contains %d atoms, should be 1 (one)\n",shisz);
750 } while(shisz != 1);
751 shatom = shidx[0];
752 fprintf(stderr,"Will calculate hydrogen bonds within a shell "
753 "of %g nm around atom %i\n",rshell,shatom+1);
754 }
755
756 /* analyze trajectory */
757
758 natoms=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
759 if ( natoms > top.atoms.nr )
760 fatal_error(0,"Topology (%d atoms) does not match trajectory (%d atoms)",
761 top.atoms.nr,natoms);
762 bBox=ir.eBox!=ebtNONE;
763 init_grid(bBox, box, rcut, ngrid, &grid);
764 max_nframes = nframes = 0;
765 max_nrhb = nrhb = 0;
766 hbexist = NULL;
767 nhb = NULL;
768 nhx = NULL;
769 time = NULL;
770 danr = NULL;
771 nabin = acut/abin;
772 nrbin = rcut/rbin;
773 snew(adist,nabin);
774 snew(rdist,nrbin);
775 if (bInsert)
776 snew(insert,natoms);
777 do {
778 build_grid(nr_a,a, x,x[shatom], bBox,box,hbox, rcut, rshell, ngrid,grid);
779 if (debug)
780 dump_grid(debug, ngrid, grid);
781
782 if (nframes >= max_nframes) {
783 max_nframes += FRINC;
784 srenew(nhb,max_nframes);
785 srenew(nhx,max_nframes);
786 srenew(time,max_nframes);
787 if (bHBMap)
788 for (i=0; i<max_nrhb; i++)
789 srenew(hbexist[i],max_nframes);
790 if (bDAnr)
791 srenew(danr,max_nframes);
792 }
793 time[nframes] = t;
794 nhb[nframes] = 0;
795 for (i=0; i<max_hx; i++)
796 nhx[nframes][i]=0;
797 if (bHBMap)
798 for (i=0; i<max_nrhb; i++)
799 hbexist[i][nframes]=HB_NO;
800 if (bDAnr)
801 count_da_grid(ngrid, grid, danr[nframes]);
802 /* loop over all gridcells (xi,yi,zi) */
803 /* Removed confusing macro, DvdS 27/12/98 */
804 for(xi=0; (xi<ngrid[XX]); xi++)
805 for(yi=0; (yi<ngrid[YY]); yi++)
806 for(zi=0; (zi<ngrid[ZZ]); zi++) {
807 /* loop over groups gr0 (always) and gr1 (if necessary) */
808 for (grp=gr0; grp<=(bTwo?gr1:gr0); grp+=grINC) {
809 icell=&grid[xi][yi][zi][grp+grH];
810 /* loop over all hydrogen atoms from group (grp)
811 * in this gridcell (icell)
812 */
813 for (ai=0; ai<icell->nr; ai++) {
814 i=icell->atoms[ai];
815 /* loop over all adjacent gridcells (xj,yj,zj) */
816 /* This is a macro!!! */
817 LOOPGRIDINNER(xj,yj,zj,xjj,yjj,zjj,xi,yi,zi,ngrid) {
818 jcell=&grid[xj][yj][zj][OGRP+grA];
819 /* loop over acceptor atoms from other group (OGRP)
820 * in this adjacent gridcell (jcell)
821 */
822 for (aj=0; aj<jcell->nr; aj++) {
823 j=jcell->atoms[aj];
824 if ( (bSelected && (j==i)) || (!bSelected) ) {
825 /* check if this once was a h-bond */
826 idx=NOTSET;
827 for (k=0; (k<donors[grp][i].nrhb) && (idx==NOTSET); k++)
828 if (j == donors[grp][i].hb[k].a)
829 idx=k;
830 if ( is_hbond(a[ grp+grD][i],a[ grp+grH][i],a[OGRP+grA][j],
831 rcut,ccut,x,bBox,hbox,&dist,&ang) ) {
832 /* add to index if not already there */
833 if (idx==NOTSET) {
834 if (donors[grp][i].nrhb>=donors[grp][i].maxnr) {
835 donors[grp][i].maxnr+=10;
836 srenew(donors[grp][i].hb,donors[grp][i].maxnr);
837 }
838 /* Add a donor atom in a hbond */
839 donors[grp][i].hb[donors[grp][i].nrhb].a=j;
840 donors[grp][i].hb[donors[grp][i].nrhb].nr=nrhb;
841 idx = donors[grp][i].nrhb;
842 donors[grp][i].nrhb++;
843 if (bHBMap && (nrhb >= max_nrhb) ) {
844 max_nrhb+=HBINC;
845 srenew(hbexist,max_nrhb);
846 for (l=max_nrhb-HBINC; l<max_nrhb; l++)
847 snew(hbexist[l],max_nframes);
848 }
849 nrhb++;
850 }
851 if (bHBMap)
852 /* update matrix */
853 hbexist[donors[grp][i].hb[idx].nr][nframes] |= HB_YES;
854
855 /* count number of hbonds per frame */
856 nhb[nframes]++;
857
858 /* make angle and distance distributions */
859 ang*=RAD2DEG;
860 adist[(int)( ang/abin)]++;
861 rdist[(int)(dist/rbin)]++;
862
863 if (!bTwo) {
864 resdist=abs(top.atoms.atom[a[ grp+grD][i]].resnr-
865 top.atoms.atom[a[OGRP+grA][j]].resnr);
866 if (resdist >= max_hx)
867 resdist = max_hx-1;
868 nhx[nframes][resdist]++;
869 }
870 }
871 if (bInsert && ( (idx!=NOTSET) || bSelected ) ) {
872 /* this has been a h-bond, or we are analyzing
873 selected bonds: check for inserted */
874 bool ins_d, ins_a;
875 real ins_d_dist, ins_d_ang, ins_a_dist, ins_a_ang;
876 int ins_d_k=0,ins_a_k=0;
877
878 ins_d=ins_a=FALSE;
879 ins_d_dist=ins_d_ang=ins_a_dist=ins_a_ang=1e6;
880
881 /* loop over gridcells adjacent to i (xk,yk,zk) */
882 LOOPGRIDINNER(xk,yk,zk,xkk,ykk,zkk,xi,yi,zi,ngrid) {
883 kcell=&grid[xk][yk][zk][grIA];
884 /* loop over acceptor atoms from ins group
885 in this adjacent gridcell (kcell) */
886 for (ak=0; ak<kcell->nr; ak++) {
887 k=kcell->atoms[ak];
888 if (is_hbond(a[grp+grD][i],
889 a[grp+grH][i],
890 a[ grIA][k],
891 rcut,ccut,x,bBox,hbox,&dist,&ang))
892 if (dist<ins_d_dist) {
893 ins_d=TRUE;
894 ins_d_dist=dist;
895 ins_d_ang =ang ;
896 ins_d_k =k ;
897 }
898 }
899 }
900 ENDLOOPGRIDINNER;
901 /* loop over gridcells adjacent to j (xk,yk,zk) */
902 LOOPGRIDINNER(xk,yk,zk,xkk,ykk,zkk,xj,yj,zj,ngrid) {
903 kcell=&grid[xk][yk][zk][grIH];
904 /* loop over hydrogen atoms from ins group
905 in this adjacent gridcell (kcell) */
906 for (ak=0; ak<kcell->nr; ak++) {
907 k=kcell->atoms[ak];
908 if (is_hbond(a[ grID][k],
909 a[ grIH][k],
910 a[OGRP+grA][j],
911 rcut,ccut,x,bBox,hbox,&dist,&ang))
912 if (dist<ins_a_dist) {
913 ins_a=TRUE;
914 ins_a_dist=dist;
915 ins_a_ang =ang ;
916 ins_a_k =k ;
917 }
918 }
919 }
920 ENDLOOPGRIDINNER;
921
922 if (ins_d && ins_a &&
923 (a[grIA][ins_d_k] == a[grID][ins_a_k])) {
924 /* set insertion index */
925 insert[a[grID][ins_a_k]]=TRUE;
926 insert[a[grIH][ins_a_k]]=TRUE;
927 insert[a[grIA][ins_d_k]]=TRUE;
928
929 /* add to hbond index if not already there */
930 if (idx==NOTSET) {
931 if (donors[grp][i].nrhb>=donors[grp][i].maxnr) {
932 donors[grp][i].maxnr+=10;
933 srenew(donors[grp][i].hb,donors[grp][i].maxnr);
934 }
935 donors[grp][i].hb[donors[grp][i].nrhb].a=j;
936 donors[grp][i].hb[donors[grp][i].nrhb].nr=nrhb;
937 idx=donors[grp][i].nrhb;
938 donors[grp][i].nrhb++;
939 if (bHBMap && (nrhb>=max_nrhb) ) {
940 max_nrhb+=HBINC;
941 srenew(hbexist,max_nrhb);
942 for (l=max_nrhb-HBINC; l<max_nrhb; l++)
943 snew(hbexist[l],max_nframes);
944 }
945 nrhb++;
946 }
947
948 if (bHBMap)
949 /* mark insertion in hbond index */
950 hbexist[donors[grp][i].hb[idx].nr][nframes] |=
951 HB_INS;
952
953 /* print insertion info to file */
954 fprintf(fpins,
955 "%4g: %4u:%3.3s%4d%3.3s -> "
956 "%4u:%3.3s%4d%3.3s (%4.2f,%2.0f) - "
957 "%4u:%3.3s%4d%3.3s (%4.2f,%2.0f)\n",t,
958 a[grIA][ins_d_k]+1,
959 *top.atoms.resname[top.atoms.atom[a[grIA][ins_d_k]].resnr],
960 top.atoms.atom[a[grIA][ins_d_k]].resnr+1,
961 *top.atoms.atomname[a[grIA][ins_d_k]],
962 a[grp+grD][i]+1,
963 *top.atoms.resname[top.atoms.atom[a[grp+grD][i]].resnr],
964 top.atoms.atom[a[grp+grD][i]].resnr+1,
965 *top.atoms.atomname[a[grp+grD][i]],
966 ins_d_dist,ins_d_ang*RAD2DEG,
967 a[OGRP+grA][j]+1,
968 *top.atoms.resname[top.atoms.atom[a[OGRP+grA][j]].resnr],
969 top.atoms.atom[a[OGRP+grA][j]].resnr+1,
970 *top.atoms.atomname[a[OGRP+grA][j]],
971 ins_a_dist,ins_a_ang*RAD2DEG);
972 }
973 }
974 }
975 } /* for aj */
976 }
977 ENDLOOPGRIDINNER;
978 } /* for ai */
979 } /* for grp */
980 } /* for xi,yi,zi */
981 nframes++;
982 } while (read_next_x(status,&t,natoms,x,box));
983
984 free_grid(ngrid,&grid);
985
986 close_trj(status);
987 if (bInsert)
988 fclose(fpins);
989
990 if (nrhb==0)
991 fprintf(stderr,"No hydrogen bonds found!!\n");
992 else {
993 fprintf(stderr,"Found %d different hydrogen bonds in trajectory\n",nrhb);
994
995 /* Dump everything to output */
996 sort_hb(nr_a,donors);
997
998 fp = xvgropen(opt2fn("-num",NFILE,fnm),"Hydrogen Bonds","Time","Number");
999 for(i=0; i<nframes; i++)
1000 fprintf(fp,"%10g %10d\n",time[i],nhb[i]);
1001 fclose(fp);
1002
1003 if (opt2bSet("-dist",NFILE,fnm)) {
1004 long sum;
1005
1006 sum=0;
1007 for(i=0; i<nrbin; i++)
1008 sum+=rdist[i];
1009
1010 fp = xvgropen(opt2fn("-dist",NFILE,fnm),
1011 "Hydrogen Bond Distribution",
1012 "Hydrogen - Acceptor Distance (nm)","");
1013 for(i=0; i<nrbin; i++)
1014 fprintf(fp,"%10g %10g\n",(i+0.5)*rbin,rdist[i]/(rbin*(real)sum));
1015 fclose(fp);
1016 }
1017
1018 if (opt2bSet("-ang",NFILE,fnm)) {
1019 long sum;
1020
1021 sum=0;
1022 for(i=0; i<nabin; i++)
1023 sum+=adist[i];
1024
1025 fp = xvgropen(opt2fn("-ang",NFILE,fnm),
1026 "Hydrogen Bond Distribution",
1027 "Donor - Hydrogen - Acceptor Angle (\\SO\\N)","");
1028 for(i=0; i<nabin; i++)
1029 fprintf(fp,"%10g %10g\n",(i+0.5)*abin,adist[i]/(abin*(real)sum));
1030 fclose(fp);
1031 }
1032
1033 if (opt2bSet("-hx",NFILE,fnm)) {
1034 fp = xvgropen(opt2fn("-hx",NFILE,fnm),
1035 "Hydrogen Bonds","Time(ps)","Count");
1036 xvgr_legend(fp,NRHXTYPES,hxtypenames);
1037 for(i=0; i<nframes; i++) {
1038 fprintf(fp,"%10g",time[i]);
1039 for (j=0; j<max_hx; j++)
1040 fprintf(fp," %6d",nhx[i][j]);
1041 fprintf(fp,"\n");
1042 }
1043 fclose(fp);
1044 }
1045
1046
1047 if (opt2bSet("-ac",NFILE,fnm)) {
1048 /* build hbexist matrix in reals for autocorr */
1049 snew(rhbex,nrhb);
1050 for(i=0; i<nrhb; i++) {
1051 snew(rhbex[i],nframes);
1052 for(j=0; j<nframes; j++)
1053 rhbex[i][j]=(hbexist[i][j] & HB_YES);
1054 }
1055 low_do_autocorr(opt2fn("-ac",NFILE,fnm), "Hydrogen Bond Autocorrelation",
1056 nframes,nrhb,-1,rhbex,time[1]-time[0],eacNormal,1,
1057 TRUE,TRUE,TRUE,FALSE,0,0,0);
1058 for(i=0; i<nrhb; i++)
1059 sfree(rhbex[i]);
1060 sfree(rhbex);
1061 }
1062
1063 if (opt2bSet("-hbm",NFILE,fnm)) {
1064 t_matrix mat;
1065 int x,y;
1066
1067 mat.nx=nframes;
1068 mat.ny=nrhb;
1069 snew(mat.matrix,mat.nx);
1070 for(x=0; x<mat.nx; x++){
1071 snew(mat.matrix[x],mat.ny);
1072 y=0;
1073 for (grp=gr0; grp<=(bTwo?gr1:gr0); grp+=grINC)
1074 for(i=0; i<nr_a[grp+grD]; i++)
1075 for(j=0; j<donors[grp][i].nrhb; j++)
1076 mat.matrix[x][y++]=hbexist[donors[grp][i].hb[j].nr][x];
1077 }
1078 mat.axis_x=time;
1079 snew(mat.axis_y,mat.ny);
1080 for(j=0; j<mat.ny; j++)
1081 mat.axis_y[j]=j;
1082 sprintf(mat.title,"Hydrogen Bond Existence Map");
1083 sprintf(mat.legend,"Hydrogen Bonds");
1084 sprintf(mat.label_x,"time (ps)");
1085 sprintf(mat.label_y,"Hydrogen Bond Index");
1086 mat.bDiscrete=TRUE;
1087 if (bInsert)
1088 mat.nmap=HB_NR;
1089 else
1090 mat.nmap=2;
1091 snew(mat.map,mat.nmap);
1092 for(i=0; i<mat.nmap; i++) {
1093 mat.map[i].code.c1=hbmap[i];
1094 mat.map[i].desc=hbdesc[i];
1095 mat.map[i].rgb=hbrgb[i];
1096 }
1097 fp = opt2FILE("-hbm",NFILE,fnm,"w");
1098 write_xpm_m(fp, mat);
1099 fclose(fp);
1100 for(x=0; x<mat.nx; x++)
1101 sfree(mat.matrix[x]);
1102 sfree(mat.axis_y);
1103 sfree(mat.matrix);
1104 sfree(mat.map);
1105 }
1106 }
1107
1108 if (opt2bSet("-hbn",NFILE,fnm)) {
1109 fp = opt2FILE("-hbn",NFILE,fnm,"w");
1110 for (grp=gr0; grp<=(bTwo?gr1:gr0); grp+=grINC) {
1111 fprintf(fp,"[ %s ]",grpnames[grp/grINC]);
1112 for (i=0; i<isize[grp/grINC]; i++) {
1113 fprintf(fp,(i%15)?" ":"\n");
1114 fprintf(fp,"%4u",index[grp/grINC][i]+1);
1115 }
1116 fprintf(fp,"\n");
1117 fprintf(fp,"[ donors_hydrogens_%s ]",grpnames[grp/grINC]);
1118 for (i=0; i<nr_a[grp+grD]; i++) {
1119 fprintf(fp,(i%6)?" ":"\n");
1120 fprintf(fp,"%4u %4u",a[grp+grD][i]+1,a[grp+grH][i]+1);
1121 }
1122 fprintf(fp,"\n");
1123 fprintf(fp,"[ acceptors_%s ]",grpnames[grp/grINC]);
1124 for (i=0; i<nr_a[grp+grA]; i++) {
1125 fprintf(fp,(i%15)?" ":"\n");
1126 fprintf(fp,"%4u",a[grp+grA][i]+1);
1127 }
1128 fprintf(fp,"\n");
1129 }
1130 if (bTwo)
1131 fprintf(fp,"[ hbonds_%s-%s ]\n",grpnames[0],grpnames[1]);
1132 else
1133 fprintf(fp,"[ hbonds_%s ]\n",grpnames[0]);
1134 for (grp=gr0; grp<=(bTwo?gr1:gr0); grp+=grINC)
1135 for(i=0; i<nr_a[grp+grD]; i++)
1136 for(j=0; j<donors[grp][i].nrhb; j++)
1137 fprintf(fp,"%6u %6u %6u\n",
1138 a[ grp+grD][i]+1,
1139 a[ grp+grH][i]+1,
1140 a[OGRP+grA][donors[grp][i].hb[j].a]+1);
1141 if (bInsert) {
1142 if (bTwo)
1143 fprintf(fp,"[ insert_%s->%s-%s ]",
1144 grpnames[2],grpnames[0],grpnames[1]);
1145 else
1146 fprintf(fp,"[ insert_%s->%s ]",grpnames[2],grpnames[0]);
1147 j=0;
1148 for(i=0; i<natoms; i++)
1149 if (insert[i]) {
1150 fprintf(fp,(j%15)?" ":"\n");
1151 fprintf(fp,"%4d",i+1);
1152 j++;
1153 }
1154 fprintf(fp,"\n");
1155 }
1156 fclose(fp);
1157 }
1158
1159 if (bDAnr) {
1160 int i,j,nleg;
1161 char **legnames;
1162 char buf[STRLEN];
1163 char *danames[grINC] = {"Donors", NULL, "Acceptors"};
1164
1165 #define USE_THIS_GROUP(j) ( ( j % grINC != grH ) && ( bTwo || ( j / grINC != 1 ) ) && ( bInsert || ( j / grINC != 2 ) ) )
1166
1167 fp = xvgropen(opt2fn("-da",NFILE,fnm),
1168 "Donors and Acceptors","Time(ps)","Count");
1169 nleg = (bTwo?2:1)*2 + (bInsert?2:0);
1170 snew(legnames,nleg);
1171 i=0;
1172 for(j=0; j<grNR; j++)
1173 if ( USE_THIS_GROUP(j) ) {
1174 sprintf(buf,"%s %s",danames[j % grINC],grpnames[j / grINC]);
1175 legnames[i++] = strdup(buf);
1176 }
1177 assert(i==nleg);
1178 xvgr_legend(fp,nleg,legnames);
1179 for(i=0; i<nframes; i++) {
1180 fprintf(fp,"%10g",time[i]);
1181 for (j=0; j<grNR; j++)
1182 if ( USE_THIS_GROUP(j) )
1183 fprintf(fp," %6d",danr[i][j]);
1184 fprintf(fp,"\n");
1185 }
1186 fclose(fp);
1187 }
1188
1189 thanx(stdout);
1190
1191 return 0;
1192 }
got merge a long time ago