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

# 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

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

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
"""

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

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")

# 必要なライブラリのインポート
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()