QSGW calculation in ecalj

Overview of QSGW calculation

(Clean version is given at https://ecalj.github.io/ecaljdoc/manual/README_tutorial)

• LDA calculations are performed with the program lmf. The initial setting file is ctrl.foobar ( foobar is user-defined). Before running lmf, it is necessary to run lmfa, which is a spherically symmetric atom calculation to determine the initial conditions for the electron density (a calculation that finishes instantaneously).

• A file sigm.foobar is the key for QSGW calculations. The file sigm.foobar contains $V_{\rm xc}^{\rm QSGW}-V_{\rm xc}^{\rm LDA}$. By adding this potential term to the usual LDA calculation performed by lmf, we can perform QSGW calculations.

• Thus the problem is how to generate $V_{\rm xc}^{\rm QSGW}({\bf r},{\bf r}')$. This is calculated from the self-energy $\Sigma({\bf r},{\bf r}',\omega)$, which is calculated in the GW approximation. Roughly speaking, we obtain $V_{\rm xc}^{\rm QSGW}({\bf r},{\bf r}')$ with fixing the omega-dependence in $\Sigma({\bf r},{\bf r}',\omega)$.

• Therefore, the calculation of $V_{\rm xc}^{\rm QSGW}$ is the major part of the QSGW cycle , and is calculated in a double-structure loop. That is, there is an inner loop of lmf, and an outer loop to calculates $V_{\rm xc}^{\rm QSGW}$ using the eigenfunctions given by lmf. This outer loop can be executed with a python script called gwsc (which runs fortran programs). The computational time for QSGW is much longer than that of LDA calculation. As a guideline, it takes about 10 hours for 20 atoms (depending on the number of electrons).

• GPU can accelerate it, https://arxiv.org/abs/2506.03477. Thus we can handle large systems.

• Since automation is being promoted, there is almost no need to tweak parameters that are difficult to understand (if something strange, let me know). I think it is probably the GW kinds of computation codes which is the easiest to use. See band database in QSGW at
  https://github.com/tkotani/DOSnpSupplement/blob/main/bandpng.md
  (this is supplement of https://arxiv.org/abs/2507.19189)

A mini database for computational tests (you can skip this)

At ecalj/MATERIALS/, you can run ./jobmaterials.py. Then you can perform LDA calculations for simple materials whose crystal structures are already contained in a mini database.

./job_materials.py

gives a help, showing a list of materials. Then

./job_materials.py Si

performs LDA calculation of Si at ecalj/MATERIALS/Si/.

• Key input files are
  ctrls.si,ctrl.si
  . See below. rst.si contains self-consistent electron density. Check iterations with save.si. The log file of lmf is in llmf. Not need to know all the console outputs.

• Before QSGW, it is better to confirm the LDA level calculations are fine. In order to do the confirmation, band plot is convenient.
  For band plot we need the symmetry line as syml.si which can
  be generated by

getsyml si

Then run

job_band si -np 8

results band plots in the gnuplot.

Calculation from scratch

step1. Get ctrls (crystal structure files)

At first, get POSCAR in VASP format. Then convert POSCAR to ctrls.foobar file. ctrls is the structure file equivalent to POSCAR.

How to get POSCAR?

You can get POSCAR from Material project (or make a query to MP and use pymatgen). It is easy to write a code to make such query in vscode with the help of copilot chat.

For example, POSCAR of mp-2534 GaAs is given as:

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

This is a poscar for ba2pdo2cl2:

POSCAR_ba2pdo2cl2
1.0
-2.06443 2.06443 8.40383 
2.06443 -2.06443 8.40383 
2.06443 2.06443 -8.40383 
Ba Pd O Cl 
2 1 2 2 
Cartesian
0.0 0.0 6.5153213224
0.0 0.0 10.2923386776
0.0 0.0 0.0
0.0 2.06443 0.0
2.06443 0.0 0.0
0.0 0.0 3.1625293056
0.0 0.0 13.6451306944

If you have cif and like to convert it to POSCAR, do
cif2cell foobar.cif -p vasp --vasp-cartesian --vasp-format=5.

Convert POSCAR to ctrls.foobar

> vasp2ctrl POSCAR

givs you the ctrls.POSCAR.vasp2ctrl. Rename this to ctrls.foobar such as ctrls.mp-2534. We can use foobar(=mp-2543 in this case) as a memo of the material. ecalj do not distinguish upper and lower cases.
Use viewvesta ctrls.foobar, which calls VESTA, to see the crystal structure.

Step2. Get ctrl.foobar

We genetate ctrl.foobar from ctrls.foobar by adding computational conditions. If you have ctrls.gaas, run

>ctrlgenM1.py gaas
>cp ctrlgenM1.ctrl.gaas ctrl.gaas

After here, we only need ctrl.foobar (ctrls.foobar is not read in the following process).

Edit ctrl.foobar if necessary. Explanations are embedded in ctrl.foobar (please let me know wrong descriptions). Possible points to rewrite in ctrl.foobar:

1. Number of k points (nk1,nk2,nk3).
2. nsp=2 if magnetic
3. SpinOrbitCoupling: so=0 (none), so=1 (LdotS), 2 (LzSz). nsp=2 is required for so=1,2. so=1 does not yet support QSGW. SOC axis can also be freely selected, but currently (0,0,1) default and (1,1,0) are supported (m_augmbl.f90). If you want to set SO=1 in QSGW, currently, run QSGW calculation with so=0 or so=2 to obtain ssig file, then set so=1
4. xcfun (choice of LDA exchange correlation term). Only =1:BH, =2:VWN, =103:PBE-GGA.
5. LDA+U settings (not explained yet).
6. ssig=1.0 (If you choose QSGW80, use ssig=0.8. Effective for QSGW calculations.
   $V^{\rm xc QSGW}-V^{\rm xc LDA}$ is stored in a file sigm.foobar. We add ssig $\times (V^{\rm xc QSGW}-V^{\rm xc LDA})$ to the potential in the lmf calculation as long as sigm.foobar file is available.

• lmchk --pr60 foobar allows you to check the recognized symmetries by lmf. Turning off --pr60 or reducing 60 will reduce the verbosity of output.

At this point, you can visually check the following check files.

• SiteInfo.chk
  MT radius Atomic positions
• PlatQlat.chk
  Primitive lattice vector (plat) Primitive reciprocal lattice vector (qlat)

Step 3. LDA band calculation

With ctrl.foobar.

>lmfa foobar

This ends instantaneously.

You can check the initial electron configuration with lmfa gaas|grep conf (there are no side effects even if lmfa is repeated). This calculates spherically symmetric atoms and generates the files required for lmf below. The initial condition of electron density for lmf is given as the superposition of spherically symmetric atomic densities. In addition, calculations are performed with a fixed logarithmic derivative of the radial wave function at the MT sphere edge (in default setting). The derivatives are contained in atmpnu.* files. So, atmpnu.* are needed for lmf.

• After that, run the LDA calculation with

mpirun -np 8 lmf foobarx

foobarx can be either ctrl.foobar or foobar. For metals, increase the k-points. (In the case of ba2pdo2cl2, with the default -np 8, one iteration takes about 60 seconds. It converges within 20 iterations.)

(This LDA step is included in the gwsc script for QSGW below. but it might be better to do things one by one at first).

Step 4. Create k-path and BZ for band plot

After the calculation converges, it might be necessary to make a band plot with job_band command explain later on. The normality of the calculation of bands can be confirmed by the band plot (for magnetic systems, check the total magnetic moment and the magnetic moment for each site).

Before job_band, run getsyml gaas. Install any missing packages with pip. It is on spglib by Togo and seekpath. After finished, view BZ.html. It shows the k-path in the BZ ashow show below for ba2pdo2cl2.
It is an interactive figure written with plotly, so you can read the coordinate values.

• Samples of BZ.html by getsyml are seen at https://ecalj.sakura.ne.jp/BZ/.
  

Step 5. band plot

（this is a case for ba2pdo2cl2 ）

>job_band ctrl.ba2pdo2cl2 -np 8

A gnuplot script can be created. Edit it if necessary. If you edit syml.ba2pdo2cl2 before job_band, you can adjust the symmetry line and mesh size.
job_pdos calculates PDOS, job_tdos calculates total DOS, and job_fermisurface draws the Fermi surface.
job_fermisurface can be used to draw the shape of the CBM bottom as
ellipsoid of Si.

• The following picture is the LDA bands for the default calculation of ba2pdo2cl2 (the names of the symmetric points can be confirmed with BZ.html. In addition, look into syml.foobar). 0 eV is the Fermi energy. Since this is metallic, we see no band gap.
  .

• The defaults are fine except for the k-mesh setting. For example, it is better to increase the k mesh for Fe. In general, for semiconductors, 4 4 4 for Si is a reasonable level, 6 6 6 is a level that can be used for a paper, and 8 8 8 is a level for checking accuracy. For metals such as Fe, 8 8 8 is a reasonable level.

Here we are talking about band energies.

• In ecalj, the k mesh for lmf (ctrl) and the k mesh for GW (n1n2n3 specified in GWinput) can be different. The former has affected little on computational time, but the latter has a large effect (thus we want to reduce n1n2n3 in GWinput).

• In ecalj's band plot mode, theoretically degenerated bands because of symmetry at the BZ edge are not degenerated. This is because there are limited numbers of APW basis functions, so run the band plot with pwemax=4, etc. (Temporary solution: We want to automate it).
  

Step.6 QSGW calculation. Example for ba2pdo2cl2

This is a case of ba2pdo2cl2.

Run

mkGWinput ba2pdo2cl2

to generate GWinput.tmp, which is a setting file for QSGW.
After copying this to GWinput, you may need to edit GWinput.
Minimum thing to edit is the number of k points for the self energy (n1n2n3).
Compared with k points in ctrl (nk1,nk2,nk3), we use small numbers.
(We usually use 1/2 or 2/3 of k points given in ctrl as nk1,nk2,nk3).

There are another setting in GWinput. However, we usually do not need to touch things except n1n2n3 if you treat non-magnetic semiconductors.

Then you can run QSGW calculation with

gwsc -np 32 1 ba2pdo2cl2

. Here 1 means the number of QSGW iterations. QSGW iteration is quite time-consuming. gwsc gives minimum help (we need to explain options elsewhere).

The iteration is kept in rst.foobar:electron density, sigm.*:vxcqsgw.
(Remove these files in addition to *run files/directories if you like to start from the beginning).

• It requires 53 minutes to run one iteration of QSGW.

• job_band ba2pdo2cl2 -np 32 gives the following picture. QSGW one-shot changes band structure around Ef from that in LDA.But still metallic, no band gaps.
  

• To continue QSGW iteration, run

gwsc -np 32 nx ba2pdo2cl2

Since you did 1 already. You will have the results of 1+nx QSGQ iteration.

• when we run 8 iterations as for ba2pdo2cl2, we had band gap 2.1 eV. We saw band gap after 4th iteration.
  

• For comparison with experiments, we recommend to use ssig=0.8
  (set in ctrl.foobar file), which is called as QSGW80.

• Spin-orbit coupling. After you obtain sigm.foobar, you set SO=1 (LdotS scheme) and run lmf. Then you can include effect of SOC. Sinc SO=1
  is not implemented in the whole gwsc cycle, we have to include SOC just at the end step (We include SOC after we fix VxcQSGW).

• If you run

gwsc -np 32 5 ba2pdo2cl2 -vssig=0.8

, this overide ssig, which is defined in ctrl.ba2pdo2cl2, in lmf calculations. (Check it in save.ba2pdo2cl2)

• Example of QSGW for KTaO3 (perovskite,mp-3614）
  

xxxxxxxxxxxxxxx under construction xxxxxxxxxxxxxxxxxxxxxx

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

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,自己エネルギー（寿命）、スピンゆらぎなどを求めることができる。