src/tools/anadih.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_anadih_c = "$Id$";
  30
  31 #include <math.h>
  32 #include <stdio.h>
  33 #include "physics.h"
  34 #include "smalloc.h"
  35 #include "macros.h"
  36 #include "txtdump.h"
  37 #include "bondf.h"
  38 #include "xvgr.h"
  39 #include "typedefs.h"
  40 #include "gstat.h"
  41 #include "lutab.h"
  42
  43 static int calc_RBbin(real phi)
  44 {
  45   static const real r30  = M_PI/6.0;
  46   static const real r90  = M_PI/2.0;
  47   static const real r150 = M_PI*5.0/6.0;
  48
  49   if ((phi < r30) && (phi > -r30))
  50     return 1;
  51   else if ((phi > -r150) && (phi < -r90))
  52     return 2;
  53   else if ((phi < r150) && (phi > r90))
  54     return 3;
  55   return 0;
  56 }
  57
  58 static int calc_Nbin(real phi)
  59 {
  60   static const real r30  =  30*DEG2RAD;
  61   static const real r90  =  90*DEG2RAD;
  62   static const real r150 = 150*DEG2RAD;
  63   static const real r210 = 210*DEG2RAD;
  64   static const real r270 = 270*DEG2RAD;
  65   static const real r330 = 330*DEG2RAD;
  66   static const real r360 = 360*DEG2RAD;
  67
  68   if (phi < 0)
  69     phi += r360;
  70
  71   if ((phi > r30) && (phi < r90))
  72     return 1;
  73   else if ((phi > r150) && (phi < r210))
  74     return 2;
  75   else if ((phi > r270) && (phi < r330))
  76     return 3;
  77   return 0;
  78 }
  79
  80 void ana_dih_trans(char *fn_trans,char *fn_histo,
  81                    real **dih,int nframes,int nangles,
  82                    char *grpname,real t0,real dt,bool bRb)
  83 {
  84   FILE *fp;
  85   int  *tr_f,*tr_h;
  86   char title[256];
  87   int  i,j,ntrans;
  88   int  cur_bin,new_bin;
  89   real ttime,tt,mind, maxd, prev;
  90   int  (*calc_bin)(real);
  91
  92   /* Analysis of dihedral transitions */
  93   fprintf(stderr,"Now calculating transitions...\n");
  94
  95   if (bRb)
  96     calc_bin=calc_RBbin;
  97   else
  98     calc_bin=calc_Nbin;
  99
 100   snew(tr_h,nangles);
 101   snew(tr_f,nframes);
 102
 103   /* dih[i][j] is the dihedral angle i in frame j  */
 104   ntrans = 0;
 105   for (i=0; (i<nangles); i++)
 106   {
 107     /*#define OLDIE*/
 108 #ifdef OLDIE
 109     mind = maxd = prev = dih[i][0];
 110 #else
 111     cur_bin = calc_bin(dih[i][0]);
 112 #endif
 113     for (j=1; (j<nframes); j++)
 114     {
 115       new_bin = calc_bin(dih[i][j]);
 116 #ifndef OLDIE
 117       if (cur_bin == 0)
 118         cur_bin=new_bin;
 119       else if ((new_bin != 0) && (cur_bin != new_bin)) {
 120         cur_bin = new_bin;
 121         tr_f[j]++;
 122         tr_h[i]++;
 123         ntrans++;
 124       }
 125 #else
 126       /* why is all this md rubbish periodic? Remove 360 degree periodicity */
 127       if ( (dih[i][j] - prev) > M_PI)
 128         dih[i][j] -= 2*M_PI;
 129       else if ( (dih[i][j] - prev) < -M_PI)
 130         dih[i][j] += 2*M_PI;
 131
 132       prev = dih[i][j];
 133
 134       mind = min(mind, dih[i][j]);
 135       maxd = max(maxd, dih[i][j]);
 136       if ( (maxd - mind) > 2*M_PI/3)    /* or 120 degrees, assuming       */
 137       {                                 /* multiplicity 3. Not so general.*/
 138         tr_f[j]++;
 139         tr_h[i]++;
 140         maxd = mind = dih[i][j];        /* get ready for next transition  */
 141         ntrans++;
 142       }
 143 #endif
 144     }
 145   }
 146   fprintf(stderr,"Total number of transitions: %10d\n",ntrans);
 147   if (ntrans > 0) {
 148     ttime = (dt*nframes*nangles)/ntrans;
 149     fprintf(stderr,"Time between transitions:    %10.3f ps\n",ttime);
 150   }
 151
 152   sprintf(title,"Number of transitions: %s",grpname);
 153   fp=xvgropen(fn_trans,title,"Time (ps)","# transitions/timeframe");
 154   for(j=0; (j<nframes); j++) {
 155     tt = t0+j*dt;
 156     fprintf(fp,"%10.3f  %10d\n",tt,tr_f[j]);
 157   }
 158   ffclose(fp);
 159
 160   /* Compute histogram from # transitions per dihedral */
 161   /* Use old array */
 162   for(j=0; (j<nframes); j++)
 163     tr_f[j]=0;
 164   for(i=0; (i<nangles); i++)
 165     tr_f[tr_h[i]]++;
 166   for(j=nframes; ((tr_f[j-1] == 0) && (j>0)); j--)
 167     ;
 168
 169   ttime = dt*nframes;
 170   sprintf(title,"Transition time: %s",grpname);
 171   fp=xvgropen(fn_histo,title,"Time (ps)","#");
 172   for(i=j-1; (i>0); i--) {
 173     if (tr_f[i] != 0)
 174       fprintf(fp,"%10.3f  %10d\n",ttime/i,tr_f[i]);
 175   }
 176   ffclose(fp);
 177
 178   sfree(tr_f);
 179   sfree(tr_h);
 180 }
 181
 182 void calc_distribution_props(int nh,int histo[],real start,
 183                              int nkkk, t_karplus kkk[],
 184                              real *S2)
 185 {
 186   real d,dc,ds,c1,c2,tdc,tds;
 187   real fac,ang,invth;
 188   int  i,j,th;
 189
 190   if (nh == 0)
 191     fatal_error(0,"No points in histogram (%s, %d)",__FILE__,__LINE__);
 192   fac = 2*M_PI/nh;
 193
 194   /* Compute normalisation factor */
 195   th=0;
 196   for(j=0; (j<nh); j++)
 197     th+=histo[j];
 198   invth=1.0/th;
 199
 200   for(i=0; (i<nkkk); i++)
 201     kkk[i].Jc = 0;
 202
 203   tdc=0,tds=0;
 204   for(j=0; (j<nh); j++) {
 205     d    = invth*histo[j];
 206     ang  = j*fac-start;
 207     c1   = cos(ang);
 208     c2   = c1*c1;
 209     dc   = d*c1;
 210     ds   = d*sin(ang);
 211     tdc += dc;
 212     tds += ds;
 213     for(i=0; (i<nkkk); i++) {
 214       c1   = cos(ang+kkk[i].offset);
 215       c2   = c1*c1;
 216       kkk[i].Jc += d*(kkk[i].A*c2 + kkk[i].B*c1 + kkk[i].C);
 217     }
 218   }
 219   *S2 = tdc*tdc+tds*tds;
 220 }
 221
 222 static void calc_angles(FILE *log,matrix box,
 223                         int n3,atom_id index[],real ang[],rvec x_s[])
 224 {
 225   int  i,ix;
 226   rvec r_ij,r_kj;
 227   real costh;
 228
 229   for(i=ix=0; (ix<n3); i++,ix+=3)
 230     ang[i]=bond_angle(box,x_s[index[ix]],x_s[index[ix+1]],x_s[index[ix+2]],
 231                       r_ij,r_kj,&costh);
 232   if (debug) {
 233     fprintf(debug,"Angle[0]=%g, costh=%g, index0 = %u, %u, %u\n",
 234             ang[0],costh,index[0],index[1],index[2]);
 235     pr_rvec(debug,0,"rij",r_ij,DIM);
 236     pr_rvec(debug,0,"rkj",r_kj,DIM);
 237     pr_rvecs(debug,0,"box",box,DIM);
 238   }
 239 }
 240
 241 static real calc_fraction(real angles[], int nangles)
 242 {
 243   int i;
 244   real trans = 0, gauche = 0;
 245   real angle;
 246
 247   for (i = 0; i < nangles; i++)
 248   {
 249     angle = angles[i] * RAD2DEG;
 250
 251     if (angle > 135 && angle < 225)
 252       trans += 1.0;
 253     else if (angle > 270 && angle < 330)
 254       gauche += 1.0;
 255     else if (angle < 90 && angle > 30)
 256       gauche += 1.0;
 257   }
 258   if (trans+gauche > 0)
 259     return trans/(trans+gauche);
 260   else
 261     return 0;
 262 }
 263
 264 static void calc_dihs(FILE *log,matrix box,
 265                       int n4,atom_id index[],real ang[],rvec x_s[])
 266 {
 267   int  i,ix;
 268   rvec r_ij,r_kj,r_kl,m,n;
 269   real cos_phi,sign,aaa;
 270
 271   for(i=ix=0; (ix<n4); i++,ix+=4) {
 272     aaa=dih_angle(box,
 273                   x_s[index[ix]],x_s[index[ix+1]],x_s[index[ix+2]],
 274                   x_s[index[ix+3]],
 275                   r_ij,r_kj,r_kl,m,n,
 276                   &cos_phi,&sign);
 277     ang[i]=aaa;  /* not taking into account ryckaert bellemans yet */
 278   }
 279 }
 280
 281 void make_histo(FILE *log,
 282                 int ndata,real data[],int npoints,int histo[],
 283                 real minx,real maxx)
 284 {
 285   double dx;
 286   int    i,ind;
 287
 288   if (minx == maxx) {
 289     minx=maxx=data[0];
 290     for(i=1; (i<ndata); i++) {
 291       minx=min(minx,data[i]);
 292       maxx=max(maxx,data[i]);
 293     }
 294     fprintf(log,"Min data: %10g  Max data: %10g\n",minx,maxx);
 295   }
 296   dx=(double)npoints/(maxx-minx);
 297   if (debug)
 298     fprintf(debug,
 299             "Histogramming: ndata=%d, nhisto=%d, minx=%g,maxx=%g,dx=%g\n",
 300             ndata,npoints,minx,maxx,dx);
 301   for(i=0; (i<ndata); i++) {
 302     ind=(data[i]-minx)*dx;
 303     if ((ind >= 0) && (ind < npoints))
 304       histo[ind]++;
 305     else
 306       fprintf(log,"index = %d, data[%d] = %g\n",ind,i,data[i]);
 307   }
 308 }
 309
 310 void normalize_histo(int npoints,int histo[],real dx,real normhisto[])
 311 {
 312   int    i;
 313   double d,fac;
 314
 315   d=0;
 316   for(i=0; (i<npoints); i++)
 317     d+=dx*histo[i];
 318   if (d==0) {
 319     fprintf(stderr,"Empty histogram!\n");
 320     return;
 321   }
 322   fac=1.0/d;
 323   for(i=0; (i<npoints); i++)
 324     normhisto[i]=fac*histo[i];
 325 }
 326
 327 void read_ang_dih(char *trj_fn,char *tpb_fn,
 328                   bool bAngles,bool bSaveAll,bool bRb,
 329                   int maxangstat,int angstat[],
 330                   int *nframes,real **time,
 331                   int isize,atom_id index[],
 332                   real **trans_frac,
 333                   real **aver_angle,
 334                   real *dih[])
 335 {
 336   t_topology *top;
 337   int        i,angind,status,natoms,nat,total,teller,nangles,nat_trj,n_alloc;
 338   real       t,fraction,pifac,aa,angle;
 339   real       *angles[2];
 340   matrix     box;
 341   rvec       *x,*x_s;
 342   int        cur=0;
 343 #define prev (1-cur)
 344
 345   /* Initiate lookup table for sqrt calculations */
 346   init_lookup_table(stdout);
 347
 348   /* Read topology */
 349   top     = read_top(tpb_fn);
 350   natoms  = top->atoms.nr;
 351   nat_trj = read_first_x(&status,trj_fn,&t,&x,box);
 352
 353   /* Check for consistency of topology and trajectory */
 354   if (natoms < nat_trj)
 355     fprintf(stderr,"WARNING! Topology has fewer atoms than trajectory\n");
 356
 357   snew(x_s,nat_trj);
 358
 359   if (bAngles) {
 360     nangles=isize/3;
 361     pifac=M_PI;
 362   }
 363   else {
 364     nangles=isize/4;
 365     pifac=2.0*M_PI;
 366   }
 367   snew(angles[cur],nangles);
 368   snew(angles[prev],nangles);
 369
 370   /* Start the loop over frames */
 371   total       = 0;
 372   teller      = 0;
 373   n_alloc     = 0;
 374   *time       = NULL;
 375   *trans_frac = NULL;
 376   *aver_angle = NULL;
 377
 378   do {
 379     if (teller >= n_alloc) {
 380       n_alloc+=100;
 381       if (bSaveAll)
 382         for (i=0; (i<nangles); i++)
 383           srenew(dih[i],n_alloc);
 384       srenew(*time,n_alloc);
 385       srenew(*trans_frac,n_alloc);
 386       srenew(*aver_angle,n_alloc);
 387     }
 388
 389     (*time)[teller] = t;
 390
 391     rm_pbc(&(top->idef),nat_trj,box,x,x_s);
 392
 393     if (bAngles)
 394       calc_angles(stdout,box,isize,index,angles[cur],x_s);
 395     else {
 396       calc_dihs(stdout,box,isize,index,angles[cur],x_s);
 397
 398       /* Trans fraction */
 399       fraction = calc_fraction(angles[cur], nangles);
 400       (*trans_frac)[teller] = fraction;
 401
 402       /* Change Ryckaert-Bellemans dihedrals to polymer convention
 403        * Modified 990913 by Erik:
 404        * We actually shouldn't change the convention, since it's
 405        * calculated from polymer above, but we change the intervall
 406        * from [-180,180] to [0,360].
 407        */
 408          if (bRb) {
 409         for(i=0; (i<nangles); i++)
 410           if (angles[cur][i] <= 0.0)
 411             angles[cur][i] += 2*M_PI;
 412         }
 413
 414       /* Periodicity in dihedral space... */
 415       if (teller > 1) {
 416         for(i=0; (i<nangles); i++) {
 417           while (angles[cur][i] <= angles[prev][i] - M_PI)
 418             angles[cur][i]+=2*M_PI;
 419           while (angles[cur][i] > angles[prev][i] + M_PI)
 420             angles[cur][i]-=2*M_PI;
 421         }
 422       }
 423     }
 424
 425     /* Average angles */
 426     aa=0;
 427     for(i=0; (i<nangles); i++) {
 428       aa=aa+angles[cur][i];
 429
 430       /* angle in rad / 2Pi * max determines bin. bins go from 0 to maxangstat,
 431          even though scale goes from -pi to pi (dihedral) or -pi/2 to pi/2
 432          (angle) Basically: translate the x-axis by Pi. Translate it back by
 433          -Pi when plotting.
 434        */
 435
 436       angle = angles[cur][i];
 437       if (!bAngles) {
 438         while (angle < -M_PI)
 439           angle += 2*M_PI;
 440         while (angle >= M_PI)
 441           angle -= 2*M_PI;
 442
 443         angle+=M_PI;
 444       }
 445
 446       /* Update the distribution histogram */
 447       angind = (int) ((angle*maxangstat)/pifac + 0.5);
 448       if (angind==maxangstat)
 449         angind=0;
 450       if ( (angind < 0) || (angind >= maxangstat) )
 451         /* this will never happen */
 452         fatal_error(0,"angle (%f) index out of range (0..%d) : %d\n",
 453                     angle,maxangstat,angind);
 454
 455       angstat[angind]++;
 456       if (angind==maxangstat)
 457         fprintf(stderr,"angle %d fr %d = %g\n",i,cur,angle);
 458
 459       total++;
 460     }
 461
 462     /* average over all angles */
 463     (*aver_angle)[teller] = (aa/nangles);
 464
 465     /* this copies all current dih. angles to dih[i], teller is frame */
 466     if (bSaveAll)
 467       for (i = 0; i < nangles; i++)
 468         dih[i][teller] = angles[cur][i];
 469
 470     /* Swap buffers */
 471     cur=prev;
 472
 473     /* Increment loop counter */
 474     teller++;
 475   } while (read_next_x(status,&t,nat_trj,x,box));
 476   close_trj(status);
 477
 478   sfree(x);
 479   sfree(x_s);
 480   sfree(angles[cur]);
 481   sfree(angles[prev]);
 482
 483   *nframes=teller;
 484 }
 485