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