# 導入 気象庁に先んじて[全国的な梅雨明けが宣言された](https://www.youtube.com/watch?v=Lz24PqZkF2s)ので、少し気が早いが夏休みの自由研究に取り組んでみた。本稿では機械学習ポテンシャル「OrbMol-v2」の特性について扱う。 ## はじめに 近年、分子や材料のポテンシャルエネルギー面（PES; potential energy surface）を高速に予測する汎用機械学習原子間ポテンシャル（uMLIP; universal machine learning interatomic potential）が急速に発展している。従来の解析では第一原理計算が必要で多大な時間を要していた構造最適化や分子動力学シミュレーションを、より低コストに実行できる可能性が拡大している。 uMLIP の多くはニューラルネットワークモデルに基づいて構築されている。こうしたニューラルネットワークポテンシャル（NNP）においては、共有結合性の分子や局所的な構造変化であれば、近傍原子の情報のみから良い精度でPESの形状を予測できる。一方で、uMLIP はその性質上、長距離静電相互作用の記述を苦手とすることが知られている。 電解質などのイオン性化合物、イオン対の双極子相互作用、解離極限などでは、遠く離れた電荷同士のクーロン相互作用が重要になるが、局所的な環境のみを考慮するNNP（商用のものを含めて多くの局所的 message-passing 型 NNP が該当する）ではこうした長距離相互作用が記述できず、PESの形状が不正確になるだけでなく、解離極限や相互作用エネルギーの size-consistency/size-extensivity 的な振る舞いにも問題が生じる。これは既知の問題だが、汎用NNPにおいて、精度・速度・汎用性を同時に満足することは容易ではない。 > 局所近似の限界についてはPFN高本氏の記事が参考になる。 https://tech.preferred.jp/ja/blog/%E6%A9%9F%E6%A2%B0%E5%AD%A6%E7%BF%92%E3%83%9D%E3%83%86%E3%83%B3%E3%82%B7%E3%83%A3%E3%83%AB%E3%81%AE%E8%BF%91%E4%BC%BC%E3%81%AE%E9%99%90%E7%95%8C%E3%82%92%E8%B6%85%E3%81%88%E3%81%A6%E3%80%90%E3%82%A4/ ## OrbMol-v2 の "latent charge" 今年5月、Orbital Materials 社は uMLIP の一つである Orb シリーズの最新版 "OrbMol-v2" をリリースした。この OrbMol-v2 は上述の問題に対して興味深い仕組みを持つ。モデル内部で各原子に対して潜在的な電荷（文字通り "latent charge" と名付けられている）をモデル内で予測し、この値を長距離静電相互作用の計算に用いるのである。 この latent charge は Mulliken 電荷、Hirshfeld 電荷、Bader 電荷、RESP 電荷のような既存の部分電荷とは異なる。OrbMol-v2 は、あらかじめ定義された原子電荷を教師データとして学習しているわけではなく、エネルギーや力を再現する過程でモデル内部に現れる量として latent charge を用いている。 latent charge は、その総和が系の総電荷に一致するという制約の下で各原子に割り当てられる量であり、クーロン相互作用の計算に直接用いられる。この量は既知の何らかの電荷と対応付けることに意味があるとは言えないが、電荷として振る舞わざるを得ない制約条件をモデル内に与えた結果として部分電荷に似た性質を有する量が創発した、と見なすのが妥当と思われる。 https://www.orbitalindustries.com/news/blog/introducing-orbmol-v2 > OrbMol のドキュメントでは latent charge について「電荷ラベルなしで、エネルギー・力だけから創発した潜在的な特徴量」と述べられており、一般的な部分電荷とは異なる扱いをする必要があるとされている。これは何らかの物理的背景を基に導出された量ではないことを述べたものだが、解釈可能な量である可能性は否定していない。というわけで、この latent charge は学術的に大変興味深い量だと言える。 :::note info ### Orb系のモデルの名称について OrbMol-v2とOrb-v3は、バージョン番号が異なるため両者の関係が非常に分かりにくいが、二つの数字は別の系列を表している。 Orb-v3は、Orbital Materialsが開発するニューラルネットワーク型原子間ポテンシャル「Orb」の第3世代にあたる。単一の学習済みモデルだけを指すのではなく、共通するネットワーク設計や学習方法を採用した複数のモデルからなるモデルファミリーである。 OrbMolは、このOrb-v3系の設計を基盤として、分子系を扱うために開発されたニューラルネットワークポテンシャルである。OrbMol-v2は、その分子向けモデルであるOrbMolの第2世代にあたる。 したがって、OrbMol-v2の「v2」はOrbMol系列のバージョンを、Orb-v3の「v3」は基盤となるOrb系列の技術世代を表している。つまり、OrbMol-v2は「Orb-v3より一世代古いモデル」という意味ではなく、「第3世代のOrb技術を基盤とした、第2世代の分子向けOrbMolモデル」という理解が適切である。 ```:Orbファミリーの系統図 Orb └─ 基本ブランド／NNPファミリー ├─ Orb-v2 ├─ Orb-v3 │ ├─ orb-v3-conservative-inf-omat │ ├─ orb-v3-direct-inf-omat │ └─ OrbMol系の基盤 │ └─ OrbMol ├─ OrbMol └─ OrbMol-v2 ``` ::: :::note info ### モデルのアーキテクチャについて 余談であるが、モデルのアーキテクチャについて少し触れておく。回転対称性のある系では、原子構造全体を回転行列 $R$ を作用させたとき、エネルギーと力は$$E(Rr)=E(r)$$ $$F(Rr)=RF(r)$$をそれぞれ満足するべきである。 [UMA](https://arxiv.org/abs/2506.23971)は回転等変性をアーキテクチャに明示的に組み込んでいるため、上記の要請を自然に満たす形でモデルが構築される。UMAの基盤である[EquiformerV2](https://arxiv.org/abs/2306.12059)は、SO(3)畳み込みを効率化したeSCNを用い、高次の等変特徴量を扱う設計となっている。このアーキテクチャは少量データでの学習効率、ファインチューニング時の安定性、未知構造に対する対称性の保持において有利である。 一方、Orbは通常のGraph Network Simulator型GNNを基盤とし、相対位置ベクトルを一般的なニューラルネットワークに入力する。 これはSO(3)等変演算による計算量とメモリ消費の増大を避け、通常のGNN演算、大規模学習データ、回転データ拡張、equigrad正則化によって回転対称性を近似的に学習させることで、推論速度と大規模系へのスケーラビリティを優先する設計思想である。 したがって、Orb-v3が明示的な回転等変ニューラルネットワークを採用していないのは、単に実装が未成熟だからという理由ではない。 （参考：https://arxiv.org/html/2410.22570v1） したがって、両者のアーキテクチャの違いは技術的な新旧というより、物理的要請と計算効率のどちらを重視するかという設計上のトレードオフと理解するのが適切である。 > NNPの設計についてはPFN岡野原氏の記事が分かりやすい。 > https://hillbig.github.io/ai_compchem/ai_compchem_2.pdf ::: # 環境構築 モデルの重みファイルはHugging Faceから入手できる。 https://huggingface.co/Orbital-Materials/orbmol-v2/blob/main/README.md > オンライン環境の場合はスクリプトからモデルを呼ぶと、重みファイルが自動的にダウンロードされてキャッシュに保存されるため、わざわざ手動でダウンロードする必要はない。 ### 仮想環境を作成 ```sh conda create -n orb312 python=3.12 -y ``` ### 仮想環境を起動＆pip等を更新 必須ではないが念のため更新しておく。 ```sh python -m pip install --upgrade pip setuptools wheel ``` ### GPU版PyTorchを導入 ```sh python -m pip install torch --index-url https://download.pytorch.org/whl/cu130 ``` 上手くインストールできない場合はCUDAのバージョンが古い可能性があるので確認する。Blackwell GPU（RTX 50シリーズなど）の場合はorb-modelsパッケージはCUDA 13.xが必須なので注意。 GPUがBlackwell世代だと古いPyTorchやCUDAのバージョンが対応していないので、 `torch 2.12.1+cu130` などをインストールして対応すれば良い。 :::note warn #### 筆者がスタックした点 WSL環境でCUDAを利用する場合、PythonのモジュールがCUDA driver libraryを正しく読めないことがある。筆者の場合、WSL 側の /usr/lib/wsl/lib/libcuda.so を見つけられていないことが原因でGPUが認識されなかった（`ldconfig -p | grep -E 'libcuda|nvidia' || true` などのコマンドで確認）。 そこで、以下のようなコマンドによって `LD_LIBRARY_PATH` 内で先に `/usr/lib/wsl/lib` が読み込まれるように環境変数を更新して対応した。 ```sh export LD_LIBRARY_PATH=/usr/lib/wsl/lib${LD_LIBRARY_PATH:+:$LD_LIBRARY_PATH} ``` これにより `warp-lang` が正しい `libcuda.so` を見つけられるようになり、GPUを認識できないエラーが解消した。 ::: ### orb-modelsパッケージを導入 rdkitは必須ではない。必要なら入れておく。 ```sh python -m pip install orb-models pandas numpy matplotlib ase scipy tqdm rdkit ``` ### CUDAの実行テスト ```sh python - <<'PY' import torch print("torch:", torch.__version__) print("torch.version.cuda:", torch.version.cuda) print("cuda available:", torch.cuda.is_available()) print("arch list:", torch.cuda.get_arch_list()) x = torch.randn(64, 64, device="cuda") y = x @ x torch.cuda.synchronize() print("matmul ok:", y.shape) PY ``` ```sh:出力結果 torch: 2.12.1+cu130 torch.version.cuda: 13.0 cuda available: True arch list: ['sm_75', 'sm_80', 'sm_86', 'sm_90', 'sm_100', 'sm_120'] matmul ok: torch.Size([64, 64]) ``` > このように出力されれば、問題なく導入できている。 # OrbMolで計算する ## テスト用スクリプト ChatGPT（GPT-5.6 Sol）にテスト計算用のスクリプトとデータセットを作成させた。ファイル構成のイメージは以下の通り。 ``` orbmol_latent_charge/ ├── README.md ├── requirements.txt ├── systems.csv ├── scripts/ │ ├── 00_check_env.py │ ├── 01_run_orbmol_latent_charges.py │ ├── 02_plot_charge_summaries.py │ ├── 03_make_alkali_halide_scans.py │ └── 04_xtb_reference_charges_optional.py └── structures/ ├── sanity/ ├── polar/ ├── ions/ ├── scans/ ├── electrolytes/ └── metal_complexes ``` スクリプトや構造データはGitHub上で公開したので、以下から取得できる。 https://github.com/h-nabata/MLIP/tree/main/test/orbmol/latent_charge_test ここでは、以下の54個の化学種を対象に、テスト計算を実施する。[^Zr_crd] ``` 1. sanity: H2, CH4, C6H6, CO2, N2 2. polar: HF, HCl, H2O, NH3, CH3Cl, CH3OH, CH3CN 3. ions: NH4+, H3O+, OH-, acetate-, nitrate-, carbonate2- 4. scans: LiF, NaCl, KBr の原子間距離でスキャン 5. electrolytes: Li+(H2O)n, Cl-(H2O)n, Li+(DME) 6. metal_complexes: ferrocene, ferrocenium, TiCl4, ZrCl4, Cp2ZrCl2, Cp2ZrMe+, Cp2ZrMe(ethylene)+ ``` 上から順に、1. sanityチェックを兼ねた中性分子、2. 分子内で分極した極性分子、3. イオン、4. 二原子ポテンシャル曲線を取得するための塩、5. イオン性の電解質、6. 金属錯体、である。 ここで、各スクリプトの基本的な出力ファイルは以下の通り。 ``` results/systems_summary.csv # energy, force RMS, latent charge 総和誤差 results/atom_latent_charges.csv # 原子ごとの latent charge results/group_charges.csv # Cp, Me, ethylene, solvent などの group charge results/run_metadata.json # Python/Torch/ASE/orb-models/device/precision ``` > 内部的には、以下のように `pretrained.orbmol_v2()` によって「推論用に重みを読み込んだ力場モデル本体」および「ASE の Atoms オブジェクトを OrbMol が読めるグラフ形式へ変換するためのアダプタ」を作成している。 > ```py:モデルの呼び出し > orbff, atoms_adapter = pretrained.orbmol_v2( > device="cuda", > precision="float32-high", > compile=False, > ) > ``` > > ```:OrbMol-v2 の学習済みモデル本体（orbff）の構造 > orbff > ├── model > │ └── MoleculeGNS backbone > ├── heads > │ ├── energy > │ │ └── ChargeConditionedEnergyHead > │ ├── confidence > │ │ └── ConfidenceHead > │ └── latent_charges > │ └── LatentChargeHead > ├── coulomb_module > │ └── CoulombModule > └── pair_repulsion_fn > └── ZBL 型の短距離反発項 > ``` ## latent charge の確認 早速、`1. sanity` のセットから順にOrbMolで計算してみる。 ### 無極性分子 ```sh:実行例 python scripts/01_run_orbmol_latent_charges.py \ --systems systems.csv \ --outdir results_sanity \ --device cuda \ --limit 5 ``` > スクリプトの内容は参考資料に掲載しているので、興味のある方は参照されたい。また、以下のように `device` を `cpu` にすれば orbmol_v2 の呼び出し時にCPUが指定されるようにしている。 > ``` > python scripts/01_run_orbmol_latent_charges.py \ > --systems systems.csv \ > --outdir results_smoke_cpu \ > --device cpu \ > --limit 5 > ``` 筆者の環境での実行時間は、モデルの読み込みを含めて約10秒程度であった。結果を集約した `atom_latent_charges.csv` の抜粋を以下に示す。
原子ごとの結果
|system_id|energy_eV|element|x|y|z|latent_charge| |:----|----:|:----|----:|----:|----:|----:| |H2|-31.588|H|-0.3700|0.0000|0.0000|0.0008| |H2|-31.588|H|0.3700|0.0000|0.0000|-0.0008| |CH4|-1101.818|C|0.0000|0.0000|0.0000|-0.1193| |CH4|-1101.818|H|0.6293|0.6293|0.6293|0.0298| |CH4|-1101.818|H|0.6293|-0.6293|-0.6293|0.0301| |CH4|-1101.818|H|-0.6293|0.6293|-0.6293|0.0298| |CH4|-1101.818|H|-0.6293|-0.6293|0.6293|0.0297| |C6H6|-6319.211|C|1.3970|0.0000|0.0000|-0.0782| |C6H6|-6319.211|C|0.6985|1.2098|0.0000|-0.0782| |C6H6|-6319.211|C|-0.6985|1.2098|0.0000|-0.0778| |C6H6|-6319.211|C|-1.3970|0.0000|0.0000|-0.0782| |C6H6|-6319.211|C|-0.6985|-1.2098|0.0000|-0.0783| |C6H6|-6319.211|C|0.6985|-1.2098|0.0000|-0.0785| |C6H6|-6319.211|H|2.4790|0.0000|0.0000|0.0783| |C6H6|-6319.211|H|1.2395|2.1469|0.0000|0.0783| |C6H6|-6319.211|H|-1.2395|2.1469|0.0000|0.0781| |C6H6|-6319.211|H|-2.4790|0.0000|0.0000|0.0781| |C6H6|-6319.211|H|-1.2395|-2.1469|0.0000|0.0781| |C6H6|-6319.211|H|1.2395|-2.1469|0.0000|0.0783| |CO2|-5132.730|O|-1.1600|0.0000|0.0000|-0.1533| |CO2|-5132.730|C|0.0000|0.0000|0.0000|0.3067| |CO2|-5132.730|O|1.1600|0.0000|0.0000|-0.1534| |N2|-2981.070|N|-0.5500|0.0000|0.0000|0.0014| |N2|-2981.070|N|0.5500|0.0000|0.0000|-0.0014|
以上の結果は、OrbMol-v2 の latent charge が少なくとも単純な中性分子に対して、総電荷制約・対称性・局所極性をおおむね自然に反映していることを示している。 ### 極性分子
原子ごとの結果
|system_id|energy_eV|element|x|y|z|latent_charge| |:----|----:|:----|----:|----:|----:|----:| |HF|-2733.711|H|0.0000|0.0000|0.0000|0.2641| |HF|-2733.711|F|0.9200|0.0000|0.0000|-0.2641| |HCl|-12538.756|H|0.0000|0.0000|0.0000|0.1501| |HCl|-12538.756|Cl|1.2700|0.0000|0.0000|-0.1501| |H₂O|-2079.864|O|0.0000|0.0000|0.0000|-0.4180| |H₂O|-2079.864|H|0.7575|0.0000|0.5865|0.2089| |H₂O|-2079.864|H|-0.7575|0.0000|0.5865|0.2091| |NH₃|-1538.811|N|0.0000|0.0000|0.0000|-0.4071| |NH₃|-1538.811|H|0.9365|0.0000|0.3784|0.1356| |NH₃|-1538.811|H|-0.4682|0.8110|0.3784|0.1355| |NH₃|-1538.811|H|-0.4682|-0.8110|0.3784|0.1360| |CH₃Cl|-13608.026|Cl|0.0000|0.0000|0.0000|-0.1761| |CH₃Cl|-13608.026|C|1.7800|0.0000|0.0000|-0.0290| |CH₃Cl|-13608.026|H|2.1433|-1.0277|0.0000|0.0685| |CH₃Cl|-13608.026|H|2.1433|0.5138|-0.8900|0.0683| |CH₃Cl|-13608.026|H|2.1433|0.5138|0.8900|0.0683| |CH₃OH|-3148.560|C|0.0000|0.0000|0.0000|0.0051| |CH₃OH|-3148.560|H|-0.3633|1.0277|0.0000|0.0666| |CH₃OH|-3148.560|H|-0.3633|-0.5138|-0.8900|0.0387| |CH₃OH|-3148.560|H|-0.3633|-0.5138|0.8900|0.0389| |CH₃OH|-3148.560|O|1.4300|0.0000|0.0000|-0.3698| |CH₃OH|-3148.560|H|1.9500|0.7500|0.0000|0.2205| |CH₃CN|-3612.348|C|-1.4500|0.0000|0.0000|-0.0380| |CH₃CN|-3612.348|H|-1.8133|1.0277|0.0000|0.0833| |CH₃CN|-3612.348|H|-1.8133|-0.5138|-0.8900|0.0836| |CH₃CN|-3612.348|H|-1.8133|-0.5138|0.8900|0.0833| |CH₃CN|-3612.348|C|0.0000|0.0000|0.0000|0.1178| |CH₃CN|-3612.348|N|1.1600|0.0000|0.0000|-0.3300|
H が正、ハロゲン原子が負となり、結合極性と整合する。HF の電荷分離は ±0.264、HCl では ±0.150 であり、電気陰性度差の大きい HF の方が強く分極している。電気陰性度の高い中心原子側が負となり、環境の等価な H にはほぼ同じ値が割り当てられている。 CH₃OHではメチル基の H にも正電荷が割り当てられているが、3つの H は完全には等価ではない。ごく僅かな差ではあるが、これは OH 基に対する空間的な相対位置の違いを反映している可能性がある。事実、メチル基の水素の環境が完全に等価となるCH₃Clや直線状の置換基CNを有するCH₃CNでは H のlatent chargeがいずれも同じ値となっている。 ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697094/a56b6b9e-e457-4f6d-9875-bdc41296a4c4.png) ### イオン性分子
原子ごとの結果
|system_id|energy_eV|element|x|y|z|latent_charge| |:----|----:|:----|----:|----:|----:|----:| |NH₄⁺|-1547.854|N|0.0000|0.0000|0.0000|-0.0489| |NH₄⁺|-1547.854|H|0.6004|0.6004|0.6004|0.2618| |NH₄⁺|-1547.854|H|0.6004|-0.6004|-0.6004|0.2626| |NH₄⁺|-1547.854|H|-0.6004|0.6004|-0.6004|0.2621| |NH₄⁺|-1547.854|H|-0.6004|-0.6004|0.6004|0.2623| |H₃O⁺|-2087.098|O|0.0000|0.0000|0.0000|-0.0169| |H₃O⁺|-2087.098|H|0.8882|0.0000|0.4142|0.3387| |H₃O⁺|-2087.098|H|-0.4441|0.7692|0.4142|0.3391| |H₃O⁺|-2087.098|H|-0.4441|-0.7692|0.4142|0.3391| |OH⁻|-2062.795|O|0.0000|0.0000|0.0000|-0.9796| |OH⁻|-2062.795|H|0.9700|0.0000|0.0000|-0.0204| |CH₃COO⁻|-6218.877|C|-1.5000|0.0000|0.0000|-0.0860| |CH₃COO⁻|-6218.877|H|-1.8633|1.0277|0.0000|0.0192| |CH₃COO⁻|-6218.877|H|-1.8633|-0.5138|-0.8900|0.0201| |CH₃COO⁻|-6218.877|H|-1.8633|-0.5138|0.8900|0.0203| |CH₃COO⁻|-6218.877|C|0.0000|0.0000|0.0000|0.1748| |CH₃COO⁻|-6218.877|O|0.6200|1.0800|0.0000|-0.5710| |CH₃COO⁻|-6218.877|O|0.6200|-1.0800|0.0000|-0.5775| |NO₃⁻|-7630.682|N|0.0000|0.0000|0.0000|0.2311| |NO₃⁻|-7630.682|O|1.2400|0.0000|0.0000|-0.4108| |NO₃⁻|-7630.682|O|-0.6200|1.0739|0.0000|-0.4102| |NO₃⁻|-7630.682|O|-0.6200|-1.0739|0.0000|-0.4100| |CO₃²⁻|-7176.325|C|0.0000|0.0000|0.0000|-0.0723| |CO₃²⁻|-7176.325|O|1.2900|0.0000|0.0000|-0.6454| |CO₃²⁻|-7176.325|O|-0.6450|1.1172|0.0000|-0.6435| |CO₃²⁻|-7176.325|O|-0.6450|-1.1172|0.0000|-0.6387|
NH₄⁺・H₃O⁺では総電荷 +1 が主に H 原子へ分配されており、各 H の値はほぼ等しい。中心の N や O はわずかに負となった。OH⁻では負電荷のほぼ全てが O に局在しており、化学的直観とよく整合する結果である。CH₃COO⁻の2つの O はそれぞれ約 −0.57 であり、負電荷がほぼ均等に分配されている。これは酢酸イオンにおける共鳴・電荷非局在化を反映した結果と解釈できる。 NO₃⁻とCO₃²⁻では3つの O の latent charge はほぼ同じ値をとっている。ここで、CO₃²⁻では中心の C が −0.072 となっており、炭素に負の latent charge が割り当てられている。負電荷の密度が高い分子の場合には、一般によく知られる原子電荷とは異なる挙動を示す可能性がある。 > 参考までに、ORCA（ωB97M-V/def2-TZVPD）で得られる原子電荷は以下のようになる。数値上は latent charge はHirshfeld電荷に近い値を与えているように見える。 > ``` > ----------------------- > MULLIKEN ATOMIC CHARGES > ----------------------- > 0 C : 0.542142 > 1 O : -0.847378 > 2 O : -0.847382 > 3 O : -0.847382 > > ---------------------- > LOEWDIN ATOMIC CHARGES > ---------------------- > 0 C : -1.094417 > 1 O : -0.301859 > 2 O : -0.301862 > 3 O : -0.301862 > > ------------------ > HIRSHFELD ANALYSIS > ------------------ > 0 C : 0.018844 > 1 O : -0.672017 > 2 O : -0.672019 > 3 O : -0.672020 > ``` ### ハロゲン化アルカリ金属 次の節で詳細を述べるので、ここでは後回しにする。 ### 電解質 イオンへの水、DME（1,2-ジメトキシエタン）の配位構造で計算してみる。 > ここでのDMEはジメチルエーテルではない。
原子ごとの結果
|system_id|energy_eV|element|x|y|z|latent_charge| |:----|----:|:----|----:|----:|----:|----:| |Li(H₂O)⁺|-2279.350|Li|0.1235|0.0001|0.0002|0.7823| |Li(H₂O)⁺|-2279.350|O|1.9602|-0.0001|-0.0005|-0.2724| |Li(H₂O)⁺|-2279.350|H|2.5448|-0.7690|0.0000|0.2453| |Li(H₂O)⁺|-2279.350|H|2.5445|0.7689|0.0003|0.2449| |Li(H₂O)₂⁺|-4360.505|Li|0.0000|-0.0003|-0.0006|0.6308| |Li(H₂O)₂⁺|-4360.505|O|1.8665|-0.0004|-0.0002|-0.2886| |Li(H₂O)₂⁺|-4360.505|H|2.4500|-0.7691|-0.0007|0.2367| |Li(H₂O)₂⁺|-4360.505|H|2.4489|0.7695|0.0008|0.2365| |Li(H₂O)₂⁺|-4360.505|O|-1.8665|0.0005|0.0021|-0.2879| |Li(H₂O)₂⁺|-4360.505|H|-2.4499|0.7691|-0.0007|0.2360| |Li(H₂O)₂⁺|-4360.505|H|-2.4490|-0.7693|-0.0007|0.2364| |Li(H₂O)₄⁺|-8521.974|Li|-0.0004|-0.0013|-0.0003|0.4826| |Li(H₂O)₄⁺|-8521.974|O|1.1716|1.1703|1.0341|-0.3070| |Li(H₂O)₄⁺|-8521.974|H|1.9833|0.8982|1.4759|0.2181| |Li(H₂O)₄⁺|-8521.974|H|0.8984|1.9832|1.4737|0.2180| |Li(H₂O)₄⁺|-8521.974|O|1.1715|-1.1705|-1.0342|-0.3066| |Li(H₂O)₄⁺|-8521.974|H|0.8987|-1.9831|-1.4741|0.2181| |Li(H₂O)₄⁺|-8521.974|H|1.9838|-0.8982|-1.4748|0.2181| |Li(H₂O)₄⁺|-8521.974|O|-1.1719|1.1713|-1.0348|-0.3070| |Li(H₂O)₄⁺|-8521.974|H|-0.8996|1.9841|-1.4747|0.2182| |Li(H₂O)₄⁺|-8521.974|H|-1.9833|0.8982|-1.4768|0.2181| |Li(H₂O)₄⁺|-8521.974|O|-1.1697|-1.1701|1.0363|-0.3072| |Li(H₂O)₄⁺|-8521.974|H|-1.9832|-0.8980|1.4747|0.2182| |Li(H₂O)₄⁺|-8521.974|H|-0.8990|-1.9841|1.4751|0.2184| |Cl(H₂O)⁻|-14604.756|Cl|0.0682|0.0000|-0.0003|-0.8555| |Cl(H₂O)⁻|-14604.756|O|3.2165|-0.0002|0.0001|-0.4458| |Cl(H₂O)⁻|-14604.756|H|2.5715|0.7227|0.0002|0.1506| |Cl(H₂O)⁻|-14604.756|H|2.5708|-0.7224|0.0000|0.1508| |Cl(H₂O)₂⁻|-16685.134|Cl|0.0000|0.0000|0.0000|-0.7608| |Cl(H₂O)₂⁻|-16685.134|O|3.2321|0.0000|0.0001|-0.4440| |Cl(H₂O)₂⁻|-16685.134|H|2.5958|0.7285|0.0000|0.1621| |Cl(H₂O)₂⁻|-16685.134|H|2.5956|-0.7285|0.0000|0.1623| |Cl(H₂O)₂⁻|-16685.134|O|-3.2322|0.0000|-0.0001|-0.4441| |Cl(H₂O)₂⁻|-16685.134|H|-2.5957|-0.7285|-0.0001|0.1622| |Cl(H₂O)₂⁻|-16685.134|H|-2.5957|0.7285|0.0001|0.1622| |Cl(H₂O)₄⁻|-20845.740|Cl|-0.0001|0.0001|0.0000|-0.6144| |Cl(H₂O)₄⁻|-20845.740|O|1.8658|1.8657|1.8615|-0.4285| |Cl(H₂O)₄⁻|-20845.740|H|0.9817|2.0188|1.5033|0.1659| |Cl(H₂O)₄⁻|-20845.740|H|2.0188|0.9817|1.5031|0.1656| |Cl(H₂O)₄⁻|-20845.740|O|1.8658|-1.8657|-1.8616|-0.4277| |Cl(H₂O)₄⁻|-20845.740|H|2.0187|-0.9817|-1.5031|0.1656| |Cl(H₂O)₄⁻|-20845.740|H|0.9817|-2.0187|-1.5032|0.1660| |Cl(H₂O)₄⁻|-20845.740|O|-1.8657|1.8657|-1.8615|-0.4275| |Cl(H₂O)₄⁻|-20845.740|H|-2.0188|0.9817|-1.5032|0.1657| |Cl(H₂O)₄⁻|-20845.740|H|-0.9817|2.0187|-1.5032|0.1657| |Cl(H₂O)₄⁻|-20845.740|O|-1.8657|-1.8658|1.8614|-0.4284| |Cl(H₂O)₄⁻|-20845.740|H|-0.9816|-2.0188|1.5032|0.1659| |Cl(H₂O)₄⁻|-20845.740|H|-2.0188|-0.9816|1.5033|0.1662| |Li(DME)⁺|-8604.637|Li|0.0000|0.0000|0.0000|0.6430| |Li(DME)⁺|-8604.637|C|1.3303|-0.3188|-2.2951|0.0231| |Li(DME)⁺|-8604.637|H|1.4446|-1.4056|-2.3155|0.0715| |Li(DME)⁺|-8604.637|H|1.5006|0.0749|-3.2990|0.0868| |Li(DME)⁺|-8604.637|O|0.0000|0.0000|-1.8637|-0.2339| |Li(DME)⁺|-8604.637|C|-0.9833|-0.2564|-2.8775|0.0139| |Li(DME)⁺|-8604.637|H|-0.8178|0.4089|-3.7242|0.0707| |Li(DME)⁺|-8604.637|H|-0.9232|-1.2943|-3.2058|0.0711| |Li(DME)⁺|-8604.637|H|-1.9608|-0.0608|-2.4464|0.0751| |Li(DME)⁺|-8604.637|C|2.3008|0.3137|-1.3278|0.0231| |Li(DME)⁺|-8604.637|H|2.3324|1.4002|-1.4465|0.0719| |Li(DME)⁺|-8604.637|H|3.3014|-0.0897|-1.4971|0.0874| |Li(DME)⁺|-8604.637|O|1.8607|0.0000|0.0012|-0.2345| |Li(DME)⁺|-8604.637|C|2.8638|0.2410|0.9994|0.0135| |Li(DME)⁺|-8604.637|H|3.6946|-0.4507|0.8625|0.0714| |Li(DME)⁺|-8604.637|H|3.2203|1.2698|0.9389|0.0716| |Li(DME)⁺|-8604.637|H|2.4058|0.0693|1.9691|0.0745| |Li(DME)₂⁺|-17010.387|Li|0.0000|0.0000|0.0000|0.4844| |Li(DME)₂⁺|-17010.387|C|1.2921|-0.3811|-2.4175|0.0268| |Li(DME)₂⁺|-17010.387|H|1.3834|-1.4726|-2.3835|0.0593| |Li(DME)₂⁺|-17010.387|H|1.4364|-0.0510|-3.4498|0.0764| |Li(DME)₂⁺|-17010.387|O|0.0000|0.0000|-1.9527|-0.2299| |Li(DME)₂⁺|-17010.387|C|-1.0398|-0.4263|-2.8312|0.0120| |Li(DME)₂⁺|-17010.387|H|-0.9408|0.0655|-3.8011|0.0625| |Li(DME)₂⁺|-17010.387|H|-1.0029|-1.5098|-2.9708|0.0590| |Li(DME)₂⁺|-17010.387|H|-1.9876|-0.1469|-2.3780|0.0602| |Li(DME)₂⁺|-17010.387|C|2.3084|0.2776|-1.5202|0.0267| |Li(DME)₂⁺|-17010.387|H|2.3257|1.3603|-1.6817|0.0626| |Li(DME)₂⁺|-17010.387|H|3.3026|-0.1271|-1.7263|0.0761| |Li(DME)₂⁺|-17010.387|O|1.9337|0.0000|-0.1735|-0.2292| |Li(DME)₂⁺|-17010.387|C|2.9442|0.3501|0.7709|0.0119| |Li(DME)₂⁺|-17010.387|H|3.8467|-0.2362|0.5874|0.0616| |Li(DME)₂⁺|-17010.387|H|3.1828|1.4145|0.7091|0.0618| |Li(DME)₂⁺|-17010.387|H|2.5581|0.1211|1.7611|0.0604| |Li(DME)₂⁺|-17010.387|C|-1.9351|0.2522|1.9762|0.0265| |Li(DME)₂⁺|-17010.387|H|-1.3356|0.5904|2.8275|0.0617| |Li(DME)₂⁺|-17010.387|H|-2.9877|0.4640|2.1814|0.0760| |Li(DME)₂⁺|-17010.387|O|-1.5138|0.9369|0.7982|-0.2302| |Li(DME)₂⁺|-17010.387|C|-1.8213|2.3306|0.8402|0.0116| |Li(DME)₂⁺|-17010.387|H|-2.9026|2.4773|0.8823|0.0623| |Li(DME)₂⁺|-17010.387|H|-1.3560|2.8011|1.7103|0.0616| |Li(DME)₂⁺|-17010.387|H|-1.4318|2.7802|-0.0695|0.0624| |Li(DME)₂⁺|-17010.387|C|-1.7475|-1.2244|1.7418|0.0263| |Li(DME)₂⁺|-17010.387|H|-2.4553|-1.5966|0.9930|0.0599| |Li(DME)₂⁺|-17010.387|H|-1.9061|-1.7713|2.6750|0.0767| |Li(DME)₂⁺|-17010.387|O|-0.4161|-1.4135|1.2734|-0.2289| |Li(DME)₂⁺|-17010.387|C|-0.0470|-2.7882|1.1892|0.0118| |Li(DME)₂⁺|-17010.387|H|-0.0511|-3.2430|2.1820|0.0623| |Li(DME)₂⁺|-17010.387|H|-0.7362|-3.3325|0.5383|0.0589| |Li(DME)₂⁺|-17010.387|H|0.9579|-2.8297|0.7758|0.0583|
溶媒和殻への latent charge の分配についてまとめると以下のようになる。この結果は、溶媒和数の増加に伴って、系の総電荷が中心イオンから溶媒分子へ段階的に分散する様子を、latent charge が表現していることを示している。 | 系 | 中心イオン | 溶媒和殻全体 | 溶媒1分子当たり | | --------- | ------: | ------: | ------: | | Li(H₂O)⁺ | +0.7823 | +0.2177 | +0.2177 | | Li(H₂O)₂⁺ | +0.6308 | +0.3692 | +0.1846 | | Li(H₂O)₄⁺ | +0.4826 | +0.5174 | +0.1294 | | Cl(H₂O)⁻ | −0.8555 | −0.1445 | −0.1445 | | Cl(H₂O)₂⁻ | −0.7608 | −0.2392 | −0.1196 | | Cl(H₂O)₄⁻ | −0.6144 | −0.3856 | −0.0964 | | Li(DME)⁺ | +0.6430 | +0.3570 | +0.3570 | | Li(DME)₂⁺ | +0.4844 | +0.5156 | +0.2578 | 例えば、Li(H₂O)⁺の系の水分子に着目すると、孤立水分子では O が約 −0.418、H が約 +0.209であったのに対し、Li⁺配位下では O の負電荷が小さく、H の正電荷が大きくなる。その結果、水分子全体が正の電荷を帯びている。これはLi⁺への配位によって、水分子内のlatent charge分布が再構成されていることを示している。 Cl(H₂O)⁻の系でも逆符号で同様の傾向を示している。特筆すべき点として、4水和系ではClには −0.614 が残り、Li(H₂O)₄⁺のLiの +0.483 よりも中心イオンへの局在性が高くなっている。これは、Li⁺系がOの孤立電子対との直接的な配位結合を形成するのに対し、Cl⁻系では主に水分子のHを介した水素結合性相互作用になるという局所環境の違いと整合する結果だと言える。 > あくまでも相対的な値の関係がそうなっているように見えるというだけで、絶対値から電荷分離の強度を評価してよいかについては議論の余地がある。 また、水とDME（1,2-ジメトキシエタン）を比較すると、2個のOドナーを持つ系同士、4個のOドナーを持つ系同士で、Liの latent charge がほぼ一致している。 | Oドナー数 | 水系 | Liの電荷 | DME系 | Liの電荷 | | ----: | --------- | ------: | --------- | ------: | | 2 | Li(H₂O)₂⁺ | +0.6308 | Li(DME)⁺ | +0.6430 | | 4 | Li(H₂O)₄⁺ | +0.4826 | Li(DME)₂⁺ | +0.4844 | この結果は、Liの latent charge が溶媒分子の種類や原子数よりも、第一配位圏にあるOドナー数とLi–O距離に強く支配されている可能性を示しており、OrbMol-v2の局所環境の表現として自然な挙動に思われる。 > DMEは1分子当たり2個のO原子でLi⁺に配位している。Li(DME)₂⁺の場合は以下のような配位構造をとる。 > ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697094/d8b8b439-808d-4088-9679-4f7d1dd7a496.png) > > Li–O距離も配位数に対応していた。配位数が少ないほどより強く分極して接近する、という描像が適切に表現できていると考えられる。 > - Li(H₂O)₂⁺：約 1.87 Å > - Li(DME)⁺：約 1.86 Å > - Li(H₂O)₄⁺：約 1.95 Å > - Li(DME)₂⁺：おおむね 1.94–1.95 Å の範囲 > ちなみに、ChatGPTの用意した初期構造とOrbMolでの最適化後の構造では、latent charge の変化は 0.01–0.03e 程度に留まり、基本的な傾向は変わらなかった。latent charge の振る舞いは構造に対して一定の堅牢性を有しているようである。ただし、ChatGPTの作成した初期構造がある程度妥当なものだったという点も考慮すべきである。 ### 金属錯体
原子ごとの結果
|system_id|energy_eV|element|x|y|z|latent_charge|charge|spin| |:----|:----|:----|:----|:----|:----|:----|:----|:----| |ferrocene|-44918.855|Fe|0.0005|-0.0005|-0.0012|0.3527|0|1| |ferrocene|-44918.855|C|1.2081|-0.0002|1.6717|-0.0949|0|1| |ferrocene|-44918.855|C|0.3731|1.1487|1.6700|-0.0954|0|1| |ferrocene|-44918.855|C|-0.9776|0.7098|1.6692|-0.0955|0|1| |ferrocene|-44918.855|C|-0.9774|-0.7103|1.6702|-0.0963|0|1| |ferrocene|-44918.855|C|0.3732|-1.1490|1.6715|-0.0957|0|1| |ferrocene|-44918.855|H|2.2856|0.0000|1.6417|0.0604|0|1| |ferrocene|-44918.855|H|0.7066|2.1732|1.6418|0.0603|0|1| |ferrocene|-44918.855|H|-1.8495|1.3429|1.6413|0.0605|0|1| |ferrocene|-44918.855|H|-1.8490|-1.3435|1.6412|0.0606|0|1| |ferrocene|-44918.855|H|0.7063|-2.1737|1.6412|0.0603|0|1| |ferrocene|-44918.855|C|0.9773|0.7104|-1.6709|-0.0957|0|1| |ferrocene|-44918.855|C|-0.3733|1.1493|-1.6709|-0.0963|0|1| |ferrocene|-44918.855|C|-1.2082|0.0005|-1.6703|-0.0969|0|1| |ferrocene|-44918.855|C|-0.3736|-1.1487|-1.6701|-0.0964|0|1| |ferrocene|-44918.855|C|0.9774|-0.7097|-1.6698|-0.0953|0|1| |ferrocene|-44918.855|H|1.8493|1.3430|-1.6422|0.0606|0|1| |ferrocene|-44918.855|H|-0.7062|2.1739|-1.6419|0.0605|0|1| |ferrocene|-44918.855|H|-2.2857|0.0000|-1.6408|0.0609|0|1| |ferrocene|-44918.855|H|-0.7061|-2.1734|-1.6405|0.0609|0|1| |ferrocene|-44918.855|H|1.8492|-1.3428|-1.6411|0.0607|0|1| |ferrocenium|-44912.554|Fe|0.0022|-0.0004|0.0026|0.5190|1|2| |ferrocenium|-44912.554|C|1.2074|0.0006|1.7101|-0.0526|1|2| |ferrocenium|-44912.554|C|0.3730|1.1491|1.7119|-0.0528|1|2| |ferrocenium|-44912.554|C|-0.9771|0.7102|1.7077|-0.0526|1|2| |ferrocenium|-44912.554|C|-0.9772|-0.7094|1.7046|-0.0530|1|2| |ferrocenium|-44912.554|C|0.3732|-1.1484|1.7053|-0.0534|1|2| |ferrocenium|-44912.554|H|2.2855|0.0014|1.6906|0.1011|1|2| |ferrocenium|-44912.554|H|0.7072|2.1741|1.6941|0.1014|1|2| |ferrocenium|-44912.554|H|-1.8496|1.3435|1.6879|0.1014|1|2| |ferrocenium|-44912.554|H|-1.8500|-1.3424|1.6852|0.1015|1|2| |ferrocenium|-44912.554|H|0.7070|-2.1735|1.6859|0.1012|1|2| |ferrocenium|-44912.554|C|0.9771|0.7096|-1.7057|-0.0540|1|2| |ferrocenium|-44912.554|C|-0.3734|1.1481|-1.7049|-0.0549|1|2| |ferrocenium|-44912.554|C|-1.2079|-0.0006|-1.7065|-0.0543|1|2| |ferrocenium|-44912.554|C|-0.3737|-1.1487|-1.7084|-0.0538|1|2| |ferrocenium|-44912.554|C|0.9763|-0.7099|-1.7066|-0.0530|1|2| |ferrocenium|-44912.554|H|1.8503|1.3419|-1.6907|0.1015|1|2| |ferrocenium|-44912.554|H|-0.7069|2.1735|-1.6890|0.1013|1|2| |ferrocenium|-44912.554|H|-2.2862|-0.0010|-1.6898|0.1020|1|2| |ferrocenium|-44912.554|H|-0.7062|-2.1745|-1.6935|0.1023|1|2| |ferrocenium|-44912.554|H|1.8490|-1.3432|-1.6909|0.1017|1|2| |TiCl₄|-73210.739|Ti|-0.0005|0.0000|-0.0007|0.9591|0|1| |TiCl₄|-73210.739|Cl|1.2496|1.2493|1.2494|-0.2411|0|1| |TiCl₄|-73210.739|Cl|1.2495|-1.2492|-1.2490|-0.2387|0|1| |TiCl₄|-73210.739|Cl|-1.2492|1.2493|-1.2491|-0.2402|0|1| |TiCl₄|-73210.739|Cl|-1.2495|-1.2495|1.2495|-0.2391|0|1| |ZrCl₄|-51377.029|Zr|-0.0002|-0.0001|-0.0001|1.1914|0|1| |ZrCl₄|-51377.029|Cl|1.3402|1.3401|1.3400|-0.2992|0|1| |ZrCl₄|-51377.029|Cl|1.3402|-1.3400|-1.3401|-0.2967|0|1| |ZrCl₄|-51377.029|Cl|-1.3399|1.3402|-1.3400|-0.2983|0|1| |ZrCl₄|-51377.029|Cl|-1.3403|-1.3402|1.3403|-0.2972|0|1| |Cp₂ZrCl₂|-36862.206|Zr|0.8887|0.0013|-0.0011|1.0140|0|1| |Cp₂ZrCl₂|-36862.206|C|1.0229|0.0004|2.5288|-0.1006|0|1| |Cp₂ZrCl₂|-36862.206|C|0.2723|1.1417|2.1739|-0.0908|0|1| |Cp₂ZrCl₂|-36862.206|C|-0.9224|0.7105|1.5662|-0.1102|0|1| |Cp₂ZrCl₂|-36862.206|C|-0.9221|-0.7110|1.5665|-0.1090|0|1| |Cp₂ZrCl₂|-36862.206|C|0.2731|-1.1412|2.1739|-0.0908|0|1| |Cp₂ZrCl₂|-36862.206|H|2.0049|0.0011|2.9730|0.0759|0|1| |Cp₂ZrCl₂|-36862.206|H|0.5862|2.1647|2.2961|0.0751|0|1| |Cp₂ZrCl₂|-36862.206|H|-1.7084|1.3459|1.1891|0.0724|0|1| |Cp₂ZrCl₂|-36862.206|H|-1.7071|-1.3472|1.1888|0.0724|0|1| |Cp₂ZrCl₂|-36862.206|H|0.5871|-2.1641|2.2969|0.0754|0|1| |Cp₂ZrCl₂|-36862.206|C|0.8094|0.6992|-2.4404|-0.0923|0|1| |Cp₂ZrCl₂|-36862.206|C|-0.3907|1.1403|-1.8271|-0.1076|0|1| |Cp₂ZrCl₂|-36862.206|C|-1.1429|-0.0012|-1.4780|-0.1001|0|1| |Cp₂ZrCl₂|-36862.206|C|-0.3907|-1.1428|-1.8265|-0.1073|0|1| |Cp₂ZrCl₂|-36862.206|C|0.8093|-0.7022|-2.4400|-0.0927|0|1| |Cp₂ZrCl₂|-36862.206|H|1.6033|1.3338|-2.7962|0.0757|0|1| |Cp₂ZrCl₂|-36862.206|H|-0.6801|2.1673|-1.6682|0.0738|0|1| |Cp₂ZrCl₂|-36862.206|H|-2.1152|-0.0009|-1.0119|0.0705|0|1| |Cp₂ZrCl₂|-36862.206|H|-0.6806|-2.1698|-1.6679|0.0743|0|1| |Cp₂ZrCl₂|-36862.206|H|1.6037|-1.3367|-2.7952|0.0760|0|1| |Cp₂ZrCl₂|-36862.206|Cl|2.4445|1.8619|-0.0014|-0.3783|0|1| |Cp₂ZrCl₂|-36862.206|Cl|2.4549|-1.8509|0.0006|-0.3756|0|1| |Cp₂ZrMe⁺|-12891.928|Zr|0.6010|-0.0185|0.0027|1.0670|1|1| |Cp₂ZrMe⁺|-12891.928|C|0.9749|-0.0308|2.4674|-0.0736|1|1| |Cp₂ZrMe⁺|-12891.928|C|0.2250|1.1312|2.1808|-0.0802|1|1| |Cp₂ZrMe⁺|-12891.928|C|-1.0465|0.7317|1.6989|-0.0903|1|1| |Cp₂ZrMe⁺|-12891.928|C|-1.0858|-0.6843|1.6973|-0.0859|1|1| |Cp₂ZrMe⁺|-12891.928|C|0.1633|-1.1532|2.1670|-0.0817|1|1| |Cp₂ZrMe⁺|-12891.928|H|1.9763|-0.0611|2.8652|0.0937|1|1| |Cp₂ZrMe⁺|-12891.928|H|0.5641|2.1487|2.3129|0.0939|1|1| |Cp₂ZrMe⁺|-12891.928|H|-1.8568|1.3890|1.4217|0.0919|1|1| |Cp₂ZrMe⁺|-12891.928|H|-1.9247|-1.2982|1.4065|0.0916|1|1| |Cp₂ZrMe⁺|-12891.928|H|0.4431|-2.1895|2.2961|0.0947|1|1| |Cp₂ZrMe⁺|-12891.928|C|0.7436|0.7138|-2.3802|-0.0746|1|1| |Cp₂ZrMe⁺|-12891.928|C|-0.5259|1.1163|-1.8850|-0.0890|1|1| |Cp₂ZrMe⁺|-12891.928|C|-1.2817|-0.0504|-1.6219|-0.0848|1|1| |Cp₂ZrMe⁺|-12891.928|C|-0.4726|-1.1696|-1.9131|-0.0873|1|1| |Cp₂ZrMe⁺|-12891.928|C|0.7781|-0.6927|-2.3921|-0.0780|1|1| |Cp₂ZrMe⁺|-12891.928|H|1.5400|1.3698|-2.6953|0.0941|1|1| |Cp₂ZrMe⁺|-12891.928|H|-0.8758|2.1334|-1.7841|0.0928|1|1| |Cp₂ZrMe⁺|-12891.928|H|-2.2988|-0.0796|-1.2649|0.0901|1|1| |Cp₂ZrMe⁺|-12891.928|H|-0.7686|-2.2052|-1.8292|0.0940|1|1| |Cp₂ZrMe⁺|-12891.928|H|1.6091|-1.3018|-2.7147|0.0944|1|1| |Cp₂ZrMe⁺|-12891.928|C|2.8065|0.1346|-0.0096|-0.2941|1|1| |Cp₂ZrMe⁺|-12891.928|H|2.9819|-0.9645|0.0692|0.0349|1|1| |Cp₂ZrMe⁺|-12891.928|H|3.2927|0.4496|-0.9318|0.0429|1|1| |Cp₂ZrMe⁺|-12891.928|H|3.3275|0.5813|0.8362|0.0435|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|Zr|0.9678|0.2302|-0.0195|1.0660|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|C|1.3003|-0.0665|2.4336|-0.0729|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|C|0.5221|1.1000|2.2685|-0.0832|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|C|-0.7315|0.7239|1.7208|-0.0912|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|C|-0.7288|-0.6799|1.5544|-0.0862|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|C|0.5286|-1.1670|1.9839|-0.0826|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|H|2.2960|-0.1151|2.8411|0.0925|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|H|0.8282|2.1036|2.5287|0.0931|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|H|-1.5555|1.3822|1.4907|0.0837|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|H|-1.5480|-1.2695|1.1735|0.0816|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|H|0.8377|-2.2022|1.9942|0.0926|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|C|1.1498|0.6873|-2.4668|-0.0791|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|C|-0.1517|1.0921|-2.0591|-0.0883|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|C|-0.8658|-0.0574|-1.6623|-0.0860|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|C|0.0007|-1.1716|-1.7911|-0.0912|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|C|1.2381|-0.7077|-2.3125|-0.0743|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|H|1.9322|1.3326|-2.8356|0.0931|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|H|-0.5389|2.1010|-2.0783|0.0903|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|H|-1.8826|-0.0791|-1.2992|0.0764|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|H|-0.2550|-2.2005|-1.5866|0.0890|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|H|2.1005|-1.3142|-2.5410|0.0933|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|C|3.1880|0.1039|0.0127|-0.2982|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|H|3.3501|-0.9834|0.1349|0.0357|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|H|3.7060|0.3895|-0.9035|0.0408|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|H|3.6872|0.5812|0.8579|0.0413|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|C|-4.2220|0.6923|0.0854|-0.1019|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|C|-4.2742|-0.6314|0.0983|-0.1013|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|H|-4.2164|1.2641|1.0052|0.0658|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|H|-4.2320|1.2564|-0.8400|0.0666|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|H|-4.3114|-1.1853|1.0284|0.0673|1|1| |Cp₂ZrMe⁺+ethylene|-15030.016|H|-4.3297|-1.2095|-0.8166|0.0675|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|Zr|0.7793|0.4615|-0.0601|0.9767|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|C|1.2847|-0.2313|2.2913|-0.0711|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|C|0.4611|0.9044|2.3871|-0.0947|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|C|-0.8197|0.5690|1.8705|-0.0916|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|C|-0.7796|-0.7728|1.4518|-0.0854|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|C|0.5302|-1.2644|1.6832|-0.0918|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|H|2.3096|-0.2998|2.6154|0.0926|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|H|0.7458|1.8514|2.8191|0.0881|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|H|-1.6842|1.2150|1.8346|0.0894|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|H|-1.6022|-1.3314|1.0332|0.0890|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|H|0.8783|-2.2668|1.4834|0.0910|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|C|1.1558|0.4395|-2.5533|-0.0877|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|C|-0.1801|0.8565|-2.3487|-0.0944|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|C|-0.8911|-0.2125|-1.7620|-0.0950|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|C|0.0091|-1.2978|-1.6118|-0.0910|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|C|1.2691|-0.8936|-2.1070|-0.0761|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|H|1.9457|1.0271|-2.9946|0.0916|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|H|-0.5944|1.8155|-2.6233|0.0878|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|H|-1.9394|-0.2166|-1.5056|0.0894|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|H|-0.2403|-2.2734|-1.2266|0.0892|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|H|2.1620|-1.4963|-2.1393|0.0935|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|C|3.0275|0.4677|0.0083|-0.2981|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|H|3.4007|-0.5489|0.1722|0.0380|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|H|3.4800|0.8233|-0.9203|0.0361|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|H|3.4144|1.0745|0.8323|0.0358|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|C|1.4029|3.1408|-0.2565|-0.0528|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|C|0.1129|3.1889|0.0960|-0.0564|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|H|2.1959|3.2389|0.4747|0.0780|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|H|1.7000|3.0724|-1.2964|0.0724|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|H|-0.1843|3.3528|1.1232|0.0743|1|1| |Cp₂ZrMe(ethylene)⁺|-15030.989|H|-0.6772|3.1667|-0.6467|0.0734|1|1|
等価なCp環、Cl原子および配位子にはほぼ等しい値が割り当てられ、構造対称性が反映されていることが観察された。 まとめると以下のようになる。 | 系 | 金属 | Cp環1 | Cp環2 | その他の配位子 | | ------------------- | ---------: | ------: | ------: | --------------------------: | | ferrocene | Fe +0.3527 | −0.1757 | −0.1770 | — | | ferrocenium | Fe +0.5190 | +0.2422 | +0.2388 | — | | TiCl₄ | Ti +0.9591 | — | — | Cl合計 −0.9591 | | ZrCl₄ | Zr +1.1914 | — | — | Cl合計 −1.1914 | | Cp₂ZrCl₂ | Zr +1.0140 | −0.1302 | −0.1297 | Cl合計 −0.7539 | | Cp₂ZrMe⁺ | Zr +1.0670 | +0.0541 | +0.0517 | Me −0.1728 | | Cp₂ZrMe⁺ + ethylene | Zr +1.0660 | +0.0274 | +0.0232 | Me −0.1804、ethylene +0.0640 | | Cp₂ZrMe(ethylene)⁺ | Zr +0.9767 | +0.0155 | +0.0073 | Me −0.1882、ethylene +0.1889 | ジルコノセン錯体とエチレンの系について、Zr中心からエチレンが約5.2 Å離れた配置（Cp₂ZrMe⁺ + ethylene）では、Zrのlatent chargeは (+1.0660) であり、単離したCp₂ZrMe⁺の値とほぼ同じであった。この構造においてエチレンは化学的には解離しているので、妥当な結果である。 また、エチレンがZrに完全に配位したCp₂ZrMe(ethylene)⁺では、エチレンの2つの炭素とZrの距離はそれぞれ約2.76 Åおよび2.81 Åであり、C=C中点とZrの距離は約2.70 Åであった。二つのZr–C距離が近いことから、構造は明瞭な $\eta^2$-ethylene型配位を示している。このとき、エチレン全体のgroup chargeは (+0.0640) から (+0.1889) へ増加した。この変化は、電子不足なZr中心へエチレンの電子密度が供与されるという描像と定性的に整合する。[^C=C] ### 全体的な傾向 結果をまとめると以下の通り。結論として、 latent charge は電荷に相当する量として自然な振る舞いを示していると言える。 - latent charge の符号と相対的な大きさは一般的な極性・電気陰性度の大小と整合的である。 - latent charge には分子構造の対称性が反映されている。これは無極性分子だけでなく共鳴構造を有するイオンでも同様であり、分子間であっても成立している。 - 溶媒和数が増えるほど中心イオンの latent charge の絶対値は小さくなる。溶媒分子の分極や電荷の非局在化をモデル内部で適切に表現できている可能性が高い。 - 錯形成によるPESの安定化と金属–配位子間のlatent charge再分配が連動していることが確認できた。 物理的な背景を有する電荷ではないにせよ、latent charge は物理的にかなり妥当な振る舞いを示しているように思われる。 ## 原子間距離と latent charge の関係 続いて、ハロゲン化アルカリ金属の解離において latent charge がどのように振舞うのかを確認する。 ### NaCl ![scan_NaCl_latent_charge.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697094/0a815e4b-7751-41d1-bdb9-f5fc561a2c42.png) ![scan_NaCl_energy_gradient.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697094/1b515971-c4bb-44c3-bcf4-5572c5e8f88c.png) ### KBr ![scan_KBr_latent_charge.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697094/bdb95e2e-9177-4ee9-8976-1831857aa2c3.png) ![scan_KBr_energy_gradient.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697094/104261c1-ed7e-4b54-ad0c-c5bfe59e312e.png) ### LiF ![scan_LiF_latent_charge.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697094/d27d6e5a-4f80-4cc7-9d72-106ab4b853b1.png) ![scan_LiF_energy_gradient.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697094/ee8a4a82-1e88-4df4-bc8f-ae47e4731f6c.png) ### 全体的な傾向 これらの結果は重要な示唆を含んでいる。 アルカリハライドの latent charge は原子間距離とともに孤立イオンの値へ近づき、カットオフ半径である 6 Å以上では一定値となった。結合伸長に伴ってアルカリ金属とハロゲンの電荷分離が増大するような内部表現を形成している点は、電荷の振る舞いとして自然である。 一方、その極限値は NaCl、KBr、LiF でそれぞれ約 ±0.97、±0.95、±0.89 であり、電気素量（±1.0）には一致しなかった。このような単純なイオンの場合でも、latent charge によって真空中の完全な電離を議論するには不十分であることが示唆された。[^latent] また、LiFではlatent chargeがNaClやKBrより短い 4 Å程度の距離で一定となる傾向が見られた。これはカットオフ半径だけでは説明できず、元素ペアごとに学習された実効的な距離依存性が反映されたものと考えられる。計算した3種類のイオンのうち、LiFはサイズが小さく平衡核間距離の短いイオンである。平衡核間距離で規格化すればNaClやKBrと同程度の相対的解離距離で飽和している可能性がある。 > 実際、原子間距離を絶対値ではなく、おおよその平衡結合長 $r_e$ で規格化して比較すると、latent chargeがsaturateする距離は > 　$\mathrm{NaCl}:6 \ \mathring{\mathrm{A}} ≈2.5r_e$ > 　$\mathrm{KBr}:7 \ \mathring{\mathrm{A}} ≈2.5r_e$ > 　$\mathrm{LiF}:4 \ \mathring{\mathrm{A}} ≈2.5r_e$​ > となっている。仮説の範疇だが、相対的な結合伸長率に対する latent charge の振る舞いとしては、どのイオン分子でもそれほど異なっていないのかもしれない。 一方、極限値が約±0.89に留まる点は原子間距離のスケーリングだけでは説明できない。LiFの解離極限におけるlatent chargeは、学習の過程でポテンシャル面の再現性が優先されたために中途半端な値を示すようになったものと思われる。 また、ポテンシャルエネルギーがカットオフ半径である 6 Åを超えても増大している点は興味深い。これは明らかにlatent chargeを用いた 1/R 型の相互作用によってPESの尾部（Coulomb tail）が形成されていると考えられる。このような性質は他の多くの汎用機械学習ポテンシャルでは正確に記述できないため、より正確に物理的性質を反映しているという点でOrbMolの優位性の一つと言えるだろう。 ## UMA、DFTとの比較 アンモニウムイオンのプロトンを窒素原子から引き離していって得られる構造に対して、OrbMol-v2、UMA-s-1.2、DFT（ORCA；ωB97M-V/def2-TZVPD）でエネルギー、勾配、電荷（UMAを除く）を比較した。 > 参考までに、N–H間距離が 4 Åの構造は下図のようになっている。 > ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697094/5b435221-815b-4f9b-8e02-4d15143aa4e1.png) ![comparison_relative_energy.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697094/38d1bbd1-66d7-4bcc-95fe-b2e131c62f68.png) ![comparison_gradient.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697094/c655b989-aafa-421e-a850-4ceb1f704f29.png) ![comparison_group_charges.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697094/ba42c158-d4fa-40e2-9c8e-70d6b75bb23e.png) まず注目すべき点として、UMA-s-1.2のエネルギー値は2Å程度の距離までしかDFT計算の結果を再現しないのに対し、OrbMol-v2では3Å程度までDFT計算を再現していることが挙げられる。NH4+という非常に単純な分子であっても、このような違いが見られるのは興味深い。 電荷の振る舞いについては、なお不明な点が多い。脱離するH原子のlatent chargeは、N–H間距離の増加に伴って約2.5–3 Åまでは増加するものの、それ以降は減少へ転じ、6 Å以降では約 (+0.08) の一定値を示した。これに対し、ORCAから得られたHirshfeld電荷およびMulliken電荷は長距離側でも増加を続け、約 (+0.64) に達した。したがって、OrbMol-v2のlatent chargeは、通常の電子密度分割に基づく部分電荷とは定性的にも異なる解離挙動を示した。 latent chargeが極大を示して減少へ転じる約3 Å付近から、OrbMol-v2とORCAの相対エネルギー曲線にも差が現れ始めた。OrbMol-v2は短距離側ではORCAのPESを良好に再現する一方、長距離側では解離エネルギーを約1 eV以上高く評価した。このDFTとの乖離の一致は、latent chargeの異常な再分配とPESのずれが、結合解離領域における共通のモデル誤差に由来する可能性を示唆している。 :::note warn #### NH₄⁺の解離極限について NH₄⁺の物理的な解離極限はヘテロリティック（中性NH₃とH⁺への解離）であると思われがちだが、実は気相中ではホモリティック解離（NH₃⁺ラジカルとHラジカルへの解離）の方が約3.5 eV程度、熱力学的に安定である。 OrbMol-v2では長距離領域でも正電荷の大部分がNH₃フラグメント側に残っており、ホモリティックな解離過程を反映しているようにも思われる。しかしながら、ホモリティック解離の生成物は2つのラジカル種となるため、全体としては開殻一重項または三重項として記述する必要がある。 今回、ORCAの計算は閉殻一重項として計算しているため、ヘテロリティックな解離過程が得られていると考えられる。実際、[アンモニアのプロトン親和力](https://webbook.nist.gov/cgi/cbook.cgi?ID=C7664417)は約 8.85 eV であり、先ほどのエネルギープロファイルから推定される解離エネルギーにほぼ一致することから、計算上はヘテロリティック解離として記述されていると推察される。 つまり、OrbMol-v2のlatent chargeは、物理的な原子電荷やフラグメント電荷とは必ずしも対応しない内部変数として振る舞っている可能性が高い。上で見たように、OrbMol-v2でのNH₄⁺解離では、遠方の水素に割り当てられるlatent chargeは中性水素に近づく一方、予測エネルギーはイオン性解離極限に近い値を示した。この結果は、latent chargeによるフラグメント間の電荷配分と、モデルが予測するポテンシャルエネルギーとの間に、物理的に整合しない挙動が生じ得ることを示唆している。 > 解離エネルギーの値は、イオン化エネルギーやプロトン親和力などの実験値を用いて「ヘスの法則」から求められる。因みに、NH₄⁺のホモリティック解離エネルギーは約 5.32 eV と計算できる。 ::: ## サイズに対する無矛盾性・示量性 冒頭でも述べた通り、一般にNNPは局所的な環境のみを入力として推論を行うため、「サイズに対する無矛盾性（size consistency）」や「サイズに対する示量性（size extensivity）」が正確に表現できないという問題がある。latent chargeは系全体にわたって総電荷が保たれるような制約を満たすため、この問題を克服している可能性がある。[^cutoff] そこで、NH₃、NH₄ラジカル、NH₄⁺カチオンを100 Å間隔で配置した構造で1分子あたりのエネルギーを求めた。エネルギーの単位は eV である。
OrbMol の結果（数値表）
|species|n_fragments|charge|multiplicity|formal_interfragment_coulomb_eV|total_energy_eV| |:----|:----|:----|:----|:----|:----| |NH₃|1|0|1|0.000|-1538.812| |NH₃|2|0|1|0.000|-3077.625| |NH₃|3|0|1|0.000|-4616.437| |NH₃|4|0|1|0.000|-6155.249| |NH₃|5|0|1|0.000|-7694.061| |NH₃|6|0|1|0.000|-9232.875| |NH₃|8|0|1|0.000|-12310.499| |NH₃|10|0|1|0.000|-15388.124| |NH₃|12|0|1|0.000|-18465.747| |NH₃|16|0|1|0.000|-24620.998| |NH₄⁺|1|1|1|0.000|-1547.874| |NH₄⁺|2|2|1|0.144|-3095.580| |NH₄⁺|3|3|1|0.360|-4643.156| |NH₄⁺|4|4|1|0.624|-6190.658| |NH₄⁺|5|5|1|0.924|-7738.229| |NH₄⁺|6|6|1|1.253|-9285.853| |NH₄⁺|8|8|1|1.979|-12380.835| |NH₄⁺|10|10|1|2.778|-15475.577| |NH₄⁺|12|12|1|3.634|-18570.326| |NH₄⁺|16|16|1|5.485|-24759.517| |NH₄|1|0|2|0.000|-1552.136| |NH₄|2|0|3|0.000|-3105.481| |NH₄|3|0|4|0.000|-4657.118| |NH₄|4|0|5|0.000|-6209.479| |NH₄|5|0|6|0.000|-7762.029| |NH₄|6|0|7|0.000|-9313.970| |NH₄|8|0|9|0.000|-12419.343| |NH₄|10|0|11|0.000|-15525.100| |NH₄|12|0|13|0.000|-18628.100| |NH₄|16|0|17|0.000|-24837.616|
UMA の結果（数値表）
|species|分子数|電荷|スピン|クーロン成分|総エネルギー| |:----|:----|:----|:----|:----|:----| |NH₃|1|0|1|0.000|-1538.811| |NH₃|2|0|1|0.000|-3077.623| |NH₃|3|0|1|0.000|-4616.434| |NH₃|4|0|1|0.000|-6155.245| |NH₃|5|0|1|0.000|-7694.056| |NH₃|6|0|1|0.000|-9232.868| |NH₃|8|0|1|0.000|-12310.498| |NH₃|10|0|1|0.000|-15388.132| |NH₃|12|0|1|0.000|-18465.766| |NH₃|16|0|1|0.000|-24621.032| |NH₄⁺|1|1|1|0.000|-1547.872| |NH₄⁺|2|2|1|0.144|-3094.409| |NH₄⁺|3|3|1|0.360|-4640.358| |NH₄⁺|4|4|1|0.624|-6186.138| |NH₄⁺|5|5|1|0.924|-7728.091| |NH₄⁺|6|6|1|1.253|-9267.410| |NH₄⁺|8|8|1|1.979|-12342.177| |NH₄⁺|10|10|1|2.778|-15409.942| |NH₄⁺|12|12|1|3.634|-18594.755| |NH₄⁺|16|16|1|5.485|-24825.354| |NH₄|1|0|2|0.000|-1552.132| |NH₄|2|0|3|0.000|-3100.393| |NH₄|3|0|4|0.000|-4654.551| |NH₄|4|0|5|0.000|-6207.085| |NH₄|5|0|6|0.000|-7759.962| |NH₄|6|0|7|0.000|-9312.245| |NH₄|8|0|9|0.000|-12415.663| |NH₄|10|0|11|0.000|-15517.568| |NH₄|12|0|13|0.000|-18624.691| |NH₄|16|0|17|0.000|-24841.737|
以下に、横軸に分子数、1分子当たりに規格化したエネルギー誤差を縦軸に取ったプロットを示す。中性NH₃では両モデルともほぼ正常にエネルギーを求められている。一方で、ラジカル種やカチオン種では「大きさについての示量性（size extensivity）」が破綻する挙動を見せた。特に、UMAでは誤差が顕著である。 ![energy_per_fragment_deviation.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697094/d258d17d-b847-45c0-8226-7c7118bee3eb.png) NH₄⁺については、形式電荷によるフラグメント間クーロン反発項を差し引いたプロット（下図）を見ると、OrbMol-v2の偏差はほぼゼロ近傍まで低下していることが分かる。つまり、NH₄⁺については$$E_N \simeq N E_1+E_{\text {Coul }}^{\text {inter }}$$の関係が成立しており、OrbMol-v2はクーロン項を除けばほぼ加算的なエネルギーを与えている（＝1分子あたりの誤差がほぼゼロに張り付いている）。 ![nh4plus_coulomb_corrected_deviation.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697094/8c41c7b3-870c-4cac-9ca7-909dcb0f69b9.png) 一方で、ラジカル種に関してはOrbMol-v2でも示量性が崩れていることが分かった。 - OrbMol-v2：概ね −0.2～−0.6 eV/分子 - UMA-s-1p2：およそ +1.9～−0.5 eV/分子 OrbMolはUMAよりスピン方向の示量性が良い傾向を示すが、NH₄ラジカルの最大スピン系列に対してサイズ無矛盾ではない。したがって、複数のラジカル種が遠距離に存在する系ではエネルギーの定量性が損なわれている可能性に注意する必要がある。 | 系 | OrbMol-v2 | UMA-s-1p2 | 評価 | | --------- | ---------------- | -------------- | ---------------------- | | NH₃ | ほぼ示量的 | ほぼ示量的 | 両モデルで妥当 | | 中性NH₄ラジカル | 明確な非示量性 | 非常に大きく非単調な非示量性 | サイズに対する示量性が崩れており要改善 | | NH₄⁺ | クーロン項の補正後はほぼ示量的 | 1分子あたり数eV規模で非示量的 | OrbMolは良好、UMAは適用範囲外の挙動 | OrbMol-v2は中性閉殻系に加えて、等価な複数荷電フラグメントに対しても、クーロン相互作用を分離すれば比較的良好な示量性を示すことが分かった。一方、複数ラジカルの高スピン状態では、両モデルとも明確な非示量性が残る。UMA-s-1p2は中性閉殻系では示量的であるが、総電荷または総スピンを系サイズとともに増加させる条件では、エネルギーが大きく非加算的となる。 > 上記の計算では、あえて各多分子構造を最適化せず、共通の固定構造に対する1点計算を行っている。これは、 > - OrbMolとUMAで異なる最適化構造になること > - 構造変化とサイズ依存性の混同 > - NH₄⁺間の長距離反発による分子全体の並進 > > を避け、純粋にモデルの原子数・総電荷・スピン条件への応答を比較するためである。 :::note info ### モデルの破綻の不連続性について UMA-s-1p2では、NH₄⁺の分子数がN=10までの間は、1分子あたりのエネルギーはある程度線形に変化しているように見えるが、N=12以降になると急変して符号まで逆転してしまった。このOrbMol-v2との挙動の差異は、両者のアーキテクチャ設計の違いに由来すると推測される。 UMAでは各原子のエネルギーについて - その原子の周りの局所的な環境（近くにどんな原子がどう並んでいるか） - 系全体の総電荷Qという「グローバルな条件」 の両方を入力情報として計算する。このグローバルな条件（総電荷Q）は、"embedding table"という一種の対応表を使って、整数の電荷を数値ベクトルに変換し、それを各原子のエネルギー計算に混ぜ込む形で使われている。 普通の物理量として考えれば「電荷8」と「電荷10」や、「電荷10」と「電荷12」は数値としても近く、影響も連続的となるはずだが、"embedding table" から得られるベクトルが連続的である保証は無く、特に学習データの少ない多電荷の系では非線形かつ予測不能な挙動を招きやすいと考えられる。 プロットに現れているUMAポテンシャルの破綻は、元々のアーキテクチャ設計と、学習データの疎な領域への外挿の掛け合わせによって顕在化したアーティファクトである可能性が高い。 > そもそもDFTでも片手で収まらないほどの多電荷系を計算することはまず無いので、実運用上はあまりお目にかからないタイプのエラーではある。 ::: ## 原子間距離と confidence の関係 OrbMol-v2 の特長として、勾配の予測誤差を推論できる点が挙げられる。これは紛らわしいことに "confidence" と命名されているが、「信頼度」というよりは「勾配予測値の誤差の区間ごとの確率分布」であり、この値が大きいほど**モデルの予測値がDFT計算の結果から外れている可能性が高い**ことを意味する。各原子の勾配にどの程度の誤差が見込まれるかを示すモデル内部の推定値、と解釈するのが適切である。 ### NH₄⁺からのプロトン解離経路 先ほどのNH₄⁺の脱プロトン経路について、"confidence" に相当する量をプロットすると、以下のようになる。ここで、横軸はエチレンのC=Cの中点とZr原子との距離、縦軸は相対エネルギー（eV）および勾配（eV/Å）である。OrbMol-v2のプロットには「OrbMolが見積もった勾配誤差の平均の推定値」をカラーマップで示した。 ![comparison_relative_energy_with_confidence.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697094/77b6d88c-96c8-4c6d-8b45-e14d336904b1.png) ![comparison_gradient_with_confidence.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697094/e3149e2e-df59-401c-986b-e3396621b73c.png) 勾配の予測誤差が大きくなっている区間で、OrbMol-v2とDFTのプロットが乖離していることが観察される。この区間においてDFTよりも勾配が大きめになったために、OrbMol-v2では解離エネルギーが過大に評価されている。（勾配誤差が積分されて蓄積するため） ### 錯体からのエチレン解離経路 エチレンが配位したジルコノセンに対して、エチレンの解離経路上の12点で、OrbMol-v2、UMA-s-1.2、DFT（ORCA；ωB97M-V/def2-TZVPD）でそれぞれ勾配計算を実施した。[^Zr] > ここでの解離経路はエチレンを剛体として平行移動させただけの構造の列であり、各フラグメントのエネルギーは最小化されていないことには注意。 > 参考までに、Zr-(C=C)間距離が 4.2 Åの構造は下図のようになっている。 > ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697094/14730d69-8383-4624-aef2-d908b336d296.png) 以下に、各計算レベルでのエネルギー・勾配のプロファイルを示す。横軸はエチレンのC=Cの中点とZr原子との距離、縦軸は相対エネルギー（eV）および勾配（eV/Å）である。OrbMol-v2のプロットには「OrbMolが見積もった勾配誤差の平均の推定値」をカラーマップで示した。 ![comparison_relative_energy_with_confidence.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697094/c5baf615-f955-4b7c-901d-8d2938ca8641.png) ![comparison_gradient_with_confidence.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/697094/857e3eca-5a93-429c-a43f-d9e6a3751e05.png) プロットを見る限り、エチレンの解離エネルギーは概ね次のようになっている。 | 手法 | 解離側の相対エネルギー | | --------- | ------------: | | ORCA-DFT | 約1.10–1.11 eV | | UMA-s-1p2 | 約1.07 eV | | OrbMol-v2 | 約1.21–1.22 eV | ORCAに対する誤差はUMA-s-1.2の方が小さく、ORCAのポテンシャル曲線とほぼ重なっている。一方、OrbMol-v2は3 Å付近からORCAより高い相対エネルギーを与え始め、4–7 Åでは約0.1–0.15 eV高くなっている。この結果からは、OrbMolはエチレンの配位構造をやや過剰に安定化しており、解離エネルギーを過大評価しているようである。 ポテンシャル曲線が類似していることを考えると、勾配もOrbMolよりUMAの方がORCAに近くなる。実際、3.1–3.5 Å付近の勾配極大では、 - ORCA：約0.69–0.70 eV Å$^{−1}$ - UMA：約0.69–0.70 eV Å$^{−1}$ - OrbMol：約0.77–0.78 eV Å$^{−1}$ 5.5 Å付近では、 - ORCA：約0.08 eV Å$^{−1}$ - UMA：約0.10 eV Å$^{−1}$ - OrbMol：約0.17 eV Å$^{−1}$ となっており、OrbMolは中距離領域でも引力の大きさを過大評価していた。これは、金属錯体の系ではUMAポテンシャルの方が記述力が高い可能性を示唆する結果である。 見た目の印象としては、エネルギープロファイルの曲率がDFTのそれと乖離する区間で "confidence" の値が増加しているように見える。例えば、勾配予測値の推定誤差が 0.03 eV/Å を上回るような点を記録しておき、DFTとの誤差を検算する、という使い方はできそうである。 今回、「"confidence" の値」と「DFTとのエネルギーの乖離」が嚙み合わないケースは見つけられなかった。今後このモデルが様々な系で利用される中で、"confidence" が高い（＝予測値としては不安定）にもかかわらずエネルギーの再現度が高いケースや、反対に "confidence" が低い（＝予測値としては安定的）にもかかわらずエネルギーの再現度が低いケースも出てくる可能性がある。 # 総括 OrbMol-v2 の "latent charge" は、何らかの結合組み換えの際に、モデル内部で電荷がどのように割り当てられるかを可視化する指標の一つとして扱える場面はありそうに思われる。ただし、これはあくまでモデル内部で学習された数値であり、量子化学的に評価可能な部分電荷そのものではない。実際にNH₄⁺の解離を調べると、latent charge が示唆する解離後の電荷状態と、エネルギーの予測値から読み取れる解離極限との間に食い違いが見られた。 また、"confidence" から予測される勾配誤差は、モデルにとって予測が難しい構造を特定する手がかりにはなるが、エネルギーの正確性や、予測された電子状態の物理的な妥当性まで保証するものではない。 latent charge と confidence は、「モデルが内部でどのように振る舞っているかを診断するための手がかり」として利用し、量子化学計算の結果や物理・化学的な知見と照らし合わせながら解釈する必要がある。 # 参考資料 ## ジルコノセンにエチレンが配位した構造 ```:Cp2ZrMe_ethylene_crd.xyz 31 Cp2ZrMe_ethylene_crd Zr 0.779296294862450 0.461540020244534 -0.060089391426586 C 1.284663373296015 -0.231318371990548 2.291268217388600 C 0.461118230301621 0.904406319220694 2.387081887252628 C -0.819739751616609 0.568977757003585 1.870521102037246 C -0.779593267224746 -0.772763180177518 1.451812450116817 C 0.530236026777381 -1.264436770792146 1.683161329965894 H 2.309636851491787 -0.299811235936479 2.615356101437578 H 0.745819452733803 1.851389545102458 2.819061926237970 H -1.684236683060149 1.215021289125552 1.834609743240058 H -1.602171690033230 -1.331352301061905 1.033160070990444 H 0.878314482067817 -2.266819438188470 1.483361454190772 C 1.155792924770431 0.439515322301455 -2.553294909645758 C -0.180131079518398 0.856549011065022 -2.348736458813501 C -0.891052643936091 -0.212536336953277 -1.762024271772002 C 0.009070753384614 -1.297785307245149 -1.611801174966658 C 1.269112366049425 -0.893584233676866 -2.107027255330991 H 1.945722264367253 1.027109960491772 -2.994643360109799 H -0.594376986389860 1.815456171995227 -2.623270944106561 H -1.939445528935323 -0.216570809048458 -1.505580970088253 H -0.240325753194718 -2.273362487062104 -1.226632926770974 H 2.162017138325702 -1.496316838936970 -2.139295155717484 C 3.027504497259511 0.467657115176977 0.008251336032641 H 3.400747725677556 -0.548900813467845 0.172241273908601 H 3.480046315567146 0.823290377659096 -0.920273288952848 H 3.414419805393677 1.074500619941054 0.832289958041418 C 1.402880878829465 3.140790642870058 -0.256459203045283 C 0.112931474344231 3.188850181451401 0.095969677643536 H 2.195872468100097 3.238897311811067 0.474696576366770 H 1.700027907037965 3.072393118493063 -1.296413346780839 H -0.184267910981800 3.352803389896210 1.123176285395068 H -0.677214735551277 3.166679156988618 -0.646728586328544 ``` ## ORCA用にOpenMPI 4.1.8を別途導入する場合 ORCAのバージョン6.1.1はOpenMPI 4.1.8に動的にリンクするため、並列計算を実行するにはOpenMPI 4.1.8を導入する必要がある。 以下にWSL環境に導入する手順の例を示す。基本的にはこれらを上から順に実行すればよい。筆者の環境では5分程度で完了した。 ```sh:ビルドに必要なものを導入する sudo apt update sudo apt install -y \ build-essential \ wget \ tar \ perl \ gfortran \ libhwloc-dev ``` ```sh:OpenMPI 4.1.8をダウンロードする mkdir -p ~/src ~/opt cd ~/src wget https://download.open-mpi.org/release/open-mpi/v4.1/openmpi-4.1.8.tar.gz tar -xzf openmpi-4.1.8.tar.gz cd openmpi-4.1.8 ``` ```sh:ORCA専用の場所へインストールする ./configure \ --prefix="$HOME/opt/openmpi-4.1.8" \ --disable-mpi-fortran ``` ```sh:ビルドする make -j"$(nproc)" make install ``` ```sh:.bashrcに設定する cat >> ~/.bashrc <<'EOF' # OpenMPI 4.1.8 for ORCA export OPENMPI_ROOT="$HOME/opt/openmpi-4.1.8" export PATH="$OPENMPI_ROOT/bin:$PATH" export LD_LIBRARY_PATH="$OPENMPI_ROOT/lib:${LD_LIBRARY_PATH:-}" EOF ``` ```sh:反映する source ~/.bashrc ``` ```sh:確認する mpirun --version which mpirun which orca ``` 以下のように出力されればOK。 ``` mpirun (Open MPI) 4.1.8 Report bugs to http://www.open-mpi.org/community/help/ /home/hogehoge/opt/openmpi-4.1.8/bin/mpirun /home/hogehoge/orca_6_1_1_avx2/orca ``` ```sh:共有ライブラリの不足を確認する ldd "$ORCA_ROOT/orca" | grep "not found" ``` 何も表示されなければOK。 [^Zr_crd]: ジルコノセンにエチレンが配位した構造もここに追加している。xyz構造は補足資料を参照のこと。 [^C=C]: エチレン配位構造におけるC=C結合長は、解離した分子状態における平衡位置の1.3248 Åから 1.3382 Åへ約0.013 Å伸長していた。この伸長量は比較的小さいため、強いメタラサイクル型の活性化ではなく、主として通常の $\pi$ 配位に対応する構造変化と見なせる。この辺りはOrbMolの学習データであるOMolデータセットに多数の金属錯体のデータセットが含まれているため、良い精度でDFTレベルの計算を再現できているものと思われる。 [^cutoff]: 一般に、原子間距離がカットオフ半径（多くの場合 6 Å）を超えると分子グラフが連結でなくなるため、グラフニューラルネットワーク内のmessage passingが断絶して相互作用を表現できなくなる。このアーティファクトにより、サイズに対する無矛盾性・示量性を満たせなくなることがある。特に、電荷を有する系や、電気的に中性であっても離れた部位で電荷分離している系ではこの問題が顕著となる。 [^latent]: これはOrbMolの開発チームが「電荷ラベルなしで、エネルギー・力だけから創発した潜在的な特徴量」と述べていることとも対応する。厳密な電荷として取り扱うのではなく、あくまでもポテンシャルエネルギーを上手く再現するために現れてくる量、と見なす立場である。 [^Zr]: これは全くの余談だが、ORCAでのエネルギー計算をCPU 8コアを利用する設定で走らせたところ、1点あたり約12分を要した。CPUはAMD Ryzen 7 5800X3Dを使用。計算中の温度は91-92℃に張り付いていた。