From 8246a6a1abdaac316d98c89d01ec018895fc0074 Mon Sep 17 00:00:00 2001 From: Daniel Pozsar Date: Mon, 4 Nov 2024 12:07:11 +0100 Subject: [PATCH] updated standard output and return array --- src/grogu_magn/useful.py | 97 +---- test.ipynb | 883 ++++++++++++++++++++++++++++++++------- 2 files changed, 743 insertions(+), 237 deletions(-) diff --git a/src/grogu_magn/useful.py b/src/grogu_magn/useful.py index 07aecc1..85b1486 100644 --- a/src/grogu_magn/useful.py +++ b/src/grogu_magn/useful.py @@ -23,7 +23,6 @@ from pprint import pprint 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]]) @@ -228,38 +227,7 @@ 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_pair_atomic_indices(pair, magnetic_entities, dh): - atomic_indices = "" - # 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 #############################################################################" - ) +def print_parameters(simulation_parameters): print( "================================================================================================================================================================" ) @@ -272,11 +240,8 @@ def print_output(simulation_parameters, magnetic_entities, pairs, dh, times): print( "================================================================================================================================================================" ) - try: - print("Cell [Ang]: ") - print(simulation_parameters["geom"].cell) - except: - print("Geometry could not be read.") + print("Cell [Ang]: ") + print(simulation_parameters["cell"]) print( "================================================================================================================================================================" ) @@ -288,19 +253,19 @@ def print_output(simulation_parameters, magnetic_entities, pairs, dh, times): print( "================================================================================================================================================================" ) - print("number of k points: ", simulation_parameters["kset"]) - print("k point directions: ", simulation_parameters["kdirs"]) - print( - "================================================================================================================================================================" - ) print("Parameters for the contour integral:") + print("Number of k points: ", simulation_parameters["kset"]) + print("k point directions: ", simulation_parameters["kdirs"]) print("Ebot: ", simulation_parameters["ebot"]) print("Eset: ", simulation_parameters["eset"]) print("Esetp: ", simulation_parameters["esetp"]) print( "================================================================================================================================================================" ) - print("Atomic informations: ") + + +def print_atoms_and_pairs(magnetic_entities, pairs): + print("Atomic information: ") print( "----------------------------------------------------------------------------------------------------------------------------------------------------------------" ) @@ -312,30 +277,10 @@ def print_output(simulation_parameters, magnetic_entities, pairs, dh, times): ) # 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): + # iterate over atoms + for tag, xyz in zip(mag_ent["tags"], mag_ent["xyz"]): # 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(f"{tag} {xyz[0]} {xyz[1]} {xyz[2]}") print("") print( @@ -351,21 +296,14 @@ def print_output(simulation_parameters, magnetic_entities, pairs, dh, times): ) # 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.", + f"{pair['tags'][0]} {pair['tags'][1]} {pair['Ruc']} d [Ang] Not yet." ) # print magnetic parameters - print("Isotropic: ", J_iso) - print("DMI: ", D) - print("Symmetric-anisotropy: ", J_S) + print("Isotropic: ", pair["J_iso"]) + print("DMI: ", pair["D"]) + print("Symmetric-anisotropy: ", pair["J_S"]) print("Energies for debugging: ") pprint(np.array(pair["energies"])) print( @@ -380,6 +318,9 @@ def print_output(simulation_parameters, magnetic_entities, pairs, dh, times): print( "================================================================================================================================================================" ) + + +def print_runtime_informations(times): print("Runtime information: ") print(f"Total runtime: {times['end_time'] - times['start_time']} s") print( diff --git a/test.ipynb b/test.ipynb index e3b237e..360b8e2 100644 --- a/test.ipynb +++ b/test.ipynb @@ -2,12 +2,29 @@ "cells": [ { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "[Mac:23435] shmem: mmap: an error occurred while determining whether or not /var/folders/yh/dx7xl94n3g52ts3td8qcxjcc0000gn/T//ompi.Mac.501/jf.0/46137344/sm_segment.Mac.501.2c00000.0 could be created.\n" + ] + }, + { + "data": { + "text/plain": [ + "'0.14.3'" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "import os\n", - "from sys import stdout\n", "from tqdm import tqdm\n", "from timeit import default_timer as timer\n", "\n", @@ -35,12 +52,45 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "================================================================================================================================================================\n", + "Input file: \n", + "Not yet specified.\n", + "Number of nodes in the parallel cluster: 1\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", + "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, 0]), array([0, 0, 1])]\n", + "[0 0 1] --» [array([1, 0, 0]), array([0, 1, 0])]\n", + "================================================================================================================================================================\n", + "Parameters for the contour integral:\n", + "Number of k points: 10\n", + "k point directions: xy\n", + "Ebot: -15\n", + "Eset: 50\n", + "Esetp: 1000\n", + "================================================================================================================================================================\n" + ] + } + ], "source": [ "# this cell mimicks an input file\n", - "fdf = sisl.get_sile(\"./data/lat3_791/Fe3GeTe2.fdf\") # ./Jij_for_Marci_6p45ang/CrBr.fdf\n", + "fdf = sisl.get_sile(\n", + " \"/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf\"\n", + ") # ./Jij_for_Marci_6p45ang/CrBr.fdf\n", "# this information needs to be given at the input!!\n", "scf_xcf_orientation = np.array([0, 0, 1]) # z\n", "# list of reference directions for around which we calculate the derivatives\n", @@ -57,7 +107,7 @@ "# human readable definition of magnetic entities ./lat3_791/Fe3GeTe2.fdf\n", "magnetic_entities = [\n", " dict(atom=3, l=2),\n", - " # dict(atom=4, l=2),\n", + " dict(atom=4, l=2),\n", " dict(atom=5, l=2),\n", " # dict(atom=[3, 4],),\n", "]\n", @@ -65,23 +115,11 @@ "pairs = [\n", " # isotropic should be -82 meV\n", " dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])),\n", - " dict(ai=0, aj=1, Ruc=np.array([1, 0, 0])),\n", - " dict(ai=0, aj=1, Ruc=np.array([-1, 0, 0])),\n", - " dict(ai=0, aj=1, Ruc=np.array([2, 0, 0])),\n", - " dict(ai=0, aj=1, Ruc=np.array([-2, 0, 0])),\n", - " dict(ai=0, aj=1, Ruc=np.array([0, 1, 0])),\n", - " dict(ai=0, aj=1, Ruc=np.array([0, -1, 0])),\n", - " dict(ai=0, aj=1, Ruc=np.array([0, 2, 0])),\n", - " dict(ai=0, aj=1, Ruc=np.array([0, -2, 0])),\n", - " dict(ai=1, aj=0, Ruc=np.array([0, 0, 0])),\n", - " dict(ai=1, aj=0, Ruc=np.array([1, 0, 0])),\n", - " dict(ai=1, aj=0, Ruc=np.array([-1, 0, 0])),\n", - " dict(ai=1, aj=0, Ruc=np.array([2, 0, 0])),\n", - " dict(ai=1, aj=0, Ruc=np.array([-2, 0, 0])),\n", - " dict(ai=1, aj=0, Ruc=np.array([0, 1, 0])),\n", - " dict(ai=1, aj=0, Ruc=np.array([0, -1, 0])),\n", - " dict(ai=1, aj=0, Ruc=np.array([0, 2, 0])),\n", - " dict(ai=1, aj=0, Ruc=np.array([0, -2, 0])),\n", + " # dict(ai=1, aj=0, Ruc=np.array([0, 0, 0])),\n", + " dict(ai=1, aj=2, Ruc=np.array([0, 0, 0])),\n", + " # dict(ai=2, aj=1, Ruc=np.array([0, 0, 0])),\n", + " dict(ai=0, aj=2, Ruc=np.array([0, 0, 0])),\n", + " # dict(ai=2, aj=0, Ruc=np.array([0, 0, 0])),\n", " # these should all be around -41.9 in the isotropic part\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", @@ -113,11 +151,11 @@ "\"\"\"\n", "\n", "# Brilloun zone sampling and Green function contour integral\n", - "kset = 20\n", + "kset = 10\n", "kdirs = \"xy\"\n", - "ebot = -30\n", + "ebot = -15\n", "eset = 50\n", - "esetp = 10000\n", + "esetp = 1000\n", "\n", "\n", "# MPI parameters\n", @@ -125,8 +163,6 @@ "size = comm.Get_size()\n", "rank = comm.Get_rank()\n", "root_node = 0\n", - "if rank == root_node:\n", - " print(\"Number of nodes in the parallel cluster: \", size)\n", "\n", "simulation_parameters = dict(\n", " path=\"Not yet specified.\",\n", @@ -143,20 +179,62 @@ "# digestion of the input\n", "# read in hamiltonian\n", "dh = fdf.read_hamiltonian()\n", - "try:\n", - " simulation_parameters[\"geom\"] = fdf.read_geometry()\n", - "except:\n", - " print(\"Error reading geometry.\")\n", + "simulation_parameters[\"cell\"] = fdf.read_geometry().cell\n", "\n", "# unit cell index\n", "uc_in_sc_idx = dh.lattice.sc_index([0, 0, 0])\n", "\n", + "print_parameters(simulation_parameters)\n", "times[\"setup_time\"] = timer()" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/danielpozsar/Documents/oktatás/elte/phd/grogu_project/.venv/lib/python3.9/site-packages/matplotlib/cbook.py:1762: ComplexWarning: Casting complex values to real discards the imaginary part\n", + " return math.isfinite(val)\n", + "/Users/danielpozsar/Documents/oktatás/elte/phd/grogu_project/.venv/lib/python3.9/site-packages/matplotlib/cbook.py:1398: ComplexWarning: Casting complex values to real discards the imaginary part\n", + " return np.asarray(x, float)\n" + ] + }, + { + "data": { + "text/plain": [ + "-12.806739" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGiCAYAAAAba+fDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA3T0lEQVR4nO3deXhU9aH/8c9km5BtQvYEEkgAQfZFjHErKrJorSi1bq1g3YWqoK3SX0XBalx6tddeqn1uRbRWqfa6VNuqiCxawhaIC0sEDCRAFkhIJguZJDPn9wdmdCRAApk5M5P363nmeZgzZ5JPjieZj9/zPedYDMMwBAAAECRCzA4AAADQnSg3AAAgqFBuAABAUKHcAACAoEK5AQAAQYVyAwAAggrlBgAABBXKDQAACCqUGwAAEFQoNwAAIKh4tdzk5+dr/Pjxio2NVUpKiqZNm6bi4mKPdZqbmzVr1iwlJiYqJiZG06dPV2VlpTdjAQCAIObVcrNq1SrNmjVLa9eu1bJly9Ta2qpJkyapsbHRvc6cOXP07rvv6o033tCqVau0f/9+XXnlld6MBQAAgpjFlzfOPHDggFJSUrRq1Sqdf/75qqurU3Jysl599VX9+Mc/liRt375dp59+ugoKCnTWWWf5KhoAAAgSYb78ZnV1dZKkhIQESVJhYaFaW1s1ceJE9zpDhgxRVlbWMcuNw+GQw+FwP3e5XKqpqVFiYqIsFouXfwIAANAdDMNQfX29MjIyFBLSvQeSfFZuXC6X7rnnHp1zzjkaPny4JKmiokIRERGKj4/3WDc1NVUVFRUdfp38/HwtWLDA23EBAIAPlJWVqW/fvt36NX1WbmbNmqUvv/xSn3766Sl9nXnz5mnu3Lnu53V1dcrKylJZWZni4uJONSYAdAuXy1De48vV6HDqrTvP1qDUWLMjAX7FbrcrMzNTsbHd/7vhk3Ize/Zsvffee1q9erVHO0tLS1NLS4tqa2s9Rm8qKyuVlpbW4deyWq2yWq1HLY+Li6PcAPAbpdVNOiyrIqNCNDInXeGhXHkD6Ig3ppR49bfNMAzNnj1bb731lj7++GNlZ2d7vD5u3DiFh4dr+fLl7mXFxcUqLS1VXl6eN6MBgFdtLbdLkgalxlBsAB/z6sjNrFmz9Oqrr+qdd95RbGysex6NzWZTr169ZLPZdNNNN2nu3LlKSEhQXFycfvGLXygvL48zpQAEtG3flJvT0xlRBnzNq+XmueeekyRNmDDBY/mLL76omTNnSpKeeeYZhYSEaPr06XI4HJo8ebL++Mc/ejMWAHiVo82p9SU1kig3gBl8ep0bb7Db7bLZbKqrq2PODQBTtTpd+tcX5frdh8UqqzksSfrH7HM0sm+8ucEAP+TNz2+fXucGAIJNq9Ol9z7frw+3VOrTHQdV72iTJKXEWvXA1CEUG8AElBsAOAkul6F3P9+vZ5Z9pd3VTe7lidERmnF2f918XraiIvgTC5iB3zwAOAm3v1KoD7ceuclvYnSErj+rny4ckqKRfWwKCeFq6YCZKDcA0EWtTpc+3l4lSbr7okG69fwcRVv5cwr4C34bAaCL9tceVpvLUERYiO6+aBAjNYCf4cpSANBFe76ZY5OVEEWxAfwQ5QYAumhPzZFy0z8xyuQkADpCuQGALiqtbpQkZSVEm5wEQEcoNwDQRe2HpfoxcgP4JcoNAHRR6TeHpbIoN4BfotwAQBcYhvHtyE0C5QbwR5QbAOiCA/UOHW51KsQi9e1NuQH8EeUGALqg/UypjPheigjjTyjgj/jNBIAuYDIx4P8oNwDQBZwGDvg/yg0AdEH7YSlGbgD/RbkBgC7gTCnA/1FuAKAL9nxzWKpfIoelAH9FuQGATrI3t+pQU6skLuAH+DPKDQB0Uuk3h6SSYiIUYw0zOQ2AY6HcAEAntc+3yWK+DeDXKDcA0El7aphvAwQCyg0AdEJVfbPe2rRPEqeBA/6Og8YAcAJlNU366QvrtKe6ScmxVv14XF+zIwE4DsoNABzH3kNNuur5AlXYm5WZ0Et/veksbpgJ+DnKDQAcQ5vTpXuWFqnC3qxBKTF65eZcpcZFmh0LwAkw5wYAjmHRil3auOeQYqxhWjxzPMUGCBCUGwDoQOGeQ3r24x2SpEemDVMmp38DAYNyAwDfU1HXrLte2yyny9DlozN0xRgmEAOBhHIDAN9R3eDQ9X9eq321h9U/MUqPTBtudiQAXUS5AYBv1B1u1c9eWK9dBxqVbovUKzfnKi4y3OxYALqIcgMAkgzD0F2vbdbWcruSYiL015tzOeUbCFCUGwCQtHRDmVZ9dUDWsBC9/PNc5STHmB0JwEmi3ADo8fYeatJv39sqSfrl5MEamhFnciIAp4JyA6BHc7kM/ervn6uxxakz+vXWjedkmx0JwCmi3ADo0f5euFdrdlUrMjxET101SqEhFrMjAThFlBsAPVaDo01PflAsSZp78WnKToo2ORGA7kC5AdBjPbdypw42ONQ/MUozz+ZwFBAsvFpuVq9ercsuu0wZGRmyWCx6++23PV6fOXOmLBaLx2PKlCnejAQAko5MIv7fT0okSfMuOV0RYfy/HhAsvPrb3NjYqFGjRmnRokXHXGfKlCkqLy93P1577TVvRgIASdKT7xerpc2l3OwETRqaanYcAN0ozJtffOrUqZo6depx17FarUpLS/NmDADwsOtAg/7x2X5ZLNKDPxwqi4VJxEAwMX0cduXKlUpJSdHgwYN1xx13qLq6+rjrOxwO2e12jwcAdMUbG/dKki4cnKLhfWwmpwHQ3UwtN1OmTNHLL7+s5cuX64knntCqVas0depUOZ3OY74nPz9fNpvN/cjMzPRhYgCBrs3p0pubjpSbq87gbt9AMPLqYakTueaaa9z/HjFihEaOHKkBAwZo5cqVuuiiizp8z7x58zR37lz3c7vdTsEB0Gmf7DioqnqHEqIjdOEQ5toAwcj0w1LflZOTo6SkJO3cufOY61itVsXFxXk8AKCz3igskyRdPjqDM6SAIOVXv9l79+5VdXW10tPTzY4CIAjVNLZo2dZKSdJV4xjxBYKVVw9LNTQ0eIzClJSUqKioSAkJCUpISNCCBQs0ffp0paWladeuXfrVr36lgQMHavLkyd6MBaCH+r/CvWp1GhreJ46bYwJBzKvlZuPGjbrgggvcz9vnysyYMUPPPfecPv/8c7300kuqra1VRkaGJk2apEceeURWq9WbsQAEuQZHmz7dcVBV9c1qbnWqyu7Qyq8OaGdVgyRGbYBg59VyM2HCBBmGcczXP/jgA29+ewBBrL65VfXNbTIkHW5xqrSmUSUHm1Swq1qrdxxQS5vrqPeEhVg0eViafnIG5QYIZqaeLQUAXVVlb9Z/L9+hpRvK5HQd+3+e+idGaUhanKIiQhVtDdOZ2Qn6weBkxUWG+zAtADNQbgD4tS/31emVtXvU0uZSi9Ol5duqdLj1yLWwIkJDZLFI4aEhykyIUnZSlE5Pi9OkYWk6LTWGKw8DPRTlBoBfW/DuFm3Yfchj2ZiseD0wZYhycxJNSgXAn1FuAPitwy1OFZXVSpLmTDxN0dZQDUiJ0YTTkhmVAXBMlBsAfmtT6SG1Og1l2CJ110UDKTQAOsWvLuIHAN+17usjN9LNzUmk2ADoNMoNAL+1tqRGkpSbnWByEgCBhHIDwC81t3473+ZMyg2ALqDcAPBLRWW1amlzKTnWquykaLPjAAgglBsAfmnd198ekmK+DYCuoNwA8EvrSr6dTAwAXUG5AeB3Wtpc2lR65MJ9ZzHfBkAXUW4A+J2Ne2rU3OpSYnSEBqbEmB0HQICh3ADwK+9/Wa7b/lIoSTpnYBLzbQB0GVcoBmCaspomrSyuUnFlvVrbDFU3OvTRtipJ0ujMeP36ktNNTgggEFFuAPhUS5tLf9tQqiVrdmvXgcYO17ntBzm6b9JghYcyuAyg6yg3AHziQL1DH26t0B9X7NK+2sOSpNAQi8b1663x/XurV3iowkJDNL5/b43rxyRiACePcgPAqz7YUqE/fLxDX+6zu5elxFo1+8KBmjamj+Iiw01MByAYUW4AeE1tU4t+8epmtThdkqThfeI0bXQfXZ/bT70iQk1OByBYUW4AeM07RfvV4nRpSFqs/nJTrpJjrWZHAtADMFsPgNe8vrFMknTN+EyKDQCfodwA8Iot++u0Zb9dEaEhunx0H7PjAOhBKDcAvOKNjXslSRcPTVXv6AiT0wDoSSg3ALqdo82pt4v2SZKuOqOvyWkA9DSUGwDd7t9fVKi2qVVpcZE6b1Cy2XEA9DCcLQWg27hchhb/p0RPvL9d0pFRm9AQ7g0FwLcoNwC6RXOrU3e8UqgVxQckSZOHpeqOCQNMTgWgJ6LcAOgWS9eXakXxAUWEhWj+D4fq+tws7ugNwBSUGwCnzDAM/e2bs6PmTR2in57Vz+REAHoyJhQDOGVb9tu1rfzINW2uGMM1bQCYi3ID4JS1X4l40rBUxUdxTRsA5qLcADglza1Ovb35yDVtfnJGpslpAIByA+AUfbClQvbmNvWJ76VzBiaZHQcAmFAMoOsOtzhVuOeQSmua9Je1eyRJ08dxTRsA/oFyA6DTGh1t+svaPfrzJ1/rYEOLe3mIRbpqHLdZAOAfKDcAjuuDLRVa/dUB7TrQoC377apvbpMkpcVFamhGnDJ799J5g5KVmRBlclIAOIJyA+CYKu3Nuv2VQhnGt8tykqI164KBunx0hsJCmbYHwP9QbgAcU8nBRhmGlBRj1a8vGaKBKTEalmFjbg0Av+bV/+1avXq1LrvsMmVkZMhisejtt9/2eN0wDM2fP1/p6enq1auXJk6cqB07dngzEoAu2HfosCRpcFqMrhzbVyP7xlNsAPg9r5abxsZGjRo1SosWLerw9SeffFLPPvusnn/+ea1bt07R0dGaPHmympubvRkLQCft/abc9I1nPg2AwOHVw1JTp07V1KlTO3zNMAz9/ve/129+8xtdfvnlkqSXX35Zqampevvtt3XNNdd4MxqATthX2yRJ6tO7l8lJAKDzTJsNWFJSooqKCk2cONG9zGazKTc3VwUFBcd8n8PhkN1u93gA8I59tUdGbvrEU24ABA7Tyk1FRYUkKTU11WN5amqq+7WO5Ofny2azuR+ZmVzuHfAW92EpRm4ABJCAO49z3rx5qqurcz/KysrMjgQEJZfLUHntkflvHJYCEEhMKzdpaWmSpMrKSo/llZWV7tc6YrVaFRcX5/EA0P0ONDjU4nQpNMSitLhIs+MAQKeZVm6ys7OVlpam5cuXu5fZ7XatW7dOeXl5ZsUC8I32Q1JpcZFcrA9AQPHq2VINDQ3auXOn+3lJSYmKioqUkJCgrKws3XPPPfrtb3+rQYMGKTs7Ww8++KAyMjI0bdo0b8YC0Al7D3GmFIDA5NVys3HjRl1wwQXu53PnzpUkzZgxQ0uWLNGvfvUrNTY26tZbb1Vtba3OPfdcvf/++4qMZAgcMFv7mVJ9OVMKQIDxarmZMGGCjO/elOZ7LBaLFi5cqIULF3ozBoCT0H51YkZuAAQaDqQD6BDXuAEQqCg3ADr07TVuuPUCgMBCuQFwFMMwOCwFIGBRbgAc5VBTqw63OiVJ6TYm+AMILJQbAEdpH7VJjrUqMjzU5DQA0DWUGwBHab/GDfeUAhCIKDcAjsKZUgACGeUGwFH2MpkYQACj3AA4ivuwFCM3AAKQV69QDCAw1Da1aE91kz7fV6e/F+7VZ2W1kqS+CVzjBkDgodwAPZRhGHqnaL9+92Gx+zBUu7AQi6aOSFdeTqJJ6QDg5FFugB6ktqlFew8dVmlNk174tESFew65X0uJtap/UrQmDU3VtDF9lBRjNTEpAJw8yg0QxFqdLv35kxJt2F2jL/fVqare4fF6r/BQzb5woGac3V8xVv4cAAgO/DUDgthbm/fpife3eyxLirGqT+9eGp4Rp9kXDlS6jUnDAIIL5QYIYu8U7ZMkXTmmj64/K0tD0uIUzQgNgCDHXzkgSFXam7VmV7Ukac7FpymTM58A9BBc5wYIUu9+tl+GIY3r15tiA6BHodwAQeqdov2SpMtHZ5icBAB8i3IDBKFdBxr0xb46hYZYdOmIdLPjAIBPUW6AIPSPb0ZtzhuUpESuVwOgh2FCMRAEDMPQ3kOHtan0kDbtOaR/fHak3Ewb3cfkZADge5QbIIBtLj2k//3ka60vOaSDDZ4X6EuIjtDFQ1NNSgYA5qHcAAGorKZJT35QrHe/GaGRpPBQi4Zm2DQ2K15js3rr7AGJXNMGQI/EXz4gQLhchv6z66BeWbtHH22rktNlyGKRpo/tq2vGZ2p4H5siw0PNjgkApqPcAAFgy/46/frNL/TZ3jr3snMHJmneJUM0LMNmYjIA8D+UG8CPNbc69cxHX+nPn5TI6TIUHRGq6eP66vrcfhqcFmt2PADwS5QbwE+1OV2a/eomfbStSpJ06Yh0PXTZUKXERZqcDAD8G+UG8EOGYWj+P7boo21VsoaF6Nlrx2jysDSzYwFAQKDcAH6mzenSohW79Oq6Ulks0n9fQ7EBgK6g3AB+wOUy9MnOg/pH0X4t316p2qZWSdLDlw3TlOEUGwDoCsoNYKL65la993m5Xvi0RDurGtzLe0eF6/YfDNCMs/ubFw4AAhTlBvAxp8vQ3wvL9M8vKlSw66BanYYkKcYapulj++iSEeka16+3wkK59RsAnAzKDeBDza1O3b10sz7YUulelpMUretys3T1+EzFRoabmA4AggPlBvCR2qYW3fzSRm3cc0gRoSGafeFAXTIiXQNTYsyOBgBBhXIDeFl53WG9vmGvXltfqgp7s2Ijw/S/N5yhs3ISzY4GAEGJcgN40csFu/XwP7bIdWRajfrE99LimeO5ujAAeBHlBvCSRkebnvqgWC5DOrN/gq7LzdKU4Wnc3BIAvIxyA3jJW5v3qb65TdlJ0Vp661kKCbGYHQkAegTONQW8wDAMvVywW5L0s7P6UWwAwIdMLzcPP/ywLBaLx2PIkCFmxwJOScHX1fqqskFREaH68Rl9zY4DAD2KXxyWGjZsmD766CP387Awv4gFnLSX1uyWJF05to/iuHYNAPiUX7SIsLAwpaV17v45DodDDofD/dxut3srFtBpza1ObdpzSKU1Taqqd2jZ1iMX6ZuR19/cYADQA/lFudmxY4cyMjIUGRmpvLw85efnKysrq8N18/PztWDBAh8nBDwdamxRUVmtNpfVakNJjQpLD6mlzeWxztkDEjUolVO+AcDXLIZhGGYG+Pe//62GhgYNHjxY5eXlWrBggfbt26cvv/xSsbFHfzB0NHKTmZmpuro6xcXF+TI6eqgPt1Tojr9uktPl+auTGmfVsAybkmIilBIbqavHZyozIcqklADg3+x2u2w2m1c+v00fuZk6dar73yNHjlRubq769eun119/XTfddNNR61utVlmtVl9GBDys3nFATpeh1DirzhmQpDFZ8cobkKQBydGyWDgrCgDMZnq5+b74+Hiddtpp2rlzp9lRgA5V1B0ZOfzFhYP007P6mZwGAPB9pp8K/n0NDQ3atWuX0tPTzY4CdKjS3ixJSouLNDkJAKAjppeb++67T6tWrdLu3bu1Zs0aXXHFFQoNDdW1115rdjSgQxXt5cZGuQEAf2T6Yam9e/fq2muvVXV1tZKTk3Xuuedq7dq1Sk5ONjsacJRWp0sHG44clkpl5AYA/JLp5Wbp0qVmRwA6rareIcOQwkMtSoyOMDsOAKADph+WAgJJRd2RQ1IpsZHcLwoA/BTlBuiCSubbAIDfo9wAXdA+csOZUgDgvyg3QBe0j9wwmRgA/BflBuiC8vaRGxtXyQYAf0W5AbqggpEbAPB7lBugC9oPS6XbepmcBABwLJQboJMMw2BCMQAEAMoN0El1h1vlaHNJklLimHMDAP6KcgN0Uvt8m95R4YoMDzU5DQDgWCg3QCe1nynFZGIA8G+UG6CTKuu4OjEABALKDdBJFe4zpSg3AODPKDdAJ3F1YgAIDJQboJM4DRwAAgPlBuikCrtDkpTKYSkA8GuUG6CTKuoOS2LkBgD8HeUG6ITmVqcONbVKotwAgL8LMzsA4M8cbU4VldZqc1mtJMkaFqL4qHBzQwEAjotyAxzDocYWXfWnAu2sanAv658YLYvFYmIqAMCJUG6ADrS0uXT7K4XaWdUgW69wjc2K18CUGF0xpq/Z0QAAJ0C5Ab7HMAz95u0vtK6kRjHWML1+W54Gp8WaHQsA0EmUG/R4DY42ldce1t5Dh7VxT43+s7NaRWW1CrFIf7huDMUGAAIM5QY9UpvTpX98tl/PrdylHd+ZU9POYpEW/GiYLhicYkI6AMCpoNygxyiuqNdne2u1o7JeH2ypVGlNk/u12MgwZdh6aXgfm3JzEpSXk6jMhCgT0wIAThblBj3CH1fu1JPvF3ssS4iO0M3nZeu6M7MUHxVhUjIAQHej3CDovVO0z11scrMTdHp6nIZmxOmHI9MVFcGvAAAEG/6yI6ht2F2jX77xuSTppnOz9eAPh5qcCADgbdx+AUHr/S8rdNOSDWpxujR5WKp+fcnpZkcCAPgAIzcIOo2ONv32n1v12voySdIZ/Xrr91ePUWgIVxYGgJ6AcoOA5XIZsje3qrapVYeaWlRa06QPt1ZqxfYqNbU4ZbFIt50/QHMvPk0RYQxSAkBPQblBQFj11QHl/2ubGlva5HQaamp1qu5wqwyj4/X7JUYp/8oROntAkm+DAgBMR7lBQHh+5S5tr6jv8LXoiFDFR0UoITpCZw9M1CXD0zWyr40bXAJAD0W5gd9ramnTxj01kqQ//Wyc0m2RigwPVXxUuGy9wmUNCzU5IQDAn1Bu4PfWldSo1WmoT3wvTRqayogMAOC4mGUJv/fpjoOSpHMHJlFsAAAnRLmB32svN+edxuRgAMCJUW7g16rszSqurJfFIp3DmU8AgE7wi3KzaNEi9e/fX5GRkcrNzdX69evNjgQ/8ck3ozbDM2zqHc3NLQEAJ2Z6ufnb3/6muXPn6qGHHtKmTZs0atQoTZ48WVVVVWZHgx/4dOc3820GMWoDAOgc08vN008/rVtuuUU33nijhg4dqueff15RUVFavHix2dFgMsMw3OXmvIGUGwBA55hablpaWlRYWKiJEye6l4WEhGjixIkqKCjo8D0Oh0N2u93jgeC0cc8hHah3KDI8ROP69zY7DgAgQJhabg4ePCin06nU1FSP5ampqaqoqOjwPfn5+bLZbO5HZmamL6LCx9aX1OjnSzZIki4cksKF+gAAnWb6Yamumjdvnurq6tyPsrIysyOhm73/ZYV+9sI61Te3aXz/3sq/YqTZkQAAAcTUKxQnJSUpNDRUlZWVHssrKyuVlpbW4XusVqusVqsv4sHHquqb9eg/t+mdov2SpImnp+p/rhujyHBGbQAAnWdquYmIiNC4ceO0fPlyTZs2TZLkcrm0fPlyzZ4928xo8JED9Q5t2F2jdV9X683N+1Tf3CaLRfr5OdmaN3WIwkIDbnARAGAy0+8tNXfuXM2YMUNnnHGGzjzzTP3+979XY2OjbrzxRrOjoZu1tLm0o6pe28rrVbjnkNaVVOvrA40e64zoY9OjVwzXyL7x5oQEAAQ808vN1VdfrQMHDmj+/PmqqKjQ6NGj9f777x81yRiBq6q+WY+8t03//qJcbS7jqNeHpMXqzOwEnT0gURcPTVNoCPePAgCcPIthGEd/2gQQu90um82muro6xcXFmR2nxzMMQ9WNLWp0tKnV6VLhnkN67F/bVXe4VZIUFxmmIelxGtnHptycRI3v31vxUVx5GAB6Gm9+fps+coPA19zq1O8+KNby7VXaV3tYLW2uo9YZ3idOj10xQiP62LizNwDAqyg3OCW7Dzbqzr9u0tZyz4spRkWEKiIsRDHWMP30rH66+dxsJgcDAHyCcoOTtuqrA5r1101qcLQpITpCCy8fplF945UaF6mIMIoMAMAclBuclKKyWt3+l0IdbnXqzP4JevbaMUqzRZodCwAAyg26bvfBRv18yQYdbnXq/NOS9cKMMxTOIScAgJ/gEwldUrinRjNeXK+axhaN6GPTc9ePpdgAAPwKIzc4oTanSxt2H9Lzq3Zp1VcHJEmZCb20eOZ4RVvZhQAA/oVPJnSopc2lj7dX6p9fVGj1Vwfc16kJDbHoqnF9Nefi05Qcyz2+AAD+h3IDD4caW/SHj3fqrc17daip1b08PipcU4al6Y4JA9QvMdrEhAAAHB/lBm4bdtfortc2q7yuWZKUEmvVFWP6aNKwVI3O7M1tEQAAAYFyAxmGoedXfa3ffVgsp8tQTlK0fvPD03X+oGQuvAcACDiUmx6u0dGmX/79M/3riwpJ0hVj+ui304YzURgAELD4BOvBiivqdffSzdpeUa/wUIsWXj5c14zP5N5PAICARrnpARxtTn3y1UFV2JvlaHPpQL1Dy7dVakdVgyQpKcaqP/1srMb1SzA5KQAAp45yE8T21x7Wi/8p0f9t2qeaxpajXg8PtegHp6XokWnDlG7rZUJCAAC6H+UmSDU62nTlH9eown7kzKfUOKtG9Y1XZHiooq2hys1O1AVDUmTrFW5yUgAAuhflJkj9+ZMSVdib1Se+lxb8aJgmDObMJwBAz0C5CUIH6h360+pdkqRfX3K6Jg5NNTkRAAC+w//KB6H/Xv6VmlqcGpUZr0tGpJkdBwAAn2LkJkg0tbRp76HD2lXVoNfWl0mS5k0dwmndAIAeh3ITBLaV23XlH9focKvTvWzi6Sk6KyfRxFQAAJiDchMEXi7YrcOtTkVHhCozIUpZCVF68IdDzY4FAIApKDcBrqmlTe9+Vi5JemHmeEZrAAA9HhOKA9y/v6hQg6NNWQlRys3mCsMAAFBuAtwbhUcmD181ri+ThwEAEOUmoO2pbtTar2tksUjTx/U1Ow4AAH6BchPA/l64V5J03qBkZcRzbygAACQmFAcUp8vQ3wvL9M8vKrS/9rD2VDdKOnJICgAAHEG5CQCGYeiTHQf12L+2aXtFvcdr2UnRupjbKwAA4Ea58VMtbS5tLbfrwy0V+ucX5dpT3SRJiosM020/GKDRmfFKt0UqMyFK4dwQEwAAN8qNn1mxvUq/X75D2/bb1eJ0uZdHhofomvFZuvuiQeodHWFiQgAA/Bvlxo/UNbVq9qub1Nhy5DYKtl7hystJ1KUj03XhkBRFW/nPBQDAifBp6Uf+sna3GlucOi01Rv97wxnKSoji2jUAAHQR5cZPNLc69eJ/dkuS7pwwUP0So80NBABAgGImqp94Y2OZqhtb1Ld3L/1wZLrZcQAACFiUGz/Q5nTpT6u/liTden6Owjj7CQCAk8anqB9Ysma39h46rMToCF01LtPsOAAABDTm3JjIMAwtWrFTv/vwK0nSbT/IUa+IUJNTAQAQ2Ewduenfv78sFovH4/HHHzczktfZm1v1/pcVenVdqea+/pm72My6YIBuOS/H5HQAAAQ+00duFi5cqFtuucX9PDY21sQ03nWosUWXPvuJ9tc1u5dZLNJDPxyqmedkm5gMAIDgYXq5iY2NVVpamtkxvM4wDP2/t7/Q/rpmJcdaNaqvTQnREfrhyAydf1qy2fEAAAgaFsMwDLO+ef/+/dXc3KzW1lZlZWXpuuuu05w5cxQWduzO5XA45HA43M/tdrsyMzNVV1enuLg4X8Q+KW9u2qu5r3+msBCL3rrzHI3oazM7EgAAprHb7bLZbF75/DZ15Oauu+7S2LFjlZCQoDVr1mjevHkqLy/X008/fcz35Ofna8GCBT5MeerKapr00DtbJEn3TBxEsQEAwIu6feTmgQce0BNPPHHcdbZt26YhQ4YctXzx4sW67bbb1NDQIKvV2uF7A23kZu3X1brrtc2qqndoXL/e+tutZ3EdGwBAj+fNkZtuLzcHDhxQdXX1cdfJyclRRMTRd7besmWLhg8fru3bt2vw4MGd+n7e3DinatGKnfqvD4vlMqSBKTFacuN49e0dZXYsAABMF1CHpZKTk5WcfHITZIuKihQSEqKUlJRuTuV7G3fX6KkPiiVJV47to99OG66oCNPnbwMAEPRM+7QtKCjQunXrdMEFFyg2NlYFBQWaM2eOfvrTn6p3795mxeo2u6ubJEl5OYl6+iejzQ0DAEAPYlq5sVqtWrp0qR5++GE5HA5lZ2drzpw5mjt3rlmRulVDc6skKSHm6MNvAADAe0wrN2PHjtXatWvN+vZe1+BokyTFWjkUBQCAL3HajpfUf1NuYig3AAD4FOXGSxqavyk3kZQbAAB8iXLjJY2M3AAAYArKjZe459wwcgMAgE9Rbrykvv2wlDXc5CQAAPQslBsvaR+5Yc4NAAC+Rbnxkgbm3AAAYArKjZe4z5ai3AAA4FOUGy+p57AUAACmoNx4gaPNqZY2lyRGbgAA8DXKjRc0Opzuf1NuAADwLcqNF7TPt4mKCFVoiMXkNAAA9CyUGy+odxy5IzijNgAA+B7lxgvaD0sxmRgAAN+j3HhBwzcjN7GM3AAA4HOUGy+o547gAACYhnLjBVydGAAA81BuvKD9bKloyg0AAD5HufGC9pEb5twAAOB7lBsvYM4NAADmodx4wbdzbsJNTgIAQM9DufGCBkZuAAAwDeXGC5hzAwCAeSg3XsCp4AAAmIdy4wXucsNhKQAAfI5y4wXuOTeM3AAA4HOUGy/gsBQAAOah3HQzl8vgsBQAACai3HSzxpY2978ZuQEAwPcoN92sfdQmPNQiaxibFwAAX+PTt5t9dzKxxWIxOQ0AAD0P5aab1TPfBgAAU1Fuutm3IzfcVwoAADNQbroZt14AAMBclJtuxmngAACYi3LTzbg6MQAA5qLcdLP2kZtoyg0AAKag3HQz95wbDksBAGAKr5WbRx99VGeffbaioqIUHx/f4TqlpaW69NJLFRUVpZSUFP3yl79UW1tbh+sGinoOSwEAYCqvfQK3tLToqquuUl5enl544YWjXnc6nbr00kuVlpamNWvWqLy8XDfccIPCw8P12GOPeSuW13HTTAAAzOW1kZsFCxZozpw5GjFiRIevf/jhh9q6dateeeUVjR49WlOnTtUjjzyiRYsWqaWlxVuxvK6huVUSZ0sBAGAW0+bcFBQUaMSIEUpNTXUvmzx5sux2u7Zs2XLM9zkcDtntdo+HP+E6NwAAmMu0clNRUeFRbCS5n1dUVBzzffn5+bLZbO5HZmamV3N2lXvODSM3AACYokvl5oEHHpDFYjnuY/v27d7KKkmaN2+e6urq3I+ysjKvfr+uYs4NAADm6tIn8L333quZM2ced52cnJxOfa20tDStX7/eY1llZaX7tWOxWq2yWq2d+h5m4FRwAADM1aVP4OTkZCUnJ3fLN87Ly9Ojjz6qqqoqpaSkSJKWLVumuLg4DR06tFu+hy852pxaX1LjvkIxF/EDAMAcXvsELi0tVU1NjUpLS+V0OlVUVCRJGjhwoGJiYjRp0iQNHTpUP/vZz/Tkk0+qoqJCv/nNbzRr1iy/HpnpyOsbyvTwu1vU1OKUJEWEhah3VITJqQAA6Jm8Vm7mz5+vl156yf18zJgxkqQVK1ZowoQJCg0N1Xvvvac77rhDeXl5io6O1owZM7Rw4UJvRfKal9fuVlOLU8mxVl0wOFk/HpepyPBQs2MBANAjWQzDMMwOcSrsdrtsNpvq6uoUFxfn8+/vchka+tD7am51acV9E5SdFO3zDAAABBpvfn5zb6lTVHaoSc2tLlnDQpSVEGV2HAAAejzKzSkqrqiXJA1IjlFoiMXkNAAAgHJzinZUNUiSTkuNMTkJAACQKDen7KvKIyM3g1JjTU4CAAAkys0p+6qyfeSGcgMAgD+g3JwCp8vQrgNHys1gyg0AAH6BcnMK9lQ3qqXNpV7hoerbu5fZcQAAgCg3p6R9vs3AlBiFcKYUAAB+gXJzCtrn2wziTCkAAPwG5eYUtI/cMJkYAAD/Qbk5BTsqmUwMAIC/odycpFanS18f5LAUAAD+hnJzkvZUN6rVaSg6IlR94jlTCgAAfxFmdoBAc6ixRSu/qtK7n5VLkgamxspi4UwpAAD8BeWmCyrtzbr02U91sMHhXnZWToKJiQAAwPdRbrrgfz7eqYMNDmXYIvWj0X107sAknT0g0exYAADgOyg3nVRW06SlG0olSf/1k9HKo9QAAOCXmFDcSX/4eIdanYbOGZhIsQEAwI9Rbjrh6wMN+r9N+yRJ904abHIaAABwPByWOobNpYe0Zle1tlfUq3B3jZwuQxcNSdHYrN5mRwMAAMdBuTmGdz8r1+L/lLifx0WG6b7JjNoAAODvKDfHkDcgUQcbHDo9PU5D0mM1um+8ekdHmB0LAACcAOXmGC4emqqLh6aaHQMAAHQRE4oBAEBQodwAAICgQrkBAABBhXIDAACCCuUGAAAEFcoNAAAIKpQbAAAQVCg3AAAgqFBuAABAUKHcAACAoEK5AQAAQYVyAwAAggrlBgAABBXKDQAACCqUGwAAEFS8Vm4effRRnX322YqKilJ8fHyH61gslqMeS5cu9VYkAADQA4R56wu3tLToqquuUl5enl544YVjrvfiiy9qypQp7ufHKkIAAACd4bVys2DBAknSkiVLjrtefHy80tLSvBUDAAD0MKbPuZk1a5aSkpJ05plnavHixTIM47jrOxwO2e12jwcAAEA7r43cdMbChQt14YUXKioqSh9++KHuvPNONTQ06K677jrme/Lz892jQgAAAN9nMU40VPIdDzzwgJ544onjrrNt2zYNGTLE/XzJkiW65557VFtbe8KvP3/+fL344osqKys75joOh0MOh8P93G63KzMzU3V1dYqLizvxDwEAAExnt9tls9m88vndpZGbe++9VzNnzjzuOjk5OScdJjc3V4888ogcDoesVmuH61it1mO+BgAA0KVyk5ycrOTkZG9lUVFRkXr37k15AQAAJ81rc25KS0tVU1Oj0tJSOZ1OFRUVSZIGDhyomJgYvfvuu6qsrNRZZ52lyMhILVu2TI899pjuu+8+b0UCAAA9gNfKzfz58/XSSy+5n48ZM0aStGLFCk2YMEHh4eFatGiR5syZI8MwNHDgQD399NO65ZZbvBUJAAD0AF2aUOyPvDkhCQAAeIc3P79Nv84NAABAd6LcAACAoEK5AQAAQYVyAwAAggrlBgAABBXKDQAACCqUGwAAEFQoNwAAIKhQbgAAQFCh3AAAgKBCuQEAAEGFcgMAAIIK5QYAAAQVyg0AAAgqlBsAABBUKDcAACCoUG4AAEBQCTM7wKkyDEOSZLfbTU4CAAA6q/1zu/1zvDsFfLmpr6+XJGVmZpqcBAAAdFV9fb1sNlu3fk2L4Y3K5EMul0v79+9XbGysLBZLt35tu92uzMxMlZWVKS4urlu/dqBhW3hie3yLbeGJ7fEttoUntse32rfF1q1bNXjwYIWEdO8smYAfuQkJCVHfvn29+j3i4uJ6/I7Yjm3hie3xLbaFJ7bHt9gWntge3+rTp0+3FxuJCcUAACDIUG4AAEBQodwch9Vq1UMPPSSr1Wp2FNOxLTyxPb7FtvDE9vgW28IT2+Nb3t4WAT+hGAAA4LsYuQEAAEGFcgMAAIIK5QYAAAQVyg0AAAgqlBsAABBUKDfHsGjRIvXv31+RkZHKzc3V+vXrzY7kdfn5+Ro/frxiY2OVkpKiadOmqbi42GOdCRMmyGKxeDxuv/12kxJ718MPP3zUzzpkyBD3683NzZo1a5YSExMVExOj6dOnq7Ky0sTE3tO/f/+jtoXFYtGsWbMkBf9+sXr1al122WXKyMiQxWLR22+/7fG6YRiaP3++0tPT1atXL02cOFE7duzwWKempkbXX3+94uLiFB8fr5tuukkNDQ0+/Cm6x/G2RWtrq+6//36NGDFC0dHRysjI0A033KD9+/d7fI2O9qfHH3/cxz9J9zjRvjFz5syjftYpU6Z4rNMT9g1JHf4NsVgseuqpp9zrdNe+QbnpwN/+9jfNnTtXDz30kDZt2qRRo0Zp8uTJqqqqMjuaV61atUqzZs3S2rVrtWzZMrW2tmrSpElqbGz0WO+WW25ReXm5+/Hkk0+alNj7hg0b5vGzfvrpp+7X5syZo3fffVdvvPGGVq1apf379+vKK680Ma33bNiwwWM7LFu2TJJ01VVXudcJ5v2isbFRo0aN0qJFizp8/cknn9Szzz6r559/XuvWrVN0dLQmT56s5uZm9zrXX3+9tmzZomXLlum9997T6tWrdeutt/rqR+g2x9sWTU1N2rRpkx588EFt2rRJb775poqLi/WjH/3oqHUXLlzosb/84he/8EX8bneifUOSpkyZ4vGzvvbaax6v94R9Q5LHNigvL9fixYtlsVg0ffp0j/W6Zd8wcJQzzzzTmDVrlvu50+k0MjIyjPz8fBNT+V5VVZUhyVi1apV72Q9+8APj7rvvNi+UDz300EPGqFGjOnyttrbWCA8PN9544w33sm3bthmSjIKCAh8lNM/dd99tDBgwwHC5XIZh9Kz9QpLx1ltvuZ+7XC4jLS3NeOqpp9zLamtrDavVarz22muGYRjG1q1bDUnGhg0b3Ov8+9//NiwWi7Fv3z6fZe9u398WHVm/fr0hydizZ497Wb9+/YxnnnnGu+FM0NH2mDFjhnH55Zcf8z09ed+4/PLLjQsvvNBjWXftG4zcfE9LS4sKCws1ceJE97KQkBBNnDhRBQUFJibzvbq6OklSQkKCx/K//vWvSkpK0vDhwzVv3jw1NTWZEc8nduzYoYyMDOXk5Oj6669XaWmpJKmwsFCtra0e+8mQIUOUlZUV9PtJS0uLXnnlFf385z+XxWJxL+9J+8V3lZSUqKKiwmNfsNlsys3Nde8LBQUFio+P1xlnnOFeZ+LEiQoJCdG6det8ntmX6urqZLFYFB8f77H88ccfV2JiosaMGaOnnnpKbW1t5gT0gZUrVyolJUWDBw/WHXfcoerqavdrPXXfqKys1D//+U/ddNNNR73WHftGwN8VvLsdPHhQTqdTqampHstTU1O1fft2k1L5nsvl0j333KNzzjlHw4cPdy+/7rrr1K9fP2VkZOjzzz/X/fffr+LiYr355psmpvWO3NxcLVmyRIMHD1Z5ebkWLFig8847T19++aUqKioUERFx1B/s1NRUVVRUmBPYR95++23V1tZq5syZ7mU9ab/4vvb/3h39zWh/raKiQikpKR6vh4WFKSEhIaj3l+bmZt1///269tprPe6Cfdddd2ns2LFKSEjQmjVrNG/ePJWXl+vpp582Ma13TJkyRVdeeaWys7O1a9cu/frXv9bUqVNVUFCg0NDQHrtvvPTSS4qNjT3qUH537RuUG3Ro1qxZ+vLLLz3mmEjyOA48YsQIpaen66KLLtKuXbs0YMAAX8f0qqlTp7r/PXLkSOXm5qpfv356/fXX1atXLxOTmeuFF17Q1KlTlZGR4V7Wk/YLdE5ra6t+8pOfyDAMPffccx6vzZ071/3vkSNHKiIiQrfddpvy8/OD7r5L11xzjfvfI0aM0MiRIzVgwACtXLlSF110kYnJzLV48WJdf/31ioyM9FjeXfsGh6W+JykpSaGhoUed9VJZWam0tDSTUvnW7Nmz9d5772nFihXq27fvcdfNzc2VJO3cudMX0UwVHx+v0047TTt37lRaWppaWlpUW1vrsU6w7yd79uzRRx99pJtvvvm46/Wk/aL9v/fx/makpaUddUJCW1ubampqgnJ/aS82e/bs0bJlyzxGbTqSm5urtrY27d692zcBTZSTk6OkpCT370ZP2zck6ZNPPlFxcfEJ/45IJ79vUG6+JyIiQuPGjdPy5cvdy1wul5YvX668vDwTk3mfYRiaPXu23nrrLX388cfKzs4+4XuKiookSenp6V5OZ76Ghgbt2rVL6enpGjdunMLDwz32k+LiYpWWlgb1fvLiiy8qJSVFl1566XHX60n7RXZ2ttLS0jz2BbvdrnXr1rn3hby8PNXW1qqwsNC9zscffyyXy+UugsGivdjs2LFDH330kRITE0/4nqKiIoWEhBx1eCYY7d27V9XV1e7fjZ60b7R74YUXNG7cOI0aNeqE6570vnHKU5KD0NKlSw2r1WosWbLE2Lp1q3Hrrbca8fHxRkVFhdnRvOqOO+4wbDabsXLlSqO8vNz9aGpqMgzDMHbu3GksXLjQ2Lhxo1FSUmK88847Rk5OjnH++eebnNw77r33XmPlypVGSUmJ8Z///MeYOHGikZSUZFRVVRmGYRi33367kZWVZXz88cfGxo0bjby8PCMvL8/k1N7jdDqNrKws4/777/dY3hP2i/r6emPz5s3G5s2bDUnG008/bWzevNl9BtDjjz9uxMfHG++8847x+eefG5dffrmRnZ1tHD582P01pkyZYowZM8ZYt26d8emnnxqDBg0yrr32WrN+pJN2vG3R0tJi/OhHPzL69u1rFBUVefwdcTgchmEYxpo1a4xnnnnGKCoqMnbt2mW88sorRnJysnHDDTeY/JOdnONtj/r6euO+++4zCgoKjJKSEuOjjz4yxo4dawwaNMhobm52f42esG+0q6urM6KiooznnnvuqPd3575BuTmGP/zhD0ZWVpYRERFhnHnmmcbatWvNjuR1kjp8vPjii4ZhGEZpaalx/vnnGwkJCYbVajUGDhxo/PKXvzTq6urMDe4lV199tZGenm5EREQYffr0Ma6++mpj586d7tcPHz5s3HnnnUbv3r2NqKgo44orrjDKy8tNTOxdH3zwgSHJKC4u9ljeE/aLFStWdPi7MWPGDMMwjpwO/uCDDxqpqamG1Wo1LrrooqO2U3V1tXHttdcaMTExRlxcnHHjjTca9fX1Jvw0p+Z426KkpOSYf0dWrFhhGIZhFBYWGrm5uYbNZjMiIyON008/3Xjsscc8PuwDyfG2R1NTkzFp0iQjOTnZCA8PN/r162fccsstR/2Pck/YN9r96U9/Mnr16mXU1tYe9f7u3DcshmEYXRvrAQAA8F/MuQEAAEGFcgMAAIIK5QYAAAQVyg0AAAgqlBsAABBUKDcAACCoUG4AAEBQodwAAICgQrkBAABBhXIDAACCCuUGAAAElf8P7Y3kKRa5/nAAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "plt.plot(np.sort(dh.eig()))\n", + "plt.ylim(None, 20)\n", + "np.real(dh.eig()).min()" + ] + }, + { + "cell_type": "code", + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -249,7 +327,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -259,6 +337,24 @@ " magnetic_entities[i][\"orbital_indeces\"] = parsed\n", " # calculate spin box indexes\n", " magnetic_entities[i][\"spin_box_indeces\"] = blow_up_orbindx(parsed)\n", + " # if orbital is not set use all\n", + " if \"l\" not in mag_ent.keys():\n", + " mag_ent[\"l\"] = \"all\"\n", + " if isinstance(mag_ent[\"atom\"], int):\n", + " mag_ent[\"tags\"] = [\n", + " f\"[{mag_ent['atom']}]{dh.atoms[mag_ent['atom']].tag}({mag_ent['l']})\"\n", + " ]\n", + " mag_ent[\"xyz\"] = [dh.xyz[mag_ent[\"atom\"]]]\n", + " if isinstance(mag_ent[\"atom\"], list):\n", + " mag_ent[\"tags\"] = []\n", + " mag_ent[\"xyz\"] = []\n", + " # iterate over atoms\n", + " for atom_idx in mag_ent[\"atom\"]:\n", + " mag_ent[\"tags\"].append(\n", + " f\"[{atom_idx}]{dh.atoms[atom_idx].tag}({mag_ent['l']})\"\n", + " )\n", + " mag_ent[\"xyz\"].append(dh.xyz[atom_idx])\n", + "\n", " # calculate size for Greens function generation\n", " spin_box_shape = len(mag_ent[\"spin_box_indeces\"])\n", "\n", @@ -284,7 +380,26 @@ " # calculate size for Greens function generation\n", " spin_box_shape_i = len(magnetic_entities[pair[\"ai\"]][\"spin_box_indeces\"])\n", " spin_box_shape_j = len(magnetic_entities[pair[\"aj\"]][\"spin_box_indeces\"])\n", - "\n", + " pair[\"tags\"] = []\n", + " for mag_ent in [magnetic_entities[pair[\"ai\"]], magnetic_entities[pair[\"aj\"]]]:\n", + " tag = \"\"\n", + " # get atoms of magnetic entity\n", + " atoms_idx = mag_ent[\"atom\"]\n", + " orbitals = mag_ent[\"l\"]\n", + "\n", + " # if magnetic entity contains one atoms\n", + " if isinstance(atoms_idx, int):\n", + " tag += f\"[{atoms_idx}]{dh.atoms[atoms_idx].tag}({orbitals})\"\n", + "\n", + " # if magnetic entity contains more than one atoms\n", + " if isinstance(atoms_idx, list):\n", + " # iterate over atoms\n", + " atom_group = \"{\"\n", + " for atom_idx in atoms_idx:\n", + " atom_group += f\"[{atom_idx}]{dh.atoms[atom_idx].tag}({orbitals})--\"\n", + " # end {} of the atoms in the magnetic entity\n", + " tag += atom_group[:-2] + \"}\"\n", + " pair[\"tags\"].append(tag)\n", " pair[\"energies\"] = [] # we will store the second order energy derivations here\n", "\n", " pair[\"Gij\"] = [] # Greens function\n", @@ -311,21 +426,29 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "k loop: 0%| | 0/100 [00:00