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