diff --git a/README.md b/README.md index 0df3eef..1c68879 100644 --- a/README.md +++ b/README.md @@ -18,7 +18,7 @@ More on the theoretical background can be seen on [arXiv](https://arxiv.org/abs/ - check on the validity of magnetic entities input - create a module instead of a script - add eset-kset-directions parallelization - +- there is a sign problem in J_symm ## Building wheel diff --git a/src/grogupy/__init__.py b/src/grogupy/__init__.py index 81e521d..37eca92 100644 --- a/src/grogupy/__init__.py +++ b/src/grogupy/__init__.py @@ -21,8 +21,35 @@ """Docstring in init. """ -from grogupy.constants import * -from grogupy.core import * -from grogupy.io import * -from grogupy.magnetism import * -from grogupy.utilities import * + +__all__ = [] + +__version__ = "0.0.0" + +import warnings + +from . import constants, core, io, magnetism, utilities + +# Pre-import most submodules. (We make exceptions for submodules that have +# special dependencies or where that would take too long.) + +__all__.extend(["constants", "core", "magnetism", "io", "utilities"]) + +# Make selected functionality available directly in the root namespace. +from .magnetism import MagneticEntity, Pair, Simulation + +__all__.extend(["Simulation", "MagneticEntity", "Pair"]) + +# Importing plotter might not work, but this does not have to be a problem -- +# only no plotting will be available. +try: + from . import plotter + from .plotter import plot +except: + pass +else: + __all__.extend(["plotter", "plot"]) +try: + from tqdm import tqdm +except: + print("Please install tqdm for nice progress bar.") diff --git a/src/grogupy/constants.py b/src/grogupy/constants.py index b813c4c..13e5b12 100644 --- a/src/grogupy/constants.py +++ b/src/grogupy/constants.py @@ -29,35 +29,34 @@ TAU_Z: Final[np.array] = np.array([[1, 0], [0, -1]]) TAU_0: Final[np.array] = np.array([[1, 0], [0, 1]]) -# list of accepted input parameters -ACCEPTED_INPUTS: Final[list] = [ - "infile", - "outfile", - "scf_xcf_orientation", - "ref_xcf_orientations", - "kset", - "kdirs", - "ebot", - "eset", - "esetp", - "parallel_solver_for_Gk", - "padawan_mode", +INFILE = "" +OUTFILE = "" +SCF_XCF_ORIENTATION = np.array([0, 0, 1]) +REF_XCF_ORIENTATIONS = [ + dict(o=np.array([1, 0, 0]), vw=[np.array([0, 1, 0]), np.array([0, 0, 1])]), + dict(o=np.array([0, 1, 0]), vw=[np.array([1, 0, 0]), np.array([0, 0, 1])]), + dict(o=np.array([0, 0, 1]), vw=[np.array([1, 0, 0]), np.array([0, 1, 0])]), ] -DEFAULT_ARGUMENTS: Final[dict] = dict( - infile=None, - outfile=None, - scf_xcf_orientation=np.array([0, 0, 1]), - ref_xcf_orientations=[ - dict(o=np.array([1, 0, 0]), vw=[np.array([0, 1, 0]), np.array([0, 0, 1])]), - dict(o=np.array([0, 1, 0]), vw=[np.array([1, 0, 0]), np.array([0, 0, 1])]), - dict(o=np.array([0, 0, 1]), vw=[np.array([1, 0, 0]), np.array([0, 1, 0])]), - ], - kset=2, - kdirs="xyz", - ebot=None, - eset=42, - esetp=1000, - parallel_solver_for_Gk=False, - padawan_mode=True, -) +KSET = 10 +KDIRS = "xy" +EBOT = -15 +ESET = 300 +ESETP = 1000 +PARALLEL_SOLVER_FOR_GK = False + +DEFAULT_PARAMETERS = dict() +DEFAULT_PARAMETERS["INFILE"] = INFILE +DEFAULT_PARAMETERS["OUTFILE"] = OUTFILE +DEFAULT_PARAMETERS["SCF_XCF_ORIENTATION"] = np.array([0, 0, 1]) +DEFAULT_PARAMETERS["REF_XCF_ORIENTATIONS"] = [ + dict(o=np.array([1, 0, 0]), vw=[np.array([0, 1, 0]), np.array([0, 0, 1])]), + dict(o=np.array([0, 1, 0]), vw=[np.array([1, 0, 0]), np.array([0, 0, 1])]), + dict(o=np.array([0, 0, 1]), vw=[np.array([1, 0, 0]), np.array([0, 1, 0])]), +] +DEFAULT_PARAMETERS["KSET"] = KSET +DEFAULT_PARAMETERS["KDIRS"] = KDIRS +DEFAULT_PARAMETERS["EBOT"] = EBOT +DEFAULT_PARAMETERS["ESET"] = ESET +DEFAULT_PARAMETERS["ESETP"] = ESETP +DEFAULT_PARAMETERS["PARALLEL_SOLVER_FOR_GK"] = PARALLEL_SOLVER_FOR_GK diff --git a/src/grogupy/core.py b/src/grogupy/core.py index 727ebbc..74b47d8 100644 --- a/src/grogupy/core.py +++ b/src/grogupy/core.py @@ -21,11 +21,11 @@ """Docstring in core. """ +import warnings + import numpy as np -from numpy.linalg import inv -from sisl.physics import Hamiltonian -from grogupy.magnetism import blow_up_orbindx, parse_magnetic_entity +from grogupy.constants import TAU_Y from grogupy.utilities import commutator @@ -70,287 +70,6 @@ def calc_Vu(H: np.array, Tu: np.array) -> np.array: return Vu1, Vu2 -def build_hh_ss(dh: Hamiltonian) -> tuple[np.array, np.array]: - """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 : sisl.physics.Hamiltonian - Hamiltonian read in by sisl - - Returns: - 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 - h11 = dh.tocsr(dh.M11r) - h11 += dh.tocsr(dh.M11i) * 1.0j - h11 = h11.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype("complex128") - - h22 = dh.tocsr(dh.M22r) - h22 += dh.tocsr(dh.M22i) * 1.0j - h22 = h22.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype("complex128") - - h12 = dh.tocsr(dh.M12r) - h12 += dh.tocsr(dh.M12i) * 1.0j - h12 = h12.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype("complex128") - - h21 = dh.tocsr(dh.M21r) - h21 += dh.tocsr(dh.M21i) * 1.0j - h21 = h21.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype("complex128") - - sov = ( - dh.tocsr(dh.S_idx) - .toarray() - .reshape(NO, dh.n_s, NO) - .transpose(0, 2, 1) - .astype("complex128") - ) - - # Reorganization of Hamiltonian and overlap matrix elements to SPIN BOX representation - U = np.vstack( - [np.kron(np.eye(NO, dtype=int), [1, 0]), np.kron(np.eye(NO, dtype=int), [0, 1])] - ) - # This is the permutation that transforms ud1ud2 to u12d12 - # That is this transforms FROM SPIN BOX to ORBITAL BOX => U - # the inverse transformation is U.T u12d12 to ud1ud2 - # That is FROM ORBITAL BOX to SPIN BOX => U.T - - # From now on everything is in SPIN BOX!! - hh, ss = np.array( - [ - U.T - @ np.block([[h11[:, :, i], h12[:, :, i]], [h21[:, :, i], h22[:, :, i]]]) - @ U - for i in range(dh.lattice.nsc.prod()) - ] - ), np.array( - [ - U.T - @ np.block( - [[sov[:, :, i], sov[:, :, i] * 0], [sov[:, :, i] * 0, sov[:, :, i]]] - ) - @ U - for i in range(dh.lattice.nsc.prod()) - ] - ) - - return hh, ss - - -def setup_pairs_and_magnetic_entities( - magnetic_entities: list, pairs: list, dh, simulation_parameters: dict -) -> tuple[list, list]: - """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: - 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: - 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 - mag_ent["orbital_indices"] = parsed - mag_ent["spin_box_indices"] = blow_up_orbindx( - parsed - ) # calculate spin box indexes - # 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"] = [] - # iterate over atoms - for atom_idx in mag_ent["atom"]: - mag_ent["tags"].append( - f"[{atom_idx}]{dh.atoms[atom_idx].tag}({mag_ent['l']})" - ) - mag_ent["xyz"].append(dh.xyz[atom_idx]) - - # calculate size for Greens function generation - spin_box_shape = len(mag_ent["spin_box_indices"]) - - # 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 - mag_ent["Vu2"] = [] - - mag_ent["Gii"] = [] # Greens function - mag_ent["Gii_tmp"] = [] # Greens function for parallelization - for _ in simulation_parameters["ref_xcf_orientations"]: - # Rotations for every quantization axis - mag_ent["Vu1"].append([]) - mag_ent["Vu2"].append([]) - # Greens functions for every quantization axis - mag_ent["Gii"].append( - np.zeros( - (simulation_parameters["eset"], spin_box_shape, spin_box_shape), - dtype="complex128", - ) - ) - mag_ent["Gii_tmp"].append( - np.zeros( - (simulation_parameters["eset"], spin_box_shape, spin_box_shape), - dtype="complex128", - ) - ) - - # for every site we have to store 2x3 Greens function (and the associated _tmp-s) - # in the 3 reference directions, because G_ij and G_ji are both needed - for pair in pairs: - # calculate distance - xyz_ai = magnetic_entities[pair["ai"]]["xyz"] - xyz_aj = magnetic_entities[pair["aj"]]["xyz"] - xyz_aj = xyz_aj + pair["Ruc"] @ simulation_parameters["cell"] - pair["dist"] = np.linalg.norm(xyz_ai - xyz_aj) - - # 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 = "" - # get atoms of magnetic entity - atoms_idx = mag_ent["atom"] - orbitals = mag_ent["l"] - - # if magnetic entity contains one atoms - if isinstance(atoms_idx, int): - tag += f"[{atoms_idx}]{dh.atoms[atoms_idx].tag}({orbitals})" - - # if magnetic entity contains more than one atoms - if isinstance(atoms_idx, list): - # iterate over atoms - atom_group = "{" - for atom_idx in atoms_idx: - atom_group += f"[{atom_idx}]{dh.atoms[atom_idx].tag}({orbitals})--" - # end {} of the atoms in the magnetic entity - tag += atom_group[:-2] + "}" - pair["tags"].append(tag) - pair["energies"] = [] # we will store the second order energy derivations here - - pair["Gij"] = [] # Greens function - pair["Gji"] = [] - pair["Gij_tmp"] = [] # Greens function for parallelization - pair["Gji_tmp"] = [] - for _ in simulation_parameters["ref_xcf_orientations"]: - # Greens functions for every quantization axis - pair["Gij"].append( - np.zeros( - (simulation_parameters["eset"], spin_box_shape_i, spin_box_shape_j), - dtype="complex128", - ) - ) - pair["Gij_tmp"].append( - np.zeros( - (simulation_parameters["eset"], spin_box_shape_i, spin_box_shape_j), - dtype="complex128", - ) - ) - pair["Gji"].append( - np.zeros( - (simulation_parameters["eset"], spin_box_shape_j, spin_box_shape_i), - dtype="complex128", - ) - ) - pair["Gji_tmp"].append( - np.zeros( - (simulation_parameters["eset"], spin_box_shape_j, spin_box_shape_i), - dtype="complex128", - ) - ) - - return pairs, magnetic_entities - - -def parallel_Gk(HK: np.array, SK: np.array, eran: np.array, eset: int) -> np.array: - """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: np.array, SK: np.array, eran: np.array, eset: int) -> np.array: - """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) - - return Gk - - def remove_clutter_for_save(pairs: list, magnetic_entities: list) -> tuple[list, list]: """Removes unimportant data from the dictionaries. diff --git a/src/grogupy/grogu.py b/src/grogupy/grogu.py index ac1df4f..b358655 100644 --- a/src/grogupy/grogu.py +++ b/src/grogupy/grogu.py @@ -21,9 +21,6 @@ """Docstring in grogupy. """ -import warnings -from argparse import ArgumentParser -from sys import getsizeof from timeit import default_timer as timer # use numpy number of threads one @@ -39,260 +36,13 @@ import numpy as np import sisl from mpi4py import MPI -try: - from tqdm import tqdm - - tqdm_imported: bool = True -except: - print("Please install tqdm for nice progress bar.") - tqdm_imported: bool = False - - from grogupy import * -def parse_command_line() -> dict: - """This function can read input from the command line.""" - - parser = ArgumentParser() - - parser.add_argument( - "-i", - "--input", - dest="grogupy_infile", - default=None, - type=str, - help="Input file name for the Grogupy fdf file", - required=True, - ) - parser.add_argument( - "--siesta-input", - dest="infile", - default=None, - type=str, - help="Input file name for the Siesta fdf file", - ) - parser.add_argument( - "-o", - "--output", - dest="outfile", - default=None, - type=str, - help="Output file name", - ) - parser.add_argument( - "--kset", - dest="kset", - default=None, - type=int, - help="k-space resolution of calculation", - ) - parser.add_argument( - "--kdirs", - dest="kdirs", - default=None, - type=str, - help="Definition of k-space dimensionality", - ) - parser.add_argument( - "--ebot", - dest="ebot", - default=None, - type=float, - help="Bottom energy of the contour", - ) - parser.add_argument( - "--eset", - dest="eset", - default=None, - type=int, - help="Number of energy points on the contour", - ) - parser.add_argument( - "--eset-p", - dest="esetp", - default=None, - type=int, - help="Parameter tuning the distribution on the contour", - ) - parser.add_argument( - "--parallel-green", - dest="parallel_solver_for_Gk", - default=None, - type=bool, - help="Whether to use the parallel or sequential solver for Greens function", - ) - parser.add_argument( - "--padawan-mode", - dest="padawan_mode", - default=None, - type=bool, - help="If it is on it turns on extra helpful information for new users", - ) - - # convert to dictionary - cmd_line_args = parser.parse_args() - cmd_line_args = vars(cmd_line_args) - - return cmd_line_args - - def main(simulation_parameters: dict, magnetic_entities: list, pairs: list) -> None: - # runtime information - times: dict = dict() - times["start_time"] = timer() - - # MPI parameters - comm = MPI.COMM_WORLD - size: int = comm.Get_size() - rank: int = comm.Get_rank() - root_node: int = 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 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 - - # add third direction for off-diagonal anisotropy - for orientation in simulation_parameters["ref_xcf_orientations"]: - o1: np.array = orientation["vw"][0] - o2: np.array = orientation["vw"][1] - orientation["vw"] = np.vstack((orientation["vw"], (o1 + o2) / np.sqrt(2))) - - # 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: int = 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") - - times["setup_time"] = timer() - print(f"Setup done. Elapsed time: {times['setup_time']} s") - print( - "================================================================================================================================================================" - ) - - NO: Final[int] = dh.no # shorthand for number of orbitals in the unit cell - - # reformat Hamiltonian and Overlap matrix for manipulations - hh, ss = build_hh_ss(dh) - - # copy arrays for tests - if rank == root_node: - hh_test, ss_test = hh.copy(), ss.copy() - - # 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 - - # test symmetrization - if rank == root_node: - diff_hh = abs(hh - hh_test).max() - diff_ss = abs(ss - ss_test).max() - if diff_hh > 1e-12: - warnings.warn( - f"Hamiltonian changed after symmetrization. Largest change is {diff_hh}" - ) - if diff_ss > 1e-12: - warnings.warn( - f"Overlap matrix changed after symmetrization. Largest change is {diff_hh}" - ) - - # identifying TRS and TRB parts of the Hamiltonian - TAUY: np.array = np.kron(np.eye(NO), TAU_Y) - hTR: np.array = np.array( - [TAUY @ hh[i].conj() @ TAUY for i in range(dh.lattice.nsc.prod())] - ) - hTRS: np.array = (hh + hTR) / 2 - hTRB: np.array = (hh - hTR) / 2 - - # extracting the exchange field - traced: np.array = [ - spin_tracer(hTRB[i]) for i in range(dh.lattice.nsc.prod()) - ] # equation 77 - XCF: np.array = 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]), - ] - ) - - # check if exchange field has scalar part - max_xcfs: float = 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 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( - "================================================================================================================================================================" - ) - - # initialize pairs and magnetic entities based on input information - pairs, magnetic_entities = setup_pairs_and_magnetic_entities( - magnetic_entities, pairs, dh, simulation_parameters - ) - - 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( - "================================================================================================================================================================" - ) - # generate k space sampling kset = make_kset( dirs=simulation_parameters["kdirs"], NUMK=simulation_parameters["kset"] @@ -304,92 +54,11 @@ def main(simulation_parameters: dict, magnetic_entities: list, pairs: list) -> N # split the k points based on MPI size kpcs: np.array = 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: list = [] - - # iterate over the reference directions (quantization axes) - for i, orient in enumerate(simulation_parameters["ref_xcf_orientations"]): - # obtain rotated exchange field and Hamiltonian - R: np.array = RotMa2b(simulation_parameters["scf_xcf_orientation"], orient["o"]) - rot_XCF: np.array = np.einsum("ij,jklm->iklm", R, XCF) - rot_H_XCF: np.array = sum( - [np.kron(rot_XCF[i], tau) for i, tau in enumerate([TAU_X, TAU_Y, TAU_Z])] - ) - rot_H_XCF_uc: np.array = rot_H_XCF[uc_in_sc_idx] - - # obtain total Hamiltonian with the rotated exchange field - rot_H: np.array = 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.array = 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: int = 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( - "================================================================================================================================================================" - ) - - # 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 {2*NO}x{2*NO}={4*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: int = 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 * 32 * len(hamiltonians) * simulation_parameters['eset']} 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( - "================================================================================================================================================================" - ) + hamiltonians: list = [] # MPI barrier comm.Barrier() @@ -445,59 +114,6 @@ def main(simulation_parameters: dict, magnetic_entities: list, pairs: list) -> N 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( - "================================================================================================================================================================" - ) - - 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: list = [] - # 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.array = 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"] - ) - - # 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: int = magnetic_entities[pair["ai"]] - site_j: int = magnetic_entities[pair["aj"]] - storage: list = [] - # iterate over the first order local perturbations in all possible orientations for the two sites - # actually all possible orientations without the orientation for the off-diagonal anisotropy - # that is why we only take the first two of each Vu1 - for Vui in site_i["Vu1"][i][:2]: - for Vuj in site_j["Vu1"][i][:2]: - # The Szunyogh-Lichtenstein formula - traced: np.array = 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: K, consistency = calculate_anisotropy_tensor(mag_ent) @@ -539,20 +155,3 @@ def main(simulation_parameters: dict, magnetic_entities: list, pairs: list) -> N save_pickle(simulation_parameters["outfile"], results) print("\n\n\n\n\n") - - -if __name__ == "__main__": - # loading parameters - command_line_arguments = parse_command_line() - command_line_arguments = dict(command_line_arguments) - fdf_arguments, magnetic_entities, pairs = read_grogupy_fdf( - command_line_arguments["grogupy_infile"] - ) - - # combine the inputs to a single dictionary - simulation_parameters = process_input_args( - DEFAULT_ARGUMENTS, fdf_arguments, command_line_arguments, ACCEPTED_INPUTS - ) - - # run grogupy - main(simulation_parameters, magnetic_entities, pairs) diff --git a/src/grogupy/io.py b/src/grogupy/io.py index 282f95a..7d2cdb5 100644 --- a/src/grogupy/io.py +++ b/src/grogupy/io.py @@ -23,13 +23,10 @@ from os import environ from pickle import dump, load -from typing import Final import numpy as np from sisl.io import fdfSileSiesta -from grogupy.constants import ACCEPTED_INPUTS - def read_grogupy_fdf(path: str) -> tuple[dict, list, list]: """It reads the simulation parameters, magnetic entities and pairs from the fdf. @@ -168,7 +165,7 @@ def process_input_args( DEFAULT_ARGUMENTS: dict, fdf_arguments: dict, command_line_arguments: dict, - accepted_inputs: dict = ACCEPTED_INPUTS, + accepted_inputs: dict, ) -> dict: """It returns the final simulation parameters based on the inputs. diff --git a/src/grogupy/magnetism.py b/src/grogupy/magnetism.py index 5067d30..9c28206 100644 --- a/src/grogupy/magnetism.py +++ b/src/grogupy/magnetism.py @@ -21,234 +21,1031 @@ """Docstring in magnetism. """ +import copy +import warnings +from os import environ from typing import Union import numpy as np -from sisl.physics import Hamiltonian - - -def blow_up_orbindx(orb_indices: np.array) -> np.array: - """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: np.array) -> dict: - """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: Hamiltonian, - atom: Union[None, int, list] = None, - l: Union[None, int, list] = None, - orb: Union[None, int, list] = None, -) -> np.array: - """Function to define orbital indexes of a given magnetic entity. - - There are five possible input types: - 1. Cluster: a list of atoms - 2. Atom: must be only one atom - 3. AtomShell: one atom and a list of shells indexed in the atom - 4. AtomOrbital: one atom and a list of orbitals indexed in the atom - 5. Orbitals: a list of orbitals indexed in the Hamiltonian - Everything else raises an Exception. - - 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 or list of integers, optional - Defining the angular momentum channel. Defaults to None - orb : integer or list of integers, optional - Defining the orbital index in the Hamiltonian or on the atom. Defaults to None - - Returns: - list - The orbital indexes of the given magnetic entity - """ - - # the case for the Orbitals keyword - if atom is None and l is None: - orbital_indexes = orb - - # the case for the Cluster keyword - elif isinstance(atom, list) and l is None and orb is None: - orbital_indexes = [] - for a in atom: - a_orb_idx = dh.geometry.a2o(a, all=True) - orbital_indexes.append(a_orb_idx) - orbital_indexes = np.array(orbital_indexes).flatten() - - # the case for the Atom keyword - elif isinstance(atom, int) and l is None and orb is None: - orbital_indexes = dh.geometry.a2o(atom, all=True) - - # the case for the Atomshell keyword - elif isinstance(atom, int) and l is not None and orb is None: - orbital_indexes = dh.geometry.a2o(atom, all=True) - # make sure l is a list for the following step - if isinstance(l, int): - l = [l] - mask = [orbital.l in l for orbital in dh.geometry.atoms[atom].orbitals] - orbital_indexes = orbital_indexes[mask] - - # the case for the AtomOrbital keyword - elif isinstance(atom, int) and l is None and orb is not None: - orbital_indexes = dh.geometry.a2o(atom, all=True) - # make sure orb is a list for the following step - if isinstance(orb, int): - orb = [orb] - orbital_indexes = orbital_indexes[orb] - - return orbital_indexes # numpy array containing integers labeling orbitals associated to a magnetic entity. - - -def calculate_anisotropy_tensor(mag_ent: dict) -> tuple[np.array, float]: - """Calculates the renormalized anisotropy tensor from the energies. - - Args: - mag_ent : dict - An element from the magnetic entities - - Returns: - K : np.array_like - elements of the anisotropy tensor - consistency_check : float - Absolute value of the difference from the consistency check - """ - - # get the energies - energies = mag_ent["energies"] - - K = np.zeros((3, 3)) - - # calculate the diagonal tensor elements - K[0, 0] = energies[1, 1] - energies[1, 0] - K[1, 1] = energies[0, 1] - energies[0, 0] - K[2, 2] = 0 - - # calculate the off-diagonal tensor elements - K[0, 1] = (energies[2, 0] + energies[2, 1]) / 2 - energies[2, 2] - K[1, 0] = K[0, 1] - K[0, 2] = (energies[1, 0] + energies[1, 1]) / 2 - energies[1, 2] - K[2, 0] = K[0, 2] - K[1, 2] = (energies[0, 0] + energies[0, 1]) / 2 - energies[0, 2] - K[2, 1] = K[1, 2] - - # perform consistency check - calculated_diff = K[1, 1] - K[0, 0] - expected_diff = energies[2, 0] - energies[2, 1] - consistency_check = abs(calculated_diff - expected_diff) - - return K, consistency_check - - -def calculate_exchange_tensor(pair: dict) -> tuple[float, np.array, np.array, np.array]: - """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 : dict - An element from the pairs - - Returns: - 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)) - D = np.zeros(3) - - # J matrix calculations - # J(1,1) = mean([DEij(2,2,2), DEij(2,2,3)]) - J[0, 0] = np.mean([energies[1, 3], energies[2, 3]]) - - # J(1,2) = -mean([DEij(1,2,3), DEij(2,1,3)]) - J[0, 1] = -np.mean([energies[2, 1], energies[2, 2]]) - J[1, 0] = J[0, 1] - - # J(1,3) = -mean([DEij(1,2,2), DEij(2,1,2)]) - J[0, 2] = -np.mean([energies[1, 1], energies[1, 2]]) - J[2, 0] = J[0, 2] - - # J(2,2) = mean([DEij(2,2,1), DEij(1,1,3)]) - J[1, 1] = np.mean([energies[0, 3], energies[2, 0]]) - - # J(2,3) = -mean([DEij(1,2,1), DEij(2,1,1)]) - J[1, 2] = -np.mean([energies[0, 1], energies[0, 2]]) - J[2, 1] = J[1, 2] - - # J(3,3) = mean([DEij(1,1,1), DEij(1,1,2)]) - J[2, 2] = np.mean([energies[0, 0], energies[1, 0]]) - - # D vector calculations - # D(1) = mean([DEij(1,2,1), -DEij(2,1,1)]) - D[0] = np.mean([energies[0, 1], -energies[0, 2]]) +import sisl +from mpi4py import MPI +from scipy.special import roots_legendre +from tqdm import tqdm as tqdm - # D(2) = mean([DEij(2,1,2), -DEij(1,2,2)]) - D[1] = np.mean([energies[1, 2], -energies[1, 1]]) +import grogupy.constants as constants +from grogupy.core import calc_Vu, onsite_projection +from grogupy.utilities import RotMa2b, spin_tracer, tau_u - # D(3) = mean([DEij(1,2,3), -DEij(2,1,3)]) - D[2] = np.mean([energies[2, 1], -energies[2, 2]]) - - J_iso = np.trace(J) / 3 - # based on the grogu output pdf - # traceless symmetric exchange matrix: - # Jxx, Jyy, Jxy, Jxz, Jyz - 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 + +class Hamiltonian: + def __init__(self) -> None: + self.H = None + self.S = None + self.NO = None + self.cell = None + self.nsc = None + self.sc_off = None + self.unit_cell_index = None + self.G = None + self.orientation = None + self.hTRS = None + self.hTRB = None + self.XCF = None + self.HK = None + self.SK = None + + def copy(self): + return copy.deepcopy(self) + + def rot_H_and_XCF(self, scf_xcf, orientation): + # obtain rotated exchange field and Hamiltonian + R: np.array = RotMa2b(scf_xcf, orientation) + rot_XCF: np.array = np.einsum("ij,jklm->iklm", R, self.XCF) + rot_H_XCF: np.array = sum( + [ + np.kron(rot_XCF[i], tau) + for i, tau in enumerate( + [constants.TAU_X, constants.TAU_Y, constants.TAU_Z] + ) + ] + ) + + rot_H_XCF_uc: np.array = rot_H_XCF[self.unit_cell_index] + + # obtain total Hamiltonian with the rotated exchange field + rot_H: np.array = self.hTRS + rot_H_XCF # equation 76 + + return rot_H, rot_H_XCF_uc + + def build_custom_H_and_S( + self, dh: sisl.physics.Hamiltonian + ) -> tuple[np.array, np.array]: + """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 : sisl.physics.Hamiltonian + Hamiltonian read in by sisl + + Returns: + 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 + """ + + self.NO = dh.no + self.cell = dh.cell + self.nsc = dh.lattice.nsc + self.sc_off = dh.sc_off + self.unit_cell_index = dh.lattice.sc_index([0, 0, 0]) + self.sislHamiltonian = dh + + NO = dh.no # shorthand for number of orbitals in the unit cell + + # preprocessing Hamiltonian and overlap matrix elements + h11 = dh.tocsr(dh.M11r) + h11 += dh.tocsr(dh.M11i) * 1.0j + h11 = ( + h11.toarray() + .reshape(NO, dh.n_s, NO) + .transpose(0, 2, 1) + .astype("complex128") + ) + + h22 = dh.tocsr(dh.M22r) + h22 += dh.tocsr(dh.M22i) * 1.0j + h22 = ( + h22.toarray() + .reshape(NO, dh.n_s, NO) + .transpose(0, 2, 1) + .astype("complex128") + ) + + h12 = dh.tocsr(dh.M12r) + h12 += dh.tocsr(dh.M12i) * 1.0j + h12 = ( + h12.toarray() + .reshape(NO, dh.n_s, NO) + .transpose(0, 2, 1) + .astype("complex128") + ) + + h21 = dh.tocsr(dh.M21r) + h21 += dh.tocsr(dh.M21i) * 1.0j + h21 = ( + h21.toarray() + .reshape(NO, dh.n_s, NO) + .transpose(0, 2, 1) + .astype("complex128") + ) + + sov = ( + dh.tocsr(dh.S_idx) + .toarray() + .reshape(NO, dh.n_s, NO) + .transpose(0, 2, 1) + .astype("complex128") + ) + + # Reorganization of Hamiltonian and overlap matrix elements to SPIN BOX representation + U = np.vstack( + [ + np.kron(np.eye(NO, dtype=int), [1, 0]), + np.kron(np.eye(NO, dtype=int), [0, 1]), + ] + ) + # This is the permutation that transforms ud1ud2 to u12d12 + # That is this transforms FROM SPIN BOX to ORBITAL BOX => U + # the inverse transformation is U.T u12d12 to ud1ud2 + # That is FROM ORBITAL BOX to SPIN BOX => U.T + + # From now on everything is in SPIN BOX!! + hh, ss = np.array( + [ + U.T + @ np.block([[h11[:, :, i], h12[:, :, i]], [h21[:, :, i], h22[:, :, i]]]) + @ U + for i in range(dh.lattice.nsc.prod()) + ] + ), np.array( + [ + U.T + @ np.block( + [[sov[:, :, i], sov[:, :, i] * 0], [sov[:, :, i] * 0, sov[:, :, i]]] + ) + @ U + for i in range(dh.lattice.nsc.prod()) + ] + ) + + # copy arrays for tests + hh_test, ss_test = hh.copy(), ss.copy() + + # 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 + + # test symmetrization + diff_hh = abs(hh - hh_test).max() + diff_ss = abs(ss - ss_test).max() + if diff_hh > 1e-12: + warnings.warn( + f"Hamiltonian changed after symmetrization. Largest change is {diff_hh}" + ) + if diff_ss > 1e-12: + warnings.warn( + f"Overlap matrix changed after symmetrization. Largest change is {diff_hh}" + ) + + self.H = hh + self.S = ss + self.hTRS, self.hTRB, self.XCF = self._extract_exchange_field() + + def _extract_exchange_field(self) -> np.array: + """Extracts the exchange field from the Hamiltonian. + + Args: + hh : (NO, NO) np.array_like + Hamiltonian in spin box representation + Returns: + np.array_like: + XCF field + """ + + # identifying TRS and TRB parts of the Hamiltonian + TAUY: np.array = np.kron(np.eye(self.NO), constants.TAU_Y) + hTR: np.array = np.array( + [TAUY @ self.H[i].conj() @ TAUY for i in range(self.nsc.prod())] + ) + hTRS: np.array = (self.H + hTR) / 2 + hTRB: np.array = (self.H - hTR) / 2 + + # extracting the exchange field + traced: np.array = [ + spin_tracer(hTRB[i]) for i in range(self.nsc.prod()) + ] # equation 77 + XCF: np.array = 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]), + ] + ) + + # check if exchange field has scalar part + max_xcfs: float = 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}" + ) + + return hTRS, hTRB, XCF + + def _parallel_Gk(self, eran: np.array, eset: int) -> np.array: + """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 np.linalg.inv(self.SK * eran.reshape(eset, 1, 1) - self.HK) + + def _sequential_Gk(self, eran: np.array, eset: int) -> np.array: + """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((eset, self.HK.shape[0], self.HK.shape[1]), dtype="complex128") + # fills the holder sequentially by the Greens function on a given energy + for j in range(eset): + Gk[j] = np.linalg.inv(self.SK * eran[j] - self.HK) + + return Gk + + def setup_at_k_point(self, k: tuple = (0, 0, 0)) -> tuple[np.array, np.array]: + """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 indices 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 + """ + + # 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 @ self.sc_off.T) + + # phases applied to the hamiltonian + self.HK = np.einsum("abc,a->bc", self.H, phases) + self.SK = np.einsum("abc,a->bc", self.S, phases) + + +class Contour: + def __init__(self, emin: float, eset: int, esetp: float, emax: float = 0) -> None: + self.emin = emin + self.emax = emax + self.eset = eset + self.esetp = esetp + samples, weights = self._make_contour(emin, emax, eset, esetp) + self.samples = samples + self.weights = weights + + def _make_contour( + self, emin: float, emax: float, enum: int, p: float + ) -> tuple[np.array, float]: + """A more sophisticated contour generator. + + 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: + 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: + Container + Contains all the information for the contour integral + """ + 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 + + return ze, we + + +class KSpace: + def __init__(self, kdirs: str, kset: int, parallel_size: int = 1) -> None: + self.kdirs = kdirs + self.kset = kset + self.kpoints = self._generate_k_points(self.kdirs, self.kset) + self.weights = np.ones(len(self.kpoints)) / len(self.kpoints) + self.__split_kpoints(parallel_size) + + def _generate_k_points(self, kdirs: str, kset: int) -> np.array: + """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 returned. The total number of k points is the NUMK**(dimensions) + + + Args: + kdirs : str + Directions of the k points in the Brillouin zone. They are the three lattice vectors + kset : int + The number of k points in a direction + + Returns: + 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 kdirs for d in "xyz"])): + return np.array([[0, 0, 0]]) + + kran = len(kdirs) * [np.linspace(0, 1, kset, endpoint=False)] + mg = np.meshgrid(*kran) + directions = dict() + + for d in enumerate(kdirs): + directions[d[1]] = mg[d[0]].flatten() + for d in "xyz": + if not (d in kdirs): + directions[d] = 0 * directions[kdirs[0]] + kset = np.array([directions[d] for d in "xyz"]).T + + return kset + + def __split_kpoints(self, size): + self.parallel_weights = np.array_split(self.weights, size) + self.parallel_kpoints = np.array_split(self.kpoints, size) + self.parallel_kpoints[0] = tqdm(self.parallel_kpoints[0]) + + +class MagneticEntity: + number_of_entities = 0 + + def __init__( + self, + dh: sisl.physics.Hamiltonian, + atom: Union[None, int, list] = None, + l: Union[None, int, list] = None, + orb: Union[None, int, list] = None, + ) -> None: + self.orbital_indices = self.__parse_for_orbital_indices(dh, atom, l, orb) + self.atom = self.__parse_for_atom_indices(dh, self.orbital_indices) + self.spin_box_indices = self.__blow_up_orbital_index(self.orbital_indices) + self.SBS = len(self.spin_box_indices) + self.xyz = [dh.xyz[i] for i in self.atom] + try: + self.tag = [f"[{i}]{dh.atoms[i].tag}({self.l[i]})" for i in self.atom] + except: + self.tag = [f"[{i}]{dh.atoms[i].tag}(Unknown)" for i in self.atom] + self.Vu1 = [] + self.Vu2 = [] + self.Gii = [] + self._Gii_tmp = [] + self._energies = [] + self.K = [] + self.K_consistency = None + + MagneticEntity.number_of_entities += 1 + + def __parse_for_atom_indices(self, dh, orbital_indices): + atoms = [] + for orb in orbital_indices: + atoms.append(dh.o2a(orb)) + atoms = np.sort(np.unique(np.array(atoms))) + return atoms + + def __blow_up_orbital_index(self, orb_indices: np.array) -> np.array: + """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 __parse_for_orbital_indices( + self, + dh: Hamiltonian, + atom: Union[None, int, list] = None, + l: Union[None, int, list] = None, + orb: Union[None, int, list] = None, + ) -> np.array: + """Function to define orbital indices of a given magnetic entity. + + There are five possible input types: + 1. Cluster: a list of atoms + 2. Atom: must be only one atom + 3. AtomShell: one atom and a list of shells indexed in the atom + 4. AtomOrbital: one atom and a list of orbitals indexed in the atom + 5. Orbitals: a list of orbitals indexed in the Hamiltonian + Everything else raises an Exception. + + 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 or list of integers, optional + Defining the angular momentum channel. Defaults to None + orb : integer or list of integers, optional + Defining the orbital index in the Hamiltonian or on the atom. Defaults to None + + Returns: + list + The orbital indices of the given magnetic entity + """ + + # the case for the Orbitals keyword + if atom is None and l is None: + orbital_indices = orb + + # the case for the Cluster keyword + elif isinstance(atom, list) and l is None and orb is None: + orbital_indices = [] + for a in atom: + a_orb_idx = dh.geometry.a2o(a, all=True) + orbital_indices.append(a_orb_idx) + orbital_indices = np.array(orbital_indices).flatten() + + # the case for the Atom keyword + elif isinstance(atom, int) and l is None and orb is None: + orbital_indices = dh.geometry.a2o(atom, all=True) + + # the case for the Atomshell keyword + elif isinstance(atom, int) and l is not None and orb is None: + orbital_indices = dh.geometry.a2o(atom, all=True) + # make sure l is a list for the following step + if isinstance(l, int): + l = [l] + mask = [orbital.l in l for orbital in dh.geometry.atoms[atom].orbitals] + orbital_indices = orbital_indices[mask] + + # the case for the AtomOrbital keyword + elif isinstance(atom, int) and l is None and orb is not None: + orbital_indices = dh.geometry.a2o(atom, all=True) + # make sure orb is a list for the following step + if isinstance(orb, int): + orb = [orb] + orbital_indices = orbital_indices[orb] + + return orbital_indices # numpy array containing integers labeling orbitals associated to a magnetic entity. + + def calculate_energies(self, weights): + energies = [] + # iterate over the quantization axes + for i, Gii in enumerate(self.Gii): + storage: list = [] + # iterate over the first and second order local perturbations + for Vu1, Vu2 in zip(self.Vu1[i], self.Vu2[i]): + # The Szunyogh-Lichtenstein formula + traced: np.array = 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(np.trapz(-1 / np.pi * np.imag(traced * weights))) + # fill up the magnetic entities dictionary with the energies + energies.append(storage) + # convert to np array + self._energies = np.array(energies) + + def calculate_anisotropy(self): + """Calculates the renormalized anisotropy tensor from the energies. + + Args: + mag_ent : dict + An element from the magnetic entities + + Returns: + K : np.array_like + elements of the anisotropy tensor + consistency_check : float + Absolute value of the difference from the consistency check + """ + + # get the energies + energies = self._energies + + K = np.zeros((3, 3)) + + # calculate the diagonal tensor elements + K[0, 0] = energies[1, 1] - energies[1, 0] + K[1, 1] = energies[0, 1] - energies[0, 0] + K[2, 2] = 0 + + # calculate the off-diagonal tensor elements + K[0, 1] = (energies[2, 0] + energies[2, 1]) / 2 - energies[2, 2] + K[1, 0] = K[0, 1] + K[0, 2] = (energies[1, 0] + energies[1, 1]) / 2 - energies[1, 2] + K[2, 0] = K[0, 2] + K[1, 2] = (energies[0, 0] + energies[0, 1]) / 2 - energies[0, 2] + K[2, 1] = K[1, 2] + + # perform consistency check + calculated_diff = K[1, 1] - K[0, 0] + expected_diff = energies[2, 0] - energies[2, 1] + consistency_check = abs(calculated_diff - expected_diff) + + self.K = K * sisl.unit_convert("eV", "meV") + self.K_consistency = consistency_check * sisl.unit_convert("eV", "meV") + + def dictionary(self): + out = dict() + warnings.warn("Tag creation by orb index") + + out["orbital_indices"] = self.orbital_indices + out["spin_box_indices"] = self.spin_box_indices + out["xyz"] = self.xyz + out["tag"] = self.tag + out["energies"] = self.energies() + + return out + + +class Pair: + number_of_pairs = 0 + + def __init__( + self, + M1: MagneticEntity, + M2: MagneticEntity, + hamiltonian: Hamiltonian, + supercell_shift: Union[list, np.array] = np.array([0, 0, 0]), + ) -> None: + self.M1 = M1 + self.M2 = M2 + self.supercell_shift = np.array(supercell_shift) + self.supercell_shift_xyz = self.supercell_shift @ hamiltonian.cell + self.xyz = np.array([M1.xyz, M2.xyz + self.supercell_shift_xyz]) + self.distance = np.linalg.norm(self.xyz[0] - self.xyz[1]) + self.SBS1 = M1.SBS + self.SBS2 = M2.SBS + self.tags = np.array([M1.tag, M2.tag]) + self.Gij = [] + self.Gji = [] + self._Gij_tmp = [] + self._Gji_tmp = [] + self._energies = [] + + Pair.number_of_pairs += 1 + + def calculate_energies(self, weights): + energies = [] + # iterate over the quantization axes + for i, (Gij, Gji) in enumerate(zip(self.Gij, self.Gji)): + storage: list = [] + # iterate over the first order local perturbations in all possible orientations for the two sites + # actually all possible orientations without the orientation for the off-diagonal anisotropy + # that is why we only take the first two of each Vu1 + for Vui in self.M1.Vu1[i][:2]: + for Vuj in self.M1.Vu2[i][:2]: + # The Szunyogh-Lichtenstein formula + traced: np.array = np.trace( + (Vui @ Gij @ Vuj @ Gji), axis1=1, axis2=2 + ) # this is the on site projection + # evaluation of the contour integral + storage.append(np.trapz(-1 / np.pi * np.imag(traced * weights))) + # fill up the pairs dictionary with the energies + energies.append(storage) + # convert to np array + self._energies = np.array(energies) + + def calculate_exchange_tensor(self) -> tuple[float, np.array, np.array, np.array]: + """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 : dict + An element from the pairs + + Returns: + 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 = self._energies + # Initialize output arrays + J = np.zeros((3, 3)) + D = np.zeros(3) + + # J matrix calculations + # J(1,1) = mean([DEij(2,2,2), DEij(2,2,3)]) + J[0, 0] = np.mean([energies[1, 3], energies[2, 3]]) + + # J(1,2) = -mean([DEij(1,2,3), DEij(2,1,3)]) + J[0, 1] = -np.mean([energies[2, 1], energies[2, 2]]) + J[1, 0] = J[0, 1] + + # J(1,3) = -mean([DEij(1,2,2), DEij(2,1,2)]) + J[0, 2] = -np.mean([energies[1, 1], energies[1, 2]]) + J[2, 0] = J[0, 2] + + # J(2,2) = mean([DEij(2,2,1), DEij(1,1,3)]) + J[1, 1] = np.mean([energies[0, 3], energies[2, 0]]) + + # J(2,3) = -mean([DEij(1,2,1), DEij(2,1,1)]) + J[1, 2] = -np.mean([energies[0, 1], energies[0, 2]]) + J[2, 1] = J[1, 2] + + # J(3,3) = mean([DEij(1,1,1), DEij(1,1,2)]) + J[2, 2] = np.mean([energies[0, 0], energies[1, 0]]) + + # D vector calculations + # D(1) = mean([DEij(1,2,1), -DEij(2,1,1)]) + D[0] = np.mean([energies[0, 1], -energies[0, 2]]) + + # D(2) = mean([DEij(2,1,2), -DEij(1,2,2)]) + D[1] = np.mean([energies[1, 2], -energies[1, 1]]) + + # D(3) = mean([DEij(1,2,3), -DEij(2,1,3)]) + D[2] = np.mean([energies[2, 1], -energies[2, 2]]) + + J_iso = np.trace(J) / 3 + # based on the grogu output pdf + # traceless symmetric exchange matrix: + # Jxx, Jyy, Jxy, Jxz, Jyz + J_S = np.array([J[0, 0] - J_iso, J[1, 1] - J_iso, J[0, 1], J[0, 2], J[1, 2]]) + + self.J = J * sisl.unit_convert("eV", "meV") + self.J_S = J_S * sisl.unit_convert("eV", "meV") + self.J_iso = J_iso * sisl.unit_convert("eV", "meV") + self.D = D * sisl.unit_convert("eV", "meV") + + def dictionary(self): + out = dict() + out["supercell_shift"] = self.supercell_shift + out["supercell_shift_xyz"] = self.supercell_shift_xyz + out["xyz"] = self.xyz + out["distance"] = self.distance + out["tags"] = self.tags + out["energies"] = self.energies() + + return out + + +class Simulation: + root_node = 0 + + def __init__(self, simulation_parameters: dict) -> None: + self.INFILE = simulation_parameters["INFILE"] + self.EIGFILE = simulation_parameters["INFILE"][:-3] + "EIG" + self.OUTFILE = simulation_parameters["OUTFILE"] + self.SCF_XCF_ORIENTATION = simulation_parameters["SCF_XCF_ORIENTATION"] + self.REF_XCF_ORIENTATIONS = simulation_parameters["REF_XCF_ORIENTATIONS"] + self.KSET = simulation_parameters["KSET"] + self.KDIRS = simulation_parameters["KDIRS"] + self.EBOT = simulation_parameters["EBOT"] + self.ESET = simulation_parameters["ESET"] + self.ESETP = simulation_parameters["ESETP"] + self.PARALLEL_SOLVER_FOR_GK = simulation_parameters["PARALLEL_SOLVER_FOR_GK"] + self.EMAX = 0 + self.AUTOMATIC_EBOT = False + + self.magnetic_entities = [] + self.pairs = [] + + if not self.OUTFILE.endswith(".pickle"): + self.OUTFILE += ".pickle" + + if self.EBOT is None: + self.EBOT = self.read_siesta_emin(self.EIGFILE) + self.AUTOMATIC_EBOT = True + + try: + self.SLURM_ID = environ["SLURM_JOB_ID"] + except: + self.SLURM_ID = "Could not be determined." + # add third direction for off-diagonal anisotropy + for orientation in self.REF_XCF_ORIENTATIONS: + o1: np.array = orientation["vw"][0] + o2: np.array = orientation["vw"][1] + orientation["vw"] = np.vstack((orientation["vw"], (o1 + o2) / np.sqrt(2))) + + # MPI parameters + self.comm = MPI.COMM_WORLD + self.parallel_size = self.comm.Get_size() + self.rank = self.comm.Get_rank() + + def __print__(self) -> None: + """It prints the parameters and the description of the job. + + Args: + self : + It contains the simulations parameters + """ + + print( + "================================================================================================================================================================" + ) + try: + print(f"SLURM job ID: {environ['SLURM_JOB_ID']}") + except: + print("SLURM job ID not found.") + print("Input file: ") + print(self.INFILE) + print("Output file: ") + print(self.OUTFILE) + print( + "Number of nodes in the parallel cluster: ", + self.PARALLEL_SIZE, + ) + if self.PARALLEL_SOLVER_FOR_GK: + print("solver used for Greens function calculation: parallel") + else: + print("solver used for Greens function calculation: sequential") + print( + "================================================================================================================================================================" + ) + print("Cell [Ang]: ") + print("NOT YET") + print( + "================================================================================================================================================================" + ) + print("DFT axis: ") + print(self.SCF_XCF_ORIENTATION) + print("Quantization axis and perpendicular rotation directions:") + for ref in self.REF_XCF_ORIENTATIONS: + print(ref["o"], " --» ", ref["vw"]) + print( + "================================================================================================================================================================" + ) + print("Parameters for the contour integral:") + print("Number of k points: ", self.KSET) + print("k point directions: ", self.KDIRS) + if self.AUTOMATIC_EBOT: + print( + "Ebot: ", + self.EBOT, + " WARNING: This was automatically determined!", + ) + else: + print("Ebot: ", self.EBOT) + print("Eset: ", self.ESET) + print("Esetp: ", self.ESETP) + print( + "================================================================================================================================================================" + ) + + def add_KSpace(self, kspace: KSpace): + self.KSpace = kspace + + def add_Hamiltonian(self, hamiltonian: Hamiltonian): + self.Hamiltonian = hamiltonian + self.sislHamiltonian = hamiltonian.sislHamiltonian + + def add_Contour(self, contour: Contour): + self.Contour = contour + + def add_magnetic_entities(self, magnetic_entities: Union[MagneticEntity, list]): + if isinstance(magnetic_entities, MagneticEntity): + magnetic_entities = [magnetic_entities] + + self.magnetic_entities = magnetic_entities + + def add_pairs(self, pairs: Union[Pair, dict, list]): + if isinstance(pairs, Pair): + self.pairs.append(pairs) + + elif isinstance(pairs, dict): + self.pairs.append( + Pair( + self.magnetic_entities[pairs["ai"]], + self.magnetic_entities[pairs["aj"]], + self.sislHamiltonian, + pairs["Ruc"], + ) + ) + + elif isinstance(pairs, list): + for pair in pairs: + if isinstance(pair, Pair): + self.pairs.append(pair) + elif isinstance(pair, dict): + self.pairs.append( + Pair( + self.magnetic_entities[pair["ai"]], + self.magnetic_entities[pair["aj"]], + self.sislHamiltonian, + pair["Ruc"], + ) + ) + else: + raise Exception("Bad input") + else: + raise Exception("Bad input") + + def read_siesta_emin(self, eigfile: str) -> float: + """It reads the lowest energy level from the siesta run. + + It uses the .EIG file from siesta that contains the eigenvalues. + + Args: + eigfile : str + The path to the .EIG file + + Returns: + float + The energy minimum + """ + + # read the file + eigenvalues = sisl.io.eigSileSiesta(eigfile).read_data() + + return eigenvalues.min() + + def rotate_XCF(self): + self.rotated_hamiltonians = [] + for mag_ent in self.magnetic_entities: + for i in range(len(self.REF_XCF_ORIENTATIONS)): + mag_ent.Vu1.append([]) + mag_ent.Vu2.append([]) + mag_ent.Gii.append( + np.zeros((self.ESET, mag_ent.SBS, mag_ent.SBS), dtype="complex128") + ) + mag_ent._Gii_tmp.append( + np.zeros((self.ESET, mag_ent.SBS, mag_ent.SBS), dtype="complex128") + ) + for pair in self.pairs: + for i in range(len(self.REF_XCF_ORIENTATIONS)): + pair.Gij.append( + np.zeros((self.ESET, pair.SBS1, pair.SBS2), dtype="complex128") + ) + pair._Gij_tmp.append( + np.zeros((self.ESET, pair.SBS1, pair.SBS2), dtype="complex128") + ) + pair.Gji.append( + np.zeros((self.ESET, pair.SBS2, pair.SBS1), dtype="complex128") + ) + pair._Gji_tmp.append( + np.zeros((self.ESET, pair.SBS2, pair.SBS1), dtype="complex128") + ) + # iterate over the reference directions (quantization axes) + for i, orient in enumerate(self.REF_XCF_ORIENTATIONS): + rot_H, rot_H_XCF_uc = self.Hamiltonian.rot_H_and_XCF( + self.SCF_XCF_ORIENTATION, orient["o"] + ) + + # store the relevant information of the Hamiltonian + current_hamiltonian = Hamiltonian() + current_hamiltonian.H = rot_H + current_hamiltonian.S = self.Hamiltonian.S + current_hamiltonian.sc_off = self.Hamiltonian.sc_off + + self.rotated_hamiltonians.append(current_hamiltonian) + + # these are the rotations (for now) perpendicular to the quantization axis + for u in orient["vw"]: + # section 2.H + Tu: np.array = np.kron(np.eye(self.Hamiltonian.NO, dtype=int), tau_u(u)) + Vu1, Vu2 = calc_Vu(rot_H_XCF_uc, Tu) + + for mag_ent in self.magnetic_entities: + # fill up the perturbed potentials (for now) based on the on-site projections + mag_ent.Vu1[i].append( + onsite_projection( + Vu1, mag_ent.spin_box_indices, mag_ent.spin_box_indices + ) + ) + mag_ent.Vu2[i].append( + onsite_projection( + Vu2, mag_ent.spin_box_indices, mag_ent.spin_box_indices + ) + ) + + def solve(self): + # sampling the integrand on the contour and the BZ + for i, k in enumerate(self.KSpace.parallel_kpoints[self.rank]): + # weight of k point in BZ integral + wk: float = self.KSpace.parallel_weights[self.rank][i] + # iterate over reference directions + for i, hamiltonian_orientation in enumerate(self.rotated_hamiltonians): + # calculate Hamiltonian and Overlap matrix in a given k point + hamiltonian_orientation.setup_at_k_point(k) + + if self.PARALLEL_SOLVER_FOR_GK: + Gk = hamiltonian_orientation._parallel_Gk( + self.Contour.samples, self.Contour.eset + ) + else: + # solve Greens function sequentially for the energies, because of memory bound + Gk = hamiltonian_orientation._sequential_Gk( + self.Contour.samples, self.Contour.eset + ) + + # store the Greens function slice of the magnetic entities + for mag_ent in self.magnetic_entities: + mag_ent._Gii_tmp[i] += ( + onsite_projection( + Gk, mag_ent.spin_box_indices, mag_ent.spin_box_indices + ) + * wk + ) + + for pair in self.pairs: + # add phase shift based on the cell difference + phase: np.array = np.exp( + 1j * 2 * np.pi * k @ pair.supercell_shift.T + ) + + # store the Greens function slice of the magnetic entities + pair._Gij_tmp[i] += ( + onsite_projection( + Gk, pair.M1.spin_box_indices, pair.M2.spin_box_indices + ) + * phase + * wk + ) + pair._Gji_tmp[i] += ( + onsite_projection( + Gk, pair.M2.spin_box_indices, pair.M1.spin_box_indices + ) + / phase + * wk + ) + + # sum reduce partial results of mpi nodes + for i in range(len(self.rotated_hamiltonians)): + for mag_ent in self.magnetic_entities: + self.comm.Reduce( + mag_ent._Gii_tmp[i], mag_ent.Gii[i], root=self.root_node + ) + + for pair in self.pairs: + self.comm.Reduce(pair._Gij_tmp[i], pair.Gij[i], root=self.root_node) + self.comm.Reduce(pair._Gji_tmp[i], pair.Gji[i], root=self.root_node) + + for mag_ent in self.magnetic_entities: + mag_ent.calculate_energies(self.Contour.weights) + mag_ent.calculate_anisotropy() + for pair in self.pairs: + pair.calculate_energies(self.Contour.weights) + pair.calculate_exchange_tensor() diff --git a/src/grogupy/plotter.py b/src/grogupy/plotter.py new file mode 100644 index 0000000..42cc760 --- /dev/null +++ b/src/grogupy/plotter.py @@ -0,0 +1,25 @@ +# Copyright (c) [2024] [] +# +# 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. + +import numpy as np + + +def plot(): + return np.radians(2) diff --git a/src/grogupy/utilities.py b/src/grogupy/utilities.py index eebaa3b..9be8c27 100644 --- a/src/grogupy/utilities.py +++ b/src/grogupy/utilities.py @@ -24,155 +24,64 @@ from typing import Union import numpy as np -from scipy.special import roots_legendre -from sisl.io.siesta import eigSileSiesta from grogupy.constants import TAU_0, TAU_X, TAU_Y, TAU_Z -def commutator(a: np.array, b: np.array) -> np.array: - """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 spin_tracer(M: np.array) -> dict: + """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. -# define some useful functions -def hsk( - H: np.array, ss: np.array, sc_off: list, k: tuple = (0, 0, 0) -) -> tuple[np.array, np.array]: - """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) + M : np.array_like + Traceable matrix Returns: - np.array_like - Hamiltonian at the given k point - np.array_like - Overlap matrix at the given k point + 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] - # 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) + M_o = dict() + M_o["x"] = M12 + M21 + M_o["y"] = 1j * (M12 - M21) + M_o["z"] = M11 - M22 + M_o["c"] = M11 + M22 - # phases applied to the hamiltonian - HK = np.einsum("abc,a->bc", H, phases) - SK = np.einsum("abc,a->bc", ss, phases) + return M_o - return HK, SK +def commutator(a: np.array, b: np.array) -> np.array: + """Shorthand for commutator. -def make_kset(dirs: str = "xyz", NUMK: int = 20) -> np.array: - """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 returned. The total number of k points is the NUMK**(dimensions) - + Commutator of two matrices in the mathematical sense. 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 + a : np.array_like + The first matrix + b : np.array_like + The second matrix Returns: 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]]) - - kran = len(dirs) * [np.linspace(0, 1, NUMK, endpoint=False)] - mg = np.meshgrid(*kran) - directions = dict() - - for d in enumerate(dirs): - directions[d[1]] = mg[d[0]].flatten() - for d in "xyz": - if not (d in dirs): - directions[d] = 0 * directions[dirs[0]] - kset = np.array([directions[d] for d in "xyz"]).T - - return kset - - -# just an empty container class -class Container: - pass - - -def make_contour( - emin: float = -20, emax: float = 0.0, enum: int = 42, p: float = 150 -) -> Container: - """A more sophisticated contour generator. - - 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: - 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: - Container - Contains all the information for the contour integral + The commutator of a and b """ - 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 + return a @ b - b @ a - cont = Container() - cont.R = R - cont.z0 = z0 - cont.ze = ze - cont.we = we - cont.enum = enum - return cont +# define some useful functions def tau_u(u: Union[list, np.array]) -> np.array: @@ -267,44 +176,3 @@ def RotMa2b(a: np.array, b: np.array, eps: float = 1e-10) -> np.array: # kill off small numbers M[abs(M) < eps] = 0.0 return M - - -def read_siesta_emin(eigfile: str) -> float: - """It reads the lowest energy level from the siesta run. - - It uses the .EIG file from siesta that contains the eigenvalues. - - Args: - eigfile : str - The path to the .EIG file - - Returns: - float - The energy minimum - """ - - # read the file - eigenvalues = eigSileSiesta(eigfile).read_data() - - return eigenvalues.min() - - -def int_de_ke(traced: np.array, we: float) -> np.array: - """It numerically integrates the traced matrix. - - It is a wrapper from numpy.trapz and it contains the - relevant constants to calculate the energy integral from - equation 93 or 96. - - Args: - traced : np.array_like - The trace of a matrix or a matrix product - we : float - The weight of a point on the contour - - Returns: - float - The energy calculated from the integral formula - """ - - return np.trapz(-1 / np.pi * np.imag(traced * we)) diff --git a/test.ipynb b/test.ipynb index 17a8f97..cc58468 100644 --- a/test.ipynb +++ b/test.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 24, + "execution_count": 141, "metadata": {}, "outputs": [], "source": [ @@ -16,7 +16,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 142, "metadata": {}, "outputs": [ { @@ -58,29 +58,549 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": null, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(array([[-8.20799665e+01, 1.42696825e-04, -8.30905002e-04],\n", + " [ 1.41821186e-04, -8.20626555e+01, -1.88623702e-01],\n", + " [ 6.97129733e-04, 1.88607418e-01, -8.92465137e+01]]), -84.46304519743876, array([[ 2.38307873e+00, 1.42259005e-04, -6.68876347e-05],\n", + " [ 1.42259005e-04, 2.40038974e+00, -8.14185600e-06],\n", + " [-6.68876347e-05, -8.14185600e-06, -4.78346847e+00]]), array([[ 0.00000000e+00, 4.37819604e-07, -7.64017368e-04],\n", + " [-4.37819604e-07, 0.00000000e+00, -1.88615560e-01],\n", + " [ 7.64017368e-04, 1.88615560e-01, 0.00000000e+00]]))\n", + "[[-8.20799665e+01 1.42259005e-04 -6.68876347e-05]\n", + " [ 1.42259005e-04 -8.20626555e+01 -8.14185600e-06]\n", + " [-6.68876347e-05 -8.14185600e-06 -8.92465137e+01]] -84.46304519743875 [ 2.38307873e+00 2.40038974e+00 1.42259005e-04 -6.68876347e-05\n", + " -8.14185600e-06] [ 1.88615560e-01 -7.64017368e-04 -4.37819604e-07]\n" + ] + }, + { + "ename": "LinAlgError", + "evalue": "Singular matrix", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mLinAlgError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[187], line 60\u001b[0m\n\u001b[1;32m 57\u001b[0m M \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mouter(ekkprime, ekprimek)\n\u001b[1;32m 58\u001b[0m V \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m E[j] \u001b[38;5;241m*\u001b[39m ekkprime\u001b[38;5;241m.\u001b[39mflatten()\n\u001b[0;32m---> 60\u001b[0m K \u001b[38;5;241m=\u001b[39m \u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlinalg\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msolve\u001b[49m\u001b[43m(\u001b[49m\u001b[43mM\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mV\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m<__array_function__ internals>:200\u001b[0m, in \u001b[0;36msolve\u001b[0;34m(*args, **kwargs)\u001b[0m\n", + "File \u001b[0;32m~/Documents/oktatás/elte/phd/grogu_project/.venv/lib/python3.9/site-packages/numpy/linalg/linalg.py:386\u001b[0m, in \u001b[0;36msolve\u001b[0;34m(a, b)\u001b[0m\n\u001b[1;32m 384\u001b[0m signature \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mDD->D\u001b[39m\u001b[38;5;124m'\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m isComplexType(t) \u001b[38;5;28;01melse\u001b[39;00m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mdd->d\u001b[39m\u001b[38;5;124m'\u001b[39m\n\u001b[1;32m 385\u001b[0m extobj \u001b[38;5;241m=\u001b[39m get_linalg_error_extobj(_raise_linalgerror_singular)\n\u001b[0;32m--> 386\u001b[0m r \u001b[38;5;241m=\u001b[39m \u001b[43mgufunc\u001b[49m\u001b[43m(\u001b[49m\u001b[43ma\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mb\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msignature\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43msignature\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mextobj\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mextobj\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 388\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m wrap(r\u001b[38;5;241m.\u001b[39mastype(result_t, copy\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m))\n", + "File \u001b[0;32m~/Documents/oktatás/elte/phd/grogu_project/.venv/lib/python3.9/site-packages/numpy/linalg/linalg.py:89\u001b[0m, in \u001b[0;36m_raise_linalgerror_singular\u001b[0;34m(err, flag)\u001b[0m\n\u001b[1;32m 88\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_raise_linalgerror_singular\u001b[39m(err, flag):\n\u001b[0;32m---> 89\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m LinAlgError(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mSingular matrix\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n", + "\u001b[0;31mLinAlgError\u001b[0m: Singular matrix" + ] + } + ], + "source": [ + "import pickle\n", + "import numpy as np\n", + "\n", + "with open(\"./Fe3GeTe2_fdf_test.pickle\", \"rb\") as file:\n", + " out = pickle.load(file)\n", + "\n", + "ref_xcf = out[\"parameters\"][\"ref_xcf_orientations\"]\n", + "energies = out[\"pairs\"][0][\"energies\"]\n", + "\n", + "\n", + "def fit_energies(energies, ref_xcf):\n", + " M = np.zeros((9, 9))\n", + " V = np.zeros(9)\n", + "\n", + " for i in range(len(ref_xcf)):\n", + " E = energies[i]\n", + " vw = ref_xcf[i][\"vw\"][:2]\n", + " o = ref_xcf[i][\"o\"]\n", + "\n", + " vw = np.array([[vw[0], vw[0]], [vw[0], vw[1]], [vw[1], vw[0]], [vw[1], vw[1]]])\n", + "\n", + " for j in range(len(vw)):\n", + " aa = np.cross(vw[j, 0], o)\n", + " bb = np.cross(vw[j, 1], o)\n", + " ekkprime = np.outer(aa, bb)\n", + " ekprimek = np.outer(bb, aa)\n", + " M += np.outer(ekkprime, ekprimek)\n", + " V += E[j] * ekkprime.flatten()\n", + "\n", + " J = np.linalg.solve(M, V)\n", + " J = J.reshape(3, 3) * sisl.unit_convert(\"eV\", \"meV\")\n", + " J_s = 0.5 * (J + J.T)\n", + " D = 0.5 * (J - J.T)\n", + " J_iso = np.trace(J_s) / 3\n", + " J_S = J_s - np.eye(3) * J_iso\n", + "\n", + " return J, J_iso, J_S, D\n", + "\n", + "\n", + "print(fit_energies(energies, ref_xcf))\n", + "print(\n", + " out[\"pairs\"][0][\"J\"],\n", + " out[\"pairs\"][0][\"J_iso\"],\n", + " out[\"pairs\"][0][\"J_S\"],\n", + " out[\"pairs\"][0][\"D\"],\n", + ")\n", + "\n", + "\n", + "energies = out[\"magnetic_entities\"][0][\"energies\"]\n", + "M = np.zeros((9, 9))\n", + "V = np.zeros(9)\n", + "\n", + "for i in range(len(ref_xcf)):\n", + " E = energies[i]\n", + " vw = ref_xcf[i][\"vw\"]\n", + " o = ref_xcf[i][\"o\"]\n", + "\n", + " vw = np.array([[vw[0], vw[0]], [vw[1], vw[1]], [vw[2], vw[2]]])\n", + "\n", + " for j in range(len(vw)):\n", + " aa = np.cross(vw[j, 0], o)\n", + " bb = np.cross(vw[j, 1], o)\n", + " ekkprime = np.outer(aa, bb)\n", + " ekprimek = np.outer(bb, aa)\n", + " M += np.outer(ekkprime, ekprimek)\n", + " V += E[j] * ekkprime.flatten()\n", + "\n", + "K = np.linalg.solve(M, V)" + ] + }, + { + "cell_type": "code", + "execution_count": 144, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Matrix([[a1], [a2], [a3]])\n", + "Matrix([[b1], [b2], [b3]])\n", + "Matrix([[J11, J12, J13], [J21, J22, J23], [J31, J32, J33]])\n", + "J11*a1*b1 + J12*a1*b2 + J13*a1*b3 + J21*a2*b1 + J22*a2*b2 + J23*a2*b3 + J31*a3*b1 + J32*a3*b2 + J33*a3*b3\n", + "J11*a1*b1 + J12*a1*b2 + J13*a1*b3 + J21*a2*b1 + J22*a2*b2 + J23*a2*b3 + J31*a3*b1 + J32*a3*b2 + J33*a3*b3\n", + "Matrix([[J11], [J12], [J13], [J21], [J22], [J23], [J31], [J32], [J33]])\n", + "Matrix([[J11*a1*b1 + J12*a1*b2 + J13*a1*b3 + J21*a2*b1 + J22*a2*b2 + J23*a2*b3 + J31*a3*b1 + J32*a3*b2 + J33*a3*b3]])\n", + "Matrix([[a1*b1], [a1*b2], [a1*b3], [a2*b1], [a2*b2], [a2*b3], [a3*b1], [a3*b2], [a3*b3]])\n" + ] + } + ], + "source": [ + "import sympy as sy\n", + "\n", + "a = sy.symbols(\"a1:4\")\n", + "b = sy.symbols(\"b1:4\")\n", + "J = sy.symbols(\"J(1:4)(1:4)\")\n", + "\n", + "\n", + "aa = sy.Matrix(a)\n", + "bb = sy.Matrix(b)\n", + "\n", + "print(aa)\n", + "print(bb)\n", + "\n", + "JJ = sy.Matrix(J).reshape(3, 3)\n", + "\n", + "print(JJ)\n", + "\n", + "print((aa.T @ JJ @ bb)[0].expand())\n", + "print(sy.trace(JJ @ (bb @ aa.T)))\n", + "print(JJ.reshape(9, 1))\n", + "print(JJ.reshape(1, 9) @ (aa @ bb.T).reshape(9, 1))\n", + "print((aa @ bb.T).reshape(9, 1))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 151, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([[ 0., 1., 0.],\n", - " [ 0., 0., 1.],\n", - " [ 100., 1000., 100.]])" + "{'parameters': {'infile': '/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf',\n", + " 'outfile': './Fe3GeTe2_fdf_test.pickle',\n", + " 'scf_xcf_orientation': array([0., 0., 1.]),\n", + " 'ref_xcf_orientations': [{'o': array([1., 0., 0.]),\n", + " 'vw': array([[0. , 1. , 0. ],\n", + " [0. , 0. , 1. ],\n", + " [0. , 0.70710678, 0.70710678]])},\n", + " {'o': array([0., 1., 0.]),\n", + " 'vw': array([[1. , 0. , 0. ],\n", + " [0. , 0. , 1. ],\n", + " [0.70710678, 0. , 0.70710678]])},\n", + " {'o': array([0., 0., 1.]),\n", + " 'vw': array([[1. , 0. , 0. ],\n", + " [0. , 1. , 0. ],\n", + " [0.70710678, 0.70710678, 0. ]])}],\n", + " 'kset': 10,\n", + " 'kdirs': 'xy',\n", + " 'ebot': -12.906878959999998,\n", + " 'eset': 600,\n", + " 'esetp': 1000.0,\n", + " 'parallel_solver_for_Gk': False,\n", + " 'padawan_mode': True,\n", + " 'parallel_size': 1,\n", + " 'automatic_ebot': True,\n", + " 'cell': array([[ 3.79100000e+00, 0.00000000e+00, 0.00000000e+00],\n", + " [-1.89550000e+00, 3.28310231e+00, 0.00000000e+00],\n", + " [ 1.25954923e-15, 2.18160327e-15, 2.05700000e+01]])},\n", + " 'magnetic_entities': [{'atom': 3,\n", + " 'l': [2],\n", + " 'orbital_indices': array([41, 42, 43, 44, 45, 46, 47, 48, 49, 50], dtype=int32),\n", + " 'spin_box_indices': array([ 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94,\n", + " 95, 96, 97, 98, 99, 100, 101]),\n", + " 'tags': ['[3]Fe([2])'],\n", + " 'xyz': [array([-7.33915874e-06, 4.14927851e-06, 1.16575858e+01])],\n", + " 'energies': array([[1.17772527, 1.17771974, 1.17772812],\n", + " [1.17764184, 1.177567 , 1.17761357],\n", + " [1.17892446, 1.17892335, 1.17892284]]),\n", + " 'K': array([[-0.07483657, 0.00106563, -0.0091485 ],\n", + " [ 0.00106563, -0.00553059, -0.00561574],\n", + " [-0.0091485 , -0.00561574, 0. ]]),\n", + " 'K_consistency': 6.819575556149537e-05},\n", + " {'atom': 4,\n", + " 'l': [2],\n", + " 'orbital_indices': array([56, 57, 58, 59, 60, 61, 62, 63, 64, 65], dtype=int32),\n", + " 'spin_box_indices': array([112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124,\n", + " 125, 126, 127, 128, 129, 130, 131]),\n", + " 'tags': ['[4]Fe([2])'],\n", + " 'xyz': [array([-7.32698766e-06, 4.15827452e-06, 8.91242254e+00])],\n", + " 'energies': array([[1.17774 , 1.17766249, 1.17769662],\n", + " [1.1776063 , 1.17760676, 1.17759665],\n", + " [1.17893438, 1.17893546, 1.17893594]]),\n", + " 'K': array([[ 0.00045681, -0.00101949, 0.00987698],\n", + " [-0.00101949, -0.07750577, 0.00462617],\n", + " [ 0.00987698, 0.00462617, 0. ]]),\n", + " 'K_consistency': 7.687732837413641e-05},\n", + " {'atom': 5,\n", + " 'l': [2],\n", + " 'orbital_indices': array([71, 72, 73, 74, 75, 76, 77, 78, 79, 80], dtype=int32),\n", + " 'spin_box_indices': array([142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154,\n", + " 155, 156, 157, 158, 159, 160, 161]),\n", + " 'tags': ['[5]Fe([2])'],\n", + " 'xyz': [array([ 1.89546671, 1.09439132, 10.2850027 ])],\n", + " 'energies': array([[0.82604465, 0.82600318, 0.82602942],\n", + " [0.82596082, 0.82597202, 0.82596098],\n", + " [0.83074663, 0.83074663, 0.83074663]]),\n", + " 'K': array([[ 1.12053382e-02, -4.12902668e-06, 5.44290076e-03],\n", + " [-4.12902668e-06, -4.14783466e-02, -5.50288071e-03],\n", + " [ 5.44290076e-03, -5.50288071e-03, 0.00000000e+00]]),\n", + " 'K_consistency': 5.2686617336705766e-05}],\n", + " 'pairs': [{'ai': 0,\n", + " 'aj': 1,\n", + " 'Ruc': array([0, 0, 0]),\n", + " 'dist': 2.745163300331324,\n", + " 'tags': ['[3]Fe([2])', '[4]Fe([2])'],\n", + " 'energies': array([[-8.92392323e-02, 1.88623702e-04, -1.88607418e-04,\n", + " -8.91422199e-02],\n", + " [-8.92537950e-02, 8.30905002e-07, -6.97129733e-07,\n", + " -8.91766857e-02],\n", + " [-7.49830910e-02, -1.42696825e-07, -1.41821186e-07,\n", + " -7.49832472e-02]]),\n", + " 'J_iso': -84.46304519743875,\n", + " 'J_S': array([ 2.38307873e+00, 2.40038974e+00, 1.42259005e-04, -6.68876347e-05,\n", + " -8.14185600e-06]),\n", + " 'D': array([ 1.88615560e-01, -7.64017368e-04, -4.37819604e-07]),\n", + " 'J': array([[-8.20799665e+01, 1.42259005e-04, -6.68876347e-05],\n", + " [ 1.42259005e-04, -8.20626555e+01, -8.14185600e-06],\n", + " [-6.68876347e-05, -8.14185600e-06, -8.92465137e+01]])},\n", + " {'ai': 0,\n", + " 'aj': 2,\n", + " 'Ruc': array([0, 0, 0]),\n", + " 'dist': 2.5835033632437767,\n", + " 'tags': ['[3]Fe([2])', '[5]Fe([2])'],\n", + " 'energies': array([[-0.04075836, 0.00119288, -0.00112921, -0.04080086],\n", + " [-0.04056911, 0.00207638, -0.00198342, -0.04056736],\n", + " [-0.04367583, -0.00022997, -0.00023025, -0.04340975]]),\n", + " 'J_iso': -41.630212919435856,\n", + " 'J_S': array([-0.35834154, -0.60813515, 0.23011103, -0.046484 , -0.03183377]),\n", + " 'D': array([ 1.16104682e+00, -2.02989980e+00, 1.42005544e-04]),\n", + " 'J': array([[-4.19885545e+01, 2.30111030e-01, -4.64839962e-02],\n", + " [ 2.30111030e-01, -4.22383481e+01, -3.18337650e-02],\n", + " [-4.64839962e-02, -3.18337650e-02, -4.06637362e+01]])},\n", + " {'ai': 1,\n", + " 'aj': 2,\n", + " 'Ruc': array([0, 0, 0]),\n", + " 'dist': 2.583501767937866,\n", + " 'tags': ['[4]Fe([2])', '[5]Fe([2])'],\n", + " 'energies': array([[-0.04077002, -0.00117329, 0.00113811, -0.04082399],\n", + " [-0.04057372, -0.00204387, 0.00197878, -0.0405759 ],\n", + " [-0.04367653, -0.00022997, -0.00023026, -0.04341045]]),\n", + " 'J_iso': -41.63843353643032,\n", + " 'J_S': array([-0.35473711, -0.61182503, 0.23011277, 0.03254538, 0.01758745]),\n", + " 'D': array([-1.15570028e+00, 2.01132360e+00, 1.42811452e-04]),\n", + " 'J': array([[-4.19931707e+01, 2.30112766e-01, 3.25453781e-02],\n", + " [ 2.30112766e-01, -4.22502586e+01, 1.75874450e-02],\n", + " [ 3.25453781e-02, 1.75874450e-02, -4.06718714e+01]])},\n", + " {'ai': 0,\n", + " 'aj': 2,\n", + " 'Ruc': array([-1, -1, 0]),\n", + " 'dist': 2.5834973202859075,\n", + " 'tags': ['[3]Fe([2])', '[5]Fe([2])'],\n", + " 'energies': array([[-4.05164545e-02, -2.37526408e-03, 2.29370374e-03,\n", + " -4.04788963e-02],\n", + " [-4.09074301e-02, -2.05208851e-07, -8.74099273e-08,\n", + " -4.09801092e-02],\n", + " [-4.32782877e-02, 1.10167960e-07, 1.95301791e-07,\n", + " -4.38099217e-02]]),\n", + " 'J_iso': -41.66184991462266,\n", + " 'J_S': array([-7.33165557e-01, -2.16742080e-01, -1.52734875e-04, 1.46309389e-04,\n", + " 4.07801676e-02]),\n", + " 'D': array([-2.33448391e+00, 5.88994618e-05, -4.25669156e-05]),\n", + " 'J': array([[-4.23950155e+01, -1.52734875e-04, 1.46309389e-04],\n", + " [-1.52734875e-04, -4.18785920e+01, 4.07801676e-02],\n", + " [ 1.46309389e-04, 4.07801676e-02, -4.07119423e+01]])},\n", + " {'ai': 1,\n", + " 'aj': 2,\n", + " 'Ruc': array([-1, -1, 0]),\n", + " 'dist': 2.583495745338251,\n", + " 'tags': ['[4]Fe([2])', '[5]Fe([2])'],\n", + " 'energies': array([[-4.05171468e-02, 2.37529161e-03, -2.29373272e-03,\n", + " -4.04795842e-02],\n", + " [-4.08811902e-02, 1.57284973e-07, -4.45157439e-07,\n", + " -4.09421171e-02],\n", + " [-4.32789799e-02, 1.09351818e-07, 1.96420584e-07,\n", + " -4.38106181e-02]]),\n", + " 'J_iso': -41.651606039885934,\n", + " 'J_S': array([-7.24761554e-01, -2.27675997e-01, -1.52886201e-04, 1.43936233e-04,\n", + " -4.07794409e-02]),\n", + " 'D': array([ 2.33451216e+00, -3.01221206e-04, -4.35343826e-05]),\n", + " 'J': array([[-4.23763676e+01, -1.52886201e-04, 1.43936233e-04],\n", + " [-1.52886201e-04, -4.18792820e+01, -4.07794409e-02],\n", + " [ 1.43936233e-04, -4.07794409e-02, -4.06991685e+01]])},\n", + " {'ai': 0,\n", + " 'aj': 2,\n", + " 'Ruc': array([-1, 0, 0]),\n", + " 'dist': 2.583541444641373,\n", + " 'tags': ['[3]Fe([2])', '[5]Fe([2])'],\n", + " 'energies': array([[-0.04075762, 0.00117311, -0.00113797, -0.04081104],\n", + " [-0.04055754, -0.00207635, 0.00198375, -0.04055565],\n", + " [-0.04366146, 0.00022983, 0.00023003, -0.04339595]]),\n", + " 'J_iso': -41.62320943278901,\n", + " 'J_S': array([-0.35258884, -0.61303902, -0.22992743, 0.04629825, -0.01757271]),\n", + " 'D': array([ 1.15553939e+00, 2.03005138e+00, -1.00670740e-04]),\n", + " 'J': array([[-4.19757983e+01, -2.29927428e-01, 4.62982483e-02],\n", + " [-2.29927428e-01, -4.22362485e+01, -1.75727100e-02],\n", + " [ 4.62982483e-02, -1.75727100e-02, -4.06575816e+01]])},\n", + " {'ai': 1,\n", + " 'aj': 2,\n", + " 'Ruc': array([-1, 0, 0]),\n", + " 'dist': 2.5835398672184064,\n", + " 'tags': ['[4]Fe([2])', '[5]Fe([2])'],\n", + " 'energies': array([[-0.04074738, -0.00119276, 0.00112911, -0.04078932],\n", + " [-0.04056216, 0.00204384, -0.00197854, -0.04056432],\n", + " [-0.04366217, 0.00022983, 0.00023003, -0.04339665]]),\n", + " 'J_iso': -41.62033567430795,\n", + " 'J_S': array([-0.36015198, -0.6054113 , -0.22992935, -0.03265044, 0.03182307]),\n", + " 'D': array([-1.16093344e+00, -2.01119076e+00, -9.97127859e-05]),\n", + " 'J': array([[-4.19804877e+01, -2.29929354e-01, -3.26504440e-02],\n", + " [-2.29929354e-01, -4.22257470e+01, 3.18230690e-02],\n", + " [-3.26504440e-02, 3.18230690e-02, -4.06547724e+01]])},\n", + " {'ai': 1,\n", + " 'aj': 2,\n", + " 'Ruc': array([-2, 0, 0]),\n", + " 'dist': 5.951322298958084,\n", + " 'tags': ['[4]Fe([2])', '[5]Fe([2])'],\n", + " 'energies': array([[-1.72611605e-03, -4.36527901e-05, 1.26576311e-04,\n", + " -1.82722454e-03],\n", + " [-1.80427815e-03, -2.63038710e-04, 2.20144926e-04,\n", + " -1.73941292e-03],\n", + " [-1.75263911e-03, 1.36519152e-03, -1.52871309e-03,\n", + " -1.52707483e-03]]),\n", + " 'J_iso': -1.7294575999129687,\n", + " 'J_S': array([ 0.09621372, -0.06047422, 0.08176078, 0.02144689, -0.04146176]),\n", + " 'D': array([-0.08511455, 0.24159182, 1.4469523 ]),\n", + " 'J': array([[-1.63324388, 0.08176078, 0.02144689],\n", + " [ 0.08176078, -1.78993182, -0.04146176],\n", + " [ 0.02144689, -0.04146176, -1.7651971 ]])},\n", + " {'ai': 1,\n", + " 'aj': 2,\n", + " 'Ruc': array([-3, 0, 0]),\n", + " 'dist': 9.638732176310562,\n", + " 'tags': ['[4]Fe([2])', '[5]Fe([2])'],\n", + " 'energies': array([[-7.80160843e-04, 3.56600605e-05, -6.10780433e-05,\n", + " -7.42312414e-04],\n", + " [-7.94983595e-04, -5.77431729e-05, 5.08380963e-05,\n", + " -8.34625148e-04],\n", + " [ 2.17237059e-04, -1.43586609e-04, 1.74512092e-04,\n", + " 2.13496975e-04]]),\n", + " 'J_iso': -0.45355799416161025,\n", + " 'J_S': array([ 0.14299391, 0.19102032, -0.01546274, 0.00345254, 0.01270899]),\n", + " 'D': array([ 0.04836905, 0.05429063, -0.15904935]),\n", + " 'J': array([[-0.31056409, -0.01546274, 0.00345254],\n", + " [-0.01546274, -0.26253768, 0.01270899],\n", + " [ 0.00345254, 0.01270899, -0.78757222]])}],\n", + " 'runtime': {'start_time': 0.881972083,\n", + " 'setup_time': 0.979693,\n", + " 'H_and_XCF_time': 1.361691666,\n", + " 'site_and_pair_dictionaries_time': 1.40248525,\n", + " 'k_set_time': 1.432338791,\n", + " 'reference_rotations_time': 1.68477325,\n", + " 'green_function_inversion_time': 281.406242666,\n", + " 'end_time': 281.949395875}}" ] }, - "execution_count": 26, + "execution_count": 151, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "asd = read_grogupy_fdf(\"input.fdf\")[0][\"ref_xcf_orientations\"][0][\"vw\"]\n", + "out" + ] + }, + { + "cell_type": "code", + "execution_count": 171, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(array([0., 0., 1.]), array([0., 1., 0.]))" + ] + }, + "execution_count": 171, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "scf_xcf_orientation = out[\"parameters\"][\"scf_xcf_orientation\"]\n", + "\n", + "ref_xcf_orientations = out[\"parameters\"][\"ref_xcf_orientations\"]\n", + "\n", + "scf_xcf_orientation\n", + "\n", + "from grogupy import *\n", + "from grogupy.utilities import *\n", + "\n", "\n", - "asd2 = np.array([100, 1000, 100])\n", - "np.vstack((asd, asd2))" + "def generate_perpendicular(orientation, reference):\n", + " x = np.array([1, 0, 0])\n", + " y = np.array([0, 1, 0])\n", + " z = np.array([0, 0, 1])\n", + "\n", + " R_to_reference = RotMa2b(z, reference)\n", + " x = R_to_reference @ x\n", + " y = R_to_reference @ y\n", + " z = R_to_reference @ z\n", + "\n", + " R_to_orientation = RotMa2b(z, orientation)\n", + " x = R_to_orientation @ x\n", + " y = R_to_orientation @ y\n", + "\n", + " return x, y\n", + "\n", + "\n", + "generate_perpendicular(ref_xcf_orientations[0][\"o\"], scf_xcf_orientation)" + ] + }, + { + "cell_type": "code", + "execution_count": 172, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[1. 0. 0.]\n", + "[[0. 1. 0. ]\n", + " [0. 0. 1. ]\n", + " [0. 0.70710678 0.70710678]]\n", + "(array([0., 0., 1.]), array([0., 1., 0.]))\n", + "[0. 1. 0.]\n", + "[[1. 0. 0. ]\n", + " [0. 0. 1. ]\n", + " [0.70710678 0. 0.70710678]]\n", + "(array([-1., 0., 0.]), array([ 0., 0., -1.]))\n", + "[0. 0. 1.]\n", + "[[1. 0. 0. ]\n", + " [0. 1. 0. ]\n", + " [0.70710678 0.70710678 0. ]]\n", + "(array([-1., 0., 0.]), array([0., 1., 0.]))\n" + ] + } + ], + "source": [ + "for ref in ref_xcf_orientations:\n", + " print(ref[\"o\"])\n", + " print(ref[\"vw\"])\n", + " print(generate_perpendicular(ref[\"o\"], scf_xcf_orientation))" ] }, + { + "cell_type": "code", + "execution_count": 182, + "metadata": {}, + "outputs": [ + { + "ename": "LinAlgError", + "evalue": "Singular matrix", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mLinAlgError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[182], line 25\u001b[0m\n\u001b[1;32m 22\u001b[0m M \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mouter(ekkprime, ekprimek)\n\u001b[1;32m 23\u001b[0m V \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m E[j] \u001b[38;5;241m*\u001b[39m ekkprime\u001b[38;5;241m.\u001b[39mflatten()\n\u001b[0;32m---> 25\u001b[0m K \u001b[38;5;241m=\u001b[39m \u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlinalg\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msolve\u001b[49m\u001b[43m(\u001b[49m\u001b[43mM\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mV\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m<__array_function__ internals>:200\u001b[0m, in \u001b[0;36msolve\u001b[0;34m(*args, **kwargs)\u001b[0m\n", + "File \u001b[0;32m~/Documents/oktatás/elte/phd/grogu_project/.venv/lib/python3.9/site-packages/numpy/linalg/linalg.py:386\u001b[0m, in \u001b[0;36msolve\u001b[0;34m(a, b)\u001b[0m\n\u001b[1;32m 384\u001b[0m signature \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mDD->D\u001b[39m\u001b[38;5;124m'\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m isComplexType(t) \u001b[38;5;28;01melse\u001b[39;00m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mdd->d\u001b[39m\u001b[38;5;124m'\u001b[39m\n\u001b[1;32m 385\u001b[0m extobj \u001b[38;5;241m=\u001b[39m get_linalg_error_extobj(_raise_linalgerror_singular)\n\u001b[0;32m--> 386\u001b[0m r \u001b[38;5;241m=\u001b[39m \u001b[43mgufunc\u001b[49m\u001b[43m(\u001b[49m\u001b[43ma\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mb\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msignature\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43msignature\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mextobj\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mextobj\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 388\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m wrap(r\u001b[38;5;241m.\u001b[39mastype(result_t, copy\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m))\n", + "File \u001b[0;32m~/Documents/oktatás/elte/phd/grogu_project/.venv/lib/python3.9/site-packages/numpy/linalg/linalg.py:89\u001b[0m, in \u001b[0;36m_raise_linalgerror_singular\u001b[0;34m(err, flag)\u001b[0m\n\u001b[1;32m 88\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_raise_linalgerror_singular\u001b[39m(err, flag):\n\u001b[0;32m---> 89\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m LinAlgError(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mSingular matrix\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n", + "\u001b[0;31mLinAlgError\u001b[0m: Singular matrix" + ] + } + ], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "code", "execution_count": 53, diff --git a/test_module.ipynb b/test_module.ipynb new file mode 100644 index 0000000..d271842 --- /dev/null +++ b/test_module.ipynb @@ -0,0 +1,312 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "[Daniels-MacBook-Air.local:03366] shmem: mmap: an error occurred while determining whether or not /var/folders/yh/dx7xl94n3g52ts3td8qcxjcc0000gn/T//ompi.Daniels-MacBook-Air.501/jf.0/246480896/sm_segment.Daniels-MacBook-Air.501.eb10000.0 could be created.\n" + ] + } + ], + "source": [ + "from grogupy.magnetism import *\n", + "from grogupy.constants import DEFAULT_PARAMETERS" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "DEFAULT_PARAMETERS[\"INFILE\"] = (\n", + " \"/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf\"\n", + ")\n", + "DEFAULT_PARAMETERS[\"KSET\"] = 10\n", + "my_simulation = Simulation(DEFAULT_PARAMETERS)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + " 0%| | 0/100 [00:00" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "kset = KSpace(DEFAULT_PARAMETERS[\"KDIRS\"], DEFAULT_PARAMETERS[\"KSET\"])\n", + "my_simulation.add_KSpace(kset)\n", + "kset" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "contour = Contour(\n", + " DEFAULT_PARAMETERS[\"EBOT\"], DEFAULT_PARAMETERS[\"ESET\"], DEFAULT_PARAMETERS[\"ESETP\"]\n", + ")\n", + "my_simulation.add_Contour(contour)\n", + "contour" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/danielpozsar/Documents/oktatás/elte/phd/grogu_project/src/grogupy/magnetism.py:170: UserWarning: Hamiltonian changed after symmetrization. Largest change is 1.1920928955078125e-07\n", + " warnings.warn(\n", + "/Users/danielpozsar/Documents/oktatás/elte/phd/grogu_project/src/grogupy/magnetism.py:174: UserWarning: Overlap matrix changed after symmetrization. Largest change is 1.1920928955078125e-07\n", + " warnings.warn(\n" + ] + } + ], + "source": [ + "# read sile\n", + "fdf = sisl.get_sile(DEFAULT_PARAMETERS[\"INFILE\"])\n", + "\n", + "# read in hamiltonian\n", + "dh = fdf.read_hamiltonian()\n", + "hamiltonian = Hamiltonian()\n", + "hamiltonian.build_custom_H_and_S(dh)\n", + "my_simulation.add_Hamiltonian(hamiltonian)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "magnetic_entities = [\n", + " MagneticEntity(my_simulation.sislHamiltonian, atom=3, l=2),\n", + " MagneticEntity(my_simulation.sislHamiltonian, atom=4, l=2),\n", + " MagneticEntity(my_simulation.sislHamiltonian, atom=5, l=2),\n", + "]\n", + "my_simulation.add_magnetic_entities(magnetic_entities)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "my_simulation.add_pairs(\n", + " [\n", + " dict(ai=0, aj=1, Ruc=(0, 0, 0)),\n", + " dict(ai=0, aj=2, Ruc=(0, 0, 0)),\n", + " dict(ai=1, aj=2, Ruc=(0, 0, 0)),\n", + " ]\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[array([0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01,\n", + " 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01,\n", + " 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01,\n", + " 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01,\n", + " 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01,\n", + " 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01,\n", + " 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01,\n", + " 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01,\n", + " 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01,\n", + " 0.01])]" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "my_simulation.KSpace.parallel_weights" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124,\n", + " 125, 126, 127, 128, 129, 130, 131])" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "my_simulation.rotate_XCF()\n", + "my_simulation.pairs[2].M1.spin_box_indices" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 100/100 [01:55<00:00, 1.16s/it]\n" + ] + } + ], + "source": [ + "my_simulation.solve()" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[{'atom': 3, 'l': [2], 'orbital_indices': array([41, 42, 43, 44, 45, 46, 47, 48, 49, 50], dtype=int32), 'spin_box_indices': array([ 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94,\n", + " 95, 96, 97, 98, 99, 100, 101]), 'tags': ['[3]Fe([2])'], 'xyz': [array([-7.33915874e-06, 4.14927851e-06, 1.16575858e+01])], 'energies': array([[1.17772527, 1.17771974, 1.17772812],\n", + " [1.17764184, 1.177567 , 1.17761357],\n", + " [1.17892446, 1.17892335, 1.17892284]]), 'K': array([[-0.07483657, 0.00106563, -0.0091485 ],\n", + " [ 0.00106563, -0.00553059, -0.00561574],\n", + " [-0.0091485 , -0.00561574, 0. ]]), 'K_consistency': 6.819575556149537e-05}, {'atom': 4, 'l': [2], 'orbital_indices': array([56, 57, 58, 59, 60, 61, 62, 63, 64, 65], dtype=int32), 'spin_box_indices': array([112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124,\n", + " 125, 126, 127, 128, 129, 130, 131]), 'tags': ['[4]Fe([2])'], 'xyz': [array([-7.32698766e-06, 4.15827452e-06, 8.91242254e+00])], 'energies': array([[1.17774 , 1.17766249, 1.17769662],\n", + " [1.1776063 , 1.17760676, 1.17759665],\n", + " [1.17893438, 1.17893546, 1.17893594]]), 'K': array([[ 0.00045681, -0.00101949, 0.00987698],\n", + " [-0.00101949, -0.07750577, 0.00462617],\n", + " [ 0.00987698, 0.00462617, 0. ]]), 'K_consistency': 7.687732837413641e-05}, {'atom': 5, 'l': [2], 'orbital_indices': array([71, 72, 73, 74, 75, 76, 77, 78, 79, 80], dtype=int32), 'spin_box_indices': array([142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154,\n", + " 155, 156, 157, 158, 159, 160, 161]), 'tags': ['[5]Fe([2])'], 'xyz': [array([ 1.89546671, 1.09439132, 10.2850027 ])], 'energies': array([[0.82604465, 0.82600318, 0.82602942],\n", + " [0.82596082, 0.82597202, 0.82596098],\n", + " [0.83074663, 0.83074663, 0.83074663]]), 'K': array([[ 1.12053382e-02, -4.12902668e-06, 5.44290076e-03],\n", + " [-4.12902668e-06, -4.14783466e-02, -5.50288071e-03],\n", + " [ 5.44290076e-03, -5.50288071e-03, 0.00000000e+00]]), 'K_consistency': 5.2686617336705766e-05}]\n" + ] + } + ], + "source": [ + "import pickle\n", + "\n", + "with open(\"./Fe3GeTe2_fdf_test.pickle\", \"rb\") as file:\n", + " out = pickle.load(file)\n", + "\n", + "print(out[\"magnetic_entities\"])" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[-8.92392323e-02 1.88623702e-04 -1.88607418e-04 -8.91422199e-02]\n", + " [-8.92537950e-02 8.30905002e-07 -6.97129733e-07 -8.91766857e-02]\n", + " [-7.49830910e-02 -1.42696825e-07 -1.41821186e-07 -7.49832472e-02]]\n" + ] + }, + { + "data": { + "text/plain": [ + "array([[-1.03455039e-05, -3.23134167e-05, 1.51994558e-05,\n", + " 1.47670297e-05],\n", + " [-1.34560335e-04, -1.19127887e-04, 7.37198539e-07,\n", + " -2.50259864e-06],\n", + " [ 1.15312736e-06, -2.78200117e-08, -5.45059115e-07,\n", + " 2.19295071e-06]])" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "print(out[\"pairs\"][0][\"energies\"])\n", + "my_simulation.pairs[0]._energies" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.6" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}