src/tools/anadih.c

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