最終更新：2026-02-01
概要
本記事では、自作ライブラリ(MultiOptPy)で、cubaneからcuneaneの間の反応経路探索を行ってみる。計算レベルは、Grimmeらが開発したGFN2-xTBとした。
MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonライブラリである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
使用した自作ライブラリMultiOptPyのバージョン
v1.20.7Video Demo
https://www.youtube.com/watch?v=AE61iY2HZ8Y
環境
Ubuntu (WSL2) 環境構築手順(Linux)
すでにAnacondaの導入が終わっている前提で、環境の構築手段について述べる。以下を実行すると可能である。
## 1. Download and install MultiOptPy:
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.20.7.zip
unzip v1.20.7.zip
cd MultiOptPy-1.20.7
## 2. Create and activate a conda environment:
conda create -n <env-name> python=3.12 pip
conda activate <env-name>
pip install git+https://github.com/ss0832/MultiOptPy.git@v1.20.7
# or pip install multioptpy
## 3. Copy the test configuration file and run the AutoTS workflow（任意で行う）:
cp test/config_autots_run_xtb_test.json .
python run_autots.py aldol_rxn.xyz -cfg config_autots_run_xtb_test.json
手順
1. 初期構造の準備
モデル反応系として、以下の構造を用意した。今回はファイルの名前をcubane.xyzとした。
初期構造は以下のものを使用した。
16
OptimizedStructure
C -0.585085285781 1.101643739045 -0.479348548609
H -1.059906825489 1.995151768396 -0.867642210138
C 0.599769947911 1.080552692170 0.509271370016
H 1.085440806291 1.956793085961 0.922295844461
C -1.291757906203 0.003746076987 0.343496737012
H -2.338899606165 0.006917460717 0.622296322739
C -0.106724163681 -0.017511316578 1.332033128149
H -0.193317595686 -0.031477527880 2.412244082723
C 1.291668797844 -0.003954958211 -0.343356339144
H 2.338955522192 -0.007256881173 -0.622090064687
C 0.106959569422 0.017566481568 -1.332265351526
H 0.193218934420 0.031812091312 -2.412370587674
C 0.585207172762 -1.101958351805 0.479507025912
H 1.059917632635 -1.995389536038 0.867751181898
C -0.599454517650 -1.080401598693 -0.509301871205
H -1.085992482821 -1.956233225778 -0.922520719928
生成物側の初期構造を以下に示す。
cubane_prod.xyz
16
OptimizedStructure
C -0.756708276717 0.777485668475 -0.595620012418
H -1.373299977193 1.538967681171 -1.055729133524
C 0.292101264937 1.097181733576 0.492863800715
H 0.470691637331 2.101130265502 0.856093757631
C -1.396318326634 -0.244959261371 0.350686854504
H -2.440801553044 -0.493049061847 0.433311036109
C -0.376034071069 0.066116704368 1.409671235164
H -0.467484545233 0.109077366393 2.481626297904
C 1.358197668884 0.398889086748 -0.359149334008
H 2.399660923060 0.659138540487 -0.442949483332
C 0.337800371266 0.088082743673 -1.418211797236
H 0.426094971323 0.057382838860 -2.490781239163
C -0.251135068739 -1.189456183804 0.591530795552
H -0.408811074915 -2.161305942291 1.028984852227
C 0.791358422241 -0.989325915721 -0.472002509784
H 1.394687634501 -1.815356264217 -0.810325120339
2. 遷移状態構造最適化
run_autots.pyを適切に使用することで、自動的に遷移状態構造が得られる。以下にその手順を説明していく。
初期構造をMultiOptPy-1.20.7ディレクトリ内にcubane.xyzとして保存する。その後、同じディレクトリ内で、config_cubane.jsonを作成し、以下のように記述する。
config_cubane.json
{
"work_dir": "cubane",
"top_n_candidates": 5,
"multioptpy_version": "1.20.7",
"step1_settings": {
"usextb": "GFN2-xTB",
"opt_method": ["rsirfo_block_fsb"],
"calc_exact_hess": 50,
"spin_multiplicity": 1,
"electronic_charge": 0,
"tight_convergence_criteria": true,
"model_function": ["bitss", "../cubane_prod.xyz"],
"manual_AFIR": ["0.0", "1", "2"]
},
"step2_settings": {
"usextb": "GFN2-xTB",
"NSTEP": 40,
"use_model_hessian": "fischerd3",
"save_pict": true,
"node_distance": 0.30,
"align_distances_energy_predicted": 2,
"spin_multiplicity": 1,
"electronic_charge": 0
},
"step3_settings": {
"usextb": "GFN2-xTB",
"opt_method": ["rsirfo_block_bofill"],
"calc_exact_hess": 5,
"tight_convergence_criteria": true,
"max_trust_radius": 0.20,
"min_trust_radius": 0.02,
"frequency_analysis": true,
"NSTEP": 500,
"detect_negative_eigenvalues": false,
"spin_multiplicity": 1,
"electronic_charge": 0
},
"step4_settings": {
"usextb": "GFN2-xTB",
"opt_method": ["rsirfo_block_bofill"],
"spin_multiplicity": 1,
"electronic_charge": 0,
"calc_exact_hess": 10,
"tight_convergence_criteria": true,
"frequency_analysis": true,
"intrinsic_reaction_coordinates": ["0.5", "200", "lqa"],
"step4b_opt_method": ["rsirfo_block_fsb"]
}
}
その後、以下のコマンドを実行する。
python run_autots.py cubane.xyz -cfg config_cubane.jsonこれにより、これまでの似た内容の記事で行ってきたコマンドの操作をまとめ、遷移状態構造を求める処理を自動的に行う。
具体的な処理の流れは、
Step1. バイアスポテンシャルによるNEB法のための初期経路の作成
Step2. NEB法による経路の緩和
Step3. NEB法により得られた経路のエネルギー極大値を示す構造のうち、エネルギー値が上位の最大で3個
(`run_autots.py`にて、`--top_n X`で最大値を変更可能)の構造を初期構造とした遷移状態構造の算出
(Step4.得られた遷移状態構造に対するIRC計算とIRC経路の末端に存在する構造に対する構造最適化。
こちらは`--run_step4`をコマンドで追記しなければ行わない。)となっている。
run_autots.pyのオプションの説明：
-cfg YYY.jsonは、workflowを実行するためのオプションが記されたJSONファイルの読み込み先を指定する。
これらの一連の結果は、(jsonファイルの"work_dir"にて指定した名前)のディレクトリの中に存在するファイルを開いて確認できる。
以下にすべてのstepで共通のオプションに関する説明を載せる。
"opt_method": ["rsirfo_block_fsb"]は準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。(以前のHessian更新法とは細かな点で異なる方法を使用している。具体的には、複数の座標変位や勾配変位を用いてHessianの更新を行う。)"spin_multiplicity": Zはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。（デフォルトでは1が指定される。）"electronic_charge": 0は形式電荷をMとすることを示す。（デフォルトでは0が指定される。）"othersoft": "uma-s-1p1"は今回使用するNNPを指定している。これを使用する際にASEライブラリが必要である。"use_model_hessian": "fischerd3"は、計算コストが非常に低い数式を使用して、近似したHessianを生成する機能を呼び出すオプションである。デフォルトではこの機能は使用されない。"calc_exact_hess": 50: 50回の反復計算ごとに1回正確なHessianを計算する。
※オプションの説明はMultiOptPy-v1.20.7/OPTION_README.mdにて示されている。
Step 1
Step1では、omolのデータセットを使用したuma-s-1p1モデルのNNPで得たエネルギーに対して、指定した人工力ポテンシャルを加えた上で初期構造を構造最適化を行っている。
以下のJSON内で記述したバイアスポテンシャルで、次の経路緩和アルゴリズムの初期経路として用いるトラジェクトリーを生成する。
"manual_AFIR": ["yyy", "a", "b]：yyykJ/molの活性化障壁を超えうるペア同士を近づける力を原子のラベル番号aとbのペアに構造最適化時に加えることを示す。(プログラムの仕様上yyyを0.0にして指定しないとプログラムが作動しない。)"model_function": ["bitss", "../cubane_prod.xyz"],: BITSS法と呼ばれる原系と生成系間をつなぐパスを生成する方法を用いることを示す。あとのファイル名は生成系の構造の情報が入ったxyzファイルを指定している。（相対パスで指定する場合、必ず初めに、../をつけておかないと仕様上動かない。）"shape_conditions": ["yyy", "lt", "a,b"]: 構造最適化中に、ラベル番号aの原子とラベル番号bの原子の間の距離yyy(Å)よりも大きくなった時に構造最適化を途中で打ち切る。“lt”を”gt”にすると、yyy(Å)よりも小さくなった時に途中で打ち切る。"shape_conditions": ["yyy", "lt", "a,b,c"]: 構造最適化中に、角度の中心がラベル番号bの原子で、ラベル番号aの原子とラベル番号cで作る角度yyy(degrees)よりも大きくなった時に構造最適化を途中で打ち切る。“lt”を”gt”にすると、yyy(degrees)よりも小さくなった時に途中で打ち切る。"shape_conditions": ["yyy", "lt", "a,b,c,d"]: 構造最適化中に、ラベル番号bとラベル番号cを軸とした、ラベル番号aの原子とラベル番号dの原子がなす二面角yyy(degrees)よりも大きくなった時に構造最適化を途中で打ち切る。“lt”を”gt”にすると、yyy(degrees)よりも小さくなった時に途中で打ち切る。"keep_pot": ["c", "d", "a,b"]:力の定数ca.u.で、平行距離dÅの調和ポテンシャルをa番とb番の原子ペアにかける。初期経路生成時に開裂を防ぎたい結合を保持するとき等に使用する。
Step1が正常終了していれば作成されたwork_dirディレクトリ中に、cubane_step1_traj.xyzが存在する。必要に応じて確認し、目的に沿った初期経路が得られているか確認する。もし想定とは異なる場合は、プロセスをkillして再度設定を見直してやり直す。
cubane_step1_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このcubane_step1_traj.xyzはStep2のNEB計算に使用している。
※cubane_step1_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
Step 2
Step2では、NEB法を用いることで、先ほど得られたcubane_step1_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
Step2固有のオプションについて以下に示す。
-
"NSTEP": nはn回分NEB法による経路の緩和を行うことを示す。 -
"save_pict": trueは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。 -
"node_distance": yyy: 入力された経路を経路座標上でyyy(Å)間隔で線形補間により構造を再配置して初期経路とする。 -
"align_distances_energy_predicted": a:経路緩和a回に1回、経路座標上で等間隔に線形補間により構造を置きなおす。エネルギー極大値を示す構造に対しては前後のノードの情報を使って経路座標を変数とした多項式を作り、真のエネルギー極大値の場所を連続最適化により推定し、再配置する。
MultiOptPy-v1.20.7/"work_dir"と同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。
そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
経路の緩和後の各ノードのエネルギー一覧(単位
※bias_force_rms.csvにて、各Iterationごとのすべてのノードの勾配のRMS値を確認できる。
経路緩和の結果、以下の構造がstep3の初期構造として自動的に用いられた。“work_dir”内のcubane_step3_TS_Opt_Inputs内に保存されたcubane_ts_guess_X.xyzにて確認が可能である。ts_guessの番号が小さい順にエネルギー値が高い構造を示すようになっている。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
cubane_ts_guess_1.xyz
16
0 1
C -0.731310357324 0.871780550407 -0.640419940066
H -1.248573154119 1.631451087813 -1.215111640235
C 0.226489920438 1.167694797187 0.530792012907
H 0.328338526193 2.136044059466 1.008960371729
C -1.490577595362 -0.121164024306 0.194296624613
H -2.540607088289 -0.313479103416 0.288507618837
C -0.422426665430 -0.044703386757 1.266518595853
H -0.648230408144 -0.025313607488 2.332858425323
C 1.380244261202 0.551532374815 -0.214942068515
H 2.400913918140 0.875100660425 -0.318440514346
C 0.416773256146 0.025487465600 -1.267328882817
H 0.631014337737 0.043747195852 -2.336541987850
C 0.396653056613 -1.194231022123 0.747415446731
H 0.564139136980 -2.067985094749 1.367652616959
C 0.219132266773 -1.315476871271 -0.613071603705
H 0.518026588446 -2.220485081454 -1.131145075419
Step 3
step3のオプションで、追加での説明を要するものを以下に示す。
"opt_method": ["rsirfo_block_bofill"]は遷移状態構造の最適化向けのoptimizerを指定することを意味する。準ニュートン法であるRS-I-RFO法を使用する。今回は-fcで正確なHessianを計算するようにしているので、初期Hessianは正確なHessianを使用するようになっている。(Bofill法によるHessianの更新法を細かい点で変更している。具体的には、複数の座標変位や勾配変位を用いてHessianの更新を行う。)"saddle_order": 1は一次の鞍点を求めることを指定する。（step3のデフォルトでは一次の鞍点を指定する。それ以外の値の指定は、プログラムの使用目的上想定していないので、行わないことを勧める。）"calc_exact_hess": 5は5回の反復回数当たり1回正確なHessianを計算することを指定する。"frequency_analysis": trueは収束条件を満たした後に基準振動解析を行うことを示す。（自前で実装しているため、あくまで目安として使用することを推奨する。各振動モードをvibration_animation内のxyzファイルで可視化できる。）UMAモデルから算出されるHessianは数値微分により求めているため、原子数Zが多いとZの二乗オーダーで計算コストが急増する。"tight_convergence_criteria": trueは収束条件を厳しくすることを示す。（Gaussianのtightと同等）"max_trust_radius": Dは一回の反復計算ごとの計算されるステップ幅の最大値をDÅ以下にすることを示す。デフォルトでは、"saddle_order": 1を指定すると0.1Åが指定される。"detect_negative_eigenvalues": trueは、初めの計算時(ITR. 0)に、任意の次数の鞍点(遷移状態構造等)を求める際に、正確なへシアンから算出した固有値に１つも負の固有値がない場合、計算を打ち切るオプションである。
実行して得られた正確な遷移状態構造と思われる構造を以下に示す。
(実行して得られた正確な遷移状態構造は計算開始時に、MultiOptPy-1.20.7/"work_dir"ディレクトリ内に生成された新規ディレクトリ内のcubane_ts_final_X.xyzとして保存されている。)
cubane_ts_final_1.xyz
16
OptimizedStructure
C -0.792203677524 0.811539043194 -0.698270998849
H -1.294092234889 1.549947695199 -1.314995182441
C 0.131337973633 1.165991818849 0.485189080983
H 0.154184711655 2.114082196475 1.001480127851
C -1.583941072671 -0.108305276488 0.222429143145
H -2.555219001869 0.252823383806 0.547869328240
C -0.450220108134 -0.106661305273 1.215446671908
H -0.631349894508 -0.127302813678 2.290766186295
C 1.285255187948 0.522861670572 -0.183450770983
H 2.348165142057 0.547461227113 -0.069315122206
C 0.371801054840 0.002909491323 -1.294292368852
H 0.634059362086 0.070389647136 -2.347898239734
C 0.486738992421 -1.160665390074 0.678574712913
H 0.826677974603 -1.956220384513 1.336404333525
C 0.302213882506 -1.372822330574 -0.680711856728
H 0.766591707846 -2.206028673068 -1.189225045065
19回の反復計算により、停留点に収束した構造が得られた。"frequency_analysis": trueオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。
以下に"frequency_analysis": trueオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -852.1765 170.3789 345.3664
Reduced mass [au] 5.1661 2.0750 2.2050
Force const [Dyne/A] -2.2104 0.0355 0.1550
Char temp [K] 0.0000 245.1372 496.9052
Normal mode x y z x y z x y z
C -0.00699 -0.01281 0.01080 -0.01158 -0.01770 -0.02118 0.00404 -0.02979 -0.04715
H -0.00237 -0.02271 -0.00165 0.01638 -0.01654 -0.04444 -0.04039 -0.09394 -0.08678
C 0.02316 0.00385 -0.02552 -0.02859 -0.00266 -0.00815 -0.06461 0.06762 -0.02135
H -0.03213 -0.02284 0.02277 -0.04375 -0.00781 -0.00303 -0.16765 0.06895 -0.02033
C -0.01186 0.02861 0.04307 -0.07059 0.10174 0.06341 0.04827 0.00074 0.02728
H 0.00266 0.04507 0.05007 0.09015 0.43552 0.19503 0.10752 0.12768 0.08764
C 0.00035 -0.01453 -0.00466 -0.01432 -0.02524 -0.01236 0.07986 0.00040 -0.01022
H 0.02750 -0.01904 0.00142 0.04062 -0.04271 -0.00381 0.12152 0.01157 -0.00273
C 0.01790 0.14748 -0.08723 -0.04202 -0.04877 0.04054 -0.03586 0.06127 0.07070
H 0.00036 -0.14020 0.08029 -0.04933 -0.22127 0.14547 -0.05802 -0.13972 0.30566
C 0.01064 -0.00410 0.00725 -0.01544 -0.01994 -0.00837 0.05314 -0.01147 0.03340
H -0.01711 -0.04776 -0.00306 -0.00658 0.03639 -0.00074 0.12855 0.01413 0.05310
C -0.11114 -0.12532 0.00313 0.03218 0.02438 -0.03609 0.02825 -0.03704 -0.04275
H -0.23208 -0.16539 0.00992 0.14314 0.09348 -0.00797 -0.04041 -0.06999 -0.04621
C 0.09543 0.00770 0.03935 0.10973 -0.03892 -0.04088 -0.08715 -0.03287 -0.02688
H 0.04508 0.00515 0.00470 0.29323 0.04563 -0.00567 -0.35994 -0.14333 -0.08826
(...snip...)その結果、虚振動が１つであることが確認できた。つまりこの構造は遷移状態構造である。
次に、vibration_animation内の虚振動を示す分子振動が示されたxyzファイル(mode_1_XXXi_wave_number.xyz)をAvogadroで確認可能である。
終わりに
自作ライブラリで、cubaneからcuneaneの間の反応経路探索を行った。なお、化学的な妥当性（この経路を通って異性化が起こるかどうか）の検討は行っていないので注意されたし。
参考
- 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.
- J. Chem. Phys. 157, 124107 (2022) (BITSS法について)
- J. Chem. Theory Comput. 2019, 15, 3, 1652–1671 https://doi.org/10.1021/acs.jctc.8b01176