計算化学にPythonとPsi4で入門

01_計算化学

近年のコンピューターの高速化に伴い,実験化学者であっても高度な計算を手元の計算機で実行可能な環境が整いつつあります。

計算化学とは「コンピュータを使って化学の問題を解く」学問領域で,

  • 分子力学法(MM)
  • 分子軌道法(MO)
  • 密度汎関数法(DFT)
  • 分子動力学法(MD)
  • モンテカルロ法(MC)

などの方法が存在します。ここでは計算化学とは特に量子化学計算(MOとDFT)のことを意味していくことにします。

量子化学計算といえば,従来では

  • 高性能のワークステーション
  • 高価なソフトウェア

が必要でした。

前者の計算資源については,現在でもスペックが高ければ高いにこしたことはありません。

もっとも本サイトで取り上げる例は「普通のラップトップコンピュータで計算可能な程度の大きさの分子と計算レベル」で行うこととし,

  • 基本的な量子化学計算のやりかた
  • 結果の解析方法

を身につけることに重点を置いていきます。

しかし,中程度以上の系をある程度の精度で計算したいと考えた場合には,10コアのCPUをフルに使っても数日かかることが普通にありますので,それなりの計算リソースを使える環境を準備することをおすすめします。

後者の計算用のソフトウェアについては,現在でもGaussianをはじめとした有料ソフトウェアのシェアは大きいですが,状況が変わりつつある兆しがあるように見えます。

巨大なソフトウェアに最新の計算手法を新たに実装するコストが年々大きくなりつつあり,もっとモジュラー型の開発形態が好まれつつあります。中でもPsi4は

  • オープンソースのライブラリーで開発が非常に活発
  • 高価なソフトウェアに劣らない機能が無料で利用可能
  • 初学者にも優しいプログラミング言語であるpythonを使って動かすことができる

などといった理由から近年人気を集めている量子化学計算用のプログラムです.特に最近のアップデータからよりpythonicなコードが書けるようになり,pythonを用いた他の科学技術系ライブラリーとの連携が容易になりました。

本サイトでは,これからいくつかの記事を通してPsi4の使い方を基本から解説していきたいと思います。

本ブログのPsi4関係の記事は「Psi4の使い方」という記事にまとめています。参照してみてください。
【まとめ】Psi4の使い方
本ブログではPythonのオープンソースライブラリであるPsi4を用いて,量子化学計算の初心者を対象にした記事を多数用意しています。 記事の内容を一つずつ習得することで,知識・経験ゼロの状態から徐々にステップアップしていくことで,量子化学計算の基本コンセプトが習得できるように構成...

量子化学計算用のソフトウェア

量子化学計算用のソフトウェアは有償・無償のものを含めてさまざまな種類があります.ここではそのうちいくつかのものをまとめてみました。

それぞれのソフトウェアで共通してできることもありますが,ソフトごとに向き・不向きの計算方法などもあります。いずれにしても背景となる理論(量子化学)は同じですので,1つのソフトの使い方に習熟すれば,他のソフトへの移行も比較的容易だと思います。

Psi4とは

Psi4(読み方:サイフォー)は,量子化学計算用パッケージPsiの最新リリースです。GitHub上でLGPLライセンスにて,オープンソースソフトウェアとして開発されています。

Psi4の機能と特徴

機能としては(ポスト)ハートリー・フォック法などの分子軌道法や密度汎関数法などの手法を用いて

  • エネルギー計算
  • 構造最適化計算
  • 振動数計算

などの一連の代表的な計算が可能になっています。

入力ファイルの作成は上級者向けの細かい設定が可能な一方で,多くのオプションを自動で判別・類推してくれるため,最小限の入力情報で計算が実行な可能な初心者に優しい設計になっています。

速度が必要な計算の主要部分はC++でコーディングされており,APIによってpythonからの利用が可能です。

コア開発チームにDavid Sherrillの研究グループが含まれており,彼らの研究対象である非結合性相互作用の分析においては有償ソフトを凌ぐ高度な計算手法が実装されています。この機能を使うためだけでも,他のソフトのユーザーもPsi4の使用方法を覚える価値がある,と個人的には思います。実際,他のソフトで最適化計算を行った構造を用いてpsi4で利用可能な計算手法を用いている論文がかなり見られます。

Psi4の学習環境

公式ドキュメントが大変充実しており,サンプル用のインプットファイルもたくさん用意されています。さらに公式のyoutubeチャンネルにも多数の動画(スクリーンキャスト)がアップロードされていますので,英語に苦手意識がない方であれば,学習しやすい環境が用意されていると言えます。

Psi4と量子コンピューター

量子コンピューターの開発によって恩恵を受ける分野の1つと期待されているのが「量子化学計算」です。Googleの研究者を中心として発表された「OpenFermion」というプロジェクトは,量子コンピューターの専門家(開発者)と化学者(利用者)の間を取り持つもつようなプラットフォームを目指しています。

OpenFermion-Psi4」というプラグインが用意されており,量子計算最適化のための準備をpsi4を用いて実施可能になっています。

OpenFermionについては以下の記事も参照してみてください。Google AI Blog「Announcing OpenFermion: The Open Source Chemistry Package for Quantum Computers

Psi4のインストール

公式ドキュメントに記載の方法に従って行うことで容易にpsi4をインストール可能です。他の科学技術系ライブラリーと同様にcondaを使うのが便利だと思います.Windows WSL,Linux,Mac OS(Intel, Apple Slicon)の各環境で容易にインストールが可能です。

以下の例ではpython3.10の環境にインストールしています.適宜設定を変えてみてください。

conda install psi4 python=3.10 -c conda-forge
OSのバージョン違いなどによって手元の環境でPsi4のインストールがうまくいかないことも多いようです。「Google Colabでpsi4:量子化学計算用のpython環境を手軽に構築」という記事ではGoogle Colaboratory上でpsi4の環境構築を行う方法を解説しています.参照してみてください。
Google ColabでPsi4:量子化学計算用のpython環境を手軽に構築
本ブログでは「有機合成化学者のための計算化学入門」を掲げて,特にpythonを用いて量子化学計算を行う際に必要となる環境構築方法から解説してきました.例えば「計算化学にpythonとPsi4で入門」という記事では,pythonの量子化学計算用ライブラリであるPsi4のインストール...

Psi4の2つの使い方

Psi4のインストールが終わると

  • psi4という実行可能なプログラム
  • psi4というpythonライブラリー

の2つが利用可能になります。それぞれ「Psithon」と「PsiAPI」という名前で公式ドキュメントでは呼ばれています。

Psithon

実行可能プログラム「psi4」を使って計算を実行する方法になります。

インプットファイルの記述法自体は独自形式になりますが,各所にpython風の文法で追加の処理を記述することが可能になっています。後からpythonをサポートした科学技術系ソフトでよく見られる記述方法だと思います。

当然pythonそのものではないので,いわゆるpythonicではない記述方法も多くなります。

具体的には以下のように用います.

psi4 -n 2 -i molecule.in -o molecule.log

この場合は

  • 2コアのCPUを利用
  • molecule.inを入力ファイル
  • molecule.logを出力ファイル

として設定して計算を行います。入力ファイルには分子構造や計算手法の記述が含まれています。

PsiAPI

「PsiAPI」ではpythonスクリプトにPsi4ライブラリーをimportして用いることになります。Psi4で用いる全てのクラス・メソッドは「psi4.XXX」でアクセス可能です。

計算結果をそのままpythonオブジェクトとして受け取れることから,Psi4の計算結果をそのまま他のプログラムに渡しやすい設計になっています。例えば「Psi4NumPy」というプロジェクトも公式から派生しています。

もっとも数日かかるような重たい計算を行うことが多い場合には,ログファイルで計算の途中経過を確認しながら実行することが大半だと思いますので,自動化の恩恵は少ないかもしれません。

PsiAPIを用いた場合には,具体的には以下のように実行します。

python molecule.py

molecule.pyには「import psi4」からはじまる一連のpsi4関係のコードによって

  • 計算環境(コア数,メモリ,ログファイルなど)の設定
  • 分子構造の設定
  • 計算レベルの設定
  • 計算方法の設定

といった内容が記述されています。

上記のようにターミナルからpythonスクリプトを実行する以外にも,Jupyter Notebookを用いてセルごとに実行していくことも当然可能です。

本サイトでは「PsiAPI」を用いた方法を中心にして,具体的な計算のやり方などを見ていきたいと思います。

入力ファイルの構造

PsiAPIの項で見たように入力ファイルには,

  1. 計算環境の指定
  2. 分子構造の指定
  3. 計算レベルの指定
  4. 計算方法の指定

といった内容が含まれています。

上級者はこれら項目について入力ファイルで細かい設定が可能になっています。一方で,Psi4では多くのオプションを自動で判別・類推してくれます。そのため,初心者は上記項目すべてを入力する必要は必ずしもなく,最小限の入力情報で計算を開始することが可能です。

それぞれの項目について詳細にみていきましょう。

計算環境の設定

psi4.set_num_threads(nthread=n)
psi4.set_memory(‘memory’)
psi4.set_output_file(output-file)

計算環境の設定はpsi4.set_XXXで行えます。

計算に用いるCPUのスレッド数やメモリの設定は計算速度に直結する要素です。また量子化学計算において計算のログファイルは,きちんと計算が実行されたかどうかを確認するための大切な判断材料になります。デフォルトでは標準出力にログが吐き出されてしまいますので,計算の途中結果を保存するファイルをきちんと設定しておきましょう。

下記の例では

  • 使用するスレッド:2
  • メモリ:4GB
  • ログファイル:UNIX時間を基にしたファイルネーム

として設定しています.現実的にはログファイルには計算に用いる分子の名前や計算手法も入れておいた方が後から見返す際に便利だと思います。

import psi4
import datetime
import time

t = datetime.datetime.fromtimestamp(time.time())

psi4.set_num_threads(nthread=2)
psi4.set_memory('4GB')
psi4.set_output_file('{}{}{}_{}{}.log'.format(t.year,
                                              t.month,
                                              t.day,
                                              t.hour,
                                              t.minute))

分子構造の指定

psi4.geometry(structure)

分子構造の設定では三重引用符(’’’)を用いて各原子の位置情報をXYZ座標Z-マトリックスを用いて行うのが普通です。その際に指定する分子構造を変数に格納しておき,後ほど計算を実行する際にどの構造に対して計算を行うかを指定します。

構造情報は通常は1行目には分子の電荷とスピン多重度を指定し,2行目以降に構造情報を指定します。

電荷) (スピン多重度
XYZ座標またはZ-マトリックス
・・・

以下に水分子(H2O)の例を示します。

### XYZ
h2O_xyz = psi4.geometry('''
0 1
O 0 0 -0.11
H 0 -1.4 1.2
H 0 1.2 1.2
''')

### Z-matrix
h2O_zmat = psi4.geometry('''
0 1
O
H 1 0.96
H 1 0.96 2 104.5
''')
Z-マトリックス等の分子の三次元情報を扱うフォーマットについては,「XYZ形式とZ-マトリックスは分子の立体構造を表す入力フォーマット」という記事で解説しています。参照してみてください。

計算方法の設定

psi4.set_options(dict)

計算方法の設定は大きく,

  1. 計算するレベルなど設定(DFT汎関数,基底関数系など)
  2. 計算種類の設定(エネルギー計算など)

にわけられます。

psi4.set_optionsでは,様々なオプション項目を辞書型の引数で設定可能です。今後の記事で徐々にオプションの種類などを見ていきたいと思います。

エネルギー計算を実行してみよう

psi4.energy(level-of-theory, molecule)

今回は計算の種類としてエネルギー計算を選択し,いくつかの計算レベルにて実際に計算してみましょう。

下のコードでは,さきほどZ-マトリックスで構造を指定したh2O_zmatに対して,

  1. HF/STO-3G
  2. HF/3-21G
  3. B3LYP/STO-3G
  4. B3LYP/3-21G

の4つの計算レベルでエネルギー計算を行い,その結果をprint文で表示しています。このように計算結果を直接pythonオブジェクトとして扱える点がPsiAPIの利点です。

theory = ['hf', 'b3lyp']
basis_set = ['sto-3g', '3-21G']

import itertools

for th, ba in itertools.product(theory, basis_set):
    level = th + '/' + ba
    e = psi4.energy(level, molecule=h2O_zmat)
    print('energy at the {} level of theory:\t{:.4f}'.format(level, e))

計算が実行されると以下のように表示されるとおもいます。

energy at the hf/sto-3g level of theory:    -74.9634
energy at the hf/3-21G level of theory: -75.5854
energy at the b3lyp/sto-3g level of theory: -75.3135
energy at the b3lyp/3-21G level of theory:  -75.9718

この結果から,同じ構造でも計算レベルによってエネルギーの値が異なることがわかると思います。すなわち,分子のエネルギーを比べる場合には完全に同じ計算手法・レベルで行った結果を比較する必要があります。モデルケミストリーを揃える必要があるともいいます。

続いて構造を少しだけ変えた水分子について,同じ計算レベルで計算を行ってみましょう。

次のコードではO–H結合の長さを変化させ,HF/STO-3Gレベルでエネルギー計算を行っています。

geom = '''
0 1
O
H 1 {0}
H 1 {0} 2 104.5'''

for i in [0.85, 0.90, 0.95, 1.00, 1.05]:
    h2O = psi4.geometry(geom.format(i))
    e = psi4.energy('hf/sto-3g', molecule=h2O)
    print('energy for d = {}: {:.4f}'.format(i, e))

計算が実行されると以下のように表示されると思います。

energy for d = 0.85: -74.9100
energy for d = 0.9: -74.9451
energy for d = 0.95: -74.9618
energy for d = 1.0: -74.9648
energy for d = 1.05: -74.9572

このことから,計算レベルを揃えても,分子構造が異なればエネルギーが異なることがわかります。すなわち,エネルギーは分子構造,分子の自由度(3N-6),を変数とした関数になります。

終わりに

今回は「計算化学にPythonとPsi4で入門」という話題について,

  • Psi4の特長とインストール方法
  • Psi4を使った簡単なエネルギー計算の方法

について見てきました。

実際に水分子のエネルギー計算を行いながら,

  1. エネルギーの値は用いる計算レベルによって異なる
  2. エネルギーの値は構造によって異なる

ということを確認しました。これら項目はこれから量子化学計算を実行していくにあたり何度も出てくる大切なコンセプトですので,きちんと理解しておきましょう。

エネルギーが分子の自由度の関数ということは,エネルギーが小さくなるように構造を徐々に変化させていくことで安定構造を求めることができそうです。これが「構造最適化計算」とよばれるものの基本的な考え方になります。

次回はPsi4を用いて簡単な構造最適化計算を行いながら,量子化学計算についてさらに理解を深めていきたいと思います。

>>次の記事:「計算化学の構造最適化の基本をPsi4で学ぶ

コメント

  1. okuracoffee より:

    いつもためになる質問をありがとうございます!
    質問なのですが、psi4.energy(level-of-theory, molecule)で出力される値の単位は
    「a.u」なのでしょうか?

    • 管理人 より:

      okuracoffee様

      ご質問ありがとうございます.
      出力される単位はご明察の通り「a.u.」になります.

      今後ともよろしくお願いいたします.

  2. あかね より:

    非常に内容の濃い情報をたくさん提供していただきありがとうございます。
    https://future-chem.com/psi4-intro/
    「計算環境の設定」の節のサンプル・プログラムで、末尾の t.minute) の後に ) が足りません。

タイトルとURLをコピーしました