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