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