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