src/gromacs/gmxana/gmx_covar.c

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