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