Ready-to-Use Code Examples

#!/usr/bin/env python3
"""
Complete PCM Simulation Pipeline
Demonstrates all major simulation capabilities
"""

import numpy as np
import matplotlib.pyplot as plt
from pcm import JMAKModel, PCMDevice, Reliability
from pcm.utils import create_reset_pulse, create_set_pulse, save_results

# Initialize models
jmak = JMAKModel(params={
    'Ea': 1.8,   # eV
    'k0': 1e16,  # 1/s
    'n': 3.0     # Avrami exponent
})

device = PCMDevice(params={
    'R_amorphous': 1e6,    # Ohms
    'R_crystalline': 1e3,   # Ohms
    'C_thermal': 1e-15,     # J/K
    'R_thermal': 1e7,       # K/W
})

reliability = Reliability()

# 1. Crystallization Kinetics
print("=" * 50)
print("1. CRYSTALLIZATION KINETICS")
print("=" * 50)

temps = [400, 450, 500, 550, 600]  # K
time = np.logspace(-9, -3, 100)    # 1ns to 1ms

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

for T in temps:
    X = jmak.simulate_isothermal(T, time)
    axes[0].semilogx(time * 1e6, X, label=f'{T} K')
    
    # Find t50
    t50_idx = np.argmin(np.abs(X - 0.5))
    t50 = time[t50_idx]
    print(f"  T={T}K: t50 = {t50*1e6:.2f} μs")

axes[0].set_xlabel('Time (μs)')
axes[0].set_ylabel('Crystalline Fraction')
axes[0].set_title('Isothermal Crystallization')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Arrhenius plot
inv_T = 1000 / np.array(temps)
t50_values = []
for T in temps:
    X = jmak.simulate_isothermal(T, time)
    t50_idx = np.argmin(np.abs(X - 0.5))
    t50_values.append(time[t50_idx])

axes[1].semilogy(inv_T, t50_values, 'o-', color='#9333ea', linewidth=2)
axes[1].set_xlabel('1000/T (K⁻¹)')
axes[1].set_ylabel('t₅₀ (s)')
axes[1].set_title('Arrhenius Analysis')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('crystallization_kinetics.png', dpi=150, facecolor='black')

# 2. Pulse Programming
print("\n" + "=" * 50)
print("2. PULSE PROGRAMMING")
print("=" * 50)

# RESET pulse
V_reset, t_reset = create_reset_pulse(V_reset=3.0, t_width=50e-9)
result_reset = device.simulate_voltage_pulse(V_reset, t_reset, X0=1.0)

print(f"RESET Pulse:")
print(f"  Peak Temperature: {np.max(result_reset['temperature']):.0f} K")
print(f"  Final State: {'Amorphous' if result_reset['phase'][-1] < 0.5 else 'Crystalline'}")
print(f"  Energy: {np.trapz(V_reset**2/result_reset['resistance'], t_reset)*1e12:.2f} pJ")

# SET pulse
V_set, t_set = create_set_pulse(V_set=1.5, t_width=500e-9)
result_set = device.simulate_voltage_pulse(V_set, t_set, X0=0.0)

print(f"\nSET Pulse:")
print(f"  Peak Temperature: {np.max(result_set['temperature']):.0f} K")
print(f"  Final State: {'Amorphous' if result_set['phase'][-1] < 0.5 else 'Crystalline'}")
print(f"  Energy: {np.trapz(V_set**2/result_set['resistance'], t_set)*1e12:.2f} pJ")

# Plot pulse results
fig, axes = plt.subplots(2, 3, figsize=(15, 8))

# RESET plots
axes[0,0].plot(t_reset*1e9, V_reset, 'r-', linewidth=2)
axes[0,0].set_xlabel('Time (ns)')
axes[0,0].set_ylabel('Voltage (V)')
axes[0,0].set_title('RESET Pulse')

axes[0,1].plot(t_reset*1e9, result_reset['temperature'], 'r-', linewidth=2)
axes[0,1].axhline(900, color='white', linestyle='--', alpha=0.5)
axes[0,1].set_xlabel('Time (ns)')
axes[0,1].set_ylabel('Temperature (K)')

axes[0,2].semilogy(t_reset*1e9, result_reset['resistance'], 'r-', linewidth=2)
axes[0,2].set_xlabel('Time (ns)')
axes[0,2].set_ylabel('Resistance (Ω)')

# SET plots
axes[1,0].plot(t_set*1e9, V_set, 'b-', linewidth=2)
axes[1,0].set_xlabel('Time (ns)')
axes[1,0].set_ylabel('Voltage (V)')
axes[1,0].set_title('SET Pulse')

axes[1,1].plot(t_set*1e9, result_set['temperature'], 'b-', linewidth=2)
axes[1,1].axhline(450, color='white', linestyle='--', alpha=0.5)
axes[1,1].set_xlabel('Time (ns)')
axes[1,1].set_ylabel('Temperature (K)')

axes[1,2].plot(t_set*1e9, result_set['phase'], 'b-', linewidth=2)
axes[1,2].set_xlabel('Time (ns)')
axes[1,2].set_ylabel('Crystalline Fraction')

plt.tight_layout()
plt.savefig('pulse_programming.png', dpi=150, facecolor='black')

# 3. Reliability Analysis
print("\n" + "=" * 50)
print("3. RELIABILITY ANALYSIS")
print("=" * 50)

# Retention
retention_data = reliability.retention_arrhenius()
print(f"Data Retention at 85°C: {retention_data['t50_85C']:.1f} years")

# Endurance
endurance_data = reliability.endurance_weibull(n_samples=10000)
print(f"Endurance B10 Life: {endurance_data['B10']:.2e} cycles")
print(f"Endurance B50 Life: {endurance_data['B50']:.2e} cycles")

# Save all results
results = {
    'crystallization': {'temps': temps, 'time': time},
    'reset': result_reset,
    'set': result_set,
    'retention': retention_data,
    'endurance': endurance_data
}

save_results(results, 'pcm_simulation_results.json')
print("\n✓ Results saved to pcm_simulation_results.json")

plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

def monte_carlo_variability(n_samples=10000, variation=0.2):
    """
    Monte Carlo simulation for PCM variability analysis
    
    Parameters:
    -----------
    n_samples : int
        Number of Monte Carlo samples
    variation : float
        Coefficient of variation (σ/μ)
    """
    
    # Device parameters with variability
    params = {
        'R_off': 1e6,      # Amorphous resistance
        'R_on': 1e3,       # Crystalline resistance
        'V_th': 1.2,       # Threshold voltage
        'T_c': 450,        # Crystallization temperature
        'E_a': 1.8         # Activation energy
    }
    
    # Generate random variations
    results = {}
    
    for param, mean_value in params.items():
        # Log-normal distribution for resistances
        if 'R_' in param:
            sigma = np.log(1 + variation**2)**0.5
            mu = np.log(mean_value) - sigma**2/2
            values = np.random.lognormal(mu, sigma, n_samples)
        # Normal distribution for other parameters
        else:
            values = np.random.normal(mean_value, mean_value*variation, n_samples)
        
        results[param] = values
    
    # Calculate derived quantities
    results['R_ratio'] = results['R_off'] / results['R_on']
    results['retention_85C'] = 10 * np.exp(results['E_a'] * 11605 / (85 + 273))
    
    # Plotting
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    # R_off distribution
    axes[0,0].hist(results['R_off']/1e6, bins=50, color='#9333ea', alpha=0.7, edgecolor='white')
    axes[0,0].set_xlabel('R_off (MΩ)')
    axes[0,0].set_ylabel('Count')
    axes[0,0].set_title(f'Off Resistance (CV={variation:.0%})')
    
    # R_on distribution
    axes[0,1].hist(results['R_on']/1e3, bins=50, color='#ff6b6b', alpha=0.7, edgecolor='white')
    axes[0,1].set_xlabel('R_on (kΩ)')
    axes[0,1].set_ylabel('Count')
    axes[0,1].set_title('On Resistance')
    
    # R_ratio distribution
    axes[0,2].hist(np.log10(results['R_ratio']), bins=50, color='#4ecdc4', alpha=0.7, edgecolor='white')
    axes[0,2].set_xlabel('log₁₀(R_off/R_on)')
    axes[0,2].set_ylabel('Count')
    axes[0,2].set_title('Resistance Ratio')
    
    # V_th distribution
    axes[1,0].hist(results['V_th'], bins=50, color='#ffc300', alpha=0.7, edgecolor='white')
    axes[1,0].set_xlabel('V_th (V)')
    axes[1,0].set_ylabel('Count')
    axes[1,0].set_title('Threshold Voltage')
    
    # Retention scatter plot
    axes[1,1].scatter(results['E_a'], np.log10(results['retention_85C']), 
                     alpha=0.3, color='#9333ea', s=1)
    axes[1,1].set_xlabel('Activation Energy (eV)')
    axes[1,1].set_ylabel('log₁₀(Retention @ 85°C)')
    axes[1,1].set_title('Retention Correlation')
    
    # Cumulative distribution
    sorted_ratio = np.sort(results['R_ratio'])
    cdf = np.arange(1, len(sorted_ratio)+1) / len(sorted_ratio)
    axes[1,2].semilogx(sorted_ratio, cdf*100, linewidth=2, color='#9333ea')
    axes[1,2].grid(True, alpha=0.3)
    axes[1,2].set_xlabel('R_off/R_on')
    axes[1,2].set_ylabel('Cumulative Probability (%)')
    axes[1,2].set_title('CDF of Resistance Ratio')
    
    # Add percentile lines
    percentiles = [10, 50, 90]
    for p in percentiles:
        val = np.percentile(sorted_ratio, p)
        axes[1,2].axvline(val, linestyle='--', alpha=0.5, color='white')
        axes[1,2].text(val, p, f'P{p}', ha='right', color='white')
    
    plt.tight_layout()
    
    # Print statistics
    print("=" * 50)
    print("MONTE CARLO STATISTICS")
    print("=" * 50)
    
    for param in ['R_off', 'R_on', 'R_ratio', 'V_th']:
        mean = np.mean(results[param])
        std = np.std(results[param])
        cv = std/mean
        
        print(f"\n{param}:")
        print(f"  Mean: {mean:.2e}")
        print(f"  Std:  {std:.2e}")
        print(f"  CV:   {cv:.2%}")
        print(f"  P10:  {np.percentile(results[param], 10):.2e}")
        print(f"  P50:  {np.percentile(results[param], 50):.2e}")
        print(f"  P90:  {np.percentile(results[param], 90):.2e}")
    
    return results

# Run simulation
results = monte_carlo_variability(n_samples=10000, variation=0.2)
plt.show()

// PCM Web Simulator Class
class PCMSimulator {
    constructor() {
        this.params = {
            Ea: 1.8,  // eV
            k0: 1e16, // 1/s
            n: 3.0,   // Avrami exponent
            kb: 8.617333262e-5  // eV/K
        };
        
        this.device = {
            R_amorphous: 1e6,  // Ohms
            R_crystalline: 1e3, // Ohms
            C_thermal: 1e-15,   // J/K
            R_thermal: 1e7,     // K/W
            T_ambient: 300      // K
        };
    }
    
    // JMAK crystallization model
    jmakRate(T) {
        return this.params.k0 * Math.exp(-this.params.Ea / (this.params.kb * T));
    }
    
    simulateIsothermal(T, timeArray) {
        const k = this.jmakRate(T);
        return timeArray.map(t => 
            1 - Math.exp(-Math.pow(k * t, this.params.n))
        );
    }
    
    // Pulse simulation
    simulatePulse(voltage, duration, initialState = 0) {
        const dt = duration / 1000; // 1000 time steps
        const time = [];
        const temperature = [];
        const resistance = [];
        const phase = [];
        
        let T = this.device.T_ambient;
        let X = initialState;
        
        for (let i = 0; i < 1000; i++) {
            const t = i * dt;
            time.push(t);
            
            // Calculate resistance
            const R = this.device.R_crystalline * X + 
                     this.device.R_amorphous * (1 - X);
            resistance.push(R);
            
            // Joule heating
            const P = voltage * voltage / R;
            
            // Temperature evolution
            const dT = (P - (T - this.device.T_ambient) / this.device.R_thermal) * 
                      dt / this.device.C_thermal;
            T += dT;
            temperature.push(T);
            
            // Phase change
            if (T > 900 && voltage > 2.5) {
                // Amorphization (RESET)
                X = Math.max(0, X - 0.1);
            } else if (T > 450) {
                // Crystallization (SET)
                const k = this.jmakRate(T);
                const dX = this.params.n * k * dt * Math.pow(1 - X, this.params.n);
                X = Math.min(1, X + dX);
            }
            phase.push(X);
        }
        
        return { time, temperature, resistance, phase };
    }
    
    // Create interactive plot
    createPlot(elementId, data) {
        const traces = [
            {
                x: data.time.map(t => t * 1e9), // Convert to ns
                y: data.temperature,
                name: 'Temperature',
                yaxis: 'y',
                line: { color: '#ff6b6b' }
            },
            {
                x: data.time.map(t => t * 1e9),
                y: data.phase,
                name: 'Crystalline Fraction',
                yaxis: 'y2',
                line: { color: '#9333ea' }
            }
        ];
        
        const layout = {
            title: 'PCM Device Simulation',
            xaxis: { 
                title: 'Time (ns)',
                gridcolor: '#333'
            },
            yaxis: { 
                title: 'Temperature (K)',
                gridcolor: '#333',
                titlefont: { color: '#ff6b6b' },
                tickfont: { color: '#ff6b6b' }
            },
            yaxis2: {
                title: 'Crystalline Fraction',
                titlefont: { color: '#9333ea' },
                tickfont: { color: '#9333ea' },
                overlaying: 'y',
                side: 'right',
                range: [0, 1]
            },
            paper_bgcolor: 'rgba(0,0,0,0)',
            plot_bgcolor: 'rgba(0,0,0,0.1)',
            font: { color: '#e0e0e0' },
            showlegend: true
        };
        
        Plotly.newPlot(elementId, traces, layout, {responsive: true});
    }
}

// Usage Example
const simulator = new PCMSimulator();

// Simulate RESET pulse
const resetData = simulator.simulatePulse(3.0, 50e-9, 1.0);
simulator.createPlot('reset-plot', resetData);

// Simulate SET pulse
const setData = simulator.simulatePulse(1.5, 500e-9, 0.0);
simulator.createPlot('set-plot', setData);

// Interactive controls
document.getElementById('voltage-slider').addEventListener('input', function() {
    const voltage = parseFloat(this.value);
    const duration = parseFloat(document.getElementById('duration-slider').value) * 1e-9;
    
    const data = simulator.simulatePulse(voltage, duration);
    simulator.createPlot('interactive-plot', data);
    
    // Update display
    document.getElementById('voltage-display').textContent = voltage.toFixed(1) + ' V';
    document.getElementById('peak-temp').textContent = 
        Math.max(...data.temperature).toFixed(0) + ' K';
    document.getElementById('final-state').textContent = 
        data.phase[data.phase.length - 1] > 0.5 ? 'Crystalline' : 'Amorphous';
});