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