"""
Code to calculate spin-orbit spillage.
Modified from JARVIS-Tools
https://www.nature.com/articles/s41598-019-45028-y
https://www.nature.com/articles/s41524-020-0319-4
"""
import numpy as np
from pymatgen.io.vasp.outputs import Wavecar
# from jarvis.io.vasp.outputs import Wavecar
[docs]class SOCSpillage(object):
"""
Spin-orbit spillage criteria to predict whether a material is topologically non-trival.
The spillage criteria physically signifies number of band-inverted electrons.
A non-zero, high value (generally >0.5) suggests non-trivial behavior
"""
def __init__(self, wf_noso="", wf_so=""):
"""
Requires path to WAVECAR files with and without LSORBIT = .TRUE.
Args:
wf_noso : WAVECAR without spin-orbit coupling
wf_so : WAVECAR with spin-orbit coupling
"""
self.wf_noso = wf_noso
self.wf_so = wf_so
[docs] def isclose(self, n1, n2, rel_tol=1e-7):
"""
Checking if the numbers are close enoung
"""
if abs(n1 - n2) < rel_tol:
return True
else:
return False
[docs] def orth(self, A):
"""
Helper function to create orthonormal basis
"""
u, s, vh = np.linalg.svd(A, full_matrices=False)
M, N = A.shape
eps = np.finfo(float).eps
tol = max(M, N) * np.amax(s) * eps
num = np.sum(s > tol, dtype=int)
Q = u[:, :num]
return Q, num
[docs] def overlap_so_spinpol(self):
"""
Main function to calculate SOC spillage
"""
noso = Wavecar(self.wf_noso)
so = Wavecar(self.wf_so)
bcell = np.linalg.inv(noso.a).T
tmp = np.linalg.norm(np.dot(np.diff(noso.kpoints, axis=0), bcell), axis=1)
noso_k = np.concatenate(([0], np.cumsum(tmp)))
noso_bands = np.array(noso.band_energy)[:, :, :, 0]
noso_kvecs = np.array(noso.kpoints)
noso_occs = np.array(noso.band_energy)[:, :, :, 2]
noso_nkpts = len(noso_k)
bcell = np.linalg.inv(so.a).T
tmp = np.linalg.norm(np.dot(np.diff(so.kpoints, axis=0), bcell), axis=1)
so_k = np.concatenate(([0], np.cumsum(tmp)))
so_bands = np.array([np.array(so.band_energy)[:, :, 0]])
so_kvecs = np.array(so.kpoints)
# so_occs = np.array([np.array(so.band_energy)[:, :, 2]])
so_nkpts = len(so_k)
nelec_list = []
for nk1 in range(1, noso_nkpts + 1): # no spin orbit kpoints loop
knoso = noso_kvecs[nk1 - 1, :]
for nk2 in range(1, so_nkpts + 1): # spin orbit
kso = so_kvecs[nk2 - 1, :]
if (
self.isclose(kso[0], knoso[0])
and self.isclose(kso[1], knoso[1])
and self.isclose(kso[2], knoso[2])
): # do kpoints match?
for c, e in enumerate(noso_occs[0, nk1 - 1, :]):
if e < 0.5:
cup = c
break
for c, e in enumerate(noso_occs[1, nk1 - 1, :]):
if e < 0.5:
cdn = c
break
nelec_list.append([cup, cdn, cup + cdn])
n_arr = np.array(nelec_list)
n_up = int(round(np.mean(n_arr[:, 0])))
n_dn = int(round(np.mean(n_arr[:, 1])))
n_tot = int(round(np.mean(n_arr[:, 2])))
nelec = int(n_tot)
# noso_homo_up = np.max(noso_bands[0, :, n_up - 1])
# noso_lumo_up = np.min(noso_bands[0, :, n_up])
# noso_homo_dn = np.max(noso_bands[1, :, n_dn - 1])
# noso_lumo_dn = np.min(noso_bands[1, :, n_dn])
so_homo = np.max(so_bands[0, :, nelec - 1])
so_lumo = np.min(so_bands[0, :, nelec])
# noso_direct_up = np.min(noso_bands[0, :, n_up] - noso_bands[0, :, n_up - 1])
# noso_direct_dn = np.min(noso_bands[1, :, n_dn] - noso_bands[1, :, n_dn - 1])
so_direct = np.min(so_bands[0, :, nelec] - so_bands[0, :, nelec - 1])
noso_direct = 1000000.0
noso_homo = -10000000.0
noso_lumo = 100000000.0
for i in range(noso_bands.shape[1]):
homo_k = max(noso_bands[0, i, n_up - 1], noso_bands[1, i, n_dn - 1])
lumo_k = min(noso_bands[0, i, n_up], noso_bands[1, i, n_dn])
noso_direct = min(noso_direct, lumo_k - homo_k)
noso_homo = max(noso_homo, homo_k)
noso_lumo = min(noso_lumo, lumo_k)
gamma_k = []
kpoints = []
np.set_printoptions(precision=4)
x = []
y = []
nelec_tot = 0.0
for nk1 in range(1, noso_nkpts + 1): # no spin orbit kpoints loop
knoso = noso_kvecs[nk1 - 1, :]
for nk2 in range(1, so_nkpts + 1): # spin orbit
kso = so_kvecs[nk2 - 1, :]
if (
self.isclose(kso[0], knoso[0])
and self.isclose(kso[1], knoso[1])
and self.isclose(kso[2], knoso[2])
): # do kpoints match?
# changes section 2
nelec_up = n_arr[nk1 - 1, 0]
nelec_dn = n_arr[nk1 - 1, 1]
nelec_tot = n_arr[nk1 - 1, 2]
kpoints.append(kso)
Mmn = 0.0
vnoso = np.array(
noso.coeffs[0][nk1 - 1][0]
) # noso.readBandCoeff(ispin=1, ikpt=nk1, iband=1, norm=False)
n_noso1 = vnoso.shape[0]
vnoso = np.array(
noso.coeffs[1][nk1 - 1][0]
) # noso.readBandCoeff(ispin=2, ikpt=nk1, iband=1, norm=False)
# n_noso2 = vnoso.shape[0]
vso = so.coeffs[nk1 - 1][
0
].flatten() # so.readBandCoeff(ispin=1, ikpt=nk2, iband=1, norm=False)
n_so = vso.shape[0]
vs = min(n_noso1 * 2, n_so)
Vnoso = np.zeros((vs, nelec_tot), dtype=complex)
Vso = np.zeros((vs, nelec_tot), dtype=complex)
if np.array(noso.coeffs[1][nk1 - 1]).shape[1] == vs // 2:
# if nk1==10 and nk2==10:
# print (np.array(noso.coeffs[1][nk1-1]).shape[1], )
# prepare matricies
for n1 in range(1, nelec_up + 1):
Vnoso[0:vs // 2, n1 - 1] = np.array(
noso.coeffs[0][nk1 - 1][n1 - 1]
)[0:vs // 2]
for n1 in range(1, nelec_dn + 1):
Vnoso[vs // 2:vs, n1 - 1 + nelec_up] = np.array(
noso.coeffs[1][nk1 - 1][n1 - 1]
)[0:vs // 2]
for n1 in range(1, nelec_tot + 1):
t = so.coeffs[nk2 - 1][n1 - 1].flatten()
Vso[0:vs // 2, n1 - 1] = t[0:vs // 2]
Vso[vs // 2:vs, n1 - 1] = t[
n_so // 2:n_so // 2 + vs // 2
]
Qnoso, num_noso = self.orth(Vnoso) # make orthonormal basis?
Qso, num_so = self.orth(Vso)
gamma_k.append(nelec_tot)
a = []
for n1 in range(0, nelec_tot): # noso occupied bands
v1 = Qnoso[:, n1]
aa = 0.0
for n2 in range(0, nelec_tot): # so occupied bands
v2 = Qso[:, n2]
t = np.dot(np.conj(v1), v2)
Mmn += t * t.conj()
aa += (t * t.conj()).real
a.append(aa)
gamma_k[-1] -= Mmn # eq 4 in prb 90 125133
x.append(kso)
y.append(np.real(gamma_k[-1]))
if gamma_k[-1] > 0.5:
print(
"nk1 nk2 kpoint gamma_k ",
nk1,
nk2,
kso,
knoso,
np.real(gamma_k[-1]),
"!!!!!!!!!!",
)
gmax = max(np.real(gamma_k))
nkmax = np.argmax(np.real(gamma_k))
kmax = kpoints[nkmax]
print("------------------------------------")
print()
print(" INDIRECT DIRECT HOMO/LUMO (eV)")
print(
"no spin-orbit gaps",
"{:+.3f}".format(float(noso_lumo - noso_homo)),
"{:+.3f}".format(noso_direct),
" ",
[noso_homo, noso_lumo],
)
print(
"spin-orbit gaps ",
"{:+.3f}".format(float(so_lumo - so_homo)),
"{:+.3f}".format(so_direct),
" ",
[so_homo, so_lumo],
)
print("gamma max", np.real(gmax), " at k = ", kmax)
return gmax