Home
3157 words
16 minutes
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる(Malaprade反応, NNP使用）
2025-08-11
Computational Chemistry
MultiOptPy
/
python

最終更新：2025-08-11

概要#

本記事では、自作モジュール(MultiOptPy)で、Malaprade反応の遷移状態構造を算出してみる。以前の記事でCriegee反応と呼ばれる同じ種類の生成物を与える反応を扱った。使用する計算レベルは、Meta社のFAIR Chemistryが開発したニューラルネットワークポテンシャル(NNP)であるUMA(Meta’s Universal Model for Atoms)を使用する。

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

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

Malaprade反応について：

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

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

v1.10c

環境#

Windows 11

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

Source codeのダウンロード#

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

環境構築手順#

今回は、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-1p1.ptをダウンロードする。（使用許諾が下りていれば可能である。）

https://huggingface.co/facebook/UMA

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

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

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

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

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

  • UMAのModel Checkpointはuma-s-1p1を使用した。
IMPORTANT

本記事の執筆当時(2025/8/9)はuma-s-1に関して、https://huggingface.co/facebook/UMA より、

” *Note uma-s-1 is now deprecated and will be removed by the next release. uma-s-1p1 should be used instead for a replacement.

とあり、今後のアップデートで削除されるため、使用しなかった。

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

※自作モジュールでの具体的な使用の仕方に関しては、https://github.com/ss0832/MultiOptPy/blob/main/multioptpy/ase_calculation_tools.py を参照

手順#

1. 初期構造の準備#

モデル反応系として、以下の構造を用意した。今回はファイルの名前をcriegee_rxn_2.xyz(Malaprade反応であるため、malaprade_rxnとするべきだが、今回はファイルの名前の付け方を誤っている。)とした。 初期構造は以下のものを使用した。

13

I     -0.319017518659      0.882564446429     -0.960881075624
O      0.930389162729      1.088565695032     -2.386880626457
H      1.528111083311      1.824423299579     -2.181158243652
O      1.049539243166     -0.258278745298     -0.237980140764
O     -1.201724722335      0.429892206736      0.709344949523
O     -0.343479995207      2.567797679245     -0.561919988547
O     -1.426333920487      0.029418947679     -1.972888863045
C      0.708032151212     -0.878413068314      1.013451317248
H      1.108106407246     -1.889664655640      0.981128845232
H      1.169077483324     -0.305355827305      1.816919182576
C     -0.808226097908     -0.864671310663      1.126969563124
H     -1.268090153128     -1.628615749485      0.495248586033
H     -1.126383123263     -0.997662917996      2.158646494352
初期パスを求めるための初期構造

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

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

python .\optmain.py criegee_rxn_2.xyz -os uma-s-1p1 -opt rsirfo_fsb -ma -400 8 11 -ns 20
  • -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-1p1は今回使用するNNPを指定している。これを使用する際にASEライブラリが必要である。
  • -ns Xは構造最適化のための反復計算回数を最大X回に制限する。今回は解離反応を扱うため、計算コスト削減のために指定した。

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

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

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

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

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

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

13
OptimizedStructure
I     -0.314934339591      0.882651549642     -1.116651790012
O      0.963372052484      1.168388497868     -2.545462687784
H      1.549722415442      1.871817291363     -2.224151046234
O      2.414815517611     -0.245948898055     -0.056748203652
O     -2.236670052475      0.331291426374      0.869325400341
O     -0.570748174324      2.570316785988     -0.748222707349
O     -1.572236999291      0.222707743847     -2.121747720674
C      2.484519149674     -0.962285306751      0.906108782517
H      2.168010778728     -2.022611429321      0.811480456947
H      2.311778468833     -0.517101511389      1.909283633039
C     -2.496682358041     -0.756491539530      1.315824199783
H     -2.499132649810     -1.625325043503      0.631161746341
H     -2.201813809239     -0.917409566534      2.369799936736
初期パスを求めるための構造最適化の結果（安定構造ではない）

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

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

python .\nebmain.py criegee_rxn_2_traj.xyz -ns 10 -modelhess -spng -nd 0.25 -os uma-s-1p1
  • -nd Nはノード間の距離をN Åとして初期パスを作成することを示す。
  • -ns nはn回分NEB法による経路の緩和を行うことを示す。
  • -fc MはM回あたりの経路緩和回数に対して1回だけ正確なHessianを計算し、経路緩和に使用する。これを使用すると、Hessianを使用しない場合の経路緩和アルゴリズムとは別のものを使用して、経路緩和を行う。
  • -spngは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。
  • -os uma-s-1p1は今回使用するNNPを指定している。これを使用する際にASEライブラリが必要である。

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

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

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

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

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

私が実行した環境では、4番のノード(グラフでは5番)がエネルギー極大値を示していた。

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

criegee_rxn_2_traj_4.xyz

13
0 1
I      -0.314645923875      0.888820557269     -0.984218838529
O       0.914245763776      1.104471791890     -2.425013878025
H       1.533183552087      1.814976181936     -2.195006952739
O       1.210992643064     -0.288064032044     -0.170323697875
O      -1.308136327616      0.367913143559      0.725790099745
O      -0.361792996545      2.573051288478     -0.563981398442
O      -1.476339287555      0.089595508020     -1.992750093494
C       0.996773956295     -0.857795125752      1.010625529807
H       1.193682884330     -1.923448089206      1.005795048935
H       1.405381735084     -0.248416044155      1.825480462569
C      -1.078105251537     -0.867984237998      1.124519784534
H      -1.482214593930     -1.641345046420      0.452129267505
H      -1.233026153578     -1.011775895577      2.186954666007
NEB法により緩和したパスのエネルギー極大値を示す構造

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

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

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

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

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

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

13
OptimizedStructure
I     -0.203800550398      1.017718482307     -1.005502592817
O      1.329586778262      1.340164314798     -2.105573198305
H      1.131105688408      2.149834462540     -2.604368876779
O      1.274548625852     -0.403624338466     -0.036823652524
O     -1.190223820419      0.473112174450      0.895734056194
O     -0.814250895513      2.645431327242     -1.111423057367
O     -1.125458479055     -0.022794692083     -2.035581163515
C      0.840774087274     -1.080234201257      0.984372915768
H      0.735473383783     -2.153070499383      0.812280083920
H      1.277522148333     -0.788246154840      1.940394108325
C     -0.927119105022     -0.749442393770      1.255708872103
H     -1.420959178723     -1.532461997926      0.675875197644
H     -0.907198682783     -0.896386483611      2.334907307353
遷移状態構造

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

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

Mode                                 0                   1                   2
Freq [cm^-1]                     -603.6403             54.7739             80.7870
Reduced mass [au]                   9.1907              5.0370              5.6560
Force const [Dyne/A]               -1.9731              0.0089              0.0217
Char temp [K]                       0.0000             78.8074            116.2345
Normal mode                   x         y         z            x         y         z            x         y         z   
       I                 0.00167    0.01413   -0.01138   -0.01325   -0.01035   -0.01614   -0.00342    0.00057   -0.00052
       O                 0.00233    0.00678   -0.01931    0.03601   -0.00719    0.05484    0.04478    0.05704    0.08450
       H                 0.01081   -0.01096   -0.04706    0.03646   -0.03768    0.00507    0.03208    0.04069    0.06266
       O                 0.04808   -0.06820    0.04831   -0.02112   -0.03274   -0.02094   -0.04728   -0.05329   -0.06972
       O                -0.04335   -0.04674    0.07485   -0.02012   -0.03969   -0.00190    0.03653    0.04989    0.04284
       O                -0.01509    0.01004   -0.02311    0.03728    0.01217    0.05803   -0.06827   -0.02710   -0.05319
       O                -0.00564   -0.00387   -0.00260   -0.02223    0.08953   -0.10916    0.06328   -0.03251   -0.02486
       C                 0.15008   -0.03638   -0.02387    0.04932    0.05956    0.06809   -0.02299   -0.05769   -0.06035
       H                -0.11365   -0.01709    0.01945    0.08573    0.03759    0.18298   -0.10132   -0.04607   -0.08361
       H                -0.02308    0.00230    0.04219    0.07646    0.17231    0.02167    0.06840   -0.11697   -0.08435
       C                -0.15003    0.03177    0.03970    0.04820    0.00417    0.08800    0.01794    0.05626    0.08071
       H                 0.01069   -0.01808   -0.03026    0.06581   -0.06669    0.16842   -0.07242    0.05141    0.16579
       H                 0.12110   -0.06207    0.02066    0.08745    0.08891    0.09968    0.10287    0.11151    0.08683

(...snip...)

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

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

終わりに#

　 　自作モジュールで、UMAモデルのニューラルネットワークポテンシャル(NNP)を用いて、Malaprade反応の遷移状態構造を算出する手順を説明した。

参考#

【計算化学】自作pythonモジュールで遷移状態構造を求めてみる(Malaprade反応, NNP使用）
https://ss0832.github.io/posts/20250811_mop_usage_40/
Author
ss0832
Published at
2025-08-11
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる(Tris(trifluoromethyl)phosphineが2つ配位したPd触媒によるtrifliorotoleneの還元的脱離, NNP使用）