--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#include "gmxpre.h"
+
+#include "entropy.h"
+
+#include <cmath>
+#include <cstdio>
+
+#include "gromacs/math/units.h"
+#include "gromacs/utility/fatalerror.h"
+
+static real eigval_to_frequency(real 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)
+{
+ int nskip = bLinear ? 5 : 6;
+ double S = 0;
+ double hbar = PLANCK1/(2*M_PI);
+ for (int i = nskip; (i < n); i++)
+ {
+ if (eigval[i] > 0)
+ {
+ double omega = 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);
+ }
+ }
+ else
+ {
+ fprintf(stderr, "eigval[%d] = %g\n", i, eigval[i]);
+ }
+ }
+ return S*RGAS;
+}
+
+real calc_entropy_schlitter(int n,
+ real eigval[],
+ real temperature,
+ gmx_bool bLinear)
+{
+ int 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);
+ }
+ double deter = 0;
+ for (int i = nskip; (i < n); i++)
+ {
+ double dd = 1+kteh*eigval[i]*evcorr;
+ deter += std::log(dd);
+ }
+ return 0.5*RGAS*deter;
+}
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Code for computing entropy from eigenvalues
+ *
+ * \author David van der Spoel <david.vanderspoel@icm.uu.se>
+ */
+#ifndef GMXANA_ENTROPY_H
+#define GMXANA_ENTROPY_H
+
+#include "gromacs/math/vec.h"
+#include "gromacs/utility/basedefinitions.h"
+
+/*! \brief Compute entropy using Schlitter formula
+ *
+ * Computes entropy for a molecule / molecular system using the
+ * algorithm due to Schlitter (Chem. Phys. Lett. 215 (1993)
+ * 617-621).
+ * 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);
+
+/*! \brief Compute entropy using Quasi-Harmonic formula
+ *
+ * Computes entropy for a molecule / molecular system using the
+ * Quasi-harmonic algorithm (Macromolecules 1984, 17, 1370).
+ * 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).
+ * \return the entropy (J/mol K)
+ */
+real calc_entropy_quasi_harmonic(int n,
+ real eigval[],
+ real temperature,
+ gmx_bool linear);
+
+#endif
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,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.
#include "gromacs/gmxana/gmx_ana.h"
#include "gromacs/math/do_fit.h"
#include "gromacs/math/functions.h"
-#include "gromacs/math/units.h"
#include "gromacs/math/vec.h"
#include "gromacs/pbcutil/rmpbc.h"
#include "gromacs/topology/index.h"
#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.h"
-static void calc_entropy_qh(FILE *fp, int n, real eigval[], real temp, int nskip)
-{
- int i;
- double hwkT, dS, S = 0;
- double hbar, lambda;
-
- hbar = PLANCK1/(2*M_PI);
- for (i = 0; (i < n-nskip); i++)
- {
- if (eigval[i] > 0)
- {
- double w;
- lambda = eigval[i]*AMU;
- w = std::sqrt(BOLTZMANN*temp/lambda)/NANO;
- hwkT = (hbar*w)/(BOLTZMANN*temp);
- dS = (hwkT/std::expm1(hwkT) - std::log1p(-std::exp(-hwkT)));
- S += dS;
- if (debug)
- {
- fprintf(debug, "i = %5d w = %10g lam = %10g hwkT = %10g dS = %10g\n",
- i, w, lambda, hwkT, dS);
- }
- }
- else
- {
- fprintf(stderr, "eigval[%d] = %g\n", i, eigval[i]);
- }
- }
- fprintf(fp, "The Entropy due to the Quasi Harmonic approximation is %g J/mol K\n",
- S*RGAS);
-}
-
-static void calc_entropy_schlitter(FILE *fp, int n, int nskip,
- real *eigval, real temp)
-{
- double dd, deter;
- int i;
- double hbar, kt, kteh, S;
-
- hbar = PLANCK1/(2*M_PI);
- kt = BOLTZMANN*temp;
- kteh = kt*std::exp(2.0)/(hbar*hbar)*AMU*(NANO*NANO);
- if (debug)
- {
- fprintf(debug, "n = %d, nskip = %d kteh = %g\n", n, nskip, kteh);
- }
-
- deter = 0;
- for (i = 0; (i < n-nskip); i++)
- {
- dd = 1+kteh*eigval[i];
- deter += std::log(dd);
- }
- S = 0.5*RGAS*deter;
-
- fprintf(fp, "The Entropy due to the Schlitter formula is %g J/mol K\n", S);
-}
+#include "entropy.h"
const char *proj_unit;
{
gmx_fatal(FARGS, "Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting");
}
- calc_entropy_qh(stdout, neig1, eigval1, temp, nskip);
- calc_entropy_schlitter(stdout, neig1, nskip, eigval1, temp);
+ printf("The Entropy due to the Schlitter formula is %g J/mol K\n",
+ calc_entropy_schlitter(neig1, eigval1, temp, FALSE));
}
if (bVec2)
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,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.
#include "gromacs/utility/pleasecite.h"
#include "gromacs/utility/smalloc.h"
+#include "entropy.h"
+
static double cv_corr(double nu, double T)
{
double x = PLANCK*nu/(BOLTZ*T);
"output."
};
- static gmx_bool bM = TRUE, bCons = FALSE;
+ static gmx_bool bM = TRUE, bCons = FALSE, bLinear = FALSE;
static int begin = 1, end = 50, maxspec = 4000;
static real T = 298.15, width = 1;
t_pargs pa[] =
"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" },
+ "Last eigenvector to write away. -1 (default) is use all dimensions." },
{ "-maxspec", FALSE, etINT, {&maxspec},
"Highest frequency (1/cm) to consider in the spectrum" },
{ "-T", FALSE, etREAL, {&T},
{
begin = 1;
}
- if (end > ndim)
+ if (end == -1 || end > ndim)
{
end = ndim;
}
write_eigenvectors(opt2fn("-v", NFILE, fnm), atom_index.size(), eigenvectorPtr, FALSE, begin, end,
eWXR_NO, nullptr, FALSE, top_x, bM, eigenvalues);
+ 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));
+ }
+ else
+ {
+ printf("Cannot compute entropy when -first = %d\n", begin);
+ }
+
+
return 0;
}
#
# This file is part of the GROMACS molecular simulation package.
#
-# Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
+# Copyright (c) 2013,2014,2015,2016,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.
gmx_add_gtest_executable(
${exename}
+ entropy.cpp
gmx_traj.cpp
gmx_trjconv.cpp
)
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2018,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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Tests for entropy code
+ *
+ * \author David van der Spoel <david.vanderspoel@icm.uu.se>
+ */
+#include "gmxpre.h"
+
+#include "gromacs/gmxana/entropy.h"
+
+#include "testutils/refdata.h"
+#include "testutils/testasserts.h"
+
+namespace gmx
+{
+
+namespace
+{
+
+class Entropy : public ::testing::Test
+{
+ protected:
+ test::TestReferenceData refData_;
+ test::TestReferenceChecker checker_;
+ Entropy( ) : checker_(refData_.rootChecker())
+ {
+ }
+ public:
+ void runSchlitter(real temperature, gmx_bool bLinear)
+ {
+ std::vector<double> ev = {
+ 0.00197861, 0.000389439, 0.000316043,
+ 0.000150392, 0.000110254, 8.99659e-05,
+ 8.06572e-05, 5.14339e-05, 4.34268e-05,
+ 2.16063e-05, 1.65182e-05, 1.3965e-05,
+ 1.01937e-05, 5.83113e-06, 4.59067e-06,
+ 3.39577e-06, 2.14837e-06, 1.20468e-06,
+ 9.60817e-18, 1.48941e-19, 1.13939e-19,
+ 5.02332e-20, -6.90708e-21, -1.91354e-18
+ };
+ 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);
+ checker_.setDefaultTolerance(test::relativeToleranceAsFloatingPoint(1, 1e-7));
+ checker_.checkReal(S, "entropy");
+ }
+
+ void runQuasiHarmonic(real temperature, gmx_bool bLinear)
+ {
+ std::vector<double> ev = {
+ -31.403, -7.73169, -3.80315, -2.15659e-06, -1.70991e-07,
+ 236, 4609.83, 6718.07, 6966.27, 8587.85,
+ 10736.3, 13543.7, 17721.3, 22868, 35517.8,
+ 44118.1, 75827.9, 106277, 115132, 120782,
+ 445118, 451149, 481058, 484576
+ };
+ 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);
+ checker_.setDefaultTolerance(test::relativeToleranceAsFloatingPoint(1, 1e-7));
+ checker_.checkReal(S, "entropy");
+ }
+};
+
+TEST_F(Entropy, Schlitter_300_NoLinear)
+{
+ runSchlitter(300.0, FALSE);
+}
+
+TEST_F(Entropy, Schlitter_300_Linear)
+{
+ runSchlitter(300.0, TRUE);
+}
+
+TEST_F(Entropy, QuasiHarmonic_300_NoLinear)
+{
+ runQuasiHarmonic(300.0, FALSE);
+}
+
+TEST_F(Entropy, QuasiHarmonic_200_NoLinear)
+{
+ runQuasiHarmonic(200.0, FALSE);
+}
+
+TEST_F(Entropy, QuasiHarmonic_200_Linear)
+{
+ runQuasiHarmonic(200.0, TRUE);
+}
+
+} // namespace
+
+} // namespace gmx
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Real Name="entropy">21.341497503554653</Real>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Real Name="entropy">8.4752380507583425</Real>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Real Name="entropy">21.602209746228603</Real>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Real Name="entropy">5.8298021123436206</Real>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Real Name="entropy">4.3981206131323711</Real>
+</ReferenceData>
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff