et209d
/
grogu
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

moved useful to utils, commented a lot and added anisotropy calculation

Browse Source
class-solution
Daniel Pozsar 4 months ago
parent 4b3a6de669
commit 1e8cea8729
2 changed files with 192 additions and 328 deletions
Show Stats Download Patch File Download Diff File

294
src/grogu_magn/jij.py
Unescape View File

@ -1,294 +0,0 @@
# 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.
if __name__ == "__main__":
import argparse
import sys
from itertools import permutations, product
from timeit import default_timer as timer
import numpy as np
import numpy.linalg as nl
import sisl
import tqdm
from mpi4py import MPI
from scipy.special import roots_legendre
from grogu_magn.useful import hsk, make_atran, make_contour, make_kset
start = timer()
# Some input parsing
parser = argparse.ArgumentParser()
parser.add_argument(
"--kset",
dest="kset",
default=2,
type=int,
help="k-space resolution of Jij calculation",
)
parser.add_argument(
"--kdirs",
dest="kdirs",
default="xyz",
help="Definition of k-space dimensionality",
)
parser.add_argument(
"--eset",
dest="eset",
default=42,
type=int,
help="Number of energy points on the contour",
)
parser.add_argument(
"--eset-p",
dest="esetp",
default=10,
type=int,
help="Parameter tuning the distribution on the contour",
)
parser.add_argument("--input", dest="infile", required=True, help="Input file name")
parser.add_argument(
"--output", dest="outfile", required=True, help="Output file name"
)
parser.add_argument(
"--Ebot",
dest="Ebot",
default=-20.0,
type=float,
help="Bottom energy of the contour",
)
parser.add_argument(
"--npairs",
dest="npairs",
default=1,
type=int,
help="Number of unitcell pairs in each direction for Jij calculation",
)
parser.add_argument(
"--adirs", dest="adirs", default=False, help="Definition of pair directions"
)
parser.add_argument(
"--use-tqdm",
dest="usetqdm",
default="not",
help="Use tqdm for progressbars or not",
)
parser.add_argument(
"--cutoff",
dest="cutoff",
default=100.0,
type=float,
help="Real space cutoff for pair generation in Angs",
)
parser.add_argument(
"--pairfile",
dest="pairfile",
default=False,
help="File to read pair information",
)
args = parser.parse_args()
# ----------------------------------------------------------------------
# MPI init
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
root_node = 0
if rank == root_node:
print("Number of nodes in the parallel cluster: ", size)
# ----------------------------------------------------------------------
# importing the necessary structures from SIESTA output
dat = sisl.get_sile(args.infile)
dh = dat.read_hamiltonian()
# update datastructure of the hamiltonian
# this is needed for quick Hk building
dh.hup = (
dh.tocsr(0)
.toarray()
.reshape(dh.no, dh.n_s, dh.no)
.transpose(0, 2, 1)
.astype("complex128")
)
dh.hdo = (
dh.tocsr(1)
.toarray()
.reshape(dh.no, dh.n_s, dh.no)
.transpose(0, 2, 1)
.astype("complex128")
)
dh.sov = (
dh.tocsr(2)
.toarray()
.reshape(dh.no, dh.n_s, dh.no)
.transpose(0, 2, 1)
.astype("complex128")
)
# ----------------------------------------------------------------------
# generate k space sampling
kset = make_kset(dirs=args.kdirs, NUMK=args.kset)
wk = 1 / len(kset) # weight of a kpoint in BZ integral
kpcs = np.array_split(kset, size)
if "k" in args.usetqdm:
kpcs[root_node] = tqdm.tqdm(kpcs[root_node], desc="k loop")
# ----------------------------------------------------------------------
# define pairs
if args.pairfile:
# if pair file is specified read in pair information from pairfile
if rank == root_node:
# read in pair on root node file in a format of five integer columns
# first two integers define sublattice
# second three define distance vectors of unitcells
dummy = np.loadtxt(args.pairfile, dtype=int)
atran = [
(dummy[p, 0], dummy[p, 1], [dummy[p, 2], dummy[p, 3], dummy[p, 4]])
for p in range(len(dummy))
]
else:
atran = None
# syncronize atran over all nodes
atran = comm.bcast(atran, root=root_node)
else:
# if pairfile is not specified generate pair as defined by npairs and adirs
# Take pair directions for k directions if adirs is not set
args.adirs = args.adirs if args.adirs else args.kdirs
# definition of pairs in terms of integer coordinates refering
# to unicell distances and atomic positions
atran = make_atran(len(dh.atoms), args.adirs, dist=args.npairs)
pairs = []
for i, j, uc in atran:
if nl.norm(np.dot(uc, dh.cell)) < args.cutoff:
pairs.append(
dict(
offset=uc, # lattice vector offset between the unitcells the two atoms are
aiij=[i, j], # indecies of the atoms in the unitcell
noij=[
dh.atoms.orbitals[i],
dh.atoms.orbitals[j],
], # number of orbitals on the appropriate atoms
slij=[
slice(
*(lambda x: [min(x), max(x) + 1])(dh.a2o(i, all=True))
), # slices for
slice(*(lambda x: [min(x), max(x) + 1])(dh.a2o(j, all=True))),
], # appropriate orbitals
rirj=[
dh.axyz()[i],
dh.axyz()[j],
], # real space vectors of atoms in the unit cell
Rij=np.dot(
uc, dh.cell
), # real space distance vector between unit cells
rij=np.dot(uc, dh.cell)
- dh.axyz()[i]
+ dh.axyz()[j], # real space vector between atoms
Jijz=[], # in this empty list are we going to gather the integrad of the energy integral
Jij=0, # the final results of the calculation are going to be here on the root node
)
)
if rank == root_node:
print("Number of pairs beeing calculated: ", len(pairs))
comm.Barrier()
# ----------------------------------------------------------------------
# make energy contour
# we are working in eV now !
# and sisil shifts E_F to 0 !
cont = make_contour(emin=args.Ebot, enum=args.eset, p=args.esetp)
eran = cont.ze
# ----------------------------------------------------------------------
# generating onsite matrix and overalp elements of all the atoms in the unitcell
# onsite of the origin supercell
orig_indx = np.arange(0, dh.no) + dh.sc_index([0, 0, 0]) * dh.no
# spin up
uc_up = dh.tocsr(dh.UP)[:, orig_indx].toarray()
# spin down
uc_down = dh.tocsr(dh.DOWN)[:, orig_indx].toarray()
Hs = []
# get number of atoms in the unit cell
for i in range(len(dh.atoms)):
at_indx = dh.a2o(i, all=True)
Hs.append(uc_up[:, at_indx][at_indx, :] - uc_down[:, at_indx][at_indx, :])
# ----------------------------------------------------------------------
# sampling the integrand on the contour and the BZ
for pair in pairs:
noi, noj = pair["noij"]
pair["Guij"] = np.zeros((args.eset, noi, noj), dtype="complex128")
pair["Gdji"] = np.zeros((args.eset, noj, noi), dtype="complex128")
pair["Guij_tmp"] = np.zeros((args.eset, noi, noj), dtype="complex128")
pair["Gdji_tmp"] = np.zeros((args.eset, noj, noi), dtype="complex128")
for k in kpcs[rank]:
HKU, HKD, SK = hsk(dh, k)
Gku = nl.inv(SK * eran.reshape(args.eset, 1, 1) - HKU)
Gkd = nl.inv(SK * eran.reshape(args.eset, 1, 1) - HKD)
for pair in pairs:
phase = np.exp(1j * np.dot(np.dot(k, dh.rcell), pair["Rij"]))
si, sj = pair["slij"]
pair["Guij_tmp"] += Gku[:, si, sj] * phase * wk
pair["Gdji_tmp"] += Gkd[:, sj, si] / phase * wk
# summ reduce partial results of mpi nodes
for pair in pairs:
comm.Reduce(pair["Guij_tmp"], pair["Guij"], root=root_node)
comm.Reduce(pair["Gdji_tmp"], pair["Gdji"], root=root_node)
if rank == root_node:
for pair in pairs:
i, j = pair["aiij"]
# The Szunyogh-Lichtenstein formula
pair["Jijz"] = np.trace(
(Hs[i] @ pair["Guij"]) @ (Hs[j] @ pair["Gdji"]), axis1=1, axis2=2
)
# evaluation of the contour integral
pair["Jij"] = np.trapz(np.imag(pair["Jijz"] * cont.we) / (2 * np.pi))
end = timer()
# ----------------------------------------------------------------------
# and saveing output of the calculation
np.savetxt(
args.outfile,
np.array(
[
[nl.norm(p["rij"]), p["Jij"] * sisl.unit_convert("eV", "Ry") * 1000]
+ p["aiij"]
+ list(p["offset"])
+ list(p["rij"])
for p in pairs
],
dtype=object,
),
header=str(args)
+ "\nnumber of cores = "
+ str(size)
+ "\ntime of calculation = "
+ str(end - start)
+ "\nnorm(rij),Jij[mRy],aiij,offset,rij",
fmt="%s",
)

226
src/grogu_magn/useful.py → src/grogu_magn/utils.py
Unescape View File

@ -33,15 +33,29 @@ tau_0 = np.array([[1, 0], [0, 1]])
# define some useful functions
def hsk(H, ss, sc_off, k=(0, 0, 0)):
"""Speed up Hk and Sk generation.
Calculates the Hamiltonian and the Overlap matrix at a given k point. It is faster that the sisl version.
Args:
H (np.array_like): Hamiltonian in spin box form
ss (np.array_like): Overlap matrix in spin box form
sc_off (list): supercell indexes of the Hamiltonian
k (tuple, optional): The k point where the matrices are set up. Defaults to (0, 0, 0).
Returns:
np.array_like: Hamiltonian at the given k point
np.array_like: Overlap matrix at the given k point
"""
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 two conversion lines are from the sisl source
k = np.asarray(k, np.float64)
k.shape = (-1,)
# this generates the list of phases
phases = np.exp(-1j * 2 * np.pi * k @ sc_off.T)
# phases applied to the hamiltonian
HK = np.einsum("abc,a->bc", H, phases)
SK = np.einsum("abc,a->bc", ss, phases)
@ -49,10 +63,20 @@ def hsk(H, ss, sc_off, k=(0, 0, 0)):
def make_contour(emin=-20, emax=0.0, enum=42, p=150):
"""
A more sophisticated contour generator
"""
"""A more sophisticated contour generator.
Calculates the parameters for the complex contour integral.
Args:
emin (int, optional): Energy minimum of the contour. Defaults to -20.
emax (float, optional): Energy maximum of the contour. Defaults to 0.0, so the Fermi level.
enum (int, optional): Number of sample points along the contour. Defaults to 42.
p (int, optional): Shape parameter that describes the distribution of the sample points. Defaults to 150.
Returns:
ccont: ccont
Contains all the information for the contour integral. Should clarify later...
"""
x, wl = roots_legendre(enum)
R = (emax - emin) / 2
z0 = (emax + emin) / 2
@ -64,8 +88,8 @@ def make_contour(emin=-20, emax=0.0, enum=42, p=150):
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
class ccont:
pass
cont = ccont()
@ -79,11 +103,20 @@ def make_contour(emin=-20, emax=0.0, enum=42, p=150):
def make_kset(dirs="xyz", NUMK=20):
"""
Simple k-grid generator. Depending on the value of the dirs
"""Simple k-grid generator to sample the Brillouin zone.
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 dirs argument does not contain either of x, y or z
a kset of a single k-pont at the origin is returned. The total number of k points is the NUMK**(dimensions)
Args:
dirs (str, optional): Directions of the k points in the Brillouin zone. They are the three lattice vectors. Defaults to "xyz".
NUMK (int, optional): The number of k points in a direction. Defaults to 20.
Returns:
np.array_like: An array of k points that uniformly sample the Brillouin zone in the given directions.
"""
if not (sum([d in dirs for d in "xyz"])):
return np.array([[0, 0, 0]])
@ -103,31 +136,66 @@ def make_kset(dirs="xyz", NUMK=20):
def commutator(a, b):
"Shorthand for commutator"
"""Shorthand for commutator.
Commutator of two matrices in the mathematical sense.
Args:
a (np.array_like): The first matrix.
b (np.array_like): The second matrix
Returns:
np.array_like: The commutator of a and b.
"""
return a @ b - b @ a
def tau_u(u):
"""Pauli matrix in direction u.
Returns the vector u in the basis of the Pauli matrices.
Args:
u (list or np.array_like): The direction
Returns:
np.array_like: Arbitrary direction in the base of the Pauli matrices.
"""
Pauli matrix in direction u.
"""
u = u / np.linalg.norm(u) # u is force to be of unit length
# u is force to be of unit length
u = u / np.linalg.norm(u)
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.
"""Definition for the cross-product matrix.
It acts as a cross product with vector u.
Args:
u (list or np.array_like): The second vector in the cross product
Returns:
np.array_like: The matrix that represents teh cross product with a vector.
"""
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.
Args:
theta (float): The angle of rotation.
u (np.array_like): The rotation axis.
eps (float, optional): Cutoff for small elements in the resulting matrix. Defaults to 1e-10.
Returns:
np.array_like: The rotation matrix.
"""
Definition of rotation matrix with angle theta around direction u.
"""
u = u / np.linalg.norm(u)
M = (
@ -135,25 +203,38 @@ def RotM(theta, u, eps=1e-10):
+ np.sin(theta) * crossM(u)
+ (1 - np.cos(theta)) * np.outer(u, u)
)
M[abs(M) < eps] = 0.0 # kill off small numbers
# kill off small numbers
M[abs(M) < eps] = 0.0
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.
Args:
a (np.array_like): First vector.
b (np.array_like): Second vector.
eps (float, optional): Cutoff for small elements in the resulting matrix. Defaults to 1e-10.
Returns:
np.array_like: The rotation matrix with the above property.
"""
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
# kill off small numbers
M[abs(M) < eps] = 0.0
return M
def spin_tracer(M):
"""
Spin tracer utility.
"""Spin tracer utility.
This takes an operator with the orbital-spin sequence:
orbital 1 up,
orbital 1 down,
@ -161,8 +242,14 @@ def spin_tracer(M):
orbital 2 down,
that is in the SPIN-BOX representation,
and extracts orbital dependent Pauli traces.
"""
Args:
M (np.array_like): Traceble matrix.
Returns:
dict: It contains the traced matrix with "x", "y", "z" and "c".
"""
M11 = M[0::2, 0::2]
M12 = M[0::2, 1::2]
M21 = M[1::2, 0::2]
@ -178,11 +265,15 @@ def spin_tracer(M):
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
"""Function to define orbital indexes of a given magnetic entity.
Args:
dh (sisl.physics.Hamiltonian): Hamiltonian
atom (integer or list of integers, optional): Defining atom (or atoms) in the unit cell forming the magnetic entity. Defaults to None.
l (integer, optional): Defining the angular momentum channel. Defaults to None.
Returns:
list: The orbital indexes of the given magnetic entity.
"""
# case where we deal with more than one atom defining the magnetic entity
if type(atom) == list:
@ -215,7 +306,37 @@ def blow_up_orbindx(orb_indices):
return np.array([[2 * o, 2 * o + 1] for o in orb_indices]).flatten()
def calculate_anisotropy_tensor(mag_ent, eps=1e-10):
"""_summary_
Args:
mag_ent (_type_): _description_
Returns:
_type_: _description_
"""
energies = mag_ent["energies"]
Kxx = energies[1, 1] - energies[1, 0]
Kyy = energies[0, 1] - energies[0, 0]
Kzz = 0
calculated_diff = Kyy - Kxx
expected_diff = energies[2, 0] - energies[2, 1]
consistency_check = abs(calculated_diff - expected_diff) < eps
return Kxx, Kyy, Kzz, consistency_check
def calculate_exchange_tensor(pair):
"""_summary_
Args:
pair (_type_): _description_
Returns:
_type_: _description_
"""
energies = pair["energies"]
# Initialize output arrays
J = np.zeros((3, 3))
@ -260,6 +381,11 @@ def calculate_exchange_tensor(pair):
def print_parameters(simulation_parameters):
"""_summary_
Args:
simulation_parameters (_type_): _description_
"""
print(
"================================================================================================================================================================"
)
@ -299,6 +425,12 @@ def print_parameters(simulation_parameters):
def print_atoms_and_pairs(magnetic_entities, pairs):
"""_summary_
Args:
magnetic_entities (_type_): _description_
pairs (_type_): _description_
"""
print("Atomic information: ")
print(
"----------------------------------------------------------------------------------------------------------------------------------------------------------------"
@ -317,6 +449,27 @@ def print_atoms_and_pairs(magnetic_entities, pairs):
print(f"{tag} {xyz[0]} {xyz[1]} {xyz[2]}")
print("")
print(
"================================================================================================================================================================"
)
print("Anisotropy [meV]")
print(
"----------------------------------------------------------------------------------------------------------------------------------------------------------------"
)
print("Magnetic entity x [Ang] y [Ang] z [Ang]")
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("Consistency check passed: ", mag_ent["K_consistency"])
print("Anisotropy diag: ", mag_ent["K"])
print("")
print(
"================================================================================================================================================================"
)
@ -356,6 +509,11 @@ def print_atoms_and_pairs(magnetic_entities, pairs):
def print_runtime_information(times):
"""_summary_
Args:
times (_type_): _description_
"""
print("Runtime information: ")
print(f"Total runtime: {times['end_time'] - times['start_time']} s")
print(
Write Preview
Loading…
Cancel
Save