src/gromacs/gmxana/thermochemistry.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2017,2018,2019, 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/(AVOGADRO*AMU);
  50     return std::sqrt(std::max(0.0, eigval)*factor_gmx_to_omega2);
  51 }
  52
  53 double calcZeroPointEnergy(gmx::ArrayRef<const real> eigval,
  54                            real                      scale_factor)
  55 {
  56     // Convert frequency (ps^-1) to energy (kJ/mol)
  57     double factor = PLANCK*PICO/(2.0*M_PI);
  58     double zpe    = 0;
  59     for (auto &r : eigval)
  60     {
  61         double omega = eigval_to_frequency(r);
  62         zpe         += 0.5*factor*scale_factor*omega;
  63     }
  64     return zpe;
  65 }
  66
  67 double calcVibrationalInternalEnergy(gmx::ArrayRef<const real> eigval,
  68                                      real                      temperature,
  69                                      gmx_bool                  linear,
  70                                      real                      scale_factor)
  71 {
  72     size_t nskip = linear ? 5 : 6;
  73     double Evib  = 0;
  74     double hbar  = PLANCK1/(2*M_PI);
  75     for (gmx::index i = nskip; i < eigval.ssize(); i++)
  76     {
  77         if (eigval[i] > 0)
  78         {
  79             double omega = scale_factor*eigval_to_frequency(eigval[i]);
  80             double hwkT  = (hbar*omega)/(BOLTZMANN*temperature);
  81             // Prevent overflow by checking for unreasonably large numbers.
  82             if (hwkT < 100)
  83             {
  84                 double dEvib = hwkT*(0.5 + 1.0/(std::expm1(hwkT)));
  85                 if (debug)
  86                 {
  87                     fprintf(debug, "i %d eigval %g omega %g hwkT %g dEvib %g\n",
  88                             static_cast<int>(i+1), static_cast<double>(eigval[i]),
  89                             omega, hwkT, dEvib);
  90                 }
  91                 Evib += dEvib;
  92             }
  93         }
  94     }
  95     return temperature*BOLTZ*Evib;
  96
  97 }
  98
  99 double calcVibrationalHeatCapacity(gmx::ArrayRef<const real> eigval,
 100                                    real                      temperature,
 101                                    gmx_bool                  linear,
 102                                    real                      scale_factor)
 103 {
 104     size_t nskip = linear ? 5 : 6;
 105     double cv    = 0;
 106     double hbar  = PLANCK1/(2*M_PI);
 107     for (gmx::index i = nskip; i < eigval.ssize(); i++)
 108     {
 109         if (eigval[i] > 0)
 110         {
 111             double omega = scale_factor*eigval_to_frequency(eigval[i]);
 112             double hwkT  = (hbar*omega)/(BOLTZMANN*temperature);
 113             // Prevent overflow by checking for unreasonably large numbers.
 114             if (hwkT < 100)
 115             {
 116                 double dcv   = std::exp(hwkT)*gmx::square(hwkT/std::expm1(hwkT));
 117                 if (debug)
 118                 {
 119                     fprintf(debug, "i %d eigval %g omega %g hwkT %g dcv %g\n",
 120                             static_cast<int>(i+1), static_cast<double>(eigval[i]), omega, hwkT, dcv);
 121                 }
 122                 cv += dcv;
 123             }
 124         }
 125     }
 126     return RGAS * cv;
 127 }
 128
 129 double calcTranslationalEntropy(real mass,
 130                                 real temperature,
 131                                 real pressure)
 132 {
 133     double kT = BOLTZ*temperature;
 134
 135     GMX_RELEASE_ASSERT(mass > 0, "Molecular mass should be larger than zero");
 136     GMX_RELEASE_ASSERT(pressure > 0, "Pressure should be larger than zero");
 137     GMX_RELEASE_ASSERT(temperature > 0, "Temperature should be larger than zero");
 138     // Convert bar to Pascal
 139     double P  = pressure*1e5;
 140     double qT = (std::pow(2*M_PI*mass*kT/gmx::square(PLANCK), 1.5) *
 141                  (kT/P) * (1e30/AVOGADRO));
 142     return RGAS*(std::log(qT) + 2.5);
 143 }
 144
 145 double calcRotationalEntropy(real       temperature,
 146                              int        natom,
 147                              gmx_bool   linear,
 148                              const rvec theta,
 149                              real       sigma_r)
 150 {
 151     GMX_RELEASE_ASSERT(sigma_r > 0, "Symmetry factor should be larger than zero");
 152     GMX_RELEASE_ASSERT(temperature > 0, "Temperature should be larger than zero");
 153
 154     double sR = 0;
 155     if (natom > 1)
 156     {
 157         if (linear)
 158         {
 159             GMX_RELEASE_ASSERT(theta[0] > 0, "Theta should be larger than zero");
 160             double qR = temperature/(sigma_r * theta[0]);
 161             sR        = RGAS * (std::log(qR) + 1);
 162         }
 163         else
 164         {
 165             double Q  = theta[XX]*theta[YY]*theta[ZZ];
 166             GMX_RELEASE_ASSERT(Q > 0, "Q should be larger than zero");
 167             double qR = std::sqrt(M_PI * std::pow(temperature, 3)/Q)/sigma_r;
 168             sR        = RGAS * (std::log(qR) + 1.5);
 169         }
 170     }
 171     return sR;
 172 }
 173
 174 double calcQuasiHarmonicEntropy(gmx::ArrayRef<const real> eigval,
 175                                 real                      temperature,
 176                                 gmx_bool                  bLinear,
 177                                 real                      scale_factor)
 178 {
 179     size_t nskip = bLinear ? 5 : 6;
 180     double S     = 0;
 181     double hbar  = PLANCK1/(2*M_PI);
 182     for (gmx::index i = nskip; (i < eigval.ssize()); i++)
 183     {
 184         if (eigval[i] > 0)
 185         {
 186             double omega  = scale_factor*eigval_to_frequency(eigval[i]);
 187             double hwkT   = (hbar*omega)/(BOLTZMANN*temperature);
 188             double dS     = (hwkT/std::expm1(hwkT) - std::log1p(-std::exp(-hwkT)));
 189             S            += dS;
 190             if (debug)
 191             {
 192                 fprintf(debug, "i = %5d eigval = %10g w = %10g hwkT = %10g dS = %10g\n",
 193                         static_cast<int>(i+1), static_cast<double>(eigval[i]),
 194                         omega, hwkT, dS);
 195             }
 196         }
 197         else if (debug)
 198         {
 199             fprintf(debug, "eigval[%d] = %g\n", static_cast<int>(i+1),
 200                     static_cast<double>(eigval[i]));
 201         }
 202     }
 203     return S*RGAS;
 204 }
 205
 206 double calcSchlitterEntropy(gmx::ArrayRef<const real> eigval,
 207                             real                      temperature,
 208                             gmx_bool                  bLinear)
 209 {
 210     size_t nskip  = bLinear ? 5 : 6;
 211     double hbar   = PLANCK1/(2*M_PI);             // J s
 212     double kt     = BOLTZMANN*temperature;        // J
 213     double kteh   = kt*std::exp(2.0)/(hbar*hbar); // 1/(J s^2) = 1/(kg m^2)
 214     double evcorr = NANO*NANO*AMU;
 215     if (debug)
 216     {
 217         fprintf(debug, "n = %td, kteh = %g evcorr = %g\n",
 218                 ssize(eigval), kteh, evcorr);
 219     }
 220     double deter = 0;
 221     for (gmx::index i = nskip; i < eigval.ssize(); i++)
 222     {
 223         double dd    = 1+kteh*eigval[i]*evcorr;
 224         deter       += std::log(dd);
 225     }
 226     return 0.5*RGAS*deter;
 227 }