Don't use POSIX fnmatch() for pattern matching.
[gromacs/qmmm-gamess-us.git] / src / tools / anadih.c
blob0e0b1da7d14ae30e718fec8f2129013f6b9ac69f
1 /*
2 *
3 * This source code is part of
4 *
5 * G R O M A C S
6 *
7 * GROningen MAchine for Chemical Simulations
8 *
9 * VERSION 3.2.0
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
32 * And Hey:
33 * Green Red Orange Magenta Azure Cyan Skyblue
35 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
39 #include <math.h>
40 #include <stdio.h>
41 #include <string.h>
42 #include "physics.h"
43 #include "smalloc.h"
44 #include "macros.h"
45 #include "txtdump.h"
46 #include "bondf.h"
47 #include "xvgr.h"
48 #include "typedefs.h"
49 #include "vec.h"
50 #include "gstat.h"
51 #include "confio.h"
53 void print_one(const output_env_t oenv,const char *base,const char *name,
54 const char *title, const char *ylabel,int nf,real time[],
55 real data[])
57 FILE *fp;
58 char buf[256],t2[256];
59 int k;
61 sprintf(buf,"%s%s.xvg",base,name);
62 fprintf(stderr,"\rPrinting %s ",buf);
63 sprintf(t2,"%s %s",title,name);
64 fp=xvgropen(buf,t2,"Time (ps)",ylabel,oenv);
65 for(k=0; (k<nf); k++)
66 fprintf(fp,"%10g %10g\n",time[k],data[k]);
67 ffclose(fp);
70 static int calc_RBbin(real phi, int multiplicity, real core_frac)
72 /* multiplicity and core_frac NOT used,
73 * just given to enable use of pt-to-fn in caller low_ana_dih_trans*/
74 static const real r30 = M_PI/6.0;
75 static const real r90 = M_PI/2.0;
76 static const real r150 = M_PI*5.0/6.0;
78 if ((phi < r30) && (phi > -r30))
79 return 1;
80 else if ((phi > -r150) && (phi < -r90))
81 return 2;
82 else if ((phi < r150) && (phi > r90))
83 return 3;
84 return 0;
87 static int calc_Nbin(real phi, int multiplicity, real core_frac)
89 static const real r360 = 360*DEG2RAD;
90 real rot_width, core_width, core_offset, low, hi;
91 int bin ;
92 /* with multiplicity 3 and core_frac 0.5
93 * 0<g(-)<120, 120<t<240, 240<g(+)<360
94 * 0< bin0 < 30, 30<bin1<90, 90<bin0<150, 150<bin2<210, 210<bin0<270, 270<bin3<330, 330<bin0<360
95 * so with multiplicity 3, bin1 is core g(-), bin2 is core t, bin3 is
96 core g(+), bin0 is between rotamers */
97 if (phi < 0)
98 phi += r360;
100 rot_width = 360/multiplicity ;
101 core_width = core_frac * rot_width ;
102 core_offset = (rot_width - core_width)/2.0 ;
103 for(bin = 1 ; bin <= multiplicity ; bin ++ ) {
104 low = ((bin - 1) * rot_width ) + core_offset ;
105 hi = ((bin - 1) * rot_width ) + core_offset + core_width;
106 low *= DEG2RAD ;
107 hi *= DEG2RAD ;
108 if ((phi > low) && (phi < hi))
109 return bin ;
111 return 0;
114 void ana_dih_trans(const char *fn_trans,const char *fn_histo,
115 real **dih,int nframes,int nangles,
116 const char *grpname,real t0,real dt,bool bRb,
117 const output_env_t oenv)
119 /* just a wrapper; declare extra args, then chuck away at end. */
120 int maxchi = 0 ;
121 t_dlist *dlist ;
122 int *xity;
123 int nlist = nangles ;
124 int k ;
126 snew(dlist,nlist);
127 snew(xity,nangles);
128 for(k=0; (k<nangles); k++) {
129 xity[k]=3 ;
132 low_ana_dih_trans(TRUE, fn_trans,TRUE, fn_histo, maxchi,
133 dih, nlist, dlist, nframes,
134 nangles, grpname, xity, t0, dt, bRb, 0.5,oenv);
135 sfree(dlist);
136 sfree(xity);
140 void low_ana_dih_trans(bool bTrans, const char *fn_trans,
141 bool bHisto, const char *fn_histo, int maxchi,
142 real **dih, int nlist, t_dlist dlist[], int nframes,
143 int nangles, const char *grpname, int xity[],
144 real t0, real dt, bool bRb, real core_frac,
145 const output_env_t oenv)
147 FILE *fp;
148 int *tr_f,*tr_h;
149 char title[256];
150 int i,j,k,Dih,ntrans;
151 int cur_bin,new_bin;
152 real ttime,tt;
153 real *rot_occ[NROT] ;
154 int (*calc_bin)(real,int,real);
156 /* Analysis of dihedral transitions */
157 fprintf(stderr,"Now calculating transitions...\n");
159 if (bRb)
160 calc_bin=calc_RBbin;
161 else
162 calc_bin=calc_Nbin;
164 for(k=0;k<NROT;k++) {
165 snew(rot_occ[k],nangles);
166 for (i=0; (i<nangles); i++)
167 rot_occ[k][i]=0;
169 snew(tr_h,nangles);
170 snew(tr_f,nframes);
172 /* dih[i][j] is the dihedral angle i in frame j */
173 ntrans = 0;
174 for (i=0; (i<nangles); i++)
177 /*#define OLDIE*/
178 #ifdef OLDIE
179 mind = maxd = prev = dih[i][0];
180 #else
181 cur_bin = calc_bin(dih[i][0],xity[i],core_frac);
182 rot_occ[cur_bin][i]++ ;
183 #endif
184 for (j=1; (j<nframes); j++)
186 new_bin = calc_bin(dih[i][j],xity[i],core_frac);
187 rot_occ[new_bin][i]++ ;
188 #ifndef OLDIE
189 if (cur_bin == 0)
190 cur_bin=new_bin;
191 else if ((new_bin != 0) && (cur_bin != new_bin)) {
192 cur_bin = new_bin;
193 tr_f[j]++;
194 tr_h[i]++;
195 ntrans++;
197 #else
198 /* why is all this md rubbish periodic? Remove 360 degree periodicity */
199 if ( (dih[i][j] - prev) > M_PI)
200 dih[i][j] -= 2*M_PI;
201 else if ( (dih[i][j] - prev) < -M_PI)
202 dih[i][j] += 2*M_PI;
204 prev = dih[i][j];
206 mind = min(mind, dih[i][j]);
207 maxd = max(maxd, dih[i][j]);
208 if ( (maxd - mind) > 2*M_PI/3) /* or 120 degrees, assuming */
209 { /* multiplicity 3. Not so general.*/
210 tr_f[j]++;
211 tr_h[i]++;
212 maxd = mind = dih[i][j]; /* get ready for next transition */
213 ntrans++;
215 #endif
216 } /* end j */
217 for(k=0;k<NROT;k++)
218 rot_occ[k][i] /= nframes ;
219 } /* end i */
220 fprintf(stderr,"Total number of transitions: %10d\n",ntrans);
221 if (ntrans > 0) {
222 ttime = (dt*nframes*nangles)/ntrans;
223 fprintf(stderr,"Time between transitions: %10.3f ps\n",ttime);
226 /* new by grs - copy transitions from tr_h[] to dlist->ntr[]
227 * and rotamer populations from rot_occ to dlist->rot_occ[]
228 * based on fn histogramming in g_chi. diff roles for i and j here */
230 j=0;
231 for (Dih=0; (Dih<NONCHI+maxchi); Dih++) {
232 for(i=0; (i<nlist); i++) {
233 if (((Dih < edOmega) ) ||
234 ((Dih == edOmega) && (has_dihedral(edOmega,&(dlist[i])))) ||
235 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1))) {
236 /* grs debug printf("Not OK? i %d j %d Dih %d \n", i, j, Dih) ; */
237 dlist[i].ntr[Dih] = tr_h[j] ;
238 for(k=0;k<NROT;k++)
239 dlist[i].rot_occ[Dih][k] = rot_occ[k][j] ;
240 j++ ;
245 /* end addition by grs */
247 if (bTrans) {
248 sprintf(title,"Number of transitions: %s",grpname);
249 fp=xvgropen(fn_trans,title,"Time (ps)","# transitions/timeframe",oenv);
250 for(j=0; (j<nframes); j++) {
251 tt = t0+j*dt;
252 fprintf(fp,"%10.3f %10d\n",tt,tr_f[j]);
254 ffclose(fp);
257 /* Compute histogram from # transitions per dihedral */
258 /* Use old array */
259 for(j=0; (j<nframes); j++)
260 tr_f[j]=0;
261 for(i=0; (i<nangles); i++)
262 tr_f[tr_h[i]]++;
263 for(j=nframes; ((tr_f[j-1] == 0) && (j>0)); j--)
266 ttime = dt*nframes;
267 if (bHisto) {
268 sprintf(title,"Transition time: %s",grpname);
269 fp=xvgropen(fn_histo,title,"Time (ps)","#",oenv);
270 for(i=j-1; (i>0); i--) {
271 if (tr_f[i] != 0)
272 fprintf(fp,"%10.3f %10d\n",ttime/i,tr_f[i]);
274 ffclose(fp);
277 sfree(tr_f);
278 sfree(tr_h);
279 for(k=0;k<NROT;k++)
280 sfree(rot_occ[k]);
284 void mk_multiplicity_lookup (int *xity, int maxchi, real **dih,
285 int nlist, t_dlist dlist[],int nangles)
287 /* new by grs - for dihedral j (as in dih[j]) get multiplicity from dlist
288 * and store in xity[j]
291 int j, Dih, i ;
292 char name[4];
294 j=0;
295 for (Dih=0; (Dih<NONCHI+maxchi); Dih++) {
296 for(i=0; (i<nlist); i++) {
297 strncpy(name, dlist[i].name,3) ;
298 name[3]='\0' ;
299 if (((Dih < edOmega) ) ||
300 ((Dih == edOmega) && (has_dihedral(edOmega,&(dlist[i])))) ||
301 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1))) {
302 /* default - we will correct the rest below */
303 xity[j] = 3 ;
305 /* make omegas 2fold, though doesn't make much more sense than 3 */
306 if (Dih == edOmega && (has_dihedral(edOmega,&(dlist[i])))) {
307 xity[j] = 2 ;
310 /* dihedrals to aromatic rings, COO, CONH2 or guanidinium are 2fold*/
311 if (Dih > edOmega && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)) {
312 if ( ((strstr(name,"PHE") != NULL) && (Dih == edChi2)) ||
313 ((strstr(name,"TYR") != NULL) && (Dih == edChi2)) ||
314 ((strstr(name,"PTR") != NULL) && (Dih == edChi2)) ||
315 ((strstr(name,"TRP") != NULL) && (Dih == edChi2)) ||
316 ((strstr(name,"HIS") != NULL) && (Dih == edChi2)) ||
317 ((strstr(name,"GLU") != NULL) && (Dih == edChi3)) ||
318 ((strstr(name,"ASP") != NULL) && (Dih == edChi2)) ||
319 ((strstr(name,"GLN") != NULL) && (Dih == edChi3)) ||
320 ((strstr(name,"ASN") != NULL) && (Dih == edChi2)) ||
321 ((strstr(name,"ARG") != NULL) && (Dih == edChi4)) ) {
322 xity[j] = 2;
325 j++ ;
329 if (j<nangles)
330 fprintf(stderr,"WARNING: not all dihedrals found in topology (only %d out of %d)!\n",
331 j,nangles);
332 /* Check for remaining dihedrals */
333 for(;(j < nangles); j++)
334 xity[j] = 3;
338 void mk_chi_lookup (int **lookup, int maxchi, real **dih,
339 int nlist, t_dlist dlist[])
342 /* by grs. should rewrite everything to use this. (but haven't,
343 * and at mmt only used in get_chi_product_traj
344 * returns the dihed number given the residue number (from-0)
345 * and chi (from-0) nr. -1 for chi undefined for that res (eg gly, ala..)*/
347 int i,j, Dih, Chi;
349 j=0;
350 for (Dih=0; (Dih<NONCHI+maxchi); Dih++) {
351 for(i=0; (i<nlist); i++) {
352 Chi = Dih - NONCHI ;
353 if (((Dih < edOmega) ) ||
354 ((Dih == edOmega) && (has_dihedral(edOmega,&(dlist[i])))) ||
355 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1))) {
356 /* grs debug printf("Not OK? i %d j %d Dih %d \n", i, j, Dih) ; */
357 if (Dih > edOmega ) {
358 lookup[i][Chi] = j ;
360 j++ ;
361 } else {
362 lookup[i][Chi] = -1 ;
370 void get_chi_product_traj (real **dih,int nframes,int nangles, int nlist,
371 int maxchi, t_dlist dlist[], real time[],
372 int **lookup, int *xity,bool bRb, bool bNormalize,
373 real core_frac, bool bAll, const char *fnall,
374 const output_env_t oenv)
377 bool bRotZero, bHaveChi=FALSE;
378 int accum=0, index, i,j,k,Xi,n,b ;
379 real *chi_prtrj;
380 int *chi_prhist;
381 int nbin ;
382 FILE *fp, *fpall ;
383 char hisfile[256],histitle[256], *namept;
385 int (*calc_bin)(real,int,real);
387 /* Analysis of dihedral transitions */
388 fprintf(stderr,"Now calculating Chi product trajectories...\n");
390 if (bRb)
391 calc_bin=calc_RBbin;
392 else
393 calc_bin=calc_Nbin;
395 snew(chi_prtrj, nframes) ;
397 /* file for info on all residues */
398 if (bNormalize)
399 fpall=xvgropen(fnall,"Cumulative Rotamers","Residue","Probability",oenv);
400 else
401 fpall=xvgropen(fnall,"Cumulative Rotamers","Residue","# Counts",oenv);
403 for(i=0; (i<nlist); i++) {
405 /* get nbin, the nr. of cumulative rotamers that need to be considered */
406 nbin = 1 ;
407 for (Xi = 0 ; Xi < maxchi ; Xi ++ ) {
408 index = lookup[i][Xi] ; /* chi_(Xi+1) of res i (-1 if off end) */
409 if ( index >= 0 ) {
410 n = xity[index];
411 nbin = n*nbin ;
414 nbin += 1 ; /* for the "zero rotamer", outside the core region */
416 for (j=0; (j<nframes); j++) {
418 bRotZero = FALSE ;
419 bHaveChi = TRUE ;
420 index = lookup[i][0] ; /* index into dih of chi1 of res i */
421 if (index == -1 ) {
422 b = 0 ;
423 bRotZero = TRUE ;
424 bHaveChi = FALSE ;
425 } else {
426 b = calc_bin(dih[index][j],xity[index],core_frac) ;
427 accum = b - 1 ;
428 if (b == 0 )
429 bRotZero = TRUE ;
430 for (Xi = 1 ; Xi < maxchi ; Xi ++ ) {
431 index = lookup[i][Xi] ; /* chi_(Xi+1) of res i (-1 if off end) */
432 if ( index >= 0 ) {
433 n = xity[index];
434 b = calc_bin(dih[index][j],n,core_frac);
435 accum = n * accum + b - 1 ;
436 if (b == 0 )
437 bRotZero = TRUE ;
440 accum ++ ;
442 if (bRotZero) {
443 chi_prtrj[j] = 0.0 ;
444 } else {
445 chi_prtrj[j] = accum ;
446 if (accum+1 > nbin )
447 nbin = accum+1 ;
450 if (bHaveChi){
452 if (bAll) {
453 /* print cuml rotamer vs time */
454 print_one(oenv,"chiproduct", dlist[i].name, "chi product for",
455 "cumulative rotamer", nframes,time,chi_prtrj);
458 /* make a histogram pf culm. rotamer occupancy too */
459 snew(chi_prhist, nbin) ;
460 make_histo(NULL,nframes,chi_prtrj,nbin,chi_prhist,0,nbin);
461 if (bAll) {
462 sprintf(hisfile,"histo-chiprod%s.xvg",dlist[i].name);
463 sprintf(histitle,"cumulative rotamer distribution for %s",dlist[i].name);
464 fprintf(stderr," and %s ",hisfile);
465 fp=xvgropen(hisfile,histitle,"number","",oenv);
466 fprintf(fp,"@ xaxis tick on\n");
467 fprintf(fp,"@ xaxis tick major 1\n");
468 fprintf(fp,"@ type xy\n");
469 for(k=0; (k<nbin); k++) {
470 if (bNormalize)
471 fprintf(fp,"%5d %10g\n",k,(1.0*chi_prhist[k])/nframes);
472 else
473 fprintf(fp,"%5d %10d\n",k,chi_prhist[k]);
475 fprintf(fp,"&\n");
476 ffclose(fp);
479 /* and finally print out occupancies to a single file */
480 /* get the gmx from-1 res nr by setting a ptr to the number part
481 * of dlist[i].name - potential bug for 4-letter res names... */
482 namept = dlist[i].name + 3 ;
483 fprintf(fpall, "%5s ", namept);
484 for(k=0; (k<nbin); k++) {
485 if (bNormalize)
486 fprintf(fpall," %10g",(1.0*chi_prhist[k])/nframes);
487 else
488 fprintf(fpall," %10d",chi_prhist[k]);
490 fprintf(fpall, "\n") ;
492 sfree(chi_prhist);
493 /* histogram done */
497 sfree(chi_prtrj);
498 ffclose(fpall);
499 fprintf(stderr,"\n") ;
503 void calc_distribution_props(int nh,int histo[],real start,
504 int nkkk, t_karplus kkk[],
505 real *S2)
507 real d,dc,ds,c1,c2,tdc,tds;
508 real fac,ang,invth,Jc;
509 int i,j,th;
511 if (nh == 0)
512 gmx_fatal(FARGS,"No points in histogram (%s, %d)",__FILE__,__LINE__);
513 fac = 2*M_PI/nh;
515 /* Compute normalisation factor */
516 th=0;
517 for(j=0; (j<nh); j++)
518 th+=histo[j];
519 invth=1.0/th;
521 for(i=0; (i<nkkk); i++) {
522 kkk[i].Jc = 0;
523 kkk[i].Jcsig = 0;
525 tdc=0,tds=0;
526 for(j=0; (j<nh); j++) {
527 d = invth*histo[j];
528 ang = j*fac-start;
529 c1 = cos(ang);
530 c2 = c1*c1;
531 dc = d*c1;
532 ds = d*sin(ang);
533 tdc += dc;
534 tds += ds;
535 for(i=0; (i<nkkk); i++) {
536 c1 = cos(ang+kkk[i].offset);
537 c2 = c1*c1;
538 Jc = (kkk[i].A*c2 + kkk[i].B*c1 + kkk[i].C);
539 kkk[i].Jc += histo[j]*Jc;
540 kkk[i].Jcsig += histo[j]*sqr(Jc);
543 for(i=0; (i<nkkk); i++) {
544 kkk[i].Jc /= th;
545 kkk[i].Jcsig = sqrt(kkk[i].Jcsig/th-sqr(kkk[i].Jc));
547 *S2 = tdc*tdc+tds*tds;
550 static void calc_angles(FILE *log,t_pbc *pbc,
551 int n3,atom_id index[],real ang[],rvec x_s[])
553 int i,ix,t1,t2;
554 rvec r_ij,r_kj;
555 real costh;
557 for(i=ix=0; (ix<n3); i++,ix+=3)
558 ang[i]=bond_angle(x_s[index[ix]],x_s[index[ix+1]],x_s[index[ix+2]],
559 pbc,r_ij,r_kj,&costh,&t1,&t2);
560 if (debug) {
561 fprintf(debug,"Angle[0]=%g, costh=%g, index0 = %d, %d, %d\n",
562 ang[0],costh,index[0],index[1],index[2]);
563 pr_rvec(debug,0,"rij",r_ij,DIM,TRUE);
564 pr_rvec(debug,0,"rkj",r_kj,DIM,TRUE);
568 static real calc_fraction(real angles[], int nangles)
570 int i;
571 real trans = 0, gauche = 0;
572 real angle;
574 for (i = 0; i < nangles; i++)
576 angle = angles[i] * RAD2DEG;
578 if (angle > 135 && angle < 225)
579 trans += 1.0;
580 else if (angle > 270 && angle < 330)
581 gauche += 1.0;
582 else if (angle < 90 && angle > 30)
583 gauche += 1.0;
585 if (trans+gauche > 0)
586 return trans/(trans+gauche);
587 else
588 return 0;
591 static void calc_dihs(FILE *log,t_pbc *pbc,
592 int n4,atom_id index[],real ang[],rvec x_s[])
594 int i,ix,t1,t2,t3;
595 rvec r_ij,r_kj,r_kl,m,n;
596 real cos_phi,sign,aaa;
598 for(i=ix=0; (ix<n4); i++,ix+=4) {
599 aaa=dih_angle(x_s[index[ix]],x_s[index[ix+1]],x_s[index[ix+2]],
600 x_s[index[ix+3]],pbc,
601 r_ij,r_kj,r_kl,m,n,
602 &cos_phi,&sign,&t1,&t2,&t3);
604 ang[i]=aaa; /* not taking into account ryckaert bellemans yet */
608 void make_histo(FILE *log,
609 int ndata,real data[],int npoints,int histo[],
610 real minx,real maxx)
612 double dx;
613 int i,ind;
615 if (minx == maxx) {
616 minx=maxx=data[0];
617 for(i=1; (i<ndata); i++) {
618 minx=min(minx,data[i]);
619 maxx=max(maxx,data[i]);
621 fprintf(log,"Min data: %10g Max data: %10g\n",minx,maxx);
623 dx=(double)npoints/(maxx-minx);
624 if (debug)
625 fprintf(debug,
626 "Histogramming: ndata=%d, nhisto=%d, minx=%g,maxx=%g,dx=%g\n",
627 ndata,npoints,minx,maxx,dx);
628 for(i=0; (i<ndata); i++) {
629 ind=(data[i]-minx)*dx;
630 if ((ind >= 0) && (ind < npoints))
631 histo[ind]++;
632 else
633 fprintf(log,"index = %d, data[%d] = %g\n",ind,i,data[i]);
637 void normalize_histo(int npoints,int histo[],real dx,real normhisto[])
639 int i;
640 double d,fac;
642 d=0;
643 for(i=0; (i<npoints); i++)
644 d+=dx*histo[i];
645 if (d==0) {
646 fprintf(stderr,"Empty histogram!\n");
647 return;
649 fac=1.0/d;
650 for(i=0; (i<npoints); i++)
651 normhisto[i]=fac*histo[i];
654 void read_ang_dih(const char *trj_fn,
655 bool bAngles,bool bSaveAll,bool bRb,bool bPBC,
656 int maxangstat,int angstat[],
657 int *nframes,real **time,
658 int isize,atom_id index[],
659 real **trans_frac,
660 real **aver_angle,
661 real *dih[],
662 const output_env_t oenv)
664 t_pbc *pbc;
665 int i,angind,status,natoms,total,teller;
666 int nangles,n_alloc;
667 real t,fraction,pifac,aa,angle;
668 real *angles[2];
669 matrix box;
670 rvec *x;
671 int cur=0;
672 #define prev (1-cur)
674 snew(pbc,1);
675 natoms = read_first_x(oenv,&status,trj_fn,&t,&x,box);
677 if (bAngles) {
678 nangles=isize/3;
679 pifac=M_PI;
681 else {
682 nangles=isize/4;
683 pifac=2.0*M_PI;
685 snew(angles[cur],nangles);
686 snew(angles[prev],nangles);
688 /* Start the loop over frames */
689 total = 0;
690 teller = 0;
691 n_alloc = 0;
692 *time = NULL;
693 *trans_frac = NULL;
694 *aver_angle = NULL;
696 do {
697 if (teller >= n_alloc) {
698 n_alloc+=100;
699 if (bSaveAll)
700 for (i=0; (i<nangles); i++)
701 srenew(dih[i],n_alloc);
702 srenew(*time,n_alloc);
703 srenew(*trans_frac,n_alloc);
704 srenew(*aver_angle,n_alloc);
707 (*time)[teller] = t;
709 if (pbc)
710 set_pbc(pbc,-1,box);
712 if (bAngles) {
713 calc_angles(stdout,pbc,isize,index,angles[cur],x);
715 else {
716 calc_dihs(stdout,pbc,isize,index,angles[cur],x);
718 /* Trans fraction */
719 fraction = calc_fraction(angles[cur], nangles);
720 (*trans_frac)[teller] = fraction;
722 /* Change Ryckaert-Bellemans dihedrals to polymer convention
723 * Modified 990913 by Erik:
724 * We actually shouldn't change the convention, since it's
725 * calculated from polymer above, but we change the intervall
726 * from [-180,180] to [0,360].
728 if (bRb) {
729 for(i=0; (i<nangles); i++)
730 if (angles[cur][i] <= 0.0)
731 angles[cur][i] += 2*M_PI;
734 /* Periodicity in dihedral space... */
735 if (bPBC) {
736 for(i=0; (i<nangles); i++) {
737 real dd = angles[cur][i];
738 angles[cur][i] = atan2(sin(dd),cos(dd));
741 else {
742 if (teller > 1) {
743 for(i=0; (i<nangles); i++) {
744 while (angles[cur][i] <= angles[prev][i] - M_PI)
745 angles[cur][i]+=2*M_PI;
746 while (angles[cur][i] > angles[prev][i] + M_PI)
747 angles[cur][i]-=2*M_PI;
753 /* Average angles */
754 aa=0;
755 for(i=0; (i<nangles); i++) {
756 aa=aa+angles[cur][i];
758 /* angle in rad / 2Pi * max determines bin. bins go from 0 to maxangstat,
759 even though scale goes from -pi to pi (dihedral) or -pi/2 to pi/2
760 (angle) Basically: translate the x-axis by Pi. Translate it back by
761 -Pi when plotting.
764 angle = angles[cur][i];
765 if (!bAngles) {
766 while (angle < -M_PI)
767 angle += 2*M_PI;
768 while (angle >= M_PI)
769 angle -= 2*M_PI;
771 angle+=M_PI;
774 /* Update the distribution histogram */
775 angind = (int) ((angle*maxangstat)/pifac + 0.5);
776 if (angind==maxangstat)
777 angind=0;
778 if ( (angind < 0) || (angind >= maxangstat) )
779 /* this will never happen */
780 gmx_fatal(FARGS,"angle (%f) index out of range (0..%d) : %d\n",
781 angle,maxangstat,angind);
783 angstat[angind]++;
784 if (angind==maxangstat)
785 fprintf(stderr,"angle %d fr %d = %g\n",i,cur,angle);
787 total++;
790 /* average over all angles */
791 (*aver_angle)[teller] = (aa/nangles);
793 /* this copies all current dih. angles to dih[i], teller is frame */
794 if (bSaveAll)
795 for (i = 0; i < nangles; i++)
796 dih[i][teller] = angles[cur][i];
798 /* Swap buffers */
799 cur=prev;
801 /* Increment loop counter */
802 teller++;
803 } while (read_next_x(oenv,status,&t,natoms,x,box));
804 close_trj(status);
806 sfree(x);
807 sfree(angles[cur]);
808 sfree(angles[prev]);
810 *nframes=teller;