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, 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] largestOkayIndex to be accessed
  77  * \param[in] indexUsagePurpose to be more explicit in the error message
  78  *
  79  * \throws RangeError if largestOkayIndex is larger than any element in indices
  80  *
  81  */
  82 void throwErrorIfIndexOutOfBounds(ArrayRef<const int> indices,
  83                                   const int           largestOkayIndex,
  84                                   const std::string&  indexUsagePurpose)
  85 {
  86     // do nothing if index is empty
  87     if (indices.ssize() == 0)
  88     {
  89         return;
  90     }
  91     const int largestIndex = *std::max_element(indices.begin(), indices.end());
  92     if (largestIndex < largestOkayIndex)
  93     {
  94         GMX_THROW(RangeError("The provided structure file only contains "
  95                              + std::to_string(largestOkayIndex) + " coordinates, but coordinate index "
  96                              + std::to_string(largestIndex) + " was requested for " + indexUsagePurpose
  97                              + ". Make sure to update structure files "
  98                                "and index files if you store only a part of your system."));
  99     }
 100 };
 101
 102 } // namespace
 103
 104 } // namespace gmx
 105
 106 int gmx_covar(int argc, char* argv[])
 107 {
 108     const char* desc[] = {
 109         "[THISMODULE] calculates and diagonalizes the (mass-weighted)",
 110         "covariance matrix.",
 111         "All structures are fitted to the structure in the structure file.",
 112         "When this is not a run input file periodicity will not be taken into",
 113         "account. When the fit and analysis groups are identical and the analysis",
 114         "is non mass-weighted, the fit will also be non mass-weighted.",
 115         "[PAR]",
 116         "The eigenvectors are written to a trajectory file ([TT]-v[tt]).",
 117         "When the same atoms are used for the fit and the covariance analysis,",
 118         "the reference structure for the fit is written first with t=-1.",
 119         "The average (or reference when [TT]-ref[tt] is used) structure is",
 120         "written with t=0, the eigenvectors",
 121         "are written as frames with the eigenvector number and eigenvalue",
 122         "as step number and timestamp, respectively.",
 123         "[PAR]",
 124         "The eigenvectors can be analyzed with [gmx-anaeig].",
 125         "[PAR]",
 126         "Option [TT]-ascii[tt] writes the whole covariance matrix to",
 127         "an ASCII file. The order of the elements is: x1x1, x1y1, x1z1, x1x2, ...",
 128         "[PAR]",
 129         "Option [TT]-xpm[tt] writes the whole covariance matrix to an [REF].xpm[ref] file.",
 130         "[PAR]",
 131         "Option [TT]-xpma[tt] writes the atomic covariance matrix to an [REF].xpm[ref] file,",
 132         "i.e. for each atom pair the sum of the xx, yy and zz covariances is",
 133         "written.",
 134         "[PAR]",
 135         "Note that the diagonalization of a matrix requires memory and time",
 136         "that will increase at least as fast as than the square of the number",
 137         "of atoms involved. It is easy to run out of memory, in which",
 138         "case this tool will probably exit with a 'Segmentation fault'. You",
 139         "should consider carefully whether a reduced set of atoms will meet",
 140         "your needs for lower costs."
 141     };
 142     static gmx_bool bFit = TRUE, bRef = FALSE, bM = FALSE, bPBC = TRUE;
 143     static int      end  = -1;
 144     t_pargs         pa[] = {
 145         { "-fit", FALSE, etBOOL, { &bFit }, "Fit to a reference structure" },
 146         { "-ref",
 147           FALSE,
 148           etBOOL,
 149           { &bRef },
 150           "Use the deviation from the conformation in the structure file instead of from the "
 151           "average" },
 152         { "-mwa", FALSE, etBOOL, { &bM }, "Mass-weighted covariance analysis" },
 153         { "-last", FALSE, etINT, { &end }, "Last eigenvector to write away (-1 is till the last)" },
 154         { "-pbc", FALSE, etBOOL, { &bPBC }, "Apply corrections for periodic boundary conditions" }
 155     };
 156     FILE*             out = nullptr; /* initialization makes all compilers happy */
 157     t_trxstatus*      status;
 158     t_topology        top;
 159     PbcType           pbcType;
 160     t_atoms*          atoms;
 161     rvec *            x, *xread, *xref, *xav, *xproj;
 162     matrix            box, zerobox;
 163     real *            sqrtm, *mat, *eigenvalues, sum, trace, inv_nframes;
 164     real              t, tstart, tend, **mat2;
 165     real              xj, *w_rls = nullptr;
 166     real              min, max, *axis;
 167     int               natoms, nat, nframes0, nframes, nlevels;
 168     int64_t           ndim, i, j, k, l;
 169     int               WriteXref;
 170     const char *      fitfile, *trxfile, *ndxfile;
 171     const char *      eigvalfile, *eigvecfile, *averfile, *logfile;
 172     const char *      asciifile, *xpmfile, *xpmafile;
 173     char              str[STRLEN], *fitname, *ananame;
 174     int               d, dj, nfit;
 175     int *             index, *ifit;
 176     gmx_bool          bDiffMass1, bDiffMass2;
 177     t_rgb             rlo, rmi, rhi;
 178     real*             eigenvectors;
 179     gmx_output_env_t* oenv;
 180     gmx_rmpbc_t       gpbc = nullptr;
 181
 182     t_filenm fnm[] = {
 183         { efTRX, "-f", nullptr, ffREAD },     { efTPS, nullptr, nullptr, ffREAD },
 184         { efNDX, nullptr, nullptr, ffOPTRD }, { efXVG, nullptr, "eigenval", ffWRITE },
 185         { efTRN, "-v", "eigenvec", ffWRITE }, { efSTO, "-av", "average.pdb", ffWRITE },
 186         { efLOG, nullptr, "covar", ffWRITE }, { efDAT, "-ascii", "covar", ffOPTWR },
 187         { efXPM, "-xpm", "covar", ffOPTWR },  { efXPM, "-xpma", "covara", ffOPTWR }
 188     };
 189 #define NFILE asize(fnm)
 190
 191     if (!parse_common_args(
 192                 &argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
 193     {
 194         return 0;
 195     }
 196
 197     clear_mat(zerobox);
 198
 199     fitfile    = ftp2fn(efTPS, NFILE, fnm);
 200     trxfile    = ftp2fn(efTRX, NFILE, fnm);
 201     ndxfile    = ftp2fn_null(efNDX, NFILE, fnm);
 202     eigvalfile = ftp2fn(efXVG, NFILE, fnm);
 203     eigvecfile = ftp2fn(efTRN, NFILE, fnm);
 204     averfile   = ftp2fn(efSTO, NFILE, fnm);
 205     logfile    = ftp2fn(efLOG, NFILE, fnm);
 206     asciifile  = opt2fn_null("-ascii", NFILE, fnm);
 207     xpmfile    = opt2fn_null("-xpm", NFILE, fnm);
 208     xpmafile   = opt2fn_null("-xpma", NFILE, fnm);
 209
 210     read_tps_conf(fitfile, &top, &pbcType, &xref, nullptr, box, TRUE);
 211     atoms = &top.atoms;
 212
 213     if (bFit)
 214     {
 215         printf("\nChoose a group for the least squares fit\n");
 216         get_index(atoms, ndxfile, 1, &nfit, &ifit, &fitname);
 217         // Make sure that we never attempt to access a coordinate out of range
 218         gmx::throwErrorIfIndexOutOfBounds({ ifit, ifit + nfit }, atoms->nr, "fitting");
 219         if (nfit < 3)
 220         {
 221             gmx_fatal(FARGS, "Need >= 3 points to fit!\n");
 222         }
 223     }
 224     else
 225     {
 226         nfit = 0;
 227     }
 228     printf("\nChoose a group for the covariance analysis\n");
 229     get_index(atoms, ndxfile, 1, &natoms, &index, &ananame);
 230     gmx::throwErrorIfIndexOutOfBounds({ index, index + natoms }, atoms->nr, "analysis");
 231
 232     bDiffMass1 = FALSE;
 233     if (bFit)
 234     {
 235         snew(w_rls, atoms->nr);
 236         for (i = 0; (i < nfit); i++)
 237         {
 238             w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
 239             if (i)
 240             {
 241                 bDiffMass1 = bDiffMass1 || (w_rls[ifit[i]] != w_rls[ifit[i - 1]]);
 242             }
 243         }
 244     }
 245     bDiffMass2 = FALSE;
 246     snew(sqrtm, natoms);
 247     for (i = 0; (i < natoms); i++)
 248     {
 249         if (bM)
 250         {
 251             sqrtm[i] = std::sqrt(atoms->atom[index[i]].m);
 252             if (i)
 253             {
 254                 bDiffMass2 = bDiffMass2 || (sqrtm[i] != sqrtm[i - 1]);
 255             }
 256         }
 257         else
 258         {
 259             sqrtm[i] = 1.0;
 260         }
 261     }
 262
 263     if (bFit && bDiffMass1 && !bDiffMass2)
 264     {
 265         bDiffMass1 = natoms != nfit;
 266         for (i = 0; (i < natoms) && !bDiffMass1; i++)
 267         {
 268             bDiffMass1 = index[i] != ifit[i];
 269         }
 270         if (!bDiffMass1)
 271         {
 272             fprintf(stderr,
 273                     "\n"
 274                     "Note: the fit and analysis group are identical,\n"
 275                     "      while the fit is mass weighted and the analysis is not.\n"
 276                     "      Making the fit non mass weighted.\n\n");
 277             for (i = 0; (i < nfit); i++)
 278             {
 279                 w_rls[ifit[i]] = 1.0;
 280             }
 281         }
 282     }
 283
 284     /* Prepare reference frame */
 285     if (bPBC)
 286     {
 287         gpbc = gmx_rmpbc_init(&top.idef, pbcType, atoms->nr);
 288         gmx_rmpbc(gpbc, atoms->nr, box, xref);
 289     }
 290     if (bFit)
 291     {
 292         reset_x(nfit, ifit, atoms->nr, nullptr, xref, w_rls);
 293     }
 294
 295     snew(x, natoms);
 296     snew(xav, natoms);
 297     ndim = natoms * DIM;
 298     if (std::sqrt(static_cast<real>(INT64_MAX)) < static_cast<real>(ndim))
 299     {
 300         gmx_fatal(FARGS, "Number of degrees of freedoms to large for matrix.\n");
 301     }
 302     snew(mat, ndim * ndim);
 303
 304     fprintf(stderr, "Calculating the average structure ...\n");
 305     nframes0 = 0;
 306     nat      = read_first_x(oenv, &status, trxfile, &t, &xread, box);
 307     if (nat != atoms->nr)
 308     {
 309         fprintf(stderr,
 310                 "\nWARNING: number of atoms in structure file (%d) and trajectory (%d) do not "
 311                 "match\n",
 312                 natoms,
 313                 nat);
 314     }
 315     gmx::throwErrorIfIndexOutOfBounds({ ifit, ifit + nfit }, nat, "fitting");
 316     gmx::throwErrorIfIndexOutOfBounds({ index, index + natoms }, nat, "analysis");
 317
 318     do
 319     {
 320         nframes0++;
 321         /* calculate x: a fitted struture of the selected atoms */
 322         if (bPBC)
 323         {
 324             gmx_rmpbc(gpbc, nat, box, xread);
 325         }
 326         if (bFit)
 327         {
 328             reset_x(nfit, ifit, nat, nullptr, xread, w_rls);
 329             do_fit(nat, w_rls, xref, xread);
 330         }
 331         for (i = 0; i < natoms; i++)
 332         {
 333             rvec_inc(xav[i], xread[index[i]]);
 334         }
 335     } while (read_next_x(oenv, status, &t, xread, box));
 336     close_trx(status);
 337
 338     inv_nframes = 1.0 / nframes0;
 339     for (i = 0; i < natoms; i++)
 340     {
 341         for (d = 0; d < DIM; d++)
 342         {
 343             xav[i][d] *= inv_nframes;
 344             xread[index[i]][d] = xav[i][d];
 345         }
 346     }
 347     write_sto_conf_indexed(
 348             opt2fn("-av", NFILE, fnm), "Average structure", atoms, xread, nullptr, PbcType::No, zerobox, natoms, index);
 349     sfree(xread);
 350
 351     fprintf(stderr,
 352             "Constructing covariance matrix (%dx%d) ...\n",
 353             static_cast<int>(ndim),
 354             static_cast<int>(ndim));
 355     nframes = 0;
 356     nat     = read_first_x(oenv, &status, trxfile, &t, &xread, box);
 357     tstart  = t;
 358     do
 359     {
 360         nframes++;
 361         tend = t;
 362         /* calculate x: a (fitted) structure of the selected atoms */
 363         if (bPBC)
 364         {
 365             gmx_rmpbc(gpbc, nat, box, xread);
 366         }
 367         if (bFit)
 368         {
 369             reset_x(nfit, ifit, nat, nullptr, xread, w_rls);
 370             do_fit(nat, w_rls, xref, xread);
 371         }
 372         if (bRef)
 373         {
 374             for (i = 0; i < natoms; i++)
 375             {
 376                 rvec_sub(xread[index[i]], xref[index[i]], x[i]);
 377             }
 378         }
 379         else
 380         {
 381             for (i = 0; i < natoms; i++)
 382             {
 383                 rvec_sub(xread[index[i]], xav[i], x[i]);
 384             }
 385         }
 386
 387         for (j = 0; j < natoms; j++)
 388         {
 389             for (dj = 0; dj < DIM; dj++)
 390             {
 391                 k  = ndim * (DIM * j + dj);
 392                 xj = x[j][dj];
 393                 for (i = j; i < natoms; i++)
 394                 {
 395                     l = k + DIM * i;
 396                     for (d = 0; d < DIM; d++)
 397                     {
 398                         mat[l + d] += x[i][d] * xj;
 399                     }
 400                 }
 401             }
 402         }
 403     } while (read_next_x(oenv, status, &t, xread, box) && (bRef || nframes < nframes0));
 404     close_trx(status);
 405     gmx_rmpbc_done(gpbc);
 406
 407     fprintf(stderr, "Read %d frames\n", nframes);
 408
 409     if (bRef)
 410     {
 411         /* copy the reference structure to the ouput array x */
 412         snew(xproj, natoms);
 413         for (i = 0; i < natoms; i++)
 414         {
 415             copy_rvec(xref[index[i]], xproj[i]);
 416         }
 417     }
 418     else
 419     {
 420         xproj = xav;
 421     }
 422
 423     /* correct the covariance matrix for the mass */
 424     inv_nframes = 1.0 / nframes;
 425     for (j = 0; j < natoms; j++)
 426     {
 427         for (dj = 0; dj < DIM; dj++)
 428         {
 429             for (i = j; i < natoms; i++)
 430             {
 431                 k = ndim * (DIM * j + dj) + DIM * i;
 432                 for (d = 0; d < DIM; d++)
 433                 {
 434                     mat[k + d] = mat[k + d] * inv_nframes * sqrtm[i] * sqrtm[j];
 435                 }
 436             }
 437         }
 438     }
 439
 440     /* symmetrize the matrix */
 441     for (j = 0; j < ndim; j++)
 442     {
 443         for (i = j; i < ndim; i++)
 444         {
 445             mat[ndim * i + j] = mat[ndim * j + i];
 446         }
 447     }
 448
 449     trace = 0;
 450     for (i = 0; i < ndim; i++)
 451     {
 452         trace += mat[i * ndim + i];
 453     }
 454     fprintf(stderr, "\nTrace of the covariance matrix: %g (%snm^2)\n", trace, bM ? "u " : "");
 455
 456     if (asciifile)
 457     {
 458         out = gmx_ffopen(asciifile, "w");
 459         for (j = 0; j < ndim; j++)
 460         {
 461             for (i = 0; i < ndim; i += 3)
 462             {
 463                 fprintf(out, "%g %g %g\n", mat[ndim * j + i], mat[ndim * j + i + 1], mat[ndim * j + i + 2]);
 464             }
 465         }
 466         gmx_ffclose(out);
 467     }
 468
 469     if (xpmfile)
 470     {
 471         min = 0;
 472         max = 0;
 473         snew(mat2, ndim);
 474         for (j = 0; j < ndim; j++)
 475         {
 476             mat2[j] = &(mat[ndim * j]);
 477             for (i = 0; i <= j; i++)
 478             {
 479                 if (mat2[j][i] < min)
 480                 {
 481                     min = mat2[j][i];
 482                 }
 483                 if (mat2[j][j] > max)
 484                 {
 485                     max = mat2[j][i];
 486                 }
 487             }
 488         }
 489         snew(axis, ndim);
 490         for (i = 0; i < ndim; i++)
 491         {
 492             axis[i] = i + 1;
 493         }
 494         rlo.r   = 0;
 495         rlo.g   = 0;
 496         rlo.b   = 1;
 497         rmi.r   = 1;
 498         rmi.g   = 1;
 499         rmi.b   = 1;
 500         rhi.r   = 1;
 501         rhi.g   = 0;
 502         rhi.b   = 0;
 503         out     = gmx_ffopen(xpmfile, "w");
 504         nlevels = 80;
 505         write_xpm3(out,
 506                    0,
 507                    "Covariance",
 508                    bM ? "u nm^2" : "nm^2",
 509                    "dim",
 510                    "dim",
 511                    ndim,
 512                    ndim,
 513                    axis,
 514                    axis,
 515                    mat2,
 516                    min,
 517                    0.0,
 518                    max,
 519                    rlo,
 520                    rmi,
 521                    rhi,
 522                    &nlevels);
 523         gmx_ffclose(out);
 524         sfree(axis);
 525         sfree(mat2);
 526     }
 527
 528     if (xpmafile)
 529     {
 530         min = 0;
 531         max = 0;
 532         snew(mat2, ndim / DIM);
 533         for (i = 0; i < ndim / DIM; i++)
 534         {
 535             snew(mat2[i], ndim / DIM);
 536         }
 537         for (j = 0; j < ndim / DIM; j++)
 538         {
 539             for (i = 0; i <= j; i++)
 540             {
 541                 mat2[j][i] = 0;
 542                 for (d = 0; d < DIM; d++)
 543                 {
 544                     mat2[j][i] += mat[ndim * (DIM * j + d) + DIM * i + d];
 545                 }
 546                 if (mat2[j][i] < min)
 547                 {
 548                     min = mat2[j][i];
 549                 }
 550                 if (mat2[j][j] > max)
 551                 {
 552                     max = mat2[j][i];
 553                 }
 554                 mat2[i][j] = mat2[j][i];
 555             }
 556         }
 557         snew(axis, ndim / DIM);
 558         for (i = 0; i < ndim / DIM; i++)
 559         {
 560             axis[i] = i + 1;
 561         }
 562         rlo.r   = 0;
 563         rlo.g   = 0;
 564         rlo.b   = 1;
 565         rmi.r   = 1;
 566         rmi.g   = 1;
 567         rmi.b   = 1;
 568         rhi.r   = 1;
 569         rhi.g   = 0;
 570         rhi.b   = 0;
 571         out     = gmx_ffopen(xpmafile, "w");
 572         nlevels = 80;
 573         write_xpm3(out,
 574                    0,
 575                    "Covariance",
 576                    bM ? "u nm^2" : "nm^2",
 577                    "atom",
 578                    "atom",
 579                    ndim / DIM,
 580                    ndim / DIM,
 581                    axis,
 582                    axis,
 583                    mat2,
 584                    min,
 585                    0.0,
 586                    max,
 587                    rlo,
 588                    rmi,
 589                    rhi,
 590                    &nlevels);
 591         gmx_ffclose(out);
 592         sfree(axis);
 593         for (i = 0; i < ndim / DIM; i++)
 594         {
 595             sfree(mat2[i]);
 596         }
 597         sfree(mat2);
 598     }
 599
 600
 601     /* call diagonalization routine */
 602
 603     snew(eigenvalues, ndim);
 604     snew(eigenvectors, ndim * ndim);
 605
 606     std::memcpy(eigenvectors, mat, ndim * ndim * sizeof(real));
 607     fprintf(stderr, "\nDiagonalizing ...\n");
 608     fflush(stderr);
 609     eigensolver(eigenvectors, ndim, 0, ndim, eigenvalues, mat);
 610     sfree(eigenvectors);
 611
 612     /* now write the output */
 613
 614     sum = 0;
 615     for (i = 0; i < ndim; i++)
 616     {
 617         sum += eigenvalues[i];
 618     }
 619     fprintf(stderr, "\nSum of the eigenvalues: %g (%snm^2)\n", sum, bM ? "u " : "");
 620     if (std::abs(trace - sum) > 0.01 * trace)
 621     {
 622         fprintf(stderr,
 623                 "\nWARNING: eigenvalue sum deviates from the trace of the covariance matrix\n");
 624     }
 625
 626     /* Set 'end', the maximum eigenvector and -value index used for output */
 627     if (end == -1)
 628     {
 629         if (nframes - 1 < ndim)
 630         {
 631             end = nframes - 1;
 632             fprintf(stderr, "\nWARNING: there are fewer frames in your trajectory than there are\n");
 633             fprintf(stderr, "degrees of freedom in your system. Only generating the first\n");
 634             fprintf(stderr, "%d out of %d eigenvectors and eigenvalues.\n", end, static_cast<int>(ndim));
 635         }
 636         else
 637         {
 638             end = ndim;
 639         }
 640     }
 641
 642     fprintf(stderr, "\nWriting eigenvalues to %s\n", eigvalfile);
 643
 644     sprintf(str, "(%snm\\S2\\N)", bM ? "u " : "");
 645     out = xvgropen(eigvalfile, "Eigenvalues of the covariance matrix", "Eigenvector index", str, oenv);
 646     for (i = 0; (i < end); i++)
 647     {
 648         fprintf(out, "%10d %g\n", static_cast<int>(i + 1), eigenvalues[ndim - 1 - i]);
 649     }
 650     xvgrclose(out);
 651
 652     if (bFit)
 653     {
 654         /* misuse lambda: 0/1 mass weighted analysis no/yes */
 655         if (nfit == natoms)
 656         {
 657             WriteXref = eWXR_YES;
 658             for (i = 0; i < nfit; i++)
 659             {
 660                 copy_rvec(xref[ifit[i]], x[i]);
 661             }
 662         }
 663         else
 664         {
 665             WriteXref = eWXR_NO;
 666         }
 667     }
 668     else
 669     {
 670         /* misuse lambda: -1 for no fit */
 671         WriteXref = eWXR_NOFIT;
 672     }
 673
 674     write_eigenvectors(
 675             eigvecfile, natoms, mat, TRUE, 1, end, WriteXref, x, bDiffMass1, xproj, bM, eigenvalues);
 676
 677     out = gmx_ffopen(logfile, "w");
 678
 679     fprintf(out, "Covariance analysis log, written %s\n", gmx_format_current_time().c_str());
 680
 681     fprintf(out, "Program: %s\n", argv[0]);
 682     gmx_getcwd(str, STRLEN);
 683
 684     fprintf(out, "Working directory: %s\n\n", str);
 685
 686     fprintf(out,
 687             "Read %d frames from %s (time %g to %g %s)\n",
 688             nframes,
 689             trxfile,
 690             output_env_conv_time(oenv, tstart),
 691             output_env_conv_time(oenv, tend),
 692             output_env_get_time_unit(oenv).c_str());
 693     if (bFit)
 694     {
 695         fprintf(out, "Read reference structure for fit from %s\n", fitfile);
 696     }
 697     if (ndxfile)
 698     {
 699         fprintf(out, "Read index groups from %s\n", ndxfile);
 700     }
 701     fprintf(out, "\n");
 702
 703     fprintf(out, "Analysis group is '%s' (%d atoms)\n", ananame, natoms);
 704     if (bFit)
 705     {
 706         fprintf(out, "Fit group is '%s' (%d atoms)\n", fitname, nfit);
 707     }
 708     else
 709     {
 710         fprintf(out, "No fit was used\n");
 711     }
 712     fprintf(out, "Analysis is %smass weighted\n", bDiffMass2 ? "" : "non-");
 713     if (bFit)
 714     {
 715         fprintf(out, "Fit is %smass weighted\n", bDiffMass1 ? "" : "non-");
 716     }
 717     fprintf(out, "Diagonalized the %dx%d covariance matrix\n", static_cast<int>(ndim), static_cast<int>(ndim));
 718     fprintf(out, "Trace of the covariance matrix before diagonalizing: %g\n", trace);
 719     fprintf(out, "Trace of the covariance matrix after diagonalizing: %g\n\n", sum);
 720
 721     fprintf(out, "Wrote %d eigenvalues to %s\n", static_cast<int>(end), eigvalfile);
 722     if (WriteXref == eWXR_YES)
 723     {
 724         fprintf(out, "Wrote reference structure to %s\n", eigvecfile);
 725     }
 726     fprintf(out, "Wrote average structure to %s and %s\n", averfile, eigvecfile);
 727     fprintf(out, "Wrote eigenvectors %d to %d to %s\n", 1, end, eigvecfile);
 728
 729     gmx_ffclose(out);
 730
 731     fprintf(stderr, "Wrote the log to %s\n", logfile);
 732
 733     return 0;
 734 }