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,2019, 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].size()/(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                  gmx::ArrayRef<const int> 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 (int i = 0; (i < atom_index.ssize()); i++)
 144         {
 145             size_t ai = atom_index[i];
 146             for (size_t j = 0; (j < DIM); j++)
 147             {
 148                 for (int k = 0; (k < atom_index.ssize()); 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 (int j = 0; j < atom_index.ssize(); 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                    gmx::ArrayRef<const int> 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 (int iatom = 0; (iatom < atom_index.ssize()); 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 (int j = 0; j < atom_index.ssize(); 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                                    gmx::ArrayRef<const int> 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     double ZPE    = calcZeroPointEnergy(eFreq, scale_factor);
 410     fprintf(fp, "Zero-point energy     %g kJ/mol\n", ZPE);
 411 }
 412
 413
 414 int gmx_nmeig(int argc, char *argv[])
 415 {
 416     const char            *desc[] = {
 417         "[THISMODULE] calculates the eigenvectors/values of a (Hessian) matrix,",
 418         "which can be calculated with [gmx-mdrun].",
 419         "The eigenvectors are written to a trajectory file ([TT]-v[tt]).",
 420         "The structure is written first with t=0. The eigenvectors",
 421         "are written as frames with the eigenvector number and eigenvalue",
 422         "written as step number and timestamp, respectively.",
 423         "The eigenvectors can be analyzed with [gmx-anaeig].",
 424         "An ensemble of structures can be generated from the eigenvectors with",
 425         "[gmx-nmens]. When mass weighting is used, the generated eigenvectors",
 426         "will be scaled back to plain Cartesian coordinates before generating the",
 427         "output. In this case, they will no longer be exactly orthogonal in the",
 428         "standard Cartesian norm, but in the mass-weighted norm they would be.[PAR]",
 429         "This program can be optionally used to compute quantum corrections to heat capacity",
 430         "and enthalpy by providing an extra file argument [TT]-qcorr[tt]. See the GROMACS",
 431         "manual, Chapter 1, for details. The result includes subtracting a harmonic",
 432         "degree of freedom at the given temperature.",
 433         "The total correction is printed on the terminal screen.",
 434         "The recommended way of getting the corrections out is:[PAR]",
 435         "[TT]gmx nmeig -s topol.tpr -f nm.mtx -first 7 -last 10000 -T 300 -qc [-constr][tt][PAR]",
 436         "The [TT]-constr[tt] option should be used when bond constraints were used during the",
 437         "simulation [BB]for all the covalent bonds[bb]. If this is not the case, ",
 438         "you need to analyze the [TT]quant_corr.xvg[tt] file yourself.[PAR]",
 439         "To make things more flexible, the program can also take virtual sites into account",
 440         "when computing quantum corrections. When selecting [TT]-constr[tt] and",
 441         "[TT]-qc[tt], the [TT]-begin[tt] and [TT]-end[tt] options will be set automatically as well.[PAR]",
 442         "Based on a harmonic analysis of the normal mode frequencies,",
 443         "thermochemical properties S0 (Standard Entropy),",
 444         "Cv (Heat capacity at constant volume), Zero-point energy and the internal energy are",
 445         "computed, much in the same manner as popular quantum chemistry",
 446         "programs."
 447     };
 448
 449     static gmx_bool        bM           = TRUE, bCons = FALSE;
 450     static int             begin        = 1, end = 50, maxspec = 4000, sigma_r = 1;
 451     static real            T            = 298.15, width = 1, P = 1, scale_factor = 1;
 452     static real            linear_toler = 1e-5;
 453
 454     t_pargs                pa[]  =
 455     {
 456         { "-m",  FALSE, etBOOL, {&bM},
 457           "Divide elements of Hessian by product of sqrt(mass) of involved "
 458           "atoms prior to diagonalization. This should be used for 'Normal Modes' "
 459           "analysis" },
 460         { "-first", FALSE, etINT, {&begin},
 461           "First eigenvector to write away" },
 462         { "-last",  FALSE, etINT, {&end},
 463           "Last eigenvector to write away. -1 is use all dimensions." },
 464         { "-maxspec", FALSE, etINT, {&maxspec},
 465           "Highest frequency (1/cm) to consider in the spectrum" },
 466         { "-T",     FALSE, etREAL, {&T},
 467           "Temperature for computing entropy, quantum heat capacity and enthalpy when using normal mode calculations to correct classical simulations" },
 468         { "-P",     FALSE, etREAL, {&P},
 469           "Pressure (bar) when computing entropy" },
 470         { "-sigma", FALSE, etINT, {&sigma_r},
 471           "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." },
 472         { "-scale", FALSE, etREAL, {&scale_factor},
 473           "Factor to scale frequencies before computing thermochemistry values" },
 474         { "-linear_toler", FALSE, etREAL, {&linear_toler},
 475           "Tolerance for determining whether a compound is linear as determined from the ration of the moments inertion Ix/Iy and Ix/Iz." },
 476         { "-constr", FALSE, etBOOL, {&bCons},
 477           "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." },
 478         { "-width",  FALSE, etREAL, {&width},
 479           "Width (sigma) of the gaussian peaks (1/cm) when generating a spectrum" }
 480     };
 481     FILE                  *out, *qc, *spec;
 482     t_topology             top;
 483     gmx_mtop_t             mtop;
 484     rvec                  *top_x;
 485     matrix                 box;
 486     real                  *eigenvalues;
 487     real                  *eigenvectors;
 488     real                   qcvtot, qutot, qcv, qu;
 489     int                    i, j, k;
 490     t_tpxheader            tpx;
 491     real                   value, omega, nu;
 492     real                   factor_gmx_to_omega2;
 493     real                   factor_omega_to_wavenumber;
 494     real                  *spectrum = nullptr;
 495     real                   wfac;
 496     gmx_output_env_t      *oenv;
 497     const char            *qcleg[] = {
 498         "Heat Capacity cV (J/mol K)",
 499         "Enthalpy H (kJ/mol)"
 500     };
 501     real *                 full_hessian   = nullptr;
 502     gmx_sparsematrix_t *   sparse_hessian = nullptr;
 503
 504     t_filenm               fnm[] = {
 505         { efMTX, "-f", "hessian",    ffREAD  },
 506         { efTPR, nullptr, nullptr,         ffREAD  },
 507         { efXVG, "-of", "eigenfreq", ffWRITE },
 508         { efXVG, "-ol", "eigenval",  ffWRITE },
 509         { efXVG, "-os", "spectrum",  ffOPTWR },
 510         { efXVG, "-qc", "quant_corr",  ffOPTWR },
 511         { efTRN, "-v", "eigenvec",  ffWRITE }
 512     };
 513 #define NFILE asize(fnm)
 514
 515     if (!parse_common_args(&argc, argv, 0,
 516                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
 517     {
 518         return 0;
 519     }
 520
 521     /* Read tpr file for volume and number of harmonic terms */
 522     read_tpxheader(ftp2fn(efTPR, NFILE, fnm), &tpx, TRUE);
 523     snew(top_x, tpx.natoms);
 524
 525     int natoms_tpx;
 526     read_tpx(ftp2fn(efTPR, NFILE, fnm), nullptr, box, &natoms_tpx,
 527              top_x, nullptr, &mtop);
 528     int nharm = 0;
 529     if (bCons)
 530     {
 531         nharm = get_nharm(&mtop);
 532     }
 533     std::vector<int> atom_index = get_atom_index(&mtop);
 534
 535     top = gmx_mtop_t_to_t_topology(&mtop, true);
 536
 537     bM       = TRUE;
 538     int ndim = DIM*atom_index.size();
 539
 540     if (opt2bSet("-qc", NFILE, fnm))
 541     {
 542         begin = 7;
 543         end   = ndim;
 544     }
 545     if (begin < 1)
 546     {
 547         begin = 1;
 548     }
 549     if (end == -1 || end > ndim)
 550     {
 551         end = ndim;
 552     }
 553     printf("Using begin = %d and end = %d\n", begin, end);
 554
 555     /*open Hessian matrix */
 556     int nrow, ncol;
 557     gmx_mtxio_read(ftp2fn(efMTX, NFILE, fnm), &nrow, &ncol, &full_hessian, &sparse_hessian);
 558
 559     /* If the Hessian is in sparse format we can calculate max (ndim-1) eigenvectors,
 560      * If this is not valid we convert to full matrix storage,
 561      * but warn the user that we might run out of memory...
 562      */
 563     if ((sparse_hessian != nullptr) && (end == ndim))
 564     {
 565         fprintf(stderr, "Cannot use sparse Hessian to calculate all eigenvectors.\n");
 566
 567         fprintf(stderr, "Will try to allocate memory and convert to full matrix representation...\n");
 568
 569         size_t hessianSize = static_cast<size_t>(nrow)*static_cast<size_t>(ncol);
 570         /* Allowing  Hessians larger than INT_MAX probably only makes sense
 571          * with (OpenMP) parallel diagonalization routines, since with a single
 572          * thread it will takes months.
 573          */
 574         if (hessianSize > INT_MAX)
 575         {
 576             gmx_fatal(FARGS, "Hessian size is %d x %d, which is larger than the maximum allowed %d elements.",
 577                       nrow, ncol, INT_MAX);
 578         }
 579         snew(full_hessian, hessianSize);
 580         for (i = 0; i < nrow*ncol; i++)
 581         {
 582             full_hessian[i] = 0;
 583         }
 584
 585         for (i = 0; i < sparse_hessian->nrow; i++)
 586         {
 587             for (j = 0; j < sparse_hessian->ndata[i]; j++)
 588             {
 589                 k                      = sparse_hessian->data[i][j].col;
 590                 value                  = sparse_hessian->data[i][j].value;
 591                 full_hessian[i*ndim+k] = value;
 592                 full_hessian[k*ndim+i] = value;
 593             }
 594         }
 595         gmx_sparsematrix_destroy(sparse_hessian);
 596         sparse_hessian = nullptr;
 597         fprintf(stderr, "Converted sparse to full matrix storage.\n");
 598     }
 599
 600     snew(eigenvalues, nrow);
 601
 602     if (full_hessian != nullptr)
 603     {
 604         /* Using full matrix storage */
 605         eigenvectors = allocateEigenvectors(nrow, begin, end, false);
 606
 607         nma_full_hessian(full_hessian, nrow, bM, &top, atom_index, begin, end,
 608                          eigenvalues, eigenvectors);
 609     }
 610     else
 611     {
 612         assert(sparse_hessian);
 613         /* Sparse memory storage, allocate memory for eigenvectors */
 614         eigenvectors = allocateEigenvectors(nrow, begin, end, true);
 615
 616         nma_sparse_hessian(sparse_hessian, bM, &top, atom_index, end, eigenvalues, eigenvectors);
 617     }
 618
 619     /* check the output, first 6 eigenvalues should be reasonably small */
 620     gmx_bool bSuck = FALSE;
 621     for (i = begin-1; (i < 6); i++)
 622     {
 623         if (std::abs(eigenvalues[i]) > 1.0e-3)
 624         {
 625             bSuck = TRUE;
 626         }
 627     }
 628     if (bSuck)
 629     {
 630         fprintf(stderr, "\nOne of the lowest 6 eigenvalues has a non-zero value.\n");
 631         fprintf(stderr, "This could mean that the reference structure was not\n");
 632         fprintf(stderr, "properly energy minimized.\n");
 633     }
 634
 635     /* now write the output */
 636     fprintf (stderr, "Writing eigenvalues...\n");
 637     out = xvgropen(opt2fn("-ol", NFILE, fnm),
 638                    "Eigenvalues", "Eigenvalue index", "Eigenvalue [Gromacs units]",
 639                    oenv);
 640     if (output_env_get_print_xvgr_codes(oenv))
 641     {
 642         if (bM)
 643         {
 644             fprintf(out, "@ subtitle \"mass weighted\"\n");
 645         }
 646         else
 647         {
 648             fprintf(out, "@ subtitle \"not mass weighted\"\n");
 649         }
 650     }
 651
 652     for (i = 0; i <= (end-begin); i++)
 653     {
 654         fprintf (out, "%6d %15g\n", begin+i, eigenvalues[i]);
 655     }
 656     xvgrclose(out);
 657
 658
 659     if (opt2bSet("-qc", NFILE, fnm))
 660     {
 661         qc = xvgropen(opt2fn("-qc", NFILE, fnm), "Quantum Corrections", "Eigenvector index", "", oenv);
 662         xvgr_legend(qc, asize(qcleg), qcleg, oenv);
 663         qcvtot = qutot = 0;
 664     }
 665     else
 666     {
 667         qc = nullptr;
 668     }
 669     printf("Writing eigenfrequencies - negative eigenvalues will be set to zero.\n");
 670
 671     out = xvgropen(opt2fn("-of", NFILE, fnm),
 672                    "Eigenfrequencies", "Eigenvector index", "Wavenumber [cm\\S-1\\N]",
 673                    oenv);
 674     if (output_env_get_print_xvgr_codes(oenv))
 675     {
 676         if (bM)
 677         {
 678             fprintf(out, "@ subtitle \"mass weighted\"\n");
 679         }
 680         else
 681         {
 682             fprintf(out, "@ subtitle \"not mass weighted\"\n");
 683         }
 684     }
 685     /* Spectrum ? */
 686     spec = nullptr;
 687     if (opt2bSet("-os", NFILE, fnm) && (maxspec > 0))
 688     {
 689         snew(spectrum, maxspec);
 690         spec = xvgropen(opt2fn("-os", NFILE, fnm),
 691                         "Vibrational spectrum based on harmonic approximation",
 692                         "\\f{12}w\\f{4} (cm\\S-1\\N)",
 693                         "Intensity [Gromacs units]",
 694                         oenv);
 695         for (i = 0; (i < maxspec); i++)
 696         {
 697             spectrum[i] = 0;
 698         }
 699     }
 700
 701     /* Gromacs units are kJ/(mol*nm*nm*amu),
 702      * where amu is the atomic mass unit.
 703      *
 704      * For the eigenfrequencies we want to convert this to spectroscopic absorption
 705      * wavenumbers given in cm^(-1), which is the frequency divided by the speed of
 706      * light. Do this by first converting to omega^2 (units 1/s), take the square
 707      * root, and finally divide by the speed of light (nm/ps in gromacs).
 708      */
 709     factor_gmx_to_omega2       = 1.0E21/(AVOGADRO*AMU);
 710     factor_omega_to_wavenumber = 1.0E-5/(2.0*M_PI*SPEED_OF_LIGHT);
 711
 712     value = 0;
 713     for (i = begin; (i <= end); i++)
 714     {
 715         value = eigenvalues[i-begin];
 716         if (value < 0)
 717         {
 718             value = 0;
 719         }
 720         omega = std::sqrt(value*factor_gmx_to_omega2);
 721         nu    = 1e-12*omega/(2*M_PI);
 722         value = omega*factor_omega_to_wavenumber;
 723         fprintf (out, "%6d %15g\n", i, value);
 724         if (nullptr != spec)
 725         {
 726             wfac = eigenvalues[i-begin]/(width*std::sqrt(2*M_PI));
 727             for (j = 0; (j < maxspec); j++)
 728             {
 729                 spectrum[j] += wfac*std::exp(-gmx::square(j-value)/(2*gmx::square(width)));
 730             }
 731         }
 732         if (nullptr != qc)
 733         {
 734             qcv = cv_corr(nu, T);
 735             qu  = u_corr(nu, T);
 736             if (i > end-nharm)
 737             {
 738                 qcv += BOLTZ*KILO;
 739                 qu  += BOLTZ*T;
 740             }
 741             fprintf (qc, "%6d %15g %15g\n", i, qcv, qu);
 742             qcvtot += qcv;
 743             qutot  += qu;
 744         }
 745     }
 746     xvgrclose(out);
 747     if (value >= maxspec)
 748     {
 749         printf("WARNING: high frequencies encountered (%g cm^-1).\n", value);
 750         printf("Your calculations may be incorrect due to e.g. improper minimization of\n");
 751         printf("your starting structure or due to issues in your topology.\n");
 752     }
 753     if (nullptr != spec)
 754     {
 755         for (j = 0; (j < maxspec); j++)
 756         {
 757             fprintf(spec, "%10g  %10g\n", 1.0*j, spectrum[j]);
 758         }
 759         xvgrclose(spec);
 760     }
 761     if (nullptr != qc)
 762     {
 763         printf("Quantum corrections for harmonic degrees of freedom\n");
 764         printf("Use appropriate -first and -last options to get reliable results.\n");
 765         printf("There were %d constraints in the simulation\n", nharm);
 766         printf("Total correction to cV = %g J/mol K\n", qcvtot);
 767         printf("Total correction to  H = %g kJ/mol\n", qutot);
 768         xvgrclose(qc);
 769         please_cite(stdout, "Caleman2011b");
 770     }
 771     /* Writing eigenvectors. Note that if mass scaling was used, the eigenvectors
 772      * were scaled back from mass weighted cartesian to plain cartesian in the
 773      * nma_full_hessian() or nma_sparse_hessian() routines. Mass scaled vectors
 774      * will not be strictly orthogonal in plain cartesian scalar products.
 775      */
 776     const real *eigenvectorPtr;
 777     if (full_hessian != nullptr)
 778     {
 779         eigenvectorPtr = eigenvectors;
 780     }
 781     else
 782     {
 783         /* The sparse matrix diagonalization store all eigenvectors up to end */
 784         eigenvectorPtr = eigenvectors + (begin - 1)*atom_index.size();
 785     }
 786     write_eigenvectors(opt2fn("-v", NFILE, fnm), atom_index.size(), eigenvectorPtr, FALSE, begin, end,
 787                        eWXR_NO, nullptr, FALSE, top_x, bM, eigenvalues);
 788
 789     if (begin == 1)
 790     {
 791         analyzeThermochemistry(stdout, top, top_x, atom_index, eigenvalues,
 792                                T, P, sigma_r, scale_factor, linear_toler);
 793         please_cite(stdout, "Spoel2018a");
 794     }
 795     else
 796     {
 797         printf("Cannot compute entropy when -first = %d\n", begin);
 798     }
 799
 800     return 0;
 801 }