src/gromacs/gmxana/gmx_covar.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include <cmath>
  41 #include <cstring>
  42
  43 #include "gromacs/commandline/pargs.h"
  44 #include "gromacs/fileio/confio.h"
  45 #include "gromacs/fileio/matio.h"
  46 #include "gromacs/fileio/trxio.h"
  47 #include "gromacs/fileio/xvgr.h"
  48 #include "gromacs/gmxana/eigio.h"
  49 #include "gromacs/gmxana/gmx_ana.h"
  50 #include "gromacs/linearalgebra/eigensolver.h"
  51 #include "gromacs/math/do_fit.h"
  52 #include "gromacs/math/vec.h"
  53 #include "gromacs/pbcutil/pbc.h"
  54 #include "gromacs/pbcutil/rmpbc.h"
  55 #include "gromacs/topology/index.h"
  56 #include "gromacs/topology/topology.h"
  57 #include "gromacs/utility/arrayref.h"
  58 #include "gromacs/utility/arraysize.h"
  59 #include "gromacs/utility/cstringutil.h"
  60 #include "gromacs/utility/exceptions.h"
  61 #include "gromacs/utility/fatalerror.h"
  62 #include "gromacs/utility/futil.h"
  63 #include "gromacs/utility/smalloc.h"
  64 #include "gromacs/utility/sysinfo.h"
  65
  66
  67 namespace gmx
  68 {
  69
  70 namespace
  71 {
  72
  73 /*! \brief Throw an error if any element in index exceeds a given number.
  74  *
  75  * \param[in] indices to be acessed
  76  * \param[in] numAtoms to be accessed
  77  * \param[in] indexUsagePurpose to be more explicit in the error message
  78  *
  79  * \throws RangeError if any element in indices is larger than or equal numAtoms
  80  *
  81  */
  82 void throwErrorIfIndexOutOfBounds(ArrayRef<const int> indices, const int numAtoms, const std::string& indexUsagePurpose)
  83 {
  84     // do nothing if index is empty
  85     if (indices.empty())
  86     {
  87         return;
  88     }
  89     const int largestIndex = *std::max_element(indices.begin(), indices.end());
  90     if (largestIndex >= numAtoms)
  91     {
  92         GMX_THROW(RangeError("The provided structure file only contains " + std::to_string(numAtoms)
  93                              + " coordinates, but coordinate index " + std::to_string(largestIndex + 1)
  94                              + " was requested for " + indexUsagePurpose
  95                              + ". Make sure to update structure files "
  96                                "and index files if you store only a part of your system."));
  97     }
  98 };
  99
 100 } // namespace
 101
 102 } // namespace gmx
 103
 104 int gmx_covar(int argc, char* argv[])
 105 {
 106     const char* desc[] = {
 107         "[THISMODULE] calculates and diagonalizes the (mass-weighted)",
 108         "covariance matrix.",
 109         "All structures are fitted to the structure in the structure file.",
 110         "When this is not a run input file periodicity will not be taken into",
 111         "account. When the fit and analysis groups are identical and the analysis",
 112         "is non mass-weighted, the fit will also be non mass-weighted.",
 113         "[PAR]",
 114         "The eigenvectors are written to a trajectory file ([TT]-v[tt]).",
 115         "When the same atoms are used for the fit and the covariance analysis,",
 116         "the reference structure for the fit is written first with t=-1.",
 117         "The average (or reference when [TT]-ref[tt] is used) structure is",
 118         "written with t=0, the eigenvectors",
 119         "are written as frames with the eigenvector number and eigenvalue",
 120         "as step number and timestamp, respectively.",
 121         "[PAR]",
 122         "The eigenvectors can be analyzed with [gmx-anaeig].",
 123         "[PAR]",
 124         "Option [TT]-ascii[tt] writes the whole covariance matrix to",
 125         "an ASCII file. The order of the elements is: x1x1, x1y1, x1z1, x1x2, ...",
 126         "[PAR]",
 127         "Option [TT]-xpm[tt] writes the whole covariance matrix to an [REF].xpm[ref] file.",
 128         "[PAR]",
 129         "Option [TT]-xpma[tt] writes the atomic covariance matrix to an [REF].xpm[ref] file,",
 130         "i.e. for each atom pair the sum of the xx, yy and zz covariances is",
 131         "written.",
 132         "[PAR]",
 133         "Note that the diagonalization of a matrix requires memory and time",
 134         "that will increase at least as fast as than the square of the number",
 135         "of atoms involved. It is easy to run out of memory, in which",
 136         "case this tool will probably exit with a 'Segmentation fault'. You",
 137         "should consider carefully whether a reduced set of atoms will meet",
 138         "your needs for lower costs."
 139     };
 140     static gmx_bool bFit = TRUE, bRef = FALSE, bM = FALSE, bPBC = TRUE;
 141     static int      end  = -1;
 142     t_pargs         pa[] = {
 143         { "-fit", FALSE, etBOOL, { &bFit }, "Fit to a reference structure" },
 144         { "-ref",
 145           FALSE,
 146           etBOOL,
 147           { &bRef },
 148           "Use the deviation from the conformation in the structure file instead of from the "
 149           "average" },
 150         { "-mwa", FALSE, etBOOL, { &bM }, "Mass-weighted covariance analysis" },
 151         { "-last", FALSE, etINT, { &end }, "Last eigenvector to write away (-1 is till the last)" },
 152         { "-pbc", FALSE, etBOOL, { &bPBC }, "Apply corrections for periodic boundary conditions" }
 153     };
 154     FILE*             out = nullptr; /* initialization makes all compilers happy */
 155     t_trxstatus*      status;
 156     t_topology        top;
 157     PbcType           pbcType;
 158     t_atoms*          atoms;
 159     rvec *            x, *xread, *xref, *xav, *xproj;
 160     matrix            box, zerobox;
 161     real *            sqrtm, *mat, *eigenvalues, sum, trace, inv_nframes;
 162     real              t, tstart, tend, **mat2;
 163     real              xj, *w_rls = nullptr;
 164     real              min, max, *axis;
 165     int               natoms, nat, nframes0, nframes, nlevels;
 166     int64_t           ndim, i, j, k, l;
 167     int               WriteXref;
 168     const char *      fitfile, *trxfile, *ndxfile;
 169     const char *      eigvalfile, *eigvecfile, *averfile, *logfile;
 170     const char *      asciifile, *xpmfile, *xpmafile;
 171     char              str[STRLEN], *fitname, *ananame;
 172     int               d, dj, nfit;
 173     int *             index, *ifit;
 174     gmx_bool          bDiffMass1, bDiffMass2;
 175     t_rgb             rlo, rmi, rhi;
 176     real*             eigenvectors;
 177     gmx_output_env_t* oenv;
 178     gmx_rmpbc_t       gpbc = nullptr;
 179
 180     t_filenm fnm[] = {
 181         { efTRX, "-f", nullptr, ffREAD },     { efTPS, nullptr, nullptr, ffREAD },
 182         { efNDX, nullptr, nullptr, ffOPTRD }, { efXVG, nullptr, "eigenval", ffWRITE },
 183         { efTRN, "-v", "eigenvec", ffWRITE }, { efSTO, "-av", "average.pdb", ffWRITE },
 184         { efLOG, nullptr, "covar", ffWRITE }, { efDAT, "-ascii", "covar", ffOPTWR },
 185         { efXPM, "-xpm", "covar", ffOPTWR },  { efXPM, "-xpma", "covara", ffOPTWR }
 186     };
 187 #define NFILE asize(fnm)
 188
 189     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT, NFILE, fnm, asize(pa), pa,
 190                            asize(desc), desc, 0, nullptr, &oenv))
 191     {
 192         return 0;
 193     }
 194
 195     clear_mat(zerobox);
 196
 197     fitfile    = ftp2fn(efTPS, NFILE, fnm);
 198     trxfile    = ftp2fn(efTRX, NFILE, fnm);
 199     ndxfile    = ftp2fn_null(efNDX, NFILE, fnm);
 200     eigvalfile = ftp2fn(efXVG, NFILE, fnm);
 201     eigvecfile = ftp2fn(efTRN, NFILE, fnm);
 202     averfile   = ftp2fn(efSTO, NFILE, fnm);
 203     logfile    = ftp2fn(efLOG, NFILE, fnm);
 204     asciifile  = opt2fn_null("-ascii", NFILE, fnm);
 205     xpmfile    = opt2fn_null("-xpm", NFILE, fnm);
 206     xpmafile   = opt2fn_null("-xpma", NFILE, fnm);
 207
 208     read_tps_conf(fitfile, &top, &pbcType, &xref, nullptr, box, TRUE);
 209     atoms = &top.atoms;
 210
 211     if (bFit)
 212     {
 213         printf("\nChoose a group for the least squares fit\n");
 214         get_index(atoms, ndxfile, 1, &nfit, &ifit, &fitname);
 215         // Make sure that we never attempt to access a coordinate out of range
 216         gmx::throwErrorIfIndexOutOfBounds({ ifit, ifit + nfit }, atoms->nr, "fitting");
 217         if (nfit < 3)
 218         {
 219             gmx_fatal(FARGS, "Need >= 3 points to fit!\n");
 220         }
 221     }
 222     else
 223     {
 224         nfit = 0;
 225     }
 226     printf("\nChoose a group for the covariance analysis\n");
 227     get_index(atoms, ndxfile, 1, &natoms, &index, &ananame);
 228     gmx::throwErrorIfIndexOutOfBounds({ index, index + natoms }, atoms->nr, "analysis");
 229
 230     bDiffMass1 = FALSE;
 231     if (bFit)
 232     {
 233         snew(w_rls, atoms->nr);
 234         for (i = 0; (i < nfit); i++)
 235         {
 236             w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
 237             if (i)
 238             {
 239                 bDiffMass1 = bDiffMass1 || (w_rls[ifit[i]] != w_rls[ifit[i - 1]]);
 240             }
 241         }
 242     }
 243     bDiffMass2 = FALSE;
 244     snew(sqrtm, natoms);
 245     for (i = 0; (i < natoms); i++)
 246     {
 247         if (bM)
 248         {
 249             sqrtm[i] = std::sqrt(atoms->atom[index[i]].m);
 250             if (i)
 251             {
 252                 bDiffMass2 = bDiffMass2 || (sqrtm[i] != sqrtm[i - 1]);
 253             }
 254         }
 255         else
 256         {
 257             sqrtm[i] = 1.0;
 258         }
 259     }
 260
 261     if (bFit && bDiffMass1 && !bDiffMass2)
 262     {
 263         bDiffMass1 = natoms != nfit;
 264         for (i = 0; (i < natoms) && !bDiffMass1; i++)
 265         {
 266             bDiffMass1 = index[i] != ifit[i];
 267         }
 268         if (!bDiffMass1)
 269         {
 270             fprintf(stderr,
 271                     "\n"
 272                     "Note: the fit and analysis group are identical,\n"
 273                     "      while the fit is mass weighted and the analysis is not.\n"
 274                     "      Making the fit non mass weighted.\n\n");
 275             for (i = 0; (i < nfit); i++)
 276             {
 277                 w_rls[ifit[i]] = 1.0;
 278             }
 279         }
 280     }
 281
 282     /* Prepare reference frame */
 283     if (bPBC)
 284     {
 285         gpbc = gmx_rmpbc_init(&top.idef, pbcType, atoms->nr);
 286         gmx_rmpbc(gpbc, atoms->nr, box, xref);
 287     }
 288     if (bFit)
 289     {
 290         reset_x(nfit, ifit, atoms->nr, nullptr, xref, w_rls);
 291     }
 292
 293     snew(x, natoms);
 294     snew(xav, natoms);
 295     ndim = natoms * DIM;
 296     if (std::sqrt(static_cast<real>(INT64_MAX)) < static_cast<real>(ndim))
 297     {
 298         gmx_fatal(FARGS, "Number of degrees of freedoms to large for matrix.\n");
 299     }
 300     snew(mat, ndim * ndim);
 301
 302     fprintf(stderr, "Calculating the average structure ...\n");
 303     nframes0 = 0;
 304     nat      = read_first_x(oenv, &status, trxfile, &t, &xread, box);
 305     if (nat != atoms->nr)
 306     {
 307         fprintf(stderr,
 308                 "\nWARNING: number of atoms in structure file (%d) and trajectory (%d) do not "
 309                 "match\n",
 310                 natoms, nat);
 311     }
 312     gmx::throwErrorIfIndexOutOfBounds({ ifit, ifit + nfit }, nat, "fitting");
 313     gmx::throwErrorIfIndexOutOfBounds({ index, index + natoms }, nat, "analysis");
 314
 315     do
 316     {
 317         nframes0++;
 318         /* calculate x: a fitted struture of the selected atoms */
 319         if (bPBC)
 320         {
 321             gmx_rmpbc(gpbc, nat, box, xread);
 322         }
 323         if (bFit)
 324         {
 325             reset_x(nfit, ifit, nat, nullptr, xread, w_rls);
 326             do_fit(nat, w_rls, xref, xread);
 327         }
 328         for (i = 0; i < natoms; i++)
 329         {
 330             rvec_inc(xav[i], xread[index[i]]);
 331         }
 332     } while (read_next_x(oenv, status, &t, xread, box));
 333     close_trx(status);
 334
 335     inv_nframes = 1.0 / nframes0;
 336     for (i = 0; i < natoms; i++)
 337     {
 338         for (d = 0; d < DIM; d++)
 339         {
 340             xav[i][d] *= inv_nframes;
 341             xread[index[i]][d] = xav[i][d];
 342         }
 343     }
 344     write_sto_conf_indexed(opt2fn("-av", NFILE, fnm), "Average structure", atoms, xread, nullptr,
 345                            PbcType::No, zerobox, natoms, index);
 346     sfree(xread);
 347
 348     fprintf(stderr, "Constructing covariance matrix (%dx%d) ...\n", static_cast<int>(ndim),
 349             static_cast<int>(ndim));
 350     nframes = 0;
 351     nat     = read_first_x(oenv, &status, trxfile, &t, &xread, box);
 352     tstart  = t;
 353     do
 354     {
 355         nframes++;
 356         tend = t;
 357         /* calculate x: a (fitted) structure of the selected atoms */
 358         if (bPBC)
 359         {
 360             gmx_rmpbc(gpbc, nat, box, xread);
 361         }
 362         if (bFit)
 363         {
 364             reset_x(nfit, ifit, nat, nullptr, xread, w_rls);
 365             do_fit(nat, w_rls, xref, xread);
 366         }
 367         if (bRef)
 368         {
 369             for (i = 0; i < natoms; i++)
 370             {
 371                 rvec_sub(xread[index[i]], xref[index[i]], x[i]);
 372             }
 373         }
 374         else
 375         {
 376             for (i = 0; i < natoms; i++)
 377             {
 378                 rvec_sub(xread[index[i]], xav[i], x[i]);
 379             }
 380         }
 381
 382         for (j = 0; j < natoms; j++)
 383         {
 384             for (dj = 0; dj < DIM; dj++)
 385             {
 386                 k  = ndim * (DIM * j + dj);
 387                 xj = x[j][dj];
 388                 for (i = j; i < natoms; i++)
 389                 {
 390                     l = k + DIM * i;
 391                     for (d = 0; d < DIM; d++)
 392                     {
 393                         mat[l + d] += x[i][d] * xj;
 394                     }
 395                 }
 396             }
 397         }
 398     } while (read_next_x(oenv, status, &t, xread, box) && (bRef || nframes < nframes0));
 399     close_trx(status);
 400     gmx_rmpbc_done(gpbc);
 401
 402     fprintf(stderr, "Read %d frames\n", nframes);
 403
 404     if (bRef)
 405     {
 406         /* copy the reference structure to the ouput array x */
 407         snew(xproj, natoms);
 408         for (i = 0; i < natoms; i++)
 409         {
 410             copy_rvec(xref[index[i]], xproj[i]);
 411         }
 412     }
 413     else
 414     {
 415         xproj = xav;
 416     }
 417
 418     /* correct the covariance matrix for the mass */
 419     inv_nframes = 1.0 / nframes;
 420     for (j = 0; j < natoms; j++)
 421     {
 422         for (dj = 0; dj < DIM; dj++)
 423         {
 424             for (i = j; i < natoms; i++)
 425             {
 426                 k = ndim * (DIM * j + dj) + DIM * i;
 427                 for (d = 0; d < DIM; d++)
 428                 {
 429                     mat[k + d] = mat[k + d] * inv_nframes * sqrtm[i] * sqrtm[j];
 430                 }
 431             }
 432         }
 433     }
 434
 435     /* symmetrize the matrix */
 436     for (j = 0; j < ndim; j++)
 437     {
 438         for (i = j; i < ndim; i++)
 439         {
 440             mat[ndim * i + j] = mat[ndim * j + i];
 441         }
 442     }
 443
 444     trace = 0;
 445     for (i = 0; i < ndim; i++)
 446     {
 447         trace += mat[i * ndim + i];
 448     }
 449     fprintf(stderr, "\nTrace of the covariance matrix: %g (%snm^2)\n", trace, bM ? "u " : "");
 450
 451     if (asciifile)
 452     {
 453         out = gmx_ffopen(asciifile, "w");
 454         for (j = 0; j < ndim; j++)
 455         {
 456             for (i = 0; i < ndim; i += 3)
 457             {
 458                 fprintf(out, "%g %g %g\n", mat[ndim * j + i], mat[ndim * j + i + 1],
 459                         mat[ndim * j + i + 2]);
 460             }
 461         }
 462         gmx_ffclose(out);
 463     }
 464
 465     if (xpmfile)
 466     {
 467         min = 0;
 468         max = 0;
 469         snew(mat2, ndim);
 470         for (j = 0; j < ndim; j++)
 471         {
 472             mat2[j] = &(mat[ndim * j]);
 473             for (i = 0; i <= j; i++)
 474             {
 475                 if (mat2[j][i] < min)
 476                 {
 477                     min = mat2[j][i];
 478                 }
 479                 if (mat2[j][j] > max)
 480                 {
 481                     max = mat2[j][i];
 482                 }
 483             }
 484         }
 485         snew(axis, ndim);
 486         for (i = 0; i < ndim; i++)
 487         {
 488             axis[i] = i + 1;
 489         }
 490         rlo.r   = 0;
 491         rlo.g   = 0;
 492         rlo.b   = 1;
 493         rmi.r   = 1;
 494         rmi.g   = 1;
 495         rmi.b   = 1;
 496         rhi.r   = 1;
 497         rhi.g   = 0;
 498         rhi.b   = 0;
 499         out     = gmx_ffopen(xpmfile, "w");
 500         nlevels = 80;
 501         write_xpm3(out, 0, "Covariance", bM ? "u nm^2" : "nm^2", "dim", "dim", ndim, ndim, axis,
 502                    axis, mat2, min, 0.0, max, rlo, rmi, rhi, &nlevels);
 503         gmx_ffclose(out);
 504         sfree(axis);
 505         sfree(mat2);
 506     }
 507
 508     if (xpmafile)
 509     {
 510         min = 0;
 511         max = 0;
 512         snew(mat2, ndim / DIM);
 513         for (i = 0; i < ndim / DIM; i++)
 514         {
 515             snew(mat2[i], ndim / DIM);
 516         }
 517         for (j = 0; j < ndim / DIM; j++)
 518         {
 519             for (i = 0; i <= j; i++)
 520             {
 521                 mat2[j][i] = 0;
 522                 for (d = 0; d < DIM; d++)
 523                 {
 524                     mat2[j][i] += mat[ndim * (DIM * j + d) + DIM * i + d];
 525                 }
 526                 if (mat2[j][i] < min)
 527                 {
 528                     min = mat2[j][i];
 529                 }
 530                 if (mat2[j][j] > max)
 531                 {
 532                     max = mat2[j][i];
 533                 }
 534                 mat2[i][j] = mat2[j][i];
 535             }
 536         }
 537         snew(axis, ndim / DIM);
 538         for (i = 0; i < ndim / DIM; i++)
 539         {
 540             axis[i] = i + 1;
 541         }
 542         rlo.r   = 0;
 543         rlo.g   = 0;
 544         rlo.b   = 1;
 545         rmi.r   = 1;
 546         rmi.g   = 1;
 547         rmi.b   = 1;
 548         rhi.r   = 1;
 549         rhi.g   = 0;
 550         rhi.b   = 0;
 551         out     = gmx_ffopen(xpmafile, "w");
 552         nlevels = 80;
 553         write_xpm3(out, 0, "Covariance", bM ? "u nm^2" : "nm^2", "atom", "atom", ndim / DIM,
 554                    ndim / DIM, axis, axis, mat2, min, 0.0, max, rlo, rmi, rhi, &nlevels);
 555         gmx_ffclose(out);
 556         sfree(axis);
 557         for (i = 0; i < ndim / DIM; i++)
 558         {
 559             sfree(mat2[i]);
 560         }
 561         sfree(mat2);
 562     }
 563
 564
 565     /* call diagonalization routine */
 566
 567     snew(eigenvalues, ndim);
 568     snew(eigenvectors, ndim * ndim);
 569
 570     std::memcpy(eigenvectors, mat, ndim * ndim * sizeof(real));
 571     fprintf(stderr, "\nDiagonalizing ...\n");
 572     fflush(stderr);
 573     eigensolver(eigenvectors, ndim, 0, ndim, eigenvalues, mat);
 574     sfree(eigenvectors);
 575
 576     /* now write the output */
 577
 578     sum = 0;
 579     for (i = 0; i < ndim; i++)
 580     {
 581         sum += eigenvalues[i];
 582     }
 583     fprintf(stderr, "\nSum of the eigenvalues: %g (%snm^2)\n", sum, bM ? "u " : "");
 584     if (std::abs(trace - sum) > 0.01 * trace)
 585     {
 586         fprintf(stderr,
 587                 "\nWARNING: eigenvalue sum deviates from the trace of the covariance matrix\n");
 588     }
 589
 590     /* Set 'end', the maximum eigenvector and -value index used for output */
 591     if (end == -1)
 592     {
 593         if (nframes - 1 < ndim)
 594         {
 595             end = nframes - 1;
 596             fprintf(stderr,
 597                     "\nWARNING: there are fewer frames in your trajectory than there are\n");
 598             fprintf(stderr, "degrees of freedom in your system. Only generating the first\n");
 599             fprintf(stderr, "%d out of %d eigenvectors and eigenvalues.\n", end, static_cast<int>(ndim));
 600         }
 601         else
 602         {
 603             end = ndim;
 604         }
 605     }
 606
 607     fprintf(stderr, "\nWriting eigenvalues to %s\n", eigvalfile);
 608
 609     sprintf(str, "(%snm\\S2\\N)", bM ? "u " : "");
 610     out = xvgropen(eigvalfile, "Eigenvalues of the covariance matrix", "Eigenvector index", str, oenv);
 611     for (i = 0; (i < end); i++)
 612     {
 613         fprintf(out, "%10d %g\n", static_cast<int>(i + 1), eigenvalues[ndim - 1 - i]);
 614     }
 615     xvgrclose(out);
 616
 617     if (bFit)
 618     {
 619         /* misuse lambda: 0/1 mass weighted analysis no/yes */
 620         if (nfit == natoms)
 621         {
 622             WriteXref = eWXR_YES;
 623             for (i = 0; i < nfit; i++)
 624             {
 625                 copy_rvec(xref[ifit[i]], x[i]);
 626             }
 627         }
 628         else
 629         {
 630             WriteXref = eWXR_NO;
 631         }
 632     }
 633     else
 634     {
 635         /* misuse lambda: -1 for no fit */
 636         WriteXref = eWXR_NOFIT;
 637     }
 638
 639     write_eigenvectors(eigvecfile, natoms, mat, TRUE, 1, end, WriteXref, x, bDiffMass1, xproj, bM,
 640                        eigenvalues);
 641
 642     out = gmx_ffopen(logfile, "w");
 643
 644     fprintf(out, "Covariance analysis log, written %s\n", gmx_format_current_time().c_str());
 645
 646     fprintf(out, "Program: %s\n", argv[0]);
 647     gmx_getcwd(str, STRLEN);
 648
 649     fprintf(out, "Working directory: %s\n\n", str);
 650
 651     fprintf(out, "Read %d frames from %s (time %g to %g %s)\n", nframes, trxfile,
 652             output_env_conv_time(oenv, tstart), output_env_conv_time(oenv, tend),
 653             output_env_get_time_unit(oenv).c_str());
 654     if (bFit)
 655     {
 656         fprintf(out, "Read reference structure for fit from %s\n", fitfile);
 657     }
 658     if (ndxfile)
 659     {
 660         fprintf(out, "Read index groups from %s\n", ndxfile);
 661     }
 662     fprintf(out, "\n");
 663
 664     fprintf(out, "Analysis group is '%s' (%d atoms)\n", ananame, natoms);
 665     if (bFit)
 666     {
 667         fprintf(out, "Fit group is '%s' (%d atoms)\n", fitname, nfit);
 668     }
 669     else
 670     {
 671         fprintf(out, "No fit was used\n");
 672     }
 673     fprintf(out, "Analysis is %smass weighted\n", bDiffMass2 ? "" : "non-");
 674     if (bFit)
 675     {
 676         fprintf(out, "Fit is %smass weighted\n", bDiffMass1 ? "" : "non-");
 677     }
 678     fprintf(out, "Diagonalized the %dx%d covariance matrix\n", static_cast<int>(ndim),
 679             static_cast<int>(ndim));
 680     fprintf(out, "Trace of the covariance matrix before diagonalizing: %g\n", trace);
 681     fprintf(out, "Trace of the covariance matrix after diagonalizing: %g\n\n", sum);
 682
 683     fprintf(out, "Wrote %d eigenvalues to %s\n", static_cast<int>(end), eigvalfile);
 684     if (WriteXref == eWXR_YES)
 685     {
 686         fprintf(out, "Wrote reference structure to %s\n", eigvecfile);
 687     }
 688     fprintf(out, "Wrote average structure to %s and %s\n", averfile, eigvecfile);
 689     fprintf(out, "Wrote eigenvectors %d to %d to %s\n", 1, end, eigvecfile);
 690
 691     gmx_ffclose(out);
 692
 693     fprintf(stderr, "Wrote the log to %s\n", logfile);
 694
 695     return 0;
 696 }