【PySCF】TD-DFT計算を使ってUVスペクトル可視化【Pythonで始める量子化学計算】

本記事では、Pythonを用いた量子化学計算ツールであるPySCFを使って分子の励起状態を解析する方法を解説します。PySCFを使うことで、分子の光吸収や蛍光物性の予測が可能になります。具体的には、基底状態の構造最適化から、TD-DFT計算による励起状態の計算、UVスペクトルの可視化までを扱います。この記事を読むことで、簡単に対象分子の励起状態とUVスペクトルの結果を得ることができ、蛍光分子の設計にも役立てることができます。

Gaussianを使った量子化学計算の初心者向け技術書を販売中

しばしば出くわすエラーへの対処法をはじめ
Gaussianと無料ソフトウェア Avogadro を組み合わせた物性解析手法が学べます！

技術書ページへ

1. 量子化学計算用のソフトウェア: PySCFとは

参考文献[1]から引用

PySCF（Python Simulations of Chemistry Framework）は、Pythonベースのオープンソースの量子化学計算フレームワークです。PySCFは、高度な電子構造計算をシンプルかつ効率的に実行するためのツールを提供し、研究者や学生が量子化学の計算を容易に行えるように設計されています。

PySCFの主な特徴は以下の通りです：

• モジュラー設計：PySCFは、計算の各ステップ（分子の構造定義、基底状態の計算、励起状態の計算など）をモジュールとして提供しており、必要な部分だけを組み合わせて使用することができます。
• 高性能計算：内部に高度なアルゴリズムを使用しており、大規模な計算でも高速に処理を行うことができます。
• 柔軟性：Pythonで記述されているため、カスタマイズや他のPythonライブラリとの連携が容易です。

本記事では、PySCFを用いたTD-DFT（時間依存密度汎関数理論）計算を通じて、UVスペクトルを可視化する方法について詳しく解説します。

2. 吸光と発光と分子構造の関係

吸光と発光は、分子の電子遷移に関する現象であり、化学や物理学の研究において重要な役割を果たします。

• 吸光：分子が光子を吸収し、電子が低いエネルギー準位から高いエネルギー準位へと遷移する現象です。この現象はUV-Vis（紫外-可視）スペクトルとして観測されます。
• 発光：励起された電子が元の低エネルギー準位に戻る際に光子を放出する現象です。これにより、蛍光やリン光が観測されます。

分子構造はこれらの現象に大きな影響を与えます。特に、電子の配置や分子軌道のエネルギー準位は、吸光・発光の波長や強度に直接関係します。分子構造の変化により、吸収スペクトルや発光スペクトルが変化するため、分子構造を最適化し、その結果として得られるスペクトルを解析することは非常に重要です。

3. UVスペクトルとは

UV（紫外線）スペクトルは、分子が紫外線領域の光を吸収する際の吸収強度を波長（またはエネルギー）の関数として表したものです。UVスペクトルは、分子の電子構造に関する情報を提供し、化学物質の特性評価や分析に広く利用されます。

励起状態計算の一点計算はUVスペクトルを求めることができる

励起状態計算を行うことで、分子の電子遷移エネルギーを求めることができます。このエネルギー情報を基に、吸収波長を計算し、UVスペクトルを作成することが可能です。具体的には、TD-DFT計算を使用して、各遷移のエネルギーと振動子強度（吸収強度）を求め、それを基にスペクトルを描画します。

UVスペクトルは基底状態の構造最適化構造での垂直遷移

UVスペクトルは、基底状態の最適化された構造における垂直遷移に基づいて計算されます。垂直遷移とは、基底状態の構造から電子が励起される際に、分子構造が変化しないと仮定する遷移のことです。この仮定により、計算の複雑さが軽減され、実用的な時間で計算を行うことができます[2]。

基底・励起状態ポテンシャル曲線と電子遷移の関係

4. TD-DFTとは？

TD-DFT（Time-Dependent Density Functional Theory、時間依存密度汎関数理論）は、分子の励起状態の性質を研究するための計算手法です。通常の密度汎関数理論（DFT）は基底状態のエネルギーと電子分布を求めるのに対し、TD-DFTは時間依存する系の応答を解析し、分子の光学的性質を評価するために使用されます。

TD-DFTの基本概念

TD-DFTは、電子密度の時間発展を追跡することで、分子の励起状態を解析します。これにより、電子遷移エネルギー（吸収波長）や振動子強度（吸収強度）を計算し、UV-Visスペクトルを予測することができます。TD-DFTは、以下のような特徴を持っています：

• 効率性：従来の励起状態計算手法（例えば、配置間相互作用法や多配置自己無撞着場法）に比べて計算コストが低く、大規模な分子系に対しても適用可能です。
• 精度：適切な交換相関汎関数を選択することで、高い精度で励起状態のエネルギーを予測できます。

TD-DFTの応用

TD-DFTは、分子の光学的性質や電子遷移の解析に広く利用されています。具体的には、UV-Visスペクトルの予測、蛍光・リン光の解析、光化学反応の研究などに応用されています。

TD-DFTの理論背景

TD-DFTは、電子密度汎関数理論（DFT）を拡張したものであり、時間依存する外部場に対する電子系の応答を記述します。Kohn-Sham軌道を用いて、基底状態から励起状態への遷移を記述し、これを時間依存するSchrödinger方程式に基づいて解きます。TD-DFTは、擬ポテンシャル法や長距離補正法などを組み合わせることで、高い精度での計算を可能にしています。

5. TD-DFTを行う時の注意点

TD-DFT計算を行う際には、いくつかの注意点があります。これらのポイントを押さえることで、より信頼性の高い計算結果を得ることができます[2]。

5.1 長距離補正を入れること

TD-DFT計算では、交換相関汎関数の選択が重要です。特に、長距離補正（LC、Long-Range Correction）を含む汎関数を使用することが推奨されます。長距離補正を含む汎関数は、電子間の長距離相互作用を適切に取り扱うことができ、特に電荷移動（CT）状態やπ-π*遷移などの長距離励起状態のエネルギーを正確に予測します。

例えば、CAM-B3LYP、ωB97XD、LC-ωHPBE、M06-2Xなどの汎関数は、長距離補正を含んでおり、これらの汎関数を用いることで、励起状態のエネルギー計算の精度が向上します。

5.2 基底関数に分散関数を加えること

TD-DFT計算に使用する基底関数セットも重要です。分散関数（Diffuse Functions）を含む基底関数セットを選択することで、電子が広がった分子系（例えば、共役系や大きな分子）に対する計算精度を向上させることができます。分散関数は、電子密度の広がりを正確に再現し、励起状態のエネルギー計算における精度を高めます。

5.3 励起状態の計算は基底状態の計算と同程度の精度で行うこと

励起状態計算を行う前に、基底状態の構造最適化を高精度で行うことが重要です。基底状態の計算が不正確である場合、励起状態の計算結果も不正確になる可能性が高いため、基底状態の計算と同程度の精度で励起状態の計算を行うことが求められます。

5.4 HOMO-LUMOギャップと第一励起状態の違い

HOMO-LUMOギャップは、分子のHOMO（最高被占軌道）とLUMO（最低空軌道）間のエネルギー差を示します。しかし、これは必ずしも第一励起状態のエネルギーに対応するわけではありません。第一励起状態は、実際には多くの場合、HOMO-LUMOギャップよりも高いエネルギーを持つことがあり、異なる遷移が関与することもあります。そのため、TD-DFT計算を行う際には、HOMO-LUMOギャップだけではなく、具体的な遷移エネルギーを詳しく解析することが重要です。

5.5 小さいHOMO-LUMOギャップの計算にはCASSCFが必要

小さいHOMO-LUMOギャップや多重励起状態の解析には、TD-DFTではなくCASSCF（Complete Active Space Self-Consistent Field）法が必要となります。CASSCFは、活性空間内の電子相互作用を厳密に取り扱うことができ、特に多重励起状態や強い相関を持つ系の計算に適しています。TD-DFTが適用できない場合や精度が不足する場合は、CASSCFを用いることでより正確な計算結果を得ることができます。PyscfでもCASSCFを使って計算することが出来ます。

以上の注意点を踏まえて、次の章では、実際にTD-DFTを用いたUVスペクトルの可視化方法について、具体的なPythonコードを用いて解説します。

6. TD-DFTを用いたUVスペクトル可視化実践

それでは、PySCFを用いてTD-DFT計算を実行し、得られたデータを基にUVスペクトルを可視化する具体的な手順を示します。以下のコードを実行することで、分子の励起状態を解析し、UVスペクトルを生成することができます[3-8]。

分子構造の定義と基底状態の計算

まず、解析対象となる分子の構造を定義し、基底状態の計算を行います。ここでは、1,2-oxazoleを対象にPySCFを用いて分子の基底状態の構造最適化を実行します。構造最適化の詳細はこちらの記事で解説しています。

from pyscf import gto, scf, dft, tddft
from pyscf.geomopt import geometric_solver
import py3Dmol
import matplotlib.pyplot as plt
import numpy as np
from scipy.constants import physical_constants
import pandas as pd
import time
import plotly.express as px

# 分子構造の定義
mol = gto.M(atom="""
  C      1.0701      0.4341     -0.0336
  C      0.8115     -0.9049     -0.1725
  C     -0.6249     -1.0279     -0.0726
  H      1.9842      1.0231     -0.0364
  H      1.5156     -1.7176     -0.3255
  H     -1.2289     -1.9326     -0.1322
  O     -0.0873      1.1351      0.1422
  N     -1.1414      0.1776      0.1122""",
  basis='6-31+g',
  verbose=3)

# SCF計算のセットアップ
mf = dft.RKS(mol)  # 分子オブジェクトをRKSクラスに渡す
mf.chkfile = path + '.chk'  # 計算のチェックポイントファイルの保存先
mf.xc = 'CAMB3LYP'  # 交換相関汎関数としてCAM-B3LYPを使用
mf.max_cycle = 150  # SCF計算の最大反復回数

# コールバック関数と収束パラメータの定義
geometries = []
energies = []

def cb(envs):
    mf = envs['g_scanner'].base  # 現在の計算オブジェクトを取得
    energies.append(mf.e_tot)  # エネルギーをリストに追加
    geometries.append(mol.atom_coords(unit="ANG"))  # ジオメトリーをリストに追加

conv_params = {
    'maxsteps': 50,
    'max_force': 1e-4,
    'max_energy': 1e-6
}

# 構造最適化を実行
opt = geometric_solver.optimize(mf, callback=cb, **conv_params)

# エネルギー収束のプロット
plt.figure(figsize=(10, 5))
plt.plot(energies, marker='o')
plt.title('Energy Convergence During Geometry Optimization')
plt.xlabel('Optimization Step')
plt.ylabel('Total Energy (Hartree)')
plt.grid(True)
plt.show()

TD-DFTの計算

基底状態の最適化が完了したら、次にTD-DFT計算を実行してUVスペクトルを予測します。以下のコードでは、CAM-B3LYP汎関数を使用して励起状態を計算し、得られた遷移エネルギーと振動子強度を基にUVスペクトルを生成します。具体的なステップを説明します。

import numpy as np
from scipy.constants import physical_constants
import pandas as pd
import plotly.graph_objects as go

# ハートリーからeVへの変換係数
ha_2_ev = 1 / physical_constants["electron volt-hartree relationship"][0]

# ガウス関数の定義
def gaussian(x, mu, sig):
    return np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))

まず、必要なライブラリをインポートし、ハートリー単位から電子ボルト（eV）への変換係数を定義します。次に、ガウス関数を定義します。ガウス関数は、UVスペクトルのピークをスムーズにするために使用されます。

# スペクトル解析を行う関数
def run_spectral_analysis(mol, xc="CAMB3LYP"):
    n_states = 15
    spectral_width = 0.1

    # 基底状態DFT計算
    mf = dft.RKS(mol, xc=xc).run()

run_spectral_analysis関数では、まず基底状態のDFT計算を行います。ここでは、CAM-B3LYP汎関数を使用しています。mf = dft.RKS(mol, xc=xc).run()は、Restricted Kohn-Sham（RKS）計算を実行し、基底状態のエネルギーと電子密度を求めます。

    # 励起状態DFT計算
    mytd = tddft.TDDFT(mf)
    mytd.nstates = n_states
    mytd.kernel()
    mytd.analyze()
    osc_strengths = mytd.oscillator_strength()[:n_states-5]

次に、TD-DFT計算を実行します。mytd = tddft.TDDFT(mf)でTD-DFT計算オブジェクトを作成し、mytd.nstates = n_statesで計算する励起状態の数を設定します。mytd.kernel()で計算を実行し、mytd.analyze()で計算結果を解析します。mytd.oscillator_strength()で得られる振動子強度（吸収強度）は、スペクトルのピークの高さを決定します。

    # NaNを0に置き換える
    osc_strengths = np.nan_to_num(osc_strengths)

    # スペクトルを作成
    energies_ev = mytd.e[:n_states-5] * ha_2_ev
    x_range = np.linspace(energies_ev.min() * 0.9, energies_ev.max() * 1.1, num=1000)
    intensity = np.zeros(x_range.size)

得られた振動子強度の中にNaNが含まれている場合、それを0に置き換えます。次に、励起状態のエネルギーをeV単位に変換し、スペクトルのx軸の範囲を設定します。x_rangeは、エネルギーの範囲をカバーする1000ポイントの配列です。intensityは、この範囲内でのスペクトル強度を格納するための配列です。

    for e, f in zip(energies_ev, osc_strengths):
        intensity += gaussian(x_range, e, spectral_width) * f

    # 大まかな正規化
    dx = (x_range[-1] - x_range[0]) / x_range.size
    area = (intensity * dx).sum()
    if area != 0:  # Check to prevent division by zero
        intensity /= area

    return x_range, intensity

各励起エネルギーと対応する振動子強度に対して、ガウス関数を用いてスペクトルを生成します。intensity配列にガウス分布を加えることで、スペクトルの各ピークが形成されます。最後に、スペクトルを正規化して総面積が1になるようにします。

# データ保存用の辞書
data = {"Excitation Energy (eV)": [], "Intensity": [], "Exchange-Correlation Functional": [], "Excitation Energy (nm)": [], "Excitation Energy (cm-1)":[]}

# CAM-B3LYPで計算
xcs = ["CAMB3LYP"]

for xc in xcs:
    x_range, intensity = run_spectral_analysis(mol, xc=xc)
    data["Excitation Energy (eV)"] += x_range.tolist()
    data["Intensity"] += intensity.tolist()
    data["Exchange-Correlation Functional"] += [xc] * x_range.size
    data["Excitation Energy (nm)"] += (1239.841984 / x_range).tolist()
    data["Excitation Energy (cm-1)"] += (8065.54429 / x_range).tolist()

# データフレームに変換
df = pd.DataFrame(data)

# データが正しく保存されているか確認
df.head()

次に、得られたデータを辞書に保存します。計算された励起エネルギー、強度、および対応する波長（nm）と波数（cm-1）を格納します。最後に、これらのデータをPandasのデータフレームに変換し、最初の数行を表示してデータが正しく保存されているか確認します。

UVスペクトルの可視化

最後に、得られたスペクトルデータを可視化します。ここでは、Plotlyを用いてインタラクティブなスペクトルプロットを作成します [9]。

# グラフの作成
fig = go.Figure()

# 各単位のトレースを作成
trace_ev = go.Scatter(x=df["Excitation Energy (eV)"], y=df["Intensity"], mode='lines+markers', name='eV')
trace_nm = go.Scatter(x=df["Excitation Energy (nm)"], y=df["Intensity"], mode='lines+markers', name='nm')
trace_cm = go.Scatter(x=df["Excitation Energy (cm-1)"], y=df["Intensity"], mode='lines+markers', name='cm-1')

# トレースをグラフに追加
fig.add_trace(trace_ev)
fig.add_trace(trace_nm)
fig.add_trace(trace_cm)

# ドロップダウンメニューを作成
fig.update_layout(
    updatemenus=[
        dict(
            active=0,
            buttons=list([
                dict(label="eV", method="update", args=[{"visible": [True, False, False]}, {"title": "Excitation Energy (eV)"}]),
                dict(label="nm", method="update", args=[{"visible": [False, True, False]}, {"title": "Excitation Energy (nm)"}]),
                dict(label="cm-1", method="update", args=[{"visible": [False, False, True]}, {"title": "Excitation Energy (cm-1)"}]),
            ]),
        )
    ]
)

# 初期表示の設定
fig.update_traces(visible=False)
fig.data[0].visible = True

# レンジスライダーを追加
fig.update_layout(
    title="Spectral Analysis",
    xaxis_title="Excitation Energy",
    yaxis_title="Intensity",
    xaxis=dict(
        rangeslider=dict(
            visible=True
        ),
        type="linear"
    )
)

fig.show()

結果は以下のようなUVスペクトルを得ることが出来ます。横軸の単位の変更は左のプルダウンで行うことが出来ます。ガウス関数の幅は適度に調整してください。

終わりに

いかがでしたでしょうか。PySCFはGoogle Colab上でTD-DFTの計算、UVスペクトルの可視化ができるため、非常に手軽に計算を行うことができます。さらに、Google Colab有料版などを使うと自宅PCではできないような大きな分子も計算可能です。是非気軽に試してみてください！

参考文献

Gaussianを使った量子化学計算の初心者向け技術書を販売中

しばしば出くわすエラーへの対処法をはじめ
Gaussianと無料ソフトウェア Avogadro を組み合わせた物性解析手法が学べます！

技術書ページへ