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