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