【計算化学】自作pythonモジュールで遷移状態構造を求めてみる（有機リチウム試薬）

最終更新：2025-06-10

概要#

本記事では、自作モジュール(MultiOptPy)で、有機リチウム試薬のアセトンのカルボニル炭素に対する求核付加反応の遷移状態構造を算出してみる。電子状態計算ソフトウェアにPySCFを使用する。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。

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

有機リチウム試薬について：

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

1.9.12

環境#

WSL2 (Ubuntu-22.04)

Source codeのダウンロード#

wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.9.12.zip
unzip v1.9.12.zip
cd MultiOptPy-1.9.12

手順#

1. 初期構造の準備#

モデル反応系として、以下の構造を用意した。今回は名前をorganolithium_reagents.xyzとした。有機リチウム試薬は反応溶液中で会合状態として存在するが、簡単のため今回のような構造を使用した。

15
OptimizedStructure
C     -0.774262822012     -1.735208971796      3.258242758686
H     -0.328530640598     -1.257335711643      4.138238015554
H     -0.276082350356     -2.704294163349      3.140952921602
H     -1.818935761366     -1.948524387409      3.512935148488
Li    -0.624337029674     -0.572122553251      1.604602252704
C      0.141520371074      0.662031806403     -1.021348342039
O     -0.322002405212      0.318425364195      0.057889784828
C      1.329823732299     -0.051081058666     -1.616308086662
H      2.127917609784     -0.095664081341     -0.884488299070
H      1.047227495389     -1.074265925608     -1.843869946994
H      1.684121175704      0.427558416576     -2.515981469136
C     -0.447636277991      1.798242942383     -1.807019696631
H      0.293096826355      2.584813620625     -1.913470681991
H     -0.710849280116      1.459622701414     -2.803445161040
H     -1.321070643279      2.187802001468     -1.306929198301

初期パスを求めるための初期構造

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

※あらかじめpythonで利用できる電子状態計算ソフトウェアのPySCF(導入方法はpipを使えば可能である。)が使える環境を用意しておく。

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

python optmain.py organolithium_reagents.xyz -opt rsirfo_fsb -func hf -bs 3-21g -pyscf -spin 0 -elec 0 -ma 100 1 6

-opt rsirfo_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。
-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。
-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。
-ma 100 1 6は100kJ/molの活性化障壁を超えうるペア同士を近づける力を原子のラベル番号1と6のペアに、構造最適化時に加えることを示す。
-func hfはHartree-Fock法を用いることを示す。
-bs 3-21Gは基底関数にPople系である3-21Gを用いることを示す。

これを実行すると、計算レベルHF/3-21Gで指定した人工力ポテンシャルをかけた上で初期構造を構造最適化することができる。

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

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

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

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

（この人工力ポテンシャルを加えて行った構造最適化の結果はorganolithium_reagents_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。今回の場合以下の構造が得られた。安定構造に近いものが得られていると判断できるため、次のNEB計算を行うことが可能である。）

初期パスを求めるための構造最適化の結果（安定構造ではない）

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

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

python nebmain.py organolithium_reagents_traj.xyz -ns 10 -modelhess -func hf -bs 3-21g -pyscf -spin 0 -elec 0 -spng -nd 0.5

-nd 0.5はノード間の距離を0.5Åとして初期パスを作成することを示す。
-ns 10は10回分NEB法による経路の緩和を行うことを示す。
-modelhessは少ない計算コストで算出された近似へシアンを用いることを示す。これを指定すると経路緩和のアルゴリズムが-fcを1以上に指定した場合と同様のものになる。
-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。
-spngは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。
-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。

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

MultiOptPy-1.9.12と同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。

パスの緩和後の各ノードのエネルギー一覧(単位) (energy_plot.csvに保存されている。)

NEB計算の結果の可視化

私が実行した環境では、28番のノード(グラフでは29番)がエネルギー極大値を示していた。 path_ITR_10_organolithium_reagentsディレクトリ内のorganolithium_reagents_traj_28.xyzを可視化する。

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

organolithium_reagents_traj_28.xyz

15
0 0
C      -0.620056138617     -1.183809570857      0.995650640930
H       0.187938574974     -1.900114688046      0.842647485768
H      -1.341508601136     -1.293046877450      0.202436552524
H      -1.151017867770     -1.552465211151      1.902602014038
Li     -0.623578449821      0.607593575547      2.045812350081
C       0.128444320860      0.711778264831     -0.010920333029
O      -0.047785164673      1.680024498705      0.826989627234
C       1.500399489047      0.096169118659     -0.173773385644
H       2.005093027910      0.106969161643      0.782215175653
H       1.425229288637     -0.910413854172     -0.543398340504
H       2.057470128021      0.698100296898     -0.890367638490
C      -0.738266498642      0.701987972413     -1.252280711351
H      -0.209879429934      1.326911548399     -1.968209855784
H      -0.884732154196     -0.270107420516     -1.695269543799
H      -1.687750524660      1.180423185098     -1.064134037626

NEB法により緩和したパスのエネルギー極大値を示す構造

構造が壊れていないので、これを遷移状態を求めるための初期構造とする。

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

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

python optmain.py organolithium_reagents_traj_28.xyz -opt rsirfo_bofill -func hf -bs 3-21g -pyscf -spin 0 -elec 0 -fc 5 -order 1 -freq -tcc

-opt rsirfo_bofillは遷移状態構造の最適化向けのoptimizerを指定することを意味する。準ニュートン法であるRS-I-RFO法を使用する。今回は-fcで正確なへシアンを計算するようにしているので、初期へシアンは正確なへシアンを使用するようになっている。
-order 1は一次の鞍点を求めることを指定する。（デフォルトだと極小値を求めるようになっている。）
-fc 5は5回の反復回数当たり1回正確なへシアンを計算することを指定する。
-freqは収束条件を満たした後に基準振動解析を行うことを示す。（自前で実装しているため、あくまで目安として使用することを推奨する。各振動モードをvibration_animation内のxyzファイルで可視化できる。）
-tccは収束条件を厳しくすることを示す。（Gaussianのtightと同等）

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

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

15
OptimizedStructure
C     -0.719357541669     -1.208497008698      1.061783810388
H      0.137461420263     -1.772434605511      1.441519163351
H     -0.953960905020     -1.593377899658      0.073397282697
H     -1.580596810912     -1.485354438733      1.694717078382
Li    -0.686862828180      0.630420064814      1.999743291740
C      0.209717088592      0.759851640455     -0.043980442930
O      0.104296895571      1.664267221437      0.848153447531
C      1.535662126820      0.079223785740     -0.281694545268
H      2.092327695643      0.009000019337      0.640444167929
H      1.405999676394     -0.901418831568     -0.711078190289
H      2.101306547036      0.695032141084     -0.978749184215
C     -0.744283339989      0.761738214047     -1.216649479648
H     -0.408098052627      1.535174781113     -1.904191507826
H     -0.750608565880     -0.184563056567     -1.732396005331
H     -1.743003406040      1.010937972709     -0.891018886512

遷移状態構造

30回程度の反復計算で遷移状態構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。

以下に-freqオプションで生成されたnormal_modes.txtの一部を示す。

Mode                                 0                   1                   2
Freq [cm^-1]                     -335.1051            116.6662            160.4254
Reduced mass [au]                   3.8287              1.0360              2.8651
Force const [Dyne/A]               -0.2533              0.0083              0.0434
Char temp [K]                       0.0000            167.8566            230.8164
Normal mode                   x         y         z            x         y         z            x         y         z
       C                 0.07683    0.11077   -0.12629    0.00788   -0.01207    0.01776    0.06838   -0.01072    0.02141
       H                 0.02157    0.12444    0.00525   -0.16840    0.05613    0.53886    0.19340    0.16583   -0.00306
       H                 0.03930   -0.07222   -0.03843    0.56866   -0.15287   -0.05927    0.06911   -0.10005    0.05529
       H                 0.03368    0.36149   -0.06729   -0.32529    0.04750   -0.41769    0.16323   -0.13604    0.09625
      Li                 0.01662    0.06870    0.01363   -0.00064   -0.01082   -0.01010   -0.25534    0.06341   -0.10907
       C                -0.06650   -0.13249    0.07545    0.00350    0.01111   -0.01038    0.01435   -0.00586    0.00982
       O                -0.00669   -0.00349    0.01862    0.01512   -0.00300    0.00126    0.04711   -0.01531    0.02362
       C                -0.00935   -0.02047    0.00782   -0.00604   -0.00280   -0.01894   -0.00940   -0.03124   -0.03845
       H                 0.00701   -0.03855   -0.00191   -0.00633   -0.01518   -0.01971   -0.01143   -0.11633   -0.04451
       H                 0.06387   -0.00370   -0.05524   -0.01797    0.00109   -0.02343   -0.04159    0.00461   -0.11139
       H                -0.05890    0.06881    0.04757    0.00208   -0.00596   -0.01521    0.00831    0.00216    0.00577
       C                -0.01113   -0.02931    0.01745   -0.02204    0.01579    0.01119   -0.01413    0.03520    0.03105
       H                 0.04224    0.00591    0.08457   -0.02944    0.00546   -0.00410    0.07328   -0.06327   -0.03780
       H                -0.01577    0.00213   -0.04325   -0.04950    0.01119    0.02063   -0.15665   -0.00265    0.10391
       H                -0.02165   -0.01990   -0.01720   -0.01049    0.03231    0.03461    0.02740    0.19768    0.03634

(...snip...)

その結果、虚振動が１つであることが確認できた。

次に、vibration_animation内のmode_1_335i_wave_number.xyzをAvogadroで確認すると、想定される反応系と生成系をつなぐ方向に振動していることを確認できた。

終わりに#

　　自作モジュールで、電子状態計算モジュールであるPySCFを用いて、計算レベルHF/3-21Gで有機リチウム試薬のアセトンのカルボニル炭素に対する求核付加反応の遷移状態構造を算出する手順を説明した。

参考#

https://github.com/pyscf/pyscf (PySCFのgithubのレポジトリ)
https://github.com/ss0832/MultiOptPy (自作モジュールMultiOptPyのレポジトリ)
https://avogadro.cc/ (Avogadro、分子構造可視化ツール)
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.

概要#