|
|
|
|
@ -23,23 +23,7 @@ from itertools import permutations, product
|
|
|
|
|
import numpy as np
|
|
|
|
|
from scipy.special import roots_legendre
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
# define some useful functions
|
|
|
|
|
def hsk(dh, k=(0, 0, 0)):
|
|
|
|
|
"""
|
|
|
|
|
One way to speed up Hk and Sk generation
|
|
|
|
|
"""
|
|
|
|
|
k = np.asarray(k, np.float64) # this two conversion lines
|
|
|
|
|
k.shape = (-1,) # are from the sisl source
|
|
|
|
|
|
|
|
|
|
# this generates the list of phases
|
|
|
|
|
phases = np.exp(-1j * np.dot(np.dot(np.dot(dh.rcell, k), dh.cell), dh.sc.sc_off.T))
|
|
|
|
|
|
|
|
|
|
HKU = np.einsum("abc,c->ab", dh.hup, phases)
|
|
|
|
|
HKD = np.einsum("abc,c->ab", dh.hdo, phases)
|
|
|
|
|
SK = np.einsum("abc,c->ab", dh.sov, phases)
|
|
|
|
|
|
|
|
|
|
return HKU, HKD, SK
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def make_contour(emin=-20, emax=0.0, enum=42, p=150):
|
|
|
|
|
@ -96,72 +80,120 @@ def make_kset(dirs="xyz", NUMK=20):
|
|
|
|
|
return kset
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def make_atran(nauc, dirs="xyz", dist=1):
|
|
|
|
|
"""
|
|
|
|
|
Simple pair generator. Depending on the value of the dirs
|
|
|
|
|
argument sampling in 1,2 or 3 dimensions is generated.
|
|
|
|
|
If dirs argument does not contain either of x,y or z
|
|
|
|
|
a single pair is returend.
|
|
|
|
|
"""
|
|
|
|
|
if not (sum([d in dirs for d in "xyz"])):
|
|
|
|
|
return (0, 0, [1, 0, 0])
|
|
|
|
|
# Pauli matrices
|
|
|
|
|
tau_x = np.array([[0, 1], [1, 0]])
|
|
|
|
|
tau_y = np.array([[0, -1j], [1j, 0]])
|
|
|
|
|
tau_z = np.array([[1, 0], [0, -1]])
|
|
|
|
|
tau_0 = np.array([[1, 0], [0, 1]])
|
|
|
|
|
|
|
|
|
|
dran = len(dirs) * [np.arange(-dist, dist + 1)]
|
|
|
|
|
mg = np.meshgrid(*dran)
|
|
|
|
|
dirsdict = dict()
|
|
|
|
|
|
|
|
|
|
for d in enumerate(dirs):
|
|
|
|
|
dirsdict[d[1]] = mg[d[0]].flatten()
|
|
|
|
|
for d in "xyz":
|
|
|
|
|
if not (d in dirs):
|
|
|
|
|
dirsdict[d] = 0 * dirsdict[dirs[0]]
|
|
|
|
|
def comm(a, b):
|
|
|
|
|
"Shorthand for commutator"
|
|
|
|
|
return a @ b - b @ a
|
|
|
|
|
|
|
|
|
|
ucran = np.array([dirsdict[d] for d in "xyz"]).T
|
|
|
|
|
atran = []
|
|
|
|
|
for i, j in list(product(range(nauc), repeat=2)):
|
|
|
|
|
for u in ucran:
|
|
|
|
|
if (abs(i - j) + sum(abs(u))) > 0:
|
|
|
|
|
atran.append((i, j, list(u)))
|
|
|
|
|
|
|
|
|
|
return atran
|
|
|
|
|
def tau_u(u):
|
|
|
|
|
"""
|
|
|
|
|
Pauli matrix in direction u.
|
|
|
|
|
"""
|
|
|
|
|
u = u / np.linalg.norm(u) # u is force to be of unit length
|
|
|
|
|
return u[0] * tau_x + u[1] * tau_y + u[2] * tau_z
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def add(x, y):
|
|
|
|
|
"""The sum of two numbers for testing.
|
|
|
|
|
#
|
|
|
|
|
def crossM(u):
|
|
|
|
|
"""
|
|
|
|
|
Definition for the cross-product matrix.
|
|
|
|
|
Acting as a crossproduct with vector u.
|
|
|
|
|
"""
|
|
|
|
|
return np.array([[0, -u[2], u[1]], [u[2], 0, -u[0]], [-u[1], u[0], 0]])
|
|
|
|
|
|
|
|
|
|
This function adds to numbers together. I only created this for testing documentation and examples.
|
|
|
|
|
|
|
|
|
|
Parameters
|
|
|
|
|
----------
|
|
|
|
|
x : float
|
|
|
|
|
First number
|
|
|
|
|
y : float
|
|
|
|
|
Second number added to `x`
|
|
|
|
|
def RotM(theta, u, eps=1e-10):
|
|
|
|
|
"""
|
|
|
|
|
Definition of rotation matrix with angle theta around direction u.
|
|
|
|
|
"""
|
|
|
|
|
u = u / np.linalg.norm(u)
|
|
|
|
|
|
|
|
|
|
Returns
|
|
|
|
|
-------
|
|
|
|
|
sum : int
|
|
|
|
|
The sum of the inputs
|
|
|
|
|
M = (
|
|
|
|
|
np.cos(theta) * np.eye(3)
|
|
|
|
|
+ np.sin(theta) * crossM(u)
|
|
|
|
|
+ (1 - np.cos(theta)) * np.outer(u, u)
|
|
|
|
|
)
|
|
|
|
|
M[abs(M) < eps] = 0.0 # kill off small numbers
|
|
|
|
|
return M
|
|
|
|
|
|
|
|
|
|
See Also
|
|
|
|
|
--------
|
|
|
|
|
numpy.add : Adds more than two numbers.
|
|
|
|
|
|
|
|
|
|
Notes
|
|
|
|
|
-----
|
|
|
|
|
We can create some latex notes here [1]_ :
|
|
|
|
|
def RotMa2b(a, b, eps=1e-10):
|
|
|
|
|
"""
|
|
|
|
|
Definition of rotation matrix rotating unit vector a to unit vector b.
|
|
|
|
|
Function returns array R such that R@a = b holds.
|
|
|
|
|
"""
|
|
|
|
|
v = np.cross(a, b)
|
|
|
|
|
c = a @ b
|
|
|
|
|
M = np.eye(3) + crossM(v) + crossM(v) @ crossM(v) / (1 + c)
|
|
|
|
|
M[abs(M) < eps] = 0.0 # kill off small numbers
|
|
|
|
|
return M
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
.. math:: a + b = c
|
|
|
|
|
def spin_tracer(M):
|
|
|
|
|
"""
|
|
|
|
|
Spin tracer utility.
|
|
|
|
|
This akes an operator with the orbital-spin sequence:
|
|
|
|
|
orbital 1 up,
|
|
|
|
|
orbital 1 down,
|
|
|
|
|
orbital 2 up,
|
|
|
|
|
orbital 2 down,
|
|
|
|
|
that is in the SPIN-BOX representation,
|
|
|
|
|
and extracts orbital dependent Pauli traces.
|
|
|
|
|
"""
|
|
|
|
|
|
|
|
|
|
References
|
|
|
|
|
----------
|
|
|
|
|
.. [1] https://numpydoc.readthedocs.io/en/latest/format.html
|
|
|
|
|
M11 = M[0::2, 0::2]
|
|
|
|
|
M12 = M[0::2, 1::2]
|
|
|
|
|
M21 = M[1::2, 0::2]
|
|
|
|
|
M22 = M[1::2, 1::2]
|
|
|
|
|
M_o = dict()
|
|
|
|
|
M_o["x"] = M12 + M21
|
|
|
|
|
M_o["y"] = 1j * (M12 - M21)
|
|
|
|
|
M_o["z"] = M11 - M22
|
|
|
|
|
M_o["c"] = M11 + M22
|
|
|
|
|
|
|
|
|
|
return M_o
|
|
|
|
|
|
|
|
|
|
Examples
|
|
|
|
|
--------
|
|
|
|
|
>>> add(1, 2)
|
|
|
|
|
3
|
|
|
|
|
|
|
|
|
|
def parse_magnetic_entity(dh, atom=None, l=None, **kwargs):
|
|
|
|
|
"""
|
|
|
|
|
Function to define orbital indeces of a given magnetic entity.
|
|
|
|
|
dh: a sisl Hamiltonian object
|
|
|
|
|
atom: an integer or list of integers, defining atom (or atoms) in the unicell forming the magnetic entity
|
|
|
|
|
l: integer, defining the angular momentum channel
|
|
|
|
|
"""
|
|
|
|
|
# case where we deal with more than one atom defining the magnetic entity
|
|
|
|
|
if type(atom) == list:
|
|
|
|
|
dat = []
|
|
|
|
|
for a in atom:
|
|
|
|
|
a_orb_idx = dh.geometry.a2o(a, all=True)
|
|
|
|
|
if (
|
|
|
|
|
type(l) == int
|
|
|
|
|
): # if specified we restrict to given l angular momentum channel inside each atom
|
|
|
|
|
a_orb_idx = a_orb_idx[[o.l == l for o in dh.geometry.atoms[a].orbitals]]
|
|
|
|
|
dat.append(a_orb_idx)
|
|
|
|
|
orbital_indeces = np.hstack(dat)
|
|
|
|
|
# case where we deal with a singel atom magnetic entity
|
|
|
|
|
elif type(atom) == int:
|
|
|
|
|
orbital_indeces = dh.geometry.a2o(atom, all=True)
|
|
|
|
|
if (
|
|
|
|
|
type(l) == int
|
|
|
|
|
): # if specified we restrict to given l angular momentum channel
|
|
|
|
|
orbital_indeces = orbital_indeces[
|
|
|
|
|
[o.l == l for o in dh.geometry.atoms[atom].orbitals]
|
|
|
|
|
]
|
|
|
|
|
|
|
|
|
|
return orbital_indeces # numpy array containing integers labeling orbitals associated to a magnetic entity.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def blow_up_orbindx(orb_indices):
|
|
|
|
|
"""
|
|
|
|
|
Function to blow up orbital indeces to make SPIN BOX indices.
|
|
|
|
|
"""
|
|
|
|
|
return x + y
|
|
|
|
|
return np.array([[2 * o, 2 * o + 1] for o in orb_indices]).flatten()
|
|
|
|
|
|
