B-1. 計算物理学のための Julia 入門

金賀 穂

千葉大学

自己紹介

担当

  • 金賀 穂
  • 千葉大学 大学院理学研究院
  • 後半 2 コマを担当

研究キーワード

  • 物性理論
  • 光物性
  • 磁性・スピントロニクス
  • 非平衡ダイナミクス

研究紹介

進め方

今日の立ち位置

  • 前半で扱った Julia / コーディングエージェントを、さらに実践的な研究コードの読解と穴埋めへ接続する
  • 数式をそのまま Julia 実装へ落とし、結果を見て物理解釈する

今日の到達点

  • 前半で 01_bands.png02_timeevol_current.png を再現する
  • 後半で 03_hhg_fft.png04_selection_rule.png を再現する

事前準備

  • 配布リポジトリを clone して作業ディレクトリへ入る
  • 依存関係は最初に一度だけ入れる
julia --project=. -e 'using Pkg; Pkg.instantiate()'

checkpoint の使い方

  • 例えばハンズオン1 で詰まったら checkpoint-1-band
git switch --detach checkpoint-1-band

第1コマの流れ

0

セットアップ確認、checkpoint の使い方、AI の使い方

1

高次高調波発生やGKSL方程式等の基本概念の導入

2

ハンズオン1: バンド図を作る

3

ハンズオン2: 電流を計算して時間発展をプロットする

高次高調波発生(HHG)

HHG → 最も単純な非線形光学効果の一つ

HHG time evolution

最初に持ってほしいイメージ

  • 高強度レーザーを入れると、固体電子系が非線形に応答する
  • その応答は電流や分極の時間変化として現れる
  • 高調波は、周波数空間へ写した結果として観測される

HHG time evolution

HHG time evolution

HHG time evolution

HHGの数値解析を実践してみましょう!

グラフェン

graphene lattice

実空間で持つ自由度

  • 黒丸と白丸が A/B 副格子

電子の生成消滅演算子

  • \(\hat a_{\bm r}\) / \(\hat b_{\bm r}\)
    → A/B 副格子の電子の消滅演算子

タイトバインディング模型

\[ \hat H_0 = -t_0 \sum_{\bm r}\sum_{j=1}^3 (\hat b_{\bm r+\bm d_j}^\dagger \hat a_{\bm r} + \text{h.c.}) + \Delta\sum_{\bm r} (\hat a_{\bm r}^\dagger \hat a_{\bm r} - \hat b_{\bm r}^\dagger \hat b_{\bm r}) \]

\[\begin{align} &\acomm{\hat c_\br}{\hat c'_{\br'}}=\delta_{cc'}\delta_{\br\br'}\\ &\comm{A}{B}=AB-BA\\ &\acomm{A}{B}=AB+BA \end{align}\]

逆格子と Dirac cone

graphene BZ and gap

波数空間で何をするか

  • 計算は逆格子原胞(ひし形)で進める
  • 図では 六角形 BZ に切り出して表示する

Δ の役割

  • Δ = 0 では Dirac cone が閉じる(ギャップレス)
  • Δ ≠ 0 では gap が開く

\[ {\small \hat a_\br = \frac{1}{\sqrt{N}}\sum_{\bk}\tilde a_\bk e^{-i\bk\cdot\br}, \quad \hat b_\br = \frac{1}{\sqrt{N}}\sum_{\bk}\tilde b_\bk e^{-i\bk\cdot\br} } \]

\[ \hat{H}_0=\sum_\bk\bm{C}^\dag_\bk \mathcal H(\bk)\bm{C}_\bk, \quad \bm{C}_\bk=\mqty(\tilde{a}_\bk & \tilde{b}_\bk)^\top \]

価電子帯と伝導帯

\[ \mathcal H(\bm{k})= \begin{pmatrix} \Delta & -t f(\bm{k})\\ -t f^*(\bm{k}) & -\Delta \end{pmatrix} \]

\[ {\small f(\bm{k})=\sum_{j=1}^3 e^{i\bm{k}\cdot\bm{d}_j} } \]

対角化

\[ E_\pm(\bm{k})=\pm\sqrt{\Delta^2+t^2|f(\bm{k})|^2} \] \[ \mqty(\xi_\bk & \zeta_\bk)=U^\dag_\bk\bm{C}_\bk \]

\[\begin{align*} \ket{v_\bk} =\zeta^\dag_\bk\ket{0_\bk}\\ \ket{c_\bk} =\xi^\dag_\bk\ket{0_\bk} \end{align*}\]

読み方

  • ±Δ が A/B 副格子の非対称性
  • -t f(k) が最近接ホッピングの総和

このあと重要

  • H(k) が分かれば dH/dk も解析的に書ける
  • その dH/dk が、あとで current_traces へ入る

ハンズオン1: バンド構造を描く

対象

  • src/tb.jl
  • examples/01_bands.jl

到達目標

  • f(k), dfdk, H(k), dHdk を埋める
  • examples/out/01_bands.png を出す
  • 余裕があれば test/test_tb.jl に TB テストを 2 本書く

六角形BZ上のバンド図

図: examples/01_bands.jl

詰まったら

  • checkpoint-1-band

ハンズオン1 詳細: f(k)H(k) の穴埋め

\[ f(\bm{k})=\sum_{j=1}^3 e^{i\bm{k}\cdot\bm{d}_j} \]

\[ \mathcal H(\bm{k})= \begin{pmatrix} \Delta & -t f(\bm{k})\\ -t f^*(\bm{k}) & -\Delta \end{pmatrix} \]

function f(k, dvecs = NN_VECTORS)::ComplexF64
    # ここを実装
end

function H(k, p::TBParams)::CMat2S
    # ここを実装
end

手元のコードでの見方

  • 空欄の前後にある型注釈と TBParams を先に読む
  • 数式を Julia へ写す

ハンズオン1 詳細: dfdk, dHdk

\[ \frac{\partial f}{\partial k_\alpha} = i\sum_j d_{j,\alpha} e^{i\bm{k}\cdot\bm{d}_j} \] \[ \pdv{\mathcal H}{k_\alpha}= \begin{pmatrix} 0 & -t \pdv{f}{k_\alpha}\\ -t \pdv{f^*}{k_\alpha} & 0 \end{pmatrix} \]

function dfdk(k, dvecs::NTuple{3,Vec2}=NN_VECTORS)::SVector{2,ComplexF64}
    # ここを実装
end
function dHdk(k, p::TBParams)::Tuple{CMat2S,CMat2S}
    # ここを実装
end

最初から埋まっている周辺コード

  • NN_VECTORS, b1, b2: 格子定義
  • hex_vertices(b1, b2): 六角形 BZ の頂点
  • hex_surface_coordinates(...): 表示用の座標生成
  • band_surfaces(kx, ky, tb): 各点で固有値を取って上下バンドを作る

実装時の注意

  • dfdk[1] が x, dfdk[2] が y
  • dHdk は電流用に使う

ハンズオン1 追加: ユニットテストを書く

バンド図が出た人向け

  • examples/01_bands.jl で図が出たら取り組む
  • checkpoint 到達条件ではなく、早い人向けの追加ミニ課題
  • test/test_tb.jl に 3 本書けばよい
julia --project=. -e 'using GrapheneHHG; include("test/test_tb.jl")'

進め方

  • バンド図完成
  • test/test_tb.jl を読む
  • ユニットテストを 3 本書く
  • 個別実行で確認する

ここでの位置づけ

  • この段階ではハミルトニアンだけを確かめる
  • repo 全体の Pkg.test() はハンズオン2で初めて回す

ハンズオン1 詳細: test/test_tb.jl で見ること

1 本目: エルミート性

  • H(k) がエルミート性を保っているか

2 本目: ±E 対称

  • グラフェン 2 バンド模型の基本的なバンド対称性を守れているか
  • 数式をコードへ写した結果が、固有値にも反映されているかを見る

3 本目: \(\bm K=\frac{2\pi}{3\sqrt 3}(1,\sqrt 3)\) でギャップレス

  • \(\Delta=0\)のとき、\(\bm K\)点でギャップレスとなる

AI活用 Prompt: ハンズオン1 の穴埋め前確認

GitHub Copilot Agent

利用文脈: 手元の src/tb.jl を埋める前に、式とコードの対応を確認する。

src/tb.jl の空欄を埋める前に理解を整理したいです。
`NN_VECTORS`, `TBParams`, 型注釈から確実に言えることを先にまとめ、
その後 f, dfdk, H, dHdk について、数式との対応、満たすべき制約
(複素共役、x/y順、エルミート性、`Δ = 0` の固有値対称性)、
候補コード、自分で確認する問いとテスト観点の順に説明してください。
確実に言えることと補助的な提案は分けて書いてください。

大事なのは、AIの実装の物理骨格をちゃんと理解し、自分のものとすること。(≠1行1行をチェック)
数式との対応を意識して、コードの意味をちゃんと理解することが重要です。
また、関数が満たすべき物理的制約をきちんと満たしているか確認するユニットテストを追加すべし。

FB1: バンド図から何を読むか

六角形BZ上のバンド図

図: examples/01_bands.jl

まず確認すること

  • 六角形 BZ の角が \(\bm K\)/\(\bm K'\)1に対応しているか
  • Δ = 0 なら\(\bm K\)/\(\bm K'\)点付近で cone が閉じているか
  • Δ ≠ 0 に変えたとき gap が開くか

この段階で掴めればよいこと

  • H(k) の形がバンド図へどう写るか

次に追加: 静的なバンド図から J(t)

ここまで

  • H(k)dH/dk が書けた
  • 六角形 BZ 上でバンド構造を描けた

次に足すもの

  • A(t) による時間依存
  • L_k による緩和
  • ρ_k(t) の時間発展
  • current_traces による Jx(t), Jy(t) の観測

パイエルス置換

レーザーの効果をどう入れるか?

→ パイエルス置換 \[ t_0 \rightarrow t_{\br,\br+\bm d_j}(t) = t_0\exp\left(-ie\int_\br^{\br+\bm d_j} \bm{A}(t)\cdot d\br\right) \]

\[\begin{align} \hat{H}(t)= & -\sum_{\bm{r},j}(t_{\br,{\br+\bm d_j}}(t)\hat b_{\br+\bm d_j}^\dag\hat a_\br + t_{{\br+\bm d_j},\br}(t)\hat a_\br^\dag\hat b_{\br+\bm d_j}) + \Delta\sum_\br(\hat a^\dag_\br\hat a_\br - \hat b^\dag_{\br+\bm d_1}\hat b_{\br+\bm d_1}) \\ =&\sum_\bk\bm{C}^\dag_\bk \mathcal H(\bk+e\bm{A}(t))\bm{C}_\bk \\ \mathcal H_\bk(t) &:= \mathcal H(\bk+e\bm{A}(t)) \end{align}\]

\[ \bm{k}\rightarrow \bm{k}+\bm{A}(t), \qquad \bm{E}(t)=-\partial_t\bm{A}(t) \]

pulse = default_pulse(; A0 = 0.5, ω0 = 0.7, fwhm_cycles = 5.0)
At = A(t, pulse)

本講義では\(x\)軸方向直線偏光として実装してある。

量子マスター(GKSL)方程式

\[ \dot{\rho}_\bk(t)= -i[\mathcal H_\bk(t), \rho_\bk(t)] +\gamma\left( L_\bk \rho_\bk(t) L_\bk^\dagger -\frac{1}{2}\{L_\bk^\dagger L_\bk,\rho_\bk(t)\} \right) \]

なぜGKSL型のマスター方程式を使うのか?

  • 自然に緩和を取り込める
  • マルコフ性/正値性・トレース保存により確率解釈が可能

このモデルで入れるもの

  • 各時刻でハミルトニアンの波数 kk + A(t) へずらして評価する
  • ρ_k(t) は 2×2 密度行列
  • \(L_\bk\)は緩和過程を表す演算子であり、本講義では\(L_\bk = |v_\bk\rangle\langle c_\bk|\)と設定する
  • \(L_\bk\)A=0 の固有状態から前計算する

ハンズオン2: J(t) を計算・プロットする

必須で編集

  • src/rhs.jl
  • src/observables.jlcurrent_traces

読むだけ

  • src/lindblad.jl
  • src/workflows.jl
  • examples/02_timeevol_current.jl

到達目標

  • rhs!current_traces を埋めて ρ_k(t) の右辺を組み立てる
  • 02_timeevol_current.png を出す

current Jx(t) and Jy(t)

図: examples/02_timeevol_current.jl

詰まったら

  • checkpoint-2-rhs

ハンズオン2 詳細: GKSL の各項と rhs!

\[ {\small \dot{\rho}_\bk= -i[\mathcal H_\bk(t), \rho_\bk] +\gamma\left( L_\bk \rho_\bk L_\bk^\dagger -\frac{1}{2}\{L_\bk^\dagger L_\bk,\rho_\bk\} \right) } \]

function commutator!(dρ, Hk, ρk)::Nothing
    # ここを実装
end

function dissipator!(dρ, ρk, Lk, LdLk, γ)::Nothing
    # ここを実装
end

function rhs!(dρ, ρ, p, t)::Nothing
    # ここを実装
end

ここで押さえること

  • \(-i[\mathcal H, \rho]\) がユニタリな回転
  • \(L \rho L^\dagger\)\(-\{L^\dagger L,\rho\}/2\) が散逸項
  • \(L_\bk\)\(L_\bk^\dagger L_\bk\)A=0 の固有状態から前計算する
  • 毎時刻に固有分解しないことで、実装を軽く保つ

ハンズオン2 詳細: current_traces

@inline function _current_component(ρk, dHdk)::Float64
    # ここを実装
end

function current_traces(ρ, dHdk::NamedTuple{(:x, :y)})
    # ここを実装(_current_componentを使って書く)
end

ここで押さえること

  • A 微分、したがって k 微分したハミルトニアンが応答演算子になる
  • この教材では \(\expval{\bm J_\bk(t)}=\Tr[\rho_\bk(t)\pdv{\mathcal H_\bk(t)}{\bk}]\)real(tr(ρ * dHdk))) を各 k で評価して平均する

実装時の注意

  • dHdk.x[i]dHdk.y[i] の取り違えだけ注意する

ハンズオン2 詳細: Pkg.test()

julia --project=. -e 'using Pkg; Pkg.test()'
for i in eachindex(dρ)
    # ここを実装(トレース保存・エルミート性)
end

for LdL in cache.LdL
    # ここを実装(L^\dagger Lは射影演算子?)
end

test/test_gksl.jl で見ること

  • tr(dρ)=0: 規格化を壊していないか
  • エルミート性: 密度行列らしさを保てているか
  • \(L^\dagger L\) が projector 的な性質を持つか

ハンズオン2での位置づけ

  • rhs! を埋めて波形が出たら、そのまま Pkg.test() で確認する
  • 図が出るだけでなく、GKSL の基本性質も崩れていないかを見る

AI活用: 実装前に確認すること

GitHub Copilot Agent

利用文脈: 手元の src/rhs.jlsrc/observables.jl を埋める前に、GKSL の各項と rhs!current_traces の役割を整理する。

`src/rhs.jl` と `src/observables.jl` を理解しながら実装したいです。
まず `commutator!`, `dissipator!`, `rhs!`, `current_traces` の役割とデータの流れを整理し、
その後 GKSL 各項との対応、`A(t)`, `γ`, `L_k`, `kshift` が効く場所、
候補コード、確認観点(`tr(dρ)=0`, エルミート性, `L^\dagger L`, `J(t)`)
の順に説明してください。最後に `current_traces` の集計内容を 2〜3 行で補足し、
確実に読めることと推測は分けてください。

FB2: J(t) の波形をどう読むか

current Jx(t) and Jy(t)

図: examples/02_timeevol_current.jl

ここで観察すること

  • Jx(t) は包絡中心で振幅が大きくなる
  • Jx(t) には包絡の上に高速振動が重なっている
  • \(\Delta=0\)では Jy(t) は強く抑制される
  • この図の Jy10^{-16} 程度なので、ほぼ数値ノイズとして読む

発展課題

  • \(\Delta\neq0\)でどう変わるか?

後半への橋渡し

前半でできたこと

  • H(k)dH/dk を式とコードで結んだ
  • k + A(t) による光駆動を入れた
  • GKSL で ρ_k(t) を開放系として更新した
  • Jx(t), Jy(t) という観測対象を得た

後半でやること

  • J(t) と動的対称性の関係を確認する
  • J(t) を FFT して HHG スペクトルへ変換する
  • Δ = 0 / Δ ≠ 0 の比較から選択則を読む

第2コマの流れ

ここまでに作ったもの

  • examples/out/01_bands.png
  • examples/out/02_timeevol_current.png

ここまでで分かっていること

  • H(k)dH/dk はコードで書ける
  • A(t) を入れて GKSL方程式を書くと J(t) が出る
  • \(\Delta=0\)では JxJy の振る舞いに差がある

第2コマの出力

  • examples/out/03_hhg_fft.png
  • examples/out/04_selection_rule.png

流れ

  • J(t) の再確認
  • ハンズオン3
  • FB3
  • ハンズオン4
  • FB4
  • 余裕があれば CI / README / 公開の話

J(t) をもう一度見る: 動的対称性の入口

current in time domain

図: examples/02_timeevol_current.jl

いま見ること

  • Jx がどのような周波数成分を含んでいるか
  • フーリエ変換によって周波数空間に移ると、どんな構造が見えるか

ハンズオン3: FFT で HHG スペクトルを出す

必須で編集

  • src/fft.jl

読むだけ

  • examples/03_hhg_fft.jl

到達目標

  • hhg_spectrum を埋めて n = ω/ω0 軸のスペクトルを作る
  • examples/03_hhg_fft.jl で返ってきたスペクトルを描いて保存する
  • examples/out/03_hhg_fft.png を保存する

hhg spectrum from fft

図: examples/03_hhg_fft.jl

詰まったら

  • checkpoint-3-fft

ハンズオン3 詳細: hhg_spectrum

function hhg_spectrum(Jt, dt, ω0)
    # ここを実装
end

ここですること

  • Jtから平均値を引いて 0 次成分を落とす
  • FFTW.jlrfft関数で正の周波数成分だけフーリエ変換する
  • dt とサンプル数から周波数軸を作る
  • ω0 で割って横軸を次数 n に変換する
  • power=amplitude^2 をプロット用の強度として使う
  • (ω, n, power) を返す

注意

  • 軸の細かさは入力長 Ndt で決まる

AI活用 Prompt: ハンズオン3 の前確認

GitHub Copilot Agent

利用文脈: src/fft.jlhhg_spectrum を埋める前に、周波数軸の意味と確認項目を整理する。

`src/fft.jl` の `hhg_spectrum` を、入出力の意味を理解しながら実装したいです。
`Jt`, `dt`, `ω0` と戻り値が何を表すかを先に整理し、その後
処理の順序(平均値除去, `rfft`, 周波数軸, 次数変換, power)、
`dt` とサンプル数が分解能や見える次数範囲にどう効くか、
候補コード、`examples/03_hhg_fft.jl` で確認する問い
の順に説明してください。確実に言えることと予想は分けてください。

FB3: HHGスペクトル

hhg spectrum from fft

examples/03_hhg_fft.jl

ここで見ること

  • Jx のピークが n=整数 で見える
  • 奇数次と偶数次の出方の違いが見える
  • Jy がほぼ 0 なら、周波数空間でもほぼゼロ

問い

  • 何が Jx,Jyのスペクトルの形状を支配しているのか?

動的対称性と選択則: 一般論

定義

\[ U^\dagger H(t + \Delta t) U = H(t) \]

\[ U^\dagger \hat{\bm O}(t + \Delta t) U = R \hat{\bm O}(t) \]

  • 空間操作 \(U\) と時間並進 \(0<\Delta t<T\) を組にした対称性
  • ここで、レーザー周期を \(T = 2\pi/\omega_0\) とする

HHG への言い換え

  • 半周期ずらしたときに電流が符号反転するなら、

\[ U^\dagger \hat J_\alpha(t + T/2) U = -\hat J_\alpha(t) \]

  • フーリエ成分は偶数次で打ち消し合う

\[ J_\alpha(2n\omega_0) = 0 \]

  • つまり禁制になる次数を対称性だけで読める

x方向直線偏光 + \(\Delta=0\)の場合

動的対称性と選択則まとめ

\(\Delta = 0\)

  • 残る対称性: \((U_{I_x}, t + T/2)\), \((U_{I_y}, t)\)

\[ J_x(2n\omega_0)=0,\qquad J_y(t)=0 \]

  • 観測: \(J_x\) は奇数次中心、\(J_y\) はほぼ 0

\(\Delta \neq 0\)

  • 残る対称性: \((U_{I_x}, t + T/2)\) のみ

\[ J_x(2n\omega_0)=0,\qquad J_y((2n+1)\omega_0)=0 \]

  • 観測: \(J_y\) の偶数次が立ち上がる

連続波で厳密に成り立つ議論だが、今回の有限パルスでも定性的にはこの選択則で図を読める。

ハンズオン4: Δ=0 / Δ≠0 を比べる

対象

  • examples/04_selection_rule.jl

主題

  • Δ=0 / Δ≠0 の比較から相対変化を読み取る

到達目標

  • Δ=0Δ=0.2 のスペクトルを並べて比較する
  • examples/out/04_selection_rule.png を保存する

selection comparison from example

図: examples/04_selection_rule.jl

ハンズオン4 詳細: Δ は何を変えるか

\[ \mathcal H(\bm{k})= \begin{pmatrix} \Delta & -t f(\bm{k})\\ -t f^*(\bm{k}) & -\Delta \end{pmatrix} \]

Δ = 0

  • 副格子が対等
  • Dirac cone が閉じる

Δ ≠ 0

  • A/B サイトが非対等
  • gap が開く
  • それまで抑えられていた応答成分が現れうる

ハンズオン4 詳細: コード

# ここは固定済み: 時間発展から FFT までは変えない
function run_spectrum(Δ; Nk = 20, γ = 0.1, A0 = 0.5, ω0 = 0.7, dt = 0.02)
    result = simulate_currents(TBParams(t = 1.0, Δ = Δ); Nk = Nk, γ = γ, A0 = A0, ω0 = ω0, dt = dt)
    return hhg_spectra(result.currents, result.dt, result.pulse.ω0)
end
# handson3 と同じ仕様に固定し、
# "Δ=0.0"と"Δ=0.2"を重ねてプロットしてみよう

周辺コードの役割

  • run_spectrum が「時間発展から FFT まで」を一気につなぐ
  • plot_selection_ruleJx/Jy を上下 2 段で描く
  • 比較条件は Δ 以外固定
  • handson3 と同じ軸仕様で比較し、違いは Δ による変化へ寄せる

AI活用 Prompt: ハンズオン4 の前確認

GitHub Copilot Agent

利用文脈: 比較図に基づく観測事実と、対称性の説明を紐づける。

`examples/out/04_selection_rule.png` を理解したいです。
まず `run_spectrum` など比較図まわりのコードを読み、`Δ` 以外の条件が固定であることを確認してください。
その後、先に読む周辺コードと比較条件、図から直接言える観測事実
(`Δ=0` と `Δ≠0` の Jx/Jy の違い)、各事実に対応する説明候補
(動的対称性・副格子非対称性)、図だけでは断定できない点の順に整理してください。
観測事実と解釈は分け、推測には推測と明記してください。

FB4: Δ=0 / Δ≠0 比較から何を言うか

selection comparison from example

図: examples/04_selection_rule.jl

この図から言えること

  • Δ を変えても Jx はほとんど変わらない
  • Jy は偶数次応答 n=2,4,6,... が現れる
  • 動的対称性の帰結と整合している

動的対称性: \(\Delta \neq 0\) の場合

壊れる対称性

  • \((U_{I_y}, t)\)\(\Delta=0\) のときだけ成り立つ
  • \(\Delta\neq0\) では A/B 副格子が非対等になり、\(y\to -y\) で元に戻らない
  • したがって \(J_y(t)=0\) はもう要請されない

残る対称性

  • \((U_{I_x}, t + T/2)\)\(\Delta\) に依らず残る
  • \(J_x(2n\omega_0)=0\), \(J_y((2n+1)\omega_0)=0\)
  • 図では \(J_x\) はほぼ重なり、\(J_y\) の偶数次が立ち上がる

周辺話題いくつか

PkgTemplates.jl: 雛形は入口をそろえる

先にそろえておきたい器

  • Project.toml と package 名
  • src/ と公開入口
  • test/runtests.jl
  • README と CI の土台

手で始めると後回しになりやすい入口を、最初にまとめて固定する。

PkgTemplates の役割

  • 最低限の構成を最初に固定する
  • src/, test/, README, CI の迷いを減らす
  • 実装そのものの代わりではない

雛形だけでは足りない

  • 中身の妥当性は自分で検証する
  • 物理の説明や例題の読み方は別に要る
  • 雛形は入口であり、研究コードそのものではない

可視化依存は ext に分ける

将来 package 化すると起こること

  • plot_hhg のような標準可視化関数を足したくなる
  • ただし MakiePlots を本体 [deps] に入れると、計算だけ使いたい人にも重い依存が乗る
  • 本体では関数名だけ公開し、描画実装は ext/ に逃がす
  • そうすると backend を後から増やしやすく、コア計算も軽く保ちやすい
[weakdeps]
CairoMakie = "..."

[extensions]
GrapheneHHGCairoMakieExt = "CairoMakie"
# src/GrapheneHHG.jl
function plot_hhg end

# ext/GrapheneHHGPlotExt.jl
GrapheneHHG.plot_hhg(spec) = ...

ドキュメント 1: README

README で最初に知りたいこと

  • これは何を計算する repo か
  • どう環境構築し、どう実行するか
  • 何が出力され、どの図が再現できるか

README が入口になる理由

  • 受講者はまず「何を打てば動くか」を見る
  • 査読者はまず「第三者が再現できるか」を見る
  • 将来の自分も、最初に README を読む

ドキュメント 2: 詳細 docs

役割の分け方

置く場所 主な役割
README 導入、インストール、最短の実行手順
詳細 docs 背景、API、図の読み方、設計の理由
examples/ 実際に再現する最短経路

分ける利点

  • README を短く保てる
  • 背景説明を削らずに置ける
  • API の説明を source と同期しやすい

この repo で対応する場所

  • README.md: セットアップ、checkpoint、examples の入口
  • docs/src/index.md: API 一覧と @docs / @autodocs
  • docs/make.jl: HTML 生成と公開設定
  • examples/02_timeevol_current.jl: simulate_currents の利用例

読む順番

  • README で「何を動かすか」をつかむ
  • docs で関数や型の意味を確認する
  • examples/ で図を再現する

ドキュメント 3: docstring から API docs へ

"""
授業で使う既定のガウシアンパルスを作る。
# キーワード引数
- `A0`, `ω0`, `t0`, `fwhm_cycles`
# 戻り値
- `PulseParams`
"""
function default_pulse(; A0=0.5, ω0=0.7, t0=0.0, fwhm_cycles=5.0)
    ...
end
## 講義で使う主要 API
```@docs
PulseParams
default_pulse
A
```

この repo での使い方

  • src/pulse.jl の説明を """ ... """ にする
  • docs/src/index.md@docs に並べる
  • ?default_pulse や HTML docs から同じ説明を読める

docstring に入れるもの

  • 役割
  • 主要引数
  • 戻り値
  • 必要なら補足

ドキュメント 4: Documenter.jl + GitHub Pages

# docs/make.jl
makedocs(;
    modules=[GrapheneHHG],
    pages=["Home" => "index.md"],
)

deploydocs(; repo="github.com/xxxxxxxx/GrapheneHHG.jl", devbranch="main")
# .github/workflows/ci.yml
- uses: julia-actions/julia-docdeploy@v1
  env:
    GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
    DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
1

docs/src/ を読む

2

makedocs で HTML 化する

3

CI の docs job が deploy を担当する

4

GitHub Pages で公開する

この repo で見る場所

  • docs/make.jl が module とページ構成を決める
  • .github/workflows/ci.yml が build と公開の入口になる
  • docstring を増やすほど @docs の中身も増やせる

CIテスト: 手元の確認を自動化する

name: CI
on:
  push:
    branches:
      - main
    tags: ["*"]
  pull_request:
  workflow_dispatch:
jobs:
  test:
    ...
    steps:
      - uses: actions/checkout@v6
      - uses: julia-actions/setup-julia@v2
      - uses: julia-actions/cache@v3
      - uses: julia-actions/julia-buildpkg@v1
      - uses: julia-actions/julia-runtest@v1

CIテストをする理由

  • ビルドが通ることを保証できる
  • 壊れていれば早く気がつく

カバレッジ 1: 未検証箇所の発見

1

テストを実行する

2

通った行と通っていない行を記録する

3

未確認の分岐や異常系を見つける

4

必要なテストを追加する

coverage で分かること

  • どの分岐がまだ通っていないか
  • 異常系や端の条件を見落としていないか
  • examplestests の間に穴がないか

coverage で分からないこと

  • 物理的な妥当性そのもの
  • 比較条件が適切かどうか
  • 図の解釈が正しいかどうか

高スコア = 十分 ではない

カバレッジ 2: どこから増やすか

この repo で優先したい箇所

  • hhg_spectrum の異常入力
  • current_traces の長さ不一致
  • TODO から完成版へ変わる箇所
  • examples でしか触らない経路

増やし方の原則

  • 壊れたとき困る経路から埋める
  • 物理的に意味のある条件を優先する

将来の運用例

  • coverage レポートを CI で生成する
  • Codecov のような可視化サービスへ送る
  • PR ごとに「どこが未確認か」を見る

ただし

  • ここで大事なのは数値より運用の意味

論文投稿 1: JOSS が見る研究ソフトの観点

README

  • 何のソフトか
  • どう導入し、どう実行するか

docs / examples

  • どんな入力で何が出るか
  • 代表的な使い方を追えるか

tests / CI

  • 第三者が壊れていないと確認できるか
  • 継続的に検証されているか

software としての価値

  • 論文本文の補助ではなく
  • 再利用可能な道具として見てもらえるか

JOSS 的な見方

  • software そのものがレビュー対象
  • 他人が install して main example を追えることが重要
  • 説明、検証、実行例がそろって初めて再利用しやすい

この講義との接続

  • README を書く
  • テストと CI をそろえる
  • docs の層を分ける
  • その積み重ねが投稿準備になる

論文投稿 2: 授業 repo から投稿準備へ

授業後に足すもの

  • README を利用者向けに磨く
  • 詳細 docs を追加する
  • tests / CI を継続運用する
  • coverage の抜けを埋める
  • release と投稿用の説明文を整える

コメント

  • 再現できるだけではまだ半分
  • 他人が読めて、試せて、直せる形にして初めて共有できる

まとめ

できるようになったこと

  • Juliaを使って Jx(t), Jy(t) を FFT して HHG スペクトルへ変換できる
  • 動的対称性と選択則からスペクトル形状を説明できる

selection comparison from example

実装と物理の接続

  • hhg_spectrum が何をしているか読める
  • current_traces が何を集計しているか説明できる
  • 図から読めることと、理論的な説明を分けて話せる
  • README / docs / CI / テストが、再現結果の共有を支えると分かる

Appendix

Additional remark

  • 汎用コード化に向けて。各パラメータはハードコードしない。
  • ArgParse.jlを使うと便利。

パラメータ比較の原則

物理パラメータ

  • γ: コヒーレンスの減衰を強める
  • A0: 非線形性を強める
  • ω0: 横軸 n = ω/ω0 の基準を決める

数値条件

  • dt を変えると周波数軸が変わる
  • 入力長を変えると分解能が変わる
  • 比較では、同時に変える量を 1 個へ絞る