最終更新:2025-08-27
概要
本記事では、自作モジュール(MultiOptPy)で、アルキンへのヒドロホウ素化 (ホウ素化試薬
MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
ヒドロホウ素化反応について:
- https://www.chem-station.com/odos/2009/06/brown-brown-hydroboration.html
- https://en.wikipedia.org/wiki/Hydroboration
今回使用したニューラルネットワークポテンシャルについて:
- 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.15.1b
環境
Windows 11
※Windows 11環境下でAnaconda PowerShell Promptを使用した。
Source codeのダウンロード
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.15.1b.zip
unzip v1.15.1b.zip
cd MultiOptPy-v1.15.1b
環境構築手順
今回は、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.15.1b
ディレクトリ内に存在するsoftware_path.conf
に対して、uma-s-1.pt
の絶対パスを用いて以下を追記する。
uma-s-1::(uma-s-1.ptの絶対パス)
これで、MultiOptPy-v1.15.1b
がNNPuma-s-1
を使用できるようになる。
使用するNNPに関する具体的な説明
今回使用するNNPについて具体的に説明する。
- UMAのModel Checkpointは
uma-s-1
を使用した。 - 小分子系のトレーニングセットである
Omol25
(omol
)を使用して学習したニューラルネットワークポテンシャルを使用する。
※自作モジュールでの具体的な使用の仕方に関しては、ase_calculation_tools.py
を参照
手順
1. 初期構造の準備
モデル反応系として、以下の構造を用意した。今回はファイルの名前をalkyne_hydroboration_BHCy2.xyz
とした。 初期構造は以下のものを使用した。
43
C -1.329604699625 -2.198441329275 -2.370961262034
C -0.386178954888 -1.790364592569 -2.987649762238
H 0.446180873789 -1.430017185906 -3.541757557139
C -2.480628088959 -2.693213172454 -1.624185226351
H -2.301536030721 -2.606920712210 -0.552601301048
H -3.371839791651 -2.115769200542 -1.869928145255
H -2.671350761149 -3.739340341570 -1.862583140198
B -0.120271935272 -0.269561861578 0.376777972333
H -0.744063179813 -1.274147403539 0.592513093818
C -0.888280790070 1.037993480494 0.005460987589
C -2.225031216239 0.874571191471 -0.721401831121
C -1.121351049360 1.719765492914 1.379348484079
H -0.241946621401 1.713769028257 -0.573234415765
C -2.963637387237 2.201871153380 -0.871552969295
H -2.850312953393 0.177234982599 -0.149382484044
H -2.062802572537 0.417744890560 -1.700087278117
C -1.859538776017 3.048671655706 1.225191832107
H -1.718172662916 1.054324999130 2.015997902340
H -0.169869140083 1.878740427963 1.894423283009
C -3.180706806620 2.868220358677 0.483890058333
H -3.920601689491 2.050012668499 -1.377032256846
H -2.373806686348 2.869726726690 -1.509861559247
H -2.031856934170 3.496317622448 2.207213145512
H -1.222975714538 3.743183609271 0.666348518920
H -3.679867552729 3.832292037225 0.361596159166
H -3.847531556542 2.242522287958 1.088334068358
C 1.429968371513 -0.282391265641 0.596397620887
C 1.905388465888 -1.205205161845 1.724259116759
C 2.089057311590 -0.707264652577 -0.735412451558
H 1.781120736403 0.739207062868 0.806227624890
C 3.426988149939 -1.231931837689 1.839486558025
H 1.537717806975 -2.218750274799 1.525311043015
H 1.461673318964 -0.897639457644 2.675697909595
C 3.611964808561 -0.728051267309 -0.620469551491
H 1.732594454937 -1.706190165230 -1.010714912322
H 1.781111772151 -0.033893675902 -1.540842477221
C 4.069672062029 -1.637955133272 0.515989931531
H 3.736863702364 -1.912736924811 2.636090636567
H 3.783270080883 -0.234359337713 2.120950388053
H 4.055324506645 -1.050747917133 -1.566109819048
H 3.969148854532 0.290775351056 -0.433866689403
H 5.158845120839 -1.621286377353 0.599550433771
H 3.786873153764 -2.670765778606 0.282578321082
2. 遷移状態構造最適化のための初期構造の算出
a. 初期構造をカレントディレクトリにalkyne_hydroboration_BHCy2.xyz
として保存し、以下を実行する。
python .\optmain.py .\alkyne_hydroboration_BHCy2.xyz -opt rsirfo_fsb -os uma-s-1 -modelhess -ma 200 1 9 200 2 8
-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ライブラリが必要である。
これを実行すると、omol
のデータセットを使用したuma-s-1
モデルのNNPで得たエネルギーに対して、指定した人工力ポテンシャルを加えた上で初期構造を構造最適化することができる。
結果はyyyy_mm_dd
(今日の年月日)ディレクトリの中に存在するディレクトリを開いて確認できる。
正常終了していれば、このディレクトリ中に、alkyne_hydroboration_BHCy2_traj.xyz
が存在するので、これをコピーして、MultiOptPy-v1.15.1b
ディレクトリに置く。
alkyne_hydroboration_BHCy2_traj.xyz
は構造最適化の過程をAvogadro(公式ページ:https://avogadro.cc/ )等で可視化して確認できるようにしている。このalkyne_hydroboration_BHCy2_traj.xyz
は次のNEB計算に使用する。
※alkyne_hydroboration_BHCy2_traj.xyz
をアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
(この人工力ポテンシャルを加えて行った構造最適化の結果はalkyne_hydroboration_BHCy2_optimized.xyz
で確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-ma
の設定を見直してb.をやり直す。今回の場合以下の構造が得られた。生成系の安定構造に近いものが得られていると判断できるため、次のNEB計算を行うことが可能である。)
b. alkyne_hydroboration_BHCy2_traj.xyz
を初期パスとして、NEB法で経路の緩和を行う。
NEB法を用いることで、先ほど得られたalkyne_hydroboration_BHCy2_traj.xyz
全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。(この時点ではまだ正確な遷移状態構造は求められていない。)
python .\nebmain.py alkyne_hydroboration_BHCy2_traj.xyz -os uma-s-1 -modelhess -ns 10 -spng -nd 0.3
-nd N
はノード間の距離をN Åとして初期パスを作成することを示す。-ns n
はn回分NEB法による経路の緩和を行うことを示す。-fc M
はM回あたりの経路緩和回数に対して1回だけ正確なHessianを計算し、経路緩和に使用する。これを使用すると、Hessianを使用しない場合の経路緩和アルゴリズムとは別のものを使用して、経路緩和を行う。-spng
は緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。-modelhess
は経路の緩和に計算コストの低いHessianを用いることを指定する。これに伴い、Hessianを使う経路緩和アルゴリズムに変更となる。-os uma-s-1
は今回使用するNNPを指定している。これを使用する際にASEライブラリが必要である。
c. 初期構造の決定及び遷移状態構造の計算
MultiOptPy-v1.15.1b
と同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。 そのディレクトリ内のenergy_plot.csv
を確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
パスの緩和後の各ノードのエネルギー一覧(単位
※bias_force_rms.csv
にて、各Iterationごとのすべてのノードの勾配のRMS値を確認できる。
経路緩和の結果得られた経路のエネルギー極大値を示す構造を遷移状態構造を求める初期構造とした。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
alkyne_hydroboration_BHCy2_traj_10.xyz
(グラフでは11番)
43
0 1
C -1.130529243315 -1.819070140764 -1.432651257611
C -0.243412591750 -1.074783966579 -1.890185711971
H 0.283172254928 -0.677490584153 -2.729339657130
C -2.386193495551 -2.570076928964 -1.212310208778
H -2.555852561633 -2.681936062836 -0.148268253905
H -3.219363930129 -2.006981022266 -1.647348634124
H -2.368490080867 -3.550157644606 -1.689669962177
B -0.197233366581 -0.505735676521 -0.292469822627
H -0.925564392740 -1.418047626526 0.217697488353
C -0.967984937997 0.918786107016 -0.214116191049
C -2.352721213658 0.884460561951 -0.855457435817
C -1.093315876264 1.530429088441 1.194381053637
H -0.339260688760 1.620930153940 -0.787952736047
C -3.031636483219 2.245261984806 -0.906631747324
H -2.975125203617 0.195196844840 -0.265347003029
H -2.286633154566 0.467612727218 -1.861888690941
C -1.784891925789 2.901900222115 1.166400203083
H -1.678949995641 0.865660828510 1.842964279262
H -0.117283264329 1.642149623802 1.678537926048
C -3.152884451952 2.840956098206 0.489739602678
H -4.012052438754 2.159417175152 -1.381763456965
H -2.429660886706 2.914301144506 -1.532242479599
H -1.892454184393 3.297126904420 2.180792633063
H -1.156696455887 3.618853274749 0.623896978072
H -3.607254644778 3.834158928851 0.442342622193
H -3.825665235914 2.216291032513 1.089021962973
C 1.357641959707 -0.482929203485 0.276167286462
C 1.772028237915 -1.292619373273 1.516135072749
C 2.290726802992 -0.781140656491 -0.901857793756
H 1.514849686385 0.573855515680 0.547200530844
C 3.285810407119 -1.235661606191 1.787384712680
H 1.481622912318 -2.342536277548 1.395619911732
H 1.248553021783 -0.928601012815 2.407523896679
C 3.773743728937 -0.728943200484 -0.604988858226
H 2.041296018357 -1.786666953300 -1.264683104319
H 2.058427392428 -0.083466770897 -1.709563954094
C 4.101215548065 -1.643749391702 0.563432275346
H 3.539378455391 -1.878289960022 2.634914406075
H 3.578178896005 -0.218468136438 2.076779206991
H 4.322197641337 -1.018234654110 -1.505194314217
H 4.064305908802 0.299332476540 -0.362391115328
H 5.172249538851 -1.626540291514 0.775938714398
H 3.845712293472 -2.674553551770 0.289451625716
構造が壊れていないので、これを遷移状態を求めるための初期構造とする。
alkyne_hydroboration_BHCy2_traj_10.xyz
をMultiOptPy-v1.15.1b
と同じディレクトリ内にコピーする。
そして、以下を実行する。
python .\optmain.py .\alkyne_hydroboration_BHCy2_traj_10.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と同等)
実行して得られた正確な遷移状態構造と思われる構造を以下に示す。
(実行して得られた正確な遷移状態構造は計算開始時に、yyyy_mm_dd
ディレクトリ内に生成された新規ディレクトリ内のYYY_optimized.xyz
として保存されている。)
43
OptimizedStructure
C -1.195295510927 -2.035449631874 -0.897293958562
C -0.426756034019 -1.429759607064 -1.659739784565
H 0.071669006214 -1.316521638950 -2.600664188903
C -2.175146433759 -2.993282040106 -0.369061915183
H -1.776583889807 -3.513029947550 0.500746747784
H -3.092050692176 -2.480543970082 -0.082402296651
H -2.403206961394 -3.722262954486 -1.147324924502
B -0.318803501455 -0.499841118953 -0.212256615098
H -1.078338967360 -1.125296312511 0.569016426285
C -0.990907441765 0.950968940421 -0.451749323699
C -2.479743882375 0.885268351639 -0.801655074796
C -0.790804756599 1.849880068828 0.774752716606
H -0.477888737006 1.439469689311 -1.295268579926
C -3.101170124092 2.266661588072 -0.994254374184
H -3.007902101346 0.371837004158 0.014389720756
H -2.634746196017 0.280198285205 -1.699264486001
C -1.407399862781 3.234755451119 0.588266414080
H -1.249740091567 1.364933085406 1.646893387217
H 0.273371769730 1.949544651950 1.003912740034
C -2.889364079159 3.140392536536 0.238530575942
H -4.167476418898 2.177678434934 -1.219091091551
H -2.635365457953 2.750589679075 -1.860446415556
H -1.266657784776 3.838715879666 1.488576351987
H -0.881639872197 3.751932544956 -0.222815858216
H -3.309956323606 4.136294059278 0.078886928244
H -3.430133871980 2.701421355963 1.085191672214
C 1.171923416879 -0.631725603843 0.396323471501
C 1.609940860912 -2.068293589435 0.684777152893
C 2.218081170966 0.066625204243 -0.478377125720
H 1.158332290382 -0.105142438459 1.361800911587
C 3.000890540035 -2.131771384979 1.311939006744
H 1.620448338939 -2.640393438690 -0.251542887137
H 0.879178392993 -2.557822971279 1.337562762606
C 3.614426184887 0.010681915437 0.138213736418
H 2.246971603091 -0.422729395214 -1.461839362208
H 1.928144946292 1.105895008708 -0.660212842704
C 4.029976964238 -1.428986073607 0.429555117004
H 3.297214717756 -3.169856162389 1.485211440535
H 2.975303796567 -1.642195733365 2.292045315638
H 4.344722911800 0.490887173924 -0.518469008608
H 3.611508512394 0.578707394841 1.075539971616
H 5.016204810044 -1.454343036287 0.899359292364
H 4.118768758895 -1.974091254546 -0.517761746290
停留点に収束した分子構造が得られた。-freq
オプションにより生成されたnormal_modes.txt
やvibration_animation
ディレクトリ内の振動モードのアニメーションを確認した。
以下に-freq
オプションで生成されたnormal_modes.txt
の一部を示す。
Mode 0 1 2
Freq [cm^-1] -439.5525 36.6339 45.0449
Reduced mass [au] 2.3568 2.9757 3.0368
Force const [Dyne/A] -0.2683 0.0024 0.0036
Char temp [K] 0.0000 52.7079 64.8096
Normal mode x y z x y z x y z
C 0.00102 -0.06803 -0.12338 0.02265 0.00853 -0.04770 0.02889 0.02356 0.00245
C -0.01238 -0.02420 -0.01995 0.02350 0.02150 -0.03708 0.01034 0.03703 -0.00522
H 0.18668 0.35117 0.13724 0.02859 0.03436 -0.03273 -0.00021 0.04607 -0.00964
C 0.00984 0.02262 0.01671 0.02611 -0.00673 -0.06862 0.06035 -0.00041 0.01689
H 0.04130 0.06230 0.02375 0.02321 -0.01501 -0.07221 0.08710 0.01020 0.01103
H 0.03866 0.06560 0.02589 0.01959 -0.01787 -0.06947 0.05150 -0.02400 0.03060
H -0.04006 -0.01377 0.06369 0.03864 0.00174 -0.08022 0.06749 -0.00668 0.02067
B -0.06552 -0.00273 0.16679 0.00422 0.00757 -0.02841 0.01533 0.03434 -0.00294
H 0.04279 0.24227 0.33182 0.00462 -0.00920 -0.04065 0.01564 0.03450 -0.00263
C 0.00645 -0.00157 -0.01370 -0.01466 0.00112 -0.02117 0.00672 0.03199 0.00315
C 0.00290 0.00240 -0.00340 -0.03746 -0.02289 0.07811 0.01913 0.01847 -0.04782
C 0.00385 0.00312 -0.01008 0.05548 0.03465 -0.05729 -0.04270 0.00610 0.03035
H 0.01768 -0.01665 -0.01506 -0.07351 -0.01423 -0.06588 0.02867 0.05505 0.03002
C 0.00066 0.00124 -0.00182 -0.06667 -0.03542 0.08317 -0.00103 0.01277 -0.02246
H 0.01198 0.00395 0.00179 0.02321 -0.00624 0.12756 0.00076 -0.01842 -0.08273
H -0.00401 0.00080 -0.00210 -0.09004 -0.04945 0.10477 0.05940 0.04544 -0.07265
C 0.00028 0.00191 -0.00327 0.02851 0.02306 -0.05481 -0.06979 -0.00226 0.05682
H 0.00306 -0.00548 -0.01664 0.11791 0.05345 -0.01390 -0.05031 -0.02790 0.00747
H 0.00311 0.00551 -0.00743 0.06924 0.05287 -0.12958 -0.04959 0.02396 0.05537
C 0.00003 -0.00018 -0.00061 0.00618 -0.00146 0.04667 -0.05805 -0.02483 0.01360
H 0.00053 -0.00030 -0.00091 -0.08084 -0.05417 0.15776 0.00817 -0.00113 -0.06059
H -0.00179 0.00471 -0.00146 -0.13124 -0.05197 0.03922 0.01785 0.05112 0.00911
H 0.00039 -0.00290 0.00014 0.08258 0.04960 -0.08108 -0.10736 -0.02790 0.07990
H -0.00270 0.00818 -0.00106 -0.03078 0.00565 -0.10436 -0.06038 0.03593 0.08731
H -0.00152 -0.00045 0.00111 -0.01499 -0.01000 0.04919 -0.07702 -0.02922 0.03626
H 0.00201 -0.00254 -0.00067 0.06747 0.01457 0.09410 -0.07143 -0.06719 -0.01687
C 0.00951 0.00250 -0.01747 -0.00202 0.01129 -0.01428 0.01735 0.02798 -0.00895
C 0.00105 0.00228 -0.00311 -0.02237 0.01613 0.04070 -0.02591 0.02873 0.06331
C 0.00515 -0.00015 -0.01027 0.01511 -0.03239 -0.02827 0.03021 -0.04755 -0.05368
H 0.01844 0.00436 -0.01607 -0.00512 0.04435 -0.03240 0.03994 0.07715 -0.03560
C -0.00001 0.00007 -0.00118 -0.03028 0.02336 0.05914 -0.01754 0.01354 0.04277
H 0.00059 -0.00600 0.00120 -0.01924 -0.01729 0.06103 -0.06204 -0.02420 0.09514
H -0.00478 0.01529 -0.00172 -0.03463 0.04865 0.05113 -0.03147 0.09188 0.10417
C 0.00225 0.00003 -0.00319 0.00703 -0.02375 -0.00916 0.03956 -0.07034 -0.07662
H 0.01980 0.00460 -0.01377 0.02258 -0.07134 -0.00852 -0.00228 -0.09058 -0.03312
H 0.00279 -0.00045 -0.00812 0.02708 -0.03654 -0.07262 0.06296 -0.04490 -0.09293
C 0.00105 -0.00002 -0.00026 -0.01183 -0.01784 0.04787 -0.00884 -0.07174 -0.01442
H 0.00023 0.00006 -0.00078 -0.04295 0.02640 0.09900 -0.05203 0.01239 0.09486
H -0.00133 0.00058 -0.00150 -0.03758 0.05760 0.04185 0.01756 0.06865 0.01618
H 0.00718 0.00040 0.00256 0.02024 -0.05557 -0.01772 0.04482 -0.12943 -0.11398
H -0.00519 -0.00041 -0.00293 0.00156 0.01037 -0.02987 0.07939 -0.02521 -0.10390
H 0.00082 0.00055 0.00042 -0.01830 -0.01115 0.06182 0.00027 -0.08584 -0.03432
H 0.00232 -0.00094 0.00017 -0.00470 -0.05228 0.06836 -0.05033 -0.12291 0.01110
(...snip...)
その結果、虚振動が1つであることが確認できた。つまりこの構造は遷移状態構造である。
次に、vibration_animation
内の虚振動を示す分子振動が示されたxyzファイル(mode_1_XXXi_wave_number.xyz
)をAvogadroで確認すると、想定される反応系と生成系をつなぐ方向に振動していることを確認できた。
終わりに
自作モジュールで、UMAモデルのニューラルネットワークポテンシャル(NNP)を用いて、アルキンへのヒドロホウ素化 (ホウ素化試薬
参考
- 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へのアカウント登録と配布元の使用許諾が必要である。)
- 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.