# Spin Hall conductivity

author: Jae-Mo Lihm (jaemo.lihm@gmail.com) and Minsu Ghim (minsu.ghim.physics@gmail.com)

In this tutorial, we calculate the spin Berry curvature and the spin Hall conductivity of bcc Platinum. We compare the two methods for calculating the spin velocity matrix, which we call the “Qiao” method  and the “Ryoo” method .

Both methods use the Kubo formula to calculate spin Hall conductivity under time-reversal symmetry:

\begin{equation} \sigma^{{\rm SHC}, \gamma}_{\alpha\beta} = \frac{-e\hbar}{N_k V_c}\sum_{\bf k}\sum_{n,m}\left(f_{n{\bf k}}-f_{m{\bf k}}\right)\frac{\textrm{Im}\left[\langle\psi_{n{\bf k}}\vert \frac{1}{2}\{ s^{\gamma}, v_\alpha \} \vert\psi_{m{\bf k}}\rangle\langle\psi_{m{\bf k}}\vert v_\beta\vert\psi_{n{\bf k}}\rangle\right]}{(\varepsilon_{n{\bf k}}-\varepsilon_{m{\bf k}})^2-(\hbar\omega+i\eta)^2}\,, \label{eq:shc}\tag{1} \end{equation}

where $$\alpha$$, $$\beta$$, $$\gamma$$ are respectively the direction of spin current, applied electric field, and spin polarisation.

The “Ryoo” method requires .chk, .eig, .mmn, .spn, .sHu, and .sIu files to calculate the spin velocity matrix in (\ref{eq:shc}), $$\langle\psi_{n{\bf k}}\vert \frac{1}{2}\{ s^{\gamma}, v_\alpha \} \vert\psi_{m{\bf k}}\rangle$$,from pw2wannier90.x, while the “Qiao” method does not use the last two files, and instead applies an approximation $$\mathbf{1}=\sum_{l\in \it{ab\,initio}} \vert u_{l{\bf q}}\rangle\langle u_{l{\bf q}}\vert$$. The sHu and sIu files are calculated by setting write_sHu = .true. and write_sIu = .true. to the pw2wannier90.x input file: see data_Pt/pw2wan.in.

:

# Preliminary (Do only once)

# Set environment variables - not mandatory but recommended
import os

import wannierberri as wberri
import numpy as np
import scipy
import matplotlib.pyplot as plt

#  This block is needed if you run this cell for a second time
#  because one cannot initiate two parallel environments at a time
try:
parallel.shutdown()
except NameError:
pass

# parallel = wberri.Parallel(num_cpus=1, progress_step_percent=10)
parallel = wberri.Serial(progress_step_percent=10)

The autoreload extension is already loaded. To reload it, use:
No need to shutdown Serial()


## Model, band structure

We load the system from a Wannier90 output. Note the arguments SHCryoo=True and SHCqiao=True which are required to compute spin Hall conductivity using the Ryoo and Qiao methods, respectively.

We set symmetry using the set_symmetry_from_structure method, which calls spglib to automatically determine the symmetry of the system.

:

system = wberri.System_w90("data_Pt/Pt", berry=True, SHCryoo=True, SHCqiao=True)
system.set_structure([[0., 0., 0.]], ["Pt"])
system.set_symmetry_from_structure()

efermi = 18.1605

using fortio to read
Reading restart information from file data_Pt/Pt.chk :
Time to read .chk : 0.15277957916259766
Time for MMN.__init__() : 0.8216688632965088 , read : 0.8122344017028809 , headstring 0.00943446159362793
----------
SPN
---------

reading data_Pt/Pt.spn : Created on 13May2022 at 15:23:23
----------
SPN OK
---------

----------
sIu
---------
reading data_Pt/Pt.sIu : <Created on 13May2022 at 15:23:23>
----------
sIu OK
---------

----------
sHu
---------
reading data_Pt/Pt.sHu : <Created on 13May2022 at 15:23:23>
----------
sHu OK
---------

time for FFT_q_to_R : 1.6445608139038086 s
using ws_distance
irvec_new_all shape (93,)
using ws_dist for Ham_R
using ws_dist for AA_R
using ws_dist for SS_R
using ws_dist for SA_R
using ws_dist for SHA_R
using ws_dist for SR_R
using ws_dist for SH_R
using ws_dist for SHR_R
Number of wannier functions: 18
Number of R points: 93
Recommended size of FFT grid [4 4 4]
Real-space lattice:
[[-1.95599772  0.          1.95599772]
[ 0.          1.95599772  1.95599772]
[-1.95599772  1.95599772  0.        ]]

:

path = wberri.Path(
system,
k_nodes=[
[0.25, 0.75, 0.50], # W
[0.50, 0.50, 0.50], # L
[0.00, 0.00, 0.00], # Gamma
[0.50, 0.00, 0.50], # X
[0.50, 0.25, 0.75], # W
[0.00, 0.00, 0.00], # Gamma
],
labels=["W", "L", "$\Gamma$", "X", "W", "$\Gamma$"],
length=300,
)

from wannierberri import calculators as calc
calculators = {}
calculators["tabulate"] = calc.TabulatorAll(
{"Energy": calc.tabulate.Energy()},
ibands=np.arange(system.num_wann),
mode="path",
)

path_result = wberri.run(
system,
grid=path,
calculators=calculators,
parallel=parallel,
print_Kpoints = False,
)

calculator not described

Calculation along a path - checking calculators for compatibility
tabulate <wannierberri.calculators.TabulatorAll object at 0x7fbc89a4b130>
All calculators are compatible
Symmetrization switched off for Path
The set of k points is a Path() with 322 points and labels {0: 'W', 54: 'L', 120: '$\\Gamma$', 197: 'X', 235: 'W', 321: '$\\Gamma$'}
WARNING : symmetry is not used for a tabulation along path
generating K_list
Done
Done, sum of weights:322.0
symgroup : None
processing 322 K points : in serial.
# K-points calculated  Wall time (sec)  Est. remaining (sec)
32              0.1                   1.0
64              0.2                   0.8
96              0.3                   0.7
128              0.4                   0.6
160              0.5                   0.5
192              0.6                   0.4
224              0.7                   0.3
256              0.8                   0.2
288              0.9                   0.1
320              1.0                   0.0
time for processing    322 K-points in serial:     1.0007 ; per K-point          0.0031 ; proc-sec per K-point          0.0031
time1 =  0.052553415298461914
Totally processed 322 K-points

:

fig = path_result.results["tabulate"].plot_path_fat(path, close_fig=False, show_fig=False)

ax = fig.get_axes()
ax.axhline(efermi, c="r", ls="--")
plt.show(fig) ## Static spin Hall conductivity

We calculate the static (i.e. DC) spin Hall conductivity. We fix $$\omega$$ to 0 and scan the Fermi energy.

:

from wannierberri import calculators as calc

efermi_list = np.linspace(efermi - 1.0, efermi + 1.0, 101, True)

kwargs = dict(
Efermi=efermi_list,
omega=np.array([0.]),
smr_fixed_width = 0.1, # Smearing for frequency in eV
kBT = 0.026, # Smearing for Fermi level (Fermi-Dirac factor) in eV (not Kelvin)
)

calculators = dict(
SHC_ryoo = calc.dynamic.SHC(SHC_type="ryoo", **kwargs),
SHC_qiao = calc.dynamic.SHC(SHC_type="qiao", **kwargs),
)

a more laconic implementation of the energy factor

a more laconic implementation of the energy factor


:

nk = 30
grid = wberri.Grid(system, NK=nk)
result = wberri.run(
system,
grid=grid,
calculators=calculators,
parallel=parallel,
print_Kpoints = False,
)

determining grids from NK=30 (<class 'int'>), NKdiv=None (<class 'NoneType'>), NKFFT=None (<class 'NoneType'>)
Minimal symmetric FFT grid :  [4 4 4]
The grids were set to NKdiv=[6 6 6], NKFFT=[5 5 5], NKtot=[30 30 30]
Grid is regular
The set of k points is a Grid() with NKdiv=[6 6 6], NKFFT=[5 5 5], NKtot=[30 30 30]
generating K_list
Done in 0.005416393280029297 s
excluding symmetry-equivalent K-points from initial grid
Done in 0.1416611671447754 s
Done in 0.14187836647033691 s
K_list contains 16 Irreducible points(7.41%) out of initial 6x6x6=216 grid
Done, sum of weights:0.9999999999999997
symgroup : <wannierberri.symmetry.Group object at 0x7ff9bc556490>
processing 16 K points : in serial.
# K-points calculated  Wall time (sec)  Est. remaining (sec)
2             13.7                  95.7
4             26.2                  78.7
6             41.2                  68.7
8             57.4                  57.4
10             74.9                  45.0
12             97.6                  32.5
14            121.8                  17.4
16            138.2                   0.0
time for processing     16 K-points in serial:   139.6814 ; per K-point          8.7301 ; proc-sec per K-point          8.7301
time1 =  0.008065938949584961
Totally processed 16 K-points


The SHC data has 5 indices: 1. The Fermi level index, 2. The frequency index, 3. The spin current direction index, 4. The electric field direction index, and 5. The spin polarization index.

:

print("result.results[\"SHC_ryoo\"].data.shape = ", result.results["SHC_ryoo"].data.shape)

result.results["SHC_ryoo"].data.shape =  (101, 1, 3, 3, 3)

:

shc_ryoo = result.results["SHC_ryoo"].data[:, 0, 0, 1, 2]
shc_qiao = result.results["SHC_qiao"].data[:, 0, 0, 1, 2]

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes.plot(efermi_list, shc_ryoo.real, label="Ryoo")
axes.plot(efermi_list, shc_qiao.real, label="Qiao")
axes.plot(efermi_list, shc_ryoo.imag)
axes.plot(efermi_list, shc_qiao.imag)
for ax in axes:
ax.set_xlabel("Efermi (eV)")
ax.axhline(0, c="k")
axes.set_ylabel("Re(SHC)")
axes.set_ylabel("Im(SHC)")
axes.legend()
plt.show() ## Dynamic spin Hall conductivity

We calculate the dynamic (i.e. frequency-dependent, AC) spin Hall conductivity. We fix the Fermi energy to the value efermi and scan the frequency in the range omega.

The smr_fixed_width parameter controls the smearing of the frequency-dependent terms (delta functions and principal values).

:

from wannierberri import calculators as calc

omega = np.linspace(0, 4, 101, True)

kwargs = dict(
Efermi=np.array([efermi]),
omega=omega,
smr_fixed_width = 0.1, # Smearing for frequency in eV
kBT = 0.026, # Smearing for Fermi level (Fermi-Dirac factor)
)

calculators = dict(
SHC_ryoo = calc.dynamic.SHC(SHC_type="ryoo", **kwargs),
SHC_qiao = calc.dynamic.SHC(SHC_type="qiao", **kwargs),
)

a more laconic implementation of the energy factor

a more laconic implementation of the energy factor


:

nk = 30
grid = wberri.Grid(system, NK=nk)
result = wberri.run(
system,
grid=grid,
calculators=calculators,
parallel=parallel,
print_Kpoints = False,
)

determining grids from NK=30 (<class 'int'>), NKdiv=None (<class 'NoneType'>), NKFFT=None (<class 'NoneType'>)
Minimal symmetric FFT grid :  [4 4 4]
The grids were set to NKdiv=[6 6 6], NKFFT=[5 5 5], NKtot=[30 30 30]
Grid is regular
The set of k points is a Grid() with NKdiv=[6 6 6], NKFFT=[5 5 5], NKtot=[30 30 30]
generating K_list
Done in 0.003348112106323242 s
excluding symmetry-equivalent K-points from initial grid
Done in 0.1131124496459961 s
Done in 0.11329960823059082 s
K_list contains 16 Irreducible points(7.41%) out of initial 6x6x6=216 grid
Done, sum of weights:0.9999999999999997
symgroup : <wannierberri.symmetry.Group object at 0x7ff9bc556490>
processing 16 K points : in serial.
# K-points calculated  Wall time (sec)  Est. remaining (sec)
2              9.5                  66.4
4             19.2                  57.6
6             29.1                  48.5
8             39.3                  39.3
10             48.4                  29.1
12             57.9                  19.3
14             67.5                   9.6
16             76.2                   0.0
time for processing     16 K-points in serial:    77.2919 ; per K-point          4.8307 ; proc-sec per K-point          4.8307
time1 =  0.008670330047607422
Totally processed 16 K-points

:

shc_ryoo = result.results["SHC_ryoo"].data[0, :, 0, 1, 2]
shc_qiao = result.results["SHC_qiao"].data[0, :, 0, 1, 2]

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes.plot(omega, shc_ryoo.real, label="Ryoo")
axes.plot(omega, shc_qiao.real, label="Qiao")
axes.plot(omega, shc_ryoo.imag)
axes.plot(omega, shc_qiao.imag)
for ax in axes:
ax.set_xlabel("omega (eV)")
ax.axhline(0, c="k")
axes.set_ylabel("Re(SHC)")
axes.set_ylabel("Im(SHC)")
axes.legend()
plt.show() ## Spin Berry curvature

To understand the microscopic origin of the spin Hall conductivity, one may inspect the k-resolved spin Berry curvature. Eq. (\ref{eq:shc}) is recast into the sum of a Berry-curvature-like term, the spin Berry curvature.

The spin Berry curvature is

:nbsphinx-math:begin{equation}

Omega^{n, gamma}_{alphabeta}({bf k}) = -sum_{m neq n}frac{2textrm{Im}left[langlepsi_{n{bf k}}vert frac{1}{2}{ s^{gamma}, v_alpha } vertpsi_{m{bf k}}ranglelanglepsi_{m{bf k}}vert v_betavertpsi_{n{bf k}}rangleright]}{(varepsilon_{n{bf k}}-varepsilon_{m{bf k}})^2-(ieta)^2},, label{eq:sbc}tag{2}

end{equation}

and the k-resolved spin Berry curvature summed over the band index is

:nbsphinx-math:begin{equation}

Omega^{gamma}_{alphabeta}({bf k}) = sum_{n}f_{n{bf k}}Omega^{n, gamma}_{alphabeta}({bf k}),. label{eq:sbc_k_resolved}tag{3}

end{equation} :nbsphinx-math:begin{equation}

sigma^{{rm SHC}, gamma}_{alphabeta} = frac{-ehbar}{N_kOmega_c}sum_{bf k}Omega^{gamma}_{alphabeta}({bf k}) label{eq:shc_sbc}tag{4}

end{equation}

Therefore, where in the k-space contributes to the total SHC can be investigated using the k-resolved spin Berry curvature.

Here, we compute the spin Berry curvature again using the Ryoo method and the Qiao method. Note that we pass the spin curren type as a kwargs_formula to the calculator, e.g. kwargs_formula=dict(spin_current_type="ryoo").

:

from wannierberri import calculators as calc
calculators = {}
calculators["tabulate"] = calc.TabulatorAll(
{
"Energy": calc.tabulate.Energy(),
'spin_berry_ryoo': calc.tabulate.SpinBerry(kwargs_formula=dict(spin_current_type="ryoo"), degen_thresh=1e-2),
'spin_berry_qiao': calc.tabulate.SpinBerry(kwargs_formula=dict(spin_current_type="qiao"), degen_thresh=1e-2),
},
ibands=np.arange(system.num_wann),
mode="path",
)

calculator not described

calculator not described

calculator not described


:

path = wberri.Path(
system,
k_nodes=[
[0.25, 0.75, 0.50], # W
[0.50, 0.50, 0.50], # L
[0.00, 0.00, 0.00], # Gamma
[0.50, 0.00, 0.50], # X
[0.50, 0.25, 0.75], # W
[0.00, 0.00, 0.00], # Gamma
],
labels=["W", "L", "$\Gamma$", "X", "W", "$\Gamma$"],
length=600,
)

result_spin_berry = wberri.run(
system,
grid=path,
calculators=calculators,
parallel = parallel,
print_Kpoints = False,
)

Calculation along a path - checking calculators for compatibility
tabulate <wannierberri.calculators.TabulatorAll object at 0x7ff9bb6850a0>
All calculators are compatible
Symmetrization switched off for Path
The set of k points is a Path() with 643 points and labels {0: 'W', 108: 'L', 241: '$\\Gamma$', 394: 'X', 471: 'W', 642: '$\\Gamma$'}
WARNING : symmetry is not used for a tabulation along path
generating K_list
Done
Done, sum of weights:643.0
symgroup : None
processing 643 K points : in serial.
# K-points calculated  Wall time (sec)  Est. remaining (sec)
64              2.6                  23.6
128              5.3                  21.5
192              8.5                  19.9
256             11.2                  16.9
320             13.5                  13.6
384             16.0                  10.8
448             18.4                   8.0
512             20.7                   5.3
576             23.2                   2.7
640             25.5                   0.1
time for processing    643 K-points in serial:    25.6413 ; per K-point          0.0399 ; proc-sec per K-point          0.0399
time1 =  0.17302536964416504
Totally processed 643 K-points


Now we sum over bands to compute the k-resolved spin Berry curvature:

$\Omega^{\gamma}_{\alpha\beta}({\bf k}) = \sum_{n \in {\rm occ.}} \Omega^{n, \gamma}_{\alpha\beta}({\bf k})$
:

nk = path.K_list.shape
spin_berry_ryoo = np.zeros((nk, 3, 3, 3))
spin_berry_qiao = np.zeros((nk, 3, 3, 3))

for iband in range(system.num_wann):
# Get the data for iband-th band
e = result_spin_berry.results["tabulate"].get_data("Energy", iband)
spin_berry_ryoo_nk = result_spin_berry.results["tabulate"].get_data("spin_berry_ryoo", iband)
spin_berry_qiao_nk = result_spin_berry.results["tabulate"].get_data("spin_berry_qiao", iband)

# Select k-point indices where the iband-th band is occupied
inds_occupied = e < efermi

# Add the spin Berry curvature of those bands
spin_berry_ryoo[inds_occupied] += spin_berry_ryoo_nk[inds_occupied]
spin_berry_qiao[inds_occupied] += spin_berry_qiao_nk[inds_occupied]

:

def get_signed_log10(x):
return np.log10(abs(x)) * np.sign(x)

kline = path.getKline()
plt.plot(kline, get_signed_log10(spin_berry_ryoo[:, 0, 1, 2]), "k-", label='Ryoo')
plt.plot(kline, get_signed_log10(spin_berry_qiao[:, 0, 1, 2]), "r--", label='Qiao')

for i in path.labels.keys():
plt.axvline(kline[i], c="k", lw=1)
plt.xticks([kline[i] for i in path.labels.keys()], path.labels.values())
plt.xlim([min(kline), max(kline)])
plt.axhline(0, c="k", lw=1)
plt.legend()
plt.title("$\mathrm{log}_{10} \Omega_\mathbf{k}$")

plt.show() You can find that the spin Berry curvature calculated using the Qiao method shows more “wiggles” than the Ryoo method. This numerical difference has been first reported in T. Ng et al, PRB 104 014412 (2021):

It is worth noting that there is jittering along Γ-Z, which occurs in the same path in WTe2 using the same method  (Qiao et al). However, such jittering disappears and the spin Berry curvature along Γ-Z becomes a smooth function using the method in Ref.  (Ryoo et al).

## Generating .sHu and .sIu from .mmn and .spn: mmn2uHu

Even if you have not obtained .shu and .sIu from an ab initio code, you can make them from the overlap matrix and the spin matrix. Wannierberri provides the utility wannierberri.utils.mmn2uHu, which calculated the matrices .uHu, .uIu, .sHu, and/or .sIu from the .mmn, .spn, .eig matrices, and also reduces the number of bands in .amn, .mmn, .eig and .spn files, by means of the sum-over-states formula

:nbsphinx-math:begin{equation}

langle u_{m{bf q}}verthat{s}hat{H}_{bf q}vert u_{n{bf q}+mathbf{b}}rangle approx sum_l^{l_{rm max}} left(s_{lm}({bf q})right)^* E_{l{bf q}} M_{ln}^{mathbf{b}}({bf q}),.

label{eq:sHu}tag{5} end{equation}

:nbsphinx-math:begin{equation}

langle u_{m{bf q}}verthat{s}vert u_{n{bf q}+mathbf{b}}rangle approx sum_l^{l_{rm max}} left(s_{lm}({bf q})right)^* M_{ln}^{mathbf{b}}({bf q}),.

label{eq:sIu}tag{6} end{equation}

Here, $$l_{\rm max}$$ cannot exceed the number of bands included in the Wannier90 calculation (i.e. the num_bands parameter in Pt.win).

The mmn2uHu utility can be particularly useful when the calculation of sHu and sIu files are not implemented in the DFT code you are using.

:

from wannierberri.utils import mmn2uHu as mmn2uHu
os.chdir("data_Pt")
mmn2uHu.run_mmn2uHu(PREFIX="Pt", writeSHU=True, writeSIU=True, NBout=18, NBsum=24)
#In case of direct execution of mmn2uHu module,
#python3 -m wannierberri.utils.mmn2uHu Pt NBout=18,NBsum=24,targets=sHu,sIu
os.chdir("..")

# Rename sHu and sIu files
import shutil
shutil.move("reduced_NB=24/Pt_nbs=24.sHu", "reduced_NB=24/Pt.sHu")
shutil.move("reduced_NB=24/Pt_nbs=24.sIu", "reduced_NB=24/Pt.sIu")

# Copy chk and spn files
shutil.copyfile("Pt.chk", "reduced_NB=24/Pt.chk")
shutil.copyfile("Pt.spn", "reduced_NB=24/Pt.spn")

os.chdir("..")

----------
---------

k-point 1 of 64
k-point 2 of 64
k-point 3 of 64
k-point 4 of 64
k-point 5 of 64
k-point 6 of 64
k-point 7 of 64
k-point 8 of 64
k-point 9 of 64
k-point 10 of 64
k-point 11 of 64
k-point 12 of 64
k-point 13 of 64
k-point 14 of 64
k-point 15 of 64
k-point 16 of 64
k-point 17 of 64
k-point 18 of 64
k-point 19 of 64
k-point 20 of 64
k-point 21 of 64
k-point 22 of 64
k-point 23 of 64
k-point 24 of 64
k-point 25 of 64
k-point 26 of 64
k-point 27 of 64
k-point 28 of 64
k-point 29 of 64
k-point 30 of 64
k-point 31 of 64
k-point 32 of 64
k-point 33 of 64
k-point 34 of 64
k-point 35 of 64
k-point 36 of 64
k-point 37 of 64
k-point 38 of 64
k-point 39 of 64
k-point 40 of 64
k-point 41 of 64
k-point 42 of 64
k-point 43 of 64
k-point 44 of 64
k-point 45 of 64
k-point 46 of 64
k-point 47 of 64
k-point 48 of 64
k-point 49 of 64
k-point 50 of 64
k-point 51 of 64
k-point 52 of 64
k-point 53 of 64
k-point 54 of 64
k-point 55 of 64
k-point 56 of 64
k-point 57 of 64
k-point 58 of 64
k-point 59 of 64
k-point 60 of 64
k-point 61 of 64
k-point 62 of 64
k-point 63 of 64
k-point 64 of 64
----------
---------

k-point 0 of 64
k-point 1 of 64
k-point 2 of 64
k-point 3 of 64
k-point 4 of 64
k-point 5 of 64
k-point 6 of 64
k-point 7 of 64
k-point 8 of 64
k-point 9 of 64
k-point 10 of 64
k-point 11 of 64
k-point 12 of 64
k-point 13 of 64
k-point 14 of 64
k-point 15 of 64
k-point 16 of 64
k-point 17 of 64
k-point 18 of 64
k-point 19 of 64
k-point 20 of 64
k-point 21 of 64
k-point 22 of 64
k-point 23 of 64
k-point 24 of 64
k-point 25 of 64
k-point 26 of 64
k-point 27 of 64
k-point 28 of 64
k-point 29 of 64
k-point 30 of 64
k-point 31 of 64
k-point 32 of 64
k-point 33 of 64
k-point 34 of 64
k-point 35 of 64
k-point 36 of 64
k-point 37 of 64
k-point 38 of 64
k-point 39 of 64
k-point 40 of 64
k-point 41 of 64
k-point 42 of 64
k-point 43 of 64
k-point 44 of 64
k-point 45 of 64
k-point 46 of 64
k-point 47 of 64
k-point 48 of 64
k-point 49 of 64
k-point 50 of 64
k-point 51 of 64
k-point 52 of 64
k-point 53 of 64
k-point 54 of 64
k-point 55 of 64
k-point 56 of 64
k-point 57 of 64
k-point 58 of 64
k-point 59 of 64
k-point 60 of 64
k-point 61 of 64
k-point 62 of 64
k-point 63 of 64
----------
MMN OK
---------

----------
AMN
---------

AMN size= (27648, 2)
24 18 64
----------
AMN  - OK
---------

[('uHu', False)]
----------
uHu  NBsum=24
---------
uHu from mmn red to 24 sum 24 bnd 2022-05-16T23:29:17.629221
60
using scipy.io to write
k-point 1 of 64
k-point 2 of 64
k-point 3 of 64
k-point 4 of 64
k-point 5 of 64
k-point 6 of 64
k-point 7 of 64
k-point 8 of 64
k-point 9 of 64
k-point 10 of 64
k-point 11 of 64
k-point 12 of 64
k-point 13 of 64
k-point 14 of 64
k-point 15 of 64
k-point 16 of 64
k-point 17 of 64
k-point 18 of 64
k-point 19 of 64
k-point 20 of 64
k-point 21 of 64
k-point 22 of 64
k-point 23 of 64
k-point 24 of 64
k-point 25 of 64
k-point 26 of 64
k-point 27 of 64
k-point 28 of 64
k-point 29 of 64
k-point 30 of 64
k-point 31 of 64
k-point 32 of 64
k-point 33 of 64
k-point 34 of 64
k-point 35 of 64
k-point 36 of 64
k-point 37 of 64
k-point 38 of 64
k-point 39 of 64
k-point 40 of 64
k-point 41 of 64
k-point 42 of 64
k-point 43 of 64
k-point 44 of 64
k-point 45 of 64
k-point 46 of 64
k-point 47 of 64
k-point 48 of 64
k-point 49 of 64
k-point 50 of 64
k-point 51 of 64
k-point 52 of 64
k-point 53 of 64
k-point 54 of 64
k-point 55 of 64
k-point 56 of 64
k-point 57 of 64
k-point 58 of 64
k-point 59 of 64
k-point 60 of 64
k-point 61 of 64
k-point 62 of 64
k-point 63 of 64
k-point 64 of 64
----------
uHu OK
---------

----------
SPN
---------

Created on 13May2022 at 15:23:23
using scipy.io to write
----------
SPN OK
---------

[('sHu', False), ('sIu', False)]
----------
sHu  NBsum=24
---------
sHu from mmn red to 24 sum 24 bnd 2022-05-16T23:29:18.531152
60
using scipy.io to write
k-point 1 of 64
k-point 2 of 64
k-point 3 of 64
k-point 4 of 64
k-point 5 of 64
k-point 6 of 64
k-point 7 of 64
k-point 8 of 64
k-point 9 of 64
k-point 10 of 64
k-point 11 of 64
k-point 12 of 64
k-point 13 of 64
k-point 14 of 64
k-point 15 of 64
k-point 16 of 64
k-point 17 of 64
k-point 18 of 64
k-point 19 of 64
k-point 20 of 64
k-point 21 of 64
k-point 22 of 64
k-point 23 of 64
k-point 24 of 64
k-point 25 of 64
k-point 26 of 64
k-point 27 of 64
k-point 28 of 64
k-point 29 of 64
k-point 30 of 64
k-point 31 of 64
k-point 32 of 64
k-point 33 of 64
k-point 34 of 64
k-point 35 of 64
k-point 36 of 64
k-point 37 of 64
k-point 38 of 64
k-point 39 of 64
k-point 40 of 64
k-point 41 of 64
k-point 42 of 64
k-point 43 of 64
k-point 44 of 64
k-point 45 of 64
k-point 46 of 64
k-point 47 of 64
k-point 48 of 64
k-point 49 of 64
k-point 50 of 64
k-point 51 of 64
k-point 52 of 64
k-point 53 of 64
k-point 54 of 64
k-point 55 of 64
k-point 56 of 64
k-point 57 of 64
k-point 58 of 64
k-point 59 of 64
k-point 60 of 64
k-point 61 of 64
k-point 62 of 64
k-point 63 of 64
k-point 64 of 64
----------
sHu OK
---------

----------
sIu  NBsum=24
---------
sIu from mmn red to 24 sum 24 bnd 2022-05-16T23:29:18.929873
60
using scipy.io to write
k-point 1 of 64
k-point 2 of 64
k-point 3 of 64
k-point 4 of 64
k-point 5 of 64
k-point 6 of 64
k-point 7 of 64
k-point 8 of 64
k-point 9 of 64
k-point 10 of 64
k-point 11 of 64
k-point 12 of 64
k-point 13 of 64
k-point 14 of 64
k-point 15 of 64
k-point 16 of 64
k-point 17 of 64
k-point 18 of 64
k-point 19 of 64
k-point 20 of 64
k-point 21 of 64
k-point 22 of 64
k-point 23 of 64
k-point 24 of 64
k-point 25 of 64
k-point 26 of 64
k-point 27 of 64
k-point 28 of 64
k-point 29 of 64
k-point 30 of 64
k-point 31 of 64
k-point 32 of 64
k-point 33 of 64
k-point 34 of 64
k-point 35 of 64
k-point 36 of 64
k-point 37 of 64
k-point 38 of 64
k-point 39 of 64
k-point 40 of 64
k-point 41 of 64
k-point 42 of 64
k-point 43 of 64
k-point 44 of 64
k-point 45 of 64
k-point 46 of 64
k-point 47 of 64
k-point 48 of 64
k-point 49 of 64
k-point 50 of 64
k-point 51 of 64
k-point 52 of 64
k-point 53 of 64
k-point 54 of 64
k-point 55 of 64
k-point 56 of 64
k-point 57 of 64
k-point 58 of 64
k-point 59 of 64
k-point 60 of 64
k-point 61 of 64
k-point 62 of 64
k-point 63 of 64
k-point 64 of 64
----------
sIu OK
---------


---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
/tmp/ipykernel_1923/239562535.py in <module>
4 #In case of direct execution of mmn2uHu module,
5 #python3 -m wannierberri.utils.mmn2uHu Pt NBout=18,NBsum=24,targets=sHu,sIu
----> 6 os.chdir("-")

FileNotFoundError: [Errno 2] No such file or directory: '-'

:

system_mmn2uhu = wberri.System_w90("data_Pt/reduced_NB=24/Pt", berry=True, SHCryoo=True, SHCqiao=True)
system_mmn2uhu.set_structure([[0., 0., 0.]], ["Pt"])
system_mmn2uhu.set_symmetry_from_structure()

using fortio to read
Reading restart information from file data_Pt/reduced_NB=24/Pt.chk :
Time to read .chk : 0.1561284065246582
Time for MMN.__init__() : 0.909125566482544 , read : 0.8991467952728271 , headstring 0.009978771209716797
----------
SPN
---------

reading data_Pt/reduced_NB=24/Pt.spn : Created on 13May2022 at 15:23:23
----------
SPN OK
---------

----------
sIu
---------
reading data_Pt/reduced_NB=24/Pt.sIu : <sIu from mmn red to 24 sum 24 bnd 2022-05-16T23:29:18.929873>
----------
sIu OK
---------

----------
sHu
---------
reading data_Pt/reduced_NB=24/Pt.sHu : <sHu from mmn red to 24 sum 24 bnd 2022-05-16T23:29:18.531152>
----------
sHu OK
---------

time for FFT_q_to_R : 1.7439723014831543 s
using ws_distance
irvec_new_all shape (93,)
using ws_dist for Ham_R
using ws_dist for AA_R
using ws_dist for SS_R
using ws_dist for SA_R
using ws_dist for SHA_R
using ws_dist for SR_R
using ws_dist for SH_R
using ws_dist for SHR_R
Number of wannier functions: 18
Number of R points: 93
Recommended size of FFT grid [4 4 4]
Real-space lattice:
[[-1.95599772  0.          1.95599772]
[ 0.          1.95599772  1.95599772]
[-1.95599772  1.95599772  0.        ]]

:

from wannierberri import calculators as calc

efermi_list = np.linspace(efermi - 1.0, efermi + 1.0, 101, True)

kwargs = dict(
Efermi=efermi_list,
omega=np.array([0.]),
smr_fixed_width = 0.1, # Smearing for frequency in eV
kBT = 0.026, # Smearing for Fermi level (Fermi-Dirac factor) in eV (not Kelvin)
)

calculators = dict(
SHC_ryoo = calc.dynamic.SHC(SHC_type="ryoo", **kwargs),
)

nk = 30
result_pw2w90 = wberri.run(
system,
grid=wberri.Grid(system, NK=nk),
calculators=calculators,
parallel=parallel,
print_Kpoints = False,
)

result_mmn2uhu = wberri.run(
system_mmn2uhu,
grid=wberri.Grid(system, NK=nk),
calculators=calculators,
parallel=parallel,
print_Kpoints = False,
)

a more laconic implementation of the energy factor

determining grids from NK=30 (<class 'int'>), NKdiv=None (<class 'NoneType'>), NKFFT=None (<class 'NoneType'>)
Minimal symmetric FFT grid :  [4 4 4]
The grids were set to NKdiv=[6 6 6], NKFFT=[5 5 5], NKtot=[30 30 30]
Grid is regular
The set of k points is a Grid() with NKdiv=[6 6 6], NKFFT=[5 5 5], NKtot=[30 30 30]
generating K_list
Done in 0.0029904842376708984 s
excluding symmetry-equivalent K-points from initial grid
Done in 0.13321137428283691 s
Done in 0.13338375091552734 s
K_list contains 16 Irreducible points(7.41%) out of initial 6x6x6=216 grid
Done, sum of weights:0.9999999999999997
symgroup : <wannierberri.symmetry.Group object at 0x7f676be449a0>
processing 16 K points : in serial.
# K-points calculated  Wall time (sec)  Est. remaining (sec)
2              8.5                  59.2
4             19.2                  57.5
6             26.7                  44.5
8             34.9                  34.9
10             42.3                  25.4
12             49.7                  16.6
14             57.2                   8.2
16             64.8                   0.0
time for processing     16 K-points in serial:    65.4276 ; per K-point          4.0892 ; proc-sec per K-point          4.0892
time1 =  0.003787517547607422
Totally processed 16 K-points
determining grids from NK=30 (<class 'int'>), NKdiv=None (<class 'NoneType'>), NKFFT=None (<class 'NoneType'>)
Minimal symmetric FFT grid :  [4 4 4]
The grids were set to NKdiv=[6 6 6], NKFFT=[5 5 5], NKtot=[30 30 30]
Grid is regular
The set of k points is a Grid() with NKdiv=[6 6 6], NKFFT=[5 5 5], NKtot=[30 30 30]
generating K_list
Done in 0.0033483505249023438 s
excluding symmetry-equivalent K-points from initial grid
Done in 0.12982869148254395 s
Done in 0.12998056411743164 s
K_list contains 16 Irreducible points(7.41%) out of initial 6x6x6=216 grid
Done, sum of weights:0.9999999999999997
symgroup : <wannierberri.symmetry.Group object at 0x7f676c2c1a60>
processing 16 K points : in serial.
# K-points calculated  Wall time (sec)  Est. remaining (sec)
2              7.5                  52.8
4             15.1                  45.2
6             23.9                  39.8
8             31.9                  31.9
10             39.4                  23.7
12             47.0                  15.7
14             54.5                   7.8
16             61.9                   0.0
time for processing     16 K-points in serial:    62.4541 ; per K-point          3.9034 ; proc-sec per K-point          3.9034
time1 =  0.004681825637817383
Totally processed 16 K-points

:

shc_pw2w90 = result_pw2w90.results["SHC_ryoo"].data[:, 0, 0, 1, 2]
shc_mmn2uHu = result_mmn2uhu.results["SHC_ryoo"].data[:, 0, 0, 1, 2]

plt.plot(efermi_list, shc_pw2w90.real, label="sHu from pw2wannier90")
plt.plot(efermi_list, shc_mmn2uHu.real, label="sHu from mmn2uHu")
plt.xlabel("Efermi (eV)")
plt.axhline(0, c="k")
plt.ylabel("Re(SHC)")
plt.legend()
plt.show() ## Further questions

If you are interested, try to answer the following questions: - Try to converge the calculation using a different value of smr_fixed_width. In principle, to achieve an ideal convergence to the zero-smearing limit, one needs to first converge SHC increasing the grid size for a fixed smr_fixed_width, and then repeat the procedure with smaller smr_fixed_width until convergence. - What happens if one include more bands in the NSCF calculation? Does the two methods converge to the same result? (To answer this question, one needs to perform additional DFT calculations.)