Hull

Dependencies

/// script requires-python = ">=3.10" dependencies = [ "a2c-ase @ git+https://github.com/abhijeetgangan/a2c_ase.git", "numpy", "pymatgen", "tqdm", "plotly", "kaleido", ] ///

Kob-Andersen Binary System - a2c Workflow

Hull exploration with a2c workflow for the classic Kob-Andersen binary Lennard-Jones glass former. Uses reduced units (LJ natural units: sigma, epsilon).

Setup and Imports¶

The Kob-Andersen model is a classic binary Lennard-Jones system that forms metallic glasses. We use reduced units: distances in sigma, energies in epsilon. We will demonstrate how we can explore the hull for a selected compositions.

In [1]:

Copied!





import os
from collections import Counter

import numpy as np
from ase import Atoms
from ase.build import bulk
from pymatgen.analysis.phase_diagram import PDEntry, PDPlotter, PhaseDiagram
from pymatgen.core.composition import Composition
from tqdm import tqdm

from a2c_ase.potentials.mlj import MultiLennardJones
from a2c_ase.runner import melt_quench_md, relax_unit_cell
from a2c_ase.utils import extract_crystallizable_subcells, random_packed_structure

IS_CI = os.getenv("CI") is not None
import os
from collections import Counter

import numpy as np
from ase import Atoms
from ase.build import bulk
from pymatgen.analysis.phase_diagram import PDEntry, PDPlotter, PhaseDiagram
from pymatgen.core.composition import Composition
from tqdm import tqdm

from a2c_ase.potentials.mlj import MultiLennardJones
from a2c_ase.runner import melt_quench_md, relax_unit_cell
from a2c_ase.utils import extract_crystallizable_subcells, random_packed_structure

IS_CI = os.getenv("CI") is not None

Kob-Andersen Parameters¶

Classic binary LJ glass former with reduced units (dimensionless):

A particles (Ni): 80%, σ=1.0, ε=1.0 (reference)
B particles (P): 20%, σ=0.88, ε=0.5
Cross: σ_AB=0.8, ε_AB=1.5
Glass transition: T_g ≈ 0.435 (in reduced units)

In [2]:

Copied!





# System configuration (80:20 A:B composition)
comp = Composition("Ni80P20")

# Cell in LJ units (sigma as length unit)
cell_size = 4.4
cell = np.array([[cell_size, 0.0, 0.0], [0.0, cell_size, 0.0], [0.0, 0.0, cell_size]])

# Kob-Andersen calculator (reduced/LJ units)
calculator = MultiLennardJones(
    sigma={"Ni": 1.0, "P": 0.88},  # LJ sigma in natural units
    epsilon={"Ni": 1.0, "P": 0.5},  # LJ epsilon in natural units
    cross_interactions={("Ni", "P"): {"sigma": 0.8, "epsilon": 1.5}},
    rc=2.5,  # Cutoff in units of sigma
    smooth=True,
)
# System configuration (80:20 A:B composition)
comp = Composition("Ni80P20")

# Cell in LJ units (sigma as length unit)
cell_size = 4.4
cell = np.array([[cell_size, 0.0, 0.0], [0.0, cell_size, 0.0], [0.0, 0.0, cell_size]])

# Kob-Andersen calculator (reduced/LJ units)
calculator = MultiLennardJones(
    sigma={"Ni": 1.0, "P": 0.88},  # LJ sigma in natural units
    epsilon={"Ni": 1.0, "P": 0.5},  # LJ epsilon in natural units
    cross_interactions={("Ni", "P"): {"sigma": 0.8, "epsilon": 1.5}},
    rc=2.5,  # Cutoff in units of sigma
    smooth=True,
)

Simulation Parameters¶

All parameters in reduced LJ units (dimensionless):

Energy: ε (epsilon) = 1.0
Distance: σ (sigma) = 1.0
Temperature: T* = kT/ε
Time: τ = √(mσ²/ε) (dimensionless)

In [3]:

Copied!





global_seed = 42
fmax = 0.01  # Force convergence in reduced units

# Reduce parameters for CI testing
max_iter = 20 if IS_CI else 200

# MD parameters (LJ units)
md_log_interval = 50
md_equi_steps = 100 if IS_CI else 2500
md_cool_steps = 100 if IS_CI else 2500
md_final_steps = 100 if IS_CI else 2500
md_T_high = 4.0  # High T* (reduced units, above glass transition ~0.8)
md_T_low = 0.4  # Low T* (reduced units, below glass transition)
md_time_step = 0.005  # Timestep in reduced units
md_friction = 1 / (100 * md_time_step)  # Friction coefficient

if IS_CI:
    print("Running in CI mode with reduced parameters")
global_seed = 42
fmax = 0.01  # Force convergence in reduced units

# Reduce parameters for CI testing
max_iter = 20 if IS_CI else 200

# MD parameters (LJ units)
md_log_interval = 50
md_equi_steps = 100 if IS_CI else 2500
md_cool_steps = 100 if IS_CI else 2500
md_final_steps = 100 if IS_CI else 2500
md_T_high = 4.0  # High T* (reduced units, above glass transition ~0.8)
md_T_low = 0.4  # Low T* (reduced units, below glass transition)
md_time_step = 0.005  # Timestep in reduced units
md_friction = 1 / (100 * md_time_step)  # Friction coefficient

if IS_CI:
    print("Running in CI mode with reduced parameters")

Running in CI mode with reduced parameters

Step 1: Generate Random Packed Structure¶

Create initial random configuration with A and B particles.

In [4]:

Copied!





packed_atoms, log_data = random_packed_structure(
    composition=comp,
    cell=cell,
    seed=global_seed,
    diameter=2.5,
    max_iter=max_iter,
    fmax=fmax,
    verbose=True,
    auto_diameter=False,
)
print(f"Generated packed structure: {packed_atoms}")
print(f"Number of Ni (A) atoms: {sum(1 for s in packed_atoms.symbols if s == 'Ni')}")
print(f"Number of P (B) atoms: {sum(1 for s in packed_atoms.symbols if s == 'P')}")
packed_atoms, log_data = random_packed_structure(
    composition=comp,
    cell=cell,
    seed=global_seed,
    diameter=2.5,
    max_iter=max_iter,
    fmax=fmax,
    verbose=True,
    auto_diameter=False,
)
print(f"Generated packed structure: {packed_atoms}")
print(f"Number of Ni (A) atoms: {sum(1 for s in packed_atoms.symbols if s == 'Ni')}")
print(f"Number of P (B) atoms: {sum(1 for s in packed_atoms.symbols if s == 'P')}")

Reduce atom overlap using the soft_sphere potential
Initial energy: 186.9468
Step: 0, E: 186.9468, Fmax: 1.7763, Min dist: 0.1037
Step: 1, E: 186.1019, Fmax: 1.7140, Min dist: 0.1120
Step: 2, E: 184.5282, Fmax: 1.5892, Min dist: 0.1259

Step: 3, E: 182.9326, Fmax: 1.4488, Min dist: 0.1413
Step: 4, E: 181.4750, Fmax: 1.3043, Min dist: 0.1570
Step: 5, E: 180.1550, Fmax: 1.1980, Min dist: 0.1732
Step: 6, E: 178.9719, Fmax: 1.0905, Min dist: 0.1900

Step: 7, E: 177.9247, Fmax: 0.9773, Min dist: 0.2075
Step: 8, E: 177.0115, Fmax: 0.8635, Min dist: 0.2259
Step: 9, E: 176.2302, Fmax: 0.7401, Min dist: 0.2453
Step: 10, E: 175.5772, Fmax: 0.6216, Min dist: 0.2661
Step: 11, E: 175.0457, Fmax: 0.5197, Min dist: 0.2883

Step: 12, E: 174.6261, Fmax: 0.4208, Min dist: 0.3123
Step: 13, E: 174.3075, Fmax: 0.3560, Min dist: 0.3384
Step: 14, E: 174.0753, Fmax: 0.3206, Min dist: 0.3669
Step: 15, E: 173.9128, Fmax: 0.3017, Min dist: 0.3980
Step: 16, E: 173.8007, Fmax: 0.3088, Min dist: 0.4317
Step: 17, E: 173.7181, Fmax: 0.3481, Min dist: 0.4679

Step: 18, E: 173.6440, Fmax: 0.3657, Min dist: 0.5062
Step: 19, E: 173.5592, Fmax: 0.3609, Min dist: 0.5459
Step: 20, E: 173.4504, Fmax: 0.3319, Min dist: 0.5582
Final energy: 173.4504
Generated packed structure: Atoms(symbols='Ni80P20', pbc=True, cell=[4.4, 4.4, 4.4], calculator=SoftSphere(...))
Number of Ni (A) atoms: 80
Number of P (B) atoms: 20

In [5]:

Copied!





# Relax the packed structure so that the initial structure doesn't have
# large forces for melt-quench MD.
packed_atoms, logger = relax_unit_cell(
    atoms=packed_atoms, calculator=calculator, max_iter=max_iter, fmax=fmax, verbose=True
)
# Relax the packed structure so that the initial structure doesn't have
# large forces for melt-quench MD.
packed_atoms, logger = relax_unit_cell(
    atoms=packed_atoms, calculator=calculator, max_iter=max_iter, fmax=fmax, verbose=True
)

Initial energy: 24103.049587 eV
Initial volume: 85.184 Å³
Initial pressure: 1208.881066 eV/Å³
Step 0: E = 24103.049587 eV, Fmax = 92475.904440 eV/Å, P = 1208.881066 eV/Å³, V = 85.184 Å³

Step 1: E = 9201.490457 eV, Fmax = 9653.251133 eV/Å, P = 500.139940 eV/Å³, V = 85.186 Å³

Step 2: E = 6933.804645 eV, Fmax = 13343.286963 eV/Å, P = 390.315734 eV/Å³, V = 85.190 Å³

Step 3: E = 5566.112354 eV, Fmax = 8726.639546 eV/Å, P = 323.785760 eV/Å³, V = 85.195 Å³

Step 4: E = 4520.803711 eV, Fmax = 5452.583236 eV/Å, P = 272.783127 eV/Å³, V = 85.200 Å³

Step 5: E = 3890.780941 eV, Fmax = 5077.364589 eV/Å, P = 241.627518 eV/Å³, V = 85.207 Å³
Step 6: E = 3381.857741 eV, Fmax = 4181.581220 eV/Å, P = 216.357357 eV/Å³, V = 85.215 Å³

Step 7: E = 2907.676084 eV, Fmax = 4455.688123 eV/Å, P = 192.861648 eV/Å³, V = 85.224 Å³

Step 8: E = 2517.345085 eV, Fmax = 4636.406730 eV/Å, P = 173.471604 eV/Å³, V = 85.234 Å³

Step 9: E = 2204.754508 eV, Fmax = 4495.762339 eV/Å, P = 157.920712 eV/Å³, V = 85.245 Å³

Step 10: E = 1969.702478 eV, Fmax = 2682.887341 eV/Å, P = 146.213460 eV/Å³, V = 85.258 Å³

Step 11: E = 1842.936639 eV, Fmax = 2818.194326 eV/Å, P = 139.733149 eV/Å³, V = 85.271 Å³

Step 12: E = 1145.771025 eV, Fmax = 730.623441 eV/Å, P = 104.316168 eV/Å³, V = 85.281 Å³

Step 13: E = 833.955456 eV, Fmax = 1207.098799 eV/Å, P = 87.827478 eV/Å³, V = 85.294 Å³
Step 14: E = 595.840582 eV, Fmax = 516.955791 eV/Å, P = 75.143666 eV/Å³, V = 85.310 Å³

Step 15: E = 416.414893 eV, Fmax = 384.055368 eV/Å, P = 65.443928 eV/Å³, V = 85.329 Å³
Step 16: E = 256.310322 eV, Fmax = 322.465777 eV/Å, P = 56.774559 eV/Å³, V = 85.349 Å³

Step 17: E = 123.041348 eV, Fmax = 232.079330 eV/Å, P = 49.503373 eV/Å³, V = 85.371 Å³
Step 18: E = 25.253097 eV, Fmax = 309.653331 eV/Å, P = 44.049215 eV/Å³, V = 85.395 Å³

Step 19: E = -52.095182 eV, Fmax = 452.306409 eV/Å, P = 39.653908 eV/Å³, V = 85.421 Å³
Step 20: E = -127.192693 eV, Fmax = 321.594908 eV/Å, P = 35.421277 eV/Å³, V = 85.448 Å³

Optimization completed:
Final energy: -127.192693 eV
Final volume: 85.448 Å³
Final pressure: 35.421277 eV/Å³
Steps taken: 20

Step 2: Melt-Quench MD Simulation¶

Heat to T=4.0, quench to T=0.4 (near T_g≈0.435).

In [6]:

Copied!





amorphous_atoms, md_log = melt_quench_md(
    atoms=packed_atoms,
    calculator=calculator,
    equi_steps=md_equi_steps,
    cool_steps=md_cool_steps,
    final_steps=md_final_steps,
    T_high=md_T_high,
    T_low=md_T_low,
    time_step=md_time_step,
    friction=md_friction,
    seed=global_seed,
    verbose=True,
    log_interval=md_log_interval,
)
print(f"Amorphous structure ready: {amorphous_atoms}")
amorphous_atoms, md_log = melt_quench_md(
    atoms=packed_atoms,
    calculator=calculator,
    equi_steps=md_equi_steps,
    cool_steps=md_cool_steps,
    final_steps=md_final_steps,
    T_high=md_T_high,
    T_low=md_T_low,
    time_step=md_time_step,
    friction=md_friction,
    seed=global_seed,
    verbose=True,
    log_interval=md_log_interval,
)
print(f"Amorphous structure ready: {amorphous_atoms}")

Step 0/300: T = 4.0 K, E_pot = -127.193 eV, E_kin = 0.052 eV

Step 50/300: T = 387.7 K, E_pot = -132.315 eV, E_kin = 5.012 eV

Step 100/300: T = 1372.4 K, E_pot = -146.128 eV, E_kin = 17.739 eV

Step 150/300: T = 2665.8 K, E_pot = -165.371 eV, E_kin = 34.458 eV

Step 200/300: T = 4000.5 K, E_pot = -186.808 eV, E_kin = 51.710 eV

Step 250/300: T = 5196.7 K, E_pot = -208.140 eV, E_kin = 67.173 eV

Melt-quench simulation completed:
Final temperature: 6187.7 K
Final energy: -228.141 eV
Amorphous structure ready: Atoms(symbols='Ni80P20', pbc=True, cell=[[4.404302085217416, 0.0004704354711196714, -0.00011483479060379136], [0.0004703656985870155, 4.405225004480768, 0.0005501805715166627], [-0.00011486374903597653, 0.0005501986925379699, 4.404104872625505]], momenta=..., calculator=MultiLennardJones(...))

Step 3: Extract Crystallizable Subcells¶

Divide the glass into overlapping subcells to search for local crystalline order. In this example, we will only consider the following compositions: Ni, NiP, Ni2P, Ni4P, NiP2, NiP4, P.

In [7]:

Copied!





crystallizable_cells = extract_crystallizable_subcells(
    atoms=amorphous_atoms,
    d_frac=0.25,  # Grid spacing (larger for binary system)
    n_min=2,
    n_max=8,  # Allow larger subcells for binary compounds
    cubic_only=False,  # Allow non-cubic structures
    allowed_atom_counts=None,  # Don't restrict by count
    restrict_to_compositions=[
        "Ni",
        "NiP",
        "Ni2P",
        "Ni4P",
        "NiP2",
        "NiP4",
        "P",
    ],  # Only consider these compositions
    max_coeff=None,
    elements=None,
)
print(f"Extracted {len(crystallizable_cells)} candidate subcells")
crystallizable_cells = extract_crystallizable_subcells(
    atoms=amorphous_atoms,
    d_frac=0.25,  # Grid spacing (larger for binary system)
    n_min=2,
    n_max=8,  # Allow larger subcells for binary compounds
    cubic_only=False,  # Allow non-cubic structures
    allowed_atom_counts=None,  # Don't restrict by count
    restrict_to_compositions=[
        "Ni",
        "NiP",
        "Ni2P",
        "Ni4P",
        "NiP2",
        "NiP4",
        "P",
    ],  # Only consider these compositions
    max_coeff=None,
    elements=None,
)
print(f"Extracted {len(crystallizable_cells)} candidate subcells")

Created 261 subcells from amorphous structure
Subcells kept after filtering: 261
Extracted 261 candidate subcells

Step 4: Optimize Subcells¶

Relax each subcell to find stable crystalline phases.

In [8]:

Copied!





relaxed_structures = []
print("Optimizing candidate structures...")

for atoms in tqdm(crystallizable_cells[:20] if IS_CI else crystallizable_cells):
    try:
        relaxed, logger = relax_unit_cell(
            atoms=atoms, calculator=calculator, max_iter=max_iter, fmax=fmax, verbose=False
        )

        final_energy = relaxed.get_potential_energy()
        energy_per_atom = final_energy / len(relaxed)

        relaxed_structures.append((relaxed, energy_per_atom, final_energy))
    except Exception as e:
        print(f"Optimization failed: {e}")
        continue

print(f"Successfully optimized {len(relaxed_structures)} structures")
relaxed_structures = []
print("Optimizing candidate structures...")

for atoms in tqdm(crystallizable_cells[:20] if IS_CI else crystallizable_cells):
    try:
        relaxed, logger = relax_unit_cell(
            atoms=atoms, calculator=calculator, max_iter=max_iter, fmax=fmax, verbose=False
        )

        final_energy = relaxed.get_potential_energy()
        energy_per_atom = final_energy / len(relaxed)

        relaxed_structures.append((relaxed, energy_per_atom, final_energy))
    except Exception as e:
        print(f"Optimization failed: {e}")
        continue

print(f"Successfully optimized {len(relaxed_structures)} structures")

Optimizing candidate structures...

  0%|          | 0/20 [00:00<?, ?it/s]

  5%|▌         | 1/20 [00:01<00:23,  1.22s/it]

 10%|█         | 2/20 [00:02<00:24,  1.36s/it]

 15%|█▌        | 3/20 [00:04<00:23,  1.35s/it]

 20%|██        | 4/20 [00:05<00:20,  1.31s/it]

 25%|██▌       | 5/20 [00:06<00:19,  1.31s/it]

 30%|███       | 6/20 [00:07<00:16,  1.16s/it]

 35%|███▌      | 7/20 [00:08<00:16,  1.24s/it]

 40%|████      | 8/20 [00:10<00:15,  1.28s/it]

 45%|████▌     | 9/20 [00:11<00:13,  1.22s/it]

 50%|█████     | 10/20 [00:12<00:13,  1.31s/it]

 55%|█████▌    | 11/20 [00:14<00:11,  1.29s/it]

 60%|██████    | 12/20 [00:15<00:10,  1.33s/it]

 65%|██████▌   | 13/20 [00:17<00:10,  1.44s/it]

 70%|███████   | 14/20 [00:18<00:08,  1.34s/it]

 75%|███████▌  | 15/20 [00:20<00:07,  1.48s/it]

 80%|████████  | 16/20 [00:21<00:05,  1.49s/it]

 85%|████████▌ | 17/20 [00:23<00:04,  1.51s/it]

 90%|█████████ | 18/20 [00:24<00:03,  1.55s/it]

 95%|█████████▌| 19/20 [00:25<00:01,  1.35s/it]

100%|██████████| 20/20 [00:27<00:00,  1.51s/it]

100%|██████████| 20/20 [00:27<00:00,  1.38s/it]

Successfully optimized 20 structures

Step 5: Construct Convex Hull¶

Determine thermodynamic stability using pymatgen's phase diagram.

In [9]:

Copied!





# Build convex hull with proper composition handling
entries = []
for atoms, _, total_e in relaxed_structures:
    # Get composition from symbol counts (more explicit and correct)
    symbol_counts = Counter(atoms.get_chemical_symbols())
    comp_obj = Composition(symbol_counts)  # e.g., {'Ni': 2, 'P': 1} -> Ni2P
    entries.append(PDEntry(comp_obj, total_e))

# Compute reference energies using proper crystal structures
print("\nComputing reference energies...")

# Pure Ni: FCC structure
ni_fcc = bulk("Ni", "fcc", a=1.5)
ni_relaxed, _ = relax_unit_cell(ni_fcc, calculator, max_iter=max_iter, fmax=fmax, verbose=False)
e_ni_total = ni_relaxed.get_potential_energy()
print(f"Ni (FCC): {e_ni_total / len(ni_relaxed):.4f} ε/atom")

# Pure P: FCC structure
p_fcc = bulk("P", "fcc", a=1.3)
p_relaxed, _ = relax_unit_cell(p_fcc, calculator, max_iter=max_iter, fmax=fmax, verbose=False)
e_p_total = p_relaxed.get_potential_energy()
print(f"P (FCC): {e_p_total / len(p_relaxed):.4f} ε/atom")

# Ni4P structure
alat = 4.0
ni4p = Atoms(
    "Ni4P",
    positions=[
        [0.0, 0.0, 0.0],
        [0.0, alat / 2, alat / 2],
        [alat / 2, 0.0, alat / 2],
        [alat / 2, alat / 2, 0.0],
        [alat / 2, alat / 2, alat / 2],
    ],
    cell=[alat, alat, alat],
    pbc=True,
)
ni4p_relaxed, _ = relax_unit_cell(ni4p, calculator, max_iter=max_iter, fmax=fmax, verbose=False)
e_ni4p_total = ni4p_relaxed.get_potential_energy()
print(f"Ni4P (Rocksalt): {e_ni4p_total / len(ni4p_relaxed):.4f} ε/atom")

# Add references with TOTAL energies
entries.extend(
    [
        PDEntry(Composition("Ni"), e_ni_total),
        PDEntry(Composition("P"), e_p_total),
        PDEntry(Composition("Ni4P"), e_ni4p_total),
    ]
)
pd = PhaseDiagram(entries)

for i, (atoms, e_per_atom, total_e) in enumerate(relaxed_structures):
    symbol_counts = Counter(atoms.get_chemical_symbols())
    comp_obj = Composition(symbol_counts)
    entry = PDEntry(comp_obj, total_e)
    e_above_hull = pd.get_e_above_hull(entry)

    print(
        f"{i + 1:2d}. {comp_obj.reduced_formula:10s} | "
        f"E/atom: {e_per_atom:8.4f} ε | "
        f"E_hull: {e_above_hull:8.4f} ε/atom"
    )

print(f"\nTotal structures analyzed: {len(relaxed_structures)}")
print(f"Total entries in phase diagram: {len(entries)}")

# Plot and save phase diagram
plotter = PDPlotter(pd, show_unstable=1000.0, backend="plotly")  # Very large value to show all
plotter.get_plot()
plotter.show()
# Build convex hull with proper composition handling
entries = []
for atoms, _, total_e in relaxed_structures:
    # Get composition from symbol counts (more explicit and correct)
    symbol_counts = Counter(atoms.get_chemical_symbols())
    comp_obj = Composition(symbol_counts)  # e.g., {'Ni': 2, 'P': 1} -> Ni2P
    entries.append(PDEntry(comp_obj, total_e))

# Compute reference energies using proper crystal structures
print("\nComputing reference energies...")

# Pure Ni: FCC structure
ni_fcc = bulk("Ni", "fcc", a=1.5)
ni_relaxed, _ = relax_unit_cell(ni_fcc, calculator, max_iter=max_iter, fmax=fmax, verbose=False)
e_ni_total = ni_relaxed.get_potential_energy()
print(f"Ni (FCC): {e_ni_total / len(ni_relaxed):.4f} ε/atom")

# Pure P: FCC structure
p_fcc = bulk("P", "fcc", a=1.3)
p_relaxed, _ = relax_unit_cell(p_fcc, calculator, max_iter=max_iter, fmax=fmax, verbose=False)
e_p_total = p_relaxed.get_potential_energy()
print(f"P (FCC): {e_p_total / len(p_relaxed):.4f} ε/atom")

# Ni4P structure
alat = 4.0
ni4p = Atoms(
    "Ni4P",
    positions=[
        [0.0, 0.0, 0.0],
        [0.0, alat / 2, alat / 2],
        [alat / 2, 0.0, alat / 2],
        [alat / 2, alat / 2, 0.0],
        [alat / 2, alat / 2, alat / 2],
    ],
    cell=[alat, alat, alat],
    pbc=True,
)
ni4p_relaxed, _ = relax_unit_cell(ni4p, calculator, max_iter=max_iter, fmax=fmax, verbose=False)
e_ni4p_total = ni4p_relaxed.get_potential_energy()
print(f"Ni4P (Rocksalt): {e_ni4p_total / len(ni4p_relaxed):.4f} ε/atom")

# Add references with TOTAL energies
entries.extend(
    [
        PDEntry(Composition("Ni"), e_ni_total),
        PDEntry(Composition("P"), e_p_total),
        PDEntry(Composition("Ni4P"), e_ni4p_total),
    ]
)
pd = PhaseDiagram(entries)

for i, (atoms, e_per_atom, total_e) in enumerate(relaxed_structures):
    symbol_counts = Counter(atoms.get_chemical_symbols())
    comp_obj = Composition(symbol_counts)
    entry = PDEntry(comp_obj, total_e)
    e_above_hull = pd.get_e_above_hull(entry)

    print(
        f"{i + 1:2d}. {comp_obj.reduced_formula:10s} | "
        f"E/atom: {e_per_atom:8.4f} ε | "
        f"E_hull: {e_above_hull:8.4f} ε/atom"
    )

print(f"\nTotal structures analyzed: {len(relaxed_structures)}")
print(f"Total entries in phase diagram: {len(entries)}")

# Plot and save phase diagram
plotter = PDPlotter(pd, show_unstable=1000.0, backend="plotly")  # Very large value to show all
plotter.get_plot()
plotter.show()

Computing reference energies...

Ni (FCC): -7.5898 ε/atom

P (FCC): -3.9678 ε/atom

Ni4P (Rocksalt): -0.0487 ε/atom
 1. Ni         | E/atom:  -6.2349 ε | E_hull:   1.3549 ε/atom
 2. Ni         | E/atom:  -6.0107 ε | E_hull:   1.5791 ε/atom
 3. Ni         | E/atom:  -7.5310 ε | E_hull:   0.0588 ε/atom
 4. Ni         | E/atom:  -7.5515 ε | E_hull:   0.0383 ε/atom
 5. Ni         | E/atom:  -5.2685 ε | E_hull:   2.3213 ε/atom
 6. Ni2P       | E/atom:  -4.9320 ε | E_hull:   1.4505 ε/atom
 7. Ni4P       | E/atom:  -6.0789 ε | E_hull:   0.7865 ε/atom
 8. Ni         | E/atom:  -6.2571 ε | E_hull:   1.3327 ε/atom
 9. Ni         | E/atom:  -3.5867 ε | E_hull:   4.0031 ε/atom
10. Ni         | E/atom:  -5.8764 ε | E_hull:   1.7134 ε/atom
11. Ni         | E/atom:  10.7466 ε | E_hull:  18.3364 ε/atom
12. Ni         | E/atom:  -2.9921 ε | E_hull:   4.5977 ε/atom
13. Ni         | E/atom:  -6.7571 ε | E_hull:   0.8327 ε/atom
14. Ni         | E/atom:  -6.4479 ε | E_hull:   1.1419 ε/atom
15. Ni         | E/atom:  -6.8673 ε | E_hull:   0.7225 ε/atom
16. Ni         | E/atom:  -5.4651 ε | E_hull:   2.1247 ε/atom
17. Ni         | E/atom:   3.4467 ε | E_hull:  11.0365 ε/atom
18. Ni         | E/atom:  -4.6192 ε | E_hull:   2.9706 ε/atom
19. Ni         | E/atom:  -5.7809 ε | E_hull:   1.8089 ε/atom
20. Ni         | E/atom:   2.2377 ε | E_hull:   9.8275 ε/atom

Total structures analyzed: 20
Total entries in phase diagram: 23