src/tools/anadih.c

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