第 4 章　場の高速化と二軸周期境界

4.1 Direct、treecode、FMM

$N$ 個の source に対する $M$ 点の Direct 評価は $O(MN)$ の panel 和である。遠方近似を含まないので、小規模な計算では他の手法の基準になる。

treecode は source の octree を評価点ごとに走査し、十分に遠い node を一つの monopole で置き換える[2]。 FMM は target 側にも局所展開を作り、遠方群の寄与を多くの評価点で共有する[4]。どちらも近傍は第3章の解析 panel 積分で評価する。

4.2 treecode の遠方判定と電荷相殺

node 内の総電荷、絶対電荷和、電荷中心を

\begin{equation*} Q_n=\sum_{i\in n}q_i,\quad S_n=\sum_{i\in n}|q_i|,\quad \bm{c}_{Q,n}=\frac{\sum_{i\in n}q_i\bm{c}_i}{Q_n}\tag{4.1} \end{equation*}

とする。 $Q_n$ が極小のときは幾何中心を用いる。 node 半径 $R_n$ は全三角形の頂点を覆うようにとり、評価点との距離 $d_n$ に対して

\begin{equation*} R_n<\theta(d_n-R_n)\tag{4.2} \end{equation*}

を遠方の条件とする。重心だけから半径を決めると、評価点に近い大きな panel を誤って monopole にまとめるためである。

さらに、正負の電荷が相殺する node では電荷中心が不安定になるので、 node 内の電荷が同符号（丸め誤差の範囲で $S_n=|Q_n|$ ）のときだけ遠方近似を受理し、異符号の電荷が混在する node では子へ降りる。この規則は相殺時の精度を守る一方、細かく正負が混在する分布では高速化の効果を弱める。

4.3 Cartesian FMM の展開

多重指数 $\alpha=(\alpha_x,\alpha_y,\alpha_z)$ を用い、 $|\alpha|=\alpha_x+\alpha_y+\alpha_z$ 、 $\alpha!=\alpha_x!\alpha_y!\alpha_z!$ とする。中心 $\bm{c}$ の panel moment を、階乗を含む規約で

\begin{equation*} M_\alpha(\bm{c}) =\sum_i\frac{q_i}{A_i\alpha!} \int_{T_i}(\bm{y}-\bm{c})^\alpha\mathop{}\!\mathrm{d}S_{\bm{y}}\tag{4.3} \end{equation*}

と定義する。重心点の monomial ではなく、三角形上の厳密な面積平均を用いる。

FMM の手順は、leaf の moment を作る P2M、親へ集約する M2M、遠方 source を target の局所展開へ移す M2L、子へ伝える L2L、評価点で計算する L2P からなる。近傍は解析 panel の Direct 和を加える。例えば、中心間変位を $\bm{d}=\bm{c}_{\mathrm{child}}-\bm{c}_{\mathrm{parent}}$ とすると

\begin{equation*} M_\beta(\bm{c}_{\mathrm{parent}}) =\sum_{\mathrm{child}}\sum_{\alpha\le\beta} M_\alpha(\bm{c}_{\mathrm{child}}) \frac{\bm{d}^{\beta-\alpha}}{(\beta-\alpha)!}\tag{4.4} \end{equation*}

で集約する。 $G(\bm{r})=1/\lVert\bm{r}\rVert$ 、 $\bm{R}=\bm{c}_t-\bm{c}_s$ に対して M2L は

\begin{equation*} L_\alpha(\bm{c}_t)\mathrel{+}= \sum_\beta(-1)^{|\beta|}M_\beta(\bm{c}_s) \partial^{\alpha+\beta}G(\bm{R})\tag{4.5} \end{equation*}

となる。 $L_\alpha$ は局所電位の微分係数であり、

\begin{equation*} \phi_{\mathrm{far}}(\bm{x})\simeq k_{\mathrm c}\sum_{|\alpha|\le p}L_\alpha \frac{(\bm{x}-\bm{c}_t)^\alpha}{\alpha!},\qquad \bm{E}_{\mathrm{far}}=-\nabla\phi_{\mathrm{far}}\tag{4.6} \end{equation*}

で評価する。

木、相互作用 pair、平行移動の係数は geometry だけで決まるので初期化時に一度だけ作り、バッチごとには電荷に依存する moment と局所展開だけを更新する。展開次数は $p=4$ とし、精度は開き角 $\theta$ と leaf の要素数で調整する。図4.1に、展開の集約と共有、近傍 Direct 和の関係を示す。

FMM の展開係数の集約と共有 — **図 4.1**　FMM の上向き・相互作用・下向きの流れ。矢印は展開係数の変換であり、粒子の移動ではない。木の一枝と一つの遠方 pair を模式的に示す。近傍の解析 panel 和は遠方展開とは別に加える。

4.4 二軸周期場の分割

$x,y$ 方向に周期 $L_x,L_y$ 、 $z$ 方向に非周期の slab を考える。表面 mesh は一つの周期セルの三角形だけを持ち、場の評価と衝突判定で周期画像を扱う。セルの横面積を $A_{\mathrm c}=L_xL_y$ とし、電位と電場の評価作用素を $K$ と書くと、

\begin{equation*} K_{\mathrm s}= \underbrace{K_{\mathrm{shell}}(N) +\mathcal F\bigl[R_{\mathrm{Ewald}}+s\bigr]-K_0^{\mathrm{sym}}-S}_{k\ne0} +K_0^{\mathrm{phys}}\tag{4.7} \end{equation*}

へ分ける。 $K_{\mathrm{shell}}(N)$ はセル自身と、各周期軸に $N$ 層までの画像の直接評価であり、近傍は panel 積分、それ以外は画像の平行移動を含めた M2L で評価する。 $R_{\mathrm{Ewald}}=K_{\mathrm{Ewald}}-K_{\mathrm{shell}}(N)$ はその外側の画像からの滑らかな残差、 $\mathcal F$ は次節の作用素による近似である。

$K_{\mathrm{shell}}(N)$ をセル上で面平均すると、無限面電荷の平均場 $K_0^{\mathrm{sym}}$ （ $k_{\mathrm c}$ を除いた単位で、高さ $h$ の単位電荷あたり $-2\pi|h|/A_{\mathrm c}$ ）と滑らかな残り $s(h)$ の和になる。 $s$ は $1/r$ を画像ブロック上で積分した閉じた式から $K_0^{\mathrm{sym}}$ を引いたもので、 $h=0$ でも解析的である。残差 $R_{\mathrm{Ewald}}$ の面平均はちょうど $-s$ なので、作用素には面平均が 0 の $R_{\mathrm{Ewald}}+s$ だけを近似させ、 $s$ は source ごとの高さの関数として $S(z)=\sum_iq_is_i(z)$ を実行時に引く。これにより $k\ne0$ 部分は面平均を持たず、面平均成分は Gauss 則で閉じる $K_0^{\mathrm{phys}}$ として一度だけ加わる。図4.2は、画像の範囲による分割と、面平均による分割を区別して示す。

周期画像と零モードの分担 — **図 4.2**　有限画像と遠方残差、および非零モードと零モードの分担。左は $N=1$ の画像集合の模式図で、外側の画像は無限に続く。右の省略した遠方作用素の内訳は式(4.7)に示す。 $z$ は非周期であり、面平均場を一度だけ加えることが重要である。

4.5 Ewald 和から生成する遠方作用素

参照の 2 軸周期 Coulomb 和は、kernel を

\begin{equation*} \frac1r=\frac{\operatorname{erfc}(a r)}r +\frac{\operatorname{erf}(a r)}r\tag{4.8} \end{equation*}

に分け、第 1 項を実空間の画像和、第 2 項を周期 2 軸の逆空間和として評価し、非周期軸は実座標のまま残す[7]。 $a$ は実空間と逆空間の計算量を配分する数値 parameter であり、Debye 遮蔽長の逆数ではない。参照は十分な層数で打ち切った有限和とし、作用素の生成時にだけ評価する。

geometry を固定すれば遠方残差は電荷に対して線形である。 source 木の深さ 2 の節点 $s$ の多重極係数を $\bm{M}_s$ 、target 木の深さ 2 の節点 $t$ の遠方局所係数を $\bm{L}_t^{\mathrm{far}}$ とすると、

\begin{equation*} \bm{L}_t^{\mathrm{far}}=\sum_s\mathsf A_{t,s}\bm{M}_s\tag{4.9} \end{equation*}

と書ける。source 全体を根の多重極一つで表すと、source の高さ方向の広がりが、除いた最も近い画像までの距離（ $N=1$ で周期長の約 2 倍）に近い場合に展開がほとんど収束しないので、深さ 2 の節点に分けて展開の収束を確保する。 $\mathsf A_{t,s}$ は、source 節点を囲む代理点の電荷から多重極への写像（最小ノルムの擬似逆）と、 target 節点の周囲に置いた検査点で残差の電場を局所展開で再現する正則化 QR の最小二乗解を合成して求める。擬似逆を求める前に、長さの次数が異なる多重極の行を正規化する。電場の fit では決まらない局所展開の定数項は、同じ検査点の電位の残差から別に求める。

作用素は geometry だけで決まるので一度生成すれば再利用できる。バッチごとの計算では、通常の M2L の後、L2L の前に $\bm{L}_t\leftarrow\bm{L}_t+\sum_s\mathsf A_{t,s}\bm{M}_s$ を加え、その後は通常の L2L と L2P で評価する。粒子の評価ごとに Ewald 和を計算することはない。

大きな有限画像場と零モードを差し引いて小さな非零場を得る条件では、作用素の fit 残差が誤差を支配し得る。弱い外部場と比較する研究では、非零場の絶対誤差を Fourier 和などの独立な参照で確認する。

4.6 Gauss 則で閉じる物理的零モード

三角形 $i$ のうち高さ $z$ 以下の面積割合を $F_i(z)$ とおくと、累積電荷は

\begin{equation*} C(z)=\sum_iq_iF_i(z)\tag{4.10} \end{equation*}

である。傾いた三角形の $F_i$ は頂点高さの間で区分二次、水平三角形は sheet の跳びになる。 $F_i$ は geometry だけで決まるので区分多項式の係数を前計算し、電荷の更新ごとに $C(z)$ とその積分を組み立てる。下側遠方電場 $E_{\mathrm b}$ と gauge $(z_g,\phi_g)$ を指定すると

\begin{align*} E_{0,z}(z)&=E_{\mathrm b}+\frac{C(z)}{\epsilon_0A_{\mathrm c}},\tag{4.11}\\ \phi_0(z)&=\phi_g-E_{\mathrm b}(z-z_g) -\frac1{\epsilon_0A_{\mathrm c}}\int_{z_g}^{z}C(\zeta)\mathop{}\!\mathrm{d}\zeta\tag{4.12} \end{align*}

となる。全電荷 $Q=\sum_iq_i$ に対する二つの closure は

\begin{equation*} E_{\mathrm b}= \begin{cases} -Q/(2\epsilon_0A_{\mathrm c}) & \text{対称な上下真空},\\ 0 & \text{下側電束を零に指定} \end{cases}\tag{4.13} \end{equation*}

である。どちらも材料内部の誘電分極や遮蔽を解くモデルではない。

非中性 cell の遠方には一定電場と線形電位が残り得る。また $Q=0$ でも、高さによって正負の電荷が分離すれば $C(z)$ は零とは限らない。零モードを単に消すと、この高さ分布と Gauss 則を失う。外部シースを接続する場合は、gauge を上端面の電位の基準に用い、表面電荷の零モードと外部空間電荷の役割を分けて二重に数えない（第8章）。