src/gromacs/gmxana/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     {
  67         fprintf(fp, "%10g  %10g\n", time[k], data[k]);
  68     }
  69     ffclose(fp);
  70 }
  71
  72 static int calc_RBbin(real phi, int gmx_unused multiplicity, real gmx_unused core_frac)
  73 {
  74     /* multiplicity and core_frac NOT used,
  75      * just given to enable use of pt-to-fn in caller low_ana_dih_trans*/
  76     static const real r30  = M_PI/6.0;
  77     static const real r90  = M_PI/2.0;
  78     static const real r150 = M_PI*5.0/6.0;
  79
  80     if ((phi < r30) && (phi > -r30))
  81     {
  82         return 1;
  83     }
  84     else if ((phi > -r150) && (phi < -r90))
  85     {
  86         return 2;
  87     }
  88     else if ((phi < r150) && (phi > r90))
  89     {
  90         return 3;
  91     }
  92     return 0;
  93 }
  94
  95 static int calc_Nbin(real phi, int multiplicity, real core_frac)
  96 {
  97     static const real r360 = 360*DEG2RAD;
  98     real              rot_width, core_width, core_offset, low, hi;
  99     int               bin;
 100     /* with multiplicity 3 and core_frac 0.5
 101      * 0<g(-)<120, 120<t<240, 240<g(+)<360
 102      * 0< bin0 < 30, 30<bin1<90, 90<bin0<150, 150<bin2<210, 210<bin0<270, 270<bin3<330, 330<bin0<360
 103      * so with multiplicity 3, bin1 is core g(-), bin2 is core t, bin3 is
 104        core g(+), bin0 is between rotamers */
 105     if (phi < 0)
 106     {
 107         phi += r360;
 108     }
 109
 110     rot_width   = 360/multiplicity;
 111     core_width  = core_frac * rot_width;
 112     core_offset = (rot_width - core_width)/2.0;
 113     for (bin = 1; bin <= multiplicity; bin++)
 114     {
 115         low  = ((bin - 1) * rot_width ) + core_offset;
 116         hi   = ((bin - 1) * rot_width ) + core_offset + core_width;
 117         low *= DEG2RAD;
 118         hi  *= DEG2RAD;
 119         if ((phi > low) && (phi < hi))
 120         {
 121             return bin;
 122         }
 123     }
 124     return 0;
 125 }
 126
 127 void ana_dih_trans(const char *fn_trans, const char *fn_histo,
 128                    real **dih, int nframes, int nangles,
 129                    const char *grpname, real *time, gmx_bool bRb,
 130                    const output_env_t oenv)
 131 {
 132     /* just a wrapper; declare extra args, then chuck away at end. */
 133     int      maxchi = 0;
 134     t_dlist *dlist;
 135     int     *multiplicity;
 136     int      nlist = nangles;
 137     int      k;
 138
 139     snew(dlist, nlist);
 140     snew(multiplicity, nangles);
 141     for (k = 0; (k < nangles); k++)
 142     {
 143         multiplicity[k] = 3;
 144     }
 145
 146     low_ana_dih_trans(TRUE, fn_trans, TRUE, fn_histo, maxchi,
 147                       dih, nlist, dlist, nframes,
 148                       nangles, grpname, multiplicity, time, bRb, 0.5, oenv);
 149     sfree(dlist);
 150     sfree(multiplicity);
 151
 152 }
 153
 154 void low_ana_dih_trans(gmx_bool bTrans, const char *fn_trans,
 155                        gmx_bool bHisto, const char *fn_histo, int maxchi,
 156                        real **dih, int nlist, t_dlist dlist[], int nframes,
 157                        int nangles, const char *grpname, int multiplicity[],
 158                        real *time, gmx_bool bRb, real core_frac,
 159                        const output_env_t oenv)
 160 {
 161     FILE *fp;
 162     int  *tr_f, *tr_h;
 163     char  title[256];
 164     int   i, j, k, Dih, ntrans;
 165     int   cur_bin, new_bin;
 166     real  ttime, tt;
 167     real *rot_occ[NROT];
 168     int   (*calc_bin)(real, int, real);
 169     real  dt;
 170
 171     if (1 <= nframes)
 172     {
 173         return;
 174     }
 175     /* Assumes the frames are equally spaced in time */
 176     dt = (time[nframes-1]-time[0])/(nframes-1);
 177
 178     /* Analysis of dihedral transitions */
 179     fprintf(stderr, "Now calculating transitions...\n");
 180
 181     if (bRb)
 182     {
 183         calc_bin = calc_RBbin;
 184     }
 185     else
 186     {
 187         calc_bin = calc_Nbin;
 188     }
 189
 190     for (k = 0; k < NROT; k++)
 191     {
 192         snew(rot_occ[k], nangles);
 193         for (i = 0; (i < nangles); i++)
 194         {
 195             rot_occ[k][i] = 0;
 196         }
 197     }
 198     snew(tr_h, nangles);
 199     snew(tr_f, nframes);
 200
 201     /* dih[i][j] is the dihedral angle i in frame j  */
 202     ntrans = 0;
 203     for (i = 0; (i < nangles); i++)
 204     {
 205
 206         /*#define OLDIE*/
 207 #ifdef OLDIE
 208         mind = maxd = prev = dih[i][0];
 209 #else
 210         cur_bin = calc_bin(dih[i][0], multiplicity[i], core_frac);
 211         rot_occ[cur_bin][i]++;
 212 #endif
 213         for (j = 1; (j < nframes); j++)
 214         {
 215             new_bin = calc_bin(dih[i][j], multiplicity[i], core_frac);
 216             rot_occ[new_bin][i]++;
 217 #ifndef OLDIE
 218             if (cur_bin == 0)
 219             {
 220                 cur_bin = new_bin;
 221             }
 222             else if ((new_bin != 0) && (cur_bin != new_bin))
 223             {
 224                 cur_bin = new_bin;
 225                 tr_f[j]++;
 226                 tr_h[i]++;
 227                 ntrans++;
 228             }
 229 #else
 230             /* why is all this md rubbish periodic? Remove 360 degree periodicity */
 231             if ( (dih[i][j] - prev) > M_PI)
 232             {
 233                 dih[i][j] -= 2*M_PI;
 234             }
 235             else if ( (dih[i][j] - prev) < -M_PI)
 236             {
 237                 dih[i][j] += 2*M_PI;
 238             }
 239
 240             prev = dih[i][j];
 241
 242             mind = min(mind, dih[i][j]);
 243             maxd = max(maxd, dih[i][j]);
 244             if ( (maxd - mind) > 2*M_PI/3) /* or 120 degrees, assuming       */
 245             {                              /* multiplicity 3. Not so general.*/
 246                 tr_f[j]++;
 247                 tr_h[i]++;
 248                 maxd = mind = dih[i][j]; /* get ready for next transition  */
 249                 ntrans++;
 250             }
 251 #endif
 252         } /* end j */
 253         for (k = 0; k < NROT; k++)
 254         {
 255             rot_occ[k][i] /= nframes;
 256         }
 257     } /* end i */
 258     fprintf(stderr, "Total number of transitions: %10d\n", ntrans);
 259     if (ntrans > 0)
 260     {
 261         ttime = (dt*nframes*nangles)/ntrans;
 262         fprintf(stderr, "Time between transitions:    %10.3f ps\n", ttime);
 263     }
 264
 265     /* new by grs - copy transitions from tr_h[] to dlist->ntr[]
 266      * and rotamer populations from rot_occ to dlist->rot_occ[]
 267      * based on fn histogramming in g_chi. diff roles for i and j here */
 268
 269     j = 0;
 270     for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
 271     {
 272         for (i = 0; (i < nlist); i++)
 273         {
 274             if (((Dih  < edOmega) ) ||
 275                 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
 276                 ((Dih  > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)))
 277             {
 278                 /* grs debug  printf("Not OK? i %d j %d Dih %d \n", i, j, Dih) ; */
 279                 dlist[i].ntr[Dih] = tr_h[j];
 280                 for (k = 0; k < NROT; k++)
 281                 {
 282                     dlist[i].rot_occ[Dih][k] = rot_occ[k][j];
 283                 }
 284                 j++;
 285             }
 286         }
 287     }
 288
 289     /* end addition by grs */
 290
 291     if (bTrans)
 292     {
 293         sprintf(title, "Number of transitions: %s", grpname);
 294         fp = xvgropen(fn_trans, title, "Time (ps)", "# transitions/timeframe", oenv);
 295         for (j = 0; (j < nframes); j++)
 296         {
 297             fprintf(fp, "%10.3f  %10d\n", time[j], tr_f[j]);
 298         }
 299         ffclose(fp);
 300     }
 301
 302     /* Compute histogram from # transitions per dihedral */
 303     /* Use old array */
 304     for (j = 0; (j < nframes); j++)
 305     {
 306         tr_f[j] = 0;
 307     }
 308     for (i = 0; (i < nangles); i++)
 309     {
 310         tr_f[tr_h[i]]++;
 311     }
 312     for (j = nframes; ((tr_f[j-1] == 0) && (j > 0)); j--)
 313     {
 314         ;
 315     }
 316
 317     ttime = dt*nframes;
 318     if (bHisto)
 319     {
 320         sprintf(title, "Transition time: %s", grpname);
 321         fp = xvgropen(fn_histo, title, "Time (ps)", "#", oenv);
 322         for (i = j-1; (i > 0); i--)
 323         {
 324             if (tr_f[i] != 0)
 325             {
 326                 fprintf(fp, "%10.3f  %10d\n", ttime/i, tr_f[i]);
 327             }
 328         }
 329         ffclose(fp);
 330     }
 331
 332     sfree(tr_f);
 333     sfree(tr_h);
 334     for (k = 0; k < NROT; k++)
 335     {
 336         sfree(rot_occ[k]);
 337     }
 338
 339 }
 340
 341 void mk_multiplicity_lookup (int *multiplicity, int maxchi,
 342                              int nlist, t_dlist dlist[], int nangles)
 343 {
 344     /* new by grs - for dihedral j (as in dih[j]) get multiplicity from dlist
 345      * and store in multiplicity[j]
 346      */
 347
 348     int  j, Dih, i;
 349     char name[4];
 350
 351     j = 0;
 352     for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
 353     {
 354         for (i = 0; (i < nlist); i++)
 355         {
 356             strncpy(name, dlist[i].name, 3);
 357             name[3] = '\0';
 358             if (((Dih  < edOmega) ) ||
 359                 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
 360                 ((Dih  > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)))
 361             {
 362                 /* default - we will correct the rest below */
 363                 multiplicity[j] = 3;
 364
 365                 /* make omegas 2fold, though doesn't make much more sense than 3 */
 366                 if (Dih == edOmega && (has_dihedral(edOmega, &(dlist[i]))))
 367                 {
 368                     multiplicity[j] = 2;
 369                 }
 370
 371                 /* dihedrals to aromatic rings, COO, CONH2 or guanidinium are 2fold*/
 372                 if (Dih > edOmega && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1))
 373                 {
 374                     if ( ((strstr(name, "PHE") != NULL) && (Dih == edChi2))  ||
 375                          ((strstr(name, "TYR") != NULL) && (Dih == edChi2))  ||
 376                          ((strstr(name, "PTR") != NULL) && (Dih == edChi2))  ||
 377                          ((strstr(name, "TRP") != NULL) && (Dih == edChi2))  ||
 378                          ((strstr(name, "HIS") != NULL) && (Dih == edChi2))  ||
 379                          ((strstr(name, "GLU") != NULL) && (Dih == edChi3))  ||
 380                          ((strstr(name, "ASP") != NULL) && (Dih == edChi2))  ||
 381                          ((strstr(name, "GLN") != NULL) && (Dih == edChi3))  ||
 382                          ((strstr(name, "ASN") != NULL) && (Dih == edChi2))  ||
 383                          ((strstr(name, "ARG") != NULL) && (Dih == edChi4))  )
 384                     {
 385                         multiplicity[j] = 2;
 386                     }
 387                 }
 388                 j++;
 389             }
 390         }
 391     }
 392     if (j < nangles)
 393     {
 394         fprintf(stderr, "WARNING: not all dihedrals found in topology (only %d out of %d)!\n",
 395                 j, nangles);
 396     }
 397     /* Check for remaining dihedrals */
 398     for (; (j < nangles); j++)
 399     {
 400         multiplicity[j] = 3;
 401     }
 402
 403 }
 404
 405 void mk_chi_lookup (int **lookup, int maxchi,
 406                     int nlist, t_dlist dlist[])
 407 {
 408
 409     /* by grs. should rewrite everything to use this. (but haven't,
 410      * and at mmt only used in get_chi_product_traj
 411      * returns the dihed number given the residue number (from-0)
 412      * and chi (from-0) nr. -1 for chi undefined for that res (eg gly, ala..)*/
 413
 414     int i, j, Dih, Chi;
 415
 416     j = 0;
 417     for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
 418     {
 419         for (i = 0; (i < nlist); i++)
 420         {
 421             Chi = Dih - NONCHI;
 422             if (((Dih  < edOmega) ) ||
 423                 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
 424                 ((Dih  > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)))
 425             {
 426                 /* grs debug  printf("Not OK? i %d j %d Dih %d \n", i, j, Dih) ; */
 427                 if (Dih > edOmega)
 428                 {
 429                     lookup[i][Chi] = j;
 430                 }
 431                 j++;
 432             }
 433             else
 434             {
 435                 lookup[i][Chi] = -1;
 436             }
 437         }
 438     }
 439
 440 }
 441
 442
 443 void get_chi_product_traj (real **dih, int nframes, int nlist,
 444                            int maxchi, t_dlist dlist[], real time[],
 445                            int **lookup, int *multiplicity, gmx_bool bRb, gmx_bool bNormalize,
 446                            real core_frac, gmx_bool bAll, const char *fnall,
 447                            const output_env_t oenv)
 448 {
 449
 450     gmx_bool bRotZero, bHaveChi = FALSE;
 451     int      accum = 0, index, i, j, k, Xi, n, b;
 452     real    *chi_prtrj;
 453     int     *chi_prhist;
 454     int      nbin;
 455     FILE    *fp, *fpall;
 456     char     hisfile[256], histitle[256], *namept;
 457
 458     int      (*calc_bin)(real, int, real);
 459
 460     /* Analysis of dihedral transitions */
 461     fprintf(stderr, "Now calculating Chi product trajectories...\n");
 462
 463     if (bRb)
 464     {
 465         calc_bin = calc_RBbin;
 466     }
 467     else
 468     {
 469         calc_bin = calc_Nbin;
 470     }
 471
 472     snew(chi_prtrj, nframes);
 473
 474     /* file for info on all residues */
 475     if (bNormalize)
 476     {
 477         fpall = xvgropen(fnall, "Cumulative Rotamers", "Residue", "Probability", oenv);
 478     }
 479     else
 480     {
 481         fpall = xvgropen(fnall, "Cumulative Rotamers", "Residue", "# Counts", oenv);
 482     }
 483
 484     for (i = 0; (i < nlist); i++)
 485     {
 486
 487         /* get nbin, the nr. of cumulative rotamers that need to be considered */
 488         nbin = 1;
 489         for (Xi = 0; Xi < maxchi; Xi++)
 490         {
 491             index = lookup[i][Xi]; /* chi_(Xi+1) of res i (-1 if off end) */
 492             if (index >= 0)
 493             {
 494                 n    = multiplicity[index];
 495                 nbin = n*nbin;
 496             }
 497         }
 498         nbin += 1; /* for the "zero rotamer", outside the core region */
 499
 500         for (j = 0; (j < nframes); j++)
 501         {
 502
 503             bRotZero = FALSE;
 504             bHaveChi = TRUE;
 505             index    = lookup[i][0]; /* index into dih of chi1 of res i */
 506             if (index == -1)
 507             {
 508                 b        = 0;
 509                 bRotZero = TRUE;
 510                 bHaveChi = FALSE;
 511             }
 512             else
 513             {
 514                 b     = calc_bin(dih[index][j], multiplicity[index], core_frac);
 515                 accum = b - 1;
 516                 if (b == 0)
 517                 {
 518                     bRotZero = TRUE;
 519                 }
 520                 for (Xi = 1; Xi < maxchi; Xi++)
 521                 {
 522                     index = lookup[i][Xi]; /* chi_(Xi+1) of res i (-1 if off end) */
 523                     if (index >= 0)
 524                     {
 525                         n     = multiplicity[index];
 526                         b     = calc_bin(dih[index][j], n, core_frac);
 527                         accum = n * accum + b - 1;
 528                         if (b == 0)
 529                         {
 530                             bRotZero = TRUE;
 531                         }
 532                     }
 533                 }
 534                 accum++;
 535             }
 536             if (bRotZero)
 537             {
 538                 chi_prtrj[j] = 0.0;
 539             }
 540             else
 541             {
 542                 chi_prtrj[j] = accum;
 543                 if (accum+1 > nbin)
 544                 {
 545                     nbin = accum+1;
 546                 }
 547             }
 548         }
 549         if (bHaveChi)
 550         {
 551
 552             if (bAll)
 553             {
 554                 /* print cuml rotamer vs time */
 555                 print_one(oenv, "chiproduct", dlist[i].name, "chi product for",
 556                           "cumulative rotamer", nframes, time, chi_prtrj);
 557             }
 558
 559             /* make a histogram pf culm. rotamer occupancy too */
 560             snew(chi_prhist, nbin);
 561             make_histo(NULL, nframes, chi_prtrj, nbin, chi_prhist, 0, nbin);
 562             if (bAll)
 563             {
 564                 sprintf(hisfile, "histo-chiprod%s.xvg", dlist[i].name);
 565                 sprintf(histitle, "cumulative rotamer distribution for %s", dlist[i].name);
 566                 fprintf(stderr, "  and %s  ", hisfile);
 567                 fp = xvgropen(hisfile, histitle, "number", "", oenv);
 568                 fprintf(fp, "@ xaxis tick on\n");
 569                 fprintf(fp, "@ xaxis tick major 1\n");
 570                 fprintf(fp, "@ type xy\n");
 571                 for (k = 0; (k < nbin); k++)
 572                 {
 573                     if (bNormalize)
 574                     {
 575                         fprintf(fp, "%5d  %10g\n", k, (1.0*chi_prhist[k])/nframes);
 576                     }
 577                     else
 578                     {
 579                         fprintf(fp, "%5d  %10d\n", k, chi_prhist[k]);
 580                     }
 581                 }
 582                 fprintf(fp, "&\n");
 583                 ffclose(fp);
 584             }
 585
 586             /* and finally print out occupancies to a single file */
 587             /* get the gmx from-1 res nr by setting a ptr to the number part
 588              * of dlist[i].name - potential bug for 4-letter res names... */
 589             namept = dlist[i].name + 3;
 590             fprintf(fpall, "%5s ", namept);
 591             for (k = 0; (k < nbin); k++)
 592             {
 593                 if (bNormalize)
 594                 {
 595                     fprintf(fpall, "  %10g", (1.0*chi_prhist[k])/nframes);
 596                 }
 597                 else
 598                 {
 599                     fprintf(fpall, "  %10d", chi_prhist[k]);
 600                 }
 601             }
 602             fprintf(fpall, "\n");
 603
 604             sfree(chi_prhist);
 605             /* histogram done */
 606         }
 607     }
 608
 609     sfree(chi_prtrj);
 610     ffclose(fpall);
 611     fprintf(stderr, "\n");
 612
 613 }
 614
 615 void calc_distribution_props(int nh, int histo[], real start,
 616                              int nkkk, t_karplus kkk[],
 617                              real *S2)
 618 {
 619     real d, dc, ds, c1, c2, tdc, tds;
 620     real fac, ang, invth, Jc;
 621     int  i, j, th;
 622
 623     if (nh == 0)
 624     {
 625         gmx_fatal(FARGS, "No points in histogram (%s, %d)", __FILE__, __LINE__);
 626     }
 627     fac = 2*M_PI/nh;
 628
 629     /* Compute normalisation factor */
 630     th = 0;
 631     for (j = 0; (j < nh); j++)
 632     {
 633         th += histo[j];
 634     }
 635     invth = 1.0/th;
 636
 637     for (i = 0; (i < nkkk); i++)
 638     {
 639         kkk[i].Jc    = 0;
 640         kkk[i].Jcsig = 0;
 641     }
 642     tdc = 0, tds = 0;
 643     for (j = 0; (j < nh); j++)
 644     {
 645         d    = invth*histo[j];
 646         ang  = j*fac-start;
 647         c1   = cos(ang);
 648         c2   = c1*c1;
 649         dc   = d*c1;
 650         ds   = d*sin(ang);
 651         tdc += dc;
 652         tds += ds;
 653         for (i = 0; (i < nkkk); i++)
 654         {
 655             c1            = cos(ang+kkk[i].offset);
 656             c2            = c1*c1;
 657             Jc            = (kkk[i].A*c2 + kkk[i].B*c1 + kkk[i].C);
 658             kkk[i].Jc    += histo[j]*Jc;
 659             kkk[i].Jcsig += histo[j]*sqr(Jc);
 660         }
 661     }
 662     for (i = 0; (i < nkkk); i++)
 663     {
 664         kkk[i].Jc    /= th;
 665         kkk[i].Jcsig  = sqrt(kkk[i].Jcsig/th-sqr(kkk[i].Jc));
 666     }
 667     *S2 = tdc*tdc+tds*tds;
 668 }
 669
 670 static void calc_angles(t_pbc *pbc,
 671                         int n3, atom_id index[], real ang[], rvec x_s[])
 672 {
 673     int  i, ix, t1, t2;
 674     rvec r_ij, r_kj;
 675     real costh;
 676
 677     for (i = ix = 0; (ix < n3); i++, ix += 3)
 678     {
 679         ang[i] = bond_angle(x_s[index[ix]], x_s[index[ix+1]], x_s[index[ix+2]],
 680                             pbc, r_ij, r_kj, &costh, &t1, &t2);
 681     }
 682     if (debug)
 683     {
 684         fprintf(debug, "Angle[0]=%g, costh=%g, index0 = %d, %d, %d\n",
 685                 ang[0], costh, index[0], index[1], index[2]);
 686         pr_rvec(debug, 0, "rij", r_ij, DIM, TRUE);
 687         pr_rvec(debug, 0, "rkj", r_kj, DIM, TRUE);
 688     }
 689 }
 690
 691 static real calc_fraction(real angles[], int nangles)
 692 {
 693     int  i;
 694     real trans = 0, gauche = 0;
 695     real angle;
 696
 697     for (i = 0; i < nangles; i++)
 698     {
 699         angle = angles[i] * RAD2DEG;
 700
 701         if (angle > 135 && angle < 225)
 702         {
 703             trans += 1.0;
 704         }
 705         else if (angle > 270 && angle < 330)
 706         {
 707             gauche += 1.0;
 708         }
 709         else if (angle < 90 && angle > 30)
 710         {
 711             gauche += 1.0;
 712         }
 713     }
 714     if (trans+gauche > 0)
 715     {
 716         return trans/(trans+gauche);
 717     }
 718     else
 719     {
 720         return 0;
 721     }
 722 }
 723
 724 static void calc_dihs(t_pbc *pbc,
 725                       int n4, atom_id index[], real ang[], rvec x_s[])
 726 {
 727     int  i, ix, t1, t2, t3;
 728     rvec r_ij, r_kj, r_kl, m, n;
 729     real sign, aaa;
 730
 731     for (i = ix = 0; (ix < n4); i++, ix += 4)
 732     {
 733         aaa = dih_angle(x_s[index[ix]], x_s[index[ix+1]], x_s[index[ix+2]],
 734                         x_s[index[ix+3]], pbc,
 735                         r_ij, r_kj, r_kl, m, n,
 736                         &sign, &t1, &t2, &t3);
 737
 738         ang[i] = aaa; /* not taking into account ryckaert bellemans yet */
 739     }
 740 }
 741
 742 void make_histo(FILE *log,
 743                 int ndata, real data[], int npoints, int histo[],
 744                 real minx, real maxx)
 745 {
 746     double dx;
 747     int    i, ind;
 748
 749     if (minx == maxx)
 750     {
 751         minx = maxx = data[0];
 752         for (i = 1; (i < ndata); i++)
 753         {
 754             minx = min(minx, data[i]);
 755             maxx = max(maxx, data[i]);
 756         }
 757         fprintf(log, "Min data: %10g  Max data: %10g\n", minx, maxx);
 758     }
 759     dx = (double)npoints/(maxx-minx);
 760     if (debug)
 761     {
 762         fprintf(debug,
 763                 "Histogramming: ndata=%d, nhisto=%d, minx=%g,maxx=%g,dx=%g\n",
 764                 ndata, npoints, minx, maxx, dx);
 765     }
 766     for (i = 0; (i < ndata); i++)
 767     {
 768         ind = (data[i]-minx)*dx;
 769         if ((ind >= 0) && (ind < npoints))
 770         {
 771             histo[ind]++;
 772         }
 773         else
 774         {
 775             fprintf(log, "index = %d, data[%d] = %g\n", ind, i, data[i]);
 776         }
 777     }
 778 }
 779
 780 void normalize_histo(int npoints, int histo[], real dx, real normhisto[])
 781 {
 782     int    i;
 783     double d, fac;
 784
 785     d = 0;
 786     for (i = 0; (i < npoints); i++)
 787     {
 788         d += dx*histo[i];
 789     }
 790     if (d == 0)
 791     {
 792         fprintf(stderr, "Empty histogram!\n");
 793         return;
 794     }
 795     fac = 1.0/d;
 796     for (i = 0; (i < npoints); i++)
 797     {
 798         normhisto[i] = fac*histo[i];
 799     }
 800 }
 801
 802 void read_ang_dih(const char *trj_fn,
 803                   gmx_bool bAngles, gmx_bool bSaveAll, gmx_bool bRb, gmx_bool bPBC,
 804                   int maxangstat, int angstat[],
 805                   int *nframes, real **time,
 806                   int isize, atom_id index[],
 807                   real **trans_frac,
 808                   real **aver_angle,
 809                   real *dih[],
 810                   const output_env_t oenv)
 811 {
 812     t_pbc       *pbc;
 813     t_trxstatus *status;
 814     int          i, angind, natoms, total, teller;
 815     int          nangles, n_alloc;
 816     real         t, fraction, pifac, aa, angle;
 817     real        *angles[2];
 818     matrix       box;
 819     rvec        *x;
 820     int          cur = 0;
 821 #define prev (1-cur)
 822
 823     snew(pbc, 1);
 824     natoms = read_first_x(oenv, &status, trj_fn, &t, &x, box);
 825
 826     if (bAngles)
 827     {
 828         nangles = isize/3;
 829         pifac   = M_PI;
 830     }
 831     else
 832     {
 833         nangles = isize/4;
 834         pifac   = 2.0*M_PI;
 835     }
 836     snew(angles[cur], nangles);
 837     snew(angles[prev], nangles);
 838
 839     /* Start the loop over frames */
 840     total       = 0;
 841     teller      = 0;
 842     n_alloc     = 0;
 843     *time       = NULL;
 844     *trans_frac = NULL;
 845     *aver_angle = NULL;
 846
 847     do
 848     {
 849         if (teller >= n_alloc)
 850         {
 851             n_alloc += 100;
 852             if (bSaveAll)
 853             {
 854                 for (i = 0; (i < nangles); i++)
 855                 {
 856                     srenew(dih[i], n_alloc);
 857                 }
 858             }
 859             srenew(*time, n_alloc);
 860             srenew(*trans_frac, n_alloc);
 861             srenew(*aver_angle, n_alloc);
 862         }
 863
 864         (*time)[teller] = t;
 865
 866         if (pbc)
 867         {
 868             set_pbc(pbc, -1, box);
 869         }
 870
 871         if (bAngles)
 872         {
 873             calc_angles(pbc, isize, index, angles[cur], x);
 874         }
 875         else
 876         {
 877             calc_dihs(pbc, isize, index, angles[cur], x);
 878
 879             /* Trans fraction */
 880             fraction              = calc_fraction(angles[cur], nangles);
 881             (*trans_frac)[teller] = fraction;
 882
 883             /* Change Ryckaert-Bellemans dihedrals to polymer convention
 884              * Modified 990913 by Erik:
 885              * We actually shouldn't change the convention, since it's
 886              * calculated from polymer above, but we change the intervall
 887              * from [-180,180] to [0,360].
 888              */
 889             if (bRb)
 890             {
 891                 for (i = 0; (i < nangles); i++)
 892                 {
 893                     if (angles[cur][i] <= 0.0)
 894                     {
 895                         angles[cur][i] += 2*M_PI;
 896                     }
 897                 }
 898             }
 899
 900             /* Periodicity in dihedral space... */
 901             if (bPBC)
 902             {
 903                 for (i = 0; (i < nangles); i++)
 904                 {
 905                     real dd = angles[cur][i];
 906                     angles[cur][i] = atan2(sin(dd), cos(dd));
 907                 }
 908             }
 909             else
 910             {
 911                 if (teller > 1)
 912                 {
 913                     for (i = 0; (i < nangles); i++)
 914                     {
 915                         while (angles[cur][i] <= angles[prev][i] - M_PI)
 916                         {
 917                             angles[cur][i] += 2*M_PI;
 918                         }
 919                         while (angles[cur][i] > angles[prev][i] + M_PI)
 920                         {
 921                             angles[cur][i] -= 2*M_PI;
 922                         }
 923                     }
 924                 }
 925             }
 926         }
 927
 928         /* Average angles */
 929         aa = 0;
 930         for (i = 0; (i < nangles); i++)
 931         {
 932             aa = aa+angles[cur][i];
 933
 934             /* angle in rad / 2Pi * max determines bin. bins go from 0 to maxangstat,
 935                even though scale goes from -pi to pi (dihedral) or -pi/2 to pi/2
 936                (angle) Basically: translate the x-axis by Pi. Translate it back by
 937                -Pi when plotting.
 938              */
 939
 940             angle = angles[cur][i];
 941             if (!bAngles)
 942             {
 943                 while (angle < -M_PI)
 944                 {
 945                     angle += 2*M_PI;
 946                 }
 947                 while (angle >= M_PI)
 948                 {
 949                     angle -= 2*M_PI;
 950                 }
 951
 952                 angle += M_PI;
 953             }
 954
 955             /* Update the distribution histogram */
 956             angind = (int) ((angle*maxangstat)/pifac + 0.5);
 957             if (angind == maxangstat)
 958             {
 959                 angind = 0;
 960             }
 961             if ( (angind < 0) || (angind >= maxangstat) )
 962             {
 963                 /* this will never happen */
 964                 gmx_fatal(FARGS, "angle (%f) index out of range (0..%d) : %d\n",
 965                           angle, maxangstat, angind);
 966             }
 967
 968             angstat[angind]++;
 969             if (angind == maxangstat)
 970             {
 971                 fprintf(stderr, "angle %d fr %d = %g\n", i, cur, angle);
 972             }
 973
 974             total++;
 975         }
 976
 977         /* average over all angles */
 978         (*aver_angle)[teller] = (aa/nangles);
 979
 980         /* this copies all current dih. angles to dih[i], teller is frame */
 981         if (bSaveAll)
 982         {
 983             for (i = 0; i < nangles; i++)
 984             {
 985                 dih[i][teller] = angles[cur][i];
 986             }
 987         }
 988
 989         /* Swap buffers */
 990         cur = prev;
 991
 992         /* Increment loop counter */
 993         teller++;
 994     }
 995     while (read_next_x(oenv, status, &t, x, box));
 996     close_trj(status);
 997
 998     sfree(x);
 999     sfree(angles[cur]);
1000     sfree(angles[prev]);
1001
1002     *nframes = teller;
1003 }