From 457432683bd102561bf1d8606ffcfda084f2e304 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Sat, 17 Jan 2026 18:28:17 -0800 Subject: [PATCH 1/6] update: add python code to get VBO polar --- .../workflows/requirements-with-made.txt | 5 + .../workflows/vbo_polar.pyi | 143 ++++++++++++++++++ 2 files changed, 148 insertions(+) create mode 100644 other/materials_designer/workflows/requirements-with-made.txt create mode 100644 other/materials_designer/workflows/vbo_polar.pyi 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..4c4aee16 --- /dev/null +++ b/other/materials_designer/workflows/requirements-with-made.txt @@ -0,0 +1,5 @@ +munch == 2.5.0 +numpy >= 1.19.5 +scipy >= 1.5.4 +matplotlib >= 3.0.0 +mat3ra-made >= 2024.11.12.0 diff --git a/other/materials_designer/workflows/vbo_polar.pyi b/other/materials_designer/workflows/vbo_polar.pyi new file mode 100644 index 00000000..cbcb5346 --- /dev/null +++ b/other/materials_designer/workflows/vbo_polar.pyi @@ -0,0 +1,143 @@ +# ------------------------------------------------------------------ # +# 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 bulk ESP reference. # +# # +# Input: Coordinates defining slab regions (from structure data) # +# - z_min_1, z_max_1: z-coordinates of first material region # +# - z_min_2, z_max_2: z-coordinates of second material region # +# # +# Output: Va, Vb - average ESP values for each slab region # +# ------------------------------------------------------------------ # +import json + +import matplotlib +import ase.io +from mat3ra.made.material import Material +from mat3ra.made.tools.convert import from_ase + +matplotlib.use('Agg') # Non-interactive backend for headless environments +import matplotlib.pyplot as plt +import numpy as np +from munch import Munch +from scipy.stats import linregress + +# Material index: 0=Interface, 1=Left, 2=Right + +# Read structure from pw_scf.out (generated by previous pw_scf calculation) +# Files are named: pw_scf.out (index 0), pw_scf.out-1 (index 1), pw_scf.out-2 (index 2) +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_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_1 = Material.create(from_ase(atoms_1)) +material_2 = Material.create(from_ase(atoms_2)) + +z_min_1 = material_1.basis.coordinates.get_min_value_along_axis("z") +z_max_1 = material_1.basis.coordinates.get_max_value_along_axis("z") +z_min_2 = material_2.basis.coordinates.get_min_value_along_axis("z") +z_max_2 = material_2.basis.coordinates.get_max_value_along_axis("z") + + + +# Data from context: macroscopic average potential profile +{ % raw %}profile = {{average_potential_profile}}{ % endraw %} + +X = np.array(profile.xDataArray) # z-coordinates (angstrom) +Y = np.array(profile.yDataSeries[1]) # Macroscopic average V̄(z) + + +# Slab region coordinates (passed from workflow, set by user based on structure) +# These define z-coordinate ranges for fitting in each material region +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 +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 and get average ESP for each slab +Va, slope_a, intercept_a = fit_and_average(X, Y, slab1_start, slab1_end) +Vb, slope_b, intercept_b = fit_and_average(X, Y, slab2_start, slab2_end) + +# Output compatible with VBO workflow (same format as find_extrema.pyi) +result = { + "minima": [Va, Vb], + "maxima": [], +} + +print(json.dumps(result, indent=4)) + +# 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, color='red', linestyle=':', linewidth=2, label=f'Avg ESP Slab 1 = {Va:.3f} eV') +plt.axhline(Vb, color='blue', linestyle=':', linewidth=2, label=f'Avg ESP Slab 2 = {Vb:.3f} eV') + +plt.xlabel('z-coordinate (Å)', fontsize=12) +plt.ylabel('Macroscopic Average Potential (eV)', fontsize=12) +plt.title(f'Polar Interface VBO Calculation', fontsize=14, fontweight='bold') +plt.legend(loc='best', fontsize=10) +plt.grid(True, alpha=0.3) +plt.tight_layout() + +# Save plot to file with interface part in filename +filename = f'polar_vbo_fit_interface.png' +plt.savefig(filename, dpi=150, bbox_inches='tight') +plt.close() From 1832d1b24ed0803c05ba633c7cbe419e95c0e344 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Mon, 19 Jan 2026 10:54:58 -0800 Subject: [PATCH 2/6] update: get 4 values --- .../workflows/vbo_polar.pyi | 126 +++++++++++------- 1 file changed, 81 insertions(+), 45 deletions(-) diff --git a/other/materials_designer/workflows/vbo_polar.pyi b/other/materials_designer/workflows/vbo_polar.pyi index cbcb5346..8355d969 100644 --- a/other/materials_designer/workflows/vbo_polar.pyi +++ b/other/materials_designer/workflows/vbo_polar.pyi @@ -6,13 +6,24 @@ # # # 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 bulk ESP reference. # +# the average value of the fit as the ESP reference. # # # -# Input: Coordinates defining slab regions (from structure data) # -# - z_min_1, z_max_1: z-coordinates of first material region # -# - z_min_2, z_max_2: z-coordinates of second material region # +# 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 # # # -# Output: Va, Vb - average ESP values for each slab region # +# Input: # +# - profile_left, profile_right: ESP profiles for bulk materials # +# - profile_interface: ESP profile for interface structure # +# - z_min_1, z_max_1: z-coordinates of slab 1 (left material) # +# - z_min_2, z_max_2: z-coordinates of slab 2 (right material) # +# # +# Output: VBO (Valence Band Offset) # # ------------------------------------------------------------------ # import json @@ -21,16 +32,16 @@ import ase.io from mat3ra.made.material import Material from mat3ra.made.tools.convert import from_ase -matplotlib.use('Agg') # Non-interactive backend for headless environments +# 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 munch import Munch +from types import SimpleNamespace from scipy.stats import linregress -# Material index: 0=Interface, 1=Left, 2=Right - # Read structure from pw_scf.out (generated by previous pw_scf calculation) -# Files are named: pw_scf.out (index 0), pw_scf.out-1 (index 1), pw_scf.out-2 (index 2) +# 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_1 = f"./pw_scf.out-1" pw_scf_output_2 = f"./pw_scf.out-2" @@ -42,22 +53,37 @@ atoms_2 = ase.io.read(pw_scf_output_2, format="espresso-out") material_1 = Material.create(from_ase(atoms_1)) material_2 = Material.create(from_ase(atoms_2)) +material_1.to_cartesian() +material_2.to_cartesian() + z_min_1 = material_1.basis.coordinates.get_min_value_along_axis("z") z_max_1 = material_1.basis.coordinates.get_max_value_along_axis("z") z_min_2 = material_2.basis.coordinates.get_min_value_along_axis("z") z_max_2 = material_2.basis.coordinates.get_max_value_along_axis("z") - - # Data from context: macroscopic average potential profile -{ % raw %}profile = {{average_potential_profile}}{ % endraw %} +with open("./.mat3ra/checkpoint", "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]) -X = np.array(profile.xDataArray) # z-coordinates (angstrom) -Y = np.array(profile.yDataSeries[1]) # Macroscopic average V̄(z) - - -# Slab region coordinates (passed from workflow, set by user based on structure) -# These define z-coordinate ranges for fitting in each material region 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) @@ -66,7 +92,6 @@ def get_region_indices(x_data, x_min, x_max): 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. @@ -90,54 +115,65 @@ def fit_and_average(x_data, y_data, start_idx, end_idx): return float(avg_value), float(slope), float(intercept) - -# Get indices for each slab region +# 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 and get average ESP for each slab -Va, slope_a, intercept_a = fit_and_average(X, Y, slab1_start, slab1_end) -Vb, slope_b, intercept_b = fit_and_average(X, Y, slab2_start, slab2_end) +# 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) -# Output compatible with VBO workflow (same format as find_extrema.pyi) -result = { - "minima": [Va, Vb], - "maxima": [], -} +# 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(json.dumps(result, indent=4)) +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) +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') +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') + 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') + plt.plot(x_fit2, y_fit2, color="darkblue", linestyle="--", linewidth=2, label="Fit Slab 2") # Plot average ESP values -plt.axhline(Va, color='red', linestyle=':', linewidth=2, label=f'Avg ESP Slab 1 = {Va:.3f} eV') -plt.axhline(Vb, color='blue', linestyle=':', linewidth=2, label=f'Avg ESP Slab 2 = {Vb:.3f} eV') +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 Calculation', fontsize=14, fontweight='bold') -plt.legend(loc='best', fontsize=10) +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() -# Save plot to file with interface part in filename -filename = f'polar_vbo_fit_interface.png' -plt.savefig(filename, dpi=150, bbox_inches='tight') +filename = f"polar_vbo_fit_interface.png" +plt.savefig(filename, dpi=150, bbox_inches="tight") plt.close() From 050c4604929ddfd8174ccbc74064413e7f885b47 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Mon, 19 Jan 2026 12:48:02 -0800 Subject: [PATCH 3/6] update: use simple coords match --- .../workflows/vbo_polar.pyi | 28 ++++++++++++++----- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/other/materials_designer/workflows/vbo_polar.pyi b/other/materials_designer/workflows/vbo_polar.pyi index 8355d969..94e8314b 100644 --- a/other/materials_designer/workflows/vbo_polar.pyi +++ b/other/materials_designer/workflows/vbo_polar.pyi @@ -20,10 +20,10 @@ # Input: # # - profile_left, profile_right: ESP profiles for bulk materials # # - profile_interface: ESP profile for interface structure # -# - z_min_1, z_max_1: z-coordinates of slab 1 (left material) # -# - z_min_2, z_max_2: z-coordinates of slab 2 (right material) # # # # Output: VBO (Valence Band Offset) # +# # +# NEW: Slab boundaries auto-detected using fingerprint matching # # ------------------------------------------------------------------ # import json @@ -42,27 +42,41 @@ 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() -z_min_1 = material_1.basis.coordinates.get_min_value_along_axis("z") -z_max_1 = material_1.basis.coordinates.get_max_value_along_axis("z") -z_min_2 = material_2.basis.coordinates.get_min_value_along_axis("z") -z_max_2 = material_2.basis.coordinates.get_max_value_along_axis("z") +# 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 -with open("./.mat3ra/checkpoint", "r") as f: +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"] From d1a33cc7d15dad1e7f4603f5fcfdd74e92924fde Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Mon, 19 Jan 2026 16:27:44 -0800 Subject: [PATCH 4/6] update: reqs --- .../workflows/requirements-with-made.txt | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/other/materials_designer/workflows/requirements-with-made.txt b/other/materials_designer/workflows/requirements-with-made.txt index 4c4aee16..80a75e4e 100644 --- a/other/materials_designer/workflows/requirements-with-made.txt +++ b/other/materials_designer/workflows/requirements-with-made.txt @@ -1,5 +1,6 @@ -munch == 2.5.0 -numpy >= 1.19.5 -scipy >= 1.5.4 -matplotlib >= 3.0.0 -mat3ra-made >= 2024.11.12.0 +munch==2.5.0 +numpy>=1.19.5 +scipy>=1.5.4 +matplotlib>=3.0.0 +ase>=3.22.1 +mat3ra-made[tools]>=2024.11.12.0 From e522970f77a06ebd312fe9de85791cc9db832dc1 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Mon, 19 Jan 2026 18:18:06 -0800 Subject: [PATCH 5/6] update: reqs --- other/materials_designer/workflows/requirements-with-made.txt | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/other/materials_designer/workflows/requirements-with-made.txt b/other/materials_designer/workflows/requirements-with-made.txt index 80a75e4e..0b9105a3 100644 --- a/other/materials_designer/workflows/requirements-with-made.txt +++ b/other/materials_designer/workflows/requirements-with-made.txt @@ -1,5 +1,4 @@ -munch==2.5.0 -numpy>=1.19.5 +numpy<2.0 scipy>=1.5.4 matplotlib>=3.0.0 ase>=3.22.1 From 878de8c0e78965c45318ba91371ea2ba18015c96 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Mon, 19 Jan 2026 18:18:20 -0800 Subject: [PATCH 6/6] update: reqs --- other/materials_designer/workflows/requirements-with-made.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/other/materials_designer/workflows/requirements-with-made.txt b/other/materials_designer/workflows/requirements-with-made.txt index 0b9105a3..c4a71fac 100644 --- a/other/materials_designer/workflows/requirements-with-made.txt +++ b/other/materials_designer/workflows/requirements-with-made.txt @@ -2,4 +2,4 @@ numpy<2.0 scipy>=1.5.4 matplotlib>=3.0.0 ase>=3.22.1 -mat3ra-made[tools]>=2024.11.12.0 +mat3ra-made[tools]>=2024.11.12.post0