src/gromacs/gmxana/gmx_rms.c

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