diff --git a/.coveragerc b/.coveragerc index c0d8a60..9d0de59 100644 --- a/.coveragerc +++ b/.coveragerc @@ -8,5 +8,5 @@ skip_empty = true exclude_also = raise NotImplementedError - ; don't complain abput abstract methods, they aren't run: - @(abc\.)?abstractmethod \ No newline at end of file + ; don't complain about abstract methods, they aren't run: + @(abc\.)?abstractmethod diff --git a/README.md b/README.md index 0220a40..ef0bcd1 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,7 @@ CVsim -------- -`cvsim` is a Python package for cyclic voltammogram (CV) simulation via a semi-analytical method developed by Oldham and Myland1. It is valid for CV experiments performed on disk macroelectrodes, and uses a semiintegration algorithm. In the limit of infinitely small potential steps, this algorithm is an exact solution. +`cvsim` is a Python package for cyclic voltammogram (CV) simulation via a semi-analytical method developed by Oldham and Myland1. It is valid for CV experiments performed on disk macroelectrodes, and uses a semi-integration algorithm. In the limit of infinitely small potential steps, this algorithm is an exact solution. Due to the precision of standard potentiostats, simulations that use a potential step of 1-5 mV typically provide a reasonable accuracy-computing time trade-off, where accuracy sanity checks (e.g. Randles-Sevcik relationship for Er and Eq mechanisms) have been performed. @@ -31,7 +31,7 @@ git clone git@github.com:ericfell/CVsim.git ## Package Structure -- `mechanisms.py`: Specifies one-/two-electron process mechanism to be simulated. +- `mechanisms.py`: Specifies one-/two-electron process mechanisms to be simulated. - `fit_curve.py`: Performs fitting of experimental CV data according to a desired mechanism. diff --git a/cv_fitting.py b/cv_fitting.py deleted file mode 100644 index 1e8d36c..0000000 --- a/cv_fitting.py +++ /dev/null @@ -1,85 +0,0 @@ - -import numpy as np -import matplotlib.pyplot as plt -from mechanisms import E_rev, E_q, E_qC, EE, SquareScheme -import scipy.constants as spc -from scipy.optimize import curve_fit - -F = spc.physical_constants['Faraday constant'][0] -R = spc.R - -# Part A: how to simulate schemes - -# example 1 electron schemes -potential1, current1 = E_q(0.3, -0.3, 0.0, 1.0, 1.0, 1e-5, 1e-5, 0.3, 1e-3).simulate() -potential2, current2 = E_qC(0.3, -0.3, 0.0, 1.0, 1.0, 1e-5, 1e-5, 0.3, 1e-3, 0.1, 1).simulate() - -# example 2 electron scheme -potential3, current3 = EE(-0.5, 0.4, -0.1, 0.1, 1.0, 1, 1e-5, 1e-5, 1e-5, 0.5, 0.5, 1e-2, 5e-3).simulate() - -fig, ax = plt.subplots() -ax.plot(potential1, [x*1000 for x in current1], label='$E_{q}$') -ax.plot(potential2, [x*1000 for x in current2], label='$E_{q}C$') -ax.plot(potential3, [x*1000 for x in current3], label='$E_{q}E_{q}$') -ax.xlabel('Potential (V)') -ax.ylabel('Current (mA)') -ax.tight_layout() -ax.legend() -plt.show() -############################################################################## -############################################################################## -# Part B: how to fit experimental data - -#for simplicity an 'experiment' dataset is simulated and then fit. Ideally you -#would input real data to fit - -################## 1) 'experiment' data - -fig, ax = plt.subplots() -ax.plot(potential1, [x*1000 for x in current1], label='experiment') -plt.show() -################## 2) turn test data into real function -new_p = np.linspace(1, 1200, 1200) -fig, ax = plt.subplots() -ax.plot(new_p, [x*1000 for x in current1], label='experiment') -plt.show() -################## 3) fit curve to data of real function -def func(x, a, b, c): #modify variables as needed - E_start = x[-1] - E_reverse = np.round(np.split(x, 2)[0][-1], 4) - print(E_reverse) - #here is where you choose what is known and unknown - setup1 = OneElectronCV(E_start, E_reverse, 0, 1, 1, 1, a, b, 5, 298) - potential_fit, current_fit = setup1.quasireversible(0.5, c) - return current_fit - -popt, pcov = curve_fit(func, potential, current, bounds=(0, [1e-5, 1e-5, 1e-3]))#p0=(1e-9, 2e-9, 0.09) -fitted_current = func(potential, *popt) -sigma = np.sqrt(np.diag(pcov)) # one standard deviation of the parameters -# print fits and std dev. for a,b,c etc. -print([f"{x:.2E} +/- {y:.0E}" for x,y in zip(popt, sigma)]) - -plt.figure(3) -plt.plot(new_p, [x*1000 for x in fitted_current], '--', label='Fit') -plt.xlabel('Timestep') -plt.ylabel('Current (mA)') -plt.tight_layout() -plt.legend() -plt.show() - -############## 4) turn fit back into CV form and overlay experiment data -plt.figure(2) -plt.plot(potential, [x*1000 for x in fitted_current], '--', label='Fit') -plt.xlabel('Potential (V)') -plt.ylabel('Current (mA)') -plt.tight_layout() -plt.legend() -plt.show() - -plt.figure(4) -plt.plot(new_p, [((x - y)*1e6) for x,y in zip(current, fitted_current)], '--') -plt.xlabel('Timestep') -plt.ylabel("Residuals ($\mu$A) ") -plt.tight_layout() -plt.show() - diff --git a/pyproject.toml b/pyproject.toml index aae43d0..9c32128 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,7 +5,7 @@ build-backend = "hatchling.build" [project] name = "cvsim" -version = "0.0.1" +version = "1.0.0" authors = [ { name="Eric Fell", email="efell@g.harvard.edu" }, { name="Jeremy Fell", email="jfell@sfu.ca"}, diff --git a/src/cvsim/fit_curve.py b/src/cvsim/fit_curve.py index 51acd51..84c20e2 100644 --- a/src/cvsim/fit_curve.py +++ b/src/cvsim/fit_curve.py @@ -107,7 +107,7 @@ def __init__( thetas = [round((i - delta_theta)) for i in [start_potential_mv, switch_potential_mv]] forward_scan = np.arange(thetas[0], thetas[1], step=delta_theta * -1) reverse_scan = np.append(forward_scan[-2::-1], start_potential_mv) - self.voltage_to_fit = np.concatenate([forward_scan, reverse_scan]) / 1000 + self.voltage = np.concatenate([forward_scan, reverse_scan]) / 1000 # Contains only variables with a user-specified fixed value. # These params are shared by all CVsim mechanisms @@ -121,8 +121,8 @@ def __init__( # These params are shared by all CVsim mechanisms self.default_vars = { 'reduction_potential': [ - round((self.voltage_to_fit[np.argmax(self.current_to_fit)] - + self.voltage_to_fit[np.argmin(self.current_to_fit)]) / 2, 3), + round((self.voltage[np.argmax(self.current_to_fit)] + + self.voltage[np.argmin(self.current_to_fit)]) / 2, 3), min(self.start_potential, self.switch_potential), max(self.start_potential, self.switch_potential), ], @@ -236,14 +236,14 @@ def fit_function( # fit raw data but exclude first data point, as semi-analytical method skips time=0 fit_results = curve_fit( f=fit_function, - xdata=self.voltage_to_fit, + xdata=self.voltage, ydata=self.current_to_fit[1:], p0=initial_guesses, bounds=[lower_bounds, upper_bounds], ) popt, pcov = list(fit_results) - current_fit = fit_function(self.voltage_to_fit, *popt) + current_fit = fit_function(self.voltage, *popt) sigma = np.sqrt(np.diag(pcov)) # one standard deviation of the parameters final_fit = {} @@ -253,13 +253,13 @@ def fit_function( # Semi-analytical method does not compute the first point (i.e. time=0) # so the starting voltage data point with a zero current is reinserted - self.voltage_to_fit = np.insert(self.voltage_to_fit, 0, self.start_potential) + self.voltage = np.insert(self.voltage, 0, self.start_potential) current_fit = np.insert(current_fit, 0, 0) - return self.voltage_to_fit, current_fit, final_fit + return self.voltage, current_fit, final_fit class FitE_rev(FitMechanism): - """Scheme for fitting a CV for a reversible (Nernstian) one electron transfer mechanism.""" + """Scheme for fitting a CV for a reversible (Nernstian) one-electron transfer mechanism.""" def _scheme(self, get_var: Callable[[str], float]) -> CyclicVoltammetryScheme: return E_rev( @@ -282,7 +282,7 @@ def fit( diffusion_product: _ParamGuess = None, ) -> tuple[np.ndarray, np.ndarray, dict]: """ - Fits the CV for a reversible (Nernstian) one electron transfer mechanism. + Fits the CV for a reversible (Nernstian) one-electron transfer mechanism. If a parameter is given, it must be a: float for initial guess of parameter; tuple[float, float] for (lower bound, upper bound) of the initial guess; or tuple[float, float, float] for (initial guess, lower bound, upper bound). @@ -301,7 +301,7 @@ def fit( Returns ------- - voltage_to_fit : np.ndarray + voltage : np.ndarray Array of potential (V) values of the CV fit. current_fit : np.ndarray Array of current (A) values of the CV fit. @@ -319,7 +319,7 @@ def fit( class FitE_q(FitMechanism): """ - Scheme for fitting a CV for a quasi-reversible one electron transfer mechanism. + Scheme for fitting a CV for a quasi-reversible one-electron transfer mechanism. Parameters ---------- @@ -426,7 +426,7 @@ def fit( k0: _ParamGuess = None, ) -> tuple[np.ndarray, np.ndarray, dict]: """ - Fits the CV for a quasi-reversible one electron transfer mechanism. + Fits the CV for a quasi-reversible one-electron transfer mechanism. If a parameter is given, it must be a: float for initial guess of parameter; tuple[float, float] for (lower bound, upper bound) of the initial guess; or tuple[float, float, float] for (initial guess, lower bound, upper bound). @@ -451,7 +451,7 @@ def fit( Returns ------- - voltage_to_fit : np.ndarray + voltage : np.ndarray Array of potential (V) values of the CV fit. current_fit : np.ndarray Array of current (A) values of the CV fit. @@ -470,7 +470,7 @@ def fit( class FitE_qC(FitMechanism): """ - Scheme for fitting a CV for a quasi-reversible one electron transfer, followed by + Scheme for fitting a CV for a quasi-reversible one-electron transfer, followed by a reversible first order homogeneous chemical transformation mechanism. Parameters @@ -598,7 +598,7 @@ def fit( k_backward: _ParamGuess = None, ) -> tuple[np.ndarray, np.ndarray, dict]: """ - Fits the CV for a quasi-reversible one electron transfer, followed by a reversible first + Fits the CV for a quasi-reversible one-electron transfer, followed by a reversible first order homogeneous chemical transformation mechanism. If a parameter is given, it must be a: float for initial guess of parameter; tuple[float, float] for (lower bound, upper bound) of the initial guess; or tuple[float, float, float] for @@ -630,7 +630,7 @@ def fit( Returns ------- - voltage_to_fit : np.ndarray + voltage : np.ndarray Array of potential (V) values of the CV fit. current_fit : np.ndarray Array of current (A) values of the CV fit. @@ -756,8 +756,8 @@ def __init__( # default [initial guess, lower bound, upper bound] self.default_vars |= { 'reduction_potential2': [ - round((self.voltage_to_fit[np.argmax(self.current_to_fit)] - + self.voltage_to_fit[np.argmin(self.current_to_fit)]) / 2, 3), + round((self.voltage[np.argmax(self.current_to_fit)] + + self.voltage[np.argmin(self.current_to_fit)]) / 2, 3), min(self.start_potential, self.switch_potential), max(self.start_potential, self.switch_potential), ], @@ -838,7 +838,7 @@ def fit( Returns ------- - voltage_to_fit : np.ndarray + voltage : np.ndarray Array of potential (V) values of the CV fit. current_fit : np.ndarray Array of current (A) values of the CV fit. @@ -990,8 +990,8 @@ def __init__( # default [initial guess, lower bound, upper bound] self.default_vars |= { 'reduction_potential2': [ - round((self.voltage_to_fit[np.argmax(self.current_to_fit)] - + self.voltage_to_fit[np.argmin(self.current_to_fit)]) / 2, 3), + round((self.voltage[np.argmax(self.current_to_fit)] + + self.voltage[np.argmin(self.current_to_fit)]) / 2, 3), min(self.start_potential, self.switch_potential), max(self.start_potential, self.switch_potential), ], @@ -1090,7 +1090,7 @@ def fit( Returns ------- - voltage_to_fit : np.ndarray + voltage : np.ndarray Array of potential (V) values of the CV fit. current_fit : np.ndarray Array of current (A) values of the CV fit. diff --git a/src/cvsim/mechanisms.py b/src/cvsim/mechanisms.py index 885092f..76edf49 100644 --- a/src/cvsim/mechanisms.py +++ b/src/cvsim/mechanisms.py @@ -161,7 +161,7 @@ def _voltage_profile_setup(self, reduction_potential2: float | None = None) -> t @abstractmethod def simulate(self) -> tuple[np.ndarray, np.ndarray]: """ - Simulates current-potential profile for desired mechanism + Simulates current-potential profile for desired mechanism. Returns ------- @@ -176,14 +176,14 @@ def simulate(self) -> tuple[np.ndarray, np.ndarray]: class E_rev(CyclicVoltammetryScheme): """ - Provides a current-potential profile for a reversible (Nernstian) one electron transfer mechanism. + Provides a current-potential profile for a reversible (Nernstian) one-electron transfer mechanism. This is equation (7:7) from [1]. """ def simulate(self) -> tuple[np.ndarray, np.ndarray]: """ - Simulates the CV for a reversible one electron transfer mechanism. + Simulates the CV for a reversible one-electron transfer mechanism. Returns ------- @@ -208,7 +208,7 @@ def simulate(self) -> tuple[np.ndarray, np.ndarray]: class E_q(CyclicVoltammetryScheme): """ - Provides a current-potential profile for a quasi-reversible one electron transfer mechanism. + Provides a current-potential profile for a quasi-reversible one-electron transfer mechanism. This is equation (8:3) from [1]. Parameters @@ -278,7 +278,7 @@ def __init__( def simulate(self) -> tuple[np.ndarray, np.ndarray]: """ - Simulates the CV for a quasi-reversible one electron transfer mechanism. + Simulates the CV for a quasi-reversible one-electron transfer mechanism. Returns ------- @@ -307,7 +307,7 @@ def simulate(self) -> tuple[np.ndarray, np.ndarray]: class E_qC(CyclicVoltammetryScheme): """ - Provides a current-potential profile for a quasi-reversible one electron transfer, followed by a reversible first + Provides a current-potential profile for a quasi-reversible one-electron transfer, followed by a reversible first order homogeneous chemical transformation mechanism. This is equation (10:4) from [1]. @@ -388,7 +388,7 @@ def __init__( def simulate(self) -> tuple[np.ndarray, np.ndarray]: """ - Simulates the CV for a quasi-reversible one electron transfer followed by a reversible first order homogeneous + Simulates the CV for a quasi-reversible one-electron transfer followed by a reversible first order homogeneous chemical transformation mechanism. Returns @@ -580,7 +580,7 @@ def simulate(self) -> tuple[np.ndarray, np.ndarray]: class SquareScheme(CyclicVoltammetryScheme): """ Provides a current-potential profile for two quasi-reversible, one-electron transfers of homogeneously - interconverting reactants (square scheme). + interconverting reactants (Square Scheme). This is equations (14:14 and 14:15) from [1]. Parameters @@ -685,7 +685,7 @@ def __init__( def simulate(self) -> tuple[np.ndarray, np.ndarray]: """ Simulates the CV for two quasi-reversible, one-electron transfers of homogeneously interconverting - reactants (square scheme). + reactants (Square Scheme). Returns -------