第5章: フォノンエンジニアリング

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh

def phononic_band_structure_1D(N_layers=100, N_k=200, a=1.0,
                                rho_A=1.0, rho_B=2.0,
                                c_A=1.0, c_B=0.5):
    """
    1次元フォノニック結晶（2層周期）のバンド構造計算

    Parameters:
    -----------
    N_layers : int
        計算に使う周期の数
    N_k : int
        k点の数
    a : float
        格子定数（1周期の長さ）
    rho_A, rho_B : float
        材料AとBの密度
    c_A, c_B : float
        材料AとBの弾性定数

    Returns:
    --------
    k_points : array
        波数ベクトル
    omega : array
        角周波数（バンド）
    """
    # 各層の厚さ（等しいと仮定）
    d_A = a / 2
    d_B = a / 2

    # 波数範囲（第1ブリルアンゾーン）
    k_points = np.linspace(-np.pi/a, np.pi/a, N_k)

    # 平面波基底の数
    N_G = 2 * N_layers + 1
    G_indices = np.arange(-N_layers, N_layers + 1)

    # バンドを格納する配列
    omega_bands = np.zeros((N_k, N_G))

    for i_k, k in enumerate(k_points):
        # 動的行列を構築
        D = np.zeros((N_G, N_G), dtype=complex)

        for i, G_i in enumerate(G_indices):
            for j, G_j in enumerate(G_indices):
                G_diff = G_i - G_j

                if G_diff == 0:
                    # 対角成分
                    # フーリエ係数 <1/rho> と 
                    inv_rho_avg = (d_A/rho_A + d_B/rho_B) / a
                    c_avg = (d_A*c_A + d_B*c_B) / a

                    D[i, j] = c_avg * (k + 2*np.pi*G_i/a)**2 / inv_rho_avg
                else:
                    # 非対角成分（フーリエ係数）
                    # 矩形波のフーリエ級数
                    inv_rho_G = (1/rho_A - 1/rho_B) * (d_A/a) * np.sinc(G_diff * d_A/a)
                    c_G = (c_A - c_B) * (d_A/a) * np.sinc(G_diff * d_A/a)

                    k_i = k + 2*np.pi*G_i/a
                    k_j = k + 2*np.pi*G_j/a

                    D[i, j] = c_G * k_i * k_j / inv_rho_avg

        # 固有値問題を解く
        eigenvalues = eigh(D, eigvals_only=True)
        omega_bands[i_k, :] = np.sqrt(np.abs(eigenvalues))

    return k_points, omega_bands

# 計算実行
k, omega = phononic_band_structure_1D(
    N_layers=50,
    N_k=300,
    rho_A=1.0, rho_B=4.0,  # 大きな密度対比
    c_A=1.0, c_B=0.3       # 弾性定数対比
)

# プロット
plt.figure(figsize=(10, 6))
for i in range(omega.shape[1]):
    plt.plot(k, omega[:, i], 'b-', linewidth=0.5)

plt.xlabel('Wave vector k (π/a)', fontsize=12)
plt.ylabel('Frequency ω (arb. units)', fontsize=12)
plt.title('1D Phononic Crystal Band Structure', fontsize=14)
plt.axvline(x=-1, color='r', linestyle='--', alpha=0.3, label='Brillouin zone')
plt.axvline(x=1, color='r', linestyle='--', alpha=0.3)
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.savefig('phononic_band_1d.png', dpi=300)
plt.show()

# バンドギャップの検出
def find_band_gaps(k, omega, k_point=0):
    """
    特定のk点でのバンドギャップを検出

    Parameters:
    -----------
    k : array
        波数ベクトル
    omega : array
        周波数バンド
    k_point : float
        ギャップを調べるk点（デフォルトはΓ点）

    Returns:
    --------
    gaps : list of tuples
        (下端周波数, 上端周波数, ギャップ幅)のリスト
    """
    # 指定k点に最も近いインデックス
    idx = np.argmin(np.abs(k - k_point))

    # その点での周波数を昇順ソート
    freqs = np.sort(omega[idx, :])

    # バンド間の隙間を検出
    gaps = []
    for i in range(len(freqs) - 1):
        gap_width = freqs[i+1] - freqs[i]
        if gap_width > 0.01:  # 閾値（数値誤差を除外）
            gaps.append((freqs[i], freqs[i+1], gap_width))

    return gaps

# Γ点（k=0）でのバンドギャップ
gaps_gamma = find_band_gaps(k, omega, k_point=0)
print("Band gaps at Γ point (k=0):")
for i, (lower, upper, width) in enumerate(gaps_gamma):
    print(f"  Gap {i+1}: [{lower:.3f}, {upper:.3f}], width = {width:.3f}")

# ブリルアンゾーン境界（k=π/a）でのバンドギャップ
gaps_edge = find_band_gaps(k, omega, k_point=np.pi)
print("\nBand gaps at zone edge (k=π/a):")
for i, (lower, upper, width) in enumerate(gaps_edge):
    print(f"  Gap {i+1}: [{lower:.3f}, {upper:.3f}], width = {width:.3f}")

import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import mean_absolute_error, r2_score
import matplotlib.pyplot as plt

# サンプルデータセット作成（実際はMaterials ProjectやCOD等から取得）
def create_sample_dataset(n_samples=500):
    """
    熱伝導率予測用のサンプルデータセット

    実際の応用では、Materials Project APIや文献データから構築
    """
    np.random.seed(42)

    # 記述子（特徴量）
    data = {
        # 元素特性
        'avg_atomic_mass': np.random.uniform(20, 200, n_samples),
        'mass_variance': np.random.uniform(0, 50, n_samples),
        'avg_atomic_radius': np.random.uniform(0.5, 2.0, n_samples),

        # 結合特性
        'avg_bond_strength': np.random.uniform(100, 1000, n_samples),  # kJ/mol
        'electronegativity_diff': np.random.uniform(0, 2.5, n_samples),

        # 構造特性
        'density': np.random.uniform(2, 20, n_samples),  # g/cm³
        'packing_fraction': np.random.uniform(0.4, 0.74, n_samples),
        'space_group_number': np.random.randint(1, 230, n_samples),

        # フォノン関連（簡易計算または既知）
        'debye_temperature': np.random.uniform(200, 1000, n_samples),  # K
        'gruneisen_parameter': np.random.uniform(1, 3, n_samples),
    }

    df = pd.DataFrame(data)

    # 熱伝導率の簡易的な生成（実際は実測値または第一原理計算値）
    # 物理的に妥当な依存性を模擬
    kappa_lattice = (
        300 / df['avg_atomic_mass']**0.5 *  # 軽い元素ほど高い
        df['avg_bond_strength']**0.3 /      # 強い結合ほど高い
        (df['mass_variance'] + 1) *         # 質量揺らぎで低下
        df['debye_temperature']**0.5 /      # Debye温度と相関
        (df['gruneisen_parameter']**2 + 1)  # 非調和性で低下
    )

    # ノイズ追加（実験誤差を模擬）
    kappa_lattice *= np.random.lognormal(0, 0.3, n_samples)

    df['kappa_lattice'] = kappa_lattice

    return df

# データセット作成
df = create_sample_dataset(n_samples=1000)

print("Dataset shape:", df.shape)
print("\nFeature statistics:")
print(df.describe())

# 特徴量と目的変数の分離
X = df.drop('kappa_lattice', axis=1)
y = df['kappa_lattice']

# 訓練・テストセット分割
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# ランダムフォレストモデルの訓練
rf_model = RandomForestRegressor(
    n_estimators=200,
    max_depth=15,
    min_samples_split=5,
    min_samples_leaf=2,
    random_state=42,
    n_jobs=-1
)

rf_model.fit(X_train, y_train)

# 予測
y_pred_train = rf_model.predict(X_train)
y_pred_test = rf_model.predict(X_test)

# 評価指標
mae_train = mean_absolute_error(y_train, y_pred_train)
mae_test = mean_absolute_error(y_test, y_pred_test)
r2_train = r2_score(y_train, y_pred_train)
r2_test = r2_score(y_test, y_pred_test)

print(f"\n=== Model Performance ===")
print(f"Training MAE: {mae_train:.2f} W/(m·K)")
print(f"Test MAE: {mae_test:.2f} W/(m·K)")
print(f"Training R²: {r2_train:.3f}")
print(f"Test R²: {r2_test:.3f}")

# クロスバリデーション
cv_scores = cross_val_score(
    rf_model, X, y, cv=5,
    scoring='neg_mean_absolute_error'
)
print(f"\n5-fold CV MAE: {-cv_scores.mean():.2f} ± {cv_scores.std():.2f} W/(m·K)")

# 特徴量重要度
feature_importance = pd.DataFrame({
    'feature': X.columns,
    'importance': rf_model.feature_importances_
}).sort_values('importance', ascending=False)

print("\n=== Feature Importance ===")
print(feature_importance)

# プロット1: 予測 vs 実測
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Training set
axes[0].scatter(y_train, y_pred_train, alpha=0.5, s=10)
axes[0].plot([y_train.min(), y_train.max()],
             [y_train.min(), y_train.max()],
             'r--', lw=2)
axes[0].set_xlabel('Actual κ (W/(m·K))', fontsize=12)
axes[0].set_ylabel('Predicted κ (W/(m·K))', fontsize=12)
axes[0].set_title(f'Training Set (R² = {r2_train:.3f})', fontsize=14)
axes[0].grid(True, alpha=0.3)

# Test set
axes[1].scatter(y_test, y_pred_test, alpha=0.5, s=10, color='green')
axes[1].plot([y_test.min(), y_test.max()],
             [y_test.min(), y_test.max()],
             'r--', lw=2)
axes[1].set_xlabel('Actual κ (W/(m·K))', fontsize=12)
axes[1].set_ylabel('Predicted κ (W/(m·K))', fontsize=12)
axes[1].set_title(f'Test Set (R² = {r2_test:.3f})', fontsize=14)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('kappa_prediction.png', dpi=300)
plt.show()

# プロット2: 特徴量重要度
plt.figure(figsize=(10, 6))
plt.barh(feature_importance['feature'], feature_importance['importance'])
plt.xlabel('Importance', fontsize=12)
plt.ylabel('Feature', fontsize=12)
plt.title('Feature Importance for Thermal Conductivity Prediction', fontsize=14)
plt.tight_layout()
plt.savefig('feature_importance.png', dpi=300)
plt.show()

# 新材料の予測例
def predict_new_material(model, features_dict):
    """
    新材料の熱伝導率を予測

    Parameters:
    -----------
    model : trained model
        訓練済みモデル
    features_dict : dict
        記述子の辞書

    Returns:
    --------
    predicted_kappa : float
        予測された格子熱伝導率
    """
    # DataFrameに変換
    features_df = pd.DataFrame([features_dict])

    # 予測
    kappa_pred = model.predict(features_df)[0]

    return kappa_pred

# 例：ダイヤモンド様の材料
diamond_like = {
    'avg_atomic_mass': 12.0,  # 軽い（C）
    'mass_variance': 0.0,     # 単元素
    'avg_atomic_radius': 0.77,
    'avg_bond_strength': 711, # 強い共有結合
    'electronegativity_diff': 0.0,
    'density': 3.5,
    'packing_fraction': 0.34,
    'space_group_number': 227,  # Fd-3m
    'debye_temperature': 2200,
    'gruneisen_parameter': 1.0,  # 低非調和性
}

kappa_diamond_pred = predict_new_material(rf_model, diamond_like)
print(f"\n予測: ダイヤモンド様材料の κ_L = {kappa_diamond_pred:.0f} W/(m·K)")
print(f"参考: 実際のダイヤモンドは ~2000 W/(m·K)")

# 例：重元素合金（低κ材料）
heavy_alloy = {
    'avg_atomic_mass': 150.0,  # 重い
    'mass_variance': 30.0,     # 大きな質量揺らぎ
    'avg_atomic_radius': 1.5,
    'avg_bond_strength': 200,  # 弱い結合
    'electronegativity_diff': 1.2,
    'density': 15.0,
    'packing_fraction': 0.68,
    'space_group_number': 225,
    'debye_temperature': 300,
    'gruneisen_parameter': 2.5,  # 高非調和性
}

kappa_heavy_pred = predict_new_material(rf_model, heavy_alloy)
print(f"\n予測: 重元素合金の κ_L = {kappa_heavy_pred:.1f} W/(m·K)")
print(f"参考: Bi₂Te₃は ~1.5 W/(m·K)")