PySCFを用いた溶媒モデルの計算を解説！ 気相だけでなく溶媒を考慮したエネルギー計算を行い、実験データと比較する方法を紹介します。PCMの基本、PySCFの実装手順、構造最適化や振動数解析の応用までを網羅。様々な溶媒モデルの比較や実験データの解析に役立つ内容を詳しく解説しています。
Google colab, pyscf 2.7.0
Gaussianを使った量子化学計算の初心者向け技術書を販売中
しばしば出くわすエラーへの対処法をはじめ
Gaussianと無料ソフトウェア Avogadro を組み合わせた物性解析手法が学べます！
量子化学計算用のソフトウェア: PySCFとは
参考文献[1]から引用
PySCF（Python Simulations of Chemistry Framework）は、Pythonベースのオープンソースの量子化学計算フレームワークです。PySCFは、高度な電子構造計算をシンプルかつ効率的に実行するためのツールを提供し、研究者や学生が量子化学の計算を容易に行えるように設計されています。
PySCFの主な特徴は以下の通りです：
- モジュラー設計：PySCFは、計算の各ステップ（分子の構造定義、基底状態の計算、励起状態の計算など）をモジュールとして提供しており、必要な部分だけを組み合わせて使用することができます。
- 高性能計算：内部に高度なアルゴリズムを使用しており、大規模な計算でも高速に処理を行うことができます。
- 柔軟性: Pythonベースであるため、スクリプトを通じて簡単に操作や拡張が可能です。ユーザーは、計算手法やパラメータをカスタマイズして特定のニーズに合わせることができます。
- 拡張性: オープンソースであるため、研究者や開発者が自分の研究に合わせてPySCFを拡張することが容易です。新しいアルゴリズムや方法論を追加することで、独自の研究を進めることができます。また、拡張機能が開発され続けている事も特徴です。
本記事では、PySCFを用いて溶媒効果を取り入れた量子化学計算を実行し、その結果を詳細に解析します。
1. PCMとは？
1.1 溶媒効果の取り扱い
化学計算において、分子はしばしば溶媒中で存在し、その影響を受けます。溶媒効果は、分子のエネルギー状態、振動数、電子構造、励起状態などに影響を与えます。
しかし、分子軌道法や密度汎関数法（DFT）などの量子化学計算は通常気相を仮定しており、溶媒効果を考慮しないことが一般的です。
溶媒効果を計算に取り入れる方法として、主に以下の2種類があります：
- 連続溶媒モデル（Continuum Solvent Model）
- 溶媒を誘電率を持つ連続的な媒質として扱う。
- 個々の溶媒分子を考慮せず、平均的な溶媒効果を計算する。
- 計算コストが低く、大規模系にも適用可能。
- 代表的な手法：PCM（Polarizable Continuum Model）
- 明示的溶媒モデル（Explicit Solvent Model）
- 溶媒分子を個々に考慮し、分子間相互作用を直接計算する。
- 計算コストが高く、特に分子動力学（MD）シミュレーションと併用される。
- 代表的な手法：QM/MM（量子力学/分子力学）、クラスターモデル
ここでは、連続溶媒モデルの代表的な手法であるPCM（Polarizable Continuum Model） に焦点を当て、PySCFでの実装方法を説明します。
1.2 PCMの概要[3]
Polarizable Continuum Model（PCM） は、分子を誘電率を持つ連続溶媒中に埋め込むことで、溶媒の極性による分子のエネルギー変化を取り入れる方法です。
PCMは溶媒を誘電率によって特徴づけ、分子の周囲に仮想的な電荷分布（反応場）を作ることで、溶媒和の影響を考慮します。
1.3 PySCFで利用可能なPCMモデル
PySCFでは、以下のPCM手法が利用可能です：
- C-PCM（Conducting PCM）
- 最も一般的なPCMモデルで、溶媒を理想導体と仮定する。
- 高い誘電率の溶媒（水など）に適している。
- IEF-PCM（Integral Equation Formalism PCM）
- C-PCMの改良版で、誘電率の影響をより厳密に取り扱う。
- GaussianのPCMではデフォルト設定
- ほぼすべての誘電率の溶媒に適用可能。
- SS(V)PE（Surface and Volume Polarization for Electrostatics）
- 分子の形状に基づいて反応場を決定する方法。
- 高精度な計算が必要な場合に適用。
- COSMO（COnductor-like Screening MOdel）
- C-PCMの近似モデルであり、計算コストを削減しながら溶媒効果を再現する。
2. PySCFにおけるPCMの実装
PySCFでは、簡単にPCMを適用することが可能です。
以下のコードでは、アセトン（C₃H₆O）のSCF計算にPCMを適用し、水中でのエネルギーを計算します。
コード①：PCMを適用したSCF計算
# 必要なライブラリのインストール
!pip install pyscf numpy from pyscf import gto, scf, solvent
# 分子定義（アセトン）
mol = gto.M(atom="""
C 0.0000 0.0000 0.0000
O 1.2000 0.0000 0.0000
H -0.6000 0.9000 0.0000
H -0.6000 -0.9000 0.0000
""", basis='6-31G(d)')
# SCF計算 + PCM（溶媒：水）
mf = scf.RHF(mol).PCM()
mf.with_solvent.method = 'IEF-PCM'
# 使用するPCMモデル（IEF-PCMを指定）
mf.with_solvent.eps = 78.3553 # 水の誘電率
mf.kernel()
3. 励起状態とPCM
PCMは、励起状態のエネルギー計算（TDDFT）にも適用可能です。
励起状態では、溶媒の影響によって遷移エネルギーが変化するため、溶媒和を考慮した解析が必要になります[4]。
3.1 溶媒中の励起状態計算（TDDFT）
- 時間依存密度汎関数理論（TDDFT）は、分子の励起状態エネルギーを求める計算手法です。
溶媒効果を考慮する場合、TDDFTとPCMを組み合わせることで、より実験に近いエネルギーが得られます。
コード②：PCMを適用したTDDFT計算
from pyscf import tdscf
# TDDFT計算（PCMあり）
mf = scf.RHF(mol).PCM()
mf.with_solvent.method = 'IEF-PCM'
mf.with_solvent.eps = 78.3553 # 水の誘電率
mf.kernel()
td = tdscf.TDA(mf).run() # TDDFT（TDA近似）で励起状態計算3.2 溶媒によるエネルギー変化の比較
異なる溶媒に対する分子のエネルギーの変化を計算し、溶媒の影響を評価します。
コード③：異なる溶媒でのSCFエネルギー比較
solvents = {'Water': 78.3553, 'Methanol': 32.613, 'Hexane': 1.89}
for name, eps in solvents.items():
mf = scf.RHF(mol).PCM()
mf.with_solvent.method = 'IEF-PCM'
mf.with_solvent.eps = eps
energy = mf.kernel()
print(f"溶媒: {name}, 誘電率: {eps}, SCFエネルギー: {energy:.6f} Hartree")→ 出力結果
溶媒: Water, 誘電率: 78.3553, SCFエネルギー: -113.873253 Hartree
溶媒: Methanol, 誘電率: 32.613, SCFエネルギー: -113.873003 Hartree
溶媒: Hexane, 誘電率: 1.89, SCFエネルギー: -113.867753 Hartreeこの結果から、誘電率の高い溶媒ほどエネルギーが低下し、分子がより安定化することが分かります。
その他の溶媒に関しては以下のような値となります必要に応じて参照してください。さらに情報が必要な場合は下記のサイトにアクセスしてください。
参考：https://gaussian.com/scrf/
solvent_dielectric_constants = {
"Water": 78.3553,
"Dimethylsulfoxide": 46.826,
"N,N-Dimethylformamide": 37.219,
"Nitromethane": 36.562,
"Methanol": 32.613,
"Ethanol": 24.852,
"Acetone": 20.493,
"Dichloroethane": 10.125,
"Dichloromethane": 8.93,
"Tetrahydrofuran": 7.4297,
"Chlorobenzene": 5.6968,
"Chloroform": 4.7113,
"Diethylether": 4.2400,
"Toluene": 2.3741,
"Benzene": 2.2706,
"1,4-Dioxane": 2.2099,
"Cyclohexane": 2.0160
}4. 溶媒環境での分子振動数解析
溶媒の誘電率は、分子の振動数（特に伸縮振動） に影響を与えることが知られています。本節では、ホルムアルデヒド（H₂CO）のC=O伸縮振動に着目し、異なる溶媒環境での振動数を比較します[5]。
4.1 振動数解析とは？
振動数解析（Frequency Analysis）は、分子の振動モード（伸縮振動、曲げ振動など）を計算する手法です。
計算により得られる振動数は、実験的な赤外（IR）スペクトルと比較することで、分子構造の同定や反応機構の理解に役立ちます。
本節では、異なる溶媒環境でホルムアルデヒドのC=O伸縮振動数がどのように変化するかを解析します。
4.2 構造最適化と振動数計算
PySCF properties モジュール
PySCFでは、分子のヘシアン行列（Hessian Matrix）の計算や熱力学解析（Thermodynamics Analysis）を行うための拡張モジュールとして、pyscf/properties が提供されています。
このモジュールを利用することで、振動数解析やエネルギーの補正計算がより簡単に実施可能になります。
例えば、ヘシアン計算に基づく振動数解析（Hessian Analysis）を行う際に、thermo.harmonic_analysis() を用いて振動数情報の抽出と解析が可能になります。
インストール方法：
!pip install git+https://github.com/pyscf/propertiesこのコマンドにより、最新のPySCF properties モジュールをGitHubから直接インストールできます。
以下のコードでは、異なる溶媒中でのホルムアルデヒドの構造最適化と振動数計算を行い、溶媒による影響を比較します。
コード④：異なる溶媒中でのホルムアルデヒドの構造最適化と振動数計算
import numpy as np
from pyscf import gto, dft, lib
from pyscf.geomopt import geometric_solver
from pyscf.hessian import rks, thermo
import matplotlib.pyplot as plt
# ホルムアルデヒドの分子構造（初期座標）
h2co = gto.M(atom="""
C 0.0000 0.0000 0.0000
O 1.2000 0.0000 0.0000
H -0.6000 0.9000 0.0000
H -0.6000 -0.9000 0.0000
""", basis='6-31G(d)')
# 溶媒リスト（誘電率） https://gaussian.com/scrf/
solvents = {
'Gas phase': None,
'Water': 78.3553,
'Acetonitrile': 35.688,
'Methanol': 32.613,
'Chloroform': 4.7113,
'Hexane': 1.89
}
# 各溶媒での振動数を格納
vibration_data = {}
for solvent, eps in solvents.items():
print(f"\n=== 計算開始：{solvent} ===")
# DFT計算の設定
mf = dft.RKS(h2co)
mf.xc = "B3LYP"
# PCMを適用する場合
if eps:
mf = mf.PCM()
mf.with_solvent.method = 'IEF-PCM'
mf.with_solvent.eps = eps
# 構造最適化
print(f"{solvent}: 構造最適化中...")
optimized_mf = geometric_solver.optimize(mf)
# DFT計算の設定
mf = dft.RKS(optimized_mf)
mf.xc = "B3LYP"
# DFT計算を実行
mf.kernel()
# 振動数計算のためのHessianオブジェクトを生成
hessian = rks.Hessian(mf)
# 振動数計算を実行
print(f"{solvent}: 振動数計算中...")
freqs = hessian.kernel()
freq_info = thermo.harmonic_analysis(mf.mol, freqs)
# 振動数を昇順にソートし、小さい方から3つ目をC=O伸縮振動と仮定
sorted_freqs = np.sort(freq_info["freq_wavenumber"])
co_stretch_freq = sorted_freqs[3] # 0-based index: 3番目（4つ目）
vibration_data[solvent] = co_stretch_freq
print(f"{solvent}: C=O 伸縮振動数（小さい方から3番目）: {co_stretch_freq:.2f} cm⁻¹")
import matplotlib.pyplot as plt
import numpy as np
# プロット設定
plt.figure(figsize=(8, 5))
# バーの色を変更（強調）
colors = ['skyblue' if v != max(vibration_data.values()) else 'red' for v in vibration_data.values()]
# バープロットの作成
bars = plt.bar(vibration_data.keys(), vibration_data.values(), color=colors)
# 軸ラベル
plt.xlabel("Solvent")
plt.ylabel("C=O Stretching Frequency (cm⁻¹)")
plt.title("Solvent Dependence of Formaldehyde C=O Stretching")
# Y軸のスケール調整
min_freq = min(vibration_data.values())
max_freq = max(vibration_data.values())
margin = (max_freq - min_freq) * 0.1 # マージンを10%追加
plt.ylim(min_freq - margin, max_freq + margin)
# 各バーに数値を表示
for bar in bars:
height = bar.get_height()
plt.text(bar.get_x() + bar.get_width()/2, height, f'{height:.1f} cm⁻¹', ha='center', va='bottom', fontsize=10)
# X軸の回転
plt.xticks(rotation=45)
# グリッド追加
plt.grid(axis='y', linestyle='--', alpha=0.7)
# 表示
plt.show()計算結果
異なる溶媒環境でのホルムアルデヒドのC=O伸縮振動数の結果を以下に示します。
| 溶媒 | 計算値 (cm⁻¹) |
|---|---|
| 気相 | 1851.6 |
| 水 | 1817.2 |
| アセトニトリル | 1818.1 |
| メタノール | 1818.3 |
| クロロホルム | 1827.4 |
| ヘキサン | 1839.4 |
実験では、アセトニトリル中と気相のC=O伸縮振動数の差が約23 cm⁻¹ であることが報告されています。今回の計算結果では 33.5 cm⁻¹（1851.6 – 1818.1 cm⁻¹） の差となりました。計算結果は実験データとおおよそ一致する傾向を示しているが、この差は、基底関数の選択、構造最適化による影響と考えられます。
グラフの説明
以下のグラフは、異なる溶媒環境でのC=O伸縮振動数の変化を示しています。
- 気相（Gas Phase） の振動数（1851.6 cm⁻¹）が最も高く、赤色で強調されています。
- 極性溶媒（水、アセトニトリル、メタノール）では、振動数が低下し、1817–1818 cm⁻¹の範囲に収まっています。
- 低極性溶媒（クロロホルム、ヘキサン）では、振動数が比較的高くなる傾向を示しています。
この結果から、溶媒の誘電率が大きいほど、C=O伸縮振動数が低下する傾向が確認できました。
5. 他の溶媒モデル（SMD, ddCOSMOなど）との比較
Polarizable Continuum Model（PCM） は溶媒効果を簡単に考慮できる手法ですが、他にもさまざまな溶媒モデルが提案されています。ここでは、SMD（Solvation Model Density） および ddCOSMO（domain-decomposition COSMO） をPCMと比較しながら紹介します。
5.1.1 SMD（Solvation Model Density）
SMDは、PCMの発展型であり、溶媒和自由エネルギーを高精度に計算するために設計された手法です。
このモデルは、誘電率以外の溶媒パラメータ（分子間相互作用や分散力など） を考慮し、より現実的な溶媒効果を再現します[6]。
SMDの特徴
✅ 溶媒パラメータを考慮（静電相互作用、誘起双極子相互作用、分散相互作用など）
✅ 非極性溶媒にも適用可能（PCMは極性溶媒向け）
✅ 溶媒和自由エネルギーをより精密に計算
コード⑤：SMDをPySCFで使用する方法
PySCFでは、SMDはIEF-PCMと組み合わせて利用できます。
from pyscf import gto, dft
# 分子定義（ホルムアルデヒド）
mol = gto.M(atom="""
C 0.0000 0.0000 0.0000
O 1.2000 0.0000 0.0000
H -0.6000 0.9000 0.0000
H -0.6000 -0.9000 0.0000
""", basis='6-31G(d)')
# SMDモデルの適用
mf = dft.RKS(mol).SMD()
mf.with_solvent.solvent = 'water' # 水を溶媒として指定
mf.kernel()PCMと異なり、SMDでは分子の分極性や分散相互作用も考慮されるため、溶媒効果の精度が向上する可能性があります。
5.1.2 ddCOSMO（Domain-Decomposition COSMO）
ddCOSMO（ドメイン分解COSMO）は、PCMの計算をより効率的に行うための改良手法です。
特に、大規模な系（タンパク質や高分子）に適用しやすいモデルとして知られています[7]。
ddCOSMOの特徴
✅ 大規模系に適用しやすい（計算負荷が低い）
✅ PCMと類似したアプローチだが、高精度な計算が可能
✅ 水などの高誘電率の溶媒に適用しやすい
コード⑥：ddCOSMOをPySCFで使用する方法
PySCFでは、SCF計算にddCOSMOを適用することが可能です。
from pyscf import gto, dft
# 分子定義（ホルムアルデヒド）
mol = gto.M(atom="""
C 0.0000 0.0000 0.0000
O 1.2000 0.0000 0.0000
H -0.6000 0.9000 0.0000
H -0.6000 -0.9000 0.0000
""", basis='6-31G(d)')
# ddCOSMOモデルの適用
mf = dft.RKS(mol).DDCOSMO()
mf.kernel()結論：どの溶媒モデルを使うべきか？
- PCM： 一般的な溶媒効果の考慮に適している（特に極性溶媒向け）。
- SMD： より現実的な溶媒効果を考慮し、有機溶媒や非極性溶媒でも精度の高い計算が可能。
- ddCOSMO： 大規模系で計算負荷を抑えながら溶媒効果を取り入れる場合に適している。
6. GPUを活用した高速な構造最適化と振動数解析
近年、量子化学計算においてGPU（Graphics Processing Unit）を活用することで、計算速度を大幅に向上させることが可能になっています[8, 9]。
PySCFでは、GPUを活用した密度汎関数理論（DFT）計算をサポートしており、溶媒効果を考慮した計算（PCM）と組み合わせることで、より高速かつ正確な結果を得ることができます。
本節では、GPUを用いてSCF計算を高速化し、その後PCMを適用、そのまま構造最適化・振動数解析までを一貫して行う方法について解説します。
6.1 PySCFでのGPU計算の設定
6.1.1 PySCFのGPU対応ライブラリ
PySCFでGPUを利用するためには、以下の拡張ライブラリが提供されています：
- gpu4pyscf（CUDA 12.x対応）
- GPUを活用した密度汎関数理論（DFT）計算の高速化が可能
gpu4pyscf.dft.rksモジュールを使用することで、GPU上での計算が実行可能
6.1.2 GPU対応版PySCFのインストール
Google Colabやローカル環境でGPUを利用する場合、以下のコマンドでGPU対応ライブラリをインストールします。
!pip install pyscf
!pip install gpu4pyscf-cuda12xこのセットアップにより、GPUアクセラレーションを活用したDFT計算が可能になります。
6.2 GPUを利用したPCM計算と構造最適化
GPUを利用したPCM計算を行い、そのまま構造最適化・振動数解析までを実行します。
実装の流れ
- 分子構造の定義（ホルムアルデヒド）
- GPUを用いたPCM計算（IEF-PCM）
- 構造最適化
- ヘシアン行列の計算
- 振動数解析
- C=O伸縮振動の特定
コード⑤：GPUを利用したPCM計算と構造最適化・振動数計算
import pyscf
from pyscf import gto, dft
from gpu4pyscf.dft import rks
from pyscf.hessian import thermo
from pyscf.geomopt import geometric_solver
import numpy as np
# 1. 分子の定義
# ホルムアルデヒドの構造を座標で指定（単位：オングストローム）
mol = gto.M(
atom="""
C -0.005636 0.000021 0.000000
O 1.202978 0.000056 -0.000000
H -0.598110 0.937848 -0.000000
H -0.598072 -0.937837 -0.000000
""",
basis='6-31G(d)', # 基底関数として6-31G(d)を選択
verbose=1 # 計算の詳細出力を抑制（必要に応じて調整）
)
# 2. PCMの設定（GPU対応）
# 水を溶媒としてPCMモデル（IEF-PCM）を適用
mf = rks.RKS(mol, xc='B3LYP').density_fit() # GPU対応のRKS DFT計算、B3LYP汎関数を使用
mf.with_solvent = pyscf.solvent.PCM(mol) # PCMモデルを有効化
mf.with_solvent.method = 'IEF-PCM' # IEF-PCM法を選択
mf.with_solvent.eps = 78.3553 # 水の誘電率
# 3. SCF計算（GPUを利用）
# 自己無撞着場（SCF）計算を実行してエネルギーを取得
e_dft = mf.kernel()
print(f"SCF Energy: {e_dft:.6f} Hartree")
# 4. 構造最適化（GPUを活用）
# 分子の構造を最適化して安定な形状を取得
optimized_mol = geometric_solver.optimize(mf)
print("Optimized geometry (Angstrom):")
print(optimized_mol.atom_coords())
# 5. ヘシアン行列の計算（GPU利用）
# 最適化された分子構造でヘシアン行列を計算
mf.mol = optimized_mol # 最適化された分子を計算に反映
h = mf.Hessian()
h_dft = h.kernel()
# 6. 振動解析
# ヘシアン行列から振動周波数を計算
freq_info = thermo.harmonic_analysis(optimized_mol, h_dft)
sorted_freqs = np.sort(freq_info["freq_wavenumber"]) # 周波数を昇順にソート
# C=O伸縮振動の特定
# ホルムアルデヒドの振動モードは6つ（3N-6=6, N=4）
# C=O伸縮振動は通常、高い周波数（例：1700 cm⁻¹付近）に現れる
co_stretch_freq = sorted_freqs[3] # 最大周波数をC=O伸縮振動と仮定
print(f"C=O stretch frequency: {co_stretch_freq:.2f} cm⁻¹")6.3 GPU計算のメリットと課題
GPUを活用するメリット
計算速度の向上
- SCF計算、構造最適化、振動数解析を高速化できる。
- 特に、溶媒効果を考慮した計算（PCM）ではGPUの恩恵が大きい。
メモリ効率の向上
- GPUは行列演算を並列計算できるため、メモリの負荷を抑えながら効率的に計算が可能。
- 計算時間が短縮されることで、複数の溶媒条件での比較も容易。
大規模系にも適用可能
- 有機分子や金属錯体などの大きな分子系の計算に適用可能。
GPU計算の課題
⚠ 一部のポスト-HF計算（MP2, CCSDなど）には制限あり
- 現時点では完全なポスト-HF法のGPUサポートは限定的。
⚠ GPUメモリ（VRAM）の制約
- 計算系が大きくなると、VRAMが不足する可能性がある。
- メモリが足りない場合は、並列計算の工夫が必要。
まとめ
いかがでしょうか！量子化学計算は溶液中と解離があると言われますが、PCMを使うことでより近づけることが出来ます。今回のコードを使って、大規模な分子系の溶媒効果解析をPySCFで計算してみて下さい！
参考文献
- PySCF 公式サイト: https://pyscf.org/
- Q. Sun et al., “Recent developments in the PySCF program package”, J. Chem. Phys., 153, 024109 (2020). DOI: 10.1063/5.0006074
- B. Mennucci, “Polarizable continuum model”, WIREs Comput. Mol. Sci., 2, 386 (2012). DOI: 10.1002/wcms.1086
- R. Improta et al., “Geometries and properties of excited states in the gas phase and in solution: Theory and application of a time-dependent density functional theory polarizable continuum model”, J. Chem. Phys., 124, 124504 (2006). DOI: 10.1063/1.2222361
- J. A. Barnes et al., “Solvent effect on vibrational frequencies: cryosolution experiments and density functional calculations”, J. Mol. Spectrosc., 192, 86 (1998). DOI: 10.1006/jmsp.1998.7624
- A. V. Marenich et al., “Universal solvation model based on solute electron density and on a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions”, J. Phys. Chem. B, 113, 6378 (2009). DOI: 10.1021/jp810292n
- F. Lipparini et al., “Quantum, classical, and hybrid QM/MM calculations in solution: General implementation of the ddCOSMO linear scaling strategy”, J. Chem. Phys., 141, 184108 (2014). DOI: 10.1063/1.4901304
- PySCF GPU4PySCFリポジトリ: https://github.com/pyscf/gpu4pyscf/blob/master/README.md
- J. Wu et al., “Enhancing GPU-acceleration in the Python-based Simulations of Chemistry Framework”, arXiv preprint (2024). DOI: 10.48550/arXiv.2404.09452
- 西長亨、本田康, 『有機化学のための 量子化学計算入門: Gaussianの基本と有効活用のヒント』, 裳華房 (2022), P. 114-119.
Gaussianを使った量子化学計算の初心者向け技術書を販売中
しばしば出くわすエラーへの対処法をはじめ
Gaussianと無料ソフトウェア Avogadro を組み合わせた物性解析手法が学べます！
