Tutorial 13: Machine Learning Potentials in Practice

cd $SCRATCH
mkdir -p tutorial13
cd tutorial13

$SCRATCH/tutorial13/
├── part1_mace_intro.ipynb      # You create this
├── part2_comparison.ipynb       # You create this
├── part3_validation.ipynb       # You create this
├── cu_bulk_energies.npz         # Generated by Part 1
├── cu_surface_energy.txt        # Generated by Part 1
└── cu_vacancy_energy.txt        # Generated by Part 1

/anvil/projects/x-chm250117/class_examples/tutorial13_completed/part2_mlp_vs_eam/
├── cu_melt_eam.lammpstrj
├── cu_md_eam.log
└── msd_eam.dat

import numpy as np
import matplotlib.pyplot as plt
from ase import Atoms
from ase.build import bulk, fcc111

# Test MACE import
from mace.calculators import mace_mp

print("="*60)
print("Tutorial 13: Machine Learning Potentials with MACE")
print("="*60)
print("\n✓ All imports successful!")

# Create bulk FCC copper with experimental lattice constant
a_exp = 3.615  # Angstroms (experimental value)
cu_bulk = bulk('Cu', crystalstructure='fcc', a=a_exp, cubic=True)

print("Created FCC copper unit cell:")
print(f"  Lattice constant: {a_exp} Å")
print(f"  Atoms: {len(cu_bulk)}")
print(f"  Cell volume: {cu_bulk.get_volume():.3f} ų")

# Load MACE-MP-0 foundation model
# 'medium' is a good balance of speed and accuracy
# First run downloads ~400 MB of model weights (cached for future use)

print("Loading MACE-MP-0 foundation model...")
print("(First run downloads model weights - please wait ~1 min)")

calc = mace_mp(model="medium", dispersion=False, default_dtype="float64")

print("✓ MACE calculator loaded!")

# Attach calculator and compute properties
cu_bulk.calc = calc

energy = cu_bulk.get_potential_energy()
forces = cu_bulk.get_forces()
stress = cu_bulk.get_stress()  # Voigt notation: [xx, yy, zz, yz, xz, xy]

print("--- Single Point Calculation ---")
print(f"Total energy: {energy:.6f} eV")
print(f"Energy per atom: {energy/len(cu_bulk):.6f} eV/atom")
print(f"Max force component: {np.max(np.abs(forces)):.6f} eV/Å")
print(f"  (Should be ~0 for perfect crystal)")

# Scan different lattice constants
print("Scanning lattice constants...")
print("(This takes ~2 minutes)\n")

a_values = np.linspace(3.50, 3.75, 15)
energies = []

for i, a in enumerate(a_values):
    cu_test = bulk('Cu', crystalstructure='fcc', a=a, cubic=True)
    cu_test.calc = calc
    e = cu_test.get_potential_energy() / len(cu_test)
    energies.append(e)
    print(f"  [{i+1:2d}/15] a = {a:.3f} Å  →  E = {e:.6f} eV/atom")

energies = np.array(energies)

# Find minimum
idx_min = np.argmin(energies)
a_optimal = a_values[idx_min]
e_cohesive = energies[idx_min]

print("\n" + "="*50)
print("RESULTS:")
print("="*50)
print(f"MACE optimal lattice constant: {a_optimal:.3f} Å")
print(f"Experimental lattice constant: {a_exp:.3f} Å")
print(f"Error: {(a_optimal - a_exp)/a_exp * 100:.2f}%")
print(f"\nMACE cohesive energy: {e_cohesive:.4f} eV/atom")
print(f"Experimental cohesive energy: -3.49 eV/atom")
print(f"Error: {(e_cohesive - (-3.49))/3.49 * 100:.1f}%")

# Save results
np.savez('cu_bulk_energies.npz', 
         a_values=a_values, 
         energies=energies,
         a_optimal=a_optimal,
         e_cohesive=e_cohesive)
print("✓ Results saved to cu_bulk_energies.npz")

# Plot
plt.figure(figsize=(10, 6))
plt.plot(a_values, energies, 'bo-', linewidth=2, markersize=8, label='MACE-MP-0')
plt.axvline(a_exp, color='red', linestyle='--', linewidth=2, label=f'Exp. a = {a_exp} Å')
plt.axvline(a_optimal, color='blue', linestyle=':', linewidth=2, label=f'MACE a = {a_optimal:.3f} Å')
plt.scatter([a_optimal], [e_cohesive], s=150, c='blue', zorder=5, 
            edgecolors='black', linewidths=2)

plt.xlabel('Lattice Constant (Å)', fontsize=12)
plt.ylabel('Energy (eV/atom)', fontsize=12)
plt.title('Copper Equation of State: MACE vs Experiment', fontsize=14, fontweight='bold')
plt.legend(fontsize=11)
plt.grid(alpha=0.3)

# Annotate error
error_pct = (a_optimal - a_exp)/a_exp * 100
plt.annotate(f'Error: {error_pct:.2f}%', 
             xy=(a_optimal, e_cohesive), xytext=(a_optimal + 0.05, e_cohesive + 0.02),
             fontsize=11, arrowprops=dict(arrowstyle='->', color='black'))

plt.tight_layout()
plt.savefig('cu_equation_of_state.png', dpi=150, bbox_inches='tight')
plt.show()

print("="*60)
print("Calculating Cu(111) Surface Energy")
print("="*60)

# Bulk reference energy (use our optimal lattice constant)
a = 3.615  # Use experimental for consistency
cu_bulk_ref = bulk('Cu', crystalstructure='fcc', a=a, cubic=True)
cu_bulk_ref.calc = calc
e_bulk_per_atom = cu_bulk_ref.get_potential_energy() / len(cu_bulk_ref)
print(f"\nBulk energy per atom: {e_bulk_per_atom:.6f} eV")

# Create Cu(111) slab: 3x3 supercell, 5 layers, 12 Å vacuum
# (Smaller than production to keep runtime reasonable)
print("\nCreating Cu(111) slab...")
slab = fcc111('Cu', size=(3, 3, 5), a=a, vacuum=12.0, periodic=True)
slab.calc = calc

n_atoms = len(slab)
print(f"  Atoms in slab: {n_atoms}")
print(f"  Cell dimensions: {slab.cell[0,0]:.2f} × {slab.cell[1,1]:.2f} × {slab.cell[2,2]:.2f} Å")

print("\nCalculating slab energy (this takes ~30 seconds)...")
e_slab = slab.get_potential_energy()
print(f"  Slab total energy: {e_slab:.4f} eV")

# Surface area (two surfaces in a slab!)
cell = slab.get_cell()
surface_area = np.linalg.norm(np.cross(cell[0], cell[1]))
print(f"  Surface area: {surface_area:.2f} Ų")

# Surface energy = (E_slab - N * E_bulk) / (2 * A)
gamma = (e_slab - n_atoms * e_bulk_per_atom) / (2 * surface_area)
gamma_Jm2 = gamma * 16.0218  # Convert eV/Ų to J/m²

print("\n" + "="*50)
print("SURFACE ENERGY RESULTS:")
print("="*50)
print(f"MACE γ(111): {gamma:.6f} eV/Ų = {gamma_Jm2:.3f} J/m²")
print(f"Experimental: ~1.79 J/m²")
print(f"DFT (PBE): ~1.24 J/m²")
print("\nNote: MACE is trained on DFT-PBE, so it matches DFT, not experiment.")

# Save result
with open('cu_surface_energy.txt', 'w') as f:
    f.write("Cu(111) Surface Energy\n")
    f.write("="*40 + "\n")
    f.write(f"MACE-MP-0 (medium): {gamma_Jm2:.4f} J/m²\n")
    f.write(f"Experimental: ~1.79 J/m²\n")
    f.write(f"DFT-PBE: ~1.24 J/m²\n")
print("\n✓ Saved to cu_surface_energy.txt")

print("="*60)
print("Calculating Vacancy Formation Energy")
print("="*60)

# Perfect bulk supercell (3x3x3 = 108 atoms)
print("\nCreating 3×3×3 supercell...")
supercell = bulk('Cu', crystalstructure='fcc', a=a, cubic=True) * (3, 3, 3)
supercell.calc = calc

e_perfect = supercell.get_potential_energy()
n_perfect = len(supercell)
e_per_atom = e_perfect / n_perfect

print(f"Perfect supercell:")
print(f"  Atoms: {n_perfect}")
print(f"  Energy: {e_perfect:.4f} eV")
print(f"  Energy/atom: {e_per_atom:.6f} eV")

# Create vacancy (remove one atom)
print("\nCreating vacancy (removing 1 atom)...")
vacancy = supercell.copy()
del vacancy[0]
vacancy.calc = calc

print("Calculating vacancy energy (this takes ~30 seconds)...")
e_vacancy = vacancy.get_potential_energy()
n_vacancy = len(vacancy)

print(f"Vacancy supercell:")
print(f"  Atoms: {n_vacancy}")
print(f"  Energy: {e_vacancy:.4f} eV")

# Vacancy formation energy
e_form = e_vacancy - (n_perfect - 1) * e_per_atom

print("\n" + "="*50)
print("VACANCY FORMATION ENERGY:")
print("="*50)
print(f"MACE E_f: {e_form:.4f} eV")
print(f"Experimental: 1.28 eV")
print(f"DFT (PBE): ~1.05 eV")

with open('cu_vacancy_energy.txt', 'w') as f:
    f.write("Cu Vacancy Formation Energy\n")
    f.write("="*40 + "\n")
    f.write(f"MACE-MP-0: {e_form:.4f} eV\n")
    f.write(f"Experimental: 1.28 eV\n")
    f.write(f"DFT-PBE: ~1.05 eV\n")
print("\n✓ Saved to cu_vacancy_energy.txt")

import numpy as np
import matplotlib.pyplot as plt
from ase.build import bulk
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md.langevin import Langevin
from ase import units
from mace.calculators import mace_mp
import time

print("="*60)
print("Part 2: MACE vs EAM Comparison")
print("="*60)

# Create small system for quick demo: 2x2x2 FCC = 32 atoms
# (Production would use 4x4x4 = 256 atoms, but that's too slow for class)
a = 3.615
atoms = bulk('Cu', crystalstructure='fcc', a=a, cubic=True) * (2, 2, 2)

print(f"System: {len(atoms)} atoms (2×2×2 supercell)")
print("Note: Small system for demo. Production runs use larger cells.")

# Load MACE
print("\nLoading MACE calculator...")
calc = mace_mp(model="medium", dispersion=False, default_dtype="float32")  # float32 for speed
atoms.calc = calc

# Initialize velocities at 500 K
temperature_K = 500
MaxwellBoltzmannDistribution(atoms, temperature_K=temperature_K)

# Remove center-of-mass motion
momenta = atoms.get_momenta()
momenta -= momenta.mean(axis=0)
atoms.set_momenta(momenta)

print(f"Target temperature: {temperature_K} K")

# Setup Langevin thermostat
timestep = 2.0 * units.fs
friction = 0.02 / units.fs

dyn = Langevin(atoms, timestep, temperature_K=temperature_K, friction=friction)

# Collect data
temps_mace = []
energies_mace = []
times_mace = []

def collect_data():
    t = dyn.get_time() / units.fs
    T = atoms.get_kinetic_energy() / (1.5 * len(atoms) * units.kB)
    E = atoms.get_potential_energy()
    temps_mace.append(T)
    energies_mace.append(E)
    times_mace.append(t)

dyn.attach(collect_data, interval=5)

# Run SHORT simulation: 100 steps = 200 fs
# This is just a performance demo - not enough for production statistics!
n_steps = 100

print(f"\nRunning {n_steps} MD steps ({n_steps * 2} fs)...")
print("(This takes ~1-2 minutes)\n")

start_time = time.time()
dyn.run(n_steps)
elapsed = time.time() - start_time

print(f"✓ Completed in {elapsed:.1f} seconds")
print(f"  Performance: {n_steps/elapsed:.1f} steps/second")
print(f"  Time per step: {elapsed/n_steps*1000:.1f} ms")

# Extrapolate to estimate practical timescales
ns_per_day = (n_steps * 2e-6) * (86400 / elapsed)
print(f"\n  Estimated throughput: {ns_per_day:.3f} ns/day")
print(f"  To simulate 1 ns would take: {1.0/ns_per_day:.1f} days")

# Convert to arrays
temps_mace = np.array(temps_mace)
energies_mace = np.array(energies_mace)
times_mace = np.array(times_mace)

print(f"\nMACE MD Statistics (short run):")
print(f"  Mean temperature: {np.mean(temps_mace):.1f} K")
print(f"  Std temperature: {np.std(temps_mace):.1f} K")

# Path to pre-run EAM results
eam_dir = "/anvil/projects/x-chm250117/class_examples/tutorial13_completed/part2_mlp_vs_eam"

def parse_lammps_log(logfile):
    """Extract thermo data from LAMMPS log file"""
    data = {'Step': [], 'Time': [], 'Temp': [], 'PotEng': [], 'KinEng': [], 'TotEng': [], 'Press': []}
    reading = False
    
    with open(logfile, 'r') as f:
        for line in f:
            if 'Step' in line and 'Temp' in line and 'PotEng' in line:
                reading = True
                continue
            if reading:
                if 'Loop time' in line or line.strip() == '':
                    reading = False
                    continue
                try:
                    parts = line.split()
                    if len(parts) >= 7:
                        data['Step'].append(int(parts[0]))
                        data['Time'].append(float(parts[1]))
                        data['Temp'].append(float(parts[2]))
                        data['PotEng'].append(float(parts[3]))
                        data['KinEng'].append(float(parts[4]))
                        data['TotEng'].append(float(parts[5]))
                        data['Press'].append(float(parts[6]))
                except (ValueError, IndexError):
                    pass
    
    return {k: np.array(v) for k, v in data.items()}

# Load EAM results
eam_data = parse_lammps_log(f'{eam_dir}/cu_md_eam.log')

print("EAM Simulation (pre-run with LAMMPS):")
print(f"  System: 256 atoms (4×4×4 supercell)")
print(f"  Steps recorded: {len(eam_data['Temp'])}")
print(f"  Simulation time: {eam_data['Time'][-1]*1000:.0f} fs")
print(f"  Mean temperature: {np.mean(eam_data['Temp']):.1f} K")
print(f"  Std temperature: {np.std(eam_data['Temp']):.1f} K")

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# --- Temperature Evolution ---
ax1 = axes[0, 0]
ax1.plot(times_mace, temps_mace, 'b-', alpha=0.8, linewidth=1.5, label=f'MACE (32 atoms, {len(times_mace)} pts)')
ax1.plot(eam_data['Time']*1000, eam_data['Temp'], 'r-', alpha=0.8, linewidth=1.5, label=f'EAM (256 atoms, {len(eam_data["Temp"])} pts)')
ax1.axhline(500, color='k', linestyle='--', alpha=0.5, label='Target (500 K)')
ax1.set_xlabel('Time (fs)', fontsize=11)
ax1.set_ylabel('Temperature (K)', fontsize=11)
ax1.set_title('Temperature Evolution', fontsize=12, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(alpha=0.3)

# --- Energy per Atom ---
ax2 = axes[0, 1]
ax2.plot(times_mace, energies_mace/32, 'b-', alpha=0.8, linewidth=1.5, label='MACE')
ax2.plot(eam_data['Time']*1000, eam_data['PotEng']/256, 'r-', alpha=0.8, linewidth=1.5, label='EAM')
ax2.set_xlabel('Time (fs)', fontsize=11)
ax2.set_ylabel('Potential Energy (eV/atom)', fontsize=11)
ax2.set_title('Potential Energy Evolution', fontsize=12, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(alpha=0.3)

# --- Temperature Distribution ---
ax3 = axes[1, 0]
ax3.hist(temps_mace, bins=15, alpha=0.6, color='blue', density=True, 
         label=f'MACE: {np.mean(temps_mace):.0f}±{np.std(temps_mace):.0f} K')
ax3.hist(eam_data['Temp'], bins=30, alpha=0.6, color='red', density=True,
         label=f'EAM: {np.mean(eam_data["Temp"]):.0f}±{np.std(eam_data["Temp"]):.0f} K')
ax3.axvline(500, color='k', linestyle='--', alpha=0.5)
ax3.set_xlabel('Temperature (K)', fontsize=11)
ax3.set_ylabel('Probability Density', fontsize=11)
ax3.set_title('Temperature Distribution', fontsize=12, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(alpha=0.3)

# --- Summary Table ---
ax4 = axes[1, 1]
ax4.axis('off')

summary_text = """
COMPARISON SUMMARY
══════════════════════════════════════════════════

                        MACE            EAM
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
System size            32 atoms       256 atoms
Simulation time        {} fs         {} fs
Data points            {}             {}
Mean Temperature       {:.0f} K          {:.0f} K

PERFORMANCE COMPARISON:
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
MACE:  {:.1f} steps/second (32 atoms)
EAM:   ~10,000+ steps/second (256 atoms)

→ EAM is roughly 100-1000× faster!

WHEN TO USE EACH:
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
EAM:   Large systems, long simulations, bulk metals
MACE:  Accuracy-critical work, validation, 
       multi-element systems, defects/surfaces
""".format(
    int(times_mace[-1]), int(eam_data['Time'][-1]*1000),
    len(temps_mace), len(eam_data['Temp']),
    np.mean(temps_mace), np.mean(eam_data['Temp']),
    n_steps/elapsed
)

ax4.text(0.05, 0.95, summary_text, transform=ax4.transAxes, fontsize=10,
         verticalalignment='top', fontfamily='monospace',
         bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8))

plt.tight_layout()
plt.savefig('mace_vs_eam_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

print("\n✓ Comparison plot saved to mace_vs_eam_comparison.png")

print("""
═══════════════════════════════════════════════════════════════════
                    KEY OBSERVATIONS
═══════════════════════════════════════════════════════════════════

1. TEMPERATURE CONTROL
   Both thermostats maintain ~500 K, but fluctuations differ:
   - Smaller systems (MACE demo) have larger T fluctuations
   - This is physics, not a bug: σ(T) ∝ 1/√N
   
2. POTENTIAL ENERGY
   MACE and EAM predict different absolute energies:
   - This is expected: different reference states
   - What matters: energy DIFFERENCES (barriers, formation energies)

3. COMPUTATIONAL COST — THE KEY DIFFERENCE
   ┌─────────────────────────────────────────────────────────┐
   │  MACE: ~50-100 steps/second (CPU)                       │
   │  EAM:  ~10,000+ steps/second (CPU)                      │
   │                                                          │
   │  For 1 ns of 256 atoms:                                  │
   │    MACE: ~1-3 days                                       │
   │    EAM:  ~5 minutes                                      │
   └─────────────────────────────────────────────────────────┘

4. WHEN IS THE MACE COST JUSTIFIED?
   
   ✓ USE MACE when:
     • You need DFT accuracy but can't afford full DFT
     • Your element doesn't have a good classical potential
     • You're studying defects, surfaces, or interfaces
     • You're validating a cheaper method
     • Multi-element systems (alloys, compounds)
   
   ✗ USE EAM/CLASSICAL when:
     • You need long timescales (ns-μs)
     • You need large systems (>10,000 atoms)
     • A validated classical potential exists
     • You're studying bulk thermodynamic properties

═══════════════════════════════════════════════════════════════════
""")

import numpy as np
import matplotlib.pyplot as plt

print("="*60)
print("Part 3: Validating MACE Predictions")
print("="*60)

# Load Part 1 results
data = np.load('cu_bulk_energies.npz')
a_values = data['a_values']
energies = data['energies']
a_optimal = float(data['a_optimal'])
e_cohesive = float(data['e_cohesive'])

print(f"\nLoaded lattice constant scan from Part 1")
print(f"  MACE optimal a: {a_optimal:.4f} Å")
print(f"  MACE cohesive energy: {e_cohesive:.4f} eV/atom")

# Experimental and DFT reference values for copper
experimental = {
    'a_lattice': 3.615,      # Å
    'E_cohesive': -3.49,     # eV/atom
    'E_vacancy': 1.28,       # eV
    'gamma_111': 1.79,       # J/m²
}

dft_pbe = {
    'a_lattice': 3.63,
    'E_cohesive': -3.48,
    'E_vacancy': 1.05,
    'gamma_111': 1.24,
}

# Load our MACE results
gamma_mace = None
with open('cu_surface_energy.txt', 'r') as f:
    for line in f:
        if 'MACE' in line and 'J/m²' in line:
            # Parse "MACE-MP-0 (medium): 1.2345 J/m²"
            gamma_mace = float(line.split(':')[1].strip().split()[0])

e_vac_mace = None
with open('cu_vacancy_energy.txt', 'r') as f:
    for line in f:
        if 'MACE' in line and 'eV' in line:
            e_vac_mace = float(line.split(':')[1].strip().split()[0])

print(f"Loaded MACE results:")
print(f"  Surface energy: {gamma_mace} J/m²")
print(f"  Vacancy energy: {e_vac_mace} eV")

mace_results = {
    'a_lattice': a_optimal,
    'E_cohesive': e_cohesive,
    'E_vacancy': e_vac_mace,
    'gamma_111': gamma_mace,
}

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

properties = ['a_lattice', 'E_cohesive', 'E_vacancy', 'gamma_111']
titles = ['Lattice Constant (Å)', 'Cohesive Energy (eV/atom)', 
          'Vacancy Formation Energy (eV)', 'Surface Energy γ(111) (J/m²)']
exp_vals = [3.615, 3.49, 1.28, 1.79]  # Note: using |E_coh| for plotting

for idx, (prop, title, exp_val) in enumerate(zip(properties, titles, exp_vals)):
    ax = axes[idx // 2, idx % 2]
    
    # Get values (handle sign for cohesive energy)
    if prop == 'E_cohesive':
        mace_val = abs(mace_results[prop])
        dft_val = abs(dft_pbe[prop])
    else:
        mace_val = mace_results[prop]
        dft_val = dft_pbe.get(prop, np.nan)
    
    x = np.arange(3)
    values = [mace_val, dft_val, exp_val]
    colors = ['steelblue', 'forestgreen', 'firebrick']
    labels = ['MACE-MP-0', 'DFT-PBE', 'Experiment']
    
    bars = ax.bar(x, values, color=colors, edgecolor='black', linewidth=1.5)
    ax.set_xticks(x)
    ax.set_xticklabels(labels, fontsize=11)
    ax.set_ylabel(title, fontsize=11)
    ax.set_title(title, fontsize=12, fontweight='bold')
    ax.grid(alpha=0.3, axis='y')
    
    # Add value labels on bars
    for bar, val in zip(bars, values):
        if val is not None and not np.isnan(val):
            ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.02*exp_val,
                   f'{val:.3f}', ha='center', fontsize=10, fontweight='bold')
    
    # Add error annotation for MACE vs experiment
    if mace_val is not None and not np.isnan(mace_val):
        error = (mace_val - exp_val) / exp_val * 100
        ax.annotate(f'Error: {error:+.1f}%', xy=(0, mace_val), 
                   xytext=(0.5, mace_val * 1.1),
                   fontsize=9, color='steelblue')

plt.suptitle('MACE Validation: Comparison to DFT and Experiment', 
             fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig('mace_validation_summary.png', dpi=150, bbox_inches='tight')
plt.show()

# Calculate status indicators
def get_status(mace_val, exp_val, threshold_good=0.05, threshold_ok=0.15):
    if mace_val is None:
        return "?"
    error = abs(mace_val - exp_val) / abs(exp_val)
    if error < threshold_good:
        return "✓"
    elif error < threshold_ok:
        return "~"
    else:
        return "✗"

print("""
╔══════════════════════════════════════════════════════════════════╗
║              MACE VALIDATION CHECKLIST FOR COPPER                 ║
╠══════════════════════════════════════════════════════════════════╣
║                                                                   ║
║  Property              MACE      DFT-PBE    Exp.     Status      ║
║  ─────────────────────────────────────────────────────────────── ║""")

print(f"║  Lattice constant (Å)  {a_optimal:.3f}     3.63      3.615    {get_status(a_optimal, 3.615, 0.01, 0.02)}          ║")
print(f"║  Cohesive energy (eV)  {e_cohesive:.3f}    -3.48     -3.49     {get_status(e_cohesive, -3.49, 0.03, 0.10)}          ║")
print(f"║  Vacancy E_f (eV)      {e_vac_mace:.3f}     1.05      1.28     {get_status(e_vac_mace, 1.28, 0.10, 0.25)}          ║")
print(f"║  Surface γ(111) (J/m²) {gamma_mace:.3f}     1.24      1.79     {get_status(gamma_mace, 1.79, 0.15, 0.35)}          ║")

print("""║                                                                   ║
╠══════════════════════════════════════════════════════════════════╣
║  INTERPRETATION:                                                  ║
║  • MACE matches DFT-PBE well (it was trained on PBE data!)       ║
║  • Both underestimate surface energy (known DFT-PBE limitation)  ║
║  • Vacancy energy close to DFT, but lower than experiment        ║
║  • For Cu bulk properties: MACE is reliable                      ║
║                                                                   ║
╠══════════════════════════════════════════════════════════════════╣
║  RECOMMENDATION:                                                  ║
║  ✓ Use MACE for Cu if you trust DFT-PBE for your application    ║
║  ✗ Don't expect better-than-DFT accuracy                        ║
║  ✓ Validate for YOUR specific property before production        ║
╚══════════════════════════════════════════════════════════════════╝
""")

print("""
═══════════════════════════════════════════════════════════════════
                    MLP DECISION FRAMEWORK
═══════════════════════════════════════════════════════════════════

QUESTION 1: Is my element in the training data?
├── Yes → Foundation model might work (validate!)
└── No  → You need custom training data

QUESTION 2: Is my system similar to training data?
├── Bulk crystal → Usually good
├── Surface/interface → Check carefully
├── Defects → Validate against DFT
└── Extreme conditions (high T/P) → Likely extrapolating!

QUESTION 3: What accuracy do I need?
├── Qualitative trends → Foundation model probably OK
├── Quantitative comparison to experiment → Validate first
└── Better than DFT → NOT POSSIBLE with MLP trained on DFT!

QUESTION 4: What timescale do I need?
├── < 100 ps → Maybe just use DFT
├── 100 ps - 10 ns → MLP sweet spot
└── > 10 ns → Consider classical potential

═══════════════════════════════════════════════════════════════════
                      THE GOLDEN RULE
═══════════════════════════════════════════════════════════════════

  "An MLP is only as good as its training data.
   MACE-MP-0 is trained on DFT-PBE, so it gives you DFT-PBE
   accuracy — no better, no worse, but 1000× faster."

═══════════════════════════════════════════════════════════════════
""")

Need to simulate copper?
        │
        ├── Is accuracy critical?
        │     ├── Yes → Use MACE, validate carefully
        │     └── No  → Use EAM (faster)
        │
        ├── Need > 1 ns simulation?
        │     └── Use EAM (MACE too slow)
        │
        └── Studying defects/surfaces?
              ├── Yes → MACE likely better
              └── No  → EAM is fine for bulk