最終更新：2025-11-08
概要
本記事では、自作ライブラリ(MultiOptPy)で、BH9データセットの4. Hydrogen atom transfer, No. 40の素過程の遷移状態構造を算出してみる。計算レベルは、Meta社のFAIR Chemistryが開発したニューラルネットワークポテンシャル(NNP)であるUMA(Meta’s Universal Model for Atoms)とした。
MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonライブラリである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
BH9のデータセットについて:
- J. Chem. Theory Comput. 2022, 18, 1, 151–166
https://doi.org/10.1021/acs.jctc.1c00694
この文献のSupporting Informationから、データセットの詳細を確認できる。
有機金属錯体が関わる反応を除いたさまざまなカテゴリの反応がまとめられたデータセットである。DFTの汎関数の電子エネルギーの精度の比較などのベンチマークに主に使われる。
今回使用したニューラルネットワークポテンシャルについて：
- 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). (プレプリント)
使用した自作ライブラリMultiOptPyのバージョン
v1.19.4b環境
Windows 11※Windows 11環境下でAnaconda PowerShell Promptを使用した。
Source codeのダウンロード(Unixコマンド)
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.19.4b.zip
unzip v1.19.4b.zip
cd MultiOptPy-v1.19.4bhttps://github.com/ss0832/MultiOptPy/releases/tag/v1.19.4b にアクセスしてzipファイルをダウンロードする。Unixコマンドの場合とはディレクトリ名が異なるので都度読み替えていただけると良い。
移動先のディレクトリでrequirements.txtを参照することで、本ソースコードで必要なモジュールを把握することが出来る。導入方法は各自の状況に合わせて適宜LLMとの対話などで調べると良い。
次に述べる環境構築手順を使用する場合は、環境構築が終わった後、pip install -r requirements.txtで本自作モジュールが動作させるために最低限必要なモジュールを導入することが可能である。
環境構築手順
今回は、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.74, 先ほど作成した仮想環境をconda activate (仮想環境名)で起動させる。
5, 以下のコマンドを実行し、必要なライブラリを導入する。
pip install ase==3.26.0 fairchem-core==2.7.1 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/UMA2, ダウンロード後、MultiOptPy-v1.19.4bディレクトリ内に存在するsoftware_path.confに対して、uma-s-1p1.ptの絶対パスを用いて以下を追記する。
uma-s-1p1::(uma-s-1p1.ptの絶対パス)これで、MultiOptPy-v1.19.4bがNNPuma-s-1p1を使用できるようになる。
使用するNNPに関する具体的な説明
今回使用するNNPについて具体的に説明する。
- UMAのModel Checkpointは
uma-s-1p1を使用した。 - 小分子系のトレーニングセットである
Omol25(omol)を使用して学習したニューラルネットワークポテンシャルを使用する。
※自作ライブラリでの具体的な使用の仕方に関しては、ase_calculation_tools.py を参照すると良い。omol以外のモデルを使用したい場合は、現バージョンでは、multioptpy/Calculator/ase_tools/firechem.py内の、self.task_nameを編集することで対応可能である。
手順
1. 初期構造の準備
モデル反応系として、以下の構造を用意した。今回はファイルの名前をbh9_4_40.xyzとした。 初期構造は以下のものを使用した。
27
OptimizedStructure
C -1.636233149698 -0.958910629225 -1.013054860345
C -2.071847689607 0.291195865734 -1.427470926590
C -2.132996952041 1.331004973678 -0.513790555818
C -1.779501143419 1.133969302156 0.817408876798
C -1.396976199707 -0.138593160604 1.229510331670
C -1.315692526426 -1.180420656211 0.319740851000
S -1.849955628364 2.492726004879 1.967081208705
H -2.337659362684 0.441278917471 -2.465821174241
H -2.442668413564 2.316146176219 -0.838161341644
H -1.147539627400 -0.308893833264 2.269121568363
H -0.987476599972 -2.162285877125 0.644321120258
H -0.533282254413 2.516935936800 2.244370723092
O 1.648716281011 1.889595511380 2.172833360909
C 1.741198608765 1.047079471879 1.259559826931
C 1.337692624767 1.350004589151 -0.098043383223
C 1.347316161986 0.393742516583 -1.067294983966
C 1.758852992801 -0.914575259472 -0.756372027678
C 2.236408216291 -1.238253209260 0.526893235800
C 2.232865314914 -0.293052994475 1.506785519425
H 2.569953303457 -0.515052384483 2.510614313348
H 1.007728964632 0.619407115678 -2.071939792333
H 0.999607340288 2.357805558681 -0.299837282512
H 2.573369239852 -2.249910238645 0.712423181441
O -1.514657183607 -1.935391506926 -1.965926736928
H -1.254682378428 -2.761211178162 -1.543225501422
O 1.697599002522 -1.914612446399 -1.660188845665
H 1.249861058044 -1.609728566041 -2.4595367053762. 遷移状態構造最適化
run_autots.pyを適切に使用することで、自動的に遷移状態構造が得られる。以下にその手順を説明していく。
初期構造をカレントディレクトリにbh9_4_40.xyzとして保存する。その後、同じディレクトリ内で、config_bh9_4_40.jsonを作成し、以下のように記述する。
config_bh9_4_40.json
{
"work_dir": "bh9_4_40",
"top_n_candidates": 3,
"multioptpy_version": "v1.19.4b",
"step1_settings": {
"othersoft": "uma-s-1p1",
"opt_method": ["rsirfo_block_fsb"],
"use_model_hessian": "fischerd3",
"spin_multiplicity": 2,
"electronic_charge": 0,
"manual_AFIR": ["300", "12", "13"]
},
"step2_settings": {
"othersoft": "uma-s-1p1",
"NSTEP": 15,
"ANEB": [3, 5],
"QSM": true,
"use_model_hessian": "fischerd3",
"save_pict": true,
"node_distance_bernstein": 0.80,
"align_distances": 9999,
"spin_multiplicity": 2,
"electronic_charge": 0
},
"step3_settings": {
"othersoft": "uma-s-1p1",
"opt_method": ["rsirfo_block_bofill"],
"calc_exact_hess": 5,
"tight_convergence_criteria": true,
"max_trust_radius": 0.2,
"frequency_analysis": true,
"spin_multiplicity": 2,
"electronic_charge": 0
},
"step4_settings": {
"othersoft": "uma-s-1p1",
"opt_method": ["rsirfo_block_bofill"],
"spin_multiplicity": 2,
"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 bh9_4_40.xyz -cfg config_bh9_4_40.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を生成する機能を呼び出すオプションである。デフォルトではこの機能は使用されない。
※オプションの説明はMultiOptPy-v1.19.4b/OPTION_README.mdにて示されている。
Step 1
Step1では、omolのデータセットを使用したuma-s-1p1モデルのNNPで得たエネルギーに対して、指定した人工力ポテンシャルを加えた上で初期構造を構造最適化を行っている。
以下のJSON内で記述したバイアスポテンシャルで、次の経路緩和アルゴリズムの初期経路として用いるトラジェクトリーを生成する。
"manual_AFIR": ["yyy", "a", "b]：yyykJ/molの活性化障壁を超えうるペア同士を近づける力を原子のラベル番号aとbのペアに構造最適化時に加えることを示す。
Step1が正常終了していれば作成されたwork_dirディレクトリ中に、bh9_4_40_step1_traj.xyzが存在する。必要に応じて確認し、目的に沿った初期経路が得られているか確認する。もし想定とは異なる場合は、プロセスをkillして再度設定を見直してやり直す。
bh9_4_40_step1_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このbh9_4_40_step1_traj.xyzはStep2のNEB計算に使用している。
※bh9_4_40_step1_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
Step 2
Step2では、NEB法を用いることで、先ほど得られたbh9_4_40_step1_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
Step2固有のオプションについて以下に示す。
"NSTEP": nはn回分NEB法による経路の緩和を行うことを示す。"align_distances": Xは線形補間で、各ノード間の距離を全て等しくするための処理である。X回の反復計算ごとに本処理を行う。Xを"NSTEP": nよりも大きな数値を指定することで、初期経路に対してのみ処理を行うことが出来る。"node_distance_bernstein": Nはノード間の距離をN Åとして初期経路を作成することを示す。経路作成時に元のノードをベルンシュタイン多項式を用いてがたついた経路を滑らかにする。
→プログラムの仕様上"align_distances": Xの処理を行った後に、"node_distance_bernstein": Nの処理を行うようになっている。
"save_pict": trueは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。"ANEB": [A, B]これを指定すると、(B+1)回の緩和ごとに、エネルギー極大値を示すノードと前後のノードの間に線形補間でA個の新規ノードを内挿するようにできる。デフォルトではこのような操作は行われない。このオプションを使用するとノードの数が徐々に増えるため、計算コストが使用しない場合と比べて増加する。一方で、エネルギー極大値を示すノード周辺にノードを増加させるため、緩和している経路中のノードが遷移状態構造付近に存在する可能性が高くすることが出来る。
MultiOptPy-v1.19.4b/"work_dir"と同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。 そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
経路の緩和後の各ノードのエネルギー一覧(単位
※bias_force_rms.csvにて、各Iterationごとのすべてのノードの勾配のRMS値を確認できる。
経路緩和の結果、以下の構造がstep3の初期構造として用いられることとなった。“work_dir”内のbh9_4_40_step3_TS_Opt_Inputs内に保存されたbh9_4_40_ts_guess_X.xyzにて確認が可能である。ts_guessの番号が小さい順にエネルギー値が高い構造を示すようになっている。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
bh9_4_40_ts_guess_1.xyz
27
0 2
C -1.640844977352 -0.965666776652 -1.009880903324
C -2.062217899644 0.301384412507 -1.413683488809
C -2.074507184028 1.335870942701 -0.505650190285
C -1.688996980690 1.135373593694 0.833878883624
C -1.353261598726 -0.168401192019 1.230358129603
C -1.319702351770 -1.207313508928 0.318929358695
S -1.540788162388 2.461755772208 1.945512005519
H -2.343076638594 0.456927447523 -2.446999380342
H -2.358513251832 2.331644518188 -0.819753556195
H -1.100247795321 -0.345816207896 2.267096957391
H -1.017249136716 -2.199860122101 0.634119226687
H 0.038942916471 2.239358001384 2.140934631648
O 1.224328823168 1.886132616075 2.195496123382
C 1.513537034836 1.046717201195 1.243932995004
C 1.227007879865 1.362317276212 -0.106747603371
C 1.327886993247 0.404352360508 -1.085726824074
C 1.753186358445 -0.879921873141 -0.749387087854
C 2.212517994624 -1.164160990740 0.544701049440
C 2.096895018019 -0.221054539702 1.529894339762
H 2.419575396345 -0.431540242764 2.542071503702
H 1.033786081078 0.627447949916 -2.105336470605
H 0.881208835243 2.363734604045 -0.329369299710
H 2.605849888723 -2.152027367847 0.750914787166
O -1.542712447602 -1.930289391868 -1.970845703931
H -1.284564405868 -2.763453569184 -1.560691844530
O 1.730018921566 -1.901793088130 -1.638519929232
H 1.261940688900 -1.621717825184 -2.435247709362Step 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Åが指定される。
実行して得られた正確な遷移状態構造と思われる構造を以下に示す。
(実行して得られた正確な遷移状態構造は計算開始時に、MultiOptPy-v1.19.4b/"work_dir"ディレクトリ内に生成された新規ディレクトリ内のbh9_4_40_ts_final_X.xyzとして保存されている。)
bh9_4_40_ts_final_1.xyz
27
OptimizedStructure
C -1.734278233437 -0.910562101733 -1.017072690674
C -2.118721034669 0.398048767963 -1.311086833878
C -2.055992599778 1.359978140726 -0.329622153407
C -1.638534941034 1.044090869904 0.977848881323
C -1.298221255884 -0.288574894116 1.254955990008
C -1.335967805051 -1.255648931096 0.267748972799
S -1.521159116104 2.267251034230 2.204043484636
H -2.439144491470 0.637289978663 -2.316581747999
H -2.320869886728 2.386208506786 -0.548314538828
H -1.003749952562 -0.548054428453 2.263163242851
H -1.045335517728 -2.277043089098 0.489614405725
H 0.056396983731 2.291549847964 2.140650341582
O 1.243297274774 2.110217980551 1.942916692647
C 1.445499648796 1.145152976547 1.089250097763
C 1.198180765265 1.330448685016 -0.292435443925
C 1.396085796809 0.302306411565 -1.182372144066
C 1.857733351229 -0.931627491626 -0.722399864454
C 2.150136855229 -1.130326018318 0.633127540488
C 1.930415774301 -0.116078259167 1.526668229135
H 2.124205170708 -0.243327143130 2.583582591359
H 1.178764593216 0.437753914740 -2.236386015391
H 0.838689587883 2.296891569696 -0.619576873914
H 2.522079317270 -2.095573158662 0.950946956606
O -1.759070160031 -1.811920242465 -2.036544306744
H -1.508444063213 -2.680691529705 -1.703328657656
O 2.033066851404 -1.986759925537 -1.553724105863
H 1.804937087075 -1.731001471248 -2.455072050121停留点に収束した構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。
以下に-freqオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -2346.0779 26.1165 49.7694
Reduced mass [au] 1.1762 7.5308 7.5920
Force const [Dyne/A] -3.8142 0.0030 0.0111
Char temp [K] 0.0000 37.5758 71.6070
Normal mode x y z x y z x y z
C 0.00080 -0.00669 -0.00556 -0.02106 0.02833 -0.02927 -0.07320 -0.00419 0.00570
C -0.00195 0.00913 0.00114 -0.01138 0.03820 0.00424 -0.05103 0.00128 0.00098
C -0.00426 0.00097 -0.01232 -0.00037 0.01330 0.02844 -0.02010 0.00765 -0.00719
C -0.00047 0.01931 0.01834 0.00266 -0.02234 0.01847 -0.00990 0.00961 -0.01043
C 0.00518 -0.01191 0.00086 -0.00227 -0.03061 -0.01678 -0.01673 0.00862 -0.00808
C 0.00129 0.00031 0.00771 -0.01459 -0.00605 -0.04014 -0.04522 0.00206 -0.00110
S -0.00790 -0.00681 -0.00724 0.00719 -0.05469 0.05136 -0.00060 0.01389 -0.01534
H 0.00298 0.00418 0.00076 -0.01506 0.06531 0.01190 -0.06726 -0.00217 0.00538
H -0.00007 0.00045 -0.00448 0.00331 0.02020 0.05639 -0.00916 0.00983 -0.01027
H 0.00593 -0.00247 -0.00246 0.00326 -0.05691 -0.02524 -0.00445 0.01131 -0.01093
H 0.00609 0.00570 0.00334 -0.01994 -0.01353 -0.06732 -0.05701 -0.00055 0.00223
H 0.89831 -0.14376 -0.10339 -0.00032 0.02200 -0.02465 -0.00175 -0.01040 0.00113
O -0.04423 0.03486 0.03125 -0.00356 0.07232 -0.08093 -0.00216 -0.02509 0.01699
C 0.00469 -0.04251 -0.03954 -0.00104 0.03812 -0.04316 0.00496 -0.01913 0.01205
C 0.00229 0.00203 0.02282 0.00349 -0.01415 -0.05035 0.01252 -0.01299 0.01150
C 0.00409 -0.01499 0.00060 0.01488 -0.04523 -0.01222 0.04621 -0.00044 0.00401
C -0.00060 0.01185 0.01133 0.01792 -0.02633 0.03553 0.07719 0.00804 -0.00344
C -0.00196 -0.00222 -0.01651 0.00636 0.02283 0.04464 0.04904 -0.00311 0.00112
C -0.00917 0.02210 -0.00037 -0.00186 0.05561 0.00604 0.01343 -0.01706 0.00816
H -0.00151 0.00588 0.00000 -0.00806 0.09498 0.01193 -0.00093 -0.02399 0.00996
H 0.00547 -0.00448 -0.00327 0.02150 -0.08435 -0.01867 0.06063 0.00671 0.00196
H 0.00554 0.00295 0.00544 0.00137 -0.02771 -0.08789 -0.00218 -0.01699 0.01589
H 0.00232 0.00049 -0.00452 0.00656 0.03545 0.08261 0.06604 0.00198 -0.00347
O -0.00080 0.00617 0.00683 -0.03571 0.05215 -0.05008 -0.13002 -0.01590 0.01789
H 0.01497 -0.00039 -0.00288 -0.05299 0.03988 -0.06917 -0.14394 -0.01900 0.01971
O 0.00095 -0.01001 -0.00826 0.03109 -0.05465 0.07363 0.14267 0.02853 -0.01500
H 0.01108 0.00570 0.00467 0.04876 -0.08361 0.06100 0.16508 0.03447 -0.01879
(...snip...)その結果、虚振動が１つであることが確認できた。つまりこの構造は遷移状態構造である。
次に、vibration_animation内の虚振動を示す分子振動が示されたxyzファイル(mode_1_XXXi_wave_number.xyz)をAvogadroで確認すると、求められた遷移状態構造の中に、想定される反応系と生成系をつなぐ方向に振動している構造が存在することを確認できた。
終わりに
自作ライブラリで、UMAモデルのニューラルネットワークポテンシャル(NNP, uma-s-1p1)を用いて、BH9データセットの4. Hydrogen atom transfer, No. 40の反応のある1つの遷移状態構造を算出する手順を説明した。
参考
- 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. Phys. Chem. A 2019, 123, 3761−3781 http://dx.doi.org/10.1021/acs.jpca.9b01546 (MOBH35のデータセットについて)