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