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