注意 この記事は永遠に完成しないと思います。
はじめに
固体物理の第一原理計算コードのうち、Quantum Espresso(以下、QE)のコードを読んでいく。
どこで何をしているかを把握して、理屈と実装をつなげていくことが目標。
なお、必要な前提知識としては、固体の電子状態計算に関する本である「物質の電子状態」が有名。
ここでは、Quantum EspressoのVersion 6.1(qe-6.1)を対象とする。1
コードリーディング
コールグラフ
プログラムから簡単なコードグラフを作成するスクリプト2を作成し、pwscfから呼ばれる関数を書いた。
コールグラフ(簡約版)
なお、vimで、set foldmethod=markerとすると折りたためるように、マーカーを挿入してある。
より深いレベルまで出力したものは、ここにある。
主なコードは、主にソースツリーのPW/src/の中、もしくは共通性の高いもの（ポテンシャル関連や、交換・相関汎関数のパラメータなど）はModules/の中に置かれている。
- pwscf.f90: エントリポイントはここ。
- run_pwscf.f90 : ここで、入力ファイルを読んだり、MPIのセットアップをしたりする。そして、構造最適化など、scf計算をする場合は、electrons()でscf計算を行い、その後、forces()という関数をつかって力の計算を行う。DOSやバンドなどを求めるためのnon-scf計算の場合は、non_scf()をコールする。
変数の一覧
全体で使いまわされる変数リスト（作成中）はここ。
主要なデータ構造
波動関数とポテンシャルをどういうデータ構造で格納しているかについて。詳細は、PW/src/scf_mod.f90の中にある構造体をみる。
SCF計算
electrons.f90というコードがSCFループを記述している。
ここにはおもに、electrons()とelectrons_scf()というサブルーチンが定義されており、前者が後者をコールする形になっている。
electrons_scf()のなかで、SCFループを回す。
electrons()は何をするかというと、electrons_scf() で得られた密度に対して、Hartree Fock交換積分を計算する。ただし、この必要があるのはハイブリッド汎関数を対象とする時のみなので、そうでない汎関数(LDAやGGA)の時は、事実上、electrons_scf()だけで終了する。
SCFループの中身： electrons_scf()の流れ
electrons_scf()の中でSCFをとく。SCFループの中で呼ばれる関数は主に以下のようになる。
- c_bands() : 実際にハミルトニアン（フォック行列）を各k点に対して求め、対角化を行う。
- sum_band(): c_bands()で得られた波動関数から、密度を求める。
- mix_rho(): Fock行列を作るときに使用した入力密度(ρ_in)と、sum_band()で計算して得られた新しい密度(ρ_out)を混ぜる。
- v_of_rho(): これは、scfが収束していないときのみ呼ばれる。次回のループで使われるc_bands()で使うための、入力ポテンシャル(Vin)を、さっき得た密度から計算する。
なお、この途中で、エネルギーの計算を行う。収束のチェックは、Harris-Weinert-Foukles Energyと実際にKohn Sham方程式から得られるエネルギーの値の差が十分に小さくなったかどうかを見ている。
Harris-Weinert-Foukles Energyについては、物質の電子状態の「9.2 全エネルギー汎関数」の節のなかの「密度の陽関数」という項目を参照。簡単にいえば、
-
Kohn-Sham方程式から得られるエネルギー: 入力の密度から、各k点ごとにフォック行列を求めて、固有値（エネルギー）を求める。このときに同時に得られる固有ベクトルがその時の電子密度を記述するのだが、この密度を汎関数の引数とする。これは、各バンドのエネルギーを計算する際のエネルギーは前の密度を使って求めているものであるため、自己無撞着が達成された時以外は実際上の意味を持たない。（と私は認識しているのだが、ちょっとあやしい）
-
Harris-Foukles Energyは、入力密度のみから得られるエネルギーである。KSエネルギーでは、フォック行列の固有ベクトルから計算される密度を交換・相関汎関数に渡してエネルギーを求めたが、こちらはその代わりに入力密度をそのまま使う。
-
これらは、自己無撞着が達成されるとほぼ同一のエネルギーを示すことになるはずであり、そのことを利用して、自己無撞着への達成具合の評価としている。
これらを次々に読んでいけばよい。
c_bands()
c_bands.f90というファイルには、c_bands()とc_bands_nscf()という二つのサブルーチンがある。後者は名前の通り、バンド構造やDOSを求めるためのnon-scf計算の時に呼ばれる。
ここでは、c_bands()について読む。
g2_kin()で、運動エネルギーの計算
sum_band()
mix_rho()
電荷密度から、次のSCFの初期電荷密度を計算する
v_of_rho()
電子密度からポテンシャルを計算する
最後に
また書きます。