Pycalphad

Use pycalphad when:

• Calculating phase diagrams (binary, ternary, or multicomponent systems)
• Determining phase equilibria at specific temperature/composition/pressure conditions
• Computing thermodynamic properties from CALPHAD databases
• Analyzing metastability and driving forces for phase transformations
• Modeling systems with ordered phases or charged species
• Automating thermodynamic calculations for materials design
• Creating custom thermodynamic models or property frameworks
• Visualizing phase stability regions and property maps

Don't use when:

• You need first-principles calculations (use VASP, Quantum ESPRESSO, GPAW instead)
• You need molecular dynamics simulations (use ASE, LAMMPS instead)
• Working with kinetic models (pycalphad is for equilibrium thermodynamics)

Core Concepts

1. Database (TDB Files)

Thermodynamic databases in Thermo-Calc format containing:

• Phase definitions and constituents
• Gibbs energy models for each phase
• Model parameters fitted to experimental data
• Temperature, composition, and pressure ranges

2. Components and Species

• Components: Independent chemical elements (e.g., Al, Fe, Ni)
• Species: Charged or molecular entities in ionic systems

3. Phases

Distinct thermodynamic phases with specific crystal structures and energy models:

• Liquid, FCC, BCC, HCP, intermetallic compounds, etc.
• Each phase has a sublattice model defining site occupancy

4. State Variables

Conditions that define the system state:

• Temperature (T)
• Pressure (P)
• Composition (mole fractions: X_i, Y_i, etc.)
• Chemical potential (MU)

5. Calculate vs Equilibrium

• calculate(): Evaluate properties at specified conditions (no equilibrium minimization)
• equilibrium(): Find stable phase assemblage by minimizing Gibbs energy

Quick Reference

Installation

# via pip
pip install pycalphad

# via conda
conda install -c conda-forge pycalphad

# development version
pip install git+https://github.com/pycalphad/pycalphad.git

Dependencies: Python 3.9+, numpy, scipy, xarray, sympy, matplotlib

Common Workflows

1. Binary Phase Diagram

from pycalphad import Database, binplot
import matplotlib.pyplot as plt

# Load thermodynamic database
db = Database('alzn_mey.tdb')

# Define components and phases
comps = ['AL', 'ZN', 'VA']  # VA = vacancy for gas phase
phases = ['LIQUID', 'FCC_A1', 'HCP_A3']

# Plot isobaric (constant pressure) binary diagram
fig = plt.figure(figsize=(8, 6))
binplot(db, comps, phases,
        conditions={'P': 101325, 'T': (300, 1000, 10), 'X(ZN)': (0, 1, 0.01)})
plt.xlabel('X(ZN)')
plt.ylabel('Temperature (K)')
plt.title('Al-Zn Binary Phase Diagram')
plt.savefig('al-zn-diagram.png')

2. Equilibrium Calculation at Fixed Conditions

from pycalphad import Database, equilibrium, variables as v
import numpy as np

# Load database
db = Database('nist_ni_al.tdb')

# Define system
comps = ['NI', 'AL', 'VA']
phases = ['LIQUID', 'FCC_L12', 'BCC_B2']

# Calculate equilibrium at 1500K, 1 atm, 50 at% Al
result = equilibrium(db, comps, phases,
                    {v.T: 1500, v.P: 101325, v.X('AL'): 0.5})

# Access results
print("Stable phases:", result.Phase.values.squeeze())
print("Phase fractions:", result.NP.values.squeeze())
print("Compositions:", result.X.values.squeeze())
print("Gibbs energy:", result.GM.values.squeeze())

3. Property Map (T-X Diagram with Property Overlay)

from pycalphad import Database, equilibrium, variables as v
import numpy as np
import matplotlib.pyplot as plt

db = Database('crfe.tdb')
comps = ['CR', 'FE', 'VA']
phases = ['LIQUID', 'BCC_A2', 'SIGMA']

# Create T-X grid
temps = np.linspace(1000, 2000, 50)
x_cr = np.linspace(0, 1, 50)
T_grid, X_grid = np.meshgrid(temps, x_cr)

# Calculate equilibrium at each point
result = equilibrium(db, comps, phases,
                    {v.T: T_grid.flatten(),
                     v.P: 101325,
                     v.X('CR'): X_grid.flatten()},
                    pdens=500)

# Extract heat capacity
cp = result.CPM.values.reshape(T_grid.shape)

# Plot
plt.contourf(X_grid, T_grid, cp, levels=20, cmap='viridis')
plt.colorbar(label='Heat Capacity (J/mol-atom-K)')
plt.xlabel('X(CR)')
plt.ylabel('Temperature (K)')
plt.title('Cr-Fe Heat Capacity Map')
plt.savefig('crfe_cp_map.png')

4. Activity Calculation

from pycalphad import Database, equilibrium, variables as v

db = Database('feni.tdb')
comps = ['FE', 'NI', 'VA']
phases = ['FCC_A1']

# Calculate at specific conditions
result = equilibrium(db, comps, phases,
                    {v.T: 1200, v.P: 101325, v.X('NI'): 0.3})

# Extract activities (relative to pure element reference state)
activity_fe = result['ACR_FE'].values.squeeze()
activity_ni = result['ACR_NI'].values.squeeze()

print(f"Activity of Fe: {activity_fe:.4f}")
print(f"Activity of Ni: {activity_ni:.4f}")

# Chemical potentials
mu_fe = result['MU_FE'].values.squeeze()
mu_ni = result['MU_NI'].values.squeeze()

print(f"Chemical potential Fe: {mu_fe:.2f} J/mol")
print(f"Chemical potential Ni: {mu_ni:.2f} J/mol")

5. Ternary Phase Diagram

from pycalphad import Database, equilibrium, variables as v
from pycalphad.plot.eqplot import eqplot
import numpy as np
import matplotlib.pyplot as plt

db = Database('ternary.tdb')
comps = ['AL', 'CU', 'ZN', 'VA']
phases = ['LIQUID', 'FCC_A1', 'HCP_A3', 'THETA']

# Calculate equilibrium at constant T
result = equilibrium(db, comps, phases,
                    {v.T: 800, v.P: 101325,
                     v.X('CU'): (0, 1, 0.02),
                     v.X('ZN'): (0, 1, 0.02)},
                    pdens=2000)

# Plot using ternary coordinates
fig = plt.figure(figsize=(8, 8))
eqplot(result, x=v.X('CU'), y=v.X('ZN'))
plt.title('Al-Cu-Zn Ternary at 800K')
plt.savefig('alcuzn_ternary.png')

6. Driving Force for Precipitation

from pycalphad import Database, equilibrium, variables as v
import numpy as np

db = Database('alcu.tdb')
comps = ['AL', 'CU', 'VA']

# Supersaturated parent phase
parent_phases = ['FCC_A1']

# Calculate parent phase energy (metastable)
parent_eq = equilibrium(db, comps, parent_phases,
                       {v.T: 500, v.P: 101325, v.X('CU'): 0.02})
gm_parent = parent_eq.GM.values.squeeze()

# Equilibrium with precipitate phase allowed
all_phases = ['FCC_A1', 'THETA']
eq = equilibrium(db, comps, all_phases,
                {v.T: 500, v.P: 101325, v.X('CU'): 0.02})
gm_stable = eq.GM.values.squeeze()

# Driving force for precipitation
driving_force = gm_parent - gm_stable
print(f"Driving force: {driving_force:.2f} J/mol-atom")

7. T0 Temperature Calculation

from pycalphad import Database, equilibrium, variables as v
import numpy as np
from scipy.optimize import brentq

db = Database('steel.tdb')
comps = ['FE', 'C', 'VA']

def t0_condition(temp, composition):
    """Calculate GM difference between phases at equal composition"""
    # FCC energy
    fcc_eq = equilibrium(db, comps, ['FCC_A1'],
                        {v.T: temp, v.P: 101325, v.X('C'): composition})
    gm_fcc = fcc_eq.GM.values.squeeze()

    # BCC energy
    bcc_eq = equilibrium(db, comps, ['BCC_A2'],
                        {v.T: temp, v.P: 101325, v.X('C'): composition})
    gm_bcc = bcc_eq.GM.values.squeeze()

    return gm_fcc - gm_bcc

# Find T0 temperature where FCC and BCC have equal energy
composition = 0.01  # 1 at% C