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/legacyheaders/typedefs.h"
  43 #include "gromacs/legacyheaders/macros.h"
  44 #include "gromacs/legacyheaders/copyrite.h"
  45 #include "gromacs/utility/futil.h"
  46 #include "gstat.h"
  47 #include "gromacs/legacyheaders/txtdump.h"
  48 #include "eigio.h"
  49 #include "gromacs/topology/mtop_util.h"
  50 #include "gromacs/math/units.h"
  51 #include "gmx_ana.h"
  52
  53 #include "gromacs/commandline/pargs.h"
  54 #include "gromacs/fileio/mtxio.h"
  55 #include "gromacs/fileio/xvgr.h"
  56 #include "gromacs/linearalgebra/eigensolver.h"
  57 #include "gromacs/linearalgebra/sparsematrix.h"
  58 #include "gromacs/math/vec.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 = 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 = 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, rdum;
 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     if (bM)
 222     {
 223         for (iatom = 0; (iatom < natoms); iatom++)
 224         {
 225             for (j = 0; (j < DIM); j++)
 226             {
 227                 row = DIM*iatom+j;
 228                 for (k = 0; k < sparse_hessian->ndata[row]; k++)
 229                 {
 230                     col      = sparse_hessian->data[row][k].col;
 231                     katom    = col/3;
 232                     mass_fac = gmx_invsqrt(top->atoms.atom[iatom].m*top->atoms.atom[katom].m);
 233                     sparse_hessian->data[row][k].value *= mass_fac;
 234                 }
 235             }
 236         }
 237     }
 238     fprintf(stderr, "\nDiagonalizing to find eigenvectors 1 through %d...\n", neig);
 239     fflush(stderr);
 240
 241     sparse_eigensolver(sparse_hessian, neig, eigenvalues, eigenvectors, 10000000);
 242
 243     /* Scale output eigenvectors */
 244     if (bM && eigenvectors != NULL)
 245     {
 246         for (i = 0; i < neig; i++)
 247         {
 248             for (j = 0; j < natoms; j++)
 249             {
 250                 mass_fac = gmx_invsqrt(top->atoms.atom[j].m);
 251                 for (k = 0; (k < DIM); k++)
 252                 {
 253                     eigenvectors[i*ndim+j*DIM+k] *= mass_fac;
 254                 }
 255             }
 256         }
 257     }
 258 }
 259
 260
 261
 262 int gmx_nmeig(int argc, char *argv[])
 263 {
 264     const char            *desc[] = {
 265         "[THISMODULE] calculates the eigenvectors/values of a (Hessian) matrix,",
 266         "which can be calculated with [gmx-mdrun].",
 267         "The eigenvectors are written to a trajectory file ([TT]-v[tt]).",
 268         "The structure is written first with t=0. The eigenvectors",
 269         "are written as frames with the eigenvector number as timestamp.",
 270         "The eigenvectors can be analyzed with [gmx-anaeig].",
 271         "An ensemble of structures can be generated from the eigenvectors with",
 272         "[gmx-nmens]. When mass weighting is used, the generated eigenvectors",
 273         "will be scaled back to plain Cartesian coordinates before generating the",
 274         "output. In this case, they will no longer be exactly orthogonal in the",
 275         "standard Cartesian norm, but in the mass-weighted norm they would be.[PAR]",
 276         "This program can be optionally used to compute quantum corrections to heat capacity",
 277         "and enthalpy by providing an extra file argument [TT]-qcorr[tt]. See the GROMACS",
 278         "manual, Chapter 1, for details. The result includes subtracting a harmonic",
 279         "degree of freedom at the given temperature.",
 280         "The total correction is printed on the terminal screen.",
 281         "The recommended way of getting the corrections out is:[PAR]",
 282         "[TT]gmx nmeig -s topol.tpr -f nm.mtx -first 7 -last 10000 -T 300 -qc [-constr][tt][PAR]",
 283         "The [TT]-constr[tt] option should be used when bond constraints were used during the",
 284         "simulation [BB]for all the covalent bonds[bb]. If this is not the case, ",
 285         "you need to analyze the [TT]quant_corr.xvg[tt] file yourself.[PAR]",
 286         "To make things more flexible, the program can also take virtual sites into account",
 287         "when computing quantum corrections. When selecting [TT]-constr[tt] and",
 288         "[TT]-qc[tt], the [TT]-begin[tt] and [TT]-end[tt] options will be set automatically as well.",
 289         "Again, if you think you know it better, please check the [TT]eigenfreq.xvg[tt]",
 290         "output."
 291     };
 292
 293     static gmx_bool        bM    = TRUE, bCons = FALSE;
 294     static int             begin = 1, end = 50, maxspec = 4000;
 295     static real            T     = 298.15, width = 1;
 296     t_pargs                pa[]  =
 297     {
 298         { "-m",  FALSE, etBOOL, {&bM},
 299           "Divide elements of Hessian by product of sqrt(mass) of involved "
 300           "atoms prior to diagonalization. This should be used for 'Normal Modes' "
 301           "analysis" },
 302         { "-first", FALSE, etINT, {&begin},
 303           "First eigenvector to write away" },
 304         { "-last",  FALSE, etINT, {&end},
 305           "Last eigenvector to write away" },
 306         { "-maxspec", FALSE, etINT, {&maxspec},
 307           "Highest frequency (1/cm) to consider in the spectrum" },
 308         { "-T",     FALSE, etREAL, {&T},
 309           "Temperature for computing quantum heat capacity and enthalpy when using normal mode calculations to correct classical simulations" },
 310         { "-constr", FALSE, etBOOL, {&bCons},
 311           "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." },
 312         { "-width",  FALSE, etREAL, {&width},
 313           "Width (sigma) of the gaussian peaks (1/cm) when generating a spectrum" }
 314     };
 315     FILE                  *out, *qc, *spec;
 316     int                    status, trjout;
 317     t_topology             top;
 318     gmx_mtop_t             mtop;
 319     int                    ePBC;
 320     rvec                  *top_x;
 321     matrix                 box;
 322     real                  *eigenvalues;
 323     real                  *eigenvectors;
 324     real                   rdum, mass_fac, qcvtot, qutot, qcv, qu;
 325     int                    natoms, ndim, nrow, ncol, count, nharm, nvsite;
 326     char                  *grpname;
 327     int                    i, j, k, l, d, gnx;
 328     gmx_bool               bSuck;
 329     atom_id               *index;
 330     t_tpxheader            tpx;
 331     int                    version, generation;
 332     real                   value, omega, nu;
 333     real                   factor_gmx_to_omega2;
 334     real                   factor_omega_to_wavenumber;
 335     real                  *spectrum = NULL;
 336     real                   wfac;
 337     output_env_t           oenv;
 338     const char            *qcleg[] = {
 339         "Heat Capacity cV (J/mol K)",
 340         "Enthalpy H (kJ/mol)"
 341     };
 342     real *                 full_hessian   = NULL;
 343     gmx_sparsematrix_t *   sparse_hessian = NULL;
 344
 345     t_filenm               fnm[] = {
 346         { efMTX, "-f", "hessian",    ffREAD  },
 347         { efTPR, NULL, NULL,         ffREAD  },
 348         { efXVG, "-of", "eigenfreq", ffWRITE },
 349         { efXVG, "-ol", "eigenval",  ffWRITE },
 350         { efXVG, "-os", "spectrum",  ffOPTWR },
 351         { efXVG, "-qc", "quant_corr",  ffOPTWR },
 352         { efTRN, "-v", "eigenvec",  ffWRITE }
 353     };
 354 #define NFILE asize(fnm)
 355
 356     if (!parse_common_args(&argc, argv, 0,
 357                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 358     {
 359         return 0;
 360     }
 361
 362     /* Read tpr file for volume and number of harmonic terms */
 363     read_tpxheader(ftp2fn(efTPR, NFILE, fnm), &tpx, TRUE, &version, &generation);
 364     snew(top_x, tpx.natoms);
 365
 366     read_tpx(ftp2fn(efTPR, NFILE, fnm), NULL, box, &natoms,
 367              top_x, NULL, NULL, &mtop);
 368     if (bCons)
 369     {
 370         nharm = get_nharm(&mtop, &nvsite);
 371     }
 372     else
 373     {
 374         nharm  = 0;
 375         nvsite = 0;
 376     }
 377     top = gmx_mtop_t_to_t_topology(&mtop);
 378
 379     bM   = TRUE;
 380     ndim = DIM*natoms;
 381
 382     if (opt2bSet("-qc", NFILE, fnm))
 383     {
 384         begin = 7+DIM*nvsite;
 385         end   = DIM*natoms;
 386     }
 387     if (begin < 1)
 388     {
 389         begin = 1;
 390     }
 391     if (end > ndim)
 392     {
 393         end = ndim;
 394     }
 395     printf("Using begin = %d and end = %d\n", begin, end);
 396
 397     /*open Hessian matrix */
 398     gmx_mtxio_read(ftp2fn(efMTX, NFILE, fnm), &nrow, &ncol, &full_hessian, &sparse_hessian);
 399
 400     /* Memory for eigenvalues and eigenvectors (begin..end) */
 401     snew(eigenvalues, nrow);
 402     snew(eigenvectors, nrow*(end-begin+1));
 403
 404     /* If the Hessian is in sparse format we can calculate max (ndim-1) eigenvectors,
 405      * and they must start at the first one. If this is not valid we convert to full matrix
 406      * storage, but warn the user that we might run out of memory...
 407      */
 408     if ((sparse_hessian != NULL) && (begin != 1 || end == ndim))
 409     {
 410         if (begin != 1)
 411         {
 412             fprintf(stderr, "Cannot use sparse Hessian with first eigenvector != 1.\n");
 413         }
 414         else if (end == ndim)
 415         {
 416             fprintf(stderr, "Cannot use sparse Hessian to calculate all eigenvectors.\n");
 417         }
 418
 419         fprintf(stderr, "Will try to allocate memory and convert to full matrix representation...\n");
 420
 421         snew(full_hessian, nrow*ncol);
 422         for (i = 0; i < nrow*ncol; i++)
 423         {
 424             full_hessian[i] = 0;
 425         }
 426
 427         for (i = 0; i < sparse_hessian->nrow; i++)
 428         {
 429             for (j = 0; j < sparse_hessian->ndata[i]; j++)
 430             {
 431                 k                      = sparse_hessian->data[i][j].col;
 432                 value                  = sparse_hessian->data[i][j].value;
 433                 full_hessian[i*ndim+k] = value;
 434                 full_hessian[k*ndim+i] = value;
 435             }
 436         }
 437         gmx_sparsematrix_destroy(sparse_hessian);
 438         sparse_hessian = NULL;
 439         fprintf(stderr, "Converted sparse to full matrix storage.\n");
 440     }
 441
 442     if (full_hessian != NULL)
 443     {
 444         /* Using full matrix storage */
 445         nma_full_hessian(full_hessian, nrow, bM, &top, begin, end,
 446                          eigenvalues, eigenvectors);
 447     }
 448     else
 449     {
 450         /* Sparse memory storage, allocate memory for eigenvectors */
 451         snew(eigenvectors, ncol*end);
 452         nma_sparse_hessian(sparse_hessian, bM, &top, end, eigenvalues, eigenvectors);
 453     }
 454
 455     /* check the output, first 6 eigenvalues should be reasonably small */
 456     bSuck = FALSE;
 457     for (i = begin-1; (i < 6); i++)
 458     {
 459         if (fabs(eigenvalues[i]) > 1.0e-3)
 460         {
 461             bSuck = TRUE;
 462         }
 463     }
 464     if (bSuck)
 465     {
 466         fprintf(stderr, "\nOne of the lowest 6 eigenvalues has a non-zero value.\n");
 467         fprintf(stderr, "This could mean that the reference structure was not\n");
 468         fprintf(stderr, "properly energy minimized.\n");
 469     }
 470
 471     /* now write the output */
 472     fprintf (stderr, "Writing eigenvalues...\n");
 473     out = xvgropen(opt2fn("-ol", NFILE, fnm),
 474                    "Eigenvalues", "Eigenvalue index", "Eigenvalue [Gromacs units]",
 475                    oenv);
 476     if (output_env_get_print_xvgr_codes(oenv))
 477     {
 478         if (bM)
 479         {
 480             fprintf(out, "@ subtitle \"mass weighted\"\n");
 481         }
 482         else
 483         {
 484             fprintf(out, "@ subtitle \"not mass weighted\"\n");
 485         }
 486     }
 487
 488     for (i = 0; i <= (end-begin); i++)
 489     {
 490         fprintf (out, "%6d %15g\n", begin+i, eigenvalues[i]);
 491     }
 492     gmx_ffclose(out);
 493
 494
 495     if (opt2bSet("-qc", NFILE, fnm))
 496     {
 497         qc = xvgropen(opt2fn("-qc", NFILE, fnm), "Quantum Corrections", "Eigenvector index", "", oenv);
 498         xvgr_legend(qc, asize(qcleg), qcleg, oenv);
 499         qcvtot = qutot = 0;
 500     }
 501     else
 502     {
 503         qc = NULL;
 504     }
 505     printf("Writing eigenfrequencies - negative eigenvalues will be set to zero.\n");
 506
 507     out = xvgropen(opt2fn("-of", NFILE, fnm),
 508                    "Eigenfrequencies", "Eigenvector index", "Wavenumber [cm\\S-1\\N]",
 509                    oenv);
 510     if (output_env_get_print_xvgr_codes(oenv))
 511     {
 512         if (bM)
 513         {
 514             fprintf(out, "@ subtitle \"mass weighted\"\n");
 515         }
 516         else
 517         {
 518             fprintf(out, "@ subtitle \"not mass weighted\"\n");
 519         }
 520     }
 521     /* Spectrum ? */
 522     spec = NULL;
 523     if (opt2bSet("-os", NFILE, fnm) && (maxspec > 0))
 524     {
 525         snew(spectrum, maxspec);
 526         spec = xvgropen(opt2fn("-os", NFILE, fnm),
 527                         "Vibrational spectrum based on harmonic approximation",
 528                         "\\f{12}w\\f{4} (cm\\S-1\\N)",
 529                         "Intensity [Gromacs units]",
 530                         oenv);
 531         for (i = 0; (i < maxspec); i++)
 532         {
 533             spectrum[i] = 0;
 534         }
 535     }
 536
 537     /* Gromacs units are kJ/(mol*nm*nm*amu),
 538      * where amu is the atomic mass unit.
 539      *
 540      * For the eigenfrequencies we want to convert this to spectroscopic absorption
 541      * wavenumbers given in cm^(-1), which is the frequency divided by the speed of
 542      * light. Do this by first converting to omega^2 (units 1/s), take the square
 543      * root, and finally divide by the speed of light (nm/ps in gromacs).
 544      */
 545     factor_gmx_to_omega2       = 1.0E21/(AVOGADRO*AMU);
 546     factor_omega_to_wavenumber = 1.0E-5/(2.0*M_PI*SPEED_OF_LIGHT);
 547
 548     for (i = begin; (i <= end); i++)
 549     {
 550         value = eigenvalues[i-begin];
 551         if (value < 0)
 552         {
 553             value = 0;
 554         }
 555         omega = sqrt(value*factor_gmx_to_omega2);
 556         nu    = 1e-12*omega/(2*M_PI);
 557         value = omega*factor_omega_to_wavenumber;
 558         fprintf (out, "%6d %15g\n", i, value);
 559         if (NULL != spec)
 560         {
 561             wfac = eigenvalues[i-begin]/(width*sqrt(2*M_PI));
 562             for (j = 0; (j < maxspec); j++)
 563             {
 564                 spectrum[j] += wfac*exp(-sqr(j-value)/(2*sqr(width)));
 565             }
 566         }
 567         if (NULL != qc)
 568         {
 569             qcv = cv_corr(nu, T);
 570             qu  = u_corr(nu, T);
 571             if (i > end-nharm)
 572             {
 573                 qcv += BOLTZ*KILO;
 574                 qu  += BOLTZ*T;
 575             }
 576             fprintf (qc, "%6d %15g %15g\n", i, qcv, qu);
 577             qcvtot += qcv;
 578             qutot  += qu;
 579         }
 580     }
 581     gmx_ffclose(out);
 582     if (NULL != spec)
 583     {
 584         for (j = 0; (j < maxspec); j++)
 585         {
 586             fprintf(spec, "%10g  %10g\n", 1.0*j, spectrum[j]);
 587         }
 588         gmx_ffclose(spec);
 589     }
 590     if (NULL != qc)
 591     {
 592         printf("Quantum corrections for harmonic degrees of freedom\n");
 593         printf("Use appropriate -first and -last options to get reliable results.\n");
 594         printf("There were %d constraints and %d vsites in the simulation\n",
 595                nharm, nvsite);
 596         printf("Total correction to cV = %g J/mol K\n", qcvtot);
 597         printf("Total correction to  H = %g kJ/mol\n", qutot);
 598         gmx_ffclose(qc);
 599         please_cite(stdout, "Caleman2011b");
 600     }
 601     /* Writing eigenvectors. Note that if mass scaling was used, the eigenvectors
 602      * were scaled back from mass weighted cartesian to plain cartesian in the
 603      * nma_full_hessian() or nma_sparse_hessian() routines. Mass scaled vectors
 604      * will not be strictly orthogonal in plain cartesian scalar products.
 605      */
 606     write_eigenvectors(opt2fn("-v", NFILE, fnm), natoms, eigenvectors, FALSE, begin, end,
 607                        eWXR_NO, NULL, FALSE, top_x, bM, eigenvalues);
 608
 609     return 0;
 610 }