src/gromacs/gmxana/gmx_rms.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include <cmath>
  41 #include <cstdlib>
  42
  43 #include <algorithm>
  44
  45 #include "gromacs/commandline/pargs.h"
  46 #include "gromacs/commandline/viewit.h"
  47 #include "gromacs/fileio/confio.h"
  48 #include "gromacs/fileio/matio.h"
  49 #include "gromacs/fileio/trxio.h"
  50 #include "gromacs/fileio/xvgr.h"
  51 #include "gromacs/gmxana/cmat.h"
  52 #include "gromacs/gmxana/gmx_ana.h"
  53 #include "gromacs/gmxana/princ.h"
  54 #include "gromacs/math/do_fit.h"
  55 #include "gromacs/math/functions.h"
  56 #include "gromacs/math/units.h"
  57 #include "gromacs/math/utilities.h"
  58 #include "gromacs/math/vec.h"
  59 #include "gromacs/pbcutil/rmpbc.h"
  60 #include "gromacs/topology/ifunc.h"
  61 #include "gromacs/topology/index.h"
  62 #include "gromacs/topology/topology.h"
  63 #include "gromacs/utility/arraysize.h"
  64 #include "gromacs/utility/fatalerror.h"
  65 #include "gromacs/utility/futil.h"
  66 #include "gromacs/utility/pleasecite.h"
  67 #include "gromacs/utility/smalloc.h"
  68
  69 static void norm_princ(const t_atoms* atoms, int isize, int* index, int natoms, rvec* x)
  70 {
  71     int  i, m;
  72     rvec princ, vec;
  73
  74     /* equalize principal components: */
  75     /* orient principal axes, get principal components */
  76     orient_princ(atoms, isize, index, natoms, x, nullptr, princ);
  77
  78     /* calc our own principal components */
  79     clear_rvec(vec);
  80     for (m = 0; m < DIM; m++)
  81     {
  82         for (i = 0; i < isize; i++)
  83         {
  84             vec[m] += gmx::square(x[index[i]][m]);
  85         }
  86         vec[m] = std::sqrt(vec[m] / static_cast<real>(isize));
  87         /* calculate scaling constants */
  88         vec[m] = 1.0 / (std::sqrt(3.0) * vec[m]);
  89     }
  90
  91     /* scale coordinates */
  92     for (i = 0; i < natoms; i++)
  93     {
  94         for (m = 0; m < DIM; m++)
  95         {
  96             x[i][m] *= vec[m];
  97         }
  98     }
  99 }
 100
 101 int gmx_rms(int argc, char* argv[])
 102 {
 103     const char* desc[] = {
 104         "[THISMODULE] compares two structures by computing the root mean square",
 105         "deviation (RMSD), the size-independent [GRK]rho[grk] similarity parameter",
 106         "([TT]rho[tt]) or the scaled [GRK]rho[grk] ([TT]rhosc[tt]), ",
 107         "see Maiorov & Crippen, Proteins [BB]22[bb], 273 (1995).",
 108         "This is selected by [TT]-what[tt].[PAR]",
 109
 110         "Each structure from a trajectory ([TT]-f[tt]) is compared to a",
 111         "reference structure. The reference structure",
 112         "is taken from the structure file ([TT]-s[tt]).[PAR]",
 113
 114         "With option [TT]-mir[tt] also a comparison with the mirror image of",
 115         "the reference structure is calculated.",
 116         "This is useful as a reference for 'significant' values, see",
 117         "Maiorov & Crippen, Proteins [BB]22[bb], 273 (1995).[PAR]",
 118
 119         "Option [TT]-prev[tt] produces the comparison with a previous frame",
 120         "the specified number of frames ago.[PAR]",
 121
 122         "Option [TT]-m[tt] produces a matrix in [REF].xpm[ref] format of",
 123         "comparison values of each structure in the trajectory with respect to",
 124         "each other structure. This file can be visualized with for instance",
 125         "[TT]xv[tt] and can be converted to postscript with [gmx-xpm2ps].[PAR]",
 126
 127         "Option [TT]-fit[tt] controls the least-squares fitting of",
 128         "the structures on top of each other: complete fit (rotation and",
 129         "translation), translation only, or no fitting at all.[PAR]",
 130
 131         "Option [TT]-mw[tt] controls whether mass weighting is done or not.",
 132         "If you select the option (default) and ",
 133         "supply a valid [REF].tpr[ref] file masses will be taken from there, ",
 134         "otherwise the masses will be deduced from the [TT]atommass.dat[tt] file in",
 135         "[TT]GMXLIB[tt]. This is fine for proteins, but not",
 136         "necessarily for other molecules. A default mass of 12.011 amu (carbon)",
 137         "is assigned to unknown atoms. You can check whether this happened by",
 138         "turning on the [TT]-debug[tt] flag and inspecting the log file.[PAR]",
 139
 140         "With [TT]-f2[tt], the 'other structures' are taken from a second",
 141         "trajectory, this generates a comparison matrix of one trajectory",
 142         "versus the other.[PAR]",
 143
 144         "Option [TT]-bin[tt] does a binary dump of the comparison matrix.[PAR]",
 145
 146         "Option [TT]-bm[tt] produces a matrix of average bond angle deviations",
 147         "analogously to the [TT]-m[tt] option. Only bonds between atoms in the",
 148         "comparison group are considered."
 149     };
 150     static gmx_bool bPBC = TRUE, bFitAll = TRUE, bSplit = FALSE;
 151     static gmx_bool bDeltaLog = FALSE;
 152     static int      prev = 0, freq = 1, freq2 = 1, nlevels = 80, avl = 0;
 153     static real     rmsd_user_max = -1, rmsd_user_min = -1, bond_user_max = -1, bond_user_min = -1,
 154                 delta_maxy = 0.0;
 155     /* strings and things for selecting difference method */
 156     enum
 157     {
 158         ewSel,
 159         ewRMSD,
 160         ewRho,
 161         ewRhoSc,
 162         ewNR
 163     };
 164     int         ewhat;
 165     const char* what[ewNR + 1]     = { nullptr, "rmsd", "rho", "rhosc", nullptr };
 166     const char* whatname[ewNR]     = { nullptr, "RMSD", "Rho", "Rho sc" };
 167     const char* whatlabel[ewNR]    = { nullptr, "RMSD (nm)", "Rho", "Rho sc" };
 168     const char* whatxvgname[ewNR]  = { nullptr, "RMSD", "\\8r\\4", "\\8r\\4\\ssc\\N" };
 169     const char* whatxvglabel[ewNR] = { nullptr, "RMSD (nm)", "\\8r\\4", "\\8r\\4\\ssc\\N" };
 170     /* strings and things for fitting methods */
 171     enum
 172     {
 173         efSel,
 174         efFit,
 175         efReset,
 176         efNone,
 177         efNR
 178     };
 179     int             efit;
 180     const char*     fit[efNR + 1] = { nullptr, "rot+trans", "translation", "none", nullptr };
 181     const char*     fitgraphlabel[efNR + 1] = { nullptr, "lsq fit", "translational fit", "no fit" };
 182     static int      nrms                    = 1;
 183     static gmx_bool bMassWeighted           = TRUE;
 184     t_pargs         pa[]                    = {
 185         { "-what", FALSE, etENUM, { what }, "Structural difference measure" },
 186         { "-pbc", FALSE, etBOOL, { &bPBC }, "PBC check" },
 187         { "-fit", FALSE, etENUM, { fit }, "Fit to reference structure" },
 188         { "-prev", FALSE, etINT, { &prev }, "Compare with previous frame" },
 189         { "-split", FALSE, etBOOL, { &bSplit }, "Split graph where time is zero" },
 190         { "-fitall", FALSE, etBOOL, { &bFitAll }, "HIDDENFit all pairs of structures in matrix" },
 191         { "-skip", FALSE, etINT, { &freq }, "Only write every nr-th frame to matrix" },
 192         { "-skip2", FALSE, etINT, { &freq2 }, "Only write every nr-th frame to matrix" },
 193         { "-max", FALSE, etREAL, { &rmsd_user_max }, "Maximum level in comparison matrix" },
 194         { "-min", FALSE, etREAL, { &rmsd_user_min }, "Minimum level in comparison matrix" },
 195         { "-bmax", FALSE, etREAL, { &bond_user_max }, "Maximum level in bond angle matrix" },
 196         { "-bmin", FALSE, etREAL, { &bond_user_min }, "Minimum level in bond angle matrix" },
 197         { "-mw", FALSE, etBOOL, { &bMassWeighted }, "Use mass weighting for superposition" },
 198         { "-nlevels", FALSE, etINT, { &nlevels }, "Number of levels in the matrices" },
 199         { "-ng", FALSE, etINT, { &nrms }, "Number of groups to compute RMS between" },
 200         { "-dlog", FALSE, etBOOL, { &bDeltaLog }, "HIDDENUse a log x-axis in the delta t matrix" },
 201         { "-dmax", FALSE, etREAL, { &delta_maxy }, "HIDDENMaximum level in delta matrix" },
 202         { "-aver", FALSE, etINT, { &avl }, "HIDDENAverage over this distance in the RMSD matrix" }
 203     };
 204     int natoms_trx, natoms_trx2, natoms;
 205     int i, j, k, m;
 206 #define NFRAME 5000
 207     int        maxframe = NFRAME, maxframe2 = NFRAME;
 208     real       t, *w_rls, *w_rms, *w_rls_m = nullptr, *w_rms_m = nullptr;
 209     gmx_bool   bNorm, bAv, bFreq2, bFile2, bMat, bBond, bDelta, bMirror, bMass;
 210     gmx_bool   bFit, bReset;
 211     t_topology top;
 212     PbcType    pbcType;
 213     t_iatom*   iatom = nullptr;
 214
 215     matrix box = { { 0 } };
 216     rvec * x, *xp, *xm = nullptr, **mat_x = nullptr, **mat_x2, *mat_x2_j = nullptr, vec1, vec2;
 217     t_trxstatus* status;
 218     char         buf[256], buf2[256];
 219     int          ncons = 0;
 220     FILE*        fp;
 221     real         rlstot = 0, **rls, **rlsm = nullptr, *time, *time2, *rlsnorm = nullptr,
 222          **rmsd_mat = nullptr, **bond_mat = nullptr, *axis, *axis2, *del_xaxis, *del_yaxis,
 223          rmsd_max, rmsd_min, rmsd_avg, bond_max, bond_min, ang;
 224     real **  rmsdav_mat = nullptr, av_tot, weight, weight_tot;
 225     real **  delta = nullptr, delta_max, delta_scalex = 0, delta_scaley = 0, *delta_tot;
 226     int      delta_xsize = 0, del_lev = 100, mx, my, abs_my;
 227     gmx_bool bA1, bA2, bPrev, bTop, *bInMat = nullptr;
 228     int      ifit, *irms, ibond = 0, *ind_bond1 = nullptr, *ind_bond2 = nullptr, n_ind_m = 0;
 229     int *    ind_fit, **ind_rms, *ind_m = nullptr, *rev_ind_m = nullptr, *ind_rms_m = nullptr;
 230     char *   gn_fit, **gn_rms;
 231     t_rgb    rlo, rhi;
 232     gmx_output_env_t* oenv;
 233     gmx_rmpbc_t       gpbc = nullptr;
 234
 235     t_filenm fnm[] = {
 236         { efTPS, nullptr, nullptr, ffREAD }, { efTRX, "-f", nullptr, ffREAD },
 237         { efTRX, "-f2", nullptr, ffOPTRD },  { efNDX, nullptr, nullptr, ffOPTRD },
 238         { efXVG, nullptr, "rmsd", ffWRITE }, { efXVG, "-mir", "rmsdmir", ffOPTWR },
 239         { efXVG, "-a", "avgrp", ffOPTWR },   { efXVG, "-dist", "rmsd-dist", ffOPTWR },
 240         { efXPM, "-m", "rmsd", ffOPTWR },    { efDAT, "-bin", "rmsd", ffOPTWR },
 241         { efXPM, "-bm", "bond", ffOPTWR }
 242     };
 243 #define NFILE asize(fnm)
 244
 245     if (!parse_common_args(&argc,
 246                            argv,
 247                            PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW,
 248                            NFILE,
 249                            fnm,
 250                            asize(pa),
 251                            pa,
 252                            asize(desc),
 253                            desc,
 254                            0,
 255                            nullptr,
 256                            &oenv))
 257     {
 258         return 0;
 259     }
 260     /* parse enumerated options: */
 261     ewhat = nenum(what);
 262     if (ewhat == ewRho || ewhat == ewRhoSc)
 263     {
 264         please_cite(stdout, "Maiorov95");
 265     }
 266     efit   = nenum(fit);
 267     bFit   = efit == efFit;
 268     bReset = efit == efReset;
 269     if (bFit)
 270     {
 271         bReset = TRUE; /* for fit, reset *must* be set */
 272     }
 273     else
 274     {
 275         bFitAll = FALSE;
 276     }
 277
 278     /* mark active cmdline options */
 279     bMirror = opt2bSet("-mir", NFILE, fnm); /* calc RMSD vs mirror of ref. */
 280     bFile2  = opt2bSet("-f2", NFILE, fnm);
 281     bMat    = opt2bSet("-m", NFILE, fnm);
 282     bBond   = opt2bSet("-bm", NFILE, fnm);
 283     bDelta  = (delta_maxy > 0); /* calculate rmsd vs delta t matrix from *
 284                                  *      your RMSD matrix (hidden option       */
 285     bNorm  = opt2bSet("-a", NFILE, fnm);
 286     bFreq2 = opt2parg_bSet("-skip2", asize(pa), pa);
 287     if (freq <= 0)
 288     {
 289         fprintf(stderr,
 290                 "The number of frames to skip is <= 0. "
 291                 "Writing out all frames.\n\n");
 292         freq = 1;
 293     }
 294     if (!bFreq2)
 295     {
 296         freq2 = freq;
 297     }
 298     else if (bFile2 && freq2 <= 0)
 299     {
 300         fprintf(stderr,
 301                 "The number of frames to skip in second trajectory is <= 0.\n"
 302                 "  Writing out all frames.\n\n");
 303         freq2 = 1;
 304     }
 305
 306     bPrev = (prev > 0);
 307     if (bPrev)
 308     {
 309         fprintf(stderr,
 310                 "WARNING: using option -prev with large trajectories will\n"
 311                 "         require a lot of memory and could lead to crashes\n");
 312         prev = abs(prev);
 313         if (freq != 1)
 314         {
 315             fprintf(stderr, "WARNING: option -skip also applies to -prev\n");
 316         }
 317     }
 318
 319     if (bFile2 && !bMat && !bBond)
 320     {
 321         fprintf(stderr,
 322                 "WARNING: second trajectory (-f2) useless when not calculating matrix (-m/-bm),\n"
 323                 "         will not read from %s\n",
 324                 opt2fn("-f2", NFILE, fnm));
 325         bFile2 = FALSE;
 326     }
 327
 328     if (bDelta)
 329     {
 330         bMat = TRUE;
 331         if (bFile2)
 332         {
 333             fprintf(stderr,
 334                     "WARNING: second trajectory (-f2) useless when making delta matrix,\n"
 335                     "         will not read from %s\n",
 336                     opt2fn("-f2", NFILE, fnm));
 337             bFile2 = FALSE;
 338         }
 339     }
 340
 341     bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, &xp, nullptr, box, TRUE);
 342     snew(w_rls, top.atoms.nr);
 343     snew(w_rms, top.atoms.nr);
 344
 345     if (!bTop && bBond)
 346     {
 347         fprintf(stderr,
 348                 "WARNING: Need a run input file for bond angle matrix,\n"
 349                 "         will not calculate bond angle matrix.\n");
 350         bBond = FALSE;
 351     }
 352
 353     if (bReset)
 354     {
 355         fprintf(stderr, "Select group for %s fit\n", bFit ? "least squares" : "translational");
 356         get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &ifit, &ind_fit, &gn_fit);
 357     }
 358     else
 359     {
 360         ifit = 0;
 361     }
 362
 363     if (bReset)
 364     {
 365         if (bFit && ifit < 3)
 366         {
 367             gmx_fatal(FARGS, "Need >= 3 points to fit!\n");
 368         }
 369
 370         bMass = FALSE;
 371         for (i = 0; i < ifit; i++)
 372         {
 373             if (bMassWeighted)
 374             {
 375                 w_rls[ind_fit[i]] = top.atoms.atom[ind_fit[i]].m;
 376             }
 377             else
 378             {
 379                 w_rls[ind_fit[i]] = 1;
 380             }
 381             bMass = bMass || (top.atoms.atom[ind_fit[i]].m != 0);
 382         }
 383         if (!bMass)
 384         {
 385             fprintf(stderr, "All masses in the fit group are 0, using masses of 1\n");
 386             for (i = 0; i < ifit; i++)
 387             {
 388                 w_rls[ind_fit[i]] = 1;
 389             }
 390         }
 391     }
 392
 393     if (bMat || bBond)
 394     {
 395         nrms = 1;
 396     }
 397
 398     snew(gn_rms, nrms);
 399     snew(ind_rms, nrms);
 400     snew(irms, nrms);
 401
 402     fprintf(stderr, "Select group%s for %s calculation\n", (nrms > 1) ? "s" : "", whatname[ewhat]);
 403     get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm), nrms, irms, ind_rms, gn_rms);
 404
 405     if (bNorm)
 406     {
 407         snew(rlsnorm, irms[0]);
 408     }
 409     snew(rls, nrms);
 410     for (j = 0; j < nrms; j++)
 411     {
 412         snew(rls[j], maxframe);
 413     }
 414     if (bMirror)
 415     {
 416         snew(rlsm, nrms);
 417         for (j = 0; j < nrms; j++)
 418         {
 419             snew(rlsm[j], maxframe);
 420         }
 421     }
 422     snew(time, maxframe);
 423     for (j = 0; j < nrms; j++)
 424     {
 425         bMass = FALSE;
 426         for (i = 0; i < irms[j]; i++)
 427         {
 428             if (bMassWeighted)
 429             {
 430                 w_rms[ind_rms[j][i]] = top.atoms.atom[ind_rms[j][i]].m;
 431             }
 432             else
 433             {
 434                 w_rms[ind_rms[j][i]] = 1.0;
 435             }
 436             bMass = bMass || (top.atoms.atom[ind_rms[j][i]].m != 0);
 437         }
 438         if (!bMass)
 439         {
 440             fprintf(stderr, "All masses in group %d are 0, using masses of 1\n", j);
 441             for (i = 0; i < irms[j]; i++)
 442             {
 443                 w_rms[ind_rms[j][i]] = 1;
 444             }
 445         }
 446     }
 447     /* Prepare reference frame */
 448     if (bPBC)
 449     {
 450         gpbc = gmx_rmpbc_init(&top.idef, pbcType, top.atoms.nr);
 451         gmx_rmpbc(gpbc, top.atoms.nr, box, xp);
 452     }
 453     if (bReset)
 454     {
 455         reset_x(ifit, ind_fit, top.atoms.nr, nullptr, xp, w_rls);
 456     }
 457     if (bMirror)
 458     {
 459         /* generate reference structure mirror image: */
 460         snew(xm, top.atoms.nr);
 461         for (i = 0; i < top.atoms.nr; i++)
 462         {
 463             copy_rvec(xp[i], xm[i]);
 464             xm[i][XX] = -xm[i][XX];
 465         }
 466     }
 467     if (ewhat == ewRhoSc)
 468     {
 469         norm_princ(&top.atoms, ifit, ind_fit, top.atoms.nr, xp);
 470     }
 471
 472     /* read first frame */
 473     natoms_trx = read_first_x(oenv, &status, opt2fn("-f", NFILE, fnm), &t, &x, box);
 474     if (natoms_trx != top.atoms.nr)
 475     {
 476         fprintf(stderr, "\nWARNING: topology has %d atoms, whereas trajectory has %d\n", top.atoms.nr, natoms_trx);
 477     }
 478     natoms = std::min(top.atoms.nr, natoms_trx);
 479     if (bMat || bBond || bPrev)
 480     {
 481         snew(mat_x, NFRAME);
 482
 483         if (bPrev)
 484         {
 485             /* With -prev we use all atoms for simplicity */
 486             n_ind_m = natoms;
 487         }
 488         else
 489         {
 490             /* Check which atoms we need (fit/rms) */
 491             snew(bInMat, natoms);
 492             for (i = 0; i < ifit; i++)
 493             {
 494                 bInMat[ind_fit[i]] = TRUE;
 495             }
 496             n_ind_m = ifit;
 497             for (i = 0; i < irms[0]; i++)
 498             {
 499                 if (!bInMat[ind_rms[0][i]])
 500                 {
 501                     bInMat[ind_rms[0][i]] = TRUE;
 502                     n_ind_m++;
 503                 }
 504             }
 505         }
 506         /* Make an index of needed atoms */
 507         snew(ind_m, n_ind_m);
 508         snew(rev_ind_m, natoms);
 509         j = 0;
 510         for (i = 0; i < natoms; i++)
 511         {
 512             if (bPrev || bInMat[i])
 513             {
 514                 ind_m[j]     = i;
 515                 rev_ind_m[i] = j;
 516                 j++;
 517             }
 518         }
 519         snew(w_rls_m, n_ind_m);
 520         snew(ind_rms_m, irms[0]);
 521         snew(w_rms_m, n_ind_m);
 522         for (i = 0; i < ifit; i++)
 523         {
 524             w_rls_m[rev_ind_m[ind_fit[i]]] = w_rls[ind_fit[i]];
 525         }
 526         for (i = 0; i < irms[0]; i++)
 527         {
 528             ind_rms_m[i]          = rev_ind_m[ind_rms[0][i]];
 529             w_rms_m[ind_rms_m[i]] = w_rms[ind_rms[0][i]];
 530         }
 531         sfree(bInMat);
 532     }
 533
 534     if (bBond)
 535     {
 536         ncons = 0;
 537         for (k = 0; k < F_NRE; k++)
 538         {
 539             if (IS_CHEMBOND(k))
 540             {
 541                 ncons += top.idef.il[k].nr / 3;
 542             }
 543         }
 544         fprintf(stderr, "Found %d bonds in topology\n", ncons);
 545         snew(ind_bond1, ncons);
 546         snew(ind_bond2, ncons);
 547         ibond = 0;
 548         for (k = 0; k < F_NRE; k++)
 549         {
 550             if (IS_CHEMBOND(k))
 551             {
 552                 iatom = top.idef.il[k].iatoms;
 553                 ncons = top.idef.il[k].nr / 3;
 554                 for (i = 0; i < ncons; i++)
 555                 {
 556                     bA1 = FALSE;
 557                     bA2 = FALSE;
 558                     for (j = 0; j < irms[0]; j++)
 559                     {
 560                         if (iatom[3 * i + 1] == ind_rms[0][j])
 561                         {
 562                             bA1 = TRUE;
 563                         }
 564                         if (iatom[3 * i + 2] == ind_rms[0][j])
 565                         {
 566                             bA2 = TRUE;
 567                         }
 568                     }
 569                     if (bA1 && bA2)
 570                     {
 571                         ind_bond1[ibond] = rev_ind_m[iatom[3 * i + 1]];
 572                         ind_bond2[ibond] = rev_ind_m[iatom[3 * i + 2]];
 573                         ibond++;
 574                     }
 575                 }
 576             }
 577         }
 578         fprintf(stderr, "Using %d bonds for bond angle matrix\n", ibond);
 579         if (ibond == 0)
 580         {
 581             gmx_fatal(FARGS, "0 bonds found");
 582         }
 583     }
 584
 585     /* start looping over frames: */
 586     int tel_mat = 0;
 587     int teller  = 0;
 588     int frame   = 0;
 589     do
 590     {
 591         if (bPBC)
 592         {
 593             gmx_rmpbc(gpbc, natoms, box, x);
 594         }
 595
 596         if (bReset)
 597         {
 598             reset_x(ifit, ind_fit, natoms, nullptr, x, w_rls);
 599         }
 600         if (ewhat == ewRhoSc)
 601         {
 602             norm_princ(&top.atoms, ifit, ind_fit, natoms, x);
 603         }
 604
 605         if (bFit)
 606         {
 607             /*do the least squares fit to original structure*/
 608             do_fit(natoms, w_rls, xp, x);
 609         }
 610
 611         if (frame % freq == 0)
 612         {
 613             /* keep frame for matrix calculation */
 614             if (bMat || bBond || bPrev)
 615             {
 616                 if (tel_mat >= NFRAME)
 617                 {
 618                     srenew(mat_x, tel_mat + 1);
 619                 }
 620                 snew(mat_x[tel_mat], n_ind_m);
 621                 for (i = 0; i < n_ind_m; i++)
 622                 {
 623                     copy_rvec(x[ind_m[i]], mat_x[tel_mat][i]);
 624                 }
 625             }
 626             tel_mat++;
 627
 628             /*calculate energy of root_least_squares*/
 629             if (bPrev)
 630             {
 631                 j = tel_mat - prev - 1;
 632                 if (j < 0)
 633                 {
 634                     j = 0;
 635                 }
 636                 for (i = 0; i < n_ind_m; i++)
 637                 {
 638                     copy_rvec(mat_x[j][i], xp[ind_m[i]]);
 639                 }
 640                 if (bReset)
 641                 {
 642                     reset_x(ifit, ind_fit, natoms, nullptr, xp, w_rls);
 643                 }
 644                 if (bFit)
 645                 {
 646                     do_fit(natoms, w_rls, x, xp);
 647                 }
 648             }
 649             for (j = 0; (j < nrms); j++)
 650             {
 651                 rls[j][teller] = calc_similar_ind(ewhat != ewRMSD, irms[j], ind_rms[j], w_rms, x, xp);
 652             }
 653             if (bNorm)
 654             {
 655                 for (j = 0; (j < irms[0]); j++)
 656                 {
 657                     rlsnorm[j] += calc_similar_ind(ewhat != ewRMSD, 1, &(ind_rms[0][j]), w_rms, x, xp);
 658                 }
 659             }
 660
 661             if (bMirror)
 662             {
 663                 if (bFit)
 664                 {
 665                     /*do the least squares fit to mirror of original structure*/
 666                     do_fit(natoms, w_rls, xm, x);
 667                 }
 668
 669                 for (j = 0; j < nrms; j++)
 670                 {
 671                     rlsm[j][teller] =
 672                             calc_similar_ind(ewhat != ewRMSD, irms[j], ind_rms[j], w_rms, x, xm);
 673                 }
 674             }
 675             time[teller] = output_env_conv_time(oenv, t);
 676
 677             teller++;
 678         }
 679         frame++;
 680         if (teller >= maxframe)
 681         {
 682             maxframe += NFRAME;
 683             srenew(time, maxframe);
 684             for (j = 0; (j < nrms); j++)
 685             {
 686                 srenew(rls[j], maxframe);
 687             }
 688             if (bMirror)
 689             {
 690                 for (j = 0; (j < nrms); j++)
 691                 {
 692                     srenew(rlsm[j], maxframe);
 693                 }
 694             }
 695         }
 696     } while (read_next_x(oenv, status, &t, x, box));
 697     close_trx(status);
 698
 699     int tel_mat2 = 0;
 700     int teller2  = 0;
 701     int frame2   = 0;
 702     if (bFile2)
 703     {
 704         snew(time2, maxframe2);
 705
 706         fprintf(stderr, "\nWill read second trajectory file\n");
 707         snew(mat_x2, NFRAME);
 708         natoms_trx2 = read_first_x(oenv, &status, opt2fn("-f2", NFILE, fnm), &t, &x, box);
 709         if (natoms_trx2 != natoms_trx)
 710         {
 711             gmx_fatal(FARGS,
 712                       "Second trajectory (%d atoms) does not match the first one"
 713                       " (%d atoms)",
 714                       natoms_trx2,
 715                       natoms_trx);
 716         }
 717         frame2 = 0;
 718         do
 719         {
 720             if (bPBC)
 721             {
 722                 gmx_rmpbc(gpbc, natoms, box, x);
 723             }
 724
 725             if (bReset)
 726             {
 727                 reset_x(ifit, ind_fit, natoms, nullptr, x, w_rls);
 728             }
 729             if (ewhat == ewRhoSc)
 730             {
 731                 norm_princ(&top.atoms, ifit, ind_fit, natoms, x);
 732             }
 733
 734             if (bFit)
 735             {
 736                 /*do the least squares fit to original structure*/
 737                 do_fit(natoms, w_rls, xp, x);
 738             }
 739
 740             if (frame2 % freq2 == 0)
 741             {
 742                 /* keep frame for matrix calculation */
 743                 if (bMat)
 744                 {
 745                     if (tel_mat2 >= NFRAME)
 746                     {
 747                         srenew(mat_x2, tel_mat2 + 1);
 748                     }
 749                     snew(mat_x2[tel_mat2], n_ind_m);
 750                     for (i = 0; i < n_ind_m; i++)
 751                     {
 752                         copy_rvec(x[ind_m[i]], mat_x2[tel_mat2][i]);
 753                     }
 754                 }
 755                 tel_mat2++;
 756
 757                 time2[teller2] = output_env_conv_time(oenv, t);
 758
 759                 teller2++;
 760             }
 761             frame2++;
 762             if (teller2 >= maxframe2)
 763             {
 764                 maxframe2 += NFRAME;
 765                 srenew(time2, maxframe2);
 766             }
 767         } while (read_next_x(oenv, status, &t, x, box));
 768         close_trx(status);
 769     }
 770     else
 771     {
 772         mat_x2   = mat_x;
 773         time2    = time;
 774         tel_mat2 = tel_mat;
 775         freq2    = freq;
 776     }
 777     gmx_rmpbc_done(gpbc);
 778
 779     if (bMat || bBond)
 780     {
 781         /* calculate RMS matrix */
 782         fprintf(stderr, "\n");
 783         if (bMat)
 784         {
 785             fprintf(stderr, "Building %s matrix, %dx%d elements\n", whatname[ewhat], tel_mat, tel_mat2);
 786             snew(rmsd_mat, tel_mat);
 787         }
 788         if (bBond)
 789         {
 790             fprintf(stderr, "Building bond angle matrix, %dx%d elements\n", tel_mat, tel_mat2);
 791             snew(bond_mat, tel_mat);
 792         }
 793         snew(axis, tel_mat);
 794         snew(axis2, tel_mat2);
 795         rmsd_max = 0;
 796         if (bFile2)
 797         {
 798             rmsd_min = 1e10;
 799         }
 800         else
 801         {
 802             rmsd_min = 0;
 803         }
 804         rmsd_avg = 0;
 805         bond_max = 0;
 806         bond_min = 1e10;
 807         for (j = 0; j < tel_mat2; j++)
 808         {
 809             axis2[j] = time2[freq2 * j];
 810         }
 811         if (bDelta)
 812         {
 813             if (bDeltaLog)
 814             {
 815                 delta_scalex = 8.0 / std::log(2.0);
 816                 delta_xsize  = gmx::roundToInt(std::log(tel_mat / 2.) * delta_scalex) + 1;
 817             }
 818             else
 819             {
 820                 delta_xsize = tel_mat / 2;
 821             }
 822             delta_scaley = 1.0 / delta_maxy;
 823             snew(delta, delta_xsize);
 824             for (j = 0; j < delta_xsize; j++)
 825             {
 826                 snew(delta[j], del_lev + 1);
 827             }
 828             if (avl > 0)
 829             {
 830                 snew(rmsdav_mat, tel_mat);
 831                 for (j = 0; j < tel_mat; j++)
 832                 {
 833                     snew(rmsdav_mat[j], tel_mat);
 834                 }
 835             }
 836         }
 837
 838         if (bFitAll)
 839         {
 840             snew(mat_x2_j, natoms);
 841         }
 842         for (i = 0; i < tel_mat; i++)
 843         {
 844             axis[i] = time[freq * i];
 845             fprintf(stderr, "\r element %5d; time %5.2f  ", i, axis[i]);
 846             fflush(stderr);
 847             if (bMat)
 848             {
 849                 snew(rmsd_mat[i], tel_mat2);
 850             }
 851             if (bBond)
 852             {
 853                 snew(bond_mat[i], tel_mat2);
 854             }
 855             for (j = 0; j < tel_mat2; j++)
 856             {
 857                 if (bFitAll)
 858                 {
 859                     for (k = 0; k < n_ind_m; k++)
 860                     {
 861                         copy_rvec(mat_x2[j][k], mat_x2_j[k]);
 862                     }
 863                     do_fit(n_ind_m, w_rls_m, mat_x[i], mat_x2_j);
 864                 }
 865                 else
 866                 {
 867                     mat_x2_j = mat_x2[j];
 868                 }
 869                 if (bMat)
 870                 {
 871                     if (bFile2 || (i < j))
 872                     {
 873                         rmsd_mat[i][j] = calc_similar_ind(
 874                                 ewhat != ewRMSD, irms[0], ind_rms_m, w_rms_m, mat_x[i], mat_x2_j);
 875                         if (rmsd_mat[i][j] > rmsd_max)
 876                         {
 877                             rmsd_max = rmsd_mat[i][j];
 878                         }
 879                         if (rmsd_mat[i][j] < rmsd_min)
 880                         {
 881                             rmsd_min = rmsd_mat[i][j];
 882                         }
 883                         rmsd_avg += rmsd_mat[i][j];
 884                     }
 885                     else
 886                     {
 887                         rmsd_mat[i][j] = rmsd_mat[j][i];
 888                     }
 889                 }
 890                 if (bBond)
 891                 {
 892                     if (bFile2 || (i <= j))
 893                     {
 894                         ang = 0.0;
 895                         for (m = 0; m < ibond; m++)
 896                         {
 897                             rvec_sub(mat_x[i][ind_bond1[m]], mat_x[i][ind_bond2[m]], vec1);
 898                             rvec_sub(mat_x2_j[ind_bond1[m]], mat_x2_j[ind_bond2[m]], vec2);
 899                             ang += std::acos(cos_angle(vec1, vec2));
 900                         }
 901                         bond_mat[i][j] = ang * 180.0 / (M_PI * ibond);
 902                         if (bond_mat[i][j] > bond_max)
 903                         {
 904                             bond_max = bond_mat[i][j];
 905                         }
 906                         if (bond_mat[i][j] < bond_min)
 907                         {
 908                             bond_min = bond_mat[i][j];
 909                         }
 910                     }
 911                     else
 912                     {
 913                         bond_mat[i][j] = bond_mat[j][i];
 914                     }
 915                 }
 916             }
 917         }
 918         if (bFile2)
 919         {
 920             rmsd_avg /= static_cast<real>(tel_mat) * static_cast<real>(tel_mat2);
 921         }
 922         else
 923         {
 924             rmsd_avg /= tel_mat * (tel_mat - 1) / 2.;
 925         }
 926         if (bMat && (avl > 0))
 927         {
 928             rmsd_max = 0.0;
 929             rmsd_min = 0.0;
 930             rmsd_avg = 0.0;
 931             for (j = 0; j < tel_mat - 1; j++)
 932             {
 933                 for (i = j + 1; i < tel_mat; i++)
 934                 {
 935                     av_tot     = 0;
 936                     weight_tot = 0;
 937                     for (my = -avl; my <= avl; my++)
 938                     {
 939                         if ((j + my >= 0) && (j + my < tel_mat))
 940                         {
 941                             abs_my = std::abs(my);
 942                             for (mx = -avl; mx <= avl; mx++)
 943                             {
 944                                 if ((i + mx >= 0) && (i + mx < tel_mat))
 945                                 {
 946                                     weight = avl + 1.0 - std::max(std::abs(mx), abs_my);
 947                                     av_tot += weight * rmsd_mat[i + mx][j + my];
 948                                     weight_tot += weight;
 949                                 }
 950                             }
 951                         }
 952                     }
 953                     rmsdav_mat[i][j] = av_tot / weight_tot;
 954                     rmsdav_mat[j][i] = rmsdav_mat[i][j];
 955                     if (rmsdav_mat[i][j] > rmsd_max)
 956                     {
 957                         rmsd_max = rmsdav_mat[i][j];
 958                     }
 959                 }
 960             }
 961             rmsd_mat = rmsdav_mat;
 962         }
 963
 964         if (bMat)
 965         {
 966             fprintf(stderr, "\n%s: Min %f, Max %f, Avg %f\n", whatname[ewhat], rmsd_min, rmsd_max, rmsd_avg);
 967             rlo.r = 1;
 968             rlo.g = 1;
 969             rlo.b = 1;
 970             rhi.r = 0;
 971             rhi.g = 0;
 972             rhi.b = 0;
 973             if (rmsd_user_max != -1)
 974             {
 975                 rmsd_max = rmsd_user_max;
 976             }
 977             if (rmsd_user_min != -1)
 978             {
 979                 rmsd_min = rmsd_user_min;
 980             }
 981             if ((rmsd_user_max != -1) || (rmsd_user_min != -1))
 982             {
 983                 fprintf(stderr, "Min and Max value set to resp. %f and %f\n", rmsd_min, rmsd_max);
 984             }
 985             sprintf(buf, "%s %s matrix", gn_rms[0], whatname[ewhat]);
 986             write_xpm(opt2FILE("-m", NFILE, fnm, "w"),
 987                       0,
 988                       buf,
 989                       whatlabel[ewhat],
 990                       output_env_get_time_label(oenv),
 991                       output_env_get_time_label(oenv),
 992                       tel_mat,
 993                       tel_mat2,
 994                       axis,
 995                       axis2,
 996                       rmsd_mat,
 997                       rmsd_min,
 998                       rmsd_max,
 999                       rlo,
1000                       rhi,
1001                       &nlevels);
1002             /* Print the distribution of RMSD values */
1003             if (opt2bSet("-dist", NFILE, fnm))
1004             {
1005                 low_rmsd_dist(opt2fn("-dist", NFILE, fnm), rmsd_max, tel_mat, rmsd_mat, oenv);
1006             }
1007
1008             if (bDelta)
1009             {
1010                 snew(delta_tot, delta_xsize);
1011                 for (j = 0; j < tel_mat - 1; j++)
1012                 {
1013                     for (i = j + 1; i < tel_mat; i++)
1014                     {
1015                         mx = i - j;
1016                         if (mx < tel_mat / 2)
1017                         {
1018                             if (bDeltaLog)
1019                             {
1020                                 mx = gmx::roundToInt(std::log(static_cast<real>(mx)) * delta_scalex);
1021                             }
1022                             my = gmx::roundToInt(rmsd_mat[i][j] * delta_scaley * static_cast<real>(del_lev));
1023                             delta_tot[mx] += 1.0;
1024                             if ((rmsd_mat[i][j] >= 0) && (rmsd_mat[i][j] <= delta_maxy))
1025                             {
1026                                 delta[mx][my] += 1.0;
1027                             }
1028                         }
1029                     }
1030                 }
1031                 delta_max = 0;
1032                 for (i = 0; i < delta_xsize; i++)
1033                 {
1034                     if (delta_tot[i] > 0.0)
1035                     {
1036                         delta_tot[i] = 1.0 / delta_tot[i];
1037                         for (j = 0; j <= del_lev; j++)
1038                         {
1039                             delta[i][j] *= delta_tot[i];
1040                             if (delta[i][j] > delta_max)
1041                             {
1042                                 delta_max = delta[i][j];
1043                             }
1044                         }
1045                     }
1046                 }
1047                 fprintf(stderr, "Maximum in delta matrix: %f\n", delta_max);
1048                 snew(del_xaxis, delta_xsize);
1049                 snew(del_yaxis, del_lev + 1);
1050                 for (i = 0; i < delta_xsize; i++)
1051                 {
1052                     del_xaxis[i] = axis[i] - axis[0];
1053                 }
1054                 for (i = 0; i < del_lev + 1; i++)
1055                 {
1056                     del_yaxis[i] = delta_maxy * static_cast<real>(i) / static_cast<real>(del_lev);
1057                 }
1058                 sprintf(buf, "%s %s vs. delta t", gn_rms[0], whatname[ewhat]);
1059                 fp = gmx_ffopen("delta.xpm", "w");
1060                 write_xpm(fp,
1061                           0,
1062                           buf,
1063                           "density",
1064                           output_env_get_time_label(oenv),
1065                           whatlabel[ewhat],
1066                           delta_xsize,
1067                           del_lev + 1,
1068                           del_xaxis,
1069                           del_yaxis,
1070                           delta,
1071                           0.0,
1072                           delta_max,
1073                           rlo,
1074                           rhi,
1075                           &nlevels);
1076                 gmx_ffclose(fp);
1077             }
1078             if (opt2bSet("-bin", NFILE, fnm))
1079             {
1080                 /* NB: File must be binary if we use fwrite */
1081                 fp = ftp2FILE(efDAT, NFILE, fnm, "wb");
1082                 for (i = 0; i < tel_mat; i++)
1083                 {
1084                     if (static_cast<int>(fwrite(rmsd_mat[i], sizeof(**rmsd_mat), tel_mat2, fp)) != tel_mat2)
1085                     {
1086                         gmx_fatal(FARGS, "Error writing to output file");
1087                     }
1088                 }
1089                 gmx_ffclose(fp);
1090             }
1091         }
1092         if (bBond)
1093         {
1094             fprintf(stderr, "\nMin. angle: %f, Max. angle: %f\n", bond_min, bond_max);
1095             if (bond_user_max != -1)
1096             {
1097                 bond_max = bond_user_max;
1098             }
1099             if (bond_user_min != -1)
1100             {
1101                 bond_min = bond_user_min;
1102             }
1103             if ((bond_user_max != -1) || (bond_user_min != -1))
1104             {
1105                 fprintf(stderr,
1106                         "Bond angle Min and Max set to:\n"
1107                         "Min. angle: %f, Max. angle: %f\n",
1108                         bond_min,
1109                         bond_max);
1110             }
1111             rlo.r = 1;
1112             rlo.g = 1;
1113             rlo.b = 1;
1114             rhi.r = 0;
1115             rhi.g = 0;
1116             rhi.b = 0;
1117             sprintf(buf, "%s av. bond angle deviation", gn_rms[0]);
1118             write_xpm(opt2FILE("-bm", NFILE, fnm, "w"),
1119                       0,
1120                       buf,
1121                       "degrees",
1122                       output_env_get_time_label(oenv),
1123                       output_env_get_time_label(oenv),
1124                       tel_mat,
1125                       tel_mat2,
1126                       axis,
1127                       axis2,
1128                       bond_mat,
1129                       bond_min,
1130                       bond_max,
1131                       rlo,
1132                       rhi,
1133                       &nlevels);
1134         }
1135     }
1136
1137     bAv = opt2bSet("-a", NFILE, fnm);
1138
1139     /* Write the RMSD's to file */
1140     if (!bPrev)
1141     {
1142         sprintf(buf, "%s", whatxvgname[ewhat]);
1143     }
1144     else
1145     {
1146         sprintf(buf,
1147                 "%s with frame %g %s ago",
1148                 whatxvgname[ewhat],
1149                 time[prev * freq] - time[0],
1150                 output_env_get_time_label(oenv).c_str());
1151     }
1152     fp = xvgropen(opt2fn("-o", NFILE, fnm), buf, output_env_get_xvgr_tlabel(oenv), whatxvglabel[ewhat], oenv);
1153     if (output_env_get_print_xvgr_codes(oenv))
1154     {
1155         fprintf(fp,
1156                 "@ subtitle \"%s%s after %s%s%s\"\n",
1157                 (nrms == 1) ? "" : "of ",
1158                 gn_rms[0],
1159                 fitgraphlabel[efit],
1160                 bFit ? " to " : "",
1161                 bFit ? gn_fit : "");
1162     }
1163     if (nrms != 1)
1164     {
1165         xvgr_legend(fp, nrms, gn_rms, oenv);
1166     }
1167     for (i = 0; (i < teller); i++)
1168     {
1169         if (bSplit && i > 0 && std::abs(time[bPrev ? freq * i : i] / output_env_get_time_factor(oenv)) < 1e-5)
1170         {
1171             fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
1172         }
1173         fprintf(fp, "%12.7f", time[bPrev ? freq * i : i]);
1174         for (j = 0; (j < nrms); j++)
1175         {
1176             fprintf(fp, " %12.7f", rls[j][i]);
1177             if (bAv)
1178             {
1179                 rlstot += rls[j][i];
1180             }
1181         }
1182         fprintf(fp, "\n");
1183     }
1184     xvgrclose(fp);
1185
1186     if (bMirror)
1187     {
1188         /* Write the mirror RMSD's to file */
1189         sprintf(buf, "%s with Mirror", whatxvgname[ewhat]);
1190         sprintf(buf2, "Mirror %s", whatxvglabel[ewhat]);
1191         fp = xvgropen(opt2fn("-mir", NFILE, fnm), buf, output_env_get_xvgr_tlabel(oenv), buf2, oenv);
1192         if (nrms == 1)
1193         {
1194             if (output_env_get_print_xvgr_codes(oenv))
1195             {
1196                 fprintf(fp, "@ subtitle \"of %s after lsq fit to mirror of %s\"\n", gn_rms[0], bFit ? gn_fit : "");
1197             }
1198         }
1199         else
1200         {
1201             if (output_env_get_print_xvgr_codes(oenv))
1202             {
1203                 fprintf(fp, "@ subtitle \"after lsq fit to mirror %s\"\n", bFit ? gn_fit : "");
1204             }
1205             xvgr_legend(fp, nrms, gn_rms, oenv);
1206         }
1207         for (i = 0; (i < teller); i++)
1208         {
1209             if (bSplit && i > 0 && std::abs(time[i]) < 1e-5)
1210             {
1211                 fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
1212             }
1213             fprintf(fp, "%12.7f", time[i]);
1214             for (j = 0; (j < nrms); j++)
1215             {
1216                 fprintf(fp, " %12.7f", rlsm[j][i]);
1217             }
1218             fprintf(fp, "\n");
1219         }
1220         xvgrclose(fp);
1221     }
1222
1223     if (bAv)
1224     {
1225         sprintf(buf, "Average %s", whatxvgname[ewhat]);
1226         sprintf(buf2, "Average %s", whatxvglabel[ewhat]);
1227         fp = xvgropen(opt2fn("-a", NFILE, fnm), buf, "Residue", buf2, oenv);
1228         for (j = 0; (j < nrms); j++)
1229         {
1230             fprintf(fp, "%10d  %10g\n", j, rlstot / static_cast<real>(teller));
1231         }
1232         xvgrclose(fp);
1233     }
1234
1235     if (bNorm)
1236     {
1237         fp = xvgropen("aver.xvg", gn_rms[0], "Residue", whatxvglabel[ewhat], oenv);
1238         for (j = 0; (j < irms[0]); j++)
1239         {
1240             fprintf(fp, "%10d  %10g\n", j, rlsnorm[j] / static_cast<real>(teller));
1241         }
1242         xvgrclose(fp);
1243     }
1244     do_view(oenv, opt2fn_null("-a", NFILE, fnm), "-graphtype bar");
1245     do_view(oenv, opt2fn("-o", NFILE, fnm), nullptr);
1246     do_view(oenv, opt2fn_null("-mir", NFILE, fnm), nullptr);
1247     do_view(oenv, opt2fn_null("-m", NFILE, fnm), nullptr);
1248     do_view(oenv, opt2fn_null("-bm", NFILE, fnm), nullptr);
1249     do_view(oenv, opt2fn_null("-dist", NFILE, fnm), nullptr);
1250
1251     return 0;
1252 }