GaAs: Construction of MLWFs with SOC and calculation of the shift current

In this example we describe how to obtain Wannier functions for GaAs with spin-orbit coupling. After generating the Wannier functions we use them to obtain the interpolated band structure and to calculate the shift current. This tutorial example follows closely the example on shift current calculation in the Wannier90 tutorial.

We use the following inpgen.inp file as input for the input generator:

 GaAs with SOC
 &input  film =f /

&lattice latsys='cF' a= 10.734/

           2
31 -0.125 -0.125 -0.125
33    0.125 0.125 0.125
&soc 0.0 0.0 /

Run inpgen2 -f inpgen.inp

Run fleur_MPI in order to converge the charge density. Check that the charge density is converged.

After converging the charge density we start with the construction of MLWFs. For this purpose, generate a new subdirectory Wann in your working directory. Copy inp.xml, cdn.hdf and sym.xml into the Wann subdirectory. Change into the Wann subdirectory and obtain from the inp.xml file information on the local orbitals. In this example, the input generator has defined the following local orbitals:

<lo type="SCLO" l="2" n="3" eDeriv="0"/>

for the Ga-1 species and

<lo type="SCLO" l="2" n="3" eDeriv="0"/>

for the As-1 species. The local orbitals are lower in energy than the valence states. We will therefore skip them. There are 20 local orbitals in total. The index of the first valence band is therefore 21. Since we would like to prepare MLWS for the calculation of the shift current, it is not sufficient to consider only the 8 MLWFs corresponding to the bonding states below the gap. We need the antibonding states above the gap as well. Therefore we will disentangle 16 MLWFs from 36 bands. The index of the first band that we consider in the MLWF calculation is 21. The index of the last band that we consider in the MLWF calculation is 20+36=56. Therefore we add the section wannier to the output section of inp.xml with minSpinUp="21" and maxSpinUp="56" in the bandSelection:

   <output dos="F" band="F" vacdos="F" slice="F" mcd="F" wannier="T">
      ...
      <wannier>
        <bandSelection minSpinUp="21" maxSpinUp="56"/>
        <jobList>  projgen prepwan90 stopopt  </jobList>
      </wannier>
   </output>

Additionally, we added wannier=T to the output-line. And we added a jobList with key projgen, which will generate a list of initial projections. Additionally we specified the key prepwan90 in jobList, which will generate the WF1.win and bkpts files. The key stopopt stops fleur_MPI after executing these two tasks. Since the highest band we need is maxSpinUp=56, we need to specify in inp.xml that we wish to calculate eigenfunctions up to the index 56. Therefore, we edit the line cutoffs in inp.xml by setting numbands to 56:

  <calculationSetup>
      <cutoffs Kmax="3.50000000" Gmax="10.50000000" GmaxXC="10.50000000" numbands="56"/>

Generate the projgen_inp file with the definition of initial projections:

Ga -3 0 0
As -3 0 0

Generate a kpts.xml file containing a uniform mesh of 8x8x8 k-points. Run fleur_MPI. Since the tasks projgen and prepwan90 are not parallelized, it is sufficient to run fleur_MPI on a single node. Check that the code has generated the proj file with the definitions of the projections for 16 MLWFs. Modify the jobList in inp.xml as follows:

   <output dos="F" band="F" vacdos="F" slice="F" mcd="F" wannier="T">
      ...
      <wannier>
        <bandSelection minSpinUp="21" maxSpinUp="56"/>
        <jobList>  matrixamn matrixmmn    </jobList>
      </wannier>
   </output>

Run fleur_MPI. The tasks matrixamn and matrixmmn are parallelized and you may run fleur_MPi in parallel. Here, we disentangle 16 MLWFs from 36 bands. Therefore, we specify the upper bound of the frozen window in WF1.win. Thus, edit the WF1.win file and set dis_froz_max=9.5 then run wannier90.x and converge the MLWFs. You may run wannier90.x in parallel. The converged MLWFs should look similar to this:

 Final State
  WF centre and spread    1  ( -2.251565, -2.251565, -2.251565 )     7.69807124
  WF centre and spread    2  ( -2.251565, -0.431935, -0.431935 )     7.69807124
  WF centre and spread    3  ( -0.431935, -2.251565, -0.431935 )     7.69807124
  WF centre and spread    4  ( -0.431935, -0.431935, -2.251565 )     7.69807124
  WF centre and spread    5  (  0.919214,  0.919214,  0.919214 )     7.13296569
  WF centre and spread    6  (  0.919214,  1.764286,  1.764286 )     7.13296569
  WF centre and spread    7  (  1.764286,  0.919214,  1.764286 )     7.13296569
  WF centre and spread    8  (  1.764286,  1.764286,  0.919214 )     7.13296569
  WF centre and spread    9  ( -2.251565, -2.251565, -2.251565 )     7.69807124
  WF centre and spread   10  ( -2.251565, -0.431935, -0.431935 )     7.69807124
  WF centre and spread   11  ( -0.431935, -2.251565, -0.431935 )     7.69807124
  WF centre and spread   12  ( -0.431935, -0.431935, -2.251565 )     7.69807124
  WF centre and spread   13  (  0.919214,  0.919214,  0.919214 )     7.13296569
  WF centre and spread   14  (  0.919214,  1.764286,  1.764286 )     7.13296569
  WF centre and spread   15  (  1.764286,  0.919214,  1.764286 )     7.13296569
  WF centre and spread   16  (  1.764286,  1.764286,  0.919214 )     7.13296569
  Sum of centres and spreads ( -0.000000, -0.000000, -0.000000 )   118.64829544

Now, you may use postw90.x in order to compute the shift current in GaAs according to the Wannier90 tutorial.