import catmap
from catmap import ReactionModelWrapper
from catmap.model import ReactionModel
from catmap.functions import smooth_piecewise_linear
from catmap.functions import parse_constraint
import pylab as plt
import numpy as np
try:
from scipy.optimize import fmin
from scipy import integrate
except ImportError:
fmin = None
integrate = None
from itertools import combinations_with_replacement as combinations
[docs]class FirstOrderInteractions(ReactionModelWrapper):
"""Class for implementing 'first-order adsorbate interaction model.
Should be sub-classed by scaler."""
[docs] def __init__(self,reaction_model=None):
if reaction_model is None:
reaction_model = ReactionModel()
self._rxm = reaction_model
defaults = dict(
cross_interaction_mode='geometric_mean',
transition_state_cross_interaction_mode='intermediate_state',
default_self_interaction_parameter = 0,
interaction_response_function = 'linear',
interaction_response_parameters = {'slope':1,'cutoff':0.25,
'smoothing':0.05},
interaction_fitting_mode=None,
input_overrides_fit = True, #user inputs override fitted values
interaction_strength = 1,
#weight interaction parameters by this
numerical_delta_theta = None
)
self._rxm.update(defaults)
self._required = {'cross_interaction_mode':str,
'transition_state_cross_interaction_mode':str,
'interaction_fitting_mode':None
}
self._log_strings = {'interaction_fitting_success':
"interaction parameter ${param}@${surface} = ${paramval} (avg_error=${error})"}
[docs] def parameterize_interactions(self):
self._parameterized = True
self.get_interaction_transition_state_scaling_matrix()
if self.interaction_fitting_mode is not None:
if '+' in self.interaction_fitting_mode:
self.interaction_fitting_mode = self.interaction_fitting_mode.split('+')
try:
self.fit()
except ValueError as inst:
#This is due to poor error handling elsewhere... if self.fit() raises
# a ValueError it looks like a convergence issue. While this hack
# is not ideal, it at least isolates the problem to within self.fit()
raise UserWarning('Parameter fitting failed with ValueError: '+inst.message)
[docs] def get_energy_error(self, epsilon_ij, theta, Ediff, Eint,
parameter_name, surface_name):
i_surf = self.surface_names.index(surface_name)
all_ads = self.adsorbate_names + self.transition_state_names
E0_list = [self.species_definitions[sp]['formation_energy'][i_surf]
for sp in all_ads]
if None in E0_list:
no_Es = [a for a,E in zip(all_ads,E0_list) if E is None]
raise UserWarning('No adsorption energy specified for: '+str(no_Es))
Ediff,Ediff_ads = Ediff
ads_a, ads_b = parameter_name.split('&')
if Ediff_ads == ads_a:
i = all_ads.index(ads_a)
j = all_ads.index(ads_b)
elif Ediff_ads == ads_b:
i = all_ads.index(ads_b)
j = all_ads.index(ads_a)
else:
raise UserWarning('Differential adsorbate not in parameter name')
if self.interaction_cross_term_names:
param_names = self.adsorbate_names + self.interaction_cross_term_names
else:
param_names = self.adsorbate_names
info = self.thermodynamics.adsorbate_interactions.get_interaction_info()
params = [info[pi][i_surf] for pi in param_names]
int_strength = self.interaction_strength
self.interaction_strength = 1.0
eps_matrix = self.params_to_matrix(E0_list+params)
self.interaction_strength = int_strength
eps_matrix[i,j] = eps_matrix[j,i] = epsilon_ij
eps_list = list(eps_matrix.ravel())
if 'differential' in self.interaction_fitting_mode:
include_derivative = True
else:
include_derivative = False
include_integral = False
if 'integral' in self.interaction_fitting_mode:
include_integral = True
if 'numerical_differential' in self.interaction_fitting_mode:
include_integral = True
Eint_m,Ediff_m,Ejac = self.interaction_function(theta,E0_list,eps_list,
self.interaction_response_function,include_derivative,
include_integral)
if 'numerical_differential' in self.interaction_fitting_mode:
if 'differential' in self.interaction_fitting_mode:
raise UserWarning('interaction_fitting_mode cannot contain both '
'numerical_differential and differential.')
if 'delta_theta' in self.species_definitions[all_ads[i]]:
delta_theta = self.species_definitions[all_ads[i]]['delta_theta']
elif self.numerical_delta_theta is None:
raise UserWarning('The change in coverage used to numerically compute '
'differential binding energies (numerical_delta_theta) '
'must be specified.')
else:
delta_theta = self.numerical_delta_theta
new_theta = [x for x in theta]
new_theta[i] -= delta_theta
Eint_i = self.interaction_function(new_theta,E0_list,eps_list,
self.interaction_response_function,False,True)[0]
Ediff_model = (Eint_m - Eint_i)/(delta_theta)
diff_err = (Ediff-Ediff_model)
elif 'differential' in self.interaction_fitting_mode:
Ediff_model = Ediff_m[i]
diff_err = (Ediff-Ediff_model)
else:
diff_err = None
if 'integral' in self.interaction_fitting_mode:
if Eint_m is None:
raise AttributeError('Interaction model has no function for ' +\
'computing integral interacting energies. Fit to differential '+\
'energies instead.')
int_err = (Eint-Eint_m)
else:
int_err = None
return diff_err, int_err
[docs] def error_norm(self,diff_err, int_err):
if diff_err is None:
return abs(int_err)
elif int_err is None:
return abs(diff_err)
elif self.integral_error_fitting_weight:
w = self.integral_error_fitting_weight
return abs(diff_err) + w*abs(int_err)
else:
return abs(diff_err) + abs(int_err)
[docs] def fit_interaction_parameter(self,theta_list,E_diffs,E_ints,param_name,surf_name):
def target_function(eps_ij, theta_list, E_diffs,E_ints,param_name,surf_name):
error = 0
for theta,Ed,Ei in zip(theta_list,E_diffs,E_ints):
d_err, i_err = self.get_energy_error(eps_ij, theta,Ed,Ei,
param_name,surf_name)
err_norm = self.error_norm(d_err, i_err)
error += err_norm
return error/len(theta_list)
minimize = lambda x: target_function(x,theta_list,E_diffs,E_ints,
param_name,surf_name)
x0 = getattr(self,'cross_interaction_initial_guess',0)
if x0 is None:
x0 = 0
verbose_fitting = getattr(self,'verbose_interaction_parameter_fitting',False)
eps_ij = fmin(minimize,[x0],disp=verbose_fitting)[0]
err = minimize(eps_ij)
self.log('interaction_fitting_success',
param = param_name,
paramval = str(eps_ij),
error = str(round(err,2)),
surface = surf_name
)
return eps_ij
[docs] def fit(self):
def required_params_to_fitting_info(param, key, required_params, fitting_info):
#helper function to transfer data from "required_params_dict" to "fitting_info"
theta_list,Ediffs,Eints = zip(*required_params_dict[key])
eps_ij = self.fit_interaction_parameter(theta_list,Ediffs,Eints,param,surf)
fitting_info[param] = [required_params_dict[key],eps_ij]
del required_params_dict[key]
def sync_with_species_defs(fitting_info,surf_id):
#helper function to dump information into species_definitions
for key in [k for k in fitting_info.keys() if '&' in k]:
eps_ij = fitting_info[key][-1]
ads_a,ads_b = sorted(key.split('&'))
if ads_a == ads_b:
old_param = self.species_definitions[ads_a].get('self_interaction_parameter',None)
if (old_param is None):
if old_param is None:
old_param = [None]*len(self.surface_names)
old_param[surf_id] = eps_ij
self.species_definitions[ads_a]['self_interaction_parameter'] = old_param
else:
cross_ints_a = self.species_definitions[ads_a].get('cross_interaction_parameters',None)
cross_ints_b = self.species_definitions[ads_b].get('cross_interaction_parameters',None)
if cross_ints_a == cross_ints_b == None:
cross_ints = {}
cross_key = ads_a
elif cross_ints_b is None:
cross_ints = cross_ints_a
cross_key = ads_a
elif cross_ints_a is None:
cross_ints = cross_ints_b
cross_key = ads_b
else:
if ads_a in cross_ints_b:
cross_ints = cross_ints_b
cross_key = ads_b
else:
cross_ints = cross_ints_a
cross_key = ads_a
if cross_key == ads_a:
other_key = ads_b
else:
other_key = ads_a
none_vec = [None]*len(self.surface_names)
if cross_key == ads_a:
cross_ints_vec = cross_ints.get(ads_b,none_vec)
elif cross_key == ads_b:
cross_ints_vec = cross_ints.get(ads_a,none_vec)
cross_ints_vec[surf_id] = eps_ij
cross_ints[other_key] = cross_ints_vec
self.species_definitions[cross_key]['cross_interaction_parameters'] = cross_ints
all_ads = self.adsorbate_names + self.transition_state_names
all_fitting_info = {}
for surf_id,surf in enumerate(self.surface_names):
#collect all coverage/energy pairs
all_cvgs = []
for key in self.species_definitions:
info = self.species_definitions[key]
if info.get('coverage_dependent_energy',[None]*(surf_id+1))[surf_id]:
for cvg_i in info['coverage_dependent_energy'][surf_id]:
cvg,Ediff,Eint= cvg_i
Ediff = [Ediff,key]
all_cvgs.append([cvg,Ediff,Eint])
required_params_dict = {}
fitting_info = {}
#organize coverages by which parameters the correspond to
for cvg_Ed_Eint_i in all_cvgs:
cvg_i,Ed_i,Eint_i = cvg_Ed_Eint_i
reqd = str(self.required_interaction_parameters(cvg_i))
if reqd not in required_params_dict:
required_params_dict[reqd] = [cvg_Ed_Eint_i]
else:
if cvg_Ed_Eint_i not in required_params_dict[reqd]:
required_params_dict[reqd] += [cvg_Ed_Eint_i]
for key in list(required_params_dict.keys()):
print(key)
keylist = eval(key)
cvgs = required_params_dict[key]
if len(keylist) == 1:
ads = keylist[0].split('&')[0]
#pull out low-coverage limit calcs since they are already used
if len(cvgs) == 1:
fitting_info[ads] = [cvgs,cvgs[0][1][0]]
del required_params_dict[key]
else:
for cvg in cvgs:
theta,Ed,Ei = cvg
Ed,Ed_ads = Ed
Ef = self.species_definitions[Ed_ads]['formation_energy'][surf_id]
if Ed == Ef:
fitting_info[ads] = [[cvg],Ed]
cvgs.remove(cvg)
self_keys = [k for k in required_params_dict if len(eval(k)) == 1]
cross_keys_3 = [k for k in required_params_dict if len(eval(k)) == 3]
cross_keys_6 = [k for k in required_params_dict if len(eval(k)) == 6]
user_inputs = []
#get user-input values
for ads in all_ads:
param = '&'.join([ads,ads])
default_vec = [self.default_self_interaction_parameter]*len(self.surface_names)
if param not in fitting_info:
val = self.species_definitions[ads].get('self_interaction_parameter',
None)
if val is not None:
user_inputs.append(param)
val = val[surf_id]
else:
val = default_vec[surf_id]
#ensure that all self-interactions have a value
fitting_info[param] = [[],val]
#get user-input cross-interactions
for ads in all_ads:
cross_ints = self.species_definitions[ads].get('cross_interaction_parameters',None)
if cross_ints:
for key in cross_ints:
param1 = '&'.join([ads,key])
param2 = '&'.join([key,ads])
val = cross_ints[key][surf_id]
if val is not None:
user_inputs += [param1,param2]
if not self.input_overrides_fit:
user_inputs = []
#fit known cross-interactions
for key in self_keys:
param = eval(key)[0]
if param not in user_inputs:
required_params_to_fitting_info(param,key,required_params_dict,fitting_info)
sync_with_species_defs(fitting_info,surf_id)
for key in cross_keys_3:
params = eval(key)
params = [pi for pi in params if pi not in list(fitting_info.keys())+user_inputs]
assert len(params) <= 1 #if not then something is very weird.
if params:
param = params[0]
if param not in user_inputs:
required_params_to_fitting_info(param,key,required_params_dict,fitting_info)
else:
print('Using user-input value for '+eval(key)[1])
sync_with_species_defs(fitting_info,surf_id)
underdetermined = cross_keys_6
new_underdetermined = []
initialized = False
while underdetermined != new_underdetermined:
#hacky way of back-substituting... this can be cleaned up.
if not initialized:
new_underdetermined = [k for k in underdetermined]
initialized = True
underdetermined = [k for k in new_underdetermined]
for key in underdetermined:
params = eval(key)
params = [pi for pi in params if pi not in fitting_info.keys()+user_inputs]
if len(params) == 1: #then we can do fitting
param = params[0]
if param not in user_inputs:
required_params_to_fitting_info(param,key,required_params_dict,fitting_info)
new_underdetermined.remove(key)
elif len(params) == 0:
cvg_states = []
for cvg,Ed,Ei in required_params_dict[key]:
cvg_state = ','.join(['='.join([a,str(c)]) for a,c in zip(all_ads,cvg) if c])
cvg_state += ' (E_'+Ed[-1]+'='+str(Ed[0])+',E_int='+str(Ei)+')'
cvg_states.append(cvg_state)
print('Warning: Not fitting to coverage-depedent energies for :'+'\n'.join(cvg_states)+'. '
'Parameters are already fully determined.')
new_underdetermined.remove(key)
else:
pass
sync_with_species_defs(fitting_info,surf_id)
for key in underdetermined:
missing = [k for k in eval(key) if k not in fitting_info]
cvg_states = []
for cvg,Ed,Ei in required_params_dict[key]:
cvg_state = ','.join(['='.join([a,str(c)]) for a,c in zip(all_ads,cvg) if c])
cvg_state += ' (E_'+Ed[-1]+'='+str(Ed[0])+',E_int='+str(Ei)+')'
cvg_states.append(cvg_state)
state_str = '\n'.join(cvg_states)
missing_str = ' and/or '.join(missing)
print('Warning: Not fitting to coverage-depedent energies for:'+state_str+'. '
'Parameters are under-determined; add calculations to determine ' + \
missing_str +' to make this data useful.')
all_fitting_info[surf] = fitting_info
self.get_interaction_info()
self.fitting_info = all_fitting_info #store as attribute so that it can be double-checked externally
[docs] def required_interaction_parameters(self,cvg):
nonzeros = [i for i,c in enumerate(cvg) if c]
all_combos = []
for combo in combinations(nonzeros,2):
all_combos.append('&'.join([self.parameter_names[i] for i in combo]))
return all_combos
[docs] def get_interaction_info(self):
interaction_dict = {}
n_ads = len(self.adsorbate_names)
cross_term_names = []
for a in self.adsorbate_names:
interaction_dict[a] = self.species_definitions[a].get('self_interaction_parameter',None)
if not interaction_dict[a]:
eps_ii0 = self.default_self_interaction_parameter
if a not in getattr(self,'_self_interaction_warned',[]):
print('Warning: No self-interaction parameter specified for '+ \
a+'. Assuming default ('+str(eps_ii0)+')')
self._self_interaction_warned = getattr(self,
'_self_interaction_warned',[]) + [a]
interaction_dict[a] = [eps_ii0]*len(self.surface_names)
cross_params = self.species_definitions[a].get(
'cross_interaction_parameters',{})
for cp in cross_params:
if cp not in self.adsorbate_names+self.transition_state_names:
raise ValueError(
'Cross parameter name must be in adsorbate names. The '+\
'name ' + cp + ' is not in ' +\
str(self.adsorbate_names+self.transition_state_names))
params = cross_params[cp]
if len(params) != len(self.surface_names):
raise ValueError('Cross parameters must be specified for each surface. '+\
str(len(params)) + ' parameters were specified, but there are '+\
str(len(self.surface_names)) + 'surfaces for '+ a +','+cp+'.' )
ads_a, ads_b = sorted([a,cp]) #sort to avoid duplicates
name = '&'.join([ads_a,ads_b])
if name not in cross_term_names:
cross_term_names.append(name)
interaction_dict[name] = params
if cross_term_names:
cross_term_names = tuple(cross_term_names)
self.interaction_cross_term_names = cross_term_names
return interaction_dict
[docs] def get_interaction_scaling_matrix(self):
interaction_dict = self.get_interaction_info()
cross_names = self.interaction_cross_term_names
if cross_names:
param_names = self.adsorbate_names + cross_names
else:
param_names = self.adsorbate_names
constraint_dict = {}
if not self.interaction_scaling_constraint_dict:
self.interaction_scaling_constraint_dict = self.scaling_constraint_dict.copy()
for ads in self.interaction_scaling_constraint_dict:
if '-' not in ads:
constr = self.interaction_scaling_constraint_dict[ads]
new_constr = []
#preserve only 0 constraints
for ci in constr:
if ci != 0:
new_constr.append(None)
else:
new_constr.append(0)
constraint_dict[ads] = new_constr
for ads in self.interaction_scaling_constraint_dict:
if '&' in ads:
a,b = ads.split('&')
a,b = sorted([a,b])
new_ads = '&'.join([a,b])
constraint_dict[new_ads] = self.interaction_scaling_constraint_dict[ads]
else:
constraint_dict[ads] = self.interaction_scaling_constraint_dict[ads]
#get mins/maxs
interaction_mins = []
interaction_maxs = []
if self.default_interaction_constraints is None:
self.default_interaction_constraints = [None]*(len(self.descriptor_names)+1)
for p in param_names:
if p not in constraint_dict:
constr = self.default_interaction_constraints
else:
constr = constraint_dict[p]
minvs,maxvs = parse_constraint(constr,p)
interaction_mins.append(minvs)
interaction_maxs.append(maxvs)
#get coefficient matrix
C = catmap.functions.scaling_coefficient_matrix(
interaction_dict, self.descriptor_dict,
self.surface_names,
param_names,
interaction_mins,
interaction_maxs)
self.interaction_coefficient_matrix = C.T
return self.interaction_coefficient_matrix
[docs] def get_interaction_transition_state_scaling_matrix(self):
#get TS scaling matrix equivalent
if self.transition_state_cross_interaction_mode == 'transition_state_scaling':
if self.transition_state_scaling_matrix is None:
raise AttributeError('Transition state scaling can only be used '+\
' for interactions if the transition_state_scaling_matrix is defined.')
TS_weight_matrix = []
for params in list(self.transition_state_scaling_matrix):
TS_weight_matrix.append(params[:-1])
else:
#take TS interaction parameters as some linear combination of initial/final state
if self.transition_state_cross_interaction_mode == 'initial_state':
weight = 0
elif self.transition_state_cross_interaction_mode == 'final_state':
weight = 1
elif self.transition_state_cross_interaction_mode.startswith('intermediate_state'):
if self.transition_state_cross_interaction_mode == 'intermediate_state':
weight = 0.5
else:
mname = self.transition_state_cross_interaction_mode
junk,w = mname.rsplit('(',1)
w,junk = w.rsplit(')',1)
try:
weight = float(w)
except:
raise AttributeError('Intermediate state weight is not floatable:'+w)
else:
raise AttributeError('Undefined transition_state_cross_interaction_mode: ' +\
self.transition_state_cross_interaction_mode)
TS_weight_matrix = self.get_TS_weight_matrix(weight)
self.interaction_transition_state_scaling_matrix = TS_weight_matrix
[docs] def get_interaction_matrix(self,descriptors):
full_descriptors = list(descriptors) + [1.]
param_vector = np.dot(self.coefficient_matrix,full_descriptors)
eps_vector = self.params_to_matrix(param_vector)
return eps_vector
[docs] def params_to_matrix(self,param_vector):
n_ads = len(self.adsorbate_names)
n_TS = len(self.transition_state_names)
all_names = self.adsorbate_names + self.transition_state_names
n_tot = len(all_names)
epsilon_matrix = np.zeros((n_tot,n_tot))
self_interactions = param_vector[n_tot:n_tot+n_ads]
cross_interactions = param_vector[n_tot+n_ads:]
for ads,param in zip(self.adsorbate_names,self_interactions):
max_val = self.species_definitions[ads].get('max_self_interaction',None)
if max_val is not None:
param = min(param,max_val)
self_interactions[self.adsorbate_names.index(ads)] = param
for i,e_ii in enumerate(self_interactions):
epsilon_matrix[i,i] = e_ii
for i, e_i in enumerate(self_interactions):
for j, e_j in enumerate(self_interactions):
if not epsilon_matrix[i,j]:
if self.cross_interaction_mode == 'geometric_mean':
epsilon_matrix[i,j] = np.sqrt(abs(e_i)*abs(e_j))
elif self.cross_interaction_mode == 'arithmetic_mean':
epsilon_matrix[i,j] = (e_i+e_j)/2.
elif self.cross_interaction_mode == 'neglect':
epsilon_matrix[i,j] = 0
if self.non_interacting_site_pairs:
for site_1, site_2 in self.non_interacting_site_pairs:
idxs1 = [self.adsorbate_names.index(sp) for
sp in self.adsorbate_names if
self.species_definitions[sp]['site']==site_1]
idxs2 = [self.adsorbate_names.index(sp) for
sp in self.adsorbate_names if
self.species_definitions[sp]['site']==site_2]
for i in idxs1:
for j in idxs2:
epsilon_matrix[i,j] = epsilon_matrix[j,i] = 0
for i,TS_params in enumerate(list(self.interaction_transition_state_scaling_matrix)):
i += n_ads
TS_params = list(TS_params)
for j, epsilon_params in enumerate(list(epsilon_matrix[0:n_ads,0:n_ads])):
e_TS = np.dot(TS_params,epsilon_params)
epsilon_matrix[i,j] = epsilon_matrix[j,i] = e_TS
if len(cross_interactions)>0:
for name, param in zip(self.interaction_cross_term_names,cross_interactions):
a,b = name.split('&')
i = all_names.index(a)
j = all_names.index(b)
epsilon_matrix[i,j] = epsilon_matrix[j,i] = param
epsilon_matrix *= self.interaction_strength
self._interaction_matrix = epsilon_matrix.tolist()
return epsilon_matrix
[docs] def get_TS_weight_matrix(self,weight):
"""Helper function to get `weights' of how
to distribute TS-cross interactions between IS/FS.
Should not be called externally."""
weight_matrix = []
for TS in self.transition_state_names:
weight_TS = [0]*len(self.adsorbate_names)
IS = None
FS = None
for rxn in self.elementary_rxns:
if TS in rxn[1]:
if IS or FS:
print('Warning: Ambiguous initial/final state for '+TS)
IS = rxn[0]
FS = rxn[-1]
for ads in IS:
if ads in self.adsorbate_names:
idx = self.adsorbate_names.index(ads)
weight_TS[idx] += (1-weight)
for ads in FS:
if ads in self.adsorbate_names:
idx = self.adsorbate_names.index(ads)
weight_TS[idx] += weight
weight_matrix.append(weight_TS)
return weight_matrix
[docs] @staticmethod
def smooth_piecewise_linear_response(*args,**kwargs):
return smooth_piecewise_linear(*args,**kwargs)[:2]
[docs] @staticmethod
def piecewise_linear_response(*args,**kwargs):
kwargs['smoothing'] = 0
return smooth_piecewise_linear(*args,**kwargs)[:2]
[docs] @staticmethod
def linear_response(*args,**kwargs):
kwargs['smoothing'] = 0
kwargs['cutoff'] = 0
return smooth_piecewise_linear(*args,**kwargs)[:2]
[docs] def fit_old(self): ##REMOVE AFTER WORKING
all_adsnames = self.adsorbate_names+self.transition_state_names
E0_list = []
for i, surf in enumerate(self.surface_names):
fitting_info = {}
for sp in self.species_definitions:
cvg_dep = self.species_definitions[sp].get('coverage_dependent_energy',
[None]*len(self.surface_names))[i]
if cvg_dep:
for theta_E in cvg_dep:
theta,Ed,Ei = theta_E
n_theta = sum(ti != 0 for ti in theta)
ads_idx = all_adsnames.index(sp)
if n_theta == 1:
if sp in fitting_info:
if theta_E not in fitting_info[sp]:
fitting_info[sp].append(theta_E)
else:
fitting_info[sp] = [theta_E]
elif n_theta == 2:
coads_idx = [j for j,cvg in enumerate(theta) if cvg and j != ads_idx][0]
coads = all_adsnames[coads_idx]
ads_a, ads_b = sorted([sp,coads])
name = '&'.join([ads_a,ads_b])
if name in fitting_info:
if theta_E not in fitting_info[name]:
fitting_info[name].append(theta_E)
else:
fitting_info[name] = [theta_E]
else:
print('Warning: Ignoring coverage dependent entry for '+sp+\
'. Only 1-2 adsorbates can have non-zero coverages')
for key in fitting_info.keys():
#clean up fitting_info
if '&' not in key and len(fitting_info[key]) == 1:
#move into cross-parameter info if possible
#This may not be totally general, but the goal
#is to ensure that cross-parameters are
#also fit to the low-coverage limit
cross_keys = [k for k in fitting_info.keys() if (
'&' in k and key in k.split('&'))]
for k in cross_keys:
fitting_info[k] = fitting_info[key] + fitting_info[k]
del fitting_info[key] #a single entry is the same one used for Ef
elif self.input_overrides_fit:
#user inputs should over-ride fit, therefore delete any fitting
#data for parameters the user already input
if '&' not in key:
if self.species_definitions[key].get(
'self_interaction_parameter',
[None]*len(self.surface_names))[i] is not None:
del fitting_info[key]
else:
ads_a, ads_b = key.split('&')
#have to check for cross interaction specified either way
ab_crossint = self.species_definitions[ads_a].get(
'cross_interaction_parameters',{})
if ads_b in ab_crossint:
del fitting_info[key]
ba_crossint = self.species_definitions[ads_b].get(
'cross_interaction_parameters',{})
if ads_a in ba_crossint:
del fitting_info[key]
self_params = [k for k in fitting_info if '&' not in k]
cross_params = [k for k in fitting_info if '&' in k]
E0_list = [self.species_definitions[sp]['formation_energy'][i] \
for sp in all_adsnames]
for key in self_params:
thetas,Ediffs,Eints = zip(*fitting_info[key])
eps_ii = self.fit_interaction_parameter(thetas,Ediffs,Eints,key,surf)
eps = self.species_definitions[key].get('self_interaction_parameter',
[None]*len(self.surface_names))
eps[i] = eps_ii
self.species_definitions[key]['self_interaction_parameter'] = eps
for key in cross_params:
thetas,Ediffs,Eints = zip(*fitting_info[key])
ads_a, ads_b = key.split('&')
idx_a = all_adsnames.index(ads_a)
idx_b = all_adsnames.index(ads_b)
thetas_a = [ti[idx_a] for ti in thetas]
thetas_b = [ti[idx_b] for ti in thetas]
stdev_a = np.std(thetas_a)
stdev_b = np.std(thetas_b)
if 0 not in [stdev_a,stdev_b]:
if 'differential' in self.interaction_fitting_mode:
print('Warning: no constant coverage for fitting '+\
'cross-interaction parameter '+key+'. Using '+\
'adsorbate with minimum standard deviation '+\
'of coverages as differential adsorbate')
if stdev_b < stdev_a:
#switch adsorbate so that the one at constant coverage
#is the adsorbate used for differential adsorption energies.
ads_a, ads_b = ads_b, ads_a
fit_param = '&'.join([ads_a,ads_b])
eps_ij = self.fit_interaction_parameter(thetas,Ediffs,Eints,
fit_param,surf)
eps_dict_b = self.species_definitions[ads_b].get(
'cross_interaction_parameters',
{})
if ads_a not in eps_dict_b:
eps_dict_b[ads_a] = [None]*len(self.surface_names)
eps_dict_b[ads_a][i] = eps_ij
self.species_definitions[ads_b][
'cross_interaction_parameters'] = eps_dict_b
#make sure that new parameters get incorporated into interaction matrix
self.get_interaction_info()
self._fitting_info = fitting_info