You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
383 lines
13 KiB
383 lines
13 KiB
# Copyright (c) [2024] [Daniel Pozsar]
|
|
#
|
|
# Permission is hereby granted, free of charge, to any person obtaining a copy
|
|
# of this software and associated documentation files (the "Software"), to deal
|
|
# in the Software without restriction, including without limitation the rights
|
|
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|
# copies of the Software, and to permit persons to whom the Software is
|
|
# furnished to do so, subject to the following conditions:
|
|
#
|
|
# The above copyright notice and this permission notice shall be included in all
|
|
# copies or substantial portions of the Software.
|
|
#
|
|
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
|
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
|
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
|
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
|
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
|
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
|
|
# SOFTWARE.
|
|
|
|
from itertools import permutations, product
|
|
from pprint import pprint
|
|
|
|
import numpy as np
|
|
from scipy.special import roots_legendre
|
|
|
|
# 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]])
|
|
|
|
|
|
# define some useful functions
|
|
def hsk(H, ss, sc_off, 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 * 2 * np.pi * k @ sc_off.T)
|
|
|
|
HK = np.einsum("abc,a->bc", H, phases)
|
|
SK = np.einsum("abc,a->bc", ss, phases)
|
|
|
|
return HK, SK
|
|
|
|
|
|
def make_contour(emin=-20, emax=0.0, enum=42, p=150):
|
|
"""
|
|
A more sophisticated contour generator
|
|
"""
|
|
|
|
x, wl = roots_legendre(enum)
|
|
R = (emax - emin) / 2
|
|
z0 = (emax + emin) / 2
|
|
y1 = -np.log(1 + np.pi * p)
|
|
y2 = 0
|
|
|
|
y = (y2 - y1) / 2 * x + (y2 + y1) / 2
|
|
phi = (np.exp(-y) - 1) / p
|
|
ze = z0 + R * np.exp(1j * phi)
|
|
we = -(y2 - y1) / 2 * np.exp(-y) / p * 1j * (ze - z0) * wl
|
|
|
|
class ccont:
|
|
# just an empty container class
|
|
pass
|
|
|
|
cont = ccont()
|
|
cont.R = R
|
|
cont.z0 = z0
|
|
cont.ze = ze
|
|
cont.we = we
|
|
cont.enum = enum
|
|
|
|
return cont
|
|
|
|
|
|
def make_kset(dirs="xyz", NUMK=20):
|
|
"""
|
|
Simple k-grid generator. Depending on the value of the dirs
|
|
argument k sampling in 1,2 or 3 dimensions is generated.
|
|
If dirs argument does not contain either of x,y or z
|
|
a kset of a single k-pont at the origin is returend.
|
|
"""
|
|
if not (sum([d in dirs for d in "xyz"])):
|
|
return np.array([[0, 0, 0]])
|
|
|
|
kran = len(dirs) * [np.linspace(0, 1, NUMK, endpoint=False)]
|
|
mg = np.meshgrid(*kran)
|
|
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]]
|
|
kset = np.array([dirsdict[d] for d in "xyz"]).T
|
|
|
|
return kset
|
|
|
|
|
|
def commutator(a, b):
|
|
"Shorthand for commutator"
|
|
return a @ b - b @ a
|
|
|
|
|
|
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 crossM(u):
|
|
"""
|
|
Definition for the cross-product matrix.
|
|
Acting as a cross product with vector u.
|
|
"""
|
|
return np.array([[0, -u[2], u[1]], [u[2], 0, -u[0]], [-u[1], u[0], 0]])
|
|
|
|
|
|
def RotM(theta, u, eps=1e-10):
|
|
"""
|
|
Definition of rotation matrix with angle theta around direction u.
|
|
"""
|
|
u = u / np.linalg.norm(u)
|
|
|
|
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
|
|
|
|
|
|
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
|
|
|
|
|
|
def spin_tracer(M):
|
|
"""
|
|
Spin tracer utility.
|
|
This takes 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.
|
|
"""
|
|
|
|
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
|
|
|
|
|
|
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 np.array([[2 * o, 2 * o + 1] for o in orb_indices]).flatten()
|
|
|
|
|
|
def calculate_exchange_tensor(pair):
|
|
energies = pair["energies"]
|
|
# Initialize output arrays
|
|
J = np.zeros((3, 3))
|
|
D = np.zeros(3)
|
|
|
|
# J matrix calculations
|
|
# J(1,1) = mean([DEij(2,2,2), DEij(2,2,3)])
|
|
J[0, 0] = np.mean([energies[1, 3], energies[2, 3]])
|
|
|
|
# J(1,2) = -mean([DEij(1,2,3), DEij(2,1,3)])
|
|
J[0, 1] = -np.mean([energies[2, 1], energies[2, 2]])
|
|
J[1, 0] = J[0, 1]
|
|
|
|
# J(1,3) = -mean([DEij(1,2,2), DEij(2,1,2)])
|
|
J[0, 2] = -np.mean([energies[1, 1], energies[1, 2]])
|
|
J[2, 0] = J[0, 2]
|
|
|
|
# J(2,2) = mean([DEij(2,2,1), DEij(1,1,3)])
|
|
J[1, 1] = np.mean([energies[0, 3], energies[2, 0]])
|
|
|
|
# J(2,3) = -mean([DEij(1,2,1), DEij(2,1,1)])
|
|
J[1, 2] = -np.mean([energies[0, 1], energies[0, 2]])
|
|
J[2, 1] = J[1, 2]
|
|
|
|
# J(3,3) = mean([DEij(1,1,1), DEij(1,1,2)])
|
|
J[2, 2] = np.mean([energies[0, 0], energies[1, 0]])
|
|
|
|
# D vector calculations
|
|
# D(1) = mean([DEij(1,2,1), -DEij(2,1,1)])
|
|
D[0] = np.mean([energies[0, 1], -energies[0, 2]])
|
|
|
|
# D(2) = mean([DEij(2,1,2), -DEij(1,2,2)])
|
|
D[1] = np.mean([energies[1, 2], -energies[1, 1]])
|
|
|
|
# D(3) = mean([DEij(1,2,3), -DEij(2,1,3)])
|
|
D[2] = np.mean([energies[2, 1], -energies[2, 2]])
|
|
|
|
J_iso = np.trace(J) / 3
|
|
J_S = (J - J_iso * np.eye(3)).flatten()
|
|
|
|
return J_iso, J_S, D, J
|
|
|
|
|
|
def print_parameters(simulation_parameters):
|
|
print(
|
|
"================================================================================================================================================================"
|
|
)
|
|
print("Input file: ")
|
|
print(simulation_parameters["path"])
|
|
print("Output file: ")
|
|
print(simulation_parameters["outpath"])
|
|
print(
|
|
"Number of nodes in the parallel cluster: ",
|
|
simulation_parameters["parallel_size"],
|
|
)
|
|
print(
|
|
"================================================================================================================================================================"
|
|
)
|
|
print("Cell [Ang]: ")
|
|
print(simulation_parameters["cell"])
|
|
print(
|
|
"================================================================================================================================================================"
|
|
)
|
|
print("DFT axis: ")
|
|
print(simulation_parameters["scf_xcf_orientation"])
|
|
print("Quantization axis and perpendicular rotation directions:")
|
|
for ref in simulation_parameters["ref_xcf_orientations"]:
|
|
print(ref["o"], " --» ", ref["vw"])
|
|
print(
|
|
"================================================================================================================================================================"
|
|
)
|
|
print("Parameters for the contour integral:")
|
|
print("Number of k points: ", simulation_parameters["kset"])
|
|
print("k point directions: ", simulation_parameters["kdirs"])
|
|
print("Ebot: ", simulation_parameters["ebot"])
|
|
print("Eset: ", simulation_parameters["eset"])
|
|
print("Esetp: ", simulation_parameters["esetp"])
|
|
print(
|
|
"================================================================================================================================================================"
|
|
)
|
|
|
|
|
|
def print_atoms_and_pairs(magnetic_entities, pairs):
|
|
print("Atomic information: ")
|
|
print(
|
|
"----------------------------------------------------------------------------------------------------------------------------------------------------------------"
|
|
)
|
|
print(
|
|
"[atom index]Element(orbitals) x [Ang] y [Ang] z [Ang] Sx Sy Sz Q Lx Ly Lz Jx Jy Jz"
|
|
)
|
|
print(
|
|
"----------------------------------------------------------------------------------------------------------------------------------------------------------------"
|
|
)
|
|
# iterate over magnetic entities
|
|
for mag_ent in magnetic_entities:
|
|
# iterate over atoms
|
|
for tag, xyz in zip(mag_ent["tags"], mag_ent["xyz"]):
|
|
# coordinates and tag
|
|
print(f"{tag} {xyz[0]} {xyz[1]} {xyz[2]}")
|
|
print("")
|
|
|
|
print(
|
|
"================================================================================================================================================================"
|
|
)
|
|
print("Exchange [meV]")
|
|
print(
|
|
"----------------------------------------------------------------------------------------------------------------------------------------------------------------"
|
|
)
|
|
print("Magnetic entity1 Magnetic entity2 [i j k] d [Ang]")
|
|
print(
|
|
"----------------------------------------------------------------------------------------------------------------------------------------------------------------"
|
|
)
|
|
# iterate over pairs
|
|
for pair in pairs:
|
|
# print pair parameters
|
|
print(
|
|
f"{pair['tags'][0]} {pair['tags'][1]} {pair['Ruc']} d [Ang] {pair['dist']}"
|
|
)
|
|
# print magnetic parameters
|
|
print("Isotropic: ", pair["J_iso"])
|
|
print("DMI: ", pair["D"])
|
|
print("Symmetric-anisotropy: ", pair["J_S"])
|
|
print("J: ", pair["J"].flatten())
|
|
print("Energies for debugging: ")
|
|
pprint(np.array(pair["energies"]))
|
|
print(
|
|
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)"
|
|
)
|
|
o1, o2, o3 = pair["energies"]
|
|
pprint(np.array([o2[-1], o3[0], o1[0]]))
|
|
print("Test J_xx = E(y,z) = E(z,y)")
|
|
print(o2[-1], o3[-1])
|
|
print("")
|
|
|
|
print(
|
|
"================================================================================================================================================================"
|
|
)
|
|
|
|
|
|
def print_runtime_information(times):
|
|
print("Runtime information: ")
|
|
print(f"Total runtime: {times['end_time'] - times['start_time']} s")
|
|
print(
|
|
"----------------------------------------------------------------------------------------------------------------------------------------------------------------"
|
|
)
|
|
print(f"Initial setup: {times['setup_time'] - times['start_time']} s")
|
|
print(
|
|
f"Hamiltonian conversion and XC field extraction: {times['H_and_XCF_time'] - times['setup_time']:.3f} s"
|
|
)
|
|
print(
|
|
f"Pair and site datastructure creatrions: {times['site_and_pair_dictionaries_time'] - times['H_and_XCF_time']:.3f} s"
|
|
)
|
|
print(
|
|
f"k set cration and distribution: {times['k_set_time'] - times['site_and_pair_dictionaries_time']:.3f} s"
|
|
)
|
|
print(
|
|
f"Rotating XC potential: {times['reference_rotations_time'] - times['k_set_time']:.3f} s"
|
|
)
|
|
print(
|
|
f"Greens function inversion: {times['green_function_inversion_time'] - times['reference_rotations_time']:.3f} s"
|
|
)
|
|
print(
|
|
f"Calculate energies and magnetic components: {times['end_time'] - times['green_function_inversion_time']:.3f} s"
|
|
)
|