【PySCF】NMRスペクトルを描写【Pythonで始める量子化学計算】

この記事では、PySCFを用いたNMRスペクトルの計算と描写について解説します。分子構造の解析や実験データとの比較を行う際に役立つ内容を提供します。NMRの基本的な理論から始まり、PySCFを使った計算手順、得られた結果の解釈方法と可視化の実践例、さらに応用例まで幅広くカバーしています。この記事を読むことで、量子化学計算のスキルを向上させ、NMRスペクトルの計算を通じて分子の構造や性質を深く理解することができるようになります。

動作検証済み環境

Google colab, pyscf 2.7.0, py3Dmol, PubChemPy 1.0.4., rdkit 2024.3.5, Avogadro 1.2.0

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


参考文献[1]から引用

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

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

  • モジュラー設計:PySCFは、計算の各ステップ(分子の構造定義、基底状態の計算、励起状態の計算など)をモジュールとして提供しており、必要な部分だけを組み合わせて使用することができます。
  • 高性能計算:内部に高度なアルゴリズムを使用しており、大規模な計算でも高速に処理を行うことができます。
  • 柔軟性: Pythonベースであるため、スクリプトを通じて簡単に操作や拡張が可能です。ユーザーは、計算手法やパラメータをカスタマイズして特定のニーズに合わせることができます。
  • 拡張性: オープンソースであるため、研究者や開発者が自分の研究に合わせてPySCFを拡張することが容易です。新しいアルゴリズムや方法論を追加することで、独自の研究を進めることができます。また、拡張機能が開発され続けている事も特徴です。

1. NMRスペクトルとは


NMR(核磁気共鳴)スペクトルは、有機化合物をはじめとした分子の構造解析に広く用いられる手法です。分子中の特定核種(¹H、¹³Cなど)が外部磁場からどのような磁気的影響(遮蔽あるいは反遮蔽)を受けるかを調べることで、その核種の化学環境を知ることができます。実験では、分子ごとに特徴的な化学シフト値を示すため、分子の同定や構造解析に不可欠な情報となります。

NMRスペクトルの中でも、¹H NMR¹³C NMRは特に有機化合物の解析で重要です。たとえば¹H NMRは、水素原子の置かれている環境のわずかな違いを示すため、分子構造や立体配置の推定に大きく貢献します。計算化学の手法を組み合わせることで、実験スペクトルとの比較や未知分子の予測など、研究・開発に役立つ多角的なアプローチが可能です。

2. 必要な準備


2.1 必要なPython環境とライブラリ

  • Python 3.7以降(推奨)
  • NumPy, SciPy, Matplotlib, Plotlyなどの科学技術計算・可視化用ライブラリ
  • RDKitやpy3Dmolなど、分子構造を扱うためのライブラリ(必要に応じて)
  • その他、本書で利用しているライブラリとして pubchempy, ipywidgets, plotly など

Google Colabを用いる場合は、上記ライブラリをpipで導入すればすぐに実行可能です。Jupyter Notebook等でも同様ですが、GPU利用時は環境(CUDAのバージョンなど)に注意が必要です。

下記の流れで環境設定をしてください。

2.2 PySCFのインストール方法

PySCFはPyPIからインストール可能です。そのほか必要なコマンドで導入します。

# Google Colabでのインストールコマンド
!pip install pyscf
!pip install numpy
!pip install geometric
!pip install rdkit-pypi
!pip install pubchempy
!pip install plotly
!pip install matplotlib
!pip install py3Dmol
!pip install gpu4pyscf-cuda12x

GPU計算をサポートするgpu4pyscfを利用したい場合は、GPU対応の環境が整っているか(ドライバやCUDAのバージョンなど)を確認した上で、同様にインストールしてください。Google Colabであれば、ランタイムランタイムのタイプを変更 で「T4」を選択することでCUDAデバイスが使用可能になります。

NMRスペクトルの計算には以下のコマンドでpyscf.propertiesを追加インストールする必要があります。

!pip install git+https://github.com/pyscf/properties

これはPySCFの拡張モジュールで、NMR以外にもいくつかの分子物性計算が行えるようになります。

3. NMR計算の基本


3.1 NMR計算の理論背景

NMR化学シフトの計算は、外部磁場中の核に対する電子の局所的な“遮蔽”の程度を理論的に評価することで求められます。分子軌道法を拡張して磁場の影響を取り込むアプローチとしては、以下のような方法が知られています。

  • CSGT法(Continuous Set of Gauge Transformations)かつてはGaussianなどの一部計算ソフトでも採用されていましたが、精度に課題があり、現在ではあまり用いられていません。

  • GIAO法(Gauge-Including Atomic Orbital)化学シフトの計算において、安定した結果が得られるため、多くのプログラムで標準的に使用されています。PySCFもGIAO法を標準でサポートしており、高精度のNMR計算が可能です。

3.2 化学シフトとは

NMR実験で一般的に用いられる基準物質がテトラメチルシラン(TMS)です。テトラメチルシランの¹Hまたは¹³C核種の「遮蔽定数(isotropic shielding constant)」をσ(TMS)とすると、目的分子の対応核の遮蔽定数σ(target)との差分によって化学シフトδ(ppm単位)が定義されます。

$$ δ(target)= \sigma(\text{TMS}) – \sigma(\text{target}) $$

計算手順としては、

  1. TMSの構造最適化・遮蔽定数の計算
  2. 目的分子の構造最適化・遮蔽定数の計算
  3. 遮蔽定数の差をとる

ことで理論的な化学シフトを導出します。

汎関数の選び方

遮蔽定数は使う汎関数によって値が異なり、レベル(HFやDFTなど)や基底関数系がシフト値に影響を与えます。

  • DFT計算でも計算は可能ですが、HFによる計算を推奨している論文も多く見られます。
  • 一方、B3LYP/6-31G(d)での構造最適化を行い、その後HF/6-31G(d)でNMRを計算する、いわゆるHF/6-31G(d)//B3LYP/6-31G(d)という方法もよく採用されます。

著者の経験上、HF/6-31G(d) が実験値に近いシフトを再現しやすい場合が多くあります。さらに、大型の基底関数(aug-cc-pVDZなど)を用いると、より高精度が期待できますが、計算コストが増加する点は留意が必要です。

いずれのレベルを選択するにせよ、化学シフトの相対値を比較する上では、一定の再現性が期待できます。研究の目的や分子の大きさ、計算コストとのバランスを見ながら手法を選んでください。

4. PySCFを用いたNMR計算の手順


以下の手順で、PySCFを用いてNMRスペクトルを理論的に算出します。

4.1 全体の流れ

  1. 分子構造の定義TMSの座標データを基に分子構造を指定します(TMSはその高い対称性によりNMR基準物質として最適です)。

  2. SCF(自己無撞着場)計算HFまたはDFTを用いて電子状態を解き、分子軌道を計算します。

  3. NMR遮蔽定数(Shielding Constant)の計算GIAO法(Gauge-Including Atomic Orbital)を用いて各原子核の遮蔽定数を求めます。

  4. 化学シフトの算出TMSを基準とし、各原子核の遮蔽定数から化学シフトを計算します。

  5. 結果の可視化計算結果を基に、NMRスペクトルをプロットして解析します。

4.2 TMSについて

NMR計算で基準物質として使用されるTMS(テトラメチルシラン)は、以下の特徴を持っています:

  • 正四面体対称性Si原子を中心に4つのメチル基が等価に結合しており、すべての¹Hおよび¹³C原子が同等の化学環境にあります。

  • 単一ピーク

    ¹H NMR:水素原子が等価であるため、1種類のピークが得られます。
    ¹³C NMR:炭素原子も等価で、単一のピークとして現れます。
  • 安定性化学的に安定で反応性が低いため、実験と計算の基準物質として適しています。

4.3 分子構造データの準備(TMSを例に)

今回はTMSを例に計算します。以下では、B3LYP/6-31G(d)法でTd対称性を保ったまま構造最適化した座標を、直接Pythonの文字列で定義しています。分子構造ファイル(XYZ, Mol, SDFなど)からの読み込みや、PubChem/RDKitを介したSMILES指定でも問題ありません。

nmr_coords = """
Si  0.00000000  0.00000000  0.00000000
C   1.09477108  1.09477108  1.09477108
H   1.74260012  1.74260012  0.49163604
H   0.49163604  1.74260012  1.74260012
H   1.74260012  0.49163604  1.74260012
C  -1.09477108 -1.09477108  1.09477108
H  -1.74260012 -1.74260012  0.49163604
H  -0.49163604 -1.74260012  1.74260012
H  -1.74260012 -0.49163604  1.74260012
C  -1.09477108  1.09477108 -1.09477108
H  -1.74260012  1.74260012 -0.49163604
H  -1.74260012  0.49163604 -1.74260012
H  -0.49163604  1.74260012 -1.74260012
C   1.09477108 -1.09477108 -1.09477108
H   1.74260012 -1.74260012 -0.49163604
H   0.49163604 -1.74260012 -1.09477108
H   1.74260012 -0.49163604 -1.74260012
"""

4.4 NMR遮蔽定数と化学シフトの計算

以下は、遮蔽定数を計算し、化学シフトを算出するための関数です。

def calculate_nmr_shielding(mol, method='HF', xc='B3LYP', selected_nmr_types=['1H', '13C'],
                            scaling_factors={'1H':0.966, '13C':0.966},
                            TMS_shifts={'1H':31.61, '13C':190.0},
                            basis_nmr="aug-cc-pVDZ"):
    """
    NMR遮蔽定数を計算し、化学シフトを算出する関数。
    """
    # 分子オブジェクトをコピーし、基底関数を変更
    mol_bkp = mol.copy()
    mol_bkp.basis = basis_nmr

    # 計算手法の選択
    if method.upper() == 'HF':
        mf = scf.RHF(mol_bkp)
    elif method.upper() == 'DFT':
        mf = dft.RKS(mol_bkp)
        mf.xc = xc
    else:
        raise ValueError("Invalid method. Choose 'HF' or 'DFT'.")

    # SCF計算
    mf.conv_tol = 1e-9
    mf.max_cycle = 300
    mf.kernel()

    if not mf.converged:
        raise RuntimeError("SCF calculation did not converge.")

    # NMR遮蔽定数を計算
    nmr_calc = nmr.RHF(mf) if method.upper() == 'HF' else nmr.RKS(mf)
    shielding_tensors = nmr_calc.kernel()

    # 化学シフトを計算
    chemical_shifts = []
    for idx, tensor in enumerate(shielding_tensors):
        atom_symbol = mol.atom_symbol(idx)

        # 核種ごとの計算
        if atom_symbol == 'H' and '1H' in selected_nmr_types:
            ref = TMS_shifts['1H']
            scaling = scaling_factors['1H']
        elif atom_symbol == 'C' and '13C' in selected_nmr_types:
            ref = TMS_shifts['13C']
            scaling = scaling_factors['13C']
        else:
            continue

        # 化学シフトを算出
        iso_shielding = np.trace(tensor) / 3.0
        delta = ref - (iso_shielding * scaling)
        chemical_shifts.append({
            'atom_index': idx,
            'element': atom_symbol,
            'iso_shielding': iso_shielding,
            'chemical_shift': delta,
            'type': '1H' if atom_symbol == 'H' else '13C'
        })

    return chemical_shifts

4.4 実行の流れ

main_nmr関数は、TMS分子のNMRスペクトルを計算し、結果を出力する一連のプロセスを管理するメイン関数です。以下にそのコードと各ステップの詳細を示します。


コード

def main_nmr():
    # 1. TMS分子の定義
    mol = define_molecule_from_coords(nmr_coords, basis='6-31G(d)', charge=0, spin=0)

    # 2. NMR遮蔽定数の計算
    chemical_shifts = calculate_nmr_shielding(mol, method='HF', xc='B3LYP')

    # 3. 結果の出力
    for shift in chemical_shifts:
        print(f"Atom {shift['atom_index']} ({shift['element']}): "
              f"Chemical Shift = {shift['chemical_shift']:.2f} ppm")

解説

ステップ 1: TMS分子の定義

mol = define_molecule_from_coords(nmr_coords, basis=’6-31G(d)’, charge=0, spin=0)

この部分では、TMS分子をPySCFの分子オブジェクトに変換しています。

  1. define_molecule_from_coords 関数
    • 入力: TMSの座標データ(nmr_coords)、基底関数(basis)、分子の電荷(charge)、スピン多重度(spin)。
    • 出力: PySCFで計算に使用できる分子オブジェクト(mol)。
    • 詳細:
      TMSは正四面体対称性を持つため、各メチル基の原子は等価な位置に配置されています。これにより計算結果も対称性を反映します
      基底関数は 6-31G(d) を使用しており、計算の精度と効率のバランスを考慮しています。
    • ポイント
      高対称性分子のため、計算の収束が比較的早い。
      電荷が0、スピンが0であるため、閉殻系の計算が可能。

ステップ 2: NMR遮蔽定数の計算

chemical_shifts = calculate_nmr_shielding(mol, method=’HF’, xc=’B3LYP’)

このステップでは、分子の遮蔽定数を計算し、TMSを基準とした化学シフトを算出します。

  1. calculate_nmr_shielding 関数
    • 入力: 分子オブジェクト(mol)、計算手法(method)、汎関数(xc)。
    • 出力: 各原子の化学シフトを含むリスト(chemical_shifts)。
    • 詳細:
      • method='HF' はHartree-Fock法を指定。高速でありながら適度な精度を持つ計算手法です。
      • 遮蔽定数の計算にはGIAO法が使用され、外部磁場による影響を分子軌道法で正確に反映します。
      • 化学シフトは、遮蔽定数(isotropic shielding)を基にTMSとの差分を計算。
  2. 計算結果の構造
    • chemical_shifts はリスト形式で、各原子の情報を辞書として格納しています。辞書の内容:
      • atom_index: 原子のインデックス。
      • element: 原子の種類(例: H, C)。
      • iso_shielding: 遮蔽定数。
      • chemical_shift: 化学シフト(ppm)。
  3. ポイント
    • 高対称性分子では、遮蔽定数や化学シフトの結果が等価な原子で一致することが期待されます。
    • 結果は実験値との比較により検証可能。

ステップ 3: 結果の出力

for shift in chemical_shifts: print(f”Atom {shift[‘atom_index’]} ({shift[‘element’]}): ” f”Chemical Shift = {shift[‘chemical_shift’]:.2f} ppm”)

計算結果を整形してコンソールに出力します。

  1. 内容
    • 各原子について、インデックス、元素記号、化学シフト(ppm単位)を表示します。
  2. 出力

    以下のような形式で出力されます:


    Atom 1 (C) - 13C: Isotropic Shielding = 200.31 ppm, 
    Atom 2 (H) - 1H: Isotropic Shielding = 32.62 ppm,
    Atom 3 (H) - 1H: Isotropic Shielding = 32.62 ppm,
    Atom 4 (H) - 1H: Isotropic Shielding = 32.62 ppm,
    Atom 5 (C) - 13C: Isotropic Shielding = 200.31 ppm,
    Atom 6 (H) - 1H: Isotropic Shielding = 32.62 ppm,
    Atom 7 (H) - 1H: Isotropic Shielding = 32.62 ppm,
    Atom 8 (H) - 1H: Isotropic Shielding = 32.62 ppm,
    Atom 9 (C) - 13C: Isotropic Shielding = 200.31 ppm,
    Atom 10 (H) - 1H: Isotropic Shielding = 32.62 ppm,
    Atom 11 (H) - 1H: Isotropic Shielding = 32.62 ppm,
    Atom 12 (H) - 1H: Isotropic Shielding = 32.62 ppm,
    Atom 13 (C) - 13C: Isotropic Shielding = 200.31 ppm,
    Atom 14 (H) - 1H: Isotropic Shielding = 32.62 ppm,
    Atom 15 (H) - 1H: Isotropic Shielding = 32.62 ppm,
    Atom 16 (H) - 1H: Isotropic Shielding = 32.62 ppm,
    • TMSの¹H NMRでは、すべての水素原子が等価なため、化学シフトは同一です。

  3. ポイント
    • HF/6-31G(d)におけるTMSの値として、後続の化学シフト解析に使用できます。
    • 実験データと比較して、計算手法や基底関数の適切さを確認できます。

5. 化学シフトとNMRスペクトルの可視化


5.1 構造最適化~振動数計算~NMRスペクトル計算(アセトンの解析)

NMR計算の精度を上げるには、対象分子の構造を十分に最適化しておくことが重要です。以下に示すサンプルコードでは、アセトンの構造最適化(DFT/B3LYPレベル)→振動数計算→NMRスペクトル計算を順次行っています。Google ColabなどのGPU環境が使える場合はgpu4pyscfを利用すると高速化が期待できます。まずは、コードの解説から始めます。

アセトン¹H NMRの場合、計算手順に従いTMSの¹H遮蔽定数(約32.6 ppm)との差をとることで、理論的な化学シフトを求めます。具体となり、

$δH(Acetone)=σ(TMS)−σ(Acetone)×(スケーリングファクター)$

で算出され、計算条件(HFやB3LYP、基底関数系など)次第で2.1 ppm前後の値が得られます。これは実験値と概ね同程度のオーダーであるため、相対値としては十分に参考になると考えられます。絶対値を合わせる場合は、スケーリングファクターを調整してください。

5.2 化学シフトの結果の確認方法

計算の最後に出力される例として、下記のような遮蔽定数および化学シフト値が得られます。

Atom 1 (C) - 13C: Isotropic Shielding = 172.20 ppm, Chemical Shift = 28.11 ppm
Atom 2 (C) - 13C: Isotropic Shielding = 172.20 ppm, Chemical Shift = 28.11 ppm
Atom 3 (C) - 13C: Isotropic Shielding = -10.59 ppm, Chemical Shift = 210.90 ppm
Atom 4 (H) - 1H: Isotropic Shielding = 30.41 ppm, Chemical Shift = 2.21 ppm
Atom 5 (H) - 1H: Isotropic Shielding = 30.63 ppm, Chemical Shift = 1.99 ppm
Atom 6 (H) - 1H: Isotropic Shielding = 30.63 ppm, Chemical Shift = 1.99 ppm
Atom 7 (H) - 1H: Isotropic Shielding = 30.63 ppm, Chemical Shift = 1.99 ppm
Atom 8 (H) - 1H: Isotropic Shielding = 30.41 ppm, Chemical Shift = 2.21 ppm
Atom 9 (H) - 1H: Isotropic Shielding = 30.63 ppm, Chemical Shift = 1.99 ppm

たとえば水素原子(1H)に注目すると、出力された化学シフトは約2.0~2.2 ppm程度の範囲にあります。アセトンのH原子はカルボニル基側に平行に配置されたもの2個と反対側に向いた4個の2種類が存在しますが、NMRの時間スケールでは十分早く自由回転しているので遮蔽定数も平均を取ってください。平均するとおよそ2.1 ppm程度になります。

参考までに、産業技術総合研究所(産総研)のデータベースに記載されているアセトンの¹H NMR実験値が約2.16 ppmであることを踏まえると、計算結果は絶対値としても実験と非常に近い値を示していることが分かります。

このように、理論計算と実験値が1 ppm未満の差に収まる場合は、NMR計算としては相当良い一致と考えられます。一方で、分子の周囲環境(溶媒効果や温度など)をより正確に再現したい場合には、PCM(Polarizable Continuum Model)などの溶媒モデルを導入してみるのも一つの方法です。

5.3 可視化ツールを用いたNMRスペクトルの描写

コードでは、Plotlyを用いてガウシアン型ピークを仮定した簡単なNMRスペクトルをプロットしています。peak_widthを変更すればピークの幅を調整可能です。また、実際の実験スペクトルに近づけたい場合には、積分値の表示やベースライン処理など、追加の処理を行うことも考えられます。

コード

以上を得られるコードは以下です。Google Colabで実行すると上記の結果が得られます。

# 必要なライブラリのインポート
import os
import numpy as np
import matplotlib.pyplot as plt
import plotly.graph_objects as go
from pyscf import gto, scf, dft, lib
from pyscf.geomopt import geometric_solver
from pyscf.prop import nmr
from gpu4pyscf import dft as gpu_dft
import pubchempy as pcp
from google.colab import drive
import py3Dmol
from pyscf.hessian import rks, thermo

# RDKitのインポート
from rdkit import Chem
from rdkit.Chem import AllChem

# Google Driveをマウント
drive.mount('/content/drive')

# フォルダの作成とパスの設定
folder_path = "/content/drive/MyDrive/pyscf_NMR"
base_filename = "Aceton"
path = os.path.join(folder_path, base_filename)

if not os.path.exists(folder_path):
    os.makedirs(folder_path)
    print(f"Created folder: {folder_path}")

# 並列計算のためにスレッド数を設定(必要に応じて調整)
lib.num_threads(4)  # 利用可能なCPUコア数に応じて設定

def define_molecule(input_method='name', molecule_input=None, unit='Angstrom', basis='aug-cc-pVDZ', charge=0, spin=0):
    """
    分子をSMILES、名前、または直接座標から定義する関数。

    parameters:
        input_method (str): 'smiles', 'name', または 'coords' のいずれか。
        molecule_input (str):
            - 'smiles'の場合はSMILES文字列。
            - 'name'の場合は分子名。
            - 'coords'の場合はPySCFのatom定義文字列。
        unit (str): 座標の単位 ('Angstrom' または 'Bohr')。
        basis (str): 基底関数。
        charge (int): 分子の電荷。
        spin (int): 分子のスピン。

    returns:
        mol_pyscf (pyscf.gto.Mole): PySCFで定義された分子オブジェクト。
    """
    if input_method not in ['smiles', 'name', 'coords']:
        raise ValueError("input_method は 'smiles', 'name', または 'coords' のいずれかでなければなりません。")

    if input_method == 'coords':
        if molecule_input is None:
            raise ValueError("input_method が 'coords' の場合、molecule_input にatom定義文字列を指定してください。")
        # PySCFのMoleオブジェクトを直接作成
        mol_pyscf = gto.M(
            atom=molecule_input,
            basis=basis,
            charge=charge,
            spin=spin,
            symmetry=False,
            verbose=4
        )
        print("座標から分子オブジェクトが定義されました。")
        return mol_pyscf

    elif input_method == 'smiles':
        if molecule_input is None:
            raise ValueError("input_method が 'smiles' の場合、molecule_input にSMILES文字列を指定してください。")
        # RDKitで分子を定義
        mol_rdkit = Chem.MolFromSmiles(molecule_input)
        if mol_rdkit is None:
            raise ValueError("SMILES文字列からRDKit分子オブジェクトの作成に失敗しました。")

        # 水素原子を追加
        mol_rdkit = Chem.AddHs(mol_rdkit)

        # 3D構造の生成
        result = AllChem.EmbedMolecule(mol_rdkit, AllChem.ETKDG())
        if result != 0:
            raise ValueError("3D構造の埋め込みに失敗しました。")

        # UFFによる構造最適化
        result = AllChem.UFFOptimizeMolecule(mol_rdkit)
        if result != 0:
            print("UFFによる構造最適化が完全に収束しませんでした。")

        # XYZ形式で分子構造を取得
        xyz_coords = Chem.MolToXYZBlock(mol_rdkit)

        # 原子の数とヘッダー情報を追加
        num_atoms = mol_rdkit.GetNumAtoms()
        xyz_full = f"{num_atoms}\\nMolecule defined by SMILES\\n" + xyz_coords

        # 必要な行のみを抽出(5行目以降)
        xyz_lines = xyz_full.splitlines()
        xyz_trimmed = "\\n".join(xyz_lines[4:])

        # XYZデータから原子シンボルと座標をパース
        symbols = []
        coords = []
        for line in xyz_trimmed.strip().split('\\n'):
            parts = line.strip().split()
            if len(parts) < 4:
                continue
            symbol = parts[0]
            try:
                x, y, z = map(float, parts[1:4])
            except ValueError:
                print(f"無効な座標値が含まれています: {line}")
                continue
            symbols.append(symbol)
            coords.append([x, y, z])

        if unit.lower() in ['bohr', 'bohr_unit', 'bohr_units']:
            # ボーア単位からアンガストローム単位に変換
            coords = [[x * 0.529177, y * 0.529177, z * 0.529177] for x, y, z in coords]

        # PySCF用の原子定義文字列を作成
        atoms = [f"{sym} {x:.6f} {y:.6f} {z:.6f}" for sym, (x, y, z) in zip(symbols, coords)]
        atom_str = "\\n".join(atoms)

        # PySCFのMoleオブジェクトを作成
        mol_pyscf = gto.M(
            atom=atom_str,
            basis=basis,
            charge=charge,
            spin=spin,
            symmetry=False,
            verbose=4
        )

        print("SMILESから分子オブジェクトが定義されました。")
        return mol_pyscf

    elif input_method == 'name':
        if molecule_input is None:
            raise ValueError("input_method が 'name' の場合、molecule_input に分子名を指定してください。")
        # PubChemから化合物を取得
        compounds = pcp.get_compounds(molecule_input, 'name')
        if not compounds:
            raise ValueError("指定された分子名に一致する化合物が見つかりません。")

        # 最初の化合物を選択
        compound = compounds[0]

        # SMILESを取得
        smiles = compound.isomeric_smiles
        if not smiles:
            raise ValueError("化合物のSMILESが取得できませんでした。")

        # SMILESからRDKit分子オブジェクトを作成
        mol_rdkit = Chem.MolFromSmiles(smiles)
        if mol_rdkit is None:
            raise ValueError("SMILESからRDKit分子オブジェクトの作成に失敗しました。")

        # 水素原子を追加
        mol_rdkit = Chem.AddHs(mol_rdkit)

        # 3D構造の生成
        result = AllChem.EmbedMolecule(mol_rdkit, AllChem.ETKDG())
        if result != 0:
            raise ValueError("3D構造の埋め込みに失敗しました。")

        # UFFによる構造最適化
        result = AllChem.UFFOptimizeMolecule(mol_rdkit)
        if result != 0:
            print("UFFによる構造最適化が完全に収束しませんでした。")

        # XYZ形式で分子構造を取得
        xyz_coords = Chem.MolToXYZBlock(mol_rdkit)

        # 原子の数とヘッダー情報を追加
        num_atoms = mol_rdkit.GetNumAtoms()
        xyz_full = f"{num_atoms}\\nMolecule defined by name\\n" + xyz_coords

        # 必要な行のみを抽出(5行目以降)
        xyz_lines = xyz_full.splitlines()
        xyz_trimmed = "\\n".join(xyz_lines[4:])

        # XYZデータから原子シンボルと座標をパース
        symbols = []
        coords = []
        for line in xyz_trimmed.strip().split('\\n'):
            parts = line.strip().split()
            if len(parts) < 4:
                continue
            symbol = parts[0]
            try:
                x, y, z = map(float, parts[1:4])
            except ValueError:
                print(f"無効な座標値が含まれています: {line}")
                continue
            symbols.append(symbol)
            coords.append([x, y, z])

        if unit.lower() in ['bohr', 'bohr_unit', 'bohr_units']:
            # ボーア単位からアンガストローム単位に変換
            coords = [[x * 0.529177, y * 0.529177, z * 0.529177] for x, y, z in coords]

        # PySCF用の原子定義文字列を作成
        atoms = [f"{sym} {x:.6f} {y:.6f} {z:.6f}" for sym, (x, y, z) in zip(symbols, coords)]
        atom_str = "\\n".join(atoms)

        # PySCFのMoleオブジェクトを作成
        mol_pyscf = gto.M(
            atom=atom_str,
            basis=basis,
            charge=charge,
            spin=spin,
            symmetry=False,
            verbose=1
        )

        print("名前から分子オブジェクトが定義されました。")
        return mol_pyscf

def optimize_geometry(mol, path, xc='B3LYP', max_cycle=300, conv_params=None, basis_opt="aug-cc-pVDZ"):
    """
    GPU対応PySCFを用いて構造最適化を行う関数。

    parameters:
        mol (pyscf.gto.Mole): 最適化対象の分子オブジェクト。
        path (str): 出力ファイルのパス。
        xc (str): DFT汎関数。
        max_cycle (int): 最適化の最大サイクル数。
        conv_params (dict): 収束条件の辞書。
        basis_opt (str): 構造最適化に使用する基底関数。

    returns:
        mf (gpu4pyscf.dft.RKS): 最適化後のDFTオブジェクト。
    """
    if conv_params is None:
        conv_params = {
            'convergence_energy': 1e-6,  # エネルギーの収束閾値 (Eh)
            'convergence_grms': 3e-4,    # 勾配のRMS収束閾値 (Eh/Bohr)
            'convergence_gmax': 4.5e-4,  # 最大勾配の収束閾値 (Eh/Bohr)
            'convergence_drms': 1.2e-3,  # 変位のRMS収束閾値 (Angstrom)
            'convergence_dmax': 1.8e-3   # 最大変位の収束閾値 (Angstrom)
        }

    # GPU対応DFT計算の設定
    mf = gpu_dft.RKS(mol)
    mf.xc = xc
    mf.max_cycle = max_cycle
    mf.chkfile = path + '.chk'

    energies = []
    geometries = []

    # コールバック関数の定義
    def callback(envs):
        current_mf = envs['g_scanner'].base
        energies.append(current_mf.e_tot)
        geometries.append(mol.atom_coords())  # PySCFのメソッド
        print(f"Step {len(energies)}: Energy = {current_mf.e_tot:.6f} Eh")

    # 構造最適化の実行
    try:
        geometric_solver.optimize(mf, callback=callback, **conv_params)
    except Exception as e:
        print(f"構造最適化中にエラーが発生しました: {e}")
        return None

    # エネルギー収束のプロット
    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()

    # ジオメトリーの保存
    for i, coords in enumerate(geometries):
        filename = f"{path}_geometry_step_{i+1}.xyz"
        try:
            with open(filename, 'w') as f:
                f.write(f"{mol.natm}\\n")
                f.write(f"Step {i+1} Energy: {energies[i]} Eh\\n")
                for j in range(mol.natm):
                    symbol = mol.atom_symbol(j)
                    x, y, z = coords[j]
                    f.write(f"{symbol} {x:.6f} {y:.6f} {z:.6f}\\n")
        except Exception as e:
            print(f"ジオメトリーの保存中にエラーが発生しました: {e}")

    print("All geometries have been saved.")

    # SCF収束の確認
    if not energies:
        print("構造最適化が正常に完了しませんでした。")
        return None
    else:
        print("構造最適化が正常に完了しました。")
        return mf

def calculate_nmr_shielding(mol, method='HF', xc='B3LYP', selected_nmr_types=['1H', '13C'],
                            scaling_factors={'1H':0.966, '13C':0.966},
                            TMS_shifts={'1H':31.61, '13C':200.31},
                            basis_nmr="aug-cc-pVDZ"):
    """
    NMR遮蔽定数を計算し、化学シフトを算出する関数。HFまたはDFTを選択可能。

    parameters:
        mol (pyscf.gto.Mole): 分子オブジェクト。
        method (str): 計算手法 ('HF' または 'DFT')。
        xc (str): DFT汎関数(DFTを選択した場合に使用)。
        selected_nmr_types (list): 計算対象のNMR核種(例: ['1H', '13C'])。
        scaling_factors (dict): スケーリングファクター。
        TMS_shifts (dict): TMS参照シフト。
        basis_nmr (str): NMR計算に使用する基底関数。

    returns:
        chemical_shifts (list): 計算された化学シフトのリスト。
    """
    # 基底関数をNMR計算用に変更
    mol_bkp = mol.copy()
    mol_bkp.basis = basis_nmr

    # 計算手法の選択
    if method.upper() == 'HF':
        print("Using Hartree-Fock (HF) method for NMR shielding calculation.")
        mf = scf.RHF(mol_bkp)
    elif method.upper() == 'DFT':
        print(f"Using Density Functional Theory (DFT) method for NMR shielding calculation with functional {xc}.")
        mf = dft.RKS(mol_bkp)
        mf.xc = xc
    else:
        print("Invalid method selected. Choose 'HF' or 'DFT'.")
        return []

    # SCF計算の設定
    mf.conv_tol = 1e-9  # 収束基準を厳しく設定
    mf.max_cycle = 300  # 最大反復回数を増やす

    # SCF計算の実行
    try:
        mf.kernel()
    except Exception as e:
        print(f"SCF計算中にエラーが発生しました: {e}")
        return []

    # SCF計算の収束確認
    if not mf.converged:
        print("SCF計算が収束しませんでした。")
        return []
    else:
        print("SCF計算が正常に収束しました。")

    # NMR計算の実行
    try:
        nmr_calc = nmr.RHF(mf)
        shielding_tensors = nmr_calc.kernel()
    except Exception as e:
        print(f"NMR遮蔽定数の計算中にエラーが発生しました: {e}")
        return []

    # 化学シフトの計算
    chemical_shifts = []
    for idx, tensor in enumerate(shielding_tensors):
        atom_symbol = mol.atom_symbol(idx)
        if atom_symbol == 'H' and '1H' in selected_nmr_types:
            nmr_type = '1H'
            sigma_ref = TMS_shifts.get('1H', 31.61)  # デフォルト値を設定
            scaling = scaling_factors.get('1H', 1.0)
        elif atom_symbol == 'C' and '13C' in selected_nmr_types:
            nmr_type = '13C'
            sigma_ref = TMS_shifts.get('13C', 190.0)  # 実際には適切な参照値を設定
            scaling = scaling_factors.get('13C', 1.0)
        else:
            continue  # 対応していない核種はスキップ

        iso_shielding = np.trace(tensor) / 3.0
        delta = sigma_ref - (iso_shielding * scaling)

        chemical_shifts.append({
            'atom_index': idx,
            'element': atom_symbol,
            'iso_shielding': iso_shielding,
            'chemical_shift': delta,
            'type': nmr_type
        })
        print(f"Atom {idx} ({atom_symbol}) - {nmr_type}: Isotropic Shielding = {iso_shielding:.2f} ppm, Chemical Shift = {delta:.2f} ppm")

    if not chemical_shifts:
        print("指定されたNMR核種に対応するデータがありませんでした。")

    return chemical_shifts

def calculate_vibrational_frequencies(path, xc='B3LYP'):
    """
    振動数計算を行い、結果を表示する関数。

    parameters:
        mf (gpu4pyscf.dft.RKS): DFTオブジェクト。
        path (str): 出力ファイルのパス。

    returns:
        None
    """
    try:
        # チェックポイントファイルから分子構造を読み込み
        mol = lib.chkfile.load_mol(path + ".chk")
        mf = dft.RKS(mol)
        mf.chkfile = path + ".chk"
        mf.xc = xc
        # DFT計算を実行
        mf.kernel()   
        # 振動数計算のためのHessianオブジェクトを生成
        hessian = rks.Hessian(mf)
        # 振動数計算を実行
        freqs = hessian.kernel()

        if freqs is None:
            raise ValueError("Hessian計算が失敗し、freqsがNoneになりました。")

        # 振動数計算の結果を解析
        freq_info = thermo.harmonic_analysis(mf.mol, freqs)
    except Exception as e:
        print(f"振動数計算中にエラーが発生しました: {e}")
        return

    # 振動数計算の結果を表示する関数
    def display_vibrational_frequencies(freq_info):
        """
        振動数計算の結果を確認し、表示します。
        """
        if freq_info["freq_error"] == 0:
            print("振動数計算は正常に完了しました。")
        else:
            print("振動数計算結果に虚振動がありました。")

        print("{:>10} {:>20} {:>15}".format("Index", "Freq (cm^-1)", "Type"))
        for i, (freq, wavenumber) in enumerate(zip(freq_info["freq_au"], freq_info["freq_wavenumber"])):
            if freq.imag != 0:  # 虚数部分がある場合
                freq_type = "Imaginary"
                formatted_freq = f"{freq.real:.2f}+{freq.imag:.2f}j"
            else:  # 虚数部分がない場合
                freq_type = "Real"
                formatted_freq = f"{freq.real:.2f}+0.00j"
            
            print("{:>10} {:>20} {:>15}".format(i + 1, formatted_freq, freq_type))

    # 結果を表示
    display_vibrational_frequencies(freq_info)

def plot_nmr_spectrum(chemical_shifts, peak_width=0.2):
    """
    化学シフトデータを基にNMRスペクトルをプロットする関数。

    parameters:
        chemical_shifts (list): 計算された化学シフトのリスト。
        peak_width (float): ピークの幅。
    """
    # Plotlyがインストールされていない場合は以下の行のコメントを外してインストールしてください
    # !pip install plotly

    # 色の設定
    nmr_colors = {
        '1H': 'blue',
        '13C': 'green',
        '17O': 'red',
        # 必要に応じて追加
    }

    fig = go.Figure()

    for shift_data in chemical_shifts:
        nmr_type = shift_data['type']
        shift = shift_data['chemical_shift']
        color = nmr_colors.get(nmr_type, 'black')

        # ガウシアンピークの生成
        x = np.linspace(shift - 5*peak_width, shift + 5*peak_width, 1000)
        y = np.exp(-((x - shift)**2) / (2 * (peak_width / 2)**2))

        fig.add_trace(go.Scatter(
            x=x,
            y=y,
            mode='lines',
            line=dict(color=color, width=2),
            name=f'{nmr_type}: Atom {shift_data["atom_index"]}'
        ))

        # ピーク中心に縦線を追加
        fig.add_trace(go.Scatter(
            x=[shift, shift],
            y=[0, max(y)],
            mode='lines',
            line=dict(color=color, width=1),
            showlegend=False
        ))

    # レイアウト設定
    fig.update_layout(
        title='NMR Spectrum',
        xaxis=dict(title='Chemical Shift (ppm)', autorange='reversed'),
        yaxis=dict(title='Intensity (a.u.)'),
        legend_title='NMR Types',
        template='plotly_white'
    )

    fig.show()

def plot_molecule_py3Dmol(mol):
    """
    py3Dmolを使用して分子構造を3Dで可視化し、各原子に番号を表示する関数。
    """
    # XYZフォーマットの生成
    symbols = [mol.atom_symbol(i) for i in range(mol.natm)]
    coords = mol.atom_coords(unit="Angstrom")
    xyz = f"{mol.natm}\\n"
    xyz += "Generated by PySCF and py3Dmol\\n"
    for i, (sym, coord) in enumerate(zip(symbols, coords), start=1):
        x, y, z = coord
        xyz += f"{sym} {x:.6f} {y:.6f} {z:.6f}\\n"

    # py3Dmolビューの作成
    view = py3Dmol.view(width=800, height=600)
    view.addModel(xyz, 'xyz')
    view.setStyle({"stick": {"radius": 0.05}, "sphere": {"radius": 0.5}})  # スタイルをスティックモデルに設定

    # 各原子に番号をラベルとして追加
    for i, (sym, coord) in enumerate(zip(symbols, coords), start=1):
        x, y, z = coord
        view.addLabel(
            str(i),
            {
                'position': {'x': x, 'y': y, 'z': z},
                'backgroundColor': 'white',
                'fontColor': 'black',
                'fontSize': 14,
                'fontWeight': 'bold',
                'font': 'Arial',
                'showBackground': True
            }
        )

    view.zoomTo()
    view.show()

def main():
    try:
        # 入力方法の選択
        # 'coords', 'smiles', 'name'
        input_method = 'coords'  # ここで 'coords', 'smiles', 'name' を選択
        # - 'smiles' の場合: SMILES文字列
        # - 'name' の場合: 分子名
        # molecule_input = "Aceton"     
        # - 'coords' の場合: PySCFのatom定義文字列
        molecule_input = '''
          O        0.0000070000      1.4152910000      0.0000000000                 
          C       -0.0000030000     -0.6180500000     -1.2920750000                 
          C       -0.0000030000     -0.6180500000      1.2920750000                 
          C       -0.0000030000      0.1738500000      0.0000000000                 
          H       -0.0000270000      0.0632760000     -2.1447380000                 
          H       -0.8820880000     -1.2688680000     -1.3458520000                 
          H        0.8821140000     -1.2688230000     -1.3458740000                 
          H       -0.8820880000     -1.2688680000      1.3458520000                 
          H       -0.0000270000      0.0632760000      2.1447380000                 
          H        0.8821140000     -1.2688230000      1.3458740000                 
         '''

        # 分子の定義
        mol = define_molecule(
            input_method=input_method,
            molecule_input=molecule_input,
            unit='Angstrom',
            basis='6-31G(d)',
            charge=0,
            spin=0
        )

        # 構造最適化(GPU対応)
        mf = optimize_geometry(mol, path, xc='B3LYP', max_cycle=300, basis_opt='6-31G(d)')

        if mf is None:
            print("構造最適化が正常に完了しませんでした。NMR計算を中止します。")
            return

        # 振動数計算
        calculate_vibrational_frequencies(path, xc='B3LYP')

        # NMR遮蔽定数の計算
        chemical_shifts = calculate_nmr_shielding(
            mol,
            method='HF',              # 'HF' または 'DFT' を選択
            xc='B3LYP',                # DFTを選択した場合に使用
            selected_nmr_types=['1H', '13C'],
            scaling_factors={'1H': 1.0, '13C': 1.0},  # 必要に応じて値を調整
            TMS_shifts={'1H': 32.62, '13C': 200.31},      # 実際の値に置き換えてください
            basis_nmr="6-31G(d)"
        )

        if not chemical_shifts:
            print("NMR遮蔽定数の計算が完了しませんでした。")
            return

        # NMRスペクトルのプロット
        plot_nmr_spectrum(chemical_shifts, peak_width=0.2)

        # 分子構造の可視化
        plot_molecule_py3Dmol(mol)

    except Exception as e:
        print(f"スクリプトの実行中にエラーが発生しました: {e}")

if __name__ == "__main__":
    main()

以下は、コード内の各ステップの簡単な解説です。

1. 入力方法の選択

  • input_method: 分子の定義方法を指定。
    • 'coords': 直接座標(XYZ形式など)を指定。
    • 'smiles': SMILES形式で分子を定義。
    • 'name': 分子名を指定して、データベースから取得。
  • molecule_input: 分子データそのもの。
    • ここでは座標データを用いてアセトン分子を定義しています。

2. 分子の定義

  • define_molecule関数: 分子オブジェクトをPySCF形式で生成。
    • basis: 計算で使用する基底関数系(例: 6-31G(d))。
    • charge: 分子の電荷。
    • spin: 未対電子数(例: 中性分子では通常0)。

3. 構造最適化

  • 目的: 分子の構造をエネルギー的に最安定化する。
  • optimize_geometry関数:
    • xc: 汎関数(例: B3LYP)。
    • max_cycle: 構造最適化の最大反復回数。
    • basis_opt: 最適化で使用する基底関数。
  • 最適化が失敗した場合、スクリプトは処理を中断します。

4. 振動数計算

  • 目的: 虚振動の有無を確認し、分子が安定構造にあることを確かめる。
  • calculate_vibrational_frequencies関数:
    • 振動モードの計算とエネルギー状態の確認を実施。

5. NMR遮蔽定数の計算

  • 目的: 各原子における遮蔽定数(Shielding Constant)と化学シフト(Chemical Shift)を計算。
  • 主要パラメータ:
    • method: HFまたはDFTを選択。
    • selected_nmr_types: 計算する核種を指定(例: 1H, 13C)。
    • scaling_factors: 実験値に合わせるためのスケーリング係数。
    • TMS_shifts: 基準物質TMSの遮蔽定数。
    • basis_nmr: NMR計算で使用する基底関数。
  • 計算が失敗した場合、処理を中断。

6. NMRスペクトルのプロット

  • 目的: 計算された化学シフトを基にNMRスペクトルを描画。
  • peak_width: ピークの幅を調整するためのパラメータ。

7. 分子構造の可視化

  • 目的: 分子の構造を3Dで可視化。
  • ツール: py3Dmolライブラリを使用

6. 実際の分子へのNMR計算適用例


6.1.1 光学異性体のスペクトル比較

キラル分子の異性体を分子構造で与え、NMR計算を行うと、理想的には同じ化学シフト結果が得られます。しかし、溶媒や立体効果により微妙に化学シフトが異なる場合もあります。計算上は同じ構造エネルギー面上にある場合でも、相対的な化学シフトを比較することで、異性体の差を検討することが可能です。

たとえば(R)-体、(S)-体のアミノ酸などを計算した際に、特定のプロトンの化学シフトがわずかに異なるケースも報告されています。そのような小さな差分も、計算を駆使することで確認できます。

6.1.2 絶対値が合わなくても相対値で分別が可能

計算で得られる化学シフト値は、基底関数や計算レベル、スケーリングファクターにより誤差が生じることがあります。しかし、研究現場では絶対的なppm値よりも相対的な位置関係が重要視されるケースが多々あります。実験値との比較には適切なスケーリングを施す、あるいは差分を評価するなどの工夫によって、構造推定や異性体の同定に有用な情報が得られるでしょう。

特に、複数の分子や官能基が類似の環境にある場合、隣り合うピーク同士の相対的なずれが重要になることも多いです。

おわりに


本記事では、PySCFを用いたNMRスペクトル計算の手順とその理論背景について解説しました。NMRスペクトルは、有機合成や化学物質の同定において欠かせない手法であり、実験データと理論計算を組み合わせることで、分子構造の詳細な解析が可能となります。是非試してみて下さい!

参考文献


  1. PySCF 公式サイト: https://pyscf.org/

  2. Q. Sun et al., “Recent developments in the PySCF program package”, J. Chem. Phys., 153, 024109 (2020). DOI: 10.1063/5.0006074

  3. Q. Sun et al., “PySCF: the Python-based simulations of chemistry framework”, WIREs Comput. Mol. Sci., 8, e1340 (2018). DOI: 10.1002/wcms.1340

  4. U. Ekström et al., “Arbitrary-order density functional response theory from automatic differentiation”, J. Chem. Theory Comput., 6, 1971 (2010). DOI: 10.1021/ct100117s

  5. S. Lehtola et al., “Recent developments in libxc — A comprehensive library of functionals for density functional theory”, SoftwareX, 7, 1 (2018). DOI: 10.1016/j.softx.2017.11.002

  6. PySCF GPU4PySCFリポジトリ: https://github.com/pyscf/gpu4pyscf/blob/master/README.md

  7. M. Xie and A. Yu. Sokolov, “Exploring GPU Acceleration for PySCF with GPU4PySCF”, arXiv preprint (2024). DOI: 10.48550/arXiv.2404.09452

  8. M. Xie and A. Yu. Sokolov, “GPU-accelerated Hartree-Fock calculations in PySCF”, arXiv preprint (2024). DOI: 10.48550/arXiv.2407.09700

  9. 産業技術総合研究所 NMRスペクトルデータベース: https://sdbs.db.aist.go.jp/CNmrSpectralView.aspx?imgdir=cds&fname=CDS03412&sdbsno=319

  10. 西長亨、本田康, 『有機化学のための 量子化学計算入門: Gaussianの基本と有効活用のヒント』, 裳華房 (2022) P. 114-119.


コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です