第6章: pycalphadによる相図計算実践

学習目標

この最終章を完了することで、以下の実践的スキルを習得できます：

• ✅ pycalphadライブラリの基本操作と主要機能を理解する
• ✅ TDBファイル（熱力学データベース）の読み込みと構造を理解する
• ✅ 二元系相図（Al-Cu、Fe-C等）を計算し可視化できる
• ✅ 三元系等温断面図（Al-Cu-Mg等）を計算できる
• ✅ 平衡相組成と相分率を計算できる
• ✅ 駆動力（driving force）を計算し相変態を解析できる
• ✅ 実材料系での応用例を理解し、研究ワークフローを確立できる
• ✅ 相図計算のベストプラクティスを実践できる

1. pycalphadの概要と環境構築

1.1 pycalphadとは

pycalphadは、CALPHAD法（CALculation of PHAse Diagrams）に基づいた相図計算を Pythonで実行するためのオープンソースライブラリです。CALPHADデータベース（TDBファイル）を読み込み、 多成分系の相平衡、熱力学量、駆動力などを計算できます。

1.2 インストールとセットアップ

# ターミナル/コマンドプロンプトで実行

# 基本的なインストール（condaを推奨）
conda install -c conda-forge pycalphad

# または pip でのインストール
pip install pycalphad

# 可視化ライブラリも併せてインストール
pip install matplotlib numpy scipy xarray

# インストール確認
python -c "import pycalphad; print(pycalphad.__version__)"
# 出力例: 0.10.3

# インタラクティブな作業にはJupyter Notebookがおすすめ
pip install jupyter
jupyter notebook

1.3 TDBファイルの取得方法

pycalphadはTDB（ThermoCalc DataBase）形式の熱力学データベースを使用します。 TDBファイルには相の熱力学的記述（ギブスエネルギー関数）が含まれています。

2. pycalphadの基本操作

2.1 データベースの読み込みと基本情報の取得

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

# デモ用の簡易Al-Cu TDBファイル（文字列として定義）
# 実際の研究では外部ファイルから読み込む
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 !
"""

# TDBファイルの読み込み
# 実際のファイルの場合: db = Database('alcu.tdb')
from io import StringIO
db = Database(StringIO(demo_tdb))

print("=== データベース基本情報 ===\n")

# 含まれる元素
print(f"元素: {sorted(db.elements)}")

# 定義されている相
print(f"相: {sorted(db.phases.keys())}")

# 各相の構成情報
print("\n=== 相の詳細情報 ===")
for phase_name in sorted(db.phases.keys()):
    phase = db.phases[phase_name]
    print(f"\n{phase_name}相:")
    print(f"  サブラティス数: {len(phase.constituents)}")
    print(f"  構成元素: {phase.constituents}")

# 元素の基本情報
print("\n=== 元素情報 ===")
for element in sorted(db.elements):
    if element != 'VA':  # 空孔を除く
        print(f"{element}:")
        # 元素のエンタルピー参照状態
        # 実際のTDBにはより詳細な情報が含まれる

print("\n=== pycalphad 変数の確認 ===")
print(f"温度変数: {v.T}")
print(f"圧力変数: {v.P}")
print(f"組成変数（例）: {v.X('AL')}")

print("\nデータベースの読み込みに成功しました。")
print("次は相図計算に進みます。")

2.2 変数と条件の設定

pycalphadでは、variablesモジュールを使って計算条件を設定します。 主要な変数は以下の通りです：

3. 二元系相図の計算

3.1 Al-Cu二元系相図

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

# 前のコード例で定義したdemo_tdbを使用
from io import StringIO
db = Database(StringIO(demo_tdb))

# 計算条件の設定
components = ['AL', 'CU', 'VA']  # VA（空孔）は必須
phases = list(db.phases.keys())  # ['FCC_A1', 'LIQUID']

# 温度範囲：600 K 〜 1400 K（Al融点付近をカバー）
temperatures = np.linspace(600, 1400, 50)

# 組成範囲：Cu 0 〜 100 at.%（Alの残りを自動計算）
cu_compositions = np.linspace(0, 1, 50)

print("=== Al-Cu二元系相図の計算 ===\n")
print(f"温度範囲: {temperatures[0]:.0f} - {temperatures[-1]:.0f} K")
print(f"組成範囲: Cu {cu_compositions[0]:.1%} - {cu_compositions[-1]:.1%}")
print(f"計算点数: {len(temperatures)} × {len(cu_compositions)} = {len(temperatures)*len(cu_compositions)}")

# 平衡計算の実行
try:
    eq_result = equilibrium(
        db,
        components,
        phases,
        {
            v.X('CU'): cu_compositions,
            v.T: temperatures,
            v.P: 101325,  # 標準圧力（Pa）
            v.N: 1.0      # 全モル数（任意の値でOK）
        },
        output='GM'  # モルギブスエネルギーを出力
    )

    print("\n計算完了！")

    # 結果の可視化
    fig, ax = plt.subplots(figsize=(10, 7))

    # 相境界の抽出と描画
    # eq_result.Phase には各点での安定相情報が格納されている
    # NP（安定相の数）が変化する点が相境界

    # 簡易的な可視化（実際にはより洗練された方法を使用）
    for phase_name in phases:
        # 各相が安定な領域を色分け表示
        phase_fraction = eq_result.NP.where(
            eq_result.Phase == phase_name
        ).values

        # 相図の描画（カラーマップ）
        im = ax.contourf(
            cu_compositions * 100,  # 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)

    # 凡例を追加
    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("\n相図を保存しました: alcu_phase_diagram_demo.png")

except Exception as e:
    print(f"\n計算中にエラーが発生しました: {e}")
    print("デモ用簡易TDBでは完全な相図計算ができない場合があります。")
    print("実際の研究では高精度なTDBファイルを使用してください。")

# 実際の研究でのベストプラクティス
print("\n=== 実際の研究での推奨事項 ===")
print("1. 高精度な商用TDB（TCAL等）またはpycalphad-dataを使用")
print("2. 計算点数を増やして相境界を滑らかに（100×100以上）")
print("3. binplot()関数を使った専用の相図描画")
print("4. 実験データと比較して精度を検証")

3.2 Fe-C二元系相図（鉄鋼材料の基礎）

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

# Fe-C系の簡易TDB（実際にはTCAL5等の精密データベースを使用）
# ここではデモ用の概念的なコードを示します

print("=== Fe-C二元系相図の計算例 ===\n")

# 実際の計算には適切なTDBファイルが必要
# 例: db = Database('tcfe9.tdb')  # Thermo-Calc鉄鋼データベース

# デモモード：計算手順の説明
print("【計算手順】")
print("1. Fe-C系TDBファイルの読み込み")
print("   db = Database('tcfe9.tdb')")
print()
print("2. 計算条件の設定")
print("   components = ['FE', 'C', 'VA']")
print("   phases = ['LIQUID', 'BCC_A2', 'FCC_A1', 'CEMENTITE']")
print("   # BCC_A2: フェライト（α鉄）")
print("   # FCC_A1: オーステナイト（γ鉄）")
print("   # CEMENTITE: セメンタイト（Fe3C）")
print()
print("3. 温度・組成範囲の設定")
print("   temperatures = np.linspace(800, 1800, 100)  # K")
print("   carbon_content = np.linspace(0, 0.08, 100)  # C濃度 0-8 wt.%")
print()
print("4. 平衡計算")
print("   eq = equilibrium(db, components, phases, {")
print("       v.X('C'): carbon_content,")
print("       v.T: temperatures,")
print("       v.P: 101325")
print("   })")
print()
print("5. 相図の可視化")
print("   # 液相線、固相線、共析温度（727°C）などが可視化される")
print()

# 重要な相変態温度の理論値
print("【Fe-C系の重要な相変態温度】")
phase_transformations = {
    "共析温度（Eutectoid）": "727°C (1000 K)",
    "共析組成": "0.76 wt.% C",
    "包析温度（Peritectic）": "1493°C (1766 K)",
    "共晶温度（Eutectic）": "1147°C (1420 K)",
    "α-γ変態温度（純Fe）": "910°C (1183 K)"
}

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

# 実際の相図の特徴を可視化（模式図）
fig, ax = plt.subplots(figsize=(11, 8))

# 模式的な相図の描画（実データではなく概念図）
# 実際の計算結果を表示するには上記のequilibrium()が必要

# 温度軸（縦軸）
T_liquid = np.array([1768, 1493, 1147])  # K
T_alpha = np.array([1183, 1000, 1000])
T_gamma = np.array([1493, 1147, 1000])

# 組成軸（横軸、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')

# 重要な点をマーク
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')

# 相領域のラベル
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("\n模式図を保存しました: fec_phase_diagram_schematic.png")
print("\n【注意】上記は模式図です。実際の研究ではTDBファイルからの計算が必要です。")

4. 三元系相図の計算

4.1 等温断面図（Isothermal Section）

三元系では、温度を固定した等温断面図を計算することが一般的です。 三角図（Gibbs三角形）上に相領域を表示します。

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("=== Al-Cu-Mg三元系等温断面図の計算例 ===\n")

# 実際の計算には三元系TDBファイルが必要
# 例: db = Database('alzn_mey.tdb')  # Al-Zn-Mg系などのTDB

# 計算手順の説明（デモモード）
print("【計算手順】")
print("1. 三元系TDBファイルの読み込み")
print("   db = Database('al-cu-mg.tdb')")
print()
print("2. 等温断面の温度を設定（例：600 K）")
print("   temperature = 600  # K (約327°C)")
print()
print("3. 三角図上の組成グリッドを生成")
print("   # Al-Cu-Mg三角図上の点を均等に配置")
print("   # X(Al) + X(Cu) + X(Mg) = 1 の制約")
print()
print("4. 各組成点で平衡計算")
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. 三角図上に相領域を可視化")
print()

# 三角図の座標変換関数（Gibbs三角形用）
def ternary_to_cartesian(a, b, c):
    """
    三成分（a, b, c）の組成を三角図のデカルト座標（x, y）に変換

    Parameters:
    -----------
    a, b, c : float or array
        3成分の組成（a + b + c = 1）

    Returns:
    --------
    x, y : float or array
        三角図上の座標
    """
    x = 0.5 * (2 * b + c) / (a + b + c)
    y = (np.sqrt(3) / 2) * c / (a + b + c)
    return x, y

# 三角図の枠を描画
fig, ax = plt.subplots(figsize=(10, 9))

# 三角形の頂点（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)

# 頂点のラベル
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')

# グリッド線（10%刻み）
for i in range(1, 10):
    frac = i / 10

    # Al-Cu軸に平行な線（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)

    # Cu-Mg軸に平行な線（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)

    # Mg-Al軸に平行な線（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)

# 模式的な相領域を表示（実際の計算結果の代わり）
# 実際にはequilibrium()の結果からcontour plotを作成

# FCC_A1相領域（Al-richな固溶体）
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)')

# その他の相領域（概念的）
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)

# 凡例
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("\n三角図を保存しました: al-cu-mg_ternary_schematic.png")
print("\n【注意】上記は模式図です。実際の研究ではTDBファイルからの計算が必要です。")

print("\n=== 三元系計算の実践的ポイント ===")
print("1. 計算時間：三元系は二元系より計算コストが高い（数分〜数時間）")
print("2. グリッド解像度：50×50程度が標準、高精度なら100×100")
print("3. 相の選択：計算対象の相を適切に選ぶ（全相だと重くなる）")
print("4. 可視化：三角図プロット、等高線図、相分率マップ等")

4.2 垂直断面図（Vertical Section）

三元系の特定の組成比を固定し、温度を変化させた垂直断面図も重要です。 例えば、Al-Cu-Mg系でMg濃度を固定し、Al-Cu比と温度の関係を調べることができます。

5. 平衡相組成と相分率の計算

5.1 指定条件での平衡計算

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

print("=== 平衡相組成と相分率の計算 ===\n")

# 前のコード例のdemo_tdb（Al-Cu系）を使用
from io import StringIO
db = Database(StringIO(demo_tdb))

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

# 計算条件：特定の温度と組成
temperature = 900  # K （約627°C）
cu_content = 0.3   # Cu 30 at.%

print(f"計算条件:")
print(f"  温度: {temperature} K ({temperature - 273.15:.1f}°C)")
print(f"  組成: Al-{cu_content*100:.0f}at.%Cu")
print(f"  圧力: 101325 Pa (1 atm)")
print()

try:
    # 平衡計算
    eq = equilibrium(
        db, components, phases,
        {
            v.X('CU'): cu_content,
            v.T: temperature,
            v.P: 101325,
            v.N: 1.0
        },
        output=['GM', 'NP', 'X']  # モルギブスエネルギー、相分率、組成
    )

    print("平衡計算完了！\n")
    print("=== 計算結果 ===\n")

    # 安定相の抽出
    stable_phases = []
    for phase_name in phases:
        # 相分率（NP）が0より大きい相が安定相
        phase_fraction = float(eq.NP.sel(Phase=phase_name).values)
        if phase_fraction > 1e-6:  # 微小な数値誤差を除外
            stable_phases.append({
                'Phase': phase_name,
                'Fraction': phase_fraction
            })

    if stable_phases:
        print("安定相:")
        for phase_data in stable_phases:
            print(f"  {phase_data['Phase']}: {phase_data['Fraction']:.4f} ({phase_data['Fraction']*100:.2f}%)")
    else:
        print("安定相が見つかりませんでした（計算エラーの可能性）")

    print()

    # 各相の組成
    print("各相の組成:")
    for phase_data in stable_phases:
        phase_name = phase_data['Phase']
        print(f"\n  {phase_name}相:")

        # Al, 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("    組成データの取得に失敗")

    # ギブスエネルギー
    print(f"\nシステム全体のモルギブスエネルギー:")
    gm_total = float(eq.GM.values)
    print(f"  {gm_total:.2f} J/mol")

    print("\n=== 実際の材料設計への応用 ===")
    print("この情報から以下が予測できます：")
    print("1. どの相が析出するか（相分率から判断）")
    print("2. 各相の組成（偏析の程度）")
    print("3. 熱処理後の組織（相の安定性）")
    print("4. 機械的性質への影響（相の種類と分率）")

except Exception as e:
    print(f"計算エラー: {e}")
    print("デモ用TDBでは詳細な計算ができない場合があります。")

# 温度を変えて相分率の変化を追跡
print("\n\n=== 温度を変化させた相分率の変化 ===\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))

    # 可視化
    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("\n相分率グラフを保存しました: phase_fraction_vs_temperature.png")

5.2 レバー則による相分率計算

二相共存領域では、レバー則（lever rule）を使って各相の量を求めることができます。 pycalphadは自動的にこれを計算しますが、原理を理解することが重要です。

6. 駆動力（Driving Force）の計算

6.1 駆動力とは

駆動力（driving force）は、ある相から別の相への変態がどれだけ進行しやすいかを 示す熱力学的指標です。化学ポテンシャルの差として定義されます。

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

print("=== 相変態駆動力の計算 ===\n")

# demo_tdbを使用
from io import StringIO
db = Database(StringIO(demo_tdb))

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

# 計算条件
cu_composition = 0.2  # Cu 20 at.%
temperatures = np.linspace(700, 1200, 30)

print(f"組成: Al-{cu_composition*100:.0f}at.%Cu")
print(f"温度範囲: {temperatures[0]:.0f} - {temperatures[-1]:.0f} K")
print()

# 各温度での各相のギブスエネルギーを計算
fcc_energies = []
liquid_energies = []
driving_forces = []

for T in temperatures:
    try:
        # FCC相のギブスエネルギー
        eq_fcc = equilibrium(
            db, components, ['FCC_A1'],  # FCC相のみ
            {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)

        # LIQUID相のギブスエネルギー
        eq_liquid = equilibrium(
            db, components, ['LIQUID'],  # LIQUID相のみ
            {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)

        # 駆動力：LIQUID - FCC（正なら液相が安定）
        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("計算完了！\n")

# 結果の可視化
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 10))

# サブプロット1：ギブスエネルギー
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)

# サブプロット2：駆動力
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.png\n")

# 融点の推定（駆動力がゼロになる温度）
try:
    # ゼロクロス点を探す
    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:  # 符号が変わる
                T_melting = (temperatures[i] + temperatures[i+1]) / 2
                print(f"推定融点: {T_melting:.1f} K ({T_melting - 273.15:.1f}°C)\n")
                break
except:
    pass

print("=== 駆動力計算の応用 ===")
print("1. 相変態の熱力学的実行可能性の評価")
print("2. 過冷却・過加熱の程度の定量化")
print("3. 核生成・成長速度の推定（速度論と組み合わせ）")
print("4. 熱処理条件の最適化")
print("5. 準安定相の安定性評価")

6.2 駆動力と相変態速度の関係

駆動力は相変態の「熱力学的な推進力」を表しますが、実際の変態速度は 速度論（拡散、界面移動など）にも依存します。 両者を組み合わせた解析が、実際の材料プロセス設計には不可欠です。

7. 実践的ワークフローと応用例

7.1 完全な研究ワークフロー

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 実践的ワークフロー統合スクリプト")
print("=" * 60)
print()

# ============================================
# ステップ1: プロジェクト設定
# ============================================
print("【ステップ1】プロジェクト設定\n")

project_name = "Al-Cu Alloy Design"
target_composition = "Al-20at.%Cu"
target_temperature_range = (600, 1200)  # K
objective = "最適熱処理条件の決定"

print(f"プロジェクト名: {project_name}")
print(f"対象合金: {target_composition}")
print(f"温度範囲: {target_temperature_range[0]} - {target_temperature_range[1]} K")
print(f"目的: {objective}")
print()

# ============================================
# ステップ2: データベース準備
# ============================================
print("【ステップ2】TDBファイルの読み込み\n")

# 実際の研究では外部TDBファイルを使用
# db = Database('path/to/your/database.tdb')

# デモ用
from io import StringIO
db = Database(StringIO(demo_tdb))

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

print(f"使用するTDB: demo_al_cu.tdb")
print(f"元素: {[c for c in components if c != 'VA']}")
print(f"相: {phases}")
print()

# ============================================
# ステップ3: 相図計算
# ============================================
print("【ステップ3】相図計算\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"計算中... (温度{len(T_range)}点 × 組成{len(X_Cu_range)}点)")

try:
    eq_phase_diagram = equilibrium(
        db, components, phases,
        {v.X('CU'): X_Cu_range, v.T: T_range, v.P: 101325},
        output='GM'
    )
    print("✓ 相図計算完了")
except Exception as e:
    print(f"✗ 相図計算エラー: {e}")
    eq_phase_diagram = None

print()

# ============================================
# ステップ4: 特定条件での平衡計算
# ============================================
print("【ステップ4】特定条件での平衡計算\n")

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

print(f"条件: {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("\n安定相:")
    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("✓ 平衡計算完了")
except Exception as e:
    print(f"✗ 平衡計算エラー: {e}")
    eq_specific = None

print()

# ============================================
# ステップ5: 駆動力解析
# ============================================
print("【ステップ5】駆動力解析\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("✓ 駆動力解析完了")
else:
    print("✗ 駆動力解析エラー")

print()

# ============================================
# ステップ6: 結果の可視化
# ============================================
print("【ステップ6】結果の可視化\n")

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

# (1) 相図（簡易版）
ax1 = axes[0, 0]
if eq_phase_diagram is not None:
    # 相図の描画（実際にはより詳細な処理が必要）
    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) 相分率 vs 温度
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:
    # ダミーデータ
    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) 駆動力
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) 材料設計指針
ax4 = axes[1, 1]
ax4.axis('off')
design_guidelines = f"""
材料設計指針

【最適熱処理条件】
• 溶体化処理: 950-1000 K
• 時効処理: 450-500 K
• 冷却速度: 急冷推奨

【予想される析出相】
• FCC_A1 (α-Al matrix)
• θ相 (Al₂Cu, 析出強化)

【機械的性質予測】
• 高強度化: θ相析出により
• 延性: マトリックスの性質に依存

【次のステップ】
1. 実験での検証
2. 微細組織観察（SEM/TEM）
3. 機械試験（引張、硬さ）
"""
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("✓ 総合解析図を保存: workflow_comprehensive_analysis.png")
print()

# ============================================
# ステップ7: レポート生成
# ============================================
print("【ステップ7】レポート生成\n")

report = f"""
{'='*60}
pycalphad 計算レポート
{'='*60}

プロジェクト名: {project_name}
実行日時: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}

【計算条件】
- 対象系: {target_composition}
- 温度範囲: {target_temperature_range[0]} - {target_temperature_range[1]} K
- 圧力: 101325 Pa (1 atm)

【使用データベース】
- TDB: demo_al_cu.tdb
- 相: {', '.join(phases)}

【計算結果サマリー】
- 相図計算: {'成功' if eq_phase_diagram is not None else '失敗'}
- 平衡計算: {'成功' if eq_specific is not None else '失敗'}
- 駆動力解析: {'成功' if driving_forces_analysis else '失敗'}

【材料設計への提言】
{objective}に向けて、以下の条件を推奨します：
1. 溶体化処理: 950-1000 K で完全溶解
2. 急冷: 室温付近の準安定状態を保持
3. 時効処理: 450-500 K で析出強化相を生成

【参考文献】
[1] Pycalphad Documentation: https://pycalphad.org/
[2] CALPHAD methodology review papers
[3] Al-Cu binary phase diagram experimental data

{'='*60}
"""

print(report)

# レポートをファイルに保存
with open('pycalphad_calculation_report.txt', 'w', encoding='utf-8') as f:
    f.write(report)

print("✓ レポートを保存: pycalphad_calculation_report.txt")
print()

print("=" * 60)
print("ワークフロー完了！")
print("=" * 60)

7.2 ベストプラクティス

8. 演習問題

まとめと次のステップ

次に学ぶべき内容

学習の継続のために

最終メッセージ

材料熱力学は、材料科学の最も基盤的な学問領域の一つです。 「なぜこの相が安定なのか」「どの条件でどの相が共存するのか」という 本質的な問いに答えるための強力なツールです。

本シリーズでは、基礎理論からpycalphadを使った実践的計算まで、 体系的に学習しました。特に最終章では、TDBファイルを使った リアルな研究ワークフローを体験しました。

ここで習得した知識とスキルは、合金設計、プロセス最適化、 新材料探索、Materials Informaticsなど、あらゆる材料研究の基盤となります。

これからは、興味のある材料系や研究テーマに対して、 このシリーズで学んだ技術を積極的に応用してください。 pycalphadは強力なツールであり、あなたの研究を大きく加速させる可能性を秘めています。

継続的な学習と実践が、真のスキル習得への道です。
本シリーズが、あなたの材料科学研究の新たなスタートとなることを願っています。

🎓 材料熱力学入門シリーズ 完 🎓