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