Home
3435 words
17 minutes
【計算化学】自作pythonライブラリで遷移状態構造を求めてみる(3-butyl-1-ylラジカルの閉環反応(3-exo-trig閉環体の生成), NNP使用）
2025-09-08
Computational Chemistry
MultiOptPy
/
python

最終更新：2025-09-08

概要#

本記事では、自作ライブラリ(MultiOptPy)で、3-butyl-1-ylラジカルの閉環反応(3-exo-trig閉環体の生成)の遷移状態構造を算出してみる。使用する計算レベルは、Meta社のFAIR Chemistryが開発したニューラルネットワークポテンシャル(NNP)であるUMA(Meta’s Universal Model for Atoms)を使用する。

MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonライブラリである。

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

ラジカル反応について

Baldwin則について

分子内閉環反応の起こりやすさをまとめた経験則である。 本記事で取り上げたラジカルは、3-exo-trig閉環体か4-endo-trig閉環体かどちらかが生成する。Baldwin則により3-exo-trig閉環体が速度論支配的に生成されることが予想できる。

今回使用したニューラルネットワークポテンシャルについて：

使用した自作ライブラリMultiOptPyのバージョン#

v1.15.4

環境#

Windows 11

※Windows 11環境下でAnaconda PowerShell Promptを使用した。

Source codeのダウンロード#

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

環境構築手順#

今回は、Windows 11のPower Shellを使用した。初めに、NNPを使用できる環境が整ったAnaconda PowerShell Promptを用意する手順を説明する。

1, https://repo.anaconda.com/archive/ より、Anaconda3-2025.06-1-Windows-x86_64.exeでAnacondaをインストールする。

2, 検索機能を使い、スタートからAnaconda PowerShell Promptを開く。

3, 以下のコマンドを実行し、仮想環境を作成する。

conda create -n (任意の仮想環境名) python=3.12.7

4, 先ほど作成した仮想環境をconda activate (仮想環境名)で起動させる。

5, 以下のコマンドを実行し、必要なライブラリを導入する。

pip install fairchem-core==2.2.0 ase==3.25.0 torch==2.6.0
  • fairchem-coreは、FAIR Chemistryが管理しているNNPを動作させるために必要なライブラリである。
  • aseはNNPに電子エネルギーを算出したい分子構造を渡すために必要なインターフェイスの役割を果たすために必要なライブラリである。
  • torchはPyTorchライブラリを指す。これはニューラルネットワークなどの機械学習を行ったり、学習結果を扱ったりするために必須なライブラリである。

これで、Anaconda PowerShell Promptから仮想環境を立ち上げることで、NNPを使用する準備が整えることが出来る。

次に、NNPを使用するために必要なModelの情報が保存されている.ptファイルのダウンロードおよびNNPの自作ライブラリへの導入方法について説明する。

1, 以下のサイトにアクセスして、uma-s-1.ptをダウンロードする。（使用許諾が下りていれば可能である。）

https://huggingface.co/facebook/UMA

2, ダウンロード後、MultiOptPy-v1.15.4ディレクトリ内に存在するsoftware_path.confに対して、uma-s-1.ptの絶対パスを用いて以下を追記する。

uma-s-1::(uma-s-1.ptの絶対パス)

これで、MultiOptPy-v1.15.4がNNPuma-s-1を使用できるようになる。

使用するNNPに関する具体的な説明#

今回使用するNNPについて具体的に説明する。

  • UMAのModel Checkpointはuma-s-1を使用した。
  • 小分子系のトレーニングセットであるOmol25(omol)を使用して学習したニューラルネットワークポテンシャルを使用する。

※自作ライブラリでの具体的な使用の仕方に関しては、ase_calculation_tools.py を参照

手順#

1. 初期構造の準備#

モデル反応系として、以下の構造を用意した。今回はファイルの名前を3_buten_1_yl_radical_3_exo_trig.xyzとした。 初期構造は以下のものを使用した。

11
OptimizedStructure
C     -1.606445179746     -0.564693010552     -0.686844598902
H     -1.232298495709     -1.394878690635     -1.275400900647
H     -2.679109119500     -0.487013586411     -0.567780570147
C     -0.782947695256      0.316349531868     -0.140896858545
H     -1.189724670640      1.132612444034      0.449347827977
C      0.714295157272      0.267956995457     -0.250716423472
H      0.990070968352     -0.556384714048     -0.925290773246
H      1.077324981397      1.179156773784     -0.737604931983
C      1.389927201787      0.108666945041      1.064716729727
H      0.908942535392     -0.436710216751      1.863960527766
H      2.409964316651      0.434937528213      1.206509971474
初期パスを求めるための初期構造

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

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

python optmain.py .\3_buten_1_yl_radical_3_exo_trig.xyz -spin 2 -elec 0 -os uma-s-1 -opt rsirfo_fsb -modelhess -ma 200 4 9
  • -opt rsirfo_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。
  • -spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。（デフォルトでは1が指定される。）今回はUMAモデルのNNPを使用し、ラジカル反応を扱うためスピンだ重度を2とした。
  • -elec Mは形式電荷をMとすることを示す。（デフォルトでは0が指定される。）
  • -ma yyy a bはyyykJ/molの活性化障壁を超えうるペア同士を近づける力を原子のラベル番号aとbのペアに構造最適化時に加えることを示す。
  • -os uma-s-1は今回使用するNNPを指定している。これを使用する際にASEライブラリが必要である。

これを実行すると、omolのデータセットを使用したuma-s-1モデルのNNPで得たエネルギーに対して、指定した人工力ポテンシャルを加えた上で初期構造を構造最適化することができる。

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

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

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

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

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

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

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

python .\nebmain.py 3_buten_1_yl_radical_3_exo_trig_traj.xyz -os uma-s-1 -ns 15 -aneb 4 5 -modelhess -spng -nd 0.5 -elec 0 -spin 2
  • -nd Nはノード間の距離をN Åとして初期パスを作成することを示す。
  • -ns nはn回分NEB法による経路の緩和を行うことを示す。
  • -fc MはM回あたりの経路緩和回数に対して1回だけ正確なHessianを計算し、経路緩和に使用する。これを使用すると、Hessianを使用しない場合の経路緩和アルゴリズムとは別のものを使用して、経路緩和を行う。
  • -spngは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。
  • -os uma-s-1は今回使用するNNPを指定している。これを使用する際にASEライブラリが必要である。
  • -aneb A Bこれを指定すると、(B+1)回の緩和ごとに、エネルギー極大値を示すノードと前後のノードの間に線形補間でA個の新規ノードを内挿するようにできる。デフォルトではこのような操作は行われない。このオプションを使用するとノードの数が徐々に増えるため、計算コストが使用しない場合と比べて増加する。一方で、エネルギー極大値を示すノード周辺にノードを増加させるため、緩和している経路中のノードが遷移状態構造付近に存在する可能性が高くすることが出来る。

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

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

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

NEB計算の結果の可視化
NEB計算の結果の可視化

※bias_force_rms.csvにて、各Iterationごとのすべてのノードの勾配のRMS値を確認できる。

経路緩和の結果、経路のエネルギー極大値を示す構造の中から目視で、緩和後に得られた経路の14番(グラフ上では15番)の構造を遷移状態構造を求める初期構造として採用した。

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

3_buten_1_yl_radical_3_exo_trig_traj_14.xyz

11
0 2
C      -1.462602360314     -0.557266538171     -0.376611754907
H      -1.176193027852     -1.460939083321     -0.902160158616
H      -2.493817315563     -0.455646643584     -0.073277618759
C      -0.547126253744      0.423591335725     -0.128204110687
H      -0.872148147028      1.367658596861      0.285855715107
C       0.883484167417      0.304152085042     -0.520329018893
H       1.075151513443     -0.507816575349     -1.215310112772
H       1.348451635385      1.231593655019     -0.832748308612
C       1.016448811738     -0.023015985553      0.909421638384
H       0.858644321878     -1.034788999108      1.234803144299
H       1.369706654640      0.712478152439      1.618560585458
NEB法により緩和した経路から得られた遷移状態構造を求めるための初期構造

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

3_buten_1_yl_radical_3_exo_trig_traj_14.xyzをMultiOptPy-v1.15.4と同じディレクトリ内にコピーする。

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

python .\optmain.py .\3_buten_1_yl_radical_3_exo_trig_traj_14.xyz -os uma-s-1 -fc 5 -order 1 -opt rsirfo_bofill -freq -tcc -elec 0 -spin 2
  • -opt rsirfo_bofillは遷移状態構造の最適化向けのoptimizerを指定することを意味する。準ニュートン法であるRS-I-RFO法を使用する。今回は-fcで正確なへシアンを計算するようにしているので、初期へシアンは正確なへシアンを使用するようになっている。
  • -order 1は一次の鞍点を求めることを指定する。（デフォルトだと極小値を求めるようになっている。）
  • -fc 5は5回の反復回数当たり1回正確なへシアンを計算することを指定する。
  • -freqは収束条件を満たした後に基準振動解析を行うことを示す。（自前で実装しているため、あくまで目安として使用することを推奨する。各振動モードをvibration_animation内のxyzファイルで可視化できる。）UMAモデルから算出されるHessianは数値微分により求めているため、原子数Zが多いとZの二乗オーダーで計算コストが急増する。
  • -tccは収束条件を厳しくすることを示す。（Gaussianのtightと同等）

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

(実行して得られた正確な遷移状態構造は計算開始時に、yyyy_mm_ddディレクトリ内に生成された新規ディレクトリ内のYYY_optimized.xyzとして保存されている。)

11
OptimizedStructure
C     -1.468875856954     -0.548813927354     -0.363531856404
H     -1.204040982408     -1.481700821280     -0.844183007456
H     -2.486797742042     -0.430758851012     -0.022807880740
C     -0.531987532719      0.427694855861     -0.149112303704
H     -0.843727419408      1.380088181312      0.259426282701
C      0.903122130453      0.313855111011     -0.536432343714
H      1.112503672617     -0.489064055062     -1.236425010822
H      1.381483925100      1.240951196062     -0.834392999433
C      0.997952515561     -0.042948025314      0.892482330164
H      0.936782597150     -1.075107290607      1.194151785597
H      1.203584692650      0.705803626381      1.640825003811
遷移状態構造

停留点に収束した分子構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。

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

Mode                                 0                   1                   2
Freq [cm^-1]                     -617.5517            276.2505            391.6846
Reduced mass [au]                   9.1610              2.1474              1.8457
Force const [Dyne/A]               -2.0585              0.0966              0.1668
Char temp [K]                       0.0000            397.4629            563.5468
Normal mode                   x         y         z            x         y         z            x         y         z   
       C                 0.00491   -0.02059    0.00690   -0.05672    0.02967   -0.10893    0.11718   -0.00707   -0.02325
       H                -0.00298   -0.02365    0.01480   -0.12450    0.13492   -0.34900    0.36293    0.09688   -0.09313
       H                 0.00364    0.01355   -0.01441   -0.05194    0.01790   -0.08867    0.07538   -0.27821   -0.05878
       C                 0.11884   -0.02568    0.13829    0.00588   -0.06900    0.10882   -0.03511    0.11088    0.05571
       H                 0.07368    0.02688   -0.02496    0.03643   -0.10342    0.21155   -0.04681    0.09727    0.07064
       C                 0.05118    0.01302   -0.06486   -0.02422   -0.05201   -0.00144   -0.07615   -0.06208   -0.01223
       H                 0.02651    0.00439   -0.05505   -0.08409   -0.12989    0.07101   -0.29567   -0.13575    0.00583
       H                 0.05704    0.00973   -0.04706   -0.06130   -0.08028   -0.15033    0.07314   -0.16876   -0.09974
       C                -0.18810    0.02984   -0.05739    0.07242    0.08148    0.02233   -0.01997   -0.01095   -0.00507
       H                 0.00336    0.01117   -0.07376    0.25776    0.10050    0.12684   -0.20913    0.02331    0.07006
       H                -0.00441   -0.00164   -0.07289    0.05917    0.17773   -0.06886    0.20742   -0.00125   -0.07537
       
(...snip...)

その結果、虚振動が１つであることが確認できた。つまりこの構造は遷移状態構造である。

次に、vibration_animation内の虚振動を示す分子振動が示されたxyzファイル(mode_1_XXXi_wave_number.xyz)をAvogadroで確認すると、想定される反応系と生成系をつなぐ方向に振動していることを確認できた。

終わりに#

　 　自作ライブラリで、UMAモデルのニューラルネットワークポテンシャル(NNP)を用いて、3-butyl-1-ylラジカルの閉環反応(3-exo-trig閉環体の生成)のある1つの遷移状態構造を算出する手順を説明した。

(追記)　NEB計算時のスピン多重度の指定に誤りがあった。しかしながら、目的とする遷移状態構造は求められているため、このまま公開する。

参考#

【計算化学】自作pythonライブラリで遷移状態構造を求めてみる(3-butyl-1-ylラジカルの閉環反応(3-exo-trig閉環体の生成), NNP使用）
https://ss0832.github.io/posts/20250908_mop_usage_105/
Author
ss0832
Published at
2025-09-08
【計算化学】自作pythonライブラリで遷移状態構造を求めてみる(Co錯体触媒が関わる酸化的環化反応(配位子:Tris(2-thienyl)arsine), NNP使用）
【計算化学】自作pythonライブラリで遷移状態構造を求めてみる(3-butyl-1-ylラジカルの閉環反応(4-endo-trig閉環体の生成), NNP使用）