529 lines
22 KiB
Python
529 lines
22 KiB
Python
|
import numpy as np
|
||
|
|
import scipy.integrate as integrate
|
||
|
|
from scipy.optimize import minimize
|
||
|
|
import numbers
|
||
|
|
import typing
|
||
|
|
from tqdm import tqdm
|
||
|
|
from echem.core.useful_funcs import nearest_array_index, ClassMethods
|
||
|
|
E_F_SHE_VAC = -4.5 # Fermi Energy of Standard Hydrogen Electrode with respect to vacuum
|
||
|
|
|
||
|
|
|
||
|
|
class GM(ClassMethods):
|
||
|
|
"""This class calculates the final Fermi and Redox species distributions according
|
||
|
|
to the Gerischer-Marcus formalism.
|
||
|
|
|
||
|
|
Parameters:
|
||
|
|
-----------
|
||
|
|
DOS: np.ndarray, optional
|
||
|
|
The values of DOS in 1D numpy array. If not specified values will be taken from saved data.
|
||
|
|
|
||
|
|
E: np.ndarray, optional
|
||
|
|
The corresponding to the DOS energy mesh. If not specified values will be taken from saved data.
|
||
|
|
|
||
|
|
efermi: np.ndarray. optional
|
||
|
|
System Fermi level. If not specified values will be taken from saved data.
|
||
|
|
|
||
|
|
vacuum_lvl: np.ndarray, optional
|
||
|
|
System vacuum level. If not specified values will be taken from saved data.
|
||
|
|
"""
|
||
|
|
def __init__(self, path_to_data='Saved_data', DOS=None, E=None, efermi=None, vacuum_lvl=None):
|
||
|
|
# variables that might be defined through __init__ function
|
||
|
|
self.E = E
|
||
|
|
self.DOS = DOS
|
||
|
|
self.efermi = efermi
|
||
|
|
self.vacuum_lvl = vacuum_lvl
|
||
|
|
|
||
|
|
# variables that should be defined through set_params function
|
||
|
|
self.C_EDL = None
|
||
|
|
self.T = None
|
||
|
|
self.l = None
|
||
|
|
self.sheet_area = None
|
||
|
|
|
||
|
|
# variables that will be created during calculations
|
||
|
|
self.sigma_Q_arr = None
|
||
|
|
|
||
|
|
# variable that define numerical parameters of quantum charge calculation
|
||
|
|
self.__SIGMA_0 = 0.5
|
||
|
|
self.__SIGMA_ACCURACY = 1e-3
|
||
|
|
self.__SIGMA_RANGE = 4
|
||
|
|
|
||
|
|
if DOS is None:
|
||
|
|
try:
|
||
|
|
self.DOS = np.load(path_to_data + '/DOS.npy')
|
||
|
|
except OSError:
|
||
|
|
print('File DOS.npy does not exist')
|
||
|
|
|
||
|
|
if E is None:
|
||
|
|
try:
|
||
|
|
self.E = np.load(path_to_data + '/E.npy')
|
||
|
|
except OSError:
|
||
|
|
print('File E_DOS.npy does not exist')
|
||
|
|
|
||
|
|
if efermi is None:
|
||
|
|
try:
|
||
|
|
self.efermi = np.load(path_to_data + '/efermi.npy')
|
||
|
|
except OSError:
|
||
|
|
print('File efermi.npy does not exist')
|
||
|
|
|
||
|
|
if vacuum_lvl is None:
|
||
|
|
try:
|
||
|
|
self.vacuum_lvl = np.load(path_to_data + '/vacuum_lvl.npy')
|
||
|
|
except OSError:
|
||
|
|
print('File vacuum_lvl.npy does not exist')
|
||
|
|
|
||
|
|
def set_params(self, C_EDL, T, l, sheet_area):
|
||
|
|
"""Sets parameters of calculation
|
||
|
|
|
||
|
|
Parameters:
|
||
|
|
----------
|
||
|
|
C_EDL: float, str
|
||
|
|
float: Capacitance of electric double layer (microF/cm^2)
|
||
|
|
str: 'Q' calculating in the Quantum Capacitance Dominating limit (C_Q << C_EDL)
|
||
|
|
str: 'Cl' calculating in the Classical limit (C_Q >> C_EDL)
|
||
|
|
|
||
|
|
T: int, float
|
||
|
|
Temperature. It is used in computing Fermi function and distribution function of redox system states
|
||
|
|
|
||
|
|
l: float
|
||
|
|
Reorganization energy in eV
|
||
|
|
"""
|
||
|
|
self.C_EDL = C_EDL
|
||
|
|
self.T = T
|
||
|
|
self.l = l
|
||
|
|
|
||
|
|
self.sheet_area = sheet_area
|
||
|
|
|
||
|
|
def set_params_advance(self, SIGMA_0=0.5, ACCURACY_SIGMA=1e-3, SIGMA_RANGE=4):
|
||
|
|
"""
|
||
|
|
Sets numerical parameters that are used in quantum charge density calculations. Delete cashed
|
||
|
|
results of charge calculations.
|
||
|
|
Args:
|
||
|
|
SIGMA_0: float, optional
|
||
|
|
Initial guess for charge at equilibrium
|
||
|
|
ACCURACY_SIGMA: float, optional
|
||
|
|
Accuracy of charge calculation
|
||
|
|
SIGMA_RANGE: float, optional
|
||
|
|
It defines the minimum and maximum calculated charge
|
||
|
|
"""
|
||
|
|
self.__SIGMA_0 = SIGMA_0
|
||
|
|
self.__SIGMA_ACCURACY = ACCURACY_SIGMA
|
||
|
|
self.__SIGMA_RANGE = SIGMA_RANGE
|
||
|
|
self.sigma_Q_arr = None
|
||
|
|
|
||
|
|
@staticmethod
|
||
|
|
def fermi_func(E, T):
|
||
|
|
"""
|
||
|
|
Calculates Fermi-Dirac Distribution
|
||
|
|
Args:
|
||
|
|
E: Energies
|
||
|
|
T: Temperature in K
|
||
|
|
"""
|
||
|
|
k = 8.617e-5 # eV/K
|
||
|
|
return 1 / (1 + np.exp(E / (k * T)))
|
||
|
|
|
||
|
|
@staticmethod
|
||
|
|
def W_ox(E, T, l):
|
||
|
|
"""
|
||
|
|
Distribution of oxidized states
|
||
|
|
Args:
|
||
|
|
E (np.array): Energies
|
||
|
|
T (float): Temperature
|
||
|
|
l (float): Reorganization energy
|
||
|
|
"""
|
||
|
|
k = 8.617e-5 # eV/K
|
||
|
|
W_0 = (1 / np.sqrt(4 * k * T * l))
|
||
|
|
return W_0 * np.exp(- (E - l) ** 2 / (4 * k * T * l))
|
||
|
|
|
||
|
|
@staticmethod
|
||
|
|
def W_red(E, T, l):
|
||
|
|
"""
|
||
|
|
Distribution of reduced states
|
||
|
|
Args:
|
||
|
|
E (np.array): Energies
|
||
|
|
T (float): Temperature
|
||
|
|
l (float): Reorganization energy
|
||
|
|
"""
|
||
|
|
k = 8.617e-5 # eV/K
|
||
|
|
W_0 = (1 / np.sqrt(4 * k * T * l))
|
||
|
|
return W_0 * np.exp(- (E + l) ** 2 / (4 * k * T * l))
|
||
|
|
|
||
|
|
def compute_C_quantum(self, dE_Q_arr):
|
||
|
|
"""
|
||
|
|
Calculates differential quantum capacitance
|
||
|
|
Q = e * int{DOS(E) * [f(E) - f(E + deltaE)] dE}
|
||
|
|
C_Q = - dQ/d(deltaE) = - (e / (4*k*T)) * int{DOS(E) * sech^2[(E+deltaE)/(2*k*T)] dE}
|
||
|
|
Args:
|
||
|
|
dE_Q_arr (np.array, float): Energy shift at which C_Q is calculated
|
||
|
|
Returns:
|
||
|
|
Quantum capacitance in accordance with energy displacement(s)
|
||
|
|
TODO check constants
|
||
|
|
"""
|
||
|
|
self.check_existence('T')
|
||
|
|
self.check_existence('sheet_area')
|
||
|
|
|
||
|
|
k = 8.617e-5 # eV/K
|
||
|
|
|
||
|
|
elementary_charge = 1.6e-19 # C
|
||
|
|
k_1 = 1.38e-23 # J/K
|
||
|
|
const = (1e6 * elementary_charge ** 2) / (4 * k_1 * self.sheet_area) # micro F / cm^2
|
||
|
|
|
||
|
|
if isinstance(dE_Q_arr, typing.Iterable):
|
||
|
|
|
||
|
|
C_q_arr = np.zeros_like(dE_Q_arr)
|
||
|
|
|
||
|
|
for i, dE_Q in enumerate(dE_Q_arr):
|
||
|
|
E_2 = self.E - dE_Q # energy range for cosh function
|
||
|
|
cosh = np.cosh(E_2 / (2 * k * self.T))
|
||
|
|
integrand = (self.DOS / cosh) / cosh
|
||
|
|
C_q = (const / self.T) * integrate.simps(integrand, self.E)
|
||
|
|
C_q_arr[i] = C_q
|
||
|
|
|
||
|
|
return C_q_arr
|
||
|
|
|
||
|
|
def compute_C_total(self, E_diff_arr, add_info=False):
|
||
|
|
|
||
|
|
sigma_arr = np.zeros_like(E_diff_arr)
|
||
|
|
for i, E_diff in tqdm(enumerate(E_diff_arr), total=len(E_diff_arr)):
|
||
|
|
sigma_arr[i] = self.compute_sigma(E_diff, sigma_0=sigma_arr[i-1])
|
||
|
|
|
||
|
|
C_tot_arr = np.zeros_like(E_diff_arr)
|
||
|
|
C_Q_arr = np.zeros_like(E_diff_arr)
|
||
|
|
|
||
|
|
for i, (E_diff, sigma) in enumerate(zip(E_diff_arr, sigma_arr)):
|
||
|
|
ind = nearest_array_index(self.sigma_Q_arr, sigma)
|
||
|
|
E_step = self.__SIGMA_ACCURACY
|
||
|
|
E_start = - self.__SIGMA_RANGE
|
||
|
|
dE_Q = E_start + E_step * ind
|
||
|
|
C_Q = self.compute_C_quantum([dE_Q])
|
||
|
|
C_Q_arr[i] = C_Q[0]
|
||
|
|
C_tot = C_Q * self.C_EDL / (C_Q + self.C_EDL)
|
||
|
|
C_tot_arr[i] = C_tot
|
||
|
|
|
||
|
|
if add_info is False:
|
||
|
|
return C_tot_arr
|
||
|
|
else:
|
||
|
|
return C_tot_arr, C_Q_arr, sigma_arr
|
||
|
|
|
||
|
|
def compute_sigma_EDL(self, dE_EDL):
|
||
|
|
"""
|
||
|
|
Calculates charge corresponding to the potential drop of -dE_EDL/|e|.
|
||
|
|
Takes into account integral capacitance C_EDL
|
||
|
|
Args:
|
||
|
|
dE_EDL (float, np.array): Electron energy shift due to potential drop
|
||
|
|
Returns:
|
||
|
|
Charge or Sequence of charges
|
||
|
|
"""
|
||
|
|
self.check_existence('C_EDL')
|
||
|
|
return - self.C_EDL * dE_EDL
|
||
|
|
|
||
|
|
def compute_sigma_quantum(self, dE_Q_arr):
|
||
|
|
"""
|
||
|
|
Computes surface charge density induced by depletion or excess of electrons
|
||
|
|
|
||
|
|
Parameters:
|
||
|
|
----------
|
||
|
|
dE_Q_arr: np.ndarray, float
|
||
|
|
Shift in Fermi level due to quantum capacitance
|
||
|
|
|
||
|
|
Returns:
|
||
|
|
-------
|
||
|
|
sigmas: np.ndarray, float
|
||
|
|
Computed values (or one value) of surface charge densities
|
||
|
|
"""
|
||
|
|
|
||
|
|
self.check_existence('T')
|
||
|
|
self.check_existence('sheet_area')
|
||
|
|
|
||
|
|
elementary_charge = 1.6e-13 # micro coulomb
|
||
|
|
|
||
|
|
if isinstance(dE_Q_arr, typing.Iterable):
|
||
|
|
y_fermi = self.fermi_func(self.E, self.T)
|
||
|
|
|
||
|
|
sigmas = []
|
||
|
|
|
||
|
|
for dE_Q in dE_Q_arr:
|
||
|
|
E_2 = self.E - dE_Q # energy range for shifted Fermi_Dirac function
|
||
|
|
y_fermi_shifted = self.fermi_func(E_2, self.T)
|
||
|
|
integrand = self.DOS * (y_fermi - y_fermi_shifted)
|
||
|
|
sigma = (elementary_charge / self.sheet_area) * integrate.simps(integrand, self.E)
|
||
|
|
sigmas.append(sigma)
|
||
|
|
|
||
|
|
return sigmas
|
||
|
|
|
||
|
|
elif isinstance(dE_Q_arr, numbers.Real):
|
||
|
|
y_fermi = self.fermi_func(self.E, self.T)
|
||
|
|
|
||
|
|
E_2 = self.E - dE_Q_arr # energy range for shifted Fermi_Dirac function
|
||
|
|
y_fermi_shifted = self.fermi_func(E_2, self.T)
|
||
|
|
integrand = self.DOS * (y_fermi - y_fermi_shifted)
|
||
|
|
sigma = (elementary_charge / self.sheet_area) * integrate.simps(integrand, self.E)
|
||
|
|
|
||
|
|
return sigma
|
||
|
|
else:
|
||
|
|
raise TypeError(f'Invalid type of dE_Q_arr: {type(dE_Q_arr)}')
|
||
|
|
|
||
|
|
def compute_sigma(self, E_diff, sigma_0=None):
|
||
|
|
|
||
|
|
def error_E_diff(sigma, E_diff, sigma_Q_arr):
|
||
|
|
ind = nearest_array_index(sigma_Q_arr, sigma)
|
||
|
|
dE_Q = E_start + E_step * ind
|
||
|
|
dE_EDL = - sigma / self.C_EDL
|
||
|
|
dE_total = dE_Q + dE_EDL
|
||
|
|
|
||
|
|
return (dE_total - E_diff) ** 2
|
||
|
|
|
||
|
|
for var in ['T', 'l', 'C_EDL']:
|
||
|
|
self.check_existence(var)
|
||
|
|
|
||
|
|
E_step = self.__SIGMA_ACCURACY
|
||
|
|
E_start = - self.__SIGMA_RANGE
|
||
|
|
if sigma_0 is None:
|
||
|
|
sigma_0 = self.__SIGMA_0
|
||
|
|
# check if we've already calculated sigma_Q_arr in another run
|
||
|
|
if self.sigma_Q_arr is None:
|
||
|
|
E_range = np.arange(E_start, -E_start, E_step)
|
||
|
|
sigma_Q_arr = self.compute_sigma_quantum(E_range)
|
||
|
|
self.sigma_Q_arr = sigma_Q_arr
|
||
|
|
else:
|
||
|
|
sigma_Q_arr = self.sigma_Q_arr
|
||
|
|
|
||
|
|
result = minimize(error_E_diff, np.array([sigma_0]), args=(E_diff, sigma_Q_arr))
|
||
|
|
sigma = result.x[0]
|
||
|
|
|
||
|
|
return sigma
|
||
|
|
|
||
|
|
def compute_distributions(self, V_std, overpot=0, reverse=False, add_info=False):
|
||
|
|
"""Computes Fermi-Dirac and Redox species distributions according to Gerischer-Markus formalism
|
||
|
|
with Quantum Capacitance
|
||
|
|
|
||
|
|
Parameters:
|
||
|
|
----------
|
||
|
|
V_std: float
|
||
|
|
Standard potential (vs SHE) of a redox couple (Volts)
|
||
|
|
overpot: float, optional
|
||
|
|
Overpotential (Volts). It shifts the electrode Fermi energy to -|e|*overpot
|
||
|
|
reverse: bool, optional
|
||
|
|
If reverse is False the process of electron transfer from electrode to the oxidized state of the
|
||
|
|
redox species is considered and vice versa
|
||
|
|
add_info: bool, optional
|
||
|
|
If False the func returns Fermi-Dirac and Redox species distributions
|
||
|
|
If True additionally returns dE_Q (Fermi energy shift due to the quantum capacitance),
|
||
|
|
sigma (surface charge) and E_diff (the whole energy shift with respect to the original Fermi level)
|
||
|
|
|
||
|
|
Returns:
|
||
|
|
-------
|
||
|
|
y_fermi: np.array
|
||
|
|
Fermi-Dirac distribution
|
||
|
|
y_redox: np.array
|
||
|
|
Redox species distributions
|
||
|
|
dE_Q: np.array, optional (if add_info == True)
|
||
|
|
Total shift of the Fermi energy due to the Quantum Capacitance
|
||
|
|
sigma: np.array, optional (if add_info == True)
|
||
|
|
surface charge in microF/cm^2
|
||
|
|
E_F_redox: np.array, optional (if add_info == True)
|
||
|
|
The sum of two energy displacement of the electrode due to the difference in Fermi level of Redox couple
|
||
|
|
and the electrode and overpotential. It splits into dE_Q and dE_EDL
|
||
|
|
"""
|
||
|
|
|
||
|
|
E_F_redox = E_F_SHE_VAC - self.efermi - V_std + self.vacuum_lvl - overpot
|
||
|
|
sigma = self.compute_sigma(E_F_redox)
|
||
|
|
|
||
|
|
ind = nearest_array_index(self.sigma_Q_arr, sigma)
|
||
|
|
E_step = self.__SIGMA_ACCURACY
|
||
|
|
E_start = - self.__SIGMA_RANGE
|
||
|
|
dE_Q = E_start + E_step * ind
|
||
|
|
|
||
|
|
E_fermi = self.E - dE_Q
|
||
|
|
E_DOS_redox = self.E - dE_Q - overpot
|
||
|
|
|
||
|
|
if reverse:
|
||
|
|
y_fermi = 1 - self.fermi_func(E_fermi, self.T)
|
||
|
|
y_redox = self.W_red(E_DOS_redox, self.T, self.l)
|
||
|
|
else:
|
||
|
|
y_fermi = self.fermi_func(E_fermi, self.T)
|
||
|
|
y_redox = self.W_ox(E_DOS_redox, self.T, self.l)
|
||
|
|
|
||
|
|
if not add_info:
|
||
|
|
return y_fermi, y_redox
|
||
|
|
else:
|
||
|
|
return y_fermi, y_redox, dE_Q, sigma, E_F_redox
|
||
|
|
|
||
|
|
def compute_k_HET(self, V_std_pot_arr, overpot_arr, reverse=False, add_info=False):
|
||
|
|
"""Computes integral k_HET using Gerischer-Markus formalism with quantum capacitance
|
||
|
|
|
||
|
|
Parameters:
|
||
|
|
----------
|
||
|
|
V_std_pot_arr: float, np.ndarray
|
||
|
|
A range of varying a standard potential
|
||
|
|
overpot_arr: float, np.ndarray
|
||
|
|
A range of varying an overpotential
|
||
|
|
reverse: bool, optional
|
||
|
|
if reverse is False the process of electron transfer from electrode to the oxidized state of the
|
||
|
|
redox mediator is considered and vice versa
|
||
|
|
|
||
|
|
Returns:
|
||
|
|
-------
|
||
|
|
k_HET: np.array
|
||
|
|
Calculated heterogeneous electron transfer rate constant according to Gerischer-Marcus model with quantum
|
||
|
|
capacitance
|
||
|
|
dE_Q_arr: np.ndarray, optional (if add_info == True)
|
||
|
|
Total shift of the Fermi energy due to the Quantum Capacitance for all calculated redox potentials or
|
||
|
|
overpotentials
|
||
|
|
sigma_arr: np.ndarray, optional (if add_info == True)
|
||
|
|
surface charge in microF/cm^2 for all calculated redox potentials or overpotentials
|
||
|
|
E_F_redox_arr: np.ndarray, optional (if add_info == True)
|
||
|
|
The sum of two energy displacement of the electrode due to the difference in Fermi level of Redox couple
|
||
|
|
and the electrode and overpotential. It splits into dE_Q and dE_EDL. For all calculated redox potentials
|
||
|
|
or overpotentials
|
||
|
|
y_fermi_arr: 2D np.ndarray, optional (if add_info == True)
|
||
|
|
Fermi-Dirac distribution for all calculated redox potentials or overpotentials
|
||
|
|
y_redox_arr: 2D np.ndarray, optional (if add_info == True)
|
||
|
|
Redox species distributions for all calculated redox potentials or overpotentials
|
||
|
|
"""
|
||
|
|
|
||
|
|
if isinstance(self.C_EDL, numbers.Real):
|
||
|
|
if isinstance(V_std_pot_arr, typing.Iterable) and isinstance(overpot_arr, numbers.Real):
|
||
|
|
k_HET = np.zeros_like(V_std_pot_arr)
|
||
|
|
if not add_info:
|
||
|
|
for i, V_std in tqdm(enumerate(V_std_pot_arr), total=len(V_std_pot_arr)):
|
||
|
|
y_fermi, y_redox = self.compute_distributions(V_std, reverse=reverse, overpot=overpot_arr)
|
||
|
|
integrand = self.DOS * y_fermi * y_redox
|
||
|
|
k_HET[i] = integrate.simps(integrand, self.E) / self.sheet_area
|
||
|
|
return k_HET
|
||
|
|
else:
|
||
|
|
dE_Q_arr = np.zeros_like(V_std_pot_arr)
|
||
|
|
sigma_arr = np.zeros_like(V_std_pot_arr)
|
||
|
|
E_F_redox_arr = np.zeros_like(V_std_pot_arr)
|
||
|
|
y_fermi_arr = np.zeros((len(V_std_pot_arr), len(self.E)))
|
||
|
|
y_redox_arr = np.zeros((len(V_std_pot_arr), len(self.E)))
|
||
|
|
for i, V_std in tqdm(enumerate(V_std_pot_arr), total=len(V_std_pot_arr)):
|
||
|
|
y_fermi, y_redox, dE_Q, sigma, E_F_redox = self.compute_distributions(V_std, reverse=reverse,
|
||
|
|
overpot=overpot_arr,
|
||
|
|
add_info=add_info)
|
||
|
|
integrand = self.DOS * y_fermi * y_redox
|
||
|
|
k_HET[i] = integrate.simps(integrand, self.E) / self.sheet_area
|
||
|
|
dE_Q_arr[i] = dE_Q
|
||
|
|
sigma_arr[i] = sigma
|
||
|
|
E_F_redox_arr[i] = E_F_redox
|
||
|
|
y_fermi_arr[i] = y_fermi
|
||
|
|
y_redox_arr[i] = y_redox
|
||
|
|
return k_HET, dE_Q_arr, sigma_arr, E_F_redox_arr, y_fermi_arr, y_redox_arr
|
||
|
|
|
||
|
|
elif isinstance(overpot_arr, typing.Iterable) and isinstance(V_std_pot_arr, numbers.Real):
|
||
|
|
k_HET = np.zeros_like(overpot_arr)
|
||
|
|
if not add_info:
|
||
|
|
for i, overpot in tqdm(enumerate(overpot_arr), total=len(overpot_arr)):
|
||
|
|
y_fermi, y_redox = self.compute_distributions(V_std_pot_arr, reverse=reverse, overpot=overpot)
|
||
|
|
integrand = self.DOS * y_fermi * y_redox
|
||
|
|
k_HET[i] = integrate.simps(integrand, self.E) / self.sheet_area
|
||
|
|
|
||
|
|
return k_HET
|
||
|
|
else:
|
||
|
|
dE_Q_arr = np.zeros_like(overpot_arr)
|
||
|
|
sigma_arr = np.zeros_like(overpot_arr)
|
||
|
|
E_F_redox_arr = np.zeros_like(overpot_arr)
|
||
|
|
y_fermi_arr = np.zeros((len(overpot_arr), len(self.E)))
|
||
|
|
y_redox_arr = np.zeros((len(overpot_arr), len(self.E)))
|
||
|
|
for i, overpot in tqdm(enumerate(overpot_arr), total=len(overpot_arr)):
|
||
|
|
y_fermi, y_redox, dE_Q, sigma, E_F_redox = self.compute_distributions(V_std_pot_arr,
|
||
|
|
reverse=reverse,
|
||
|
|
overpot=overpot,
|
||
|
|
add_info=add_info)
|
||
|
|
integrand = self.DOS * y_fermi * y_redox
|
||
|
|
k_HET[i] = integrate.simps(integrand, self.E) / self.sheet_area
|
||
|
|
dE_Q_arr[i] = dE_Q
|
||
|
|
sigma_arr[i] = sigma
|
||
|
|
E_F_redox_arr[i] = E_F_redox
|
||
|
|
y_fermi_arr[i] = y_fermi
|
||
|
|
y_redox_arr[i] = y_redox
|
||
|
|
return k_HET, dE_Q_arr, sigma_arr, E_F_redox_arr, y_fermi_arr, y_redox_arr
|
||
|
|
|
||
|
|
else:
|
||
|
|
raise ValueError('One and only one type of V_std_pot_arr and overpot arr must be Sequence. The other \
|
||
|
|
must be a Real number')
|
||
|
|
|
||
|
|
elif self.C_EDL == 'Cl':
|
||
|
|
if isinstance(V_std_pot_arr, typing.Iterable) and isinstance(overpot_arr, numbers.Real):
|
||
|
|
E_fermi = self.E
|
||
|
|
E_DOS_redox = self.E - overpot_arr
|
||
|
|
|
||
|
|
if reverse:
|
||
|
|
y_fermi = 1 - self.fermi_func(E_fermi, self.T)
|
||
|
|
y_redox = self.W_red(E_DOS_redox, self.T, self.l)
|
||
|
|
else:
|
||
|
|
y_fermi = self.fermi_func(E_fermi, self.T)
|
||
|
|
y_redox = self.W_ox(E_DOS_redox, self.T, self.l)
|
||
|
|
|
||
|
|
integrand = self.DOS * y_fermi * y_redox
|
||
|
|
k_HET = np.ones_like(V_std_pot_arr) * integrate.simps(integrand, self.E)
|
||
|
|
|
||
|
|
return k_HET
|
||
|
|
|
||
|
|
elif isinstance(overpot_arr, typing.Sequence) and isinstance(V_std_pot_arr, numbers.Real):
|
||
|
|
k_HET = np.zeros_like(overpot_arr)
|
||
|
|
|
||
|
|
for i, overpot in tqdm(enumerate(overpot_arr), total=len(overpot_arr)):
|
||
|
|
E_fermi = self.E
|
||
|
|
E_DOS_redox = self.E - overpot
|
||
|
|
|
||
|
|
if reverse:
|
||
|
|
y_fermi = 1 - self.fermi_func(E_fermi, self.T)
|
||
|
|
y_redox = self.W_red(E_DOS_redox, self.T, self.l)
|
||
|
|
else:
|
||
|
|
y_fermi = self.fermi_func(E_fermi, self.T)
|
||
|
|
y_redox = self.W_ox(E_DOS_redox, self.T, self.l)
|
||
|
|
|
||
|
|
integrand = self.DOS * y_fermi * y_redox
|
||
|
|
k_HET[i] = integrate.simps(integrand, self.E)
|
||
|
|
|
||
|
|
return k_HET
|
||
|
|
|
||
|
|
else:
|
||
|
|
raise ValueError('One and only one type of V_std_pot_arr and overpot arr must be Sequence. The other \
|
||
|
|
must be Real number')
|
||
|
|
|
||
|
|
elif self.C_EDL == 'Q':
|
||
|
|
if isinstance(V_std_pot_arr, typing.Iterable) and isinstance(overpot_arr, numbers.Real):
|
||
|
|
k_HET = np.zeros_like(V_std_pot_arr)
|
||
|
|
|
||
|
|
for i, V_std in tqdm(enumerate(V_std_pot_arr), total=len(V_std_pot_arr)):
|
||
|
|
E_F_redox = E_F_SHE_VAC - self.efermi - V_std + self.vacuum_lvl
|
||
|
|
E_DOS_redox = self.E - E_F_redox
|
||
|
|
E_fermi = E_DOS_redox - overpot_arr
|
||
|
|
|
||
|
|
if reverse:
|
||
|
|
y_fermi = 1 - self.fermi_func(E_fermi, self.T)
|
||
|
|
y_redox = self.W_red(E_DOS_redox, self.T, self.l)
|
||
|
|
else:
|
||
|
|
y_fermi = self.fermi_func(E_fermi, self.T)
|
||
|
|
y_redox = self.W_ox(E_DOS_redox, self.T, self.l)
|
||
|
|
|
||
|
|
integrand = self.DOS * y_fermi * y_redox
|
||
|
|
k_HET[i] = integrate.simps(integrand, self.E)
|
||
|
|
|
||
|
|
return k_HET
|
||
|
|
|
||
|
|
elif isinstance(overpot_arr, typing.Iterable) and isinstance(V_std_pot_arr, numbers.Real):
|
||
|
|
k_HET = np.zeros_like(overpot_arr)
|
||
|
|
|
||
|
|
for i, overpot in tqdm(enumerate(overpot_arr), total=len(overpot_arr)):
|
||
|
|
E_F_redox = E_F_SHE_VAC - self.efermi - V_std_pot_arr + self.vacuum_lvl - overpot
|
||
|
|
E_fermi = self.E - E_F_redox
|
||
|
|
E_DOS_redox = self.E - E_F_redox - overpot
|
||
|
|
|
||
|
|
if reverse:
|
||
|
|
y_fermi = 1 - self.fermi_func(E_fermi, self.T)
|
||
|
|
y_redox = self.W_red(E_DOS_redox, self.T, self.l)
|
||
|
|
else:
|
||
|
|
y_fermi = self.fermi_func(E_fermi, self.T)
|
||
|
|
y_redox = self.W_ox(E_DOS_redox, self.T, self.l)
|
||
|
|
|
||
|
|
integrand = self.DOS * y_fermi * y_redox
|
||
|
|
k_HET[i] = integrate.simps(integrand, self.E)
|
||
|
|
|
||
|
|
return k_HET
|
||
|
|
|
||
|
|
else:
|
||
|
|
raise ValueError('One and only one type of V_std_pot_arr and overpot arr must be Sequence. The other \
|
||
|
|
must be Real number')
|