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