3636 words
18 minutes
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる（クライゼン転位）
最終更新：2025-05-03
概要
本記事では、自作モジュール(MultiOptPy)で、クライゼン転位の遷移状態構造を算出してみる。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy クライゼン転位：https://en.wikipedia.org/wiki/Claisen_rearrangement
使用した自作モジュールMultiOptPyのバージョン
1.8.6
環境
WSL2 (Ubuntu-22.04)
手順
1. 初期構造の準備
モデル反応系として、以下の構造をxyzファイル形式で用意した。今回は名前をclaisen_rearrangement.xyzとした。
(事前に構造最適化したものを使用すると失敗しにくい。)
14
C -1.897439922880 -0.613803649185 -1.412116540079
H -1.872193454868 -1.671977245468 -1.222919220472
H -2.850918253474 -0.203267508112 -1.684120813893
C -0.819063557695 0.134355553194 -1.321695098598
H -0.862171525022 1.191862303747 -1.508180908564
C 0.548953558862 -0.396337492618 -0.984226220678
H 0.493132775086 -1.432142442997 -0.669287742820
H 1.205971909925 -0.326448551480 -1.838925001221
O 1.192366041333 0.411155137871 0.027873340469
C 0.645474815019 0.384016111213 1.297624450891
H -0.389078365079 0.102802158635 1.349521531919
C 1.342363326529 0.711641413440 2.359675380422
H 2.372199194871 0.992795834752 2.276465141235
H 0.890403457393 0.715348377010 3.330311701388
2. 遷移状態構造最適化のための初期構造の算出
※あらかじめpythonで利用できる電子状態計算ソフトウェアのpsi4(導入方法は、psi4公式ページ[https://psicode.org/installs/v191/] を参照)が使える環境を用意しておく。
a. 以下のURLにアクセスし、Source codeをダウンロードする。
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.8.6.zip
unzip v1.8.6.zip
cd MultiOptPy-1.8.6
b. 初期構造をカレントディレクトリにclaisen_rearrangement.xyzとして保存し、以下を実行する。
python optmain.py claisen_rearrangement.xyz -ma 300 1 12 -bs 3-21g -func hf
-ma 300 1 12は原子のラベル番号1と12のペアに、300kJ/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法を使用。
これを実行すると、計算レベルHF/3-21Gで先ほどの人工力ポテンシャルをかけた上で初期構造を構造最適化する。
結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。
正常終了していれば、このディレクトリ中に、claisen_rearrangement_traj.xyzが存在するので、これをコピーして、MultiOptPy-1.8.6ディレクトリに置く。
claisen_rearrangement_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このclaisen_rearrangement_traj.xyzは次のNEB計算に使用する。
※claisen_rearrangement_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
（この人工力ポテンシャルを加えて行った構造最適化の結果はclaisen_rearrangement_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。）
c. claisen_rearrangement_traj.xyzを初期パスとして、NEB法で経路の緩和を行う。
NEB法を用いることで、先ほど得られたclaisen_rearrangement_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
LUP法と呼ばれるNEB法とよく似た方法があるが、本モジュールではNEB法を用いている。
python nebmain.py claisen_rearrangement_traj.xyz -ns 10 -spng -nd 0.20 -bs 3-21g -func hf
-spngは経路の緩和途中の各ノードのエネルギーや勾配のRMSを可視化する。-nd 0.20はノード間の距離を0.20Åとして初期パスを作成することを示す。-ns 10は10回分NEB法による経路の緩和を行うことを示す。-bsは基底関数の指定（特に指定しなければ6-31G(d)が使われる。）今回はPople系の基底関数である3-21Gを使用。-funcは計算手法の指定（Hartree-Fock(HF), post-HF(MP2等), DFTの汎関数）が可能。（特に指定しなければB3LYPが使われる。）今回はHartree-Fock法を使用。
(参考)
d. 初期構造の決定及び遷移状態構造の計算
MultiOptPy-1.8.6と同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。このディレクトリ内のplot_energy_9.pngを開き、エネルギーが極大値を示すノード番号を確認する。
私が実行した環境では、42番である。 これらの２つの構造を確認する。 path_ITR_10_claisen_rearrangementディレクトリ内のclaisen_rearrangement_traj_42.xyzを可視化する。
claisen_rearrangement_traj_42.xyz
14
0 1
C -1.168608436583 0.165657054855 0.025227450231
H -1.506653211082 -0.836798947662 0.209674814793
H -1.908428446912 0.881971058610 0.318756914605
C -0.427779922894 0.417142332234 -1.181245925517
H -0.347077772058 1.435769406913 -1.523655250009
C 0.603817692389 -0.467844537370 -1.448121769838
H 0.385626785707 -1.505556275161 -1.344828380077
H 1.303702400233 -0.241057339793 -2.226912921490
O 1.639818113250 -0.334160514164 -0.069874170466
C 0.969566828738 -0.438833249975 1.022219786055
H 0.927024262727 -1.389036267812 1.545183616128
C -0.036077190348 0.492469962408 1.340303756265
H 0.102244300964 1.520882573091 1.057691370749
H -0.537175404132 0.299394743824 2.275580708572
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
claisen_rearrangement_traj_42.xyzをMultiOptPy-1.8.6と同じディレクトリ内にコピーする。
そして、以下を実行する。
python optmain.py claisen_rearrangement_traj_42.xyz -opt ersprfo_bofill -order 1 -fc 5 -bs '6-31G(d)' -func B3LYP
-opt ersprfo_bofillは遷移状態構造の最適化向けのoptimizerを指定することを意味する。-order 1は一次の鞍点を求めることを指定する。（デフォルトだと極小値を求めるようになっている。）-fc 5は5回の反復回数当たり1回正確なへシアンを計算することを指定する。
実行して得られた正確な遷移状態構造を以下に示す。
(実行して得られた正確な遷移状態構造はclaisen_rearrangement_traj_42_optimized.xyzとして保存されている。)
14
OptimizedStructure
C -1.284333758764 0.159636856641 -0.122362479678
H -1.526352468674 -0.849151519024 0.152035696970
H -2.014758660597 0.886890139244 0.178183446013
C -0.464264924553 0.405219579374 -1.207685024376
H -0.378747831441 1.412480277379 -1.568245239189
C 0.559917493427 -0.487089497540 -1.493893284402
H 0.405496618009 -1.537553564631 -1.343637333875
H 1.284686751492 -0.246398827488 -2.243691608262
O 1.688723805355 -0.281092677304 -0.009561729332
C 0.981899472158 -0.428759059005 1.059912117017
H 0.931794632423 -1.412438216759 1.506772659987
C 0.077378650844 0.512318972482 1.488446615495
H 0.217637134265 1.537624381646 1.215135192695
H -0.479076913945 0.328313154986 2.388590970937
検証
今回求めた遷移状態構造が適当なものか確かめるために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.284333758764 0.159636856641 -0.122362479678
H -1.526352468674 -0.849151519024 0.152035696970
H -2.014758660597 0.886890139244 0.178183446013
C -0.464264924553 0.405219579374 -1.207685024376
H -0.378747831441 1.412480277379 -1.568245239189
C 0.559917493427 -0.487089497540 -1.493893284402
H 0.405496618009 -1.537553564631 -1.343637333875
H 1.284686751492 -0.246398827488 -2.243691608262
O 1.688723805355 -0.281092677304 -0.009561729332
C 0.981899472158 -0.428759059005 1.059912117017
H 0.931794632423 -1.412438216759 1.506772659987
C 0.077378650844 0.512318972482 1.488446615495
H 0.217637134265 1.537624381646 1.215135192695
H -0.479076913945 0.328313154986 2.388590970937
""")
# Set calculation parameters
psi4.set_options({
'basis': '3-21g',
'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 -267.23858631 -5.11e-06 o 7.00e-05 * 2.11e-05 * 6.43e-04 * 2.22e-04 * ~
----------------------------------------------------------------------------------------------
残念ながら1回の計算で収束条件は満たされてはいなかったが、今回自作モジュールで求めた構造で収束条件（電子状態計算ソフトウェアGaussianのデフォルトの基準と同等）を現実的な反復回数で満たせることがわかる。
次に、基準振動数解析の結果を確認する。以下に結果の一例を示す。 frequency_analysis.txtで以下の内容を確認できる。
=== Vibrational Frequency Analysis ===
Extracted frequency data:
Vibration 1 8 9
Freq [cm^-1] 763.3506i 204.0826 380.9595
Irrep A A A
Reduced mass [u] 10.4938 2.8146 3.9605
Force const [mDyne/A] -3.6027 0.0691 0.3387
Turning point v=0 [a0] 0.0000 0.4578 0.2825
RMS dev v=0 [a0 u^1/2] 0.0000 0.5431 0.3975
IR activ [km/mol] 85.6557 3.7257 24.7738
Char temp [K] 0.0000 293.6294 548.1159
----------------------------------------------------------------------------------
1 C 0.22 0.03 0.39 -0.01 0.16 -0.03 0.25 -0.01 0.20
2 H -0.17 -0.02 -0.08 -0.22 0.22 -0.00 0.22 0.01 0.25
3 H 0.01 -0.03 0.04 0.12 0.32 -0.09 0.23 -0.01 0.15
4 C 0.10 -0.03 -0.07 0.05 -0.03 -0.03 0.07 -0.09 0.02
5 H 0.04 -0.00 -0.01 0.17 -0.08 -0.14 0.00 -0.11 -0.04
6 C -0.33 0.01 -0.37 -0.05 -0.17 0.06 0.14 -0.03 -0.01
7 H 0.07 0.05 0.17 -0.13 -0.13 0.30 0.21 -0.06 -0.10
8 H -0.09 0.01 -0.14 -0.07 -0.39 -0.02 0.08 0.06 -0.05
9 O 0.17 0.01 0.32 0.02 0.20 -0.03 -0.19 0.09 -0.09
10 C 0.11 -0.04 -0.08 0.05 -0.01 -0.04 -0.21 -0.04 -0.11
11 H -0.00 0.00 0.04 0.16 -0.08 -0.18 -0.48 -0.14 -0.39
12 C -0.33 -0.00 -0.29 -0.06 -0.17 0.07 -0.02 0.06 0.04
13 H 0.15 0.06 0.13 -0.11 -0.10 0.27 0.18 0.07 0.18
14 H -0.12 0.01 -0.15 -0.07 -0.37 0.02 -0.10 0.03 -0.02
...(snip)...
Vibration 40 41 42
Freq [cm^-1] 3402.0034 3427.7793 3429.8397
Irrep A A A
Reduced mass [u] 1.1149 1.1128 1.1143
Force const [mDyne/A] 7.6024 7.7038 7.7232
Turning point v=0 [a0] 0.1782 0.1777 0.1775
RMS dev v=0 [a0 u^1/2] 0.1330 0.1325 0.1325
IR activ [km/mol] 4.3058 6.7614 6.0380
Char temp [K] 4894.7254 4931.8113 4934.7758
----------------------------------------------------------------------------------
1 C 0.03 -0.09 -0.00 0.00 -0.01 -0.00 0.00 -0.01 -0.00
2 H 0.15 0.61 -0.18 0.01 0.05 -0.01 0.01 0.04 -0.01
3 H -0.49 0.49 0.21 -0.04 0.04 0.02 -0.03 0.03 0.01
4 C -0.00 -0.02 0.01 0.00 -0.00 0.00 -0.00 -0.01 0.00
5 H 0.02 0.16 -0.06 0.00 0.02 -0.01 0.01 0.13 -0.04
6 C 0.00 0.00 -0.00 -0.00 -0.00 0.00 -0.06 -0.06 0.06
7 H -0.01 -0.03 0.00 0.00 0.02 -0.00 0.06 0.50 -0.07
8 H -0.05 -0.02 0.05 0.02 0.01 -0.02 0.58 0.19 -0.58
9 O -0.00 -0.00 0.00 0.00 0.00 -0.00 0.00 0.00 0.00
10 C 0.00 0.00 -0.00 -0.00 -0.00 0.00 0.00 0.00 -0.00
11 H 0.00 -0.00 0.00 0.00 0.06 -0.02 -0.00 -0.00 -0.00
12 C 0.00 0.01 -0.00 -0.03 -0.08 0.05 0.00 0.00 -0.00
13 H -0.00 -0.05 0.02 0.10 0.81 -0.22 -0.00 -0.03 0.01
14 H -0.03 -0.01 0.05 0.27 0.08 -0.43 -0.01 -0.00 0.02
Found 1 imaginary frequencies:
1: 763.3506i cm^-1
Primary imaginary frequency: 763.3506i cm^-1
Primary imaginary frequency: 763.3506i cm^-1より、１個の虚振動が確認された。よって、今回自作モジュール(MultiOptPy)で求めた遷移状態構造を用いて、Claisen転位の正確な遷移状態構造を求められたことがわかる。
そして、虚振動方向を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、分子構造可視化ツール)
- https://en.wikipedia.org/wiki/Claisen_rearrangement (クライゼン転位)
- https://www.chem-station.com/odos/2009/06/claisen-claisen-rearrangement.html (クライゼン転位, Chem-Stationの記事)
- 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モジュールで遷移状態構造を求めてみる（クライゼン転位）
https://ss0832.github.io/posts/20250503_mop_usage_4/