vbo_polar.pyi

# ------------------------------------------------------------------ #
# 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"
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()