全部「粒子」で、全部それぞれの差分方程式で解いています。ずるはしていません。3つの階層は 105 倍ずつスケールが違うので、
同時に見せるためにパネルごとに倍率を変えてあります（各パネルに実スケールバーと倍率を表示）。
① 電子：点電子をクーロン力で回しても s/p 軌道は絶対に出ません（古典電子は放射して落ちるだけ）。
なので変分モンテカルロ (VMC) で、多体波動関数 |Ψ|² を粒子（ウォーカー）が厳密にサンプリングします。
スレーター行列式の反対称性が殻構造・パウリ排他・p軌道のローブ・フントの規則を作ります。
② 原子核：陽子と中性子を量子分子動力学 (QMD) のガウス波束で。Skyrme + 対称エネルギー + クーロン + パウリポテンシャル。
③ 核子：3つのクォークを相対論的コーネル型ポテンシャル + Yストリング（格子QCDが支持するシュタイナー木）で。
電子の描画は2モードあります（スライダー「電子の描画」）。どちらも分布は厳密に |Ψ|² です
— メトロポリスの採択判定のおかげで、時間刻み τ を変えても定常分布は変わりません。変わるのは歩き方だけ:
0 = 雲 … τ を大きく取るので相関が速く切れ、密度とエネルギーが速く収束します。
連続する位置はほぼ独立標本なので、線で結んだら嘘になります（だから点でしか描いていません）。
1 = 粒子 … τ を小さく取ると、位置の列が |Ψ|² を定常分布に持つ拡散過程の本物の標本経路になるので、
軌跡付きの動く粒子として正直に描けます。核に引き込む力はドリフト項 D·τ·∇ln|Ψ|² そのものです。
「電子を外から落とす」ボタンを押すと、外側にばらまかれた電子がこの力で引き込まれ、
約150フレームで平衡（酸素なら ⟨r⟩ = 1.00 bohr）に戻ります。代償は相関時間が長くなること（エネルギーの収束が遅い）。
ただし、どちらも「電子の実時間の運動」ではありません。定常状態には古典的な軌道が存在しません
（だから軌跡はギザギザのブラウン運動になります。滑らかな楕円軌道に描いたらそれは嘘です）。
これは正しい分布を運ぶ拡散過程であり、束縛電子の正直な粒子描像です。
入っています。local_energy() の中にカットオフ無しの厳密な 1/r で全ペア:
ただし「力 → 加速度」ではなく、次の2つの経路で分布を決めています。 ① 変分原理：最適化の目的関数 ⟨E_L⟩ にこの項が入っているので、雲の形はクーロンエネルギーが選んでいる。 ② ドリフト項：∇ln|Ψ|² の Ψ は、そのクーロンエネルギーが選んだ Ψ。 ジャストロウ因子の係数 1/2 と 1/4 も、1/r₁₂ の特異性を打ち消すカスプ条件から出た値で、電子間反発そのものが起源です。
検証1：虚位定理。純クーロン系の厳密な基底状態は −⟨V⟩/⟨T⟩ = 2 ちょうど。 運動エネルギーは波動関数のラプラシアンから、ポテンシャルは 1/r の和から、まったく別のコードで計算しているので、 これが 2 の近くに来ること自体が検証になります（酸素で 2.14。2 からのずれが試行関数の粗さを測っています）。
| 酸素 | Hartree |
|---|---|
| ⟨T⟩ 運動（ラプラシアンから） | 64.405 |
| ⟨Ven⟩ 電子‑核引力 | -163.759 |
| ⟨Vee⟩ 電子間反発 | 25.827 |
| −⟨V⟩/⟨T⟩（厳密なら 2） | 2.1416 |
検証2：対相関関数 g(r₁₂) = P(同じウォーカー内) / P(別ウォーカー間)。 分母は同じ1体密度を持ちながら相関だけ無い参照分布なので、1 からのずれは相関そのものです。
| r₁₂ (bohr) | 平行スピン | 反平行スピン | 反平行（ジャストロウOFF） |
|---|---|---|---|
| 0.1 | 0.022 | 0.937 | 0.986 |
| 0.3 | 0.186 | 0.987 | 1.006 |
| 0.5 | 0.655 | 0.970 | 1.001 |
| 0.9 | 1.004 | 0.991 | 1.002 |
平行スピンは r₁₂ → 0 で 0 に落ちます（フェルミホール。これは力ではなく行列式の反対称性が起源）。 反平行は 0 にはならないが抑制されます（クーロンホール。1/r₁₂ の反発が起源）。 ジャストロウを切ると反平行のホールだけが埋まり、平行のホールは残ります。 同時に ⟨Vee⟩ が 25.83 → 28.82 Ha に増えます — 電子が近づけるようになった分、反発エネルギーを余計に払っているということです。
ちなみにこの質問がきっかけで最適化がひとつ良くなりました。虚位定理から全体の長さスケールは s* = −⟨V⟩/(2⟨T⟩) と測定から直接読めるので、粗いグリッド探索の後にこれを入れました （精緻化の前後を同じ統計で測って低い方を採るので、探索より悪くなることはありません）。 平均誤差 3.38% → 3.05%、H₂ は 5.74% → 2.37%、He は 1.60% → 0.75% に改善しました。
正しいです。そしてこれは測って確かめられます。決定的な実験は「経過時間を一定に保ったまま
時間刻み τ を細かくして、経路長が収束するか」です。
滑らかな曲線を粗くサンプルしているだけなら経路長は一定値に収束します。
本物のブラウン運動（どこでも微分不可能）なら τ-1/2 で発散します。
./atom_os walk 7 40 の実測値（酸素、経過時間 T = 40 に固定）:
| τ₀ | ステップ数 | 経路長 L | L·√τ | 曲がり角の cos |
|---|---|---|---|---|
| 0.16 | 250 | 139.3 | 55.7 | -0.109 |
| 0.04 | 1,000 | 278.8 | 55.8 | -0.072 |
| 0.01 | 4,000 | 577.1 | 57.7 | -0.046 |
| 0.0025 | 16,000 | 1,199.2 | 60.0 | +0.003 |
| 0.000625 | 64,000 | 2,396.0 | 59.9 | -0.004 |
時間刻みを256倍細かくしても経路長は収束せず、τ を 1/4 にするたびにきれいに 2 倍になり続けます （不変量は L·√τ）。連続するステップの方向の相関 cos も 0、つまり無相関で、1（滑らか）ではありません。 拡大すると滑らかな曲線が現れるのではなく、もっと細かいジグザグが現れます （ハウスドルフ次元 2 のフラクタル）。だからカクカクは描画の粗さではなく、過程そのものの性質です。
物理的な理由は慣性が無いことです。この差分方程式は時間について 1 階の確率微分方程式で、 ニュートンの運動方程式ではありません。古典粒子の軌道が滑らかなのは速度が連続だからですが、 拡散過程では各瞬間の速度が定義できません（無限大になる）。
ちなみに、本当に滑らかな量子的軌道を与える解釈もあります：ボーム力学の v = ∇S/m（S は波動関数の位相）です。ところがここで扱っている基底状態の波動関数は すべて実数なので S = 0、つまりボーム速度は厳密に 0 — 電子は完全に静止してしまいます。 2p の m = ±1 のような複素軌道なら滑らかな円運動が出ます。 「滑らかに回る電子」を見たいなら、そこまで含めて選ぶ必要がある、ということです。
なお、画面に描いている折れ線は1 本の線分がちょうど 1 メトロポリスステップです （1フレームに複数ステップ進む設定でも、途中の位置を飛ばして直線で結んでしまわないようにしてあります。 そうすると実際より滑らかに見えてしまうので）。
① 電子：メトロポリス・ドリフト拡散差分方程式
試行関数（原子単位 ħ = me = e = 1）：
aij の値は電子-電子カスプ条件そのもの（勝手な値ではありません）。 f(R) = |ΨT(R)|² を定常解にもつフォッカー・プランク方程式
df/dt = D grad·( grad f − F f ), F = grad ln|Psi_T|^2 = 2 grad Psi / Psiに対応するランジュバン方程式を伊藤離散化すると（1電子ずつ動かす）：
r_i^{n+1} = r_i^n + D tau_i F_i(R^n) + sqrt(2 D tau_i) xi , xi ~ N(0,1)^3, D = 1/2τ が有限だとこれはバイアスを持つので、提案密度
G(r->r') = (4 pi D tau)^{-3/2} exp( −|r' − r − D tau F(r)|^2 / (4 D tau) ) A = min( 1, |Psi(R')|^2 G(r'->r) / ( |Psi(R)|^2 G(r->r') ) )でメトロポリス・ヘイスティングス補正します。これでどんな τ でも分布は厳密に |ΨT|²（時間刻み誤差ゼロ）。 τ は原子核の近くで小さく、価電子領域で大きくしてありますが、逆向きの遷移で τ(r′) を使うので厳密性は保たれます （だから両方の規格化因子をきちんと残しています）。1電子移動の Ψ 比は逆行列を使って O(N) で：
ratio_det = sum_j Phi_j(r_i') · Ainv[j][i] grad ln det|r_i' = ( sum_j Ainv[j][i] grad Phi_j(r_i') ) / ratio_det観測量（＝これが正しさの証拠）：
E_L(R) = −1/2 sum_i (lap Psi/Psi)_i + V(R) (lap Psi/Psi)_i = (lap det/det)_i + lap_i J + |grad_i J|^2 + 2 (grad_i ln det)·(grad_i J) V = −sum_A Z_A/|r_i−R_A| + sum_{i<j} 1/r_ij + sum_{A<B} Z_A Z_B/R_AB<EL> は変分エネルギー、つまり真の基底状態の上限です。画面には厳密な非相対論的値と並べて表示しています。 誤差棒はブロック平均法（連続サンプルは相関するので √(var/N) では過小評価になる）。
② 核子：QMD 差分方程式
各核子は幅 L の固定ガウス波束。積状態での期待値をとると古典ハミルトニアンになります：
パウリポテンシャルが古典核子MDにおける反対称性の代役で、同種核子を同じ位相空間セルから追い出しフェルミ運動を供給します （Wilets, Dorso–Randrup の方法）。H が P に運動項以外で依存するので dR/dt = ∂H/∂P ≠ P/M、 したがって普通の leapfrog は使えません。導出される更新式は陽的中点法（2次ルンゲ・クッタ）：
R* = R^n + (dt/2) ∂H/∂P(R^n,P^n) , P* = P^n − (dt/2) ∂H/∂R(R^n,P^n) R^{n+1} = R^n + dt ∂H/∂P(R*,P*) , P^{n+1} = P^n − dt ∂H/∂R(R*,P*)使う解析的な微分（d rho_ij/dR_i = −rho_ij (R_i−R_j)/(2L)）：
dE_sk2/dR_i = (alpha/rho0) sum_j d rho_ij/dR_i dE_sk3/dR_i = (beta g/((g+1) rho0^g)) sum_j (rb_i^{g−1} + rb_j^{g−1}) d rho_ij/dR_i d/dr [erf(r/a)/r] = 2 exp(−r^2/a^2)/(a sqrt(pi) r) − erf(r/a)/r^2 dVp/dR_i = −Vp (R_i−R_j)/q0^2 , dVp/dP_i = −Vp (P_i−P_j)/p0^2基底状態は同じ H の位相空間最急降下で作ります（R と P の両方を降下させる：パウリ項があるので H の最小点は P ≠ 0）。ここで大事な事実：完全に緩和した古典核は静止します（最小点では dR/dt = ∂H/∂P = 0）。 実際の原子核は核子あたり約20 MeV の量子的ゼロ点フェルミ運動を持っていますが、古典モデルはそれを「運動」として持てません。 したがって画面で見えている運動は、明示的に与えた励起エネルギー E* と集団回転 J です（両方とも実在の物理量で、数値を表示しています）。
③ クォーク：相対論的シンプレクティック差分方程式
H = sum_i sqrt(p_i^2 + m_i^2) + V(x) V = −kappa sum_{i<j} erf(r_ij/r0)/r_ij + sigma · L_Y(x) + V0 kappa = (2/3) alpha_s hbar c … 3クォーク色一重項の色因子 L_Y = 最小全ストリング長 = 3点のシュタイナー（フェルマー）木H = T(p) + V(x) と分離できるので、厳密にシンプレクティックな kick–drift–kick が使えます （エネルギーはドリフトなしで保存。画面に相対ドリフトを表示しています）：
p^{n+1/2} = p^n − (dt/2) grad V(x^n) x^{n+1} = x^n + dt · p^{n+1/2} / sqrt((p^{n+1/2})^2 + m^2) p^{n+1} = p^{n+1/2} − (dt/2) grad V(x^{n+1})dx/dt = p/E はまさに相対論的速度 v/c。ストリングの勾配は包絡線定理で、接合点 S が最小点にあるので
dL_Y/dx_i = (x_i − S)/|x_i − S| … 大きさ1の単位ベクトルつまりどれだけ引き離しても内向きの力は σ で一定＝クォークは絶対に外に出られない。 これは押しつけた仕様ではなく σ·L という項の帰結です（角度が120°以上になると接合点はその頂点に退化するので、 その場合の勾配も別に導出して実装しています。ここを手抜きするとエネルギー保存が壊れます）。
核子の大きさはどこから来るか：円軌道の力の釣り合い
(A) p^2 / sqrt(p^2 + m^2) = (R/3) dV_eq/dRだけでは決まりません（どの R でも解があり、古典的には R → 0 に潰れる）。核子の大きさは量子的なので、 水素原子のボーア半径を決めるのと同じボーア＝ゾンマーフェルト条件で閉じます：
(B) p R = n hbar （n = 1 が基底のピンホイール）(B) を (A) に代入すると R について1本の方程式になり、二分法で解いています。R も p も手で入れていません。 V0 は構成子クォーク模型の標準的な加算定数で、n = 1 のときに核子質量に合わせて一度だけ決め、以後触りません。 だから n = 2, 3 … の軌道は本当に重くなります（実在の N* 共鳴と同じ）。
① 電子 VMC 全エネルギー（Hartree） 厳密値は非相対論的極限の文献値 [Chakravorty et al., PRA 47, 3649 (1993)]。EVMC ≥ E厳密 が変分原理。
| 系 | 電子数 | EVMC | E厳密 | 誤差 |
|---|---|---|---|---|
| H | 1 | -0.5000 | -0.5000 | 0.00% |
| He | 2 | -2.8819 | -2.9037 | 0.75% |
| O | 8 | -73.765 | -75.067 | 1.73% |
| Ne | 10 | -126.795 | -128.938 | 1.66% |
| Ar | 18 | -514.673 | -527.540 | 2.44% |
| H₂ | 2 | -1.1466 | -1.1745 | 2.37% |
| H₂O | 10 | -72.476 | -76.438 | 5.18% |
| O₂ | 16 | -147.593 | -150.327 | 1.82% |
| N₂ | 14 | -107.411 | -109.542 | 1.95% |
| CH₄ | 10 | -38.464 | -40.515 | 5.06% |
| NH₃ | 10 | -52.251 | -56.564 | 7.62% |
| HF | 10 | -95.274 | -100.459 | 5.16% |
| CO | 14 | -108.842 | -113.317 | 3.95% |
H₂ は分離した水素原子2個（-1.0000）より低い＝化学結合がモデルから出てきている（結合エネルギー 0.147 Ha = 3.99 eV、実験値 4.75 eV）。
② 原子核 QMD 基底状態 パラメータは α, β, γ（核物質の飽和条件で決まる）を固定したまま、 L, Cp, q₀, p₀, Cz, az, Cs を座標降下で下の2列だけに合わせたもの（相対誤差RMS 6.5%）。
| 核種 | E/A モデル | E/A 実験 | 誤差 | rrms モデル | rrms 実験 | 誤差 |
|---|---|---|---|---|---|---|
| He-4 | -5.597 | -7.074 | +20.9% | 1.232 | 1.428 | -13.7% |
| Be-9 | -6.906 | -6.463 | -6.9% | 2.264 | 2.361 | -4.1% |
| C-12 | -7.693 | -7.680 | -0.2% | 2.229 | 2.309 | -3.5% |
| O-16 | -8.097 | -7.976 | -1.5% | 2.415 | 2.553 | -5.4% |
| Ne-20 | -8.261 | -8.032 | -2.9% | 2.727 | 2.875 | -5.1% |
| Mg-24 | -8.284 | -8.261 | -0.3% | 3.010 | 2.929 | +2.8% |
| Si-28 | -8.413 | -8.448 | +0.4% | 3.131 | 2.997 | +4.5% |
| S-32 | -8.480 | -8.493 | +0.2% | 3.262 | 3.141 | +3.9% |
| Ar-40 | -8.640 | -8.595 | -0.5% | 3.507 | 3.313 | +5.9% |
単位は MeV と fm。実験の rrms は電荷半径から陽子自身の大きさを差し引いた点核子半径 √(Rch² − 0.769 fm²)。
③ クォーク（基底のピンホイール、n = 1）
| 量 | 値 | 備考 |
|---|---|---|
| 核子質量 | 938.918 MeV | V₀ をこれに合わせて1回だけ決定 |
| V₀ | -1581.5 MeV | 構成子クォーク模型の加算定数 |
| 軌道半径 R | 0.3677 fm | (A)+(B) から導出、手入力なし |
| 運動量 p | 536.7 MeV/c | 同上 |
| 速度 v/c | 0.848 | だから相対論が必要 |
| ストリング長 LY | 1.103 fm | = 3R |
| エネルギー保存 | 1.0 × 10-13 | 40万ステップ後の相対ドリフト |
核の時間発展のエネルギー保存も同様に測っていて、Δt = 0.20 fm/c で4000ステップ後の相対ドリフトは 3.9 × 10-5
（画面の dE がその実測値です）。