src/gromacs/gmxana/gmx_nmeig.c

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