Estimating Irving-Kirkwood pressure profiles from per-atom virial and density profiles: a closure for the contour-dependent term in planar liquid films

技術要約：原子ごとのビリアルおよび密度プロファイルからのアーヴィング–カークウッド圧力プロファイルの推定

問題提起
不均一な流体の局所的な圧力テンソルは一意的ではなく、原子間に働く力が空間内にどのように分布するかという「輪郭（contour）」の選択に依存する。平面分子動力学（MD）シミュレーションでは、主に2つの接線方向圧力プロファイルが用いられる：

1. アーヴィング–カークウッド（IK）プロファイル： 標準的な選択肢であり、対ビリアルを2原子を結ぶ直線に沿って一様に分布させるものである。これの計算には、すべてのスラブに対して全原子対に対する専用の計算コストの高いパスが必要であり、著者のテストでは、シミュレーションの実行時間の大部分（例：55%）を占めることもある。
2. 原子ごとのビリアル（PAV）プロファイル： 標準的なMDコード（例：LAMMPS）によって、低コストで計算される（原子ごとのビリアルをビン分割する）。接線成分については、これはハラシマ（Harasima）の輪郭と一致する。

表面張力（$\gamma$）のような積分量は輪郭の選択によらず不変であるが、局所的なプロファイルは、対量 $Q(z)$ の微分によって異なる。すなわち $P_{IK}(z) - P_{PAV}(z) = dQ/dz$ である。この差は、非等価な界面（例：非対称な吸着）を持つ薄膜において重要となる。なぜなら、輪郭の選択によって、両方の面間の張力の分配が変わるからである。本論文は、高コストなIK計算を行うことなく、安価なPAVプロファイルから正確なIKプロファイルを復元する必要性に対処している。

手法
著者は、密度プロファイル $\rho(z)$ と共存液相のバルク径分布関数 $g(r)$ のみを用いて、ゲージ項 $Q(z)$ およびその微分 $dQ/dz$ を予測する、閉鎖近似に基づく推定法を提案している。

1. 理論的定式化：

   • IKとPAVの重みの差は、ランプ関数 $q(z)$ の全微分として表される。
   • 量 $Q(z)$ は、これらのランプ関数に対するビリアル重み付きの和として定義される。
   • 閉鎖近似 $\rho^{(2)}(\mathbf{r}_i, \mathbf{r}_j) \approx \rho(z_i)\rho(z_j)g(r_{ij})$ の下で、$Q(z)$ は、対ポテンシャルの微分 $u'(r)$ と $g(r)$ から構成されるカーネル $K(\zeta)$ に対する密度プロファイルの二重積分へと簡約される。

2. 改良：

   • 対局所補間（Pair-Local Interpolation）： 界面において環境が変化することを考慮するため、カーネル $K$ は「液体」カーネル（バルク $g(r)$ を使用）と「希薄」カーネル（低密度極限 $g(r) = \exp[-u(r)/k_BT]$ を使用）の間で補間される。補間ウェイトは、対の平均密度に依存する。
   • 混合系への拡張： 多成分系の場合、閉鎖は各成分のペアに対して適用される。バルクの参照値が存在しないクロスペア（例：希薄な吸着種）については、主要成分の $g(r)$ をペアの直径でスケーリングする共形近似が用いられる。
   • 垂直成分のテスト： 同じ閉鎖を垂直成分（$zz$）にも適用する。機械的平衡においてIK垂直圧力が一定（蒸気圧に等しい）である必要があるため、補正項$dQ_1/dz$ は自己整合性テストとして機能する。つまり、IKの実行なしに、補正されたPAVプロファイルが平坦でない場合、その閉鎖誤差を定量化できる。

3. 推定器：
   推定されたIK接線プロファイルは、以下のように計算される：
   $$P_{IK, est}(z) = P_{PAV}(z) + \frac{dQ}{dz}$$
   ここで $P_{PAV}$ は測定されたプロファイルであり、$dQ/dz$は密度プロファイルとバルク$g(r)$ から導出される。

主な結果
本手法は、厚さが7から25$\sigma$の範囲にある20種類のレナード・ジョーンズ自由浮遊膜、3つのカットオフ半径（3.5, 5.0, 7.5$\sigma$）、および2つの温度（0.72および0.90）に対してテストされた。

• 接線プロファイルの精度：

  • 推定されたIKプロファイルと測定されたIKプロファイルの間の平方根平均二乗（rms）残差は 0.001–0.002 $\epsilon/\sigma^3$ であり、補正前のPAVの差（0.010–0.021）と比較して10分の1に減少した。
  • 液体内部において、残差は $T=0.72$ において統計ノイズの2倍以内であり、$T=0.90$ ではノイズレベルと一致した。
  • 界面において、推定器は輪郭の差の 86–93% を捉えている。残りの誤差は、閉鎖の精度が低下する急峻な密度勾配に跨る短い対に起因する、液体の側面に現れる系統的なダブルット（二重構造）である。
  • 中心面異方性（$\Delta P = P_{IK} - P_{PAV}$、機械的安定性基準に不可欠）は、0.03の範囲内で 0.0012 以内の精度で再現された。補正前のPAV値は、薄膜においてこの異方性を22–28%過小評価している。

• 垂直成分：

  • 垂直成分に閉鎖を適用することで、PAV $zz$ プロファイルに見られる偽の構造の 83–92% が除去され、機械的平衡に必要な平坦性が回復した。これは、IKシミュレーションを必要としない検証指標となる。

• 表面張力と分配：

  • 表面張力はすべてのプロファイルにおいて不変であり、予想通りであった。
  • 非等価な界面を持つ膜（表面活性吸着種を用いてシミュレート）について、本手法は各面の張力を解明することに成功した。クリーンな面では、純粋な膜の張力が統計誤差内で回収された。
  • 輪郭の選択によって張力が転送される量を示すゲージ項 $Q(z)$ は、全表面張力の 0.6% 以内の精度で回収された。

意義と主張
本論文は、提案された推定器により、標準的で低コストなPAVシミュレーションから、物理的に動機付けられたアーヴィング–カークウッド接線圧力プロファイルが得られることを主張している。

• 計算効率： 本手法は、高コストでCPU負荷の高いIKのペアパスを排除する。研究対象となったシステムでは、サンプリング段階の壁時間（wall time）を約2.2倍（3.84hから1.73hへ）短縮した。これは、以前はIK計算が実行時間の55%を消費していたためである。
• 物理的意味： 著者は、バルクの内部を持たない薄膜においては、局所的な表面張力の分配は輪郭に依存すると論じている。IKの輪郭は、応力を担う原子の間に応力を配置するため、物理的に意味のある選択として特定されている。この推定器は、追加の計算コストなしに、非対称な膜（例：優先的な吸着がある場合）における局所的な分配を可能にする。
• 検証： 本手法は、直接的なIK測定だけでなく、「平坦性テスト」を通じた垂直成分の検証によっても妥当性が示されている。これはIKデータなしで行える。
• 限界： 著者は、本手法が閉鎖近似に依存しており、界面が広く（高温）、液体が構造化されていない場合に最も高い性能を発揮することを指摘している。界面における残留誤差（約10%の差）は、単純な閉鎖では捉えきれない短距離相関に起因するが、この誤差は機械的安定性に関する結論に影響を与えないほど小さいとみなされている。

要約すると、本論文は、遍在するPAV圧力プロファイルを補正するための実用的かつ低コストな計算ツールを提供し、薄膜や複雑な界面における正確な局所応力解析および機械的安定性基準を可能にするものである。