src/gromacs/gmxana/gmx_rms.cpp

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