【PySCF】IRスペクトル・熱力学的特性の解析【Pythonで始める量子化学計算】

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

from google.colab import drive
drive.mount('/content/drive')

import os

# 計算結果を保存するパスを設定
folder_path = "/content/drive/MyDrive/pyscf_IR_thermo"
base_filename = "anticlinal_conformation"

# フルパスの設定
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}")

# 例としてチェックポイントファイル名の設定
checkpoint_file = path + ".chk"
print(f"Checkpoint file will be saved as: {checkpoint_file}")



from pyscf import gto

mol = gto.M(
    atom='''
    F       -2.9513875514     -1.0928382404      0.2476531628
    F       -5.6179262375      1.2415516978     -0.2476528977
    C       -3.5272115798      0.1385618493     -0.0055750514
    C       -5.0421022105      0.0101516046      0.0055753029
    H       -3.1869047581      0.5023934144     -0.9980363055
    H       -3.2040823737      0.8590234395      0.7752476523
    H       -5.3652314210     -0.7103099785     -0.7752474058
    H       -5.3824090296     -0.3536799681      0.9980365550
    ''',
    basis='6-31G(d)',
    charge=0,
    spin=0,
    symmetry=True,
    verbose=4
)

from pyscf import dft

mf = dft.RKS(mol)
mf.chkfile = path + '.chk'
mf.xc = 'B3LYP'
mf.max_cycle = 150  # SCF計算の最大反復回数

energies = []  # 各ステップのエネルギーを保存するリスト
geometries = []  # 各ステップの分子構造を保存するリスト

# コールバック関数の定義
def cb(envs):
    mf = envs['g_scanner'].base
    if not hasattr(mf, 'e_tot'):
        print("エネルギーが更新されていません。")
        return
    energies.append(mf.e_tot)
    geometries.append(mol.atom_coords(unit="ANG"))
    print(f"ステップ {len(energies)}: エネルギー = {mf.e_tot} Hartree")

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

from pyscf.geomopt import geometric_solver
import matplotlib.pyplot as plt

# 構造最適化
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()

# ジオメトリーの保存
for i, coords in enumerate(geometries):
    filename = f"{path}_geometry_step_{i+1}.xyz"
    with open(filename, 'w') as f:
        f.write(f"{mol.natm}\\\\n")
        f.write(f"Step {i+1} Energy: {energies[i]}\\\\n")
        for j in range(mol.natm):
            symbol = mol.atom_pure_symbol(j)
            x, y, z = coords[j]
            f.write(f"{symbol} {x:.6f} {y:.6f} {z:.6f}\\\\n")

print("All geometries have been saved.")



import ipywidgets
from ipywidgets import interact
import py3Dmol

# モデルのXYZ座標をファイルから読み込む
def read_xyz(file_path):
    with open(file_path, 'r') as file:
        return file.read()

# py3Dmolでのビューアを設定する関数
def show_structure(xyz_content, style):
    view = py3Dmol.view(width=800, height=400)
    view.addModel(xyz_content, 'xyz')
    if style == 'ball and stick':
        view.setStyle({'stick':{}, 'sphere':{'radius': 0.5}})
    elif style == 'stick':
        view.setStyle({'stick': {}})
    elif style == 'line':
        view.setStyle({'line': {'linewidth': 5}})
    elif style == 'sphere':
        view.setStyle({'sphere': {}})
    view.setBackgroundColor('white')
    view.zoomTo()
    return view.show()

# ファイルパスを設定
filename_template = f"{path}_geometry_step_"

# 最適化ステップの数に応じたインタラクション
def style_selector(idx, s):
    filename = f"{filename_template}{idx}.xyz"
    xyz_content = read_xyz(filename)
    show_structure(xyz_content, s)

# インタラクションウィジェットの設定
max_steps = i + 1  # 最適化ステップの総数
interact(style_selector,
         idx=ipywidgets.IntSlider(min=1, max=max_steps, step=1, value=max_steps, description='Step'),
         s=ipywidgets.Dropdown(
             options=['ball and stick', 'stick','line', 'sphere'],
             value='ball and stick',
             description='Style:'))



import matplotlib.pyplot as plt
import sys
from pyscf import gto, dft, lib
from pyscf.hessian import rks  # 振動数計算のためのモジュール
from pyscf.prop import infrared
from pyscf.hessian import thermo

# チェックポイントファイルから分子情報をロード
# 'path'はチェックポイントファイルのパスです。チェックポイントファイルには、前回の計算結果が保存されています。
mol = lib.chkfile.load_mol(path + '.chk')

# DFT計算の設定
# チェックポイントファイルを指定して、前回の計算結果を引き継ぎます。
mf = dft.RKS(mol)
mf.chkfile = path + '.chk'  # チェックポイントファイルの指定
mf.xc = 'b3lyp'  # 使用する交換・相関汎関数の指定

# DFT計算の実行
# 分子の電子構造を計算します。
mf.kernel()

# IRスペクトル解析の実行
# infraredモジュールを使用して赤外スペクトルを計算します。
mf_ir = infrared.rks.Infrared(mf).run()

# IRスペクトルの概要を表示
# 計算結果のサマリーを表示します。
mf_ir.summary()

# 熱力学的計算の実行
# thermoモジュールを使用して、計算した振動周波数を基に熱力学的性質を評価します。
thermo.dump_thermo(mol, thermo.thermo(mf, mf_ir.vib_dict["freq_au"], 298.15, 101325))

# IRスペクトルのプロット
# 振動数のスケーリングファクターとして0.960を適用
# 詳細はCCCBDBの推奨を参照（<https://cccbdb.nist.gov/vibscalejust.asp）>
fig, ax, ax2 = mf_ir.plot_ir(w=100, scale=0.960)

# プロットのタイトルを設定
# プロットの見出しを設定し、詳細を含めます。
ax.set_title(r"Infrared Spectra of Methylformate Molecule ($\\\\\\\\mathrm{CH2FCH2F}$)" "\\\\\\\\n"
              r"B3LYP/6-31G(d) with scaling factor 0.960 and Lorentz broadening with $\\\\\\\\mathrm{FWHW = 100 cm^{-1}}")

# レイアウトの調整とプロットの表示
# プロットのレイアウトを調整して表示します。
fig.tight_layout()
fig.show()



import os

# 計算結果を保存するパスを設定
folder_path = "/content/drive/MyDrive/pyscf_IR_thermo"
base_filename2 = "gauche_conformation"

# フルパスの設定
path2 = os.path.join(folder_path, base_filename2)

# フォルダが存在しない場合は作成
if not os.path.exists(folder_path):
    os.makedirs(folder_path)
    print(f"Created folder: {folder_path}")

# 例としてチェックポイントファイル名の設定
checkpoint_file2 = path2 + ".chk"
print(f"Checkpoint file will be saved as: {checkpoint_file2}")



import matplotlib.pyplot as plt
import sys
from pyscf import gto, dft, lib
from pyscf.geomopt import geometric_solver
from pyscf.hessian import rks  # 振動数計算のためのモジュール
from pyscf.hessian import thermo

# 分子の定義
mol = gto.M(
    atom='''
    F       -2.9646670035      1.3527250959      0.1781513936
    F       -4.6216795179     -0.7305791160     -0.6508177073
    C       -4.2423816872      1.5905363784     -0.2946175127
    C       -5.1302460519      0.3878960539     -0.0160002927
    H       -4.1984219564      1.7820726269     -1.3875645139
    H       -4.6566700057      2.4872569027      0.2125627625
    H       -5.1771329542      0.1985823306      1.0772122356
    H       -6.1561483328      0.5846824999     -0.3925288715
    ''',
    basis='6-31G(d)',
    verbose=4
)

# DFT計算の設定
mf = dft.RKS(mol)
mf.chkfile = checkpoint_file2
mf.xc = 'B3LYP'
mf.max_cycle = 150  # SCF計算の最大反復回数

energies = []  # 各ステップのエネルギーを保存するリスト
geometries = []  # 各ステップの分子構造を保存するリスト

# コールバック関数の定義
def cb(envs):
    mf = envs['g_scanner'].base
    if not hasattr(mf, 'e_tot'):
        print("エネルギーが更新されていません。")
        return
    energies.append(mf.e_tot)
    geometries.append(mol.atom_coords(unit="ANG"))
    print(f"ステップ {len(energies)}: エネルギー = {mf.e_tot} Hartree")

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

# 構造最適化
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()

# ジオメトリーの保存
for i, coords in enumerate(geometries):
    filename = f"{path2}_geometry_step_{i+1}.xyz"
    with open(filename, 'w') as f:
        f.write(f"{mol.natm}\\\\n")
        f.write(f"Step {i+1} Energy: {energies[i]}\\\\n")
        for j in range(mol.natm):
            symbol = mol.atom_pure_symbol(j)
            x, y, z = coords[j]
            f.write(f"{symbol} {x:.6f} {y:.6f} {z:.6f}\\\\n")

print("All geometries have been saved.")



import ipywidgets
from ipywidgets import interact
import py3Dmol

# モデルのXYZ座標をファイルから読み込む
def read_xyz(file_path):
    with open(file_path, 'r') as file:
        return file.read()

# py3Dmolでのビューアを設定する関数
def show_structure(xyz_content, style):
    view = py3Dmol.view(width=800, height=400)
    view.addModel(xyz_content, 'xyz')
    if style == 'ball and stick':
        view.setStyle({'stick':{}, 'sphere':{'radius': 0.5}})
    elif style == 'stick':
        view.setStyle({'stick': {}})
    elif style == 'line':
        view.setStyle({'line': {'linewidth': 5}})
    elif style == 'sphere':
        view.setStyle({'sphere': {}})
    view.setBackgroundColor('white')
    view.zoomTo()
    return view.show()

# ファイルパスを設定
filename_template = f"{path2}_geometry_step_"

# 最適化ステップの数に応じたインタラクション
def style_selector(idx, s):
    filename = f"{filename_template}{idx}.xyz"
    xyz_content = read_xyz(filename)
    show_structure(xyz_content, s)

# インタラクションウィジェットの設定
max_steps = i + 1  # 最適化ステップの総数
interact(style_selector,
         idx=ipywidgets.IntSlider(min=1, max=max_steps, step=1, value=max_steps, description='Step'),
         s=ipywidgets.Dropdown(
             options=['ball and stick', 'stick','line', 'sphere'],
             value='ball and stick',
             description='Style:'))



import pandas as pd
import os
from pyscf import gto, dft, lib
from pyscf.hessian import rks, thermo

# Hartree to eV conversion factor
HARTREE_TO_EV = 27.2114

def calculate_thermo_properties(path):
    # Checkpointファイルから分子情報をロード
    mol = lib.chkfile.load_mol(path + '.chk')

    # DFT計算の設定
    mf = dft.RKS(mol)
    mf.chkfile = path + '.chk'
    mf.xc = 'b3lyp'

    # DFT計算
    mf.kernel()

    # 振動数計算の実行
    hessian = rks.Hessian(mf)
    freqs = hessian.kernel()

    # Thermochemistry analysis at 298.15 K and 1 atmospheric pressure
    freq_info = thermo.harmonic_analysis(mf.mol, freqs)
    thermo_info = thermo.thermo(mf, freq_info['freq_au'], 298.15, 101325)

    # 熱化学特性を辞書形式で返す
    return {
        "File Name": os.path.basename(path),
        "Imaginary Freq (cm^-1)": freq_info['freq_error'],
        "Temperature (K)": thermo_info['temperature'][0],
        "Pressure (atm)": thermo_info['pressure'][0] / 101325,  # Pa to atmに変換
        "Frequencies scaled by": 1.0000,
        "Electronic Energy (EE) (Hartree)": thermo_info['E0'][0],
        "Zero-point Energy Correction (Hartree)": thermo_info['ZPE'][0],
        "Thermal Correction to Energy (Hartree)": thermo_info['E_vib'][0],
        "Thermal Correction to Enthalpy (Hartree)": thermo_info['H_vib'][0],
        "Thermal Correction to Free Energy (Hartree)": thermo_info['G_vib'][0],
        "EE + Zero-point Energy (Hartree)": thermo_info['E_0K'][0],
        "EE + Thermal Energy Correction (Hartree)": thermo_info['E_tot'][0],
        "EE + Thermal Enthalpy Correction (Hartree)": thermo_info['H_tot'][0],
        "EE + Thermal Free Energy Correction (Hartree)": thermo_info['G_tot'][0],
        "E (Thermal) (Hartree)": thermo_info['E_0K'][0],  # Thermal energyとして仮定
        "Heat Capacity (Cv) (J/mol·K)": thermo_info['Cv_tot'][0],
        "Entropy (S) (J/mol·K)": thermo_info['S_tot'][0]
    }

# 最初の構造の計算結果
thermo_properties_1 = calculate_thermo_properties(path)

# 2番目の構造の計算結果
thermo_properties_2 = calculate_thermo_properties(path2)

# データフレームに追加
df = pd.DataFrame([thermo_properties_1, thermo_properties_2])

# 最小エネルギーを持つ構造を見つける
min_energy = df["EE + Zero-point Energy (Hartree)"].min()

# 相対エネルギーを計算し、eV単位に変換
df["Relative Energy (eV)"] = (df["EE + Zero-point Energy (Hartree)"] - min_energy) * HARTREE_TO_EV

# 行と列を入れ替える
transposed_df_with_units = df.set_index("File Name").T

# 結果を表示
transposed_df_with_units