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