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