Basis Set Exchangeは基底関数系に関するデータベース:PythonとPsi4を用いてECPを自在に扱う

01_計算化学

計算化学では分子軌道を,構成原子の原子軌道の線形結合で表現します.これまでさまざまな原子軌道を表す関数が提唱されています.このような原子軌道を表す関数を基底関数(basis function)と呼び,その集まりを基底関数系(basis set)といいます.

基底関数系は量子化学計算の精度を決める大切な構成要素であり,精度の改善を目指して多くの人によって色々な形のものが報告されています.このような今まで報告されている基底関系を一カ所に集めたデータベースが「Basis Set Exchange(BSE)」になります.

2007年の公開以来,ソフトを用いて計算化学を利用する多くのエンドユーザーに支持されるばかりでなく,理論開発研究者の間でも論文の著者自身が新たなデータを配布するプラットフォームとして用いられています.基底関数系の入手先を一元化することで,異なるユーザー・ソフトウェア間においても計算の再現性を担保しやすくなりました

この記事では2019年になって一新された新生BSEについて,WEBインターフェースやPython APIの使い方について触れ,実際にBSEから取得した基底関数系をPsi4を用いた量子化学計算に適用してみます.

新しくなったBasis Set Exchangeの技術的な詳細については以下の論文を参照してみてください.
New Basis Set Exchange: An Open, Up-to-Date Resource for the Molecular Sciences CommunityJ. Chem. Inf. Model. 2019, 59, 4814-4820.
今回の記事でも量子化学計算の計算エンジンとしてPsi4を用いています。本ブログのPsi4関係の記事は「Psi4の使い方」という記事にまとめています。参照してみてください。
【まとめ】Psi4の使い方
本ブログではPythonのオープンソースライブラリであるPsi4を用いて,量子化学計算の初心者を対象にした記事を多数用意しています。 記事の内容を一つずつ習得することで,これから計算化学を始めたい人が知識・経験ゼロの状態から徐々にステップアップしていくことで,量子化学計算の基本コ...

Basis Set Exchangeとは

Basis Set Exchangeは基底関数に関する包括的なデータベースです.下のようにWEBサイト上で,

  • 左側:基底関数系の一覧から選択
  • 右側:周期表から興味のある元素を選択

どちらかを選択することで,総数600近くから絞り込んでいきます.

BSE

また選択した基底関数系は,各種計算化学用ソフトウェアに対応した形式にてダウンロードが可能になっています.

BSE format

ここではヨウ素を含む化合物のハロゲン結合の計算でたびたび用いられるLanL2DZ(d,p)をダウンロードしてみます.こちらはGaussian16やpsi4にはデフォルトでは組み込まれていませんので,BSEから入手する必要があります.

今回は「psi4」形式で試してみます.

  • LANL2DZdp
  • ヨウ素
  • psi4

が選択されていることを確認します.

BSE lanl2dzdp

Get Basis Set」をクリックすることでテキストファル形式のウィンドウが現れます.「Download」ボタンを押すことで望みの基底関数系が入手できました.

Get bs

Basis Set Exchangeをpythonから使う

新しくなったBSEではpythonから操作可能なライブラリが用意されており,pipでインストール可能です.

pip install basis_set_exchange

コマンドラインから,またはpythonスクリプト中で機械的に任意の基底関数系を取得可能です.

この記事ではpythonスクリプト中での使い方をみてみます.

pythonライブラリの利用法

basis_set_exchange.get_basis(basis-set, elements, fmt)

pythonライブラリの「basis_set_exchange」でもget_basis関数を用いて,

  • どの基底関数系を
  • どの原子に対して(原子番号で指定)
  • どのフォーマットで

ダウンロードするかを指定します.

先ほどと同じく,ヨウ素用のLanL2DZ(d,p)をpsi4用フォーマットで入手します.

import basis_set_exchange as bse
lanl2dzdp = bse.get_basis('lanl2dzdp', elements=[53], fmt='psi4')
print(lanl2dzdp)

得られたデータはpythonの文字列オブジェクトになります.このオブジェクトを後ほど,psi4を用いた独自の基底関数系設定のところで利用していきます.

spherical

!----------------------------------------------------------------------
! Basis Set Exchange
! Version v0.8.12
! https://www.basissetexchange.org
!----------------------------------------------------------------------
!   Basis set: LANL2DZdp
! Description: DZP  Double Zeta + Polarization + Diffuse ECP
!        Role: orbital
!     Version: 0  (Data from the original Basis Set Exchange)
!----------------------------------------------------------------------

...略

****

I     0
I-ECP     3     46
f potential
  5
0      1.0715702             -0.0747621
1     44.1936028            -30.0811224
2     12.9367609            -75.3722721
2      3.1956412            -22.0563758
2      0.8589806             -1.6979585
s-f potential
  5
0    127.9202670              2.9380036
1     78.6211465             41.2471267
2     36.5146237            287.8680095
2      9.9065681            114.3758506
2      1.9420086             37.6547714
p-f potential
  5
0     13.0035304              2.2222630
1     76.0331404             39.4090831

...略

Psi4で基底関数系の設定方法

psi4.set_options({‘basis’: basis-set})

計算化学における電荷:psi4を用いた電子密度解析」などの記事ではPsi4を用いた計算では例えば,

psi4.energy('mp2/sto-3g', molecule=mol)

のようにエネルギー計算を実行する関数の中で基底関数系を指定していました.計算で用いる基底関数系はbasisという(グローバル)オプション変数で指定できますので,上記コードは以下のものと同じです.

psi4.set_options({'basis': 'sto-3g'})
psi4.energy('mp2', molecule=mol)

それでは原子の種類ごとに異なる基底関数系を用いて計算を行いたい場合はどうすればよいでしょうか?

basisとして独自に定義した基底関数系を用いればよいです.

基底関数系を独自に定義

psi4.basis_helper(block, name, set_option=True)

原子毎に異なる基底関数系を用いて計算する際には,「どの原子にどの基底関数系を用いるか」を以下のようなフォーマットで記していきます.

assign 原子(省略可) basis-set

PsiAPI以前から使われている記法をそのまま移植しているためpythonicではないフォーマットですが,この形式のブロックを連ねてstringオブジェクトとして作成します.

このようにして得た文字列をbasis_helper関数に渡すことで,作成した独自の定義をbasisオプションとして設定できるようになります.

この説明ではややわかりにくいと思いますので,以下でメタンを用いてコード例を示します.

ここではmybs1mybs2という2つを定義しています.

  • mybs1:全ての原子にSTO-3Gを設定(通常のSTO-3Gと実質同じ)
  • mybs2:炭素に3-21Gを,残りの原子にSTO-3Gを設定

bs2では最初全ての原子にSTO-3Gを割り当てていますが,そのあとで炭素を3-21Gに上書きしています.

basis_helperではset_optionFalseにしない場合は,basisオプションに作成した定義をそのまま設定しますので,すぐに利用可能です.

以下のコードでは確認と学習のためにpsi4.core.BasisSet.buildを用いて,ログファイルに独自の定義の詳細を描き込んでいます.

import psi4
psi4.set_output_file('CH4_bs.txt')

mol1 = '''
0 1
 C                  0.12632979   -0.69148935    0.00000000
 H                  0.48298421   -1.70029935    0.00000000
 H                  0.48300263   -0.18709116    0.87365150
 H                  0.48300263   -0.18709116   -0.87365150
 H                 -0.94367021   -0.69147617    0.00000000
 '''

mybs1 = '''
assign sto-3g
'''

mybs2 = '''
assign sto-3g
assign C 3-21G
'''

psi4mol = psi4.geometry(mol1)
psi4.basis_helper(mybs1, name='mybs1')
f1 = psi4.core.BasisSet.build(psi4mol)
f1.print_detail_out()

psi4.basis_helper(mybs2, name='mybs2')
f2 = psi4.core.BasisSet.build(psi4mol)
f2.print_detail_out()

ログファイルの内容を一部抜粋してみます.

    Name: MYBS1
    Role: ORBITAL
    Keyword: None
    atoms 1   entry C          line    61 file /Users/.../psi4/basis/sto-3g.gbs 
    atoms 2-5 entry H          line    19 file /Users/.../psi4/basis/sto-3g.gbs 
...
    Name: MYBS2
    Role: ORBITAL
    Keyword: None
    atoms 1   entry C          line    68 file /Users/.../psi4/basis/3-21g.gbs 
    atoms 2-5 entry H          line    19 file /Users/.../psi4/basis/sto-3g.gbs 

さらに分子の入力座標を工夫することで,同じ元素の異なる原子毎に異なる基底関数系を指定することも可能です.

以下の例ではメタンの水素を3種類(H, H2, H3)にわけ,それぞれに

  • H:sto-3g
  • H2:cc-pVDZ
  • H3:pp-pVTZ

のように異なる基底関数系を割り当てたmybs3を定義しています.

mol2 = '''
0 1
 C                  0.12632979   -0.69148935    0.00000000
 H                  0.48298421   -1.70029935    0.00000000
 H2                  0.48300263   -0.18709116    0.87365150
 H2                  0.48300263   -0.18709116   -0.87365150
 H3                 -0.94367021   -0.69147617    0.00000000
'''

mybs3 = '''
assign sto-3g
assign H2 cc-pVDZ
assign H3 cc-pVTZ
'''

psi4mol2 = psi4.geometry(mol2)
psi4.basis_helper(mybs3, name='mybs3')
f3 = psi4.core.BasisSet.build(psi4mol2)
f3.print_detail_out()

同様にログファイルを抜粋します.

    Name: MYBS3
    Role: ORBITAL
    Keyword: None
    atoms 1   entry C          line    61 file /Users/.../psi4/basis/sto-3g.gbs 
    atoms 2   entry H          line    19 file /Users/.../psi4/basis/sto-3g.gbs 
    atoms 3-4 entry H          line    22 file /Users/.../psi4/basis/cc-pvdz.gbs 
    atoms 5   entry H          line    23 file /Users/.../psi4/basis/cc-pvtz.gbs 

もっともこの例のように分子内の異なる水素にそれぞれ違った基底関数をわりあてて計算することはまずありません.こういった指定法が最も使われるのが,次節で説明する「ECP(有効内殻ポテンシャル)」の指定です.

psi4を用いたECPの設定

注意)psi4でのECPの利用は未だβ版の扱いです.

ここではヨウ化メチルと塩化物イオンからなる系の計算を考えてみます.その際,

  • 水素・炭素・塩素
  • ヨウ素

の2種類について別々の基底関数系を割り当てていき,ヨウ素についてはECPを用いることとします.ECPとしては元々psi4で利用可能な「LanL2DZ」に加えて,先ほどBSEから取得した「LanL2DZ(d,p)」を用いることで以下のように5つの形式を定義します.

名前 水素・炭素・塩素 ヨウ素
bs1 aug-cc-pVDZ LanL2DZ
bs2 aug-cc-pVTZ LanL2DZ
bs3 aug-cc-pVDZ LanL2DZ(d,p)
bs4 aug-cc-pVTZ LanL2DZ(d,p)
bs5 aug-cc-pVQZ LanL2DZ(d,p)

BSEからダウンロードしたデータを用いて,指定の原子に割り当てるには以下のように[original-bs]を定義して,それをassignしていきます.「中身」と書いてある部分にダウンロードした文字列オブジェクトを用います.

assign 原子 original-bs
[original-bs]
中身・・・

具体例を見た方がわかりやすいと思います.

以下では[ecp]という名前で定義しています.basis_helper関数set_optionFalseにしていますので,この段階ではbasisオプションには何も設定されていません.

import psi4
import basis_set_exchange as bse
lanl2dzdp = bse.get_basis('lanl2dzdp', elements=[53], fmt='psi4')

bs1 = '''
assign aug-cc-pVDZ
assign I lanl2dz
'''

bs2 = '''
assign aug-cc-pVTZ
assign I lanl2dz
'''

bs3 = '''
assign aug-cc-pVDZ
assign I ecp

[ecp]
{}'''.format(lanl2dzdp)

bs4 = '''
assign aug-cc-pVTZ
assign I ecp

[ecp]
{}
'''.format(lanl2dzdp)

bs5 = '''
assign aug-cc-pVQZ
assign I ecp

[ecp]
{}
'''.format(lanl2dzdp)

psi4.basis_helper(bs1, name='bs1', set_option=False)
psi4.basis_helper(bs2, name='bs2', set_option=False)
psi4.basis_helper(bs3, name='bs3', set_option=False)
psi4.basis_helper(bs4, name='bs4', set_option=False)
psi4.basis_helper(bs5, name='bs5', set_option=False)

Psi4でECPを用いてハロゲン結合のエネルギーを計算

今回の記事ではヨウ素ー塩素距離を変えるごとに構造最適化は行っていません.構造が変化する部分はメチル基程度ですので,結果にはほぼ影響ないと考えますが一応念頭に置いてください.

続いてヨウ化メチルと塩化物イオンを用いて,前項で定義した5つの手法を用いてヨウ素ー塩素距離を変えながら相互作用エネルギーを求めていきます.その際,ハロゲン結合の角度特異性を考えて塩化物イオンは炭素ーヨウ素結合の延長線上に位置することとします

まずはライブラリのインポートと,塩化物イオンの位置を決めるためのC-I結合に平行な単位ベクトルを求めます.

import psi4
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
%matplotlib inline
import time

au2kcal = psi4.constants.hartree2kcalmol
now = time.strftime('%Y%m%d-%H%M')
psi4.set_num_threads(10)
psi4.set_memory('14GB')
psi4.set_output_file(now + '_BSE.log')

geom = '''
0 1
 C                 -2.32047871   -0.58510637    0.00000000
 H                 -1.96382429   -1.59391638    0.00000000
 H                 -1.96380587   -0.08070818    0.87365150
 H                 -3.39047871   -0.58509319    0.00000000
 I                 -1.62046660    0.40483400   -1.71464314
--
-1 1
 Cl                {}            {}            {}
'''

r0 = np.array([-2.32047871, -0.58510637, 0.00000000])
r1 = np.array([-1.62046660, 0.40483400, -1.71464314])
r = (r1 - r0)/np.sqrt(np.sum((r1 - r0) ** 2))

以下のコードで実際にエネルギーを計算しています.なおDFTの汎関数にはwB97X-Dを用いています.

calc_energy関数では,錯体のエネルギーに加えて,「Psi4で分子間相互作用エネルギーの計算:超分子法と基底関数重なり誤差」で解説したCP法を用いた相互作用エネルギーを計算します.set_options({‘basis’: bs})の部分で用いる基底関数系をbasisオプションに設定しています.

この関数を用いて原子間距離を1.5から4.4Åまで0.1Åずつ変えながら,先ほど定義したbs1からbs5を用いてエネルギー計算を行います.

basis_set = ['bs1', 'bs2', 'bs3', 'bs4', 'bs5']

def calc_energy(molecule, bs):
    psi4.set_options({'basis': bs})
    int_E = psi4.energy('wb97x-d', bsse_type='cp', molecule=molecule)
    calc_E = psi4.energy('wb97x-d', molecule=molecule)
    psi4.core.clean()

    return int_E * au2kcal, calc_E


electronic_energy = []
int_energy = []

for i in np.arange(1.5, 4.5, 0.1):
    x, y, z = r * i + r1
    psi4mol = psi4.geometry(geom.format(x, y, z))
    calced_e = []
    calced_i = []
    for bs in basis_set:
        int_e, calc_e = calc_energy(psi4mol, bs)
        calced_i.append(int_e)
        calced_e.append(calc_e)

    electronic_energy.append(calced_e)
    int_energy.append(calced_i)


df1 = pd.DataFrame(electronic_energy, index=[0.1*i for i in range(15,45)])
df2 = pd.DataFrame(int_energy, index=[0.1*i for i in range(15,45)])
df1.columns = ['bs1', 'bs2', 'bs3', 'bs4', 'bs5']
df2.columns = df1.columns

得られた結果の全体図と拡大図が以下のようになります.bs5を用いて最小エネルギーを与える点も併せて示しておきました.この距離(約3.13Å)での相互作用エネルギーは約-7.7 kcal/molでした.

ヨウ素はサイズの大きい原子ですので,近距離(2.5Å以下)では反発が非常に大きいことが見てとれます.また相互作用エネルギーの大きさに関しては,塩化物イオンの基底関数系を大きくする(bs1 vs bs2)よりもヨウ素上のECPを変える(bs1 vs bs3)の方が影響が大きいこともわかります.

Energy plot

おわりに

今回は「Basis Set Exchangeは基底関数系に関するデータベース:PythonとPsi4を用いてECPを自在に扱う」という話題について,

  • Basis Set Exchange(BSE)とはなにか
  • BSEのWEBインターフェースの使い方
  • BSEをpythonスクリプトからの利用方法
  • Psi4で独自の基底関数系を定義する方法
  • ヨウ化メチルと塩化物イオンを例としたエネルギー計算

などについてみてきました.今回の内容を消化することで,分子中の原子について思い通りの方法で基底関数系を割り当てて計算可能になったことと思います.

コメント

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