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