第6章：実践：物性計算ワークフロー

from ase import Atoms
from ase.build import bulk
from ase.visualize import view
import numpy as np

# Si結晶構造作成（ダイヤモンド構造）
si = bulk('Si', 'diamond', a=5.43)

print("Si結晶情報:")
print(f"  格子定数: {si.cell.cellpar()[0]:.3f} Å")
print(f"  原子数: {len(si)} atoms")
print(f"  空間群: Fd-3m (227)")
print(f"  原子座標:")
for i, pos in enumerate(si.get_scaled_positions()):
    print(f"    Si{i+1}: ({pos[0]:.3f}, {pos[1]:.3f}, {pos[2]:.3f})")

# 2x2x2スーパーセル作成
si_supercell = si.repeat((2, 2, 2))
print(f"\n2x2x2スーパーセル: {len(si_supercell)} atoms")

# 構造の可視化（ASE GUIが利用可能な場合）
# view(si)

# CIF形式で保存
from ase.io import write
write('si_unitcell.cif', si)
write('si_supercell.cif', si_supercell)
print("\n✅ CIFファイル保存完了")

from ase.build import bulk
from ase.calculators.espresso import Espresso
import numpy as np
import matplotlib.pyplot as plt

# Si結晶作成
si = bulk('Si', 'diamond', a=5.43)

# 擬ポテンシャル設定
pseudopotentials = {'Si': 'Si.pbe-n-rrkjus_psl.1.0.0.UPF'}

# =====================================
# 1. k点収束テスト（エネルギーカットオフ固定）
# =====================================
kpts_list = [2, 4, 6, 8, 10, 12]
energies_k = []
ecutwfc_fixed = 30.0  # Ry

print("k点収束テスト実行中...")
for k in kpts_list:
    calc = Espresso(
        pseudopotentials=pseudopotentials,
        input_data={
            'control': {'calculation': 'scf'},
            'system': {'ecutwfc': ecutwfc_fixed, 'occupations': 'smearing'},
        },
        kpts=(k, k, k),
    )
    si.calc = calc
    energy = si.get_potential_energy()
    energies_k.append(energy)
    print(f"  k={k:2d}: E = {energy:.6f} eV")

# =====================================
# 2. エネルギーカットオフ収束テスト（k点固定）
# =====================================
ecutwfc_list = [20, 25, 30, 35, 40, 45, 50]
energies_ecut = []
kpts_fixed = (8, 8, 8)

print("\nエネルギーカットオフ収束テスト実行中...")
for ecut in ecutwfc_list:
    calc = Espresso(
        pseudopotentials=pseudopotentials,
        input_data={
            'control': {'calculation': 'scf'},
            'system': {'ecutwfc': ecut, 'occupations': 'smearing'},
        },
        kpts=kpts_fixed,
    )
    si.calc = calc
    energy = si.get_potential_energy()
    energies_ecut.append(energy)
    print(f"  ecutwfc={ecut:2d} Ry: E = {energy:.6f} eV")

# =====================================
# 3. 収束グラフ作成
# =====================================
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))

# k点収束
ax1.plot(kpts_list, energies_k, 'o-', linewidth=2, markersize=8)
ax1.axhline(y=energies_k[-1], color='r', linestyle='--',
            label=f'Converged: {energies_k[-1]:.4f} eV')
ax1.set_xlabel('k-point mesh (k×k×k)', fontsize=12)
ax1.set_ylabel('Total Energy (eV)', fontsize=12)
ax1.set_title('k-point Convergence Test', fontsize=14)
ax1.grid(True, alpha=0.3)
ax1.legend()

# エネルギーカットオフ収束
ax2.plot(ecutwfc_list, energies_ecut, 's-', linewidth=2, markersize=8)
ax2.axhline(y=energies_ecut[-1], color='r', linestyle='--',
            label=f'Converged: {energies_ecut[-1]:.4f} eV')
ax2.set_xlabel('Energy Cutoff (Ry)', fontsize=12)
ax2.set_ylabel('Total Energy (eV)', fontsize=12)
ax2.set_title('Energy Cutoff Convergence Test', fontsize=14)
ax2.grid(True, alpha=0.3)
ax2.legend()

plt.tight_layout()
plt.savefig('si_convergence_tests.png', dpi=300)
print("\n✅ 収束テストグラフ保存: si_convergence_tests.png")

# =====================================
# 4. 収束判定（エネルギー変化 < 1 meV/atom）
# =====================================
threshold = 0.001  # eV/atom

delta_k = abs(energies_k[-1] - energies_k[-2]) / len(si)
delta_ecut = abs(energies_ecut[-1] - energies_ecut[-2]) / len(si)

print("\n📊 収束判定:")
print(f"  k点: ΔE = {delta_k*1000:.3f} meV/atom → {'✅ 収束' if delta_k < threshold else '❌ 未収束'}")
print(f"  ecutwfc: ΔE = {delta_ecut*1000:.3f} meV/atom → {'✅ 収束' if delta_ecut < threshold else '❌ 未収束'}")

print(f"\n🎯 推奨設定:")
print(f"  k点メッシュ: {kpts_list[-2]}×{kpts_list[-2]}×{kpts_list[-2]}")
print(f"  エネルギーカットオフ: {ecutwfc_list[-2]} Ry")

from ase.build import bulk
from ase.calculators.espresso import Espresso
from ase.dft.kpoints import special_points, bandpath
import matplotlib.pyplot as plt
import numpy as np

# Si結晶作成
si = bulk('Si', 'diamond', a=5.43)

# 収束テストから得た最適パラメータ
optimal_kpts = (10, 10, 10)
optimal_ecutwfc = 40.0  # Ry

pseudopotentials = {'Si': 'Si.pbe-n-rrkjus_psl.1.0.0.UPF'}

# =====================================
# ステップ1：SCF計算
# =====================================
print("ステップ1：SCF計算実行中...")
calc_scf = Espresso(
    pseudopotentials=pseudopotentials,
    input_data={
        'control': {
            'calculation': 'scf',
            'restart_mode': 'from_scratch',
            'prefix': 'si',
        },
        'system': {
            'ecutwfc': optimal_ecutwfc,
            'occupations': 'smearing',
            'smearing': 'gaussian',
            'degauss': 0.01,
        },
        'electrons': {
            'conv_thr': 1.0e-8,
        }
    },
    kpts=optimal_kpts,
)
si.calc = calc_scf
energy_scf = si.get_potential_energy()
print(f"✅ SCF完了: E = {energy_scf:.6f} eV")

# =====================================
# ステップ2：バンド構造計算
# =====================================
print("\nステップ2：バンド構造計算実行中...")

# 高対称点パス定義（FCC: L-Γ-X-W-K-Γ）
points = special_points['fcc']
path = bandpath(['L', 'Gamma', 'X', 'W', 'K', 'Gamma'], si.cell, npoints=100)

calc_bands = Espresso(
    pseudopotentials=pseudopotentials,
    input_data={
        'control': {
            'calculation': 'bands',
            'restart_mode': 'restart',
            'prefix': 'si',
        },
        'system': {
            'ecutwfc': optimal_ecutwfc,
            'occupations': 'smearing',
            'smearing': 'gaussian',
            'degauss': 0.01,
        },
    },
    kpts=path,
)
si.calc = calc_bands
bands = calc_bands.band_structure()

print("✅ バンド構造計算完了")

# =====================================
# ステップ3：バンドギャップ計算
# =====================================
energies = bands.energies
kpoints = bands.kpts
fermi_level = calc_bands.get_fermi_level()

# 価電子帯頂上（VBM）と伝導帯底（CBM）を探索
vbm_energy = np.max(energies[energies < fermi_level])
cbm_energy = np.min(energies[energies > fermi_level])
bandgap = cbm_energy - vbm_energy

print(f"\n📊 バンド構造解析結果:")
print(f"  フェルミレベル: {fermi_level:.4f} eV")
print(f"  VBM: {vbm_energy:.4f} eV")
print(f"  CBM: {cbm_energy:.4f} eV")
print(f"  バンドギャップ: {bandgap:.4f} eV")
print(f"  実験値: 1.12 eV (300K)")
print(f"  誤差: {abs(bandgap - 1.12)/1.12*100:.1f}%")

# =====================================
# ステップ4：可視化
# =====================================
fig, ax = plt.subplots(figsize=(8, 6))
bands.plot(ax=ax, emin=-10, emax=10)
ax.axhline(y=0, color='k', linestyle='--', linewidth=1, label='Fermi level')
ax.set_ylabel('Energy (eV)', fontsize=12)
ax.set_title(f'Si Band Structure (Eg = {bandgap:.3f} eV)', fontsize=14)
ax.legend()
plt.tight_layout()
plt.savefig('si_bandstructure.png', dpi=300)
print("\n✅ バンド構造図保存: si_bandstructure.png")

from ase import Atoms
from ase.io import write
import numpy as np

# GaN ウルツ鉱構造（六方晶, P6_3mc）
a = 3.189  # Å
c = 5.185  # Å
u = 0.377  # 内部パラメータ

# 六方晶格子ベクトル
cell = [
    [a, 0, 0],
    [-a/2, a*np.sqrt(3)/2, 0],
    [0, 0, c]
]

# 原子座標（分率座標）
positions = [
    [1/3, 2/3, 0],      # Ga1
    [2/3, 1/3, 1/2],    # Ga2
    [1/3, 2/3, u],      # N1
    [2/3, 1/3, 1/2+u],  # N2
]

# Atoms作成
gan = Atoms(
    symbols='Ga2N2',
    scaled_positions=positions,
    cell=cell,
    pbc=True
)

print("GaN結晶情報:")
print(f"  格子定数: a = {a:.3f} Å, c = {c:.3f} Å")
print(f"  c/a比: {c/a:.3f}")
print(f"  空間群: P6_3mc (186)")
print(f"  原子数: {len(gan)} atoms")
print(f"\n原子座標（カルテシアン, Å）:")
for i, (sym, pos) in enumerate(zip(gan.get_chemical_symbols(), gan.positions)):
    print(f"  {sym}{i+1}: ({pos[0]:7.3f}, {pos[1]:7.3f}, {pos[2]:7.3f})")

# CIF保存
write('gan_wurtzite.cif', gan)
print("\n✅ GaN構造保存完了: gan_wurtzite.cif")

# =====================================
# 光学計算のための設定例
# =====================================
print("\n📋 GaN光学計算推奨設定:")
print("  1. 構造最適化:")
print("     - 格子定数最適化: 必須")
print("     - 内部座標最適化: u パラメータ")
print("  2. 収束設定:")
print("     - k点: 12×12×8 以上（六方晶系）")
print("     - ecutwfc: 50 Ry 以上")
print("  3. 光学計算:")
print("     - 空バンド数: 伝導帯 30本以上")
print("     - スメアリング: 0.01 Ry (Gaussian)")
print("     - 偏光方向: ⊥c軸（E_⊥）と ∥c軸（E_∥）")

from ase.calculators.espresso import Espresso
import numpy as np
import matplotlib.pyplot as plt

# GaN構造（前のコードで作成済み）
# gan = [GaN Atoms object]

pseudopotentials = {
    'Ga': 'Ga.pbe-dnl-rrkjus_psl.1.0.0.UPF',
    'N': 'N.pbe-n-radius_5.UPF'
}

# =====================================
# ステップ1：SCF計算（高密度k点）
# =====================================
print("ステップ1：SCF計算（光学用）...")
calc_scf = Espresso(
    pseudopotentials=pseudopotentials,
    input_data={
        'control': {
            'calculation': 'scf',
            'prefix': 'gan',
        },
        'system': {
            'ecutwfc': 50.0,
            'nbnd': 40,  # 空バンド含む
            'occupations': 'smearing',
            'smearing': 'gaussian',
            'degauss': 0.01,
        },
        'electrons': {
            'conv_thr': 1.0e-9,
        }
    },
    kpts=(12, 12, 8),  # 六方晶系
)
gan.calc = calc_scf
energy_scf = gan.get_potential_energy()
print(f"✅ SCF完了: E = {energy_scf:.6f} eV")

# =====================================
# ステップ2：NSCF計算（より密なk点）
# =====================================
print("\nステップ2：NSCF計算（光学遷移用）...")
calc_nscf = Espresso(
    pseudopotentials=pseudopotentials,
    input_data={
        'control': {
            'calculation': 'nscf',
            'restart_mode': 'restart',
            'prefix': 'gan',
        },
        'system': {
            'ecutwfc': 50.0,
            'nbnd': 60,  # さらに多くの空バンド
            'occupations': 'smearing',
            'smearing': 'gaussian',
            'degauss': 0.01,
        },
    },
    kpts=(24, 24, 16),  # SCFの2倍密度
)
gan.calc = calc_nscf
print("✅ NSCF完了")

# =====================================
# ステップ3：誘電関数計算（epsilon.x）
# =====================================
print("\nステップ3：誘電関数計算...")

# epsilon.xの入力ファイル作成（実際の計算ではQuantum ESPRESSOのepsilon.xを使用）
epsilon_input = """
 &inputpp
   prefix = 'gan'
   outdir = './tmp/'
   calculation = 'eps'
 /
 &energy_grid
   smeartype = 'gauss'
   intersmear = 0.1
   wmin = 0.0
   wmax = 30.0
   nw = 500
 /
"""

with open('epsilon.in', 'w') as f:
    f.write(epsilon_input)

print("✅ epsilon.x入力ファイル作成: epsilon.in")
print("  次のコマンドで実行: epsilon.x < epsilon.in > epsilon.out")

# =====================================
# ステップ4：結果の可視化（例示）
# =====================================
# 実際の計算結果からデータを読み込む想定

# 仮想的な誘電関数データ（実際はepsilonr.datから読み込み）
omega = np.linspace(0, 10, 200)  # eV
epsilon_real = 5.0 + 2.0 * np.exp(-((omega - 3.4)/0.5)**2)  # ε₁(ω)
epsilon_imag = 3.0 * np.exp(-((omega - 3.4)/0.3)**2)  # ε₂(ω)

# 吸収係数 α(ω) = (2ω/c) * sqrt((sqrt(ε₁² + ε₂²) - ε₁)/2)
c = 2.998e8  # m/s
alpha = (2 * omega / c) * np.sqrt((np.sqrt(epsilon_real**2 + epsilon_imag**2) - epsilon_real) / 2)

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(12, 10))

# ε₁(ω)
ax1.plot(omega, epsilon_real, 'b-', linewidth=2)
ax1.axvline(x=3.4, color='r', linestyle='--', label='Direct gap (exp: 3.4 eV)')
ax1.set_xlabel('Energy (eV)', fontsize=12)
ax1.set_ylabel('Re[ε(ω)]', fontsize=12)
ax1.set_title('Real Part of Dielectric Function', fontsize=14)
ax1.grid(True, alpha=0.3)
ax1.legend()

# ε₂(ω)
ax2.plot(omega, epsilon_imag, 'r-', linewidth=2)
ax2.axvline(x=3.4, color='r', linestyle='--', label='Direct gap (exp: 3.4 eV)')
ax2.set_xlabel('Energy (eV)', fontsize=12)
ax2.set_ylabel('Im[ε(ω)]', fontsize=12)
ax2.set_title('Imaginary Part of Dielectric Function', fontsize=14)
ax2.grid(True, alpha=0.3)
ax2.legend()

# 吸収スペクトル
ax3.plot(omega, alpha, 'g-', linewidth=2)
ax3.set_xlabel('Energy (eV)', fontsize=12)
ax3.set_ylabel('Absorption Coefficient (cm⁻¹)', fontsize=12)
ax3.set_title('Optical Absorption Spectrum', fontsize=14)
ax3.set_yscale('log')
ax3.grid(True, alpha=0.3)

# 反射率 R(ω)
n = np.sqrt(epsilon_real + 1j * epsilon_imag)  # 複素屈折率
reflectance = np.abs((n - 1) / (n + 1))**2
ax4.plot(omega, reflectance, 'm-', linewidth=2)
ax4.set_xlabel('Energy (eV)', fontsize=12)
ax4.set_ylabel('Reflectance', fontsize=12)
ax4.set_title('Normal Incidence Reflectance', fontsize=14)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('gan_optical_properties.png', dpi=300)
print("\n✅ 光学特性グラフ保存: gan_optical_properties.png")

print("\n📊 GaN光学特性（代表値）:")
print(f"  直接遷移バンドギャップ: 3.4 eV (実験値)")
print(f"  静誘電率 ε(0): {epsilon_real[0]:.2f}")
print(f"  可視光領域屈折率: ~2.3")
print(f"  吸収端: 約365 nm（紫外域）")

from ase.build import bulk
from ase.calculators.espresso import Espresso
from ase.io import write
import numpy as np
import matplotlib.pyplot as plt

# Fe bcc構造
fe = bulk('Fe', 'bcc', a=2.87)

print("Fe結晶情報:")
print(f"  構造: bcc (body-centered cubic)")
print(f"  格子定数: {fe.cell.cellpar()[0]:.3f} Å")
print(f"  空間群: Im-3m (229)")
print(f"  磁性: 強磁性体（室温）")

pseudopotentials = {'Fe': 'Fe.pbe-spn-rrkjus_psl.1.0.0.UPF'}

# =====================================
# ステップ1：スピン分極SCF計算
# =====================================
print("\nステップ1：スピン分極SCF計算...")

calc_scf = Espresso(
    pseudopotentials=pseudopotentials,
    input_data={
        'control': {
            'calculation': 'scf',
            'prefix': 'fe',
        },
        'system': {
            'ecutwfc': 60.0,
            'ecutrho': 480.0,  # 金属は高いカットオフが必要
            'occupations': 'smearing',
            'smearing': 'mv',  # Marzari-Vanderbilt (金属推奨)
            'degauss': 0.02,
            'nspin': 2,  # スピン分極計算を有効化
            'starting_magnetization(1)': 0.5,  # 初期磁化（Fe）
        },
        'electrons': {
            'conv_thr': 1.0e-8,
            'mixing_beta': 0.3,  # 磁性体は収束が遅い
        }
    },
    kpts=(16, 16, 16),  # 金属は密なk点必要
)

fe.calc = calc_scf
energy_scf = fe.get_potential_energy()
print(f"✅ SCF完了: E = {energy_scf:.6f} eV")

# =====================================
# ステップ2：磁気モーメント解析
# =====================================
# 実際の計算では、Quantum ESPRESSOの出力から抽出
print("\n📊 磁気特性解析結果:")

# 仮想的な結果（実際はSCF出力から読み取り）
total_magnetization = 2.15  # μB/atom
absolute_magnetization = 2.20  # μB/atom

print(f"  全磁気モーメント: {total_magnetization:.3f} μB/atom")
print(f"  絶対磁気モーメント: {absolute_magnetization:.3f} μB/atom")
print(f"  実験値: 2.22 μB/atom")
print(f"  誤差: {abs(total_magnetization - 2.22)/2.22*100:.1f}%")

# =====================================
# ステップ3：スピン状態密度計算
# =====================================
print("\nステップ3：スピン状態密度（DOS）計算...")

calc_dos = Espresso(
    pseudopotentials=pseudopotentials,
    input_data={
        'control': {
            'calculation': 'nscf',
            'restart_mode': 'restart',
            'prefix': 'fe',
        },
        'system': {
            'ecutwfc': 60.0,
            'ecutrho': 480.0,
            'occupations': 'tetrahedra',  # DOS計算推奨
            'nspin': 2,
        },
    },
    kpts=(24, 24, 24),  # NSCFではより密なk点
)

fe.calc = calc_dos
print("✅ NSCF完了")

# projwfc.x または dos.x で状態密度計算（実際の手順）
print("  次のステップ: projwfc.x でスピン分極DOSを計算")

# =====================================
# ステップ4：スピン分極DOS可視化
# =====================================
# 仮想的なDOSデータ（実際は計算結果から読み込み）
energy = np.linspace(-10, 10, 500)
dos_up = 1.5 * np.exp(-((energy - 2)/1.5)**2) + 0.5
dos_down = -1.0 * np.exp(-((energy - 2)/1.5)**2) - 0.3
fermi_level = 0.0

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# スピン分極DOS
ax1.fill_between(energy, dos_up, 0, alpha=0.5, color='b', label='Spin up')
ax1.fill_between(energy, dos_down, 0, alpha=0.5, color='r', label='Spin down')
ax1.axvline(x=fermi_level, color='k', linestyle='--', linewidth=2, label='Fermi level')
ax1.axhline(y=0, color='k', linewidth=0.5)
ax1.set_xlabel('Energy (eV)', fontsize=12)
ax1.set_ylabel('DOS (states/eV/atom)', fontsize=12)
ax1.set_title('Fe Spin-Polarized DOS', fontsize=14)
ax1.legend()
ax1.grid(True, alpha=0.3)

# 磁気モーメント密度
mag_density = dos_up + dos_down
ax2.plot(energy, mag_density, 'g-', linewidth=2)
ax2.axvline(x=fermi_level, color='k', linestyle='--', linewidth=2, label='Fermi level')
ax2.axhline(y=0, color='k', linewidth=0.5)
ax2.set_xlabel('Energy (eV)', fontsize=12)
ax2.set_ylabel('Magnetization Density (μB/eV/atom)', fontsize=12)
ax2.set_title('Fe Magnetic Moment Density', fontsize=14)
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('fe_magnetic_properties.png', dpi=300)
print("\n✅ 磁気特性グラフ保存: fe_magnetic_properties.png")

# =====================================
# ステップ5：交換スピン分裂解析
# =====================================
print("\n📊 交換スピン分裂解析:")
exchange_splitting = 2.1  # eV（d軌道）
print(f"  d軌道交換分裂: {exchange_splitting:.2f} eV")
print(f"  スピン分極率: {(dos_up[250] - abs(dos_down[250]))/(dos_up[250] + abs(dos_down[250]))*100:.1f}%")

from ase import Atoms
from ase.io import write
import numpy as np

# BaTiO₃ ペロブスカイト構造（立方晶相, 常誘電性）
a = 4.00  # Å（立方晶格子定数）

# 立方晶相（高温、常誘電性）
bto_cubic = Atoms(
    symbols='BaTiO3',
    scaled_positions=[
        [0.0, 0.0, 0.0],  # Ba at corner
        [0.5, 0.5, 0.5],  # Ti at body center
        [0.5, 0.5, 0.0],  # O1 at face center
        [0.5, 0.0, 0.5],  # O2
        [0.0, 0.5, 0.5],  # O3
    ],
    cell=[a, a, a],
    pbc=True
)

print("BaTiO₃ 立方晶相（常誘電性、> 120°C）:")
print(f"  格子定数: a = {a:.3f} Å")
print(f"  空間群: Pm-3m (221)")
print(f"  中心対称性: あり")

# 正方晶相（室温、強誘電性）
a_tetra = 3.99  # Å
c_tetra = 4.03  # Å（c軸方向に伸びる）
delta_z_Ti = 0.015  # Tiの変位
delta_z_O = 0.020  # 軸方向Oの変位

bto_tetra = Atoms(
    symbols='BaTiO3',
    scaled_positions=[
        [0.0, 0.0, 0.0],           # Ba
        [0.5, 0.5, 0.5 + delta_z_Ti],  # Ti（下方変位）
        [0.5, 0.5, 0.0 - delta_z_O],   # O1（上方変位）
        [0.5, 0.0, 0.5],           # O2
        [0.0, 0.5, 0.5],           # O3
    ],
    cell=[a_tetra, a_tetra, c_tetra],
    pbc=True
)

print(f"\nBaTiO₃ 正方晶相（強誘電性、室温）:")
print(f"  格子定数: a = {a_tetra:.3f} Å, c = {c_tetra:.3f} Å")
print(f"  c/a比: {c_tetra/a_tetra:.4f}")
print(f"  空間群: P4mm (99)")
print(f"  中心対称性: なし（極性構造）")
print(f"  Ti変位: {delta_z_Ti * c_tetra:.4f} Å")

# CIF保存
write('bto_cubic.cif', bto_cubic)
write('bto_tetragonal.cif', bto_tetra)
print("\n✅ BaTiO₃構造保存完了")

# =====================================
# 自発分極計算の準備
# =====================================
print("\n📋 自発分極計算手順:")
print("  1. 立方晶相の構造最適化（参照状態）")
print("  2. 正方晶相の構造最適化（分極状態）")
print("  3. Berry位相法で分極計算:")
print("     P = (P_tetra - P_cubic) / V")
print("  4. 実験値との比較:")
print("     - 実験値: 26 μC/cm² (室温)")
print("     - 計算値: 通常20-30 μC/cm²（DFT精度）")

from ase.calculators.espresso import Espresso
from ase.phonons import Phonons
import numpy as np
import matplotlib.pyplot as plt

# BaTiO₃正方晶相（前のコードで作成）
# bto_tetra = [BaTiO3 Atoms object]

pseudopotentials = {
    'Ba': 'Ba.pbe-spn-rrkjus_psl.1.0.0.UPF',
    'Ti': 'Ti.pbe-spn-rrkjus_psl.1.0.0.UPF',
    'O': 'O.pbe-n-rrkjus_psl.1.0.0.UPF'
}

# =====================================
# ステップ1：構造最適化
# =====================================
print("ステップ1：構造最適化...")

calc_relax = Espresso(
    pseudopotentials=pseudopotentials,
    input_data={
        'control': {
            'calculation': 'vc-relax',  # 格子定数も最適化
            'prefix': 'bto',
        },
        'system': {
            'ecutwfc': 60.0,
            'ecutrho': 600.0,  # 酸化物は高いカットオフ必要
            'occupations': 'fixed',
        },
        'electrons': {
            'conv_thr': 1.0e-9,
        },
        'ions': {
            'ion_dynamics': 'bfgs',
        },
        'cell': {
            'cell_dynamics': 'bfgs',
        }
    },
    kpts=(8, 8, 8),
)

bto_tetra.calc = calc_relax
energy_relax = bto_tetra.get_potential_energy()
optimized_cell = bto_tetra.get_cell()

print(f"✅ 構造最適化完了")
print(f"  最適化後格子定数: a = {optimized_cell[0][0]:.4f} Å, c = {optimized_cell[2][2]:.4f} Å")

# =====================================
# ステップ2：誘電定数計算（DFPT）
# =====================================
print("\nステップ2：誘電定数計算（DFPT: 密度汎関数摂動論）...")

# ph.xの入力ファイル作成
ph_input = """
 &inputph
   prefix = 'bto'
   tr2_ph = 1.0d-14
   epsil = .true.
   zeu = .true.
   ldisp = .false.
   qplot = .false.
 /
 0.0 0.0 0.0
"""

with open('ph_epsilon.in', 'w') as f:
    f.write(ph_input)

print("✅ ph.x入力ファイル作成: ph_epsilon.in")
print("  次のコマンドで実行: ph.x < ph_epsilon.in > ph_epsilon.out")
print("  出力から誘電テンソルとBorn有効電荷を取得")

# =====================================
# ステップ3：結果解析（例示）
# =====================================
# 実際の計算結果から読み取る想定

print("\n📊 BaTiO₃誘電特性（計算結果例）:")

# 電子誘電率（高周波）
epsilon_inf = np.array([
    [6.5, 0.0, 0.0],
    [0.0, 6.5, 0.0],
    [0.0, 0.0, 6.8]
])

print("\n1. 電子誘電率 ε∞（光学誘電率）:")
print(f"  ε∞_⊥（a軸方向）: {epsilon_inf[0,0]:.2f}")
print(f"  ε∞_∥（c軸方向）: {epsilon_inf[2,2]:.2f}")

# イオン誘電率寄与（格子振動）
epsilon_ionic = np.array([
    [150, 0, 0],
    [0, 150, 0],
    [0, 0, 200]
])

print("\n2. イオン誘電率寄与 Δε（格子振動）:")
print(f"  Δε_⊥: {epsilon_ionic[0,0]:.0f}")
print(f"  Δε_∥: {epsilon_ionic[2,2]:.0f}")

# 静誘電率
epsilon_static = epsilon_inf + epsilon_ionic
print("\n3. 静誘電率 ε₀:")
print(f"  ε₀_⊥: {epsilon_static[0,0]:.0f}")
print(f"  ε₀_∥: {epsilon_static[2,2]:.0f}")
print(f"  実験値（室温）: ε₀_⊥ ≈ 170, ε₀_∥ ≈ 90")

# Born有効電荷
Z_Ti_star = 7.2  # e（Ti原子）
Z_O_star = -5.8  # e（軸方向O）

print(f"\n4. Born有効電荷:")
print(f"  Z*(Ti): +{Z_Ti_star:.2f} e（名目価数: +4）")
print(f"  Z*(O, axial): {Z_O_star:.2f} e（名目価数: -2）")
print(f"  → 強い共有結合性を反映")

# =====================================
# ステップ4：ソフトモード解析
# =====================================
print("\n5. ソフトモード（強誘電不安定性）:")
omega_soft = 50  # cm⁻¹（Γ点のA1モード）
print(f"  A1モード振動数: {omega_soft:.0f} cm⁻¹")
print(f"  → 立方晶では虚数振動数（構造不安定）")
print(f"  → 正方晶では正の振動数（安定）")

# =====================================
# ステップ5：誘電率の温度依存性（Curie-Weiss則）
# =====================================
T = np.linspace(130, 400, 100)  # K
T_c = 393  # K（キュリー温度）
C = 1.7e5  # K（Curie定数）

epsilon_para = epsilon_inf[0,0] + C / (T - T_c)  # 常誘電相（T > Tc）

fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(T, epsilon_para, 'r-', linewidth=2, label='Paraelectric phase (T > Tc)')
ax.axvline(x=T_c, color='k', linestyle='--', linewidth=2, label=f'Tc = {T_c} K')
ax.set_xlabel('Temperature (K)', fontsize=12)
ax.set_ylabel('Dielectric Constant ε', fontsize=12)
ax.set_title('BaTiO₃ Temperature Dependence (Curie-Weiss Law)', fontsize=14)
ax.set_ylim(0, 10000)
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('bto_curie_weiss.png', dpi=300)
print("\n✅ 誘電率温度依存性グラフ保存: bto_curie_weiss.png")

"""
convergence_tester.py
任意の材料に対してk点・エネルギーカットオフ収束テストを実行
"""

import numpy as np
import matplotlib.pyplot as plt
from ase.calculators.espresso import Espresso
from typing import List, Tuple, Dict
import json

class ConvergenceTester:
    """DFT計算パラメータの収束テストを実行するクラス"""

    def __init__(self, atoms, pseudopotentials, threshold_meV=1.0):
        """
        Parameters
        ----------
        atoms : ase.Atoms
            テスト対象の原子構造
        pseudopotentials : dict
            擬ポテンシャル辞書 {'元素': 'ファイル名'}
        threshold_meV : float
            収束判定閾値（meV/atom）
        """
        self.atoms = atoms
        self.pseudopotentials = pseudopotentials
        self.threshold = threshold_meV / 1000  # eV/atomに変換
        self.results = {}

    def test_kpoints(self,
                     kpts_list: List[int],
                     ecutwfc: float,
                     base_input: Dict = None) -> Tuple[List[float], Dict]:
        """
        k点メッシュの収束テスト

        Parameters
        ----------
        kpts_list : List[int]
            テストするk点リスト（例: [4, 6, 8, 10, 12]）
        ecutwfc : float
            固定するエネルギーカットオフ（Ry）
        base_input : Dict
            追加の計算設定

        Returns
        -------
        energies : List[float]
            各k点での全エネルギー
        convergence_info : Dict
            収束情報
        """
        print(f"\n{'='*60}")
        print(f"k点収束テスト開始（ecutwfc = {ecutwfc} Ry 固定）")
        print(f"{'='*60}")

        if base_input is None:
            base_input = {
                'control': {'calculation': 'scf'},
                'system': {'ecutwfc': ecutwfc},
                'electrons': {'conv_thr': 1.0e-8},
            }
        else:
            base_input['system']['ecutwfc'] = ecutwfc

        energies = []
        for k in kpts_list:
            calc = Espresso(
                pseudopotentials=self.pseudopotentials,
                input_data=base_input,
                kpts=(k, k, k),
            )
            self.atoms.calc = calc
            energy = self.atoms.get_potential_energy()
            energies.append(energy)
            print(f"  k = {k:2d}×{k:2d}×{k:2d}: E = {energy:12.6f} eV")

        # 収束判定
        converged_idx = self._check_convergence(energies)
        optimal_k = kpts_list[converged_idx] if converged_idx is not None else kpts_list[-1]

        convergence_info = {
            'converged': converged_idx is not None,
            'optimal_k': optimal_k,
            'threshold_meV': self.threshold * 1000,
            'final_energy': energies[-1],
        }

        print(f"\n📊 収束判定:")
        if converged_idx is not None:
            print(f"  ✅ 収束達成: k = {optimal_k}")
        else:
            print(f"  ⚠️ 未収束: より密なk点が必要")

        self.results['kpoints'] = {
            'kpts_list': kpts_list,
            'energies': energies,
            'convergence_info': convergence_info
        }

        return energies, convergence_info

    def test_ecutwfc(self,
                     ecutwfc_list: List[float],
                     kpts: Tuple[int, int, int],
                     base_input: Dict = None) -> Tuple[List[float], Dict]:
        """
        エネルギーカットオフの収束テスト

        Parameters
        ----------
        ecutwfc_list : List[float]
            テストするカットオフリスト（Ry）
        kpts : Tuple[int, int, int]
            固定するk点メッシュ
        base_input : Dict
            追加の計算設定

        Returns
        -------
        energies : List[float]
            各カットオフでの全エネルギー
        convergence_info : Dict
            収束情報
        """
        print(f"\n{'='*60}")
        print(f"エネルギーカットオフ収束テスト開始（kpts = {kpts} 固定）")
        print(f"{'='*60}")

        if base_input is None:
            base_input = {
                'control': {'calculation': 'scf'},
                'system': {},
                'electrons': {'conv_thr': 1.0e-8},
            }

        energies = []
        for ecut in ecutwfc_list:
            base_input['system']['ecutwfc'] = ecut
            calc = Espresso(
                pseudopotentials=self.pseudopotentials,
                input_data=base_input,
                kpts=kpts,
            )
            self.atoms.calc = calc
            energy = self.atoms.get_potential_energy()
            energies.append(energy)
            print(f"  ecutwfc = {ecut:5.1f} Ry: E = {energy:12.6f} eV")

        # 収束判定
        converged_idx = self._check_convergence(energies)
        optimal_ecut = ecutwfc_list[converged_idx] if converged_idx is not None else ecutwfc_list[-1]

        convergence_info = {
            'converged': converged_idx is not None,
            'optimal_ecutwfc': optimal_ecut,
            'threshold_meV': self.threshold * 1000,
            'final_energy': energies[-1],
        }

        print(f"\n📊 収束判定:")
        if converged_idx is not None:
            print(f"  ✅ 収束達成: ecutwfc = {optimal_ecut} Ry")
        else:
            print(f"  ⚠️ 未収束: より高いカットオフが必要")

        self.results['ecutwfc'] = {
            'ecutwfc_list': ecutwfc_list,
            'energies': energies,
            'convergence_info': convergence_info
        }

        return energies, convergence_info

    def _check_convergence(self, energies: List[float]) -> int:
        """
        エネルギー変化から収束を判定

        Returns
        -------
        converged_idx : int or None
            収束したインデックス（Noneは未収束）
        """
        n_atoms = len(self.atoms)
        for i in range(1, len(energies)):
            delta = abs(energies[i] - energies[i-1]) / n_atoms
            if delta < self.threshold:
                return i
        return None

    def plot_results(self, filename='convergence_tests.png'):
        """収束テスト結果のプロット"""
        fig, axes = plt.subplots(1, 2, figsize=(14, 5))

        # k点収束
        if 'kpoints' in self.results:
            data = self.results['kpoints']
            ax = axes[0]
            ax.plot(data['kpts_list'], data['energies'], 'o-',
                   linewidth=2, markersize=8, label='Calculated')
            if data['convergence_info']['converged']:
                opt_k = data['convergence_info']['optimal_k']
                ax.axvline(x=opt_k, color='r', linestyle='--',
                          label=f'Converged: {opt_k}')
            ax.set_xlabel('k-points (k×k×k)', fontsize=12)
            ax.set_ylabel('Total Energy (eV)', fontsize=12)
            ax.set_title('k-point Convergence', fontsize=14)
            ax.grid(True, alpha=0.3)
            ax.legend()

        # エネルギーカットオフ収束
        if 'ecutwfc' in self.results:
            data = self.results['ecutwfc']
            ax = axes[1]
            ax.plot(data['ecutwfc_list'], data['energies'], 's-',
                   linewidth=2, markersize=8, label='Calculated')
            if data['convergence_info']['converged']:
                opt_ecut = data['convergence_info']['optimal_ecutwfc']
                ax.axvline(x=opt_ecut, color='r', linestyle='--',
                          label=f'Converged: {opt_ecut} Ry')
            ax.set_xlabel('Energy Cutoff (Ry)', fontsize=12)
            ax.set_ylabel('Total Energy (eV)', fontsize=12)
            ax.set_title('Energy Cutoff Convergence', fontsize=14)
            ax.grid(True, alpha=0.3)
            ax.legend()

        plt.tight_layout()
        plt.savefig(filename, dpi=300)
        print(f"\n✅ 収束テストグラフ保存: {filename}")

    def save_results(self, filename='convergence_results.json'):
        """結果をJSONファイルに保存"""
        with open(filename, 'w') as f:
            json.dump(self.results, f, indent=2)
        print(f"✅ 収束テスト結果保存: {filename}")


# =====================================
# 使用例
# =====================================
if __name__ == "__main__":
    from ase.build import bulk

    # Si結晶
    si = bulk('Si', 'diamond', a=5.43)
    pseudopotentials = {'Si': 'Si.pbe-n-rrkjus_psl.1.0.0.UPF'}

    # テスター初期化
    tester = ConvergenceTester(si, pseudopotentials, threshold_meV=1.0)

    # k点テスト
    tester.test_kpoints(
        kpts_list=[4, 6, 8, 10, 12],
        ecutwfc=30.0
    )

    # エネルギーカットオフテスト
    tester.test_ecutwfc(
        ecutwfc_list=[20, 25, 30, 35, 40, 45],
        kpts=(8, 8, 8)
    )

    # 結果の可視化と保存
    tester.plot_results('si_convergence.png')
    tester.save_results('si_convergence.json')

    print("\n🎉 収束テスト完了！")

"""
materials_workflow_manager.py
複数材料の物性計算を自動化するワークフロー管理システム
"""

import os
import json
import time
from pathlib import Path
from typing import List, Dict, Any
from datetime import datetime
from ase import Atoms
from ase.io import read, write
from ase.calculators.espresso import Espresso
import yaml

class MaterialsWorkflow:
    """物性計算ワークフローの自動化システム"""

    def __init__(self, config_file: str):
        """
        Parameters
        ----------
        config_file : str
            ワークフロー設定ファイル（YAML形式）
        """
        with open(config_file, 'r') as f:
            self.config = yaml.safe_load(f)

        self.results_dir = Path(self.config['output_dir'])
        self.results_dir.mkdir(exist_ok=True)

        self.log_file = self.results_dir / 'workflow.log'
        self._log("ワークフローマネージャー初期化完了")

    def run_workflow(self):
        """全材料のワークフローを実行"""
        materials = self.config['materials']

        self._log(f"\n{'='*70}")
        self._log(f"ワークフロー実行開始: {len(materials)}材料")
        self._log(f"{'='*70}")

        results_summary = []

        for i, mat_config in enumerate(materials, 1):
            self._log(f"\n[{i}/{len(materials)}] {mat_config['name']} 処理開始...")

            try:
                result = self._process_material(mat_config)
                results_summary.append(result)
                self._log(f"✅ {mat_config['name']} 完了")
            except Exception as e:
                self._log(f"❌ {mat_config['name']} エラー: {str(e)}")
                results_summary.append({
                    'name': mat_config['name'],
                    'status': 'failed',
                    'error': str(e)
                })

        # 最終レポート作成
        self._generate_report(results_summary)
        self._log("\n🎉 全ワークフロー完了！")

    def _process_material(self, mat_config: Dict) -> Dict:
        """単一材料の完全ワークフロー実行"""
        mat_name = mat_config['name']
        mat_dir = self.results_dir / mat_name
        mat_dir.mkdir(exist_ok=True)

        result = {
            'name': mat_name,
            'status': 'running',
            'start_time': datetime.now().isoformat()
        }

        # ステップ1：構造作成
        atoms = self._create_structure(mat_config)
        write(mat_dir / f'{mat_name}.cif', atoms)
        result['structure_created'] = True

        # ステップ2：収束テスト（設定で有効な場合）
        if self.config.get('run_convergence_tests', True):
            conv_result = self._run_convergence_tests(atoms, mat_config, mat_dir)
            result['convergence'] = conv_result
            optimal_kpts = conv_result['optimal_kpts']
            optimal_ecutwfc = conv_result['optimal_ecutwfc']
        else:
            optimal_kpts = tuple(mat_config.get('kpts', [8, 8, 8]))
            optimal_ecutwfc = mat_config.get('ecutwfc', 40.0)

        # ステップ3：物性計算
        properties = mat_config.get('properties', ['scf', 'bands', 'dos'])
        result['properties'] = {}

        for prop in properties:
            self._log(f"  - {prop}計算実行中...")
            prop_result = self._calculate_property(
                atoms, prop, mat_config,
                optimal_kpts, optimal_ecutwfc, mat_dir
            )
            result['properties'][prop] = prop_result

        # ステップ4：後処理と可視化
        self._postprocess(atoms, result, mat_dir)

        result['status'] = 'completed'
        result['end_time'] = datetime.now().isoformat()

        # 材料ごとの結果保存
        with open(mat_dir / 'results.json', 'w') as f:
            json.dump(result, f, indent=2)

        return result

    def _create_structure(self, mat_config: Dict) -> Atoms:
        """結晶構造作成"""
        if 'structure_file' in mat_config:
            return read(mat_config['structure_file'])

        # 組み込み構造生成
        from ase.build import bulk
        name = mat_config['formula']
        crystal_type = mat_config.get('crystal_type', 'fcc')
        lattice_constant = mat_config.get('lattice_constant', 4.0)

        return bulk(name, crystal_type, a=lattice_constant)

    def _run_convergence_tests(self, atoms: Atoms,
                                mat_config: Dict,
                                output_dir: Path) -> Dict:
        """収束テスト実行"""
        from convergence_tester import ConvergenceTester

        pseudopotentials = mat_config['pseudopotentials']
        tester = ConvergenceTester(atoms, pseudopotentials, threshold_meV=1.0)

        # k点テスト
        kpts_range = mat_config.get('kpts_test_range', [4, 6, 8, 10, 12])
        ecutwfc_test = mat_config.get('ecutwfc_test', 30.0)
        tester.test_kpoints(kpts_range, ecutwfc_test)

        # エネルギーカットオフテスト
        ecutwfc_range = mat_config.get('ecutwfc_test_range', [20, 30, 40, 50])
        kpts_test = tuple(mat_config.get('kpts_test', [8, 8, 8]))
        tester.test_ecutwfc(ecutwfc_range, kpts_test)

        # 結果保存
        tester.plot_results(str(output_dir / 'convergence.png'))
        tester.save_results(str(output_dir / 'convergence.json'))

        return {
            'optimal_kpts': tuple([tester.results['kpoints']['convergence_info']['optimal_k']] * 3),
            'optimal_ecutwfc': tester.results['ecutwfc']['convergence_info']['optimal_ecutwfc']
        }

    def _calculate_property(self, atoms: Atoms,
                           property_type: str,
                           mat_config: Dict,
                           kpts: tuple,
                           ecutwfc: float,
                           output_dir: Path) -> Dict:
        """指定された物性を計算"""
        pseudopotentials = mat_config['pseudopotentials']

        # 物性タイプに応じた計算設定
        if property_type == 'scf':
            return self._run_scf(atoms, pseudopotentials, kpts, ecutwfc, output_dir)

        elif property_type == 'bands':
            return self._run_bands(atoms, pseudopotentials, kpts, ecutwfc, output_dir)

        elif property_type == 'dos':
            return self._run_dos(atoms, pseudopotentials, kpts, ecutwfc, output_dir)

        elif property_type == 'optical':
            return self._run_optical(atoms, pseudopotentials, kpts, ecutwfc, output_dir)

        elif property_type == 'magnetic':
            return self._run_magnetic(atoms, pseudopotentials, kpts, ecutwfc, output_dir)

        else:
            raise ValueError(f"未対応の物性タイプ: {property_type}")

    def _run_scf(self, atoms, pseudopotentials, kpts, ecutwfc, output_dir):
        """SCF計算"""
        calc = Espresso(
            pseudopotentials=pseudopotentials,
            input_data={
                'control': {'calculation': 'scf', 'prefix': 'scf'},
                'system': {'ecutwfc': ecutwfc},
                'electrons': {'conv_thr': 1.0e-8},
            },
            kpts=kpts,
        )
        atoms.calc = calc
        energy = atoms.get_potential_energy()
        forces = atoms.get_forces()

        return {
            'energy': energy,
            'energy_per_atom': energy / len(atoms),
            'max_force': float(np.max(np.abs(forces))),
        }

    def _run_bands(self, atoms, pseudopotentials, kpts, ecutwfc, output_dir):
        """バンド構造計算"""
        # 実装は簡略化（実際はコード例3に準拠）
        return {'status': 'completed', 'output': 'bands.dat'}

    def _run_dos(self, atoms, pseudopotentials, kpts, ecutwfc, output_dir):
        """状態密度計算"""
        return {'status': 'completed', 'output': 'dos.dat'}

    def _run_optical(self, atoms, pseudopotentials, kpts, ecutwfc, output_dir):
        """光学特性計算"""
        return {'status': 'completed', 'output': 'epsilon.dat'}

    def _run_magnetic(self, atoms, pseudopotentials, kpts, ecutwfc, output_dir):
        """磁気特性計算"""
        return {'status': 'completed', 'magnetization': 2.15}

    def _postprocess(self, atoms: Atoms, result: Dict, output_dir: Path):
        """結果の後処理と可視化"""
        # プロット生成などの後処理
        pass

    def _generate_report(self, results_summary: List[Dict]):
        """最終レポート生成"""
        report_file = self.results_dir / 'workflow_report.md'

        with open(report_file, 'w') as f:
            f.write("# 物性計算ワークフロー実行レポート\n\n")
            f.write(f"実行日時: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}\n\n")

            f.write("## 実行サマリー\n\n")
            total = len(results_summary)
            completed = sum(1 for r in results_summary if r['status'] == 'completed')
            failed = total - completed

            f.write(f"- 総材料数: {total}\n")
            f.write(f"- 成功: {completed}\n")
            f.write(f"- 失敗: {failed}\n\n")

            f.write("## 材料別詳細\n\n")
            for result in results_summary:
                f.write(f"### {result['name']}\n\n")
                f.write(f"- ステータス: {result['status']}\n")
                if 'properties' in result:
                    f.write(f"- 計算物性: {', '.join(result['properties'].keys())}\n")
                if 'error' in result:
                    f.write(f"- エラー: {result['error']}\n")
                f.write("\n")

        self._log(f"\n✅ レポート生成完了: {report_file}")

    def _log(self, message: str):
        """ログ出力"""
        timestamp = datetime.now().strftime('%Y-%m-%d %H:%M:%S')
        log_message = f"[{timestamp}] {message}"
        print(log_message)

        with open(self.log_file, 'a') as f:
            f.write(log_message + '\n')


# =====================================
# 設定ファイル例（workflow_config.yaml）
# =====================================
example_config = """
output_dir: ./workflow_results
run_convergence_tests: true

materials:
  - name: Si
    formula: Si
    crystal_type: diamond
    lattice_constant: 5.43
    pseudopotentials:
      Si: Si.pbe-n-rrkjus_psl.1.0.0.UPF
    properties: [scf, bands, dos]
    kpts_test_range: [4, 6, 8, 10]
    ecutwfc_test_range: [25, 30, 35, 40]

  - name: GaN
    formula: GaN
    structure_file: gan_wurtzite.cif
    pseudopotentials:
      Ga: Ga.pbe-dnl-rrkjus_psl.1.0.0.UPF
      N: N.pbe-n-radius_5.UPF
    properties: [scf, bands, optical]
    kpts: [12, 12, 8]
    ecutwfc: 50.0

  - name: Fe
    formula: Fe
    crystal_type: bcc
    lattice_constant: 2.87
    pseudopotentials:
      Fe: Fe.pbe-spn-rrkjus_psl.1.0.0.UPF
    properties: [scf, magnetic]
    kpts: [16, 16, 16]
    ecutwfc: 60.0
"""

# 使用例
if __name__ == "__main__":
    # 設定ファイル作成
    with open('workflow_config.yaml', 'w') as f:
        f.write(example_config)

    # ワークフロー実行
    workflow = MaterialsWorkflow('workflow_config.yaml')
    workflow.run_workflow()