diff --git a/other/materials_designer/workflows/requirements-with-made.txt b/other/materials_designer/workflows/requirements-with-made.txt
new file mode 100644
index 00000000..c4a71fac
--- /dev/null
+++ b/other/materials_designer/workflows/requirements-with-made.txt
@@ -0,0 +1,5 @@
+numpy<2.0
+scipy>=1.5.4
+matplotlib>=3.0.0
+ase>=3.22.1
+mat3ra-made[tools]>=2024.11.12.post0
diff --git a/other/materials_designer/workflows/vbo_polar.pyi b/other/materials_designer/workflows/vbo_polar.pyi
new file mode 100644
index 00000000..94e8314b
--- /dev/null
+++ b/other/materials_designer/workflows/vbo_polar.pyi
@@ -0,0 +1,193 @@
+# ------------------------------------------------------------------ #
+# Linear Fit of ESP for Polar Interface VBO Calculation              #
+# ------------------------------------------------------------------ #
+#
+# Reference: Choudhary & Garrity, arXiv:2401.02021 (InterMat)        #
+#                                                                    #
+# For polar interfaces, ESP shows linear gradient in bulk regions    #
+# due to internal electric field. We fit each slab region and use    #
+# the average value of the fit as the ESP reference.                 #
+#                                                                    #
+# VBO Calculation:                                                    #
+#   1. Fit interface profile over slab 1 region → Va_interface       #
+#   2. Fit interface profile over slab 2 region → Vb_interface       #
+#   3. Fit bulk left profile over slab 1 region → Va_bulk_left       #
+#   4. Fit bulk right profile over slab 2 region → Vb_bulk_right     #
+#   5. VBO = (∆V_interface) - (∆V_bulk)                              #
+#      where ∆V_interface = Vb_interface - Va_interface              #
+#            ∆V_bulk = Vb_bulk_right - Va_bulk_left                  #
+#                                                                    #
+# Input:                                                              #
+#   - profile_left, profile_right: ESP profiles for bulk materials   #
+#   - profile_interface: ESP profile for interface structure         #
+#                                                                    #
+# Output: VBO (Valence Band Offset)                                  #
+#                                                                    #
+# NEW: Slab boundaries auto-detected using fingerprint matching      #
+# ------------------------------------------------------------------ #
+import json
+
+import matplotlib
+import ase.io
+from mat3ra.made.material import Material
+from mat3ra.made.tools.convert import from_ase
+
+# Non-interactive backend for running the script on the server, if working in Jupyter, comment out the next line
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+from types import SimpleNamespace
+from scipy.stats import linregress
+
+# Read structure from pw_scf.out (generated by previous pw_scf calculation)
+# Material index: 0=Interface, 1=Left, 2=Right
+# Files are named: pw_scf.out, pw_scf.out-1, pw_scf.out-2
+pw_scf_output = f"./pw_scf.out"
+pw_scf_output_1 = f"./pw_scf.out-1"
+pw_scf_output_2 = f"./pw_scf.out-2"
+
+# Read atomic structure from espresso output
+atoms = ase.io.read(pw_scf_output, format="espresso-out")
+atoms_1 = ase.io.read(pw_scf_output_1, format="espresso-out")
+atoms_2 = ase.io.read(pw_scf_output_2, format="espresso-out")
+
+# Convert ASE Atoms to Material using mat3ra-made
+material = Material.create(from_ase(atoms))
+material_1 = Material.create(from_ase(atoms_1))
+material_2 = Material.create(from_ase(atoms_2))
+
+material.to_cartesian()
+material_1.to_cartesian()
+material_2.to_cartesian()
+
+# Get the z-coordinate boundaries of each slab using element-based matching
+coords = material.basis.coordinates.values
+elements = material.basis.elements.values
+z_elements = sorted(zip([c[2] for c in coords], elements))
+n_left = len(material_1.basis.elements.values)
+
+z_max_1 = z_elements[n_left - 1][0]  # Last atom of left slab
+z_min_2 = z_elements[n_left][0]  # First atom of right slab
+z_min_1 = z_elements[0][0]
+z_max_2 = z_elements[-1][0]
+
+print(f"Detected Slab 1 (left) boundaries: z = {z_min_1:.3f} to {z_max_1:.3f} Å")
+print(f"Detected Slab 2 (right) boundaries: z = {z_min_2:.3f} to {z_max_2:.3f} Å")
+
+# Data from context: macroscopic average potential profile
+CHECKPOINT_FILE = "./.mat3ra/checkpoint.json"
+with open(CHECKPOINT_FILE, "r") as f:
+    checkpoint_data = json.load(f)
+    profile_interface = SimpleNamespace(
+        **checkpoint_data["scope"]["local"]["average-electrostatic-potential"]["average_potential_profile"]
+    )
+    profile_left = SimpleNamespace(
+        **checkpoint_data["scope"]["local"]["average-electrostatic-potential-left"]["average_potential_profile"]
+    )
+    profile_right = SimpleNamespace(
+        **checkpoint_data["scope"]["local"]["average-electrostatic-potential-right"]["average_potential_profile"]
+    )
+
+# Interface ESP profile
+X = np.array(profile_interface.xDataArray)  # z-coordinates (angstrom)
+Y = np.array(profile_interface.yDataSeries[1])  # Macroscopic average V̄(z)
+
+# Bulk material ESP profiles
+X_left = np.array(profile_left.xDataArray)
+Y_left = np.array(profile_left.yDataSeries[1])
+X_right = np.array(profile_right.xDataArray)
+Y_right = np.array(profile_right.yDataSeries[1])
+
+def get_region_indices(x_data, x_min, x_max):
+    """Get array indices corresponding to coordinate range."""
+    mask = (x_data >= x_min) & (x_data <= x_max)
+    indices = np.where(mask)[0]
+    if len(indices) == 0:
+        return 0, len(x_data)
+    return indices[0], indices[-1] + 1
+
+def fit_and_average(x_data, y_data, start_idx, end_idx):
+    """
+    Fit linear regression to region and return average value, slope, and intercept.
+
+    The average of the fitted line equals the mean of y-values,
+    but fitting helps smooth out oscillations and validates linearity.
+    """
+    x_region = x_data[start_idx:end_idx]
+    y_region = y_data[start_idx:end_idx]
+
+    if len(x_region) < 2:
+        avg = float(np.mean(y_region)) if len(y_region) > 0 else 0.0
+        return avg, 0.0, avg
+
+    slope, intercept, r_value, _, _ = linregress(x_region, y_region)
+
+    # Average value of linear fit over the region
+    # V_avg = (1/L) * integral(slope*x + intercept) = slope*x_mid + intercept
+    x_mid = (x_region[0] + x_region[-1]) / 2.0
+    avg_value = slope * x_mid + intercept
+
+    return float(avg_value), float(slope), float(intercept)
+
+# Get indices for each slab region in interface profile
+slab1_start, slab1_end = get_region_indices(X, z_min_1, z_max_1)
+slab2_start, slab2_end = get_region_indices(X, z_min_2, z_max_2)
+
+# Fit interface regions to get average ESP at interface
+Va_interface, slope_a, intercept_a = fit_and_average(X, Y, slab1_start, slab1_end)
+Vb_interface, slope_b, intercept_b = fit_and_average(X, Y, slab2_start, slab2_end)
+
+# Get indices for slab regions in bulk profiles
+slab1_start_left, slab1_end_left = get_region_indices(X_left, z_min_1, z_max_1)
+slab2_start_right, slab2_end_right = get_region_indices(X_right, z_min_2, z_max_2)
+
+# Fit bulk material profiles over the same z-ranges as their slabs
+Va_bulk_left, _, _ = fit_and_average(X_left, Y_left, slab1_start_left, slab1_end_left)
+Vb_bulk_right, _, _ = fit_and_average(X_right, Y_right, slab2_start_right, slab2_end_right)
+
+# Calculate VBO using interface and bulk ESP values
+# VBO = (interface difference) - (bulk difference)
+VBO = (Vb_interface - Va_interface) - (Vb_bulk_right - Va_bulk_left)
+
+print(f"Interface ESP Slab 1 (Va_interface): {Va_interface:.3f} eV")
+print(f"Interface ESP Slab 2 (Vb_interface): {Vb_interface:.3f} eV")
+print(f"Bulk ESP Left (Va_bulk): {Va_bulk_left:.3f} eV")
+print(f"Bulk ESP Right (Vb_bulk): {Vb_bulk_right:.3f} eV")
+print(f"Interface ∆V: {Vb_interface - Va_interface:.3f} eV")
+print(f"Bulk ∆V: {Vb_bulk_right - Va_bulk_left:.3f} eV")
+print(f"Valence Band Offset (VBO): {VBO:.3f} eV")
+
+# Generate visualization plot
+plt.figure(figsize=(10, 6))
+plt.plot(X, Y, label="Macroscopic Average Potential", linewidth=2)
+
+# Highlight fitting regions
+plt.axvspan(z_min_1, z_max_1, color="red", alpha=0.2, label="Slab 1 Region")
+plt.axvspan(z_min_2, z_max_2, color="blue", alpha=0.2, label="Slab 2 Region")
+
+# Plot fitted lines
+if slab1_end > slab1_start:
+    x_fit1 = X[slab1_start:slab1_end]
+    y_fit1 = slope_a * x_fit1 + intercept_a
+    plt.plot(x_fit1, y_fit1, color="darkred", linestyle="--", linewidth=2, label="Fit Slab 1")
+
+if slab2_end > slab2_start:
+    x_fit2 = X[slab2_start:slab2_end]
+    y_fit2 = slope_b * x_fit2 + intercept_b
+    plt.plot(x_fit2, y_fit2, color="darkblue", linestyle="--", linewidth=2, label="Fit Slab 2")
+
+# Plot average ESP values
+plt.axhline(Va_interface, color="red", linestyle=":", linewidth=2, label=f"Avg ESP Slab 1 = {Va_interface:.3f} eV")
+plt.axhline(Vb_interface, color="blue", linestyle=":", linewidth=2, label=f"Avg ESP Slab 2 = {Vb_interface:.3f} eV")
+
+plt.xlabel("z-coordinate (Å)", fontsize=12)
+plt.ylabel("Macroscopic Average Potential (eV)", fontsize=12)
+plt.title(f"Polar Interface VBO = {VBO:.3f} eV", fontsize=14, fontweight="bold")
+plt.legend(loc="best", fontsize=10)
+plt.grid(True, alpha=0.3)
+plt.tight_layout()
+# plt.show()
+
+filename = f"polar_vbo_fit_interface.png"
+plt.savefig(filename, dpi=150, bbox_inches="tight")
+plt.close()