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