【計算化学】自作pythonライブラリで遷移状態構造を求めてみる(BH9データセットの9. Nucleophilic addition, No. 6, NNP(UMA)使用)

最終更新：2025-12-22

概要#

本記事では、自作ライブラリ(MultiOptPy)で、BH9データセットの9. Nucleophilic addition, No. 6の素過程の遷移状態構造を算出してみる。計算レベルは、Meta社のFAIR Chemistryが開発したニューラルネットワークポテンシャル(NNP)であるUMA(Meta’s Universal Model for Atoms)とした。

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

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

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

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). (プレプリント)

BH9のデータセットについて:

J. Chem. Theory Comput. 2022, 18, 1, 151–166

https://doi.org/10.1021/acs.jctc.1c00694

この文献のSupporting Informationから、データセットの詳細を確認できる。

有機金属錯体が関わる反応を除いたさまざまなカテゴリの反応がまとめられたデータセットである。DFTの汎関数の電子エネルギーの精度の比較などのベンチマークに主に使われることを想定していると考えられる。

使用した自作ライブラリMultiOptPyのバージョン#

v1.20.4

Video Demo#

https://www.youtube.com/watch?v=AE61iY2HZ8Y

環境#

Windows 11

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

Source codeのダウンロード(Unixコマンド)#

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

https://github.com/ss0832/MultiOptPy/releases/tag/v1.20.4 にアクセスして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.7

4, 先ほど作成した仮想環境を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の自作ライブラリへの導入方法について説明する。

UMAを使用可能にするための手順#

1, 以下のサイトにアクセスして、uma-s-1p1.ptをダウンロードする。（使用許諾が下りていれば可能である。）

https://huggingface.co/facebook/UMA

2, ダウンロード後、MultiOptPy-1.20.4ディレクトリ内に存在するsoftware_path.confに対して、uma-s-1p1.ptの絶対パスを用いて以下を追記する。

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

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

※Linuxの場合#

すでにAnacondaの導入が終わっている前提で、環境の構築手段について述べる。以下を実行すると可能である。

## 1. Download and install MultiOptPy:
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.20.4.zip
unzip v1.20.4.zip
cd MultiOptPy-1.20.4

## 2. Create and activate a conda environment:
conda env create -f environment.yml
conda activate test_mop

## 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


#### Installation via pip (Linux)
conda create -n <env-name> python=3.12 pip
conda activate <env-name>
pip install git+https://github.com/ss0832/MultiOptPy.git@v1.20.4 
# or pip install multioptpy
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.20.4.zip
unzip v1.20.4.zip
cd MultiOptPy-1.20.4

あとは、「UMAを使用可能にするための手順」と同じように進めれば問題なく導入可能である。

使用する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_9_6.xyzとした。初期構造は以下のものを使用した。

16
OptimizedStructure
C      2.215855740596     -0.747273864659     -0.627373706423
C      1.453690653488      0.196444770110      0.304728412069
C      0.436240158150     -0.325896821815      1.061585001682
O      1.789271976761      1.427567290020      0.240193794991
C     -1.766575989307      0.432442488294     -0.110954174508
C     -2.583152224315     -0.577073182601      0.091609620527
C     -1.202498901921      1.562522903990     -0.431218270948
H      3.283194622010     -0.694666968677     -0.395860302716
H      1.882852363492     -1.784507058018     -0.546922026179
H      2.092064850926     -0.407762724190     -1.659946679744
H     -0.085666562780      0.295090616212      1.778285869447
H      0.198187527536     -1.381010129986      1.029913683402
H     -0.155314458058      1.787215921301     -0.220752686804
H     -1.818030594599      2.301674035356     -0.938288696532
H     -3.118392729556     -0.677892159803      1.028561221526
H     -2.621726432423     -1.406875115537     -0.603561059788

初期経路を求めるための初期構造

2. 遷移状態構造最適化#

　run_autots.pyを適切に使用することで、自動的に遷移状態構造が得られる。以下にその手順を説明していく。

初期構造をMultiOptPy-1.20.4ディレクトリ内にbh9_9_6.xyzとして保存する。その後、同じディレクトリ内で、config_bh9_9_6.jsonを作成し、以下のように記述する。

config_bh9_9_6.json

{
  "work_dir": "bh9_9_6",
  "top_n_candidates": 3,
  "multioptpy_version": "1.20.4",
  
  "step1_settings": {
    "othersoft": "uma-s-1p1",
    "opt_method": ["rsirfo_block_fsb"],
    "use_model_hessian": "fischerd3",
    "spin_multiplicity": 1,
    "electronic_charge": -1, 
	"manual_AFIR": ["300", "3", "5"]
  },
  
  "step2_settings": {
    "othersoft": "uma-s-1p1",
    "NSTEP": 18,
    "use_model_hessian": "fischerd3",
    "save_pict": true,
	"QSMv2": true,
    "node_distance": 0.40,
    "align_distances_energy_predicted": 1,
    "spin_multiplicity": 1,
    "electronic_charge": -1
  },
  
  "step3_settings": {
    "othersoft": "uma-s-1p1",
    "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": -1
  },

  "step4_settings": {
    "othersoft": "uma-s-1p1",
	"opt_method": ["rsirfo_block_bofill"],
    "spin_multiplicity": 1,
    "electronic_charge": -1,
	"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_9_6.xyz -cfg config_bh9_9_6.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.20.4/OPTION_README.mdにて示されている。

Step 1#

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

以下のJSON内で記述したバイアスポテンシャルで、次の経路緩和アルゴリズムの初期経路として用いるトラジェクトリーを生成する。

"manual_AFIR": ["yyy", "a", "b]：yyykJ/molの活性化障壁を超えうるペア同士を近づける力を原子のラベル番号aとbのペアに構造最適化時に加えることを示す。
"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"]:力の定数c a.u.で、平行距離dÅの調和ポテンシャルをa番とb番の原子ペアにかける。初期経路生成時に開裂を防ぎたい結合を保持するとき等に使用する。

Step1が正常終了していれば作成されたwork_dirディレクトリ中に、bh9_9_6_step1_traj.xyzが存在する。必要に応じて確認し、目的に沿った初期経路が得られているか確認する。もし想定とは異なる場合は、プロセスをkillして再度設定を見直してやり直す。

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

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

Step 2#

Step2では、NEB法を用いることで、先ほど得られたbh9_9_6_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.4/"work_dir"と同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。

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

NEB計算の結果の可視化

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

経路緩和の結果、以下の構造がstep3の初期構造として自動的に用いられた。“work_dir”内のbh9_9_6_step3_TS_Opt_Inputs内に保存されたbh9_9_6_ts_guess_X.xyzにて確認が可能である。ts_guessの番号が小さい順にエネルギー値が高い構造を示すようになっている。

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

bh9_9_6_ts_guess_1.xyz

16
-1 1
C       2.171017606891     -0.749246417776     -0.622901840982
C       1.382311984330      0.209082433841      0.266959290660
C       0.186952697338     -0.255288044549      0.886013551752
O       1.794221830782      1.383262535091      0.317245656485
C      -1.449816054944      0.328707833182     -0.003519592538
C      -2.512563688177     -0.552191159298      0.084633438404
C      -1.119284860293      1.542563829792     -0.405729096166
H       3.228814207621     -0.679209284922     -0.358315269285
H       1.844731741161     -1.786264495207     -0.545503810913
H       2.076979819839     -0.417922718518     -1.658912396777
H      -0.123693596911      0.286437844899      1.768138194386
H       0.046750549106     -1.324678507966      0.938252309994
H      -0.116919343294      1.921685679702     -0.247247867914
H      -1.844870733974      2.180973721863     -0.894005514496
H      -3.040977452426     -0.658714037335      1.030484618032
H      -2.523654707049     -1.429199212797     -0.555591670643

NEB法により緩和した経路から得られた遷移状態構造を求めるための初期構造 (No.1)

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.4/"work_dir"ディレクトリ内に生成された新規ディレクトリ内のbh9_9_6_ts_final_X.xyzとして保存されている。)

bh9_9_6_ts_final_1.xyz

16
OptimizedStructure
C      2.009469879096     -0.704821747108     -0.722372302787
C      1.418003067867      0.168446587918      0.380767873820
C      0.339831555923     -0.347494205816      1.113850706422
O      1.830752053287      1.347160276130      0.456719773053
C     -1.420354885104      0.350470078080     -0.017917057822
C     -2.358219539095     -0.606208128185      0.052439252440
C     -1.099125150795      1.579964373384     -0.397358719913
H      3.098950965198     -0.645812538784     -0.688033600888
H      1.694035872118     -1.747070390219     -0.654170068300
H      1.684974310186     -0.299936468579     -1.685178938406
H      0.013375455484      0.193827433802      1.993361725685
H      0.132421450980     -1.408940813344      1.085030755196
H     -0.285114517777      2.129186713350      0.060129628772
H     -1.575462196312      1.986928679939     -1.282737324980
H     -3.361632150741     -0.339102306687      0.362690188341
H     -2.121906170315     -1.656597543880     -0.057221890633

遷移状態構造 (No.1)

24回の反復計算により、停留点に収束した構造が得られた。"frequency_analysis": trueオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。

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

Mode                                 0                   1                   2
Freq [cm^-1]                     -428.8111             59.3677             81.5807
Reduced mass [au]                   4.8465              3.3385              2.1046
Force const [Dyne/A]               -0.5251              0.0069              0.0083
Char temp [K]                       0.0000             85.4169            117.3764
Normal mode                   x         y         z            x         y         z            x         y         z     
       C                -0.00447   -0.00022   -0.00146    0.02026    0.07490   -0.04157    0.14768   -0.02697    0.05097
       C                -0.01306    0.01757   -0.02502    0.00931   -0.00249    0.01250    0.01263   -0.00617   -0.03623
       C                -0.14418    0.05923   -0.09030   -0.00411   -0.04931   -0.03957   -0.01812   -0.01677   -0.08495
       O                 0.00464   -0.02104    0.01734    0.00599   -0.00611    0.09292   -0.02007    0.00459   -0.03190
       C                 0.15634   -0.05725    0.09245   -0.00877   -0.01283   -0.01456   -0.03932    0.00765    0.01668
       C                 0.03087    0.00234    0.00125   -0.06035    0.04687    0.11048   -0.05957    0.02568    0.03899
       C                -0.01211   -0.00595    0.00520    0.03779   -0.05984   -0.12696   -0.04350    0.01467    0.03543
       H                -0.00649    0.00999   -0.00164    0.01826    0.12181   -0.06786    0.14855   -0.17607    0.31659
       H                -0.00445    0.00262    0.01268    0.06555    0.05825   -0.08516   -0.00077    0.01035   -0.07112
       H                -0.01429   -0.01542   -0.00304   -0.02561    0.10812   -0.01264    0.44162    0.06331   -0.00867
       H                 0.02814    0.00084    0.01082   -0.02074   -0.10204   -0.01377   -0.03692   -0.02947   -0.08436
       H                -0.02495    0.03097    0.00219    0.00254   -0.04917   -0.10647   -0.00145   -0.02007   -0.09182
       H                -0.03503    0.08707   -0.07328    0.07685   -0.05976   -0.19451    0.00063   -0.01730   -0.00418
       H                -0.14047   -0.05193    0.05528    0.01491   -0.10451   -0.13512   -0.11053    0.05945    0.09180
       H                -0.00983    0.11071   -0.22722   -0.02889    0.13292    0.13799   -0.02915    0.03512    0.12897
       H                -0.02570   -0.02807    0.16182   -0.12798    0.02345    0.18971   -0.09117    0.02459   -0.01948
       
(...snip...)

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

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

終わりに#

　　自作ライブラリで、UMAモデルのニューラルネットワークポテンシャル(NNP, uma-s-1p1)を用いて、BH9データセットの9. Nucleophilic addition, No. 6の反応のある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.