From 4a0849e87bf3675e1b55be4f7239e8aa124eab85 Mon Sep 17 00:00:00 2001 From: mieand Date: Wed, 13 Jan 2016 16:55:37 +0100 Subject: [PATCH 1/5] Added thermochemistry corrections and species data --- kmos/species.py | 98 +++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 79 insertions(+), 19 deletions(-) diff --git a/kmos/species.py b/kmos/species.py index cb6b33de..c7d4c876 100644 --- a/kmos/species.py +++ b/kmos/species.py @@ -1,4 +1,8 @@ #!/usr/bin/env python +from ase.structure import molecule +from ase.thermochemistry import IdealGasThermo +from ase.thermochemistry import HarmonicThermo + """Very simple module that keeps several species commonly needed to model heterogeneous catalyst kinetics """ @@ -58,9 +62,19 @@ in the tab-delimited text format to the `janaf_data` directory.""") +def GibbsAds(energy, frequencies, T): + #Expecting T in Kelvin + #note that frequencies should have a lower bound, e.g. 56 cm-1, in order to bound entropic contributions. + cm_to_eV=1.23981e-4 + vib_energies=list(frequencies) + for i in range(len(vib_energies)): + vib_energies[i]=vib_energies[i]*cm_to_eV + thermo_ads=HarmonicThermo(vib_energies=vib_energies, electronicenergy=energy) + val=thermo_ads.get_gibbs_energy(temperature=T,verbose=False) + return val class Species(object): - def __init__(self, atoms, gas=False, janaf_file='', name=''): + def __init__(self, atoms, gas=False, janaf_file='', name='', symmetrynumber=0, frequencies=[], geometry='', spin=0): self.atoms = atoms self.gas = gas if name: @@ -71,6 +85,10 @@ def __init__(self, atoms, gas=False, janaf_file='', name=''): else: self.atoms.get_name() self.janaf_file = janaf_file + self.symmetrynumber = symmetrynumber + self.frequencies = frequencies + self.geometry = geometry + self.spin = spin # prepare chemical potential if self.gas and self.janaf_file and janaf_data is not None: @@ -99,6 +117,28 @@ def mu(self, T, p): else: raise UserWarning('%s is no gas-phase species.' % self.name) + def GibbsGas(self, energy, T, p): + #Expecting T in Kelvin, p in bar + #note that frequencies should have a lower bound, e.g. 56 cm-1, in order to bound entropic contributions. + if self.gas: + cm_to_eV=1.23981e-4 + vib_energies=list(self.frequencies) + for i in range(len(vib_energies)): + vib_energies[i]=vib_energies[i]*cm_to_eV + thermo_gas=IdealGasThermo(vib_energies=vib_energies, + electronicenergy=energy, + atoms=self.atoms, + geometry=self.geometry, + symmetrynumber=self.symmetrynumber, + spin=self.spin) + bar_to_Pa=1e5 + pressure=p*bar_to_Pa + val=thermo_gas.get_gibbs_energy(temperature=T, pressure=pressure, verbose=False) + return val + + else: + raise UserWarning('%s is no gas-phase species.' % self.name) + def _prepare_G_p0(self, filename): # from CODATA 2010 Jmol_in_eV = 1.03642E-5 @@ -125,22 +165,32 @@ def __hash__(self): # prepare all required species -H2gas = Species(ase.atoms.Atoms('H2', [[0, 0, 0], [0, 0, 1.2]],), + +#frequencies are experimental data from NIST CCCBDB and are in cm-1 +#geometry should be 'monatomic', 'linear', or 'nonlinear' +#symmetry numbers can be found in Table 10.1 and Appendix B of C. Cramer "Essentials of Computational Chemistry", 2nd Ed. +#see ase.thermochemistry module + +H2gas = Species(molecule('H2'), gas=True, janaf_file='H-050.txt', - name='H2gas') + name='H2gas', + frequencies=[4401], + geometry='linear', + symmetrynumber=2, + spin=0) H = Species(ase.atoms.Atoms('H')) -CH4gas = Species(ase.atoms.Atoms('CH4', - [[-2.14262, 3.03116, 0.00000], - [-1.07262, 3.03116, 0.00000], - [-2.49979, 4.03979, 0.00000], - [-2.51306, 2.50700, 0.85611], - [-2.49435, 2.53348, -0.87948]], - ), gas=True, - janaf_file='C-067.txt', - name='CH4gas') +CH4gas = Species(molecule('CH4'), + gas=True, + janaf_file='C-067.txt', + name='CH4gas', + frequencies=[2917,1534,1534,3019,3019,3019,1306,1306,1306], + geometry='nonlinear', + symmetrynumber=12, + spin=0) + CH4 = Species(ase.atoms.Atoms('CH4', [[-2.14262, 3.03116, 0.00000], [-1.07262, 3.03116, 0.00000], @@ -149,11 +199,13 @@ def __hash__(self): [-2.49435, 2.53348, -0.87948]], ), name='CH4') + O = Species(ase.atoms.Atoms('O', [[0, 0, 0]], cell=[10, 10, 10], ), name='O') + O2gas = Species(ase.atoms.Atoms('O2', [[0, 0, 0], [0, 0, 1.2]], @@ -188,11 +240,15 @@ def __hash__(self): janaf_file='N-009.txt', name='NO3gas') -COgas = Species(ase.atoms.Atoms('CO', [[0, 0, 0], [0, 0, 1.2]], - cell=[10, 10, 10],), +COgas = Species(molecule('CO'), gas=True, janaf_file='C-093.txt', - name='COgas') + name='COgas', + frequencies=[2170], + geometry='linear', + symmetrynumber=1, + spin=0) + CO = Species(ase.atoms.Atoms('CO', [[0, 0, 0], [0, 0, 1.2]], cell=[10, 10, 10], ), name='CO') @@ -237,10 +293,14 @@ def __hash__(self): janaf_file='Cl-073.txt', name='Cl2gas',) -H2Ogas = Species(ase.atoms.Atoms(), - gas=True, - janaf_file='H-064.txt', - name='H2Ogas',) +H2Ogas = Species(molecule('H2O'), + gas=True, + janaf_file='H-064.txt', + name='H2Ogas', + frequencies=[3657,1595,3756], + geometry='nonlinear', + symmetrynumber=2, + spin=0) H2Oliquid = Species(ase.atoms.Atoms(), gas=False, From 530dad2034fa4bb19e0413564bd0a2cf5b139fe0 Mon Sep 17 00:00:00 2001 From: mieand Date: Thu, 21 Jan 2016 15:23:30 +0100 Subject: [PATCH 2/5] Allow for scaling expressions in parameters --- kmos/__init__.py | 75 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) diff --git a/kmos/__init__.py b/kmos/__init__.py index b3cdb440..87f64c2a 100644 --- a/kmos/__init__.py +++ b/kmos/__init__.py @@ -124,6 +124,81 @@ def evaluate_rate_expression(rate_expr, parameters={}): print('No JANAF table assigned for %s' % species_name) print('Setting chemical potential to zero') replaced_tokens.append((i, '0')) + + elif token.startswith('GibbsGas_'): + #evaluate gas phase gibbs free energy using ase thermochemistry module, + #experimental data from NIST CCCBDB, the electronic energy + #and current temperature and partial pressure + from kmos import species + species_name = '_'.join(token.split('_')[1:]) + if species_name in dir(species): + if not 'T' in parameters: + raise Exception('Need "T" in parameters to evaluate gas phase gibbs free energy.') + + if not ('p_%s' % species_name) in parameters: + raise Exception('Need "p_%s" in parameters to evaluate gas phase gibbs free energy.' % species_name) + + replaced_tokens.append((i, 'species.%s.GibbsGas(%s,%s,%s)' % ( + species_name, + parameters['E_'+species_name]['value'], + parameters['T']['value'], + parameters['p_%s' % species_name]['value'], + ))) + else: + print('No NIST data assigned for %s' % species_name) + print('Setting chemical potential to zero') + replaced_tokens.append((i, '0')) + + #gibbs=eval(replaced_tokens[-1][-1]) + #print species_name+': %.3f'%gibbs + + elif token.startswith('GibbsAds_'): + #evaluate gibbs free energy of adsorbate using ase thermochemistry module, + #calculated frequencies and electronic energy and current temperature + from kmos import species + species_name = '_'.join(token.split('_')[1:]) + if not 'T' in parameters: + raise Exception('Need "T" in parameters to evaluate adsorbate gibbs free energy.') + energy=parameters['E_'+species_name]['value'] + try: + eval(energy) + except: + try: + replaced_tokens2=[] + input = StringIO.StringIO(energy).readline + tokens2 = list(tokenize.generate_tokens(input)) + except: + raise Exception('Could not tokenize expression: %s' % input) + for j, token2, _, _, _ in tokens2: + if token2 in parameters: + parameter_str = str(parameters[token2]['value']) + try: + eval(parameter_str) + replaced_tokens2.append((j, parameter_str)) + except: + try: + input = StringIO.StringIO(parameter_str).readline + tokens3 = list(tokenize.generate_tokens(input)) + except: + raise Exception('Could not tokenize expression: %s' % input) + for k, token3, _, _, _ in tokens3: + if token3 in parameters: + parameter_str = str(parameters[token3]['value']) + replaced_tokens2.append((k, parameter_str)) + else: + replaced_tokens2.append((k, token3)) + else: + replaced_tokens2.append((j, token2)) + energy = tokenize.untokenize(replaced_tokens2) + replaced_tokens.append((i, 'species.GibbsAds(%s,%s,%s)' % ( + energy, + parameters['f_'+species_name]['value'], + parameters['T']['value'], + ))) + + #gibbs=eval(replaced_tokens[-1][-1]) + #print species_name+': %.3f'%gibbs + elif token in parameters: parameter_str = str(parameters[token]['value']) # replace units used in parameters From e2d85c08cb74527db5a17875b87589ec7af0dac0 Mon Sep 17 00:00:00 2001 From: mieand Date: Fri, 22 Jan 2016 11:10:30 +0100 Subject: [PATCH 3/5] Added poss. to use analyt. exp. in param. --- kmos/__init__.py | 100 ++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 99 insertions(+), 1 deletion(-) diff --git a/kmos/__init__.py b/kmos/__init__.py index 87f64c2a..d1b8b00d 100644 --- a/kmos/__init__.py +++ b/kmos/__init__.py @@ -201,11 +201,109 @@ def evaluate_rate_expression(rate_expr, parameters={}): elif token in parameters: parameter_str = str(parameters[token]['value']) + # replace some aliases + parameter_str = parameter_str.replace('beta', '(1./(kboltzmann*T))') # replace units used in parameters for unit in units.keys: parameter_str = parameter_str.replace( unit, '%s' % eval('units.%s' % unit)) - replaced_tokens.append((i, parameter_str)) + try: + eval(parameter_str) + replaced_tokens.append((i, parameter_str)) + except: + try: + replaced_tokens2=[] + input = StringIO.StringIO(parameter_str).readline + tokens2 = list(tokenize.generate_tokens(input)) + except: + raise Exception('Could not tokenize expression: %s' % input) + for i, token2, _, _, _ in tokens2: + if token2 in ['sqrt', 'exp', 'sin', 'cos', 'pi', 'pow', 'log']: + replaced_tokens2.append((i, 'math.' + token2)) + elif token2.startswith('GibbsGas_'): + #evaluate gas phase gibbs free energy using ase thermochemistry module, + #experimental data from NIST CCCBDB, the electronic energy + #and current temperature and partial pressure + from kmos import species + species_name = '_'.join(token2.split('_')[1:]) + if species_name in dir(species): + if not 'T' in parameters: + raise Exception('Need "T" in parameters to evaluate gas phase gibbs free energy.') + + if not ('p_%s' % species_name) in parameters: + raise Exception('Need "p_%s" in parameters to evaluate gas phase gibbs free energy.' % species_name) + + replaced_tokens2.append((i, 'species.%s.GibbsGas(%s,%s,%s)' % ( + species_name, + parameters['E_'+species_name]['value'], + parameters['T']['value'], + parameters['p_%s' % species_name]['value'], + ))) + else: + print('No NIST data assigned for %s' % species_name) + print('Setting chemical potential to zero') + replaced_tokens2.append((i, '0')) + + #gibbs=eval(replaced_tokens2[-1][-1]) + #print species_name+': %.3f'%gibbs + + elif token2.startswith('GibbsAds_'): + #evaluate gibbs free energy of adsorbate using ase thermochemistry module, + #calculated frequencies and electronic energy and current temperature + from kmos import species + species_name = '_'.join(token2.split('_')[1:]) + if not 'T' in parameters: + raise Exception('Need "T" in parameters to evaluate adsorbate gibbs free energy.') + energy=parameters['E_'+species_name]['value'] + try: + eval(energy) + except: + try: + replaced_tokens3=[] + input = StringIO.StringIO(energy).readline + tokens3 = list(tokenize.generate_tokens(input)) + except: + raise Exception('Could not tokenize expression: %s' % input) + for j, token3, _, _, _ in tokens3: + if token3 in parameters: + parameter_str = str(parameters[token3]['value']) + try: + eval(parameter_str) + replaced_tokens3.append((j, parameter_str)) + except: + try: + input = StringIO.StringIO(parameter_str).readline + tokens4 = list(tokenize.generate_tokens(input)) + except: + raise Exception('Could not tokenize expression: %s' % input) + for k, token4, _, _, _ in tokens4: + if token4 in parameters: + parameter_str = str(parameters[token4]['value']) + replaced_tokens3.append((k, parameter_str)) + else: + replaced_tokens3.append((k, token4)) + else: + replaced_tokens3.append((j, token3)) + energy = tokenize.untokenize(replaced_tokens3) + replaced_tokens2.append((i, 'species.GibbsAds(%s,%s,%s)' % ( + energy, + parameters['f_'+species_name]['value'], + parameters['T']['value'], + ))) + + #gibbs=eval(replaced_tokens2[-1][-1]) + #print species_name+': %.3f'%gibbs + + elif token2 in parameters: + replaced_tokens2.append((i, str(parameters[token2]['value']))) + else: + replaced_tokens2.append((i, token2)) + parameter_str = tokenize.untokenize(replaced_tokens2) + + #print parameter_str + #print token+': %.7f'%eval(parameter_str) + + replaced_tokens.append((i, parameter_str)) else: replaced_tokens.append((i, token)) From 4aa9af7e61adff21497c0f510358b6175d3e30ec Mon Sep 17 00:00:00 2001 From: mieand Date: Mon, 25 Jan 2016 14:43:05 +0100 Subject: [PATCH 4/5] Added function to evaluate param. exp. --- kmos/__init__.py | 115 +++++++++++++++++++++++++++++++++++++++++++ kmos/run/__init__.py | 5 ++ 2 files changed, 120 insertions(+) diff --git a/kmos/__init__.py b/kmos/__init__.py index d1b8b00d..cf2614bf 100644 --- a/kmos/__init__.py +++ b/kmos/__init__.py @@ -53,6 +53,121 @@ __version__ = "0.3.16" VERSION = __version__ +def evaluate_param_expression(param, parameters={}): + import tokenize + import StringIO + import math + from kmos import units + + # convert parameters to dict if passed as list of Parameters() + if type(parameters) is list: + param_dict = {} + for parameter in parameters: + param_dict[parameter.name] = {'value': parameter.value} + parameters = param_dict + + parameter_str = str(parameters[param]['value']) + + # replace some aliases + parameter_str = parameter_str.replace('beta', '(1./(kboltzmann*T))') + # replace units used in parameters + for unit in units.keys: + parameter_str = parameter_str.replace( + unit, '%s' % eval('units.%s' % unit)) + try: + return eval(parameter_str) + except: + try: + replaced_tokens=[] + input = StringIO.StringIO(parameter_str).readline + tokens = list(tokenize.generate_tokens(input)) + except: + raise Exception('Could not tokenize expression: %s' % input) + for i, token, _, _, _ in tokens: + if token in ['sqrt', 'exp', 'sin', 'cos', 'pi', 'pow', 'log']: + replaced_tokens.append((i, 'math.' + token)) + elif token.startswith('GibbsGas_'): + #evaluate gas phase gibbs free energy using ase thermochemistry module, + #experimental data from NIST CCCBDB, the electronic energy + #and current temperature and partial pressure + from kmos import species + species_name = '_'.join(token.split('_')[1:]) + if species_name in dir(species): + if not 'T' in parameters: + raise Exception('Need "T" in parameters to evaluate gas phase gibbs free energy.') + + if not ('p_%s' % species_name) in parameters: + raise Exception('Need "p_%s" in parameters to evaluate gas phase gibbs free energy.' % species_name) + + replaced_tokens.append((i, 'species.%s.GibbsGas(%s,%s,%s)' % ( + species_name, + parameters['E_'+species_name]['value'], + parameters['T']['value'], + parameters['p_%s' % species_name]['value'], + ))) + else: + print('No NIST data assigned for %s' % species_name) + print('Setting chemical potential to zero') + replaced_tokens.append((i, '0')) + + #gibbs=eval(replaced_tokens2[-1][-1]) + #print species_name+': %.3f'%gibbs + + elif token.startswith('GibbsAds_'): + #evaluate gibbs free energy of adsorbate using ase thermochemistry module, + #calculated frequencies and electronic energy and current temperature + from kmos import species + species_name = '_'.join(token.split('_')[1:]) + if not 'T' in parameters: + raise Exception('Need "T" in parameters to evaluate adsorbate gibbs free energy.') + energy=parameters['E_'+species_name]['value'] + try: + eval(energy) + except: + try: + replaced_tokens2=[] + input = StringIO.StringIO(energy).readline + tokens2 = list(tokenize.generate_tokens(input)) + except: + raise Exception('Could not tokenize expression: %s' % input) + for j, token2, _, _, _ in tokens2: + if token2 in parameters: + parameter_str2 = str(parameters[token2]['value']) + try: + eval(parameter_str2) + replaced_tokens2.append((j, parameter_str2)) + except: + try: + input = StringIO.StringIO(parameter_str2).readline + tokens3 = list(tokenize.generate_tokens(input)) + except: + raise Exception('Could not tokenize expression: %s' % input) + for k, token3, _, _, _ in tokens3: + if token3 in parameters: + parameter_str3 = str(parameters[token3]['value']) + replaced_tokens2.append((k, parameter_str3)) + else: + replaced_tokens2.append((k, token3)) + else: + replaced_tokens2.append((j, token2)) + energy = tokenize.untokenize(replaced_tokens2) + replaced_tokens.append((i, 'species.GibbsAds(%s,%s,%s)' % ( + energy, + parameters['f_'+species_name]['value'], + parameters['T']['value'], + ))) + + #gibbs=eval(replaced_tokens2[-1][-1]) + #print species_name+': %.3f'%gibbs + + elif token in parameters: + replaced_tokens.append((i, str(parameters[token]['value']))) + else: + replaced_tokens.append((i, token)) + parameter_str = tokenize.untokenize(replaced_tokens) + return eval(parameter_str) + #print parameter_str + #print token+': %.7f'%eval(parameter_str) def evaluate_rate_expression(rate_expr, parameters={}): """Evaluates an expression for a typical kMC rate constant. diff --git a/kmos/run/__init__.py b/kmos/run/__init__.py index 5962fc0a..65e3572e 100644 --- a/kmos/run/__init__.py +++ b/kmos/run/__init__.py @@ -43,6 +43,7 @@ from copy import deepcopy from fnmatch import fnmatch from kmos import evaluate_rate_expression +from kmos import evaluate_param_expression from kmos.utils import OrderedDict import kmos.run.png from math import log @@ -279,6 +280,10 @@ def get_param_header(self): for param_name in sorted(self.settings.parameters) if self.settings.parameters[param_name].get('adjustable', False)) + def get_param_value(self,param): + """Return the evaluated value of a parameter""" + return evaluate_param_expression(param, settings.parameters) + def get_occupation_header(self): """Return the names of the fields returned by self.get_atoms().occupation. From c71768bdc2279087ddc6fc315ebbb39e645bfc9d Mon Sep 17 00:00:00 2001 From: mieand Date: Thu, 21 Jul 2016 13:52:40 +0200 Subject: [PATCH 5/5] Added data for O2gas --- kmos/species.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/kmos/species.py b/kmos/species.py index c7d4c876..b90db352 100644 --- a/kmos/species.py +++ b/kmos/species.py @@ -206,14 +206,14 @@ def __hash__(self): ), name='O') -O2gas = Species(ase.atoms.Atoms('O2', - [[0, 0, 0], - [0, 0, 1.2]], - cell=[10, 10, 10], - ), +O2gas = Species(molecule('O2'), gas=True, janaf_file='O-029.txt', - name='O2gas') + name='O2gas', + frequencies=[1580], + geometry='linear', + symmetrynumber=2, + spin=1) NOgas = Species(ase.atoms.Atoms('NO', [[0, 0, 0], [0, 0, 1.2]],