計算化学における電荷:psi4を用いた電子密度解析

化学

計算手法とエネルギー・最適化構造の関係:コンフォメーション探索における注意点」という記事では,psi4を用いて様々な計算手法を用いてエネルギー計算を行うことで,

  • 計算手法によって考慮できる相互作用が異なる

ことを説明しました.その結果

  • 構造最適化できないことがある
  • コンフォマー間のエネルギー差が逆転する場合がある

ことなどを具体例を通して見てきました.

ところで,エネルギー計算により得られるのは分子のエネルギーだけではありません.今回はシュレディンガー方程式の解である波動関数から得られる,様々な情報について見ていきたいと思います.

今回も量子化学計算のエンジンとしてpsi4を用いています.本ブログではpsi4を用いたエネルギーと構造最適化の基本について「計算化学にpythonとpsi4で入門」や「計算化学の構造最適化の基本をpsi4で学ぶ」という記事で解説しています.適宜参照してみてください.

計算化学にpythonとpsi4で入門
近年のコンピューターの高速化に伴い,実験化学者であっても高度な計算を手元の計算機で実行可能な環境が整いつつあります.計算化学とは「コンピュータを使って化学の問題を解く」学問領域で,分子力学法(MM)分子軌道法(MO)密度汎関数法(DFT)分子動力学法(MD)モン...
計算化学の構造最適化の基本をpsi4で学ぶ
「計算化学にpythonとpsi4で入門」という記事ではpythonで扱える量子化学計算ソフトウェア「psi4」の紹介と,簡単なエネルギー計算のやり方を扱いました.その際に得られるエネルギーは用いる計算レベルによってエネルギーの値が異なる計算する構造によってエネルギー...

電子密度解析とは何か

有機化学反応の多くは,「カチオン」と「アニオン」,「正電荷」と「負電荷」によって理解されています.例えばカルボニル基の求電子性は「C=O」の結合の電荷の偏りによって説明されます.

このような形式的な「電荷」を,波動関数の絶対値の二乗として得られる電子密度を用いて各原子に割り当てることを密度解析(population analysis)といいます.

割り当てるルールによっていくつかの方法が知られていますが,よく用いられるのは

  • Mulliken(マリケン)密度解析
  • Löwdin(ローディン)密度解析

の2つで,どちらもpsi4で利用可能です.

他の電荷としては

  • 自然密度解析(NPA)

や静電ポテンシャルへのフィッティングによって電荷を決める

  • Merz–Kollman電荷

などが有名です.

日本語で読める電子密度解析に関する解説としては,以下の資料(PDF)があげられます.Gaussian03のログファイルの読み方を中心に解説していますが,他のソフトを利用する際にも役に立つ本質的な説明がなされています.
名古屋大学情報連携基盤センターニュース Vol.7, No.1-2008.2-

Mulliken電荷とLöwdin電荷

イオン結合性化合物の場合に,Mulliken電荷の方がよい結果を与えることが知られています.一方で,共有結合性化合物ではどちらもあまり変わらないことから,一般的にはMulliken電荷が用いられています.

とはいえ,先ほど述べたように電荷とは観測できない概念的なものですので,興味のある化合物群で傾向を見て判断するのが重要です.

それではpsi4とRDKitを用いて計算を行っていきましょう.まずは必要なライブラリーをimportしておきます.またSMILESを用いてRDKitのMolオブジェクトを生成し,xyz形式へと書き出す関数smi2xyzを定義しておきます.

## ライブラリのインポート
import psi4
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

import time

from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Draw, PandasTools
from rdkit.Chem.Draw import IPythonConsole

print('rdkit version: ', rdBase.rdkitVersion) # rdkit version:  2019.09.3
print('psi4 version: ', psi4.__version__) # psi4 version:  1.3.2
print('numpy version: ', np.__version__) # numpy version:  1.17.4

## SMILESをxyz形式に変換
def smi2xyz(smiles):
    mol = Chem.AddHs(Chem.MolFromSmiles(smiles))
    AllChem.EmbedMolecule(mol, AllChem.ETKDGv2())
    AllChem.UFFOptimizeMolecule(mol)
    conf = mol.GetConformer(-1)
    
    xyz = '0 1'
    for atom, (x,y,z) in zip(mol.GetAtoms(), conf.GetPositions()):
        xyz += '\n'
        xyz += '{}\t{}\t{}\t{}'.format(atom.GetSymbol(), x, y, z)
        
    return xyz

## psi4の設定:各自の環境に合わせてください
psi4.set_num_threads(4)
psi4.set_memory('10GB')

Mulliken電荷

psi4.oeprop(wfn, ‘MULLIKEN_CHARGES’)
wfn.atomic_point_charges()

エネルギー計算や構造最適化計算における戻り値は通常エネルギーのみですが,return_wfnオプションTrueにすることでWavefunctionオブジェクトが得られます.

波動関数の解析にはpsi4.oeprop(One-Electron properties)を用います.解析を行うと,Wavefunctionオブジェクトの所定のメソッドによって計算結果がアクセス可能になります.

Mulliken電荷の問題点

Mulliken密度解析では基底関数系にSTO-3Gなどの最小基底系を利用することを想定しています.そのため,大きな基底関数系を用いた場合には正しくない電荷分布が得られることがあります.このような基底関数系による一貫性の無さがMulliken電荷の問題点として広く認識されています.

以下のコードでは塩化メチルについて,基底関数系の大きさを変えていくことでMulliken電荷がどのように変化していくかを見ていきます.

MeCl = psi4.geometry(smi2xyz('CCl'))
basis_set = ['sto-3g', '3-21G', '6-31G(d)', '6-31+G(d,p)', '6-311++G(2d,p)']

import pandas as pd
df = pd.DataFrame()

psi4.set_output_file('MeCl.txt')
psi4.optimize('hf/sto-3g', molecule=MeCl)
for basis in basis_set:
    _, wfn_cl = psi4.energy('hf/' + basis, molecule=MeCl, return_wfn=True)
    psi4.oeprop(wfn_cl, 'MULLIKEN_CHARGES')
    df['hf/'+basis] = np.array(wfn_cl.atomic_point_charges())

df.index = ['C', 'Cl', 'H', 'H', 'H']

結果は以下の表のようになりました.電気陰性度から考えると炭素の方が塩素よりもポジティブになりそうですが,多くの場合に塩素原子の方がポジティブになっていることがわかります.

hf/sto-3g hf/3-21G hf/6-31G(d) hf/6-31+G(d,p) hf/6-311++G(2d,p)
C -0.13547 -0.75487 -0.53096 -0.43272 -0.40001
Cl -0.17187 -0.03742 -0.11365 -0.09194 -0.05562
H 0.10245 0.26410 0.21487 0.17489 0.15188
H 0.10245 0.26410 0.21487 0.17489 0.15188
H 0.10245 0.26410 0.21487 0.17489 0.15188

Löwdin電荷

psi4.oeprop(wfn, ‘LOWDIN_CHARGES’)
wfn.atomic_point_charges()

Löwdin電荷についても同様にpsi4.oepropを用いて求めることができます.

Mulliken電荷をRDKitの類似度マップで可視化

Draw.SimilarityMaps.GetSimilarityMapFromWeights(mol, weights)

RDKitの類似度マップを用いて原子ごとの寄与を可視化する」という記事でGasteiger電荷について行った電荷マップの可視化をMulliken電荷でも行ってみます.

分子にはアスピリンを用いることとし,HF/STO-3Gで計算します.

## 分子の構築
psi4.set_output_file('aspirin.txt')
aspirin = psi4.geometry(smi2xyz('CC(OC1=C(C(=O)O)C=CC=C1)=O CC(=O)OC1C=CC=CC=1C(O)=O'))
_, wfn_asp = psi4.optimize('hf/sto-3g', molecule=aspirin, return_wfn=True)
rdkit_aspirin = Chem.AddHs(Chem.MolFromSmiles('CC(OC1=C(C(=O)O)C=CC=C1)=O CC(=O)OC1C=CC=CC=1C(O)=O'))
## Mulliken電荷の計算
psi4.oeprop(wfn_asp, 'MULLIKEN_CHARGES')
mulliken = np.array(wfn_asp.atomic_point_charges())
## 類似度マップで可視化
from rdkit.Chem.Draw import SimilarityMaps
fig = SimilarityMaps.GetSimilarityMapFromWeights(rdkit_aspirin, mulliken, colorMap='RdBu')

青い部分が正電荷を帯びている部分になります.

Aspirin

双極子モーメント

psi4.oeprop(wfn, ‘DIPOLE’)
psi4.variable(‘SCF DIPOLE X’)
psi4.variable(‘SCF DIPOLE Y’)
psi4.variable(‘SCF DIPOLE Z’)

双極子モーメントの各座標成分はpsi4.variableで取得可能です.分子全体の双極子モーメントの大きさを求めたい場合は自分でベクトルの大きさを計算する必要があります.

$$ dipole = \sqrt{x^2 + y^2 + z^2} $$

単位はデバイ(D, debye)です.

先ほどのアスピリンの双極子モーメントは1.41Dでした.

dipole_x, dipole_y, dipole_z = psi4.variable('SCF DIPOLE X'), psi4.variable('SCF DIPOLE Y'), psi4.variable('SCF DIPOLE Z')
print('{:.3f}'.format(np.sqrt(dipole_x ** 2 + dipole_y ** 2 + dipole_z ** 2)))
## 1.410

Wiberg結合次数

psi4.oeprop(wfn, ‘WIBERG_LOWDIN_INDICES’)

よく構造化学の分野の論文で使われているWiberg結合次数も求めることが可能です.ログファイル中に結合次数行列として出力されます.

以下の例ではスチレンを用いて,通常の二重結合と芳香族結合の結合次数が異なることを確かめてみましょう.

styrene = psi4.geometry(smi2xyz('C=Cc1ccccc1'))
psi4.set_output_file('Sty_Wiberg.txt')
_, wfn_sty = psi4.optimize('hf/sto-3g', molecule=styrene, return_wfn=True)
psi4.oeprop(wfn_sty, 'WIBERG_LOWDIN_INDICES')

出力は以下のように16原子 x 16原子の行列になります.下のログファイルで太字にしたように,

  • 二重結合:1.95
  • 芳香族結合:1.40

と違いが出ています.

Wiberg Bond Indices using Orthogonal Lowdin Orbitals:

  Irrep: 1 Size: 16 x 16

              1                   2                   3                   4                   5
    1     0.00000000000     1.95324109401     0.00073818541     0.03136594674     0.00033188156
    2     1.95324109401     0.00000000000     1.04581620233     0.00159509968     0.00630944500
    3     0.00073818541     1.04581620233     0.00000000000     1.40797156479     0.00079115441
    4     0.03136594674     0.00159509968     1.40797156479     0.00000000000     1.45142188611
    5     0.00033188156     0.00630944500     0.00079115441     1.45142188611     0.00000000000
    6     0.01628633247     0.00025892466     0.11073590864     0.00089592724     1.44069290921
    7     0.00020670447     0.00679507725     0.00085229099     0.11094786487     0.00093765801
    8     0.02894944073     0.00075567731     1.40471817351059     0.00090681862     0.11142786162
    9     0.98120867928     0.00023437704     0.00410605550     0.00011580928     0.00002081635
   10     0.98394387746     0.00036838891     0.00987370036     0.00034254791     0.00014266115
   11     0.00023335137     0.97589819640     0.00045570928     0.00167940083     0.00005966414
   12     0.00007998189     0.00268424198     0.00007474515     0.97872144722     0.00003084802
   13     0.00008273165     0.00049303434     0.00641896335     0.00002751182     0.98002610760
   14     0.00001230556     0.00008864824     0.00022518555     0.00636499688     0.00003852454
   15     0.00001742722     0.00050475695     0.00639030217     0.00022730512     0.00634128375
   16     0.00029904504     0.00283767782     0.00006327527     0.00634345703     0.00024429452
・・・省略・・・

終わりに

今回は「計算化学における電荷:psi4を用いた電子密度解析」という話題について,

  • 電子密度解析とは何か
  • Mulliken電荷とLöwdin電荷
  • psi4を用いた密度解析のやりかた
  • RDKitの類似度マップを用いた可視化
  • 双極子モーメントの求め方
  • Wiberg結合次数の求め方

などについて説明していきました.

今回取り扱った以外にもpsi4.Wavefunctionオブジェクトからは,HOMOやLUMOといった分子軌道に関する情報も取得可能で,外部プログラムと組み合わせることで軌道の可視化も可能になります.こちらも重要なトピックですので,いずれ扱いたいと思います.

これまでエネルギー計算構造最適化計算について説明を重ねてきました.次回は,調和振動子近似の説明からpsi4における振動数計算について説明することで計算化学について理解を深めていきたいと思います.

コメント

  1. 化学系メーカ勤務者 より:

    たいへん勉強になります。
    ありがとうございます。

    一つ教えて頂きたいことがございます。
    psi4 を用いて、下記のような励起準位の情報取得は可能でしょうか?
    (psi4のHPと睨めっこしていますが、量子計算初心者ということもあり行き詰まっております)
    もし可能であれば、参考になるサイト等あれば教えて頂ければ幸いです。

    ####
    Excited State 2: Singlet-A 3.8723 eV 320.19 nm f=0.1017 =0.000
    109 ->111 0.66378
    110 ->113 -0.20418
    ####

    • 管理人 より:

      ご質問ありがとうございます.

      GaussianのCISまたはTD-DFTのアウトプットファイルからの抜粋とお見受けしました.

      TD-DFTについては近々リリース予定(?)のpsi4 version 1.4で利用可能になるようです.ChemRxivに上がっているプレプリントのFig6にも簡単な例が載っていまして,
      手許の環境ではベンゼンを対象としてこちらのコードは動かすことができました.

      しかし汎関数をCAM-B3LYPに変えるとTDに対応していない旨のエラーが出るなど,何が実装されているか細かくコードレベルで見ていかないとわからないのが現状です.第1,第2励起エネルギーなども取得できましたが,今の段階ではマニュアルが揃っておらず学習コストが高そうです.

      他の(無償)ソフトですと,やはりFireflyかGamessを用いるのがよいと思います.日本語ですと例えば以下の記事が参考になると思います.

      – 「紫外可視吸収スペクトルを計算してみよう」:こちらの記事では,TD-DFTを用いたUV/Visスペクトルの描画方法が丁寧に解説されています
      – 「ホルムアルデヒドの励起状態計算」:こちらの記事ではホルムアルデヒドを例にCIS,TD-DFT,EOM-CCSDを用いた垂直遷移と,CISを用いた励起状態での構造最適化の解説がなされています

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

      • 化学系メーカ勤務者 より:

        早速のご回答ありがとうございます。

        TD-DFT の抜粋でございます。
        どこを探しても説明が見当たらなかったので、TD-DFTがまだ実装されていないとのこと納得しました。
        Gamess を用いるか、psi4 1.4 のリリースを待つか少し検討してみます。

        理論計算、インフォオマティクス関連の記事楽しみにしております。
        今後ともよろしくお願いいたします。

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