Chapter 1: 金属結合と結晶構造

import numpy as np

# 物理定数
e = 1.602e-19  # 電子の電荷 (C)
m_e = 9.109e-31  # 電子の質量 (kg)

# 銅のパラメータ
n = 8.45e28  # 電子密度 (m^-3)
tau = 2.5e-14  # 緩和時間 (s)

# 電気伝導度の計算
sigma = (n * e**2 * tau) / m_e

# 電気抵抗率の計算
rho = 1 / sigma

print("=" * 60)
print("銅の電気伝導度と電気抵抗率（Drudeモデル）")
print("=" * 60)
print(f"\n与えられた値:")
print(f"  電子密度 n:    {n:.2e} m^-3")
print(f"  緩和時間 τ:    {tau:.2e} s")
print(f"\n計算結果:")
print(f"  電気伝導度 σ:  {sigma:.3e} S/m")
print(f"  電気抵抗率 ρ:  {rho:.3e} Ω·m")
print(f"                 {rho * 1e8:.2f} μΩ·cm")
print(f"\n参考: 実測値は約 1.68 μΩ·cm（室温）")
print("=" * 60)

# 説明
print("\n【計算の説明】")
print("-" * 60)
print("Drudeモデルの電気伝導度の式:")
print("  σ = (n × e² × τ) / m")
print("\nここで、")
print("  n  : 電子密度")
print("  e  : 電子の電荷 (1.602 × 10⁻¹⁹ C)")
print("  τ  : 緩和時間（電子が散乱されるまでの平均時間）")
print("  m  : 電子の質量 (9.109 × 10⁻³¹ kg)")
print("\n電気抵抗率は電気伝導度の逆数:")
print("  ρ = 1 / σ")
print("-" * 60)

# 検証
print("\n【結果の検証】")
print("-" * 60)
theoretical = rho * 1e8
experimental = 1.68
error = abs(theoretical - experimental) / experimental * 100
print(f"理論値（Drudeモデル）: {theoretical:.2f} μΩ·cm")
print(f"実測値（室温）:         {experimental:.2f} μΩ·cm")
print(f"誤差:                   {error:.1f}%")
print("\nDrudeモデルは古典的な近似ですが、")
print("オーダーは合っており、金属の電気伝導性を")
print("定性的に説明できることがわかります。")
print("=" * 60)

import numpy as np

# 格子定数
a = 4.05  # Å

print("=" * 70)
print("FCC構造の計算問題")
print("=" * 70)
print(f"\n与えられた値:")
print(f"  格子定数 a = {a} Å")

# (1) 原子半径の計算
# FCC では、面対角線上に原子が接触している
# 面対角線の長さ = 4r = a√2
r = a / (2 * np.sqrt(2))

print(f"\n(1) 原子半径の計算:")
print(f"    FCC構造では、面対角線上に原子が接触しています。")
print(f"    面対角線の長さ = 4r = a√2")
print(f"    したがって、r = a / (2√2)")
print(f"    ")
print(f"    計算: r = {a} / (2√2)")
print(f"           = {a} / {2 * np.sqrt(2):.4f}")
print(f"           = {r:.4f} Å")

# (2) 充填率の計算
# 単位格子あたりの原子数
n_atoms = 4

# 単位格子の体積
V_cell = a**3

# 原子の体積（球として）
V_atom = (4/3) * np.pi * r**3

# 充填率
packing_fraction = (n_atoms * V_atom) / V_cell

print(f"\n(2) 充填率の計算:")
print(f"    単位格子あたりの原子数: n = {n_atoms}")
print(f"    単位格子の体積: V_cell = a³ = {a}³ = {V_cell:.3f} ų")
print(f"    原子1個の体積: V_atom = (4/3)πr³ = {V_atom:.3f} ų")
print(f"    ")
print(f"    充填率 = (n × V_atom) / V_cell")
print(f"           = ({n_atoms} × {V_atom:.3f}) / {V_cell:.3f}")
print(f"           = {packing_fraction:.4f}")
print(f"           = {packing_fraction * 100:.2f}%")

# 理論値との比較
theoretical_packing = np.pi / (3 * np.sqrt(2))
print(f"    ")
print(f"    理論値: π / (3√2) = {theoretical_packing:.4f} = {theoretical_packing * 100:.2f}%")

# (3) 単位格子あたりの原子数
print(f"\n(3) 単位格子あたりの原子数:")
print(f"    頂点の原子: 8個の頂点 × 1/8 = 1個")
print(f"    面心の原子: 6個の面 × 1/2 = 3個")
print(f"    合計:       1 + 3 = {n_atoms}個")

# まとめ
print("\n" + "=" * 70)
print("【解答まとめ】")
print("=" * 70)
print(f"(1) 原子半径 r:              {r:.4f} Å = {r * 100:.2f} pm")
print(f"(2) 充填率:                  {packing_fraction:.4f} = {packing_fraction * 100:.2f}%")
print(f"(3) 単位格子あたりの原子数:  {n_atoms}個")
print("=" * 70)

# 補足情報
print("\n【補足】")
print("-" * 70)
print("FCC構造の特徴:")
print("  - 最密充填構造の一つ（充填率 74%）")
print("  - 配位数: 12")
print("  - すべり系が多く、延性が高い")
print("  - 代表的な金属: Cu, Al, Au, Ag, Ni, Pb")
print("-" * 70)

import numpy as np
import matplotlib.pyplot as plt

def calculate_d_spacing(a, h, k, l):
    """
    立方晶系の面間隔を計算

    Parameters:
    -----------
    a : float
        格子定数 (Å)
    h, k, l : int
        ミラー指数

    Returns:
    --------
    d : float
        面間隔 (Å)
    """
    d = a / np.sqrt(h**2 + k**2 + l**2)
    return d

# 格子定数
a = 3.5  # Å

# 計算する面
planes = [
    (1, 0, 0),
    (1, 1, 0),
    (1, 1, 1),
    (2, 0, 0),
    (2, 2, 0)
]

print("=" * 80)
print("立方晶系における面間隔の計算")
print("=" * 80)
print(f"\n与えられた値:")
print(f"  格子定数 a = {a} Å")
print(f"\n立方晶系の面間隔の公式:")
print(f"  d_hkl = a / √(h² + k² + l²)")

print(f"\n計算結果:")
print("=" * 80)
print(f"{'(hkl)':<8} {'h²+k²+l²':<12} {'√(h²+k²+l²)':<18} {'d_hkl (Å)':<15}")
print("-" * 80)

d_values = []
plane_labels = []

for h, k, l in planes:
    sum_squares = h**2 + k**2 + l**2
    sqrt_sum = np.sqrt(sum_squares)
    d = calculate_d_spacing(a, h, k, l)
    d_values.append(d)
    plane_labels.append(f"({h}{k}{l})")

    print(f"({h}{k}{l}){'':<5} {sum_squares:<12} {sqrt_sum:<18.4f} {d:.4f}")

print("=" * 80)

# 面間隔の可視化
fig, ax = plt.subplots(figsize=(12, 7))

colors = ['#FF6B6B', '#4ECDC4', '#45B7D1', '#FFA07A', '#98D8C8']
bars = ax.bar(plane_labels, d_values, color=colors, alpha=0.8, edgecolor='black', linewidth=2)

ax.set_ylabel('面間隔 d (Å)', fontsize=13, fontweight='bold')
ax.set_xlabel('ミラー指数 (hkl)', fontsize=13, fontweight='bold')
ax.set_title(f'立方晶系の面間隔（a = {a} Å）', fontsize=15, fontweight='bold', pad=15)
ax.grid(axis='y', alpha=0.3, linestyle='--')

# バーに数値を表示
for bar, d in zip(bars, d_values):
    height = bar.get_height()
    ax.text(bar.get_x() + bar.get_width()/2., height + 0.05,
            f'{d:.3f} Å', ha='center', va='bottom', fontsize=11, fontweight='bold')

plt.tight_layout()
plt.savefig('d_spacing_comparison.png', dpi=300, bbox_inches='tight')
plt.show()

print("\n✓ グラフを 'd_spacing_comparison.png' として保存しました")

# X線回折への応用
print("\n【X線回折への応用】")
print("=" * 80)
print("Braggの法則: 2d sinθ = nλ")
print(f"\nCu Kα線（λ = 1.5418 Å）を使った場合の回折角 2θ:")
print("-" * 80)

lambda_cu = 1.5418  # Å

for (h, k, l), d in zip(planes, d_values):
    if d > lambda_cu / 2:  # Bragg条件を満たす場合のみ
        sin_theta = lambda_cu / (2 * d)
        theta = np.arcsin(sin_theta) * 180 / np.pi
        two_theta = 2 * theta
        print(f"({h}{k}{l}):  d = {d:.4f} Å  →  2θ = {two_theta:.2f}°")
    else:
        print(f"({h}{k}{l}):  d = {d:.4f} Å  →  回折なし（λ/2 より小さい）")

print("=" * 80)

# 面間隔の順序
print("\n【面間隔の順序】")
print("-" * 80)
sorted_indices = np.argsort(d_values)[::-1]  # 降順
print("面間隔の大きい順:")
for i, idx in enumerate(sorted_indices):
    h, k, l = planes[idx]
    print(f"  {i+1}. ({h}{k}{l}): {d_values[idx]:.4f} Å")
print("-" * 80)

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

# 物理定数
k_B = constants.k  # ボルツマン定数 (J/K)
e = constants.e  # 電子の電荷 (C)
m_e = constants.m_e  # 電子の質量 (kg)

# 銅のフェルミエネルギー
E_F_eV = 7.05  # eV
E_F = E_F_eV * e  # J に変換

print("=" * 80)
print("銅のフェルミ-ディラック分布とフェルミ速度の計算")
print("=" * 80)
print(f"\n与えられた値:")
print(f"  フェルミエネルギー E_F = {E_F_eV} eV = {E_F:.3e} J")

# (1) フェルミ-ディラック分布の計算
def fermi_dirac(E, E_F, T):
    """
    フェルミ-ディラック分布関数

    Parameters:
    -----------
    E : array-like
        エネルギー (J)
    E_F : float
        フェルミエネルギー (J)
    T : float
        温度 (K)

    Returns:
    --------
    f : array-like
        占有確率
    """
    if T == 0:
        return np.where(E <= E_F, 1.0, 0.0)
    else:
        # オーバーフロー対策
        x = (E - E_F) / (k_B * T)
        x = np.clip(x, -100, 100)
        return 1 / (1 + np.exp(x))

# エネルギー範囲（eV単位で設定）
E_range_eV = np.linspace(0, 15, 1000)
E_range = E_range_eV * e  # J に変換

# 0 K と 300 K での分布
f_0K = fermi_dirac(E_range, E_F, 0)
f_300K = fermi_dirac(E_range, E_F, 300)

# プロット
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# 左図: 分布関数の比較
ax1.plot(E_range_eV, f_0K, 'b-', linewidth=2.5, label='T = 0 K')
ax1.plot(E_range_eV, f_300K, 'r-', linewidth=2.5, label='T = 300 K')
ax1.axvline(E_F_eV, color='green', linestyle='--', linewidth=2, label=f'$E_F$ = {E_F_eV} eV')
ax1.fill_between(E_range_eV, 0, f_300K, alpha=0.2, color='red')
ax1.set_xlabel('エネルギー (eV)', fontsize=12, fontweight='bold')
ax1.set_ylabel('占有確率 f(E)', fontsize=12, fontweight='bold')
ax1.set_title('フェルミ-ディラック分布の温度依存性', fontsize=13, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, 15)
ax1.set_ylim(-0.05, 1.05)

# 右図: フェルミエネルギー付近の拡大図
E_zoom_eV = np.linspace(E_F_eV - 1, E_F_eV + 1, 500)
E_zoom = E_zoom_eV * e
f_zoom_0K = fermi_dirac(E_zoom, E_F, 0)
f_zoom_300K = fermi_dirac(E_zoom, E_F, 300)

ax2.plot(E_zoom_eV, f_zoom_0K, 'b-', linewidth=2.5, label='T = 0 K')
ax2.plot(E_zoom_eV, f_zoom_300K, 'r-', linewidth=2.5, label='T = 300 K')
ax2.axvline(E_F_eV, color='green', linestyle='--', linewidth=2, label=f'$E_F$ = {E_F_eV} eV')
ax2.axhline(0.5, color='gray', linestyle=':', linewidth=1.5, alpha=0.7)
ax2.set_xlabel('エネルギー (eV)', fontsize=12, fontweight='bold')
ax2.set_ylabel('占有確率 f(E)', fontsize=12, fontweight='bold')
ax2.set_title('フェルミエネルギー付近の拡大図', fontsize=13, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)
ax2.set_xlim(E_F_eV - 1, E_F_eV + 1)
ax2.set_ylim(-0.05, 1.05)

plt.tight_layout()
plt.savefig('fermi_dirac_temperature_comparison.png', dpi=300, bbox_inches='tight')
plt.show()

print("\n(1) フェルミ-ディラック分布:")
print("-" * 80)
print(f"  0 K:   E < E_F ではすべての状態が占有 (f = 1)")
print(f"         E > E_F ではすべての状態が空 (f = 0)")
print(f"         ステップ関数として振る舞う")
print(f"  ")
print(f"  300 K: E = E_F で f = 0.5")
print(f"         熱エネルギー k_B T = {k_B * 300 / e * 1000:.2f} meV 程度の幅で")
print(f"         分布が「ぼける」")
print(f"         しかし、E_F = {E_F_eV} eV >> k_B T なので、")
print(f"         ステップ関数に近い形を保つ")

# (2) フェルミ速度の計算
# E_F = (1/2) m v_F^2 より v_F = √(2 E_F / m)
v_F = np.sqrt(2 * E_F / m_e)

print(f"\n(2) フェルミ速度の計算:")
print("-" * 80)
print(f"  フェルミエネルギーと電子の運動エネルギーの関係:")
print(f"    E_F = (1/2) m v_F²")
print(f"  ")
print(f"  したがって、フェルミ速度は:")
print(f"    v_F = √(2 E_F / m)")
print(f"  ")
print(f"  計算:")
print(f"    v_F = √(2 × {E_F:.3e} J / {m_e:.3e} kg)")
print(f"        = {v_F:.3e} m/s")
print(f"        = {v_F / 1e6:.3f} × 10⁶ m/s")
print(f"  ")
print(f"  光速 c = {constants.c:.3e} m/s との比:")
print(f"    v_F / c = {v_F / constants.c:.4f} = {v_F / constants.c * 100:.2f}%")
print("-" * 80)

print("\n" + "=" * 80)
print("【解答まとめ】")
print("=" * 80)
print(f"(1) フェルミ-ディラック分布:")
print(f"    - 0 K ではステップ関数")
print(f"    - 300 K では E_F 付近で k_B T 程度の幅でぼける")
print(f"    - グラフを 'fermi_dirac_temperature_comparison.png' に保存")
print(f"")
print(f"(2) フェルミ速度:")
print(f"    v_F = {v_F / 1e6:.3f} × 10⁶ m/s （光速の約 {v_F / constants.c * 100:.1f}%）")
print("=" * 80)

print("\n【物理的意義】")
print("-" * 80)
print("フェルミ速度は、絶対零度においてフェルミエネルギーを持つ")
print("電子の速度です。この速度は非常に大きく（約 10⁶ m/s）、")
print("金属中の電子が高速で運動していることを示しています。")
print("これが金属の高い電気伝導性と熱伝導性の原因です。")
print("-" * 80)

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

print("=" * 80)
print("FCC金属のすべり系解析：(111)面 vs (100)面")
print("=" * 80)

# FCC構造のすべり系
print("\n【FCC構造のすべり系】")
print("-" * 80)
print("FCC金属（Cu, Al, Au, Ag, Niなど）の主要なすべり系:")
print("  すべり面: {111} 面（最密充填面）")
print("  すべり方向: <110> 方向（最密充填方向）")
print("  ")
print("(111)面のすべり系: 4つの{111}面 × 3つの<110>方向 = 12個")
print("(100)面のすべり系: 通常はすべらないが、高温では活性化")

# (111)面と(100)面の原子配列密度の計算
print("\n【原子配列密度の比較】")
print("-" * 80)

# 格子定数
a = 1.0  # 任意単位

# (111)面の原子配列密度
# (111)面は正三角形の格子
# 三角形の面積 = (√3/2) × (a/√2)² = (√3/4) × a²
area_111_triangle = (np.sqrt(3) / 4) * a**2

# (111)面の単位三角形には1個の原子が含まれる
# （頂点の3個 × 1/3 = 1個）
atoms_per_triangle_111 = 1

# 原子配列密度 = 原子数 / 面積
density_111 = atoms_per_triangle_111 / area_111_triangle

print(f"\n(111)面の解析:")
print(f"  格子構造: 正三角形格子（最密充填）")
print(f"  単位三角形の面積: (√3/4) × a² = {area_111_triangle:.4f} a²")
print(f"  単位三角形あたりの原子数: {atoms_per_triangle_111}")
print(f"  原子配列密度: {density_111:.4f} / a²")
print(f"                = {density_111:.4f} × (1/a²)")

# (100)面の原子配列密度
# (100)面は正方格子
# 正方形の面積 = (a/√2)² = a²/2
area_100_square = a**2 / 2

# (100)面の単位正方形には1個の原子が含まれる
# （頂点の4個 × 1/4 = 1個）
atoms_per_square_100 = 1

# 原子配列密度
density_100 = atoms_per_square_100 / area_100_square

print(f"\n(100)面の解析:")
print(f"  格子構造: 正方格子")
print(f"  単位正方形の面積: a²/2 = {area_100_square:.4f} a²")
print(f"  単位正方形あたりの原子数: {atoms_per_square_100}")
print(f"  原子配列密度: {density_100:.4f} / a²")
print(f"                = {density_100:.4f} × (1/a²)")

# 比較
ratio = density_111 / density_100
print(f"\n原子配列密度の比:")
print(f"  (111)面 / (100)面 = {ratio:.4f}")
print(f"  (111)面は(100)面の {ratio:.2f}倍 密に原子が配列")

# すべりやすさの評価
print("\n" + "=" * 80)
print("【すべりやすさの評価】")
print("=" * 80)

print("\n(111)面が主要なすべり面である理由:")
print("-" * 80)
print(f"1. 最密充填面である")
print(f"   → 原子配列密度が最も高い ({density_111:.4f} / a²)")
print(f"   → 原子間結合が強く、面全体として安定")
print(f"")
print(f"2. 面間隔が最も大きい")
print(f"   → d_111 = a / √3 = {a / np.sqrt(3):.4f} a")
print(f"   → d_100 = a / √1 = {a / np.sqrt(1):.4f} a")
print(f"   → 面間隔が大きいほど、すべりに必要なせん断応力が小さい")
print(f"")
print(f"3. すべり方向 <110> が最密充填方向")
print(f"   → 原子間距離が最短 (a/√2 = {a / np.sqrt(2):.4f} a)")
print(f"   → すべりの際の原子の移動距離が最小")
print(f"")
print(f"4. すべり系の数が多い")
print(f"   → {111} 面: 4面 × 3方向 = 12個のすべり系")
print(f"   → 多様な応力方向に対応可能 → 延性が高い")

print("\n(100)面がすべりにくい理由:")
print("-" * 80)
print(f"1. 原子配列密度が低い ({density_100:.4f} / a²)")
print(f"   → (111)面の {1/ratio:.2f}倍")
print(f"   → 面全体としての結合が弱い")
print(f"")
print(f"2. 面間隔が小さい")
print(f"   → d_100 = {a / np.sqrt(1):.4f} a （(111)面の {np.sqrt(3):.2f}倍密）")
print(f"   → すべりに必要なせん断応力が大きい")
print(f"")
print(f"3. 常温では活性化しない")
print(f"   → 高温（>0.5T_m）で初めて活性化")

# 臨界分解せん断応力（CRSS）の概念
print("\n" + "=" * 80)
print("【臨界分解せん断応力（CRSS）】")
print("=" * 80)

print("\nSchmidの法則:")
print("  τ_R = σ cos φ cos λ")
print("  ")
print("  ここで、")
print("  τ_R : 分解せん断応力")
print("  σ  : 引張応力")
print("  φ  : 応力軸とすべり面法線のなす角")
print("  λ  : 応力軸とすべり方向のなす角")
print("  ")
print("すべりが開始する条件:")
print("  τ_R ≥ τ_CRSS")
print("  ")
print("FCC金属の典型的なCRSS値（室温）:")
print("  Cu:  約 0.5 - 1.0 MPa")
print("  Al:  約 0.2 - 0.5 MPa")
print("  Au:  約 0.3 - 0.8 MPa")

# 可視化：すべり系の模式図
fig = plt.figure(figsize=(14, 6))

# 左図: (111)面のすべり系
ax1 = fig.add_subplot(121, projection='3d')
ax1.set_title('(111)面のすべり系\n（最密充填面、12個のすべり系）',
              fontsize=12, fontweight='bold', pad=15)

# (111)面を描画（簡略化）
plane_vertices_111 = np.array([
    [0, 0, 1],
    [1, 0, 0],
    [0, 1, 0]
])

# 面を描画
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
poly_111 = Poly3DCollection([plane_vertices_111], alpha=0.5,
                             facecolor='cyan', edgecolor='blue', linewidth=2)
ax1.add_collection3d(poly_111)

# すべり方向のベクトルを描画
slip_directions_111 = [
    ([0, 0, 1], [1, 0, 0]),
    ([1, 0, 0], [0, 1, 0]),
    ([0, 1, 0], [0, 0, 1])
]

for start, end in slip_directions_111:
    start = np.array(start)
    end = np.array(end)
    ax1.quiver(start[0], start[1], start[2],
               end[0] - start[0], end[1] - start[1], end[2] - start[2],
               color='red', arrow_length_ratio=0.2, linewidth=2)

ax1.set_xlabel('X', fontsize=10, fontweight='bold')
ax1.set_ylabel('Y', fontsize=10, fontweight='bold')
ax1.set_zlabel('Z', fontsize=10, fontweight='bold')
ax1.set_xlim([0, 1])
ax1.set_ylim([0, 1])
ax1.set_zlim([0, 1])

# 右図: (100)面のすべり系
ax2 = fig.add_subplot(122, projection='3d')
ax2.set_title('(100)面のすべり系\n（非最密充填面、常温では不活性）',
              fontsize=12, fontweight='bold', pad=15)

# (100)面を描画
plane_vertices_100 = np.array([
    [1, 0, 0],
    [1, 1, 0],
    [1, 1, 1],
    [1, 0, 1]
])

poly_100 = Poly3DCollection([plane_vertices_100], alpha=0.5,
                             facecolor='lightcoral', edgecolor='red', linewidth=2)
ax2.add_collection3d(poly_100)

ax2.set_xlabel('X', fontsize=10, fontweight='bold')
ax2.set_ylabel('Y', fontsize=10, fontweight='bold')
ax2.set_zlabel('Z', fontsize=10, fontweight='bold')
ax2.set_xlim([0, 1])
ax2.set_ylim([0, 1])
ax2.set_zlim([0, 1])

plt.tight_layout()
plt.savefig('fcc_slip_systems_comparison.png', dpi=300, bbox_inches='tight')
plt.show()

print("\n" + "=" * 80)
print("【結論】")
print("=" * 80)
print("(111)面が(100)面よりすべりやすい理由:")
print("  1. 最密充填面であり、原子配列密度が高い")
print(f"     → (111): {density_111:.4f} / a²")
print(f"     → (100): {density_100:.4f} / a² （{ratio:.2f}倍の差）")
print("  2. 面間隔が大きく、すべりに必要なせん断応力が小さい")
print("  3. すべり系が12個と多く、多様な応力に対応可能")
print("  4. FCC金属の高い延性の主要因")
print("=" * 80)

print("\n✓ グラフを 'fcc_slip_systems_comparison.png' として保存しました")