4457 words
22 minutes
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる（ヘテロDiels-Alder反応）
最終更新：2025-05-03
概要
本記事では、自作モジュール(MultiOptPy)で、1,3-ブタジエンと二酸化硫黄が関わるヘテロDiels-Alder反応の遷移状態構造を算出してみる。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。
今回は、初めに低い計算レベル(HF/3-21G)で計算した後に、高い計算レベル(HF/6-31G)で遷移状態構造を算出する。この手順のメリットの1つは、全ての過程を高い計算レベルで行うよりも計算コストを削減することができる点にある。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
ヘテロDiels-Alder反応:
- J. Am. Chem. Soc. 1998, 120, 50, 13276–13277
- J. Am. Chem. Soc. 1992, 114, 23, 9210–9211
使用した自作モジュールMultiOptPyのバージョン
1.8.7e
環境
WSL2 (Ubuntu-22.04)
手順
1. 初期構造の準備
モデル反応系として、1,3-ブタジエン1分子と二酸化硫黄1分子を含む以下の構造をxyzファイル形式で用意した。今回は名前をhetero_diels_alder_rxn.xyzとした。
以下の構造を次のような条件で構造最適化を行い、初期構造を用意した。
python optmain.py hetero_diels_alder_rxn.xyz -ma 50 1-10 11-13 -func hf -bs 3-21g
-ma 50 1-10 11-13は原子のラベル番号1から10までで構成されるフラグメントと11-13までで構成されるフラグメントのペアに、50kJ/molの活性化障壁を超えうる人工力ポテンシャルをそれぞれ付加することを示す。(詳しくはAFIR法[https://doi.org/10.1063/1.3457903]を参照)→目的とする反応に該当する結合の組み換えが起らない程度に1,3-ブタジエンと二酸化硫黄を接近させ、後に行う構造最適化の時間短縮をするために行った。
-bsは基底関数の指定（特に指定しなければ6-31G(d)が使われる。）今回はPople系の基底関数である3-21Gを使用。-funcは計算手法の指定（Hartree-Fock(HF), post-HF(MP2等), DFTの汎関数）が可能。（特に指定しなければB3LYPが使われる。）今回はHartree-Fock法を使用。
<初期構造を求めるために用意した構造>
13
C -1.16894511 0.01272347 0.02114536
H -0.64133082 -0.91802344 0.03642152
H -2.23887748 0.02138948 0.02948964
C -0.48207648 1.18057761 -0.00877504
H -1.00969078 2.11132448 -0.02405296
C 1.05782616 1.16810497 -0.02078763
C 1.72607186 -0.01072508 -0.00143726
H 2.79600423 -0.01939110 -0.00978145
H 1.60014422 2.09018580 -0.04441092
H 1.18375382 -0.93280586 0.02218874
S 0.40256407 0.97428323 2.78675893
O -0.15489339 2.54849452 2.78675893
O 2.07256407 0.97428323 2.78675893
実行した結果、以下の構造が得られた(hetero_diels_alder_rxn_optimized.xyzより得られる)。こちらをhetero_diels_alder_rxn.xyzとリネームして初期構造として用いた。
13
OptimizedStructure
C -1.748629973135 -0.706675950511 -1.048075491396
H -1.309588006219 -1.524779806432 -1.587069211946
H -2.819038214545 -0.699281599984 -0.975181437192
C -1.018331499811 0.258355972136 -0.522847220341
H -1.478899819117 1.091322659580 -0.026224530667
C 0.455743162198 0.312656137245 -0.577331925397
C 1.260579232935 -0.727541507514 -0.421550110625
H 2.326683573678 -0.624763852271 -0.478802790143
H 0.876592379982 1.290456903924 -0.725313734954
H 0.878424777411 -1.711683893465 -0.222993513915
S 0.557700444608 0.543130140841 2.160513691443
O 0.009061560446 1.925110520229 1.791038884053
O 2.009702381571 0.573694276221 2.633837391079
2. 遷移状態構造最適化のための初期構造の算出
※あらかじめpythonで利用できる電子状態計算ソフトウェアのpsi4(導入方法は、psi4公式ページ[https://psicode.org/installs/v191/] を参照)が使える環境を用意しておく。
a. 以下のURLにアクセスし、Source codeをダウンロードする。
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.8.7e.zip
unzip v1.8.7e.zip
cd MultiOptPy-1.8.7e
b. 初期構造をカレントディレクトリにhetero_diels_alder_rxn_rearrangement.xyzとして保存し、以下を実行する。
python optmain.py hetero_diels_alder_rxn.xyz -ma 150 1 11 150 7 11 -bs 3-21g -func hf -opt rsirfo_fsb
-ma 150 1 11 150 7 11は原子のラベル番号1と11のペアと原子のラベル番号7と11のペアに、150kJ/molの活性化障壁を超えうる人工力ポテンシャルをそれぞれ付加することを示す。(詳しくはAFIR法[https://doi.org/10.1063/1.3457903]を参照)-bsは基底関数の指定（特に指定しなければ6-31G(d)が使われる。）今回はPople系の基底関数である3-21Gを使用。-funcは計算手法の指定（Hartree-Fock(HF), post-HF(MP2等), DFTの汎関数）が可能。（特に指定しなければB3LYPが使われる。）今回はHartree-Fock法を使用。-opt rsirfo_fsbは構造最適化につかう最適化アルゴリズムの変更を示す。デフォルトで指定しているものよりも効率が良いので今回はこちらを指定した。
これを実行すると、計算レベルHF/3-21Gで先ほどの人工力ポテンシャルをかけた上で初期構造を構造最適化する。
結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。
正常終了していれば、このディレクトリ中に、hetero_diels_alder_rxn_traj.xyzが存在するので、これをコピーして、MultiOptPy-1.8.7eディレクトリに置く。
hetero_diels_alder_rxn_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このhetero_diels_alder_rxn_traj.xyzは次のNEB計算に使用する。
※hetero_diels_alder_rxn_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
（この人工力ポテンシャルを加えて行った構造最適化の結果はhetero_diels_alder_rxn_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。）
c. hetero_diels_alder_rxn_traj.xyzを初期パスとして、NEB法で経路の緩和を行う。
NEB法を用いることで、先ほど得られたhetero_diels_alder_rxn_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
LUP法と呼ばれるNEB法とよく似た方法があるが、本モジュールではNEB法を用いている。
python nebmain.py hetero_diels_alder_rxn_traj.xyz -ns 10 -spng -nd 0.30 -bs 6-31g -func hf
-spngは経路の緩和途中の各ノードのエネルギーや勾配のRMSを可視化する。-nd 0.30はノード間の距離を0.30Åとして初期パスを作成することを示す。-ns 10は10回分NEB法による経路の緩和を行うことを示す。-bsは基底関数の指定（特に指定しなければ6-31G(d)が使われる。）今回はPople系の基底関数である6-31Gを使用。→ここで、低い計算レベル(HF/3-21G)でを用意したパスを、高い計算レベルでパスの緩和を行う。これにより、高い計算レベル向けの遷移状態構造最適化のための初期構造の算出を行うことができる。これは、前の過程にて、目的とする遷移状態付近を通る近似反応経路を得るために試行錯誤することとなった際に、効率性の面から都合が良い。
-funcは計算手法の指定（Hartree-Fock(HF), post-HF(MP2等), DFTの汎関数）が可能。（特に指定しなければB3LYPが使われる。）今回はHartree-Fock法を使用。
(参考)
d. 初期構造の決定及び遷移状態構造の計算
MultiOptPy-1.8.7eと同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。このディレクトリ内のplot_energy_9.pngを開き、エネルギーが極大値を示すノード番号を確認する。
<パスの緩和後の各ノードのエネルギー一覧(単位
-701.797547885856,-701.7962139171517,-701.7812875070348,-701.778949456086,-701.7761906562428,-701.7687589944584,-701.7624849696018,-701.7626488278896,-701.7636197096543,-701.7615922742301,-701.7616209350683,-701.7591564699749,-701.7542015618326,-701.7541846157952,-701.7592045202412,-701.7581350738144,-701.7592578704257,-701.760697378149,-701.7482687301608,-701.7383954577864,-701.7382089481032,-701.7294447225203,-701.7411870019764,-701.7538177224916,-701.7588472307489,-701.7703996015691,-701.775198885782,-701.782172031723,-701.7865698193236,-701.7912147056113,-701.7953583903235,-701.7980482738085,-701.8007612242224,-701.8034370664958,-701.8052066438805,-701.8062875606162,-701.8074924839188,-701.8088366452839
私が実行した環境では、21番である。 path_ITR_10_hetero_diels_alder_rxnディレクトリ内のhetero_diels_alder_rxn_traj_21.xyzを可視化する。
hetero_diels_alder_rxn_traj_21.xyz
13
0 1
C -1.281262998076 -0.524772566354 -0.041539347252
H -0.840345804380 -1.500484168573 -0.110454874245
H -2.365106775185 -0.562120359860 0.057310510593
C -0.841966178638 0.521423841934 -1.007305207182
H -1.534397602631 1.265077659429 -1.314305190669
C 0.474892683059 0.679009573428 -1.073112577047
C 1.174863163256 -0.255690851795 -0.111121532355
H 2.226833168601 -0.125772252746 -0.018767958085
H 0.940810975420 1.554388521010 -1.458044551080
H 0.948364102967 -1.267025881353 -0.432965288777
S 0.133691649484 -0.110247587860 1.286806397457
O 0.033265006107 1.415422903620 1.841177576334
O 0.930358610016 -1.089208830880 2.382322042307
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
ただし、hetero_diels_alder_rxn_traj_21.xyzで遷移状態構造最適化を行ったところ、今回はうまくいかなかったため、hetero_diels_alder_rxn_traj_21.xyz付近の勾配のRMS値が極小値を示す構造（hetero_diels_alder_rxn_traj_19.xyz）を遷移状態構造を求めるための初期構造とした。
<パスの緩和後の各ノードの勾配のRMS値一覧(単位
0.0018959512303263358,0.0024597119751811262,0.012351803287516429,0.010985035200775449,0.00887262514899807,0.010154719898884447,0.010713841393081177,0.00992097901262472,0.00919272281049582,0.010838339296992702,0.011020652555659644,0.010717692983918847,0.011468308119941248,0.0115086852161414,0.008458537395259296,0.008121651739174694,0.00952329993226271,0.00843045495279052,0.008817656241721445,0.006930959121922963,0.007385172802435283,0.019592813510850633,0.012457251355933499,0.011573952748862325,0.013518754011860418,0.007919019173458983,0.00935517106965961,0.005633199584006061,0.005193020816799143,0.004901919966772464,0.004438440378183392,0.004708262684103144,0.004645054067670604,0.0037448520697715803,0.0036734584061398623,0.0033657023018593612,0.0031889913050963934,0.0016949065027171194
hetero_diels_alder_rxn_traj_19.xyz
13
0 1
C -1.425463645832 -0.557053789613 -0.238316739281
H -0.941188042636 -1.504835588518 -0.135168890051
H -2.487122553972 -0.536816382943 -0.066710070354
C -0.883205358249 0.481194005408 -0.982482319686
H -1.504507972913 1.311700100760 -1.251763182033
C 0.496506425477 0.618911775141 -1.034931211354
C 1.270916714997 -0.300442274058 -0.306669063542
H 2.308509102096 -0.118315091812 -0.125086203924
H 0.924378112104 1.557151193928 -1.339459672439
H 1.004009279748 -1.329491080018 -0.423216309000
S 0.221034931886 -0.106266277446 1.473650839120
O -0.044346066854 1.425327818177 1.884038404601
O 1.060479074149 -0.941064409006 2.546114417943
hetero_diels_alder_rxn_traj_19.xyzをMultiOptPy-1.8.7eと同じディレクトリ内にコピーする。
そして、以下を実行する。
python optmain.py hetero_diels_alder_rxn_traj_19.xyz -opt rsirfo_bofill -order 1 -fc 5 -bs 6-31g -func hf -tcc
-opt rsirfo_bofillは遷移状態構造の最適化向けのoptimizerを指定することを意味する。-order 1は一次の鞍点を求めることを指定する。（デフォルトだと極小値を求めるようになっている。）-fc 5は5回の反復回数当たり1回正確なへシアンを計算することを指定する。-tccは構造最適化の収束条件をデフォルトよりも厳しくすることを示す。
実行して得られた正確な遷移状態構造を以下に示す。
(実行して得られた正確な遷移状態構造はhetero_diels_alder_rxn_traj_19_optimized.xyzとして保存されている。)
13
OptimizedStructure
C -1.506538518986 -0.548097638401 -0.277384342262
H -1.067341513411 -1.513851443962 -0.135306136723
H -2.539145978388 -0.469958017929 0.003592281796
C -0.920575753540 0.436857838281 -1.068751661369
H -1.541803032297 1.201155612283 -1.491309963461
C 0.446872351938 0.604667649086 -1.072243411080
C 1.280034732310 -0.186524992342 -0.216347478453
H 2.288026023548 0.156985642945 -0.067968219981
H 0.858024291374 1.482722770488 -1.530416300319
H 1.225437957431 -1.259017555165 -0.316316394179
S 0.246732504424 -0.098770906325 1.555466345385
O 0.356240156513 1.420417023568 2.098713758054
O 0.874036779085 -1.226585982528 2.518271522592
検証
今回求めた遷移状態構造が適当なものか確かめるためにPsi4に実装されている機能を用いて検証する。使用したPsi4のバージョンは1.9.1である。
具体的には、今回求めた正確な遷移状態構造を初期構造としてもう一度遷移状態構造最適化を行う。最適化が反復回数が１回だけであれば、自作モジュールで正確に停留点が求められていることになる。
後に、得られた遷移状態構造に対して基準振動数解析を行う。そして、虚振動数が１つだけであることを確認する。これにより、正確な遷移状態構造が求められていることが確認する。
最後に、虚振動方向を見て、遷移状態が反応系と生成系をつないでいるかを推測する。
以下のpythonコードを実行する。
こちらをクリックしてコードの詳細を表示
import psi4
import numpy as np
import os
import re
bohr2ang = 0.52918
# Set memory allocation and number of threads
psi4.set_memory('2 GB')
psi4.set_num_threads(nthread=8)
# Initialize output file
psi4.core.set_output_file('ts_vib_output.dat', False)
# Create output directory
if not os.path.exists('vib_results'):
os.makedirs('vib_results')
# Define molecule
mol = psi4.geometry("""
0 1
C -1.506538518986 -0.548097638401 -0.277384342262
H -1.067341513411 -1.513851443962 -0.135306136723
H -2.539145978388 -0.469958017929 0.003592281796
C -0.920575753540 0.436857838281 -1.068751661369
H -1.541803032297 1.201155612283 -1.491309963461
C 0.446872351938 0.604667649086 -1.072243411080
C 1.280034732310 -0.186524992342 -0.216347478453
H 2.288026023548 0.156985642945 -0.067968219981
H 0.858024291374 1.482722770488 -1.530416300319
H 1.225437957431 -1.259017555165 -0.316316394179
S 0.246732504424 -0.098770906325 1.555466345385
O 0.356240156513 1.420417023568 2.098713758054
O 0.874036779085 -1.226585982528 2.518271522592
""")
# Set calculation parameters
psi4.set_options({
'basis': '6-31g',
'reference': 'rhf',
'scf_type': 'direct',
'G_CONVERGENCE': 'gau',
'opt_type': 'ts'
})
# Perform transition state optimization
print("\n\n=== Starting Transition State Optimization ===\n\n")
ts_energy, ts_wfn = psi4.optimize('hf', return_wfn=True)
# Save the optimized TS structure
mol.save_xyz_file("vib_results/ts_optimized.xyz", True)
# Frequency calculation to get vibrational modes
print("\n\n=== Running Frequency Analysis ===\n\n")
freq_energy, freq_wfn = psi4.frequency('hf', return_wfn=True)
# Extract the atom symbols
atom_symbols = []
for atom in range(mol.natom()):
atom_symbols.append(mol.symbol(atom))
# Get the equilibrium coordinates
coords = np.array([[mol.x(i), mol.y(i), mol.z(i)] for i in range(mol.natom())])
# Function to extract vibration modes from the output file using the improved regex
def extract_vibration_modes():
with open('ts_vib_output.dat', 'r') as f:
content = f.read()
# Find all frequency blocks (each block contains 3 columns of data)
freq_blocks = re.findall(
r'Freq \[cm\^-1\]\s+([\d\.]+i?)\s+([\d\.]+i?)\s+([\d\.]+i?)',
content
)
# Find all displacement vector blocks
vector_blocks = re.findall(
r'-{10,}.*?\n((?:\s+\d+\s+[A-Z][a-z]?\s+(?:-?\d+\.\d+\s+){9,}.*?\n)+)',
content, re.DOTALL
)
# Process frequencies to find imaginary ones
imag_freqs = []
all_freqs = []
# Flatten frequency data from blocks
for block in freq_blocks:
for freq in block:
all_freqs.append(freq)
if 'i' in freq:
imag_freqs.append(freq)
if len(imag_freqs) == 1:
print("Transition state structure was found!")
elif len(imag_freq_value) == 0:
print("Munimum 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.")
# If no data found, return None
if not all_freqs or not vector_blocks:
return None, [], False
# Find the first imaginary frequency (if any)
imag_freq_idx = -1
for i, freq in enumerate(all_freqs):
if 'i' in freq:
imag_freq_idx = i
break
# If no imaginary frequency, use the first normal mode
if imag_freq_idx == -1:
imag_freq_idx = 0
is_imaginary = False
else:
is_imaginary = True
# Determine which block and column in block contains our target frequency
block_idx = imag_freq_idx // 3
col_idx = imag_freq_idx % 3
if block_idx >= len(vector_blocks):
return None, [], is_imaginary
# Parse the displacement vectors for our target frequency
vectors = []
lines = vector_blocks[block_idx].strip().split('\n')
for line in lines:
parts = line.strip().split()
if len(parts) >= 10: # Atom index, symbol, and 3 sets of x,y,z
# Extract the correct column based on which frequency we want
start_idx = 2 + (col_idx * 3)
x = float(parts[start_idx])
y = float(parts[start_idx + 1])
z = float(parts[start_idx + 2])
vectors.append([x, y, z])
# Extract just the numeric part of the frequency
freq_value = float(all_freqs[imag_freq_idx].rstrip('i'))
return freq_value, vectors, is_imaginary
# Extract the vibration mode
imag_freq_value, mode_vectors, is_imaginary = extract_vibration_modes()
# If extraction failed, use the example data
if not mode_vectors or len(mode_vectors) != mol.natom():
print("Using example vibration data because extraction failed or was incomplete.")
exit()
else:
print(f"Successfully extracted vibration data: frequency = {imag_freq_value}{' (imaginary)' if is_imaginary else ''}")
# Save frequency information to file
with open('vib_results/frequency_analysis.txt', 'w') as f:
f.write("=== Vibrational Frequency Analysis ===\n\n")
# Extract the frequency section from the output file
with open('ts_vib_output.dat', 'r') as vib_file:
vib_content = vib_file.read()
# Find and extract all frequency blocks with their data
freq_sections = re.findall(
r'(Vibration\s+\d+.*?\n\s*Freq \[cm\^-1\].*?\n.*?-{10,}.*?\n(?:\s+\d+\s+[A-Z][a-z]?.*?\n)+)',
vib_content, re.DOTALL
)
if freq_sections:
f.write("Extracted frequency data:\n\n")
for section in freq_sections:
f.write(section + "\n")
# Find imaginary frequencies
imag_freqs = re.findall(r'Freq \[cm\^-1\].*?([\d\.]+i)', vib_content)
if imag_freqs:
f.write(f"\nFound {len(imag_freqs)} imaginary frequencies:\n")
for i, freq in enumerate(imag_freqs):
f.write(f"{i+1}: {freq} cm^-1\n")
f.write(f"\nPrimary imaginary frequency: {imag_freqs[0]} cm^-1\n")
else:
f.write("\nNo imaginary frequencies found in output. This may not be a transition state.\n")
# Frequency string representation (with 'i' if imaginary)
freq_str = f"{imag_freq_value:.4f}{'i' if is_imaginary else ''}"
print(f"Creating animation for {'imaginary ' if is_imaginary else ''}frequency: {freq_str} cm^-1")
# Create animation of the normal mode
num_frames = 20 # Number of frames for animation
amplitude = 1.0 # Amplitude factor for visualization
# 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 = coords + amplitude * phase * np.array(mode_vectors)
# Write frame to XYZ file
f.write(f"{mol.natom()}\n")
f.write(f"{'Imaginary' if is_imaginary 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)):
coord = coord * bohr2ang
f.write(f"{symbol} {coord[0]:12.8f} {coord[1]:12.8f} {coord[2]:12.8f}\n")
# Create additional animation with larger amplitude to visualize the motion better
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 = coords + enhanced_amplitude * phase * np.array(mode_vectors)
f.write(f"{mol.natom()}\n")
f.write(f"Amplified {'imaginary' if is_imaginary 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)):
coord = coord * bohr2ang
f.write(f"{symbol} {coord[0]:12.8f} {coord[1]:12.8f} {coord[2]:12.8f}\n")
print("\n\n=== Calculation Complete ===")
print("Results saved to vib_results/ directory:")
print("- Transition state structure: vib_results/ts_optimized.xyz")
print("- Frequency analysis: vib_results/frequency_analysis.txt")
print(f"- {'Imaginary' if 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")
実行して得られたts_irc_output.datを見る。
初めに遷移状態構造最適化についてである。収束条件の確認を行う。以下は素の出力例である。
==> Convergence Check <==
Measures of convergence in internal coordinates in au.
Criteria marked as inactive (o), active & met (*), and active & unmet ( ).
----------------------------------------------------------------------------------------------
Step Total Energy Delta E Max Force RMS Force Max Disp RMS Disp
----------------------------------------------------------------------------------------------
Convergence Criteria o 4.50e-04 * 3.00e-04 * 1.80e-03 * 1.20e-03 *
----------------------------------------------------------------------------------------------
5 -701.75041455 1.22e-06 o 4.40e-05 * 1.49e-05 * 5.49e-04 * 1.77e-04 * ~
----------------------------------------------------------------------------------------------
残念ながら1回の計算で収束条件は満たされてはいなかったが、今回自作モジュールで求めた構造で収束条件（電子状態計算ソフトウェアGaussianのデフォルトの基準と同等）を現実的な反復回数で満たせることがわかる。
次に、基準振動数解析の結果を確認する。以下に結果の一例を示す。 frequency_analysis.txtで以下の内容を確認できる。
=== Vibrational Frequency Analysis ===
Extracted frequency data:
Vibration 1 8 9
Freq [cm^-1] 471.6399i 88.3200 150.7091
Irrep A A A
Reduced mass [u] 4.3679 7.0103 5.8702
Force const [mDyne/A] -0.5725 0.0322 0.0786
Turning point v=0 [a0] 0.0000 0.4410 0.3689
RMS dev v=0 [a0 u^1/2] 0.0000 0.8256 0.6320
IR activ [km/mol] 46.7642 2.5789 72.1157
Char temp [K] 0.0000 127.0729 216.8369
----------------------------------------------------------------------------------
1 C 0.16 0.04 0.31 0.04 -0.12 0.03 -0.09 -0.08 -0.21
2 H -0.02 -0.12 -0.15 0.14 -0.05 0.24 -0.01 -0.02 -0.08
3 H 0.24 0.15 0.55 0.02 -0.18 -0.02 -0.13 -0.15 -0.35
4 C 0.07 0.03 -0.00 -0.03 -0.21 -0.13 -0.03 -0.02 -0.11
5 H -0.03 -0.14 -0.15 -0.10 -0.36 -0.30 -0.02 -0.09 -0.25
6 C -0.05 0.02 -0.04 -0.04 -0.07 -0.13 -0.06 0.06 0.16
7 C -0.09 0.04 0.26 0.03 0.16 -0.01 -0.16 0.08 0.36
8 H -0.12 0.13 0.28 -0.01 0.28 -0.04 -0.21 0.17 0.49
9 H 0.07 -0.17 -0.31 -0.13 -0.11 -0.28 -0.00 0.05 0.19
10 H 0.18 0.05 -0.10 0.13 0.13 0.17 0.04 0.09 0.08
11 S 0.01 -0.03 -0.16 -0.01 0.16 0.03 -0.05 -0.02 0.03
12 O -0.05 0.03 -0.05 0.23 0.02 0.32 0.25 0.08 -0.17
13 O -0.06 -0.07 -0.03 -0.21 -0.13 -0.18 0.12 -0.06 -0.05
...(snip)...
Vibration 37 38 39
Freq [cm^-1] 3388.1276 3403.5680 3456.4803
Irrep A A A
Reduced mass [u] 1.1006 1.0991 1.1179
Force const [mDyne/A] 7.4437 7.5016 7.8692
Turning point v=0 [a0] 0.1797 0.1794 0.1765
RMS dev v=0 [a0 u^1/2] 0.1333 0.1330 0.1320
IR activ [km/mol] 0.8287 4.3925 0.2389
Char temp [K] 4874.7613 4896.9766 4973.1055
----------------------------------------------------------------------------------
1 C 0.00 -0.01 0.00 -0.00 0.01 -0.00 0.08 -0.06 -0.01
2 H -0.05 0.09 -0.02 0.04 -0.09 0.01 -0.31 0.66 -0.11
3 H 0.00 0.00 -0.00 0.00 -0.00 0.00 -0.63 0.05 0.18
4 C -0.03 0.04 -0.02 0.03 -0.03 0.02 0.01 -0.01 0.00
5 H 0.39 -0.48 0.27 -0.31 0.37 -0.21 -0.07 0.09 -0.05
6 C -0.00 -0.01 0.01 -0.03 -0.06 0.03 -0.00 -0.00 0.00
7 C -0.05 -0.04 -0.01 -0.03 -0.02 -0.01 -0.00 0.00 -0.00
8 H 0.59 0.20 0.09 0.29 0.10 0.04 0.00 0.00 0.00
9 H 0.05 0.12 -0.07 0.30 0.63 -0.33 0.01 0.01 -0.01
10 H 0.01 0.34 0.03 0.00 0.16 0.01 -0.00 -0.00 -0.00
11 S 0.00 -0.00 -0.00 0.00 -0.00 -0.00 -0.00 0.00 0.00
12 O 0.00 0.00 0.00 0.00 0.00 0.00 0.00 -0.00 0.00
13 O -0.00 -0.00 0.00 0.00 -0.00 0.00 -0.00 -0.00 -0.00
Found 1 imaginary frequencies:
1: 471.6399i cm^-1
Primary imaginary frequency: 471.6399i cm^-1
Primary imaginary frequency: 471.6399i cm^-1より、１個の虚振動が確認された。よって、今回自作モジュール(MultiOptPy)で求めた遷移状態構造は、正確な遷移状態構造であることがわかる。
そして、虚振動方向をvibration_mode_animation.xyzかvibration_mode_animation_amplified.xyz(振動が強調されている。)から確認すると、生成系と反応系をつなぐように振動していることが確認できる。（IRC計算をしたほうが正確だが、今回は単純な反応なので虚振動方向の確認で十分と判断した。）
※アニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
終わりに
自作モジュールで、遷移状態構造を算出する手順を説明した。さらに、求めた遷移状態構造を既存のモジュールを用いて正しい遷移状態構造なのかを検証した。
参考
- https://github.com/psi4/psi4 (Psi4のgithubのレポジトリ)
- https://github.com/ss0832/MultiOptPy (自作モジュールMultiOptPyのレポジトリ)
- https://avogadro.cc/ (Avogadro、分子構造可視化ツール)
- J. Am. Chem. Soc. 1998, 120, 50, 13276–13277 (ヘテロDiels-Alder反応)
- J. Am. Chem. Soc. 1992, 114, 23, 9210–9211 (ヘテロDiels-Alder反応)
- The Journal of Chemical Physics 2010, 132 (24), 241102.
- The Journal of Chemical Physics 1991, 94 (1), 751–760.
- In Classical and Quantum Dynamics in Condensed Phase Simulations; WORLD SCIENTIFIC: LERICI, Villa Marigola, 1998; pp 385–404.
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる（ヘテロDiels-Alder反応）
https://ss0832.github.io/posts/20250503_mop_usage_5/