1/r の発散の防止について

方針

　ポテンシャル等 $f\left({\frac {1}{r}}\right)$　で定義される関数は、$r\ll 1$　において発散を生じる。 これを防ぐため、$f\left({\frac {1}{r}}\right)\to f\left({\frac {1}{r^{mod}}}\right)$ へ置き換えることを考える。$\xi ={\frac {1}{r^{mod}}}$ とし、$f\left({\frac {1}{r^{mod}}}\right)$ の$r$ 微分を考えると、

${\begin{aligned}{\frac {df}{dr}}&={\frac {df}{d\xi }}{\frac {d\xi }{dr}}\\{\frac {\partial ^{2}f}{\partial r^{2}}}&={\frac {d^{2}f}{d\xi ^{2}}}\left({\frac {d\xi }{dr}}\right)^{2}+{\frac {df}{d\xi }}{\frac {d^{2}\xi }{dr^{2}}}\\{\frac {\partial ^{3}f}{\partial r^{3}}}&={\frac {d^{3}f}{d\xi ^{3}}}\left({\frac {d\xi }{dr}}\right)^{3}+3{\frac {d^{2}f}{d\xi ^{2}}}{\frac {d\xi }{dr}}{\frac {d^{2}\xi }{dr^{2}}}+{\frac {df}{d\xi }}{\frac {d^{3}\xi }{dr^{3}}}\\\end{aligned}}$

となっている。 つまり、例えば、LJポテンシャルなら、

${\begin{aligned}U_{LJ}(\xi )&=4\varepsilon (\sigma ^{12}\xi ^{12}-\sigma ^{6}\xi ^{6})\\\mathbf {F} _{LJ}(\mathbf {r} )&=-4\varepsilon (12\sigma ^{12}\xi ^{11}-6\sigma ^{6}\xi ^{5}){\frac {d\xi }{dr}}\mathbf {e} \end{aligned}}$

相互作用テンソルなら、

${\begin{aligned}T^{(0)}&={\frac {\xi }{4\pi \varepsilon _{0}}}\\\mathbf {T} ^{(1)}&={\frac {1}{4\pi \varepsilon _{0}}}{\frac {d\xi }{dr}}\mathbf {e} \\\mathbf {T} ^{(2)}&={\frac {1}{4\pi \varepsilon _{0}}}\left({\frac {d^{2}\xi }{dr^{2}}}\mathbf {ee} +{\frac {d\xi }{dr}}\mathbf {e} ^{(2)}\right)\\T_{ijk}^{(3)}&={\frac {1}{4\pi \varepsilon _{0}}}\left({\frac {d^{3}\xi }{dr^{3}}}e_{i}e_{j}e_{k}+{\frac {\partial ^{2}\xi }{\partial r^{2}}}(e_{ij}^{(2)}e_{k}+e_{ki}^{(2)}e_{j}+e_{jk}^{(2)}e_{i})+{\frac {\partial \xi }{\partial r}}e_{ijk}^{(3)}\right)\end{aligned}}$

である。 従って、$\xi ={\frac {1}{r^{mod}}}$　及びその微分をあらかじめ計算しておくことで、プログラムの変更を最小限にして${\frac {1}{r}}$　の発散を阻止することができる。 また、長距離部分のカットオフも同様の方法で組み込むことができる。

ξの形状

　微分の計算のしやすさから、

$\xi (r)={\frac {1}{r^{mod}}}={\begin{cases}\sum _{n=0}^{8}a_{n}r^{n}&r<r_{1}\\{\frac {1}{r}}&r_{1}<r<r_{2}\\\sum _{n=0}^{8}b_{n}(r-r_{3})^{n}&r_{2}<r<r_{3}\\b_{0}&r_{3}<r\end{cases}}$

と置くのが良い。 この関数形では、$r>r_{3}$ 　において ξ が0でないが、エネルギーの基準点を

$U^{mod}(\xi )=U(\xi )-U(b_{0})$

と補正することにより、トラジェクトリー計算には影響がない。（多項式の次数を上げることで$b_{0}=0$　とすることも可能だが、その場合、ξ の高次の微分でスイッチング関数の影響が無視できなくなる。）

　$a_{n},b_{n}$　は、4回微分が連続という条件を課すことで、

${\begin{aligned}a_{n}&=c_{n}(0,r_{1})\\b_{n}&=c_{n}(r_{3},r_{2})\\c_{0}(r_{1},r_{2})&=-{\frac {-70r_{2}^{4}-35r_{2}^{3}r_{d}-15r_{2}^{2}r_{d}^{2}-5r_{2}r_{d}^{3}-r_{d}^{4}}{70r_{2}^{5}}}\\c_{1}(r_{1},r_{2})&=0\\c_{2}(r_{1},r_{2})&=0\\c_{3}(r_{1},r_{2})&=0\\c_{4}(r_{1},r_{2})&=0\\c_{5}(r_{1},r_{2})&=-{\frac {35r_{2}^{3}+30r_{2}^{2}r_{d}+15r_{2}r_{d}^{2}+4r_{d}^{3}}{5r_{2}^{5}r_{d}^{4}}}\\c_{6}(r_{1},r_{2})&=-{\frac {-14r_{2}^{3}-13r_{2}^{2}r_{d}-7r_{2}r_{d}^{2}-2r_{d}^{3}}{r_{2}^{5}r_{d}^{5}}}\\c_{7}(r_{1},r_{2})&=-{\frac {70r_{2}^{3}+68r_{2}^{2}r_{d}+39r_{2}r_{d}^{2}+12r_{d}^{3}}{7r_{2}^{5}r_{d}^{6}}}\\c_{8}(r_{1},r_{2})&=-{\frac {-5r_{2}^{3}-5r_{2}^{2}r_{d}-3r_{2}r_{d}^{2}-r_{d}^{3}}{2r_{2}^{5}r_{d}^{7}}}\\\end{aligned}}$

となる。 なお、$r_{d}=r_{2}-r_{1}$ である。（導出はファイル:Switch.nb 参照。 関数の形状や微分形もこちらを参照。 3階微分はやや波うちが大きいが、$r_{d}$　を大きく取ることで改善される。）

switch.nb

上の $1/r$ および $1/r^{mod}$ の形をプロットする mathematica スクリプトである。スウィッチの境界は $r_{c}-r_{s}=5.5$ A, $r_{c}=11.5$ A としている。実行したプロット結果は、Kikkawa et al. J. Am. Chem. Soc. 2015, 137, 8022-8025 の Supporting Information 中の Figure S3 にも記載されている。