Tutorial Si with SOC (WIEN2k) - rubel75/mstar GitHub Wiki
We calculate the effective masses in Si with spin-orbit coupling (SOC) with GGA-PBE XC. Here is a master script to run the tutorial:
#!/bin/bash
# SCRATCH setting will ensure that vector files are in the case folder
export SCRATCH=./
# case name = directory name
case=${PWD##*/}
# change editor settings
EDITOR2=$EDITOR # save settings
export EDITOR=cat # prevent opening files for editing
# initialization
cp /path/to/mstar/test_files/W2k/Si.struct ${case}.struct
init_lapw -b -vxc 13 -rkmax 7 -numk 500
# run SCF without SOC
run_lapw -ec 0.00001 -cc 0.0001
save_lapw -d noSOC
# init SOC
echo -e "0 0 1\n\n\n\nN\n" | init_so_lapw
# run with SOC
run_lapw -ec 0.00001 -cc 0.0001 -so
# increase Emax
sed -i 's/4 -9\.0 5\.0/4 -9.0 10.0/g' ${case}.in1 # Emax=10 Ry for LAPW1
sed -i 's/^-10 1\.9/-10 7.0/g' ${case}.inso # Emax=7 Ry for LAPWSO
# WFs and eigenvalues with increased band range
x lapw1
x lapwso
# fake spin-polarized calculation for optic
rm ${case}.vspup
ln -s ${case}.vsp ${case}.vspup
rm ${case}.vspdn
ln -s ${case}.vsp ${case}.vspdn
ln -s ${case}.vectorso ${case}.vectorsoup
# optics
cp $WIENROOT/SRC_templates/case.inop ${case}.inop # get input file template
sed -i 's/^OFF/ON/g' ${case}.inop # enable writing of momentum matrix elements in *.mommat2*
sed -i 's/^-5\.0 3\.0/-5.0 7.0/g' ${case}.inop # extend Emax = 7 Ry
x optic -so -up
# run mstar
/path/to/mstar/mstar ${case}.mommat2up 1e-5
# reset the editor settings back
export EDITOR=$EDITOR2
Results for m* are quite forgiving to RKmax parameter and k-mesh density since it is not a finite difference technique. There is no need to make a dense k mesh as one would normally do for optical calculations.
When mstar calculation is executed, you will see this output (20 k points, ~180 bands each):
Detected input arguments = 2
Input mommat file = mass.mommat2up
Degeneracy tolerance dEtol = 1e-5 [Ha]
Confirming text-to-number conversion dEtol = 1.00000E-05 [Ha]
The input file mass.mommat2up was found.
number of lines in mommat file = 327804
Entering the main loop...
KP: 1 bands: 180 progress: 4%
KP: 2 bands: 176 progress: 9%
KP: 3 bands: 180 progress: 14%
KP: 4 bands: 176 progress: 19%
KP: 5 bands: 176 progress: 24%
KP: 6 bands: 176 progress: 28%
KP: 7 bands: 178 progress: 33%
KP: 8 bands: 180 progress: 38%
KP: 9 bands: 180 progress: 43%
KP: 10 bands: 182 progress: 48%
KP: 11 bands: 176 progress: 53%
KP: 12 bands: 180 progress: 58%
KP: 13 bands: 186 progress: 63%
KP: 14 bands: 180 progress: 68%
KP: 15 bands: 184 progress: 74%
KP: 16 bands: 188 progress: 79%
KP: 17 bands: 180 progress: 84%
KP: 18 bands: 180 progress: 89%
KP: 19 bands: 186 progress: 94%
KP: 20 bands: 186 progress: 100%
Summary of the output:
(1) Components of the (m0/m*_ij) tensor are stored in file minv_ij-up.dat
(2) The conductivity inverse eff. masses m0/m_c = 1/3*Tr(m0/m*_ij)
are stored in file minv_c-up.dat
(3) Principal components of the (m0/m*_ij) tensor
are stored in file minv_pr-up.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-up.dat
Let's inspect the output file minv_ij-up.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: 180 NEMAX: 180
...
3 -4.335E+00 -4.335E+00 -4.335E+00 -6.862E-09 -2.179E-07 -4.309E-07
4 -4.335E+00 -4.335E+00 -4.335E+00 5.618E-07 4.223E-07 4.312E-07
5 -5.159E+00 -5.159E+00 -5.159E+00 -2.639E+00 -2.639E+00 -2.639E+00
6 -5.159E+00 -5.159E+00 -5.159E+00 -2.639E+00 -2.639E+00 -2.639E+00
7 -3.738E+00 -3.738E+00 -3.738E+00 2.639E+00 2.639E+00 2.639E+00
8 -3.738E+00 -3.738E+00 -3.738E+00 2.639E+00 2.639E+00 2.639E+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.335) = -0.23
m_LH(xx)/m0 = 1/(-5.159) = -0.19
m_HH(xx)/m0 = 1/(-3.738) = -0.27
and well agree with results obtained from the band curvature -0.23, -0.19, and -0.26.
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-up.dat:
# KP: 15 NBCDER: 184 NEMAX: 184
...
9 1.071E+00 5.186E+00 5.186E+00 -3.791E-07 -3.873E-11 -2.610E-11
10 1.071E+00 5.186E+00 5.186E+00 -4.346E-08 3.840E-11 9.193E-12
...
We can determine longitudinal and transverse masses in the conduction band:
m_CB(xx)/m0 = 1/1.071 = 0.93
m_CB(yy)/m0 = 1/5.186 = 0.19
The corresponding band curvature masses are 0.96 and 0.20. The error of 3% is acceptable.
If you want to get a feel for convergence with respect to Emax, you can edit ${case}.inop and set Emax to a smaller value (e.g., 1 Ry). Then re-run optic and mstar. You should see fewer bands in the output minv_ij-up.dat file. Also, the accuracy of m* will drop (see Fig. 1 in arXiv:2007.03816 [cond-mat.mtrl-sci]).