src/gromacs/gmxana/anadih.cpp

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