src/gromacs/gmxana/gmx_nmeig.cpp

   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  * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include <cmath>
  40 #include <cstring>
  41
  42 #include "gromacs/commandline/pargs.h"
  43 #include "gromacs/fileio/mtxio.h"
  44 #include "gromacs/fileio/xvgr.h"
  45 #include "gromacs/gmxana/eigio.h"
  46 #include "gromacs/gmxana/gmx_ana.h"
  47 #include "gromacs/gmxana/gstat.h"
  48 #include "gromacs/legacyheaders/copyrite.h"
  49 #include "gromacs/legacyheaders/macros.h"
  50 #include "gromacs/legacyheaders/txtdump.h"
  51 #include "gromacs/legacyheaders/typedefs.h"
  52 #include "gromacs/linearalgebra/eigensolver.h"
  53 #include "gromacs/linearalgebra/sparsematrix.h"
  54 #include "gromacs/math/units.h"
  55 #include "gromacs/math/vec.h"
  56 #include "gromacs/topology/mtop_util.h"
  57 #include "gromacs/utility/futil.h"
  58 #include "gromacs/utility/gmxassert.h"
  59 #include "gromacs/utility/smalloc.h"
  60
  61 static double cv_corr(double nu, double T)
  62 {
  63     double x  = PLANCK*nu/(BOLTZ*T);
  64     double ex = std::exp(x);
  65
  66     if (nu <= 0)
  67     {
  68         return BOLTZ*KILO;
  69     }
  70     else
  71     {
  72         return BOLTZ*KILO*(ex*sqr(x)/sqr(ex-1) - 1);
  73     }
  74 }
  75
  76 static double u_corr(double nu, double T)
  77 {
  78     double x  = PLANCK*nu/(BOLTZ*T);
  79     double ex = std::exp(x);
  80
  81     if (nu <= 0)
  82     {
  83         return BOLTZ*T;
  84     }
  85     else
  86     {
  87         return BOLTZ*T*(0.5*x - 1 + x/(ex-1));
  88     }
  89 }
  90
  91 static int get_nharm_mt(gmx_moltype_t *mt)
  92 {
  93     static int harm_func[] = { F_BONDS };
  94     int        i, ft, nh;
  95
  96     nh = 0;
  97     for (i = 0; (i < asize(harm_func)); i++)
  98     {
  99         ft  = harm_func[i];
 100         nh += mt->ilist[ft].nr/(interaction_function[ft].nratoms+1);
 101     }
 102     return nh;
 103 }
 104
 105 static int get_nvsite_mt(gmx_moltype_t *mt)
 106 {
 107     static int vs_func[] = {
 108         F_VSITE2, F_VSITE3, F_VSITE3FD, F_VSITE3FAD,
 109         F_VSITE3OUT, F_VSITE4FD, F_VSITE4FDN, F_VSITEN
 110     };
 111     int        i, ft, nh;
 112
 113     nh = 0;
 114     for (i = 0; (i < asize(vs_func)); i++)
 115     {
 116         ft  = vs_func[i];
 117         nh += mt->ilist[ft].nr/(interaction_function[ft].nratoms+1);
 118     }
 119     return nh;
 120 }
 121
 122 static int get_nharm(gmx_mtop_t *mtop, int *nvsites)
 123 {
 124     int j, mt, nh, nv;
 125
 126     nh = 0;
 127     nv = 0;
 128     for (j = 0; (j < mtop->nmolblock); j++)
 129     {
 130         mt  = mtop->molblock[j].type;
 131         nh += mtop->molblock[j].nmol * get_nharm_mt(&(mtop->moltype[mt]));
 132         nv += mtop->molblock[j].nmol * get_nvsite_mt(&(mtop->moltype[mt]));
 133     }
 134     *nvsites = nv;
 135     return nh;
 136 }
 137
 138 static void
 139 nma_full_hessian(real     *           hess,
 140                  int                  ndim,
 141                  gmx_bool             bM,
 142                  t_topology     *     top,
 143                  int                  begin,
 144                  int                  end,
 145                  real     *           eigenvalues,
 146                  real     *           eigenvectors)
 147 {
 148     int  i, j, k, l;
 149     real mass_fac;
 150     int  natoms;
 151
 152     natoms = top->atoms.nr;
 153
 154     /* divide elements hess[i][j] by sqrt(mas[i])*sqrt(mas[j]) when required */
 155
 156     if (bM)
 157     {
 158         for (i = 0; (i < natoms); i++)
 159         {
 160             for (j = 0; (j < DIM); j++)
 161             {
 162                 for (k = 0; (k < natoms); k++)
 163                 {
 164                     mass_fac = gmx_invsqrt(top->atoms.atom[i].m*top->atoms.atom[k].m);
 165                     for (l = 0; (l < DIM); l++)
 166                     {
 167                         hess[(i*DIM+j)*ndim+k*DIM+l] *= mass_fac;
 168                     }
 169                 }
 170             }
 171         }
 172     }
 173
 174     /* call diagonalization routine. */
 175
 176     fprintf(stderr, "\nDiagonalizing to find vectors %d through %d...\n", begin, end);
 177     fflush(stderr);
 178
 179     eigensolver(hess, ndim, begin-1, end-1, eigenvalues, eigenvectors);
 180
 181     /* And scale the output eigenvectors */
 182     if (bM && eigenvectors != NULL)
 183     {
 184         for (i = 0; i < (end-begin+1); i++)
 185         {
 186             for (j = 0; j < natoms; j++)
 187             {
 188                 mass_fac = gmx_invsqrt(top->atoms.atom[j].m);
 189                 for (k = 0; (k < DIM); k++)
 190                 {
 191                     eigenvectors[i*ndim+j*DIM+k] *= mass_fac;
 192                 }
 193             }
 194         }
 195     }
 196 }
 197
 198
 199
 200 static void
 201 nma_sparse_hessian(gmx_sparsematrix_t     *     sparse_hessian,
 202                    gmx_bool                     bM,
 203                    t_topology     *             top,
 204                    int                          neig,
 205                    real     *                   eigenvalues,
 206                    real     *                   eigenvectors)
 207 {
 208     int  i, j, k;
 209     int  row, col;
 210     real mass_fac;
 211     int  iatom, katom;
 212     int  natoms;
 213     int  ndim;
 214
 215     natoms = top->atoms.nr;
 216     ndim   = DIM*natoms;
 217
 218     /* Cannot check symmetry since we only store half matrix */
 219     /* divide elements hess[i][j] by sqrt(mas[i])*sqrt(mas[j]) when required */
 220
 221     GMX_RELEASE_ASSERT(sparse_hessian != NULL, "NULL matrix pointer provided to nma_sparse_hessian");
 222
 223     if (bM)
 224     {
 225         for (iatom = 0; (iatom < natoms); iatom++)
 226         {
 227             for (j = 0; (j < DIM); j++)
 228             {
 229                 row = DIM*iatom+j;
 230                 for (k = 0; k < sparse_hessian->ndata[row]; k++)
 231                 {
 232                     col      = sparse_hessian->data[row][k].col;
 233                     katom    = col/3;
 234                     mass_fac = gmx_invsqrt(top->atoms.atom[iatom].m*top->atoms.atom[katom].m);
 235                     sparse_hessian->data[row][k].value *= mass_fac;
 236                 }
 237             }
 238         }
 239     }
 240     fprintf(stderr, "\nDiagonalizing to find eigenvectors 1 through %d...\n", neig);
 241     fflush(stderr);
 242
 243     sparse_eigensolver(sparse_hessian, neig, eigenvalues, eigenvectors, 10000000);
 244
 245     /* Scale output eigenvectors */
 246     if (bM && eigenvectors != NULL)
 247     {
 248         for (i = 0; i < neig; i++)
 249         {
 250             for (j = 0; j < natoms; j++)
 251             {
 252                 mass_fac = gmx_invsqrt(top->atoms.atom[j].m);
 253                 for (k = 0; (k < DIM); k++)
 254                 {
 255                     eigenvectors[i*ndim+j*DIM+k] *= mass_fac;
 256                 }
 257             }
 258         }
 259     }
 260 }
 261
 262
 263
 264 int gmx_nmeig(int argc, char *argv[])
 265 {
 266     const char            *desc[] = {
 267         "[THISMODULE] calculates the eigenvectors/values of a (Hessian) matrix,",
 268         "which can be calculated with [gmx-mdrun].",
 269         "The eigenvectors are written to a trajectory file ([TT]-v[tt]).",
 270         "The structure is written first with t=0. The eigenvectors",
 271         "are written as frames with the eigenvector number as timestamp.",
 272         "The eigenvectors can be analyzed with [gmx-anaeig].",
 273         "An ensemble of structures can be generated from the eigenvectors with",
 274         "[gmx-nmens]. When mass weighting is used, the generated eigenvectors",
 275         "will be scaled back to plain Cartesian coordinates before generating the",
 276         "output. In this case, they will no longer be exactly orthogonal in the",
 277         "standard Cartesian norm, but in the mass-weighted norm they would be.[PAR]",
 278         "This program can be optionally used to compute quantum corrections to heat capacity",
 279         "and enthalpy by providing an extra file argument [TT]-qcorr[tt]. See the GROMACS",
 280         "manual, Chapter 1, for details. The result includes subtracting a harmonic",
 281         "degree of freedom at the given temperature.",
 282         "The total correction is printed on the terminal screen.",
 283         "The recommended way of getting the corrections out is:[PAR]",
 284         "[TT]gmx nmeig -s topol.tpr -f nm.mtx -first 7 -last 10000 -T 300 -qc [-constr][tt][PAR]",
 285         "The [TT]-constr[tt] option should be used when bond constraints were used during the",
 286         "simulation [BB]for all the covalent bonds[bb]. If this is not the case, ",
 287         "you need to analyze the [TT]quant_corr.xvg[tt] file yourself.[PAR]",
 288         "To make things more flexible, the program can also take virtual sites into account",
 289         "when computing quantum corrections. When selecting [TT]-constr[tt] and",
 290         "[TT]-qc[tt], the [TT]-begin[tt] and [TT]-end[tt] options will be set automatically as well.",
 291         "Again, if you think you know it better, please check the [TT]eigenfreq.xvg[tt]",
 292         "output."
 293     };
 294
 295     static gmx_bool        bM    = TRUE, bCons = FALSE;
 296     static int             begin = 1, end = 50, maxspec = 4000;
 297     static real            T     = 298.15, width = 1;
 298     t_pargs                pa[]  =
 299     {
 300         { "-m",  FALSE, etBOOL, {&bM},
 301           "Divide elements of Hessian by product of sqrt(mass) of involved "
 302           "atoms prior to diagonalization. This should be used for 'Normal Modes' "
 303           "analysis" },
 304         { "-first", FALSE, etINT, {&begin},
 305           "First eigenvector to write away" },
 306         { "-last",  FALSE, etINT, {&end},
 307           "Last eigenvector to write away" },
 308         { "-maxspec", FALSE, etINT, {&maxspec},
 309           "Highest frequency (1/cm) to consider in the spectrum" },
 310         { "-T",     FALSE, etREAL, {&T},
 311           "Temperature for computing quantum heat capacity and enthalpy when using normal mode calculations to correct classical simulations" },
 312         { "-constr", FALSE, etBOOL, {&bCons},
 313           "If constraints were used in the simulation but not in the normal mode analysis (this is the recommended way of doing it) you will need to set this for computing the quantum corrections." },
 314         { "-width",  FALSE, etREAL, {&width},
 315           "Width (sigma) of the gaussian peaks (1/cm) when generating a spectrum" }
 316     };
 317     FILE                  *out, *qc, *spec;
 318     t_topology             top;
 319     gmx_mtop_t             mtop;
 320     rvec                  *top_x;
 321     matrix                 box;
 322     real                  *eigenvalues;
 323     real                  *eigenvectors;
 324     real                   qcvtot, qutot, qcv, qu;
 325     int                    natoms, ndim, nrow, ncol, nharm, nvsite;
 326     int                    i, j, k;
 327     gmx_bool               bSuck;
 328     t_tpxheader            tpx;
 329     int                    version, generation;
 330     real                   value, omega, nu;
 331     real                   factor_gmx_to_omega2;
 332     real                   factor_omega_to_wavenumber;
 333     real                  *spectrum = NULL;
 334     real                   wfac;
 335     output_env_t           oenv;
 336     const char            *qcleg[] = {
 337         "Heat Capacity cV (J/mol K)",
 338         "Enthalpy H (kJ/mol)"
 339     };
 340     real *                 full_hessian   = NULL;
 341     gmx_sparsematrix_t *   sparse_hessian = NULL;
 342
 343     t_filenm               fnm[] = {
 344         { efMTX, "-f", "hessian",    ffREAD  },
 345         { efTPR, NULL, NULL,         ffREAD  },
 346         { efXVG, "-of", "eigenfreq", ffWRITE },
 347         { efXVG, "-ol", "eigenval",  ffWRITE },
 348         { efXVG, "-os", "spectrum",  ffOPTWR },
 349         { efXVG, "-qc", "quant_corr",  ffOPTWR },
 350         { efTRN, "-v", "eigenvec",  ffWRITE }
 351     };
 352 #define NFILE asize(fnm)
 353
 354     if (!parse_common_args(&argc, argv, 0,
 355                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 356     {
 357         return 0;
 358     }
 359
 360     /* Read tpr file for volume and number of harmonic terms */
 361     read_tpxheader(ftp2fn(efTPR, NFILE, fnm), &tpx, TRUE, &version, &generation);
 362     snew(top_x, tpx.natoms);
 363
 364     read_tpx(ftp2fn(efTPR, NFILE, fnm), NULL, box, &natoms,
 365              top_x, NULL, NULL, &mtop);
 366     if (bCons)
 367     {
 368         nharm = get_nharm(&mtop, &nvsite);
 369     }
 370     else
 371     {
 372         nharm  = 0;
 373         nvsite = 0;
 374     }
 375     top = gmx_mtop_t_to_t_topology(&mtop);
 376
 377     bM   = TRUE;
 378     ndim = DIM*natoms;
 379
 380     if (opt2bSet("-qc", NFILE, fnm))
 381     {
 382         begin = 7+DIM*nvsite;
 383         end   = DIM*natoms;
 384     }
 385     if (begin < 1)
 386     {
 387         begin = 1;
 388     }
 389     if (end > ndim)
 390     {
 391         end = ndim;
 392     }
 393     printf("Using begin = %d and end = %d\n", begin, end);
 394
 395     /*open Hessian matrix */
 396     gmx_mtxio_read(ftp2fn(efMTX, NFILE, fnm), &nrow, &ncol, &full_hessian, &sparse_hessian);
 397
 398     /* Memory for eigenvalues and eigenvectors (begin..end) */
 399     snew(eigenvalues, nrow);
 400     snew(eigenvectors, nrow*(end-begin+1));
 401
 402     /* If the Hessian is in sparse format we can calculate max (ndim-1) eigenvectors,
 403      * and they must start at the first one. If this is not valid we convert to full matrix
 404      * storage, but warn the user that we might run out of memory...
 405      */
 406     if ((sparse_hessian != NULL) && (begin != 1 || end == ndim))
 407     {
 408         if (begin != 1)
 409         {
 410             fprintf(stderr, "Cannot use sparse Hessian with first eigenvector != 1.\n");
 411         }
 412         else if (end == ndim)
 413         {
 414             fprintf(stderr, "Cannot use sparse Hessian to calculate all eigenvectors.\n");
 415         }
 416
 417         fprintf(stderr, "Will try to allocate memory and convert to full matrix representation...\n");
 418
 419         snew(full_hessian, nrow*ncol);
 420         for (i = 0; i < nrow*ncol; i++)
 421         {
 422             full_hessian[i] = 0;
 423         }
 424
 425         for (i = 0; i < sparse_hessian->nrow; i++)
 426         {
 427             for (j = 0; j < sparse_hessian->ndata[i]; j++)
 428             {
 429                 k                      = sparse_hessian->data[i][j].col;
 430                 value                  = sparse_hessian->data[i][j].value;
 431                 full_hessian[i*ndim+k] = value;
 432                 full_hessian[k*ndim+i] = value;
 433             }
 434         }
 435         gmx_sparsematrix_destroy(sparse_hessian);
 436         sparse_hessian = NULL;
 437         fprintf(stderr, "Converted sparse to full matrix storage.\n");
 438     }
 439
 440     if (full_hessian != NULL)
 441     {
 442         /* Using full matrix storage */
 443         nma_full_hessian(full_hessian, nrow, bM, &top, begin, end,
 444                          eigenvalues, eigenvectors);
 445     }
 446     else
 447     {
 448         /* Sparse memory storage, allocate memory for eigenvectors */
 449         snew(eigenvectors, ncol*end);
 450         nma_sparse_hessian(sparse_hessian, bM, &top, end, eigenvalues, eigenvectors);
 451     }
 452
 453     /* check the output, first 6 eigenvalues should be reasonably small */
 454     bSuck = FALSE;
 455     for (i = begin-1; (i < 6); i++)
 456     {
 457         if (std::abs(eigenvalues[i]) > 1.0e-3)
 458         {
 459             bSuck = TRUE;
 460         }
 461     }
 462     if (bSuck)
 463     {
 464         fprintf(stderr, "\nOne of the lowest 6 eigenvalues has a non-zero value.\n");
 465         fprintf(stderr, "This could mean that the reference structure was not\n");
 466         fprintf(stderr, "properly energy minimized.\n");
 467     }
 468
 469     /* now write the output */
 470     fprintf (stderr, "Writing eigenvalues...\n");
 471     out = xvgropen(opt2fn("-ol", NFILE, fnm),
 472                    "Eigenvalues", "Eigenvalue index", "Eigenvalue [Gromacs units]",
 473                    oenv);
 474     if (output_env_get_print_xvgr_codes(oenv))
 475     {
 476         if (bM)
 477         {
 478             fprintf(out, "@ subtitle \"mass weighted\"\n");
 479         }
 480         else
 481         {
 482             fprintf(out, "@ subtitle \"not mass weighted\"\n");
 483         }
 484     }
 485
 486     for (i = 0; i <= (end-begin); i++)
 487     {
 488         fprintf (out, "%6d %15g\n", begin+i, eigenvalues[i]);
 489     }
 490     xvgrclose(out);
 491
 492
 493     if (opt2bSet("-qc", NFILE, fnm))
 494     {
 495         qc = xvgropen(opt2fn("-qc", NFILE, fnm), "Quantum Corrections", "Eigenvector index", "", oenv);
 496         xvgr_legend(qc, asize(qcleg), qcleg, oenv);
 497         qcvtot = qutot = 0;
 498     }
 499     else
 500     {
 501         qc = NULL;
 502     }
 503     printf("Writing eigenfrequencies - negative eigenvalues will be set to zero.\n");
 504
 505     out = xvgropen(opt2fn("-of", NFILE, fnm),
 506                    "Eigenfrequencies", "Eigenvector index", "Wavenumber [cm\\S-1\\N]",
 507                    oenv);
 508     if (output_env_get_print_xvgr_codes(oenv))
 509     {
 510         if (bM)
 511         {
 512             fprintf(out, "@ subtitle \"mass weighted\"\n");
 513         }
 514         else
 515         {
 516             fprintf(out, "@ subtitle \"not mass weighted\"\n");
 517         }
 518     }
 519     /* Spectrum ? */
 520     spec = NULL;
 521     if (opt2bSet("-os", NFILE, fnm) && (maxspec > 0))
 522     {
 523         snew(spectrum, maxspec);
 524         spec = xvgropen(opt2fn("-os", NFILE, fnm),
 525                         "Vibrational spectrum based on harmonic approximation",
 526                         "\\f{12}w\\f{4} (cm\\S-1\\N)",
 527                         "Intensity [Gromacs units]",
 528                         oenv);
 529         for (i = 0; (i < maxspec); i++)
 530         {
 531             spectrum[i] = 0;
 532         }
 533     }
 534
 535     /* Gromacs units are kJ/(mol*nm*nm*amu),
 536      * where amu is the atomic mass unit.
 537      *
 538      * For the eigenfrequencies we want to convert this to spectroscopic absorption
 539      * wavenumbers given in cm^(-1), which is the frequency divided by the speed of
 540      * light. Do this by first converting to omega^2 (units 1/s), take the square
 541      * root, and finally divide by the speed of light (nm/ps in gromacs).
 542      */
 543     factor_gmx_to_omega2       = 1.0E21/(AVOGADRO*AMU);
 544     factor_omega_to_wavenumber = 1.0E-5/(2.0*M_PI*SPEED_OF_LIGHT);
 545
 546     for (i = begin; (i <= end); i++)
 547     {
 548         value = eigenvalues[i-begin];
 549         if (value < 0)
 550         {
 551             value = 0;
 552         }
 553         omega = std::sqrt(value*factor_gmx_to_omega2);
 554         nu    = 1e-12*omega/(2*M_PI);
 555         value = omega*factor_omega_to_wavenumber;
 556         fprintf (out, "%6d %15g\n", i, value);
 557         if (NULL != spec)
 558         {
 559             wfac = eigenvalues[i-begin]/(width*std::sqrt(2*M_PI));
 560             for (j = 0; (j < maxspec); j++)
 561             {
 562                 spectrum[j] += wfac*std::exp(-sqr(j-value)/(2*sqr(width)));
 563             }
 564         }
 565         if (NULL != qc)
 566         {
 567             qcv = cv_corr(nu, T);
 568             qu  = u_corr(nu, T);
 569             if (i > end-nharm)
 570             {
 571                 qcv += BOLTZ*KILO;
 572                 qu  += BOLTZ*T;
 573             }
 574             fprintf (qc, "%6d %15g %15g\n", i, qcv, qu);
 575             qcvtot += qcv;
 576             qutot  += qu;
 577         }
 578     }
 579     xvgrclose(out);
 580     if (NULL != spec)
 581     {
 582         for (j = 0; (j < maxspec); j++)
 583         {
 584             fprintf(spec, "%10g  %10g\n", 1.0*j, spectrum[j]);
 585         }
 586         xvgrclose(spec);
 587     }
 588     if (NULL != qc)
 589     {
 590         printf("Quantum corrections for harmonic degrees of freedom\n");
 591         printf("Use appropriate -first and -last options to get reliable results.\n");
 592         printf("There were %d constraints and %d vsites in the simulation\n",
 593                nharm, nvsite);
 594         printf("Total correction to cV = %g J/mol K\n", qcvtot);
 595         printf("Total correction to  H = %g kJ/mol\n", qutot);
 596         xvgrclose(qc);
 597         please_cite(stdout, "Caleman2011b");
 598     }
 599     /* Writing eigenvectors. Note that if mass scaling was used, the eigenvectors
 600      * were scaled back from mass weighted cartesian to plain cartesian in the
 601      * nma_full_hessian() or nma_sparse_hessian() routines. Mass scaled vectors
 602      * will not be strictly orthogonal in plain cartesian scalar products.
 603      */
 604     write_eigenvectors(opt2fn("-v", NFILE, fnm), natoms, eigenvectors, FALSE, begin, end,
 605                        eWXR_NO, NULL, FALSE, top_x, bM, eigenvalues);
 606
 607     return 0;
 608 }