【計算化学】自作pythonモジュールで遷移状態構造を求めてみる(Claisen転位, NNP使用）

最終更新：2025-08-09

概要#

本記事では、自作モジュール(MultiOptPy)で、Claisen転位の遷移状態構造を算出してみる。以前同じ遷移状態構造を計算したので、今回は、別の計算レベルで算出する。使用する計算レベルは、Meta社のFAIR Chemistryが開発したニューラルネットワークポテンシャル(NNP)であるUMA(Meta’s Universal Model for Atoms)を使用する。

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

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

Claisen転位について：

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

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へのアカウント登録と配布元の使用許諾が必要である。)

使用した自作モジュール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

環境構築手順#

今回は、WSL2の代わりに、Windows 11のPower Shellを使用した。初めに、NNPを使用できる環境が整ったAnaconda PowerShell Promptを用意する手順を説明する。

1, 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

fairchem-coreは、FAIR Chemistryが管理しているNNPを動作させるために必要なライブラリである。
aseはNNPに電子エネルギーを算出したい分子構造を渡すために必要なインターフェイスの役割を果たすために必要なライブラリである。

これで、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
本記事の執筆当時は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. 初期構造の準備#

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

14

C     -1.897439922880     -0.613803649185     -1.412116540079
H     -1.872193454868     -1.671977245468     -1.222919220472
H     -2.850918253474     -0.203267508112     -1.684120813893
C     -0.819063557695      0.134355553194     -1.321695098598
H     -0.862171525022      1.191862303747     -1.508180908564
C      0.548953558862     -0.396337492618     -0.984226220678
H      0.493132775086     -1.432142442997     -0.669287742820
H      1.205971909925     -0.326448551480     -1.838925001221
O      1.192366041333      0.411155137871      0.027873340469
C      0.645474815019      0.384016111213      1.297624450891
H     -0.389078365079      0.102802158635      1.349521531919
C      1.342363326529      0.711641413440      2.359675380422
H      2.372199194871      0.992795834752      2.276465141235
H      0.890403457393      0.715348377010      3.330311701388

初期パスを求めるための初期構造

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

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

python .\optmain.py .\claisen_rxn.xyz  -ma 300 1 12 -os uma-s-1p1 -opt rsirfo_fsb

-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ライブラリが必要である。

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

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

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

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

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

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

14
OptimizedStructure
C     -0.919758147665      0.098573692958      0.161372519491
H     -1.166904662393     -0.968738312451      0.189103091221
H     -1.878604911258      0.624942742402      0.208113539120
C     -0.268055112508      0.405436485347     -1.160182921513
H      0.003653923419      1.445880553364     -1.323124161730
C      0.010415769355     -0.489222308717     -2.095214159144
H     -0.237619372072     -1.536698583486     -1.964030075049
H      0.491611536436     -0.206719356076     -3.022073871171
O      2.237085170855      0.098696574202      1.740235402163
C      1.179962381759     -0.338377830948      1.375153365747
H      1.081178136081     -1.409220558701      1.093128480435
C     -0.097078705678      0.465493254332      1.283108654462
H      0.185099615813      1.521377752704      1.283141793090
H     -0.620985622144      0.288575895072      2.231268342878

初期パスを求めるための構造最適化の結果（安定構造ではない）

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

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

python .\nebmain.py claisen_rxn_traj.xyz -ns 10 -fc 10 -spng -nd 0.5 -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計算の結果の可視化

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

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

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

claisen_rxn_traj_18.xyz

14
0 1
C      -1.079855629287      0.268036119772      0.206178309299
H      -1.470591644257     -0.671717343607      0.565674858808
H      -1.870781623346      1.026723894611      0.189999615165
C      -0.442371694252      0.352835064025     -1.096678538135
H      -0.713428182529      1.192891579303     -1.725224480187
C       0.638119420424     -0.424549088075     -1.434630803772
H       0.631425125518     -1.509133564412     -1.281656442730
H       1.205775789109      0.012696403000     -2.245612654348
O       1.677146237888     -0.117056294449     -0.075451768629
C       1.012821706507     -0.496317920387      0.919812253520
H       0.770661487387     -1.570605826761      1.025344248324
C      -0.020724897224      0.409030723890      1.379189599172
H       0.203436317966      1.450309879594      1.308261702839
H      -0.541632413905      0.076856373497      2.264794100675

NEB法により緩和したパスのエネルギー極大値を示す構造

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

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

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

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

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

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

14
OptimizedStructure
C     -1.304150236813      0.143731525802     -0.098807271799
H     -1.493961657045     -0.888767476561      0.159663372537
H     -2.041091559730      0.847849206433      0.266120350533
C     -0.511988951487      0.450669567494     -1.184047659169
H     -0.448168660388      1.480290037623     -1.514607099733
C      0.522359133824     -0.421352223961     -1.518205044794
H      0.382684939541     -1.484159974446     -1.371279230439
H      1.221112127403     -0.162483379991     -2.299627413453
O      1.680096113007     -0.210496393137     -0.075809937452
C      1.016648485351     -0.449875306082      0.995497596184
H      0.992919313953     -1.482755265312      1.367455576120
C      0.138994484614      0.459553987579      1.536618907875
H      0.236436977474      1.506206903093      1.289346060126
H     -0.391890509705      0.211588791467      2.447681793464

遷移状態構造

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

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

Mode                                 0                   1                   2
Freq [cm^-1]                     -601.9758            175.0855            312.8339
Reduced mass [au]                  10.0764              2.8492              6.1525
Force const [Dyne/A]               -2.1514              0.0515              0.3548
Char temp [K]                       0.0000            251.9090            450.0982
Normal mode                   x         y         z            x         y         z            x         y         z
       C                 0.07611    0.00922    0.10995   -0.00012    0.09984   -0.00595   -0.11040   -0.00329   -0.11275
       H                -0.03415   -0.00303   -0.00859   -0.11938    0.12941    0.02242   -0.06808   -0.00644   -0.10076
       H                 0.02473   -0.00602    0.03607    0.06467    0.18991   -0.04796   -0.10008   -0.00196   -0.09864
       C                 0.01941   -0.00826   -0.01498    0.02665   -0.02100   -0.02330   -0.03371    0.02240   -0.04041
       H                 0.01123   -0.00213    0.00081    0.09181   -0.04975   -0.09946   -0.02514    0.02218   -0.03845
       C                -0.10345   -0.00171   -0.12121   -0.02221   -0.10134    0.03815   -0.08702   -0.00548   -0.07562
       H                 0.01192    0.00766    0.03286   -0.05438   -0.07802    0.18515   -0.07094    0.00706    0.00917
       H                -0.04475    0.00504   -0.06441   -0.03680   -0.23059   -0.01732   -0.04587   -0.02088   -0.04249
       O                 0.06663    0.00512    0.09947    0.00222    0.12124   -0.01741    0.08103   -0.01034    0.07220
       C                 0.02555   -0.00670   -0.01795    0.02808   -0.00622   -0.02975    0.06651    0.00719    0.05403
       H                -0.01594   -0.00822    0.00188    0.08547   -0.04071   -0.12079    0.11695    0.02294    0.11053
       C                -0.09995   -0.00036   -0.08515   -0.03066   -0.10199    0.03689    0.06595   -0.00663    0.08108
       H                 0.03907    0.01526    0.02726   -0.05642   -0.07091    0.15743   -0.00494   -0.00943    0.04047
       H                -0.06938    0.00320   -0.06441   -0.03087   -0.21702    0.00587    0.08690   -0.01837    0.08960

(...snip...)

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

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

※参考：全振動数（単位は`cm^-1`）#

基準振動数：36個

（併進と回転成分は取り除いてある。）

-601.9758            175.0855            312.8339

327.7458            393.2901            442.8391

472.6033            533.0068            746.3639

803.0586            869.2177            906.1828

963.7304            996.3106           1000.7841

1024.5573           1034.7009           1058.3318

1106.8463           1253.7325           1268.2601

1292.0947           1349.0815           1417.3071

1451.6365           1510.6263           1535.2538

1628.2460           3016.8645           3145.7071

3148.4397           3153.7466           3177.6811

3221.4359           3238.3078           3239.6346

終わりに#

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

参考#

https://github.com/pyscf/pyscf (PySCFのgithubのレポジトリ)
https://github.com/ss0832/MultiOptPy (自作モジュールMultiOptPyのレポジトリ)
https://avogadro.cc/ (Avogadro、分子構造可視化ツール)
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.

概要#

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

環境#

Source codeのダウンロード#

環境構築手順#

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

手順#

1. 初期構造の準備#

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

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

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

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

※参考：全振動数（単位はcm^-1）#

終わりに#

参考#

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

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

※参考：全振動数（単位は`cm^-1`）#