Home
3453 words
17 minutes
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる
2025-04-19
Computational Chemistry
MultiOptPy
/
python

最終更新：2025-04-22

概要#

本記事では、自作モジュール(MultiOptPy)で、ヒドロホウ素化反応の遷移状態構造を算出してみる。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。

MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy

ヒドロホウ素化の説明：https://en.wikipedia.org/wiki/Hydroboration

R2C=CR2 + BHR2 → R2(BR2)C-CHR2

ヒドロホウ素化とは、アルケンに対してボランが付加する反応である。この反応によるボランの付加後に様々な反応を行うことが可能である。これにより、多様な分子変換が可能となる。

ヒドロホウ素化を行う代表的な意義を簡単に説明する。 ヒドロホウ素化後の酸化反応により、anti-Markovnikov型で水和体が得られる。アルケンに対して酸性条件下水和させると、多置換アルコールが多く得られる（Markovnikov則）。つまり、ヒドロホウ素化-酸化反応で、少置換数炭素アルコールを多く得られることに大きな意義がある。

使用した自作モジュールMultiOptPyのバージョン#

1.8.1c

環境#

WSL2 (Ubuntu-22.04)

手順#

1. 初期構造の準備#

モデル反応系として、以下の構造をxyzファイル形式で用意した。今回は名前をhydroboration.xyzとした。この構造はエチレンとボランの２分子で構成されている。

必要に応じてボランを9-BBN等に変えたり、エチレンに置換基を生やしてみたりしてもよい。

(事前に構造最適化したものを使用すると失敗しにくい。)

10

C      0.783536233546     -0.473210549843     -0.421681866859
C     -0.575465790810     -0.631512043438     -0.565808889130
H     -1.106598711787     -0.166899521673     -1.387899404870
H     -1.098959221364     -1.419965423361     -0.038274921716
H      1.358870827068     -1.122104737064      0.229743332825
H      1.355214803731      0.141282534944     -1.108891181524
B     -0.222297442373      0.745833387840      0.650142559883
H      0.933328222546      0.922239961540      0.993344830125
H     -0.617986677491      1.707046448923      0.053897085459
H     -0.809642243064      0.297289942133      1.595428455808

2. 遷移状態構造最適化のための初期構造の算出#

※あらかじめpythonで利用できる電子状態計算ソフトウェアのpsi4(導入方法は、psi4公式ページ[https://psicode.org/installs/v191/] を参照)が使える環境を用意しておく。

a. 以下のURLにアクセスし、Source codeをダウンロードする。#

wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.8.1c.zip
unzip v1.8.1c.zip
cd MultiOptPy-1.8.1c

b. 初期構造をカレントディレクトリにhydroboration.xyzとして保存し、以下を実行する。#

python optmain.py hydroboration.xyz -ma 50 2 7 50 8 1 -bs '6-31G(d)' -func B3LYP
  • -ma 50 2 7 50 8 1は原子のラベル番号2と7のペアと原子のラベル番号8と1のペアに、50kJ/molの活性化障壁を超えうる人工力ポテンシャルをそれぞれ付加することを示す。(詳しくはAFIR法[https://doi.org/10.1063/1.3457903]を参照)
  • -bsは基底関数の指定（特に指定しなければ6-31G(d)が使われる。）
  • -funcは計算手法の指定（Hartree-Fock(HF), post-HF(MP2等), DFTの汎関数）が可能。（特に指定しなければB3LYPが使われる。）

これを実行すると、計算レベルB3LYP/6-31G(d)で先ほどの人工力ポテンシャルをかけた上で初期構造を構造最適化する。

結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。

正常終了していれば、このディレクトリ中に、hydroboration_traj.xyzが存在するので、これをコピーして、MultiOptPy-1.8.1cディレクトリに置く。

hydroboration_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このhydroboration_traj.xyzは次のNEB計算に使用する。

※hydroboration_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。

（この人工力ポテンシャルを加えて行った構造最適化の結果はhydroboration_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。）

c. hydroboration_traj.xyzを初期パスとして、NEB法で経路の緩和を行う。#

NEB法を用いることで、先ほど得られたhydroboration_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）

LUP法と呼ばれるNEB法とよく似た方法があるが、本モジュールではNEB法を用いている。

python nebmain.py hydroboration_traj.xyz -ns 10 -spng -nd 0.15 -bs '6-31G(d)' -func B3LYP
  • -spngは経路の緩和途中の各ノードのエネルギーや勾配のRMSを可視化する。
  • -nd 0.15はノード間の距離を0.15Åとして初期パスを作成することを示す。
  • -ns 10は10回分NEB法による経路の緩和を行うことを示す。
  • -bsは基底関数の指定（特に指定しなければ6-31G(d)が使われる。）
  • -funcは計算手法の指定（Hartree-Fock(HF), post-HF(MP2等), DFTの汎関数）が可能。（特に指定しなければB3LYPが使われる。）

(参考)

d. 初期構造の決定及び遷移状態構造の計算#

MultiOptPy-1.8.1cと同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。このディレクトリ内のplot_energy_9.pngを開き、エネルギーが極大値を示すノード番号を確認する。

私が実行した環境では、3番である。 これらの２つの構造を確認する。 path_ITR_10_hydroborationディレクトリ内のhydroboration_traj_3.xyzを可視化する。

hydroboration_traj_3.xyz

10
0 1
C       0.806227839292     -0.388900330358     -0.336781636062
C      -0.594030558212     -0.552157975946     -0.497687302821
H      -1.087556722317     -0.170264966525     -1.383890351806
H      -1.081040285830     -1.408091525966     -0.041467592260
H       1.370406793246     -1.115788317229      0.239195417240
H       1.365997581863      0.150553146709     -1.094516347090
B      -0.321102945767      0.683243904726      0.586469865153
H       0.912265361621      0.782352078329      0.851822164165
H      -0.592268983609      1.713840137792      0.060355872055
H      -0.778898080287      0.305213848467      1.616499911425

※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。

hydroboration_traj_3.xyzをMultiOptPy-1.8.1cと同じディレクトリ内にコピーする。

そして、以下を実行する。

python optmain.py hydroboration_traj_3.xyz -opt rfo3_bofill -order 1 -fc 5 -bs '6-31G(d)' -func B3LYP
  • -opt rfo3_bofillは遷移状態構造の最適化向けのoptimizerを指定する。
  • -order 1は一次の鞍点を求めることを指定する。（デフォルトだと極小値を求めるようになっている。）
  • -fc 5は5回の反復回数当たり1回正確なへシアンを計算することを指定する。

実行して得られた正確な遷移状態構造を以下に示す。

(実行して得られた正確な遷移状態構造はhydroboration_traj_3_optimized.xyzとして保存されている。)

10
OptimizedStructure
C      0.784821119651     -0.448395240206     -0.377976183396
C     -0.587738861447     -0.590398727981     -0.556204691183
H     -1.067996946747     -0.164277564001     -1.429152881021
H     -1.117062031001     -1.399435764717     -0.066674000230
H      1.335157682162     -1.116164149994      0.276826282619
H      1.384178316069      0.117657945131     -1.083621475571
B     -0.270232354672      0.708431432318      0.631947919477
H      0.925646115199      0.901365797011      0.852173782016
H     -0.664091082401      1.703542038936      0.095129815649
H     -0.722681956815      0.287674233503      1.657551431640

検証#

今回求めた遷移状態構造が適当なものか確かめるために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      0.784821119651     -0.448395240206     -0.377976183396
C     -0.587738861447     -0.590398727981     -0.556204691183
H     -1.067996946747     -0.164277564001     -1.429152881021
H     -1.117062031001     -1.399435764717     -0.066674000230
H      1.335157682162     -1.116164149994      0.276826282619
H      1.384178316069      0.117657945131     -1.083621475571
B     -0.270232354672      0.708431432318      0.631947919477
H      0.925646115199      0.901365797011      0.852173782016
H     -0.664091082401      1.703542038936      0.095129815649
H     -0.722681956815      0.287674233503      1.657551431640
""")

# Set calculation parameters
psi4.set_options({
    'basis': '6-31G(d)',
    '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('b3lyp', 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('b3lyp', 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 *
	----------------------------------------------------------------------------------------------
	     1    -105.21930029   -1.05e+02 o    2.77e-05 *    5.65e-06 *    3.02e-04 *    1.23e-04 *  ~
	----------------------------------------------------------------------------------------------

一回の計算で、収束条件（電子状態計算ソフトウェアGaussianのデフォルトの基準と同等）を満たしていることがわかる。

次に、基準振動数解析の結果を確認する。以下に結果の一例を示す。

=== Vibrational Frequency Analysis ===

  non-mass-weighted Hessian:       Symmetric? True   Hermitian? True   Lin Dep Dim?  3 (0)
  projection of translations (True) and rotations (True) removed 6 degrees of freedom (6)
  total projector:                 Symmetric? True   Hermitian? True   Lin Dep Dim?  6 (6)
  mass-weighted Hessian:           Symmetric? True   Hermitian? True   Lin Dep Dim?  3 (0)
  pre-proj  low-frequency mode:    0.0000i [cm^-1]
  pre-proj  low-frequency mode:    0.0000i [cm^-1]
  pre-proj  low-frequency mode:    0.0000i [cm^-1]
  pre-proj  low-frequency mode:    0.0000  [cm^-1]
  pre-proj  low-frequency mode:    1.5737  [cm^-1]
  pre-proj  low-frequency mode:    4.0978  [cm^-1]
  pre-proj  low-frequency mode:    6.3332  [cm^-1]
  pre-proj  all modes:['270.7069i' '0.0000i' '0.0000i' '0.0000' '1.5737' '4.0978' '6.3332'
 '234.9413' '564.6771' '691.8276' '757.8373' '837.0446' '895.6937'
 '1007.2197' '1055.5235' '1064.1080' '1143.1079' '1144.9831' '1191.7997'
 '1254.8953' '1316.9096' '1498.8766' '1586.1523' '2397.0525' '2570.5145'
 '2671.2666' '3179.9184' '3189.2019' '3259.2592' '3280.7052']
  projected mass-weighted Hessian: Symmetric? True   Hermitian? True   Lin Dep Dim?  6 (6)
  post-proj low-frequency mode:  270.7068i [cm^-1] (V)
  post-proj low-frequency mode:    0.0000i [cm^-1] (TR)
  post-proj low-frequency mode:    0.0000i [cm^-1] (TR)
  post-proj low-frequency mode:    0.0000i [cm^-1] (TR)
  post-proj low-frequency mode:    0.0000  [cm^-1] (TR)
  post-proj low-frequency mode:    0.0000  [cm^-1] (TR)
  post-proj low-frequency mode:    0.0000  [cm^-1] (TR)
  post-proj  all modes:['270.7068i' '0.0000i' '0.0000i' '0.0000i' '0.0000' '0.0000' '0.0000'
 '234.9407' '564.6771' '691.8276' '757.8372' '837.0446' '895.6937'
 '1007.2197' '1055.5235' '1064.1080' '1143.1079' '1144.9831' '1191.7997'
 '1254.8953' '1316.9096' '1498.8766' '1586.1523' '2397.0525' '2570.5145'
 '2671.2666' '3179.9184' '3189.2019' '3259.2592' '3280.7052']

post-proj low-frequency mode: 270.7068i [cm^-1] (V)より、１個の虚振動が確認された。よって、今回自作モジュール(MultiOptPy)で求めた遷移状態構造は、正確な遷移状態構造であることがわかる。

そして、虚振動方向をvibration_mode_animation.xyzかvibration_mode_animation_amplified.xyz(振動が強調されている。)から確認すると、生成系と反応系をつなぐように振動していることが確認できる。（IRC計算をしたほうが正確だが、今回は単純な反応なので虚振動方向の確認で十分と判断した。）

※アニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。

終わりに#

　 　自作モジュールで、遷移状態構造を算出する手順を説明した。さらに、求めた遷移状態構造を既存のモジュールを用いて正しい遷移状態構造なのかを検証した。

参考#

【計算化学】自作pythonモジュールで遷移状態構造を求めてみる
https://ss0832.github.io/posts/20250418_mop_usage_1/
Author
ss0832
Published at
2025-04-19
【Gaussian】Gaussian 16実行ファイルの機能解説
【C】C言語の基本