Script calculates tafel slope or RDE ORR Data for Pt/C catalyst, starting from raw data text file containing iR corrected potential and capacitive current corrected current density.
import numpy as np
import pandas as pd
from pathlib import Path
# ============================================================
# 1. USER SETTINGS
# ============================================================
INPUT_FILE = "processed_lsv (1).txt"
OUTPUT_PLOTTABLE_FILE = "tafel_plottable_data.txt"
OUTPUT_SUMMARY_FILE = "tafel_summary.txt"
# Limiting current density magnitude for ORR, in mA/cm2
# Use positive value even if the measured ORR current is negative.
LIMITING_CURRENT_DENSITY_ABS = 6.0 # mA/cm2
# Potential range used for the actual Tafel linear fit
TAFEL_HIGH_V = 1.0
TAFEL_LOW_V = 0.94
# You can change this to 0.94 if needed:
# TAFEL_LOW_V = 0.94
# Potential range exported for plotting
PLOT_HIGH_V = 1.02
PLOT_LOW_V = 0.80
# ============================================================
# 2. HELPER FUNCTIONS
# ============================================================
def read_lsv_file(filename):
"""
Reads a text file containing:
column 1 = potential vs RHE, V
column 2 = capacitive-corrected current density, mA/cm2
The function accepts tab, spaces, or comma-separated data.
"""
path = Path(filename)
if not path.exists():
raise FileNotFoundError(f"Input file not found: {filename}")
df = pd.read_csv(
path,
sep=r"\s+|,|\t",
engine="python",
comment="#"
)
if df.shape[1] < 2:
raise ValueError("The input file must contain at least two columns.")
potential = pd.to_numeric(df.iloc[:, 0], errors="coerce")
current = pd.to_numeric(df.iloc[:, 1], errors="coerce")
clean = pd.DataFrame({
"Potential_V_vs_RHE": potential,
"Current_mA_per_cm2": current
}).dropna()
clean = clean.sort_values("Potential_V_vs_RHE").reset_index(drop=True)
return clean
def closest_value(array, target):
"""
Finds the closest available value in array to the requested target.
"""
array = np.asarray(array)
idx = np.argmin(np.abs(array - target))
return array[idx]
def calculate_kinetic_current_abs(current_mA_cm2, jlim_abs):
"""
Calculates absolute kinetic current density using the Koutecky-Levich equation.
For cathodic ORR data:
measured current is usually negative.
We use magnitudes:
1/j = 1/j_k + 1/j_L
Therefore:
j_k = j * j_L / (j_L - j)
where j, j_k, and j_L are positive current-density magnitudes.
Units:
j = mA/cm2
j_L = mA/cm2
j_k = mA/cm2
"""
j_abs = np.abs(current_mA_cm2)
kinetic_abs = (j_abs * jlim_abs) / (jlim_abs - j_abs)
return kinetic_abs
def linear_fit_tafel(potential, log_jk):
"""
Fits:
E = intercept + slope * log10(abs(j_k))
The slope has units of V/dec.
The reported Tafel slope is abs(slope) * 1000 in mV/dec.
"""
slope_V_dec, intercept_V = np.polyfit(log_jk, potential, 1)
fitted = intercept_V + slope_V_dec * log_jk
ss_res = np.sum((potential - fitted) ** 2)
ss_tot = np.sum((potential - np.mean(potential)) ** 2)
if ss_tot == 0:
r_squared = np.nan
else:
r_squared = 1 - ss_res / ss_tot
tafel_slope_mV_dec = abs(slope_V_dec) * 1000
return slope_V_dec, intercept_V, tafel_slope_mV_dec, r_squared
# ============================================================
# 3. MAIN CALCULATION
# ============================================================
df = read_lsv_file(INPUT_FILE)
potential = df["Potential_V_vs_RHE"].to_numpy()
current = df["Current_mA_per_cm2"].to_numpy()
# Calculate absolute kinetic current density
df["Current_abs_mA_per_cm2"] = np.abs(df["Current_mA_per_cm2"])
df["Kinetic_current_abs_mA_per_cm2"] = calculate_kinetic_current_abs(
df["Current_mA_per_cm2"],
LIMITING_CURRENT_DENSITY_ABS
)
# Keep only physically meaningful values:
# current magnitude must be positive and smaller than limiting current
valid = (
(df["Current_abs_mA_per_cm2"] > 0) &
(df["Current_abs_mA_per_cm2"] < LIMITING_CURRENT_DENSITY_ABS) &
(df["Kinetic_current_abs_mA_per_cm2"] > 0) &
np.isfinite(df["Kinetic_current_abs_mA_per_cm2"])
)
df_valid = df.loc[valid].copy()
df_valid["log10_abs_kinetic_current_mA_per_cm2"] = np.log10(
df_valid["Kinetic_current_abs_mA_per_cm2"]
)
# ============================================================
# 4. FIND CLOSEST REAL POTENTIALS FOR FIT RANGE
# ============================================================
available_potentials = df_valid["Potential_V_vs_RHE"].to_numpy()
actual_fit_high = closest_value(available_potentials, TAFEL_HIGH_V)
actual_fit_low = closest_value(available_potentials, TAFEL_LOW_V)
fit_upper = max(actual_fit_high, actual_fit_low)
fit_lower = min(actual_fit_high, actual_fit_low)
fit_mask = (
(df_valid["Potential_V_vs_RHE"] >= fit_lower) &
(df_valid["Potential_V_vs_RHE"] <= fit_upper)
)
fit_df = df_valid.loc[fit_mask].copy()
if len(fit_df) < 3:
raise ValueError(
"Not enough valid data points in the selected Tafel range. "
"Try widening the potential range."
)
slope_V_dec, intercept_V, tafel_slope_mV_dec, r_squared = linear_fit_tafel(
fit_df["Potential_V_vs_RHE"].to_numpy(),
fit_df["log10_abs_kinetic_current_mA_per_cm2"].to_numpy()
)
fit_df["Fitted_potential_V_vs_RHE"] = (
intercept_V +
slope_V_dec * fit_df["log10_abs_kinetic_current_mA_per_cm2"]
)
# ============================================================
# 5. EXPORT PLOTTABLE DATA FROM 1.02 V TO 0.80 V
# ============================================================
actual_plot_high = closest_value(available_potentials, PLOT_HIGH_V)
actual_plot_low = closest_value(available_potentials, PLOT_LOW_V)
plot_upper = max(actual_plot_high, actual_plot_low)
plot_lower = min(actual_plot_high, actual_plot_low)
plot_mask = (
(df_valid["Potential_V_vs_RHE"] >= plot_lower) &
(df_valid["Potential_V_vs_RHE"] <= plot_upper)
)
plot_df = df_valid.loc[plot_mask, [
"Potential_V_vs_RHE",
"log10_abs_kinetic_current_mA_per_cm2"
]].copy()
plot_df.to_csv(
OUTPUT_PLOTTABLE_FILE,
sep="\t",
index=False,
float_format="%.8f"
)
# ============================================================
# 6. EXPORT SUMMARY
# ============================================================
with open(OUTPUT_SUMMARY_FILE, "w", encoding="utf-8") as f:
f.write("ORR Tafel Slope Calculation Summary\n")
f.write("==================================\n\n")
f.write(f"Input file: {INPUT_FILE}\n")
f.write(f"Limiting current density magnitude used: {LIMITING_CURRENT_DENSITY_ABS:.6f} mA/cm2\n\n")
f.write("Koutecky-Levich correction used:\n")
f.write("1/j = 1/j_k + 1/j_L\n")
f.write("j_k = j * j_L / (j_L - j)\n")
f.write("where j, j_k, and j_L are absolute current-density magnitudes in mA/cm2.\n\n")
f.write("Tafel fit equation used:\n")
f.write("E = intercept + slope * log10(abs(j_k))\n\n")
f.write("Requested Tafel fit range:\n")
f.write(f"{TAFEL_HIGH_V:.6f} V to {TAFEL_LOW_V:.6f} V vs RHE\n\n")
f.write("Actual closest available fit range used:\n")
f.write(f"{actual_fit_high:.8f} V to {actual_fit_low:.8f} V vs RHE\n\n")
f.write(f"Number of points in fit: {len(fit_df)}\n")
f.write(f"Slope from E vs log10(abs(j_k)): {slope_V_dec:.8f} V/dec\n")
f.write(f"Tafel slope magnitude: {tafel_slope_mV_dec:.2f} mV/dec\n")
f.write(f"Intercept: {intercept_V:.8f} V\n")
f.write(f"R_squared: {r_squared:.6f}\n\n")
f.write("Plottable output file:\n")
f.write(f"{OUTPUT_PLOTTABLE_FILE}\n")
f.write("Columns:\n")
f.write("1. Potential_V_vs_RHE\n")
f.write("2. log10_abs_kinetic_current_mA_per_cm2\n")
# ============================================================
# 7. PRINT SUMMARY TO SCREEN
# ============================================================
print("ORR Tafel slope calculation completed.")
print(f"Input file: {INPUT_FILE}")
print(f"Limiting current density used: {LIMITING_CURRENT_DENSITY_ABS:.3f} mA/cm2")
print()
print(f"Requested fit range: {TAFEL_HIGH_V:.3f} to {TAFEL_LOW_V:.3f} V vs RHE")
print(f"Actual closest fit range used: {actual_fit_high:.8f} to {actual_fit_low:.8f} V vs RHE")
print(f"Number of fit points: {len(fit_df)}")
print()
print(f"Slope from E vs log10(abs(j_k)): {slope_V_dec:.8f} V/dec")
print(f"Tafel slope magnitude: {tafel_slope_mV_dec:.2f} mV/dec")
print(f"R_squared: {r_squared:.6f}")
print()
print(f"Plottable data saved to: {OUTPUT_PLOTTABLE_FILE}")
print(f"Summary saved to: {OUTPUT_SUMMARY_FILE}")