最終更新:2025-09-18
概要
本記事では、自作ライブラリ(MultiOptPy)で、Pyrridineが配位したCu錯体触媒が関わるCMD(Concerted Metalation-Deprotonation)反応の遷移状態構造を算出してみる。使用する計算レベルは、Meta社のFAIR Chemistryが開発したニューラルネットワークポテンシャル(NNP)であるUMA(Meta’s Universal Model for Atoms)を使用する。
MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonライブラリである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
CMD機構(Concerted metalation-deprotonation mechanism)について
今回使用したニューラルネットワークポテンシャルについて:
- 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.5b
環境
Windows 11
※Windows 11環境下でAnaconda PowerShell Promptを使用した。
Source codeのダウンロード
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.15.5b.zip
unzip v1.15.5b.zip
cd MultiOptPy-v1.15.5b
環境構築手順
今回は、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.5b
ディレクトリ内に存在するsoftware_path.conf
に対して、uma-s-1.pt
の絶対パスを用いて以下を追記する。
uma-s-1::(uma-s-1.ptの絶対パス)
これで、MultiOptPy-v1.15.5b
がNNPuma-s-1
を使用できるようになる。
使用するNNPに関する具体的な説明
今回使用するNNPについて具体的に説明する。
- UMAのModel Checkpointは
uma-s-1
を使用した。 - 小分子系のトレーニングセットである
Omol25
(omol
)を使用して学習したニューラルネットワークポテンシャルを使用する。
※自作ライブラリでの具体的な使用の仕方に関しては、ase_calculation_tools.py
を参照
手順
1. 初期構造の準備
モデル反応系として、以下の構造を用意した。今回はファイルの名前をcmd_cu.xyz
とした。 初期構造は以下のものを使用した。
40
OptimizedStructure
Cu 0.190289951178 0.686511884538 -0.871022607760
C 2.606776328857 -1.168399922725 -1.339436780657
C 0.737346967964 -1.735063725517 -2.541489481949
C 3.355563859096 -2.172330627148 -1.926878718381
H 3.045349587432 -0.517752447158 -0.591912349649
C 1.415599742906 -2.762304436881 -3.179760088008
H -0.312304808842 -1.531204946846 -2.720409724879
C 2.746629491378 -2.984228056775 -2.871546324704
H 4.390369671408 -2.311442804334 -1.646273644689
H 0.896272918843 -3.374596401973 -3.903760082711
H 3.301815163429 -3.778525449812 -3.353874881084
C 2.339257602289 1.782779875983 0.792752685285
C 0.568416095689 1.040487134854 2.070239638346
C 3.066832267035 2.159975849325 1.903110539872
H 2.731897208852 1.936352878025 -0.205406922760
C 1.252227269038 1.382046826128 3.236755285417
C 2.509275644478 1.946465899666 3.155390942911
H 4.042320357548 2.608804549570 1.783518034408
H 0.790860703564 1.189219699886 4.194915181378
H 3.049090328604 2.214838289672 4.054302962129
C -0.803432241241 0.485165435918 2.127717526453
C -1.208299576024 -0.521603061789 1.256579864922
C -1.721216087957 1.005951858929 3.039176348589
C -2.515534299790 -0.983204226002 1.275018688768
H -0.505907543582 -0.965215259987 0.559771740616
C -3.022792123418 0.534453717777 3.067530119488
H -1.426035524334 1.809650061618 3.702469205890
C -3.422734593232 -0.458061145188 2.181830674473
H -2.817993840471 -1.733642692525 0.558432994740
H -3.731584602719 0.957175161453 3.767617387911
H -4.445069621722 -0.813391481011 2.192488967622
O -1.246126826917 1.409631813386 -1.986226479582
O -2.068768377404 -0.663287083877 -1.907597241032
C -2.128840673309 0.526992260984 -2.258055876265
C -3.313467760790 1.022909825317 -3.065178527195
H -3.852447857732 1.772596250901 -2.484992441209
H -2.959425461596 1.510819364163 -3.973310015720
H -3.981440882034 0.202989961502 -3.315451888719
N 1.121521588441 1.231156124030 0.861223710735
N 1.325709955086 -0.942720954075 -1.638258423001
2. 遷移状態構造最適化のための初期構造の算出
a. 初期構造をカレントディレクトリにcmd_cu.xyz
として保存し、以下を実行する。
python optmain.py .\cmd_cu.xyz -os uma-s-1 -opt rsirfo_fsb -modelhess -ma 250 1 22 250 25 33
-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
(今日の年月日)ディレクトリの中に存在するディレクトリを開いて確認できる。
正常終了していれば、このディレクトリ中に、cmd_cu_traj.xyz
が存在するので、これをコピーして、MultiOptPy-v1.15.5b
ディレクトリに置く。
cmd_cu_traj.xyz
は構造最適化の過程をAvogadro(公式ページ:https://avogadro.cc/ )等で可視化して確認できるようにしている。このcmd_cu_traj.xyz
は次のNEB計算に使用する。
※cmd_cu_traj.xyz
をアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
(この人工力ポテンシャルを加えて行った構造最適化の結果はcmd_cu_optimized.xyz
で確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-ma
の設定を見直してb.をやり直す。今回の場合以下の構造が得られた。生成系の安定構造に近いものが得られていると判断できるため、次のNEB計算を行うことが可能である。)
b. cmd_cu_traj.xyz
を初期パスとして、NEB法で経路の緩和を行う。
NEB法を用いることで、先ほど得られたcmd_cu_traj.xyz
全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。(この時点ではまだ正確な遷移状態構造は求められていない。)
python nebmain.py cmd_cu_traj.xyz -os uma-s-1 -ns 15 -aneb 3 5 -modelhess -spng -nd 0.5
-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.5b
と同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。 そのディレクトリ内のenergy_plot.csv
を確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
パスの緩和後の各ノードのエネルギー一覧(単位
※bias_force_rms.csv
にて、各Iterationごとのすべてのノードの勾配のRMS値を確認できる。
経路緩和の結果、経路のエネルギー極大値を示す構造の中から目視で、緩和後に得られた経路の11番(グラフ上では12番)の構造を遷移状態構造を求める初期構造として採用した。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
cmd_cu_traj_11.xyz
40
0 1
Cu 0.099497085750 0.448893682523 -0.519663449360
C 2.578229371168 -1.171873096938 -1.307064426908
C 0.744176372492 -1.744469320083 -2.540545228313
C 3.348561407221 -2.156694065338 -1.902519419352
H 3.001844271590 -0.529473677442 -0.544426959203
C 1.435965641811 -2.757455543741 -3.184080448661
H -0.299588650201 -1.561523993252 -2.751142036642
C 2.764635947277 -2.969860226822 -2.860471567134
H 4.383188361268 -2.282083564769 -1.614798196771
H 0.929452577751 -3.367396409642 -3.919401146206
H 3.331832908788 -3.756178054215 -3.341344204638
C 2.442752238101 1.837464615538 0.829945636697
C 0.594179710562 1.094955792088 2.033204876681
C 3.130419947818 2.188210768962 1.972440503072
H 2.861410792647 1.980861809775 -0.160786952186
C 1.253514843251 1.383063140480 3.231103984673
C 2.518297717342 1.938621488950 3.193380830824
H 4.113011573135 2.630668349824 1.902964077050
H 0.779989329053 1.163273424333 4.176866235853
H 3.032346094706 2.177915224279 4.115037169959
C -0.805561162927 0.567176894636 1.974582690601
C -1.312759471598 -0.218278513889 0.911716090473
C -1.634091946915 0.944537746172 3.031744470909
C -2.598209797525 -0.733798125275 1.117193825518
H -1.397408692129 -0.581992601323 -0.555611822693
C -2.937285071020 0.494959988016 3.137820576112
H -1.263434292765 1.648891672540 3.764426832426
C -3.399660847229 -0.409676240505 2.200739137864
H -3.012922703554 -1.393424241618 0.364653217342
H -3.577686490788 0.852920217279 3.931511481480
H -4.391294727043 -0.835760447214 2.283605876735
O -1.203500190177 1.330248343027 -1.939106603747
O -2.040626306442 -0.704660507887 -1.615833526598
C -2.077804588495 0.473893466529 -2.140535972520
C -3.264660045352 0.827128399907 -2.988208231161
H -3.840049434033 1.567093392706 -2.430215175777
H -2.909662401239 1.304016574979 -3.897781352510
H -3.895762577444 -0.028737072357 -3.212398501502
N 1.221310659371 1.307338601999 0.865788988158
N 1.297352545774 -0.958797892228 -1.612791280544
構造が壊れていないので、これを遷移状態を求めるための初期構造とする。
cmd_cu_traj_11.xyz
をMultiOptPy-v1.15.5b
と同じディレクトリ内にコピーする。
そして、以下を実行する。
python .\optmain.py .\cmd_cu_traj_11.xyz -fc 5 -order 1 -opt rsirfo_bofill -freq -tcc -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
として保存されている。)
40
OptimizedStructure
Cu 0.078481399687 -0.337642001068 -0.431667701288
C 1.698267561445 -2.792037698150 -0.694400994490
C 1.236061624025 -1.755112545255 -2.690042831170
C 2.401501781260 -3.801151279913 -1.327364828638
H 1.579563430938 -2.780856719205 0.381464914884
C 1.922087648290 -2.726073576411 -3.399783205659
H 0.749250804127 -0.922301629671 -3.181178623152
C 2.516362045096 -3.768448359184 -2.707965591969
H 2.845313471243 -4.594679281878 -0.742820644843
H 1.982624690897 -2.658785645642 -4.476966109315
H 3.058212453841 -4.543078283660 -3.234646902286
C 1.839141805561 2.192530051019 0.303489762442
C 0.262860959745 1.451702959231 1.818695920918
C 2.119198772726 3.322121492382 1.048575363477
H 2.347195604849 2.000640946685 -0.634134507553
C 0.483930785932 2.571372169768 2.619919645258
C 1.421273314721 3.509603143258 2.233690082639
H 2.856569733792 4.034068391292 0.705525000629
H -0.097906580437 2.709010995848 3.520256611033
H 1.598153404882 4.386858609039 2.842881770568
C -0.739812298117 0.407418585962 2.143410835761
C -1.359815327883 -0.262480964109 1.069883985918
C -1.032930704456 0.107627808265 3.473015396227
C -2.296469524496 -1.241282903298 1.419581311228
H -1.792351956652 0.512676956623 -0.169294517158
C -1.965440255964 -0.872632020538 3.772266027651
H -0.516545392017 0.617627100612 4.278606897765
C -2.602404045571 -1.547636898424 2.739314079476
H -2.822353104434 -1.771962773160 0.631789877202
H -2.190800759357 -1.110834185683 4.804234321616
H -3.338544155397 -2.309909587673 2.966414802656
O -0.570050414268 0.750631204993 -2.245895814514
O -2.339247195326 1.186750348871 -0.953081667712
C -1.665938771705 1.294513379580 -2.044556479032
C -2.306322621807 2.152091817892 -3.104782007475
H -2.453845018981 3.157867211912 -2.711973528336
H -1.684761147383 2.185379061696 -3.994930554492
H -3.290289811599 1.751520857444 -3.346857126284
N 0.945555585609 1.278233498115 0.677048138757
N 1.124222207184 -1.783340237565 -1.357721110741
停留点に収束した分子構造が得られた。-freq
オプションにより生成されたnormal_modes.txt
やvibration_animation
ディレクトリ内の振動モードのアニメーションを確認した。
以下に-freq
オプションで生成されたnormal_modes.txt
の一部を示す。
Mode 0 1 2
Freq [cm^-1] -810.6862 16.1136 20.0864
Reduced mass [au] 1.2500 5.8950 4.7276
Force const [Dyne/A] -0.4840 0.0009 0.0011
Char temp [K] 0.0000 23.1839 28.8998
Normal mode x y z x y z x y z
Cu 0.00349 0.00543 -0.01060 -0.03818 -0.01772 -0.00423 -0.01022 -0.04268 0.01225
C 0.00147 -0.00302 0.00039 0.04576 0.03436 0.00347 -0.00897 -0.03976 -0.01976
C -0.00357 0.00338 0.00103 -0.00976 0.00015 -0.00144 0.02248 0.01174 -0.00043
C 0.00114 -0.00111 0.00046 0.09315 0.06473 0.00769 0.00645 -0.01995 -0.03426
H 0.00596 -0.00559 0.00088 0.04886 0.03543 0.00381 -0.02788 -0.06893 -0.02159
C -0.00029 0.00191 0.00170 0.03525 0.02901 0.00255 0.03936 0.03365 -0.01402
H -0.01444 0.00734 0.01386 -0.04923 -0.02483 -0.00468 0.02787 0.02287 0.01314
C -0.00076 0.00118 -0.00006 0.08786 0.06210 0.00718 0.03130 0.01761 -0.03128
H 0.00192 -0.00116 -0.00050 0.13358 0.08998 0.01127 -0.00076 -0.03375 -0.04752
H -0.00103 0.00318 0.00197 0.02936 0.02568 0.00201 0.05850 0.06289 -0.01111
H -0.00158 0.00106 -0.00091 0.12438 0.08547 0.01039 0.04416 0.03407 -0.04227
C -0.00176 -0.00243 -0.00260 0.02277 -0.05188 -0.01166 -0.02478 -0.01473 -0.03170
C -0.00171 0.00810 -0.00258 0.01407 -0.01825 -0.00429 -0.01945 0.00406 -0.01740
C -0.00167 -0.00334 0.00039 0.05343 -0.05979 -0.01117 -0.04017 0.00698 -0.05874
H -0.01048 -0.00207 -0.00551 0.01311 -0.06219 -0.01487 -0.02097 -0.03138 -0.02627
C -0.00159 0.00062 0.00099 0.04393 -0.02434 -0.00389 -0.03418 0.02549 -0.04346
C 0.00012 -0.00035 -0.00036 0.06377 -0.04554 -0.00734 -0.04455 0.02728 -0.06450
H -0.00221 -0.00388 -0.00117 0.06832 -0.07642 -0.01370 -0.04836 0.00781 -0.07469
H -0.00243 0.00426 -0.00061 0.05210 -0.01199 -0.00053 -0.03806 0.04037 -0.04820
H -0.00026 0.00001 -0.00161 0.08730 -0.05068 -0.00679 -0.05625 0.04428 -0.08557
C 0.00403 0.00961 -0.00570 -0.00634 0.00278 0.00114 -0.01171 0.00279 0.00206
C -0.08070 0.06211 -0.00284 -0.02463 0.01215 0.00571 -0.00606 -0.01939 0.01277
C 0.00955 -0.00848 -0.01037 -0.00749 0.01200 0.00295 -0.01294 0.02576 0.00673
C -0.00739 -0.00321 -0.00623 -0.04393 0.03296 0.01190 -0.00183 -0.01805 0.02744
H 0.39335 -0.47695 0.62798 -0.01730 0.01135 0.00999 0.00936 -0.00894 0.01997
C 0.00005 -0.00029 -0.00779 -0.02696 0.03246 0.00934 -0.00829 0.02585 0.02182
H 0.01666 -0.00909 -0.01804 0.00628 0.00314 -0.00032 -0.01796 0.04383 -0.00147
C 0.01068 -0.01089 -0.00365 -0.04537 0.04309 0.01387 -0.00239 0.00390 0.03247
H -0.04080 0.05295 -0.01635 -0.05832 0.04242 0.01523 0.00244 -0.03468 0.03592
H -0.00967 0.01211 -0.00734 -0.02785 0.04016 0.01091 -0.00961 0.04306 0.02550
H 0.00950 -0.00719 -0.00675 -0.06067 0.05937 0.01879 0.00126 0.00404 0.04468
O 0.01854 -0.01335 0.00276 -0.03647 -0.02310 -0.00537 0.01729 -0.00065 0.02163
O 0.02409 -0.00999 -0.00471 -0.01584 0.01803 0.00852 0.03003 0.02574 0.03032
C -0.03204 -0.00448 0.04979 -0.02582 -0.00303 0.00097 0.03647 0.03314 0.03522
C 0.00766 -0.00152 -0.00550 -0.02222 -0.00283 -0.00116 0.07322 0.09278 0.06111
H 0.00435 -0.00028 -0.00930 -0.00631 0.00221 -0.00818 0.16431 0.09994 0.07787
H 0.00442 0.00345 -0.01180 -0.02863 -0.01824 -0.00622 0.05333 0.05727 0.04586
H 0.00651 0.00124 -0.00585 -0.02923 0.00794 0.00940 0.03491 0.17449 0.08301
N -0.00491 -0.00163 0.00138 0.00352 -0.03150 -0.00834 -0.01449 -0.01616 -0.01110
N -0.00244 -0.00681 0.00222 -0.00484 0.00258 -0.00097 -0.00115 -0.02428 -0.00303
(...snip...)
その結果、虚振動が1つであることが確認できた。つまりこの構造は遷移状態構造である。
次に、vibration_animation
内の虚振動を示す分子振動が示されたxyzファイル(mode_1_XXXi_wave_number.xyz
)をAvogadroで確認すると、想定される反応系と生成系をつなぐ方向に振動していることを確認できた。
終わりに
自作ライブラリで、UMAモデルのニューラルネットワークポテンシャル(NNP)を用いて、Pyrridineが配位したCu錯体触媒が関わるCMD反応のある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.