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, 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] / 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, teller, teller2, tel_mat, tel_mat2;
 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     tel_mat = 0;
 614     teller  = 0;
 615     do
 616     {
 617         if (bPBC)
 618         {
 619             gmx_rmpbc(gpbc, natoms, box, x);
 620         }
 621
 622         if (bReset)
 623         {
 624             reset_x(ifit, ind_fit, natoms, nullptr, x, w_rls);
 625         }
 626         if (ewhat == ewRhoSc)
 627         {
 628             norm_princ(&top.atoms, ifit, ind_fit, natoms, x);
 629         }
 630
 631         if (bFit)
 632         {
 633             /*do the least squares fit to original structure*/
 634             do_fit(natoms, w_rls, xp, x);
 635         }
 636
 637         if (teller % freq == 0)
 638         {
 639             /* keep frame for matrix calculation */
 640             if (bMat || bBond || bPrev)
 641             {
 642                 if (tel_mat >= NFRAME)
 643                 {
 644                     srenew(mat_x, tel_mat+1);
 645                 }
 646                 snew(mat_x[tel_mat], n_ind_m);
 647                 for (i = 0; i < n_ind_m; i++)
 648                 {
 649                     copy_rvec(x[ind_m[i]], mat_x[tel_mat][i]);
 650                 }
 651             }
 652             tel_mat++;
 653         }
 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         if (teller >= maxframe)
 708         {
 709             maxframe += NFRAME;
 710             srenew(time, maxframe);
 711             for (j = 0; (j < nrms); j++)
 712             {
 713                 srenew(rls[j], maxframe);
 714             }
 715             if (bMirror)
 716             {
 717                 for (j = 0; (j < nrms); j++)
 718                 {
 719                     srenew(rlsm[j], maxframe);
 720                 }
 721             }
 722         }
 723     }
 724     while (read_next_x(oenv, status, &t, x, box));
 725     close_trj(status);
 726
 727     if (bFile2)
 728     {
 729         snew(time2, maxframe2);
 730
 731         fprintf(stderr, "\nWill read second trajectory file\n");
 732         snew(mat_x2, NFRAME);
 733         natoms_trx2 =
 734             read_first_x(oenv, &status, opt2fn("-f2", NFILE, fnm), &t, &x, box);
 735         if (natoms_trx2 != natoms_trx)
 736         {
 737             gmx_fatal(FARGS,
 738                       "Second trajectory (%d atoms) does not match the first one"
 739                       " (%d atoms)", natoms_trx2, natoms_trx);
 740         }
 741         tel_mat2 = 0;
 742         teller2  = 0;
 743         do
 744         {
 745             if (bPBC)
 746             {
 747                 gmx_rmpbc(gpbc, natoms, box, x);
 748             }
 749
 750             if (bReset)
 751             {
 752                 reset_x(ifit, ind_fit, natoms, nullptr, x, w_rls);
 753             }
 754             if (ewhat == ewRhoSc)
 755             {
 756                 norm_princ(&top.atoms, ifit, ind_fit, natoms, x);
 757             }
 758
 759             if (bFit)
 760             {
 761                 /*do the least squares fit to original structure*/
 762                 do_fit(natoms, w_rls, xp, x);
 763             }
 764
 765             if (teller2 % freq2 == 0)
 766             {
 767                 /* keep frame for matrix calculation */
 768                 if (bMat)
 769                 {
 770                     if (tel_mat2 >= NFRAME)
 771                     {
 772                         srenew(mat_x2, tel_mat2+1);
 773                     }
 774                     snew(mat_x2[tel_mat2], n_ind_m);
 775                     for (i = 0; i < n_ind_m; i++)
 776                     {
 777                         copy_rvec(x[ind_m[i]], mat_x2[tel_mat2][i]);
 778                     }
 779                 }
 780                 tel_mat2++;
 781             }
 782
 783             time2[teller2] = output_env_conv_time(oenv, t);
 784
 785             teller2++;
 786             if (teller2 >= maxframe2)
 787             {
 788                 maxframe2 += NFRAME;
 789                 srenew(time2, maxframe2);
 790             }
 791         }
 792         while (read_next_x(oenv, status, &t, x, box));
 793         close_trj(status);
 794     }
 795     else
 796     {
 797         mat_x2   = mat_x;
 798         time2    = time;
 799         tel_mat2 = tel_mat;
 800         freq2    = freq;
 801     }
 802     gmx_rmpbc_done(gpbc);
 803
 804     if (bMat || bBond)
 805     {
 806         /* calculate RMS matrix */
 807         fprintf(stderr, "\n");
 808         if (bMat)
 809         {
 810             fprintf(stderr, "Building %s matrix, %dx%d elements\n",
 811                     whatname[ewhat], tel_mat, tel_mat2);
 812             snew(rmsd_mat, tel_mat);
 813         }
 814         if (bBond)
 815         {
 816             fprintf(stderr, "Building bond angle matrix, %dx%d elements\n",
 817                     tel_mat, tel_mat2);
 818             snew(bond_mat, tel_mat);
 819         }
 820         snew(axis, tel_mat);
 821         snew(axis2, tel_mat2);
 822         rmsd_max = 0;
 823         if (bFile2)
 824         {
 825             rmsd_min = 1e10;
 826         }
 827         else
 828         {
 829             rmsd_min = 0;
 830         }
 831         rmsd_avg = 0;
 832         bond_max = 0;
 833         bond_min = 1e10;
 834         for (j = 0; j < tel_mat2; j++)
 835         {
 836             axis2[j] = time2[freq2*j];
 837         }
 838         if (bDelta)
 839         {
 840             if (bDeltaLog)
 841             {
 842                 delta_scalex = 8.0/std::log(2.0);
 843                 delta_xsize  = static_cast<int>(std::log(static_cast<real>(tel_mat/2))*delta_scalex+0.5)+1;
 844             }
 845             else
 846             {
 847                 delta_xsize = tel_mat/2;
 848             }
 849             delta_scaley = 1.0/delta_maxy;
 850             snew(delta, delta_xsize);
 851             for (j = 0; j < delta_xsize; j++)
 852             {
 853                 snew(delta[j], del_lev+1);
 854             }
 855             if (avl > 0)
 856             {
 857                 snew(rmsdav_mat, tel_mat);
 858                 for (j = 0; j < tel_mat; j++)
 859                 {
 860                     snew(rmsdav_mat[j], tel_mat);
 861                 }
 862             }
 863         }
 864
 865         if (bFitAll)
 866         {
 867             snew(mat_x2_j, natoms);
 868         }
 869         for (i = 0; i < tel_mat; i++)
 870         {
 871             axis[i] = time[freq*i];
 872             fprintf(stderr, "\r element %5d; time %5.2f  ", i, axis[i]);
 873             fflush(stderr);
 874             if (bMat)
 875             {
 876                 snew(rmsd_mat[i], tel_mat2);
 877             }
 878             if (bBond)
 879             {
 880                 snew(bond_mat[i], tel_mat2);
 881             }
 882             for (j = 0; j < tel_mat2; j++)
 883             {
 884                 if (bFitAll)
 885                 {
 886                     for (k = 0; k < n_ind_m; k++)
 887                     {
 888                         copy_rvec(mat_x2[j][k], mat_x2_j[k]);
 889                     }
 890                     do_fit(n_ind_m, w_rls_m, mat_x[i], mat_x2_j);
 891                 }
 892                 else
 893                 {
 894                     mat_x2_j = mat_x2[j];
 895                 }
 896                 if (bMat)
 897                 {
 898                     if (bFile2 || (i < j))
 899                     {
 900                         rmsd_mat[i][j] =
 901                             calc_similar_ind(ewhat != ewRMSD, irms[0], ind_rms_m,
 902                                              w_rms_m, mat_x[i], mat_x2_j);
 903                         if (rmsd_mat[i][j] > rmsd_max)
 904                         {
 905                             rmsd_max = rmsd_mat[i][j];
 906                         }
 907                         if (rmsd_mat[i][j] < rmsd_min)
 908                         {
 909                             rmsd_min = rmsd_mat[i][j];
 910                         }
 911                         rmsd_avg += rmsd_mat[i][j];
 912                     }
 913                     else
 914                     {
 915                         rmsd_mat[i][j] = rmsd_mat[j][i];
 916                     }
 917                 }
 918                 if (bBond)
 919                 {
 920                     if (bFile2 || (i <= j))
 921                     {
 922                         ang = 0.0;
 923                         for (m = 0; m < ibond; m++)
 924                         {
 925                             rvec_sub(mat_x[i][ind_bond1[m]], mat_x[i][ind_bond2[m]], vec1);
 926                             rvec_sub(mat_x2_j[ind_bond1[m]], mat_x2_j[ind_bond2[m]], vec2);
 927                             ang += std::acos(cos_angle(vec1, vec2));
 928                         }
 929                         bond_mat[i][j] = ang*180.0/(M_PI*ibond);
 930                         if (bond_mat[i][j] > bond_max)
 931                         {
 932                             bond_max = bond_mat[i][j];
 933                         }
 934                         if (bond_mat[i][j] < bond_min)
 935                         {
 936                             bond_min = bond_mat[i][j];
 937                         }
 938                     }
 939                     else
 940                     {
 941                         bond_mat[i][j] = bond_mat[j][i];
 942                     }
 943                 }
 944             }
 945         }
 946         if (bFile2)
 947         {
 948             rmsd_avg /= tel_mat*tel_mat2;
 949         }
 950         else
 951         {
 952             rmsd_avg /= tel_mat*(tel_mat - 1)/2;
 953         }
 954         if (bMat && (avl > 0))
 955         {
 956             rmsd_max = 0.0;
 957             rmsd_min = 0.0;
 958             rmsd_avg = 0.0;
 959             for (j = 0; j < tel_mat-1; j++)
 960             {
 961                 for (i = j+1; i < tel_mat; i++)
 962                 {
 963                     av_tot     = 0;
 964                     weight_tot = 0;
 965                     for (my = -avl; my <= avl; my++)
 966                     {
 967                         if ((j+my >= 0) && (j+my < tel_mat))
 968                         {
 969                             abs_my = std::abs(my);
 970                             for (mx = -avl; mx <= avl; mx++)
 971                             {
 972                                 if ((i+mx >= 0) && (i+mx < tel_mat))
 973                                 {
 974                                     weight      = avl+1.0-std::max(std::abs(mx), abs_my);
 975                                     av_tot     += weight*rmsd_mat[i+mx][j+my];
 976                                     weight_tot += weight;
 977                                 }
 978                             }
 979                         }
 980                     }
 981                     rmsdav_mat[i][j] = av_tot/weight_tot;
 982                     rmsdav_mat[j][i] = rmsdav_mat[i][j];
 983                     if (rmsdav_mat[i][j] > rmsd_max)
 984                     {
 985                         rmsd_max = rmsdav_mat[i][j];
 986                     }
 987                 }
 988             }
 989             rmsd_mat = rmsdav_mat;
 990         }
 991
 992         if (bMat)
 993         {
 994             fprintf(stderr, "\n%s: Min %f, Max %f, Avg %f\n",
 995                     whatname[ewhat], rmsd_min, rmsd_max, rmsd_avg);
 996             rlo.r = 1; rlo.g = 1; rlo.b = 1;
 997             rhi.r = 0; rhi.g = 0; rhi.b = 0;
 998             if (rmsd_user_max != -1)
 999             {
1000                 rmsd_max = rmsd_user_max;
1001             }
1002             if (rmsd_user_min != -1)
1003             {
1004                 rmsd_min = rmsd_user_min;
1005             }
1006             if ((rmsd_user_max !=  -1) || (rmsd_user_min != -1))
1007             {
1008                 fprintf(stderr, "Min and Max value set to resp. %f and %f\n",
1009                         rmsd_min, rmsd_max);
1010             }
1011             sprintf(buf, "%s %s matrix", gn_rms[0], whatname[ewhat]);
1012             write_xpm(opt2FILE("-m", NFILE, fnm, "w"), 0, buf, whatlabel[ewhat],
1013                       output_env_get_time_label(oenv), output_env_get_time_label(oenv), tel_mat, tel_mat2,
1014                       axis, axis2, rmsd_mat, rmsd_min, rmsd_max, rlo, rhi, &nlevels);
1015             /* Print the distribution of RMSD values */
1016             if (opt2bSet("-dist", NFILE, fnm))
1017             {
1018                 low_rmsd_dist(opt2fn("-dist", NFILE, fnm), rmsd_max, tel_mat, rmsd_mat, oenv);
1019             }
1020
1021             if (bDelta)
1022             {
1023                 snew(delta_tot, delta_xsize);
1024                 for (j = 0; j < tel_mat-1; j++)
1025                 {
1026                     for (i = j+1; i < tel_mat; i++)
1027                     {
1028                         mx = i-j;
1029                         if (mx < tel_mat/2)
1030                         {
1031                             if (bDeltaLog)
1032                             {
1033                                 mx = static_cast<int>(std::log(static_cast<real>(mx))*delta_scalex+0.5);
1034                             }
1035                             my             = static_cast<int>(rmsd_mat[i][j]*delta_scaley*del_lev+0.5);
1036                             delta_tot[mx] += 1.0;
1037                             if ((rmsd_mat[i][j] >= 0) && (rmsd_mat[i][j] <= delta_maxy))
1038                             {
1039                                 delta[mx][my] += 1.0;
1040                             }
1041                         }
1042                     }
1043                 }
1044                 delta_max = 0;
1045                 for (i = 0; i < delta_xsize; i++)
1046                 {
1047                     if (delta_tot[i] > 0.0)
1048                     {
1049                         delta_tot[i] = 1.0/delta_tot[i];
1050                         for (j = 0; j <= del_lev; j++)
1051                         {
1052                             delta[i][j] *= delta_tot[i];
1053                             if (delta[i][j] > delta_max)
1054                             {
1055                                 delta_max = delta[i][j];
1056                             }
1057                         }
1058                     }
1059                 }
1060                 fprintf(stderr, "Maximum in delta matrix: %f\n", delta_max);
1061                 snew(del_xaxis, delta_xsize);
1062                 snew(del_yaxis, del_lev+1);
1063                 for (i = 0; i < delta_xsize; i++)
1064                 {
1065                     del_xaxis[i] = axis[i]-axis[0];
1066                 }
1067                 for (i = 0; i < del_lev+1; i++)
1068                 {
1069                     del_yaxis[i] = delta_maxy*i/del_lev;
1070                 }
1071                 sprintf(buf, "%s %s vs. delta t", gn_rms[0], whatname[ewhat]);
1072                 fp = gmx_ffopen("delta.xpm", "w");
1073                 write_xpm(fp, 0, buf, "density", output_env_get_time_label(oenv), whatlabel[ewhat],
1074                           delta_xsize, del_lev+1, del_xaxis, del_yaxis,
1075                           delta, 0.0, delta_max, rlo, rhi, &nlevels);
1076                 gmx_ffclose(fp);
1077             }
1078             if (opt2bSet("-bin", NFILE, fnm))
1079             {
1080                 /* NB: File must be binary if we use fwrite */
1081                 fp = ftp2FILE(efDAT, NFILE, fnm, "wb");
1082                 for (i = 0; i < tel_mat; i++)
1083                 {
1084                     if (static_cast<int>(fwrite(rmsd_mat[i], sizeof(**rmsd_mat), tel_mat2, fp)) != tel_mat2)
1085                     {
1086                         gmx_fatal(FARGS, "Error writing to output file");
1087                     }
1088                 }
1089                 gmx_ffclose(fp);
1090             }
1091         }
1092         if (bBond)
1093         {
1094             fprintf(stderr, "\nMin. angle: %f, Max. angle: %f\n", bond_min, bond_max);
1095             if (bond_user_max != -1)
1096             {
1097                 bond_max = bond_user_max;
1098             }
1099             if (bond_user_min != -1)
1100             {
1101                 bond_min = bond_user_min;
1102             }
1103             if ((bond_user_max !=  -1) || (bond_user_min != -1))
1104             {
1105                 fprintf(stderr, "Bond angle Min and Max set to:\n"
1106                         "Min. angle: %f, Max. angle: %f\n", bond_min, bond_max);
1107             }
1108             rlo.r = 1; rlo.g = 1; rlo.b = 1;
1109             rhi.r = 0; rhi.g = 0; rhi.b = 0;
1110             sprintf(buf, "%s av. bond angle deviation", gn_rms[0]);
1111             write_xpm(opt2FILE("-bm", NFILE, fnm, "w"), 0, buf, "degrees",
1112                       output_env_get_time_label(oenv), output_env_get_time_label(oenv), tel_mat, tel_mat2,
1113                       axis, axis2, bond_mat, bond_min, bond_max, rlo, rhi, &nlevels);
1114         }
1115     }
1116
1117     bAv = opt2bSet("-a", NFILE, fnm);
1118
1119     /* Write the RMSD's to file */
1120     if (!bPrev)
1121     {
1122         sprintf(buf, "%s", whatxvgname[ewhat]);
1123     }
1124     else
1125     {
1126         sprintf(buf, "%s with frame %g %s ago", whatxvgname[ewhat],
1127                 time[prev*freq]-time[0], output_env_get_time_label(oenv));
1128     }
1129     fp = xvgropen(opt2fn("-o", NFILE, fnm), buf, output_env_get_xvgr_tlabel(oenv),
1130                   whatxvglabel[ewhat], oenv);
1131     if (output_env_get_print_xvgr_codes(oenv))
1132     {
1133         fprintf(fp, "@ subtitle \"%s%s after %s%s%s\"\n",
1134                 (nrms == 1) ? "" : "of ", gn_rms[0], fitgraphlabel[efit],
1135                 bFit     ? " to " : "", bFit ? gn_fit : "");
1136     }
1137     if (nrms != 1)
1138     {
1139         xvgr_legend(fp, nrms, (const char**)gn_rms, oenv);
1140     }
1141     for (i = 0; (i < teller); i++)
1142     {
1143         if (bSplit && i > 0 &&
1144             std::abs(time[bPrev ? freq*i : i]/output_env_get_time_factor(oenv)) < 1e-5)
1145         {
1146             fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
1147         }
1148         fprintf(fp, "%12.7f", time[bPrev ? freq*i : i]);
1149         for (j = 0; (j < nrms); j++)
1150         {
1151             fprintf(fp, " %12.7f", rls[j][i]);
1152             if (bAv)
1153             {
1154                 rlstot += rls[j][i];
1155             }
1156         }
1157         fprintf(fp, "\n");
1158     }
1159     xvgrclose(fp);
1160
1161     if (bMirror)
1162     {
1163         /* Write the mirror RMSD's to file */
1164         sprintf(buf, "%s with Mirror", whatxvgname[ewhat]);
1165         sprintf(buf2, "Mirror %s", whatxvglabel[ewhat]);
1166         fp = xvgropen(opt2fn("-mir", NFILE, fnm), buf, output_env_get_xvgr_tlabel(oenv),
1167                       buf2, oenv);
1168         if (nrms == 1)
1169         {
1170             if (output_env_get_print_xvgr_codes(oenv))
1171             {
1172                 fprintf(fp, "@ subtitle \"of %s after lsq fit to mirror of %s\"\n",
1173                         gn_rms[0], bFit ? gn_fit : "");
1174             }
1175         }
1176         else
1177         {
1178             if (output_env_get_print_xvgr_codes(oenv))
1179             {
1180                 fprintf(fp, "@ subtitle \"after lsq fit to mirror %s\"\n", bFit ? gn_fit : "");
1181             }
1182             xvgr_legend(fp, nrms, (const char**)gn_rms, oenv);
1183         }
1184         for (i = 0; (i < teller); i++)
1185         {
1186             if (bSplit && i > 0 && std::abs(time[i]) < 1e-5)
1187             {
1188                 fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
1189             }
1190             fprintf(fp, "%12.7f", time[i]);
1191             for (j = 0; (j < nrms); j++)
1192             {
1193                 fprintf(fp, " %12.7f", rlsm[j][i]);
1194             }
1195             fprintf(fp, "\n");
1196         }
1197         xvgrclose(fp);
1198     }
1199
1200     if (bAv)
1201     {
1202         sprintf(buf, "Average %s", whatxvgname[ewhat]);
1203         sprintf(buf2, "Average %s", whatxvglabel[ewhat]);
1204         fp = xvgropen(opt2fn("-a", NFILE, fnm), buf, "Residue", buf2, oenv);
1205         for (j = 0; (j < nrms); j++)
1206         {
1207             fprintf(fp, "%10d  %10g\n", j, rlstot/teller);
1208         }
1209         xvgrclose(fp);
1210     }
1211
1212     if (bNorm)
1213     {
1214         fp = xvgropen("aver.xvg", gn_rms[0], "Residue", whatxvglabel[ewhat], oenv);
1215         for (j = 0; (j < irms[0]); j++)
1216         {
1217             fprintf(fp, "%10d  %10g\n", j, rlsnorm[j]/teller);
1218         }
1219         xvgrclose(fp);
1220     }
1221     do_view(oenv, opt2fn_null("-a", NFILE, fnm), "-graphtype bar");
1222     do_view(oenv, opt2fn("-o", NFILE, fnm), nullptr);
1223     do_view(oenv, opt2fn_null("-mir", NFILE, fnm), nullptr);
1224     do_view(oenv, opt2fn_null("-m", NFILE, fnm), nullptr);
1225     do_view(oenv, opt2fn_null("-bm", NFILE, fnm), nullptr);
1226     do_view(oenv, opt2fn_null("-dist", NFILE, fnm), nullptr);
1227
1228     return 0;
1229 }