dHvA計算

FPLOには、de Haas-van Alphen（dHvA）効果における振動数Fを計算するツールが実装されている。本ツールは、オンサーガーの量子化条件（Onsager relation）F = (¥hbar c) / (2 ¥pi e) S_mを用いて、磁場に垂直な面におけるフェルミ面の極値断面積S_mを数値的に求めることで、dHvA振動数Fを計算している。また、FPLOでは単にE = E_Fの等エネルギー面であるフェルミ面だけでなく、一般の等エネルギー面（iso surface）に対してdHvA振動数を計算することができる。

FPLOにおけるdHvA計算は、以下の2ステップに分割されている。
(1) iso-surface stage: 収束したSCF計算の結果（=.dens）から、高精度のiso surfaceを求め、+iso_b..._p..._spin...として出力する。「Iso Surface (R)un」を選択してfploを実行することで、このステップが実行される。
(2) fdhva: iso surfaceを読み込み、様々な磁場に垂直な平面で切って極値断面積を探し、Onsarger関係式からdHvA振動数Fを計算して、+area_vs_angle_...や+mass_vs_angle...として出力する。このステップはfploではなくfdhvaで行われる（ただしパラメーターはiso-surface stageと同じfeditのsub-menuで設定する）。

ステップ(1)および(2)におけるパラメーターは、feditの「(D)HvA」sub-menuにて設定できる。各ステップの概要と、重要なパラメーターを以下に記す。

(1) iso-surface stage
Brillouin zone（BZ）のうちほぼ全ての空間はiso surfaceではない。そのため単にBZ全体を探索するのでは効率が悪い。そのため、以下のような2段階でiso surface探索を行う。まず、Initial subdivision: (nk1, nk2, nk3)でメッシュを切り、粗くiso-surfaceを探索する。続いて、iso-surfaceを発見したk点をさらにBisectionだけ探索して、iso-surfaceを滑らかにする。したがって、Initial subdivision探索で発見できなかったiso-surface sheetは、Bisectionをどんなに増加させても発見できない。そのため、まずはBisection = 0にして、Initial subdivition探索で全シートが感知できているか確認すべきである。iso-surfaceはXFPLOで可視化できる；xfplo +iso_*。これで全シートが感知できていれば、Initial subdivisionは十分である。続いて有限のBisectionを設定して、再度fploを実行する。

(2) fdhva
磁場の向きを変えながら極値断面積を探索して、dHvA振動数を計算できる。計算したdHvA振動数はxfbp area_vs_angle.cmdでプロットできる。ある磁場に対して、それに垂直な面を複数用意してiso-surfaceの断面積を計算し、極値を求める。そのため、カット面の数（No of planes）に対する収束性も確認すべきである。また、バンド計算におけるk-pathのように、磁場方向の変化のさせ方を設定できる。代表的な方向とその間の点数（An(G)le subdiv）で定義される。Angle subdivは、単にプロットの滑らかさだけではなく、branchの同定に非常に重要なパラメーターである。磁場の方向を若干変化させた時の、断面積の追跡には、それなりに多くのAn(G)le subdivが必要である。

推奨ワークフロー

上述のような計算原理を踏まえると、以下のような手順を踏んでパラメーターに対する収束性を確認すべきである。

1. SCF計算を収束させる。

2. xfploで高精度のフェルミ面を描写し、存在するシートを確認する。

3. まずはBisections = 0に固定して、全シートを捕捉できる程度に細かいInitial subdivisionを決める。iso-surfaceの描写はxfplo +iso_*で可能。

4. 続いて有限のBisectionsを設定してiso-surfaceを再計算し、xfploでiso-surfaceの滑らかさを確認。一旦、目視で十分滑らかとなるまでBisectionsを上げる。

5. 続いてfdhvaを実行しdHvA振動数を計算する。

6. Bisections, No of planes, Angle subdivを変化させながら、dHvA振動数の収束性を確認する。

1. SCF計算

ここでは、dHvA計算の題材としてNbAl3​を扱う。NbAl3​は、Dirac半金属の候補として注目されている金属間化合物である。Dirac点の詳細を議論するにはスピン軌道相互作用を含むfully relativistic計算が必要となるが、フェルミ面形状には大きな影響を及ぼさないため（時間があれば両者を比較するとよい）、ここではscalar relativistic計算を用いてdHvA振動数を計算する。結晶構造はVESTAファイルを参照。12*12*12 k-mesh, PBE-GGAを持ちいれば十分だろう。これでSCF計算を収束させる。

2. フェルミ面の描写

収束したscalar relavitistic GGA計算の結果からXFPLOを用いて描写したフェルミ面を下に記す。電子面βとホール面γの2枚のシートから成るフェルミ面を示している。[001]方向に磁場を印加した場合の極値断面積は、β面が極小（neck）β1と極大（belly）β2を示し、γ面はただ1つの極大γを示す。これらに相当する3つのブランチを持つdHvA振動数が想定される。磁場を[010]方向に傾けてゆくと、β1とβ2はやがて一致して大きな断面積を示す。逆にγはky = 0のneckとky != 0の小さなneckの2つに分裂することが想像できる。

3. iso-surfaceの計算：Initial subdivitionの決定

まずはデフォルトのdHvA設定（Inititial subdivition = 12, 12, 12, Bisections = 0）で、Iso Surface (R)unをオンにしてfploを実行し、iso-surface（フェルミ面）を計算してみる。計算されたフェルミ面を描写；xfplo +iso_* すると、以下のようになる。非常に粗いが、全てのシートが検知できていることが分かるので、Inititial subdivition = 12, 12, 12を採用しよう。

4. iso-surfaceの計算：Bisectionsの決定

feditのdHvA sub-menuでBisections = 2に設定し、再びfploを実行してiso-surfaceを計算・描画する。γ面のジャギーさが目立つものの、フェルミ面の概形は捉えられているので、まずはこの状態でdHvA計算に進んでみよう。（なお、XFPLOではフェルミ面描画時の色がデフォルトでフェルミ速度を表している。これは、計算されたフェルミ面の滑らかさを目視で確認する際にも役立つ。後述のように、この段階でフェルミ速度が自然に連続的になる程度にはBisectionsを上げるのが良い。）

5. fdvaの実行

Bisections = 2で求めたフェルミ面を使って、極値断面積およびdHvA振動数の計算を行う。
fdhva > fdva.out （念のため標準出力をファイルにパイプしておく。）
数秒で終わるので、dHvA振動数をプロットしてみる。
xfbp area_vs_angle.cmd
孤立した点や分裂してしまったブランチが見えるが、[001]方向のdHvA振動数を観察してみよう。この周りでは3つのブランチが得られている。飛び抜けて大きな振動数は、飛び抜けて大きな極値断面積であるγからの寄与と予想できる。次いでβ2、β1からの寄与と予想できる。磁場を[010]に傾けてゆくと、β1（黄線）とβ2（青線）は黒線に縮退して、逆にγ（緑線）は水色線と橙線の2つのブランチに分裂している様子が確認できる。このブランチの様子はフェルミ面から予想していた振る舞いとコンシステントである。しかし、β面由来のブランチが10° < θ < 20°で分裂していたり、γ由来のブランチが同一視されていなかったり、孤立した点が複数見られる。以降では、Bisections, angles, planesに対する収束性を見て、これらが改善されるか確認する。

6. パラメーターに対する収束性の確認

Bisectionsを変更する場合、iso-surface計算からやり直す必要がある。+isoergcacheの使い回しが出来るとはいえFPLOの再計算となるので、計算資源削減のためにディレクトリを分割するのが吉であろう。ここではBisections = 2, 4を試すので、Bisections=4用にディレクトリごとコピーしておく。まず、Bisectionsを変更するだけでも、孤立した点やブランチ分岐がだいぶ改善することが分かる。さらに、planesを200に、anglesを60まで増加させると、予想されるブランチのみの綺麗な結果となる（右下図）。やはり、+iso_b*を描画した時にフェルミ速度が滑らかに見える程度にBisectionを設定するのが肝心と考えられる。