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