From 1e8cea8729d64b4961de026c875028359e1e95e1 Mon Sep 17 00:00:00 2001 From: Daniel Pozsar Date: Fri, 8 Nov 2024 08:54:16 +0100 Subject: [PATCH] moved useful to utils, commented a lot and added anisotropy calculation --- src/grogu_magn/jij.py | 294 ------------------------- src/grogu_magn/{useful.py => utils.py} | 226 ++++++++++++++++--- 2 files changed, 192 insertions(+), 328 deletions(-) delete mode 100644 src/grogu_magn/jij.py rename src/grogu_magn/{useful.py => utils.py} (66%) diff --git a/src/grogu_magn/jij.py b/src/grogu_magn/jij.py deleted file mode 100644 index c91c289..0000000 --- a/src/grogu_magn/jij.py +++ /dev/null @@ -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", - ) diff --git a/src/grogu_magn/useful.py b/src/grogu_magn/utils.py similarity index 66% rename from src/grogu_magn/useful.py rename to src/grogu_magn/utils.py index 7bd9911..0b94e06 100644 --- a/src/grogu_magn/useful.py +++ b/src/grogu_magn/utils.py @@ -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 + # just an empty container class class ccont: - # just an empty container class 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(