この記事では、フォノン計算用のパッケージ「Phonopy」を用いたフォノン計算の方法を解説します。Phonopyではフォノンのバンド構造、状態密度の計算や、その結果を用いて熱力学的物性などの物性の計算が可能です。
なお、Phonopyを使ってフォノンバンドの描画する際、バンドを描画する波数空間内の経路を、与えられたセル形状から自動的に求めるに当たってはsumoを用いるのが便利です。ここではsumoも合わせて用いていきます。
どちらのパッケージもpip installコマンドから簡単にインストールできます。
!pip install phonopy
!pip install sumo
ここでは、ダイヤモンド構造のSiを例に解説していきます。
フォノンを計算するに当たっては、まずは構造最適化が必要になります。フォノンの計算は、エネルギー的に極小値となる構造に微小変位を与えたときの力の変化をもとに計算していきますので、最適化の収束条件は厳し目に取ったほうが良いです。
また、セルの形状にはprimitive unit cellを入力構造とした方が系の対称性を正しく反映した結果を得ることができます。
ここではget_primitive_atoms関数とoptimize_atoms関数を定義して、まず前者を使って入力構造をprimitive unit cellに変換した後、後者を使って構造最適化を行っていきます（収束条件はfmax=0.0001)。
from ase.build import bulk
from ase.optimize import FIRE, LBFGS
from ase.filters import FrechetCellFilter
from ase.constraints import FixSymmetry
import pfp_api_client
from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
from pfp_api_client.pfp.estimator import Estimator, EstimatorCalcMode
from pymatgen.io.ase import AseAtomsAdaptor
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
def get_primitive_atoms(atoms_init):
atoms = atoms_init.copy()
adaptor = AseAtomsAdaptor()
struct = adaptor.get_structure(atoms)
sga = SpacegroupAnalyzer(struct, symprec=0.01)
struct_prim = sga.get_primitive_standard_structure()
atoms_prim = adaptor.get_atoms(struct_prim)
return atoms
def optimize_atoms(atoms_init, calculator, opt_algo=LBFGS, fmax=0.0001, steps=10000):
atoms = atoms_init.copy()
atoms.calc = calculator
atoms.set_constraint(FixSymmetry(atoms))
cf = FrechetCellFilter(atoms) # セル可変最適化可能
Epot = atoms.get_potential_energy()
print(f"E_pot before opt: {Epot} eV\n")
opt = opt_algo(cf, logfile=None)
opt.run(fmax=fmax, steps=steps)
atoms.set_constraint()
cell_diff = (atoms.cell.cellpar() / atoms_init.cell.cellpar() - 1.0) * 100
print("Optimized Cell:\n", atoms.cell.cellpar(), "\n")
print("Optimized Cell diff (%):\n", cell_diff, "\n")
print("Scaled positions:\n", atoms.get_scaled_positions(), "\n")
Epot = atoms.get_potential_energy()
print(f"E_pot after opt: {Epot} eV")
return atoms
estimator = Estimator(calc_mode=EstimatorCalcMode.R2SCAN, model_version="v8.0.0")
calculator = ASECalculator(estimator)
atoms = bulk("Si")
atoms.calc = calculator
atoms_opt = optimize_atoms(atoms, calculator)
Phonopyを用いてフォノン計算を行っていく際は、PhonopyのPhonopyAtomsオブジェクトを使っていきます。ASEのAtomsオブジェクトとの相互変換が必要になってくるので、ここではそれぞれの変換の関数を定義します。
from ase import Atoms
from phonopy import Phonopy
from phonopy.structure.atoms import PhonopyAtoms
def phonopy_to_ase(phonopy_atoms: PhonopyAtoms):
atoms_ase = Atoms(symbols=phonopy_atoms.symbols,
positions=phonopy_atoms.positions,
cell=phonopy_atoms.cell,
pbc=True)
return atoms_ase
def ase_to_phonopy(atoms: Atoms):
atoms_phonopy = PhonopyAtoms(symbols=atoms.get_chemical_symbols(),
cell=atoms.cell,
positions=atoms.get_positions())
return atoms_phonopy
フォノン計算を行っていく際は、原子を平衡位置から様々なパターンで微小変位させた構造を作り、微小変位前後の力の差から変位に対する力の数値微分を計算することでフォノンの振動数を計算していきます。
微小変位は基本単位構造の周期性を満たすとは限らないので、スーパーセルを用いて計算していきます。大きなスーパーセル構造を取るほど長波長振動に対する結果が精密になっていきますが、計算量も膨大になっていきますので、ほどほどのところで打ち切りにする必要があります。ここではprimitive unit cellに対して4x4x4サイズのスーパーセルもとに、そのスーパーセルサイズのもとで可能な微小変位パターンを持ったスーパーセルを生成していきます。微小変位の値はphonon.generate_displacements()メソッドのdistanceという引数で指定します。
# スーパーセルサイズの指定
supercell_matrix = [4,4,4]
# ASEのAtomsオブジェクトをPhonopyのPhonopyAtomsオブジェクトにに変換
unitcell = ase_to_phonopy(atoms_opt)
# スーパーセルの作成
phonon = Phonopy(unitcell, supercell_matrix)
# 原子の微小変位を生成
phonon.generate_displacements(distance=0.1)
supercells = phonon.supercells_with_displacements
生成した微小変位構造に対して力を計算し、微小変位のない平衡位置における力との差異から数値微分を用いて力定数行列 (interatomic force constant matrix) を計算します。
# 力の計算
forces = []
for sc in supercells:
sc_ase = phonopy_to_ase(sc)
sc_ase.calc = calculator
forces.append(sc_ase.get_forces())
# 力定数の計算
phonon.forces = forces
phonon.produce_force_constants()
ひとたび力定数行列が得られれば、あとはそれをフーリエ変換して波数空間上の任意の点における固有振動数を求めることができます。これを波数空間上の特定の経路上にプロットしたものがフォノンバンド図となります。
一般にバンド構造は、与えられたセルの対称性を反映した対称性の高い点 (高対称点) を結んだ経路上でプロットしていきます。セルの対称性ごとに標準的に使われる経路が決まっています。
ここではsumoを用いて自動的にバンドをプロットする経路を生成していきます。
from phonopy.phonon.band_structure import get_band_qpoints_and_path_connections
from sumo.symmetry.kpoints import get_path_data
from pymatgen.io.ase import AseAtomsAdaptor
# ASEのAtomsオブジェクトをpymatgenのStructureオブジェクトに変換
structure = AseAtomsAdaptor.get_structure(atoms_opt*supercell_matrix)
kpath, band_paths, labels = get_path_data(structure, mode="pymatgen", symprec=0.0001, line_density=60, phonopy=True)
print("Generated path labels:", labels)
print(len(band_paths))
ここで得られた経路および高対称点のラベル名を用いて、バンド図を描画します。
labels = ['$'+label+'$' for label in labels]
#path_connections = [True, True, True, True, True, True, True, True, False]
# フォノンバンド構造のプロット
phonon.run_band_structure(paths=band_paths, labels=labels)
phonon.plot_band_structure().show()
一方、状態密度は特定の経路上でプロットするのではなく、波数空間上の単位構造 (第一ブリルアンゾーン) 内をメッシュで区切ったグリッド上の点でフォノン振動数を計算し、振動数あたりに存在する振動状態の数をカウントして求めてきます。ここでは20x20x20でメッシュを区切って状態密度を計算してプロットします。
q_mesh = [20,20,20]
phonon.run_mesh(q_mesh)
phonon.run_total_dos()
phonon.plot_total_dos().show()
フォノンバンドとフォノン状態密度を横に並べて描画することも可能です。状態密度のピークがどのバンドに対応するかを理解するのに役立ちます。
phonon.run_mesh(q_mesh)
phonon.run_total_dos()
phonon.plot_band_structure_and_dos().show()
また、フォノンの状態密度からは、ヘルムホルツ自由エネルギーやエントロピー、比熱などの熱力学物性の温度依存性も計算することができます。あくまでも格子振動由来の物性になるため、電子由来の寄与は入っていませんが、半導体や絶縁体などでは電子由来の寄与は無視できることが多いので格子振動由来の寄与のみでほぼ説明できます。
phonon.run_mesh(q_mesh)
phonon.run_thermal_properties(t_step=10,
t_max=1000,
t_min=0)
tp_dict = phonon.get_thermal_properties_dict()
temperatures = tp_dict['temperatures']
free_energy = tp_dict['free_energy']
entropy = tp_dict['entropy']
heat_capacity = tp_dict['heat_capacity']
for t, F, S, cv in zip(temperatures, free_energy, entropy, heat_capacity):
print(("%12.3f " + "%15.7f" * 3) % ( t, F, S, cv ))
phonon.plot_thermal_properties().show()
ただし、ここまでの計算は調和近似に基づいたものになっています。物理的には、フォノン間の相互作用などが存在しない状況を計算していることになり、フォノン間の散乱、フォノンの寿命などの効果が入っていない結果になっています。この状況は低温（具体的にはデバイ温度以下）では良い近似になっているのですが、高温になってくると現実との乖離が大きくなってきます。そのようなときには3次以上の効果を取り入れる手法が必要になってきます。