The quasiharmonic entropy analysis computes vibrational entropy.
Here we add translational and rotational terms in order to compare
directly to experimental data. In addition the heat capacity and
vibrational internal energy are computed in the same manner.
The vibrational internal energy is printed as well.
A frequency scaling factor as is customary in quantum chemistry
code is available as well.
Change-Id: Ib9a76ecfc2f2cc27c90f3f9b58a924036a912e1f
gmx_fatal(FARGS, "Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting");
}
printf("The Entropy due to the Schlitter formula is %g J/mol K\n",
- calc_entropy_schlitter(neig1, eigval1, temp, FALSE));
+ calcSchlitterEntropy(gmx::arrayRefFromArray(eigval1, neig1),
+ temp, FALSE));
}
if (bVec2)
#include "gromacs/gmxana/eigio.h"
#include "gromacs/gmxana/gmx_ana.h"
#include "gromacs/gmxana/gstat.h"
+#include "gromacs/gmxana/princ.h"
#include "gromacs/linearalgebra/eigensolver.h"
#include "gromacs/linearalgebra/sparsematrix.h"
#include "gromacs/math/functions.h"
#include "gromacs/math/units.h"
#include "gromacs/math/vec.h"
+#include "gromacs/math/vecdump.h"
#include "gromacs/topology/ifunc.h"
#include "gromacs/topology/mtop_util.h"
#include "gromacs/topology/topology.h"
return eigenvectors;
}
+/*! \brief Compute heat capacity due to translational motion
+ */
+static double calcTranslationalHeatCapacity()
+{
+ return RGAS*1.5;
+}
+
+/*! \brief Compute internal energy due to translational motion
+ */
+static double calcTranslationalInternalEnergy(double T)
+{
+ return BOLTZ*T*1.5;
+}
+
+/*! \brief Compute heat capacity due to rotational motion
+ *
+ * \param[in] linear Should be TRUE if this is a linear molecule
+ * \param[in] T Temperature
+ * \return The heat capacity at constant volume
+ */
+static double calcRotationalInternalEnergy(gmx_bool linear, double T)
+{
+ if (linear)
+ {
+ return BOLTZ*T;
+ }
+ else
+ {
+ return BOLTZ*T*1.5;
+ }
+}
+
+/*! \brief Compute heat capacity due to rotational motion
+ *
+ * \param[in] linear Should be TRUE if this is a linear molecule
+ * \return The heat capacity at constant volume
+ */
+static double calcRotationalHeatCapacity(gmx_bool linear)
+{
+ if (linear)
+ {
+ return RGAS;
+ }
+ else
+ {
+ return RGAS*1.5;
+ }
+}
+
+static void analyzeThermochemistry(FILE *fp,
+ const t_topology &top,
+ rvec top_x[],
+ const std::vector<size_t> &atom_index,
+ real eigfreq[],
+ real T,
+ real P,
+ int sigma_r,
+ real scale_factor,
+ real linear_toler)
+{
+ std::vector<int> index;
+ for (auto &ai : atom_index)
+ {
+ index.push_back(static_cast<int>(ai));
+ }
+
+ rvec xcm;
+ double tmass = calc_xcm(top_x, index.size(),
+ index.data(), top.atoms.atom, xcm, FALSE);
+ double Strans = calcTranslationalEntropy(tmass, T, P);
+ std::vector<gmx::RVec> x_com;
+ x_com.resize(top.atoms.nr);
+ for (int i = 0; i < top.atoms.nr; i++)
+ {
+ copy_rvec(top_x[i], x_com[i]);
+ }
+ (void) sub_xcm(as_rvec_array(x_com.data()), index.size(), index.data(),
+ top.atoms.atom, xcm, FALSE);
+
+ rvec inertia;
+ matrix trans;
+ principal_comp(index.size(), index.data(), top.atoms.atom,
+ as_rvec_array(x_com.data()), trans, inertia);
+ bool linear = (inertia[XX]/inertia[YY] < linear_toler &&
+ inertia[XX]/inertia[ZZ] < linear_toler);
+ // (kJ/mol ps)^2/(Dalton nm^2 kJ/mol K) =
+ // KILO kg m^2 ps^2/(s^2 mol g/mol nm^2 K) =
+ // KILO^2 10^18 / 10^24 K = 1/K
+ double rot_const = gmx::square(PLANCK)/(8*gmx::square(M_PI)*BOLTZ);
+ // Rotational temperature (1/K)
+ rvec theta = { 0, 0, 0 };
+ if (linear)
+ {
+ // For linear molecules the first element of the inertia
+ // vector is zero.
+ theta[0] = rot_const/inertia[1];
+ }
+ else
+ {
+ for (int m = 0; m < DIM; m++)
+ {
+ theta[m] = rot_const/inertia[m];
+ }
+ }
+ if (debug)
+ {
+ pr_rvec(debug, 0, "inertia", inertia, DIM, TRUE);
+ pr_rvec(debug, 0, "theta", theta, DIM, TRUE);
+ pr_rvecs(debug, 0, "trans", trans, DIM);
+ fprintf(debug, "linear molecule = %s\n", linear ? "true" : "no");
+ }
+ size_t nFreq = index.size()*DIM;
+ auto eFreq = gmx::arrayRefFromArray(eigfreq, nFreq);
+ double Svib = calcQuasiHarmonicEntropy(eFreq, T, linear, scale_factor);
+
+ double Srot = calcRotationalEntropy(T, top.atoms.nr,
+ linear, theta, sigma_r);
+ fprintf(fp, "Translational entropy %g J/mol K\n", Strans);
+ fprintf(fp, "Rotational entropy %g J/mol K\n", Srot);
+ fprintf(fp, "Vibrational entropy %g J/mol K\n", Svib);
+ fprintf(fp, "Total entropy %g J/mol K\n", Svib+Strans+Srot);
+ double HeatCapacity = (calcTranslationalHeatCapacity() +
+ calcRotationalHeatCapacity(linear) +
+ calcVibrationalHeatCapacity(eFreq, T, linear, scale_factor));
+ fprintf(fp, "Heat capacity %g J/mol K\n", HeatCapacity);
+ double Evib = (calcTranslationalInternalEnergy(T) +
+ calcRotationalInternalEnergy(linear, T) +
+ calcVibrationalInternalEnergy(eFreq, T, linear, scale_factor));
+ fprintf(fp, "Internal energy %g kJ/mol\n", Evib);
+}
+
int gmx_nmeig(int argc, char *argv[])
{
"you need to analyze the [TT]quant_corr.xvg[tt] file yourself.[PAR]",
"To make things more flexible, the program can also take virtual sites into account",
"when computing quantum corrections. When selecting [TT]-constr[tt] and",
- "[TT]-qc[tt], the [TT]-begin[tt] and [TT]-end[tt] options will be set automatically as well.",
- "Again, if you think you know it better, please check the [TT]eigenfreq.xvg[tt]",
- "output."
+ "[TT]-qc[tt], the [TT]-begin[tt] and [TT]-end[tt] options will be set automatically as well.[PAR]",
+ "Based on a harmonic analysis of the normal mode frequencies,",
+ "thermochemical properties S0 (Standard Entropy),",
+ "Cv (Heat capacity at constant volume) and the internal energy can be",
+ "computed, much in the same manner as popular quantum chemistry",
+ "programs."
};
- static gmx_bool bM = TRUE, bCons = FALSE, bLinear = FALSE;
- static int begin = 1, end = 50, maxspec = 4000;
- static real T = 298.15, width = 1;
+ static gmx_bool bM = TRUE, bCons = FALSE;
+ static int begin = 1, end = 50, maxspec = 4000, sigma_r = 1;
+ static real T = 298.15, width = 1, P = 1, scale_factor = 1;
+ static real linear_toler = 1e-5;
+
t_pargs pa[] =
{
{ "-m", FALSE, etBOOL, {&bM},
"Divide elements of Hessian by product of sqrt(mass) of involved "
"atoms prior to diagonalization. This should be used for 'Normal Modes' "
"analysis" },
- { "-linear", FALSE, etBOOL, {&bLinear},
- "This should be set in order to get correct entropies for linear molecules" },
{ "-first", FALSE, etINT, {&begin},
"First eigenvector to write away" },
{ "-last", FALSE, etINT, {&end},
- "Last eigenvector to write away. -1 (default) is use all dimensions." },
+ "Last eigenvector to write away. -1 is use all dimensions." },
{ "-maxspec", FALSE, etINT, {&maxspec},
"Highest frequency (1/cm) to consider in the spectrum" },
{ "-T", FALSE, etREAL, {&T},
- "Temperature for computing quantum heat capacity and enthalpy when using normal mode calculations to correct classical simulations" },
+ "Temperature for computing entropy, quantum heat capacity and enthalpy when using normal mode calculations to correct classical simulations" },
+ { "-P", FALSE, etREAL, {&P},
+ "Pressure (bar) when computing entropy" },
+ { "-sigma", FALSE, etINT, {&sigma_r},
+ "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." },
+ { "-scale", FALSE, etREAL, {&scale_factor},
+ "Factor to scale frequencies before computing thermochemistry values" },
+ { "-linear_toler", FALSE, etREAL, {&linear_toler},
+ "Tolerance for determining whether a compound is linear as determined from the ration of the moments inertion Ix/Iy and Ix/Iz." },
{ "-constr", FALSE, etBOOL, {&bCons},
- "If constraints were used in the simulation but not in the normal mode analysis (this is the recommended way of doing it) you will need to set this for computing the quantum corrections." },
+ "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." },
{ "-width", FALSE, etREAL, {&width},
"Width (sigma) of the gaussian peaks (1/cm) when generating a spectrum" }
};
if (begin == 1)
{
- printf("The Entropy due to the Quasi Harmonic approximation is %g J/mol K\n",
- calc_entropy_quasi_harmonic(DIM*atom_index.size(),
- eigenvalues, T, bLinear));
+ analyzeThermochemistry(stdout, top, top_x, atom_index, eigenvalues,
+ T, P, sigma_r, scale_factor, linear_toler);
}
else
{
printf("Cannot compute entropy when -first = %d\n", begin);
}
-
return 0;
}
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2018,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2017,2018, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
std::vector<real> eigenvalue;
std::copy(ev.begin(), ev.end(), std::back_inserter(eigenvalue));
- real S = calc_entropy_schlitter(eigenvalue.size(), eigenvalue.data(),
- temperature, bLinear);
+ real S = calcSchlitterEntropy(eigenvalue, temperature, bLinear);
checker_.setDefaultTolerance(test::relativeToleranceAsFloatingPoint(1, 1e-7));
checker_.checkReal(S, "entropy");
}
std::vector<real> eigenvalue;
std::copy(ev.begin(), ev.end(), std::back_inserter(eigenvalue));
- real S = calc_entropy_quasi_harmonic(eigenvalue.size(), eigenvalue.data(),
- temperature, bLinear);
+ real S = calcQuasiHarmonicEntropy(eigenvalue, temperature, bLinear, 1.0);
+
checker_.setDefaultTolerance(test::relativeToleranceAsFloatingPoint(1, 1e-7));
checker_.checkReal(S, "entropy");
}
#include <cstdio>
#include "gromacs/math/units.h"
+#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
-static real eigval_to_frequency(real eigval)
+static double eigval_to_frequency(double eigval)
{
double factor_gmx_to_omega2 = 1.0E21/(AVOGADRO*AMU);
return std::sqrt(eigval*factor_gmx_to_omega2);
}
-real calc_entropy_quasi_harmonic(int n,
- real eigval[],
- real temperature,
- gmx_bool bLinear)
+double calcVibrationalInternalEnergy(gmx::ArrayRef<const real> eigval,
+ real temperature,
+ gmx_bool linear,
+ real scale_factor)
{
- int nskip = bLinear ? 5 : 6;
+ size_t nskip = linear ? 5 : 6;
+ double Evib = 0;
+ double hbar = PLANCK1/(2*M_PI);
+ for (size_t i = nskip; i < eigval.size(); i++)
+ {
+ if (eigval[i] > 0)
+ {
+ double omega = scale_factor*eigval_to_frequency(eigval[i]);
+ double hwkT = (hbar*omega)/(BOLTZMANN*temperature);
+ // Prevent overflow by checking for unreasonably large numbers.
+ if (hwkT < 100)
+ {
+ double dEvib = hwkT*(0.5 + 1.0/(std::expm1(hwkT)));
+ if (debug)
+ {
+ fprintf(debug, "i %d eigval %g omega %g hwkT %g dEvib %g\n",
+ static_cast<int>(i+1), static_cast<double>(eigval[i]),
+ omega, hwkT, dEvib);
+ }
+ Evib += dEvib;
+ }
+ }
+ }
+ return temperature*BOLTZ*Evib;
+
+}
+
+double calcVibrationalHeatCapacity(gmx::ArrayRef<const real> eigval,
+ real temperature,
+ gmx_bool linear,
+ real scale_factor)
+{
+ size_t nskip = linear ? 5 : 6;
+ double cv = 0;
+ double hbar = PLANCK1/(2*M_PI);
+ for (size_t i = nskip; i < eigval.size(); i++)
+ {
+ if (eigval[i] > 0)
+ {
+ double omega = scale_factor*eigval_to_frequency(eigval[i]);
+ double hwkT = (hbar*omega)/(BOLTZMANN*temperature);
+ // Prevent overflow by checking for unreasonably large numbers.
+ if (hwkT < 100)
+ {
+ double dcv = std::exp(hwkT)*gmx::square(hwkT/std::expm1(hwkT));
+ if (debug)
+ {
+ fprintf(debug, "i %d eigval %g omega %g hwkT %g dcv %g\n",
+ static_cast<int>(i+1), static_cast<double>(eigval[i]), omega, hwkT, dcv);
+ }
+ cv += dcv;
+ }
+ }
+ }
+ return RGAS * cv;
+}
+
+double calcTranslationalEntropy(real mass,
+ real temperature,
+ real pressure)
+{
+ double kT = BOLTZ*temperature;
+
+ GMX_RELEASE_ASSERT(mass > 0, "Molecular mass should be larger than zero");
+ GMX_RELEASE_ASSERT(pressure > 0, "Pressure should be larger than zero");
+ GMX_RELEASE_ASSERT(temperature > 0, "Temperature should be larger than zero");
+ // Convert bar to Pascal
+ double P = pressure*1e5;
+ double qT = (std::pow(2*M_PI*mass*kT/gmx::square(PLANCK), 1.5) *
+ (kT/P) * (1e30/AVOGADRO));
+ return RGAS*(std::log(qT) + 2.5);
+}
+
+double calcRotationalEntropy(real temperature,
+ int natom,
+ gmx_bool linear,
+ rvec theta,
+ real sigma_r)
+{
+ GMX_RELEASE_ASSERT(sigma_r > 0, "Symmetry factor should be larger than zero");
+ GMX_RELEASE_ASSERT(temperature > 0, "Temperature should be larger than zero");
+
+ double sR = 0;
+ if (natom > 1)
+ {
+ if (linear)
+ {
+ GMX_RELEASE_ASSERT(theta[0] > 0, "Theta should be larger than zero");
+ double qR = temperature/(sigma_r * theta[0]);
+ sR = RGAS * (std::log(qR) + 1);
+ }
+ else
+ {
+ double Q = theta[XX]*theta[YY]*theta[ZZ];
+ GMX_RELEASE_ASSERT(Q > 0, "Q should be larger than zero");
+ double qR = std::sqrt(M_PI * std::pow(temperature, 3)/Q)/sigma_r;
+ sR = RGAS * (std::log(qR) + 1.5);
+ }
+ }
+ return sR;
+}
+
+double calcQuasiHarmonicEntropy(gmx::ArrayRef<const real> eigval,
+ real temperature,
+ gmx_bool bLinear,
+ real scale_factor)
+{
+ size_t nskip = bLinear ? 5 : 6;
double S = 0;
double hbar = PLANCK1/(2*M_PI);
- for (int i = nskip; (i < n); i++)
+ for (size_t i = nskip; (i < eigval.size()); i++)
{
if (eigval[i] > 0)
{
- double omega = eigval_to_frequency(eigval[i]);
+ double omega = scale_factor*eigval_to_frequency(eigval[i]);
double hwkT = (hbar*omega)/(BOLTZMANN*temperature);
double dS = (hwkT/std::expm1(hwkT) - std::log1p(-std::exp(-hwkT)));
S += dS;
if (debug)
{
fprintf(debug, "i = %5d eigval = %10g w = %10g hwkT = %10g dS = %10g\n",
- i, eigval[i], omega, hwkT, dS);
+ static_cast<int>(i+1), static_cast<double>(eigval[i]),
+ omega, hwkT, dS);
}
}
- else
+ else if (debug)
{
- fprintf(stderr, "eigval[%d] = %g\n", i, eigval[i]);
+ fprintf(debug, "eigval[%d] = %g\n", static_cast<int>(i+1),
+ static_cast<double>(eigval[i]));
}
}
return S*RGAS;
}
-real calc_entropy_schlitter(int n,
- real eigval[],
- real temperature,
- gmx_bool bLinear)
+double calcSchlitterEntropy(gmx::ArrayRef<const real> eigval,
+ real temperature,
+ gmx_bool bLinear)
{
- int nskip = bLinear ? 5 : 6;
+ size_t nskip = bLinear ? 5 : 6;
double hbar = PLANCK1/(2*M_PI); // J s
double kt = BOLTZMANN*temperature; // J
double kteh = kt*std::exp(2.0)/(hbar*hbar); // 1/(J s^2) = 1/(kg m^2)
double evcorr = NANO*NANO*AMU;
if (debug)
{
- fprintf(debug, "n = %d, kteh = %g evcorr = %g\n", n, kteh, evcorr);
+ fprintf(debug, "n = %d, kteh = %g evcorr = %g\n",
+ static_cast<int>(eigval.size()), kteh, evcorr);
}
double deter = 0;
- for (int i = nskip; (i < n); i++)
+ for (size_t i = nskip; i < eigval.size(); i++)
{
double dd = 1+kteh*eigval[i]*evcorr;
deter += std::log(dd);
*/
/*! \internal \file
* \brief
- * Code for computing entropy from eigenvalues
+ * Code for computing entropy and heat capacity from eigenvalues
*
* \author David van der Spoel <david.vanderspoel@icm.uu.se>
*/
-#ifndef GMXANA_ENTROPY_H
-#define GMXANA_ENTROPY_H
+#ifndef GMXANA_THERMOCHEMISTRY_H
+#define GMXANA_THERMOCHEMISTRY_H
+#include "gromacs/math/units.h"
#include "gromacs/math/vec.h"
+#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/basedefinitions.h"
+/*! \brief Compute heat capacity due to vibrational motion
+ *
+ * \param[in] eigval The eigenvalues
+ * \param[in] temperature Temperature (K)
+ * \param[in] linear TRUE if this is a linear molecule
+ * \param[in] scale_factor Factor to scale frequencies by before computing cv
+ * \return The heat capacity at constant volume
+ */
+double calcVibrationalHeatCapacity(gmx::ArrayRef<const real> eigval,
+ real temperature,
+ gmx_bool linear,
+ real scale_factor);
+
+/*! \brief Compute entropy due to translational motion
+ *
+ * Following the equations in J. W. Ochterski,
+ * Thermochemistry in Gaussian, Gaussian, Inc., 2000
+ * Pitssburg PA
+ *
+ * \param[in] mass Molecular mass (Dalton)
+ * \param[in] temperature Temperature (K)
+ * \param[in] pressure Pressure (bar) at which to compute
+ * \returns The translational entropy
+ */
+double calcTranslationalEntropy(real mass,
+ real temperature,
+ real pressure);
+
+/*! \brief Compute entropy due to rotational motion
+ *
+ * Following the equations in J. W. Ochterski,
+ * Thermochemistry in Gaussian, Gaussian, Inc., 2000
+ * Pitssburg PA
+ *
+ * \param[in] temperature Temperature (K)
+ * \param[in] natom Number of atoms
+ * \param[in] linear TRUE if this is a linear molecule
+ * \param[in] theta The principal moments of inertia (unit of Energy)
+ * \param[in] sigma_r Symmetry factor, should be >= 1
+ * \returns The rotational entropy
+ */
+double calcRotationalEntropy(real temperature,
+ int natom,
+ gmx_bool linear,
+ rvec theta,
+ real sigma_r);
+
+/*! \brief Compute internal energy due to vibrational motion
+ *
+ * \param[in] eigval The eigenvalues
+ * \param[in] temperature Temperature (K)
+ * \param[in] linear TRUE if this is a linear molecule
+ * \param[in] scale_factor Factor to scale frequencies by before computing E
+ * \return The internal energy
+ */
+double calcVibrationalInternalEnergy(gmx::ArrayRef<const real> eigval,
+ real temperature,
+ gmx_bool linear,
+ real scale_factor);
+
/*! \brief Compute entropy using Schlitter formula
*
* Computes entropy for a molecule / molecular system using the
* The input should be eigenvalues from a covariance analysis,
* the units of the eigenvalues are those of energy.
*
- * \param[in] n Number of entries in the eigval array
* \param[in] eigval The eigenvalues
* \param[in] temperature Temperature (K)
* \param[in] linear True if this is a linear molecule (typically a diatomic molecule).
* \return the entropy (J/mol K)
*/
-real calc_entropy_schlitter(int n,
- real eigval[],
- real temperature,
- gmx_bool linear);
+double calcSchlitterEntropy(gmx::ArrayRef<const real> eigval,
+ real temperature,
+ gmx_bool linear);
/*! \brief Compute entropy using Quasi-Harmonic formula
*
* The input should be eigenvalues from a normal mode analysis.
* In both cases the units of the eigenvalues are those of energy.
*
- * \param[in] n Number of entries in the eigval array
* \param[in] eigval The eigenvalues
* \param[in] temperature Temperature (K)
* \param[in] linear True if this is a linear molecule (typically a diatomic molecule).
+ * \param[in] scale_factor Factor to scale frequencies by before computing S0
* \return the entropy (J/mol K)
*/
-real calc_entropy_quasi_harmonic(int n,
- real eigval[],
- real temperature,
- gmx_bool linear);
+double calcQuasiHarmonicEntropy(gmx::ArrayRef<const real> eigval,
+ real temperature,
+ gmx_bool linear,
+ real scale_factor);
#endif
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff