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