Chapter 6: Practical Phase Diagram Calculation with pycalphad

Learning Objectives

By completing this final chapter, you will acquire the following practical skills:

• ✅ Understand the basic operations and main features of the pycalphad library
• ✅ Understand how to load and structure TDB files (thermodynamic databases)
• ✅ Calculate and visualize binary phase diagrams (Al-Cu, Fe-C, etc.)
• ✅ Calculate ternary isothermal sections (Al-Cu-Mg, etc.)
• ✅ Calculate equilibrium phase compositions and phase fractions
• ✅ Calculate driving forces and analyze phase transformations
• ✅ Understand application examples in real material systems and establish a research workflow
• ✅ Practice best practices for phase diagram calculation

1. Overview of pycalphad and Environment Setup

1.1 What Is pycalphad?

pycalphad is an open-source library for performing phase diagram calculations based on the CALPHAD method (CALculation of PHAse Diagrams) in Python. It reads CALPHAD databases (TDB files) and can compute phase equilibria, thermodynamic quantities, driving forces, and more for multicomponent systems.

1.2 Installation and Setup

# Run in a terminal/command prompt

# Basic installation (conda recommended)
conda install -c conda-forge pycalphad

# Or install with pip
pip install pycalphad

# Install visualization libraries as well
pip install matplotlib numpy scipy xarray

# Verify the installation
python -c "import pycalphad; print(pycalphad.__version__)"
# Example output: 0.10.3

# Jupyter Notebook is recommended for interactive work
pip install jupyter
jupyter notebook

1.3 How to Obtain TDB Files

pycalphad uses thermodynamic databases in TDB (ThermoCalc DataBase) format. TDB files contain the thermodynamic descriptions of phases (Gibbs energy functions).

2. Basic Operations of pycalphad

2.1 Loading a Database and Retrieving Basic Information

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

# Simplified demo Al-Cu TDB file (defined as a string)
# In actual research, load from an external file
demo_tdb = """
$ Al-Cu binary system (simplified demo)
ELEMENT AL FCC_A1 26.98 4577.3 28.3 !
ELEMENT CU FCC_A1 63.546 5004.1 33.15 !

FUNCTION GHSERAL 298.15 -7976.15+137.093038*T-24.3671976*T*LN(T)
    -.001884662*T**2-8.77664E-07*T**3+74092*T**(-1); 933.47 Y
    -11276.24+223.048446*T-38.5844296*T*LN(T)+.018531982*T**2
    -5.764227E-06*T**3+74092*T**(-1); 2900 N !

FUNCTION GHSERCU 298.15 -7770.458+130.485235*T-24.112392*T*LN(T)
    -.00265684*T**2+1.29223E-07*T**3+52478*T**(-1); 1358 Y
    -13542.026+183.803828*T-31.38*T*LN(T)+3.64167E+29*T**(-9); 3200 N !

PHASE FCC_A1 % 1 1 !
CONSTITUENT FCC_A1 : AL,CU : !
PARAMETER G(FCC_A1,AL;0) 298.15 +GHSERAL; 6000 N !
PARAMETER G(FCC_A1,CU;0) 298.15 +GHSERCU; 6000 N !
PARAMETER G(FCC_A1,AL,CU;0) 298.15 -53520+7.2*T; 6000 N !

PHASE LIQUID % 1 1 !
CONSTITUENT LIQUID : AL,CU : !
PARAMETER G(LIQUID,AL;0) 298.15 +GHSERAL+11005-11.841*T; 6000 N !
PARAMETER G(LIQUID,CU;0) 298.15 +GHSERCU+12964-9.511*T; 6000 N !
PARAMETER G(LIQUID,AL,CU;0) 298.15 -66200+40*T; 6000 N !
"""

# Load the TDB file
# For an actual file: db = Database('alcu.tdb')
from io import StringIO
db = Database(StringIO(demo_tdb))

print("=== Basic Database Information ===\n")

# Elements included
print(f"Elements: {sorted(db.elements)}")

# Defined phases
print(f"Phases: {sorted(db.phases.keys())}")

# Constituent information for each phase
print("\n=== Detailed Phase Information ===")
for phase_name in sorted(db.phases.keys()):
    phase = db.phases[phase_name]
    print(f"\n{phase_name} phase:")
    print(f"  Number of sublattices: {len(phase.constituents)}")
    print(f"  Constituents: {phase.constituents}")

# Basic element information
print("\n=== Element Information ===")
for element in sorted(db.elements):
    if element != 'VA':  # Exclude vacancies
        print(f"{element}:")
        # Enthalpy reference state of the element
        # Actual TDBs contain more detailed information

print("\n=== Checking pycalphad Variables ===")
print(f"Temperature variable: {v.T}")
print(f"Pressure variable: {v.P}")
print(f"Composition variable (example): {v.X('AL')}")

print("\nDatabase loaded successfully.")
print("Next, we proceed to phase diagram calculation.")

2.2 Setting Variables and Conditions

In pycalphad, calculation conditions are set using the variables module. The main variables are as follows:

3. Calculating Binary Phase Diagrams

3.1 Al-Cu Binary Phase Diagram

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

# Use demo_tdb defined in the previous code example
from io import StringIO
db = Database(StringIO(demo_tdb))

# Set the calculation conditions
components = ['AL', 'CU', 'VA']  # VA (vacancy) is required
phases = list(db.phases.keys())  # ['FCC_A1', 'LIQUID']

# Temperature range: 600 K to 1400 K (covering the vicinity of the Al melting point)
temperatures = np.linspace(600, 1400, 50)

# Composition range: Cu 0 to 100 at.% (Al balance computed automatically)
cu_compositions = np.linspace(0, 1, 50)

print("=== Calculation of the Al-Cu Binary Phase Diagram ===\n")
print(f"Temperature range: {temperatures[0]:.0f} - {temperatures[-1]:.0f} K")
print(f"Composition range: Cu {cu_compositions[0]:.1%} - {cu_compositions[-1]:.1%}")
print(f"Number of calculation points: {len(temperatures)} × {len(cu_compositions)} = {len(temperatures)*len(cu_compositions)}")

# Run the equilibrium calculation
try:
    eq_result = equilibrium(
        db,
        components,
        phases,
        {
            v.X('CU'): cu_compositions,
            v.T: temperatures,
            v.P: 101325,  # Standard pressure (Pa)
            v.N: 1.0      # Total moles (any value is fine)
        },
        output='GM'  # Output molar Gibbs energy
    )

    print("\nCalculation complete!")

    # Visualize the results
    fig, ax = plt.subplots(figsize=(10, 7))

    # Extract and draw phase boundaries
    # eq_result.Phase stores the stable-phase information at each point
    # Points where NP (number of stable phases) changes are phase boundaries

    # Simplified visualization (a more refined approach is used in practice)
    for phase_name in phases:
        # Color-code the region where each phase is stable
        phase_fraction = eq_result.NP.where(
            eq_result.Phase == phase_name
        ).values

        # Draw the phase diagram (colormap)
        im = ax.contourf(
            cu_compositions * 100,  # Convert to at.%
            temperatures,
            phase_fraction.T,
            levels=[0, 0.01, 1],
            colors=['white', 'lightblue'] if phase_name == 'LIQUID' else ['white', 'lightcoral'],
            alpha=0.5
        )

    ax.set_xlabel('Cu (at.%)', fontsize=12, fontweight='bold')
    ax.set_ylabel('Temperature (K)', fontsize=12, fontweight='bold')
    ax.set_title('Al-Cu Binary Phase Diagram (Simplified Demo)',
                 fontsize=14, fontweight='bold')
    ax.grid(True, alpha=0.3)

    # Add a legend
    from matplotlib.patches import Patch
    legend_elements = [
        Patch(facecolor='lightcoral', alpha=0.5, label='FCC_A1 (Solid)'),
        Patch(facecolor='lightblue', alpha=0.5, label='LIQUID')
    ]
    ax.legend(handles=legend_elements, loc='upper right', fontsize=10)

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

    print("\nPhase diagram saved: alcu_phase_diagram_demo.png")

except Exception as e:
    print(f"\nAn error occurred during the calculation: {e}")
    print("The simplified demo TDB may not allow a complete phase diagram calculation.")
    print("Use a high-accuracy TDB file for actual research.")

# Best practices for actual research
print("\n=== Recommendations for Actual Research ===")
print("1. Use a high-accuracy commercial TDB (e.g., TCAL) or pycalphad-data")
print("2. Increase the number of calculation points for smoother phase boundaries (100×100 or more)")
print("3. Use the dedicated binplot() function for phase diagram plotting")
print("4. Validate accuracy by comparison with experimental data")

3.2 Fe-C Binary Phase Diagram (Fundamentals of Steel Materials)

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

# Simplified TDB for the Fe-C system (in practice, use a precise database such as TCAL5)
# Here we show conceptual demo code

print("=== Example Calculation of the Fe-C Binary Phase Diagram ===\n")

# An appropriate TDB file is required for the actual calculation
# Example: db = Database('tcfe9.tdb')  # Thermo-Calc steel database

# Demo mode: explanation of the calculation procedure
print("[Calculation Procedure]")
print("1. Load the Fe-C TDB file")
print("   db = Database('tcfe9.tdb')")
print()
print("2. Set the calculation conditions")
print("   components = ['FE', 'C', 'VA']")
print("   phases = ['LIQUID', 'BCC_A2', 'FCC_A1', 'CEMENTITE']")
print("   # BCC_A2: ferrite (alpha iron)")
print("   # FCC_A1: austenite (gamma iron)")
print("   # CEMENTITE: cementite (Fe3C)")
print()
print("3. Set the temperature and composition ranges")
print("   temperatures = np.linspace(800, 1800, 100)  # K")
print("   carbon_content = np.linspace(0, 0.08, 100)  # C content 0-8 wt.%")
print()
print("4. Equilibrium calculation")
print("   eq = equilibrium(db, components, phases, {")
print("       v.X('C'): carbon_content,")
print("       v.T: temperatures,")
print("       v.P: 101325")
print("   })")
print()
print("5. Visualize the phase diagram")
print("   # The liquidus, solidus, eutectoid temperature (727°C), etc. are visualized")
print()

# Theoretical values of important phase transformation temperatures
print("[Important Phase Transformation Temperatures of the Fe-C System]")
phase_transformations = {
    "Eutectoid temperature": "727°C (1000 K)",
    "Eutectoid composition": "0.76 wt.% C",
    "Peritectic temperature": "1493°C (1766 K)",
    "Eutectic temperature": "1147°C (1420 K)",
    "α-γ transformation temperature (pure Fe)": "910°C (1183 K)"
}

for name, value in phase_transformations.items():
    print(f"  {name}: {value}")

# Visualize the features of the actual phase diagram (schematic)
fig, ax = plt.subplots(figsize=(11, 8))

# Schematic phase diagram drawing (conceptual diagram, not real data)
# The equilibrium() call above is required to display actual calculation results

# Temperature axis (vertical)
T_liquid = np.array([1768, 1493, 1147])  # K
T_alpha = np.array([1183, 1000, 1000])
T_gamma = np.array([1493, 1147, 1000])

# Composition axis (horizontal, wt.% C)
C_liquid = np.array([0, 0.17, 4.3])
C_alpha = np.array([0, 0.02, 0.02])
C_gamma = np.array([0.17, 2.14, 0.76])

ax.plot(C_liquid, T_liquid, 'b-', linewidth=2, label='Liquidus')
ax.plot(C_alpha, T_alpha, 'r-', linewidth=2, label='α (BCC) boundary')
ax.plot(C_gamma, T_gamma, 'g-', linewidth=2, label='γ (FCC) boundary')

# Mark the important points
ax.plot(0.76, 1000, 'ko', markersize=10, label='Eutectoid point')
ax.plot(0.17, 1766, 'ms', markersize=10, label='Peritectic point')
ax.plot(4.3, 1420, 'c^', markersize=10, label='Eutectic point')

# Labels for the phase regions
ax.text(0.1, 1400, 'LIQUID', fontsize=12, fontweight='bold', color='blue')
ax.text(0.4, 1300, 'L + γ', fontsize=11, color='green')
ax.text(0.01, 1050, 'α', fontsize=12, fontweight='bold', color='red')
ax.text(0.3, 1100, 'γ', fontsize=12, fontweight='bold', color='green')
ax.text(0.4, 950, 'α + Fe₃C', fontsize=11, color='purple')
ax.text(1.5, 950, 'γ + Fe₃C', fontsize=11, color='purple')

ax.set_xlabel('Carbon content (wt.%)', fontsize=13, fontweight='bold')
ax.set_ylabel('Temperature (K)', fontsize=13, fontweight='bold')
ax.set_title('Fe-C Binary Phase Diagram (Schematic)',
             fontsize=15, fontweight='bold')
ax.set_xlim(0, 2.5)
ax.set_ylim(900, 1900)
ax.grid(True, alpha=0.3)
ax.legend(loc='upper right', fontsize=10)

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

print("\nSchematic diagram saved: fec_phase_diagram_schematic.png")
print("\n[Note] The above is a schematic. Actual research requires calculation from a TDB file.")

4. Calculating Ternary Phase Diagrams

4.1 Isothermal Sections

For ternary systems, it is common to calculate an isothermal section at a fixed temperature. Phase regions are displayed on a triangular diagram (Gibbs triangle).

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

print("=== Example Calculation of an Al-Cu-Mg Ternary Isothermal Section ===\n")

# A ternary TDB file is required for the actual calculation
# Example: db = Database('alzn_mey.tdb')  # TDB for systems such as Al-Zn-Mg

# Explanation of the calculation procedure (demo mode)
print("[Calculation Procedure]")
print("1. Load the ternary TDB file")
print("   db = Database('al-cu-mg.tdb')")
print()
print("2. Set the temperature of the isothermal section (e.g., 600 K)")
print("   temperature = 600  # K (about 327°C)")
print()
print("3. Generate a composition grid on the triangular diagram")
print("   # Evenly distribute points on the Al-Cu-Mg triangle")
print("   # Constraint: X(Al) + X(Cu) + X(Mg) = 1")
print()
print("4. Run an equilibrium calculation at each composition point")
print("   eq = equilibrium(db, components, phases, {")
print("       v.X('CU'): cu_fractions,")
print("       v.X('MG'): mg_fractions,")
print("       v.T: temperature,")
print("       v.P: 101325")
print("   })")
print()
print("5. Visualize the phase regions on the triangular diagram")
print()

# Coordinate transformation function for the triangular diagram (for the Gibbs triangle)
def ternary_to_cartesian(a, b, c):
    """
    Convert ternary compositions (a, b, c) to Cartesian coordinates (x, y) on the triangular diagram

    Parameters:
    -----------
    a, b, c : float or array
        Compositions of the three components (a + b + c = 1)

    Returns:
    --------
    x, y : float or array
        Coordinates on the triangular diagram
    """
    x = 0.5 * (2 * b + c) / (a + b + c)
    y = (np.sqrt(3) / 2) * c / (a + b + c)
    return x, y

# Draw the frame of the triangular diagram
fig, ax = plt.subplots(figsize=(10, 9))

# Vertices of the triangle (Al, Cu, Mg)
vertices = np.array([
    [0, 0],           # Al (100%)
    [1, 0],           # Cu (100%)
    [0.5, np.sqrt(3)/2]  # Mg (100%)
])

triangle = plt.Polygon(vertices, fill=False, edgecolor='black', linewidth=2)
ax.add_patch(triangle)

# Vertex labels
ax.text(0, -0.05, 'Al', fontsize=14, fontweight='bold', ha='center', va='top')
ax.text(1, -0.05, 'Cu', fontsize=14, fontweight='bold', ha='center', va='top')
ax.text(0.5, np.sqrt(3)/2 + 0.05, 'Mg', fontsize=14, fontweight='bold',
        ha='center', va='bottom')

# Grid lines (10% intervals)
for i in range(1, 10):
    frac = i / 10

    # Lines parallel to the Al-Cu edge (constant Mg)
    x1, y1 = ternary_to_cartesian(1-frac, 0, frac)
    x2, y2 = ternary_to_cartesian(0, 1-frac, frac)
    ax.plot([x1, x2], [y1, y2], 'k-', linewidth=0.5, alpha=0.3)

    # Lines parallel to the Cu-Mg edge (constant Al)
    x1, y1 = ternary_to_cartesian(frac, 1-frac, 0)
    x2, y2 = ternary_to_cartesian(frac, 0, 1-frac)
    ax.plot([x1, x2], [y1, y2], 'k-', linewidth=0.5, alpha=0.3)

    # Lines parallel to the Mg-Al edge (constant Cu)
    x1, y1 = ternary_to_cartesian(1-frac, frac, 0)
    x2, y2 = ternary_to_cartesian(0, frac, 1-frac)
    ax.plot([x1, x2], [y1, y2], 'k-', linewidth=0.5, alpha=0.3)

# Show schematic phase regions (instead of actual calculation results)
# In practice, create a contour plot from the equilibrium() results

# FCC_A1 phase region (Al-rich solid solution)
al_rich_x = [0, 0.3, 0.15, 0]
al_rich_y = [0, 0, 0.13, 0]
ax.fill(al_rich_x, al_rich_y, color='lightcoral', alpha=0.5,
        label='FCC_A1 (Al-rich)')

# Other phase regions (conceptual)
ax.text(0.15, 0.05, 'α-Al', fontsize=11, fontweight='bold', ha='center')
ax.text(0.6, 0.3, 'Intermetallic\nphases', fontsize=10, ha='center', style='italic')

ax.set_xlim(-0.1, 1.1)
ax.set_ylim(-0.1, 1.0)
ax.set_aspect('equal')
ax.axis('off')
ax.set_title('Al-Cu-Mg Ternary Isothermal Section at 600 K (Schematic)',
             fontsize=14, fontweight='bold', pad=20)

# Legend
ax.legend(loc='upper right', fontsize=10)

plt.tight_layout()
plt.savefig('al-cu-mg_ternary_schematic.png', dpi=150, bbox_inches='tight')
plt.show()

print("\nTriangular diagram saved: al-cu-mg_ternary_schematic.png")
print("\n[Note] The above is a schematic. Actual research requires calculation from a TDB file.")

print("\n=== Practical Points for Ternary Calculations ===")
print("1. Computation time: ternary systems cost more than binary systems (minutes to hours)")
print("2. Grid resolution: about 50×50 is standard, 100×100 for high accuracy")
print("3. Phase selection: choose the phases to calculate appropriately (all phases makes it heavy)")
print("4. Visualization: triangular plots, contour maps, phase fraction maps, etc.")

4.2 Vertical Sections

Vertical sections, in which a specific composition ratio of a ternary system is fixed and the temperature is varied, are also important. For example, in the Al-Cu-Mg system, one can fix the Mg concentration and investigate the relationship between the Al-Cu ratio and temperature.

5. Calculating Equilibrium Phase Compositions and Phase Fractions

5.1 Equilibrium Calculation under Specified Conditions

import pycalphad as pyc
from pycalphad import Database, equilibrium, variables as v
import numpy as np
import pandas as pd

print("=== Calculating Equilibrium Phase Compositions and Phase Fractions ===\n")

# Use demo_tdb (Al-Cu system) from the previous code example
from io import StringIO
db = Database(StringIO(demo_tdb))

components = ['AL', 'CU', 'VA']
phases = list(db.phases.keys())

# Calculation conditions: a specific temperature and composition
temperature = 900  # K (about 627°C)
cu_content = 0.3   # Cu 30 at.%

print(f"Calculation conditions:")
print(f"  Temperature: {temperature} K ({temperature - 273.15:.1f}°C)")
print(f"  Composition: Al-{cu_content*100:.0f}at.%Cu")
print(f"  Pressure: 101325 Pa (1 atm)")
print()

try:
    # Equilibrium calculation
    eq = equilibrium(
        db, components, phases,
        {
            v.X('CU'): cu_content,
            v.T: temperature,
            v.P: 101325,
            v.N: 1.0
        },
        output=['GM', 'NP', 'X']  # Molar Gibbs energy, phase fraction, composition
    )

    print("Equilibrium calculation complete!\n")
    print("=== Calculation Results ===\n")

    # Extract the stable phases
    stable_phases = []
    for phase_name in phases:
        # Phases with a phase fraction (NP) greater than 0 are stable
        phase_fraction = float(eq.NP.sel(Phase=phase_name).values)
        if phase_fraction > 1e-6:  # Exclude tiny numerical errors
            stable_phases.append({
                'Phase': phase_name,
                'Fraction': phase_fraction
            })

    if stable_phases:
        print("Stable phases:")
        for phase_data in stable_phases:
            print(f"  {phase_data['Phase']}: {phase_data['Fraction']:.4f} ({phase_data['Fraction']*100:.2f}%)")
    else:
        print("No stable phases found (possible calculation error)")

    print()

    # Composition of each phase
    print("Composition of each phase:")
    for phase_data in stable_phases:
        phase_name = phase_data['Phase']
        print(f"\n  {phase_name} phase:")

        # Compositions of Al and Cu
        try:
            al_content = float(eq.X.sel(Phase=phase_name, component='AL').values)
            cu_content_phase = float(eq.X.sel(Phase=phase_name, component='CU').values)

            print(f"    Al: {al_content:.4f} ({al_content*100:.2f} at.%)")
            print(f"    Cu: {cu_content_phase:.4f} ({cu_content_phase*100:.2f} at.%)")
        except:
            print("    Failed to retrieve composition data")

    # Gibbs energy
    print(f"\nMolar Gibbs energy of the whole system:")
    gm_total = float(eq.GM.values)
    print(f"  {gm_total:.2f} J/mol")

    print("\n=== Application to Practical Materials Design ===")
    print("From this information, the following can be predicted:")
    print("1. Which phases precipitate (judged from phase fractions)")
    print("2. Composition of each phase (degree of segregation)")
    print("3. Microstructure after heat treatment (phase stability)")
    print("4. Effect on mechanical properties (types and fractions of phases)")

except Exception as e:
    print(f"Calculation error: {e}")
    print("The demo TDB may not allow detailed calculations.")

# Track changes in phase fractions with temperature
print("\n\n=== Change in Phase Fractions with Temperature ===\n")

temperatures_range = np.linspace(700, 1100, 20)
results = []

for T in temperatures_range:
    try:
        eq_temp = equilibrium(
            db, components, phases,
            {v.X('CU'): 0.3, v.T: T, v.P: 101325, v.N: 1.0},
            output='NP'
        )

        fcc_fraction = float(eq_temp.NP.sel(Phase='FCC_A1').values)
        liquid_fraction = float(eq_temp.NP.sel(Phase='LIQUID').values)

        results.append({
            'T (K)': T,
            'T (°C)': T - 273.15,
            'FCC_A1': fcc_fraction,
            'LIQUID': liquid_fraction
        })
    except:
        pass

if results:
    df = pd.DataFrame(results)
    print(df.to_string(index=False))

    # Visualization
    import matplotlib.pyplot as plt

    fig, ax = plt.subplots(figsize=(10, 6))
    ax.plot(df['T (K)'], df['FCC_A1'], 'ro-', linewidth=2,
            markersize=6, label='FCC_A1 (Solid)')
    ax.plot(df['T (K)'], df['LIQUID'], 'bo-', linewidth=2,
            markersize=6, label='LIQUID')

    ax.set_xlabel('Temperature (K)', fontsize=12, fontweight='bold')
    ax.set_ylabel('Phase Fraction', fontsize=12, fontweight='bold')
    ax.set_title('Phase Fraction vs Temperature (Al-30at.%Cu)',
                 fontsize=14, fontweight='bold')
    ax.legend(fontsize=11)
    ax.grid(True, alpha=0.3)
    ax.set_ylim(-0.05, 1.05)

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

    print("\nPhase fraction plot saved: phase_fraction_vs_temperature.png")

5.2 Phase Fraction Calculation with the Lever Rule

In a two-phase coexistence region, the amount of each phase can be determined using the lever rule. pycalphad calculates this automatically, but it is important to understand the principle.

6. Calculating the Driving Force

6.1 What Is a Driving Force?

The driving force is a thermodynamic measure of how readily a transformation from one phase to another proceeds. It is defined as a difference in chemical potentials.

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

print("=== Calculating the Driving Force for a Phase Transformation ===\n")

# Use demo_tdb
from io import StringIO
db = Database(StringIO(demo_tdb))

components = ['AL', 'CU', 'VA']
phases = list(db.phases.keys())

# Calculation conditions
cu_composition = 0.2  # Cu 20 at.%
temperatures = np.linspace(700, 1200, 30)

print(f"Composition: Al-{cu_composition*100:.0f}at.%Cu")
print(f"Temperature range: {temperatures[0]:.0f} - {temperatures[-1]:.0f} K")
print()

# Calculate the Gibbs energy of each phase at each temperature
fcc_energies = []
liquid_energies = []
driving_forces = []

for T in temperatures:
    try:
        # Gibbs energy of the FCC phase
        eq_fcc = equilibrium(
            db, components, ['FCC_A1'],  # FCC phase only
            {v.X('CU'): cu_composition, v.T: T, v.P: 101325},
            output='GM'
        )
        G_fcc = float(eq_fcc.GM.values)
        fcc_energies.append(G_fcc)

        # Gibbs energy of the LIQUID phase
        eq_liquid = equilibrium(
            db, components, ['LIQUID'],  # LIQUID phase only
            {v.X('CU'): cu_composition, v.T: T, v.P: 101325},
            output='GM'
        )
        G_liquid = float(eq_liquid.GM.values)
        liquid_energies.append(G_liquid)

        # Driving force: LIQUID - FCC (positive means the liquid phase is stable)
        driving_force = G_liquid - G_fcc
        driving_forces.append(driving_force)

    except Exception as e:
        fcc_energies.append(np.nan)
        liquid_energies.append(np.nan)
        driving_forces.append(np.nan)

print("Calculation complete!\n")

# Visualize the results
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 10))

# Subplot 1: Gibbs energy
ax1.plot(temperatures, fcc_energies, 'r-', linewidth=2.5, label='FCC_A1 (Solid)')
ax1.plot(temperatures, liquid_energies, 'b-', linewidth=2.5, label='LIQUID')
ax1.axhline(y=0, color='gray', linestyle='--', linewidth=1, alpha=0.5)
ax1.set_xlabel('Temperature (K)', fontsize=12, fontweight='bold')
ax1.set_ylabel('Gibbs Energy (J/mol)', fontsize=12, fontweight='bold')
ax1.set_title(f'Gibbs Energy vs Temperature (Al-{cu_composition*100:.0f}at.%Cu)',
              fontsize=13, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# Subplot 2: Driving force
ax2.plot(temperatures, driving_forces, 'g-', linewidth=2.5, label='ΔG (LIQUID - FCC)')
ax2.axhline(y=0, color='red', linestyle='--', linewidth=2, alpha=0.7, label='Equilibrium')
ax2.fill_between(temperatures, 0, driving_forces,
                 where=np.array(driving_forces) < 0, alpha=0.3, color='blue',
                 label='FCC stable region')
ax2.fill_between(temperatures, 0, driving_forces,
                 where=np.array(driving_forces) > 0, alpha=0.3, color='red',
                 label='LIQUID stable region')
ax2.set_xlabel('Temperature (K)', fontsize=12, fontweight='bold')
ax2.set_ylabel('Driving Force (J/mol)', fontsize=12, fontweight='bold')
ax2.set_title(f'Driving Force for Melting (Al-{cu_composition*100:.0f}at.%Cu)',
              fontsize=13, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

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

print("Driving force analysis plot saved: driving_force_analysis.png\n")

# Estimate the melting point (the temperature at which the driving force becomes zero)
try:
    # Find the zero-crossing point
    for i in range(len(driving_forces) - 1):
        if not np.isnan(driving_forces[i]) and not np.isnan(driving_forces[i+1]):
            if driving_forces[i] * driving_forces[i+1] < 0:  # Sign change
                T_melting = (temperatures[i] + temperatures[i+1]) / 2
                print(f"Estimated melting point: {T_melting:.1f} K ({T_melting - 273.15:.1f}°C)\n")
                break
except:
    pass

print("=== Applications of Driving Force Calculations ===")
print("1. Assessing the thermodynamic feasibility of phase transformations")
print("2. Quantifying the degree of undercooling/overheating")
print("3. Estimating nucleation and growth rates (combined with kinetics)")
print("4. Optimizing heat treatment conditions")
print("5. Evaluating the stability of metastable phases")

6.2 Relationship between Driving Force and Transformation Rate

The driving force represents the "thermodynamic impetus" for a phase transformation, but the actual transformation rate also depends on kinetics (diffusion, interface migration, etc.). An analysis combining both is indispensable for designing real materials processes.

7. Practical Workflows and Application Examples

7.1 Complete Research Workflow

import pycalphad as pyc
from pycalphad import Database, equilibrium, variables as v
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from datetime import datetime

print("=" * 60)
print("pycalphad Integrated Practical Workflow Script")
print("=" * 60)
print()

# ============================================
# Step 1: Project setup
# ============================================
print("[Step 1] Project setup\n")

project_name = "Al-Cu Alloy Design"
target_composition = "Al-20at.%Cu"
target_temperature_range = (600, 1200)  # K
objective = "Determination of optimal heat treatment conditions"

print(f"Project name: {project_name}")
print(f"Target alloy: {target_composition}")
print(f"Temperature range: {target_temperature_range[0]} - {target_temperature_range[1]} K")
print(f"Objective: {objective}")
print()

# ============================================
# Step 2: Database preparation
# ============================================
print("[Step 2] Loading the TDB file\n")

# In actual research, use an external TDB file
# db = Database('path/to/your/database.tdb')

# For the demo
from io import StringIO
db = Database(StringIO(demo_tdb))

components = ['AL', 'CU', 'VA']
phases = list(db.phases.keys())

print(f"TDB in use: demo_al_cu.tdb")
print(f"Elements: {[c for c in components if c != 'VA']}")
print(f"Phases: {phases}")
print()

# ============================================
# Step 3: Phase diagram calculation
# ============================================
print("[Step 3] Phase diagram calculation\n")

T_range = np.linspace(target_temperature_range[0],
                      target_temperature_range[1], 40)
X_Cu_range = np.linspace(0, 0.5, 40)  # 0-50 at.% Cu

print(f"Calculating... ({len(T_range)} temperature points × {len(X_Cu_range)} composition points)")

try:
    eq_phase_diagram = equilibrium(
        db, components, phases,
        {v.X('CU'): X_Cu_range, v.T: T_range, v.P: 101325},
        output='GM'
    )
    print("✓ Phase diagram calculation complete")
except Exception as e:
    print(f"✗ Phase diagram calculation error: {e}")
    eq_phase_diagram = None

print()

# ============================================
# Step 4: Equilibrium calculation under specific conditions
# ============================================
print("[Step 4] Equilibrium calculation under specific conditions\n")

specific_T = 900  # K
specific_X_Cu = 0.2  # 20 at.%

print(f"Conditions: {specific_T} K, Cu {specific_X_Cu*100:.0f} at.%")

try:
    eq_specific = equilibrium(
        db, components, phases,
        {v.X('CU'): specific_X_Cu, v.T: specific_T, v.P: 101325},
        output=['GM', 'NP', 'X']
    )

    print("\nStable phases:")
    for phase_name in phases:
        fraction = float(eq_specific.NP.sel(Phase=phase_name).values)
        if fraction > 1e-6:
            print(f"  {phase_name}: {fraction:.4f} ({fraction*100:.2f}%)")

    print("✓ Equilibrium calculation complete")
except Exception as e:
    print(f"✗ Equilibrium calculation error: {e}")
    eq_specific = None

print()

# ============================================
# Step 5: Driving force analysis
# ============================================
print("[Step 5] Driving force analysis\n")

T_df_range = np.linspace(700, 1100, 25)
driving_forces_analysis = []

for T in T_df_range:
    try:
        eq_fcc = equilibrium(db, components, ['FCC_A1'],
                            {v.X('CU'): 0.2, v.T: T, v.P: 101325}, output='GM')
        eq_liq = equilibrium(db, components, ['LIQUID'],
                            {v.X('CU'): 0.2, v.T: T, v.P: 101325}, output='GM')

        G_fcc = float(eq_fcc.GM.values)
        G_liq = float(eq_liq.GM.values)
        df_value = G_liq - G_fcc

        driving_forces_analysis.append({
            'T (K)': T,
            'ΔG (J/mol)': df_value
        })
    except:
        pass

if driving_forces_analysis:
    df_analysis = pd.DataFrame(driving_forces_analysis)
    print(df_analysis.head(10).to_string(index=False))
    print("...")
    print("✓ Driving force analysis complete")
else:
    print("✗ Driving force analysis error")

print()

# ============================================
# Step 6: Visualization of results
# ============================================
print("[Step 6] Visualization of results\n")

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

# (1) Phase diagram (simplified version)
ax1 = axes[0, 0]
if eq_phase_diagram is not None:
    # Draw the phase diagram (more detailed processing is required in practice)
    ax1.text(0.5, 0.5, 'Phase Diagram\n(Requires detailed plotting)',
             ha='center', va='center', fontsize=12, transform=ax1.transAxes)
ax1.set_title('Phase Diagram', fontweight='bold')
ax1.set_xlabel('Cu (at.%)')
ax1.set_ylabel('Temperature (K)')
ax1.grid(True, alpha=0.3)

# (2) Phase fraction vs temperature
ax2 = axes[0, 1]
if not driving_forces_analysis:
    ax2.text(0.5, 0.5, 'No data', ha='center', va='center',
             transform=ax2.transAxes)
else:
    # Dummy data
    ax2.plot([700, 900, 1100], [1, 0.7, 0], 'ro-', label='FCC', linewidth=2)
    ax2.plot([700, 900, 1100], [0, 0.3, 1], 'bo-', label='LIQUID', linewidth=2)
ax2.set_title('Phase Fraction vs Temperature', fontweight='bold')
ax2.set_xlabel('Temperature (K)')
ax2.set_ylabel('Phase Fraction')
ax2.legend()
ax2.grid(True, alpha=0.3)

# (3) Driving force
ax3 = axes[1, 0]
if driving_forces_analysis:
    ax3.plot(df_analysis['T (K)'], df_analysis['ΔG (J/mol)'],
             'g-', linewidth=2.5)
    ax3.axhline(y=0, color='red', linestyle='--', linewidth=2)
    ax3.fill_between(df_analysis['T (K)'], 0, df_analysis['ΔG (J/mol)'],
                     alpha=0.3)
ax3.set_title('Driving Force (LIQUID - FCC)', fontweight='bold')
ax3.set_xlabel('Temperature (K)')
ax3.set_ylabel('ΔG (J/mol)')
ax3.grid(True, alpha=0.3)

# (4) Material design guidelines
ax4 = axes[1, 1]
ax4.axis('off')
design_guidelines = f"""
Material Design Guidelines

[Optimal Heat Treatment Conditions]
• Solution treatment: 950-1000 K
• Aging treatment: 450-500 K
• Cooling rate: rapid quenching recommended

[Expected Precipitate Phases]
• FCC_A1 (α-Al matrix)
• θ phase (Al₂Cu, precipitation strengthening)

[Predicted Mechanical Properties]
• Strengthening: via θ phase precipitation
• Ductility: depends on matrix properties

[Next Steps]
1. Experimental validation
2. Microstructure observation (SEM/TEM)
3. Mechanical testing (tensile, hardness)
"""
ax4.text(0.05, 0.95, design_guidelines,
         fontsize=10, verticalalignment='top',
         fontfamily='monospace',
         bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8))
ax4.set_title('Material Design Guidelines', fontweight='bold', loc='left')

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

print("✓ Comprehensive analysis figure saved: workflow_comprehensive_analysis.png")
print()

# ============================================
# Step 7: Report generation
# ============================================
print("[Step 7] Report generation\n")

report = f"""
{'='*60}
pycalphad Calculation Report
{'='*60}

Project name: {project_name}
Executed at: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}

[Calculation Conditions]
- Target system: {target_composition}
- Temperature range: {target_temperature_range[0]} - {target_temperature_range[1]} K
- Pressure: 101325 Pa (1 atm)

[Database Used]
- TDB: demo_al_cu.tdb
- Phases: {', '.join(phases)}

[Summary of Calculation Results]
- Phase diagram calculation: {'Success' if eq_phase_diagram is not None else 'Failure'}
- Equilibrium calculation: {'Success' if eq_specific is not None else 'Failure'}
- Driving force analysis: {'Success' if driving_forces_analysis else 'Failure'}

[Recommendations for Materials Design]
Toward {objective}, the following conditions are recommended:
1. Solution treatment: complete dissolution at 950-1000 K
2. Rapid quenching: retain the metastable state near room temperature
3. Aging treatment: form precipitation-strengthening phases at 450-500 K

[References]
[1] Pycalphad Documentation: https://pycalphad.org/
[2] CALPHAD methodology review papers
[3] Al-Cu binary phase diagram experimental data

{'='*60}
"""

print(report)

# Save the report to a file
with open('pycalphad_calculation_report.txt', 'w', encoding='utf-8') as f:
    f.write(report)

print("✓ Report saved: pycalphad_calculation_report.txt")
print()

print("=" * 60)
print("Workflow complete!")
print("=" * 60)

7.2 Best Practices

8. Exercises

Summary and Next Steps

What to Learn Next

For Continued Learning

Final Message

Materials thermodynamics is one of the most fundamental disciplines in materials science. It is a powerful tool for answering essential questions such as "Why is this phase stable?" and "Under what conditions do which phases coexist?"

In this series, you have learned systematically, from fundamental theory to practical calculations with pycalphad. In this final chapter in particular, you experienced a realistic research workflow using TDB files.

The knowledge and skills acquired here form the foundation for all kinds of materials research, including alloy design, process optimization, new materials exploration, and Materials Informatics.

From now on, actively apply the techniques you learned in this series to the material systems and research topics that interest you. pycalphad is a powerful tool with the potential to greatly accelerate your research.

Continuous learning and practice are the path to true mastery.
We hope this series marks a new start for your materials science research.

🎓 Introduction to Materials Thermodynamics Series - The End 🎓