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