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 t0,real dt,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, t0, dt, 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 t0, real dt, 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
 156   /* Analysis of dihedral transitions */
 157   fprintf(stderr,"Now calculating transitions...\n");
 158
 159   if (bRb)
 160     calc_bin=calc_RBbin;
 161   else
 162     calc_bin=calc_Nbin;
 163
 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;
 168   }
 169   snew(tr_h,nangles);
 170   snew(tr_f,nframes);
 171
 172   /* dih[i][j] is the dihedral angle i in frame j  */
 173   ntrans = 0;
 174   for (i=0; (i<nangles); i++)
 175   {
 176
 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++)
 185     {
 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++;
 196       }
 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;
 203
 204       prev = dih[i][j];
 205
 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++;
 214       }
 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);
 224   }
 225
 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 */
 229
 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++ ;
 241       }
 242     }
 243   }
 244
 245   /* end addition by grs */
 246
 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]);
 253     }
 254     ffclose(fp);
 255   }
 256
 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--)
 264     ;
 265
 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]);
 273     }
 274     ffclose(fp);
 275   }
 276
 277   sfree(tr_f);
 278   sfree(tr_h);
 279   for(k=0;k<NROT;k++)
 280     sfree(rot_occ[k]);
 281
 282 }
 283
 284 void mk_multiplicity_lookup (int *xity, int maxchi, real **dih,
 285                              int nlist, t_dlist dlist[],int nangles)
 286 {
 287   /* new by grs - for dihedral j (as in dih[j]) get multiplicity from dlist
 288    * and store in xity[j]
 289    */
 290
 291   int j, Dih, i ;
 292   char name[4];
 293
 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 ;
 304
 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 ;
 308         }
 309
 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;
 323           }
 324         }
 325         j++ ;
 326       }
 327     }
 328   }
 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;
 335
 336 }
 337
 338 void mk_chi_lookup (int **lookup, int maxchi, real **dih,
 339                     int nlist, t_dlist dlist[])
 340 {
 341
 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..)*/
 346
 347   int i,j, Dih, Chi;
 348
 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 ;
 359         }
 360         j++ ;
 361       } else {
 362         lookup[i][Chi] = -1 ;
 363       }
 364     }
 365   }
 366
 367 }
 368
 369
 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,gmx_bool bRb, gmx_bool bNormalize,
 373                            real core_frac, gmx_bool bAll, const char *fnall,
 374                            const output_env_t oenv)
 375 {
 376
 377   gmx_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;
 384
 385   int  (*calc_bin)(real,int,real);
 386
 387   /* Analysis of dihedral transitions */
 388   fprintf(stderr,"Now calculating Chi product trajectories...\n");
 389
 390   if (bRb)
 391     calc_bin=calc_RBbin;
 392   else
 393     calc_bin=calc_Nbin;
 394
 395   snew(chi_prtrj, nframes) ;
 396
 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);
 402
 403   for(i=0; (i<nlist); i++) {
 404
 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 ;
 412       }
 413     }
 414     nbin += 1 ; /* for the "zero rotamer", outside the core region */
 415
 416     for (j=0; (j<nframes); j++) {
 417
 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 ;
 438           }
 439         }
 440         accum ++ ;
 441       }
 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 ;
 448       }
 449     }
 450     if (bHaveChi){
 451
 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);
 456       }
 457
 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]);
 474         }
 475         fprintf(fp,"&\n");
 476         ffclose(fp);
 477       }
 478
 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]);
 489       }
 490       fprintf(fpall, "\n") ;
 491
 492       sfree(chi_prhist);
 493       /* histogram done */
 494     }
 495   }
 496
 497   sfree(chi_prtrj);
 498   ffclose(fpall);
 499   fprintf(stderr,"\n") ;
 500
 501 }
 502
 503 void calc_distribution_props(int nh,int histo[],real start,
 504                              int nkkk, t_karplus kkk[],
 505                              real *S2)
 506 {
 507   real d,dc,ds,c1,c2,tdc,tds;
 508   real fac,ang,invth,Jc;
 509   int  i,j,th;
 510
 511   if (nh == 0)
 512     gmx_fatal(FARGS,"No points in histogram (%s, %d)",__FILE__,__LINE__);
 513   fac = 2*M_PI/nh;
 514
 515   /* Compute normalisation factor */
 516   th=0;
 517   for(j=0; (j<nh); j++)
 518     th+=histo[j];
 519   invth=1.0/th;
 520
 521   for(i=0; (i<nkkk); i++) {
 522     kkk[i].Jc    = 0;
 523     kkk[i].Jcsig = 0;
 524   }
 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);
 541     }
 542   }
 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));
 546   }
 547   *S2 = tdc*tdc+tds*tds;
 548 }
 549
 550 static void calc_angles(FILE *log,t_pbc *pbc,
 551                         int n3,atom_id index[],real ang[],rvec x_s[])
 552 {
 553   int  i,ix,t1,t2;
 554   rvec r_ij,r_kj;
 555   real costh;
 556
 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);
 565   }
 566 }
 567
 568 static real calc_fraction(real angles[], int nangles)
 569 {
 570   int i;
 571   real trans = 0, gauche = 0;
 572   real angle;
 573
 574   for (i = 0; i < nangles; i++)
 575   {
 576     angle = angles[i] * RAD2DEG;
 577
 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;
 584   }
 585   if (trans+gauche > 0)
 586     return trans/(trans+gauche);
 587   else
 588     return 0;
 589 }
 590
 591 static void calc_dihs(FILE *log,t_pbc *pbc,
 592                       int n4,atom_id index[],real ang[],rvec x_s[])
 593 {
 594   int  i,ix,t1,t2,t3;
 595   rvec r_ij,r_kj,r_kl,m,n;
 596   real sign,aaa;
 597
 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                   &sign,&t1,&t2,&t3);
 603
 604     ang[i]=aaa;  /* not taking into account ryckaert bellemans yet */
 605   }
 606 }
 607
 608 void make_histo(FILE *log,
 609                 int ndata,real data[],int npoints,int histo[],
 610                 real minx,real maxx)
 611 {
 612   double dx;
 613   int    i,ind;
 614
 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]);
 620     }
 621     fprintf(log,"Min data: %10g  Max data: %10g\n",minx,maxx);
 622   }
 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]);
 634   }
 635 }
 636
 637 void normalize_histo(int npoints,int histo[],real dx,real normhisto[])
 638 {
 639   int    i;
 640   double d,fac;
 641
 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;
 648   }
 649   fac=1.0/d;
 650   for(i=0; (i<npoints); i++)
 651     normhisto[i]=fac*histo[i];
 652 }
 653
 654 void read_ang_dih(const char *trj_fn,
 655                   gmx_bool bAngles,gmx_bool bSaveAll,gmx_bool bRb,gmx_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)
 663 {
 664   t_pbc      *pbc;
 665   t_trxstatus *status;
 666   int        i,angind,natoms,total,teller;
 667   int        nangles,n_alloc;
 668   real       t,fraction,pifac,aa,angle;
 669   real       *angles[2];
 670   matrix     box;
 671   rvec       *x;
 672   int        cur=0;
 673 #define prev (1-cur)
 674
 675   snew(pbc,1);
 676   natoms = read_first_x(oenv,&status,trj_fn,&t,&x,box);
 677
 678   if (bAngles) {
 679     nangles=isize/3;
 680     pifac=M_PI;
 681   }
 682   else {
 683     nangles=isize/4;
 684     pifac=2.0*M_PI;
 685   }
 686   snew(angles[cur],nangles);
 687   snew(angles[prev],nangles);
 688
 689   /* Start the loop over frames */
 690   total       = 0;
 691   teller      = 0;
 692   n_alloc     = 0;
 693   *time       = NULL;
 694   *trans_frac = NULL;
 695   *aver_angle = NULL;
 696
 697   do {
 698     if (teller >= n_alloc) {
 699       n_alloc+=100;
 700       if (bSaveAll)
 701         for (i=0; (i<nangles); i++)
 702           srenew(dih[i],n_alloc);
 703       srenew(*time,n_alloc);
 704       srenew(*trans_frac,n_alloc);
 705       srenew(*aver_angle,n_alloc);
 706     }
 707
 708     (*time)[teller] = t;
 709
 710     if (pbc)
 711       set_pbc(pbc,-1,box);
 712
 713     if (bAngles) {
 714       calc_angles(stdout,pbc,isize,index,angles[cur],x);
 715     }
 716     else {
 717       calc_dihs(stdout,pbc,isize,index,angles[cur],x);
 718
 719       /* Trans fraction */
 720       fraction = calc_fraction(angles[cur], nangles);
 721       (*trans_frac)[teller] = fraction;
 722
 723       /* Change Ryckaert-Bellemans dihedrals to polymer convention
 724        * Modified 990913 by Erik:
 725        * We actually shouldn't change the convention, since it's
 726        * calculated from polymer above, but we change the intervall
 727        * from [-180,180] to [0,360].
 728        */
 729       if (bRb) {
 730         for(i=0; (i<nangles); i++)
 731           if (angles[cur][i] <= 0.0)
 732             angles[cur][i] += 2*M_PI;
 733         }
 734
 735       /* Periodicity in dihedral space... */
 736       if (bPBC) {
 737         for(i=0; (i<nangles); i++) {
 738           real dd = angles[cur][i];
 739           angles[cur][i] = atan2(sin(dd),cos(dd));
 740         }
 741       }
 742       else {
 743         if (teller > 1) {
 744           for(i=0; (i<nangles); i++) {
 745             while (angles[cur][i] <= angles[prev][i] - M_PI)
 746               angles[cur][i]+=2*M_PI;
 747             while (angles[cur][i] > angles[prev][i] + M_PI)
 748               angles[cur][i]-=2*M_PI;
 749           }
 750         }
 751       }
 752     }
 753
 754     /* Average angles */
 755     aa=0;
 756     for(i=0; (i<nangles); i++) {
 757       aa=aa+angles[cur][i];
 758
 759       /* angle in rad / 2Pi * max determines bin. bins go from 0 to maxangstat,
 760          even though scale goes from -pi to pi (dihedral) or -pi/2 to pi/2
 761          (angle) Basically: translate the x-axis by Pi. Translate it back by
 762          -Pi when plotting.
 763        */
 764
 765       angle = angles[cur][i];
 766       if (!bAngles) {
 767         while (angle < -M_PI)
 768           angle += 2*M_PI;
 769         while (angle >= M_PI)
 770           angle -= 2*M_PI;
 771
 772         angle+=M_PI;
 773       }
 774
 775       /* Update the distribution histogram */
 776       angind = (int) ((angle*maxangstat)/pifac + 0.5);
 777       if (angind==maxangstat)
 778         angind=0;
 779       if ( (angind < 0) || (angind >= maxangstat) )
 780         /* this will never happen */
 781         gmx_fatal(FARGS,"angle (%f) index out of range (0..%d) : %d\n",
 782                     angle,maxangstat,angind);
 783
 784       angstat[angind]++;
 785       if (angind==maxangstat)
 786         fprintf(stderr,"angle %d fr %d = %g\n",i,cur,angle);
 787
 788       total++;
 789     }
 790
 791     /* average over all angles */
 792     (*aver_angle)[teller] = (aa/nangles);
 793
 794     /* this copies all current dih. angles to dih[i], teller is frame */
 795     if (bSaveAll)
 796       for (i = 0; i < nangles; i++)
 797         dih[i][teller] = angles[cur][i];
 798
 799     /* Swap buffers */
 800     cur=prev;
 801
 802     /* Increment loop counter */
 803     teller++;
 804   } while (read_next_x(oenv,status,&t,natoms,x,box));
 805   close_trj(status);
 806
 807   sfree(x);
 808   sfree(angles[cur]);
 809   sfree(angles[prev]);
 810
 811   *nframes=teller;
 812 }
 813