src/tools/anadih.c

   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.
  14
  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.
  19  *
  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.
  26  *
  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.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  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"
  52
  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[])
  56 {
  57   FILE *fp;
  58   char buf[256],t2[256];
  59   int  k;
  60
  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);
  68 }
  69
  70 static int calc_RBbin(real phi, int multiplicity, real core_frac)
  71 {
  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;
  77
  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;
  85 }
  86
  87 static int calc_Nbin(real phi, int multiplicity, real core_frac)
  88 {
  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;
  99
 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 ;
 110  }
 111  return 0;
 112 }
 113
 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 *time,gmx_bool bRb,
 117                    const output_env_t oenv)
 118 {
 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 ;
 125
 126   snew(dlist,nlist);
 127   snew(xity,nangles);
 128   for(k=0; (k<nangles); k++) {
 129     xity[k]=3 ;
 130   }
 131
 132   low_ana_dih_trans(TRUE, fn_trans,TRUE, fn_histo, maxchi,
 133                     dih, nlist, dlist, nframes,
 134                     nangles, grpname, xity, time, bRb, 0.5,oenv);
 135   sfree(dlist);
 136   sfree(xity);
 137
 138 }
 139
 140 void low_ana_dih_trans(gmx_bool bTrans, const char *fn_trans,
 141                        gmx_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 *time, gmx_bool bRb, real core_frac,
 145                        const output_env_t oenv)
 146 {
 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);
 155   real dt;
 156
 157   if (1 <= nframes) {
 158     return;
 159   }
 160   /* Assumes the frames are equally spaced in time */
 161   dt=(time[nframes-1]-time[0])/(nframes-1);
 162
 163   /* Analysis of dihedral transitions */
 164   fprintf(stderr,"Now calculating transitions...\n");
 165
 166   if (bRb)
 167     calc_bin=calc_RBbin;
 168   else
 169     calc_bin=calc_Nbin;
 170
 171   for(k=0;k<NROT;k++) {
 172     snew(rot_occ[k],nangles);
 173     for (i=0; (i<nangles); i++)
 174       rot_occ[k][i]=0;
 175   }
 176   snew(tr_h,nangles);
 177   snew(tr_f,nframes);
 178
 179   /* dih[i][j] is the dihedral angle i in frame j  */
 180   ntrans = 0;
 181   for (i=0; (i<nangles); i++)
 182   {
 183
 184     /*#define OLDIE*/
 185 #ifdef OLDIE
 186     mind = maxd = prev = dih[i][0];
 187 #else
 188     cur_bin = calc_bin(dih[i][0],xity[i],core_frac);
 189     rot_occ[cur_bin][i]++ ;
 190 #endif
 191     for (j=1; (j<nframes); j++)
 192     {
 193       new_bin = calc_bin(dih[i][j],xity[i],core_frac);
 194       rot_occ[new_bin][i]++ ;
 195 #ifndef OLDIE
 196       if (cur_bin == 0)
 197         cur_bin=new_bin;
 198       else if ((new_bin != 0) && (cur_bin != new_bin)) {
 199         cur_bin = new_bin;
 200         tr_f[j]++;
 201         tr_h[i]++;
 202         ntrans++;
 203       }
 204 #else
 205       /* why is all this md rubbish periodic? Remove 360 degree periodicity */
 206       if ( (dih[i][j] - prev) > M_PI)
 207         dih[i][j] -= 2*M_PI;
 208       else if ( (dih[i][j] - prev) < -M_PI)
 209         dih[i][j] += 2*M_PI;
 210
 211       prev = dih[i][j];
 212
 213       mind = min(mind, dih[i][j]);
 214       maxd = max(maxd, dih[i][j]);
 215       if ( (maxd - mind) > 2*M_PI/3)    /* or 120 degrees, assuming       */
 216       {                                 /* multiplicity 3. Not so general.*/
 217         tr_f[j]++;
 218         tr_h[i]++;
 219         maxd = mind = dih[i][j];        /* get ready for next transition  */
 220         ntrans++;
 221       }
 222 #endif
 223     } /* end j */
 224     for(k=0;k<NROT;k++)
 225       rot_occ[k][i] /= nframes ;
 226   } /* end i */
 227   fprintf(stderr,"Total number of transitions: %10d\n",ntrans);
 228   if (ntrans > 0) {
 229     ttime = (dt*nframes*nangles)/ntrans;
 230     fprintf(stderr,"Time between transitions:    %10.3f ps\n",ttime);
 231   }
 232
 233   /* new by grs - copy transitions from tr_h[] to dlist->ntr[]
 234    * and rotamer populations from rot_occ to dlist->rot_occ[]
 235    * based on fn histogramming in g_chi. diff roles for i and j here */
 236
 237   j=0;
 238   for (Dih=0; (Dih<NONCHI+maxchi); Dih++) {
 239     for(i=0; (i<nlist); i++) {
 240       if (((Dih  < edOmega) ) ||
 241           ((Dih == edOmega) && (has_dihedral(edOmega,&(dlist[i])))) ||
 242           ((Dih  > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1))) {
 243         /* grs debug  printf("Not OK? i %d j %d Dih %d \n", i, j, Dih) ; */
 244         dlist[i].ntr[Dih] = tr_h[j] ;
 245         for(k=0;k<NROT;k++)
 246           dlist[i].rot_occ[Dih][k] = rot_occ[k][j] ;
 247         j++ ;
 248       }
 249     }
 250   }
 251
 252   /* end addition by grs */
 253
 254   if (bTrans) {
 255     sprintf(title,"Number of transitions: %s",grpname);
 256     fp=xvgropen(fn_trans,title,"Time (ps)","# transitions/timeframe",oenv);
 257     for(j=0; (j<nframes); j++) {
 258       fprintf(fp,"%10.3f  %10d\n",time[j],tr_f[j]);
 259     }
 260     ffclose(fp);
 261   }
 262
 263   /* Compute histogram from # transitions per dihedral */
 264   /* Use old array */
 265   for(j=0; (j<nframes); j++)
 266     tr_f[j]=0;
 267   for(i=0; (i<nangles); i++)
 268     tr_f[tr_h[i]]++;
 269   for(j=nframes; ((tr_f[j-1] == 0) && (j>0)); j--)
 270     ;
 271
 272   ttime = dt*nframes;
 273   if (bHisto) {
 274     sprintf(title,"Transition time: %s",grpname);
 275     fp=xvgropen(fn_histo,title,"Time (ps)","#",oenv);
 276     for(i=j-1; (i>0); i--) {
 277       if (tr_f[i] != 0)
 278         fprintf(fp,"%10.3f  %10d\n",ttime/i,tr_f[i]);
 279     }
 280     ffclose(fp);
 281   }
 282
 283   sfree(tr_f);
 284   sfree(tr_h);
 285   for(k=0;k<NROT;k++)
 286     sfree(rot_occ[k]);
 287
 288 }
 289
 290 void mk_multiplicity_lookup (int *xity, int maxchi, real **dih,
 291                              int nlist, t_dlist dlist[],int nangles)
 292 {
 293   /* new by grs - for dihedral j (as in dih[j]) get multiplicity from dlist
 294    * and store in xity[j]
 295    */
 296
 297   int j, Dih, i ;
 298   char name[4];
 299
 300   j=0;
 301   for (Dih=0; (Dih<NONCHI+maxchi); Dih++) {
 302     for(i=0; (i<nlist); i++) {
 303       strncpy(name, dlist[i].name,3) ;
 304       name[3]='\0' ;
 305       if (((Dih  < edOmega) ) ||
 306           ((Dih == edOmega) && (has_dihedral(edOmega,&(dlist[i])))) ||
 307           ((Dih  > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1))) {
 308         /* default - we will correct the rest below */
 309         xity[j] = 3 ;
 310
 311         /* make omegas 2fold, though doesn't make much more sense than 3 */
 312         if (Dih == edOmega && (has_dihedral(edOmega,&(dlist[i])))) {
 313           xity[j] = 2 ;
 314         }
 315
 316         /* dihedrals to aromatic rings, COO, CONH2 or guanidinium are 2fold*/
 317         if (Dih > edOmega && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)) {
 318           if ( ((strstr(name,"PHE") != NULL) && (Dih == edChi2))  ||
 319                ((strstr(name,"TYR") != NULL) && (Dih == edChi2))  ||
 320                ((strstr(name,"PTR") != NULL) && (Dih == edChi2))  ||
 321                ((strstr(name,"TRP") != NULL) && (Dih == edChi2))  ||
 322                ((strstr(name,"HIS") != NULL) && (Dih == edChi2))  ||
 323                ((strstr(name,"GLU") != NULL) && (Dih == edChi3))  ||
 324                ((strstr(name,"ASP") != NULL) && (Dih == edChi2))  ||
 325                ((strstr(name,"GLN") != NULL) && (Dih == edChi3))  ||
 326                ((strstr(name,"ASN") != NULL) && (Dih == edChi2))  ||
 327                ((strstr(name,"ARG") != NULL) && (Dih == edChi4))  ) {
 328             xity[j] = 2;
 329           }
 330         }
 331         j++ ;
 332       }
 333     }
 334   }
 335   if (j<nangles)
 336     fprintf(stderr,"WARNING: not all dihedrals found in topology (only %d out of %d)!\n",
 337             j,nangles);
 338   /* Check for remaining dihedrals */
 339   for(;(j < nangles); j++)
 340     xity[j] = 3;
 341
 342 }
 343
 344 void mk_chi_lookup (int **lookup, int maxchi, real **dih,
 345                     int nlist, t_dlist dlist[])
 346 {
 347
 348   /* by grs. should rewrite everything to use this. (but haven't,
 349    * and at mmt only used in get_chi_product_traj
 350    * returns the dihed number given the residue number (from-0)
 351    * and chi (from-0) nr. -1 for chi undefined for that res (eg gly, ala..)*/
 352
 353   int i,j, Dih, Chi;
 354
 355   j=0;
 356   for (Dih=0; (Dih<NONCHI+maxchi); Dih++) {
 357     for(i=0; (i<nlist); i++) {
 358       Chi = Dih - NONCHI ;
 359       if (((Dih  < edOmega) ) ||
 360           ((Dih == edOmega) && (has_dihedral(edOmega,&(dlist[i])))) ||
 361           ((Dih  > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1))) {
 362         /* grs debug  printf("Not OK? i %d j %d Dih %d \n", i, j, Dih) ; */
 363         if (Dih > edOmega ) {
 364           lookup[i][Chi] = j ;
 365         }
 366         j++ ;
 367       } else {
 368         lookup[i][Chi] = -1 ;
 369       }
 370     }
 371   }
 372
 373 }
 374
 375
 376 void get_chi_product_traj (real **dih,int nframes,int nangles, int nlist,
 377                            int maxchi, t_dlist dlist[], real time[],
 378                            int **lookup, int *xity,gmx_bool bRb, gmx_bool bNormalize,
 379                            real core_frac, gmx_bool bAll, const char *fnall,
 380                            const output_env_t oenv)
 381 {
 382
 383   gmx_bool bRotZero, bHaveChi=FALSE;
 384   int  accum=0, index, i,j,k,Xi,n,b ;
 385   real *chi_prtrj;
 386   int  *chi_prhist;
 387   int  nbin ;
 388   FILE *fp, *fpall ;
 389   char hisfile[256],histitle[256], *namept;
 390
 391   int  (*calc_bin)(real,int,real);
 392
 393   /* Analysis of dihedral transitions */
 394   fprintf(stderr,"Now calculating Chi product trajectories...\n");
 395
 396   if (bRb)
 397     calc_bin=calc_RBbin;
 398   else
 399     calc_bin=calc_Nbin;
 400
 401   snew(chi_prtrj, nframes) ;
 402
 403   /* file for info on all residues */
 404   if (bNormalize)
 405     fpall=xvgropen(fnall,"Cumulative Rotamers","Residue","Probability",oenv);
 406   else
 407     fpall=xvgropen(fnall,"Cumulative Rotamers","Residue","# Counts",oenv);
 408
 409   for(i=0; (i<nlist); i++) {
 410
 411     /* get nbin, the nr. of cumulative rotamers that need to be considered */
 412     nbin = 1 ;
 413     for (Xi = 0 ; Xi < maxchi ; Xi ++ ) {
 414       index = lookup[i][Xi] ; /* chi_(Xi+1) of res i (-1 if off end) */
 415       if ( index >= 0 ) {
 416         n = xity[index];
 417         nbin = n*nbin ;
 418       }
 419     }
 420     nbin += 1 ; /* for the "zero rotamer", outside the core region */
 421
 422     for (j=0; (j<nframes); j++) {
 423
 424       bRotZero = FALSE ;
 425       bHaveChi = TRUE ;
 426       index = lookup[i][0] ; /* index into dih of chi1 of res i */
 427       if (index == -1 ) {
 428         b = 0 ;
 429         bRotZero = TRUE ;
 430         bHaveChi = FALSE ;
 431       } else {
 432         b = calc_bin(dih[index][j],xity[index],core_frac) ;
 433         accum = b - 1 ;
 434         if (b == 0 )
 435           bRotZero = TRUE ;
 436         for (Xi = 1 ; Xi < maxchi ; Xi ++ ) {
 437           index = lookup[i][Xi] ; /* chi_(Xi+1) of res i (-1 if off end) */
 438           if ( index >= 0 ) {
 439             n = xity[index];
 440             b = calc_bin(dih[index][j],n,core_frac);
 441             accum = n * accum + b - 1 ;
 442             if (b == 0 )
 443               bRotZero = TRUE ;
 444           }
 445         }
 446         accum ++ ;
 447       }
 448       if (bRotZero) {
 449         chi_prtrj[j] = 0.0 ;
 450       } else {
 451         chi_prtrj[j] = accum ;
 452         if (accum+1 > nbin )
 453           nbin = accum+1 ;
 454       }
 455     }
 456     if (bHaveChi){
 457
 458       if (bAll) {
 459         /* print cuml rotamer vs time */
 460         print_one(oenv,"chiproduct", dlist[i].name, "chi product for",
 461                   "cumulative rotamer", nframes,time,chi_prtrj);
 462       }
 463
 464       /* make a histogram pf culm. rotamer occupancy too */
 465       snew(chi_prhist, nbin) ;
 466       make_histo(NULL,nframes,chi_prtrj,nbin,chi_prhist,0,nbin);
 467       if (bAll) {
 468         sprintf(hisfile,"histo-chiprod%s.xvg",dlist[i].name);
 469         sprintf(histitle,"cumulative rotamer distribution for %s",dlist[i].name);
 470         fprintf(stderr,"  and %s  ",hisfile);
 471         fp=xvgropen(hisfile,histitle,"number","",oenv);
 472         fprintf(fp,"@ xaxis tick on\n");
 473         fprintf(fp,"@ xaxis tick major 1\n");
 474         fprintf(fp,"@ type xy\n");
 475         for(k=0; (k<nbin); k++) {
 476           if (bNormalize)
 477             fprintf(fp,"%5d  %10g\n",k,(1.0*chi_prhist[k])/nframes);
 478           else
 479             fprintf(fp,"%5d  %10d\n",k,chi_prhist[k]);
 480         }
 481         fprintf(fp,"&\n");
 482         ffclose(fp);
 483       }
 484
 485       /* and finally print out occupancies to a single file */
 486       /* get the gmx from-1 res nr by setting a ptr to the number part
 487        * of dlist[i].name - potential bug for 4-letter res names... */
 488       namept = dlist[i].name + 3 ;
 489       fprintf(fpall, "%5s ", namept);
 490       for(k=0; (k<nbin); k++) {
 491         if (bNormalize)
 492           fprintf(fpall,"  %10g",(1.0*chi_prhist[k])/nframes);
 493         else
 494           fprintf(fpall,"  %10d",chi_prhist[k]);
 495       }
 496       fprintf(fpall, "\n") ;
 497
 498       sfree(chi_prhist);
 499       /* histogram done */
 500     }
 501   }
 502
 503   sfree(chi_prtrj);
 504   ffclose(fpall);
 505   fprintf(stderr,"\n") ;
 506
 507 }
 508
 509 void calc_distribution_props(int nh,int histo[],real start,
 510                              int nkkk, t_karplus kkk[],
 511                              real *S2)
 512 {
 513   real d,dc,ds,c1,c2,tdc,tds;
 514   real fac,ang,invth,Jc;
 515   int  i,j,th;
 516
 517   if (nh == 0)
 518     gmx_fatal(FARGS,"No points in histogram (%s, %d)",__FILE__,__LINE__);
 519   fac = 2*M_PI/nh;
 520
 521   /* Compute normalisation factor */
 522   th=0;
 523   for(j=0; (j<nh); j++)
 524     th+=histo[j];
 525   invth=1.0/th;
 526
 527   for(i=0; (i<nkkk); i++) {
 528     kkk[i].Jc    = 0;
 529     kkk[i].Jcsig = 0;
 530   }
 531   tdc=0,tds=0;
 532   for(j=0; (j<nh); j++) {
 533     d    = invth*histo[j];
 534     ang  = j*fac-start;
 535     c1   = cos(ang);
 536     c2   = c1*c1;
 537     dc   = d*c1;
 538     ds   = d*sin(ang);
 539     tdc += dc;
 540     tds += ds;
 541     for(i=0; (i<nkkk); i++) {
 542       c1   = cos(ang+kkk[i].offset);
 543       c2   = c1*c1;
 544       Jc   = (kkk[i].A*c2 + kkk[i].B*c1 + kkk[i].C);
 545       kkk[i].Jc    += histo[j]*Jc;
 546       kkk[i].Jcsig += histo[j]*sqr(Jc);
 547     }
 548   }
 549   for(i=0; (i<nkkk); i++) {
 550     kkk[i].Jc    /= th;
 551     kkk[i].Jcsig  = sqrt(kkk[i].Jcsig/th-sqr(kkk[i].Jc));
 552   }
 553   *S2 = tdc*tdc+tds*tds;
 554 }
 555
 556 static void calc_angles(FILE *log,t_pbc *pbc,
 557                         int n3,atom_id index[],real ang[],rvec x_s[])
 558 {
 559   int  i,ix,t1,t2;
 560   rvec r_ij,r_kj;
 561   real costh;
 562
 563   for(i=ix=0; (ix<n3); i++,ix+=3)
 564     ang[i]=bond_angle(x_s[index[ix]],x_s[index[ix+1]],x_s[index[ix+2]],
 565                       pbc,r_ij,r_kj,&costh,&t1,&t2);
 566   if (debug) {
 567     fprintf(debug,"Angle[0]=%g, costh=%g, index0 = %d, %d, %d\n",
 568             ang[0],costh,index[0],index[1],index[2]);
 569     pr_rvec(debug,0,"rij",r_ij,DIM,TRUE);
 570     pr_rvec(debug,0,"rkj",r_kj,DIM,TRUE);
 571   }
 572 }
 573
 574 static real calc_fraction(real angles[], int nangles)
 575 {
 576   int i;
 577   real trans = 0, gauche = 0;
 578   real angle;
 579
 580   for (i = 0; i < nangles; i++)
 581   {
 582     angle = angles[i] * RAD2DEG;
 583
 584     if (angle > 135 && angle < 225)
 585       trans += 1.0;
 586     else if (angle > 270 && angle < 330)
 587       gauche += 1.0;
 588     else if (angle < 90 && angle > 30)
 589       gauche += 1.0;
 590   }
 591   if (trans+gauche > 0)
 592     return trans/(trans+gauche);
 593   else
 594     return 0;
 595 }
 596
 597 static void calc_dihs(FILE *log,t_pbc *pbc,
 598                       int n4,atom_id index[],real ang[],rvec x_s[])
 599 {
 600   int  i,ix,t1,t2,t3;
 601   rvec r_ij,r_kj,r_kl,m,n;
 602   real sign,aaa;
 603
 604   for(i=ix=0; (ix<n4); i++,ix+=4) {
 605     aaa=dih_angle(x_s[index[ix]],x_s[index[ix+1]],x_s[index[ix+2]],
 606                   x_s[index[ix+3]],pbc,
 607                   r_ij,r_kj,r_kl,m,n,
 608                   &sign,&t1,&t2,&t3);
 609
 610     ang[i]=aaa;  /* not taking into account ryckaert bellemans yet */
 611   }
 612 }
 613
 614 void make_histo(FILE *log,
 615                 int ndata,real data[],int npoints,int histo[],
 616                 real minx,real maxx)
 617 {
 618   double dx;
 619   int    i,ind;
 620
 621   if (minx == maxx) {
 622     minx=maxx=data[0];
 623     for(i=1; (i<ndata); i++) {
 624       minx=min(minx,data[i]);
 625       maxx=max(maxx,data[i]);
 626     }
 627     fprintf(log,"Min data: %10g  Max data: %10g\n",minx,maxx);
 628   }
 629   dx=(double)npoints/(maxx-minx);
 630   if (debug)
 631     fprintf(debug,
 632             "Histogramming: ndata=%d, nhisto=%d, minx=%g,maxx=%g,dx=%g\n",
 633             ndata,npoints,minx,maxx,dx);
 634   for(i=0; (i<ndata); i++) {
 635     ind=(data[i]-minx)*dx;
 636     if ((ind >= 0) && (ind < npoints))
 637       histo[ind]++;
 638     else
 639       fprintf(log,"index = %d, data[%d] = %g\n",ind,i,data[i]);
 640   }
 641 }
 642
 643 void normalize_histo(int npoints,int histo[],real dx,real normhisto[])
 644 {
 645   int    i;
 646   double d,fac;
 647
 648   d=0;
 649   for(i=0; (i<npoints); i++)
 650     d+=dx*histo[i];
 651   if (d==0) {
 652     fprintf(stderr,"Empty histogram!\n");
 653     return;
 654   }
 655   fac=1.0/d;
 656   for(i=0; (i<npoints); i++)
 657     normhisto[i]=fac*histo[i];
 658 }
 659
 660 void read_ang_dih(const char *trj_fn,
 661                   gmx_bool bAngles,gmx_bool bSaveAll,gmx_bool bRb,gmx_bool bPBC,
 662                   int maxangstat,int angstat[],
 663                   int *nframes,real **time,
 664                   int isize,atom_id index[],
 665                   real **trans_frac,
 666                   real **aver_angle,
 667                   real *dih[],
 668                   const output_env_t oenv)
 669 {
 670   t_pbc      *pbc;
 671   t_trxstatus *status;
 672   int        i,angind,natoms,total,teller;
 673   int        nangles,n_alloc;
 674   real       t,fraction,pifac,aa,angle;
 675   real       *angles[2];
 676   matrix     box;
 677   rvec       *x;
 678   int        cur=0;
 679 #define prev (1-cur)
 680
 681   snew(pbc,1);
 682   natoms = read_first_x(oenv,&status,trj_fn,&t,&x,box);
 683
 684   if (bAngles) {
 685     nangles=isize/3;
 686     pifac=M_PI;
 687   }
 688   else {
 689     nangles=isize/4;
 690     pifac=2.0*M_PI;
 691   }
 692   snew(angles[cur],nangles);
 693   snew(angles[prev],nangles);
 694
 695   /* Start the loop over frames */
 696   total       = 0;
 697   teller      = 0;
 698   n_alloc     = 0;
 699   *time       = NULL;
 700   *trans_frac = NULL;
 701   *aver_angle = NULL;
 702
 703   do {
 704     if (teller >= n_alloc) {
 705       n_alloc+=100;
 706       if (bSaveAll)
 707         for (i=0; (i<nangles); i++)
 708           srenew(dih[i],n_alloc);
 709       srenew(*time,n_alloc);
 710       srenew(*trans_frac,n_alloc);
 711       srenew(*aver_angle,n_alloc);
 712     }
 713
 714     (*time)[teller] = t;
 715
 716     if (pbc)
 717       set_pbc(pbc,-1,box);
 718
 719     if (bAngles) {
 720       calc_angles(stdout,pbc,isize,index,angles[cur],x);
 721     }
 722     else {
 723       calc_dihs(stdout,pbc,isize,index,angles[cur],x);
 724
 725       /* Trans fraction */
 726       fraction = calc_fraction(angles[cur], nangles);
 727       (*trans_frac)[teller] = fraction;
 728
 729       /* Change Ryckaert-Bellemans dihedrals to polymer convention
 730        * Modified 990913 by Erik:
 731        * We actually shouldn't change the convention, since it's
 732        * calculated from polymer above, but we change the intervall
 733        * from [-180,180] to [0,360].
 734        */
 735       if (bRb) {
 736         for(i=0; (i<nangles); i++)
 737           if (angles[cur][i] <= 0.0)
 738             angles[cur][i] += 2*M_PI;
 739         }
 740
 741       /* Periodicity in dihedral space... */
 742       if (bPBC) {
 743         for(i=0; (i<nangles); i++) {
 744           real dd = angles[cur][i];
 745           angles[cur][i] = atan2(sin(dd),cos(dd));
 746         }
 747       }
 748       else {
 749         if (teller > 1) {
 750           for(i=0; (i<nangles); i++) {
 751             while (angles[cur][i] <= angles[prev][i] - M_PI)
 752               angles[cur][i]+=2*M_PI;
 753             while (angles[cur][i] > angles[prev][i] + M_PI)
 754               angles[cur][i]-=2*M_PI;
 755           }
 756         }
 757       }
 758     }
 759
 760     /* Average angles */
 761     aa=0;
 762     for(i=0; (i<nangles); i++) {
 763       aa=aa+angles[cur][i];
 764
 765       /* angle in rad / 2Pi * max determines bin. bins go from 0 to maxangstat,
 766          even though scale goes from -pi to pi (dihedral) or -pi/2 to pi/2
 767          (angle) Basically: translate the x-axis by Pi. Translate it back by
 768          -Pi when plotting.
 769        */
 770
 771       angle = angles[cur][i];
 772       if (!bAngles) {
 773         while (angle < -M_PI)
 774           angle += 2*M_PI;
 775         while (angle >= M_PI)
 776           angle -= 2*M_PI;
 777
 778         angle+=M_PI;
 779       }
 780
 781       /* Update the distribution histogram */
 782       angind = (int) ((angle*maxangstat)/pifac + 0.5);
 783       if (angind==maxangstat)
 784         angind=0;
 785       if ( (angind < 0) || (angind >= maxangstat) )
 786         /* this will never happen */
 787         gmx_fatal(FARGS,"angle (%f) index out of range (0..%d) : %d\n",
 788                     angle,maxangstat,angind);
 789
 790       angstat[angind]++;
 791       if (angind==maxangstat)
 792         fprintf(stderr,"angle %d fr %d = %g\n",i,cur,angle);
 793
 794       total++;
 795     }
 796
 797     /* average over all angles */
 798     (*aver_angle)[teller] = (aa/nangles);
 799
 800     /* this copies all current dih. angles to dih[i], teller is frame */
 801     if (bSaveAll)
 802       for (i = 0; i < nangles; i++)
 803         dih[i][teller] = angles[cur][i];
 804
 805     /* Swap buffers */
 806     cur=prev;
 807
 808     /* Increment loop counter */
 809     teller++;
 810   } while (read_next_x(oenv,status,&t,natoms,x,box));
 811   close_trj(status);
 812
 813   sfree(x);
 814   sfree(angles[cur]);
 815   sfree(angles[prev]);
 816
 817   *nframes=teller;
 818 }
 819