diff --git a/src/grogupy/core.py b/src/grogupy/core.py index b8c8a9c..291a659 100644 --- a/src/grogupy/core.py +++ b/src/grogupy/core.py @@ -21,15 +21,57 @@ import numpy as np from numpy.linalg import inv -from grogupy.utils import blow_up_orbindx, commutator, parse_magnetic_entity +from grogupy.magnetism import blow_up_orbindx, parse_magnetic_entity +from grogupy.utils import commutator def parallel_Gk(HK, SK, eran, eset): + """Calculates the Greens function by inversion. + + It calculates the Greens function on all the energy levels at the same time. + + Args: + HK : (NO, NO), np.array_like + Hamiltonian at a given k point + SK : (NO, NO), np.array_like + Overlap Matrix at a given k point + eran : (eset) np.array_like + Energy sample along the contour + eset : int + Number of energy samples along the contour + + Returns: + Gk : (eset, NO, NO), np.array_like + Green's function at a given k point + """ + + # Calculates the Greens function on all the energy levels return inv(SK * eran.reshape(eset, 1, 1) - HK) def sequential_GK(HK, SK, eran, eset): + """Calculates the Greens function by inversion. + + It calculates sequentially over the energy levels. + + Args: + HK : (NO, NO), np.array_like + Hamiltonian at a given k point + SK : (NO, NO), np.array_like + Overlap Matrix at a given k point + eran : (eset) np.array_like + Energy sample along the contour + eset : int + Number of energy samples along the contour + + Returns: + Gk : (eset, NO, NO), np.array_like + Green's function at a given k point + """ + + # creates an empty holder Gk = np.zeros(shape=(eset, HK.shape[0], HK.shape[1]), dtype="complex128") + # fills the holder sequentially by the Greens function on a given energy for j in range(eset): Gk[j] = inv(SK * eran[j] - HK) @@ -37,14 +79,19 @@ def sequential_GK(HK, SK, eran, eset): def calc_Vu(H, Tu): - """_summary_ + """Calculates the local perturbation in case of a spin rotation. Args: - H (_type_): _description_ - Tu (_type_): _description_ + H : (NO, NO) np.array_like + Hamiltonian + Tu : (NO, NO) array_like + Rotation around u Returns: - _type_: _description_ + Vu1 : (NO, NO) np.array_like + First order perturbed matrix + Vu2 : (NO, NO) np.array_like + Second order perturbed matrix """ Vu1 = 1j / 2 * commutator(H, Tu) # equation 100 @@ -54,15 +101,24 @@ def calc_Vu(H, Tu): def remove_clutter_for_save(pairs, magnetic_entities): - """_summary_ + """Removes unimportant data from the dictionaries. + + It is used before saving to throw away data that + is not needed for post processing. Args: - pairs (_type_): _description_ - magnetic_entities (_type_): _description_ + pairs : dict + Contains all the pair information + magnetic_entities : dict + Contains all the magnetic entity information Returns: - _type_: _description_ + pairs : dict + Contains all the reduced pair information + magnetic_entities : dict + Contains all the reduced magnetic entity information """ + # remove clutter from magnetic entities and pair information for pair in pairs: del pair["Gij"] @@ -79,14 +135,22 @@ def remove_clutter_for_save(pairs, magnetic_entities): def build_hh_ss(dh): - """_summary_ + """It builds the Hamiltonian and Overlap matrix from the sisl.dh class. + + It restructures the data in the SPIN BOX representation, where NS is + the number of supercells and NO is the number of orbitals. Args: - dh (_type_): _description_ + dh : sisl.physics.Hamiltonian + Hamiltonian read in by sisl Returns: - _type_: _description_ + hh : (NS, NO, NO) np.array_like + Hamiltonian in SPIN BOX representation + ss : (NS, NO, NO) np.array_like + Overlap matrix in SPIN BOX representation """ + NO = dh.no # shorthand for number of orbitals in the unit cell # preprocessing Hamiltonian and overlap matrix elements @@ -142,23 +206,35 @@ def build_hh_ss(dh): ] ) - return hh, ss, NO + return hh, ss def setup_pairs_and_magnetic_entities( magnetic_entities, pairs, dh, simulation_parameters ): - """_summary_ + """It creates the complete structure of the dictionaries and fills some basic data. + + It creates orbital indexes, spin box indexes, coordinates and tags for magnetic entities. + Furthermore it creates the structures for the energies, the perturbed potentials and + the Greens function calculation. It dose the same for the pairs. Args: - magnetic_entities (_type_): _description_ - pairs (_type_): _description_ - dh (_type_): _description_ - simulation_parameters (_type_): _description_ + pairs : dict + Contains the initial pair information + magnetic_entities : dict + Contains the initial magnetic entity information + dh : sisl.physics.Hamiltonian + Hamiltonian read in by sisl + simulation_parameters : dict + A set of parameters from the simulation Returns: - _type_: _description_ + pairs : dict + Contains the initial information and the complete structure + magnetic_entities : dict + Contains the initial information and the complete structure """ + # for every site we have to store 3 Greens function (and the associated _tmp-s) in the 3 reference directions for mag_ent in magnetic_entities: parsed = parse_magnetic_entity(dh, **mag_ent) # parse orbital indexes @@ -169,11 +245,14 @@ def setup_pairs_and_magnetic_entities( # if orbital is not set use all if "l" not in mag_ent.keys(): mag_ent["l"] = "all" + + # tag creation for one atom if isinstance(mag_ent["atom"], int): mag_ent["tags"] = [ f"[{mag_ent['atom']}]{dh.atoms[mag_ent['atom']].tag}({mag_ent['l']})" ] mag_ent["xyz"] = [dh.xyz[mag_ent["atom"]]] + # tag creation for more atoms if isinstance(mag_ent["atom"], list): mag_ent["tags"] = [] mag_ent["xyz"] = [] @@ -187,9 +266,8 @@ def setup_pairs_and_magnetic_entities( # calculate size for Greens function generation spin_box_shape = len(mag_ent["spin_box_indices"]) - mag_ent["energies"] = ( - [] - ) # we will store the second order energy derivations here + # we will store the second order energy derivations here + mag_ent["energies"] = [] # These will be the perturbed potentials from eq. 100 mag_ent["Vu1"] = [] # so they are independent in memory @@ -197,7 +275,7 @@ def setup_pairs_and_magnetic_entities( mag_ent["Gii"] = [] # Greens function mag_ent["Gii_tmp"] = [] # Greens function for parallelization - for i in simulation_parameters["ref_xcf_orientations"]: + for _ in simulation_parameters["ref_xcf_orientations"]: # Rotations for every quantization axis mag_ent["Vu1"].append([]) mag_ent["Vu2"].append([]) @@ -227,6 +305,7 @@ def setup_pairs_and_magnetic_entities( # calculate size for Greens function generation spin_box_shape_i = len(magnetic_entities[pair["ai"]]["spin_box_indices"]) spin_box_shape_j = len(magnetic_entities[pair["aj"]]["spin_box_indices"]) + # tag generation pair["tags"] = [] for mag_ent in [magnetic_entities[pair["ai"]], magnetic_entities[pair["aj"]]]: tag = "" @@ -253,7 +332,7 @@ def setup_pairs_and_magnetic_entities( pair["Gji"] = [] pair["Gij_tmp"] = [] # Greens function for parallelization pair["Gji_tmp"] = [] - for i in simulation_parameters["ref_xcf_orientations"]: + for _ in simulation_parameters["ref_xcf_orientations"]: # Greens functions for every quantization axis pair["Gij"].append( np.zeros( @@ -284,14 +363,19 @@ def setup_pairs_and_magnetic_entities( def onsite_projection(matrix, idx1, idx2): - """_summary_ + """It produces the slices of a matrix for the on site projection. + + The slicing is along the last two axes as these contains the orbital indexing. Args: - matrix (_type_): _description_ - idx (_type_): _description_ + matrix : (..., :, :) np.array_like + Some matrix + idx : np.array_like + The indexes of the orbitals Returns: - _type_: _description_ + np.array_like + Reduced matrix based on the projection """ return matrix[..., idx1, :][..., idx2] diff --git a/src/grogupy/grogu.py b/src/grogupy/grogu.py index 5964f88..208a946 100644 --- a/src/grogupy/grogu.py +++ b/src/grogupy/grogu.py @@ -19,457 +19,435 @@ # SOFTWARE. import os - -os.environ["OMP_NUM_THREADS"] = "1" # export OMP_NUM_THREADS=1 -os.environ["OPENBLAS_NUM_THREADS"] = "1" # export OPENBLAS_NUM_THREADS=1 -os.environ["MKL_NUM_THREADS"] = "1" # export MKL_NUM_THREADS=1 -os.environ["VECLIB_MAXIMUM_THREADS"] = "1" # export VECLIB_MAXIMUM_THREADS=1 -os.environ["NUMEXPR_NUM_THREADS"] = "1" # export NUMEXPR_NUM_THREADS=1 - - -from timeit import default_timer as timer - -# runtime information -times = dict() -times["start_time"] = timer() - - import warnings from sys import getsizeof +from timeit import default_timer as timer +import numpy as np import sisl from mpi4py import MPI -from src.grogupy import * - -# input output stuff -###################################################################### -###################################################################### -###################################################################### - -infile = ( - "/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf" -) -outfile = "./Fe3GeTe2_notebook" - -magnetic_entities = [ - dict(atom=3, l=2), - dict(atom=4, l=2), - dict(atom=5, l=2), -] -pairs = [ - dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])), - dict(ai=0, aj=2, Ruc=np.array([0, 0, 0])), - dict(ai=1, aj=2, Ruc=np.array([0, 0, 0])), - dict(ai=0, aj=2, Ruc=np.array([-1, -1, 0])), - dict(ai=1, aj=2, Ruc=np.array([-1, -1, 0])), - dict(ai=0, aj=2, Ruc=np.array([-1, 0, 0])), - dict(ai=1, aj=2, Ruc=np.array([-1, 0, 0])), - dict(ai=1, aj=2, Ruc=np.array([-2, 0, 0])), - dict(ai=1, aj=2, Ruc=np.array([-3, 0, 0])), -] -simulation_parameters = default_args -simulation_parameters["infile"] = infile -simulation_parameters["outfile"] = outfile -simulation_parameters["kset"] = 20 -simulation_parameters["kdirs"] = "xy" -simulation_parameters["eset"] = 600 -simulation_parameters["esetp"] = 10000 - -fdf = sisl.io.fdfSileSiesta("input.fdf") -fdf.get("XCF_Rotation") - -###################################################################### -###################################################################### -###################################################################### - - -# MPI parameters -comm = MPI.COMM_WORLD -size = comm.Get_size() -rank = comm.Get_rank() -root_node = 0 - -# include parallel size in simulation parameters -simulation_parameters["parallel_size"] = size - -# check versions for debugging -if rank == root_node: - try: - print("sisl version: ", sisl.__version__) - except: - print("sisl version unknown.") - try: - print("numpy version: ", np.__version__) - except: - print("numpy version unknown.") - -# rename outfile -if not simulation_parameters["outfile"].endswith(".pickle"): - simulation_parameters["outfile"] += ".pickle" - -# if ebot is not given put it 0.1 eV under the smallest energy -if simulation_parameters["ebot"] is None: - try: - eigfile = simulation_parameters["infile"][:-3] + "EIG" - simulation_parameters["ebot"] = read_siesta_emin(eigfile) - 0.1 - simulation_parameters["automatic_ebot"] = True - except: - print("Could not determine ebot.") - print("Parameter was not given and .EIG file was not found.") - -# read sile -fdf = sisl.get_sile(simulation_parameters["infile"]) - -# read in hamiltonian -dh = fdf.read_hamiltonian() - -# read unit cell vectors -simulation_parameters["cell"] = fdf.read_geometry().cell - -# unit cell index -uc_in_sc_idx = dh.lattice.sc_index([0, 0, 0]) - -if rank == root_node: - print("\n\n\n\n\n") - print( - "#################################################################### JOB INFORMATION ###########################################################################" - ) - print_job_description(simulation_parameters) - print( - "################################################################################################################################################################" - ) - print("\n\n\n\n\n") +from grogupy import * - times["setup_time"] = timer() - print(f"Setup done. Elapsed time: {times['setup_time']} s") - print( - "================================================================================================================================================================" - ) +def main(): -# reformat Hamltonian and Overlap matrix for manipulations -hh, ss, NO = build_hh_ss(dh) - -# symmetrizing Hamiltonian and Overlap matrix to make them hermitian -for i in range(dh.lattice.sc_off.shape[0]): - j = dh.lattice.sc_index(-dh.lattice.sc_off[i]) - h1, h1d = hh[i], hh[j] - hh[i], hh[j] = (h1 + h1d.T.conj()) / 2, (h1d + h1.T.conj()) / 2 - s1, s1d = ss[i], ss[j] - ss[i], ss[j] = (s1 + s1d.T.conj()) / 2, (s1d + s1.T.conj()) / 2 - -# identifying TRS and TRB parts of the Hamiltonian -TAUY = np.kron(np.eye(NO), tau_y) -hTR = np.array([TAUY @ hh[i].conj() @ TAUY for i in range(dh.lattice.nsc.prod())]) -hTRS = (hh + hTR) / 2 -hTRB = (hh - hTR) / 2 - -# extracting the exchange field -traced = [spin_tracer(hTRB[i]) for i in range(dh.lattice.nsc.prod())] # equation 77 -XCF = np.array( - [ - np.array([f["x"] / 2 for f in traced]), - np.array([f["y"] / 2 for f in traced]), - np.array([f["z"] / 2 for f in traced]), - ] -) + # constrain numpy in parallel run + os.environ["OMP_NUM_THREADS"] = "1" # export OMP_NUM_THREADS=1 + os.environ["OPENBLAS_NUM_THREADS"] = "1" # export OPENBLAS_NUM_THREADS=1 + os.environ["MKL_NUM_THREADS"] = "1" # export MKL_NUM_THREADS=1 + os.environ["VECLIB_MAXIMUM_THREADS"] = "1" # export VECLIB_MAXIMUM_THREADS=1 + os.environ["NUMEXPR_NUM_THREADS"] = "1" # export NUMEXPR_NUM_THREADS=1 -# check if exchange field has scalar part -max_xcfs = abs(np.array(np.array([f["c"] / 2 for f in traced]))).max() -if max_xcfs > 1e-12: - warnings.warn( - f"Exchange field has non negligible scalar part. Largest value is {max_xcfs}" - ) + # runtime information + times = dict() + times["start_time"] = timer() -if rank == root_node: - times["H_and_XCF_time"] = timer() - print( - f"Hamiltonian and exchange field rotated. Elapsed time: {times['H_and_XCF_time']} s" - ) - print( - "================================================================================================================================================================" - ) + # input output stuff + ###################################################################### + ###################################################################### + ###################################################################### -# initialize pairs and magnetic entities based on input information -pairs, magnetic_entities = setup_pairs_and_magnetic_entities( - magnetic_entities, pairs, dh, simulation_parameters -) + infile = "/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf" + outfile = "./Fe3GeTe2_notebook" -if rank == root_node: - times["site_and_pair_dictionaries_time"] = timer() - print( - f"Site and pair dictionaries created. Elapsed time: {times['site_and_pair_dictionaries_time']} s" - ) - print( - "================================================================================================================================================================" - ) + magnetic_entities = [ + dict(atom=3, l=2), + dict(atom=4, l=2), + dict(atom=5, l=2), + ] + pairs = [ + dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])), + dict(ai=0, aj=2, Ruc=np.array([0, 0, 0])), + dict(ai=1, aj=2, Ruc=np.array([0, 0, 0])), + dict(ai=0, aj=2, Ruc=np.array([-1, -1, 0])), + dict(ai=1, aj=2, Ruc=np.array([-1, -1, 0])), + dict(ai=0, aj=2, Ruc=np.array([-1, 0, 0])), + dict(ai=1, aj=2, Ruc=np.array([-1, 0, 0])), + dict(ai=1, aj=2, Ruc=np.array([-2, 0, 0])), + dict(ai=1, aj=2, Ruc=np.array([-3, 0, 0])), + ] + simulation_parameters = default_args + simulation_parameters["infile"] = infile + simulation_parameters["outfile"] = outfile + simulation_parameters["kset"] = 20 + simulation_parameters["kdirs"] = "xy" + simulation_parameters["eset"] = 600 + simulation_parameters["esetp"] = 10000 + + ###################################################################### + ###################################################################### + ###################################################################### + + # MPI parameters + comm = MPI.COMM_WORLD + size = comm.Get_size() + rank = comm.Get_rank() + root_node = 0 + + if rank == root_node: + # include parallel size in simulation parameters + simulation_parameters["parallel_size"] = size + + # check versions for debugging + try: + print("sisl version: ", sisl.__version__) + except: + print("sisl version unknown.") + try: + print("numpy version: ", np.__version__) + except: + print("numpy version unknown.") + + # rename outfile + if not simulation_parameters["outfile"].endswith(".pickle"): + simulation_parameters["outfile"] += ".pickle" + + # if ebot is not given put it 0.1 eV under the smallest energy + if simulation_parameters["ebot"] is None: + try: + eigfile = simulation_parameters["infile"][:-3] + "EIG" + simulation_parameters["ebot"] = read_siesta_emin(eigfile) - 0.1 + simulation_parameters["automatic_ebot"] = True + except: + print("Could not determine ebot.") + print("Parameter was not given and .EIG file was not found.") + else: + simulation_parameters["automatic_ebot"] = False -# generate k space sampling -kset = make_kset( - dirs=simulation_parameters["kdirs"], NUMK=simulation_parameters["kset"] -) + # read sile + fdf = sisl.get_sile(simulation_parameters["infile"]) -# generate weights for k points -wkset = np.ones(len(kset)) / len(kset) + # read in hamiltonian + dh = fdf.read_hamiltonian() -# split the k points based on MPI size -kpcs = np.array_split(kset, size) + # read unit cell vectors + simulation_parameters["cell"] = fdf.read_geometry().cell -# use progress bar if available -if rank == root_node and tqdm_imported: - kpcs[root_node] = tqdm(kpcs[root_node], desc="k loop") + # unit cell index + uc_in_sc_idx = dh.lattice.sc_index([0, 0, 0]) -if rank == root_node: - times["k_set_time"] = timer() - print(f"k set created. Elapsed time: {times['k_set_time']} s") - print( - "================================================================================================================================================================" - ) + if rank == root_node: + print("\n\n\n\n\n") + print( + "#################################################################### JOB INFORMATION ###########################################################################" + ) + print_job_description(simulation_parameters) + print( + "################################################################################################################################################################" + ) + print("\n\n\n\n\n") -# this will contain the three Hamiltonian in the -# reference directions needed to calculate the energy -# variations upon rotation -hamiltonians = [] - -# iterate over the reference directions (quantization axes) -for i, orient in enumerate(simulation_parameters["ref_xcf_orientations"]): - # obtain rotated exchange field and Hamiltonian - R = RotMa2b(simulation_parameters["scf_xcf_orientation"], orient["o"]) - rot_XCF = np.einsum("ij,jklm->iklm", R, XCF) - rot_H_XCF = sum( - [np.kron(rot_XCF[i], tau) for i, tau in enumerate([tau_x, tau_y, tau_z])] - ) - rot_H_XCF_uc = rot_H_XCF[uc_in_sc_idx] + times["setup_time"] = timer() + print(f"Setup done. Elapsed time: {times['setup_time']} s") + print( + "================================================================================================================================================================" + ) - # obtain total Hamiltonian with the rotated exchange field - rot_H = hTRS + rot_H_XCF # equation 76 + NO = dh.no # shorthand for number of orbitals in the unit cell + + # reformat Hamltonian and Overlap matrix for manipulations + hh, ss = build_hh_ss(dh) + + # symmetrizing Hamiltonian and Overlap matrix to make them hermitian + for i in range(dh.lattice.sc_off.shape[0]): + j = dh.lattice.sc_index(-dh.lattice.sc_off[i]) + h1, h1d = hh[i], hh[j] + hh[i], hh[j] = (h1 + h1d.T.conj()) / 2, (h1d + h1.T.conj()) / 2 + s1, s1d = ss[i], ss[j] + ss[i], ss[j] = (s1 + s1d.T.conj()) / 2, (s1d + s1.T.conj()) / 2 + + # identifying TRS and TRB parts of the Hamiltonian + TAUY = np.kron(np.eye(NO), tau_y) + hTR = np.array([TAUY @ hh[i].conj() @ TAUY for i in range(dh.lattice.nsc.prod())]) + hTRS = (hh + hTR) / 2 + hTRB = (hh - hTR) / 2 + + # extracting the exchange field + traced = [spin_tracer(hTRB[i]) for i in range(dh.lattice.nsc.prod())] # equation 77 + XCF = np.array( + [ + np.array([f["x"] / 2 for f in traced]), + np.array([f["y"] / 2 for f in traced]), + np.array([f["z"] / 2 for f in traced]), + ] + ) - # store the relevant information of the Hamiltonian - hamiltonians.append(dict(orient=orient["o"], H=rot_H)) + # check if exchange field has scalar part + max_xcfs = abs(np.array(np.array([f["c"] / 2 for f in traced]))).max() + if max_xcfs > 1e-12: + warnings.warn( + f"Exchange field has non negligible scalar part. Largest value is {max_xcfs}" + ) - if simulation_parameters["calculate_charge"]: - hamiltonians[-1]["GS"] = np.zeros( - (simulation_parameters["eset"], rot_H.shape[1], rot_H.shape[2]), - dtype="complex128", + if rank == root_node: + times["H_and_XCF_time"] = timer() + print( + f"Hamiltonian and exchange field rotated. Elapsed time: {times['H_and_XCF_time']} s" ) - hamiltonians[-1]["GS_tmp"] = np.zeros( - (simulation_parameters["eset"], rot_H.shape[1], rot_H.shape[2]), - dtype="complex128", + print( + "================================================================================================================================================================" ) - # these are the rotations (for now) perpendicular to the quantization axis - for u in orient["vw"]: - # section 2.H - Tu = np.kron(np.eye(NO, dtype=int), tau_u(u)) - Vu1, Vu2 = calc_Vu(rot_H_XCF_uc, Tu) - - for mag_ent in magnetic_entities: - idx = mag_ent["spin_box_indices"] - # fill up the perturbed potentials (for now) based on the on-site projections - mag_ent["Vu1"][i].append(onsite_projection(Vu1, idx, idx)) - mag_ent["Vu2"][i].append(onsite_projection(Vu2, idx, idx)) - -if rank == root_node: - times["reference_rotations_time"] = timer() - print( - f"Rotations done perpendicular to quantization axis. Elapsed time: {times['reference_rotations_time']} s" - ) - print( - "================================================================================================================================================================" + # initialize pairs and magnetic entities based on input information + pairs, magnetic_entities = setup_pairs_and_magnetic_entities( + magnetic_entities, pairs, dh, simulation_parameters ) -# provide helpful information to estimate the runtime and memory -# requirements of the Greens function calculations -if rank == root_node: - print("Starting matrix inversions.") - if simulation_parameters["padawan_mode"]: - print("Padawan mode: ") - print(f"Total number of k points: {kset.shape[0]}") - print(f"Number of energy samples per k point: {simulation_parameters['eset']}") - print(f"Total number of directions: {len(hamiltonians)}") + if rank == root_node: + times["site_and_pair_dictionaries_time"] = timer() print( - f"Total number of matrix inversions: {kset.shape[0] * len(hamiltonians) * simulation_parameters['eset']}" + f"Site and pair dictionaries created. Elapsed time: {times['site_and_pair_dictionaries_time']} s" ) print( - f"The shape of the Hamiltonian and the Greens function is {NO}x{NO}={NO*NO}" + "================================================================================================================================================================" + ) + + # generate k space sampling + kset = make_kset( + dirs=simulation_parameters["kdirs"], NUMK=simulation_parameters["kset"] + ) + + # generate weights for k points + wkset = np.ones(len(kset)) / len(kset) + + # split the k points based on MPI size + kpcs = np.array_split(kset, size) + + # use progress bar if available + if rank == root_node and tqdm_imported: + kpcs[root_node] = tqdm(kpcs[root_node], desc="k loop") + + if rank == root_node: + times["k_set_time"] = timer() + print(f"k set created. Elapsed time: {times['k_set_time']} s") + print( + "================================================================================================================================================================" + ) + + # this will contain the three Hamiltonian in the + # reference directions needed to calculate the energy + # variations upon rotation + hamiltonians = [] + + # iterate over the reference directions (quantization axes) + for i, orient in enumerate(simulation_parameters["ref_xcf_orientations"]): + # obtain rotated exchange field and Hamiltonian + R = RotMa2b(simulation_parameters["scf_xcf_orientation"], orient["o"]) + rot_XCF = np.einsum("ij,jklm->iklm", R, XCF) + rot_H_XCF = sum( + [np.kron(rot_XCF[i], tau) for i, tau in enumerate([tau_x, tau_y, tau_z])] ) - # https://stackoverflow.com/questions/70746660/how-to-predict-memory-requirement-for-np-linalg-inv - # memory is O(64 n**2) for complex matrices - memory_size = getsizeof(hamiltonians[0]["H"].base) / 1024 + rot_H_XCF_uc = rot_H_XCF[uc_in_sc_idx] + + # obtain total Hamiltonian with the rotated exchange field + rot_H = hTRS + rot_H_XCF # equation 76 + + # store the relevant information of the Hamiltonian + hamiltonians.append(dict(orient=orient["o"], H=rot_H)) + + # these are the rotations (for now) perpendicular to the quantization axis + for u in orient["vw"]: + # section 2.H + Tu = np.kron(np.eye(NO, dtype=int), tau_u(u)) + Vu1, Vu2 = calc_Vu(rot_H_XCF_uc, Tu) + + for mag_ent in magnetic_entities: + idx = mag_ent["spin_box_indices"] + # fill up the perturbed potentials (for now) based on the on-site projections + mag_ent["Vu1"][i].append(onsite_projection(Vu1, idx, idx)) + mag_ent["Vu2"][i].append(onsite_projection(Vu2, idx, idx)) + + if rank == root_node: + times["reference_rotations_time"] = timer() print( - f"Memory taken by a single Hamiltonian is: {getsizeof(hamiltonians[0]['H'].base) / 1024} KB" + f"Rotations done perpendicular to quantization axis. Elapsed time: {times['reference_rotations_time']} s" ) - print(f"Expected memory usage per matrix inversion: {memory_size * 32} KB") print( - f"Expected memory usage per k point for parallel inversion: {memory_size * len(hamiltonians) * simulation_parameters['eset'] * 32} KB" + "================================================================================================================================================================" ) + + # provide helpful information to estimate the runtime and memory + # requirements of the Greens function calculations + if rank == root_node: + print("Starting matrix inversions.") + if simulation_parameters["padawan_mode"]: + print("Padawan mode: ") + print(f"Total number of k points: {kset.shape[0]}") + print( + f"Number of energy samples per k point: {simulation_parameters['eset']}" + ) + print(f"Total number of directions: {len(hamiltonians)}") + print( + f"Total number of matrix inversions: {kset.shape[0] * len(hamiltonians) * simulation_parameters['eset']}" + ) + print( + f"The shape of the Hamiltonian and the Greens function is {NO}x{NO}={NO*NO}" + ) + # https://stackoverflow.com/questions/70746660/how-to-predict-memory-requirement-for-np-linalg-inv + # memory is O(64 n**2) for complex matrices + memory_size = getsizeof(hamiltonians[0]["H"].base) / 1024 + print( + f"Memory taken by a single Hamiltonian is: {getsizeof(hamiltonians[0]['H'].base) / 1024} KB" + ) + print(f"Expected memory usage per matrix inversion: {memory_size * 32} KB") + print( + f"Expected memory usage per k point for parallel inversion: {memory_size * len(hamiltonians) * simulation_parameters['eset'] * 32} KB" + ) + print( + f"Expected memory usage on root node: {len(np.array_split(kset, size)[0]) * memory_size * len(hamiltonians) * simulation_parameters['eset'] * 32 / 1024} MB" + ) print( - f"Expected memory usage on root node: {len(np.array_split(kset, size)[0]) * memory_size * len(hamiltonians) * simulation_parameters['eset'] * 32 / 1024} MB" + "================================================================================================================================================================" ) - print( - "================================================================================================================================================================" - ) -# MPI barrier -comm.Barrier() - -# make energy contour -cont = make_contour( - emin=simulation_parameters["ebot"], - enum=simulation_parameters["eset"], - p=simulation_parameters["esetp"], -) -eran = cont.ze - -# sampling the integrand on the contour and the BZ -for k in kpcs[rank]: - # weight of k point in BZ integral - wk = wkset[rank] - - # iterate over reference directions - for i, hamiltonian_orientation in enumerate(hamiltonians): - # calculate Hamiltonian and Overlap matrix in a given k point - H = hamiltonian_orientation["H"] - HK, SK = hsk(H, ss, dh.sc_off, k) - - if simulation_parameters["parallel_solver_for_Gk"]: - Gk = parallel_Gk(HK, SK, eran, simulation_parameters["eset"]) - else: - # solve Greens function sequentially for the energies, because of memory bound - Gk = sequential_GK(HK, SK, eran, simulation_parameters["eset"]) + # MPI barrier + comm.Barrier() - # saving this for total charge - if simulation_parameters["calculate_charge"]: - hamiltonian_orientation["GS_tmp"] += Gk @ SK * wk + # make energy contour + cont = make_contour( + emin=simulation_parameters["ebot"], + enum=simulation_parameters["eset"], + p=simulation_parameters["esetp"], + ) + eran = cont.ze - # store the Greens function slice of the magnetic entities - for mag_ent in magnetic_entities: - idx = mag_ent["spin_box_indices"] - mag_ent["Gii_tmp"][i] += onsite_projection(Gk, idx, idx) * wk + # sampling the integrand on the contour and the BZ + for k in kpcs[rank]: + # weight of k point in BZ integral + wk = wkset[rank] - for pair in pairs: - # add phase shift based on the cell difference - phase = np.exp(1j * 2 * np.pi * k @ pair["Ruc"].T) + # iterate over reference directions + for i, hamiltonian_orientation in enumerate(hamiltonians): + # calculate Hamiltonian and Overlap matrix in a given k point + H = hamiltonian_orientation["H"] + HK, SK = hsk(H, ss, dh.sc_off, k) - # get the pair orbital sizes from the magnetic entities - ai = magnetic_entities[pair["ai"]]["spin_box_indices"] - aj = magnetic_entities[pair["aj"]]["spin_box_indices"] + if simulation_parameters["parallel_solver_for_Gk"]: + Gk = parallel_Gk(HK, SK, eran, simulation_parameters["eset"]) + else: + # solve Greens function sequentially for the energies, because of memory bound + Gk = sequential_GK(HK, SK, eran, simulation_parameters["eset"]) # store the Greens function slice of the magnetic entities - pair["Gij_tmp"][i] += onsite_projection(Gk, ai, aj) * phase * wk - pair["Gji_tmp"][i] += onsite_projection(Gk, aj, ai) / phase * wk - -# summ reduce partial results of mpi nodes -for i in range(len(hamiltonians)): - # for total charge calculation - if simulation_parameters["calculate_charge"]: - comm.Reduce(hamiltonians[i]["GS_tmp"], hamiltonians[i]["GS"], root=root_node) - - for mag_ent in magnetic_entities: - comm.Reduce(mag_ent["Gii_tmp"][i], mag_ent["Gii"][i], root=root_node) - - for pair in pairs: - comm.Reduce(pair["Gij_tmp"][i], pair["Gij"][i], root=root_node) - comm.Reduce(pair["Gji_tmp"][i], pair["Gji"][i], root=root_node) - -if rank == root_node: - times["green_function_inversion_time"] = timer() - print( - f"Calculated Greens functions. Elapsed time: {times['green_function_inversion_time']} s" - ) - print( - "================================================================================================================================================================" - ) + for mag_ent in magnetic_entities: + idx = mag_ent["spin_box_indices"] + mag_ent["Gii_tmp"][i] += onsite_projection(Gk, idx, idx) * wk + + for pair in pairs: + # add phase shift based on the cell difference + phase = np.exp(1j * 2 * np.pi * k @ pair["Ruc"].T) + + # get the pair orbital sizes from the magnetic entities + ai = magnetic_entities[pair["ai"]]["spin_box_indices"] + aj = magnetic_entities[pair["aj"]]["spin_box_indices"] + # store the Greens function slice of the magnetic entities + pair["Gij_tmp"][i] += onsite_projection(Gk, ai, aj) * phase * wk + pair["Gji_tmp"][i] += onsite_projection(Gk, aj, ai) / phase * wk -if rank == root_node: - # Calculate total charge - if simulation_parameters["calculate_charge"]: - for hamiltonian in hamiltonians: - GS = hamiltonian["GS"] - traced = np.trace((GS), axis1=1, axis2=2) - simulation_parameters["charges"].append(int_de_ke(traced, cont.we)) - - # iterate over the magnetic entities - for tracker, mag_ent in enumerate(magnetic_entities): - # iterate over the quantization axes - for i, Gii in enumerate(mag_ent["Gii"]): - storage = [] - # iterate over the first and second order local perturbations - for Vu1, Vu2 in zip(mag_ent["Vu1"][i], mag_ent["Vu2"][i]): - # The Szunyogh-Lichtenstein formula - traced = np.trace( - (Vu2 @ Gii + 0.5 * Gii @ Vu1 @ Gii), axis1=1, axis2=2 - ) # this is the on site projection - # evaluation of the contour integral - storage.append(int_de_ke(traced, cont.we)) - # fill up the magnetic entities dictionary with the energies - magnetic_entities[tracker]["energies"].append(storage) - # convert to np array - magnetic_entities[tracker]["energies"] = np.array( - magnetic_entities[tracker]["energies"] + # summ reduce partial results of mpi nodes + for i in range(len(hamiltonians)): + for mag_ent in magnetic_entities: + comm.Reduce(mag_ent["Gii_tmp"][i], mag_ent["Gii"][i], root=root_node) + + for pair in pairs: + comm.Reduce(pair["Gij_tmp"][i], pair["Gij"][i], root=root_node) + comm.Reduce(pair["Gji_tmp"][i], pair["Gji"][i], root=root_node) + + if rank == root_node: + times["green_function_inversion_time"] = timer() + print( + f"Calculated Greens functions. Elapsed time: {times['green_function_inversion_time']} s" + ) + print( + "================================================================================================================================================================" ) - # iterate over the pairs - for tracker, pair in enumerate(pairs): - # iterate over the quantization axes - for i, (Gij, Gji) in enumerate(zip(pair["Gij"], pair["Gji"])): - site_i = magnetic_entities[pair["ai"]] - site_j = magnetic_entities[pair["aj"]] - storage = [] - # iterate over the first order local perturbations in all possible orientations for the two sites - for Vui in site_i["Vu1"][i]: - for Vuj in site_j["Vu1"][i]: + if rank == root_node: + # iterate over the magnetic entities + for tracker, mag_ent in enumerate(magnetic_entities): + # iterate over the quantization axes + for i, Gii in enumerate(mag_ent["Gii"]): + storage = [] + # iterate over the first and second order local perturbations + for Vu1, Vu2 in zip(mag_ent["Vu1"][i], mag_ent["Vu2"][i]): # The Szunyogh-Lichtenstein formula traced = np.trace( - (Vui @ Gij @ Vuj @ Gji), axis1=1, axis2=2 + (Vu2 @ Gii + 0.5 * Gii @ Vu1 @ Gii), axis1=1, axis2=2 ) # this is the on site projection # evaluation of the contour integral storage.append(int_de_ke(traced, cont.we)) - # fill up the pairs dictionary with the energies - pairs[tracker]["energies"].append(storage) - # convert to np array - pairs[tracker]["energies"] = np.array(pairs[tracker]["energies"]) - - # calculate magnetic parameters - for mag_ent in magnetic_entities: - Kxx, Kyy, Kzz, consistency = calculate_anisotropy_tensor(mag_ent) - mag_ent["K"] = np.array([Kxx, Kyy, Kzz]) * sisl.unit_convert("eV", "meV") - mag_ent["K_consistency"] = consistency - - for pair in pairs: - J_iso, J_S, D, J = calculate_exchange_tensor(pair) - pair["J_iso"] = J_iso * sisl.unit_convert("eV", "meV") - pair["J_S"] = J_S * sisl.unit_convert("eV", "meV") - pair["D"] = D * sisl.unit_convert("eV", "meV") - pair["J"] = J * sisl.unit_convert("eV", "meV") - - times["end_time"] = timer() - print("\n\n\n\n\n") - print( - "##################################################################### GROGU OUTPUT #############################################################################" - ) + # fill up the magnetic entities dictionary with the energies + magnetic_entities[tracker]["energies"].append(storage) + # convert to np array + magnetic_entities[tracker]["energies"] = np.array( + magnetic_entities[tracker]["energies"] + ) + + # iterate over the pairs + for tracker, pair in enumerate(pairs): + # iterate over the quantization axes + for i, (Gij, Gji) in enumerate(zip(pair["Gij"], pair["Gji"])): + site_i = magnetic_entities[pair["ai"]] + site_j = magnetic_entities[pair["aj"]] + storage = [] + # iterate over the first order local perturbations in all possible orientations for the two sites + for Vui in site_i["Vu1"][i]: + for Vuj in site_j["Vu1"][i]: + # The Szunyogh-Lichtenstein formula + traced = np.trace( + (Vui @ Gij @ Vuj @ Gji), axis1=1, axis2=2 + ) # this is the on site projection + # evaluation of the contour integral + storage.append(int_de_ke(traced, cont.we)) + # fill up the pairs dictionary with the energies + pairs[tracker]["energies"].append(storage) + # convert to np array + pairs[tracker]["energies"] = np.array(pairs[tracker]["energies"]) + + # calculate magnetic parameters + for mag_ent in magnetic_entities: + Kxx, Kyy, Kzz, consistency = calculate_anisotropy_tensor(mag_ent) + mag_ent["K"] = np.array([Kxx, Kyy, Kzz]) * sisl.unit_convert("eV", "meV") + mag_ent["K_consistency"] = consistency - print_parameters(simulation_parameters) - print_atoms_and_pairs(magnetic_entities, pairs) - print( - "################################################################################################################################################################" - ) - print_runtime_information(times) - print("") - - # remove unwanted stuff before saving - pairs, magnetic_entities = remove_clutter_for_save(pairs, magnetic_entities) - # create output dictionary with all the relevant data - results = dict( - parameters=simulation_parameters, - magnetic_entities=magnetic_entities, - pairs=pairs, - runtime=times, - ) + for pair in pairs: + J_iso, J_S, D, J = calculate_exchange_tensor(pair) + pair["J_iso"] = J_iso * sisl.unit_convert("eV", "meV") + pair["J_S"] = J_S * sisl.unit_convert("eV", "meV") + pair["D"] = D * sisl.unit_convert("eV", "meV") + pair["J"] = J * sisl.unit_convert("eV", "meV") + + times["end_time"] = timer() + print("\n\n\n\n\n") + print( + "##################################################################### GROGU OUTPUT #############################################################################" + ) + + print_parameters(simulation_parameters) + print_atoms_and_pairs(magnetic_entities, pairs) + print( + "################################################################################################################################################################" + ) + print_runtime_information(times) + print("") + + # remove unwanted stuff before saving + pairs, magnetic_entities = remove_clutter_for_save(pairs, magnetic_entities) + # create output dictionary with all the relevant data + results = dict( + parameters=simulation_parameters, + magnetic_entities=magnetic_entities, + pairs=pairs, + runtime=times, + ) + + # save results + save_pickle(simulation_parameters["outfile"], results) + + print("\n\n\n\n\n") - # save results - save_pickle(simulation_parameters["outfile"], results) - print("\n\n\n\n\n") +if __name__ == "__main__": + main() diff --git a/src/grogupy/io.py b/src/grogupy/io.py index 834f232..bb5b474 100644 --- a/src/grogupy/io.py +++ b/src/grogupy/io.py @@ -35,13 +35,9 @@ default_args = dict( kset=2, kdirs="xyz", ebot=None, - automatic_ebot=False, eset=42, esetp=1000, - calculate_charge=False, - charges=[], parallel_solver_for_Gk=False, - parallel_size=None, padawan_mode=True, ) @@ -60,27 +56,33 @@ default_args = dict( def save_pickle(outfile, data): - """_summary_ + """Saves the data in the outfile with pickle. Args: - outfile (_type_): _description_ - data (_type_): _description_ + outfile : str + Path to outfile + data : dict + Contains the data """ + # save dictionary with open(outfile, "wb") as output_file: dump(data, output_file) -def load_pickle(infile, data): - """_summary_ +def load_pickle(infile): + """Loads the data from the infile with pickle. Args: - infile (_type_): _description_ - data (_type_): _description_ + infile : str + Path to infile Returns: - _type_: _description_ + data : dict + A dictionary of data """ + + # open and read file with open(infile, "wb") as input_file: data = load(data, input_file) @@ -88,11 +90,13 @@ def load_pickle(infile, data): def print_parameters(simulation_parameters): - """_summary_ + """It prints the simulation parameters for the grogu out. Args: - simulation_parameters (_type_): _description_ + simulation_parameters : dict + It contains the simulations parameters """ + print( "================================================================================================================================================================" ) @@ -138,12 +142,15 @@ def print_parameters(simulation_parameters): def print_atoms_and_pairs(magnetic_entities, pairs): - """_summary_ + """It prints the pair and magnetic entity information for the grogu out. Args: - magnetic_entities (_type_): _description_ - pairs (_type_): _description_ + magnetic_entities : dict + It contains the data on the magnetic entities + pairs : dict + It contains the data on the pairs """ + print("Atomic information: ") print( "----------------------------------------------------------------------------------------------------------------------------------------------------------------" @@ -223,11 +230,13 @@ def print_atoms_and_pairs(magnetic_entities, pairs): def print_runtime_information(times): - """_summary_ + """It prints the runtime information for the grogu out. Args: - times (_type_): _description_ + times : dict + It contains the runtime data """ + print("Runtime information: ") print(f"Total runtime: {times['end_time'] - times['start_time']} s") print( @@ -255,10 +264,12 @@ def print_runtime_information(times): def print_job_description(simulation_parameters): - """_summary_ + """It prints the parameters and the description of the job. + Args: - simulation_parameters (_type_): _description_ + simulation_parameters : dict + It contains the simulations parameters """ print( diff --git a/src/grogupy/magnetism.py b/src/grogupy/magnetism.py index e312ee4..90e01ce 100644 --- a/src/grogupy/magnetism.py +++ b/src/grogupy/magnetism.py @@ -21,21 +21,119 @@ import numpy as np +def blow_up_orbindx(orb_indices): + """Function to blow up orbital indices to make SPIN BOX indices. + + Args: + orb_indices : np.array_like + These are the indices in ORBITAL BOX + + Returns: + orb_indices : np.array_like + These are the indices in SPIN BOX + """ + + orb_indices = np.array([[2 * o, 2 * o + 1] for o in orb_indices]).flatten() + + return orb_indices + + +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. + + + Args: + M : np.array_like + Traceable 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] + 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 indexes of a given magnetic entity. + + Args: + dh : sisl.physics.Hamiltonian + Hamiltonian from sisl + 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: + 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 calculate_anisotropy_tensor(mag_ent): - """_summary_ + """Calculates the renormalized anisotropy tensor from the energies. + + It uses the grogu convention for output. Args: - mag_ent (_type_): _description_ + mag_ent : dict + An element from the magnetic entities Returns: - _type_: _description_ + K : np.array_like + elements of the anisotropy tensor """ + # get the energies energies = mag_ent["energies"] + + # calculate the diagonal tensor elements Kxx = energies[1, 1] - energies[1, 0] Kyy = energies[0, 1] - energies[0, 0] Kzz = 0 + # perform consistency check calculated_diff = Kyy - Kxx expected_diff = energies[2, 0] - energies[2, 1] consistency_check = abs(calculated_diff - expected_diff) @@ -44,14 +142,28 @@ def calculate_anisotropy_tensor(mag_ent): def calculate_exchange_tensor(pair): - """_summary_ + """Calculates the exchange tensor from the energies. + + It produces the isotropic exchange, the relevant elements + from the Dzyaloshinskii-Morilla (Dm) tensor, the symmetric-anisotropy + and the complete exchange tensor. Args: - pair (_type_): _description_ + pair : dict + An element from the pairs Returns: - _type_: _description_ + J_iso : float + Isotropic exchange (Tr[J] / 3) + J_S : np.array_like + Symmetric-anisotropy (J_S = J - J_iso * I ––> Jxx, Jyy, Jxy, Jxz, Jyz) + D : np.array_like + DM elements (Dx, Dy, Dz) + J : np.array_like + Complete exchange tensor flattened (Jxx, Jxy, Jxz, Jyx, Jyy, Jyz, Jzx, Jzy, Jzz) """ + + # energies from rotations energies = pair["energies"] # Initialize output arrays J = np.zeros((3, 3)) @@ -96,16 +208,3 @@ def calculate_exchange_tensor(pair): J_S = np.array([J[0, 0] - J_iso, J[1, 1] - J_iso, J[0, 1], J[0, 2], J[1, 2]]) return J_iso, J_S, D, J - - -def int_de_ke(traced, we): - """_summary_ - - Args: - traced (_type_): _description_ - we (_type_): _description_ - - Returns: - _type_: _description_ - """ - return np.trapz(-1 / np.pi * np.imag(traced * we)) diff --git a/src/grogupy/utils.py b/src/grogupy/utils.py index 2a124c8..0cb86f7 100644 --- a/src/grogupy/utils.py +++ b/src/grogupy/utils.py @@ -29,6 +29,25 @@ tau_z = np.array([[1, 0], [0, -1]]) tau_0 = np.array([[1, 0], [0, 1]]) +def commutator(a, b): + """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 + + # define some useful functions def hsk(H, ss, sc_off, k=(0, 0, 0)): """Speed up Hk and Sk generation. @@ -36,14 +55,20 @@ def hsk(H, ss, sc_off, k=(0, 0, 0)): 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). + 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 + np.array_like + Hamiltonian at the given k point + np.array_like + Overlap matrix at the given k point """ # this two conversion lines are from the sisl source @@ -60,46 +85,6 @@ def hsk(H, ss, sc_off, k=(0, 0, 0)): return HK, SK -def make_contour(emin=-20, emax=0.0, enum=42, p=150): - """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 - 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 - - # just an empty container class - class ccont: - 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 to sample the Brillouin zone. @@ -110,12 +95,17 @@ def make_kset(dirs="xyz", NUMK=20): 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. + 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. + np.array_like + An array of k points that uniformly sample the Brillouin zone in the given directions """ + + # if there is no xyz in dirs return the Gamma point if not (sum([d in dirs for d in "xyz"])): return np.array([[0, 0, 0]]) @@ -133,19 +123,50 @@ def make_kset(dirs="xyz", NUMK=20): return kset -def commutator(a, b): - """Shorthand for commutator. +def make_contour(emin=-20, emax=0.0, enum=42, p=150): + """A more sophisticated contour generator. - Commutator of two matrices in the mathematical sense. + Calculates the parameters for the complex contour integral. It uses the + Legendre-Gauss quadrature method. It returns a class that contains + the information for the contour integral. Args: - a (np.array_like): The first matrix. - b (np.array_like): The second matrix + 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: - np.array_like: The commutator of a and b. + ccont + Contains all the information for the contour integral """ - return a @ b - b @ a + x, wl = roots_legendre(enum) + R = (emax - emin) / 2 # radius + z0 = (emax + emin) / 2 # center point + y1 = -np.log(1 + np.pi * p) # lower bound + y2 = 0 # upper bound + + y = (y2 - y1) / 2 * x + (y2 + y1) / 2 + phi = (np.exp(-y) - 1) / p # angle parameter + ze = z0 + R * np.exp(1j * phi) # complex points for path + we = -(y2 - y1) / 2 * np.exp(-y) / p * 1j * (ze - z0) * wl + + # just an empty container class + class ccont: + pass + + cont = ccont() + cont.R = R + cont.z0 = z0 + cont.ze = ze + cont.we = we + cont.enum = enum + + return cont def tau_u(u): @@ -154,10 +175,12 @@ def tau_u(u): Returns the vector u in the basis of the Pauli matrices. Args: - u (list or np.array_like): The direction + u : list or np.array_like + The direction Returns: - np.array_like: Arbitrary direction in the base of the Pauli matrices. + np.array_like + Arbitrary direction in the base of the Pauli matrices """ # u is force to be of unit length @@ -173,10 +196,12 @@ def crossM(u): It acts as a cross product with vector u. Args: - u (list or np.array_like): The second vector in the cross product + 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. + 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]]) @@ -186,12 +211,16 @@ def RotM(theta, u, eps=1e-10): 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. + 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. + np.array_like + The rotation matrix """ u = u / np.linalg.norm(u) @@ -213,12 +242,16 @@ def RotMa2b(a, b, eps=1e-10): 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. + 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. + np.array_like + The rotation matrix with the above property """ v = np.cross(a, b) @@ -230,89 +263,42 @@ def RotMa2b(a, b, eps=1e-10): 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. +def read_siesta_emin(eigfile): + """It reads the lowest energy level from the siesta run. + It uses the .EIG file from siesta that contains the eigenvalues. Args: - M (np.array_like): Traceble matrix. + eigfile : str + The path to the .EIG file Returns: - dict: It contains the traced matrix with "x", "y", "z" and "c". + float + The energy minimum """ - 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 + # read the file + eigs = eigSileSiesta(eigfile).read_data() -def parse_magnetic_entity(dh, atom=None, l=None, **kwargs): - """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. + return eigs.min() - 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: - 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 int_de_ke(traced, we): + """It numerically integrates the traced matrix. -def read_siesta_emin(eigfile): - """_summary_ + It is a wrapper from numpy.trapz and it contains the + relevant constants to calculate the energy integral from + equation 93 or 96. Args: - eigfile (_type_): _description_ + traced : np.array_like + The trace of a matrix or a matrix product + we : float + The weight of a point on the contour Returns: - _type_: _description_ + float + The energy calculated from the integral formula """ - eigs = eigSileSiesta(eigfile).read_data() - return eigs.min() + return np.trapz(-1 / np.pi * np.imag(traced * we))