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