diff --git a/catmap/__init__.py b/catmap/__init__.py index 2c029ecd..406fa24d 100644 --- a/catmap/__init__.py +++ b/catmap/__init__.py @@ -1,26 +1,28 @@ -#Standard dependencies +# Standard dependencies import os import sys import inspect import time try: import cPickle as pickle -except: +except (ImportError, ModuleNotFoundError): import _pickle as pickle import re from copy import copy from string import Template -#Non-standard dependencies +# Non-standard dependencies import numpy as np try: from scipy.interpolate import InterpolatedUnivariateSpline as spline except ImportError: - def spline_wrapper(x_data, y_data, k=3): # input kwarg k is intentionally ignored + # input kwarg k is intentionally ignored. + def spline_wrapper(x_data, y_data, k=3): # behaves like scipy.interpolate.InterpolatedUnivariateSpline for k=1 def spline_func(x): - return np.interp(x, map(float,x_data), map(float,y_data)) # loss of precision here + # loss of precision here + return np.interp(x, map(float,x_data), map(float,y_data)) return spline_func spline = spline_wrapper diff --git a/catmap/analyze/analysis_base.py b/catmap/analyze/analysis_base.py index 19a59b06..87148524 100644 --- a/catmap/analyze/analysis_base.py +++ b/catmap/analyze/analysis_base.py @@ -111,38 +111,38 @@ class MapPlot: :param aspect: :type aspect: - :param subplots_adjust_kwargs: Dictionary of keyword arguments for adjusting matplotlib subplots + :param subplots_adjust_kwargs: Dictionary of keyword arguments for + adjusting matplotlib subplots :type subplots_adjust_kwargs: dict .. todo:: Some missing descriptions """ def __init__(self): - defaults = dict( - resolution_enhancement = 1, - min = None, - max = None, - n_ticks = 8, - plot_function = None, - colorbar = True, - colormap = plt.cm.jet, - axis_label_decimals = 2, - log_scale = False, - descriptor_labels = ['X_descriptor','Y_descriptor'], - default_descriptor_pt_args = {'marker':'o'}, - default_descriptor_label_args = {}, - descriptor_pt_args = {}, - descriptor_label_args = {}, - include_descriptors = False, - plot_size = 4, - aspect = None, - subplots_adjust_kwargs = {'hspace':0.35,'wspace':0.35, - 'bottom':0.15} - ) + defaults = dict(resolution_enhancement=1, + min=None, + max=None, + n_ticks=8, + plot_function=None, + colorbar=True, + colormap=plt.cm.YlGnBu_r, + axis_label_decimals=2, + log_scale=False, + descriptor_labels=['X_descriptor', 'Y_descriptor'], + default_descriptor_pt_args={'marker': 'o'}, + default_descriptor_label_args={}, + descriptor_pt_args={}, + descriptor_label_args={}, + include_descriptors=False, + plot_size=4, + aspect=None, + subplots_adjust_kwargs={'hspace': 0.35, + 'wspace': 0.35, + 'bottom': 0.15}) for key in defaults: val = defaults[key] - if not hasattr(self,key): - setattr(self,key,val) + if not hasattr(self, key): + setattr(self, key, val) elif getattr(self,key) is None: setattr(self,key,val) diff --git a/catmap/analyze/matrix_map.py b/catmap/analyze/matrix_map.py index 96b22542..f8c20ad2 100644 --- a/catmap/analyze/matrix_map.py +++ b/catmap/analyze/matrix_map.py @@ -1,6 +1,7 @@ from .analysis_base import * from .vector_map import * + class MatrixMap(VectorMap): """ .. todo:: __doc__ diff --git a/catmap/analyze/scaling.py b/catmap/analyze/scaling.py index 93870b8f..e80740a9 100644 --- a/catmap/analyze/scaling.py +++ b/catmap/analyze/scaling.py @@ -1,26 +1,27 @@ -from .analysis_base import * +from .analysis_base import ScalingPlot, ReactionModelWrapper -class ScalingAnalysis(ScalingPlot,ReactionModelWrapper): + +class ScalingAnalysis(ScalingPlot, ReactionModelWrapper): """ .. todo:: __doc__ """ - def __init__(self,reaction_model): + def __init__(self, reaction_model): self._rxm = reaction_model self.scaling_mode = 'linear' - ScalingPlot.__init__(self,self.descriptor_names, - self.descriptor_dict,self.surface_names, - self.parameter_dict, - self.plotter_scaling_function, - self.plotter_x_axis_function) + ScalingPlot.__init__(self, self.descriptor_names, + self.descriptor_dict, self.surface_names, + self.parameter_dict, + self.plotter_scaling_function, + self.plotter_x_axis_function) - def plotter_scaling_function(self,descriptors,**kwargs): + def plotter_scaling_function(self, descriptors, **kwargs): """ .. todo:: __doc__ """ descriptors = list(descriptors) - return self.scaler.get_electronic_energies(descriptors,**kwargs) + return self.scaler.get_electronic_energies(descriptors, **kwargs) - def plotter_x_axis_function(self,descriptors,**kwargs): + def plotter_x_axis_function(self, descriptors, **kwargs): """ .. todo:: __doc__ """ @@ -32,24 +33,24 @@ def plotter_x_axis_function(self,descriptors,**kwargs): self.coefficient_matrix = C else: C = self.coefficient_matrix - for ads,c in zip( - self.adsorbate_names+self.transition_state_names,C): + for ads, c in zip( + self.adsorbate_names+self.transition_state_names, C): x_dict[ads] = sum( - [c[i]*a for i,a in enumerate(descriptors)]) + c[-1] + [c[i]*a for i, a in enumerate(descriptors)]) + c[-1] lab = '' - for i,a in enumerate(self.descriptor_names): - c_print = round(c[i],self.descriptor_decimal_precision) + for i, a in enumerate(self.descriptor_names): + c_print = round(c[i], self.descriptor_decimal_precision) if c_print != 0: lab += str(c_print)+'*$E_{'+str(a)+'}$+' - const = round(c[-1],self.descriptor_decimal_precision) + const = round(c[-1], self.descriptor_decimal_precision) if const > 0: lab += '+' + str(const) elif const < 0: lab += '-' + str(-const) - lab = lab.replace('++','+') - lab = lab.replace('+-','-') + lab = lab.replace('++', '+') + lab = lab.replace('+-', '-') labels[ads] = lab - return x_dict,labels + return x_dict, labels def get_error(self): """ @@ -58,4 +59,3 @@ def get_error(self): if not self.scaling_error: self.plot(save=False) return self.scaling_error - diff --git a/catmap/analyze/vector_map.py b/catmap/analyze/vector_map.py index 8321c905..394a6021 100644 --- a/catmap/analyze/vector_map.py +++ b/catmap/analyze/vector_map.py @@ -1,6 +1,7 @@ from .analysis_base import * from string import Template + class VectorMap(MapPlot, ReactionModelWrapper): """ .. todo:: __doc__ diff --git a/catmap/api/ase_data.py b/catmap/api/ase_data.py index a0b55350..ae95d94e 100644 --- a/catmap/api/ase_data.py +++ b/catmap/api/ase_data.py @@ -49,8 +49,9 @@ and energies. """ import warnings -from os import listdir, mkdir +from os import mkdir import os.path +from uuid import uuid4 import numpy as np import ase.db from ase.atoms import string2symbols @@ -526,13 +527,23 @@ def _db2pes(self, fname, selection=[], site_specific=False): ens = self.bee.get_ensemble_perturbations(BEEFvdW_contribs) except AttributeError: ens = 0 # np.zeros(self.bee.size) - rxn_id = str(d.path_id) + if 'path_id' in d: + rxn_id = str(d.path_id) + else: + rxn_id = uuid4().hex + if 'distance' in d: + distance = float(d.distance) + else: + distance = np.nan + if 'step' in d: + step = int(d.step) + else: + step = np.nan if rxn_id in rxn_paths: rxn_paths[rxn_id]['pes'].append(abinitio_energy) rxn_paths[rxn_id]['dbids'].append(dbid) rxn_paths[rxn_id]['ens'].append(ens) - # try: - rxn_paths[rxn_id]['distance'].append(float(d.distance)) + rxn_paths[rxn_id]['distance'].append(distance) # except AttributeError: # d0 = c.get(rxn_paths[rxn_id]['dbids'][0]) # atoms0 = c.get_atoms(rxn_paths[rxn_id]['dbids'][0]) @@ -565,8 +576,8 @@ def _db2pes(self, fname, selection=[], site_specific=False): 'dbids': [dbid], 'ens': [ens], 'site': site, - 'images': [int(d.step)], - 'distance': [float(d.distance)]} + 'images': [step], + 'distance': [distance]} return rxn_paths def _mol2ref(self, references): @@ -785,52 +796,57 @@ def pes2ts(self, freq_path=None, rtol=1.1): key = '1_' + species + '_' + m + '_' + site images = self.rxn_paths[rxn_id]['images'] pes = np.array(self.rxn_paths[rxn_id]['pes']) - s = np.argsort(images) - # Look for local minima and maxima. - localmins = np.where(np.r_[True, pes[s][1:] < pes[s][:-1]] & - np.r_[pes[s][:-1] < pes[s][1:], True])[0] - localmaxs = np.where(np.r_[True, pes[s][1:] > pes[s][:-1]] & - np.r_[pes[s][:-1] > pes[s][1:], True])[0] - # Measure path roughness - differences = np.diff(pes[s]) - roughness = np.std(differences) - # For fixed bond length (drag) calculations - g1, g2 = species.split('-') - dbond = covalent_radii[atomic_numbers[g1[0]]] - if len(g2) > 0: - dbond += covalent_radii[atomic_numbers[g2[0]]] - else: - # Assume the bond is with the surface. + if len(images) > 1: + s = np.argsort(images) + # Look for local minima and maxima. + localmins = np.where(np.r_[True, pes[s][1:] < pes[s][:-1]] & + np.r_[pes[s][:-1] < pes[s][1:], True])[0] + localmaxs = np.where(np.r_[True, pes[s][1:] > pes[s][:-1]] & + np.r_[pes[s][:-1] > pes[s][1:], True])[0] + # Measure path roughness + differences = np.diff(pes[s]) + roughness = np.std(differences) + # For fixed bond length (drag) calculations + g1, g2 = species.split('-') + dbond = covalent_radii[atomic_numbers[g1[0]]] + if len(g2) > 0: + dbond += covalent_radii[atomic_numbers[g2[0]]] + else: + # Assume the bond is with the surface. + try: + dbond += covalent_radii[atomic_numbers[ + m.split('_')[0]]] + except KeyError: + print("Bond not defined.") + if len(np.unique(images)) != len(images): + warn = True + print('non unique image number!') + print('Warning!', species, m, roughness, + len(localmaxs), len(localmins), images) + continue + if (len(localmaxs) > 1 or + len(localmins) > 2 or + len(localmins) == 1): + warn = True try: - dbond += covalent_radii[atomic_numbers[m.split('_')[0]]] + shortest = np.nanmin(self.rxn_paths[rxn_id]['distance']) + if shortest > dbond * rtol: + warn = True + s_last = np.argmax(self.rxn_paths[rxn_id]['images']) + calculate.append( + self.rxn_paths[rxn_id]['dbids'][s_last]) + continue except KeyError: - print("Bond not defined.") - if len(np.unique(images)) != len(images): - warn = True - print('non unique image number!') - print('Warning!', species, m, roughness, - len(localmaxs), len(localmins), images) - continue - if len(localmaxs) > 1 or len(localmins) > 2 or len(localmins) == 1: - warn = True - try: - shortest = np.min(self.rxn_paths[rxn_id]['distance']) - except KeyError: - print('Distances missing in reaction path.') - if shortest > dbond * rtol: - warn = True - s_last = np.argmax(self.rxn_paths[rxn_id]['images']) - calculate.append(self.rxn_paths[rxn_id]['dbids'][s_last]) - continue + print('Distances missing in reaction path.') + if warn: + warnings.warn("Warning! " + species + "*" + m + " " + + str(round(dbond * rtol, 3)) + " AA. " + + str(round(shortest, 3)) + " AA. " + + str(roughness) + " eV. " + + str(len(localmaxs)) + " local maxima. " + + str(len(localmins)) + " local minima. " + + str(len(images)) + " images.") tst = np.argmax(pes) - if warn: - warnings.warn("Warning! " + species + "*" + m + " " + - str(round(dbond * rtol, 3)) + " AA. " + - str(round(shortest, 3)) + " AA. " + - str(roughness) + " eV. " + - str(len(localmaxs)) + " local maxima. " + - str(len(localmins)) + " local minima. " + - str(len(images)) + " images.") if key not in abinitio_energies: abinitio_energies[key] = pes[tst] dbids[key] = self.rxn_paths[rxn_id]['dbids'][tst] @@ -1066,7 +1082,7 @@ def db_attach_formation_energy(self, fname, key_name): kvp = {key_name: float(self.formation_energies[key])} c.update(int(self.dbid[key]), **kvp) - def make_input_file(self, file_name, site_specific=False, + def make_input_file(self, file_name, site_specific='facet', catalyst_specific=False, covariance=None, reference=None): """ Saves the catmap input file. @@ -1193,8 +1209,18 @@ def make_input_file(self, file_name, site_specific=False, def make_nested_folders(self, project, reactions, surfaces=None, site='site', mol_db=None, slab_db=None, ads_db=None, ts_db=None, - publication='', url=''): + publication='', url='', xc='xc', code='code'): """Saves a nested directory structure. + The folder structure for catalysis-hub.org should be + + + + + + @ + .traj or + _.traj or + 'TS'.traj Parameters ---------- @@ -1219,12 +1245,16 @@ def make_nested_folders(self, project, reactions, surfaces=None, url : str url to publication. """ + data_folder = project + '/' + code + '/' + xc + # Create a header - spreadsheet = [['chemical_composition', 'facet', 'reactants', - 'products', 'reaction_energy', - 'beef_standard_deviation', - 'activation_energy', 'DFT_code', 'DFT_functional', - 'reference', 'url']] + # spreadsheet = [['chemical_composition', 'facet', 'reactants', + # 'products', 'reaction_energy', + # 'beef_standard_deviation', + # 'activation_energy', 'DFT_code', 'DFT_functional', + # 'reference', 'url']] + + # Connect to databases. if surfaces is None: surfaces = [s for s in self.reference_epot.keys() if 'slab' in s] if mol_db is None: @@ -1243,51 +1273,55 @@ def make_nested_folders(self, project, reactions, surfaces=None, c_slab = ase.db.connect(slab_db) if ts_db is not None: c_ts = ase.db.connect(ts_db) + + # Iterate over surfaces Nsurf = 0 Nrxn = 0 - # Loop over reaction expressions. - for i in range(len(reactions)): - # Separate left and right side of reaction, and ts if applicable. - states = reactions[i].replace(' ', '').split('<->') - if len(states) == 1: - states = states[0].split('->') + for slabkey in surfaces: + totraj = {} + [n, species, name, phase, + lattice, facet, cell, slab] = slabkey.split('_') + catalyst_name = name.replace('/', '') + '_' + phase + path_surface = data_folder + '/' + catalyst_name + facet_name = facet.replace('(', '').replace(')', '') + path_facet = path_surface + '/' + facet_name + + # Loop over reaction expressions. + for i, rxn in enumerate(reactions): + # Separate left and right side of reaction, and ts. + states = rxn.replace(' ', '').split('<->') if len(states) == 1: - states = states[0].split('<-') - elif len(states) < 3: - states = [states[0]] + states[-1].split('->') - if len(states) < 3: - states = states[0].split('<-') + states[1:] - # List individual species. - reactants = states[0].split('+') - tstates = states[1].split('+') - products = states[-1].split('+') - pspecies = [] - rspecies = [] - for specie in reactants: - if '_g' in specie: - rspecies.append(specie.split('_')[0] + 'gas') - elif '*_' in specie: - rspecies.append('star') - else: - rspecies.append(specie.split('_')[0] + 'star') - for specie in products: - if '_g' in specie: - pspecies.append(specie.split('_')[0] + 'gas') - elif '*_' in specie: - pspecies.append('star') - else: - pspecies.append(specie.split('_')[0] + 'star') - rname = '_'.join(rspecies) - pname = '_'.join(pspecies) - reaction_name = '__'.join([rname, pname]) - path_reaction = project + '/' + reaction_name - for slabkey in surfaces: - totraj = {} - [n, species, name, phase, - lattice, facet, cell, slab] = slabkey.split('_') - path_surface = path_reaction + '/' + name.replace('/', '') - facet_name = facet.replace('(', '').replace(')', '') - path_facet = path_surface + '/' + facet_name + states = states[0].split('->') + if len(states) == 1: + states = states[0].split('<-') + elif len(states) < 3: + states = [states[0]] + states[-1].split('->') + if len(states) < 3: + states = states[0].split('<-') + states[1:] + # List individual species. + reactants = states[0].split('+') + tstates = states[1].split('+') + products = states[-1].split('+') + pspecies = [] + rspecies = [] + for specie in reactants: + if '_g' in specie: + rspecies.append(specie.split('_')[0] + 'gas') + elif '*_' in specie: + rspecies.append('star') + else: + rspecies.append(specie.split('_')[0] + 'star') + for specie in products: + if '_g' in specie: + pspecies.append(specie.split('_')[0] + 'gas') + elif '*_' in specie: + pspecies.append('star') + else: + pspecies.append(specie.split('_')[0] + 'star') + rname = '_'.join(rspecies) + pname = '_'.join(pspecies) + reaction_name = '__'.join([rname, pname]) + path_reaction = path_facet + '/' + reaction_name DeltaE = 0. de = np.zeros(self.bee.size) ea = 0. @@ -1302,14 +1336,14 @@ def make_nested_folders(self, project, reactions, surfaces=None, n = 1 if sitesymbol == 'g': rkey = species + '_gas' - fname = path_reaction + '/' + rkey + fname = data_folder + '/gas/' + rkey elif species == '*': rkey = slabkey - fname = path_facet + '/empty_surface' + fname = path_reaction + '/empty_slab' else: rkey = '_'.join(['1', species.replace('*', ''), name, phase, lattice, facet, cell, site]) - fname = path_facet + '/' + species + fname = path_reaction + '/' + species if rkey not in self.dbid: intermediates_exist = False break @@ -1333,7 +1367,7 @@ def make_nested_folders(self, project, reactions, surfaces=None, continue totraj.update({tskey: {'dbid': self.dbid[tskey], - 'fname': path_facet + '/' + ts}}) + 'fname': path_reaction + '/' + ts}}) ea += self.formation_energies[tskey] - DeltaE # Find product structures and energies. for product in products: @@ -1345,14 +1379,14 @@ def make_nested_folders(self, project, reactions, surfaces=None, n = 1 if sitesymbol == 'g': pkey = species + '_gas' - fname = path_reaction + '/' + pkey + fname = data_folder + '/gas/' + pkey elif species == '*': pkey = slabkey - fname = path_facet + '/empty_surface' + fname = path_reaction + '/empty_slab' else: pkey = '_'.join(['1', species.replace('*', ''), name, phase, lattice, facet, cell, site]) - fname = path_facet + '/' + species + fname = path_reaction + '/' + species if pkey not in self.dbid: intermediates_exist = False break @@ -1364,23 +1398,12 @@ def make_nested_folders(self, project, reactions, surfaces=None, de += n * self.de_dict[pkey] # If all states are found for this surface, write. if intermediates_exist: - if not os.path.isdir(project): - mkdir(project) - if reaction_name not in listdir(project): - mkdir(path_reaction) - if name.replace('/', '') not in listdir(path_reaction): - mkdir(path_surface) - if facet_name not in listdir(path_surface): - mkdir(path_facet) # Loop over atomic structures to export. for trajkey in totraj.keys(): fname = totraj[trajkey]['fname'] + '.traj' if 'gas' in trajkey: - if fname in os.listdir(path_reaction): - continue - else: - atoms = c_mol.get_atoms(self.dbid[trajkey]) - d = c_mol.get(self.dbid[trajkey]) + atoms = c_mol.get_atoms(self.dbid[trajkey]) + d = c_mol.get(self.dbid[trajkey]) elif '-' in trajkey.split('_')[1]: atoms = c_ts.get_atoms(self.dbid[trajkey]) d = c_ts.get(self.dbid[trajkey]) @@ -1397,20 +1420,25 @@ def make_nested_folders(self, project, reactions, surfaces=None, atoms.set_calculator(calc) fname = fname.replace('(', '').replace(')', '') # Save trajectory file. + folder_structure = fname.split('/') + for depth in range(1, len(folder_structure)-1): + directory = '/'.join(folder_structure[:depth+1]) + if not os.path.isdir(directory): + mkdir(directory) atoms.write(fname) # Store rows for spreadsheet. - std = np.std(de) - spreadsheet.append([name, facet, rname, pname, - DeltaE, std, ea, - 'Quantum Espresso', - 'BEEF-vdW', publication, url]) + # std = np.std(de) + # spreadsheet.append([name, facet, rname, pname, + # DeltaE, std, ea, + # code, xc, + # publication, url]) # width, height, angle, covariance - Nsurf += 1 + Nrxn += 1 else: continue - Nrxn += 1 - with open(project + '/data.csv', 'wb') as f: - writer = csv.writer(f) - writer.writerows(spreadsheet) + Nsurf += 1 + # with open(project + '/data.csv', 'wb') as f: + # writer = csv.writer(f) + # writer.writerows(spreadsheet) print(Nrxn, 'reactions imported.') print(Nsurf, 'surfaces saved.') diff --git a/catmap/model.py b/catmap/model.py index ce52ed02..dead2072 100644 --- a/catmap/model.py +++ b/catmap/model.py @@ -4,6 +4,7 @@ from copy import copy import numpy as np import catmap +from catmap import griddata from string import Template from . import functions import re @@ -11,7 +12,7 @@ string2symbols = catmap.string2symbols pickle = catmap.pickle plt = catmap.plt -from catmap import griddata + class ReactionModel: """ @@ -26,41 +27,41 @@ class ReactionModel: - external parameters (temperature, pressures) - other more technical settings related to the solver and mapper """ - def __init__(self,**kwargs): # + def __init__(self, **kwargs): """Class for managing microkinetic models. :param setup_file: Specify from which to load model. :type setup_file: str """ - #Set static utility functions + # Set static utility functions. for f in dir(functions): - if not f.startswith('_') and callable(getattr(functions,f)): - setattr(self,f,getattr(functions,f)) + if not f.startswith('_') and callable(getattr(functions, f)): + setattr(self, f, getattr(functions, f)) - self._kB = 8.617332478e-5 #eV/K (from NIST) - self._h = 4.135667516e-15 #eV s (from NIST) + self._kB = 8.617332478e-5 # eV/K (from NIST) + self._h = 4.135667516e-15 # eV s (from NIST) self._default_site = 's' self._gas_sites = ['g'] - self._required = {'adsorbate_names':tuple, - 'transition_state_names':tuple, - 'gas_names':tuple, - 'descriptor_names':tuple, - 'surface_names':tuple, - 'species_definitions':dict, - 'elementary_rxns':tuple, - 'thermodynamics':None, - 'numerical_representation':None, - 'scaler':None, - 'solver':None, - 'mapper':None, - } #These are requirements of any reaction model. - - self._classes = ['parser','scaler', 'thermodynamics', - 'solver','mapper'] + # These are requirements of any reaction model. + self._required = {'adsorbate_names': tuple, + 'transition_state_names': tuple, + 'gas_names': tuple, + 'descriptor_names': tuple, + 'surface_names': tuple, + 'species_definitions': dict, + 'elementary_rxns': tuple, + 'thermodynamics': None, + 'numerical_representation': None, + 'scaler': None, + 'solver': None, + 'mapper': None} + + self._classes = ['parser', 'scaler', 'thermodynamics', + 'solver', 'mapper'] self._solved = None - #keeps track of whether or not the model was solved + # Keeps track of whether or not the model was solved. - #attributes for logging + # Attributes for logging self._log_lines = [] self._log_dict = {} @@ -68,14 +69,14 @@ def __init__(self,**kwargs): # 'loaded all data from input file', 'header_fail': 'could not save ${save_txt}'} - #modules to import in the log file to allow opening with python -i - self._log_imports = "from numpy import array\n\n"+\ + # Modules to import in the log file to allow opening with python -i + self._log_imports = "from numpy import array\n\n" + \ "try:\n" + \ " import cPickle as pickle\n\n" + \ "except: #workaround for python3.X\n" + \ " import _pickle as pickle\n\n" - #attrs taking up more space than this will be dumpted to self.data_file - self._max_log_line_length = 8000 #100 lines at 80 cols + # attrs taking up more space will be dumpted to self.data_file + self._max_log_line_length = 8000 # 100 lines at 80 cols #Character used for loop depth in std out self._pickle_attrs = [] #attributes to pickle rather than ASCII #place to store warnings @@ -122,10 +123,10 @@ def __init__(self,**kwargs): # #This is NOT idiot proof. self.model_name = self.setup_file.rsplit('.',1)[0] self.load(self.setup_file) - + # Functions for executing the kinetic model - def run(self,**kwargs): + def run(self, **kwargs): """Run the microkinetic model. If recalculate is True then data which is re-loaded will be used as an initial guess; otherwise it will be assumed to be correct. @@ -136,7 +137,7 @@ def run(self,**kwargs): """ for key in kwargs: - setattr(self,key,kwargs[key]) + setattr(self, key, kwargs[key]) # if 'all' is in output_variables, expand it to all processed output_variables # note that this expansion needs to happen in ReactionModel.run and not @@ -151,28 +152,29 @@ def run(self,**kwargs): ) )) - #ensure resolution has the proper dimensions - if not hasattr(self.resolution,'__iter__'): + # Ensure resolution has the proper dimensions. + if not hasattr(self.resolution, '__iter__'): self.resolution = [self.resolution]*len(self.descriptor_names) - #set numerical representation + # Set numerical representation. if self.numerical_representation == 'mpmath': import mpmath as mp - mp.mp.dps = self.decimal_precision + 25 #pad decimal precision + mp.mp.dps = self.decimal_precision + 25 # pad decimal precision self._math = mp self._log_imports += "\nfrom mpmath import mpf \n\n" - self._kB = mp.mpf('8.617332478e-5') #eV/K (from NIST) - self._h = mp.mpf('4.135667516e-15') #eV s (from NIST) + self._kB = mp.mpf('8.617332478e-5') # eV/K (from NIST) + self._h = mp.mpf('4.135667516e-15') # eV s (from NIST) self._mpfloat = mp.mpf self._matrix = mp.matrix self._Axb_solver = mp.lu_solve - self._math.infnorm = lambda x: mp.norm(x,'inf') - elif self.numerical_representation in ['numpy','python']: + self._math.infnorm = lambda x: mp.norm(x, 'inf') + elif self.numerical_representation in ['numpy', 'python']: self._log_imports += "\nfrom numpy import matrix \n\n" self._math = np - self._math.infnorm = lambda x: np.linalg.norm(x,np.inf) + self._math.infnorm = lambda x: np.linalg.norm(x, np.inf) self._mpfloat = float - def matrixT(*args,**kwargs): + + def matrixT(*args, **kwargs): array = np.array(args[0]) while 1 in array.shape: sum_idx = array.shape.index(1) @@ -180,28 +182,30 @@ def matrixT(*args,**kwargs): return array.T self._matrix = matrixT - def Axb_solver(A,b): + + def Axb_solver(A, b): try: - return np.linalg.solve(A,b.T) + return np.linalg.solve(A, b.T) except np.linalg.linalg.LinAlgError: raise ZeroDivisionError + self._Axb_solver = Axb_solver if self.decimal_precision > 15: - print( - 'Warning: Max precision with numpy/python is 16 digits') + print('Warning: Max precision with numpy/python is 16 digits') self.decimal_precision = 15 else: raise AttributeError( - 'Numerical representation must be mpmath, numpy, or python.') + 'Numerical representation must be mpmath, numpy, or python.') - #set up interaction model + # Set up interaction model. if self.adsorbate_interaction_model == 'first_order': - interaction_model = catmap.thermodynamics.FirstOrderInteractions(self) + interaction_model = \ + catmap.thermodynamics.FirstOrderInteractions(self) interaction_model.get_interaction_info() response_func = interaction_model.interaction_response_function if not callable(response_func): int_function = getattr(interaction_model, - response_func+'_response') + response_func+'_response') interaction_model.interaction_response_function = int_function self.thermodynamics.__dict__['adsorbate_interactions'] = interaction_model @@ -226,74 +230,79 @@ def Axb_solver(A,b): self.compatibility_check() - #determine whether or not to (re-)solve model + # Determine whether or not to (re-)solve model. has_all = True for v in self.output_variables: - if not hasattr(self,v+'_map'): + if not hasattr(self, v+'_map'): has_all = False - if not hasattr(self,'stdout'): - #any re-loaded model will have stdout + if not hasattr(self, 'stdout'): + # Any re-loaded model will have stdout. has_all = False - if self._solved == self._token() and not getattr(self,'recalculate',None): - #Do not solve the same model twice + if self._solved == self._token() and \ + not getattr(self, 'recalculate', None): + # Do not solve the same model twice. has_all = True - elif has_all and not getattr(self,'recalculate',None): - #All data has been loaded and no verification => solved. + elif has_all and not getattr(self, 'recalculate', None): + # All data has been loaded and no verification => solved. self._solved = self._token() self.log('input_success') else: - ran_dsa = False #When no map exists, run descriptor space analysis first - if not getattr(self,'coverage_map',None): - #Make "volcano plot" - if getattr(self,'descriptor_ranges',None) and getattr(self,'resolution',None): + ran_dsa = False + # When no map exists, run descriptor space analysis first. + if not getattr(self, 'coverage_map', None): + # Make "volcano plot". + if getattr(self, 'descriptor_ranges', None) and \ + getattr(self, 'resolution', None): self.descriptor_space_analysis() ran_dsa = True - #Get rates at single points - if getattr(self,'descriptors',None): + # Get rates at single points. + if getattr(self, 'descriptors', None): self.single_point_analysis(self.descriptors) - #Get rates at multiple points - if getattr(self,'descriptor_values',None): + # Get rates at multiple points. + if getattr(self, 'descriptor_values', None): self.multi_point_analysis() - #If a map exists, run descriptor space analysis last (so that single-point guesses are - #not discarded) - if not ran_dsa and getattr(self,'descriptor_ranges',None) and getattr(self,'resolution',None): - self.descriptor_space_analysis() + # If a map exists, run descriptor space analysis last + # (so that single-point guesses are not discarded). + if not ran_dsa and (getattr(self, 'descriptor_ranges', None) and + getattr(self, 'resolution', None)): + self.descriptor_space_analysis() - #Save long attrs in data_file + # Save long attrs in data_file. for attr in dir(self): if (not attr.startswith('_') and - not callable(getattr(self,attr)) and + not callable(getattr(self, attr)) and attr not in self._classes): - if (len(repr(getattr(self,attr))) > + if (len(repr(getattr(self, attr))) > self._max_log_line_length): - #line is too long for logfile -> put into pickle + # Line is too long for logfile -> put into pickle self._pickle_attrs.append(attr) pickled_data = {} for attr in self._pickle_attrs: - pickled_data[attr] = getattr(self,attr) + pickled_data[attr] = getattr(self, attr) try: - pickle.dump(pickled_data,open(self.data_file,'w')) - except: # Fallback workaround for Py3 - pickle.dump(pickled_data,open(self.data_file,'wb')) + pickle.dump(pickled_data, open(self.data_file, 'w')) + except: + # Fallback workaround for Py3 + pickle.dump(pickled_data, open(self.data_file, 'wb')) - #Make logfile + # Make logfile log_txt = self._log_imports log_txt += self._header(exclude_outputs=self._pickle_attrs) self._stdout = '\n'.join(self._log_lines) log_txt += 'stdout = ' + '"'*3+self._stdout+'"'*3 - #this construction means that self.stdout will only be set - #for models which have been re-loaded. + # This construction means that self.stdout will only be set + # for models which have been re-loaded. self._solved = self._token() if hasattr(self,'log_file'): logfile = self.log_file else: - name,suffix = self.setup_file.rsplit('.',1) + name, suffix = self.setup_file.rsplit('.', 1) if suffix != 'log': suffix = 'log' else: @@ -349,7 +358,7 @@ def generate_static_functions(self): #File IO functions - def load(self,setup_file): # + def load(self, setup_file): # """Load a 'setup file' by importing it and assigning all local variables as attributes of the kinetic model. Special attributes mapper, parser, scaler, solver will attempt to convert strings @@ -421,23 +430,23 @@ def load(self,setup_file): # self.verify() - def load_data_file(self,overwrite=False): + def load_data_file(self, overwrite=False): """ Load in output data from external files. """ if os.path.exists(self.data_file): try: - pickled_data = pickle.load(open(self.data_file,'r')) + pickled_data = pickle.load(open(self.data_file, 'r')) except: - pickled_data = pickle.load(open(self.data_file,'rb')) + pickled_data = pickle.load(open(self.data_file, 'rb')) for attr in pickled_data: if not overwrite: - if getattr(self,attr,None) is None: #don't over-write - setattr(self,attr,pickled_data[attr]) + if getattr(self, attr, None) is None: # Don't over-write + setattr(self, attr, pickled_data[attr]) else: - setattr(self,attr,pickled_data[attr]) + setattr(self, attr, pickled_data[attr]) - def parse(self,*args, **kwargs): # + def parse(self, *args, **kwargs): """ Read in all the information from the input_file. Alias to parser.parse. diff --git a/requirements.txt b/requirements.txt index d7d525dc..278394b2 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,6 +2,6 @@ numpy mpmath matplotlib scipy -ase +ase<3.17.0 graphviz tqdm