5361 words
27 minutes
【計算化学】自作pythonライブラリで遷移状態構造を求めてみる(Cu-DPPBZ系錯体触媒が関わる有機金属反応, NNP使用）
最終更新：2025-09-02
概要
本記事では、自作ライブラリ(MultiOptPy)で、Cu-DPPBZ系錯体触媒が関わる有機金属反応の遷移状態構造を算出してみる。使用する計算レベルは、Meta社のFAIR Chemistryが開発したニューラルネットワークポテンシャル(NNP)であるUMA(Meta’s Universal Model for Atoms)を使用する。
MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonライブラリである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
移動挿入反応について：
- https://en.wikipedia.org/wiki/Migratory_insertion
- https://chem.libretexts.org/Bookshelves/Inorganic_Chemistry/Organometallic_Chemistry_(Evans)/04%3A_Fundamentals_of_Organometallic_Chemistry/4.05%3A_Migratory_Insertion-_12-Insertions
今回参考にした論文
- Org. Lett. 2025, 27, 2005-2010 (論文中の計算に使用している系の配位子の置換基を変えて計算を行った。)
今回使用したニューラルネットワークポテンシャルについて：
- https://ai.meta.com/blog/meta-fair-science-new-open-source-releases/ (UMAの公開に関する記事)
- https://github.com/facebookresearch/fairchem (FAIR Chemistryの提供するGitHubのレポジトリ)
- https://fair-chem.github.io/ (同上のレポジトリの内容に関して説明したサイト)
- https://huggingface.co/facebook/UMA (NNPの配布サイト, Hugging Faceへのアカウント登録と配布元の使用許諾が必要である。)
- arXiv preprint arXiv:2505.08762 (2025). (プレプリント)
使用した自作ライブラリMultiOptPyのバージョン
v1.15.2環境
Windows 11※Windows 11環境下でAnaconda PowerShell Promptを使用した。
Source codeのダウンロード
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.15.2.zip
unzip v1.15.2.zip
cd MultiOptPy-v1.15.2環境構築手順
今回は、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.74, 先ほど作成した仮想環境を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/UMA2, ダウンロード後、MultiOptPy-v1.15.2ディレクトリ内に存在するsoftware_path.confに対して、uma-s-1.ptの絶対パスを用いて以下を追記する。
uma-s-1::(uma-s-1.ptの絶対パス)これで、MultiOptPy-v1.15.2がNNPuma-s-1を使用できるようになる。
使用するNNPに関する具体的な説明
今回使用するNNPについて具体的に説明する。
- UMAのModel Checkpointは
uma-s-1を使用した。 - 小分子系のトレーニングセットである
Omol25(omol)を使用して学習したニューラルネットワークポテンシャルを使用する。
※自作ライブラリでの具体的な使用の仕方に関しては、ase_calculation_tools.py を参照
手順
1. 初期構造の準備
モデル反応系として、以下の構造を用意した。今回はファイルの名前をboration_Cu_dppbz.gjfとした。 初期構造は以下のものを使用した。
%chk=XXX
# hf/3-21g geom=connectivity
Title Card Required
0 1
C -2.07534734 -0.07151738 -0.01487024
H -2.14242599 -1.10549585 -0.28187033
C -3.48591645 0.54614783 0.00473047
H -3.41883769 1.58012626 0.27173049
H -4.09083583 0.03120273 0.72153166
Se -0.97857754 0.86212136 -1.31449124
C -1.78037812 0.69916811 -3.07351432
F -1.86501019 -0.60538411 -3.41038363
F -3.01691602 1.24062791 -3.05633197
F -1.01716206 1.34886520 -3.97788977
Cu -1.27354677 0.09143587 1.74415283
P -0.51368013 1.80254279 2.85007247
H -1.49645579 2.38390043 3.57026627
H 0.04583795 2.74781836 2.06528519
P -0.76847630 -1.38222422 3.26124738
H -0.37964611 -2.57045151 2.75190643
H -1.81533574 -1.60188796 4.08485589
C 0.71331455 0.95238480 3.94175535
C 1.75353746 1.64990657 4.56251894
C 0.59818083 -0.48678095 4.12755277
C 2.67373753 0.95559518 5.36586963
H 1.84852325 2.70698610 4.42664944
C 1.52881663 -1.15885965 4.92515841
C 2.56232986 -0.43702619 5.54565465
H 3.46720801 1.48986094 5.84531169
H 1.45455196 -2.21741749 5.06241763
H 3.27198384 -0.95028807 6.16035255
C -4.12239720 0.41679319 -1.39160757
H -3.51747795 0.93173830 -2.10840879
H -4.18947610 -0.61718523 -1.65860763
H -5.10246812 0.84595021 -1.37798890
B 2.15605356 0.01238532 -0.64868019
H 1.65404720 1.07297880 -0.52414356
O 1.45161196 -1.32355551 -0.32719503
O 3.59831554 -0.21060735 -1.15322349
C 2.36239495 -2.18421916 -1.01363697
C 3.73704673 -1.56714149 -0.72674405
C 2.25717677 -3.64056003 -0.52409505
H 1.27343176 -4.01270602 -0.72069640
H 2.97594009 -4.24288270 -1.03933573
H 2.44902878 -3.67796158 0.52790774
C 2.08579727 -2.15759688 -2.52832348
H 2.80792455 -2.76452417 -3.03343860
H 1.10413720 -2.53855707 -2.71838759
H 2.15322013 -1.15123345 -2.88557130
C 4.02914633 -1.62583988 0.78417667
H 4.07799683 -2.64747217 1.09860491
H 4.96387826 -1.14556232 0.98568301
H 3.24853698 -1.12526178 1.31805857
C 4.87388381 -2.29080777 -1.47193865
H 5.80382262 -1.80486860 -1.26223965
H 4.92216136 -3.30898973 -1.14657283
H 4.68553171 -2.26135697 -2.52482754
2. 遷移状態構造最適化のための初期構造の算出
a. 初期構造をカレントディレクトリにboration_Cu_dppbz.xyzとして保存し、以下を実行する。
python .\optmain.py .\boration_Cu_dppbz.gjf -opt rsirfo_fsb -os uma-s-1 -ma 300 1 32 300 11 33 - Gaussianのinputファイルである
.gjfファイルでも認識するようになっている。 -opt rsirfo_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。（デフォルトでは1が指定される。）-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（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。
正常終了していれば、このディレクトリ中に、boration_Cu_dppbz_traj.xyzが存在するので、これをコピーして、MultiOptPy-v1.15.2ディレクトリに置く。
boration_Cu_dppbz_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このboration_Cu_dppbz_traj.xyzは次のNEB計算に使用する。
※boration_Cu_dppbz_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
この人工力ポテンシャルを加えて行った構造最適化の結果はboration_Cu_dppbz_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。今回の場合以下の構造が得られた。生成系の安定構造に近いものが得られていると判断できるため、次のNEB計算を行うことが可能である。今回収束に時間がかかったため、ITR. 161の時点で、以上のコマンドを実行した際に生成されるディレクトリ内にend.txtを新規作成することで停止させている。
b. boration_Cu_dppbz_traj.xyzを初期パスとして、NEB法で経路の緩和を行う。
NEB法を用いることで、先ほど得られたboration_Cu_dppbz_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
python .\nebmain.py boration_Cu_dppbz_traj.xyz -os uma-s-1 -ns 20 -modelhess -spng -nd 0.3 -aneb 2 6-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.2と同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。 そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
パスの緩和後の各ノードのエネルギー一覧(単位
※bias_force_rms.csvにて、各Iterationごとのすべてのノードの勾配のRMS値を確認できる。
経路のエネルギー極大値を示す構造の中から、緩和後の経路の58番(グラフ上では59番)の構造を遷移状態構造を求める初期構造として採用した。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
boration_Cu_dppbz_traj_58.xyz
53
0 1
C -1.694530177919 0.642624123496 -1.463079730545
H -1.775038220338 -0.316292064838 -1.957242934790
C -3.102842773257 0.915441632813 -0.820407565520
H -3.136071813953 1.906267713203 -0.349101598509
H -3.284706106168 0.189537082083 -0.014415475882
Se -1.181890366179 2.004258629935 -2.733656114757
C -2.210667573697 1.443396658051 -4.296056320858
F -2.644942665964 0.177840577646 -4.204956724428
F -3.298068841211 2.215066722675 -4.503628568022
F -1.460000761112 1.534480523013 -5.404647681960
Cu -1.138300882827 0.977447741748 0.692379349123
P -1.002547373454 2.541747560870 2.247154070200
H -2.147094814393 2.901127982564 2.985931873142
H -0.395494254094 3.803514100760 2.219858443115
P -1.249711491787 -0.539099246951 2.432723370381
H -0.805906403909 -1.866267412536 2.551148473469
H -2.422196235290 -0.645244330563 3.202531465379
C 0.060057184822 1.682130718660 3.464561783164
C 1.065295034656 2.329641733402 4.177107905226
C -0.068199293686 0.286989092021 3.564088644684
C 1.938899753622 1.605139450712 4.973094274367
H 1.194467595386 3.401032461208 4.090048328827
C 0.794781047035 -0.423041952864 4.394177688015
C 1.798015898650 0.231240144796 5.092176497925
H 2.738465540016 2.113829916260 5.494701266306
H 0.703272170437 -1.498405818729 4.479057691667
H 2.479054432532 -0.332449078930 5.716670008505
C -4.410612413227 0.847285965427 -1.662408004803
H -4.469083265554 1.558721065943 -2.486921313178
H -4.639499934604 -0.140615636686 -2.070106255009
H -5.246063688824 1.100383195559 -0.998812377771
B 0.022730761185 0.073877563581 -0.763256940973
H 0.370702822190 0.619857804717 0.446430403303
O -0.107734724910 -1.343654169726 -0.848466101929
O 1.330274239020 0.417371022150 -1.323807453079
C 1.059986264964 -1.843999708329 -1.509168982490
C 2.124427445449 -0.772112210488 -1.138495623277
C 1.277722377755 -3.261346953081 -1.028153412365
H 0.410433174327 -3.863696785809 -1.304633154371
H 2.158338668262 -3.685007179522 -1.515018033175
H 1.385614098648 -3.314472020916 0.054958070675
C 0.896346717478 -1.807469584661 -3.032109784131
H 1.670938398560 -2.415943722330 -3.502516259732
H -0.074335639442 -2.178560561051 -3.346898523456
H 0.994701872823 -0.784606987043 -3.394641148300
C 2.669007175290 -0.934580428206 0.297068755125
H 3.151177505797 -1.910496823563 0.404784364358
H 3.425747916267 -0.173146217137 0.493028505335
H 1.902188423944 -0.847857900312 1.067717872847
C 3.346752521376 -0.683302678823 -2.052335204804
H 4.022247609495 0.077129686706 -1.657612286202
H 3.899283302473 -1.621041843607 -2.096433157147
H 3.074609763341 -0.394669553299 -3.062412373674
構造が壊れていないので、これを遷移状態を求めるための初期構造とする。
boration_Cu_dppbz_traj_58.xyzをMultiOptPy-v1.15.2と同じディレクトリ内にコピーする。
そして、以下を実行する。
python .\optmain.py .\boration_Cu_dppbz_traj_58.xyz -freq -tcc -opt rsirfo_bofill -fc 5 -order 1 -os uma-s-1-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として保存されている。)
53
OptimizedStructure
C -0.844924293711 0.698979230164 -0.902717183904
H -0.651178148697 0.389497151145 -1.925160517128
C -2.112601153277 0.008836381218 -0.394488803625
H -2.459250058890 0.480688730248 0.530728781474
H -1.820966380872 -1.015684766156 -0.135654804292
Se -0.977202495502 2.637030550544 -1.114887571776
C -1.746139363945 2.801052623023 -2.913600613497
F -1.440659120524 1.760569016576 -3.702295924585
F -3.087428683567 2.920155891544 -2.931436679351
F -1.262867179258 3.911619434765 -3.495199456442
Cu 0.296425572231 0.921805827243 0.933141367493
P 0.323955860744 2.737905935485 2.260063117203
H -0.573623805326 3.812710950897 2.132508961987
H 1.472646955077 3.501110663920 2.533155886150
P -0.333023799920 -0.360059339136 2.753405512167
H 0.480308943582 -1.377003636998 3.285544836381
H -1.563552790374 -1.019785468181 2.931394526627
C -0.091649645710 2.156529619693 3.948470364955
C -0.128737908858 3.048905998033 5.020994399248
C -0.376654617127 0.803946840319 4.168790944540
C -0.441185068067 2.610112340981 6.294219828310
H 0.089451323977 4.096641366301 4.853397473283
C -0.690388214469 0.373940815614 5.458917305443
C -0.722972226794 1.266695602901 6.514143893927
H -0.465925848385 3.312720505086 7.116653405222
H -0.910893606871 -0.672010600553 5.634075150746
H -0.967857679677 0.918604580130 7.508936858980
C -3.261136430624 -0.028687473808 -1.394827985685
H -3.630237795480 0.974243784810 -1.608852583960
H -2.940370965846 -0.475899479215 -2.338241160158
H -4.097889296097 -0.616147383401 -1.011748877154
B 0.934314468004 -0.204407749771 -0.763316769265
H 1.656327681199 0.271776472970 0.207324500290
O 0.668041936182 -1.594739498220 -0.628564535867
O 1.548097223664 0.041932988580 -2.020103401005
C 0.979537837362 -2.221504155955 -1.881524365089
C 2.002335150228 -1.225271187960 -2.509989293278
C 1.537193807455 -3.606101713641 -1.598360384252
H 0.750990394013 -4.236546812291 -1.180303912576
H 1.896663708540 -4.075131887197 -2.517096238333
H 2.354951089041 -3.561358851800 -0.881817214815
C -0.300779526384 -2.350346297813 -2.703297240300
H -0.136121900182 -2.946492355234 -3.602469993364
H -1.058361105486 -2.846234063788 -2.094499886660
H -0.687655898502 -1.378048906950 -3.007829244071
C 3.420467834219 -1.448364621084 -1.989497802740
H 3.859717100360 -2.365645399380 -2.384725834831
H 4.038182273760 -0.605142154487 -2.299405882717
H 3.429407336724 -1.495183705730 -0.899478249832
C 2.003605658036 -1.193258562481 -4.027451603062
H 2.752970079471 -0.481237915813 -4.374723848223
H 2.250527178196 -2.176094878078 -4.435658411215
H 1.036115596355 -0.881624437072 -4.416640841374
停留点に収束した分子構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。
以下に-freqオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -180.8109 9.1388 9.6447
Reduced mass [au] 6.6025 5.6060 5.4272
Force const [Dyne/A] -0.1272 0.0003 0.0003
Char temp [K] 0.0000 13.1486 13.8766
Normal mode x y z x y z x y z
C 0.05258 -0.05823 -0.08788 -0.01335 -0.00551 -0.00851 -0.00085 0.01847 -0.00982
H -0.09458 0.01312 -0.13984 0.00085 -0.00177 -0.00692 0.00031 0.00955 -0.00696
C 0.05604 -0.02844 -0.01681 -0.01870 -0.00989 -0.02773 -0.00174 0.02403 -0.00478
H 0.10046 -0.01044 -0.00940 -0.03308 -0.01245 -0.03193 0.00043 0.02934 -0.00658
H 0.04249 -0.03204 -0.02154 -0.01965 -0.00960 -0.02551 -0.00348 0.02483 0.00054
Se 0.01543 -0.02171 -0.00233 -0.01541 -0.00518 -0.00531 0.00368 0.01710 -0.02319
C 0.00760 -0.00357 0.00132 0.01293 -0.00308 -0.01720 0.00133 0.00708 -0.02318
F 0.00007 0.00439 -0.00845 0.03086 0.00089 -0.01554 -0.00426 0.00002 -0.01594
F 0.00807 -0.00032 0.01297 0.01270 -0.00950 -0.03887 0.00170 0.01222 -0.02147
F 0.00429 0.00071 0.00951 0.01724 0.00109 -0.00566 0.00455 0.00085 -0.03236
Cu 0.01489 0.00566 0.02934 -0.03791 -0.00818 0.00621 -0.00516 0.03109 -0.00856
P -0.00896 0.01024 0.00612 -0.02766 -0.00060 -0.00274 -0.05043 0.01391 0.01234
H 0.00598 0.02481 0.02389 -0.04990 -0.01740 0.01001 -0.09349 -0.01972 0.03148
H 0.00973 -0.01431 -0.00767 -0.03607 0.02187 -0.03086 -0.07788 0.05557 0.01081
P -0.00451 0.00149 0.00306 -0.05822 0.00856 0.00895 0.02689 -0.00920 -0.02340
H 0.01877 0.01123 -0.01691 -0.09361 -0.02644 -0.00426 0.04715 -0.00467 -0.04603
H 0.00474 -0.01203 0.01361 -0.08318 0.05778 0.02036 0.03837 -0.03063 -0.02294
C -0.00164 0.00220 0.00424 0.01999 0.00086 0.00943 -0.01241 -0.02118 0.00986
C -0.00088 0.00092 0.00486 0.07228 -0.00284 0.01429 -0.01557 -0.03779 0.02357
C -0.00057 0.00115 0.00293 0.00654 0.00450 0.01428 0.02161 -0.03098 -0.00554
C -0.00034 -0.00032 0.00449 0.11157 -0.00314 0.02385 0.01490 -0.06359 0.02218
H -0.00159 0.00119 0.00566 0.08334 -0.00571 0.01077 -0.04169 -0.03042 0.03561
C 0.00048 0.00019 0.00259 0.04594 0.00420 0.02376 0.05189 -0.05685 -0.00682
C 0.00013 -0.00045 0.00321 0.09841 0.00041 0.02860 0.04883 -0.07318 0.00692
H -0.00041 -0.00111 0.00512 0.15303 -0.00621 0.02771 0.01253 -0.07628 0.03294
H 0.00103 0.00001 0.00207 0.03628 0.00688 0.02759 0.07851 -0.06445 -0.01869
H 0.00031 -0.00129 0.00293 0.12954 0.00012 0.03616 0.07296 -0.09341 0.00579
C 0.01862 -0.00424 0.02368 -0.00412 -0.01104 -0.04443 -0.00314 0.02108 -0.00304
H 0.03470 0.00343 0.03145 -0.00357 -0.01159 -0.04801 -0.00121 0.02076 -0.00785
H -0.02273 -0.01730 0.01583 0.01074 -0.00858 -0.04055 -0.00531 0.01557 -0.00118
H 0.01827 0.01719 0.05477 -0.00817 -0.01387 -0.05764 -0.00390 0.02485 0.00110
B -0.14637 0.06756 -0.06500 -0.01539 -0.00405 0.01136 0.00275 0.02193 0.00193
H -0.06297 0.03304 -0.07166 -0.03059 -0.01570 0.02879 0.00460 0.05297 -0.01585
O -0.05154 0.02548 -0.02225 -0.01730 -0.00533 -0.00737 0.00852 0.02395 0.03854
O -0.05033 0.01568 -0.04931 0.00574 0.00913 0.02420 -0.00650 -0.00780 -0.00836
C -0.01333 0.01220 0.00195 0.00628 0.00721 -0.00779 -0.00178 -0.00549 0.05075
C -0.01401 0.01539 0.00331 0.01659 0.01449 0.02049 -0.00885 -0.01893 0.01815
C -0.00374 0.01864 0.02108 0.00278 0.00512 -0.01103 0.00260 0.00203 0.07867
H -0.00319 0.01334 0.01493 -0.00414 0.00004 -0.03172 0.00727 0.01095 0.10093
H 0.01335 0.01478 0.02945 0.02038 0.01460 -0.00899 -0.00527 -0.01950 0.08659
H -0.01376 0.03174 0.03189 -0.01055 -0.00085 0.00455 0.00924 0.02020 0.06997
C 0.00390 -0.01628 -0.02476 0.02154 0.01353 -0.03259 -0.00915 -0.02616 0.06564
H 0.02326 -0.02609 -0.01497 0.03859 0.02116 -0.03454 -0.01638 -0.04445 0.07647
H -0.00144 -0.01538 -0.03012 0.01056 0.00798 -0.05076 -0.00220 -0.01532 0.08311
H -0.00100 -0.02441 -0.04511 0.02661 0.01581 -0.03177 -0.01430 -0.03324 0.04955
C -0.01979 0.03384 0.02823 0.00727 0.01136 0.04461 -0.00403 -0.00416 0.01128
H -0.00907 0.03010 0.04778 0.01583 0.01578 0.04387 -0.00597 -0.01297 0.02956
H -0.01931 0.03214 0.02372 0.01187 0.01520 0.06420 -0.00806 -0.01063 -0.01431
H -0.03679 0.05009 0.02895 -0.01283 0.00078 0.04432 0.00539 0.02216 0.01233
C 0.01372 -0.01719 0.00046 0.04463 0.02894 0.02083 -0.02209 -0.05535 0.01736
H 0.01150 -0.01452 0.00030 0.05021 0.03315 0.04148 -0.02634 -0.06241 -0.00626
H 0.03106 -0.02144 0.02054 0.05337 0.03313 0.01604 -0.02398 -0.06474 0.03883
H 0.01713 -0.03325 -0.02123 0.05146 0.03150 0.00591 -0.02602 -0.06637 0.01831
(...snip...)その結果、虚振動が１つであることが確認できた。つまりこの構造は遷移状態構造である。
次に、vibration_animation内の虚振動を示す分子振動が示されたxyzファイル(mode_1_XXXi_wave_number.xyz)をAvogadroで確認すると、想定される反応系と生成系をつなぐ方向に振動していることを確認できた。
以下のコマンドを実行して、求めた遷移状態構造に対するIRC計算を行った。
python .\optmain.py .\boration_Cu_dppbz_traj_58_optimized.xyz -freq -tcc -opt rsirfo_bofill -fc 10 -order 1 -os uma-s-1 -irc 0.5 200 lqa-irc X Y lqaは構造最適化後にIRC（固有反応座標）計算を行うためのオプションである。X[ang.]で最大ステップ幅を決定し、Yは最大ステップ計算回数を指定する。lqaはLQAと呼ばれるIRC計算アルゴリズムを使用することを指定する。-fc NはN回のIRCステップ計算ごとに正確なHessianを算出することを示す。構造最適化時のへシアンの計算頻度のオプションの指定と重複している。
実行して得られたIRC経路の末端構造を以下に示す。
Forward structure
53
IRC Step 66
C -0.8607279313 -0.3175485357 -1.5289743713
H -1.2578820872 -0.5732601032 -2.5160573608
C -1.8613425154 -0.6903647684 -0.4396856679
H -1.4737735124 -0.3518419351 0.5350215606
H -1.8756086671 -1.7829169454 -0.3881824944
Se -0.4711570189 1.6029786118 -1.6007552345
C -1.6485996038 2.1179160885 -3.0629078177
F -1.4508931675 1.3752026518 -4.1629209502
F -2.9639856444 2.0385492783 -2.7898188776
F -1.3935378783 3.3978886456 -3.3856043451
Cu 0.7914556302 -0.1089223699 1.1143856514
P 0.2778254771 1.7611368270 2.3624385700
H -0.9205902303 2.4562009383 2.1170291842
H 1.0838039771 2.8837109690 2.6219033182
P 0.2800572547 -1.3891612661 2.9446878732
H 1.0679681353 -2.3526632076 3.5995741423
H -0.9220820360 -2.1192939371 2.9427017000
C -0.0110317404 1.1607612601 4.0699275266
C -0.2184107984 2.0515705990 5.1238960519
C -0.0175244563 -0.2172946185 4.3237348933
C -0.4273242312 1.5891186092 6.4100157356
H -0.2126994580 3.1176658215 4.9319952874
C -0.2341060557 -0.6711737613 5.6255002888
C -0.4365721303 0.2221135580 6.6615919522
H -0.5827503843 2.2917829216 7.2179988190
H -0.2392477531 -1.7354804881 5.8273472838
H -0.5996129290 -0.1445817510 7.6663982879
C -3.2742512131 -0.1546141185 -0.6109828000
H -3.2920718538 0.9343154155 -0.5713788241
H -3.6881497387 -0.4564104141 -1.5749750765
H -3.9329001065 -0.5354355119 0.1717773495
B 0.5581386521 -1.0787763214 -1.3758972227
H 1.3105739135 -0.2294373685 -0.3861585102
O 0.5139292853 -2.4262548998 -0.8861155834
O 1.4203328354 -1.0322064779 -2.5003491073
C 1.2540050757 -3.2413613149 -1.8030683035
C 2.2197538897 -2.2150397469 -2.4684408003
C 1.9419330739 -4.3543538563 -1.0326573984
H 1.1913392333 -5.0268519857 -0.6142613007
H 2.5925868298 -4.9372430174 -1.6887387075
H 2.5376445417 -3.9544858894 -0.2135664910
C 0.2651639533 -3.8331658400 -2.8041329536
H 0.7476345927 -4.5355855853 -3.4860552140
H -0.5124170411 -4.3638131053 -2.2529760095
H -0.2063640704 -3.0431889807 -3.3901479763
C 3.4578307935 -1.9470656090 -1.6157094280
H 4.1486609629 -2.7918931421 -1.6329694903
H 3.9688752446 -1.0701407117 -2.0143715547
H 3.1806308642 -1.7383982514 -0.5819318598
C 2.6346620512 -2.5593171200 -3.8875597322
H 3.3174555424 -1.7950618209 -4.2606339281
H 3.1488664004 -3.5228499952 -3.9200001187
H 1.7726869284 -2.5953142157 -4.5507932412
backward structure
53
IRC Step 49
C -0.8209392177 -0.0113574511 -0.7945074844
H -0.2843251868 -0.3899681361 -1.6649088457
C -2.1333026427 -0.7606813453 -0.5719899084
H -2.6696940473 -0.3069551945 0.2678076913
H -1.8719029379 -1.7747101752 -0.2427730494
Se -1.0012053753 1.8878045981 -1.1232389214
C -1.5628167246 2.0182186469 -3.0024865709
F -1.1810049889 0.9475758149 -3.7234836068
F -2.8927110155 2.1489369287 -3.1883765702
F -1.0022733539 3.1026324875 -3.5716820044
Cu 0.3446236394 -0.0487479353 0.8465061324
P 0.8348356373 1.6908354091 2.2438403583
H 0.1669673868 2.9230918157 2.1272997179
H 2.1002675379 2.2274571724 2.5538024880
P -0.1078261645 -1.3139154199 2.7324590929
H 0.7308302796 -2.3165876859 3.2610079640
H -1.3206986005 -1.9850100587 2.9773102507
C 0.2845722804 1.1865599621 3.9198628241
C 0.2585501354 2.0946492286 4.9795756550
C -0.1207721662 -0.1348185369 4.1375764833
C -0.1598279234 1.7008049609 6.2370452900
H 0.5667793800 3.1196771040 4.8131223494
C -0.5383349153 -0.5211476729 5.4120938992
C -0.5583981014 0.3866926590 6.4545994633
H -0.1779727224 2.4158028017 7.0490096055
H -0.8536176296 -1.5430653996 5.5844158852
H -0.8871612184 0.0740949668 7.4369998277
C -3.0728948489 -0.8617146506 -1.7692179323
H -3.4522482561 0.1191108930 -2.0541281333
H -2.5594832072 -1.2820242835 -2.6374895910
H -3.9325420923 -1.4978465143 -1.5453203017
B 1.9302311707 -1.4692168156 -0.4760845631
H 2.0926992353 -0.7549809971 0.4880071329
O 1.5561323980 -2.7784318104 -0.3712240999
O 2.3097091400 -1.0901139913 -1.7254161848
C 1.4970344620 -3.3128796410 -1.7243698167
C 2.3656960922 -2.2928747867 -2.5372996533
C 2.0423213774 -4.7284488508 -1.7005910976
H 1.3773165693 -5.3597260314 -1.1111830133
H 2.0887100056 -5.1354409704 -2.7125241252
H 3.0355829477 -4.7717805201 -1.2592911468
C 0.0334314614 -3.3355603437 -2.1396630283
H -0.0854662299 -3.7911926341 -3.1234532925
H -0.5267502995 -3.9274756357 -1.4153202035
H -0.3952363328 -2.3365082493 -2.1642889921
C 3.8355271028 -2.6846532438 -2.6141058751
H 3.9834114754 -3.5579687459 -3.2499472260
H 4.3963652647 -1.8502712448 -3.0343021289
H 4.2373341508 -2.8997897299 -1.6232875717
C 1.8324794492 -1.9635065231 -3.9179112107
H 2.4969907728 -1.2459243984 -4.3988046054
H 1.7910767922 -2.8612471443 -4.5379423707
H 0.8386630357 -1.5227963688 -3.8673508454
次に、エネルギープロファイルを示す。
エネルギープロファイルに関して、STEP 1はbackward structureのエネルギーを示し、最後のSTEPはforward structureのエネルギーを示す。経路の途中でエネルギーがほぼ変化しない箇所が存在しないため、単一の素過程を示すIRC経路であることがわかる。エネルギーがほぼ変化しない箇所が存在した場合、2つの素過程がIRC経路中に含まれており、片方の素過程の活性化障壁がないことを示している可能性がある。
IRC経路の各末端構造とエネルギープロファイルから、今回求めた遷移状態構造が、σ結合メタセシスの遷移状態構造であることが示唆された。 σ結合メタセシスの遷移状態構造であることを今回よりもさらに強く主張するには銅の電荷のIRC経路上での変化が１価でほとんど変わらないことを示す必要があるが、今回はここまで立ち入らない。UMAモデルのNNPで部分電荷を調べる方法がわからないからである。
終わりに
自作ライブラリで、UMAモデルのニューラルネットワークポテンシャル(NNP)を用いて、Cu-DPPBZ系錯体触媒が関わる有機金属反応のある1つの遷移状態構造を算出する手順を説明した。
参考
- https://github.com/ss0832/MultiOptPy (自作ライブラリMultiOptPyのレポジトリ)
- https://avogadro.cc/ (Avogadro、分子構造可視化ツール)
- https://ai.meta.com/blog/meta-fair-science-new-open-source-releases/ (UMAの公開に関する記事)
- https://github.com/facebookresearch/fairchem (FAIR Chemistryの提供するGitHubのレポジトリ)
- https://fair-chem.github.io/ (同上のレポジトリの内容に関して説明したサイト)
- https://huggingface.co/facebook/UMA (NNPの配布サイト, Hugging Faceへのアカウント登録と配布元の使用許諾が必要である。)
- arXiv preprint arXiv:2505.08762 (2025). (プレプリント)
- 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.
【計算化学】自作pythonライブラリで遷移状態構造を求めてみる(Cu-DPPBZ系錯体触媒が関わる有機金属反応, NNP使用）
https://ss0832.github.io/posts/20250902_mop_usage_100/