Coulomb FMMコア詳細

Coulomb FMMコアは、source geometryから作るplanと、電荷から作るstateを分離し、 多重極展開と近傍Direct和から電場を評価します。ここでは、現行Fortran実装の低水準開発API、データ構造、 展開式、periodic2補正を、処理順に沿って説明します。

通常のFMMを利用するための数式と計算フローはFMM、periodic root operatorの構成と運用は periodic2遠方補正にまとめています。このページはFortran内部配列と実装手順の 詳細を扱います。

APIの詳細はbem_coulomb_fmm_core module pageでも確認できます。

• 低水準開発 API / 境界: src/physics/field_solver/fmm/api/
• 内部共通実装: src/physics/field_solver/fmm/internal/common/
• tree / plan 実装: src/physics/field_solver/fmm/internal/tree/
• state / eval 実装: src/physics/field_solver/fmm/internal/runtime/
• periodic2 実装: src/physics/field_solver/fmm/internal/periodic/

コアの内部APIはsimulatorに依存せず、mesh_typeやsim_configを直接useしません。 BEACHのfield solver adapterが、simulator側のデータを変換してこのコアを呼び出します。

1. 目的

このFMMコアは、固定されたsource geometryと可変電荷src_q(n)から、多数の評価点における Coulomb電場を計算します。低水準の検証・periodic operator構築用にはmonopole位置 src_pos(3,n)を受けるgeneric planがあり、BEACHの表面電荷には三角形3頂点を受けるpanel planを 使います。generic planはBEACHの表面source modelではなく、TOMLやC/Python APIからは選べません。

計算の前提と責務は次のとおりです。

• kernel は 3D Coulomb のみ
• source 幾何と電荷更新を分離する
• free と periodic2 のみを対象にする
• 近傍 direct 和もコア内部に含める
• simulator からは配列 API だけが見えるようにする

2. 低水準開発 API

2.1 Fortran API

コアが提供する主な手続きは次のとおりです。

call build_plan(plan, src_pos, options)             ! generic monopole plan
call build_panel_plan(plan, v0, v1, v2, options)   ! BEACH surface plan
call update_state(plan, state, src_q)             ! optional 4th argument: mpi_context
call eval_points(plan, state, target_pos, e)
call eval_point(plan, state, r, e)

call prepare_panel_sources(plan, v0, v1, v2, options)  ! panel geometry only
call eval_direct_points(plan, src_q, target_pos, e)    ! non-periodic exact direct sum
call complete_plan(plan)                               ! = build_panel_plan from here on

入力・出力の意味は次の通りです。

• src_pos(3,n): 低水準generic planのsource点座標。build_plan後は固定とみなす。
• v0(3,n), v1(3,n), v2(3,n): panel planの三角形頂点。build_panel_plan後は固定とみなす。
• src_q(n): source 点の電荷。update_state ごとに更新できる。
• mpi（省略可）: bem_mpi の mpi_context。与えると M2L の target node を rank で分担し、局所展開を allgather で揃える。 全 rank が同じ plan と src_q で呼ぶこと。結果は省略した場合と bit 単位で同じ。
• target_pos(3,m) または r(3): 評価点。
• e(3,m) または e(3): 電場ベクトル。

呼び出し側では、次の点に注意が必要です。

• コアが返す電場にはk_coulombを掛けていません。BEACHのadapterが最後に掛けます。
• build_plan / build_panel_planは幾何依存処理、update_stateは電荷依存処理です。
• eval_point(s) は plan と state が ready な前提です。
• build_panel_plan は prepare_panel_sources と complete_plan を続けて呼ぶのと同じです。 eval_direct_points / eval_direct_potential_points は前半の panel 幾何と src_q だけを使うため、 exact-direct だけの呼び出し側は木・相互作用・演算子と update_state を省けます。

2.2 C ABI / Python 連携

src/physics/field_solver/bem_field_kernel_c.f90は、Fortran APIをiso_c_bindingの opaque handle APIとして公開します。make build-kernelを実行すると、共有ライブラリ build/libbeach_field_kernel.soが生成されます。

主な C ABI:

beach_kernel_get_abi_version(major, minor)
beach_kernel_get_build_info(buffer, capacity, length)
beach_kernel_create(handle)
beach_kernel_destroy(handle)
beach_kernel_build(handle, vertex0_xyz, vertex1_xyz, vertex2_xyz, options...)
beach_kernel_update_charges(handle, src_q)
beach_kernel_eval_e(handle, target_pos, e)
beach_kernel_eval_phi(handle, target_pos, phi)
beach_kernel_eval_e_direct(handle, target_pos, e)
beach_kernel_eval_phi_direct(handle, target_pos, phi)
beach_kernel_force_on_charges(handle, target_pos, target_q, origin, force, torque)

公開ヘッダは Python package 内の beach/include/beach_field_kernel.h にあります。現行 ABI は 2.1 です。C の呼び出し側は、ほかの関数を使う前に beach_kernel_get_abi_version を呼び、 major が 2 と一致し、library の minor が必要な minor（direct APIを使う場合は 1）以上であることを 確認してください。eval_e_direct / eval_phi_direct は非周期planだけで利用でき、同じsource geometryと chargeの非周期exact-direct値を返します。周期planに対してはinvalid argumentを返します。 非周期の handle は build で panel 幾何だけを作り、update_charges では電荷を保持するだけです。 FMM の木・演算子と展開は、最初の eval_e / eval_phi / force_on_charges で作り、以後は電荷が 変わったときだけ更新します。direct だけを使う handle はこれらを作りません。評価関数も handle を 更新することがあるため、1 つの handle への呼び出しは呼び出し側で直列化してください。座標と vector の配列は values[3 * point_index + component] の順で格納します。status code、 periodic far-correction code、handle の所有権は公開ヘッダに定義しています。

Pythonでは、beach.native.kernel.FieldKernelがこのABIをctypesで呼び出します。 Python wrapper はversion問い合わせを提供するlibraryの互換性を読込時に検査します。 問い合わせsymbolがない旧libraryも拒否し、ABI v2以上を明示したlibraryだけを受理します。 calc_object_forces_kernelは、対象object自身のsource電荷をゼロにして sum(q_i E_not_self(r_i))を評価します。これにより、自己力を除外しながら、 periodic2 + cached_kneq0を含むfield kernelをそのまま利用できます。

Beach.scene() / BeachSceneは、Python側でobjectの剛体移動や回転を一時的に適用し、 変換後の三角形3頂点を同じABIへ渡します。剛体変換の補助処理には既定でNumPyを使い、 任意依存のNumba backendも選べます。電場の評価は、どちらの場合もFortran kernelが担当します。

2.3 BEACH adapterでの使い方

BEACH の field solver adapter は、各三角形の 3 頂点を build_panel_plan に渡します。 src_q(i) は三角形全体の総電荷で、面密度は src_q(i)/area(i) です。

• 初期化時は build_panel_plan の直後に update_state を実行します。
• その後のrefreshでは、mesh geometryが変わらない限り既存のplanを再利用し、src_qを渡して update_stateだけを呼びます。
• plan を再構築するのは、planが未構築の場合、source数が変わった場合、 または要素数0としてplan/stateを破棄した場合だけです。

3. データ構造

3.1 fmm_options_type

主な内部オプション:

• theta: well-separated 判定用パラメータ
• leaf_max: source octree の葉に許す最大 source 数
• order: Cartesian 展開次数（1以上）
• softening: 低水準の汎用monopole planだけが使う内部値。BEACH adapterは常に0を渡す
• use_periodic2: 2 周期軸モードの有効化
• periodic_axes(2), periodic_len(2): 周期軸と周期長
• periodic_image_layers: 近傍画像和の層数 N
• periodic_far_correction: coreが受ける値はauto, none, cached_kneq0。periodic2有効時の autoは互換性のためnoneへ正規化される
• periodic_ewald_alpha, periodic_ewald_layers: cached_kneq0のbuild-time Ewald fitで使う 分解パラメータと打切り深さ
• target_box_min/max: dual-target tree を作るときの box

BEACH の adapter は現状 order = 4 を使いますが、コア自体は1以上の可変次数を受けられます。 order = 0 は電場のfar/local展開を持てないため、plan構築時に拒否します。 periodic2 の auto は none に正規化されます。cached_kneq0 は遠方補正を明示的に有効化します。

3.2 fmm_plan_type

幾何にだけ依存する不変データです。

• 多重指数テーブル alpha, deriv_alpha
• source octree
• optional target tree
• source 葉一覧 source_leaf_nodes
• target 葉一覧 leaf_nodes
• 近傍 list near_start/near_nodes
• M2L pair cache m2l_target_nodes/m2l_source_nodes
• periodic image shift 配列
• M2L 微分表 m2l_deriv
• P2M 基底表 source_p2m_basis
• M2M/L2L の平行移動用圧縮テーブル

3.3 fmm_state_type

電荷に依存して毎回更新されるデータです。

• src_q(n)
• multipole(ncoef, nnode)
• local(ncoef, n_target_nodes)
• multipole_active(nnode)
• local_active(n_target_nodes)

multipoleはsource tree nodeごとの多重極係数、localはtarget tree nodeごとの局所展開係数です。 *_activeはzero-nodeの計算を省略するための0/1 flagです。

4. 数学的定義

4.1 source kernel

BEACH runtimeのsource kernelはP0 triangleに固定されています。低水準FMM coreには build_planによる汎用monopole planも内部API・回帰試験用として残っていますが、 BEACH field solver、公開C ABI、Python APIから選択する経路はありません。

runtimeのq_iを三角形$T_i$の総電荷、$A_i$をその面積とし、 $\sigma_i=q_i/A_i$ を三角形上の一定面密度として扱います。

$$\phi_i(\mathbf{x}) = \frac{q_i}{A_i}\int_{T_i} \frac{1}{\lVert\mathbf{x}-\mathbf{y}\rVert}\,dA_{\mathbf{y}}$$

近傍Direct評価には、辺の対数項と立体角を使う解析的P0 panel kernelを使います。 遠方P2Mには、tree nodeの中心$\mathbf{c}$に対する三角形上のmonomialの厳密な面積平均を使います。

$$M_{i,\alpha} = q_i\left\langle(\mathbf{y}-\mathbf{c})^\alpha\right\rangle_{T_i} = \frac{q_i}{A_i}\int_{T_i} (\mathbf{y}-\mathbf{c})^\alpha\,dA_{\mathbf{y}}$$

area weightingはpanel積分とP2M基底に含まれています。q_iは要素の総電荷なので、 $A_i$を重ねて掛ける必要はありません。以後のM2M/M2L/L2Lには、このpanel momentに対する Coulomb/Laplace展開を使います。近傍評価は解析的 panel kernel で処理します。

内部build_planは、必要な低水準試験でだけ $G_\epsilon(\mathbf r)=1/\sqrt{\lVert\mathbf r\rVert^2+\epsilon^2}$を使います。 これは削除済みの公開source modelを復活させるものではなく、beach.tomlから設定できません。

4.2 多重指数

多重指数 $\alpha = (\alpha_x, \alpha_y, \alpha_z)$ を使います。

$$|\alpha| = \alpha_x + \alpha_y + \alpha_z$$ $$\alpha! = \alpha_x! \, \alpha_y! \, \alpha_z!$$ $$\mathbf{r}^\alpha = r_x^{\alpha_x} r_y^{\alpha_y} r_z^{\alpha_z}$$

4.3 P2M

node centerを$c$とすると、leaf nodeのmultipole係数は

$$M_\alpha(c) = \sum_{j \in \text{leaf}} q_j \frac{(\mathbf{x}_j - \mathbf{c})^\alpha}{\alpha!}$$

で定義します。

4.4 M2M

子ノード中心 $c_{\mathrm{child}}$ の係数を親中心 $c_{\mathrm{parent}}$ に平行移動して集約します。 $\mathbf{d} = c_{\mathrm{child}} - c_{\mathrm{parent}}$ とすると

$$M_\beta(c_{\mathrm{parent}}) = \sum_{\alpha \le \beta} M_\alpha(c_{\mathrm{child}}) \frac{\mathbf{d}^{\beta-\alpha}}{(\beta-\alpha)!}$$

現行実装では $\beta - \alpha$ に対応する index と $\mathbf{d}^{\beta-\alpha} / (\beta-\alpha)!$ を build_plan 時に前計算します。

4.5 M2L

source node 中心 $c_s$、target node 中心 $c_t$ に対して $R = c_t - c_s$ とします。

局所展開係数は

$$L_\alpha(c_t) \mathrel{+}= \sum_\beta (-1)^{|\beta|} M_\beta(c_s) D^{\alpha+\beta} G(R)$$

で更新します。

ここで $D^\gamma$ は multi-index 微分です。 現行実装では、$D^{\alpha+\beta} G(R)$をpairごとに前計算し、 m2l_deriv(:, pair)へ保存します。

4.6 L2L

親中心 $c_{\mathrm{parent}}$ の局所展開を子中心 $c_{\mathrm{child}}$ へ平行移動します。 $\mathbf{d} = c_{\mathrm{child}} - c_{\mathrm{parent}}$ とすると

$$L_\alpha(c_{\mathrm{child}}) \mathrel{+}= \sum_{\gamma \ge \alpha} L_\gamma(c_{\mathrm{parent}}) \frac{\mathbf{d}^{\gamma-\alpha}}{(\gamma-\alpha)!}$$

これも build_plan 時に shift monomial を前計算します。

4.7 L2P

評価点 $x$ が属する target leaf の中心を $c_{\mathrm{leaf}}$ とし、 $\mathbf{dr} = x - c_{\mathrm{leaf}}$ とすると

$$E_k(x) = - \sum_{|\alpha| < p} L_{\alpha + e_k}(c_{\mathrm{leaf}}) \frac{\mathbf{dr}^\alpha}{\alpha!}$$

で電場を評価します。 ここで $e_k$ は軸 $k$ の単位 multi-index です。

5. build_plan のアルゴリズム

build_plan は幾何依存処理だけを行います。

5.1 source tree

source 座標の bounding box を再帰的に 8 分割して octree を作ります。 停止条件は次のどちらかです。

• source 数 <= leaf_max
• bbox が十分に小さく、これ以上分割しても意味がない

5.2 target topology

target 側は 2 通りあります。

• target_box が無効: source tree の葉をそのまま target leaf として使う
• target_box が有効: box 全体を覆う別 target tree を作る

periodic2 では target point を box 内に wrap してから target leaf を探します。

5.3 near/farとM2L pair cache

各target leafについてsource treeを再帰的に走査し、near nodeとfar nodeへ分類します。

well-separated 判定は

$$(r_s + r_t)^2 < \theta_{\mathrm{eff}}^2 \, \lVert\mathbf{d}\rVert^2$$

です。

• $r_s$: source node 半径
• $r_t$: target node 半径
• $\mathbf{d}$: node center 間ベクトル
• $\theta_{\mathrm{eff}} = \theta$ for free と periodic2

periodic2 では $\mathbf{d}$ に minimum-image 補正を入れます。

その後、dual-tree再帰によってM2L pair cacheを作り、target nodeごとの索引配列を準備します。

5.4 build時の前計算

build_plan の最後で、refresh ごとに変わらない量を前計算します。

• source_parent_of
• parent_of
• source_p2m_basis
• m2m_term_count, m2m_alpha_list, m2m_delta_list
• l2l_term_count, l2l_gamma_list, l2l_delta_list
• source_shift_monomial
• target_shift_monomial
• shift_axis1, shift_axis2
• periodic_ewald
• periodic_root_operator
• m2l_deriv

この前計算により、update_stateでは主に電荷に依存する加算だけを実行します。

5.5 擬似コード

build_plan(src_pos, options):
  initialize_basis_tables(order)
  precompute_periodic_source_grid()  # cached_kneq0 の box 固定セルと P2M 基底
  build_source_tree(src_pos)
  precompute_source_p2m_basis()
  build_target_topology(target_box)
  build_interactions()
  precompute_translation_operators()
  precompute_periodic2_ewald_data()
  precompute_periodic_root_operator()
  precompute_m2l_derivatives()

6. update_state のアルゴリズム

update_stateはlegacy実装のrefreshに相当します。source座標は固定され、src_qだけが変わることを 前提とします。

6.1 処理順

update_state(plan, state, src_q):
  ensure_state_capacity()
  copy src_q
  clear active flags
  clear multipole/local only when the tree has no source leaves or no M2L pairs
  P2M on source leaves
  M2M bottom-up
  M2L on cached pairs (this rank's block of target nodes when mpi is given)
  allgather local expansions and active flags across ranks (mpi only)
  periodic far correction on anchors
  L2L top-down
  mark state ready

mpi を与えた場合、target node は組数がほぼ等しい連続ブロックに rank 順で分けられ（m2l_rank_blocks）、 各 node の展開は 1 つの rank が従来と同じループで求める。通常精度の同一ビルドでは計算を分担しても結果は同じで、 M2L の組の評価を rank に分配できる。通信と組数の偏りがあるため、所要時間が rank 数に反比例するとは限らない。 演算子と plan の保存領域は引き続き各 rank が保持する。

6.2 OpenMP並列化

現行実装では、次の処理をOpenMPで並列化しています。

• update_state全体を1つのparallel regionで囲み、その内側でsrc_qのコピーとactive flagの初期化を行う
• P2M: source leaf ごとのループ
• M2M: 同一 depth の node ループ
• M2L: target node ごとのループ
• L2L: 同一 depth の node ループ
• build_plan 時の translation / M2L 微分前計算

各ループは1 nodeを1 threadが担当しやすい構造とし、共有配列の更新をnode単位で独立させています。

6.3 実装上の最適化

update_state では次の無駄を避けています。

• $\beta - \alpha$ の multi-index 差分を毎回計算しない
• 親子 center 差分のべき乗を毎回作り直さない
• P2M の monomial 基底を source ごとに build 時に前計算する
• M2M/L2L の有効な (alpha, delta) 項だけを圧縮して持ち、無効項を走査しない
• M2Lではsource nodeのactive flagを調べ、zero-nodeのpairを早期にskipする
• M2Lではtarget node列を繰り返し更新せず、thread-localなlocal_accへ蓄積してから書き戻す
• panel plan では各葉の近接 panel を像シフト込みの成分配列（near_vx … near_area）として build 時に保存し、 点評価は葉の区間を 64 panel ずつ SIMD ループで評価して、リスト順に合算する（panel_field_sum_translated）
• 近接 panel の成分表には、像と葉による重複を含む近接リストの 1 エントリにつき約 200 byte を各 rank で追加保持する
• P2M で target leaf ではなく source leaf 専用 index を使う

7. eval_point(s) のアルゴリズム

評価時の処理は次の通りです。

eval_point(r):
  if plan is not built or state is not ready:
    return zero vector

  if periodic2:
    wrap r into target box

  leaf = locate_target_leaf(r)
  if leaf not found or leaf is not mapped to a leaf slot:
    use direct sum over all sources
    return

  evaluate local expansion at leaf center
  add near direct interactions
  root local already carries periodic root correction when enabled

7.1 leafの特定

• periodic2 では評価点を target box 内へ wrap してから探索する
• target tree があるときは target tree の葉を使う
• target tree が無いときは source tree の葉を使う
• leaf lookupに失敗した場合、またはleafをtreeのleaf slotへ写像できない場合はDirect fallbackに入る

7.2 近傍Direct和

near listに入ったsource indexについてDirect和を計算します。periodic2では、 [-N, N] x [-N, N]の画像シフトも陽に加算します。fallbackでも同じDirect kernelを使います。

7.3 box外fallback

dual-target treeを使う場合、target boxの外にある評価点にはtarget leafがありません。この場合は、 全sourceのDirect和へfallbackします。cached_kneq0は固定target topologyを前提とするため、 target box外の評価をrejectします。

7.4 root補正の位置

cached_kneq0のroot補正は、update_stateがbox格子のtarget anchorごとのlocal展開（state%anchor_local）として 集計します。anchorは適応木の節点ではなく、点評価では葉のlocal展開に加えて、点を含むanchorのlocal展開を評価します。 通常の評価ではroot補正を再計算せず、stateに保存された展開をそのまま使います。

8. periodic2 と遠方補正

この節はFMM core内部から見た数式とfallbackを記録します。設定の選び方、operator fit、cache lifecycle、 k=0 ownershipを通した説明はperiodic2遠方補正に分離しています。

8.1 periodic2

periodic2では、2軸を周期境界、残りの1軸を開放境界とします。

近傍画像和は

$$i, j \in [-N, N]$$

の有限画像を陽に足します。

M2L でも同じ画像シフト集合を使い、各 pair の derivative を画像和で前計算します。

8.2 periodic2 Ewald（Ewald2P）補正

bem_coulomb_fmm_periodic_ewald.f90は、2周期・1開放のCoulomb fieldに対するEwald型の補正を 実装します。ここでいうexactは、コードが実際に評価する有限和を指し、理論上の無限和そのものではありません。 real-spaceとreciprocal-spaceの打切り深さは、fields.periodic.image_layers = Nと fields.periodic.ewald_layers = Lで決まります。この有限和をbuild-time oracleとして使います。

Ewald2Pはruntime particle kernelではなく、cached_kneq0を作るための build-time teacherです。teacherはproxy monopoleに適用されます。 その結果を、box 固定セルの multipoleからlocal展開へのoperatorとしてfitします。実三角形の近傍評価には解析panel kernel、 遠方source表現にはtriangle-averaged P2Mを引き続き使います。

$\alpha$はreal-spaceとreciprocal-spaceの収束を配分する数値パラメータで、Debye遮蔽を表しません。 この分割をruntimeの場評価へ組み込む流れは、periodic2遠方補正の説明にまとめています。¹

8.2.1 記法

周期軸を a_1, a_2、開放軸を f とします。 周期長、セル面積、画像集合、逆格子集合を次のように置きます。

$$L_1 = \operatorname{periodic\_len}(1),\qquad L_2 = \operatorname{periodic\_len}(2),\qquad A = L_1 L_2$$ $$\mathcal I_N = \{(i,j)\in\mathbb Z^2 \mid |i|,|j|\le N\},\qquad \mathcal K_L = \{(m,n)\in\mathbb Z^2 \mid |m|,|n|\le L,\ (m,n)\neq(0,0)\}$$

画像シフトと逆格子ベクトルは

$$\mathbf L_{ij} = iL_1\,\mathbf e_{a_1} + jL_2\,\mathbf e_{a_2},\qquad \mathbf k_{mn} = 2\pi\left(\frac{m}{L_1}\mathbf e_{a_1} + \frac{n}{L_2}\mathbf e_{a_2}\right)$$

と書けます。ソース位置を $\mathbf s$、評価点を $\mathbf r$ とし、

$$\mathbf R_{ij} = \mathbf r - \mathbf s - \mathbf L_{ij},\qquad R_{ij} = \lVert\mathbf R_{ij}\rVert,\qquad z = (\mathbf r - \mathbf s)\cdot \mathbf e_f$$

を導入します。以下では $\alpha =$ fields.periodic.ewald_alpha とします。

8.2.2 実空間項

内部helperが実装しているscreened Coulomb fieldは

$$\mathbf E_\alpha(\mathbf R) = q\left( \frac{\operatorname{erfc}(\alpha R)}{R^3} +\frac{2\alpha}{\sqrt{\pi}}\frac{e^{-\alpha^2 R^2}}{R^2} \right)\mathbf R$$

です。これはポテンシャル

$$\Phi_\alpha(\mathbf R) = q\,\frac{\operatorname{erfc}(\alpha R)}{R}$$

の勾配に一致します。

低水準Ewald helperが内側画像和から差し引くdirect fieldは

$$\mathbf E_\epsilon(\mathbf R) = q\,\frac{\mathbf R}{(R^2+\epsilon^2)^{3/2}}$$

です。BEACHのP0 panel経路では$\epsilon=0$なので、 $\mathbf E_\epsilon=\mathbf E_0=q\mathbf R/R^3$です。

実装上の real-space 補正は

$$\mathbf E_{\mathrm{real}} = \sum_{(i,j)\in\mathcal I_{N+L}} \mathbf E_\alpha(\mathbf R_{ij}) {}- \sum_{(i,j)\in\mathcal I_N} \mathbf E_\epsilon(\mathbf R_{ij})$$

です。実装ではr2 <= tiny(1.0d0)の項をskipするため、self interactionは含まれません。 add_periodic2_exact_ewald_correction_single_sourceにDirect fallbackの $\sum_{(i,j)\in\mathcal I_N}\mathbf E_\epsilon$を加えると、inner imageのdirect部分が打ち消され、 outer shell側がscreened形式に置き換わります。

8.2.3 逆空間項

add_exact_periodic2_reciprocal_space_correction が使う逆空間項は、$(m,n)\neq(0,0)$ に対して

$$\theta_{mn} = \mathbf k_{mn}\cdot(\mathbf r-\mathbf s),\qquad k_{mn} = \lVert\mathbf k_{mn}\rVert$$ $$G^\pm_{mn}(z) = e^{\pm k_{mn} z}\operatorname{erfc}\!\left(\frac{k_{mn}}{2\alpha}\pm \alpha z\right)$$

を定義すると

$$\mathbf E_{\mathrm{rec}} = q \sum_{(m,n)\in\mathcal K_L} \frac{\pi}{A} \begin{pmatrix} \frac{(\mathbf k_{mn})_{a_1}}{k_{mn}}\sin\theta_{mn}\,\bigl(G^+_{mn}(z)+G^-_{mn}(z)\bigr) \\ \frac{(\mathbf k_{mn})_{a_2}}{k_{mn}}\sin\theta_{mn}\,\bigl(G^+_{mn}(z)+G^-_{mn}(z)\bigr) \\ \cos\theta_{mn}\,\bigl(G^-_{mn}(z)-G^+_{mn}(z)\bigr) \end{pmatrix}$$

です。コードでは term_p, term_m, pair_sum に対応します。 この式は、逆格子の k=0 を除いた高周波成分に対応します。

8.2.4 k=0項

add_exact_periodic2_k0_correction が実装しているゼロモード補正は

$$\mathbf E_0 = q\,\frac{2\pi}{A}\operatorname{erf}(\alpha z)\,\mathbf e_f$$

です。single-source の Ewald teacher では k=0 の電場寄与としてこの形を保持します。

8.2.5 実装される補正

以上をまとめると、add_periodic2_exact_ewald_correction_single_source が 1 粒子分に加える補正は

$$\mathbf E_{\mathrm{corr}} = \mathbf E_{\mathrm{real}} {}+ \mathbf E_{\mathrm{rec}} {}+ \mathbf E_0$$

です。このsingle-source補正をproxy/check点で評価し、cached operatorを生成します。

fields.periodic.ewald_alpha が <= 0 の場合、resolve_periodic2_ewald_alpha は

$$\alpha = \frac{1.2}{(N+1)\min(L_1,L_2)}$$

を自動選択します。min(L_1,L_2)\le 0 なら alpha = 0 としてoperator生成を無効化します。 また内部では kmax = max(1, fields.periodic.ewald_layers) として逆空間の有限和を作ります。

cached_kneq0のcold buildはcheck pointsでEwald residualの電場と電位を評価し、box 固定セルの multipoleから target anchorのlocal展開へのoperatorをfitします。有限画像殻の面平均は閉じた式で分離し、fitには面平均0の残差だけを 渡します。電場fitで決まらない定数potential係数はpotential residualから別にfitし、versioned cacheへ保存します。 無限周期のproduction計算にはcached_kneq0を使います（詳細はperiodic2 遠方補正）。

9. 計算量の見方

展開次数$p$を固定し、interaction listの長さに上限があると仮定すると、計算量の目安は次のとおりです。

• build_plan: $O(n \log n)$ に近い
• update_state: $O(n)$ に近い
• eval_point: $O(\log n + n_{\mathrm{near}} \, n_{\mathrm{img}}^2)$ に近い
• eval_points: 上記の点評価を target ごとに並列実行

実際の定数係数は、次の設定に強く依存します。

• order
• theta
• leaf_max
• periodic_image_layers
• periodic_ewald_layers
• target tree の有無

10. 現行実装の制約

このFMMコアは、次の条件に特化した実装です。

• kernel は Coulomb 固定
• simulator adapter の既定次数は order = 4
• source 座標は build_plan 後に不変とみなす
• 対応境界は free と periodic2
• periodic2 は正確に 2 周期軸が必要
• far correctionはnone（既定）、auto、cached_kneq0。 autoはnoneとして動作し、cached backendはproductionの非零modeに使う
• eval_point(s) の返り値には k_coulomb を含めない

10.1 cached periodic nonzero operator

何を高速化するoperatorか

periodic2では、primary cellの外に同じ電荷分布がx/y方向へ無限に続きます。 通常のFMMは近傍画像[-N,N]^2を高速に評価できますが、その外側にある無限個の画像が作る 滑らかな遠方場は含みません。cached_kneq0は、この遠方の差分だけをFMM local展開へ変換する 線形operatorとして事前計算します。これにより、particleを評価するたびにEwald和を計算する必要がなくなります。

入力	cached operator	出力
現在の電荷から作ったbox 固定セルの multipole	geometry固定の行列を適用	target anchorごとの遠方local展開

cache が保存するのは、box 固定格子の source multipole を遠方 local 展開へ写す行列です。要素を含む線源セルの分だけを作り、 別の物体集合で必要になったセルはその都度追加して同じファイルへ書き足します。 電場値、粒子位置、電荷履歴は保存しません。同じ box と設定なら、電荷更新やメッシュの変更後も同じ行列を再利用できます。

1回のfield評価で何を足すか

1回のfield評価では、次の順に各成分を加算します。

順序	成分	目的
1	primary cell + 有限近傍画像	singular/near fieldを通常のFMMとdirectで評価
2	cached Ewald residual	有限画像の外側にある滑らかな無限周期遠方場を補う
3	cached teacherに含まれた対称k=0を減算	非零モードbackendをk!=0だけにする
4	場の合成時に物理的k=0を加算	symmetric_vacuumまたはe_bottom_zeroを反映

cached_kneq0単体では全電場を構成しません。手順3までを非零mode backendが担当し、手順4を electrostatic_snapshotが担当します。exclude_k0は平均場を除外する設定ではありません。 別のboundary providerが平均場を1回だけ加えるための重複防止規則です。

数式との対応

上の1—3をまとめたruntime非零モードkernelは

$$K_{k\ne0} = K_\mathrm{shell}(N) + R_\mathrm{Ewald}^{\mathrm{full}} - K_0^\mathrm{sym}$$

項	作成時期	役割
$K_\mathrm{shell}(N)$	通常のFMM plan/runtime	primary cellと有限近傍画像
$R_\mathrm{Ewald}^{\mathrm{full}}$	cache cold build	full-periodic Ewald解と有限画像和の滑らかな差
$K_0^\mathrm{sym}$	charge state更新時	cached full-periodic kernelから対称$k=0$を除く

$K_0^\mathrm{sym}$は、source高さから作った区分多項式の累積stateを二分探索して評価します。 1評価あたりの計算量は$O(\log n)$です。 最終的なsurface fieldは

$$K_\mathrm{surface}=K_{k\ne0}+K_0^\mathrm{physical}$$

です。$K_0^\mathrm{physical}$のtriangle-height積分と下側境界条件は periodic2静電場で説明します。

field fitとpotentialの定数modeは単位が異なるため、同じleast-squaresの列には含めません。 potential gaugeは、平均residualから別途固定します。

cache lifecycle

段階	動作
identity作成	geometry、target topology、order、周期長、画像層、generator/build versionからfingerprintを作る
warm read	version、fingerprint、shape、checksumが一致したoperatorだけを受理する
miss/corruption	filesystem lockを取得し、operatorを再生成する
publish	同じdirectoryの.tmpをcloseしてからatomic renameする
checkpoint	operator本体は保存しない。再生成可能なcacheとして扱う

同時に起動したjobは、同じlock fileによってcache生成を直列化します。readerは、書込み途中のfileを cache hitとして受理しません。

MPI/OpenMPによるcold build

単位	担当
cache I/Oとlock	MPI rank 0のみ
target operator slice	MPI rankへ均等分配。rank間差は最大1 target
target内のproxy列	OpenMPで並列評価
正則化QR	targetごとに1回作り、全proxy RHSで再利用
operator集約	MPI_Allreduce(SUM)で全rankへ構築

warm runのfield evaluationとcharge refreshには、all-source Ewald和もoperator再fitもありません。

cold buildとwarm runの違い

経路	Ewald teacher	QR fit	particle evaluation
初回cache miss	実行	実行してoperatorをpublish	build完了後に開始
cache hit	実行しない	実行しない	読み込んだoperatorを使用
batch charge refresh	実行しない	実行しない	新しいmultipoleへ同じoperatorを適用

cold buildでは、再利用可能なoperator行列を初回だけ生成します。batchごとに無限周期場を再構築するわけでは ありません。warm runのhot pathでは、全sourceのEwald和を計算しません。

SysA測定値

測定条件は2026-07-12の旧レゴリス入力、order 4、64 target、280 proxy、840 checkです。 時間はcache公開までではなく、表に明記した範囲の実測値です。

構成	実測時間	範囲
旧root-only、1 rank x 1 thread	31分24秒	cold operator build
QR再利用後、1 rank x 1 thread	約25分45秒	operator公開まで
1 rank x 112 threads	47.0秒	cache prime + batch 1
2 ranks x 112 threads	36.7秒	cache prime + batch 1
4 ranks x 112 threads	31.5秒	cache prime + batch 1
6 ranks x 112 threads	30.3秒	cache prime + batch 1

全並列構成でchecksumが一致し、旧operatorとの差はFrobenius相対値1.73e-15でした。 実測時間は上記fixtureに対する値であり、geometryによって変わります。

運用指針

状況	推奨
専用cache prime	1 rank x 112 threadsをcore効率とqueue footprintの基準にする
production allocationが既にある	既存rankで生成してよい。6 x 112なら同fixtureで約30秒
cold buildだけのために増員	4—6 rankの追加利得は小さい
1 core	operator公開後の粒子batchでメモリ不足となった測定例があり、運用対象外
warm cache	fingerprintとchecksumを確認し、そのまま再利用

11. 実装との対応

主な対応箇所:

• 公開 API / ラッパ: src/physics/field_solver/fmm/api/bem_coulomb_fmm_core.f90, src/physics/field_solver/fmm/api/bem_coulomb_fmm_core_build.f90, src/physics/field_solver/fmm/api/bem_coulomb_fmm_core_state.f90, src/physics/field_solver/fmm/api/bem_coulomb_fmm_core_eval.f90
• shared 型定義: fmm_options_type, fmm_plan_type, fmm_state_type （src/physics/field_solver/fmm/internal/common/bem_coulomb_fmm_types.f90）
• plan 構築: build_plan, build_panel_plan （src/physics/field_solver/fmm/internal/tree/bem_coulomb_fmm_plan_ops.f90）
• charge refresh: update_state, p2m_leaf_moments, m2m_upward_pass, m2l_accumulate, l2l_downward_pass （src/physics/field_solver/fmm/internal/runtime/bem_coulomb_fmm_state_ops.f90）
• 評価: eval_point, eval_points （src/physics/field_solver/fmm/internal/runtime/bem_coulomb_fmm_eval_ops.f90）
• periodic2 補助: has_valid_target_box, use_periodic2_cached_kneq0, use_periodic2_root_operator, build_periodic_shift_values, add_point_charge_images_field, wrap_periodic2_point, apply_periodic2_minimum_image, distance_to_source_bbox, distance_to_source_bbox_periodic （src/physics/field_solver/fmm/internal/periodic/bem_coulomb_fmm_periodic.f90）
• periodic2 Ewald/oracle: resolve_periodic2_ewald_alpha, precompute_periodic2_ewald_data, add_periodic2_exact_ewald_correction_single_source （src/physics/field_solver/fmm/internal/periodic/bem_coulomb_fmm_periodic_ewald.f90）
• periodic2 root operator: precompute_periodic_root_operator, build_missing_cell_operators （src/physics/field_solver/fmm/internal/periodic/bem_coulomb_fmm_periodic_root_ops.f90）。 セル単位 cache の読み書きと lock は lookup_cell_cache, publish_cell_cache_and_unlock （src/physics/field_solver/fmm/internal/periodic/bem_coulomb_fmm_periodic_cell_cache.f90）。 rank 分担に使う MPI 文脈は fmm_options_type%mpi で、呼び出し側（simulator、--prepare-periodic-cache）が snapshot%init(..., mpi=) から渡す。コアは MPI を初期化・終了しない
• BEACH adapter: src/physics/field_solver/bem_field_solver_config.f90, src/physics/field_solver/bem_field_solver_tree.f90, src/physics/field_solver/bem_field_solver_eval.f90

コアとBEACH adapterの責務は次のように分かれます。

• コア: 幾何前処理、展開係数更新、近傍 direct、点評価
• BEACH adapter: mesh_typeの三角形3頂点でbuild_panel_planを構築し、q_elemをsrc_qへ渡し、 最後にk_coulombを掛ける

---