From a0f1e9ece6b3d9aba9bc968196dde7441bc9b986 Mon Sep 17 00:00:00 2001 From: Daniel Pozsar Date: Wed, 30 Oct 2024 14:21:04 +0100 Subject: [PATCH] added black:jupyter to pre-commit, made the output nicer so it is easier to debug code, experimented with no results... :'( --- .pre-commit-config.yaml | 5 + local_operator.ipynb | 183 +++++++++++++------------- src/grogu_magn/grogu.py | 218 ++++++------------------------- src/grogu_magn/useful.py | 194 ++++++++++++++++++++++++---- test.ipynb | 273 +++++++++++++++------------------------ 5 files changed, 410 insertions(+), 463 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index d065163..79885c2 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -27,6 +27,11 @@ repos: hooks: - id: black language_version: python3.9 +- repo: https://github.com/psf/black-pre-commit-mirror + rev: 24.4.2 + hooks: + - id: black-jupyter + language_version: python3.9 #- repo: local # hooks: # - id: pytest-check diff --git a/local_operator.ipynb b/local_operator.ipynb index 968674a..1176a88 100644 --- a/local_operator.ipynb +++ b/local_operator.ipynb @@ -24,48 +24,47 @@ "\n", "# Define MatrixSymbols with specific sizes\n", "# A = 2x2, B = 3x3, R = 4x4\n", - "AA = MatrixSymbol('AA', 2, 2) # 2x2\n", - "AB = MatrixSymbol('AB', 2, 3) # 2x3\n", - "BA = MatrixSymbol('BA', 3, 2) # 3x2\n", - "AR = MatrixSymbol('AR', 2, 4) # 2x4\n", - "RA = MatrixSymbol('RA', 4, 2) # 4x2\n", - "BB = MatrixSymbol('BB', 3, 3) # 3x3\n", - "BR = MatrixSymbol('BR', 3, 4) # 3x4\n", - "RB = MatrixSymbol('RB', 4, 3) # 4x3\n", - "RR = MatrixSymbol('RR', 4, 4) # 4x4\n", + "AA = MatrixSymbol(\"AA\", 2, 2) # 2x2\n", + "AB = MatrixSymbol(\"AB\", 2, 3) # 2x3\n", + "BA = MatrixSymbol(\"BA\", 3, 2) # 3x2\n", + "AR = MatrixSymbol(\"AR\", 2, 4) # 2x4\n", + "RA = MatrixSymbol(\"RA\", 4, 2) # 4x2\n", + "BB = MatrixSymbol(\"BB\", 3, 3) # 3x3\n", + "BR = MatrixSymbol(\"BR\", 3, 4) # 3x4\n", + "RB = MatrixSymbol(\"RB\", 4, 3) # 4x3\n", + "RR = MatrixSymbol(\"RR\", 4, 4) # 4x4\n", "\n", "# Create block matrices VA and VB using BlockMatrix\n", "# We need to use ZeroMatrix with appropriate dimensions for padding\n", - "VA = BlockMatrix([\n", - " [AA, 0.5*AB, 0.5*AR],\n", - " [0.5*BA, ZeroMatrix(3,3), ZeroMatrix(3,4)],\n", - " [0.5*RA, ZeroMatrix(4,3), ZeroMatrix(4,4)]\n", - "])\n", - "\n", - "VB = BlockMatrix([\n", - " [ZeroMatrix(2,2), 0.5*AB, ZeroMatrix(2,4)],\n", - " [0.5*BA, BB, 0.5*BR],\n", - " [ZeroMatrix(4,2), 0.5*RB, ZeroMatrix(4,4)]\n", - "])\n", + "VA = BlockMatrix(\n", + " [\n", + " [AA, 0.5 * AB, 0.5 * AR],\n", + " [0.5 * BA, ZeroMatrix(3, 3), ZeroMatrix(3, 4)],\n", + " [0.5 * RA, ZeroMatrix(4, 3), ZeroMatrix(4, 4)],\n", + " ]\n", + ")\n", + "\n", + "VB = BlockMatrix(\n", + " [\n", + " [ZeroMatrix(2, 2), 0.5 * AB, ZeroMatrix(2, 4)],\n", + " [0.5 * BA, BB, 0.5 * BR],\n", + " [ZeroMatrix(4, 2), 0.5 * RB, ZeroMatrix(4, 4)],\n", + " ]\n", + ")\n", "\n", "# Define G matrix symbols with matching dimensions\n", - "G_AA = MatrixSymbol('G_AA', 2, 2)\n", - "G_AB = MatrixSymbol('G_AB', 2, 3)\n", - "G_BA = MatrixSymbol('G_BA', 3, 2)\n", - "G_AR = MatrixSymbol('G_AR', 2, 4)\n", - "G_RA = MatrixSymbol('G_RA', 4, 2)\n", - "G_BB = MatrixSymbol('G_BB', 3, 3)\n", - "G_BR = MatrixSymbol('G_BR', 3, 4)\n", - "G_RB = MatrixSymbol('G_RB', 4, 3)\n", - "G_RR = MatrixSymbol('G_RR', 4, 4)\n", + "G_AA = MatrixSymbol(\"G_AA\", 2, 2)\n", + "G_AB = MatrixSymbol(\"G_AB\", 2, 3)\n", + "G_BA = MatrixSymbol(\"G_BA\", 3, 2)\n", + "G_AR = MatrixSymbol(\"G_AR\", 2, 4)\n", + "G_RA = MatrixSymbol(\"G_RA\", 4, 2)\n", + "G_BB = MatrixSymbol(\"G_BB\", 3, 3)\n", + "G_BR = MatrixSymbol(\"G_BR\", 3, 4)\n", + "G_RB = MatrixSymbol(\"G_RB\", 4, 3)\n", + "G_RR = MatrixSymbol(\"G_RR\", 4, 4)\n", "\n", "# Create G matrix as a BlockMatrix\n", - "G = BlockMatrix([\n", - " [G_AA, G_AB, G_AR],\n", - " [G_BA, G_BB, G_BR],\n", - " [G_RA, G_RB, G_RR]\n", - "])\n", - "\n", + "G = BlockMatrix([[G_AA, G_AB, G_AR], [G_BA, G_BB, G_BR], [G_RA, G_RB, G_RR]])\n", "\n", "\n", "# Calculate the trace of VA@G@VB@G\n", @@ -92,47 +91,47 @@ ], "source": [ "# Define MatrixSymbols with specific sizes\n", - "AA = MatrixSymbol('AA', 2, 2) # 2x2\n", - "AB = MatrixSymbol('AB', 2, 3) # 2x3\n", - "BA = MatrixSymbol('BA', 3, 2) # 3x2\n", - "BB = MatrixSymbol('BB', 3, 3) # 3x3\n", + "AA = MatrixSymbol(\"AA\", 2, 2) # 2x2\n", + "AB = MatrixSymbol(\"AB\", 2, 3) # 2x3\n", + "BA = MatrixSymbol(\"BA\", 3, 2) # 3x2\n", + "BB = MatrixSymbol(\"BB\", 3, 3) # 3x3\n", "\n", "# Create block matrices VA and VB using BlockMatrix\n", "# R-related terms are replaced with zero matrices\n", - "VA = BlockMatrix([\n", - " [AA, 0.5*AB, ZeroMatrix(2,4)],\n", - " [0.5*BA, ZeroMatrix(3,3), ZeroMatrix(3,4)],\n", - " [ZeroMatrix(4,2), ZeroMatrix(4,3), ZeroMatrix(4,4)]\n", - "])\n", - "\n", - "VB = BlockMatrix([\n", - " [ZeroMatrix(2,2), 0.5*AB, ZeroMatrix(2,4)],\n", - " [0.5*BA, BB, ZeroMatrix(3,4)],\n", - " [ZeroMatrix(4,2), ZeroMatrix(4,3), ZeroMatrix(4,4)]\n", - "])\n", + "VA = BlockMatrix(\n", + " [\n", + " [AA, 0.5 * AB, ZeroMatrix(2, 4)],\n", + " [0.5 * BA, ZeroMatrix(3, 3), ZeroMatrix(3, 4)],\n", + " [ZeroMatrix(4, 2), ZeroMatrix(4, 3), ZeroMatrix(4, 4)],\n", + " ]\n", + ")\n", + "\n", + "VB = BlockMatrix(\n", + " [\n", + " [ZeroMatrix(2, 2), 0.5 * AB, ZeroMatrix(2, 4)],\n", + " [0.5 * BA, BB, ZeroMatrix(3, 4)],\n", + " [ZeroMatrix(4, 2), ZeroMatrix(4, 3), ZeroMatrix(4, 4)],\n", + " ]\n", + ")\n", "\n", "# Define G matrix symbols with matching dimensions (kept original)\n", - "G_AA = MatrixSymbol('G_AA', 2, 2)\n", - "G_AB = MatrixSymbol('G_AB', 2, 3)\n", - "G_BA = MatrixSymbol('G_BA', 3, 2)\n", - "G_AR = MatrixSymbol('G_AR', 2, 4)\n", - "G_RA = MatrixSymbol('G_RA', 4, 2)\n", - "G_BB = MatrixSymbol('G_BB', 3, 3)\n", - "G_BR = MatrixSymbol('G_BR', 3, 4)\n", - "G_RB = MatrixSymbol('G_RB', 4, 3)\n", - "G_RR = MatrixSymbol('G_RR', 4, 4)\n", + "G_AA = MatrixSymbol(\"G_AA\", 2, 2)\n", + "G_AB = MatrixSymbol(\"G_AB\", 2, 3)\n", + "G_BA = MatrixSymbol(\"G_BA\", 3, 2)\n", + "G_AR = MatrixSymbol(\"G_AR\", 2, 4)\n", + "G_RA = MatrixSymbol(\"G_RA\", 4, 2)\n", + "G_BB = MatrixSymbol(\"G_BB\", 3, 3)\n", + "G_BR = MatrixSymbol(\"G_BR\", 3, 4)\n", + "G_RB = MatrixSymbol(\"G_RB\", 4, 3)\n", + "G_RR = MatrixSymbol(\"G_RR\", 4, 4)\n", "\n", "# Create G matrix as a BlockMatrix (kept original)\n", - "G = BlockMatrix([\n", - " [G_AA, G_AB, G_AR],\n", - " [G_BA, G_BB, G_BR],\n", - " [G_RA, G_RB, G_RR]\n", - "])\n", + "G = BlockMatrix([[G_AA, G_AB, G_AR], [G_BA, G_BB, G_BR], [G_RA, G_RB, G_RR]])\n", "\n", "# Calculate the product and trace\n", "product = block_collapse(VA @ G @ VB @ G)\n", "result = trace(product)\n", - "print(result)\n" + "print(result)" ] }, { @@ -152,39 +151,39 @@ "from sympy import MatrixSymbol, BlockMatrix, ZeroMatrix, trace, block_collapse\n", "\n", "# Define MatrixSymbols with specific sizes\n", - "AA = MatrixSymbol('AA', 2, 2) # 2x2\n", - "BB = MatrixSymbol('BB', 3, 3) # 3x3\n", + "AA = MatrixSymbol(\"AA\", 2, 2) # 2x2\n", + "BB = MatrixSymbol(\"BB\", 3, 3) # 3x3\n", "\n", "# Create block matrices VA and VB using BlockMatrix with only AA and BB\n", - "VA = BlockMatrix([\n", - " [AA, ZeroMatrix(2,3), ZeroMatrix(2,4)],\n", - " [ZeroMatrix(3,2), ZeroMatrix(3,3), ZeroMatrix(3,4)],\n", - " [ZeroMatrix(4,2), ZeroMatrix(4,3), ZeroMatrix(4,4)]\n", - "])\n", - "\n", - "VB = BlockMatrix([\n", - " [ZeroMatrix(2,2), ZeroMatrix(2,3), ZeroMatrix(2,4)],\n", - " [ZeroMatrix(3,2), BB, ZeroMatrix(3,4)],\n", - " [ZeroMatrix(4,2), ZeroMatrix(4,3), ZeroMatrix(4,4)]\n", - "])\n", + "VA = BlockMatrix(\n", + " [\n", + " [AA, ZeroMatrix(2, 3), ZeroMatrix(2, 4)],\n", + " [ZeroMatrix(3, 2), ZeroMatrix(3, 3), ZeroMatrix(3, 4)],\n", + " [ZeroMatrix(4, 2), ZeroMatrix(4, 3), ZeroMatrix(4, 4)],\n", + " ]\n", + ")\n", + "\n", + "VB = BlockMatrix(\n", + " [\n", + " [ZeroMatrix(2, 2), ZeroMatrix(2, 3), ZeroMatrix(2, 4)],\n", + " [ZeroMatrix(3, 2), BB, ZeroMatrix(3, 4)],\n", + " [ZeroMatrix(4, 2), ZeroMatrix(4, 3), ZeroMatrix(4, 4)],\n", + " ]\n", + ")\n", "\n", "# Define G matrix symbols with matching dimensions\n", - "G_AA = MatrixSymbol('G_AA', 2, 2)\n", - "G_AB = MatrixSymbol('G_AB', 2, 3)\n", - "G_BA = MatrixSymbol('G_BA', 3, 2)\n", - "G_AR = MatrixSymbol('G_AR', 2, 4)\n", - "G_RA = MatrixSymbol('G_RA', 4, 2)\n", - "G_BB = MatrixSymbol('G_BB', 3, 3)\n", - "G_BR = MatrixSymbol('G_BR', 3, 4)\n", - "G_RB = MatrixSymbol('G_RB', 4, 3)\n", - "G_RR = MatrixSymbol('G_RR', 4, 4)\n", + "G_AA = MatrixSymbol(\"G_AA\", 2, 2)\n", + "G_AB = MatrixSymbol(\"G_AB\", 2, 3)\n", + "G_BA = MatrixSymbol(\"G_BA\", 3, 2)\n", + "G_AR = MatrixSymbol(\"G_AR\", 2, 4)\n", + "G_RA = MatrixSymbol(\"G_RA\", 4, 2)\n", + "G_BB = MatrixSymbol(\"G_BB\", 3, 3)\n", + "G_BR = MatrixSymbol(\"G_BR\", 3, 4)\n", + "G_RB = MatrixSymbol(\"G_RB\", 4, 3)\n", + "G_RR = MatrixSymbol(\"G_RR\", 4, 4)\n", "\n", "# Create G matrix as a BlockMatrix\n", - "G = BlockMatrix([\n", - " [G_AA, G_AB, G_AR],\n", - " [G_BA, G_BB, G_BR],\n", - " [G_RA, G_RB, G_RR]\n", - "])\n", + "G = BlockMatrix([[G_AA, G_AB, G_AR], [G_BA, G_BB, G_BR], [G_RA, G_RB, G_RR]])\n", "\n", "# Calculate the product and simplify\n", "product = block_collapse(VA @ G @ VB @ G)\n", diff --git a/src/grogu_magn/grogu.py b/src/grogu_magn/grogu.py index 5fd7b30..03ee14b 100644 --- a/src/grogu_magn/grogu.py +++ b/src/grogu_magn/grogu.py @@ -7,11 +7,14 @@ import sisl from mpi4py import MPI from numpy.linalg import inv from tqdm import tqdm -from useful import * + +from src.grogu_magn.useful import * def main(): - start_time = timer() + start_time = timer() # runtime information + times = dict() + times["start_time"] = timer() # this cell mimicks an input file fdf = sisl.get_sile("./lat3_791/Fe3GeTe2.fdf") @@ -28,45 +31,6 @@ def main(): ] # human readable definition of magnetic entities - # magnetic_entities = [ - # dict(atom=0, ), - # dict(atom=1, ), - # dict(atom=2, ), - # dict(atom=3, l=2), - # dict(atom=4, l=2), - # dict(atom=5, l=2), - # ] - # pairs = [ - # dict(ai=3, aj=4, Ruc=np.array([0, 0, 0])), # isotropic should be -82 meV - # dict(ai=3, aj=5, Ruc=np.array([0, 0, 0])), # these should all be around -41.9 in the isotropic part - # dict(ai=4, aj=5, Ruc=np.array([0, 0, 0])), - # dict(ai=3, aj=0, Ruc=np.array([0, 0, 0])), - # dict(ai=3, aj=1, Ruc=np.array([0, 0, 0])), - # dict(ai=3, aj=2, Ruc=np.array([0, 0, 0])), - # ] - # magnetic_entities = [ - # dict(atom=3, l=2), - # dict(atom=4, l=2), - # dict(atom=5, l=2), - # ] - - # pair information - # pairs = [ - # dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])), # isotropic should be -82 meV - # dict( - # ai=0, aj=2, Ruc=np.array([0, 0, 0]) - # ), # these should all be around -41.9 in the isotropic part - # dict(ai=1, aj=2, Ruc=np.array([0, 0, 0])), - # dict(ai=0, aj=1, Ruc=np.array([-1, 0, 0])), - # dict(ai=0, aj=2, Ruc=np.array([-1, 0, 0])), - # dict(ai=0, aj=1, Ruc=np.array([1, 0, 0])), - # dict(ai=0, aj=2, Ruc=np.array([1, 0, 0])), - # dict(ai=0, aj=1, Ruc=np.array([0, -1, 0])), - # dict(ai=0, aj=2, Ruc=np.array([0, -1, 0])), - # dict(ai=0, aj=1, Ruc=np.array([0, 1, 0])), - # dict(ai=0, aj=2, Ruc=np.array([0, 1, 0])), - # dict(ai=1, aj=2, Ruc=np.array([-1, 0, 0])), - # ] magnetic_entities = [ dict(atom=3, l=2), dict(atom=4, l=2), @@ -77,11 +41,13 @@ def main(): ] # pair information + # these should all be around -41.9 in the isotropic part + # isotropic should be -82 meV pairs = [ - dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])), # isotropic should be -82 meV - dict( - ai=0, aj=2, Ruc=np.array([0, 0, 0]) - ), # these should all be around -41.9 in the isotropic part + dict(ai=0, aj=3, Ruc=np.array([0, 0, 0])), + dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])), + dict(ai=1, aj=0, Ruc=np.array([0, 0, 0])), + dict(ai=0, aj=2, Ruc=np.array([0, 0, 0])), dict(ai=1, aj=2, Ruc=np.array([0, 0, 0])), dict(ai=0, aj=2, Ruc=np.array([-1, 0, 0])), dict(ai=1, aj=2, Ruc=np.array([-1, 0, 0])), @@ -92,7 +58,7 @@ def main(): kdirs = "xy" ebot = -30 eset = 50 - esetp = 1000 + esetp = 10000 # MPI parameters comm = MPI.COMM_WORLD @@ -125,7 +91,7 @@ def main(): # unit cell index uc_in_sc_idx = dh.lattice.sc_index([0, 0, 0]) - setup_time = timer() + times["setup_time"] = timer() NO = dh.no # shorthand for number of orbitals in the unit cell @@ -213,18 +179,16 @@ def main(): f"Exchange field has non negligible scalar part. Largest value is {max_xcfs}" ) - H_and_XCF_time = timer() + times["H_and_XCF_time"] = timer() # for every site we have to store 3 Greens function (and the associated _tmp-s) in the 3 reference directions for i, mag_ent in enumerate(magnetic_entities): parsed = parse_magnetic_entity(dh, **mag_ent) # parse orbital indexes magnetic_entities[i]["orbital_indeces"] = parsed - magnetic_entities[i]["spin_box_indeces"] = blow_up_orbindx( - parsed - ) # calculate spin box indexes - spin_box_shape = len( - mag_ent["spin_box_indeces"] - ) # calculate size for Greens function generation + # calculate spin box indexes + magnetic_entities[i]["spin_box_indeces"] = blow_up_orbindx(parsed) + # calculate size for Greens function generation + spin_box_shape = len(mag_ent["spin_box_indeces"]) mag_ent["energies"] = ( [] @@ -232,14 +196,14 @@ def main(): mag_ent["Gii"] = [] # Greens function mag_ent["Gii_tmp"] = [] # Greens function for parallelization - mag_ent["Vu1"] = [ - list([]) for _ in range(len(ref_xcf_orientations)) - ] # These will be the perturbed potentials from eq. 100 + # These will be the perturbed potentials from eq. 100 + mag_ent["Vu1"] = [list([]) for _ in range(len(ref_xcf_orientations))] mag_ent["Vu2"] = [list([]) for _ in range(len(ref_xcf_orientations))] for i in ref_xcf_orientations: + # Greens functions for every quantization axis mag_ent["Gii"].append( np.zeros((eset, spin_box_shape, spin_box_shape), dtype="complex128") - ) # Greens functions for every quantization axis + ) mag_ent["Gii_tmp"].append( np.zeros((eset, spin_box_shape, spin_box_shape), dtype="complex128") ) @@ -247,11 +211,9 @@ def main(): # 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: - spin_box_shape_i, spin_box_shape_j = len( - magnetic_entities[pair["ai"]]["spin_box_indeces"] - ), len( - magnetic_entities[pair["aj"]]["spin_box_indeces"] - ) # calculate size for Greens function generation + # calculate size for Greens function generation + spin_box_shape_i = len(magnetic_entities[pair["ai"]]["spin_box_indeces"]) + spin_box_shape_j = len(magnetic_entities[pair["aj"]]["spin_box_indeces"]) pair["energies"] = [] # we will store the second order energy derivations here @@ -273,16 +235,16 @@ def main(): np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype="complex128") ) - site_and_pair_dictionaries_time = timer() + times["site_and_pair_dictionaries_time"] = timer() kset = make_kset(dirs=kdirs, NUMK=kset) # generate k space sampling wkset = np.ones(len(kset)) / len(kset) # generate weights for k points kpcs = np.array_split(kset, size) # split the k points based on MPI size kpcs[root_node] = tqdm(kpcs[root_node], desc="k loop", file=stdout) - k_set_time = timer() + times["k_set_time"] = timer() - # this will contain all the data needed to calculate the energy variations upon rotation + # this will contain the three hamiltonians in the reference directions needed to calculate the energy variations upon rotation hamiltonians = [] # iterate over the reference directions (quantization axes) @@ -296,15 +258,16 @@ def main(): rot_H_XCF_uc = rot_H_XCF[uc_in_sc_idx] # obtain total Hamiltonian with the rotated exchange field - rot_H = hTRS + rot_H_XCF # equation 76 + rot_H = ( + hTRS + rot_H_XCF + ) # equation 76 ####################################################################################### hamiltonians.append( - dict(orient=orient["o"], H=rot_H, rotations=[]) + dict(orient=orient["o"], H=rot_H) ) # store orientation and rotated Hamiltonian - for u in orient[ - "vw" - ]: # these are the infinitezimal rotations (for now) perpendicular to the quantization axis + # these are the infinitezimal rotations (for now) perpendicular to the quantization axis + for u in orient["vw"]: Tu = np.kron(np.eye(NO, dtype=int), tau_u(u)) # section 2.H Vu1 = 1j / 2 * commutator(rot_H_XCF_uc, Tu) # equation 100 @@ -318,7 +281,7 @@ def main(): Vu2[:, mag_ent["spin_box_indeces"]][mag_ent["spin_box_indeces"], :] ) - reference_rotations_time = timer() + times["reference_rotations_time"] = timer() if rank == root_node: print("Number of magnetic entities being calculated: ", len(magnetic_entities)) @@ -339,9 +302,8 @@ def main(): # sampling the integrand on the contour and the BZ for k in kpcs[rank]: wk = wkset[rank] # weight of k point in BZ integral - for i, hamiltonian_orientation in enumerate( - hamiltonians - ): # iterate over reference directions + # iterate over reference directions + for i, hamiltonian_orientation in enumerate(hamiltonians): # calculate Greens function H = hamiltonian_orientation["H"] HK, SK = hsk(H, ss, dh.sc_off, k) @@ -375,7 +337,7 @@ def main(): comm.Reduce(pair["Gij_tmp"][i], pair["Gij"][i], root=root_node) comm.Reduce(pair["Gji_tmp"][i], pair["Gji"][i], root=root_node) - green_function_inversion_time = timer() + times["green_function_inversion_time"] = timer() if rank == root_node: # iterate over the magnetic entities @@ -393,7 +355,7 @@ def main(): storage.append(np.trapz(-1 / np.pi * np.imag(traced * cont.we))) # fill up the magnetic entities dictionary with the energies - mag_ent["energies"].append(storage) + magnetic_entities[tracker]["energies"].append(storage) # iterate over the pairs for tracker, pair in enumerate(pairs): @@ -410,109 +372,11 @@ def main(): traced = np.trace((Vui @ Gij @ Vuj @ Gji), axis1=1, axis2=2) # evaluation of the contour integral storage.append(np.trapz(-1 / np.pi * np.imag(traced * cont.we))) - # fill up the pairs dictionary with the energies pairs[tracker]["energies"].append(storage) - end_time = timer() - - print( - "############################### GROGU OUTPUT ###################################" - ) - print( - "================================================================================" - ) - print("Input file: ") - print(simulation_parameters["path"]) - print( - "Number of nodes in the parallel cluster: ", - simulation_parameters["parallel_size"], - ) - print( - "================================================================================" - ) - try: - print("Cell [Ang]: ") - print(simulation_parameters["geom"].cell) - except: - print("Geometry could not be read.") - print( - "================================================================================" - ) - print("DFT axis: ") - print(simulation_parameters["scf_xcf_orientation"]) - print("Quantization axis and perpendicular rotation directions:") - for ref in ref_xcf_orientations: - print(ref["o"], " --» ", ref["vw"]) - print( - "================================================================================" - ) - print("number of k points: ", simulation_parameters["kset"]) - print("k point directions: ", simulation_parameters["kdirs"]) - print( - "================================================================================" - ) - print("Parameters for the contour integral:") - print("Ebot: ", simulation_parameters["ebot"]) - print("Eset: ", simulation_parameters["eset"]) - print("Esetp: ", simulation_parameters["esetp"]) - print( - "================================================================================" - ) - print("Atomic informations: ") - print("") - print("") - print("Not yet specified.") - print("") - print("") - print( - "================================================================================" - ) - print("Exchange [meV]") - print( - "--------------------------------------------------------------------------------" - ) - print("Atom1 Atom2 [i j k] d [Ang]") - print( - "--------------------------------------------------------------------------------" - ) - for pair in pairs: - J_iso, J_S, D = calculate_exchange_tensor(pair) - J_iso = J_iso * sisl.unit_convert("eV", "meV") - J_S = J_S * sisl.unit_convert("eV", "meV") - D = D * sisl.unit_convert("eV", "meV") - - print(print_atomic_indices(pair, magnetic_entities, dh)) - print("Isotropic: ", J_iso) - print("DMI: ", D) - print("Symmetric-anisotropy: ", J_S) - print("") - - print( - "================================================================================" - ) - print("Runtime information: ") - print("Total runtime: ", end_time - start_time) - print( - "--------------------------------------------------------------------------------" - ) - print("Initial setup: ", setup_time - start_time) - print( - f"Hamiltonian conversion and XC field extraction: {H_and_XCF_time - setup_time:.3f} s" - ) - print( - f"Pair and site datastructure creatrions: {site_and_pair_dictionaries_time - H_and_XCF_time:.3f} s" - ) - print( - f"k set cration and distribution: {k_set_time - site_and_pair_dictionaries_time:.3f} s" - ) - print(f"Rotating XC potential: {reference_rotations_time - k_set_time:.3f} s") - print( - f"Greens function inversion: {green_function_inversion_time - reference_rotations_time:.3f} s" - ) - print( - f"Calculate energies and magnetic components: {end_time - green_function_inversion_time:.3f} s" - ) + times["end_time"] = timer() + print_output(simulation_parameters, magnetic_entities, pairs, dh, times) if __name__ == "__main__": diff --git a/src/grogu_magn/useful.py b/src/grogu_magn/useful.py index 90a33cb..55adb34 100644 --- a/src/grogu_magn/useful.py +++ b/src/grogu_magn/useful.py @@ -22,6 +22,7 @@ from itertools import permutations, product import numpy as np from scipy.special import roots_legendre +from sisl import unit_convert # Pauli matrices tau_x = np.array([[0, 1], [1, 0]]) @@ -226,33 +227,170 @@ def calculate_exchange_tensor(pair): return J_ii.sum() / 3, np.concatenate([J_ii[:2] - J_ii.sum() / 3, J_S]).flatten(), D -def print_atomic_indices(pair, magnetic_entities, dh): +def print_pair_atomic_indices(pair, magnetic_entities, dh): atomic_indices = "" - atoms = magnetic_entities[pair["ai"]] - if "l" not in atoms.keys(): - atoms["l"] = "all" - if isinstance(atoms["atom"], int): - atomic_indices += ( - f"[{atoms['atom']}]{dh.atoms[atoms['atom']].tag}({atoms['l']})" - ) - if isinstance(atoms, list): - atomic_indices += [ - f"[{atoms['atom']}]{dh.atoms[atom['atom']].tag}({atom['l']})" - for atom in atoms["atom"] - ] - atoms = magnetic_entities[pair["aj"]] - if "l" not in atoms.keys(): - atoms["l"] = "all" - atomic_indices += " " - if isinstance(atoms["atom"], int): - atomic_indices += ( - f"[{atoms['atom']}]{dh.atoms[atoms['atom']].tag}({atoms['l']})" - ) - if isinstance(atoms, list): - atomic_indices += [ - f"[{atoms['atom']}]{dh.atoms[atom['atom']].tag}({atom['l']})" - for atom in atoms["atom"] - ] - - atomic_indices += f" {pair['Ruc']} d [Ang] Not yet." + # iterate over the two magnetic entities in a pair + for mag_ent in [magnetic_entities[pair["ai"]], magnetic_entities[pair["aj"]]]: + # get atoms of magnetic entity + atoms_idx = mag_ent["atom"] + # if orbital is not set use all + if "l" not in mag_ent.keys(): + mag_ent["l"] = "all" + orbitals = mag_ent["l"] + + # if magnetic entity contains one atoms + if isinstance(atoms_idx, int): + atomic_indices += 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 + atomic_indices += atom_group[:-2] + "}" + # separate magnetic entities + atomic_indices += " " return atomic_indices + + +def print_output(simulation_parameters, magnetic_entities, pairs, dh, times): + print( + "##################################################################### GROGU OUTPUT #############################################################################" + ) + print( + "================================================================================================================================================================" + ) + print("Input file: ") + print(simulation_parameters["path"]) + print( + "Number of nodes in the parallel cluster: ", + simulation_parameters["parallel_size"], + ) + print( + "================================================================================================================================================================" + ) + try: + print("Cell [Ang]: ") + print(simulation_parameters["geom"].cell) + except: + print("Geometry could not be read.") + print( + "================================================================================================================================================================" + ) + print("DFT axis: ") + print(simulation_parameters["scf_xcf_orientation"]) + print("Quantization axis and perpendicular rotation directions:") + for ref in simulation_parameters["ref_xcf_orientations"]: + print(ref["o"], " --» ", ref["vw"]) + print( + "================================================================================================================================================================" + ) + print("number of k points: ", simulation_parameters["kset"]) + print("k point directions: ", simulation_parameters["kdirs"]) + print( + "================================================================================================================================================================" + ) + print("Parameters for the contour integral:") + print("Ebot: ", simulation_parameters["ebot"]) + print("Eset: ", simulation_parameters["eset"]) + print("Esetp: ", simulation_parameters["esetp"]) + print( + "================================================================================================================================================================" + ) + print("Atomic informations: ") + print( + "----------------------------------------------------------------------------------------------------------------------------------------------------------------" + ) + print( + "[atom index]Element(orbitals) x [Ang] y [Ang] z [Ang] Sx Sy Sz Q Lx Ly Lz Jx Jy Jz" + ) + print( + "----------------------------------------------------------------------------------------------------------------------------------------------------------------" + ) + # iterate over magnetic entities + for mag_ent in magnetic_entities: + # get atoms of magnetic entity + atoms_idx = mag_ent["atom"] + # if orbital is not set use all + if "l" not in mag_ent.keys(): + mag_ent["l"] = "all" + orbitals = mag_ent["l"] + + # if magnetic entity contains one atom + if isinstance(atoms_idx, int): + # coordinates and tag + x, y, z = dh.xyz[atoms_idx] + print( + f"[{atoms_idx}]{dh.atoms[atoms_idx].tag}({orbitals}) {x} {y} {z}" + ) + + # if magnetic entity contains more than one atoms + if isinstance(atoms_idx, list): + # iterate over atoms + for atom_idx in atoms_idx: + # coordinates and tag + x, y, z = dh.xyz[atom_idx] + print( + f"[{atom_idx}]{dh.atoms[atom_idx].tag}({orbitals}) {x} {y} {z}" + ) + print("") + + print( + "================================================================================================================================================================" + ) + print("Exchange [meV]") + print( + "----------------------------------------------------------------------------------------------------------------------------------------------------------------" + ) + print("Magnetic entity1 Magnetic entity2 [i j k] d [Ang]") + print( + "----------------------------------------------------------------------------------------------------------------------------------------------------------------" + ) + # iterate over pairs + for pair in pairs: + # calculate magnetic parameters + J_iso, J_S, D = calculate_exchange_tensor(pair) + J_iso = J_iso * unit_convert("eV", "meV") + J_S = J_S * unit_convert("eV", "meV") + D = D * unit_convert("eV", "meV") + + # print pair parameters + print( + print_pair_atomic_indices(pair, magnetic_entities, dh), + f" {pair['Ruc']} d [Ang] Not yet.", + ) + # print magnetic parameters + print("Isotropic: ", J_iso) + print("DMI: ", D) + print("Symmetric-anisotropy: ", J_S) + print("") + + print( + "================================================================================================================================================================" + ) + print("Runtime information: ") + print(f"Total runtime: {times['end_time'] - times['start_time']} s") + print( + "----------------------------------------------------------------------------------------------------------------------------------------------------------------" + ) + print(f"Initial setup: {times['setup_time'] - times['start_time']} s") + print( + f"Hamiltonian conversion and XC field extraction: {times['H_and_XCF_time'] - times['setup_time']:.3f} s" + ) + print( + f"Pair and site datastructure creatrions: {times['site_and_pair_dictionaries_time'] - times['H_and_XCF_time']:.3f} s" + ) + print( + f"k set cration and distribution: {times['k_set_time'] - times['site_and_pair_dictionaries_time']:.3f} s" + ) + print( + f"Rotating XC potential: {times['reference_rotations_time'] - times['k_set_time']:.3f} s" + ) + print( + f"Greens function inversion: {times['green_function_inversion_time'] - times['reference_rotations_time']:.3f} s" + ) + print( + f"Calculate energies and magnetic components: {times['end_time'] - times['green_function_inversion_time']:.3f} s" + ) diff --git a/test.ipynb b/test.ipynb index dfcfacd..c2e63ce 100644 --- a/test.ipynb +++ b/test.ipynb @@ -9,7 +9,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "[Mac:47319] shmem: mmap: an error occurred while determining whether or not /var/folders/yh/dx7xl94n3g52ts3td8qcxjcc0000gn/T//ompi.Mac.501/jf.0/1258881024/sm_segment.Mac.501.4b090000.0 could be created.\n" + "[Mac:54919] shmem: mmap: an error occurred while determining whether or not /var/folders/yh/dx7xl94n3g52ts3td8qcxjcc0000gn/T//ompi.Mac.501/jf.0/778436608/sm_segment.Mac.501.2e660000.0 could be created.\n" ] } ], @@ -32,7 +32,9 @@ "from numpy.linalg import inv\n", "import warnings\n", "\n", - "start_time = timer()" + "# runtime information\n", + "times = dict()\n", + "times[\"start_time\"] = timer()" ] }, { @@ -59,25 +61,28 @@ "# we can have some default for this\n", "ref_xcf_orientations = [\n", " dict(o=np.array([1, 0, 0]), vw=[np.array([0, 1, 0]), np.array([0, 0, 1])]),\n", - " dict(o=np.array([0, 1, 0]), vw=[np.array([1, 0, 1]), np.array([0, 0, 1])]),\n", + " dict(o=np.array([0, 1, 0]), vw=[np.array([1, 0, 0]), np.array([0, 0, 1])]),\n", " dict(o=np.array([0, 0, 1]), vw=[np.array([1, 0, 0]), np.array([0, 1, 0])]),\n", "]\n", "\n", "# human readable definition of magnetic entities\n", - "magnetic_entities=[dict(atom=3 ,l=2),\n", - " dict(atom=4 ,l=2),\n", - " dict(atom=5 ,l=2),\n", - " dict(atom=[3,4],)]\n", + "magnetic_entities = [\n", + " dict(atom=3, l=2),\n", + " dict(atom=4, l=2),\n", + " dict(atom=5, l=2),\n", + " dict(\n", + " atom=[3, 4],\n", + " ),\n", + "]\n", "\n", "# pair information\n", "# these should all be around -41.9 in the isotropic part\n", "# isotropic should be -82 meV\n", "pairs = [\n", + " dict(ai=0, aj=3, Ruc=np.array([0, 0, 0])),\n", " dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])),\n", " dict(ai=1, aj=0, Ruc=np.array([0, 0, 0])),\n", - " dict(\n", - " ai=0, aj=2, Ruc=np.array([0, 0, 0])\n", - " ),\n", + " dict(ai=0, aj=2, Ruc=np.array([0, 0, 0])),\n", " dict(ai=1, aj=2, Ruc=np.array([0, 0, 0])),\n", " dict(ai=0, aj=2, Ruc=np.array([-1, 0, 0])),\n", " dict(ai=1, aj=2, Ruc=np.array([-1, 0, 0])),\n", @@ -122,7 +127,7 @@ "# unit cell index\n", "uc_in_sc_idx = dh.lattice.sc_index([0, 0, 0])\n", "\n", - "setup_time = timer()" + "times[\"setup_time\"] = timer()" ] }, { @@ -215,7 +220,7 @@ " f\"Exchange field has non negligible scalar part. Largest value is {max_xcfs}\"\n", " )\n", "\n", - "H_and_XCF_time = timer()" + "times[\"H_and_XCF_time\"] = timer()" ] }, { @@ -263,12 +268,13 @@ " pair[\"Gij_tmp\"] = [] # Greens function for parallelization\n", " pair[\"Gji_tmp\"] = []\n", " for i in ref_xcf_orientations:\n", + " # Greens functions for every quantization axis\n", " pair[\"Gij\"].append(\n", " np.zeros((eset, spin_box_shape_i, spin_box_shape_j), dtype=\"complex128\")\n", " )\n", " pair[\"Gij_tmp\"].append(\n", " np.zeros((eset, spin_box_shape_i, spin_box_shape_j), dtype=\"complex128\")\n", - " ) # Greens functions for every quantization axis\n", + " )\n", " pair[\"Gji\"].append(\n", " np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype=\"complex128\")\n", " )\n", @@ -276,7 +282,7 @@ " np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype=\"complex128\")\n", " )\n", "\n", - "site_and_pair_dictionaries_time = timer()" + "times[\"site_and_pair_dictionaries_time\"] = timer()" ] }, { @@ -298,7 +304,7 @@ "kpcs = np.array_split(kset, size) # split the k points based on MPI size\n", "kpcs[root_node] = tqdm(kpcs[root_node], desc=\"k loop\", file=stdout)\n", "\n", - "k_set_time = timer()" + "times[\"k_set_time\"] = timer()" ] }, { @@ -337,14 +343,15 @@ " Vu2 = 1 / 8 * commutator(commutator(Tu, rot_H_XCF_uc), Tu) # equation 100\n", "\n", " for mag_ent in magnetic_entities:\n", + " # fill up the perturbed potentials (for now) based on the on-site projections\n", " mag_ent[\"Vu1\"][i].append(\n", " Vu1[:, mag_ent[\"spin_box_indeces\"]][mag_ent[\"spin_box_indeces\"], :]\n", - " ) # fill up the perturbed potentials (for now) based on the on-site projections\n", + " )\n", " mag_ent[\"Vu2\"][i].append(\n", " Vu2[:, mag_ent[\"spin_box_indeces\"]][mag_ent[\"spin_box_indeces\"], :]\n", " )\n", "\n", - "reference_rotations_time = timer()" + "times[\"reference_rotations_time\"] = timer()" ] }, { @@ -359,7 +366,7 @@ "Number of magnetic entities being calculated: 4\n", "We have to calculate the Greens function for three reference direction and we are going to calculate 15 energy integrals per site.\n", "The shape of the Hamiltonian and the Greens function is 84x84.\n", - "k loop: 100%|██████████| 400/400 [01:50<00:00, 3.61it/s]\n" + "k loop: 100%|██████████| 400/400 [19:35<00:00, 2.94s/it] \n" ] } ], @@ -418,7 +425,7 @@ " comm.Reduce(pair[\"Gij_tmp\"][i], pair[\"Gij\"][i], root=root_node)\n", " comm.Reduce(pair[\"Gji_tmp\"][i], pair[\"Gji\"][i], root=root_node)\n", "\n", - "green_function_inversion_time = timer()" + "times[\"green_function_inversion_time\"] = timer()" ] }, { @@ -430,86 +437,99 @@ "name": "stdout", "output_type": "stream", "text": [ - "############################### GROGU OUTPUT ###################################\n", - "================================================================================\n", + "##################################################################### GROGU OUTPUT #############################################################################\n", + "================================================================================================================================================================\n", "Input file: \n", "Not yet specified.\n", "Number of nodes in the parallel cluster: 1\n", - "================================================================================\n", + "================================================================================================================================================================\n", "Cell [Ang]: \n", "[[ 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", - "================================================================================\n", + "================================================================================================================================================================\n", "DFT axis: \n", "[0 0 1]\n", "Quantization axis and perpendicular rotation directions:\n", "[1 0 0] --» [array([0, 1, 0]), array([0, 0, 1])]\n", - "[0 1 0] --» [array([1, 0, 1]), array([0, 0, 1])]\n", + "[0 1 0] --» [array([1, 0, 0]), array([0, 0, 1])]\n", "[0 0 1] --» [array([1, 0, 0]), array([0, 1, 0])]\n", - "================================================================================\n", + "================================================================================================================================================================\n", "number of k points: 20\n", "k point directions: xy\n", - "================================================================================\n", + "================================================================================================================================================================\n", "Parameters for the contour integral:\n", "Ebot: -30\n", "Eset: 50\n", "Esetp: 10000\n", - "================================================================================\n", + "================================================================================================================================================================\n", "Atomic informations: \n", + "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n", + "[atom index]Element(orbitals) x [Ang] y [Ang] z [Ang] Sx Sy Sz Q Lx Ly Lz Jx Jy Jz\n", + "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n", + "[3]Fe(2) -7.339158738013707e-06 4.149278510690423e-06 11.657585837928032\n", "\n", + "[4]Fe(2) -7.326987662162937e-06 4.158274523275774e-06 8.912422537596708\n", "\n", - "Not yet specified.\n", + "[5]Fe(2) 1.8954667088117545 1.0943913231921656 10.285002698393109\n", "\n", + "[3]Fe(all) -7.339158738013707e-06 4.149278510690423e-06 11.657585837928032\n", + "[4]Fe(all) -7.326987662162937e-06 4.158274523275774e-06 8.912422537596708\n", "\n", - "================================================================================\n", + "================================================================================================================================================================\n", "Exchange [meV]\n", - "--------------------------------------------------------------------------------\n", - "Atom1 Atom2 [i j k] d [Ang]\n", - "--------------------------------------------------------------------------------\n", - "[3]Fe(2) [4]Fe(2) [0 0 0] d [Ang] Not yet.\n", - "Isotropic: -61.134994398006846\n", - "DMI: [-6.44078586e+00 7.51906171e+00 6.90431275e-04 9.74101032e-04\n", + "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n", + "Magnetic entity1 Magnetic entity2 [i j k] d [Ang]\n", + "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n", + "[3]Fe(2) {[3]Fe(all)--[4]Fe(all)} [0 0 0] d [Ang] Not yet.\n", + "Isotropic: -7404.692892728154\n", + "DMI: [-9.31336133e-01 -7.86063119e-04 1.83575526e-06]\n", + "Symmetric-anisotropy: [ 2.03042884e+00 1.04458862e+00 -1.40962466e-03 -6.58426045e-04\n", + " -6.97969537e-02]\n", + "\n", + "[3]Fe(2) [4]Fe(2) [0 0 0] d [Ang] Not yet.\n", + "Isotropic: -60.89330922309355\n", + "DMI: [-9.32966923e-01 -8.92579299e-04 -2.04258659e-06]\n", + "Symmetric-anisotropy: [-5.95741551e+00 7.27737654e+00 6.90431275e-04 -8.11057566e-04\n", " -5.49031203e-06]\n", - "Symmetric-anisotropy: [-9.32966923e-01 8.92579299e-04 -2.04258659e-06]\n", "\n", - "[4]Fe(2) [3]Fe(2) [0 0 0] d [Ang] Not yet.\n", - "Isotropic: -61.134994398006846\n", - "DMI: [-6.44078586e+00 7.51906171e+00 6.90431275e-04 -8.11057566e-04\n", + "[4]Fe(2) [3]Fe(2) [0 0 0] d [Ang] Not yet.\n", + "Isotropic: -60.893309223093574\n", + "DMI: [9.32966923e-01 8.92579299e-04 2.04258658e-06]\n", + "Symmetric-anisotropy: [-5.95741551e+00 7.27737654e+00 6.90431275e-04 9.74101032e-04\n", " -5.49031203e-06]\n", - "Symmetric-anisotropy: [ 9.32966923e-01 -8.92579299e-04 2.04258658e-06]\n", "\n", - "[3]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n", - "Isotropic: -60.47084255793878\n", - "DMI: [-0.17431857 0.57352788 0.07106945 6.04526746 -0.0424978 ]\n", - "Symmetric-anisotropy: [3.78506176e+00 6.13838308e+00 3.59037036e-03]\n", + "[3]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n", + "Isotropic: -60.55651225519789\n", + "DMI: [ 3.78506176e+00 -6.13838308e+00 3.59037036e-03]\n", + "Symmetric-anisotropy: [-0.34565796 0.65919757 0.07106945 -6.23149871 -0.0424978 ]\n", "\n", - "[4]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n", - "Isotropic: -60.46575240352053\n", - "DMI: [-0.16838554 0.56839367 0.07106826 -6.06303351 0.03636701]\n", - "Symmetric-anisotropy: [-3.79945963e+00 -6.15244494e+00 3.58990840e-03]\n", + "[4]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n", + "Isotropic: -60.54974989595536\n", + "DMI: [-3.79945963e+00 6.15244494e+00 3.58990840e-03]\n", + "Symmetric-anisotropy: [-0.33638052 0.65239116 0.07106826 6.24185638 0.03636701]\n", "\n", - "[3]Fe(2) [5]Fe(2) [-1 0 0] d [Ang] Not yet.\n", - "Isotropic: -6.691813632755354\n", - "DMI: [-0.79119686 0.79044933 -0.031302 -7.59758103 0.03239586]\n", - "Symmetric-anisotropy: [ 5.95251705 -7.64859703 6.50501652]\n", + "[3]Fe(2) [5]Fe(2) [-1 0 0] d [Ang] Not yet.\n", + "Isotropic: -6.6253295899433216\n", + "DMI: [5.95251705 7.64859703 6.50501652]\n", + "Symmetric-anisotropy: [-0.65822877 0.72396528 -0.031302 7.69961304 0.03239586]\n", "\n", - "[4]Fe(2) [5]Fe(2) [-1 0 0] d [Ang] Not yet.\n", - "Isotropic: -6.190626204103167\n", - "DMI: [ 0.19185206 0.28941842 -0.03129943 -4.21429122 -0.09833472]\n", - "Symmetric-anisotropy: [ 6.19414647 -4.23019689 6.50504332]\n", + "[4]Fe(2) [5]Fe(2) [-1 0 0] d [Ang] Not yet.\n", + "Isotropic: -6.123086646494101\n", + "DMI: [6.19414647 4.23019689 6.50504332]\n", + "Symmetric-anisotropy: [ 0.32693117 0.22187887 -0.03129943 4.24610256 -0.09833472]\n", "\n", - "================================================================================\n", + "================================================================================================================================================================\n", "Runtime information: \n", - "Total runtime: 111.725750458\n", - "--------------------------------------------------------------------------------\n", - "Initial setup: 0.20263604099999988\n", - "Hamiltonian conversion and XC field extraction: 0.595 s\n", + "Total runtime: 118.107833833 s\n", + "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n", + "Initial setup: 0.12060429199999945 s\n", + "Hamiltonian conversion and XC field extraction: 0.558 s\n", "Pair and site datastructure creatrions: 0.011 s\n", - "k set cration and distribution: 0.019 s\n", - "Rotating XC potential: 0.210 s\n", - "Greens function inversion: 110.600 s\n", - "Calculate energies and magnetic components: 0.087 s\n" + "k set cration and distribution: 0.018 s\n", + "Rotating XC potential: 0.214 s\n", + "Greens function inversion: 117.061 s\n", + "Calculate energies and magnetic components: 0.127 s\n" ] } ], @@ -548,105 +568,8 @@ " # fill up the pairs dictionary with the energies\n", " pairs[tracker][\"energies\"].append(storage)\n", "\n", - " end_time = timer()\n", - "\n", - " print(\n", - " \"############################### GROGU OUTPUT ###################################\"\n", - " )\n", - " print(\n", - " \"================================================================================\"\n", - " )\n", - " print(\"Input file: \")\n", - " print(simulation_parameters[\"path\"])\n", - " print(\n", - " \"Number of nodes in the parallel cluster: \",\n", - " simulation_parameters[\"parallel_size\"],\n", - " )\n", - " print(\n", - " \"================================================================================\"\n", - " )\n", - " try:\n", - " print(\"Cell [Ang]: \")\n", - " print(simulation_parameters[\"geom\"].cell)\n", - " except:\n", - " print(\"Geometry could not be read.\")\n", - " print(\n", - " \"================================================================================\"\n", - " )\n", - " print(\"DFT axis: \")\n", - " print(simulation_parameters[\"scf_xcf_orientation\"])\n", - " print(\"Quantization axis and perpendicular rotation directions:\")\n", - " for ref in ref_xcf_orientations:\n", - " print(ref[\"o\"], \" --» \", ref[\"vw\"])\n", - " print(\n", - " \"================================================================================\"\n", - " )\n", - " print(\"number of k points: \", simulation_parameters[\"kset\"])\n", - " print(\"k point directions: \", simulation_parameters[\"kdirs\"])\n", - " print(\n", - " \"================================================================================\"\n", - " )\n", - " print(\"Parameters for the contour integral:\")\n", - " print(\"Ebot: \", simulation_parameters[\"ebot\"])\n", - " print(\"Eset: \", simulation_parameters[\"eset\"])\n", - " print(\"Esetp: \", simulation_parameters[\"esetp\"])\n", - " print(\n", - " \"================================================================================\"\n", - " )\n", - " print(\"Atomic informations: \")\n", - " print(\"\")\n", - " print(\"\")\n", - " print(\"Not yet specified.\")\n", - " print(\"\")\n", - " print(\"\")\n", - " print(\n", - " \"================================================================================\"\n", - " )\n", - " print(\"Exchange [meV]\")\n", - " print(\n", - " \"--------------------------------------------------------------------------------\"\n", - " )\n", - " print(\"Atom1 Atom2 [i j k] d [Ang]\")\n", - " print(\n", - " \"--------------------------------------------------------------------------------\"\n", - " )\n", - " for pair in pairs:\n", - " J_iso, J_S, D = calculate_exchange_tensor(pair)\n", - " J_iso = J_iso * sisl.unit_convert(\"eV\", \"meV\")\n", - " J_S = J_S * sisl.unit_convert(\"eV\", \"meV\")\n", - " D = D * sisl.unit_convert(\"eV\", \"meV\")\n", - "\n", - " print(print_atomic_indices(pair, magnetic_entities, dh))\n", - " print(\"Isotropic: \", J_iso)\n", - " print(\"DMI: \", D)\n", - " print(\"Symmetric-anisotropy: \", J_S)\n", - " print(\"\")\n", - "\n", - " print(\n", - " \"================================================================================\"\n", - " )\n", - " print(\"Runtime information: \")\n", - " print(\"Total runtime: \", end_time - start_time)\n", - " print(\n", - " \"--------------------------------------------------------------------------------\"\n", - " )\n", - " print(\"Initial setup: \", setup_time - start_time)\n", - " print(\n", - " f\"Hamiltonian conversion and XC field extraction: {H_and_XCF_time - setup_time:.3f} s\"\n", - " )\n", - " print(\n", - " f\"Pair and site datastructure creatrions: {site_and_pair_dictionaries_time - H_and_XCF_time:.3f} s\"\n", - " )\n", - " print(\n", - " f\"k set cration and distribution: {k_set_time - site_and_pair_dictionaries_time:.3f} s\"\n", - " )\n", - " print(f\"Rotating XC potential: {reference_rotations_time - k_set_time:.3f} s\")\n", - " print(\n", - " f\"Greens function inversion: {green_function_inversion_time - reference_rotations_time:.3f} s\"\n", - " )\n", - " print(\n", - " f\"Calculate energies and magnetic components: {end_time - green_function_inversion_time:.3f} s\"\n", - " )" + " times[\"end_time\"] = timer()\n", + " print_output(simulation_parameters, magnetic_entities, pairs, dh, times)" ] }, { @@ -706,9 +629,27 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 27, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "xyz[-3:]: red, green, blue\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAABNoAAAHACAYAAAB0/gUQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCsElEQVR4nO3df5hWdZ0//ucMyAwWM2rK8GMnf5b4E0iTMFu1Rln1YmU/5aLtAktqW9l+VLYfUiqWJqlpbEayaUbq9vG35ie9IGUjM1lNlNZMTRSElEH9qjOICjJzf//g47QToAOcmXtmeDyu61x5n/u87/v1PkznNfdzzn1ORalUKgUAAAAA2CqV5S4AAAAAAHoDQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEAB+pa7gO6otbU1zz//fAYMGJCKiopylwPQ45VKpaxatSpDhgxJZaW/8egzAMXSZzak1wAUq6O9RtC2Ec8//3zq6+vLXQZAr7N8+fL81V/9VbnLKDt9BqBz6DN/ptcAdI536zWCto0YMGBAkvU7r6ampszVAPR8zc3Nqa+vbzu+buv0GYBi6TMb0msAitXRXiNo24i3T62uqanRlAAK5Ksr6+kzAJ1Dn/kzvQagc7xbr3EBAwAAAAAogKANAAAAAAogaAMAAACAAgjaAAAAAKAAgjYAAAAAKICgDQAAAAAKIGgDAAAAgAII2gAAAACgAII2AAAAAChA33IX0KusWJHMnp388Y/JgAHJCSckhx2WVFSUuzIAeoM330xuvDG59971veXww5NPfSqpri53ZQAA0C2tW5fccUcyd+76//7wh5N/+If1sU1nELQV5fLLkzPPTEqlpLLyz+uOOCK5/faktrac1QHQ0/32t8lxxyUvvpj0/X/t+6qrki99KbnrruRDHypvfQAA0M0sXpyMGZM888yff4X+8Y+TL385ueWW5Oiji39PXx0twq23Jv/7fyctLUlr6/qIdN269c/9+tfJ+PHlrQ+Anm3FiuSoo5KXX17/+H/2mZdeShoakhdeKF99AADQzbz+evLxjyfPPrv+8du/QpdK658bOzb5wx+Kf19B29YqlZJvfGPTXw9taVl/fuKiRV1aFgC9yKxZyWuvre8pf6mlJWlqWn92GwAAkCS54YZk+fKN/wrd2rp++bd/K/59BW1b67nnkv/+7/WB26b07bv+66MAsCVuumnjvyG8rbV1/TYAAECS9V8+rHyH1GvduvWXPy6aoG1rvf76u29TUZG88Ubn1wJA77R6dTHbAADANmL16vV/j34nb75Z/PsK2rZWfX3ynve88zZvvZXsv3/X1ANA7zNy5J+v3roxffuu3wYAAEiSDB/+zr9CV1Z2TlQjaNta/fsnJ5+c9Omz8ecrKpIddkhOOKFLywKgF/nCF/5884ONWbcu+fznu64eAADo5v75n9/5V+jW1uSLXyz+fQVtRfjmN5NhwzYM2/r0Wb/8x38k1dXlqQ2Anu+oo5LTTlv/3//z5jtvX3RiypTkiCO6vCwAAOiuhg1LLr54/X//z2u1VVSsX/7X/0r+8R+Lf19BWxFqa5Pf/CY566xkp53Wr6usXH+v2PvvT449trz1AdCzVVQkl1+e/PjHyb77/nn9/vsn116bfOc75asNAAC6qS9/OfnZz5KPfOTP63bfff3dRm+8cdNfTtwaFaXSO90uc9vU3Nyc2traNDU1paamZvMGt7YmTU3J9tsnVVWdUyBAD7NVx9VeaKv3R3Pz+vBtwIDiiwPogfSZDdknAO2tXr3+q6Q1Ne2/JNJRHT2uvsNl4dgilZXJjjuWuwoAejMfmAAAYLO8230si+KrowAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFCAsgZt9957b8aOHZshQ4akoqIit99+e7vnb7311hx99NF53/vel4qKiixatOhdX3P27NmpqKhot1RXV3fOBADo1vQZAACgK5U1aFu9enWGDx+emTNnbvL5ww47LBdddNFmvW5NTU1WrFjRtjz77LNFlAtAD6PPAAAAXalvOd/8mGOOyTHHHLPJ5ydMmJAkWbp06Wa9bkVFRQYNGrQ1pQHQC+gzAABAV+qV12h77bXXsuuuu6a+vj7HH398HnvssXfcfs2aNWlubm63AMCm6DMAAMDG9Lqgbe+9987VV1+dn/3sZ7nuuuvS2tqaQw89NH/60582OWb69Ompra1tW+rr67uwYgB6En0GAADYlF4XtI0ePToTJ07MiBEjcvjhh+fWW2/NLrvskn//93/f5JipU6emqampbVm+fHkXVgxAT6LPAAAAm1LWa7R1he222y4jR47M4sWLN7lNVVVVqqqqurAqAHoLfQYAAHhbrzuj7S+1tLTk0UcfzeDBg8tdCgC9kD4DAAC8raxntL322mvtzgBYsmRJFi1alJ122invf//78/LLL2fZsmV5/vnnkyRPPvlkkmTQoEFtd3ubOHFihg4dmunTpydJvvnNb+YjH/lI9tprr7z66qu55JJL8uyzz+aUU07p4tkBUG76DAAA0JXKGrQ99NBDOfLII9seT5kyJUkyadKkzJ49O3fccUcmT57c9vyJJ56YJJk2bVrOO++8JMmyZctSWfnnE/NeeeWVnHrqqWlsbMyOO+6Ygw46KPfff3/23XffLpgRAN2JPgMAAHSlilKpVCp3Ed1Nc3Nzamtr09TUlJqamnKXA9DjOa62Z38AFKu7H1fvvffeXHLJJVm4cGFWrFiR2267LePGjXvHMfPnz8+UKVPy2GOPpb6+PmeffXb+6Z/+qcPv2d33CUBP09Hjaq+/RhsAAEA5rV69OsOHD8/MmTM7tP2SJUty3HHH5cgjj8yiRYtyxhln5JRTTsncuXM7uVIAtlavv+soAABAOR1zzDE55phjOrz9rFmzsvvuu+fSSy9Nkuyzzz6577778t3vfjdjxozprDIBKIAz2gAAALqRBQsWpKGhod26MWPGZMGCBZscs2bNmjQ3N7dbAOh6gjYAAIBupLGxMXV1de3W1dXVpbm5OW+88cZGx0yfPj21tbVtS319fVeUCsBfELQBAAD0cFOnTk1TU1Pbsnz58nKXBLBNco02AACAbmTQoEFZuXJlu3UrV65MTU1N+vfvv9ExVVVVqaqq6oryAHgHzmgDAADoRkaPHp158+a1W3f33Xdn9OjRZaoIgI4StAEAAHSi1157LYsWLcqiRYuSJEuWLMmiRYuybNmyJOu/9jlx4sS27T/3uc/lmWeeyVe+8pU88cQT+cEPfpAbb7wxZ555ZjnKB2AzCNoAAAA60UMPPZSRI0dm5MiRSZIpU6Zk5MiROffcc5MkK1asaAvdkmT33XfPnXfembvvvjvDhw/PpZdemquuuipjxowpS/0AdJxrtAEAAHSiI444IqVSaZPPz549e6NjHnnkkU6sCoDO4Iw2AAAAACiAoA0AAAAACiBoAwAAAIACCNoAAAAAoACCNgAAAAAogKANAAAAAAogaAMAAACAAgjaAAAAAKAAgjYAAAAAKICgDQAAAAAKIGgDAAAAgAII2gAAAACgAII2AAAAACiAoA0AAAAACiBoAwAAAIACCNoAAAAAoACCNgAAAAAogKANAAAAAAogaAMAAACAAgjaAAAAAKAAgjYAAAAAKICgDQAAAAAKIGgDAAAAgAII2gAAAACgAII2AAAAACiAoA0AAAAACiBoAwAAAIACCNoAAAAAoACCNgAAAAAogKANAAAAAAogaAMAAACAAgjaAAAAAKAAgjYAAAAAKICgDQAAAAAKIGgDAAAAgAII2gAAAACgAII2AAAAACiAoA0AAAAACiBoAwAAAIACCNoAAAAAoACCNgAAAAAogKANAAAAAAogaAMAAACAApQ1aLv33nszduzYDBkyJBUVFbn99tvbPX/rrbfm6KOPzvve975UVFRk0aJFHXrdm266KcOGDUt1dXUOOOCA3HXXXcUXD0C3p88AAABdqaxB2+rVqzN8+PDMnDlzk88fdthhueiiizr8mvfff39OOumknHzyyXnkkUcybty4jBs3Lr///e+LKhuAHkKfAQAAulJFqVQqlbuIJKmoqMhtt92WcePGbfDc0qVLs/vuu+eRRx7JiBEj3vF1xo8fn9WrV+fnP/9527qPfOQjGTFiRGbNmtWhWpqbm1NbW5umpqbU1NRszjQA2IjucFzVZwB6L8fVDdknAMXq6HG1112jbcGCBWloaGi3bsyYMVmwYMEmx6xZsybNzc3tFgDYGH0GAADYlF4XtDU2Nqaurq7durq6ujQ2Nm5yzPTp01NbW9u21NfXd3aZAPRQ+gwAALApvS5o2xJTp05NU1NT27J8+fJylwRAL6LPAADAtqFvuQso2qBBg7Jy5cp261auXJlBgwZtckxVVVWqqqo6uzQAegF9BgAA2JRed0bb6NGjM2/evHbr7r777owePbpMFQHQm+gzAADAppT1jLbXXnstixcvbnu8ZMmSLFq0KDvttFPe//735+WXX86yZcvy/PPPJ0mefPLJJOvPJnj7zIGJEydm6NChmT59epLk9NNPz+GHH55LL700xx13XK6//vo89NBD+eEPf9jFswOg3PQZAACgK5X1jLaHHnooI0eOzMiRI5MkU6ZMyciRI3PuuecmSe64446MHDkyxx13XJLkxBNPzMiRIzNr1qy211i2bFlWrFjR9vjQQw/NT3/60/zwhz/M8OHDc/PNN+f222/P/vvv34UzA6A70GcAAICuVFEqlUrlLqK7aW5uTm1tbZqamlJTU1PucgB6PMfV9uwPgGL1hOPqzJkzc8kll6SxsTHDhw/P5ZdfnkMOOWST28+YMSNXXHFFli1blp133jmf+tSnMn369FRXV3fo/XrCPgHoSTp6XO1112gDAADoTm644YZMmTIl06ZNy8MPP5zhw4dnzJgxeeGFFza6/U9/+tOcddZZmTZtWh5//PH86Ec/yg033JCvfe1rXVw5AJtL0AYAANCJLrvsspx66qmZPHly9t1338yaNSvbb799rr766o1uf//99+ejH/1oPv3pT2e33XbL0UcfnZNOOikPPvhgF1cOwOYStAEAAHSStWvXZuHChWloaGhbV1lZmYaGhixYsGCjYw499NAsXLiwLVh75plnctddd+XYY4/d5PusWbMmzc3N7RYAul5Z7zoKAADQm7300ktpaWlJXV1du/V1dXV54oknNjrm05/+dF566aUcdthhKZVKWbduXT73uc+941dHp0+fnm984xuF1g7A5nNGGwAAQDcyf/78XHjhhfnBD36Qhx9+OLfeemvuvPPOnH/++ZscM3Xq1DQ1NbUty5cv78KKAXibM9oAAAA6yc4775w+ffpk5cqV7davXLkygwYN2uiYc845JxMmTMgpp5ySJDnggAOyevXqfPazn83Xv/71VFZueL5EVVVVqqqqip8AAJvFGW0AAACdpF+/fjnooIMyb968tnWtra2ZN29eRo8evdExr7/++gZhWp8+fZIkpVKp84oFYKs5ow0AAKATTZkyJZMmTcrBBx+cQw45JDNmzMjq1aszefLkJMnEiRMzdOjQTJ8+PUkyduzYXHbZZRk5cmRGjRqVxYsX55xzzsnYsWPbAjcAuidBGwAAQCcaP358XnzxxZx77rlpbGzMiBEjMmfOnLYbJCxbtqzdGWxnn312KioqcvbZZ+e5557LLrvskrFjx+Zb3/pWuaYAQAdVlJx7vIHm5ubU1tamqakpNTU15S4HoMdzXG3P/gAoluPqhuwTgGJ19LjqGm0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFKGvQdu+992bs2LEZMmRIKioqcvvtt7d7vlQq5dxzz83gwYPTv3//NDQ05KmnnnrH1zzvvPNSUVHRbhk2bFgnzgKA7kyvAQAAukpZg7bVq1dn+PDhmTlz5kafv/jii/O9730vs2bNygMPPJD3vOc9GTNmTN588813fN399tsvK1asaFvuu+++zigfgB5ArwEAALpK33K++THHHJNjjjlmo8+VSqXMmDEjZ599do4//vgkyTXXXJO6urrcfvvtOfHEEzf5un379s2gQYM6pWYAeha9BgAA6Crd9hptS5YsSWNjYxoaGtrW1dbWZtSoUVmwYME7jn3qqacyZMiQ7LHHHvmHf/iHLFu27B23X7NmTZqbm9stAPR+XdVr9BkAANg2dNugrbGxMUlSV1fXbn1dXV3bcxszatSozJ49O3PmzMkVV1yRJUuW5GMf+1hWrVq1yTHTp09PbW1t21JfX1/MJADo1rqq1+gzAACwbei2QduWOuaYY3LCCSfkwAMPzJgxY3LXXXfl1VdfzY033rjJMVOnTk1TU1Pbsnz58i6sGICeZnN7jT4DAADbhm4btL193ZuVK1e2W79y5crNuibODjvskA9+8INZvHjxJrepqqpKTU1NuwWA3q+reo0+AwAA24ZuG7TtvvvuGTRoUObNm9e2rrm5OQ888EBGjx7d4dd57bXX8vTTT2fw4MGdUSYAPZheAwAAFKmsQdtrr72WRYsWZdGiRUnWX5R60aJFWbZsWSoqKnLGGWfkggsuyB133JFHH300EydOzJAhQzJu3Li21/jEJz6R73//+22Pv/SlL+VXv/pVli5dmvvvvz9/93d/lz59+uSkk07q4tkB0B3oNQAAQFfpW843f+ihh3LkkUe2PZ4yZUqSZNKkSZk9e3a+8pWvZPXq1fnsZz+bV199NYcddljmzJmT6urqtjFPP/10XnrppbbHf/rTn3LSSSfl//v//r/ssssuOeyww/Jf//Vf2WWXXbpuYgB0G3oNAADQVSpKpVKp3EV0N83NzamtrU1TU5Pr6AAUwHG1PfsDoFg94bg6c+bMXHLJJWlsbMzw4cNz+eWX55BDDtnk9q+++mq+/vWv59Zbb83LL7+cXXfdNTNmzMixxx7boffrCfsEoCfp6HG1rGe0AQAA9HY33HBDpkyZklmzZmXUqFGZMWNGxowZkyeffDIDBw7cYPu1a9fmqKOOysCBA3PzzTdn6NChefbZZ7PDDjt0ffEAbBZBGwAAQCe67LLLcuqpp2by5MlJklmzZuXOO+/M1VdfnbPOOmuD7a+++uq8/PLLuf/++7PddtslSXbbbbeuLBmALdRt7zoKAADQ061duzYLFy5MQ0ND27rKyso0NDRkwYIFGx1zxx13ZPTo0TnttNNSV1eX/fffPxdeeGFaWlo2+T5r1qxJc3NzuwWAridoAwAA6CQvvfRSWlpaUldX1259XV1dGhsbNzrmmWeeyc0335yWlpbcddddOeecc3LppZfmggsu2OT7TJ8+PbW1tW1LfX19ofMAoGMEbQAAAN1Ia2trBg4cmB/+8Ic56KCDMn78+Hz961/PrFmzNjlm6tSpaWpqaluWL1/ehRUD8DbXaAMAAOgkO++8c/r06ZOVK1e2W79y5coMGjRoo2MGDx6c7bbbLn369Glbt88++6SxsTFr165Nv379NhhTVVWVqqqqYosHYLM5ow0AAKCT9OvXLwcddFDmzZvXtq61tTXz5s3L6NGjNzrmox/9aBYvXpzW1ta2dX/84x8zePDgjYZsAHQfgjYAAIBONGXKlFx55ZX5yU9+kscffzyf//zns3r16ra7kE6cODFTp05t2/7zn/98Xn755Zx++un54x//mDvvvDMXXnhhTjvttHJNAYAO8tVRAACATjR+/Pi8+OKLOffcc9PY2JgRI0Zkzpw5bTdIWLZsWSor/3wORH19febOnZszzzwzBx54YIYOHZrTTz89X/3qV8s1BQA6qKJUKpXKXUR309zcnNra2jQ1NaWmpqbc5QD0eI6r7dkfAMVyXN2QfQJQrI4eV311FAAAAAAKIGgDAAAAgAII2gAAAACgAII2AAAAACiAoA0AAAAACtC33AUA0L09+WTyk58kzz+fDBqUTJiQ7LdfuasCoLd45ZXkmmuS3/0uqapKxo5NxoxJ+vQpd2UAsPkEbQBsVGtrcvrpyfe/n/Ttm5RKSUVFctFFycknJ7NmrV8PAFvqlluSf/zHZM2aPwdrs2Yl+++fzJmTDB1a3voAYHP56igAG/Wtb60P2ZJk3bqkpWX9/ybJ1VcnX/96+WoDoOd74IFk/Pj1IVuptL7HvN1nnngiOfroPz8GgJ5ii4K2j3/84/nGN76xwfpXXnklH//4x7e6KADK6/XXk0su2fTzpVLyve8lTU2d8/76DEDv9+1vrz9TulTa8Ll165I//CG5666ur+ttkyZNyr333lu+AgDokbYoaJs/f36+//3vZ9y4cVm9enXb+rVr1+ZXv/pVYcUBUB733pusWvXO27z5ZnL33Z3z/voMQO/W0pL83//7zmes9e2b3HZb19X0l5qamtLQ0JAPfOADufDCC/Pcc8+VrxgAeowt/uroPffck8bGxnzkIx/J0qVLCywJgHJ7/fWObffGG51Xgz4D0Hu9fUmCd9La2rl95t3cfvvtee655/L5z38+N9xwQ3bbbbccc8wxufnmm/PWW2+VrzAAurUtDtoGDx6cX/3qVznggAPy4Q9/OPPnzy+wLADKaf/9i91uS+gzAL1XVVWy++7rvzr6Tjqzz3TELrvskilTpuR3v/tdHnjggey1116ZMGFChgwZkjPPPDNPPfVUeQsEoNvZoqCt4v91xKqqqvz0pz/N6aefnr/5m7/JD37wg0KLA6A8PvjB5PDD/3wHuL/Up09y0EHJyJGd8/76DEDv9y//8s7PV1Ymn/lM19TyblasWJG77747d999d/r06ZNjjz02jz76aPbdd99897vfLXd5AHQjfbdkUOkvrlh69tlnZ5999smkSZMKKQqA8rvqqmT06OTVV9tfQ6dv3+S9701+8pPOe299BqD3O+205Oc/T+bPX/810bf16bP+a6VXXJEMGVK28vLWW2/ljjvuyI9//OP84he/yIEHHpgzzjgjn/70p1NTU5Mkue222/KZz3wmZ555ZvkKBaBb2aKgbcmSJdlll13arfvkJz+ZYcOG5aGHHiqkMADKa6+9kocfTqZPT2bPXn+dnKqq5B//Mfna15I99ui899ZnAHq/fv3W31X03/4tufzy5E9/Wr/+r/86mTo1Oeqo8tY3ePDgtLa25qSTTsqDDz6YESNGbLDNkUcemR122KHLawOg+6oo/eVpA6S5uTm1tbVpampq+2sVwLZs3bqkuTkZMCDZbrvNH++42p79AdBeqZQ0Na0P37bffvPHd8Zx9dprr80JJ5yQ6urqQl6vq+k1AMXq6HF1i85oA2Db0rdvstNO5a4CgN6qoiLpbieGTZgwodwlANADbfFdRwEAAACAPxO0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABShr0Hbvvfdm7NixGTJkSCoqKnL77be3e75UKuXcc8/N4MGD079//zQ0NOSpp55619edOXNmdtttt1RXV2fUqFF58MEHO2kGAHR3eg0AANBVyhq0rV69OsOHD8/MmTM3+vzFF1+c733ve5k1a1YeeOCBvOc978mYMWPy5ptvbvI1b7jhhkyZMiXTpk3Lww8/nOHDh2fMmDF54YUXOmsaAHRjeg0AANBVKkqlUqncRSRJRUVFbrvttowbNy7J+jMMhgwZkn/913/Nl770pSRJU1NT6urqMnv27Jx44okbfZ1Ro0blwx/+cL7//e8nSVpbW1NfX59/+Zd/yVlnndWhWpqbm1NbW5umpqbU1NRs/eQAtnHd5bjaXXpNd9kfAL2F4+qG7BOAYnX0uNptr9G2ZMmSNDY2pqGhoW1dbW1tRo0alQULFmx0zNq1a7Nw4cJ2YyorK9PQ0LDJMQBsu/QaAACgSH3LXcCmNDY2Jknq6urara+rq2t77i+99NJLaWlp2eiYJ554YpPvtWbNmqxZs6btcXNz85aWDUAP0lW9Rp8BAIBtQ7c9o60rTZ8+PbW1tW1LfX19uUsCoBfRZwAAYNvQbYO2QYMGJUlWrlzZbv3KlSvbnvtLO++8c/r06bNZY5Jk6tSpaWpqaluWL1++ldUD0BN0Va/RZwAAYNvQbYO23XffPYMGDcq8efPa1jU3N+eBBx7I6NGjNzqmX79+Oeigg9qNaW1tzbx58zY5JkmqqqpSU1PTbgGg9+uqXqPPAADAtqGs12h77bXXsnjx4rbHS5YsyaJFi7LTTjvl/e9/f84444xccMEF+cAHPpDdd98955xzToYMGdJ2t7gk+cQnPpG/+7u/yxe/+MUkyZQpUzJp0qQcfPDBOeSQQzJjxoysXr06kydP7urpAdAN6DUAAEBXKWvQ9tBDD+XII49sezxlypQkyaRJkzJ79ux85StfyerVq/PZz342r776ag477LDMmTMn1dXVbWOefvrpvPTSS22Px48fnxdffDHnnntuGhsbM2LEiMyZM2eDi1YDsG3QawAAgK5SUSqVSuUuortpbm5ObW1tmpqafL0HoACOq+3ZHwDFclzdkH0CUKyOHle77TXaAAAAeouZM2dmt912S3V1dUaNGpUHH3ywQ+Ouv/76VFRUtLukAQDdl6ANAACgE91www2ZMmVKpk2blocffjjDhw/PmDFj8sILL7zjuKVLl+ZLX/pSPvaxj3VRpQBsLUEbAABAJ7rsssty6qmnZvLkydl3330za9asbL/99rn66qs3OaalpSX/8A//kG984xvZY489urBaALaGoA0AAKCTrF27NgsXLkxDQ0PbusrKyjQ0NGTBggWbHPfNb34zAwcOzMknn9yh91mzZk2am5vbLQB0PUEbAABAJ3nppZfS0tKywZ2p6+rq0tjYuNEx9913X370ox/lyiuv7PD7TJ8+PbW1tW1LfX39VtUNwJYRtAEAAHQTq1atyoQJE3LllVdm55137vC4qVOnpqmpqW1Zvnx5J1YJwKb0LXcBAAAAvdXOO++cPn36ZOXKle3Wr1y5MoMGDdpg+6effjpLly7N2LFj29a1trYmSfr27Zsnn3wye+655wbjqqqqUlVVVXD1AGwuZ7QBAAB0kn79+uWggw7KvHnz2ta1trZm3rx5GT169AbbDxs2LI8++mgWLVrUtvzt3/5tjjzyyCxatMhXQgG6OWe0AQAAdKIpU6Zk0qRJOfjgg3PIIYdkxowZWb16dSZPnpwkmThxYoYOHZrp06enuro6+++/f7vxO+ywQ5JssB6A7kfQBgAA0InGjx+fF198Meeee24aGxszYsSIzJkzp+0GCcuWLUtlpS8bAfQGFaVSqVTuIrqb5ubm1NbWpqmpKTU1NeUuB6DHc1xtz/4AKJbj6obsE4BidfS46s8mAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABSg2wdtq1atyhlnnJFdd901/fv3z6GHHprf/va3m9x+/vz5qaio2GBpbGzswqoB6En0GgAAoAh9y13AuznllFPy+9//Ptdee22GDBmS6667Lg0NDfnDH/6QoUOHbnLck08+mZqamrbHAwcO7IpyAeiB9BoAAKAI3fqMtjfeeCO33HJLLr744vz1X/919tprr5x33nnZa6+9csUVV7zj2IEDB2bQoEFtS2Vlt54qAGWi1wAAAEXp1p8I1q1bl5aWllRXV7db379//9x3333vOHbEiBEZPHhwjjrqqPzmN795x23XrFmT5ubmdgsA24au6DX6DAAAbBu6ddA2YMCAjB49Oueff36ef/75tLS05LrrrsuCBQuyYsWKjY4ZPHhwZs2alVtuuSW33HJL6uvrc8QRR+Thhx/e5PtMnz49tbW1bUt9fX1nTQmAbqYreo0+AwAA24aKUqlUKncR7+Tpp5/OZz7zmdx7773p06dPPvShD+WDH/xgFi5cmMcff7xDr3H44Yfn/e9/f6699tqNPr9mzZqsWbOm7XFzc3Pq6+vT1NTU7to7AGyZ5ubm1NbWdtvjamf3Gn0GoHN19z5TDvYJQLE6elzt1me0Jcmee+6ZX/3qV3nttdeyfPnyPPjgg3nrrbeyxx57dPg1DjnkkCxevHiTz1dVVaWmpqbdAsC2o7N7jT4DAADbhm4ftL3tPe95TwYPHpxXXnklc+fOzfHHH9/hsYsWLcrgwYM7sToAegO9BgAA2Bp9y13Au5k7d25KpVL23nvvLF68OF/+8pczbNiwTJ48OUkyderUPPfcc7nmmmuSJDNmzMjuu++e/fbbL2+++Wauuuqq/Od//md+8YtflHMaAHRjeg0AAFCEbh+0NTU1ZerUqfnTn/6UnXbaKZ/85CfzrW99K9ttt12SZMWKFVm2bFnb9mvXrs2//uu/5rnnnsv222+fAw88MPfcc0+OPPLIck0BgG5OrwEAAIrQ7W+GUA4uHApQLMfV9uwPgGI5rm7IPgEoVq+5GQIAAEBPN3PmzOy2226prq7OqFGj8uCDD25y2yuvvDIf+9jHsuOOO2bHHXdMQ0PDO24PQPchaAMAAOhEN9xwQ6ZMmZJp06bl4YcfzvDhwzNmzJi88MILG91+/vz5Oemkk/LLX/4yCxYsSH19fY4++ug899xzXVw5AJtL0AYAANCJLrvsspx66qmZPHly9t1338yaNSvbb799rr766o1u/x//8R/5whe+kBEjRmTYsGG56qqr0tramnnz5nVx5QBsLkEbAABAJ1m7dm0WLlyYhoaGtnWVlZVpaGjIggULOvQar7/+et56663stNNOm9xmzZo1aW5ubrcA0PUEbQAAAJ3kpZdeSktLS+rq6tqtr6urS2NjY4de46tf/WqGDBnSLqz7S9OnT09tbW3bUl9fv1V1A7BlBG0AAADd1Le//e1cf/31ue2221JdXb3J7aZOnZqmpqa2Zfny5V1YJQBv61vuAgAAAHqrnXfeOX369MnKlSvbrV+5cmUGDRr0jmO/853v5Nvf/nbuueeeHHjgge+4bVVVVaqqqra6XgC2jjPaAAAAOkm/fv1y0EEHtbuRwds3Nhg9evQmx1188cU5//zzM2fOnBx88MFdUSoABXBGGwAAQCeaMmVKJk2alIMPPjiHHHJIZsyYkdWrV2fy5MlJkokTJ2bo0KGZPn16kuSiiy7Kueeem5/+9KfZbbfd2q7l9t73vjfvfe97yzYPAN6doA0AAKATjR8/Pi+++GLOPffcNDY2ZsSIEZkzZ07bDRKWLVuWyso/f9noiiuuyNq1a/OpT32q3etMmzYt5513XleWDsBmErQBAAB0si9+8Yv54he/uNHn5s+f3+7x0qVLO78gADqFa7QBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUIBuH7StWrUqZ5xxRnbdddf0798/hx56aH7729++45j58+fnQx/6UKqqqrLXXntl9uzZXVMsAD2SXgMAABSh2wdtp5xySu6+++5ce+21efTRR3P00UenoaEhzz333Ea3X7JkSY477rgceeSRWbRoUc4444yccsopmTt3bhdXDkBPodcAAABFqCiVSqVyF7Epb7zxRgYMGJCf/exnOe6449rWH3TQQTnmmGNywQUXbDDmq1/9au688878/ve/b1t34okn5tVXX82cOXM69L7Nzc2pra1NU1NTampqtn4iANu47nxcLUev6c77A6AnclzdkH0CUKyOHle79Rlt69atS0tLS6qrq9ut79+/f+67776NjlmwYEEaGhrarRszZkwWLFiwyfdZs2ZNmpub2y0AbBu6otfoMwAAsG3o1kHbgAEDMnr06Jx//vl5/vnn09LSkuuuuy4LFizIihUrNjqmsbExdXV17dbV1dWlubk5b7zxxkbHTJ8+PbW1tW1LfX194XMBoHvqil6jzwAAwLahWwdtSXLttdemVCpl6NChqaqqyve+972cdNJJqawsrvSpU6emqampbVm+fHlhrw1A99fZvUafAQCAbUPfchfwbvbcc8/86le/yurVq9Pc3JzBgwdn/Pjx2WOPPTa6/aBBg7Jy5cp261auXJmampr0799/o2OqqqpSVVVVeO0A9Ayd3Wv0GQAA2DZ0+zPa3vae97wngwcPziuvvJK5c+fm+OOP3+h2o0ePzrx589qtu/vuuzN69OiuKBOAHkyvAQAAtka3D9rmzp2bOXPmZMmSJbn77rtz5JFHZtiwYZk8eXKS9V/HmThxYtv2n/vc5/LMM8/kK1/5Sp544on84Ac/yI033pgzzzyzXFMAoJvTawAAgCJ0+6Ctqakpp512WoYNG5aJEyfmsMMOy9y5c7PddtslSVasWJFly5a1bb/77rvnzjvvzN13353hw4fn0ksvzVVXXZUxY8aUawoAdHN6DQAAUISKUqlUKncR3U1zc3Nqa2vT1NSUmpqaDo9bsWpFZi+anT++/McM6DcgJ+x7Qg57/2GpqKjoxGoBur8tPa72Vlu6P95c92ZufOzG3PvsvalIRQ7f7fB8at9PpbpvdSdWC9D96TMb2pJ9Umptzb3/9/LcvOBHeW3dGxm24175p09fkrrd9+/kagG6v44eV7v9GW09xeUPXJ7679bn7F+enev++7pc8dAV+evZf52PX/PxNL3ZVO7yAOjhfvvcb/P+774/k26flJ/87ieZ/bvZmXDbhOw2Y7c8vOLhcpcHwLuYOXNmdtttt1RXV2fUqFF58MEH33H7m266KcOGDUt1dXUOOOCA3HXXXZ1a3ysrluSvp+yQIxadkVn9Hs112y/O19bOyV/9+IDMuuzTnfreAL2JoK0Atz5+a/73nP+dllJLWkutWde6Luta1yVJfv3srzP+5vFlrhCAnmzFqhU56tqj8vIbLydJuz7z0usvpeGahryw+oVylgjAO7jhhhsyZcqUTJs2LQ8//HCGDx+eMWPG5IUXNn7svv/++3PSSSfl5JNPziOPPJJx48Zl3Lhx+f3vf98p9ZVaW/PJiz6UBbWrkiTr+qxfWivX/+/nV/2f/Py6czvlvQF6G0HbViqVSvnGr76Rimz866EtpZbMfXpuFjUu6trCAOg1Zj00K6+tfS0tpZYNnmsptaRpTVOueviqMlQGQEdcdtllOfXUUzN58uTsu+++mTVrVrbffvtcffXVG93+3/7t3/I3f/M3+fKXv5x99tkn559/fj70oQ/l+9//fqfU99t7fpJf7vhqWjbx6bCyNblg4WWd8t4AvY2gbSs9t+q5/PfK/04pm77UXd+Kvrn9idu7rigAepWb/nDTRkO2t7WWWnPTYzd1YUUAdNTatWuzcOHCNDQ0tK2rrKxMQ0NDFixYsNExCxYsaLd9kowZM2aT2yfJmjVr0tzc3G7pqNt/fWX6brrNpLUyeWCH1Xlh6WMdfk2AbZWgbSu9/tbr77pNRUVF3njrjS6oBoDeaPVbqwvZBoCu99JLL6WlpSV1dXXt1tfV1aWxsXGjYxobGzdr+ySZPn16amtr25b6+voO1/j6ujc28f2c9t5Y/WqHXxNgWyVo20r1NfV5z3bvecdt3mp9K/sPdKceALbMyEEj07ey7yaf71vZNyMHjezCigDobqZOnZqmpqa2Zfny5R0ee8Dg4XnrXT4Z1q5JBu8xfCurBOj9BG1bqf92/XPyyJPTp6LPRp+vSEV2qN4hJ+x3QhdXBkBv8YUPf6Ht5gcbs651XT7/4c93YUUAdNTOO++cPn36ZOXKle3Wr1y5MoMGDdromEGDBm3W9klSVVWVmpqadktHnTjh4gxYm1Rs4mo4fVqTUysOTr/+7+3wawJsqwRtBfjmkd/MsJ2HbRC29anokz6VffIf/+s/Ut23ukzVAdDTHbXHUTntw6clSbub71T+vzY+5SNTcsRuR5SjNADeRb9+/XLQQQdl3rx5betaW1szb968jB49eqNjRo8e3W77JLn77rs3uf3Wes+OA3PdPl9LZSnp8xfXauvTmuzfXJ1zp/ysU94boLcRtBWgtro2v/nMb3LWYWdlp/47JUkqKyozdu+xuf8z9+fYDxxb5goB6MkqKipy+TGX58fH/zj77rJv2/r96/bPtX93bb5z9HfKWB0A72bKlCm58sor85Of/CSPP/54Pv/5z2f16tWZPHlykmTixImZOnVq2/ann3565syZk0svvTRPPPFEzjvvvDz00EP54he/2Gk1/u2Eb+W+0Vfm2FV1qWxdv27nNyrytXwsvz776Qx435BOe2+A3qSiVCpt+naZ26jm5ubU1tamqalps065Ttbf+a3pzaZsv932qepb1UkVAvQsW3Nc7Y22dn80r2lORSoyoGpAJ1QH0PP0hD7z/e9/P5dcckkaGxszYsSIfO9738uoUaOSJEcccUR22223zJ49u237m266KWeffXaWLl2aD3zgA7n44otz7LEd/wP+1uyTN197NW+seiW1A+tT2WfT1wgF2JZ09LgqaNuIntCoAXoSx9X27A+AYjmubsg+AShWR4+rvjoKAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEAB+pa7gO6oVColSZqbm8tcCUDv8Pbx9O3j67ZOnwEolj6zIb0GoFgd7TWCto1YtWpVkqS+vr7MlQD0LqtWrUptbW25yyg7fQagc+gzf6bXAHSOd+s1FSV/9tlAa2trnn/++QwYMCAVFRWbPb65uTn19fVZvnx5ampqOqHCbYP9WAz7sRj249YplUpZtWpVhgwZkspKVy3QZ96dOfZ8vX1+iTl2J/rMhram1/SUf/fuzn4shv1YDPtx63W01zijbSMqKyvzV3/1V1v9OjU1NX6AC2A/FsN+LIb9uOWcYfBn+kzHmWPP19vnl5hjd6HPtFdEr+kJ/+49gf1YDPuxGPbj1ulIr/HnHgAAAAAogKANAAAAAAogaOsEVVVVmTZtWqqqqspdSo9mPxbDfiyG/Uh3si38PJpjz9fb55eYI72Xf/di2I/FsB+LYT92HTdDAAAAAIACOKMNAAAAAAogaAMAAACAAgjaAAAAAKAAgjYAAAAAKICgbQvNnDkzu+22W6qrqzNq1Kg8+OCD77j9TTfdlGHDhqW6ujoHHHBA7rrrri6qtHvbnP04e/bsVFRUtFuqq6u7sNru5957783YsWMzZMiQVFRU5Pbbb3/XMfPnz8+HPvShVFVVZa+99srs2bM7vc7ubnP34/z58zf4WayoqEhjY2PXFMw2YVvoM5szxyuvvDIf+9jHsuOOO2bHHXdMQ0PDu+6T7mBz/x3fdv3116eioiLjxo3r3AK30ubO79VXX81pp52WwYMHp6qqKh/84Ae7/c/q5s5xxowZ2XvvvdO/f//U19fnzDPPzJtvvtlF1W4+v0tsu7aFPtMVfJ7Zeo5DxfCZpvsQtG2BG264IVOmTMm0adPy8MMPZ/jw4RkzZkxeeOGFjW5///3356STTsrJJ5+cRx55JOPGjcu4cePy+9//vosr7142dz8mSU1NTVasWNG2PPvss11YcfezevXqDB8+PDNnzuzQ9kuWLMlxxx2XI488MosWLcoZZ5yRU045JXPnzu3kSru3zd2Pb3vyySfb/TwOHDiwkypkW7Mt9JnNneP8+fNz0kkn5Ze//GUWLFiQ+vr6HH300Xnuuee6uPKO25I+lyRLly7Nl770pXzsYx/rokq3zObOb+3atTnqqKOydOnS3HzzzXnyySdz5ZVXZujQoV1cecdt7hx/+tOf5qyzzsq0adPy+OOP50c/+lFuuOGGfO1rX+viyjvO7xLbpm2hz3QFn2eK4ThUDJ9pupESm+2QQw4pnXbaaW2PW1paSkOGDClNnz59o9v//d//fem4445rt27UqFGlf/7nf+7UOru7zd2PP/7xj0u1tbVdVF3Pk6R02223veM2X/nKV0r77bdfu3Xjx48vjRkzphMr61k6sh9/+ctflpKUXnnllS6piW3PttBnNneOf2ndunWlAQMGlH7yk590VolbbUvmuG7dutKhhx5auuqqq0qTJk0qHX/88V1Q6ZbZ3PldccUVpT322KO0du3aripxq23uHE877bTSxz/+8XbrpkyZUvroRz/aqXUWxe8S245toc90BZ9niuc4VAyfacrLGW2bae3atVm4cGEaGhra1lVWVqahoSELFizY6JgFCxa02z5JxowZs8nttwVbsh+T5LXXXsuuu+6a+vr6HH/88Xnssce6otxew89isUaMGJHBgwfnqKOOym9+85tyl0MvsS30mS3tAf/T66+/nrfeeis77bRTZ5W5VbZ0jt/85jczcODAnHzyyV1R5hbbkvndcccdGT16dE477bTU1dVl//33z4UXXpiWlpauKnuzbMkcDz300CxcuLDtq2PPPPNM7rrrrhx77LFdUnNX6GnHGza0LfSZruDzTPn4eSyWzzTFE7RtppdeeiktLS2pq6trt76urm6T32VubGzcrO23BVuyH/fee+9cffXV+dnPfpbrrrsura2tOfTQQ/OnP/2pK0ruFTb1s9jc3Jw33nijTFX1PIMHD86sWbNyyy235JZbbkl9fX2OOOKIPPzww+UujV5gW+gzWzLHv/TVr341Q4YM2eAX7e5iS+Z433335Uc/+lGuvPLKrihxq2zJ/J555pncfPPNaWlpyV133ZVzzjknl156aS644IKuKHmzbckcP/3pT+eb3/xmDjvssGy33XbZc889c8QRR3Trr45uLr9L9HzbQp/pCj7PlI/jUDF8puk8fctdAHTU6NGjM3r06LbHhx56aPbZZ5/8+7//e84///wyVsa2Zu+9987ee+/d9vjQQw/N008/ne9+97u59tpry1gZbBu+/e1v5/rrr8/8+fN7zUWkV61alQkTJuTKK6/MzjvvXO5yOkVra2sGDhyYH/7wh+nTp08OOuigPPfcc7nkkksybdq0cpdXiPnz5+fCCy/MD37wg4waNSqLFy/O6aefnvPPPz/nnHNOucsDysznGboTn2k6j6BtM+28887p06dPVq5c2W79ypUrM2jQoI2OGTRo0GZtvy3Ykv34l7bbbruMHDkyixcv7owSe6VN/SzW1NSkf//+ZaqqdzjkkENy3333lbsMeoFtoc9sTQ/4zne+k29/+9u55557cuCBB3ZmmVtlc+f49NNPZ+nSpRk7dmzbutbW1iRJ37598+STT2bPPffs3KI3w5b8Gw4ePDjbbbdd+vTp07Zun332SWNjY9auXZt+/fp1as2ba0vmeM4552TChAk55ZRTkiQHHHBAVq9enc9+9rP5+te/nsrKnv9lEr9L9HzbQp/pCj7PlI/jUOfxmaYYPb/bd7F+/frloIMOyrx589rWtba2Zt68ee3+OvE/jR49ut32SXL33XdvcvttwZbsx7/U0tKSRx99NIMHD+6sMnsdP4udZ9GiRX4WKcS20Ge2tAdcfPHFOf/88zNnzpwcfPDBXVHqFtvcOQ4bNiyPPvpoFi1a1Lb87d/+bdsd1err67uy/He1Jf+GH/3oR7N48eK2ADFJ/vjHP2bw4MHdLmRLtmyOr7/++gZh2tvBYqlU6rxiu1BPO96woW2hz3QFn2fKx89j5/GZpiDlvhtDT3T99deXqqqqSrNnzy794Q9/KH32s58t7bDDDqXGxsZSqVQqTZgwoXTWWWe1bf+b3/ym1Ldv39J3vvOd0uOPP16aNm1aabvttis9+uij5ZpCt7C5+/Eb3/hGae7cuaWnn366tHDhwtKJJ55Yqq6uLj322GPlmkLZrVq1qvTII4+UHnnkkVKS0mWXXVZ65JFHSs8++2ypVCqVzjrrrNKECRPatn/mmWdK22+/fenLX/5y6fHHHy/NnDmz1KdPn9KcOXPKNYVuYXP343e/+93S7bffXnrqqadKjz76aOn0008vVVZWlu65555yTYFeZlvoM5s7x29/+9ulfv36lW6++ebSihUr2pZVq1aVawrvanPn+Je6+11HN3d+y5YtKw0YMKD0xS9+sfTkk0+Wfv7zn5cGDhxYuuCCC8o1hXe1uXOcNm1aacCAAaX/83/+T+mZZ54p/eIXvyjtueeepb//+78v1xTeld8ltk3bQp/pCj7PFMNxqBg+03QfgrYtdPnll5fe//73l/r161c65JBDSv/1X//V9tzhhx9emjRpUrvtb7zxxtIHP/jBUr9+/Ur77bdf6c477+ziirunzdmPZ5xxRtu2dXV1pWOPPbb08MMPl6Hq7uPtWzL/5fL2fps0aVLp8MMP32DMiBEjSv369SvtsccepR//+MddXnd3s7n78aKLLirtueeeperq6tJOO+1UOuKII0r/+Z//WZ7i6bW2hT6zOXPcddddN/r/02nTpnV94Zthc/8d/6fuHrSVSps/v/vvv780atSoUlVVVWmPPfYofetb3yqtW7eui6vePJszx7feeqt03nnntfWI+vr60he+8IXSK6+80vWFd5DfJbZd20Kf6Qo+z2w9x6Fi+EzTfVSUSr3kPHYAAAAAKCPXaAMAAACAAgjaAAAAAKAAgjYAAAAAKICgDQAAAAAKIGgDAAAAgAII2gAAAACgAII2AAAAACiAoA0AAAAACiBoAwAAAIACCNoAAAAAoACCNuhhXnzxxQwaNCgXXnhh27r7778//fr1y7x588pYGQC9xTXXXJP3ve99WbNmTbv148aNy4QJE8pUFQC9hc809GYVpVKpVO4igM1z1113Zdy4cbn//vuz9957Z8SIETn++ONz2WWXlbs0AHqBN954I4MHD86VV16ZE044IUnywgsvZOjQofnFL36RI488sswVAtDT+UxDbyVogx7qtNNOyz333JODDz44jz76aH7729+mqqqq3GUB0Et84QtfyNKlS3PXXXclSS677LLMnDkzixcvTkVFRZmrA6A38JmG3kjQBj3UG2+8kf333z/Lly/PwoULc8ABB5S7JAB6kUceeSQf/vCH8+yzz2bo0KE58MADc8IJJ+Scc84pd2kA9BI+09AbuUYb9FBPP/10nn/++bS2tmbp0qXlLgeAXmbkyJEZPnx4rrnmmixcuDCPPfZY/umf/qncZQHQi/hMQ2/kjDbogdauXZtDDjkkI0aMyN57750ZM2bk0UcfzcCBA8tdGgC9yBVXXJEZM2bkqKOOylNPPZW5c+eWuyQAegmfaeitBG3QA335y1/OzTffnN/97nd573vfm8MPPzy1tbX5+c9/Xu7SAOhFmpqaMmTIkKxbty7XXHNNxo8fX+6SAOglfKaht/LVUehh5s+fnxkzZuTaa69NTU1NKisrc+211+bXv/51rrjiinKXB0AvUltbm09+8pN573vfm3HjxpW7HAB6CZ9p6M2c0QYAwCZ94hOfyH777Zfvfe975S4FAKDbE7QBALCBV155JfPnz8+nPvWp/OEPf8jee+9d7pIAALq9vuUuAACA7mfkyJF55ZVXctFFFwnZAAA6yBltAAAAAFAAN0MAAAAAgAII2gAAAACgAII2AAAAACiAoA0AAAAACiBoAwAAAIACCNoAAAAAoACCNgAAAAAogKANAAAAAAogaAMAAACAAvz/e/gHozxT/KgAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ "import matplotlib.pyplot as plt\n", "\n",