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