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     if (!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         return 0;
 169     }
 170
 171     clear_mat(zerobox);
 172
 173     fitfile    = ftp2fn(efTPS, NFILE, fnm);
 174     trxfile    = ftp2fn(efTRX, NFILE, fnm);
 175     ndxfile    = ftp2fn_null(efNDX, NFILE, fnm);
 176     eigvalfile = ftp2fn(efXVG, NFILE, fnm);
 177     eigvecfile = ftp2fn(efTRN, NFILE, fnm);
 178     averfile   = ftp2fn(efSTO, NFILE, fnm);
 179     logfile    = ftp2fn(efLOG, NFILE, fnm);
 180     asciifile  = opt2fn_null("-ascii", NFILE, fnm);
 181     xpmfile    = opt2fn_null("-xpm", NFILE, fnm);
 182     xpmafile   = opt2fn_null("-xpma", NFILE, fnm);
 183
 184     read_tps_conf(fitfile, str, &top, &ePBC, &xref, NULL, box, TRUE);
 185     atoms = &top.atoms;
 186
 187     if (bFit)
 188     {
 189         printf("\nChoose a group for the least squares fit\n");
 190         get_index(atoms, ndxfile, 1, &nfit, &ifit, &fitname);
 191         if (nfit < 3)
 192         {
 193             gmx_fatal(FARGS, "Need >= 3 points to fit!\n");
 194         }
 195     }
 196     else
 197     {
 198         nfit = 0;
 199     }
 200     printf("\nChoose a group for the covariance analysis\n");
 201     get_index(atoms, ndxfile, 1, &natoms, &index, &ananame);
 202
 203     bDiffMass1 = FALSE;
 204     if (bFit)
 205     {
 206         snew(w_rls, atoms->nr);
 207         for (i = 0; (i < nfit); i++)
 208         {
 209             w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
 210             if (i)
 211             {
 212                 bDiffMass1 = bDiffMass1 || (w_rls[ifit[i]] != w_rls[ifit[i-1]]);
 213             }
 214         }
 215     }
 216     bDiffMass2 = FALSE;
 217     snew(sqrtm, natoms);
 218     for (i = 0; (i < natoms); i++)
 219     {
 220         if (bM)
 221         {
 222             sqrtm[i] = sqrt(atoms->atom[index[i]].m);
 223             if (i)
 224             {
 225                 bDiffMass2 = bDiffMass2 || (sqrtm[i] != sqrtm[i-1]);
 226             }
 227         }
 228         else
 229         {
 230             sqrtm[i] = 1.0;
 231         }
 232     }
 233
 234     if (bFit && bDiffMass1 && !bDiffMass2)
 235     {
 236         bDiffMass1 = natoms != nfit;
 237         i          = 0;
 238         for (i = 0; (i < natoms) && !bDiffMass1; i++)
 239         {
 240             bDiffMass1 = index[i] != ifit[i];
 241         }
 242         if (!bDiffMass1)
 243         {
 244             fprintf(stderr, "\n"
 245                     "Note: the fit and analysis group are identical,\n"
 246                     "      while the fit is mass weighted and the analysis is not.\n"
 247                     "      Making the fit non mass weighted.\n\n");
 248             for (i = 0; (i < nfit); i++)
 249             {
 250                 w_rls[ifit[i]] = 1.0;
 251             }
 252         }
 253     }
 254
 255     /* Prepare reference frame */
 256     if (bPBC)
 257     {
 258         gpbc = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr);
 259         gmx_rmpbc(gpbc, atoms->nr, box, xref);
 260     }
 261     if (bFit)
 262     {
 263         reset_x(nfit, ifit, atoms->nr, NULL, xref, w_rls);
 264     }
 265
 266     snew(x, natoms);
 267     snew(xav, natoms);
 268     ndim = natoms*DIM;
 269     if (sqrt(GMX_LARGE_INT_MAX) < ndim)
 270     {
 271         gmx_fatal(FARGS, "Number of degrees of freedoms to large for matrix.\n");
 272     }
 273     snew(mat, ndim*ndim);
 274
 275     fprintf(stderr, "Calculating the average structure ...\n");
 276     nframes0 = 0;
 277     nat      = read_first_x(oenv, &status, trxfile, &t, &xread, box);
 278     if (nat != atoms->nr)
 279     {
 280         fprintf(stderr, "\nWARNING: number of atoms in tpx (%d) and trajectory (%d) do not match\n", natoms, nat);
 281     }
 282     do
 283     {
 284         nframes0++;
 285         /* calculate x: a fitted struture of the selected atoms */
 286         if (bPBC)
 287         {
 288             gmx_rmpbc(gpbc, nat, box, xread);
 289         }
 290         if (bFit)
 291         {
 292             reset_x(nfit, ifit, nat, NULL, xread, w_rls);
 293             do_fit(nat, w_rls, xref, xread);
 294         }
 295         for (i = 0; i < natoms; i++)
 296         {
 297             rvec_inc(xav[i], xread[index[i]]);
 298         }
 299     }
 300     while (read_next_x(oenv, status, &t, xread, box));
 301     close_trj(status);
 302
 303     inv_nframes = 1.0/nframes0;
 304     for (i = 0; i < natoms; i++)
 305     {
 306         for (d = 0; d < DIM; d++)
 307         {
 308             xav[i][d]         *= inv_nframes;
 309             xread[index[i]][d] = xav[i][d];
 310         }
 311     }
 312     write_sto_conf_indexed(opt2fn("-av", NFILE, fnm), "Average structure",
 313                            atoms, xread, NULL, epbcNONE, zerobox, natoms, index);
 314     sfree(xread);
 315
 316     fprintf(stderr, "Constructing covariance matrix (%dx%d) ...\n", (int)ndim, (int)ndim);
 317     nframes = 0;
 318     nat     = read_first_x(oenv, &status, trxfile, &t, &xread, box);
 319     tstart  = t;
 320     do
 321     {
 322         nframes++;
 323         tend = t;
 324         /* calculate x: a (fitted) structure of the selected atoms */
 325         if (bPBC)
 326         {
 327             gmx_rmpbc(gpbc, nat, box, xread);
 328         }
 329         if (bFit)
 330         {
 331             reset_x(nfit, ifit, nat, NULL, xread, w_rls);
 332             do_fit(nat, w_rls, xref, xread);
 333         }
 334         if (bRef)
 335         {
 336             for (i = 0; i < natoms; i++)
 337             {
 338                 rvec_sub(xread[index[i]], xref[index[i]], x[i]);
 339             }
 340         }
 341         else
 342         {
 343             for (i = 0; i < natoms; i++)
 344             {
 345                 rvec_sub(xread[index[i]], xav[i], x[i]);
 346             }
 347         }
 348
 349         for (j = 0; j < natoms; j++)
 350         {
 351             for (dj = 0; dj < DIM; dj++)
 352             {
 353                 k  = ndim*(DIM*j+dj);
 354                 xj = x[j][dj];
 355                 for (i = j; i < natoms; i++)
 356                 {
 357                     l = k+DIM*i;
 358                     for (d = 0; d < DIM; d++)
 359                     {
 360                         mat[l+d] += x[i][d]*xj;
 361                     }
 362                 }
 363             }
 364         }
 365     }
 366     while (read_next_x(oenv, status, &t, xread, box) &&
 367            (bRef || nframes < nframes0));
 368     close_trj(status);
 369     gmx_rmpbc_done(gpbc);
 370
 371     fprintf(stderr, "Read %d frames\n", nframes);
 372
 373     if (bRef)
 374     {
 375         /* copy the reference structure to the ouput array x */
 376         snew(xproj, natoms);
 377         for (i = 0; i < natoms; i++)
 378         {
 379             copy_rvec(xref[index[i]], xproj[i]);
 380         }
 381     }
 382     else
 383     {
 384         xproj = xav;
 385     }
 386
 387     /* correct the covariance matrix for the mass */
 388     inv_nframes = 1.0/nframes;
 389     for (j = 0; j < natoms; j++)
 390     {
 391         for (dj = 0; dj < DIM; dj++)
 392         {
 393             for (i = j; i < natoms; i++)
 394             {
 395                 k = ndim*(DIM*j+dj)+DIM*i;
 396                 for (d = 0; d < DIM; d++)
 397                 {
 398                     mat[k+d] = mat[k+d]*inv_nframes*sqrtm[i]*sqrtm[j];
 399                 }
 400             }
 401         }
 402     }
 403
 404     /* symmetrize the matrix */
 405     for (j = 0; j < ndim; j++)
 406     {
 407         for (i = j; i < ndim; i++)
 408         {
 409             mat[ndim*i+j] = mat[ndim*j+i];
 410         }
 411     }
 412
 413     trace = 0;
 414     for (i = 0; i < ndim; i++)
 415     {
 416         trace += mat[i*ndim+i];
 417     }
 418     fprintf(stderr, "\nTrace of the covariance matrix: %g (%snm^2)\n",
 419             trace, bM ? "u " : "");
 420
 421     if (asciifile)
 422     {
 423         out = ffopen(asciifile, "w");
 424         for (j = 0; j < ndim; j++)
 425         {
 426             for (i = 0; i < ndim; i += 3)
 427             {
 428                 fprintf(out, "%g %g %g\n",
 429                         mat[ndim*j+i], mat[ndim*j+i+1], mat[ndim*j+i+2]);
 430             }
 431         }
 432         ffclose(out);
 433     }
 434
 435     if (xpmfile)
 436     {
 437         min = 0;
 438         max = 0;
 439         snew(mat2, ndim);
 440         for (j = 0; j < ndim; j++)
 441         {
 442             mat2[j] = &(mat[ndim*j]);
 443             for (i = 0; i <= j; i++)
 444             {
 445                 if (mat2[j][i] < min)
 446                 {
 447                     min = mat2[j][i];
 448                 }
 449                 if (mat2[j][j] > max)
 450                 {
 451                     max = mat2[j][i];
 452                 }
 453             }
 454         }
 455         snew(axis, ndim);
 456         for (i = 0; i < ndim; i++)
 457         {
 458             axis[i] = i+1;
 459         }
 460         rlo.r   = 0; rlo.g = 0; rlo.b = 1;
 461         rmi.r   = 1; rmi.g = 1; rmi.b = 1;
 462         rhi.r   = 1; rhi.g = 0; rhi.b = 0;
 463         out     = ffopen(xpmfile, "w");
 464         nlevels = 80;
 465         write_xpm3(out, 0, "Covariance", bM ? "u nm^2" : "nm^2",
 466                    "dim", "dim", ndim, ndim, axis, axis,
 467                    mat2, min, 0.0, max, rlo, rmi, rhi, &nlevels);
 468         ffclose(out);
 469         sfree(axis);
 470         sfree(mat2);
 471     }
 472
 473     if (xpmafile)
 474     {
 475         min = 0;
 476         max = 0;
 477         snew(mat2, ndim/DIM);
 478         for (i = 0; i < ndim/DIM; i++)
 479         {
 480             snew(mat2[i], ndim/DIM);
 481         }
 482         for (j = 0; j < ndim/DIM; j++)
 483         {
 484             for (i = 0; i <= j; i++)
 485             {
 486                 mat2[j][i] = 0;
 487                 for (d = 0; d < DIM; d++)
 488                 {
 489                     mat2[j][i] += mat[ndim*(DIM*j+d)+DIM*i+d];
 490                 }
 491                 if (mat2[j][i] < min)
 492                 {
 493                     min = mat2[j][i];
 494                 }
 495                 if (mat2[j][j] > max)
 496                 {
 497                     max = mat2[j][i];
 498                 }
 499                 mat2[i][j] = mat2[j][i];
 500             }
 501         }
 502         snew(axis, ndim/DIM);
 503         for (i = 0; i < ndim/DIM; i++)
 504         {
 505             axis[i] = i+1;
 506         }
 507         rlo.r   = 0; rlo.g = 0; rlo.b = 1;
 508         rmi.r   = 1; rmi.g = 1; rmi.b = 1;
 509         rhi.r   = 1; rhi.g = 0; rhi.b = 0;
 510         out     = ffopen(xpmafile, "w");
 511         nlevels = 80;
 512         write_xpm3(out, 0, "Covariance", bM ? "u nm^2" : "nm^2",
 513                    "atom", "atom", ndim/DIM, ndim/DIM, axis, axis,
 514                    mat2, min, 0.0, max, rlo, rmi, rhi, &nlevels);
 515         ffclose(out);
 516         sfree(axis);
 517         for (i = 0; i < ndim/DIM; i++)
 518         {
 519             sfree(mat2[i]);
 520         }
 521         sfree(mat2);
 522     }
 523
 524
 525     /* call diagonalization routine */
 526
 527     snew(eigenvalues, ndim);
 528     snew(eigenvectors, ndim*ndim);
 529
 530     memcpy(eigenvectors, mat, ndim*ndim*sizeof(real));
 531     fprintf(stderr, "\nDiagonalizing ...\n");
 532     fflush(stderr);
 533     eigensolver(eigenvectors, ndim, 0, ndim, eigenvalues, mat);
 534     sfree(eigenvectors);
 535
 536     /* now write the output */
 537
 538     sum = 0;
 539     for (i = 0; i < ndim; i++)
 540     {
 541         sum += eigenvalues[i];
 542     }
 543     fprintf(stderr, "\nSum of the eigenvalues: %g (%snm^2)\n",
 544             sum, bM ? "u " : "");
 545     if (fabs(trace-sum) > 0.01*trace)
 546     {
 547         fprintf(stderr, "\nWARNING: eigenvalue sum deviates from the trace of the covariance matrix\n");
 548     }
 549
 550     fprintf(stderr, "\nWriting eigenvalues to %s\n", eigvalfile);
 551
 552     sprintf(str, "(%snm\\S2\\N)", bM ? "u " : "");
 553     out = xvgropen(eigvalfile,
 554                    "Eigenvalues of the covariance matrix",
 555                    "Eigenvector index", str, oenv);
 556     for (i = 0; (i < ndim); i++)
 557     {
 558         fprintf (out, "%10d %g\n", (int)i+1, eigenvalues[ndim-1-i]);
 559     }
 560     ffclose(out);
 561
 562     if (end == -1)
 563     {
 564         if (nframes-1 < ndim)
 565         {
 566             end = nframes-1;
 567         }
 568         else
 569         {
 570             end = ndim;
 571         }
 572     }
 573     if (bFit)
 574     {
 575         /* misuse lambda: 0/1 mass weighted analysis no/yes */
 576         if (nfit == natoms)
 577         {
 578             WriteXref = eWXR_YES;
 579             for (i = 0; i < nfit; i++)
 580             {
 581                 copy_rvec(xref[ifit[i]], x[i]);
 582             }
 583         }
 584         else
 585         {
 586             WriteXref = eWXR_NO;
 587         }
 588     }
 589     else
 590     {
 591         /* misuse lambda: -1 for no fit */
 592         WriteXref = eWXR_NOFIT;
 593     }
 594
 595     write_eigenvectors(eigvecfile, natoms, mat, TRUE, 1, end,
 596                        WriteXref, x, bDiffMass1, xproj, bM, eigenvalues);
 597
 598     out = ffopen(logfile, "w");
 599
 600     time(&now);
 601     gmx_ctime_r(&now, timebuf, STRLEN);
 602     fprintf(out, "Covariance analysis log, written %s\n", timebuf);
 603
 604     fprintf(out, "Program: %s\n", argv[0]);
 605     gmx_getcwd(str, STRLEN);
 606
 607     fprintf(out, "Working directory: %s\n\n", str);
 608
 609     fprintf(out, "Read %d frames from %s (time %g to %g %s)\n", nframes, trxfile,
 610             output_env_conv_time(oenv, tstart), output_env_conv_time(oenv, tend), output_env_get_time_unit(oenv));
 611     if (bFit)
 612     {
 613         fprintf(out, "Read reference structure for fit from %s\n", fitfile);
 614     }
 615     if (ndxfile)
 616     {
 617         fprintf(out, "Read index groups from %s\n", ndxfile);
 618     }
 619     fprintf(out, "\n");
 620
 621     fprintf(out, "Analysis group is '%s' (%d atoms)\n", ananame, natoms);
 622     if (bFit)
 623     {
 624         fprintf(out, "Fit group is '%s' (%d atoms)\n", fitname, nfit);
 625     }
 626     else
 627     {
 628         fprintf(out, "No fit was used\n");
 629     }
 630     fprintf(out, "Analysis is %smass weighted\n", bDiffMass2 ? "" : "non-");
 631     if (bFit)
 632     {
 633         fprintf(out, "Fit is %smass weighted\n", bDiffMass1 ? "" : "non-");
 634     }
 635     fprintf(out, "Diagonalized the %dx%d covariance matrix\n", (int)ndim, (int)ndim);
 636     fprintf(out, "Trace of the covariance matrix before diagonalizing: %g\n",
 637             trace);
 638     fprintf(out, "Trace of the covariance matrix after diagonalizing: %g\n\n",
 639             sum);
 640
 641     fprintf(out, "Wrote %d eigenvalues to %s\n", (int)ndim, eigvalfile);
 642     if (WriteXref == eWXR_YES)
 643     {
 644         fprintf(out, "Wrote reference structure to %s\n", eigvecfile);
 645     }
 646     fprintf(out, "Wrote average structure to %s and %s\n", averfile, eigvecfile);
 647     fprintf(out, "Wrote eigenvectors %d to %d to %s\n", 1, end, eigvecfile);
 648
 649     ffclose(out);
 650
 651     fprintf(stderr, "Wrote the log to %s\n", logfile);
 652
 653     return 0;
 654 }