データ分析とケモインフォ

Pythonで量子化学計算（Psi4）

オープンソースの量子化学計算ソフトウェア（Psi4）の使用方法についてのまとめ。

使い方（主にコード）は、化学の新しいカタチ、ご注文はリード化合部ですか？を参考にさせていただきました。ありがとうございました。
（詳しい使い方などの情報は、上記HPの方が充実しています）

やったこと（目的）

アスピリン（アセチルサリチル酸）に対して、構造最適化計算を行い、以下の情報を取得した。

  • 最適化構造の構造情報（各原子の座標）
  • エネルギー
  • HOMO, LUMO のエネルギー
  • Mulliken電荷
  • 双極子モーメント

前準備

モジュールのインストール

conda install psi4 psi4-rt python=3.7 -c psi4
conda install -c rdkit rdkit

コード（構造最適化まで)

モジュールのインポート

import numpy as np
import pandas as pd
from rdkit import Chem
from rdkit.Chem.rdDistGeom import ETKDGv3, EmbedMolecule
from rdkit.Chem.rdForceFieldHelpers import MMFFHasAllMoleculeParams, MMFFOptimizeMolecule
from rdkit.Chem.Draw import IPythonConsole
import psi4
import datetime
import time

計算の設定~構造最適化まで

%%time
#計算時間を見てみる
# ハードウェア側の設定（計算に用いるCPUのスレッド数とメモリ設定）
psi4.set_num_threads(nthread=1)
psi4.set_memory("1GB")

# 入力する分子（アセチルサリチル酸）
smiles = 'CC(=O)Oc1ccccc1C(=O)O'

# ファイル名を決める
t = datetime.datetime.fromtimestamp(time.time())
psi4.set_output_file("{}{}{}_{}{}_{}.log".format(t.year,
                                                t.month,
                                                t.day,
                                                t.hour,
                                                t.minute,
                                                smiles))

# SMILES から三次元構造を発生させて、粗3D構造最適化
mol = Chem.MolFromSmiles(smiles)
mol = Chem.AddHs(mol)
params = ETKDGv3()
params.randomSeed = 1
EmbedMolecule(mol, params)

# MMFF（Merck Molecular Force Field） で構造最適化する
MMFFOptimizeMolecule(mol)

#UFF（Universal Force Field）普遍力場で構造最適化したい場合は
#UFFOptimizeMolecule(mol)
conf = mol.GetConformer()

# Psi4 に入力可能な形式に変換する。
# 電荷とスピン多重度を設定（下は、電荷０、スピン多重度1)
mol_input = "0 1"

#各々の原子の座標をXYZフォーマットで記述
for atom in mol.GetAtoms():
    mol_input += "\n " + atom.GetSymbol() + " " + str(conf.GetAtomPosition(atom.GetIdx()).x)\
    + " " + str(conf.GetAtomPosition(atom.GetIdx()).y)\
    + " " + str(conf.GetAtomPosition(atom.GetIdx()).z)
molecule = psi4.geometry(mol_input)

# 計算手法（汎関数）、基底関数を設定
level = "b3lyp/6-31G*"

# 計算手法（汎関数）、基底関数の例
#theory = ['hf', 'b3lyp']
#basis_set = ['sto-3g', '3-21G', '6-31G(d)', '6-31+G(d,p)', '6-311++G(2d,p)']

# 構造最適化計算を実行
energy, wave_function = psi4.optimize(level, molecule=molecule, return_wfn=True)
'''output
Optimizer: Optimization complete!
CPU times: user 15min 50s, sys: 27.5 s, total: 16min 17s
Wall time: 17min 18s
'''

スレッド数 １、メモリ １GB 使用して、1分子計算するのに計17分かかった。

ちなみに入力するファイル(XYZ形式）はこんな感じ

0 1
C 2.974609598387524 0.6502198947117926 -1.4939659177438849
C 1.5230890246254225 0.5053737824529224 -1.1530242528515486
O 0.603330778130228 1.057538083744091 -1.7483070198697948
O 1.387914769521494 -0.30271498078404824 -0.03263846326379664
C 0.054417291055592946 -0.5712371294324196 0.2945337836815555
C -0.4224582203026945 -1.8473553693968503 -0.01601581679481116
C -1.7257359235711085 -2.203041426185137 0.3353284417639963
C -2.5433749515038557 -1.2955508196005052 1.0118502341281146
C -2.060774394293738 -0.026729946932774716 1.3386540927479436
C -0.7608195905359867 0.34046831642295106 0.9729716175553647
C -0.29830011960960084 1.6797421182462917 1.3678608910543
O -0.14779504472207144 2.009190935828314 2.530629328049564
O -0.14453966540626556 2.510886086741464 0.31763612590410073
H 3.485468003428318 1.2068580777799764 -0.7047280552402798
H 3.4260440775599967 -0.33672600434994854 -1.6236316858587398
H 3.071579924882846 1.198336303548168 -2.4353976902946775
H 0.21564436873480777 -2.565765313412947 -0.5236197408913651
H -2.1022344085524307 -3.194413587471996 0.09265503994624906
H -3.5541301454100207 -1.5811057139723368 1.2953607541830348
H -2.6935167212816706 0.6706025040826632 1.8845251711162185
H -0.28841865113707016 2.0954241879804125 -0.5663389837709468

コード（最適化構造からの情報取得）

構造情報（XYZ形式）

print(molecule.save_string_xyz())
#単位はオングストローム
'''
0 1
 C    3.098884439155    0.286778541532   -1.818971930267
 C    1.631416601050    0.255588755108   -1.492773907603
 O    0.755324426944    0.834934357982   -2.102478251457
 O    1.396466682379   -0.489685966862   -0.380657317744
 C    0.074752427021   -0.810467114907   -0.015892403575
 C   -0.343553415622   -2.104118243082   -0.317545006753
 C   -1.590020749705   -2.548730143826    0.117945091675
 C   -2.402170016296   -1.705115088782    0.878418776055
 C   -1.963835487350   -0.422144280017    1.191561137714
 C   -0.729128903046    0.059858926344    0.731238467127
 C   -0.357527633804    1.464587591984    1.139219316977
 O   -0.521525170787    1.863006769904    2.268651582597
 O    0.121812264509    2.275304342291    0.176882785876
 H    3.636482218213    0.795645992511   -1.011210543547
 H    3.498253279532   -0.729173607712   -1.889712382197
 H    3.251358031652    0.823578110904   -2.755331437125
 H    0.322818412683   -2.753052223004   -0.876772337761
 H   -1.916523030934   -3.555688162495   -0.124828509423
 H   -3.368180133191   -2.049556108518    1.235066217610
 H   -2.570453292820    0.238058474783    1.802873381996
 H    0.073192650040    1.843321147028   -0.703222197850
'''

エネルギー、HOMO、LUMOエネルギー

#単位は原子単位（a.u. もしくはHartrees)
#エネルギー
print(round(energy,3),'a.u.')
# >>> -648.689 a.u.

# HOMO を表示（単位： au = Hartree）
LUMO_idx = wave_function.nalpha()
HOMO_idx = LUMO_idx - 1
homo = wave_function.epsilon_a_subset("AO", "ALL").np[HOMO_idx]
lumo = wave_function.epsilon_a_subset("AO", "ALL").np[LUMO_idx]
print('homo:',round(homo,5),' a.u.')
print('lumo:',round(lumo,5),' a.u.')
# >>> homo: -0.26022 a.u.
# >>> lumo: -0.04309 a.u.

Mulliken電荷

psi4.oeprop(wave_function, 'MULLIKEN_CHARGES')
mulliken = np.array(wave_function.atomic_point_charges())
for i, atom in enumerate(mol.GetAtoms()):
    print(atom.GetSymbol(),round(mulliken[i],4))

# 類似度マップで可視化してみる
from rdkit.Chem.Draw import SimilarityMaps

fig = SimilarityMaps.GetSimilarityMapFromWeights(mol,
mulliken,
colorMap='RdBu')
'''output
C -0.527
C 0.6015
O -0.4756
O -0.5025
C 0.2962
C -0.147
C -0.1334
C -0.1206
C -0.1705
C 0.0273
C 0.5345
O -0.4395
O -0.5899
H 0.2051
H 0.1907
H 0.1951
H 0.1468
H 0.1476
H 0.1462
H 0.1705
H 0.4445
'''

Mulliken密度解析では基底関数系にSTO-3Gなどの最小基底系を利用することを想定しているため、大きな基底関数系を用いた場合には正しくない電荷分布が得られることがあるとのこと。

双極子モーメント

双極子モーメントの各座標成分から、分子全体の双極子モーメントの大きさを求める。

dipole_x, dipole_y, dipole_z = psi4.variable('SCF DIPOLE X'), psi4.variable('SCF DIPOLE Y'), psi4.variable('SCF DIPOLE Z')
dipole_moment = np.sqrt(dipole_x ** 2 + dipole_y ** 2 + dipole_z ** 2)

#単位はデバイ(D,debye)
print(round(dipole_moment,3),'D')
# >>> 5.175 D

計算手法について

  • HF法(SCF法)：シュレディンガー(Schrödinger)方程式の解である電子状態を近似波動関数について解く分子軌道法
  • 平均化された電子雲の中で運動している電子を扱っているので、電子間の反発を過小評価していることになる。
    (HF法の近似波動関数で解いた電子状態と実際の値との差を、電子相関と呼ぶ）
  • DFT計算:シュレディンガー方程式の代わりに、等価なコーン・シャム方程式(Kohn-Sham)という計算を解く
    • HF法では含まない交換相互作用と電子相関の相互作用が正確な解に対し部分的に考慮されている。
    • 圧倒的に速く計算ができるが、汎関数を選択するかが計算結果に大きな影響を与える
    • B3LYPがよく利用される（大きく実験値からハズレることが少ないため）

基底関数について

基底関数とは

電子がどのような大きさや形で、どの軌道にいるのか」 を数式で表現したもの
例） 6-31G：内殻軌道を6個のGauss型関数、原子価軌道を内側と外側とで3個と１個のGauss型関数に分けて表現しているという意味

基底関数の選び方（参考)

全ての系に万能な基底関数系は存在せず、計算目的に応じてその都度、目的にあった関数系を選ぶことが重要

  • 一般的な有機化合物であればまずは、6-31(d)を基準に選ぶ。（計算時間のわりにそこそこいい結果が出るらしい）
  • 分子量100以内であればB3LYPと6-31G(d)の組み合わせで、実用的な計算時間内で結果が出る。

補助関数（6-31G(d)のdの部分）について

  • d, p：分極関数(分子内での結合生成やイオン間の静電相互作用の分極の効果を考慮する)
    • d：水素原子以外の原子に分極関数を適用
    • p：水素原子に分極関数を適用
  • +, ++：diffuse関数（電子分布の広がりを考慮する）
    • +：H,Heよりも重い原子にdiffuse関数を適用
    • ++：軽い元素（H,He）にもdiffuse関数を適用

補助関数（6-31G(d)のdの部分）の使い分け

  • 一般的な有機化合物：6-31G(d)　(分極関数：dを追加)
  • 負電荷を持つ化学種(アニオン), 励起状態：6-31+G(d)　(diffuse関数：+を追加)
  • 正電荷を持つ化学種(カチオン)：6-31G(d,p)　(分極関数：dpを追加)

参考

  1. 計算化学にpythonとpsi4で入門(化学の新しいカタチ)
  2. 計算化学における電荷：psi4を用いた電子密度解析(化学の新しいカタチ)
  3. 計算化学の構造最適化の基本をpsi4で学ぶ(化学の新しいカタチ)
  4. Pythonで量子化学計算(ご注文はリード化合物ですか？)
  5. psi4公式ドキュメント
  6. psi4公式 youtube
  7. 基底関数はどれを選べばいいの?(PC CHEM BASISCS.COM)
  8. HF/3-21Gとは何なのか？(PC CHEM BASISCS.COM)
  9. B3LYP/6-31Gとは何なのか？(PC CHEM BASISCS.COM)

関連書籍