ecaljによるQSGW計算

QSGW計算の概要

• LDAの計算はプログラムlmfで行う。初期設定ファイルはctrl.gaasに書く(.gaas部分はユーザー定義）。lmf実行前に、電子密度初期条件を作るために球対称原子計算lmfaを行っておく必要がある(一瞬で終わる計算）。QSGW計算ではsigmファイルを計算する必要がある。sigmファイルには、$V_{\rm xc}^{\rm QSGW}-V_{\rm xc}^{\rm LDA}$が入っている。これを通常のlmf計算に加えることでQSGW計算が行えることになる。
• 問題となるのは、$V_{\rm xc}^{\rm QSGW}$を生成するところ。これはGW計算における自己エネルギーに対応するものである。自己エネルギーは振動数（時間）依存である。一方、$V_{\rm xc}^{\rm QSGW}$は非局所的ではあるが時間依存性はない。大雑把にいって、$V_{\rm xc}^{\rm QSGW}$は自己エネルギーから時間依存性を適当に消したものである。
• なので$V_{\rm xc}^{\rm QSGW}$を計算するのがQSGWの大部分であり二重構造のループで計算することになる。すなわち、lmf内部のループと、その結果（波動関数、固有値）を使って$V_{\rm xc}^{\rm QSGW}$を計算する外側のループである。この外側のループは,gwscというpythonスクリプト(fortran プログラムを連鎖的に呼び出す）で実行できる。QSGWの計算時間はlda計算よりかなり長くかかる。目安としては、20原子とかで10時間程度かな（電子数による）。ただGPU高速化もしているので大きい系でもそれなりにはできたりもする。自動化はすすめているので、よくわからないパラメータをいじったりはほとんど必要ない。たぶん世にあるGW計算のコードの中では最も取り扱いが簡単だと思う。

インストール

物性研マシンについては別のインストールマニュアルあり。

パッケージは,https://github.com/tkotani/ecalj/
のmasterが現状でベストなバージョン(2025-1月)なのでこれをクローンする。httpsなら

git clone https://github.com/tkotani/ecalj.git

• 必要なツール等は、たぶんbash make cmake gfortran(or ifort or nvfortran,intel-mkl,VESTAとpip系のツール

pip install mp-api spglib,seekpath,,numpy,ploty,gnuplot,cif2cell```

など。pipでローカルインストールしておけばよい。python versionが古いとmp-apiなどで文句を言われる。足りないと言われたときにインストールするのでもよい。

ecaljをインストールするには、ecalj/で以下のどれかを行う。

FC=ifort ./InstallAll
FC=gfortran ./InstallAll
FC=nvfortran ./InstallAll

fortran sourceコードからのコンパイル、リンクの後、binへのコピーが行われ、さらにecalj/SRC/TestInstallへ移動してインストールテスト（事前実行してあるものとの計算比較）が行われる。終了するのに10分ぐらいはかかるかもしれない(最新のubuntuのノートPC,gfortranで５分ぐらい）。

• スクリプト./InstallAllは事前に目視確認しておいたほうがよい。インストール場所や手順などが簡潔に書いてある。その内容は「バイナリをおくためのディレクトリBINDIRを作る(~/binがデフォルト）。getsymlやviewvesta,vasp2ctrlなどのソフトリンクを作る。フォートランをコンパイルして複数のバイナリを作る。必要なスクリプトやバイナリをBINDIRにコピーする。最後にテストを実行。」

インストールテストが終わって、すべてOKと出ることを確認。以下のような感じで終了するはず。

PASSED! TEST 1 out.lmf.copt
PASSED! TEST 1 out.lmf.te
...(中略)...
PASSED! nio_gwsc/QPU
PASSED! nio_gwsc/log.nio
PASSED! fe_gwsc/QPU
PASSED! fe_gwsc/QPD
PASSED! fe_gwsc/log.fe
PASSED! ni_crpa/Screening_W-v.h
PASSED! ni_crpa/Screening_W-v_crpa.h
PASSED! srvo3_crpa/Screening_W-v.h
PASSED! srvo3_crpa/Screening_W-v_crpa.h
OK! ALL PASSED ===
   See work/summary.txt

real	4m49.712s
user	19m38.952s
sys	3m38.583s

このように表示されていれば正常インストールはできているはず。

• テストの個別実行は、ecalj/SRC/TestInstall/で
  ./testecalj.py si_gwsc
  などとして実行する。個別テストごとにsi_gwscなどのディレクトリがある。個別テストでは、ecalj/SRC/TestInstall/work/si_gwscを作り、計算を行い、ecalj/SRC/TestInstall/si_gwscにある計算結果と照合する仕組みである。テストはtestecalj.pyに書いてある。単純な仕組みなので、これを参考に自分のテストも追加できる。

• インストールが失敗しクリーンするには、ecalj/SRC/exec/でCMakeFiles/とCMakeCache を消すこと。またmodファイルが残ってトラブルこともあり、SRC以下を検索して消したほうがいい場合もある(modが残るのはgfortranとifortを同居させたりのときにトラブル可能性があるーあまりテストできていない）。git logでバージョン確認しておく。

• たまにパラレルメイク(make -j)の効かないマシンもある。以下のようにしてmake終了後、インストール操作を再開すればインストールテストが起動される。

cd ecalj/SRC/exec/
rm -rf CMakeFiles CMakeCache.txt 
FC=nvfortran cmake .
make 

(コード開発時、makeを繰り返すとsequential makeになってしまうときがある-CMakeFiles,CMakeCacheなどを消してやり直す）。

• テストのためのサンプル
  ecalj/MATERIALS/で./jobmaterials.pyを実行するとあらかじめ用意されたテスト用構造データベースを用いた計算が実行できる。
  ./job_materials.pyを実行するとヘルプがでる。テストにある物質名のリストが表示される。./job_materials.py Siを実行するとSiについてのLDA計算が行われる。ecalj/MATERIALS/Si/で行われる。ctrls.si,ctrl.siの説明は以下。またrst.siに計算の電子密度が格納されている。イテレーションの様相はsave.siで確認。lmfのログはllmfにあるのでサクッと見ておく。LDAレベルの計算が正常にできているかを確認しておく。そのためには以下で説明しているようにgetsyml siでsyml.siを生成後（ブラウザが立ち上がりBZが表示される)、job_band si -np 8でバンドプロットを確認。以下では、この計算を手動で再現したのちQSGW計算を行うことになる。

** インストールトラブルがあれば報告してもらえるとありがたいです **

ctrlsファイル（結晶構造ファイル）の作成

POSCARあるいはcifなどからctrls.foobarを生成する。

POSCARを得る

POSCARを得るにはMPからダウンロードするとかMPに直接クエリーをかけpymatgenなどで変換すれば良い。pymatgenを勉強するよりも、vscode github copilotに聞けば、mp-foobarをダウンロードしてPOSCARに書き出すpythonスクリプトぐらいは書いてくれるはず。もっと高度なqueryも書いてくれるかも。
たとえば、mp-2534のGaAsならでPOSCARは

Ga1 As1
1.0
3.5212530000000002 0.0000000000000000 2.0329969999999999
1.1737510000000000 3.3198690000000002 2.0329969999999999
0.0000000000000000 0.0000000000000000 4.0659929999999997
Ga As
1 1
direct
0.0000000000000000 0.0000000000000000 0.0000000000000000 Ga
0.2500000000000000 0.2500000000000000 0.2500000000000000 As

cifならあらかじめPOSCARに変換する必要がある。まずcifをCrystal Open Database(COD)からダウンロードする。そのあと
cif2cell foobar.cif -p vasp --vasp-cartesian --vasp-format=5で、POSCARを得ておく.

ctrlsへの変換

vasp2ctrl POSCAR

として、ctrls.POSCAR.vasp2ctrlを得ることができる。これをrenameしてctrls.foobarとする。たとえばctrls.mp-2534など。ctrls.gaasでもよい。ドット以下の拡張子で物質名のメモとするー大文字小文字を区別しない。
viewvesta ctrls.foobarでctrlsを目視確認できる。POSCARを見たのと同じ構造のはず。

ctrlファイルを得る

ctrlsをctrlに変換する。ctrls.gaasであれば、

ctrlgenM1.py gaas

を実行後、

cp ctrlgenM1.ctrl.gaas ctrl.gaas

を実行して、ctrl.gaasを得る。以下の計算で必要なのはctrlファイルのみ。

• ctrlの編集が必要なら編集する。ctrlには説明が書き込まれている(不具合指摘してください）。ctrlにおいて書き換える可能性のある点：

1. ｋ点の数(nk1,nk2,nk3).
2. 磁性のある場合はnsp=2
3. SpinOrbitCoupling:so=0(なし）、so=1(LdotS)、2(LzSz)。so=1,2ではnsp=2が必須。so=1はまだQSGWに対応していない。SOCの軸も自由に選べるがいまは(0,0,1)デフォルト、と(1,1,0)に対応(m_augmbl.f90).QSGWでSO=1にしたい場合、現状ではQSGW計算をso=0,もしくはso=2で実行してssigファイルを得た後、so=1
4. xcfun（LDA交換相関項のチョイス）.BH,VWN,PBE-GGAのみ。
5. LDA+Uの設定(説明まだ）.
6. ssig=1.0 (QSGW80を選ぶならssig=0.8とする。QSGW計算で有効。
   $V^{\rm xc QSGW}-V^{\rm xc LDA}$はssigmというファイルに格納されるが、これをどれだけの割合で通常のLDA計算に加えるかがssig

lmchk --pr60 foobarにより、認識される対称性などがチェックできる(要説明）。--pr60を消す、もしくは60を小さくすると表示量が小さくなる。

• この時点で以下のチェックファイルを目視確認できる。
• SiteInfo.chk
  MT半径　原子位置
• PlatQlat.chk
  基本格子ベクトル(plat)　基本逆格子ベクトル(qlat)

バンドプロットのためのk-path,BZ作成

計算収束後、バンドプロットする（以下のjob_bandコマンド）のに必要になる。
バンドプロットにより計算の正常性がおおむね確認できる（磁性のある場合、加えて、全磁気モーメント、サイトごとの磁気モーメントを確認すること）。
getsyml gaasを実行.足りないパッケージはpipインストールしておくこと。主にはspglib,seekpathを使っている。ブラウザが立ち上がりk-pathが表示されるはず。以下の図はba2pdo2cl2の場合。plotlyで書いているインタラクティブ図なので座標値も読める。

• getsymlによるBZ.html (座標値を読んでPlatQlat.chkと照らし合わせて確認する）。https://ecalj.sakura.ne.jp/BZ/ にサンプルをおいた。

LDAのバンド計算：

• まずlmfa foobarを実行.
  lmfa gaas|grep confで初期電子配置が確認できる(lmfaは繰り返してもサイド・エフェクトはない).これにより球対称原子を計算し、以下のlmfに必要なファイルを生成する。バレンスの電子密度は球対称原子の重ね合わせからスタート。あと、デフォルトではMT球端での動径波動関数の対数微分を固定した計算を行うのでそれの入ったatmpnuというファイルも大事。
• その後、time mpirun -np 8 lmf foobarを実行。(time は時間が知りたいとき）。foobarはctrl.gaasでもgaasでもよい。tab補完をつかうにはctrl.gaasのほうが入力しやすい。金属のときにはk点をふやすこと。(ba2pdo2cl2の場合、デフォルト設定だと-np 8でイテレーション1回に60秒ぐらい。20回以内で収束する。)

バンドプロットなどを行ってみる。

（以下、ba2pdo2cl2の場合）

job_band ctrl.ba2pdo2cl2 -np 8

gnuplotスクリプトができる。それをエディトして描画は調整。syml.ba2pdo2cl2を事前に編集してjob_band実行すれば、シンメトリラインや刻み幅を調整できる。
job_pdosはPDOS計算、job_tdosはtotal DOS計算、job_fermisurfaceはフェルミ面描画。
job_fermisurfaceでCBMボトムの形状を描いたりもできる。

ba2pdo2cl2のデフォルト計算でのLDAのバンド（対称点の名前はgnuplotファイル*.gltとgetsymlでできたBZ.htmlを用いてfirefox BZ.htmlで確認）。0eVがフェルミエネルギーで金属になっている。バンドギャップがない。
。

• ｋ点数の設定以外はデフォルトで良い。k点数はたとえばFeだと増やしたほうが良い。一般に半導体だと、Siで4 4 4:まあまあレベル、6 6 6:論文にできるレベル、8 8 8:精度確認のレベル、ぐらいに思えば良い。Feなどの金属であれば、8 8 8がまあまあレベル、になる。
  ここではバンドについて言ってる（全エネルギー計算ほどにはシビアではない）。

• ecaljではlmfのk点数（ctrl)とGWのk点数（GWinput,sigmを求める点数）を違うように取れる。前者は大きくとっても計算時間に影響でないが、後者は大きく影響する（その意味でできるだけへらしたいところ）

• ecaljのバントプロットモードに置いてBZ端で縮退すべきところが縮退しない、というふうに見えることがある。これはAPW基底関数が少ないためであるので(バンドプロットモードでは平面波のGをq=0で決めてしまうため）、pwemax=4などとしてバンドプロットを実行すること。
  (とりあえずの対症療法：自動化したい）。

QSGW計算の例 ba2pdo2cl2の場合

まずmkGWinput ba2pdo2cl2を実行するとGWinput.tmpができる。これはGW計算部分(sigmを求める計算）を制御する設定の書かれたファイルである。
これをGWinputにコピーした後

gwsc -np 32 1 ba2pdo2cl2

をまずやってみる。１がQSGW計算の繰り返し回数。バンドプロットは同じやり方でできる。gwscと引数なしで打つとヘルプが出る。GWinputではk点(n1n2n3)はいじったほうがいい。あといくつか設定があるけど半導体ならデフォルトで問題はないと思う。

gwsc -np 32 1 ba2pdo2cl2

をゼロから再実行するときは,sigm*, (NTQXX:最新版にはない）, run, RUNを消して,rst.lda.ba2pdo2cl2をrst.ba2pdo2cl2にコピーして実行する。
もし、これらのファイルを消さずに実行すると、継続実行になる。

• QSGW one-shot (gwsc -np 32 1 ba2pdo2cl2を実行後、job_band ba2pdo2cl2 -np 32を実行.qsubでサブミットしたほうがいい。所要時間53分だった。）
  
  １回目のQSGWではフェルミエネルギー近傍でバンドが分離しているが、メタルのまま）。
  GWinput内で設定するｋ点数は注意。目安として、Siで6 6 6以上だと論文出版できる精度。4 4 4が大体の結果を得るための精度。2 2 2はちょっと荒い計算。あくまで目安。これを基準に原子あたり何個かを考えて増減する。QSGWの固有値（バンドのエネルギー）の数値精度の目安は0.1eV程度と思っておくのがよい。
  ctrlではバンド計算でのｋ点数を指定、GWinputではGW計算でのｋ点数を指定する。
  QSGW計算では、GWinputで書かれたｋ点数を用いて、
  $V_{\rm xc}^{\rm QSGW}-V_{\rm xc}^{\rm LDA}$を計算する。QSGWのイテレーションのサイクルではこれを付加してバンド計算を収束するまで行うことになる。

• ba2pdo2cl2についてQSGWのイテレーションを8回回すとバンドギャップは2.1eVだった。およそ収束しているだろう。イテレーション４回目ぐらいからバンドギャップが見られる。
  
  ただし、実験との比較にはssig=0.8で行う(QSGW80)。最後の段階でもよい。SO=1として最後のバンド計算をまわしたほうがよい。

• 磁性のある金属が一般には収束が一番厄介。イテレーションが数十回必要な場合もあるし、収束サインがでていても少しずつモーメントが動いていったりもする。

* QSGW計算例：KTaO3 (ペロフスカイト,mp-3614）。

こういう計算は、>job_mp mp-3614
と打つだけで可能である（投稿準備中）。
デフォルト計算。
ワニエ化しなくとも任意のｋ点で一体tハミルトニアンを得ることができるので線形応答などの計算も自動化できる。

自動モデル化も動くようにはなってきている（未発表）。

* その他

• 物理量を求めるいろいろな機能があるが整備中。

ecalj/Samples
├── LaGaO3_relax　構造緩和サンプル（LDA/GGAレベル）
├── BOLZTRAP  ボルツトラップ用入力生成
├── CuepsPP0  誘電関数
├── FermiSurface　フェルミ面
├── MLOsamples　自動モデル化
├── MgO_PROCAR　fatband 
├── README.md
├── SOCAXIS     SOCを入れるときの軸の向きの設定
├── TETRAHEDRON_HomoGas　一様電子ガスの誘電関数(Lindhard関数）のサンプル
├── UUmatSOC　<u_kn |u_k+b m>を求める。
└── mass_fit_test 有効質量を求めるサンプル
あと、pdos,自己エネルギー（寿命）、スピンゆらぎなどを求めることができる。

• GPU並列もかなり可能（改善の余地はあるのだが）.単純なGPU可を心がけている。

• 原因不明だがmake -j のパラレルメイクが効かなくなるときがある。./CleanAllで掃除すると復活する。

• コードはfortranだがmoduleをシングルトンとして用いるようなやり方で書かれており、classで書く人にもなじみやすいはず。また単純なメイン関数以外はライブラリ化している。

• 将来的には、上位側からpython化したい。exec/lmf.pyは実験的に書いたものだけど、main/lmf.f90を置き換えたい。ただまだpythonバージョン次第で動かない場合もあったりで移行する気になれないので一旦断念してる。pytorchとか勉強したけど、テンソル演算が充実してて行列などの添字を直接にいじらずに（loopを書かないで）計算していくからコードの階層が深くならない。将来的には数値計算はpytorchになっていくのかな？とも思えます。

---