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