見出し画像

Google Colab上でQuantum Espressoを使う

Hiroyuki Hayashi

はじめに

以前、X線回折実験のデータを使って未知の結晶構造を解析する方法を解説しました。しかし、実際の研究現場では、ある一つの「壁」にぶつかることがよくあります。

それは、「陽イオンの位置は綺麗に決まったのに、酸素原子の位置がどうしても定まらない」という問題です。X線回折のような実験手法は、重い原子はX線を強く散乱するのでデータへの寄与が大きく、見つけやすいのですが、重元素が強く存在する系では酸素のような軽い原子は位置精度が落ちることがあります。

そんなとき、実験とは全く異なるアプローチである「計算科学」が強力な助っ人になります。それが今回紹介する第一原理計算です。

計算で実験の限界を補完する

結晶構造解析において、X線回折は実験室に置いて標準的な手法です。しかし、どの方法にも限界はあるもので、例えば、水素原子や酸素原子のような比較的軽い元素の配置を決定する際、その散乱強度の弱さがネックとなり、精度不足に陥ることがあります。重い陽イオンの位置がある程度決まっていても、酸素原子の具体的な座標が実験データだけでは特定できないことは珍しくありません。

そこで、回折データだけに拠らず原子位置を知るための強力な補助手法として用いられるのが「第一原理計算」です。

第一原理計算とは、量子力学の基本原理に基づき、原子や電子の振る舞いをコンピュータ上でシミュレーションする手法のことです。物質の構造が与えられたとき、量子力学の方程式を解くことにより、その構造が持つ全エネルギー $${E_{total}}$$ を算出します。

この手法を構造決定に応用する考え方は、非常にシンプルです。

  • 自然界に存在する安定な物質は、その条件下でエネルギーが最も低くなる（極小になる）構造をとります。第一原理計算は基本的に絶対零度の計算になりますが、酸素原子の位置を推定するという目的では室温でなくても十分であることが多いです。

  • 実験で位置が決まらなかった原子をコンピュータ上で様々に動かし、全エネルギーが最も小さくなる座標を探します。

  • こうして導き出された「エネルギーが最も低い構造」を、その物質の妥当な構造として採用するのです。

第一原理計算は、実験装置による測定データに直接依存せず、理論計算のみで物質の構造を予測できます。実験データと計算による予測を組み合わせることで、従来の解析では困難だった軽元素の位置特定が可能になるのです。

本稿では、この第一原理計算をGoogle Colab上で実際に実行し、物質の安定構造を計算する手順を解説していきます。

第一原理計算の詳しいことはこのnoteでは触れませんので、例えば下記の書籍や動画を参照してください。

Google Colabで始めてみよう

実際に計算を実行するには、量子力学の方程式を解くための非常に高性能なコンピュータと、それを動かすための専門的なソフトが必要です。

通常、こうした環境を個人で構築するのはハードルが高いものです。しかし、Google Colabを利用すれば、その課題は一気に解消されます。

Google Colabとは何か

Google Colab（正式名称：Google Colaboratory）は、Googleが提供している「ブラウザ上で動くPythonの実行環境」です。

本来、プログラミングや複雑な計算を行うには、自分のコンピュータにソフトをインストールしたり、環境設定という難しい作業が必要になります。しかし、Google Colabを使えば、ブラウザを開くだけで、Googleの巨大なサーバー上にある高性能なコンピュータを、誰でも無料で借りて利用することができます。

Google Colabを使うメリットは、大きく分けて2つあります。

  1. 環境構築の手間がゼロ: Quantum Espressoなどの専門的なソフトも、コマンド一つでサーバー上に準備できます。

  2. 計算パワーの確保: Googleのサーバー上で計算が行われるため、個人のノートPCのスペックに左右されず、高速な計算が可能です。

無料版の制約

Googleアカウントがあれば使えますが、無料版では以下の制約があります。

  • リソースは「保証」されていない: 常に同じ性能のコンピュータが使えるわけではありません。その時の利用状況や混雑具合によって、割り当てられる性能が変動します。

  • 計算が中断されることがある: 長時間の計算を行っている途中で、突然接続が切断されたり、リソースの再割り当てが行われたりする可能性があります。そのため、数時間かかるような計算を一度に行うことは適していません。

  • 放置は厳禁: ブラウザを開いたまま長時間放置すると、自動的に接続が終了します。計算を実行している間は、定期的に画面を確認する必要があります。

【重要：研究を止めないためのヒント】 これらの制約があるため、実際に本格的な計算を行う際は、「計算結果をGoogleドライブに随時保存する」という工夫が欠かせません。もし計算の途中で接続が切れても、保存されたデータから続きを再開できるようにしておくことが、計算科学を長く楽しむためのコツです。

あくまで「限られた資源をみんなで分かち合っている」という意識を持ち、小分けにして計算を実行する習慣をつけておきましょう。

計算を始めるための「作法」

Google Colabでの作業は、ノートブックという形式で行われます。PythonのJupyter Notebookのようなものです。この中には「コードセル」という枠があり、そこに命令文を入力して実行ボタンを押すだけで、順番に処理が進んでいきます。

本稿では、皆さんがこのノートブックに指示を出すだけで、自動的に「量子力学の方程式を解く準備」が整うようにコードを構成しています。専門的な知識がなくても、セルを順番に実行するだけで、誰でも計算の入り口に立つことができるのです。

Quantum Espresso：物質の電子状態を計算するソフトウェア

第一原理計算を実行する際、多くの研究者が世界中で利用しているオープンソースソフトウェアが「Quantum Espresso」です。

物理学の世界では、原子や電子の振る舞いを記述する「シュレーディンガー方程式」を解く必要があります。しかし、この方程式を厳密に解くことは計算機にとっても非常に困難です。そこでQuantum Espressoは、「密度汎関数法（Density Functional Theory, DFT）」という手法を用いて、電子同士の相互作用を近似的に扱い、現実的な時間で計算を終わらせる仕組みを持っています。

計算の基本的な流れ

Quantum Espressoは、以下の手順で計算を行います。

  1. 入力ファイルの作成: 計算したい結晶の形や、原子の種類、電子の密度を計算するための精度などの条件を書き込みます。

  2. 実行（計算）: pw.x というプログラム（エンジン）が、入力ファイルに基づき、量子力学の方程式を反復計算によって解きます。

  3. 結果の出力: 計算が終わると、その構造の「全エネルギー」や「原子に働く力」が記録された出力ファイルが生成されます。

Python（ASE）の利用

Quantum Espresso単体でも計算は可能ですが、研究現場では「ASE（Atomic Simulation Environment）」というPythonライブラリを組み合わせて使うのが一般的です。

直接Quantum Espressoの入力ファイルを書くこともできますが、手作業ではミスが起こりやすく、構造の変更も手間がかかります。ASEを使うことで、Pythonプログラムの中で原子の配置を定義し、その計算設定を自動的にQuantum Espressoへ受け渡すことができます。

今回のコードで行っているのは、この「Pythonで構造を作成し、Quantum Espressoへ渡して結果を受け取る」という橋渡しです。

用語の解説：擬ポテンシャル

今回のコードで「Si.pz-vbc.UPF」というファイルを使いました。これは「擬ポテンシャル（Pseudopotential）」と呼ばれるデータです。

原子核に近い電子（内殻電子）と原子核の動きをまとめて「一つの塊」として扱うことで、化学結合に関わる外側の電子（価電子）のみを効率よく計算するための工夫です。このデータがないと、原子内のすべての電子を計算することになり、膨大な計算時間が必要になってしまいます。第一原理計算においては、計算したい原子種ごとに、このデータを用意しておくことが必須となります。

このように、Quantum Espressoは「物理的な近似」と「計算効率」のバランスを取りながら、現実的な時間で物質の性質を導き出すためのツールです。次の章では、実際にコードを入力し、この仕組みを動かしてみましょう。

シリコンの安定構造を計算

それでは、実際にGoogle Colab上でプログラムを動かし、シリコン（Si）の安定構造を計算してみましょう。ここでの目的は、量子力学の方程式に基づき、ダイヤモンド構造を持つシリコンの全エネルギー$${E_{total}}$$ を算出することです。

1. 計算環境の準備

まず、プログラムの最初のブロックで、必要なソフトウェアとライブラリをインストールします。

# Quantum ESPRESSOと関連ライブラリのインストール
!apt-get update -qq
!apt-get install -y quantum-espresso

# ASE（Atomic Simulation Environment）の導入
!pip -q install --upgrade ase

# 擬ポテンシャルファイルのダウンロード
!wget -q https://gitlab.com/QEF/q-e/-/raw/3e920635bff849e891bb3fcec99acc096c3bdcf6/test-suite/epw_pl/Si.pz-vbc.UPF

ここでは、量子力学計算を行うエンジン「Quantum Espresso」と、Pythonからそれを制御するためのツール「ASE」を準備しています。

2. 計算の設定と実行

次に、結晶構造を定義し、計算の条件を指定します。

from ase.build import bulk
from ase.calculators.espresso import Espresso, EspressoProfile

# シリコンの結晶構造（ダイヤモンド構造、格子定数 a=5.43 Å）を定義
atoms = bulk("Si", "diamond", a=5.43)

# 擬ポテンシャルの指定
pseudopotentials = {"Si": "Si.pz-vbc.UPF"}

# 計算条件の設定
profile = EspressoProfile(command="pw.x", pseudo_dir=".")

calc = Espresso(
    profile=profile,
    pseudopotentials=pseudopotentials,
    input_data={
        "control": {"calculation": "scf"},
        "system": {"ecutwfc": 30, "ecutrho": 240},
        "electrons": {"conv_thr": 1e-8},
    },
    kpts=(4, 4, 4),
)

atoms.calc = calc
energy = atoms.get_potential_energy()
print("Total energy (eV):", energy)

3. 設定値の解説

コード中の数値は、計算の精度とコストのバランスを決める重要なパラメータです。

  • ecutwfc (30 Ry) / ecutrho (240 Ry):

    1. これは「平面波基底」という関数のカットオフエネルギーです。電子の波動関数をどれだけ細かく分解するかを決定します。数値を大きくするほど精度は上がりますが、計算時間は増大します。

  • kpts (4, 4, 4):

    1. 結晶の周期性を考慮するためのK点グリッドです。物質の電子状態を空間的にサンプリングする細かさを表します。

  • calculation: "scf":

    1. 「自己無撞着場（Self-Consistent Field）計算」の略です。電子密度が変化しなくなるまで繰り返し計算を行うモードで、基底状態のエネルギーを求める際の標準的な手法です。

エネルギーの単位「Ry（リュードベリ）」とは？

1 Ryは、水素原子の基底状態におけるイオン化エネルギー（約13.6 eV）に相当します。原子や分子中の電子の運動を扱う際、この単位を用いると、方程式を簡潔な形で表すことができます。

もしRyの感覚が掴みにくい場合は、1 Ry ≒ 13.6 eVという変換式を覚えておくと便利です。

エネルギーの単位：Ry

4. 計算結果の読み方

このコードを実行すると、最後に Total energy (eV): -xxxx.xxxx という数値が出力されます。これが、入力した格子定数 $${a}$$ = 5.43 Å におけるシリコンの全エネルギーです。

第一原理計算の真髄は、この計算を格子定数 $${a}$$ を変えながら何度も繰り返し、「どの $${a}$$ の値のときにエネルギーが最小になるか」を探すことにあります。その最小値を与える格子定数が、その物質の最も安定な構造であると理論的に予測されるのです。

まずはこのコードを動かし、シリコンという代表的な物質のエネルギーが得られることを確認してください。これが構造解析における計算科学の第一歩となります。

おわりに

Google Colabを活用した第一原理計算の基本と、Siの安定構造を求めるプロセスを解説してきました。

回折実験は、物質の結晶構造を解き明かすための極めて強力な手段です。しかし、そこには「見えにくいもの」や「曖昧なもの」が存在します。実験データだけで構造を完全に決定できないとき、今回の手法のような「計算による予測」を組み合わせることは、現代の研究における標準的なアプローチとなっています。

今回皆さんが実行したシリコンの計算は、「計算科学の世界の入り口」です。

本来、酸素の位置がうまく決まらないような複雑な結晶構造解析においても、考え方は同じです。実験から得られた陽イオンの位置を固定し、酸素原子をさまざまな位置に配置して構造最適化計算を行いそれぞれの全エネルギーを計算する。そうしてエネルギーが極小になる場所を探すことで、回折データからは見えなかった酸素の正しい位置を、理論という名の根拠に基づいて導き出すことができます。

第一原理計算は、単なる計算ソフトではありません。皆さんが持っている「なぜこの原子はこの場所に並んでいるのか？」という問いに対する答えを探すための、もう一つの実験装置です。

次回以降、実際に結晶構造ファイルCIFを入力とする計算方法に関して説明していきます。

いいなと思ったら応援しよう！