4007 words
20 minutes
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる（Huisgen付加環化反応）
最終更新：2025-05-15
概要
本記事では、自作モジュール(MultiOptPy)で、Huisgen付加環化反応の遷移状態構造を算出してみる。今回は電子状態計算ソフトウェアにPySCFを使用する。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
使用した自作モジュールMultiOptPyのバージョン
1.9.3e
環境
WSL2 (Ubuntu-22.04)
手順
1. 初期構造の準備
モデル反応系として、アセチレン1分子とアジドメタン1分子を含む以下の構造をxyzファイル形式で用意した。今回は名前をhuisgen_cycloaddition.xyzとした。
以下の構造を次のような条件で構造最適化を行い、初期構造を用意した。
python optmain.py huisgen_cycloaddition.xyz -opt rsirfo_fsb -ma 10 1 7 10 3 6 -pyscf -spin 0 -func hf -bs 6-31g
-ma 10 1 7 10 3 6は原子のラベル番号1と7のペアと原子のラベル番号3と6のペアに、それぞれ10kJ/molの活性化障壁を超えうる人工力ポテンシャルをそれぞれ付加することを示す。(詳しくはAFIR法[https://doi.org/10.1063/1.3457903]を参照) →2分子同士が拡散して離れないようにするための処理-opt rsirfo_fsbは構造最適化向けの高性能なオプティマイザーのアルゴリズムを使用することを示す。-pyscfは電子状態計算ソフトウェアの指定をPsi4からPySCFに変えるためのオプションである。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。-bsは基底関数の指定（特に指定しなければ6-31G(d)が使われる。）今回はPople系の基底関数である6-31Gを使用。-funcは計算手法の指定（Hartree-Fock(HF), post-HF(MP2等), DFTの汎関数）が可能。（特に指定しなければB3LYPが使われる。）今回はHartree-Fock法を使用。
<初期構造を求めるために用意した構造>
11
0 1
C -4.56210186 -1.95063691 -0.00000000
H -5.63210186 -1.95063691 -0.00000000
C -3.36090186 -1.95063691 -0.00000000
H -2.29090186 -1.95063691 0.00000000
N -3.45445669 -0.72366265 2.71172498
N -2.29796537 -0.40967800 2.42577837
N -4.47952838 -1.00196823 2.96517722
C -1.79692216 0.93168293 2.75835399
H -2.29347099 1.29052745 3.63560531
H -0.74336515 0.88420565 2.93908455
H -1.98922542 1.59668045 1.94245096
実行した結果、以下の構造が得られた(huisgen_cycloaddition_optimized.xyzより得られる)。こちらをhuisgen_cycloaddition.xyzとリネームして初期構造として用いた。
11
OptimizedStructure
C -1.453088716194 -1.494818146477 -1.669789058563
H -2.497243006950 -1.606956942337 -1.599410668329
C -0.267739344394 -1.371533340436 -1.754986136683
H 0.777119509536 -1.262445369017 -1.805843047547
N -0.636312161622 -0.121575457176 0.966154008869
N 0.583763849698 -0.012933691774 0.740495897863
N -1.712275591303 -0.397229528696 1.075695887443
C 1.174001461666 1.324917680996 0.973085750064
H 0.916563735878 1.728202871427 1.945068103613
H 2.239855576787 1.186013693374 0.924494146154
H 0.875354686899 2.028358230117 0.205035117116
2. 遷移状態構造最適化のための初期構造の算出
※あらかじめpythonで利用できる電子状態計算ソフトウェアのPySCF(導入方法はpipを使えば可能である。)が使える環境を用意しておく。
a. 以下のURLにアクセスし、Source codeをダウンロードする。
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.9.3e.zip
unzip v1.9.3e.zip
cd MultiOptPy-1.9.3e
b. 初期構造をカレントディレクトリにhuisgen_cycloaddition.xyzとして保存し、以下を実行する。
python optmain.py huisgen_cycloaddition.xyz -opt rsirfo_fsb -ma 100 1 7 100 3 6 -pyscf -spin 0 -func hf -bs 6-31g
これを実行すると、計算レベルHF/6-31Gで指定した人工力ポテンシャルをかけた上で初期構造を構造最適化する。
結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。
正常終了していれば、このディレクトリ中に、huisgen_cycloaddition_traj.xyzが存在するので、これをコピーして、MultiOptPy-1.9.3eディレクトリに置く。
huisgen_cycloaddition_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このhuisgen_cycloaddition_traj.xyzは次のNEB計算に使用する。
※huisgen_cycloaddition_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
（この人工力ポテンシャルを加えて行った構造最適化の結果はhuisgen_cycloaddition_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。）
c. huisgen_cycloaddition_traj.xyzを初期パスとして、NEB法で経路の緩和を行う。
NEB法を用いることで、先ほど得られたhuisgen_cycloaddition_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
LUP法と呼ばれるNEB法とよく似た方法があるが、本モジュールではNEB法を用いている。
python nebmain.py huisgen_cycloaddition_traj.xyz -ns 10 -spng -nd 0.2 -pyscf -spin 0 -elec 0 -func hf -bs 6-31g
-spngは経路の緩和途中の各ノードのエネルギーや勾配のRMSをmatplotlibを用いて可視化するためのオプションである。-nd 0.20はノード間の距離を0.20Åとして初期パスを作成することを示す。-ns 10は10回分NEB法による経路の緩和を行うことを示す。-bsは基底関数の指定（特に指定しなければ6-31G(d)が使われる。）今回はPople系の基底関数である6-31Gを使用。-funcは計算手法の指定（Hartree-Fock(HF), post-HF(MP2等), DFTの汎関数）が可能。（特に指定しなければB3LYPが使われる。）今回はHartree-Fock法を使用。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。-elecは形式電荷の指定である。今回は中性なので0を指定。
(参考)
d. 初期構造の決定及び遷移状態構造の計算
MultiOptPy-1.9.3eと同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。このディレクトリ内のplot_energy_9.pngを開き、エネルギーが極大値を示すノード番号を確認する。
<パスの緩和後の各ノードのエネルギー一覧(単位
-279.530816571615,-279.53048639516646,-279.5253964706841,-279.517896430136,-279.51104884576023,-279.51113077121,-279.5073649481894,-279.5070457836142,-279.5030869007446,-279.4966162154964,-279.4873383699872,-279.48201720313546,-279.49676596454185,-279.5377457317554,-279.5835942899924,-279.6032526956598,-279.62212186470396,-279.6317879543121,-279.64055774938834,-279.64439211899,-279.64677102010637,-279.6495646354455,-279.6512268283734,-279.6517985998191,-279.65337278577897,-279.6544909787659
私が実行した環境では、11番である。 path_ITR_10_huisgen_cycloadditionディレクトリ内のhuisgen_cycloaddition_traj_11.xyzを可視化する。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
huisgen_cycloaddition_traj_11.xyz
11
0 0
C -1.294858561950 -1.392140099099 -1.174975694423
H -2.251510671553 -1.817669756533 -1.316125572791
C -0.124710991267 -1.051143825215 -1.333966066195
H 0.842591217554 -0.915154381462 -1.724814498373
N -0.898382388563 0.022066822538 0.859240646109
N 0.278767072408 0.087967233265 0.427394014195
N -1.887095447861 -0.488418011175 0.534773759291
C 1.155874344668 1.195066927147 0.860256653936
H 1.065404320521 1.413381428212 1.916026543194
H 2.160507936031 0.855208873064 0.668725968058
H 0.953413170011 2.090834789258 0.283464246998
構造が壊れていないので、これを遷移状態を求めるための初期構造とする。
huisgen_cycloaddition_traj_19.xyzをMultiOptPy-1.8.7eと同じディレクトリ内にコピーする。
そして、以下を実行する。
python optmain.py huisgen_cycloaddition_traj_11.xyz -opt rsirfo_bofill -order 1 -fc 5 -pyscf -spin 0 -func hf -bs 6-31g -freq
-opt rsirfo_bofillは遷移状態構造の最適化向けのoptimizerを指定することを意味する。-order 1は一次の鞍点を求めることを指定する。（デフォルトだと極小値を求めるようになっている。）-fc 5は5回の反復回数当たり1回正確なへシアンを計算することを指定する。-freqは収束条件を満たした後に基準振動解析を行うことを示す。（自前で実装しているため、あくまで目安として使用することを推奨する。各振動モードをvibration_animation内のxyzファイルで可視化できる。）
実行して得られた正確な遷移状態構造を以下に示す。
(実行して得られた正確な遷移状態構造はhuisgen_cycloaddition_traj_11_optimized.xyzとして保存されている。)
11
OptimizedStructure
C -1.257101249762 -1.381814617903 -1.193094358983
H -2.167038144260 -1.823615927101 -1.488692112061
C -0.068150317325 -1.081678382101 -1.228462227372
H 0.942248592575 -1.032315290648 -1.522722895500
N -0.988407842383 0.121612736133 0.874179251732
N 0.221475979432 0.078853725804 0.493601652012
N -1.986127578401 -0.357749159322 0.566903556415
C 1.136368001094 1.173451995307 0.820998098890
H 1.341526437626 1.212055417570 1.882504408653
H 2.057789448227 0.954917747235 0.304744376361
H 0.767416673178 2.136281755026 0.490040249854
検証
今回求めた遷移状態構造が適当なものか確かめるためにPySCFに実装されている機能を用いて検証する。使用したPySCFのバージョンは2.9.0である。
具体的には、今回求めた正確な遷移状態構造を初期構造としてもう一度遷移状態構造最適化を行う。最適化が反復回数が1回だけであれば、自作モジュールで正確に停留点が求められていることになる。
後に、得られた遷移状態構造に対して基準振動数解析を行う。そして、虚振動数が１つだけであることを確認する。これにより、正確な遷移状態構造が求められていることが確認する。
最後に、虚振動方向を見て、遷移状態が反応系と生成系をつないでいるかを推測する。
以下のpythonコードを実行する。
こちらをクリックしてコードの詳細を表示
import numpy as np
import os
import matplotlib.pyplot as plt
from pyscf import gto, scf, hessian
from pyscf.hessian import thermo
from pyscf.geomopt import geometric_solver
from pyscf.data.elements import _symbol
# Constants
bohr2ang = 0.52918
# Create output directory
if not os.path.exists('vib_results'):
os.makedirs('vib_results')
# Define molecule
mol = gto.Mole()
mol.atom = '''
C -1.257101249762 -1.381814617903 -1.193094358983
H -2.167038144260 -1.823615927101 -1.488692112061
C -0.068150317325 -1.081678382101 -1.228462227372
H 0.942248592575 -1.032315290648 -1.522722895500
N -0.988407842383 0.121612736133 0.874179251732
N 0.221475979432 0.078853725804 0.493601652012
N -1.986127578401 -0.357749159322 0.566903556415
C 1.136368001094 1.173451995307 0.820998098890
H 1.341526437626 1.212055417570 1.882504408653
H 2.057789448227 0.954917747235 0.304744376361
H 0.767416673178 2.136281755026 0.490040249854
'''
# Using 3-21G basis as requested
mol.basis = '6-31G'
mol.charge = 0
mol.spin = 0
mol.build()
# Function to save xyz file
def save_xyz_file(filename, mol, coords=None):
if coords is None:
coords = mol.atom_coords() * bohr2ang
atom_symbols = [_symbol(mol.atom_charge(i)) for i in range(mol.natm)]
with open(filename, 'w') as f:
f.write(f"{mol.natm}\n")
f.write("Generated by PySCF\n")
for i, (symbol, coord) in enumerate(zip(atom_symbols, coords)):
f.write(f"{symbol} {coord[0]:12.8f} {coord[1]:12.8f} {coord[2]:12.8f}\n")
# Get initial coordinates for reference
initial_coords = mol.atom_coords().copy()
save_xyz_file("vib_results/initial_structure.xyz", mol)
# Setup HF calculation with the 3-21G basis
mf = scf.RHF(mol)
mf.conv_tol = 1e-9
mf.max_cycle = 200
print("\n\n=== Starting HF Calculation ===\n\n")
mf.kernel()
# Output file for logs
with open('ts_vib_output.dat', 'w') as f:
f.write(f"SCF Energy: {mf.e_tot}\n\n")
print("\n\n=== Starting Geometry Optimization ===\n\n")
try:
# GeometRIC optimization returns the optimized molecule object
mol_opt = geometric_solver.optimize(mf, maxsteps=10)
# Get optimized coordinates
ts_coords = mol_opt.atom_coords()
save_xyz_file("vib_results/optimized_structure.xyz", mol_opt)
# Run a new SCF calculation on the optimized structure
mf_opt = scf.RHF(mol_opt)
mf_opt.kernel()
print("\n\n=== Running Vibrational Analysis ===\n\n")
# Calculate Hessian matrix
h = hessian.RHF(mf_opt)
hess_matrix = h.kernel()
# Using PySCF's built-in thermo module for vibrational analysis
# This handles all the mass-weighting and frequency calculation
thermo_results = thermo.harmonic_analysis(mol_opt, hess_matrix)
# Extract frequencies and modes
raw_frequencies = thermo_results['freq_wavenumber'] # might be complex
reduced_mass = thermo_results['reduced_mass']
modes = thermo_results['norm_mode']
# Convert complex frequencies to real numbers
# Imaginary frequencies (real part ≈ 0, imaginary part ≠ 0) become negative real values
# Regular frequencies get their real part
frequencies = np.zeros(len(raw_frequencies))
for i, freq in enumerate(raw_frequencies):
if np.iscomplexobj(freq):
real_part = np.real(freq)
imag_part = np.imag(freq)
# Check if this is a true imaginary frequency (real part ≈ 0)
if abs(real_part) < 1e-10 and abs(imag_part) > 1e-10:
# Convert pure imaginary frequency to negative real value
frequencies[i] = -abs(imag_part)
else:
# For complex numbers with significant real part, use the real part
frequencies[i] = real_part
else:
# Already a real number
frequencies[i] = freq
# Check for imaginary frequencies (now represented as negative values)
imag_freqs = [i for i in range(len(frequencies)) if frequencies[i] < 0]
is_imaginary = len(imag_freqs) > 0
# Report structure type based on imaginary frequencies
if is_imaginary:
if len(imag_freqs) == 1:
print("Transition state structure was found!")
else:
num_imag_freqs = len(imag_freqs)
if num_imag_freqs == 2:
print("2nd-order saddle point was found.")
elif num_imag_freqs == 3:
print("3rd-order saddle point was found.")
else:
print(f"{num_imag_freqs}th-order saddle point was found.")
else:
print("Minimum structure was found!")
# Select mode for visualization (first imaginary mode if available, otherwise first vibrational mode)
# No need to skip first 6 modes as thermo module already handles this
mode_idx = imag_freqs[0] if is_imaginary else 0
if mode_idx >= len(modes):
mode_idx = 0 # Safety check
selected_mode = modes[mode_idx]
selected_freq = frequencies[mode_idx]
# Save frequency analysis to file
with open('vib_results/frequency_analysis.txt', 'w') as f:
f.write("=== Vibrational Frequency Analysis ===\n\n")
f.write(f"Total number of vibrational modes: {len(frequencies)}\n")
f.write(f"Number of imaginary frequencies: {len(imag_freqs)}\n\n")
f.write("Vibrational Frequencies (cm^-1):\n")
for i, (freq, mass) in enumerate(zip(frequencies, reduced_mass)):
# No need to skip modes as thermo module already handles this
f.write(f"Mode {i+1:3d}: {abs(freq):12.4f}{' (imaginary)' if freq < 0 else ''} cm^-1, ")
f.write(f"Reduced mass: {mass:.6f} amu\n")
if is_imaginary:
f.write(f"\nPrimary imaginary frequency: {abs(selected_freq):.4f}i cm^-1\n")
f.write(f"Mode index: {mode_idx}\n\n")
f.write("Displacement vectors for imaginary mode:\n")
atom_symbols = [_symbol(mol_opt.atom_charge(i)) for i in range(mol_opt.natm)]
for i, (symbol, vec) in enumerate(zip(atom_symbols, selected_mode.reshape(-1, 3))):
f.write(f"{i+1:4d} {symbol:2s} {vec[0]:12.8f} {vec[1]:12.8f} {vec[2]:12.8f}\n")
# Frequency string representation (with 'i' if imaginary)
freq_str = f"{abs(selected_freq):.4f}{'i' if selected_freq < 0 else ''}"
print(f"Creating animation for {'imaginary ' if selected_freq < 0 else ''}frequency: {freq_str} cm^-1")
# Get atom symbols for animation
atom_symbols = [_symbol(mol_opt.atom_charge(i)) for i in range(mol_opt.natm)]
# Create animation of the normal mode
num_frames = 20 # Number of frames for animation
amplitude = 1.0 # Amplitude factor for visualization
# Reshape mode if necessary
if selected_mode.ndim > 2:
selected_mode = selected_mode.reshape(-1, 3)
# Create animation file
with open("vib_results/vibration_mode_animation.xyz", 'w') as f:
for frame in range(num_frames):
# Calculate phase factor for oscillation (-1 to 1)
phase = np.sin(2 * np.pi * frame / num_frames)
# Calculate displaced coordinates
displaced_coords = ts_coords + amplitude * phase * selected_mode
# Convert to Angstroms for XYZ file
displaced_coords_ang = displaced_coords * bohr2ang
# Write frame to XYZ file
f.write(f"{mol_opt.natm}\n")
f.write(f"{'Imaginary' if selected_freq < 0 else 'Normal'} frequency animation, frame {frame+1}/{num_frames}, " +
f"frequency = {freq_str} cm^-1, phase = {phase:.2f}\n")
for i, (symbol, coord) in enumerate(zip(atom_symbols, displaced_coords_ang)):
f.write(f"{symbol} {coord[0]:12.8f} {coord[1]:12.8f} {coord[2]:12.8f}\n")
# Create additional animation with larger amplitude for better visualization
with open("vib_results/vibration_mode_animation_amplified.xyz", 'w') as f:
enhanced_amplitude = 3.0 # Larger amplitude for better visualization
for frame in range(num_frames):
phase = np.sin(2 * np.pi * frame / num_frames)
displaced_coords = ts_coords + enhanced_amplitude * phase * selected_mode
displaced_coords_ang = displaced_coords * bohr2ang
f.write(f"{mol_opt.natm}\n")
f.write(f"Amplified {'imaginary' if selected_freq < 0 else 'normal'} frequency animation, " +
f"frame {frame+1}/{num_frames}, frequency = {freq_str} cm^-1, " +
f"phase = {phase:.2f}, amplitude factor = {enhanced_amplitude:.1f}\n")
for i, (symbol, coord) in enumerate(zip(atom_symbols, displaced_coords_ang)):
f.write(f"{symbol} {coord[0]:12.8f} {coord[1]:12.8f} {coord[2]:12.8f}\n")
except Exception as e:
print(f"Error during calculation: {str(e)}")
print("\n\n=== Calculation Complete ===")
print("Results saved to vib_results/ directory:")
print("- Optimized structure: vib_results/optimized_structure.xyz")
print("- Frequency analysis: vib_results/frequency_analysis.txt")
print(f"- {'Imaginary' if 'is_imaginary' in locals() and is_imaginary else 'Normal'} mode animation: vib_results/vibration_mode_animation.xyz")
print("- Amplified animation: vib_results/vibration_mode_animation_amplified.xyz")
print("Main output file: ts_vib_output.dat")
初めに遷移状態構造最適化についてである。収束条件の確認を行う。以下は素の出力例である。
cycle 2: E = -279.485868237 dE = 2.85657e-08 norm(grad) = 0.000312294
Step 1 : Displace = 3.583e-04/6.302e-04 (rms/max) Trust = 1.000e-01 (=) Grad = 9.416e-05/1.975e-04 (rms/max) E (change) = -279.4858682367 (+2.857e-08) Quality = -0.447
Hessian Eigenvalues: 2.30000e-02 2.30000e-02 2.30000e-02 ... 7.65773e-01 9.24503e-01 1.35440e+00
Converged! =D
#==========================================================================#
#| If this code has benefited your research, please support us by citing: |#
#| |#
#| Wang, L.-P.; Song, C.C. (2016) "Geometry optimization made simple with |#
#| translation and rotation coordinates", J. Chem, Phys. 144, 214108. |#
#| http://dx.doi.org/10.1063/1.4952956 |#
#==========================================================================#
Time elapsed since start of run_optimizer: 13.247 seconds
converged SCF energy = -279.485868237098
2回のcycleでgeometricの収束条件を満たしたことが判明した。
次に、基準振動数解析の結果を確認する。以下に結果の一例を示す。 frequency_analysis.txtで以下の内容を確認できる。
=== Vibrational Frequency Analysis ===
Total number of vibrational modes: 27
Number of imaginary frequencies: 1
Vibrational Frequencies (cm^-1):
Mode 1: 630.2511 (imaginary) cm^-1, Reduced mass: 7.848099 amu
Mode 2: 42.2555 cm^-1, Reduced mass: 3.232241 amu
Mode 3: 150.5706 cm^-1, Reduced mass: 1.094742 amu
Mode 4: 226.6007 cm^-1, Reduced mass: 3.804030 amu
Mode 5: 305.3443 cm^-1, Reduced mass: 3.262436 amu
Mode 6: 414.5251 cm^-1, Reduced mass: 4.318283 amu
Mode 7: 544.8131 cm^-1, Reduced mass: 4.685317 amu
Mode 8: 628.8389 cm^-1, Reduced mass: 8.407307 amu
Mode 9: 734.6766 cm^-1, Reduced mass: 4.080603 amu
Mode 10: 801.9040 cm^-1, Reduced mass: 1.662207 amu
Mode 11: 833.8143 cm^-1, Reduced mass: 1.175761 amu
Mode 12: 895.0725 cm^-1, Reduced mass: 1.517104 amu
Mode 13: 910.3777 cm^-1, Reduced mass: 2.388001 amu
Mode 14: 1013.6408 cm^-1, Reduced mass: 1.317149 amu
Mode 15: 1231.6900 cm^-1, Reduced mass: 1.259811 amu
Mode 16: 1279.7080 cm^-1, Reduced mass: 1.424727 amu
Mode 17: 1372.8221 cm^-1, Reduced mass: 12.353443 amu
Mode 18: 1611.7700 cm^-1, Reduced mass: 1.152234 amu
Mode 19: 1649.1428 cm^-1, Reduced mass: 1.049191 amu
Mode 20: 1659.9268 cm^-1, Reduced mass: 1.049767 amu
Mode 21: 1949.5982 cm^-1, Reduced mass: 6.138693 amu
Mode 22: 2044.8515 cm^-1, Reduced mass: 6.894354 amu
Mode 23: 3200.2696 cm^-1, Reduced mass: 1.034758 amu
Mode 24: 3271.7794 cm^-1, Reduced mass: 1.104050 amu
Mode 25: 3314.0688 cm^-1, Reduced mass: 1.099551 amu
Mode 26: 3600.2958 cm^-1, Reduced mass: 1.086480 amu
Mode 27: 3683.0029 cm^-1, Reduced mass: 1.183791 amu
Primary imaginary frequency: 630.2511i cm^-1
Mode index: 0
Displacement vectors for imaginary mode:
1 C -0.02609203 0.07431014 0.12252832
2 H 0.11347027 -0.06130903 -0.12620894
3 C 0.01792925 0.08769216 0.11769790
4 H -0.02467552 -0.06016761 -0.07784565
5 N -0.00979938 0.01607093 0.03787042
6 N -0.01697452 -0.03764779 -0.08942034
7 N 0.02994218 -0.06820185 -0.10668188
8 C -0.00047686 -0.03696346 -0.03449931
9 H -0.04769059 -0.07565963 -0.02539656
10 H 0.01961073 -0.01350505 -0.00822300
11 H -0.00179344 -0.03173081 -0.01493787
Mode 1: 630.2511 (imaginary) cm^-1, Reduced mass: 7.848099 amuより、１個の虚振動が確認された。よって、今回自作モジュール(MultiOptPy)で求めた遷移状態構造で、正確な遷移状態構造が求められることがわかる。
そして、虚振動方向をvibration_mode_animation.xyzかvibration_mode_animation_amplified.xyz(振動が強調されている。)から確認すると、生成系と反応系をつなぐように振動していることが確認できる。（IRC計算をしたほうが正確だが、今回は単純な反応なので虚振動方向の確認で十分と判断した。）
※アニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
終わりに
自作モジュールで、遷移状態構造を算出する手順を説明した。さらに、求めた遷移状態構造を既存のモジュールを用いて正しい遷移状態構造なのかを検証した。
参考
- https://github.com/pyscf/pyscf (PySCFのgithubのレポジトリ)
- https://github.com/ss0832/MultiOptPy (自作モジュールMultiOptPyのレポジトリ)
- https://avogadro.cc/ (Avogadro、分子構造可視化ツール)
- https://www.chem-station.com/odos/2010/01/-huisgen-cycloaddition.html (Huisgen付加環化反応)
- https://en.wikipedia.org/wiki/Azide-alkyne_Huisgen_cycloaddition (Huisgen付加環化反応)
- The Journal of Chemical Physics 2010, 132, 241102.
- The Journal of Chemical Physics 1991, 94, 751–760.
- In Classical and Quantum Dynamics in Condensed Phase Simulations; WORLD SCIENTIFIC: LERICI, Villa Marigola, 1998; pp 385–404.
- The Journal of Chemical Physics, 2020, 153, 024109.
- The Journal of Chemical Physics, 2022, 144, 214108.
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる（Huisgen付加環化反応）
https://ss0832.github.io/posts/20250515_mop_usage_6/