Tutorial: Si with SOC (VASP)

In this tutorial we compute Si effective masses with spin-orbit coupling (SOC) and GGA-PBE XC. The input files are

INCAR

SYSTEM=Si
ISTART=0
ICHARG=2

LOPTICS = T    # compute dipole matrix elements and write them in WAVEDER file
LPEAD = F      # do not use finite difference in k space to get derivative of the cell-periodic part of the orbitals
NBANDS = 180   # more empty bands

ALGO = Normal
EDIFF=1.0E-09

LCHARG = F     # store charge density?
LWAVE = F      # store FWs?

ISMEAR=0; SIGMA=0.01

# SOC
LNONCOLLINEAR = T
NUPDOWN=0
LSORBIT = T
SAXIS =   0 0 1
GGA_COMPAT = F

Since this is a small calculation, we can do SCF with a large number of bands in one shot. A more efficient (for large cases) is a 2-step approach: (1) get the SCF charge density with a default number of bands, and then (2) add more bands. See steps 1 and 2 on VASP Wiki. Results for m* are quite forgiving to ENCUT, so there is no need to raise it beyond ENMAX in the POTCAR file.

POSCAR

Si experimental a0
   5.431
     0.5    0.5    0.0000000000000000
     0.0000000000000000    0.5    0.5
     0.5   0.0000000000000000    0.5
   Si
     2
Direct
  0.0000000000000000  0.0000000000000000  0.0000000000000000
  0.2500000000000000  0.2500000000000000  0.2500000000000000

KPOINTS

Automatic
 0
Gamma
 7 7 7
 0 0 0

Since the conduction band min is not at Gamma, the choice of k point mesh is important. With 7x7x7 you get a k point (#15) very close to the minimum. Results for m* are quite forgiving to the k-mesh because it is not a finite difference technique. There is no need to make a dense mesh as one would normally do for optical calculations.

POTCAR: PAW_PBE Si 05Jan2001 potential was used from potpaw_PBE.54/Si/POTCAR

Once INCAR, POSCAR, KPOINTS, and POTCAR files are in place, we can run VASP calculation (noncollinear mode because of SOC)

/path/to/VASP/bin/vasp_ncl

Once VASP calculation is dine, execute mstar

~/mstar/mstar WAVEDER 1e-5

The output is

 Detected input arguments = 2
 Input mommat file = WAVEDER
 Degeneracy tolerance dEtol = 1e-5 [Ha]
 Confirming text-to-number conversion dEtol =  1.00000E-05 [Ha]
 Assume VASP calculation
 The input file WAVEDER was found.
 The file EIVENVAL is also required. It is found.
  number of spins in WAVEDER file = 1
  number of k-points in WAVEDER file = 20
  smaller number of bands in WAVEDER file = 16
  number of bands per k-point in EIGENVAL file = 180
 Memory required to store matrix elements from WAVEDER file   0.0 GB
 Memory required to store dE_ij from EIGENVAL file   0.0 GB
 It will take additional   0.0 GB overhead to run the calculation
 Overall, you will need at least   0.0 GB (+ 20% incidental) of RAM to run the calculation
 WAVEDER and EIGENVAL files were read successfully
 Entering the main loop...
 KP: 1 bands: 180 progress:   5%
 KP: 2 bands: 180 progress:  10%
 KP: 3 bands: 180 progress:  15%
 KP: 4 bands: 180 progress:  20%
 KP: 5 bands: 180 progress:  25%
 KP: 6 bands: 180 progress:  30%
 KP: 7 bands: 180 progress:  35%
 KP: 8 bands: 180 progress:  40%
 KP: 9 bands: 180 progress:  45%
 KP: 10 bands: 180 progress:  50%
 KP: 11 bands: 180 progress:  55%
 KP: 12 bands: 180 progress:  60%
 KP: 13 bands: 180 progress:  65%
 KP: 14 bands: 180 progress:  70%
 KP: 15 bands: 180 progress:  75%
 KP: 16 bands: 180 progress:  80%
 KP: 17 bands: 180 progress:  85%
 KP: 18 bands: 180 progress:  90%
 KP: 19 bands: 180 progress:  95%
 KP: 20 bands: 180 progress: 100%
 Summary of the output:
 (1) Components of the (m0/m*_ij) tensor are stored in file minv_ij.dat
 (2) The conductivity inverse eff. masses m0/m_c = 1/3*Tr(m0/m*_ij)
     are stored in file minv_c.dat
 (3) Principal components of the (m0/m*_ij) tensor
     are stored in file minv_pr.dat
 (4) Density of states inverse effective mass
     m0/m*_d = m0/(m_1*m_2*m_3)**(1/3) are stored in file minv_d.dat
 See the file header for the description

Let's inspect the output file minv_ij.dat and locate the 1st k point (Gamma point). Bands 3-8 correspond to split off (SO), light hole (LH), and heavy hole (HH):

...
# band 1=xx;     2=yy;      3=zz;     4=yz;       5=xz;      6=xy
# KP: 1 NBCDER: 16 NEMAX: 180
...
 3 -4.339E+00 -4.339E+00 -4.339E+00 -3.356E-04 -3.340E-04 -3.359E-04
 4 -4.339E+00 -4.339E+00 -4.339E+00 -3.305E-04 -3.263E-04 -3.295E-04
 5 -5.185E+00 -5.185E+00 -5.185E+00 -2.640E+00 -2.640E+00 -2.640E+00
 6 -5.185E+00 -5.185E+00 -5.185E+00 -2.640E+00 -2.640E+00 -2.640E+00
 7 -3.737E+00 -3.737E+00 -3.737E+00  2.640E+00  2.640E+00  2.640E+00
 8 -3.737E+00 -3.737E+00 -3.737E+00  2.640E+00  2.640E+00  2.640E+00
...

Each eigenstate is double degenerate (due to SOC and Gamma symmetry). The effective masses in [100] direction are

m_SO(xx)/m0 = 1/(-4.339) = -0.23
m_LH(xx)/m0 = 1/(-5.185) = -0.19
m_HH(xx)/m0 = 1/(-3.737) = -0.27

Bottom of the conduction band (bands 9-10) is located at k point #15. Here is the relevant section of the output file minv_ij.dat:

# KP: 15 NBCDER: 16 NEMAX: 180
...
 9  5.181E+00  1.062E+00  5.181E+00 -1.453E-06 -1.321E-04 -1.445E-06
10  5.181E+00  1.062E+00  5.181E+00 -1.443E-06 -1.318E-04 -1.437E-06
...

We can determine longitudinal and transverse masses in the conduction band:

m_CB(xx)/m0 = 1/5.181 = 0.19
m_CB(yy)/m0 = 1/1.062 = 0.94

Results agree almost exactly with the WIEN2k all electron calculations.