Home
3150 words
16 minutes
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる(Prilezhaevエポキシ化(過酸化物:ethaneperoxoic acid), NNP使用）
2025-08-28
Computational Chemistry
MultiOptPy
/
python

最終更新：2025-08-28

概要#

本記事では、自作モジュール(MultiOptPy)で、Prilezhaevエポキシ化 (過酸化物acid)の遷移状態構造を算出してみる。使用する計算レベルは、Meta社のFAIR Chemistryが開発したニューラルネットワークポテンシャル(NNP)であるUMA(Meta’s Universal Model for Atoms)を使用する。

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

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

Prilezhaevエポキシ化反応について：

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

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

v1.15.1b

環境#

Windows 11

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

Source codeのダウンロード#

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

環境構築手順#

今回は、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.1bディレクトリ内に存在するsoftware_path.confに対して、uma-s-1.ptの絶対パスを用いて以下を追記する。

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

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

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

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

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

※自作モジュールでの具体的な使用の仕方に関しては、ase_calculation_tools.py を参照

手順#

1. 初期構造の準備#

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

15

C     -1.726556204909     -1.233431691637      1.758616444547
H     -0.946997093090     -1.848384189984      1.325940304693
H     -2.742940614285     -1.600829172101      1.690216054990
C     -1.443538369814     -0.080534904891      2.345554819244
H     -2.217476369148      0.539657516592      2.780338129769
H     -0.424337184297      0.280504413050      2.409625086141
C      0.996056475689      0.320628773852     -1.112396999573
O     -1.245658595614      0.758912426721     -0.924895603655
H     -0.923162001404      0.282975489473     -0.126728522246
O     -0.070353288983      0.876819451205     -1.713350522778
O      0.923769615186     -0.208414572858     -0.032035049279
C      2.222798697096      0.460084209283     -1.958401587586
H      2.007403565936      0.964307045641     -2.895938603060
H      2.967125793993      1.021233647368     -1.395645945774
H      2.623865573645     -0.533528441716     -2.150898005435
初期パスを求めるための初期構造

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

経路緩和の結果得られた経路のエネルギー極大値を示す構造を遷移状態構造を求める初期構造とした。

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

epoxidation_Me_traj_10.xyz(グラフでは11番)

15
0 1
C      -1.794428599303     -1.053663636679      1.455970837119
H      -1.048982948538     -1.803221845400      1.215119098758
H      -2.805161344006     -1.233104108852      1.122333922744
C      -1.489170382649      0.037904262318      2.165409597906
H      -2.247647752086      0.764067783688      2.424184393295
H      -0.471504302150      0.252168316418      2.480026289051
C       0.906474806221      0.294744104830     -1.107516352778
O      -1.083025027069      0.264045035776     -0.103399479388
H       0.093407650003     -0.151854855256      0.367721457076
O      -0.275090591778      0.606080044947     -1.430809978389
O       1.134786602330     -0.235732780468      0.036117133613
C       2.052548969914      0.531121009788     -2.009028740758
H       1.735571592759      1.028560777831     -2.917644803435
H       2.783759790613      1.130795502448     -1.467444194501
H       2.508461535739     -0.431909611390     -2.231039180311
NEB法により緩和した経路から得られた遷移状態構造を求めるための初期構造

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

epoxidation_Me_traj_10.xyzをMultiOptPy-v1.15.1bと同じディレクトリ内にコピーする。

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

python .\optmain.py .\epoxidation_Me_traj_10.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として保存されている。)

15
OptimizedStructure
C     -1.860551909178     -1.105932757004      1.512944815453
H     -1.177730830683     -1.940503983229      1.594775177188
H     -2.812277397892     -1.280956467647      1.035738073694
C     -1.576570486477      0.088440852475      2.095547595269
H     -2.296848050433      0.892091824714      2.097124633083
H     -0.662717931451      0.232146166714      2.654768493049
C      1.067043058847      0.297017646753     -1.060230600079
O     -1.014788530736      0.115204249488      0.184787739146
H     -0.061955083434     -0.202803759052      0.418466836894
O     -0.137685698286      0.655538058120     -1.309391817468
O      1.382338202261     -0.283780039865     -0.002268458962
C      2.092914761093      0.593483774965     -2.117673627700
H      1.690673872833      1.249809680630     -2.884016654918
H      2.968669516654      1.040159284241     -1.650737151591
H      2.399486506882     -0.349914531302     -2.569835053057
遷移状態構造

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

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

Mode                                 0                   1                   2
Freq [cm^-1]                     -641.9247             29.8452             72.6761
Reduced mass [au]                  12.7821              1.0688              2.6486
Force const [Dyne/A]               -3.1033              0.0006              0.0082
Char temp [K]                       0.0000             42.9405            104.5647
Normal mode                   x         y         z            x         y         z            x         y         z
       C                -0.03268   -0.01090    0.08048   -0.01651    0.01317   -0.00434    0.11133   -0.05679    0.02568
       H                -0.01201   -0.00254    0.02437   -0.03880   -0.00724   -0.02676    0.24490    0.06004    0.10905
       H                -0.00233    0.00428    0.01085   -0.02758    0.04753    0.00512    0.16691   -0.23794   -0.01856
       C                -0.03977   -0.04270    0.06192    0.02417   -0.00443    0.01182   -0.10847    0.02652   -0.03797
       H                -0.00465   -0.00929   -0.00221    0.04663    0.01566    0.03442   -0.23927   -0.09045   -0.12452
       H                -0.01696   -0.01150    0.02108    0.03559   -0.03914    0.00213   -0.16646    0.20587    0.00984
       C                 0.02290   -0.01768    0.02854   -0.00671   -0.01452   -0.00999   -0.00309   -0.00601   -0.00325
       O                 0.07302    0.07194   -0.17918    0.00046    0.00974    0.00365    0.00619    0.04601    0.01838
       H                -0.07826   -0.01484    0.04459   -0.00345   -0.01437   -0.00979    0.00178    0.01990    0.00484
       O                -0.03729   -0.02084    0.05919   -0.00167    0.00804    0.00122    0.01681    0.08905    0.03920
       O                 0.02959   -0.00196   -0.00994   -0.01112   -0.03790   -0.02098   -0.01859   -0.06443   -0.03027
       C                -0.02187    0.00660   -0.00511    0.00746    0.02054    0.01347   -0.00542   -0.04208   -0.01548
       H                -0.02721    0.00583   -0.00296   -0.14078   -0.37911   -0.25001   -0.01140   -0.08353   -0.04761
       H                -0.02539    0.01093   -0.00403   -0.23366    0.52955   -0.01535   -0.01548   -0.00272   -0.03303
       H                -0.01952    0.00745   -0.00405    0.45748   -0.00927    0.38534    0.01625   -0.05922    0.03590

(...snip...)

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

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

終わりに#

　 　自作モジュールで、UMAモデルのニューラルネットワークポテンシャル(NNP)を用いて、Prilezhaevエポキシ化 (過酸化物acid）のある1つの遷移状態構造を算出した。

参考#

【計算化学】自作pythonモジュールで遷移状態構造を求めてみる(Prilezhaevエポキシ化(過酸化物:ethaneperoxoic acid), NNP使用）
https://ss0832.github.io/posts/20250828_mop_usage_93/
Author
ss0832
Published at
2025-08-28
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる(Prilezhaevエポキシ化(過酸化物:peroxyformic acid), NNP使用）
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる(Prilezhaevエポキシ化(過酸化物:perbenzoic acid), NNP使用）