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(
 165             TRUE, fn_trans, TRUE, fn_histo, maxchi, dih, nlist, dlist, nframes, nangles, 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", "cumulative rotamer", nframes, time, chi_prtrj);
 579             }
 580
 581             /* make a histogram pf culm. rotamer occupancy too */
 582             snew(chi_prhist, nbin);
 583             make_histo(nullptr, nframes, chi_prtrj, nbin, chi_prhist, 0, nbin);
 584             if (bAll)
 585             {
 586                 sprintf(hisfile, "histo-chiprod%s.xvg", dlist[i].name);
 587                 sprintf(histitle, "cumulative rotamer distribution for %s", dlist[i].name);
 588                 fprintf(stderr, "  and %s  ", hisfile);
 589                 fp = xvgropen(hisfile, histitle, "number", "", oenv);
 590                 if (output_env_get_print_xvgr_codes(oenv))
 591                 {
 592                     fprintf(fp, "@ xaxis tick on\n");
 593                     fprintf(fp, "@ xaxis tick major 1\n");
 594                     fprintf(fp, "@ type xy\n");
 595                 }
 596                 for (k = 0; (k < nbin); k++)
 597                 {
 598                     if (bNormalize)
 599                     {
 600                         fprintf(fp, "%5d  %10g\n", k, (1.0 * chi_prhist[k]) / nframes);
 601                     }
 602                     else
 603                     {
 604                         fprintf(fp, "%5d  %10d\n", k, chi_prhist[k]);
 605                     }
 606                 }
 607                 fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 608                 xvgrclose(fp);
 609             }
 610
 611             /* and finally print out occupancies to a single file */
 612             /* get the gmx from-1 res nr by setting a ptr to the number part
 613              * of dlist[i].name - potential bug for 4-letter res names... */
 614             namept = dlist[i].name + 3;
 615             fprintf(fpall, "%5s ", namept);
 616             for (k = 0; (k < nbin); k++)
 617             {
 618                 if (bNormalize)
 619                 {
 620                     fprintf(fpall, "  %10g", (1.0 * chi_prhist[k]) / nframes);
 621                 }
 622                 else
 623                 {
 624                     fprintf(fpall, "  %10d", chi_prhist[k]);
 625                 }
 626             }
 627             fprintf(fpall, "\n");
 628
 629             sfree(chi_prhist);
 630             /* histogram done */
 631         }
 632     }
 633
 634     sfree(chi_prtrj);
 635     xvgrclose(fpall);
 636     fprintf(stderr, "\n");
 637 }
 638
 639 void calc_distribution_props(int nh, const int histo[], real start, int nkkk, t_karplus kkk[], real* S2)
 640 {
 641     real d, dc, ds, c1, c2, tdc, tds;
 642     real fac, ang, invth, Jc;
 643     int  i, j, th;
 644
 645     if (nh == 0)
 646     {
 647         gmx_fatal(FARGS, "No points in histogram (%s, %d)", __FILE__, __LINE__);
 648     }
 649     fac = 2 * M_PI / nh;
 650
 651     /* Compute normalisation factor */
 652     th = 0;
 653     for (j = 0; (j < nh); j++)
 654     {
 655         th += histo[j];
 656     }
 657     invth = 1.0 / th;
 658
 659     for (i = 0; (i < nkkk); i++)
 660     {
 661         kkk[i].Jc    = 0;
 662         kkk[i].Jcsig = 0;
 663     }
 664     tdc = 0;
 665     tds = 0;
 666     for (j = 0; (j < nh); j++)
 667     {
 668         d   = invth * histo[j];
 669         ang = j * fac - start;
 670         c1  = std::cos(ang);
 671         dc  = d * c1;
 672         ds  = d * std::sin(ang);
 673         tdc += dc;
 674         tds += ds;
 675         for (i = 0; (i < nkkk); i++)
 676         {
 677             c1 = std::cos(ang + kkk[i].offset);
 678             c2 = c1 * c1;
 679             Jc = (kkk[i].A * c2 + kkk[i].B * c1 + kkk[i].C);
 680             kkk[i].Jc += histo[j] * Jc;
 681             kkk[i].Jcsig += histo[j] * gmx::square(Jc);
 682         }
 683     }
 684     for (i = 0; (i < nkkk); i++)
 685     {
 686         kkk[i].Jc /= th;
 687         kkk[i].Jcsig = std::sqrt(kkk[i].Jcsig / th - gmx::square(kkk[i].Jc));
 688     }
 689     *S2 = tdc * tdc + tds * tds;
 690 }
 691
 692 static void calc_angles(struct t_pbc* pbc, int n3, int index[], real ang[], rvec x_s[])
 693 {
 694     int  i, ix, t1, t2;
 695     rvec r_ij, r_kj;
 696     real costh = 0.0;
 697
 698     for (i = ix = 0; (ix < n3); i++, ix += 3)
 699     {
 700         ang[i] = bond_angle(
 701                 x_s[index[ix]], x_s[index[ix + 1]], x_s[index[ix + 2]], pbc, r_ij, r_kj, &costh, &t1, &t2);
 702     }
 703     if (debug)
 704     {
 705         fprintf(debug, "Angle[0]=%g, costh=%g, index0 = %d, %d, %d\n", ang[0], costh, index[0], index[1], index[2]);
 706         pr_rvec(debug, 0, "rij", r_ij, DIM, TRUE);
 707         pr_rvec(debug, 0, "rkj", r_kj, DIM, TRUE);
 708     }
 709 }
 710
 711 static real calc_fraction(const real angles[], int nangles)
 712 {
 713     int  i;
 714     real trans = 0, gauche = 0;
 715     real angle;
 716
 717     for (i = 0; i < nangles; i++)
 718     {
 719         angle = angles[i] * RAD2DEG;
 720
 721         if (angle > 135 && angle < 225)
 722         {
 723             trans += 1.0;
 724         }
 725         else if ((angle > 270 && angle < 330) || (angle < 90 && angle > 30))
 726         {
 727             gauche += 1.0;
 728         }
 729     }
 730     if (trans + gauche > 0)
 731     {
 732         return trans / (trans + gauche);
 733     }
 734     else
 735     {
 736         return 0;
 737     }
 738 }
 739
 740 static void calc_dihs(struct t_pbc* pbc, int n4, const int index[], real ang[], rvec x_s[])
 741 {
 742     int  i, ix, t1, t2, t3;
 743     rvec r_ij, r_kj, r_kl, m, n;
 744     real aaa;
 745
 746     for (i = ix = 0; (ix < n4); i++, ix += 4)
 747     {
 748         aaa = dih_angle(x_s[index[ix]],
 749                         x_s[index[ix + 1]],
 750                         x_s[index[ix + 2]],
 751                         x_s[index[ix + 3]],
 752                         pbc,
 753                         r_ij,
 754                         r_kj,
 755                         r_kl,
 756                         m,
 757                         n,
 758                         &t1,
 759                         &t2,
 760                         &t3);
 761
 762         ang[i] = aaa; /* not taking into account ryckaert bellemans yet */
 763     }
 764 }
 765
 766 void make_histo(FILE* log, int ndata, real data[], int npoints, int histo[], real minx, real maxx)
 767 {
 768     double dx;
 769     int    i, ind;
 770
 771     if (minx == maxx)
 772     {
 773         minx = maxx = data[0];
 774         for (i = 1; (i < ndata); i++)
 775         {
 776             minx = std::min(minx, data[i]);
 777             maxx = std::max(maxx, data[i]);
 778         }
 779         fprintf(log, "Min data: %10g  Max data: %10g\n", minx, maxx);
 780     }
 781     dx = npoints / (maxx - minx);
 782     if (debug)
 783     {
 784         fprintf(debug, "Histogramming: ndata=%d, nhisto=%d, minx=%g,maxx=%g,dx=%g\n", ndata, npoints, minx, maxx, dx);
 785     }
 786     for (i = 0; (i < ndata); i++)
 787     {
 788         ind = static_cast<int>((data[i] - minx) * dx);
 789         if ((ind >= 0) && (ind < npoints))
 790         {
 791             histo[ind]++;
 792         }
 793         else
 794         {
 795             fprintf(log, "index = %d, data[%d] = %g\n", ind, i, data[i]);
 796         }
 797     }
 798 }
 799
 800 void normalize_histo(int npoints, const int histo[], real dx, real normhisto[])
 801 {
 802     int    i;
 803     double d, fac;
 804
 805     d = 0;
 806     for (i = 0; (i < npoints); i++)
 807     {
 808         d += dx * histo[i];
 809     }
 810     if (d == 0)
 811     {
 812         fprintf(stderr, "Empty histogram!\n");
 813         return;
 814     }
 815     fac = 1.0 / d;
 816     for (i = 0; (i < npoints); i++)
 817     {
 818         normhisto[i] = fac * histo[i];
 819     }
 820 }
 821
 822 void read_ang_dih(const char*             trj_fn,
 823                   gmx_bool                bAngles,
 824                   gmx_bool                bSaveAll,
 825                   gmx_bool                bRb,
 826                   gmx_bool                bPBC,
 827                   int                     maxangstat,
 828                   int                     angstat[],
 829                   int*                    nframes,
 830                   real**                  time,
 831                   int                     isize,
 832                   int                     index[],
 833                   real**                  trans_frac,
 834                   real**                  aver_angle,
 835                   real*                   dih[],
 836                   const gmx_output_env_t* oenv)
 837 {
 838     struct t_pbc* pbc;
 839     t_trxstatus*  status;
 840     int           i, angind, total, teller;
 841     int           nangles, n_alloc;
 842     real          t, fraction, pifac, angle;
 843     real*         angles[2];
 844     matrix        box;
 845     rvec*         x;
 846     int           cur = 0;
 847 #define prev (1 - cur)
 848
 849     snew(pbc, 1);
 850     read_first_x(oenv, &status, trj_fn, &t, &x, box);
 851
 852     if (bAngles)
 853     {
 854         nangles = isize / 3;
 855         pifac   = M_PI;
 856     }
 857     else
 858     {
 859         nangles = isize / 4;
 860         pifac   = 2.0 * M_PI;
 861     }
 862     snew(angles[cur], nangles);
 863     snew(angles[prev], nangles);
 864
 865     /* Start the loop over frames */
 866     total       = 0;
 867     teller      = 0;
 868     n_alloc     = 0;
 869     *time       = nullptr;
 870     *trans_frac = nullptr;
 871     *aver_angle = nullptr;
 872
 873     do
 874     {
 875         if (teller >= n_alloc)
 876         {
 877             n_alloc += 100;
 878             if (bSaveAll)
 879             {
 880                 for (i = 0; (i < nangles); i++)
 881                 {
 882                     srenew(dih[i], n_alloc);
 883                 }
 884             }
 885             srenew(*time, n_alloc);
 886             srenew(*trans_frac, n_alloc);
 887             srenew(*aver_angle, n_alloc);
 888         }
 889
 890         (*time)[teller] = t;
 891
 892         if (pbc)
 893         {
 894             set_pbc(pbc, PbcType::Unset, box);
 895         }
 896
 897         if (bAngles)
 898         {
 899             calc_angles(pbc, isize, index, angles[cur], x);
 900         }
 901         else
 902         {
 903             calc_dihs(pbc, isize, index, angles[cur], x);
 904
 905             /* Trans fraction */
 906             fraction              = calc_fraction(angles[cur], nangles);
 907             (*trans_frac)[teller] = fraction;
 908
 909             /* Change Ryckaert-Bellemans dihedrals to polymer convention
 910              * Modified 990913 by Erik:
 911              * We actually shouldn't change the convention, since it's
 912              * calculated from polymer above, but we change the intervall
 913              * from [-180,180] to [0,360].
 914              */
 915             if (bRb)
 916             {
 917                 for (i = 0; (i < nangles); i++)
 918                 {
 919                     if (angles[cur][i] <= 0.0)
 920                     {
 921                         angles[cur][i] += 2 * M_PI;
 922                     }
 923                 }
 924             }
 925
 926             /* Periodicity in dihedral space... */
 927             if (bPBC)
 928             {
 929                 for (i = 0; (i < nangles); i++)
 930                 {
 931                     real dd        = angles[cur][i];
 932                     angles[cur][i] = std::atan2(std::sin(dd), std::cos(dd));
 933                 }
 934             }
 935             else
 936             {
 937                 if (teller > 1)
 938                 {
 939                     for (i = 0; (i < nangles); i++)
 940                     {
 941                         while (angles[cur][i] <= angles[prev][i] - M_PI)
 942                         {
 943                             angles[cur][i] += 2 * M_PI;
 944                         }
 945                         while (angles[cur][i] > angles[prev][i] + M_PI)
 946                         {
 947                             angles[cur][i] -= 2 * M_PI;
 948                         }
 949                     }
 950                 }
 951             }
 952         }
 953
 954         /* Average angles */
 955         double aa = 0;
 956         for (i = 0; (i < nangles); i++)
 957         {
 958             if (!bAngles && i > 0)
 959             {
 960                 real diffa     = angles[cur][i] - angles[cur][i - 1];
 961                 diffa          = correctRadianAngleRange(diffa);
 962                 angles[cur][i] = angles[cur][i - 1] + diffa;
 963             }
 964
 965             aa = aa + angles[cur][i];
 966
 967             /* angle in rad / 2Pi * max determines bin. bins go from 0 to maxangstat,
 968                even though scale goes from -pi to pi (dihedral) or -pi/2 to pi/2
 969                (angle) Basically: translate the x-axis by Pi. Translate it back by
 970                -Pi when plotting.
 971              */
 972
 973             angle = angles[cur][i];
 974             if (!bAngles)
 975             {
 976                 angle = correctRadianAngleRange(angle);
 977                 angle += M_PI;
 978             }
 979
 980             /* Update the distribution histogram */
 981             angind = gmx::roundToInt((angle * maxangstat) / pifac);
 982             if (angind == maxangstat)
 983             {
 984                 angind = 0;
 985             }
 986             if ((angind < 0) || (angind >= maxangstat))
 987             {
 988                 /* this will never happen */
 989                 gmx_fatal(FARGS, "angle (%f) index out of range (0..%d) : %d\n", angle, maxangstat, angind);
 990             }
 991
 992             angstat[angind]++;
 993             if (angind == maxangstat)
 994             {
 995                 fprintf(stderr, "angle %d fr %d = %g\n", i, cur, angle);
 996             }
 997
 998             total++;
 999         }
1000
1001         /* average over all angles */
1002         aa                    = correctRadianAngleRange(aa / nangles);
1003         (*aver_angle)[teller] = (aa);
1004
1005         /* this copies all current dih. angles to dih[i], teller is frame */
1006         if (bSaveAll)
1007         {
1008             for (i = 0; i < nangles; i++)
1009             {
1010                 if (!bAngles)
1011                 {
1012                     dih[i][teller] = correctRadianAngleRange(angles[cur][i]);
1013                 }
1014                 else
1015                 {
1016                     dih[i][teller] = angles[cur][i];
1017                 }
1018             }
1019         }
1020
1021         /* Swap buffers */
1022         cur = prev;
1023
1024         /* Increment loop counter */
1025         teller++;
1026     } while (read_next_x(oenv, status, &t, x, box));
1027     close_trx(status);
1028
1029     sfree(x);
1030     sfree(angles[cur]);
1031     sfree(angles[prev]);
1032
1033     *nframes = teller;
1034 }