src/gromacs/gmxana/gmx_angle.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include <cmath>
  41 #include <cstring>
  42
  43 #include <algorithm>
  44
  45 #include "gromacs/commandline/pargs.h"
  46 #include "gromacs/commandline/viewit.h"
  47 #include "gromacs/correlationfunctions/autocorr.h"
  48 #include "gromacs/fileio/trrio.h"
  49 #include "gromacs/fileio/xvgr.h"
  50 #include "gromacs/gmxana/angle_correction.h"
  51 #include "gromacs/gmxana/gmx_ana.h"
  52 #include "gromacs/gmxana/gstat.h"
  53 #include "gromacs/math/functions.h"
  54 #include "gromacs/math/units.h"
  55 #include "gromacs/math/vec.h"
  56 #include "gromacs/topology/index.h"
  57 #include "gromacs/utility/arraysize.h"
  58 #include "gromacs/utility/fatalerror.h"
  59 #include "gromacs/utility/futil.h"
  60 #include "gromacs/utility/gmxassert.h"
  61 #include "gromacs/utility/pleasecite.h"
  62 #include "gromacs/utility/smalloc.h"
  63
  64 static void dump_dih_trr(int nframes, int nangles, real** dih, const char* fn, real* time)
  65 {
  66     int              i, j, k, l, m, na;
  67     struct t_fileio* fio;
  68     rvec*            x;
  69     matrix           box = { { 2, 0, 0 }, { 0, 2, 0 }, { 0, 0, 2 } };
  70
  71     na = (nangles * 2);
  72     if ((na % 3) != 0)
  73     {
  74         na = 1 + na / 3;
  75     }
  76     else
  77     {
  78         na = na / 3;
  79     }
  80     printf("There are %d dihedrals. Will fill %d atom positions with cos/sin\n", nangles, na);
  81     snew(x, na);
  82     fio = gmx_trr_open(fn, "w");
  83     for (i = 0; (i < nframes); i++)
  84     {
  85         k = l = 0;
  86         for (j = 0; (j < nangles); j++)
  87         {
  88             for (m = 0; (m < 2); m++)
  89             {
  90                 // This is just because the compler and static-analyzer cannot
  91                 // know that dih[j][i] is always valid. Since it occurs in the innermost
  92                 // loop over angles and will only trigger on coding errors, we
  93                 // only enable it for debug builds.
  94                 GMX_ASSERT(dih != nullptr && dih[j] != nullptr, "Incorrect dihedral array data");
  95                 x[k][l] = (m == 0) ? std::cos(dih[j][i]) : std::sin(dih[j][i]);
  96                 l++;
  97                 if (l == DIM)
  98                 {
  99                     l = 0;
 100                     k++;
 101                 }
 102             }
 103         }
 104         gmx_trr_write_frame(fio, i, time[i], 0, box, na, x, nullptr, nullptr);
 105     }
 106     gmx_trr_close(fio);
 107     sfree(x);
 108 }
 109
 110 int gmx_g_angle(int argc, char* argv[])
 111 {
 112     static const char* desc[] = {
 113         "[THISMODULE] computes the angle distribution for a number of angles",
 114         "or dihedrals.[PAR]",
 115         "With option [TT]-ov[tt], you can plot the average angle of",
 116         "a group of angles as a function of time. With the [TT]-all[tt] option,",
 117         "the first graph is the average and the rest are the individual angles.[PAR]",
 118         "With the [TT]-of[tt] option, [THISMODULE] also calculates the fraction of trans",
 119         "dihedrals (only for dihedrals) as function of time, but this is",
 120         "probably only fun for a select few.[PAR]",
 121         "With option [TT]-oc[tt], a dihedral correlation function is calculated.[PAR]",
 122         "It should be noted that the index file must contain",
 123         "atom triplets for angles or atom quadruplets for dihedrals.",
 124         "If this is not the case, the program will crash.[PAR]",
 125         "With option [TT]-or[tt], a trajectory file is dumped containing cos and",
 126         "sin of selected dihedral angles, which subsequently can be used as",
 127         "input for a principal components analysis using [gmx-covar].[PAR]",
 128         "Option [TT]-ot[tt] plots when transitions occur between",
 129         "dihedral rotamers of multiplicity 3 and [TT]-oh[tt]",
 130         "records a histogram of the times between such transitions,",
 131         "assuming the input trajectory frames are equally spaced in time."
 132     };
 133     static const char* opt[] = { nullptr, "angle", "dihedral", "improper", "ryckaert-bellemans",
 134                                  nullptr };
 135     static gmx_bool    bALL = FALSE, bChandler = FALSE, bAverCorr = FALSE, bPBC = TRUE;
 136     static real        binwidth = 1;
 137     t_pargs            pa[]     = {
 138         { "-type", FALSE, etENUM, { opt }, "Type of angle to analyse" },
 139         { "-all",
 140           FALSE,
 141           etBOOL,
 142           { &bALL },
 143           "Plot all angles separately in the averages file, in the order of appearance in the "
 144           "index file." },
 145         { "-binwidth",
 146           FALSE,
 147           etREAL,
 148           { &binwidth },
 149           "binwidth (degrees) for calculating the distribution" },
 150         { "-periodic", FALSE, etBOOL, { &bPBC }, "Print dihedral angles modulo 360 degrees" },
 151         { "-chandler",
 152           FALSE,
 153           etBOOL,
 154           { &bChandler },
 155           "Use Chandler correlation function (N[trans] = 1, N[gauche] = 0) rather than cosine "
 156           "correlation function. Trans is defined as phi < -60 or phi > 60." },
 157         { "-avercorr",
 158           FALSE,
 159           etBOOL,
 160           { &bAverCorr },
 161           "Average the correlation functions for the individual angles/dihedrals" }
 162     };
 163     static const char* bugs[] = {
 164         "Counting transitions only works for dihedrals with multiplicity 3"
 165     };
 166
 167     FILE*         out;
 168     real          dt;
 169     int           isize;
 170     int*          index;
 171     char*         grpname;
 172     real          maxang, S2, norm_fac, maxstat;
 173     unsigned long mode;
 174     int           nframes, maxangstat, mult, *angstat;
 175     int           i, j, nangles, first, last;
 176     gmx_bool      bAver, bRb, bPeriodic, bFrac, /* calculate fraction too?  */
 177             bTrans,                             /* worry about transtions too? */
 178             bCorr;                              /* correlation function ? */
 179     double   tfrac = 0;
 180     char     title[256];
 181     real**   dih = nullptr; /* mega array with all dih. angles at all times*/
 182     real *   time, *trans_frac, *aver_angle;
 183     t_filenm fnm[] = { { efTRX, "-f", nullptr, ffREAD },     { efNDX, nullptr, "angle", ffREAD },
 184                        { efXVG, "-od", "angdist", ffWRITE }, { efXVG, "-ov", "angaver", ffOPTWR },
 185                        { efXVG, "-of", "dihfrac", ffOPTWR }, { efXVG, "-ot", "dihtrans", ffOPTWR },
 186                        { efXVG, "-oh", "trhisto", ffOPTWR }, { efXVG, "-oc", "dihcorr", ffOPTWR },
 187                        { efTRR, "-or", nullptr, ffOPTWR } };
 188 #define NFILE asize(fnm)
 189     int               npargs;
 190     t_pargs*          ppa;
 191     gmx_output_env_t* oenv;
 192
 193     npargs = asize(pa);
 194     ppa    = add_acf_pargs(&npargs, pa);
 195     if (!parse_common_args(
 196                 &argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs, &oenv))
 197     {
 198         sfree(ppa);
 199         return 0;
 200     }
 201
 202     mult   = 4;
 203     maxang = 360.0;
 204     bRb    = FALSE;
 205
 206     GMX_RELEASE_ASSERT(opt[0] != nullptr,
 207                        "Internal option inconsistency; opt[0]==NULL after processing");
 208
 209     switch (opt[0][0])
 210     {
 211         case 'a':
 212             mult   = 3;
 213             maxang = 180.0;
 214             break;
 215         case 'd': // Intended fall through
 216         case 'i': break;
 217         case 'r': bRb = TRUE; break;
 218     }
 219
 220     if (opt2bSet("-or", NFILE, fnm))
 221     {
 222         if (mult != 4)
 223         {
 224             gmx_fatal(FARGS, "Can not combine angles with trr dump");
 225         }
 226         else
 227         {
 228             please_cite(stdout, "Mu2005a");
 229         }
 230     }
 231
 232     /* Calculate bin size */
 233     maxangstat = gmx::roundToInt(maxang / binwidth);
 234     binwidth   = maxang / maxangstat;
 235
 236     rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
 237     nangles = isize / mult;
 238     if ((isize % mult) != 0)
 239     {
 240         gmx_fatal(FARGS,
 241                   "number of index elements not multiple of %d, "
 242                   "these can not be %s\n",
 243                   mult,
 244                   (mult == 3) ? "angle triplets" : "dihedral quadruplets");
 245     }
 246
 247
 248     /* Check whether specific analysis has to be performed */
 249     bCorr  = opt2bSet("-oc", NFILE, fnm);
 250     bAver  = opt2bSet("-ov", NFILE, fnm);
 251     bTrans = opt2bSet("-ot", NFILE, fnm);
 252     bFrac  = opt2bSet("-of", NFILE, fnm);
 253     if (bTrans && opt[0][0] != 'd')
 254     {
 255         fprintf(stderr, "Option -ot should only accompany -type dihedral. Disabling -ot.\n");
 256         bTrans = FALSE;
 257     }
 258
 259     if (bChandler && !bCorr)
 260     {
 261         bCorr = TRUE;
 262     }
 263
 264     if (bFrac && !bRb)
 265     {
 266         fprintf(stderr,
 267                 "Warning:"
 268                 " calculating fractions as defined in this program\n"
 269                 "makes sense for Ryckaert Bellemans dihs. only. Ignoring -of\n\n");
 270         bFrac = FALSE;
 271     }
 272
 273     if ((bTrans || bFrac || bCorr) && mult == 3)
 274     {
 275         gmx_fatal(FARGS,
 276                   "Can only do transition, fraction or correlation\n"
 277                   "on dihedrals. Select -d\n");
 278     }
 279
 280     /*
 281      * We need to know the nr of frames so we can allocate memory for an array
 282      * with all dihedral angles at all timesteps. Works for me.
 283      */
 284     if (bTrans || bCorr || bALL || opt2bSet("-or", NFILE, fnm))
 285     {
 286         snew(dih, nangles);
 287     }
 288
 289     snew(angstat, maxangstat);
 290
 291     read_ang_dih(ftp2fn(efTRX, NFILE, fnm),
 292                  (mult == 3),
 293                  bALL || bCorr || bTrans || opt2bSet("-or", NFILE, fnm),
 294                  bRb,
 295                  bPBC,
 296                  maxangstat,
 297                  angstat,
 298                  &nframes,
 299                  &time,
 300                  isize,
 301                  index,
 302                  &trans_frac,
 303                  &aver_angle,
 304                  dih,
 305                  oenv);
 306
 307     dt = (time[nframes - 1] - time[0]) / (nframes - 1);
 308
 309     if (bAver)
 310     {
 311         sprintf(title, "Average Angle: %s", grpname);
 312         out = xvgropen(opt2fn("-ov", NFILE, fnm), title, "Time (ps)", "Angle (degrees)", oenv);
 313         for (i = 0; (i < nframes); i++)
 314         {
 315             fprintf(out, "%10.5f  %8.3f", time[i], aver_angle[i] * RAD2DEG);
 316             if (bALL)
 317             {
 318                 for (j = 0; (j < nangles); j++)
 319                 {
 320                     if (bPBC)
 321                     {
 322                         real dd = dih[j][i];
 323                         fprintf(out, "  %8.3f", std::atan2(std::sin(dd), std::cos(dd)) * RAD2DEG);
 324                     }
 325                     else
 326                     {
 327                         fprintf(out, "  %8.3f", dih[j][i] * RAD2DEG);
 328                     }
 329                 }
 330             }
 331             fprintf(out, "\n");
 332         }
 333         xvgrclose(out);
 334     }
 335     if (opt2bSet("-or", NFILE, fnm))
 336     {
 337         dump_dih_trr(nframes, nangles, dih, opt2fn("-or", NFILE, fnm), time);
 338     }
 339
 340     if (bFrac)
 341     {
 342         sprintf(title, "Trans fraction: %s", grpname);
 343         out   = xvgropen(opt2fn("-of", NFILE, fnm), title, "Time (ps)", "Fraction", oenv);
 344         tfrac = 0.0;
 345         for (i = 0; (i < nframes); i++)
 346         {
 347             fprintf(out, "%10.5f  %10.3f\n", time[i], trans_frac[i]);
 348             tfrac += trans_frac[i];
 349         }
 350         xvgrclose(out);
 351
 352         tfrac /= nframes;
 353         fprintf(stderr, "Average trans fraction: %g\n", tfrac);
 354     }
 355     sfree(trans_frac);
 356
 357     if (bTrans)
 358     {
 359         ana_dih_trans(
 360                 opt2fn("-ot", NFILE, fnm), opt2fn("-oh", NFILE, fnm), dih, nframes, nangles, grpname, time, bRb, oenv);
 361     }
 362
 363     if (bCorr)
 364     {
 365         /* Autocorrelation function */
 366         if (nframes < 2)
 367         {
 368             fprintf(stderr, "Not enough frames for correlation function\n");
 369         }
 370         else
 371         {
 372
 373             if (bChandler)
 374             {
 375                 real     dval, sixty = DEG2RAD * 60;
 376                 gmx_bool bTest;
 377
 378                 for (i = 0; (i < nangles); i++)
 379                 {
 380                     for (j = 0; (j < nframes); j++)
 381                     {
 382                         dval = dih[i][j];
 383                         if (bRb)
 384                         {
 385                             bTest = (dval > -sixty) && (dval < sixty);
 386                         }
 387                         else
 388                         {
 389                             bTest = (dval < -sixty) || (dval > sixty);
 390                         }
 391                         if (bTest)
 392                         {
 393                             dih[i][j] = dval - tfrac;
 394                         }
 395                         else
 396                         {
 397                             dih[i][j] = -tfrac;
 398                         }
 399                     }
 400                 }
 401             }
 402             if (bChandler)
 403             {
 404                 mode = eacNormal;
 405             }
 406             else
 407             {
 408                 mode = eacCos;
 409             }
 410             do_autocorr(opt2fn("-oc", NFILE, fnm),
 411                         oenv,
 412                         "Dihedral Autocorrelation Function",
 413                         nframes,
 414                         nangles,
 415                         dih,
 416                         dt,
 417                         mode,
 418                         bAverCorr);
 419         }
 420     }
 421
 422
 423     /* Determine the non-zero part of the distribution */
 424     for (first = 0; (first < maxangstat - 1) && (angstat[first + 1] == 0); first++) {}
 425     for (last = maxangstat - 1; (last > 0) && (angstat[last - 1] == 0); last--) {}
 426
 427     double aver = 0;
 428     printf("Found points in the range from %d to %d (max %d)\n", first, last, maxangstat);
 429     if (bTrans || bCorr || bALL || opt2bSet("-or", NFILE, fnm))
 430     { /* It's better to re-calculate Std. Dev per sample */
 431         real b_aver = aver_angle[0];
 432         real b      = dih[0][0];
 433         real delta;
 434         for (int i = 0; (i < nframes); i++)
 435         {
 436             delta = correctRadianAngleRange(aver_angle[i] - b_aver);
 437             b_aver += delta;
 438             aver += b_aver;
 439             for (int j = 0; (j < nangles); j++)
 440             {
 441                 delta = correctRadianAngleRange(dih[j][i] - b);
 442                 b += delta;
 443             }
 444         }
 445     }
 446     else
 447     { /* Incorrect  for Std. Dev. */
 448         real delta, b_aver = aver_angle[0];
 449         for (i = 0; (i < nframes); i++)
 450         {
 451             delta = correctRadianAngleRange(aver_angle[i] - b_aver);
 452             b_aver += delta;
 453             aver += b_aver;
 454         }
 455     }
 456     aver /= nframes;
 457     double aversig = correctRadianAngleRange(aver);
 458     aversig *= RAD2DEG;
 459     aver *= RAD2DEG;
 460     printf(" < angle >  = %g\n", aversig);
 461
 462     if (mult == 3)
 463     {
 464         sprintf(title, "Angle Distribution: %s", grpname);
 465     }
 466     else
 467     {
 468         sprintf(title, "Dihedral Distribution: %s", grpname);
 469
 470         calc_distribution_props(maxangstat, angstat, -180.0, 0, nullptr, &S2);
 471         fprintf(stderr, "Order parameter S^2 = %g\n", S2);
 472     }
 473
 474     bPeriodic = (mult == 4) && (first == 0) && (last == maxangstat - 1);
 475
 476     out = xvgropen(opt2fn("-od", NFILE, fnm), title, "Degrees", "", oenv);
 477     if (output_env_get_print_xvgr_codes(oenv))
 478     {
 479         fprintf(out, "@    subtitle \"average angle: %g\\So\\N\"\n", aver);
 480     }
 481     norm_fac = 1.0 / (nangles * nframes * binwidth);
 482     if (bPeriodic)
 483     {
 484         maxstat = 0;
 485         for (i = first; (i <= last); i++)
 486         {
 487             maxstat = std::max(maxstat, angstat[i] * norm_fac);
 488         }
 489         if (output_env_get_print_xvgr_codes(oenv))
 490         {
 491             fprintf(out, "@with g0\n");
 492             fprintf(out, "@    world xmin -180\n");
 493             fprintf(out, "@    world xmax  180\n");
 494             fprintf(out, "@    world ymin 0\n");
 495             fprintf(out, "@    world ymax %g\n", maxstat * 1.05);
 496             fprintf(out, "@    xaxis  tick major 60\n");
 497             fprintf(out, "@    xaxis  tick minor 30\n");
 498             fprintf(out, "@    yaxis  tick major 0.005\n");
 499             fprintf(out, "@    yaxis  tick minor 0.0025\n");
 500         }
 501     }
 502     for (i = first; (i <= last); i++)
 503     {
 504         fprintf(out, "%10g  %10f\n", i * binwidth + 180.0 - maxang, angstat[i] * norm_fac);
 505     }
 506     if (bPeriodic)
 507     {
 508         /* print first bin again as last one */
 509         fprintf(out, "%10g  %10f\n", 180.0, angstat[0] * norm_fac);
 510     }
 511
 512     xvgrclose(out);
 513
 514     do_view(oenv, opt2fn("-od", NFILE, fnm), "-nxy");
 515     if (bAver)
 516     {
 517         do_view(oenv, opt2fn("-ov", NFILE, fnm), "-nxy");
 518     }
 519
 520     return 0;
 521 }