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