Chapter 5: Applications and Future

5.1 Industrial Applications

Scale of Industrial Catalysis

Catalysis drives the modern chemical industry, with 90% of all chemical processes relying on catalysts. The global catalyst market exceeded $35 billion in 2024 and continues to grow.

Key Industrial Processes

Code Example 1: Industrial Process Economics

"""
Economic analysis of industrial catalytic processes
Compares capital and operating costs across major industries
"""
import numpy as np
import matplotlib.pyplot as plt

# Industrial process data (approximate values)
processes = {
    'Ammonia\n(Haber-Bosch)': {
        'capex': 2000,  # $/ton capacity
        'opex': 300,    # $/ton product
        'catalyst_cost': 5,  # $/ton product
        'energy_frac': 0.75,  # Energy fraction of opex
        'catalyst_life': 15  # years
    },
    'FCC\n(Petroleum)': {
        'capex': 500,
        'opex': 50,
        'catalyst_cost': 2,
        'energy_frac': 0.40,
        'catalyst_life': 0.02  # Continuous replacement
    },
    'Methanol\nSynthesis': {
        'capex': 1500,
        'opex': 200,
        'catalyst_cost': 3,
        'energy_frac': 0.60,
        'catalyst_life': 4
    },
    'Polyethylene\n(Ziegler-Natta)': {
        'capex': 1200,
        'opex': 150,
        'catalyst_cost': 20,  # Higher activity required
        'energy_frac': 0.30,
        'catalyst_life': 0.001  # Single-use
    },
    'Sulfuric\nAcid': {
        'capex': 300,
        'opex': 30,
        'catalyst_cost': 0.5,
        'energy_frac': 0.50,  # Exothermic - energy producer
        'catalyst_life': 10
    }
}

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

# Plot 1: CAPEX comparison
ax1 = axes[0, 0]
names = list(processes.keys())
capex = [processes[p]['capex'] for p in names]
colors = plt.cm.viridis(np.linspace(0.2, 0.8, len(names)))
bars1 = ax1.bar(names, capex, color=colors, edgecolor='black')
ax1.set_ylabel('CAPEX ($/ton capacity)', fontsize=11)
ax1.set_title('Capital Investment Intensity', fontsize=12, fontweight='bold')
ax1.tick_params(axis='x', rotation=0)
for bar, val in zip(bars1, capex):
    ax1.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 30,
             f'${val}', ha='center', fontsize=9)

# Plot 2: Operating cost breakdown
ax2 = axes[0, 1]
x = np.arange(len(names))
width = 0.35

energy_costs = [processes[p]['opex'] * processes[p]['energy_frac'] for p in names]
other_costs = [processes[p]['opex'] * (1 - processes[p]['energy_frac']) for p in names]
catalyst_costs = [processes[p]['catalyst_cost'] for p in names]

ax2.bar(x - width/2, energy_costs, width, label='Energy', color='#ff6b6b')
ax2.bar(x + width/2, other_costs, width, label='Other OPEX', color='#4ecdc4')
ax2.bar(x + width/2, catalyst_costs, width, bottom=other_costs, label='Catalyst', color='#45b7d1')

ax2.set_ylabel('Operating Cost ($/ton product)', fontsize=11)
ax2.set_title('Operating Cost Breakdown', fontsize=12, fontweight='bold')
ax2.set_xticks(x)
ax2.set_xticklabels(names, fontsize=9)
ax2.legend()

# Plot 3: Catalyst cost vs lifetime
ax3 = axes[1, 0]
catalyst_life = [processes[p]['catalyst_life'] for p in names]
cat_costs = [processes[p]['catalyst_cost'] for p in names]

scatter = ax3.scatter(catalyst_life, cat_costs, s=200, c=colors, edgecolors='black')
for i, name in enumerate(names):
    ax3.annotate(name.replace('\n', ' '), (catalyst_life[i], cat_costs[i]),
                 xytext=(5, 5), textcoords='offset points', fontsize=8)

ax3.set_xlabel('Catalyst Lifetime (years)', fontsize=11)
ax3.set_ylabel('Catalyst Cost ($/ton product)', fontsize=11)
ax3.set_title('Catalyst Economics', fontsize=12, fontweight='bold')
ax3.set_xscale('log')
ax3.grid(alpha=0.3)

# Plot 4: Cost per ton analysis
ax4 = axes[1, 1]
# Annualized total cost (simplified: 10-year payback on CAPEX)
total_cost = []
for p in names:
    annual_capex = processes[p]['capex'] / 10
    opex = processes[p]['opex']
    total_cost.append(annual_capex + opex)

bars4 = ax4.barh(names, total_cost, color=colors, edgecolor='black')
ax4.set_xlabel('Total Annualized Cost ($/ton)', fontsize=11)
ax4.set_title('Total Cost Comparison', fontsize=12, fontweight='bold')
for bar, val in zip(bars4, total_cost):
    ax4.text(val + 5, bar.get_y() + bar.get_height()/2,
             f'${val:.0f}', va='center', fontsize=9)

plt.tight_layout()
plt.show()

print("\nIndustrial Process Economics Summary:")
for name, data in processes.items():
    print(f"\n{name.replace(chr(10), ' ')}:")
    print(f"  CAPEX: ${data['capex']}/ton capacity")
    print(f"  OPEX: ${data['opex']}/ton product")
    print(f"  Catalyst contribution: ${data['catalyst_cost']}/ton ({data['catalyst_cost']/data['opex']*100:.1f}%)")

5.2 Environmental Catalysis

Automotive Emission Control

The three-way catalyst (TWC) remains the cornerstone of automotive emission control, simultaneously converting CO, HC, and NOx:

Emerging Environmental Technologies

Code Example 2: Emission Control Analysis

"""
TWC light-off behavior and conversion efficiency analysis
Demonstrates the temperature-dependent performance of automotive catalysts
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def light_off_curve(T, T50, n):
    """
    Sigmoidal light-off curve
    T50: Temperature at 50% conversion
    n: Steepness parameter
    """
    return 100 / (1 + np.exp(-n * (T - T50) / T50))

# Temperature range
T = np.linspace(100, 500, 200)  # °C

# Different pollutants have different light-off temperatures
pollutants = {
    'CO': {'T50': 220, 'n': 15, 'color': 'red'},
    'HC (C₃H₆)': {'T50': 260, 'n': 12, 'color': 'green'},
    'NOx': {'T50': 280, 'n': 10, 'color': 'blue'}
}

# Create figure
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# Left: Light-off curves
ax1 = axes[0]
for name, params in pollutants.items():
    conversion = light_off_curve(T, params['T50'], params['n'])
    ax1.plot(T, conversion, color=params['color'], linewidth=2.5, label=name)
    ax1.axhline(y=50, color='gray', linestyle='--', alpha=0.5)
    ax1.axvline(x=params['T50'], color=params['color'], linestyle=':', alpha=0.5)

ax1.set_xlabel('Temperature (°C)', fontsize=12)
ax1.set_ylabel('Conversion (%)', fontsize=12)
ax1.set_title('TWC Light-off Curves', fontsize=13, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(alpha=0.3)
ax1.set_xlim(100, 500)
ax1.set_ylim(0, 105)
ax1.text(150, 55, 'T₅₀ line', fontsize=9, color='gray')

# Middle: Operating window (Lambda = 1 requirement)
ax2 = axes[1]
lambda_vals = np.linspace(0.95, 1.05, 100)

# Simplified conversion dependencies on lambda
def co_conv(lam):
    return 100 * np.exp(-50 * (lam - 1)**2) * (1 - 0.5 * np.maximum(0, lam - 1))

def hc_conv(lam):
    return 100 * np.exp(-40 * (lam - 1)**2)

def nox_conv(lam):
    return 100 * np.exp(-60 * (lam - 1)**2) * (1 + 0.3 * np.minimum(0, lam - 1))

ax2.plot(lambda_vals, co_conv(lambda_vals), 'r-', linewidth=2.5, label='CO')
ax2.plot(lambda_vals, hc_conv(lambda_vals), 'g-', linewidth=2.5, label='HC')
ax2.plot(lambda_vals, nox_conv(lambda_vals), 'b-', linewidth=2.5, label='NOx')

# Operating window
ax2.axvspan(0.995, 1.005, alpha=0.2, color='yellow', label='Operating window')
ax2.axvline(x=1.0, color='black', linestyle='--', alpha=0.7)

ax2.set_xlabel('Air/Fuel Ratio (λ)', fontsize=12)
ax2.set_ylabel('Conversion (%)', fontsize=12)
ax2.set_title('Lambda Window', fontsize=13, fontweight='bold')
ax2.legend(fontsize=9, loc='lower left')
ax2.grid(alpha=0.3)
ax2.set_xlim(0.95, 1.05)
ax2.set_ylim(0, 105)

# Right: Emission standards progression
ax3 = axes[2]
standards = ['Euro 1\n(1992)', 'Euro 3\n(2000)', 'Euro 5\n(2009)', 'Euro 6d\n(2020)', 'Euro 7\n(2025)']
co_limits = [2.72, 0.64, 0.50, 0.50, 0.30]  # g/km for gasoline
hc_limits = [0.97, 0.20, 0.10, 0.10, 0.05]
nox_limits = [0.97, 0.15, 0.06, 0.06, 0.03]

x = np.arange(len(standards))
width = 0.25

ax3.bar(x - width, co_limits, width, label='CO', color='red', alpha=0.8)
ax3.bar(x, [h*5 for h in hc_limits], width, label='HC (×5)', color='green', alpha=0.8)
ax3.bar(x + width, [n*5 for n in nox_limits], width, label='NOx (×5)', color='blue', alpha=0.8)

ax3.set_xlabel('Emission Standard', fontsize=12)
ax3.set_ylabel('Emission Limit (g/km)', fontsize=12)
ax3.set_title('Evolution of Emission Standards', fontsize=13, fontweight='bold')
ax3.set_xticks(x)
ax3.set_xticklabels(standards, fontsize=9)
ax3.legend(fontsize=9)
ax3.set_yscale('log')

plt.tight_layout()
plt.show()

print("\nEmission Control Analysis:")
print(f"• CO light-off (T₅₀): 220°C")
print(f"• HC light-off (T₅₀): 260°C")
print(f"• NOx light-off (T₅₀): 280°C")
print(f"• Operating window: λ = 0.995-1.005 (±0.5%)")
print(f"• Euro 7 reduction vs Euro 1: CO 89%, HC 95%, NOx 97%")

5.3 Green Hydrogen Production

Water Electrolysis Technologies

Electrocatalyst Design Principles

Efficient hydrogen evolution reaction (HER) and oxygen evolution reaction (OER) catalysts follow specific design principles:

• HER: Optimal hydrogen binding energy (Sabatier principle)
• OER: Balanced M-O and M-OH binding (scaling relations)
• Stability: Resistance to dissolution and oxidation
• Conductivity: High electronic conductivity for current flow

Code Example 3: Green Hydrogen Economics

"""
Techno-economic analysis of green hydrogen production
Compares different electrolysis technologies and cost projections
"""
import numpy as np
import matplotlib.pyplot as plt

# Current and projected costs ($/kg H2)
years = np.array([2020, 2022, 2024, 2026, 2028, 2030, 2035, 2040])

# Cost trajectories by technology
technologies = {
    'PEM Electrolysis': {
        'costs': [8.0, 6.5, 5.5, 4.5, 3.5, 2.5, 1.8, 1.5],
        'color': '#2196F3',
        'marker': 'o'
    },
    'Alkaline Electrolysis': {
        'costs': [6.0, 5.0, 4.2, 3.5, 2.8, 2.2, 1.6, 1.3],
        'color': '#4CAF50',
        'marker': 's'
    },
    'AEM Electrolysis': {
        'costs': [10.0, 7.5, 5.5, 4.0, 3.0, 2.0, 1.4, 1.1],
        'color': '#FF9800',
        'marker': '^'
    },
    'SOEC (High-T)': {
        'costs': [12.0, 9.0, 6.5, 4.5, 3.2, 2.2, 1.5, 1.0],
        'color': '#9C27B0',
        'marker': 'D'
    }
}

# Comparison benchmarks
gray_h2_cost = 1.5  # Steam methane reforming
blue_h2_cost = 2.5  # SMR + CCS

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

# Plot 1: Cost trajectory
ax1 = axes[0, 0]
for name, data in technologies.items():
    ax1.plot(years, data['costs'], marker=data['marker'], color=data['color'],
             linewidth=2, markersize=8, label=name)

ax1.axhline(y=gray_h2_cost, color='gray', linestyle='--', label='Gray H₂ (SMR)')
ax1.axhline(y=blue_h2_cost, color='steelblue', linestyle='--', label='Blue H₂ (SMR+CCS)')
ax1.axhline(y=2.0, color='green', linestyle=':', linewidth=2, label='DOE Target ($2/kg)')

ax1.set_xlabel('Year', fontsize=12)
ax1.set_ylabel('Hydrogen Cost ($/kg)', fontsize=12)
ax1.set_title('Green Hydrogen Cost Projections', fontsize=13, fontweight='bold')
ax1.legend(fontsize=9, loc='upper right')
ax1.grid(alpha=0.3)
ax1.set_xlim(2020, 2040)
ax1.set_ylim(0, 13)

# Plot 2: Cost breakdown (2024)
ax2 = axes[0, 1]
components = ['Electricity', 'CAPEX\n(Stack)', 'CAPEX\n(BoP)', 'O&M', 'Catalyst']
pem_breakdown = [2.5, 1.2, 0.8, 0.6, 0.4]  # $/kg
alk_breakdown = [2.3, 0.8, 0.6, 0.3, 0.2]

x = np.arange(len(components))
width = 0.35

bars1 = ax2.bar(x - width/2, pem_breakdown, width, label='PEM', color='#2196F3')
bars2 = ax2.bar(x + width/2, alk_breakdown, width, label='Alkaline', color='#4CAF50')

ax2.set_xlabel('Cost Component', fontsize=12)
ax2.set_ylabel('Cost Contribution ($/kg H₂)', fontsize=12)
ax2.set_title('Cost Breakdown (2024)', fontsize=13, fontweight='bold')
ax2.set_xticks(x)
ax2.set_xticklabels(components, fontsize=10)
ax2.legend()

# Plot 3: Efficiency comparison
ax3 = axes[1, 0]
techs = ['AWE', 'PEM', 'AEM', 'SOEC']
current_eff = [65, 75, 70, 85]
projected_eff = [75, 85, 82, 95]

x = np.arange(len(techs))
width = 0.35

ax3.bar(x - width/2, current_eff, width, label='2024', color='#90CAF9')
ax3.bar(x + width/2, projected_eff, width, label='2030 Target', color='#1565C0')

ax3.set_xlabel('Technology', fontsize=12)
ax3.set_ylabel('System Efficiency (%)', fontsize=12)
ax3.set_title('Electrolysis Efficiency', fontsize=13, fontweight='bold')
ax3.set_xticks(x)
ax3.set_xticklabels(techs)
ax3.legend()
ax3.set_ylim(0, 100)
ax3.axhline(y=80, color='green', linestyle='--', alpha=0.5, label='Target')

# Plot 4: Catalyst cost reduction impact
ax4 = axes[1, 1]
catalyst_reduction = np.linspace(0, 90, 50)  # % reduction in Ir/Pt loading

# Impact on system cost
def cost_impact(reduction_pct, base_catalyst_cost=0.4, total_cost=5.5):
    catalyst_savings = base_catalyst_cost * reduction_pct / 100
    return total_cost - catalyst_savings

pem_cost_traj = [cost_impact(r, 0.4, 5.5) for r in catalyst_reduction]

ax4.plot(catalyst_reduction, pem_cost_traj, 'b-', linewidth=2.5, label='PEM System Cost')
ax4.fill_between(catalyst_reduction, pem_cost_traj, 5.5, alpha=0.3, color='green')
ax4.axhline(y=2.0, color='red', linestyle='--', label='$2/kg target')

# Annotate milestones
ax4.annotate('50% Ir reduction\n(Toyota/Panasonic 2025)',
             xy=(50, cost_impact(50)), xytext=(60, 4.5),
             arrowprops=dict(arrowstyle='->', color='black'),
             fontsize=9)

ax4.set_xlabel('Catalyst Loading Reduction (%)', fontsize=12)
ax4.set_ylabel('H₂ Production Cost ($/kg)', fontsize=12)
ax4.set_title('Impact of Catalyst Innovation', fontsize=13, fontweight='bold')
ax4.legend(fontsize=9)
ax4.grid(alpha=0.3)

plt.tight_layout()
plt.show()

print("\nGreen Hydrogen Economics Summary:")
print(f"• 2024 PEM cost: $5.5/kg (Electricity 45%, CAPEX 36%, O&M 11%, Catalyst 7%)")
print(f"• 2030 target: <$2/kg")
print(f"• Efficiency gains needed: 75% → 85%")
print(f"• Catalyst cost reduction potential: 50-90% through SAC and non-PGM designs")

5.4 CO₂ Utilization

CO₂ Conversion Pathways

Electrochemical CO₂ Reduction

Code Example 4: CO₂ Conversion Economics

"""
Techno-economic analysis of CO₂ conversion pathways
Compares different products and their market potential
"""
import numpy as np
import matplotlib.pyplot as plt

# CO2 conversion products
products = {
    'Methanol': {
        'market_size': 100,  # Billion $/year
        'price': 450,  # $/ton
        'co2_utilization': 1.37,  # ton CO2/ton product
        'energy_req': 10.5,  # MWh/ton product
        'trl': 9  # Technology Readiness Level
    },
    'Formic Acid': {
        'market_size': 1,
        'price': 800,
        'co2_utilization': 0.96,
        'energy_req': 5.0,
        'trl': 7
    },
    'CO (Syngas)': {
        'market_size': 10,
        'price': 150,
        'co2_utilization': 1.57,
        'energy_req': 3.5,
        'trl': 8
    },
    'Ethylene': {
        'market_size': 200,
        'price': 1200,
        'co2_utilization': 3.14,
        'energy_req': 25.0,
        'trl': 5
    },
    'Jet Fuel': {
        'market_size': 300,
        'price': 800,
        'co2_utilization': 3.5,
        'energy_req': 30.0,
        'trl': 4
    }
}

# Calculate derived metrics
for name, data in products.items():
    # CO2 cost at $50/ton
    data['co2_cost'] = 50 * data['co2_utilization']
    # Energy cost at $50/MWh
    data['energy_cost'] = 50 * data['energy_req']
    # Production cost (simplified)
    data['prod_cost'] = data['co2_cost'] + data['energy_cost'] + 100  # +100 for CAPEX/OPEX
    # Margin
    data['margin'] = data['price'] - data['prod_cost']
    # CO2 abatement cost
    data['abatement_cost'] = max(0, -data['margin']) / data['co2_utilization'] if data['margin'] < 0 else 0

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

names = list(products.keys())
colors = plt.cm.Set2(np.linspace(0, 1, len(names)))

# Plot 1: Market size vs CO2 utilization
ax1 = axes[0, 0]
market_sizes = [products[p]['market_size'] for p in names]
co2_util = [products[p]['co2_utilization'] for p in names]
trls = [products[p]['trl'] for p in names]

scatter = ax1.scatter(co2_util, market_sizes, s=[t**2 * 15 for t in trls],
                      c=colors, edgecolors='black', alpha=0.7)

for i, name in enumerate(names):
    ax1.annotate(name, (co2_util[i], market_sizes[i]),
                 xytext=(5, 5), textcoords='offset points', fontsize=9)

ax1.set_xlabel('CO₂ Utilization (ton CO₂/ton product)', fontsize=12)
ax1.set_ylabel('Market Size (Billion $/year)', fontsize=12)
ax1.set_title('Market Opportunity vs CO₂ Intensity', fontsize=13, fontweight='bold')
ax1.set_yscale('log')
ax1.grid(alpha=0.3)
ax1.text(0.05, 0.95, 'Bubble size = TRL', transform=ax1.transAxes, fontsize=9)

# Plot 2: Cost breakdown
ax2 = axes[0, 1]
x = np.arange(len(names))
width = 0.6

co2_costs = [products[p]['co2_cost'] for p in names]
energy_costs = [products[p]['energy_cost'] for p in names]
other_costs = [100] * len(names)
prices = [products[p]['price'] for p in names]

ax2.bar(x, co2_costs, width, label='CO₂', color='#ff6b6b')
ax2.bar(x, energy_costs, width, bottom=co2_costs, label='Energy', color='#ffd93d')
ax2.bar(x, other_costs, width, bottom=[c+e for c, e in zip(co2_costs, energy_costs)],
        label='Other', color='#6bcb77')
ax2.scatter(x, prices, s=100, c='black', marker='_', linewidths=3, label='Market price', zorder=5)

ax2.set_xlabel('Product', fontsize=12)
ax2.set_ylabel('Cost / Price ($/ton)', fontsize=12)
ax2.set_title('Production Cost vs Market Price', fontsize=13, fontweight='bold')
ax2.set_xticks(x)
ax2.set_xticklabels(names, rotation=45, ha='right', fontsize=9)
ax2.legend(fontsize=9)

# Plot 3: Technology readiness
ax3 = axes[1, 0]
trls = [products[p]['trl'] for p in names]
colors_trl = ['green' if t >= 7 else 'yellow' if t >= 5 else 'red' for t in trls]

bars = ax3.barh(names, trls, color=colors_trl, edgecolor='black')
ax3.axvline(x=7, color='green', linestyle='--', alpha=0.5, label='Commercial ready')
ax3.axvline(x=5, color='orange', linestyle='--', alpha=0.5, label='Demonstration')

ax3.set_xlabel('Technology Readiness Level', fontsize=12)
ax3.set_title('Technology Maturity', fontsize=13, fontweight='bold')
ax3.set_xlim(0, 10)
ax3.legend(fontsize=9)

for bar, trl in zip(bars, trls):
    ax3.text(trl + 0.1, bar.get_y() + bar.get_height()/2,
             f'TRL {trl}', va='center', fontsize=10)

# Plot 4: CO2 abatement potential
ax4 = axes[1, 1]
# Scenario: 1 Gt CO2/year utilization by 2040
target_co2 = 1000  # Mt CO2/year

potential_utilization = {}
for name, data in products.items():
    # Max utilization based on current market × 3 growth
    max_output = data['market_size'] * 3 * 1e9 / data['price']  # tons/year
    co2_captured = max_output * data['co2_utilization'] / 1e6  # Mt CO2/year
    potential_utilization[name] = min(co2_captured, 500)  # Cap at 500 Mt

util_values = [potential_utilization[p] for p in names]
ax4.barh(names, util_values, color=colors, edgecolor='black')
ax4.axvline(x=100, color='red', linestyle='--', linewidth=2, label='100 Mt CO₂/year')
ax4.set_xlabel('CO₂ Utilization Potential (Mt/year)', fontsize=12)
ax4.set_title('CO₂ Capture Potential (2040)', fontsize=13, fontweight='bold')
ax4.legend(fontsize=9)

plt.tight_layout()
plt.show()

print("\nCO₂ Conversion Analysis:")
for name, data in products.items():
    status = "Profitable" if data['margin'] > 0 else f"Gap: ${-data['margin']:.0f}/ton"
    print(f"\n{name}:")
    print(f"  Production cost: ${data['prod_cost']:.0f}/ton")
    print(f"  Market price: ${data['price']}/ton")
    print(f"  Status: {status}")
    print(f"  TRL: {data['trl']}")

5.5 AI-Driven Catalyst Design

ML Applications in Catalysis

Code Example 5: ML-Assisted Catalyst Screening

"""
ML-assisted catalyst screening using descriptor-based approach
Demonstrates feature engineering and model training for catalyst activity prediction
"""
import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error, r2_score
from sklearn.preprocessing import StandardScaler

# Simulated catalyst dataset
np.random.seed(42)
n_samples = 500

# Features (catalyst descriptors)
features = {
    'd_band_center': np.random.uniform(-4, 0, n_samples),  # eV
    'work_function': np.random.uniform(4, 6, n_samples),   # eV
    'coordination_number': np.random.randint(4, 12, n_samples),
    'electronegativity': np.random.uniform(1.5, 2.5, n_samples),
    'atomic_radius': np.random.uniform(1.2, 2.0, n_samples),  # Angstrom
    'surface_energy': np.random.uniform(1, 3, n_samples),     # J/m²
}

X = np.column_stack(list(features.values()))
feature_names = list(features.keys())

# Target: HER activity (simulated based on Sabatier principle)
# Activity peaks at optimal d-band center (~-2 eV)
def activity_function(d_band, work_fn, coord_num, eneg, radius, surf_e):
    # Volcano-like dependence on d-band
    d_band_term = -0.5 * (d_band + 2)**2
    # Work function contribution
    wf_term = -0.1 * (work_fn - 5)**2
    # Coordination effect
    coord_term = 0.05 * coord_num
    # Noise
    noise = np.random.normal(0, 0.3)
    return 2.0 + d_band_term + wf_term + coord_term + noise

y = np.array([activity_function(*X[i]) for i in range(n_samples)])
y = np.clip(y, 0, 3)  # Log(i0) in reasonable range

# Split data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Scale features
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# Train model
model = GradientBoostingRegressor(n_estimators=100, max_depth=4, random_state=42)
model.fit(X_train_scaled, y_train)

# Predictions
y_pred_train = model.predict(X_train_scaled)
y_pred_test = model.predict(X_test_scaled)

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

# Plot 1: Parity plot
ax1 = axes[0, 0]
ax1.scatter(y_train, y_pred_train, alpha=0.5, s=30, c='blue', label='Train')
ax1.scatter(y_test, y_pred_test, alpha=0.7, s=50, c='red', label='Test')
ax1.plot([0, 3], [0, 3], 'k--', linewidth=2)

mae_test = mean_absolute_error(y_test, y_pred_test)
r2_test = r2_score(y_test, y_pred_test)

ax1.set_xlabel('Actual Activity (log i₀)', fontsize=12)
ax1.set_ylabel('Predicted Activity (log i₀)', fontsize=12)
ax1.set_title('ML Model Performance', fontsize=13, fontweight='bold')
ax1.legend()
ax1.text(0.05, 0.95, f'MAE = {mae_test:.3f}\nR² = {r2_test:.3f}',
         transform=ax1.transAxes, fontsize=11, verticalalignment='top',
         bbox=dict(boxstyle='round', facecolor='wheat'))

# Plot 2: Feature importance
ax2 = axes[0, 1]
importance = model.feature_importances_
sorted_idx = np.argsort(importance)
ax2.barh([feature_names[i] for i in sorted_idx], importance[sorted_idx],
         color='steelblue', edgecolor='black')
ax2.set_xlabel('Feature Importance', fontsize=12)
ax2.set_title('Catalyst Descriptor Importance', fontsize=13, fontweight='bold')

# Plot 3: d-band center vs activity (volcano plot)
ax3 = axes[1, 0]
d_band_all = X[:, 0]
scatter = ax3.scatter(d_band_all, y, c=X[:, 2], cmap='viridis', alpha=0.6, s=30)
plt.colorbar(scatter, ax=ax3, label='Coordination Number')

# Fit polynomial for volcano
d_band_sorted = np.linspace(-4, 0, 100)
poly_fit = np.polyfit(d_band_all, y, 2)
y_fit = np.polyval(poly_fit, d_band_sorted)
ax3.plot(d_band_sorted, y_fit, 'r-', linewidth=2, label='Volcano fit')

ax3.axvline(x=-2, color='green', linestyle='--', alpha=0.5, label='Optimal d-band')
ax3.set_xlabel('d-band Center (eV)', fontsize=12)
ax3.set_ylabel('HER Activity (log i₀)', fontsize=12)
ax3.set_title('Sabatier Volcano Relationship', fontsize=13, fontweight='bold')
ax3.legend(fontsize=9)

# Plot 4: Screening results
ax4 = axes[1, 1]
# Generate new candidates
n_candidates = 1000
X_new = np.column_stack([
    np.random.uniform(-4, 0, n_candidates),
    np.random.uniform(4, 6, n_candidates),
    np.random.randint(4, 12, n_candidates),
    np.random.uniform(1.5, 2.5, n_candidates),
    np.random.uniform(1.2, 2.0, n_candidates),
    np.random.uniform(1, 3, n_candidates),
])
X_new_scaled = scaler.transform(X_new)
y_new_pred = model.predict(X_new_scaled)

# Histogram of predicted activities
ax4.hist(y_new_pred, bins=30, color='steelblue', edgecolor='black', alpha=0.7)
ax4.axvline(x=np.percentile(y_new_pred, 90), color='red', linestyle='--',
            linewidth=2, label=f'Top 10% threshold: {np.percentile(y_new_pred, 90):.2f}')

top_10_pct = np.sum(y_new_pred > np.percentile(y_new_pred, 90))
ax4.set_xlabel('Predicted Activity (log i₀)', fontsize=12)
ax4.set_ylabel('Count', fontsize=12)
ax4.set_title('Virtual Screening Results', fontsize=13, fontweight='bold')
ax4.legend(fontsize=9)
ax4.text(0.05, 0.95, f'Candidates screened: {n_candidates}\nTop performers: {top_10_pct}',
         transform=ax4.transAxes, fontsize=11, verticalalignment='top')

plt.tight_layout()
plt.show()

print("\nML Catalyst Screening Results:")
print(f"• Training samples: {len(y_train)}")
print(f"• Test MAE: {mae_test:.3f} log(i₀)")
print(f"• Test R²: {r2_test:.3f}")
print(f"• Most important feature: {feature_names[np.argmax(importance)]}")
print(f"• Candidates screened: {n_candidates}")
print(f"• Top 10% activity threshold: {np.percentile(y_new_pred, 90):.2f}")

5.6 Artificial Enzymes and De Novo Design

Artificial Enzyme Categories

Key Achievements

• C-Si Bond Formation: Cytochrome P450 variants creating silicon-carbon bonds not found in nature
• Asymmetric Catalysis: >99% ee in cyclopropanation reactions
• Stability: Designed enzymes stable at >100°C
• Activity: kcat/KM exceeding natural enzymes by 100×

5.7 Future Directions

Emerging Technologies

Self-Healing Catalysts

Catalysts that can autonomously regenerate their active sites during operation:

• Mechanism: Dynamic atom migration under reaction conditions
• Example: Pt-Ni alloys that restructure to maintain activity
• Benefit: 10× extended catalyst lifetime

Quantum Effects in Catalysis

Exploiting quantum tunneling and coherence for enhanced selectivity:

• Proton tunneling: Accelerated H-transfer in enzymes
• Electron tunneling: Long-range charge transfer
• Quantum coherence: Enhanced selectivity pathways

Code Example 6: Technology Roadmap

"""
Technology roadmap visualization for future catalyst development
Shows timelines and expected breakthroughs
"""
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

# Technology roadmap data
technologies = {
    'Single-Atom Catalysts': {
        'start': 2020, 'current_trl': 7, 'target_trl': 9, 'target_year': 2028,
        'category': 'Materials'
    },
    'MOF Catalysts': {
        'start': 2015, 'current_trl': 6, 'target_trl': 9, 'target_year': 2030,
        'category': 'Materials'
    },
    'AI Catalyst Design': {
        'start': 2022, 'current_trl': 5, 'target_trl': 8, 'target_year': 2030,
        'category': 'Digital'
    },
    'Self-Healing Catalysts': {
        'start': 2024, 'current_trl': 3, 'target_trl': 7, 'target_year': 2035,
        'category': 'Emerging'
    },
    'Artificial Enzymes': {
        'start': 2018, 'current_trl': 4, 'target_trl': 8, 'target_year': 2032,
        'category': 'Bio-inspired'
    },
    'Green H₂ (<$2/kg)': {
        'start': 2020, 'current_trl': 6, 'target_trl': 9, 'target_year': 2030,
        'category': 'Energy'
    },
    'CO₂ to Fuels': {
        'start': 2018, 'current_trl': 5, 'target_trl': 8, 'target_year': 2035,
        'category': 'Energy'
    },
    'Quantum Catalysis': {
        'start': 2025, 'current_trl': 2, 'target_trl': 5, 'target_year': 2040,
        'category': 'Emerging'
    }
}

# Category colors
category_colors = {
    'Materials': '#2196F3',
    'Digital': '#9C27B0',
    'Emerging': '#FF5722',
    'Bio-inspired': '#4CAF50',
    'Energy': '#FFC107'
}

# Create figure
fig, axes = plt.subplots(1, 2, figsize=(16, 8))

# Left: Timeline (Gantt-like)
ax1 = axes[0]
y_positions = list(range(len(technologies)))
tech_names = list(technologies.keys())

for i, (name, data) in enumerate(technologies.items()):
    color = category_colors[data['category']]

    # Current progress
    ax1.barh(i, 2026 - data['start'], left=data['start'], height=0.4,
             color=color, alpha=0.7, edgecolor='black')

    # Future projection
    ax1.barh(i, data['target_year'] - 2026, left=2026, height=0.4,
             color=color, alpha=0.3, hatch='//', edgecolor='black')

    # TRL markers
    ax1.scatter(2026, i, s=100, c='white', edgecolors='black', zorder=5)
    ax1.text(2026, i, str(data['current_trl']), ha='center', va='center', fontsize=8)

ax1.set_yticks(y_positions)
ax1.set_yticklabels(tech_names, fontsize=10)
ax1.set_xlabel('Year', fontsize=12)
ax1.set_title('Catalyst Technology Roadmap', fontsize=14, fontweight='bold')
ax1.axvline(x=2026, color='red', linestyle='--', linewidth=2, label='Current (2026)')
ax1.set_xlim(2015, 2042)
ax1.legend(loc='lower right')

# Legend for categories
legend_patches = [mpatches.Patch(color=c, label=cat) for cat, c in category_colors.items()]
ax1.legend(handles=legend_patches, loc='upper left', fontsize=9)

# Right: TRL progression
ax2 = axes[1]

for name, data in technologies.items():
    color = category_colors[data['category']]
    years = [data['start'], 2026, data['target_year']]
    trls = [1, data['current_trl'], data['target_trl']]
    ax2.plot(years, trls, 'o-', color=color, linewidth=2, markersize=8, label=name)

ax2.set_xlabel('Year', fontsize=12)
ax2.set_ylabel('Technology Readiness Level', fontsize=12)
ax2.set_title('TRL Progression', fontsize=14, fontweight='bold')
ax2.set_xlim(2015, 2042)
ax2.set_ylim(0, 10)
ax2.axhline(y=7, color='green', linestyle='--', alpha=0.5, label='Commercial ready')
ax2.axhline(y=9, color='darkgreen', linestyle='--', alpha=0.5, label='Mature')
ax2.axvline(x=2026, color='red', linestyle='--', alpha=0.5)
ax2.legend(fontsize=8, loc='lower right', ncol=2)
ax2.grid(alpha=0.3)

plt.tight_layout()
plt.show()

# Summary table
print("\nCatalyst Technology Roadmap Summary:")
print("-" * 70)
print(f"{'Technology':<25} {'Current TRL':>12} {'Target TRL':>12} {'Target Year':>12}")
print("-" * 70)
for name, data in technologies.items():
    print(f"{name:<25} {data['current_trl']:>12} {data['target_trl']:>12} {data['target_year']:>12}")

print("\n2030 Targets:")
print("• Green hydrogen: <$2/kg (current: $4-6/kg)")
print("• SAC commercial deployment in fuel cells")
print("• AI-designed catalysts in production")
print("• MOF-based CO₂ capture at Mt/year scale")

5.8 Chapter Summary

Exercises

Exercise 1: Economic Analysis

Calculate the break-even electricity price for green hydrogen at $2/kg assuming 75% system efficiency and $500/kW electrolyzer cost (10-year lifetime, 8000 hours/year operation).

Exercise 2: CO₂ Conversion

A CO₂-to-methanol plant produces 1000 tons/day with 80% carbon efficiency. Calculate annual CO₂ utilization and compare with a 500 MW coal power plant's emissions.

Exercise 3: ML Catalyst Design

Design a feature set for predicting oxygen evolution reaction (OER) activity. What descriptors would you include beyond d-band center and work function?