CIFを入力として第一原理計算による構造最適化を行う

はじめに

前回は、Google Colab上でQuantum Espresso(QE)を導入し、QEの機能でSiのダイヤモンド型結晶構造を生成してその全エネルギーを求める方法を解説しました。

今回は、格子定数と空間群の解析はそれなりに自信があるが原子位置は自信がないという粉末X線構造解析でよくぶつかる問題を念頭に、実例を使ってQEで最適な原子位置を計算してみます。この記事で使う結晶構造は、以前実際に粉末X線回折の実例で解析したSnO₂のCIFを用い、意図的に酸素原子の位置をおかしくしたものを対象にします。

対象となる結晶構造

まず、粉末X線構造解析の結果出てきたのは以下のCIFです。

これをVESTAという結晶構造描画アプリで開くと以下のようになります。赤い丸が酸素原子、灰色の丸がSnです。各原子の間の棒は、互いが閾値以下の距離にあることを示します。この構造でのSnとOの距離(棒の長さ)は約2 Å( = 0.2 nm)です。

これは粉末X線回折の解析だけでいい感じの酸素原子位置が得られたものですが、今回は解析がうまくいっていないモデルにするために、酸素原子の位置を(0.195, 0.195, 0)から(0.1, 0.1, 0)にします。そうすると以下のようになります。

Sn-O距離が、長いものは2.6 Å、短いものは1.7 Å程度になりました。これはルチル型SnO₂として不自然な解です。この構造にnonOptという名前を付けておきます。

O位置だけを最適化する

実際にGoogle Colabで原子位置を最適化をしていきましょう。今回はO原子位置だけを最適化します。もし金属原子位置も最適化したい場合は、コード中の以下の部分をコメントアウトしてください。

# ===== 1.5 O 以外を固定（= O だけ動かす） =====
mask_fix = [atom.symbol != "O" for atom in atoms]
atoms.set_constraint(FixAtoms(mask=mask_fix))

QEのセットアップ

ではまず、前回同様QEのセットアップから行います。すでに行っている場合は飛ばしてOKです。

# Quantum ESPRESSO と必要ライブラリの導入
!apt-get update -qq
!apt-get install -y quantum-espresso
!pip -q install --upgrade ase pymatgen

# ASE（最新版）
!pip -q install --upgrade ase

# pw.x が見えるか確認
!which pw.x
!pw.x -h | head

SnとOのPAW型擬ポテンシャルをダウンロード

今回は動作例としてQEの配布表にあるPAW型擬ポテンシャルを用います。

# Sn と O の擬ポテンシャルを取得
!wget -q https://pseudopotentials.quantum-espresso.org/upf_files/Sn.pbe-dn-kjpaw_psl.1.0.0.UPF
!wget -q https://pseudopotentials.quantum-espresso.org/upf_files/O.pbe-n-kjpaw_psl.1.0.0.UPF

CIFファイルをアップロードする

Google Colabではwgetしたファイルなどは/contentフォルダに保存され、ファイル名を指定するだけでスクリプトから読み込むことが可能です。このフォルダ内のものは一定時間後に消えるので、計算時間が長いとか、試行錯誤が必要な場合は自身のGoogle Driveから直接読むようにするのが基本です。今回は計算が短時間で済むので/contentにアップロードする方法を示します。

# Colab に CIF をアップロードして使う場合
from google.colab import files
uploaded = files.upload()

実行するとファイルのアップロードボタンが出てきますので、それを押してローカル上のCIFを/contentにアップロードしてください。

計算を実行する

以下のコードを実行し、10～20分くらい待ちます。

from pymatgen.core import Structure
from pymatgen.io.ase import AseAtomsAdaptor
from ase.calculators.espresso import Espresso, EspressoProfile
from ase.constraints import FixAtoms
from ase.io import read, write

# ===== 1. CIF の読み込み =====
cif_file = "SnO2-EDMA_nonOpt.cif"

pmg_structure = Structure.from_file(cif_file)
atoms = AseAtomsAdaptor.get_atoms(pmg_structure)

print("読み込んだ構造")
print(atoms)
print("格子定数")
print(atoms.cell)
print("分率座標")
for symbol, pos in zip(atoms.get_chemical_symbols(), atoms.get_scaled_positions()):
    print(symbol, pos)

# ===== 1.5 O 以外を固定（= O だけ動かす） =====
mask_fix = [atom.symbol != "O" for atom in atoms]
atoms.set_constraint(FixAtoms(mask=mask_fix))

# ===== 2. 擬ポテンシャルの指定 =====
pseudopotentials = {
    "Sn": "Sn.pbe-dn-kjpaw_psl.1.0.0.UPF",
    "O": "O.pbe-n-kjpaw_psl.1.0.0.UPF",
}

# ===== 3. Quantum ESPRESSO の設定 =====
profile = EspressoProfile(
    command="pw.x",
    pseudo_dir="."
)

calc = Espresso(
    profile=profile,
    pseudopotentials=pseudopotentials,
    input_data={
        "control": {
            "calculation": "relax",
            "forc_conv_thr": 5.0e-4,
            "tstress": False,
            "tprnfor": True,
        },
        "system": {
            "ecutwfc": 50,
            "ecutrho": 400,
            "occupations": "fixed",
        },
        "electrons": {
            "conv_thr": 1.0e-6,
            "mixing_beta": 0.3,
        },
        "ions": {
            "ion_dynamics": "bfgs",
        },
    },
    kpts=(4, 4, 6),
)

atoms.calc = calc

# ===== 4. 構造最適化の実行 =====
energy = atoms.get_potential_energy()
print("Relax 後の全エネルギー (eV):", energy)

# ===== 5. 最終構造の読み出し =====
relaxed_atoms = read("espresso.pwo")

print("\n最適化後の分率座標")
for symbol, pos in zip(relaxed_atoms.get_chemical_symbols(), relaxed_atoms.get_scaled_positions()):
    print(symbol, pos)

write("SnO2_relaxed.cif", relaxed_atoms)
print("\n最適化後の構造を SnO2_relaxed.cif に保存しました")

以下のように出力にエラーがなければ成功です。

元のCIFでは、空間群から決まる等価位置を省いて、非等価原子だけが書かれています。そのため、CIFにおけるSnもOも代表的な一つの位置だけが書かれていました。一方、今回の最適化後のCIFは、対称性情報を付けずに全原子を明示して保存しているため、見かけ上は空間群P1の構造として出力されます。必要に応じて、あとからpymatgenで空間群を再判定し、対称性付きCIFに戻します。

保存されたファイルを取得しましょう。Google Colabでは/contentに保存されます。左側にあるフォルダマークを押し、一つ上の階層に行くと「content」というフォルダが見えると思います。その中にあるSnO2_relaxed.cifをダウンロードしてください。

構造最適化後のCIF

VESTAでSnO2_relaxed.cifを開きます。すると、nonOptから酸素原子位置が改善され、元の良い構造に戻ったような見た目になっています。

構造の空間群を復元する

計算結果として得られるCIFは、対称性情報を含まない形で全ての原子位置が書かれるため、見かけ上はP1（空間群番号1）の構造になります。しかし、この原子配置は実際には元の空間群の対称性を保っています。Rietveld解析やSimulated Annealingなどの初期値として用いる場合には、この原子配置が持っている空間群を再び判定し、対称性を反映した構造に戻したいことがあります。

また、結晶構造には軸の取り方などにより複数の等価な書き方があります。そのため、ある基準に基づいて構造の表現をそろえる「標準化」という処理を行うことがあります。結晶学ではStructure tidyと呼ばれる方法がよく用いられます。

pymatgenは計算材料学向けのライブラリであり、バンド構造計算などで扱いやすい標準的なセル表現に構造を変換する機能を持っています。結果としてStructure tidyと似た形になることもありますが、アルゴリズムとしては同一の処理ではありません。

以下ではpymatgenを用いて空間群を判定し、対称性を反映した構造を作成します。

from pymatgen.core import Structure
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from pymatgen.io.cif import CifWriter

# ===== 1. 最適化後のCIFを読み込む =====
structure = Structure.from_file("SnO2_relaxed.cif")

# ===== 2. 空間群解析 =====
symprec = 1e-2
angle_tolerance = 5

sga = SpacegroupAnalyzer(
    structure,
    symprec=symprec,
    angle_tolerance=angle_tolerance
)

print("Detected space group:", sga.get_space_group_symbol(), sga.get_space_group_number())

# ===== 3. 対称性を反映した構造を作成 =====
refined_structure = sga.get_refined_structure()

# 対称性付きCIFとして保存
writer_refined = CifWriter(
    refined_structure,
    symprec=symprec,
    angle_tolerance=angle_tolerance,
    refine_struct=True
)
writer_refined.write_file("SnO2_relaxed_symmetrized.cif")

print("対称性付きCIFを保存しました: SnO2_relaxed_symmetrized.cif")

# ===== 4. 標準化した conventional cell を作成 =====
standard_structure = SpacegroupAnalyzer(
    refined_structure,
    symprec=symprec,
    angle_tolerance=angle_tolerance
).get_conventional_standard_structure()

# 標準化CIFとして保存
writer_standard = CifWriter(
    standard_structure,
    symprec=symprec,
    angle_tolerance=angle_tolerance,
    refine_struct=True
)
writer_standard.write_file("SnO2_relaxed_standardized.cif")

print("標準化CIFを保存しました: SnO2_relaxed_standardized.cif")

# ===== 5. 確認表示 =====
print("\n--- refined structure ---")
print(refined_structure)
print("Space group:", refined_structure.get_space_group_info())

print("\n--- standardized structure ---")
print(standard_structure)
print("Space group:", standard_structure.get_space_group_info())

このようにすると、酸素原子の位置は(0.195, 0.805, 0)になります。最初の(0.195, 0.195, 0)と違うように見えますが、これは分率座標として(0.195, -0.195, 0)と同じで、空間群136のWyckoff位置4fでは(x, x, 0)と(x, -x, 0)が等価なので実際には同じです。
うまく空間群136 が出ない場合はsymprecを少し大きくするとよいです。

おわりに

今回は自分の解析対象をQEで原子位置の最適化をするための一連の流れをSnO₂を用いて説明しました。まず道具としてQEを使えるようにすることを念頭に解説したので、計算条件などには何も触れませんでした。しかし、道具のことをもっとよく知ると自分の目的に合うように使いこなせるようになっていくので、私の勉強も含めいつかまとめたいと思います。

実際問題としては、ルチル型SnO₂のような簡単な構造は粉末X線構造解析で決まることがほとんどです。もっと複雑な、例えば対称性が低いとか、多様な金属元素からなる多元系酸化物であるような粉末X線構造解析で決めにくいような物質の場合は、QEのような第一原理計算では時間がかかりすぎGoogle Colabでは時間切れまたはリソース不足になることもあります。QEの計算条件を緩和するか、第一原理計算ではない方法が使える場合もありますので、次回以降解説したいと思います。