# -*- coding: utf-8 -*-
'''
Calculate thermodynamic data (S298, H298, and Cp(T)
from ab initio DFT data (energies and frequencies)
providng input thermodynamics files for
KMC (Zacros) and MKM (Chemkin and Matlab)
Gerhard R Wittreich, P.E.
Created on Fri Mar 31 2017
author wittregr
Adopted from Matlab code written and modified by mVassili Vorotnikov and Geun Ho Gu
This program contains the class objects used to read energy, vibration and
molecular configuration data and determine the standard entropy and enthalpy
and heat capacities at various temperatures.
'''
import numpy as _np
from numpy import pi
import os
import ase.io as _ase
import re
import scipy.interpolate as _sp
import datetime
from utils import constant as c
class Particle(object):
Cp_Range = _np.linspace(100, 1500, 15)
VibScalingFactor = 1 # Vibrational Scaling Factor
def __init__(self, data, dict, Base_path, Tstp=298.15):
'''
Fill object with species name and associated thermodynamic data
'''
self.name = str(data[dict['name']]) # Species name
if 'numh' in dict:
if data[dict['numc']] != '':
self.carbon = int(data[dict['numc']]) # No. of C atoms
else:
self.carbon = int(0)
if data[dict['numh']] != '':
self.hydrogen = int(data[dict['numh']]) # No. of H atoms
else:
self.hydrogen = int(0)
if data[dict['numo']] != '':
self.oxygen = int(data[dict['numo']]) # No. of O atoms
else:
self.oxygen = int(0)
if data[dict['numn']] != '':
self.nitrogen = int(data[dict['numn']]) # No. of N atoms
else:
self.nitrogen = int(0)
if 'mw' in dict:
self.MW = float(data[dict['mw']])/c.NA/1000
else:
self.MW = (self.carbon * c.MW_carbon +
self.hydrogen * c.MW_hydorgen +
self.oxygen * c.MW_oxygen +
self.nitrogen * c.MW_nitrogen)/c.NA/1000
if not hasattr(self, 'Inertia'):
self.totengpath = os.path.join(*re.split(r'\\|/',
str(data[dict['totengpath']]).
strip('.').strip('\\')))
self.etotal = float(data[dict['etotal']]) # Total energy
if data[dict['edisp']] == '':
self.edisp = float(0)
else:
self.edisp = float(data[dict['edisp']]) # Dispersion energy
self.numvibfreq = int(data[dict['numvibfreq']]) # No. vib frequencies
if not hasattr(self, 'phase'):
self.phase = None # Phase (G=gas, S=surface)
self.vibfreq = [] # Vibration frequencies
for x in range(0, self.numvibfreq):
self.vibfreq.append(Particle.VibScalingFactor *
float(data[dict['vibfreq'] + x]))
self.vibfreq = _np.array(self.vibfreq)
self.Base_path = Base_path
self.Tstp = Tstp
self.ThermoProperties()
def ThermoProperties(self):
'''
Calculate all thermodynamic properties from input data
'''
self.Determine_Phase()
'''
Get rotational data from VASP CONTCAR for gas species using
Atomic Simulation Environment (ASE) libraries for python
'''
if self.phase == 'G':
if hasattr(self, 'Inertia'):
self.I3 = self.Inertia
else:
filepath = os.path.join(self.Base_path,
self.totengpath,
'CONTCAR')
VASP = _ase.read(filepath)
self.I3 = VASP.get_moments_of_inertia()*c.A2_to_m2*c.amu_to_kg
self.MW = sum(VASP.get_masses())/c.NA/1000.
self.T_I = c.h1**2/(8*_np.pi**2*c.kb1)
'''
Calulcate common frequency data for vibrational components
'''
self.nu = self.vibfreq * 100 * c.c2
self.theta = c.h1 * self.nu / c.kb1
'''
Call Entropy method to calculate standard state entropy
'''
self.Calc_Entropy()
'''
Call Heat Capacity method to calculate heat capacities at the
temperature range specified by Cp_Range
'''
self.Calc_HeatCapacities()
'''
Call Enthalpy method to calculate standard state enthalpy
'''
self.Calc_Enthalpy()
def Determine_Phase(self):
'''
Determine species phase if one is not provided
'''
if self.phase is not None:
pass
elif hasattr(self, 'surface'):
self.phase = 'S'
elif self.islinear is not None:
self.phase = 'G'
elif hasattr(self, 'sigma'):
self.phase = 'G'
else:
'''
This should proabbaly result in an error condition vs supplying
'S' as a default value
'''
self.phase = 'S'
def Calc_Entropy(self):
'''
Calculate the vibrational, rotational and translational components and
total entropy of a species at standard conditions
'''
'''
Calculate vibrational component of entropy for gas and surface species
'''
T = self.Tstp
self.S_Tstp_vib = c.R1 * sum((self.theta/T) /
(_np.exp(self.theta/T)-1) -
_np.log(1 - _np.exp(-self.theta/T)))
self.q_vib = _np.product(_np.divide(1, (1 - _np.exp(-self.theta/T))))
'''
Calculate rotational and translational components of entropy
'''
if self.phase == 'G':
'''
Gas phase calculation
'''
if self.islinear == 0:
'''
Non-linear species
'''
I = _np.product(self.I3)
self.S_Tstp_rot = c.R1*(3./2. + 1./2. *
_np.log(pi*T**3/self.T_I**3*I) -
_np.log(self.sigma))
self.q_rot = _np.sqrt(_np.pi*I)/self.sigma *\
(T/self.T_I)**(3./2.)
else:
'''
Linear species
'''
I = _np.max(self.I3)
self.S_Tstp_rot = c.R1*(1. + _np.log(T/self.T_I*I) -
_np.log(self.sigma))
self.q_rot = (T*I/self.sigma)/self.T_I
p = 100000 # Presure of 1 atm or 100000 Pa
self.S_Tstp_trans = c.R1*(5./2. + 3./2. *
_np.log(2.*pi*self.MW/c.h1**2) +
5./2.*_np.log(c.kb1*T) -
_np.log(p))
if hasattr(self, 'A_st'):
self.q_trans2D = self.A_st * (2*pi*self.MW*c.kb1*T)/c.h1**2
else:
'''
Surface phase calculation
'''
self.S_Tstp_rot = 0.
self.S_Tstp_trans = 0.
self.q_rot = 0.
'''
Sum of all contributions to entropy for total entropy
'''
self.S_Tstp = self.S_Tstp_vib + self.S_Tstp_rot + self.S_Tstp_trans
def Calc_HeatCapacities(self):
'''
Calculate the vibrational, rotational and translational components and
total heat capacity of a species for a range of temperatures
'''
'''
Calculate vibrational contribution to heat capacity for temperature
range specified in Cp_Range for linear and non-linear species
'''
zz = []
for x in range(0, len(self.Cp_Range)):
zz.append(_np.divide(self.theta, self.Cp_Range[x]))
zz = _np.array(zz).T
self.Cp_vib = sum(_np.divide(_np.multiply(_np.power(zz, 2),
_np.exp(-zz)),
_np.power(1-_np.exp(-zz), 2)))
if self.phase == 'G':
'''
Translational and rotational Gas phase calculation
'''
self.Cp_trans = _np.array([3./2.]*len(self.Cp_Range))
if self.islinear == 0:
'''
Non-Linear species
'''
self.Cp_rot = _np.array([3./2.]*len(self.Cp_Range))
else:
'''
Linear species
'''
self.Cp_rot = _np.array([1.]*len(self.Cp_Range))
else:
'''
Surface species
'''
self.Cp_rot = _np.array([0.]*len(self.Cp_Range))
self.Cp_trans = _np.array([0.]*len(self.Cp_Range))
'''
Sum all contribution to heat capacity for total heat capapcity
'''
self.Cp = c.R1*(self.Cp_trans + self.Cp_rot + self.Cp_vib + 1)
def Calc_Enthalpy(self):
T = self.Tstp
'''
Calculate zero-point energy
'''
self.zpe = sum(_np.multiply(c.h2, self.nu)/2)*c.NA/1000
'''
Calculate vibrational component of enthalpy
'''
self.E_Tstp_vib = c.kb2 *\
sum(_np.divide(self.theta*_np.exp(-self.theta/T),
(1 - _np.exp(-self.theta/T)))) *\
c.NA/1000
'''
Calculate translational and rotational component of enthalpy
'''
if self.phase == 'G':
'''
Gas phase calculation
'''
self.E_Tstp_trans = 3./2.*c.R1*T/1000
if self.islinear == 0:
'''
Non-linear species
'''
self.E_Tstp_rot = 3./2.*c.R1*T/1000
else:
'''
Linear species
'''
self.E_Tstp_rot = 1.*c.R1*T/1000
self.pv_Tstp = 1.*c.R1*T/1000
else:
'''
Surface phase calculation
'''
self.E_Tstp_trans = 0.
self.E_Tstp_rot = 0.
self.pv_Tstp = 0.
'''
Sum all contribution to enthalpy for total enthalpy
'''
self.dfth = self.etotal*c.R1/c.R2*c.NA/1000 + self.zpe +\
self.E_Tstp_vib + self.E_Tstp_trans + self.E_Tstp_rot +\
self.pv_Tstp
[docs]class Reference(Particle):
'''
SubClass object to add specific fields for reference species
'''
def __init__(self, data, dict, Base_path, Tstp=298.15):
if data[dict['sigma']] != '':
self.sigma = int(data[dict['sigma']]) # Sigma
else:
self.sigma = int(0)
if data[dict['islinear']] != '':
self.islinear = int(data[dict['islinear']]) # Is molecule linear?
else:
self.linear = int(-1)
if 'hf298nist' in dict:
self.hf298nist = float(data[dict['hf298nist']]) # NIST Std enthpy
if 'inertia' in dict:
if data[dict['inertia']] != '':
self.Inertia = float(data[dict['inertia']])
else:
self.Inertia = float(0)
self.phase = str.upper(data[dict['phase']]) # Phase
if 'a_st' in dict and data[dict['a_st']] != '':
self.A_st = float(data[dict['a_st']])
super(Reference, self).__init__(data, dict, Base_path, Tstp=298.15)
@staticmethod
def BasisSet(RefSpecies):
A = []
b_nist = []
b_dfth = []
for x in range(0, len(RefSpecies)):
A.append([RefSpecies[x].carbon, RefSpecies[x].hydrogen,
RefSpecies[x].oxygen, RefSpecies[x].nitrogen])
b_nist.append([RefSpecies[x].hf298nist])
b_dfth.append([RefSpecies[x].dfth])
ref = _np.linalg.lstsq(A, _np.array(b_nist) - _np.array(b_dfth))[0]
return(ref)
[docs]class Target(Particle):
'''
SubClass object to add specific fields for target surface species
'''
def __init__(self, data, dict, Base_path, Tstp=298.15):
self.surface = str(data[dict['surface']]) # Surface
self.functional = str(data[dict['functional']]) # Functional
self.kpoints = str(data[dict['kpoints']]) # k-Points
self.vibfreqpath = str(data[dict['vibfreqpath']]) # Unused
self.phase = None # Phase
super(Target, self).__init__(data, dict, Base_path, Tstp=298.15)
@staticmethod
def ReferenceDFT(Species, Surface, Basis):
for x in range(0, len(Species)):
Molecule = _np.array([Species[x].carbon, Species[x].hydrogen,
Species[x].oxygen, Species[x].nitrogen])
if Species[x].phase == 'G':
Species[x].hf_Tstp = (Species[x].dfth +
_np.dot(Molecule, Basis))[0]
if hasattr(Species[x], 'edisp'):
Species[x].convedisp = (Species[x].edisp *
c.ev_atom_2_kcal_mol)
else:
Slab = next((y for y in Surface if y.name ==
Species[x].surface),
None)
if Slab is None:
print 'Error'
else:
Species[x].hf_Tstp = (Species[x].dfth +
_np.dot(Molecule, Basis) -
Slab.etotal *
c.ev_atom_2_kcal_mol)[0]
if hasattr(Species[x], 'edisp'):
Species[x].convedisp = (Species[x].edisp *
c.ev_atom_2_kcal_mol -
Slab.edisp *
c.ev_atom_2_kcal_mol)
return(Species)
@staticmethod
[docs] def CreateThermdat(Species, Base_path, Output):
'''
Calculate the seven coefficients for the NASA polynomials for two
temperature ranges and output the results in a Chemkin format thermdat
file
'''
T_mid = 500
Tstp = Species[0].Tstp
def HS_NASA(T, a):
'''
7-coefficient NASA polynomials for enthalpy and entropy
'''
Enthalpy = a[0] + a[1]*T/2 + a[2]*T**2/3 + \
a[3]*T**3/4 + a[4]*T**4/5
Entropy = a[0]*_np.log(T) + a[1]*T + a[2]*T**2/2 + \
a[3]*T**3/3 + a[4]*T**4/4
return[Enthalpy, Entropy]
for x in range(0, len(Species)):
T_rng_low = _np.linspace(min(Species[x].Cp_Range), T_mid, 1600)
T_rng_high = _np.linspace(T_mid, max(Species[x].Cp_Range), 4000)
T_func = _sp.InterpolatedUnivariateSpline(Species[x].Cp_Range,
Species[x].Cp/c.R1, k=4)
'''
Fit coefficients A1-A5 to heat capacity data
'''
Species[x].a_low = _np.polyfit(T_rng_low,
T_func(T_rng_low), 4)[::-1]
Species[x].a_high = _np.polyfit(T_rng_high,
T_func(T_rng_high), 4)[::-1]
'''
Correct A1 high temperature range coefficient to eliminate
discontinuity between high and low temperature range polynomials
'''
Species[x].a_high[0] = Species[x].a_high[0] + \
(_np.polyval(Species[x].a_low[::-1], T_mid) -
_np.polyval(Species[x].a_high[::-1], T_mid))
'''
Determine A6 coefficient for enthalpy calculations
'''
a6_high = (Species[x].hf_Tstp/c.R1/Tstp*1000 -
HS_NASA(Tstp, Species[x].a_high)[0])*Tstp
a6_low = (Species[x].hf_Tstp/c.R1/Tstp*1000 -
HS_NASA(Tstp, Species[x].a_low)[0])*Tstp
'''
Correct A6 high temperature range coefficient to eliminate
discontinuity between high and low temperature range polynomials
'''
a6_high_delta = (HS_NASA(T_mid, Species[x].a_low)[0] +
a6_low/T_mid) - \
(HS_NASA(T_mid,
Species[x].a_high)[0] + a6_high/T_mid)
a6_high = a6_high + a6_high_delta * T_mid
Species[x].a_high = _np.append(Species[x].a_high, a6_high)
Species[x].a_low = _np.append(Species[x].a_low, a6_low)
'''
Determine A7 coefficient for entropy calculations
'''
a7_high = Species[x].S_Tstp/c.R1 - \
HS_NASA(Tstp, Species[x].a_high)[1]
a7_low = Species[x].S_Tstp/c.R1 - \
HS_NASA(Tstp, Species[x].a_low)[1]
'''
Correct A7 high temperature range coefficient to eliminate
discontinuity between high and low temperature range polynomials
'''
a7_high_delta = (HS_NASA(T_mid, Species[x].a_low)[1] +
a7_low) - (HS_NASA(T_mid,
Species[x].a_high)[1] + a7_high)
a7_high = a7_high + a7_high_delta
Species[x].a_high = _np.append(Species[x].a_high, a7_high)
Species[x].a_low = _np.append(Species[x].a_low, a7_low)
'''
Write the species name, seven NASA coefficients for both a high and
a low temperature range and other data in the Chemkin thermdat
file format
'''
if os.path.isdir(os.path.join(Base_path, Output)) is False:
os.mkdir(os.path.join(Base_path, Output))
filepath = os.path.join(Base_path, Output, 'thermdat')
fid = open(filepath, 'w')
fid.truncate()
'''
Write thermdat file header
'''
fid.write('THERMO ALL\n')
for s in range(0, _np.size(Species)):
'''
Write header line for each species on line 1
'''
fid.write('%-16s' % (Species[s].name))
fid.write('%-8s' % (datetime.date.today().strftime('%Y%m%d')))
fid.write('%1s%-4i' % ('C', Species[s].carbon))
fid.write('%1s%-4i' % ('O', Species[s].oxygen))
fid.write('%1s%-4i' % ('H', Species[s].hydrogen))
fid.write('%1s%-4i' % ('N', Species[s].nitrogen))
if Species[s].name.find('(S)'):
fid.write('S')
else:
fid.write('G')
fid.write('%10.0f%10.0f%8.0f' % (min(Species[x].Cp_Range),
max(Species[x].Cp_Range),
T_mid))
fid.write('%6s%1i\n' % ('', 1))
'''
Write first five NASA coefficients for
low temperature range on line 2
'''
for x in range(0, 5):
fid.write('%15E' % (Species[s].a_low[x]))
fid.write('%4s%1i\n' % ('', 2))
'''
Write final two NASA coefficients for
low temperature range on line 2
'''
for x in range(0, 2):
fid.write('%15E' % (Species[s].a_low[x+4]))
'''
Write first three NASA coeficients for
high temperature range on line 3
'''
for x in range(0, 3):
fid.write('%15E' % (Species[s].a_high[x]))
fid.write('%4s%1i\n' % ('', 3))
'''
Write final four NASA coefficients for
high temperature range on line 4
'''
for x in range(0, 4):
fid.write('%15E' % (Species[s].a_high[x+3]))
fid.write('%19s%1i\n' % ('', 4))
'''
Write file footer and close the file
'''
fid.write('END\n')
fid.close()
return(Species)
[docs]class Surface:
'''
Class object to populate slab energies for surfaces
'''
def __init__(self, data, dict):
self.name = str(data[dict['name']]) # Surface name
self.etotal = float(data[dict['etotal']]) # Total energy-DFT
self.edisp = float(data[dict['edisp']]) # Dispersion energy-DFT
def DFTFileRead(filepath):
fid = open(filepath, 'r')
file = fid.read()
lines = file.splitlines()
dict_array = lines[2].lower().split('\t')
dict = {}
for x in range(0, len(dict_array)):
dict[dict_array[x]] = x
return(lines, dict)