src/gromacs/gmxana/gmx_nmeig.cpp

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