/*************************************************************************************
* *
* ***** OpenThermo ***** *
* Calculation of thermodynamic functions from molecular data *
* Copyright 2007-2009 Konstantin Tokarev <annulen@users.sourceforge.net> *
* and others *
* *
*************************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program 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 General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program (see COPYING); if not, see *
* http://www.gnu.org/licenses *
* *
*************************************************************************************
* Module name : kernel.cpp *
* Author : Tokarev. K *
* Last modification : 2008/07/17 *
* Description : This module contains function lnZ, which determines *
* the partition function, and function Thermo, which calculates *
* values of thermodynamic function from lnZ and its derivatives *
* *
*************************************************************************************/
#include "thermo.h" // Build options
#include "kernel.hpp"
#include "rotation.hpp"
#ifdef NUMERIC
#include "mathmodule.hpp"
#endif
/*enum ROTATION_THEORY { NO, PITZER, CRAWFORD, S_DET};
enum INTERPOLATION_THEORY { NO, POLYNOM, B_SPLINE, FOURIER };
enum ANHARMONIC_THEORY { NO, X_DIAG, X_ALL };*/
// Fundamental constants
extern const double k=1.3806504e-23, //
h=6.62606896e-34, //
R=8.314472, //
c=2.99792458e10, // speed of light (cm/s!)
Na=6.02214179e23, // Avogadro constat
cal=4.184, //
Cels = 273.15,
atm = 101325,
HartreeJ = 4.35974394e-18, // Hartree unit in J
Hartree=Na*HartreeJ, // Hartree unit in kJ/mol 2625.49961748
pi = 3.14159265358979323846;
// Log's of fundamental constants
const double log_k=log(k),
log_h=log(h),
log_pi=log(pi);
/*namespace Control
{
// Thresholds
bool Atom::No_Spin=false;
bool Structure::Symmetry_Brutal=true;
bool Structure::Top_Brutal=true;
double Structure::Symmetry_Primary=0.5;
double Structure::Symmetry_Final=0.05;
double Structure::Symmetry_Maxit=500;
double Structure::Top_Type=0.2;
double QCutOff=1e-6;
double QPrintOut=2e-5;
}*/
// This variables are accessible from lnZ and main
namespace OpenThermo
{
Molecule * Mol = NULL;
int Anharmonic_Range = 0;
int Range = 0;
}
// Local declarations
// Local definitions
bool Initial (Molecule * M, const char * LogName)
{
using namespace OpenThermo;
Mol = M;
ofstream flog;
string pg;
flog.open(LogName, ios::app);
//flog << "Initialization...\n\n";
Mol->LogAtoms(flog);
pg = Mol->GetPointGroup();
flog << "Computing molecular mass...\n";
Mol->CalculateWeight();
flog << "Moving molecule to mass center...\n";
Mol->Move2MassCenter();
flog << "Computing inertial tensor...\n";
Mol->CalculateInertialTensor();
flog << "Computing nuclear spin states number...\n";
Mol->CalculateNuclearSpinStatesNumber();
if ((pg == "") || (pg == "auto"))
{
flog.close();
Mol->TestSymmetry(LogName);
flog.open(LogName, ios::app);
}
else
Mol->SetPointGroup(pg.c_str()); // It MUST remember :)
// Seriously, this trick is needed to get right rotor
// classification when symmetry is set explicitly
if ((Mol->GetSigma()==0)||(Mol->GetSigma()>24))
{
cerr << "## Invalid rotation factor! sigma = " << Mol->GetSigma() << "\n Exiting" << endl;
flog << "## Invalid rotation factor! sigma = " << Mol->GetSigma() << "\n Exiting" << endl;
return false;
}
Mol->LogAtoms(flog);
Mol->LogInfo(flog);
if (Mol->GetNRot()!=0)
{
flog << "\n\nInitialization of internal rotors...\n\n";
if (Mol->GetNRot() == 1)
flog << Mol->GetNRot() << " internal rotor found\n" << endl;
else
flog << Mol->GetNRot() << " internal rotors found\n" << endl;
for (int j=0; j<Mol->GetNRot(); j++)
{
flog << "\n-------------------------------------------------------------------------------\n";
flog << "Rotor #" << j+1 << endl;
Rotor * R = Mol->GetRotor(j);
R->LogAtoms(flog);
if (R->GetSigma()==0)
{
flog.close();
R->TestSymmetry(LogName);
flog.open(LogName, ios::app);
}
else if ((R->GetSigma()>60) || (R->GetSigma()<0))
{
cerr << "## Invalid rotation number for rotor " << j+1
<< "\n sigma = " << R->GetSigma() << "\n Exiting" << endl;
flog << "## Invalid rotation number for rotor " << j+1
<< "\n sigma = " << R->GetSigma() << "\n Exiting" << endl;
return false;
}
Mol->CalculateIeff(j);
R->LogInfo(flog);
}
flog << "\n-------------------------------------------------------------------------------\n";
Mol->CalculateSDet();
int n = 3+Mol->GetNRot();
flog << "Molecule rotation determinant:\t" << Mol->GetSDet()*pow(1e7,n)
<< " [g^" << n << "*cm^" << 2*n << "] \n\n" << endl;
}
flog << "Using frequencies list [cm^-1]:" << endl;
for (int i=0; i<Mol->GetNFreq(); i++)
flog << " " << i+1 << "\t" << Mol->GetFreq(i) << endl;
int k;
if(Mol->GetLinear())
k = 5;
else
k = 6;
if(Mol->GetNFreq() != 3*Mol->GetNAtoms()-k-Mol->GetNRot())
{
cerr << "# WARNING: Number of frequencies isn't equal to 3N-" << k+Mol->GetNRot()
<< ". Results may be senseless!" << endl;
flog << "# WARNING: Number of frequencies isn't equal to 3N-" << k+Mol->GetNRot()
<< ". Results may be senseless!" << endl;
}
flog.close();
return true;
}
double lnZ (double T, double V, double & zpve, ostream & flog)
{
using namespace OpenThermo;
const double Qthresh = 1e-8;
double lnQ, lnQtransl, lnQvibr, lnQrot=0, lnQel, lnQspin;
double ZPVE;
double lnQvibrPart, lnQrotGr;
double log_T=log(T), log_2=0.69314718055994530942;
double log_m = log(Mol->GetWeight())+log(1e-3/Na);
flog << "\t-------------------------------------------------------------------\n";
flog << "\t PARTITION FUNCTION ln(Q) at T = " << T << " K, V = " << V << " m^3" << endl;
flog << "\t-------------------------------------------------------------------\n" << endl;
// Partition function components
// Electronic states summ
// Qel = d1*exp(-E0/kT) + ...
double Qel = 0;
for (int i=0; i<Mol->GetNLevels(); i++)
{
int d;
double E;
Mol->GetLevel(i, d, E);
Qel += d*exp(-E/(k*T));
}
lnQel = log(Qel);
flog << "Electron transitions:\t" << lnQel << endl;
// Translations
// Qtransl = pow(2*pi*m*k*T, 1.5)/(h*h*h)*V;
lnQtransl = 1.5*(log_2+log_pi+log_m+log_k+log_T)-3*log_h+log(V);
flog << "Translations:\t\t" << lnQtransl << endl;
// Vibrations
lnQvibr = 0; //cout << Mol->GetFreq(0) << endl;
ZPVE = 0;
flog << "Vibrations:\t\t\t\tContribution" << endl;
for (int i=0; i<Mol->GetNFreq(); i++)
{
/*if(i < Anharmonic_Range)
{
double E, we, w0, Xe;
float w1, w2, w3;
int v = 0;
double Qvibr = 0;
we = Mol->GetFreq(i);
//w0 = Mol->GetFreq_Anh(i);
//Xe = (we-w0)/(2*we);
Xe = Mol->GetAnh(i,i);
do
{
E = h*we*c*(v-Xe*(v*v+v));
Qvibr += exp(-E/(k*T));
v++;
}while (E < h*w0*c/(4*Xe));
ZPVE += h*we*c*(0.5-Xe*0.25);
lnQvibrPart = log(Qvibr);
lnQvibr += lnQvibrPart;
// Overtones
w1 = 2*we*(1-3*Xe);
w2 = 3*we*(1-4*Xe);
w3 = 4*we*(1-5*Xe);
flog << "\tVibration " << i+1 << "\t" << Mol->GetFreq(i) << "\t\t" << lnQvibrPart
<< "\t\tObserved frequencies : " << w0 << " "
<< w1 << " " << w2 << " " << w3 << "; xe = " << Xe << "; v_max = " << v-1 << endl;
}
else
{ */
ZPVE += 0.5*h*Mol->GetFreq(i)*c;
// Qvibr = Qvibr/(1-exp(-h*Freq[i]*c)/(k*T));
lnQvibrPart = -log(1-exp((-h*Mol->GetFreq(i)*c)/(k*T)));
lnQvibr += lnQvibrPart;
if (lnQvibrPart > 2e-5)
flog << "\tVibration " << i+1 << "\t" << Mol->GetFreq(i) << "\t\t" << lnQvibrPart << endl;
// }
}
flog << "\nVibrational total:\t\t\t" << lnQvibr << endl;
flog << "Zero Point Vibrational Energy (ZPVE)\t" << ZPVE/HartreeJ
<< " Hartree\n \t" << ZPVE*Na/1000 << " kJ/mol\n"<< endl;
zpve = ZPVE*Na/1000; // for return
// Rotations
#ifdef S_DET
int sigma = Mol->GetSigma();
for (int j=0; j<Mol->GetNRot(); j++)
{
sigma *= (Mol->GetRotor(j))->GetSigma();
//RI *= RotationIntegral(R->GetRotType(), T, R->GetRotParams());
}
// Qrot = 8pi^2/sigma * (2pikT/h^2)^0.5(NRot+3) * sqrt(|s|) * (2pi)^NRot
lnQrot = 3*log_2-log(static_cast<double>(sigma)) + 0.5*(Mol->GetNRot()+3)*(log_2+log_pi+log_k+log_T-2*log_h) +
0.5*log(Mol->GetSDet()) + Mol->GetNRot()*(log_2+log_pi);
//lnQrot = 0.5*(log(Mol->GetSDet()))-log(static_cast<double> (sigma)))
// +Mol->GetNRot()*(1.5*(log_2+log_pi)+0.5*(log_k+log_T)-log_h);
#else /* S_DET */
// Calculate partition functions of molecule and independent rotors
double Qrot=0, RI=0, Ia=0, Ib=0;
#ifndef CLASSIC_ROTATION
double lnQrot_old=-3*Qthresh, tau1=0, tau2=0;
int J=0;
#endif
#ifndef CLASSIC_ROTATION
switch (Mol->GetTopType())
{
case SPHERICAL:
// Discreet sum
//cout << "for spherical\n";
Qrot = 0;
tau1 = h*h/(8*pi*pi*k*T*Mol->GetI(1));
J=0;
do
{
if (J!=0)
lnQrot_old = log(Qrot);
Qrot += (2*J+1)*(2*J+1)*exp(-J*(J+1)*tau1);
J++;
}while (fabs(log(Qrot)-lnQrot_old) > Qthresh);
lnQrot = log(Qrot)-log(static_cast<double> (Mol->GetSigma()));
break;
// Reminder: I(1)<=I(2)<=I(3)
case OBLONG: // Ia<Ib, I(2)=I(3)
//cout << "for sym1\n";
// Discreet sum
Ia = Mol->GetI(1), Ib = (Mol->GetI(2)+Mol->GetI(3))/2;
goto generic_symmetric;
case OBLATE: // Ia>Ib, I(1)=I(2)
//cout << "for sym2\n";
// Discreet sum
Ia = Mol->GetI(3), Ib = (Mol->GetI(1)+Mol->GetI(2))/2;
generic_symmetric:
//cout << "for sym\n";
// Discreet sum
Qrot = 0;
tau1 = h*h/(8*pi*pi*k*T*Ib);
tau2 = h*h/(8*pi*pi*k*T)*(1/Ia-1/Ib);
J=0;
do
{
if (J!=0)
lnQrot_old = log(Qrot);
for (int K=0; K<=J; K++)
if (K!=0)
Qrot += 2*(2*J+1)*exp(-J*(J+1)*tau1-K*K*tau2);
else
Qrot += (2*J+1)*exp(-J*(J+1)*tau1-K*K*tau2);
J++;
}while (fabs(log(Qrot)-lnQrot_old) > Qthresh);
lnQrot = log(Qrot)-log(static_cast<double> (Mol->GetSigma()));
break;
case ASYMMETRIC:
//cout << "for asym\n";
// Qrot = sqrt(Mol->GetDetI()*pi)*pow(8*pi*pi*k*T, 1.5)/(h*h*h*sigma[0]);
lnQrot = 0.5*(log(Mol->GetDetI())+log_pi)+1.5*(log(8.0)+2*log_pi+log_k+log_T)-3*log_h-log(static_cast<double> (Mol->GetSigma()));
break;
case LINEAR:
//cout << "for linear\n";
// Classic integral
// Qrot = Ilinear*8*pi*pi*k*T/(h*h*h*sigma);
// lnQrot = log(8.0)+2*log_pi+log(Mol->GetIlinear())+log_k+log_T-2*log_h;
// Discreet sum
Qrot = 0;
Ia = (Mol->GetI(2)+Mol->GetI(3))/2;
tau1 = h*h/(8*pi*pi*k*T*Ia);
J=0;
do
{
if (J!=0)
lnQrot_old = log(Qrot);
Qrot += (2*J+1)*exp(-J*(J+1)*tau1);
J++;
}while (fabs(log(Qrot)-lnQrot_old) > Qthresh);
lnQrot = log(Qrot)-log(static_cast<double> (Mol->GetSigma()));
break;
default:
flog << "## Illegal top type for molecule!" << endl;
cout << "## Illegal top type for molecule!" << endl;
}
#else /* CLASSIC_ROTATION */
if (Mol->GetLinear())
{
Ia = (Mol->GetI(2)+Mol->GetI(3))/2;
lnQrot = log(8.0)+2*log_pi+log(Ia)+log_k+log_T-2*log_h;
}
else
//cout << "for asym\n";
// Qrot = sqrt(Mol->GetDetI()*pi)*pow(8*pi*pi*k*T, 1.5)/(h*h*h*sigma[0]);
lnQrot = 0.5*(log(Mol->GetDetI())+log_pi)+1.5*(3*log_2+2*log_pi+log_k+log_T)-3*log_h-log(static_cast<double> (Mol->GetSigma()));
#endif /* CLASSIC_ROTATION */
flog << "Rotations:\n\tWhole system:\t" << lnQrot << " \tRotat. T = " << T/exp(lnQrot) << " K";
// Should be replaced with something like if (theory != classic)
#ifndef CLASSIC_ROTATION
if (Mol->GetTopType() != ASYMMETRIC)
flog << "\tJ_max = " << J-1;
#endif /* CLASSIC_ROTATION */
flog << endl;
// Internal rotations
if(Mol->GetNRot())
{
for (int j=0; j<Mol->GetNRot(); j++)
{
Rotor * R = Mol->GetRotor(j);
if (R->GetRotType() != FREE)
RI = RotationIntegral(T, R->GetBarrier());
else
RI = 2*pi;
#ifndef CLASSIC_ROTATION
lnQrotGr = 0; lnQrot_old=-3*Qthresh;
switch(R->GetTopType())
{
case SYMMETRIC:
//cout << "for sym rotor" << endl;
if (R->GetRotType() != FREE)
goto hindered; // if hindered rotation, evaluate as quasi-classic
tau1 = h/sqrt(2*pi*k*T*R->GetIeff()*RI);
Qrot = 0;
J = 0;
do
{
if (J!=0)
lnQrot_old = log(Qrot);
Qrot += (2*J+1)*exp(-J*(J+1)*tau1);
J++;
}while (fabs(log(Qrot)-lnQrot_old) > Qthresh);
lnQrotGr = log(Qrot)-log(static_cast<double> (R->GetSigma()));
break;
hindered:
case ASYMMETRIC:
//cout << "for asym rotor" << endl;
// QrotGr = sqrt(2*pi*Ieff*kT)*RotationIntegral/(h*sigma)
lnQrotGr = 0.5*(log_2+log_pi+log(R->GetIeff())+log_k+log_T)
+log(RI)-log_h-log(static_cast<double> (R->GetSigma()));
break;
default:
flog << "## Illegal top type for rotor!" << endl;
cout << "## Illegal top type for rotor!" << endl;
}
#else /* CLASSIC_ROTATION */
//cout << "for asym rotor" << endl;
// QrotGr = sqrt(2*pi*Ieff*kT)*RotationIntegral/(h*sigma)
lnQrotGr = 0.5*(log_2+log_pi+log(R->GetIeff())+log_k+log_T)
+log(RI)-log_h-log(static_cast<double> (R->GetSigma()));
#endif /* CLASSIC_ROTATION */
lnQrot += lnQrotGr;
flog << "\tRotor #" << j+1 << "\t" << lnQrotGr << " \tRotat. T = " << T/exp(lnQrotGr) << " K";
if (R->GetRotType() == HINDERED_HARMONIC)
{
flog << "; Barrier = " << (R->GetBarrier()).GetHeight()*Na << " J/mol";
// if one harmonic in potential!
flog << "; Vibr. frequency = "
<< sqrt((R->GetBarrier()).GetHeight()/(2*R->GetIeff()))*R->GetSigma();
//flog << "; Vibr. frequency = "
// << sqrt(R->GetRotParams()[0]/(2*R->GetIeff()))*R->GetSigma()/(2*pi*c) << " 1/cm";
flog << "; Rotat. integral = " << RI << " (" << 2*pi << " for free rotor)";
}
// Should be replaced with something like if (theory != classic)
#ifndef CLASSIC_ROTATION
if ((R->GetTopType() != ASYMMETRIC)&&(R->GetRotType()==FREE))
flog << "; J_max = " << J-1;
#endif /* CLASSIC_ROTATION */
if (T/exp(lnQrotGr) >= T)
flog << "\nWARNING! Low rotational temperature!";
flog << endl;
}
flog << "Rotational total:\t" << lnQrot << "\n" << endl;
}
#endif /* S_DET */
// Nuclear spin states
// Qspin = Mol->GetNuclearSpinStatesNumber();
lnQspin = log(static_cast<double> (Mol->GetNuclearSpinStatesNumber()));
flog << "Nuclear spin:\t" << lnQspin << "\n" << endl;
// Partition function for 1 molecule
// Q = Qel * Qtransl * Qvibr * Qrot * Qspin;
lnQ = lnQel + lnQtransl + lnQvibr + lnQrot + lnQspin;
flog << "\t*****************************************************************\n" << endl;
flog << "\t\tPARTITION FUNCTION ln(Q):\t\t" << lnQ << endl;
// Partition function for 1 mol
return Na*(lnQ - log(Na) + 1);
}
void Thermo (double T, double p, double & zpve, ostream & flog, ostream & fout)
{
// Calculates thermodinamic functions using Q and its analytical or numerical derivatives
double V, Cp, S, H, U, F, G, Cv; // H is H(T) - H(0)
double PF, PF_Der1, PF_Der2;
V = R*T/p; // Ideal gaz approximation
#ifndef NUMERIC
using MolData::Mol;
//Cp = R + k*T*(2*dlnZdT(T, p) + T*d2lnZdT2(T, p));
// Cp = R + Cv = R + kT(dlnZ/dT + T*d2lnZ/dT2)
//S = k*(lnZ(T, p) + T*dlnZdT(T, p)); // S = -dF/dT = k(lnZ + dlnZ/dT)
//H = R*T + k*T*T*dlnZdT(T, p); // H(T)-H(0) = pV + U(T)-U(0) = pV + (kT^2)*dlnZ/dT*/
S=2.5*R+R*log((pow((2*pi*(Mol->GetWeight())*(1e-3/Na)*k*T),1.5)*R*T)/(h*h*h*p*Na))
+ 1.5*R+R*log((8*pi*pi*sqrt(Mol->GetDetI())*pow(2*pi*k*T,1.5))/(Mol->GetSigma()*h*h*h));
H = 4*R*T;
Cp = 4*R;
for (int i=0; i<Mol->GetNFreq(); i++)
{
double Tv, eTv;
Tv=(h*Mol->GetFreq(i)*c)/(k*T);
eTv=exp(Tv);
S += R*(Tv/(eTv-1))-R*(log(1-1/eTv));
H += R*T*(Tv/(eTv-1));
Cp += R*((Tv*Tv*eTv)/((eTv-1)*(eTv-1)));
}
#else
ofstream FakeLog("");
double z = 0;
PF = lnZ(T, V, zpve, flog);
PF_Der1 = Diff(lnZ, T, V, z, FakeLog);
PF_Der2 = Diff2(lnZ, T, V, z, FakeLog);
flog << "\t\tPARTITION FUNCTION ln(Z):\t\t" << PF << endl;
flog << "\t\tPARTITION FUNCTION 1st DERIVATIVE:\t" << PF_Der1 << endl;
flog << "\t\tPARTITION FUNCTION 2nd DERIVATIVE:\t" << PF_Der2 << endl;
Cv = k*T*(2*PF_Der1 + T*PF_Der2); // Cp = R + Cv = R + kT(dlnZ/dT + T*d2lnZ/dT2)
S = k*(PF + T*PF_Der1); // S = -dF/dT = k(lnZ + dlnZ/dT)
U = k*T*T*PF_Der1; // U(T)-U(0) = (kT^2)*dlnZ/dT
H = U + p*V; // Ideal gas only!!!
Cp = Cv + R; // Ideal gas only!!!
F = U - T*S; // F(T) - F(0)
G = H - T*S; // G(T) - G(0)
flog << "\t\tHEAT CAPACITY Cv(" << "T): \t " << Cv << " J/(mol*K)" << endl;
flog << "\t\tHEAT CAPACITY Cp(" << "T): \t " << Cp << " J/(mol*K)" << endl;
flog << "\t\tENTROPY S(" << "T): \t " << S << " J/(mol*K)" << endl;
flog << "\t\tENERGY U(" << "T)-U(0): \t " << U/1000 << " kJ/mol" << endl;
flog << "\t\tENTHALPY H(" << "T)-H(0): \t " << H/1000 << " kJ/mol" << endl;
flog << "\t\tFREE ENERGY F(" << "T)-F(0): \t " << F/1000 << " kJ/mol" << endl;
flog << "\t\tFREE ENTHALPY G(" << "T)-G(0):\t " << G/1000 << " kJ/mol" << endl;
flog << "\n\t*****************************************************************\n" << endl;
flog << endl;
#endif
fout << resetiosflags(ios::showpoint) << T << setiosflags(ios::showpoint) << setprecision(5)
<< " " << Cv << " \t" << Cp << " \t" << setprecision(6)
<< S << " \t" << U/1000 << " \t\t" << H/1000 << " \t\t"
<< F/1000 << " \t" << G/1000 << endl;
#ifdef ECHO_ON
cout << T << " " << " " << Cv << " " << Cp << " " << S
<< " " << U/1000 << " " << H/1000 << " "
<< F/1000 << " " << G/1000 << endl;
#endif
}