Home
3647 words
18 minutes
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる(Heck反応（β水素脱離過程）, NNP使用）
2025-08-14
Computational Chemistry
MultiOptPy
/
python

最終更新：2025-08-14

概要#

本記事では、自作モジュール(MultiOptPy)で、Heck反応（β水素脱離過程）の遷移状態構造を算出してみる。使用する計算レベルは、Meta社のFAIR Chemistryが開発したニューラルネットワークポテンシャル(NNP)であるUMA(Meta’s Universal Model for Atoms)を使用する。

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

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

Heck反応について：

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

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

v1.12

環境#

Windows 11

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

Source codeのダウンロード#

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

環境構築手順#

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

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

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

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

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

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

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

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

手順#

1. 初期構造の準備#

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

33
OptimizedStructure
C     -1.032100390838     -1.035550686512     -0.499635688220
C      0.010104465189     -0.523835962485     -1.531664219333
H     -0.421146697476      0.324490136927     -2.074452289930
H      0.133157384733     -1.321404712277     -2.277978583004
C     -2.368046916793     -1.235783351996     -1.235827531200
H     -2.689654798368     -0.364862803503     -1.812232342764
H     -2.270041885058     -2.054360409780     -1.958070686459
H     -3.169859022661     -1.507402859330     -0.546346963808
C     -0.622456251148     -2.427368202395     -0.009681490950
H     -1.322680643945     -2.806914845009      0.735047813124
H     -0.629254954731     -3.117181795013     -0.864166407865
H      0.365688844644     -2.443142233142      0.437886689012
Pd    -1.125799417674      0.259873864285      0.980052526621
Br     0.306058640853     -0.970893829846      2.656318969659
C      1.289621749332     -0.146589310382     -0.981732708807
C      1.536103930694      1.169107728609     -0.587457735922
C      2.305636875568     -1.079207424244     -0.784331944077
C      2.729016678327      1.533895108198      0.016785549082
H      0.777386052512      1.924725604234     -0.768410734187
C      3.500845678891     -0.723650421275     -0.177839782580
H      2.153156634569     -2.102775671521     -1.105886365114
C      3.714550358454      0.582721654604      0.234441976611
H      2.890465704150      2.562417789908      0.314631042052
H      4.267811725835     -1.472087152295     -0.025043526271
H      4.645533767041      0.859539003848      0.711547833672
P     -2.442032999849      1.615589963510     -0.239074071374
H     -2.560992252723      2.941717560530      0.206858808601
H     -2.179222205778      1.885594211469     -1.592780295468
H     -3.809948999356      1.317959980377     -0.348137654102
P     -1.046617235566      1.787279753945      2.871118363327
H     -1.688881882573      3.047000621707      2.919831876797
H     -1.483611731193      1.299268108234      4.109803064047
H      0.237209794935      2.221830580616      3.226426508832
初期パスを求めるための初期構造

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

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

python .\optmain.py .\heck_rxn_elimination.xyz -modelhess -os uma-s-1 -opt rsirfo_fsb -ma 300 3 13 -ns 200
  • -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ライブラリが必要である。
  • -ns 200は反復計算回数を最大200回とする。構造最適化を収束させる計算コストを削減するためにつけた。

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

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

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

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

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

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

33
OptimizedStructure
C     -0.609844255750     -1.201267224164     -0.559375413494
C      0.243525410109     -0.180826721368     -0.908911727726
H      0.030424493127      1.458420747710      0.541255416295
H     -0.126317082633      0.504113199690     -1.667142456405
C     -1.893669294923     -1.392845487763     -1.326503001412
H     -2.155792187147     -0.506095176862     -1.902212734457
H     -1.756132656172     -2.227249193416     -2.022383596984
H     -2.716366129561     -1.643451126692     -0.658980814775
C     -0.219208929905     -2.373365078298      0.303235596134
H     -1.072470369765     -2.672029403216      0.910259786144
H      0.051353586150     -3.223446526446     -0.331737436111
H      0.616583381472     -2.149545285577      0.963076827744
Pd    -0.955820227587      0.549975836041      0.920635433751
Br    -2.872201104624     -0.754525673147      2.036354706649
C      1.721863847592     -0.110209626207     -0.711947763452
C      2.355479252485      1.131563931005     -0.640555463217
C      2.516321452484     -1.253578113109     -0.723103306485
C      3.731888280389      1.227056389257     -0.539712929214
H      1.758121252181      2.035662540713     -0.663961413964
C      3.897850145770     -1.160957467070     -0.619002988190
H      2.060370379689     -2.226271331180     -0.835353844908
C      4.510670471769      0.077015631903     -0.518480941078
H      4.200016587970      2.201140487732     -0.480135690625
H      4.494845954671     -2.063774772257     -0.626108119407
H      5.587267417890      0.148698299031     -0.436814481305
P     -3.897844372615      1.720472286362     -0.448134218966
H     -4.586660556052      2.004313842073      0.749974517886
H     -4.765623098857      2.509946458444     -1.243315904848
H     -4.607357310514      0.525569523112     -0.686936885528
P     -1.288966920754      1.973749534258      2.642133941212
H     -0.538996114508      3.149940303229      2.802700674459
H     -2.581299127906      2.497002125793      2.762159591261
H     -1.132012174477      1.424797070420      3.919024641015
初期パスを求めるための構造最適化の結果（安定構造ではない）

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

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

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

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

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

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

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

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

私が実行した環境では、17番と19番のノード(グラフでは18番, 20番)がエネルギー極大値を示していた。今回は基質のC-H結合長が長い19番のノードを遷移状態構造を求める初期構造とした。

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

heck_rxn_elimination_traj_19.xyz

33
0 1
C      -0.863586642196     -1.239682668729     -0.645659788714
C      -0.002276440498     -0.137573541300     -1.030603623713
H      -0.548270471959      0.967419158343     -0.457669617410
H      -0.229456783969      0.293526515730     -2.005912786890
C      -2.098841809628     -1.441456689563     -1.501588186782
H      -2.407243268646     -0.528200408440     -2.013183930872
H      -1.909103904499     -2.204098975909     -2.265309277870
H      -2.930072776297     -1.796652211005     -0.892170341038
C      -0.385045703210     -2.522769826650     -0.023394584956
H      -1.205395041788     -2.973199371234      0.534072427680
H      -0.084792993712     -3.229115657107     -0.807188119089
H       0.435142730885     -2.374659070518      0.672513308867
Pd     -1.203935867580      0.261576843487      0.827528835124
Br     -1.421722057816     -1.185323358042      2.875566791146
C       1.460472519540     -0.038034366519     -0.679683007560
C       1.971275348967      1.209899497797     -0.327703074787
C       2.329133363202     -1.124201089516     -0.708544732728
C       3.295575765879      1.366805200346      0.037320866975
H       1.305515236971      2.066190276040     -0.316225449221
C       3.656812395493     -0.972176370249     -0.328048121301
H       1.974668670897     -2.098673105359     -1.012887447577
C       4.141722601217      0.266377146663      0.058411810047
H       3.663043125218      2.344126139654      0.322152760798
H       4.312261555887     -1.833300026390     -0.323408287712
H       5.170411843941      0.375826707082      0.375553122921
P      -3.598772044310      1.332563892208     -0.015911826355
H      -3.860437801296      2.618604525910      0.503196293459
H      -3.305251520172      1.817266945922     -1.310756064147
H      -4.948374984054      1.052362212010     -0.324809721950
P      -0.787953287952      2.088659968197      2.322939563535
H      -0.980369293296      3.469554121131      2.102660881710
H      -1.468954354201      1.987827458415      3.542081186993
H       0.523821888982      2.180530127594      2.816660141417
NEB法により緩和したパスのエネルギー極大値を示す構造

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

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

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

python .\optmain.py heck_rxn_elimination_traj_19.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と同等）

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

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

33
OptimizedStructure
C     -0.602531096756     -1.097671434272     -0.258582425362
C      0.242468392100     -0.088587336794     -0.817655534602
H     -0.015432249419      1.284500499236     -0.174733180217
H     -0.078243225117      0.274737751367     -1.791998051191
C     -1.825387876415     -1.512123666666     -1.040316701786
H     -2.211448870243     -0.697588499926     -1.653754801682
H     -1.568361792771     -2.350330132495     -1.698183362010
H     -2.614207133085     -1.841831120375     -0.364796087807
C     -0.044304207996     -2.184243658121      0.624137718191
H     -0.832032061272     -2.569373864631      1.270082262471
H      0.323032372733     -3.008628870413      0.001726003338
H      0.777074499139     -1.843432126277      1.249896101796
Pd    -1.027214670266      0.665480320530      0.851483226324
Br    -2.661742664137     -0.444419386110      2.382272884732
C      1.734259335206     -0.117115972121     -0.674058421450
C      2.388778072732      0.128015919108      0.527261486730
C      2.487733300391     -0.432081474413     -1.800740500365
C      3.770543789820      0.042246607067      0.605953721986
H      1.806131537618      0.394431182956      1.401982072113
C      3.870008678366     -0.515941015609     -1.724860628225
H      1.986006960505     -0.617062369868     -2.743566888322
C      4.514950896004     -0.280802470292     -0.519445628052
H      4.266755528144      0.233623625335      1.548564171479
H      4.443061539705     -0.762831616014     -2.609144730609
H      5.593473355443     -0.343517612170     -0.458498946095
P     -4.231708256692      1.474813167895     -0.437349627720
H     -4.871665457891      1.792407101582      0.779641065656
H     -5.300472521603      1.953203893053     -1.235254031505
H     -4.695910086074      0.142842098007     -0.444262570038
P     -1.321237826062      2.666259238206      2.018683682613
H     -0.715572268100      3.877946578494      1.650385128491
H     -2.637509356362      3.121687015406      2.165053361924
H     -0.949296637644      2.655387628325      3.370079229195
遷移状態構造

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

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

Mode                                 0                   1                   2
Freq [cm^-1]                     -542.8467             22.5170             36.8040
Reduced mass [au]                   1.4121              6.5505              6.3082
Force const [Dyne/A]               -0.2452              0.0020              0.0050
Char temp [K]                       0.0000             32.3969             52.9528
Normal mode                   x         y         z            x         y         z            x         y         z     
       C                 0.03967   -0.05090   -0.05883   -0.00530    0.01343   -0.03099   -0.02022   -0.02502   -0.00975
       C                -0.01120   -0.10935   -0.00638    0.01153    0.00738   -0.01662   -0.01970   -0.02881   -0.01738
       H                -0.13720    0.60834    0.49066    0.00157    0.00992   -0.02628   -0.04512   -0.02206   -0.05024
       H                -0.07611    0.14292    0.11680    0.03044    0.01249   -0.02098   -0.01208   -0.04467   -0.02590
       C                -0.01247   -0.00817   -0.00093   -0.00865    0.04169   -0.04097   -0.02307   -0.02519   -0.00568
       H                 0.01322   -0.00457   -0.01086    0.00270    0.05647   -0.02871   -0.02058   -0.02789   -0.01102
       H                -0.04831   -0.02729    0.01059   -0.01871    0.04850   -0.05356   -0.02898   -0.03129   -0.00016
       H                -0.01471    0.03023    0.01354   -0.01422    0.04135   -0.04759   -0.02341   -0.01663   -0.00172
       C                -0.00093   -0.01560    0.00663   -0.02579   -0.00610   -0.04229   -0.02212   -0.02511   -0.00865
       H                -0.00429    0.01741    0.01987   -0.03470   -0.00348   -0.05172   -0.02202   -0.02107   -0.00619
       H                -0.01116   -0.03858    0.03208   -0.03441   -0.00272   -0.05188   -0.02694   -0.02737   -0.00854
       H                 0.00337   -0.00383   -0.00701   -0.02443   -0.02664   -0.03307   -0.01969   -0.02794   -0.01065
      Pd                 0.00711    0.01769   -0.00287    0.00576    0.00997   -0.02206   -0.01924   -0.01775   -0.02207
      Br                -0.00089   -0.00042    0.00052    0.02773    0.00276   -0.00368    0.03012   -0.01166    0.03245
       C                -0.01677    0.00997    0.01492    0.00859   -0.00540    0.00962   -0.02142   -0.00757    0.00232
       C                -0.00072    0.00225    0.00755   -0.01308   -0.02563    0.02575   -0.04584   -0.01635    0.01796
       C                -0.01010    0.01159    0.00605    0.02826    0.00196    0.02084    0.00391    0.03292    0.00832
       C                 0.00122    0.00104   -0.00024   -0.01555   -0.03920    0.05253   -0.04554    0.01423    0.03883
       H                 0.00867   -0.00784    0.01428   -0.02829   -0.03152    0.01745   -0.06585   -0.04509    0.01342
       C                -0.01015    0.00472    0.00055    0.02608   -0.01151    0.04766    0.00496    0.06468    0.02893
       H                -0.01326    0.00235    0.00914    0.04565    0.01778    0.00849    0.02383    0.04118   -0.00388
       C                -0.00228    0.00104   -0.00037    0.00403   -0.03240    0.06348   -0.01985    0.05532    0.04384
       H                 0.00301   -0.00461    0.00006   -0.03281   -0.05540    0.06490   -0.06533    0.00711    0.05069
       H                -0.01483   -0.00131   -0.00094    0.04160   -0.00592    0.05617    0.02510    0.09719    0.03290
       H                -0.00240   -0.00497   -0.00352    0.00222   -0.04327    0.08438   -0.01927    0.08037    0.05954
       P                 0.00024   -0.00076    0.00107   -0.08354   -0.03095    0.07952    0.05577    0.06574    0.00381
       H                -0.00648   -0.00157   -0.00166   -0.07512   -0.07308    0.09489    0.06592    0.05773    0.01128
       H                 0.00791    0.00642   -0.00388   -0.11120   -0.05137    0.10408    0.05932    0.09126    0.01481
       H                -0.00344    0.00135   -0.00289   -0.04076   -0.04572    0.06190    0.03746    0.07220   -0.01099
       P                -0.00415   -0.02011   -0.00024   -0.00103    0.01484   -0.03166    0.01543    0.00259   -0.04612
       H                 0.02002   -0.05542   -0.05833   -0.02015    0.01916   -0.04879   -0.03944    0.00649   -0.12295
       H                -0.00087   -0.01319   -0.00883   -0.00454   -0.00005   -0.01752    0.01723   -0.02236    0.04779
       H                 0.00216   -0.00447   -0.00428    0.01755    0.02954   -0.03661    0.12198    0.04518   -0.07484

(...snip...)

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

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

終わりに#

　 　自作モジュールで、UMAモデルのニューラルネットワークポテンシャル(NNP)を用いて、Heck反応（β水素脱離過程）のある1つの遷移状態構造を算出する手順を説明した。

参考#

【計算化学】自作pythonモジュールで遷移状態構造を求めてみる(Heck反応（β水素脱離過程）, NNP使用）
https://ss0832.github.io/posts/20250814_mop_usage_49/
Author
ss0832
Published at
2025-08-14
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる(Heck反応, 挿入素過程, NNP使用）
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる(SN2反応（基質-chloromethane, 求核剤-bromide ion）, NNP使用）