src/gromacs/gmxana/thermochemistry.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2017,2018,2019,2020,2021, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include "gmxpre.h"
  36
  37 #include "thermochemistry.h"
  38
  39 #include <cmath>
  40 #include <cstdio>
  41
  42 #include "gromacs/math/units.h"
  43 #include "gromacs/utility/arrayref.h"
  44 #include "gromacs/utility/fatalerror.h"
  45 #include "gromacs/utility/gmxassert.h"
  46
  47 static double eigval_to_frequency(double eigval)
  48 {
  49     double factor_gmx_to_omega2 = 1.0E21 / (gmx::c_avogadro * gmx::c_amu);
  50     return std::sqrt(std::max(0.0, eigval) * factor_gmx_to_omega2);
  51 }
  52
  53 double calcZeroPointEnergy(gmx::ArrayRef<const real> eigval, real scale_factor)
  54 {
  55     // Convert frequency (ps^-1) to energy (kJ/mol)
  56     double factor = gmx::c_planck * gmx::c_pico / (2.0 * M_PI);
  57     double zpe    = 0;
  58     for (auto& r : eigval)
  59     {
  60         double omega = eigval_to_frequency(r);
  61         zpe += 0.5 * factor * scale_factor * omega;
  62     }
  63     return zpe;
  64 }
  65
  66 double calcVibrationalInternalEnergy(gmx::ArrayRef<const real> eigval, real temperature, gmx_bool linear, real scale_factor)
  67 {
  68     size_t nskip = linear ? 5 : 6;
  69     double Evib  = 0;
  70     double hbar  = gmx::c_planck1 / (2 * M_PI);
  71     for (gmx::index i = nskip; i < eigval.ssize(); i++)
  72     {
  73         if (eigval[i] > 0)
  74         {
  75             double omega = scale_factor * eigval_to_frequency(eigval[i]);
  76             double hwkT  = (hbar * omega) / (gmx::c_boltzmann * temperature);
  77             // Prevent overflow by checking for unreasonably large numbers.
  78             if (hwkT < 100)
  79             {
  80                 double dEvib = hwkT * (0.5 + 1.0 / (std::expm1(hwkT)));
  81                 if (debug)
  82                 {
  83                     fprintf(debug,
  84                             "i %d eigval %g omega %g hwkT %g dEvib %g\n",
  85                             static_cast<int>(i + 1),
  86                             static_cast<double>(eigval[i]),
  87                             omega,
  88                             hwkT,
  89                             dEvib);
  90                 }
  91                 Evib += dEvib;
  92             }
  93         }
  94     }
  95     return temperature * gmx::c_boltz * Evib;
  96 }
  97
  98 double calcVibrationalHeatCapacity(gmx::ArrayRef<const real> eigval, real temperature, gmx_bool linear, real scale_factor)
  99 {
 100     size_t nskip = linear ? 5 : 6;
 101     double cv    = 0;
 102     double hbar  = gmx::c_planck1 / (2 * M_PI);
 103     for (gmx::index i = nskip; i < eigval.ssize(); i++)
 104     {
 105         if (eigval[i] > 0)
 106         {
 107             double omega = scale_factor * eigval_to_frequency(eigval[i]);
 108             double hwkT  = (hbar * omega) / (gmx::c_boltzmann * temperature);
 109             // Prevent overflow by checking for unreasonably large numbers.
 110             if (hwkT < 100)
 111             {
 112                 double dcv = std::exp(hwkT) * gmx::square(hwkT / std::expm1(hwkT));
 113                 if (debug)
 114                 {
 115                     fprintf(debug,
 116                             "i %d eigval %g omega %g hwkT %g dcv %g\n",
 117                             static_cast<int>(i + 1),
 118                             static_cast<double>(eigval[i]),
 119                             omega,
 120                             hwkT,
 121                             dcv);
 122                 }
 123                 cv += dcv;
 124             }
 125         }
 126     }
 127     return gmx::c_universalGasConstant * cv;
 128 }
 129
 130 double calcTranslationalEntropy(real mass, real temperature, real pressure)
 131 {
 132     double kT = gmx::c_boltz * temperature;
 133
 134     GMX_RELEASE_ASSERT(mass > 0, "Molecular mass should be larger than zero");
 135     GMX_RELEASE_ASSERT(pressure > 0, "Pressure should be larger than zero");
 136     GMX_RELEASE_ASSERT(temperature > 0, "Temperature should be larger than zero");
 137     // Convert bar to Pascal
 138     double P  = pressure * 1e5;
 139     double qT = (std::pow(2 * M_PI * mass * kT / gmx::square(gmx::c_planck), 1.5) * (kT / P)
 140                  * (1e30 / gmx::c_avogadro));
 141     return gmx::c_universalGasConstant * (std::log(qT) + 2.5);
 142 }
 143
 144 double calcRotationalEntropy(real temperature, int natom, gmx_bool linear, const rvec theta, real sigma_r)
 145 {
 146     GMX_RELEASE_ASSERT(sigma_r > 0, "Symmetry factor should be larger than zero");
 147     GMX_RELEASE_ASSERT(temperature > 0, "Temperature should be larger than zero");
 148
 149     double sR = 0;
 150     if (natom > 1)
 151     {
 152         if (linear)
 153         {
 154             GMX_RELEASE_ASSERT(theta[0] > 0, "Theta should be larger than zero");
 155             double qR = temperature / (sigma_r * theta[0]);
 156             sR        = gmx::c_universalGasConstant * (std::log(qR) + 1);
 157         }
 158         else
 159         {
 160             double Q = theta[XX] * theta[YY] * theta[ZZ];
 161             GMX_RELEASE_ASSERT(Q > 0, "Q should be larger than zero");
 162             double qR = std::sqrt(M_PI * std::pow(temperature, 3) / Q) / sigma_r;
 163             sR        = gmx::c_universalGasConstant * (std::log(qR) + 1.5);
 164         }
 165     }
 166     return sR;
 167 }
 168
 169 double calcQuasiHarmonicEntropy(gmx::ArrayRef<const real> eigval, real temperature, gmx_bool bLinear, real scale_factor)
 170 {
 171     size_t nskip = bLinear ? 5 : 6;
 172     double S     = 0;
 173     double hbar  = gmx::c_planck1 / (2 * M_PI);
 174     for (gmx::index i = nskip; (i < eigval.ssize()); i++)
 175     {
 176         if (eigval[i] > 0)
 177         {
 178             double omega = scale_factor * eigval_to_frequency(eigval[i]);
 179             double hwkT  = (hbar * omega) / (gmx::c_boltzmann * temperature);
 180             double dS    = (hwkT / std::expm1(hwkT) - std::log1p(-std::exp(-hwkT)));
 181             S += dS;
 182             if (debug)
 183             {
 184                 fprintf(debug,
 185                         "i = %5d eigval = %10g w = %10g hwkT = %10g dS = %10g\n",
 186                         static_cast<int>(i + 1),
 187                         static_cast<double>(eigval[i]),
 188                         omega,
 189                         hwkT,
 190                         dS);
 191             }
 192         }
 193         else if (debug)
 194         {
 195             fprintf(debug, "eigval[%d] = %g\n", static_cast<int>(i + 1), static_cast<double>(eigval[i]));
 196         }
 197     }
 198     return S * gmx::c_universalGasConstant;
 199 }
 200
 201 double calcSchlitterEntropy(gmx::ArrayRef<const real> eigval, real temperature, gmx_bool bLinear)
 202 {
 203     size_t nskip  = bLinear ? 5 : 6;
 204     double hbar   = gmx::c_planck1 / (2 * M_PI);        // J s
 205     double kt     = gmx::c_boltzmann * temperature;     // J
 206     double kteh   = kt * std::exp(2.0) / (hbar * hbar); // 1/(J s^2) = 1/(kg m^2)
 207     double evcorr = gmx::c_nano * gmx::c_nano * gmx::c_amu;
 208     if (debug)
 209     {
 210         fprintf(debug, "n = %td, kteh = %g evcorr = %g\n", ssize(eigval), kteh, evcorr);
 211     }
 212     double deter = 0;
 213     for (gmx::index i = nskip; i < eigval.ssize(); i++)
 214     {
 215         double dd = 1 + kteh * eigval[i] * evcorr;
 216         deter += std::log(dd);
 217     }
 218     return 0.5 * gmx::c_universalGasConstant * deter;
 219 }