src/gromacs/gmxana/thermochemistry.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2017,2018, 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/fatalerror.h"
  44
  45 static real eigval_to_frequency(real eigval)
  46 {
  47     double factor_gmx_to_omega2       = 1.0E21/(AVOGADRO*AMU);
  48     return std::sqrt(eigval*factor_gmx_to_omega2);
  49 }
  50
  51 real calc_entropy_quasi_harmonic(int      n,
  52                                  real     eigval[],
  53                                  real     temperature,
  54                                  gmx_bool bLinear)
  55 {
  56     int    nskip = bLinear ? 5 : 6;
  57     double S     = 0;
  58     double hbar  = PLANCK1/(2*M_PI);
  59     for (int i = nskip; (i < n); i++)
  60     {
  61         if (eigval[i] > 0)
  62         {
  63             double omega  = eigval_to_frequency(eigval[i]);
  64             double hwkT   = (hbar*omega)/(BOLTZMANN*temperature);
  65             double dS     = (hwkT/std::expm1(hwkT) - std::log1p(-std::exp(-hwkT)));
  66             S            += dS;
  67             if (debug)
  68             {
  69                 fprintf(debug, "i = %5d eigval = %10g w = %10g hwkT = %10g dS = %10g\n",
  70                         i, eigval[i], omega, hwkT, dS);
  71             }
  72         }
  73         else
  74         {
  75             fprintf(stderr, "eigval[%d] = %g\n", i, eigval[i]);
  76         }
  77     }
  78     return S*RGAS;
  79 }
  80
  81 real calc_entropy_schlitter(int      n,
  82                             real     eigval[],
  83                             real     temperature,
  84                             gmx_bool bLinear)
  85 {
  86     int    nskip  = bLinear ? 5 : 6;
  87     double hbar   = PLANCK1/(2*M_PI);             // J s
  88     double kt     = BOLTZMANN*temperature;        // J
  89     double kteh   = kt*std::exp(2.0)/(hbar*hbar); // 1/(J s^2) = 1/(kg m^2)
  90     double evcorr = NANO*NANO*AMU;
  91     if (debug)
  92     {
  93         fprintf(debug, "n = %d, kteh = %g evcorr = %g\n", n, kteh, evcorr);
  94     }
  95     double deter = 0;
  96     for (int i = nskip; (i < n); i++)
  97     {
  98         double dd    = 1+kteh*eigval[i]*evcorr;
  99         deter       += std::log(dd);
 100     }
 101     return 0.5*RGAS*deter;
 102 }