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