You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
grogu/test.ipynb

1803 lines
156 KiB

{
"cells": [
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'0.14.3'"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import os\n",
"from tqdm import tqdm\n",
"from sys import getsizeof\n",
"from timeit import default_timer as timer\n",
"\n",
"os.environ[\"OMP_NUM_THREADS\"] = \"1\" # export OMP_NUM_THREADS=4\n",
"os.environ[\"OPENBLAS_NUM_THREADS\"] = \"1\" # export OPENBLAS_NUM_THREADS=4\n",
"os.environ[\"MKL_NUM_THREADS\"] = \"1\" # export MKL_NUM_THREADS=6\n",
"os.environ[\"VECLIB_MAXIMUM_THREADS\"] = \"1\" # export VECLIB_MAXIMUM_THREADS=4\n",
"os.environ[\"NUMEXPR_NUM_THREADS\"] = \"1\" # export NUMEXPR_NUM_THREADS=6\n",
"\n",
"import numpy as np\n",
"import sisl\n",
"from src.grogu_magn.useful import *\n",
"from mpi4py import MPI\n",
"import pickle\n",
"from numpy.linalg import inv\n",
"import warnings\n",
"\n",
"# runtime information\n",
"times = dict()\n",
"times[\"start_time\"] = timer()\n",
"########################\n",
"# it works if data is in downloads folder\n",
"########################\n",
"sisl.__version__"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"dat = sisl.io.siesta.eigSileSiesta(\n",
" \"/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.EIG\"\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-12.806878959999999"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dat.read_data().min()"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"================================================================================================================================================================\n",
"Input file: \n",
"Not yet specified.\n",
"Output file: \n",
"test.pickle\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: 20\n",
"k point directions: xy\n",
"Ebot: -13\n",
"Eset: 600\n",
"Esetp: 10000\n",
"================================================================================================================================================================\n",
"Setup done. Elapsed time: 198.668371833 s\n",
"================================================================================================================================================================\n"
]
}
],
"source": [
"# this cell mimicks an input file\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",
"\n",
"outfile = \"test\"\n",
"if not outfile.endswith(\".pickle\"):\n",
" outfile += \".pickle\"\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",
"# o is the quantization axis, v and w are two axes perpendicular to it\n",
"# at this moment the user has to supply o,v,w on the input.\n",
"# 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, 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",
"\n",
"# 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=5, l=2),\n",
"]\n",
"# pair information ./lat3_791/Fe3GeTe2.fdf\n",
"pairs = [\n",
" # isotropic should be -82 meV\n",
" dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])),\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",
"]\n",
"\n",
"\"\"\"\n",
"# human readable definition of magnetic entities ./Jij_for_Marci_6p45ang/CrBr.fdf\n",
"magnetic_entities = [\n",
" dict(atom=0, l=2),\n",
" dict(atom=1, l=2),\n",
" dict(atom=2, l=2),\n",
"]\n",
"# pair information ./Jij_for_Marci_6p45ang/CrBr.fdf\n",
"pairs = [\n",
" dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])),\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=1, Ruc=np.array([1, 0, 0])),\n",
" dict(ai=0, aj=2, Ruc=np.array([1, 0, 0])),\n",
" dict(ai=0, aj=1, Ruc=np.array([-1, 0, 0])),\n",
" dict(ai=0, aj=2, Ruc=np.array([-1, 0, 0])),\n",
" dict(ai=0, aj=1, Ruc=np.array([0, 1, 0])),\n",
" dict(ai=0, aj=2, Ruc=np.array([0, 1, 0])),\n",
" dict(ai=0, aj=1, Ruc=np.array([0, 1, 0])),\n",
" dict(ai=0, aj=2, Ruc=np.array([0, 1, 0])),\n",
"]\n",
"\"\"\"\n",
"\n",
"# Brilloun zone sampling and Green function contour integral\n",
"kset = 20\n",
"kdirs = \"xy\"\n",
"ebot = -13\n",
"eset = 600\n",
"esetp = 10000\n",
"\n",
"\n",
"# MPI parameters\n",
"comm = MPI.COMM_WORLD\n",
"size = comm.Get_size()\n",
"rank = comm.Get_rank()\n",
"root_node = 0\n",
"\n",
"simulation_parameters = dict(\n",
" path=\"Not yet specified.\",\n",
" outpath=outfile,\n",
" scf_xcf_orientation=scf_xcf_orientation,\n",
" ref_xcf_orientations=ref_xcf_orientations,\n",
" kset=kset,\n",
" kdirs=kdirs,\n",
" ebot=ebot,\n",
" eset=eset,\n",
" esetp=esetp,\n",
" parallel_size=size,\n",
")\n",
"\n",
"# digestion of the input\n",
"# read in hamiltonian\n",
"dh = fdf.read_hamiltonian()\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",
"if rank == root_node:\n",
" print_parameters(simulation_parameters)\n",
" times[\"setup_time\"] = timer()\n",
" print(f\"Setup done. Elapsed time: {times['setup_time']} s\")\n",
" print(\n",
" \"================================================================================================================================================================\"\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"xyz[-3:]: red, green, blue\n"
]
},
{
"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": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABLwAAAGsCAYAAADXMb4GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB3yUlEQVR4nO3deZxbdb3/8Xf2ZPalna2d7i1l3wqlgCxSCogIF0RRvFZB0GtBoPeK1Ct4QaSCilyQS8GrqFdB8aeAG0upLAKllEJZS/eNtjPdZiazJTlJzu+PLJPMks5MM3Mymdfz8ehjkpOTk2++hJkz7/l8P8dmmqYpAAAAAAAAIE/YrR4AAAAAAAAAkE0EXgAAAAAAAMgrBF4AAAAAAADIKwReAAAAAAAAyCsEXgAAAAAAAMgrBF4AAAAAAADIKwReAAAAAAAAyCtOqwfQXTQa1c6dO1VcXCybzWb1cAAAwAhhmqZaW1tVV1cnu52/6eUizvMAAMBgDOY8L+cCr507d6q+vt7qYQAAgBFq+/btGj9+vNXDQC84zwMAAAdjIOd5ORd4FRcXS4q9iZKSkqwf3zAMPfvss5o3b55cLlfWjz/SMT+ZMT+ZMT+ZMT+ZMT+ZMT+ZGYahJ554Ql/5yleS5xLIPUN9ngcAAPKT3+9XfX39gM7zci7wSpS3l5SUDFngVVBQoJKSEn5h6AXzkxnzkxnzkxnzkxnzkxnzk1lifiSxVC6HDfV5HgAAyG8DOc+jwQUAAAAAAADyCoEXAAAAAAAA8gqBFwAAAAAAAPIKgRcAAAAAAADyCoEXAAAAAAAA8gqBFwAAAAAAAPIKgRcAAAAAAADyCoEXAAAAAAAA8gqBFwAAAAAAAPIKgRcAAAAAAADyCoEXAAAAAAAA8gqBFwAAAAAAAPIKgRcAAAAAAADyitPqAQAAgPwXMCJ656MWmaap2VMqrR4O8tSG3W1qaAloQkWBJlQWWD0cAABgIQIvAACQdfvbQ3pjy36t2tqklVv2690dLTIipuZMqdSjVxN4YWj88tXN+s1r23TdWdN1w9kzrB4OAACwEIEXAAA4aE3tIT23plFvbGnSG1v3a+Oe9h77jC32qK7MZ8HoAAAAMNoQeAEAgIP2uZ+9pg8bWtO2Ta8q0qxJFZo1sVwnTKpQfYVPNpvNohECAABgNCHwAgAAB23rvg5J0hfnTNTpM8bquAnlKi90WzwqAAAAjFYEXgAA4KCYpqlgOCJJuubMaaoq8Vo8IgAAAIx29oHsHIlEdPPNN2vy5Mny+XyaOnWqvve978k0zeQ+pmnqlltuUW1trXw+n+bOnav169dnfeAAACA3hKOmovFTAY/TYe1gAAAAAA0w8Lrzzjv1wAMP6Kc//anWrFmjO++8U3fddZfuu+++5D533XWX7r33Xi1ZskQrVqxQYWGhzjnnHAUCgawPHgAAWC8YjiZve1wDOrUAAAAAhsSAljS++uqruvDCC3X++edLkiZNmqRHH31Ur7/+uqRYddc999yj73znO7rwwgslSb/+9a9VXV2tJ554QpdddlmPYwaDQQWDweR9v98vSTIMQ4ZhDO5dZZA45lAcOx8wP5kxP5kxP5kxP5kxP5nl8vy0d4aSt23RiAwjmmHvoZGL8wIAAADrDCjwOvnkk/XQQw9p3bp1mjFjht5++229/PLLuvvuuyVJmzdvVkNDg+bOnZt8TmlpqWbPnq3ly5f3GngtXrxYt956a4/tzz77rAoKCgb6fvpt6dKlQ3bsfMD8ZMb8ZMb8ZMb8ZMb8ZJaL89MclCSnHDZTTz/9lNXDAQAAAAYWeN10003y+/2aOXOmHA6HIpGIvv/97+vyyy+XJDU0NEiSqqur055XXV2dfKy7RYsWaeHChcn7fr9f9fX1mjdvnkpKSgb0ZvrDMAwtXbpUZ599tlwuV9aPP9IxP5kxP5kxP5kxP5kxP5nl8vxs3dchvfmyvG6nPvGJcywZg2EYevLJJy15bQAAAOSeAQVejz32mH7729/qkUce0eGHH67Vq1fr+uuvV11dnebPnz+oAXg8Hnk8nh7bXS7XkJ7QD/XxRzrmJzPmJzPmJzPmJzPmJ7NcnJ9IvCWo1+nIubEBAABgdBpQ4PXNb35TN910U3Jp4pFHHqmtW7dq8eLFmj9/vmpqaiRJjY2Nqq2tTT6vsbFRxxxzTPZGDQAAckYo3rTe7aRhPQAAAHLDgM5MOzo6ZLenP8XhcCgajZ3oTp48WTU1NVq2bFnycb/frxUrVmjOnDlZGC4AAMg1wXBEkuQh8AIAAECOGFCF1wUXXKDvf//7mjBhgg4//HC99dZbuvvuu3XFFVdIkmw2m66//nrdfvvtmj59uiZPnqybb75ZdXV1uuiii4Zi/AAAwGLBeIWXx+mweCQAAABAzIACr/vuu08333yzvv71r2v37t2qq6vTV7/6Vd1yyy3JfW688Ua1t7fr6quvVnNzs0499VQ9/fTT8nq9WR88AACwXmJJo8dFhRcAAAByw4DOTIuLi3XPPfdo69at6uzs1MaNG3X77bfL7XYn97HZbLrtttvU0NCgQCCg5557TjNmzMj6wAEAQG5ILGl0Owi8RrJIJKKbb75ZkydPls/n09SpU/W9731Ppmkm9zFNU7fccotqa2vl8/k0d+5crV+/3sJRAwAA9I4zUwAAcFCCVHjlhTvvvFMPPPCAfvrTn2rNmjW68847ddddd+m+++5L7nPXXXfp3nvv1ZIlS7RixQoVFhbqnHPOUSAQsHDkAAAAPQ1oSSMAAEB3QYMeXvng1Vdf1YUXXqjzzz9fkjRp0iQ9+uijev311yXFqrvuuecefec739GFF14oSfr1r3+t6upqPfHEE8mreKcKBoMKBoPJ+36/fxjeCQAAABVeAADgIAUjscCLJY0j28knn6xly5Zp3bp1kqS3335bL7/8ss477zxJ0ubNm9XQ0KC5c+cmn1NaWqrZs2dr+fLlvR5z8eLFKi0tTf6rr68f+jcCAAAgKrwAAMBBChqxHl4saRzZbrrpJvn9fs2cOVMOh0ORSETf//73dfnll0uSGhoaJEnV1dVpz6uurk4+1t2iRYu0cOHC5H2/30/oBQAAhgWBFwAAOCjJHl5OAq+R7LHHHtNvf/tbPfLIIzr88MO1evVqXX/99aqrq9P8+fMHdUyPxyOPx5PlkQIAABwYgRcAADgoXYEXPbxGsm9+85u66aabkr24jjzySG3dulWLFy/W/PnzVVNTI0lqbGxUbW1t8nmNjY065phjrBgyAABAn/hTLAAAOCiheODlpsJrROvo6JDdnv7f0OFwKBqN/fedPHmyampqtGzZsuTjfr9fK1as0Jw5c4Z1rAAAAAdChRcAADgowXC8hxeB14h2wQUX6Pvf/74mTJigww8/XG+99ZbuvvtuXXHFFZIkm82m66+/XrfffrumT5+uyZMn6+abb1ZdXZ0uuugiawcPAADQDYEXAAA4KCxpzA/33Xefbr75Zn3961/X7t27VVdXp69+9au65ZZbkvvceOONam9v19VXX63m5madeuqpevrpp+X1ei0cOQAAQE8EXgAA4KAEjXjgxVUaR7Ti4mLdc889uueee/rcx2az6bbbbtNtt902fAMDAAAYBM5MAQDAQQlF4j28HJxWAAAAIDdwZgoAAA5K0Ij38KLCCwAAADmCM1MAAHBQ6OEFAACAXEPgBQAADgpXaQQAAECu4cwUAAAclFC8wstN4AUAAIAcwZkpAAA4KF1LGjmtAAAAQG7gzBQAABwUengBAAAg1xB4AQCAg5JY0shVGgEAAJArODMFAAAHJdG03u3gtAIAAAC5gTNTAABwUBJLGr1UeAEAACBHcGYKAAAOStCghxcAAAByC4EXAAA4KKEIV2kEAABAbuHMFAAADFo4ElUkakqS3AReAAAAyBGcmQIAgEFL9O+SWNIIAACA3EHgBQAABi018KLCCwAAALmCM1MAADBooXjg5bTb5LDbLB4NAAAAEEPgBQAABi0YjkiiYT0AAAByC2enAABg0BJLGj0u+ncBAAAgdxB4AQCAQQsa8cCLCi8AAADkEM5OAQDAoIUisSWNNKwHAABALuHsFAAADBoVXgAAAMhFnJ0CAIBBS/bwctLDCwAAALmDwAsAAAwaV2kEAABALuLsFAAADFqiwoseXgAAAMglnJ0CAIBB61rSyCkFAAAAcgdnpwAAYNDo4QUAAIBcROAFAAAGLZQIvFycUgAAACB3cHYKAAAGLdG03u3glAIAAAC5g7NTAAAwaEGDCi8AAADkHs5OAQDAoNHDCwAAALmIwAsAAAxaiKs0AgAAIAdxdgoAAAYt2cOLwAsAAAA5hLNTAAAwaCxpBAAAQC4i8AIAAIMWZEkjAAAAchBnpwAAYNBCLGkEAABADuLsFAAADBoVXgAAAMhFnJ0CAIBBCxrxwMtFDy8AAADkDgIvAAAwaImrNFLhBQAAgFzC2SkAABi0UCRW4UUPLwAAAOQSzk4BAMCgJZc0EngBAAAgh3B2CgAABq2raT09vAAAAJA7CLwAAMCg0cMLAAAAuYizUwAAMGihMEsaAQAAkHs4OwUAAIPGkkYAAADkIgIvAAAwaMnAy8UpBQAAAHIHZ6cAAGBQwpGoIlFTEksaAQAAkFs4OwUAAIMSikSTt90EXgAAAMghnJ0CAIBBCRopgZeDUwoAAADkDs5OAQDAoCT6dzntNjkJvAAAAJBDODsFAACDEkpeoZHTCQAAAOQWzlABAMCgBMMRSfTvAgAAQO7hDBUAAAxKMFnh5bB4JAAAAEC6AQdeO3bs0Be+8AVVVlbK5/PpyCOP1BtvvJF83DRN3XLLLaqtrZXP59PcuXO1fv36rA4aAABYL1Hh5XHx9zMAAADklgGdoTY1NemUU06Ry+XSU089pQ8++EA//vGPVV5entznrrvu0r333qslS5ZoxYoVKiws1DnnnKNAIJD1wQMAAOskKry4QiMAAAByjXMgO995552qr6/Xww8/nNw2efLk5G3TNHXPPffoO9/5ji688EJJ0q9//WtVV1friSee0GWXXZalYQMAAKsllzRS4QUAAIAcM6DA689//rPOOeccXXrppXrxxRc1btw4ff3rX9dVV10lSdq8ebMaGho0d+7c5HNKS0s1e/ZsLV++vNfAKxgMKhgMJu/7/X5JkmEYMgxjUG8qk8Qxh+LY+YD5yYz5yYz5yYz5yYz5ySwX56cjEJIUq/CyelxWvz4AAAByy4ACr02bNumBBx7QwoUL9e1vf1srV67UN77xDbndbs2fP18NDQ2SpOrq6rTnVVdXJx/rbvHixbr11lt7bH/22WdVUFAwkOENyNKlS4fs2PmA+cmM+cmM+cmM+cmM+cksl+Zn1V6bJIdam/fr73//u9XDAQAAAJIGFHhFo1HNmjVLd9xxhyTp2GOP1XvvvaclS5Zo/vz5gxrAokWLtHDhwuR9v9+v+vp6zZs3TyUlJYM6ZiaGYWjp0qU6++yz5XK5sn78kY75yYz5yYz5yYz5yYz5ySwX56fzzR3S+vdVV12lT3ziOEvHYhiGnnzySUvHAAAAgNwxoMCrtrZWhx12WNq2Qw89VH/84x8lSTU1NZKkxsZG1dbWJvdpbGzUMccc0+sxPR6PPB5Pj+0ul2tIT+iH+vgjHfOTGfOTGfOTGfOTGfOTWS7NT9i0SZK8bkfOjAkAAACQBniVxlNOOUVr165N27Zu3TpNnDhRUqyBfU1NjZYtW5Z83O/3a8WKFZozZ04WhgsAAHJFsmm902HxSJAtO3bs0Be+8AVVVlbK5/PpyCOP1BtvvJF83DRN3XLLLaqtrZXP59PcuXO1fv16C0cMAADQuwEFXjfccINee+013XHHHdqwYYMeeeQRPfTQQ1qwYIEkyWaz6frrr9ftt9+uP//5z3r33Xf1xS9+UXV1dbrooouGYvwAAMAiwXBEkuRxcpXGfNDU1KRTTjlFLpdLTz31lD744AP9+Mc/Vnl5eXKfu+66S/fee6+WLFmiFStWqLCwUOecc44CgYCFIwcAAOhpQEsaTzjhBD3++ONatGiRbrvtNk2ePFn33HOPLr/88uQ+N954o9rb23X11VerublZp556qp5++ml5vd6sDx4AAFgnFK/wchN45YU777xT9fX1evjhh5PbJk+enLxtmqbuuecefec739GFF14oSfr1r3+t6upqPfHEEwO6GjcAAMBQG/AZ6ic/+Um9++67CgQCWrNmja666qq0x202m2677TY1NDQoEAjoueee04wZM7I2YAAAkBtY0phf/vznP2vWrFm69NJLVVVVpWOPPVY/+9nPko9v3rxZDQ0Nmjt3bnJbaWmpZs+ereXLl/d6zMWLF6u0tDT5r76+fsjfBwAAgDSIwAsAAECSgkY88HJxOpEPNm3apAceeEDTp0/XM888o3/7t3/TN77xDf3qV7+SJDU0NEiSqqur055XXV2dfKy7RYsWqaWlJflv+/btQ/smAAAA4ga0pBEAACAhFKGHVz6JRqOaNWuW7rjjDknSscceq/fee09LlizR/PnzB3XMvq7GDQAAMNQ4QwUAAIOSqPCih1d+qK2t1WGHHZa27dBDD9W2bdskSTU1NZKkxsbGtH0aGxuTjwEAAOQKzlABAMCg0MMrv5xyyilau3Zt2rZ169Zp4sSJkmIN7GtqarRs2bLk436/XytWrNCcOXOGdawAAAAHwpJGAAAwKMEwSxrzyQ033KCTTz5Zd9xxhz7zmc/o9ddf10MPPaSHHnpIUuzCRNdff71uv/12TZ8+XZMnT9bNN9+suro6XXTRRdYOHgAAoBsCLwAAMCihZIUXgVc+OOGEE/T4449r0aJFuu222zR58mTdc889uvzyy5P73HjjjWpvb9fVV1+t5uZmnXrqqXr66afl9XotHDkAAEBPBF4AAGBQEksa6eGVPz75yU/qk5/8ZJ+P22w23XbbbbrtttuGcVQAAAADxxkqAAAYFHp4AQAAIFdR4QUAAPrlo6YO/fKVLfrHh7vV6A+oPUQPLwAAAOQmAi8AANCDaZpatbVJ25s6tKc1qLc/atHT7zUoEjXT9htf7tNR40stGiUAAADQOwIvAADQw6Ovb9e3H3+3x/ZTplXqi3MmaUZ1scYWe1Todshms1kwQgAAAKBvBF4AACBNNGrqoZc2SpKOGl+qqWOLVFPq1SePqtXhdVRzAQAAIPcReAEAgDT/+HC3tuzrUInXqUevOkmFHk4XAAAAMLLQZRYAAKT5+cubJUmfO3ECYRcAAABGJAIvAACQ9MFOv5Zv2ieH3ab5J0+yejgAAADAoBB4AQCApF+8EqvuOu+IGtWV+SweDQAAADA4rFMAAGAUW9vQqt+t3KadzZ0KhaN6ZcM+SdIVp062eGQAAADA4BF4AQAwSjS0BLRqa5OC4Yg6QhE9/V6DXt6wt8d+J0wq13ETyi0YIQAAAJAdBF4AAIwSn3lwubbt70jbZrdJ8w6r0SnTKuVxOuRx2XXqtDEWjRAAAADIDgIvAABGgf3toWTY9bHpY+Rx2jWtqliXz56g+ooCi0cHAAAAZBeBFwAAo8DGPW2SpHFlPv3flbMtHg0AAAAwtLhKIwAAo8CmeOA1ZWyhxSMBAAAAhh6BFwAAo8CmPe2SpCljCLwAAACQ/wi8AAAYBTYmAq+xRRaPBAAAABh6BF4AAIwCm/aypBEAAACjB4EXAAB5zohEtW1f7AqNVHgBAABgNCDwAgAgz23f36Fw1JTP5VBtidfq4QAAAABDjsALAIA8l2hYP3lMoex2m8WjAQAAAIYegRcAAHmO/l0AAAAYbQi8AADIc5u4QiMAAABGGQIvAADyXCLwmkqFFwAAAEYJAi8AAPLcxj3xJY1jqPACAADA6EDgBQBAHmvpMLSvPSRJmkyFF0YJ0+oBAAAAyxF4AQCQxzbGG9ZXl3hU5HFaPBpgaI0p8kiSXlq3R6ZJ7AUAwGhG4AUAQB7r6t/Fckbkv8/PniCvy67V25v1/NrdVg8HAABYiMALAIA8tinRv4vljBgFqoq9mj9nkiTp7qXrqPICAGAUI/ACACCPmKapt7Y16X//uUkPvLBRz6/dI4mG9Rg9vnr6VBW6HXpvh1/PvN9o9XAAAIBFaOYBAMAI19Qe0vs7/VqxeZ/+/PZObd3X0WOfGdXFFowMGH4VhW5dcepk3fePDfrJ0nWad1i17Hab1cMCAADDjMALAIAcEY5EtWFPm3Y1BxQwIgqEIwoYUQWMiNoDht7dbte7z6xTKGJqX1tIDf6AdjZ3aldLIO04PpdDH5s+RqU+l5wOm8aXF+jkqZUWvStg+H3l1Cn65atbtLaxVX99d5c+dXSd1UMCAADDjMALAAALRaOmfrNiq/745g59uMuvYDiaYW+79NGWXh+ZWFmgI+pKdfZh1Tr7sGoVckVGjGKlBS5d9bEpunvpOi3++xq9s71ZtWU+1ZV6VVvmU22pV2OKPHJQ+QUAQN7ibBgAAIvsawvqP/7wdrLPliQVe5yaOKZAPpdDXpdDHqdDXpddbodNjTs/0iFTJ6vA41JlkVs1JV5VlXg1vbpIJV6Xhe8EyD1fPmWSfvnqFu1qCeh/X97c43Gn3abqEq9q4yFYXalXNaVe1Zb6VFcW+1pZ6GY5JAAAIxSBFwAAFnh1w17d8NhqNfqDcjvt+ua8Q3T2YdWaUFHQ6y/YhmHo73/fpk+cd4hcLsIt4ECKvS499tWT9OK6vdoVX/q7qyX2tdEfUDhqakdzp3Y0d0pbm3o9htthV3WpR7WlsaqwRBhWU+JVXbxSrKLQLZuNUAwAgFxD4AUAwDDqCIV151Mf6lfLt0qSpo4t1E8/f5wOrS2xeGRA/plWVaxpVT0v2BCORLWnLaidzbEQrKElkLy9syWghpZO7W4NKhSJavv+Tm3f39nna7id9ngY5u0Kxsp8qi3xqrbMq7pSn8oKXIRiAAAMMwIvAACGUCRq6uUNe7VpT5t2tQT0zPsNyasofn72BH3n/ENV4ObHMTCcnA57PJzySSrvdR8jElWjPxALw1oCParEdjYHtLctqFA4qq37Onq9OmqC12VPqxKLhWKxMKymNPa1xOckFAMAIIs4wwYAYAj9buU2/efj76Vtqy316s5LjtJpM8ZaNCoAB+Jy2DW+vEDjywv63CcUjoViO5s741dN7QrEdrV0aldzQPvaQwoYUW3e267Ne9v7PFaB29GzSixeOVZXFgvG6NUHAED/EXgBADCEnv8w1pD+mPoyHT+xXBMrC3ThMeNU6uMXV2Ckczvtqq8oUH1F36FYwIjEQ7GeYVjidlOHoY5QRBv3tGvjnr5DsSKPU7Xx5vp1pT7VlnnTeovVlvq4QisAAHH8RAQAYIhEoqZe37xPknTrpw7X0fVl1g4IwLDzuhyaWFmoiZWFfe7TGYp09RJLLJ/0py6jDKil01BbMKz1u9u0fndbn8cq9jrTwrALjxmnk6ZUDsVbAwAgpxF4AQAwRD5s8MsfCKvI49ThdTSlB9A7n9uhKWOLNGVsUa+PB4yINu1p1+rtzXpzW5Pe2tbUZyVYayCstYFWrW1slSSt3t6ip6772JCNHQCAXEXgBQDAEHlt035J0qxJ5XI67BaPBkCuMU1T+9tDavAH4g3yg2poifUDa/AH1dgSUIM/Vt3VH3abNLbYo5oSr6pLYksfLzi6bojfBQAAuYnACwCAIfLapthyRpYTAaNPMBzRbn8wFl61BJJXfNzlDySDrN3+oEKRaL+O53PFmtongqzqEq9qSjzJ27WlPo0pchOuAwAQR+AFAMAQiEZNvb45VuE1e3KFxaMBkC2maaql0+gWZAXV4I/14GrwB9XoD2h/e6jfxxxT5I6HVokgy6vq0tjXRKBV4nXKZrMN4TsDACC/EHgBADAEPmxoVUunoUK3Q0eMK7V6OAD6wYhEtbs1mFaR1eiPNY1PLDts9AcUMPpXleV22mOhVTzASg20ako9qi7xqqrYK7eTqiwAALKNwAsAgCGQWM54/KQKuVhiBFjKNE21BsPJpYTJQCt+uyFepbWvPSjT7N8xywtcyeWFqT2zEvdrSrwqK3BRlQUAgEUIvAAAGAJd/btYzggMl/ZgWH9+e6e27e9IC7ca/AF1hCL9OobLYVNVcfcgy5OszKot9amqxCOvyzHE7wYAABwMAi8AALIsGjX1+pZY/y4a1gPD5xcvb9aPl64b1HNLvE5NGVuk+ooCFXudKnQ7VOB2qtAT+2qaUkunoXDU1L72oIo8ThV4uvZjWSIAALmFwAsAgCxb/VGzmjsMFbgdOpL+XcCw+cRRtVq/u01NHSG1B8PqCEXUFv/aHgwrGO6795Y/ENbq7c1avb15UK/tcthU6HGq0O1UgduhAo9TRfGwrDB+v9Dt6NrH41Ch2xm/3/V4gcepovjjLIcGAGDwCLwAAMiiVVubdNWv35AknTJtDL+wAsNo6tgi3fu5Y/t8PByJqj0UUUcorPZg19f2YFjtoa5gLPlYKKyOYETtif3S7ofVHoooFA/RjIip5g5DzR1G1t6P22FPCcbSK86KPLFgrTD+NXY/Zb9EuJa8T4gGABhdCLwAAMiSp97dpet/v1rBcFRHjivV9//lCKuHBCCF02FXqc+uUp8ra8c0ItFkUNbRazAWSdvenlJx1lWBlr5fKBIL0UKRqEId0eyGaE572nLN1Kq0REAWu58esvWoSksJ1pyEaACAHETgBQDAQdq6r10/eOpDPfVegyTprJlVuvdzx6rQw49ZIN+5hiBEC4Wj6gzFg7NQWG3BiDriFWWx+12BWldVWjitei0ZrMXDt2SIFo4qFI6qKYshmsdpT1aaJYKxom73uwdrvVWlxfqiOVTgIkQDABw8zsQBABiASNTUht1t2rC7Tdv2d2h9Y6v+8s5OGRFTdpt0xSmTddN5M/llDcCguZ12uZ12lRZkN0SLLdPsCs8SQVlqMJasOIvv15Zc3hm7n6hKaw+GFY6akqRgOKpgOKT97VkbbjJE66o469n/rK8+aH3t57DbsjdAAEDOI/ACACCDzlBEb2zdr+Ub9+mNLU16b2eLOkKRHvudPmOsvv2JQ3VITbEFowSAzGIhmltlBdk7ZiJEa+u2TDPRE61nn7Tu/dN6BnBDGaJ5XfYeQVnPqrSUiwekLOk8rLZUNaXe7A0GADDkCLwAAKOWaZqKRE1FTFPtwYhWbokFW+/uaJG/01BbMKy9bUEZETPteQVuh2bWFGtiZaHqKwp00pQKnTx1jEXvAgCs0RWiubNyPNM0FYpE05Zrpi/f7La8s5dgrXsA1x6KKBIP0QJGVAEjpH2DCNGKPU69cfNceZyOrLxXAMDQI/ACAIxKL6/fqxseW609rcED7ltX6tWcqWM0e0qFjq0v05SxRSyNAYAss9ls8jgd8jgdKi/sPUQzTVMdoYj2t4fU1BFK+WqoqT2k/R0hNSe2txva1x7SvvagTLPXwx1Qscep8kK3Tp5aKTdL1QFgRDmowOsHP/iBFi1apOuuu0733HOPJCkQCOjf//3f9bvf/U7BYFDnnHOO/ud//kfV1dXZGC8AAAdtxaZ9+sqvVypgRNO2Tx1bqDlTK3XCpAqNLfao2ONSRZFbdaVe2WwEXACQbQEjFl6lBVjtIe3v6AqwmtpDakq5HwpHD3zgXhS4HSovcKui0K3yQrcqClzxr26Vxb+WF7pUkdhW4JbbScgFACPVoAOvlStX6sEHH9RRRx2Vtv2GG27Q3/72N/3hD39QaWmprrnmGl188cV65ZVXDnqwAAAcrLe2NemKX8bCrjMOGasfXXq03E67XHa7fG6WqgDAYAXDETW1G2rqCKWFVfvj23qGWqEef3joL7fTrspCd+8BVmEsrEoNsMoL3PK6+B4PAKPJoAKvtrY2XX755frZz36m22+/Pbm9paVFP//5z/XII4/o4x//uCTp4Ycf1qGHHqrXXntNJ510UnZGDQDAIDT6A5r/i9fVHoro5KmVWvKF4/kFCAB6YUSi8eDK0P72+DLBAwRY7b1c0KM/XA5bV3CVDLBc8cCqa3t5SoDlczmovAUAZDSowGvBggU6//zzNXfu3LTAa9WqVTIMQ3Pnzk1umzlzpiZMmKDly5f3GngFg0EFg139U/x+vyTJMAwZhjGY4WWUOOZQHDsfMD+ZMT+ZMT+ZMT+ZDcf8PPzyJvkDYR1WW6z/+dzRcigqY5DVBcONz09mzEv20boif0Sippo7unpdpQZVzR29B1itgfCgXstht6m8wBUPp9wpoZWrWzVWohLLpSKPk/AKAJB1Aw68fve73+nNN9/UypUrezzW0NAgt9utsrKytO3V1dVqaGjo9XiLFy/Wrbfe2mP7s88+q4KCLF43uZulS5cO2bHzAfOTGfOTGfOTGfOT2VDNTygi/d+bDkk2nVTcrBeXPTskrzPU+PxgONC6Ird1hMJqaAn0aNbelOyFlb6ssKXTGFTTdptNKi+IhVIV/QiwygvdKvY4ZeeiHgCAHDCgwGv79u267rrrtHTpUnm93qwMYNGiRVq4cGHyvt/vV319vebNm6eSkpKsvEYqwzC0dOlSnX322XK5XFk//kjH/GTG/GTG/GTG/GQ21PPzh1Ufqf31DzS+zKsbL//YiLvKIp+fzAzD0JNPPmn1MPJCNltX9FXJj8wCRkQ7mju1fX+HPmrq1Pam2NeP4vf3tYcGddxSX6KfVUoFVnIZYc8Aq8TnGnHfKwEASBhQ4LVq1Srt3r1bxx13XHJbJBLRSy+9pJ/+9Kd65plnFAqF1NzcnFbl1djYqJqaml6P6fF45PF4emx3uVxDekI/1Mcf6ZifzJifzJifzJifzIZifkzT1K9f2y5Jmn/yZHk9vV/ufiTg84Ohls3WFX1V8o92oXBUO5s7U8KsDm3f36mP4sHW7tbgAY9R5HGmN2s/QIBV5nPJ6eCKgwCA0WNAgddZZ52ld999N23bl7/8Zc2cOVPf+ta3VF9fL5fLpWXLlumSSy6RJK1du1bbtm3TnDlzsjdqAAAGYPmmffqwoVU+l0OfmVVv9XCAnJXt1hV9VfLnu3Akql0tga7KrJTqrO1NHWrwBw64xLDQ7VB9RYHGl/s0vrzra31F7Gupj+AbAIBMBhR4FRcX64gjjkjbVlhYqMrKyuT2K6+8UgsXLlRFRYVKSkp07bXXas6cOVyhEQBgmYdf2SJJuuT4cSot4JdEoDdD0bqir0r+kS4SNbW7NZCsykp+jQdcu1oCikQzJ1pelz0ZZNUnvlZ03S8rcNHIHQCAgzCoqzRm8pOf/ER2u12XXHJJ2tV7AACwwsOvbNbSDxolSV86eZK1gwFy2FC0rhipTNPUntagtjd1LTNM9NP6qKlDO5o7ZUQyB1puh13jyn09KrMSgdaYIjeBFgAAQ+igA68XXngh7b7X69X999+v+++//2APDQDAoEWjpu585kM9+OImSdJXTp2saVXFFo8KyF2jqXWFaZra3x5Kawif2iB+R1OnguFoxmM47TbVlSUCrXiVVkWiWqtAVcUerlYIAICFsl7hBQBALvjOk+/pkRXbJEnfPOcQff2MqRaPCMht+dS6wjRN+TvDPRrCp1ZsdYQiGY9ht0m1pT6N623JYUWBqos9NIEHACCHEXgBAPLO5r3temTFNtlt0p2XHKVLaVQPZMVIaF3RGYroovtf0drG1gPuW13i6RFmjS8vUH15gWpKvXI7CbQAABipCLwAAHnnT29+JEk6bcZYwi7gIIzE1hVvf9ScDLvGFLm7+mZVpDeIryvzyetyWDxaAAAwVAi8AAB5JRo19ac3d0iSLjluvMWjATDcNu5pkySdechYPfzlEy0eDQAAsAp12gCAvPLa5n3a0dypYq9TZx9WbfVwAAyzjbvbJUlTxxZZPBIAAGAlAi8AQF7546pYddcnj6pjuRIwCiUqvKZWEXgBADCaEXgBAPJGezCsp97bJUn69PHjLB4NACskAy8qvAAAGNUIvAAAeeOp9xrUEYpo8phCHTeh3OrhABhmnaGIdjR3SpKmji20eDQAAMBKNK0HAOSsbfs6dN8/1mt3a1DBcETBcFShcFTBcFTBcKTrthFVKBJVJGpKki4+dpxsNpvFowcw3DbvbZdpSmUFLlUUuq0eDgAAsBCBFwAgJ/3l7Z1a9Kd31RYMD+h5lYVuXTqrfohGBSCXNfoDkqTx5T5CbwAARjkCLwBAzmgJSc9+0KjnPtyrx9+KNZ+fNbFcl504QR6nPfbP5ZDbYZfHZZfbYZfXZZfH6ZA7/niRxymngxX7wGjUaUQkSQUuTnEBABjtOBsAAFiuLRjW/J+/rlXbnNKqtyVJNpu04Ixpun7udAIsAP0SiAdeHhffMwAAGO0IvAAAlvvRM2u1aluzbDJ1SHWxjplQrn85dpxmT6m0emgARpBEhZfX5bB4JAAAwGoEXgAAS729vVm/Wr5FkvS1Q6Na+PmT5XK5rB0UgBEpYEQlST4CLwAARj3qvQEAlglHolr0p3dlmtKnjqrVzDLT6iEBGMECyQovTnEBABjtqPACAAyLSNTUvvag9rQG1dJhqD0U0asb9+qDXX6V+lz69nkztOKl7VYPE8AIFmRJIwAAiCPwAgBkXTRqqrnT0Oa97Xpx3R69sHa33tvRomgfBVzf/sRMVRZ5hneQAPJOMBxb0uhxUuEFAMBoR+AFAMiaVVv36xuPrtauls5ewy2bTaos9Ki8wKUCj1MFLoeOri/TpcfXKxIJD/+AAeSVqBn7xmO32SweCQAAsBqBFwAga3772jbtaO5M3q8sdOukKZU645CxOnnaGNWUeOWw9/6LaCQyXKMEkK8SQbuNwAsAgFGPwAsAkBWmaerlDXslSQ/+6/H6+MwquRwsKwIwfOIFXuojVwcAAKMIv4kAALJi45427W4NyuO06/QZYwm7AAw7ljQCAIAEfhsBAGTFy+tj1V2zJpVzhTQAljCTgZfFAwEAAJZjSSMAICte2bhPknTy1DEWjwTAaJXo4fXG1ib98pXNKvG5VOpzqcTnUonXpRKfUyVelwrcDvp8AQCQ5wi8AAAHLRyJ6rVNscDr1GkEXgCsUeCJVZe+unGfXo2H8L1x2G0q8Tp7BGFpt33dbhOYAQAwohB4AQAO2rs7WtQaCKvE69QR40qtHg6AUerKUyfL47BrT1tI/oAhf2f8XyAc/2rIiJiKRE01dRhq6jAG9TrdA7PSHuFYL2Fayn2fi8AMAIChRuAFADhoiUqKk6ZUykHzHAAWqSr2auG8Q/p83DRNBYxoVxgWMOTvDKfcD2fc3tJpKBw9+MDMabcdIBjruT11aabXZScwAwDgAAi8AAAHLdGw/tTpLGcEkLtsNpt8bod8boeqS7wDfr5pmuo0It3CsP6HZi2dhiJRU+Goqf3tIe1vDw3qfbgctn5Uk/VdZUZgBgAYDQi8AAAHJWBEtGpbkyQa1gPIbzabTQVupwrcTtWUHnxg1tLZLTTrHpT1EppFoqaMiKl97SHtG4rALHnb1W3ZZldo5nESmAEAch+BFwAgo0172vTORy2KxJfxbNvfofd2tuiDnX61BcMKR0yFIlFVl3g0dWyh1cMFgJyVjcCsIxTpEYS19NKrrK/QLGrqoAMzt8OerBor7meVWWpg5nU5BvW6AAAMBIEXAKBPwXBEFz/wqpr70afm4uPG8xd/ABhCNptNhR6nCj1O1Q7i+iCmaao9FOmjqiy9V1lvoVlrIBaYhSJR7W0LaW/bIAMzp/2ASzDHlfl03hG1cjvtg3oNAAAIvAAAfXrnoxY1dxgqcDs0a1KFHLZYU+gjxpfq8LoSjSn0yG6P/fJSVTzwagUAwPCx2Wwq8jhV5HGqTr4BPz8aNdUeCqcFY3tag/qoqVMfNXWkfQ2Go30eJxSOam9bUHvbgplf8DLpwmPGDXicAABIBF4AgAxWbIpdffHMQ6p0/+XHWTwaAEC2BYyuiq+Wfjbfb03ZHor0HWz1l8Nu61HhVVvq05yplVl4hwCA0YrACwDQpxWb90uSTpxcYfFIAAC9CRg9e3plCqu6bw9lqMTqL7tNPXt4ZVyyGHusNH6/wO1gSTwAIOsIvAAAvTIiUa3aGrv64uwpBF4AMBSC4cgBQ6lMTemzFVgVZwyr4k3nuwVWiduFBFYAgBxE4AUA6NV7O1rUEYqorMClGVXFVg8HAHJSfwKrntu77mfqddVfNptU7Em9ImLPwKrXaqv47UK3U3Y7gRUAIL8QeAEAepVYznjCpAp+EQKQt0LhaD/Cqr7Dq4CR/cCqX2FVYrvPpSICKwAAeiDwAgD06vV44DWb/l0ARoho1NQbW5u0pzXYI6xq6aW6KluBlSQVezOEUn1sL41XWRFYAQCQfQReAIAeIlFTK+OB10lTuEoWgJHhwZc26c6nPxzUcxMVVsV9hFWlGYKsIo9TDgIrAAByCoEXAKCHNbv8ag2GVexx6tDaEquHAwD9sqO5Q5I0vtynmTUlmZcCJnpdeV0q8hJYAQCQbwi8AAA9vLZpnyRp1qRyfgkEMOJcctx43XD2DKuHAQAALETgBQCQJIUjUb2zo0Wb9rTr8bd2SJJms5wRAAAAwAhE4AUAkCRd97vV+tu7u9K2nTJ1jEWjAQAAAIDBI/ACAGhtQ6v+9u4u2WzSyVMrNXlMoU6YVKEjx5daPTQAAAAAGDACLwCAHnxpoyTp3MNr9MAXjrd4NAAAAABwcOxWDwAAYK0dzZ368+qdkqSvnT7V4tEAAAAAwMEj8AKAUe5//7lJ4aipk6dW6uj6MquHAwAAAAAHjcALAEaxpvaQfvf6dknSv51BdRcAAACA/EAPLwDIc/vbQ1rb0KoNe9rUHgwrFI6qPRjWht1ten+nX51GRIfXlejUaVyREQAAAEB+IPACgDy12x/QZQ+9pk172zPu57DbtPDsGbLZbMM0MgAAAAAYWgReAJCnlry4KRl21Vf4dEh1sUp9bnlcdnmdDk0aU6BDqot1SE2xygrcFo8WAAAAALKHwAsA8lBTe0iPvr5NkvTLL5+gMw6psnhEAAAAADB8aFoPAHnoV8u3JHtznT5jrNXDAYBhtb2pQ6ZpWj0MAABgIQIvAMgzHaGwfvnqFkmxKy/SmwvAaHFsfbkk6U9v7tBX/2+VWjoNi0cEAACsQuAFAHnmd69vV3OHoYmVBTrviFqrhwMAw+bi48bp9ouOkNth17MfNOpTP31Z7+9ssXpYAADAAvTwAoARzjRNPf7WDq3csl97WkN6ffM+SdJXT5sqh53qLgCjh81m0xdOmqijxpfq337zprbu69C//M+ruu1Th+uzJ9RT8QoAwChC4AUAI1hDS0Df/H9v65/r96Ztryv16uLjxlk0KgCw1lHjy/S3b5yqf3/sbS37cLdu+tO7WrmlSbdfdIR8bofVwwMAAMOAwAsARpBo1NS2/R1a19iqDxta9fOXN6ul05DHadf8kydpQkWBxhR5dNzEMnld/FIHYPQqK3DrZ1+cpSUvbdSPnlmrP775kd7f2aL/ufw4TRlbZPXwAADAECPwAoARYOu+dv3xzR3605sf6aOmzrTHjhpfqrs/c4ymVfELHACksttt+voZ03RsfbmuffQtfdjQqk/99BXdeclROv8oehwCAJDPCLwAIAeta2zV//5zkzbuadfWfR3a2xZMPuZx2jWtqkjTq4p0/KQKXXZCvVwOrkECAH2ZM7VSf//Gqbr20be0YvN+LXjkTTX6D9MVp062emgAAGCIEHgBQA6JRk398tUt+sHTHyoUjia322zSqdPG6NPHj9e8w2roQQMAA1RV4tXnZ0/Qis37JUkrt+wn8AIAII8ReAFAjtjTGtTCx1YnG9CfcchYffr48ZpYUaiJYwpU4nVZPEIAGJmMSFQ/eOpD/fzlzZKkOVMq9b2LjrB4VAAAYCgReAFADli9vVlf+79VavAH5HHa9Z3zD9UXTpoom81m9dAAYETb3RrQNb99S69viVV2fe30qfqPeTPkZCk4AAB5jZ/0AGCxx1Zu12eWLFeDP6ApYwv112tP1b/OmUTYBWBYLV68WCeccIKKi4tVVVWliy66SGvXrk3bJxAIaMGCBaqsrFRRUZEuueQSNTY2WjTiA1u5Zb/Ov/dlvb5lv4o9Tj34r8frpvNmEnYBADAKDOinfT6eCAGAlV5ct0c3/vEdhSJRzTusWk8uOEXTq4utHhaAUejFF1/UggUL9Nprr2np0qUyDEPz5s1Te3t7cp8bbrhBf/nLX/SHP/xBL774onbu3KmLL77YwlH3zjRN/fzlzfrcQ69pT2tQh1QX68lrTtE5h9dYPTQAADBMBrSkMXEidMIJJygcDuvb3/625s2bpw8++ECFhYWSYidCf/vb3/SHP/xBpaWluuaaa3TxxRfrlVdeGZI3AAAjVTRq6q6nP5QkXXZCve74lyNlt1PVBcAaTz/9dNr9X/7yl6qqqtKqVat02mmnqaWlRT//+c/1yCOP6OMf/7gk6eGHH9ahhx6q1157TSeddFKPYwaDQQWDXVeZ9fv9Q/smJLUHw/rWH9/RX9/ZJUm68Jg6Lb74SBW46eQBAMBoMqCf/ENxIgQAo9XT7zfo/Z1+FXmcuvHcmYRdAHJKS0uLJKmiokKStGrVKhmGoblz5yb3mTlzpiZMmKDly5f3ep63ePFi3XrrrcMzYEnNHSF9eslybdjdJqfdpps/eZi+OId+iAAAjEYH9aeubJwI9fWXP8MwZBjGwQyvV4ljDsWx8wHzkxnzkxnzk1nq/IQjUf3omdiS8CtOnqhit23Uzxufn8yYn8yYl+yKRqO6/vrrdcopp+iII2JXM2xoaJDb7VZZWVnavtXV1WpoaOj1OIsWLdLChQuT9/1+v+rr64ds3H97d5c27G7TmCKPHvzX43T8xIohey0AAJDbBh14ZetEqK+//D377LMqKCgY7PAOaOnSpUN27HzA/GTG/GTG/GS2dOlSrdht06a9DhU6TY1rW6u//33tgZ84SvD5yYz5wXBYsGCB3nvvPb388ssHdRyPxyOPx5OlUR1YRzAiSTpt+hjCLgAARrlBB17ZOhHq6y9/8+bNU0lJyUEduzeGYWjp0qU6++yz5XK5sn78kY75yYz5yYz56V0kampfe0i7mtr15PMr5BgzSc80NEoK6ZqzDtHFp06yeog5gc9PZsxPZoZh6Mknn7R6GHnhmmuu0V//+le99NJLGj9+fHJ7TU2NQqGQmpub0/642djYqJqa3GgGH4pEJUkursIIAMCoN6jAK5snQn395c/lcg3pCf1QH3+kY34yY34yG+3zs31/h55cvUMf7PJrza5Wbd3XrqiZeNQhbdguSaqv8OnLp06Ry+WwbKy5aLR/fg6E+cFQMU1T1157rR5//HG98MILmjx5ctrjxx9/vFwul5YtW6ZLLrlEkrR27Vpt27ZNc+bMsWLIPRiJwMtJzy4AAEa7AQVe+XAiBABDyTRNfe03q/T+zvQrkdlt0pgijwoV0KmHTdCR48t19mHV8hJ2AcgRCxYs0COPPKInn3xSxcXFyXYUpaWl8vl8Ki0t1ZVXXqmFCxeqoqJCJSUluvbaazVnzpycuTCRQYUXAACIG1DglQ8nQgAwlN7d0aL3d/rldtr172fP0GF1JZpWVaSxRR6Z0Yj+/ve/6xOfOJQKHQA554EHHpAknXHGGWnbH374YX3pS1+SJP3kJz+R3W7XJZdcomAwqHPOOUf/8z//M8wj7ZsRiZXTugm8AAAY9QYUeOXDiRAADKXfrYwtVzzviBp99fSpaY8Z0YgVQwKAfjFN84D7eL1e3X///br//vuHYUQDFwrHKrze3NakB17YqAK3Qz63QwXxfz6Xs+u226ECd+y+x2mXzcYySAAA8smAlzQeSK6fCAHAUOkIhfXn1TslSZ89od7i0QDA6JNYJr5yS5NWbmnq9/PsNsnncsjn7grECuKBWG+BWdo2t1MFrp4hmi/5HAdhGgAAFhj0VRoBAOn+9s4utQXDmlhZoJMmV1o9HAAYda44dZI8TruaOkLqCEXUGYqoIxSO3TYiPbYF4xVhUVNqD0XUHhqaSlxfWiAWC8kKuwVmPUK0lMCsR/CWCNZcDtnthGkAAPSGwAsAsuSxN2LLGT8zq55fQADAAlXFXt1w9ox+7x+JmuoIheMhWDwQM8Jdt5Pb4/sY3UK0xONGRJ2hsNqDiWAtrIARTb5OpxHbrvbsv2eP096jGs3ncqjQE7/v6iVEO8ASz8TznPRCAwCMYAReAJAFG3a3aeWWJjnsNn36+PFWDwcA0A8Ou03FXpeKvdm/kEg0aqZXlRnh3kO0UFdIlv54t21G+v6JTiPBcFTBcFRNHUbW34PbYU+pSotXmLmcKvD0FZgdeIlngSsWxLmdhGkAgKFF4AUAB2F9Y6v+/m6Dnli9Q5J05iFVqi7xWjwqAIDV7HabCj1OFXqyf7ptmqYCRrTHcs20EC0RmBnpIVtH98d7We4ZjYdpoUhUoc6oWjqzH6Y57baUZZrOHss+C1Mr1npZ4tlXLzUfFyEAAMQReAHAILy2aZ/uXbZer27cl9xW4Hbo62dOzfAsAAAOns0WC4t8boey3THSNE0Fw9GUJZx9LPE0IrElnInH+1ju2b16LRxP08JRU62BsFoDYUnBrL4Hu009lnimLvss7G2JZ7deaWnVaSnP97oI0wBgpCDwAoB+ChgRPfN+g367Ypte37xfUuwv1KfPGKtzj6jR2YdVq6zAbfEoAQAYPJvNJq/LIa/LofIhOH4oGab1vsSzPdQtZOtluWdvz+8MRRSKdF2EoC0YVlswnPXx22zdLkLg6n5BgfiSTk969dpxE8t1TH1Z1scDAOgbgRcAHMD2/R164MWN+svqnWqNnzy7HXZ95oTx+trpUzW+vMDiEQIAMDK4nXa5nXaVKvt904xINHaBgG5LONuCYe1vC2lfe1D72kLam3J7X1tQe9tDCoWjB34BSaapZBg3EAVuh9757jwuBAAAw4jACwD60OgP6L5/rNfvV26XEYktwRhX5tMlx43T52ZPUG2pz+IRAgCQH8KRaI9+Y71dPTNtWy/LKXtr+h/sZ5h1sLwue7Kiq9CT3sB/ztRKwi4AGGYEXgDQjRGJ6uFXNusnS9fHLiMv6WPTx+jrZ0zT7MkVstvp3QEAGH0OtByxK3zqJZwyem+Y33054lDqvhwxrTF+L1ec7NlMv/crTib6fHF+AAC5hcALAFKs3t6sm/74jj5saJUkHTehTDeeO1MnTcl2W2AAALJrIA3nu1dJ9RVYtffScH4odW84X9DjCo00nAcA9A+BFwDE/fWdnbrh96tlREyVFbj07U8cqkuPH8/JMQAga0zTVMCIdjViN3pWSSUroIzel/B1D6za48/pNCKKDEMo5bTbeg+kUpbw+dy9hFOuXqqkUp7nczvkcRJKAQCyg8ALACQ9+vo2ffvxd2Wa0jmHV+uOfzlSlUUeq4cFALBANGomQ6UDLeHr7SqC3XtOtXfb3xz6TEpuhz2t4qnQ3e1qgr0s4esZWPW+hM/tpBcVACD3EXgByHvhSFS7W4Pa1dKpPa0hhaNRRaKmOkIRNbQEtHlvu/789k5J0udnT9D3LjxCDvpwAMCIE42aWrWtSfvbQz3DKSO9d1T3JXztwa7wKmAMT5Nzj9PeS5VUyvK8vpbw9RpYpYRTLgcN0gEAox6BF4ARZ/v+Dm3Y0yaZUtQ0tW1/h97f6dfahlZ1hMIylbhseOwXmPZQuF9/TV9w5lT9x7xDWEoBACPUL17ZrNv/tmbYXq/I41R5oUsVhR5VFLhUXuhWRYFbFUWxr2UFbhV5ulVWpTRD548rAAAMHQIvACPK5r3t+uS9/1R7KDKg5zntNtWUejW22COXwy6n3Savy6HqEq9qSrw6ZkKZTp8xdohGDQAYDsdPLNfR40vV3GkoYEQUMKIKGBEFw0NTsdUWDKstGNb2/Z197uN22OV12eV1OeL/Um875HXG7vtSHvPEb/tSn+N0yOt2xL6mHCPtefS/AgAgicALwIgRjkR1w+9Xqz0UUVWxR1UlHtlkU3WJR4fVleqw2hKVF7iS+xe4nSryOlXkcaqy0M3lwgEgzx07oVxPXnNqj+3RaOzqhQEjokA4tpwxYEQVCEdigZgRVacRSQvJAuGIAqGIAonnGRF1Gl23g/Hnd4bi+ybCNSOqUKQrYAtFYvf9gfCwzIHHGevd1T0YS96Ob481iE9/LD08c8SPk34MX0og53YQsAEAcheBF4AR46fPb9Dq7c0q9jr1xIJTVFfms3pIAIARwB6/qqDP7RiW14tEzWQwFghH4wFbRMGUYKwzNVyLV6El9ksEaJ1GRMFuIVwirAum7JN6ZcZgOBqvaDOG/H3abJK3WzDmcTnk6x609VqdlhqexZ7vczvSQrnuYZ2LvmQAgAEg8AIwIry1rUn3/WODJOn2i44g7AIA5CyH3aZCj1OFnuE51TYi0fTqNCO9gi2QUsEW7HY/fb9ejpES3CWCukRfTNOUOuPbhoPDbktWoXmc3SvT4tvjgZnPHQ/aeltGGg/hYgGbPa3SzZeyHz3WAGBkI/ACYDnTNNXgD6jRH1SjP6Cm9pAipinTlHa1dOrVjfv0zkctikRNXXB0nS48ZpzVQwYAIGe4HHa5HHYVe4f+tUzTVCgSjVWZ9QjP0peApoZmqRVtyeq05HLQzGFdQiRqxvumDf37lCSXw9azUi0esHlSQ7S06jR7vMqte9iWun/PsM7jtNN6AQCyjMALwLBavb1Z6xpbtac1qJ3NnVrX2KoPd7WqNXjg3ibHTSjT7RceMQyjBAAAvbHZbPI4Yz2+5HMd+AkHyTRT+q+lhWcpwVi3Pmrpy0G79utM6ckW7GVZaSAcVSjlAgdGxJQRCat1mPqvuZ32tD5qmfqw+dzx0C0ZwqXv1z1087kcGl9eQNUagFGFwAvAsDBNUz946kM9+NKmXh932G2qLvaoqsSrykK3HHab7Dabir1OnTi5QnOmVmp8ecEwjxoAAFjJZuuqshoOkaiZ1msttTotmOyjlnJhA6P3irZE77Zgt2Wl/k5DzR1G2oUNEkLxwK2l74t+HpTTZ4zVr644cWgODgA5iMALwJALR6L69uPv6rE3PpIknTKtUrWlPlWXeDStqkiH1pZoypgiuZ00owUAAMPDNE0ZETPz1TrjgVcwZZllZyh9yWWwt6qz3oKwcFf/s+Hmctg0oYI/HAIYXQi8AAyZvW1BvbZpn36/crv+uX6vHHabFl98pD4zq97qoQEAgBwUjkT77APWY0ljOB429bKksWuZYvrzgt2CrKhFAZTDbkv2/qIBPwAMDQIvAFm3rz2kf/9/q/TKhn3JbW6nXT/93LGad3iNhSMDAAADEY2aPUKj1Mqm/i7xSy4L7OXKkalVUmGLEiibTV3BkTMWKnl66Y3ldTrkdXfvr2WPN7LvrYl97324XA6q2gFgqBF4Aciq/UHpcz97XZv3dUiSZtYU6+SpY3TJ8eN0eF2pxaMDAGBk697EvbN7uNTPJu6pVVJ9NnE3or32mhounkRwlKF5u9eVWiXVvVIqHlx1C7J6hlR2uR122WxUQQFAPiHwAnBQOkJh7WoJqCMY0a7mdt3zrkMtRofGlfn0qytO0LSqYquHCADAkDFNU6FItFsvp/Qr/3VvXh7otl8wcQVBI703VI8gK37bKm6HvffgKaWyqUfY1GN5XeryvPRjpIdUBFAAgIND4AWgX/a3h/TYG9u1vz2k1kBYe1qDWr+7Vdv2d3RrwGrTtLGF+s1XTlJNqdeq4QIARjEj0ksvp7Sm4+mVTd2X1/VeJRVRZ2IZX0plVKdhXSNyp92WVqXUvXdTojeUr6+AypXeG6r7sdL6STnpAwUAGFkIvAD0yw+f+VCPvr6918eKPU4VeZ0qcDtUGm3Vg185UWMJuwAAw+yFtbt17aNvqTUQtuT17TalNQr3dGso3ntlU9c+npQqqe4VT71VUznpAwUAQJ8IvAAcUDAc0V/f2SVJ+uyseo0r96mswKVpY4s0o6ZYY4o8kiTDMPT3v/9dZQUuK4cLABilPtjl7xF2+bqFRb31hepRJZW4Ml78tqeXiqeu53Xt53LYWIYHAECOIPACcEAvrt2j1kBYNSVeLb74SNlZ0gAAyEEfmzZWd2mtvC67Xv/PuSr2OAmgAAAYpaiDBnBAf357pyTpk0fVEnYBAHLWEeNKVF3iUcCI6s2tTYRdAACMYgReADJqD4b13JpGSdIFR9dZPBoAAPpms9n08ZnVkqRla3ZbPBoAAGAlAi8AGT23plEBI6qJlQU6anyp1cMBACCjsw+rkiQtW9Mo06rLJwIAAMsReAHI6C/x5YyfOrqOpSEAgJx38tQx8rrs2tkS0JpdrVYPBwAAWITAC0CfGv0Bvbhuj6RY4AUAQK7zuhw6ddpYSdLPX96sgBGxeEQAAMAKXKURGKWiUVNtobA27WnXux81690dLdrbFlJbMKy2QFg7WzrV3GFIkmbWFGt6dbHFIwYAoH/+5dhxem5No/745kd6Y+t+ffeCw5K9vQAAwOhA4AWMAr9dsVV/XPWRWgPhZKDVFgqrP61NxhZ7dN1Z04d+kAAAZMn5R9UqHD1Gd/x9jbbu69AVv3xDZ82s0i0XHKaJlYVWDw8AAAwDAi8gz927bL3uXrquz8fLC1w6anyZjh5fqvHlBSrwOFTodqqm1KsJFQUq9PBtAgAw8lx4zDiddWi17vvHev38n5u17MPd+uf6vbr6tCn6+plTVeDm5xsAAPmMn/RAHksNu75+xlSdOm2MirxOFXtdKvI4Vex1yuO004weAJCXijxOLTrvUF16fL1u/cv7+uf6vfrp8xv0pzc/0nc+eZjOO6KGn4EAAOQpAi9gBPMHDH2w068Pd/m1ty2k/R0hNbWHtK89pH1tQW3c0y5JuvHcQ/T1M6ZZPFoAAKwxrapIv77iRD3zfqO+99cPtKO5U1//7Zs6ZVql/uuCw+lTCQBAHiLwAkaAcCSqDxtatXp7szbsbtPmve3avLdd2/Z3HPC53zp3pv7tjKnDMEoAAHKXzWbTuUfU6PQZY/XAixu15MWNemXDPp333//Ul0+ZpG+cNV3FXpfVwwQAAFlC4AXkKCMS1bI1jfr9yu1asXm/OkK9X1Z9XJlPh9WVqK7Uq/JCtyoS/wrcqq8oUH1FwTCPHACA3OVzO7Tw7Bn69HHjddtfP9Bzaxr1s39u1hOrd+qBy4/TrEkVVg8RAABkAYEXYDHTNNXoD2rjnjZt3dehhpZO7WgO6MV1e7S3LZjcr9jj1DETynRYbYkmjSnUpMpCHVJTrIpCt4WjBwBgZJpQWaCrT5uiNbv82tHcqT2tQf3prR0EXgAA5AkCL2AY+QOGnv9wt17btE87mgOxcKupU+19VG+NKfLoM7PG61PH1GlGVbHsdhrrAgBwsN7Ysl8/eW6dXtmwT5Lkctj0mVn1+o95h1g8MgAAkC0EXsAQMU1Tq7Y2ac0uv7bs69Dahlat2LxPRsTssa/DbtPEigJNGlOo2lKv6sp8mlFdrDMOGSuXw27B6AEAyD+rtjbpnufW6Z/r90qSnHabLp1VrwVnTtX4cloAAACQTwi8gCEQCkf1nSfe1WNvfNTjsaljCzX3sGpNHVuk2lKvakt9mlBRILeTYAsAgKHw1rYm/eS59Xpp3R5JiaBrvL5+xjR6XQIAkKcIvIAs29cW1L/95k29vmW/7DbpzEOqNHlMoSaOKdScKZWaVlVk9RABABgVVm9v1j3PrdMLa2NBl8Nu06ePG69rPk7QBQBAviPwAgZgT2tQb+yxad9r29QWiqql00j75+80tKO5U62BsIo9Tt33+WN1xiFVVg8bAIBR5Z2PmvWTpev0fErQdclx43TNmdM1oZKgCwCA0YDACxiAK3/9ptY0OKQNH2bcb1Jlgf53/ixNqyoeppEBAIB9bUHd+P/e0bIPdye3lfpcuvq0KTqstkQ7WzrV3BlSgdupArdDBW6HfG6H3A67bDYuDAMAQD4h8AL6acvedq1paJVdpuYdXqPyQrdKfC6V9vJvZk0JPbkAABhmv1u5PS3skqSWTkM/fGZtxuc57Tb54gFYehjmVIHLoQJP12M+l0OFnq7HkrfdjvhjzmSQVuByyMnFZwAAsASBF9BPz61plCRNLTF132VHy+VyWTwiAACQ6rMn1MuIRLW/PaSOUEQdoXD8a0SdoYjaQ2F1xu93hMLJKyeHo6ZaA2G1BsKSglkdk9tpj4VlrlgIVuiJhWYFbocKPPFALR6uFSaCspTQrcDtTIZxhSm3fS6H7Haq0gAA6AuBF9BPy9bE/mJ8RIVp8UgAAEBvxhR5dP3cGf3e34hE+wzDuoKy2O32lNuJ/RO3U5+TeCwaP10IhaMKhaNqlpH19+tzdS3L7F6dlrjdV+VaX+Fagdshj5MlngCAkY/AC+iHlg5Dr2/ZL0k6opzACwCAfOBy2FXqs6vUl92qbdM0FQxHY2GYEVFHMCVAM8JqDyaCsnD88a7HOkKR2OPx2x3BiDqM1DAuknydTiOiTiMitWd1+LLblFZZlhaWxZdtJpZsJivV4o8VuJ2xJaCuruq0Qo9DBa7YbVo+AACGC4EX0A8vrNutSNTU1LGFGuNtsXo4AAAgh9lsNnldDnldDpVn+djRqKlAOHLASrNEdVp7KCVcS9mnM1m11vVYMByNvYYptQXDaguGszx6yeWwdQVj7nh/NFcf4Vq3ZZw9K9fSK9UcLPEEAKQg8AL6IbGc8ayZY6UwgRcAALCG3W6Lhz7ZP42PRM1kGNbRY5ln+rLNHuFaauWaEU5WrSX2CcfXeBoRU0YkLH8g+2GaJ9EvLVFZlhKGJe4nK9dcXZVpPSrXeumXxhJPABh5CLyAAzAiUb2wNhZ4ffyQsWp8f4PFIwIAAMg+h92mYq9Lxd7sX5gnlFzimbJUM76kszMUUXswrE4jJVwLZngsLXjr6pcWDEcVDEfV1JH9fmmpVWcFLmfyyp0+Vyw0K/R03U6Eax6XQy6HXS6HLf71wLfdTrucdptcTrvc8ceoXAOAwSHwAg7gjS1N8gfCqih065j6Mj3zvtUjAgAAGFnczliYU6qh6ZfWfalmR58XIch0gYL0xzqNrn5p3funDSebLdZvzh0Pxpwpt10Oe/x+V3jmdNiSYZnLaZfLHn/MeeDgzR1/fvJ+/DlOe9ftns/veTwq4gDkAgIvoBsjEtXSDxr1jw93q7kjpA272yRJZxwylr+wAQAA5JDUfmkVhe6sHjsaNZOVZYnqtNQLDnQasQsMJMO11AsUGBEFjajC0aiMSFRG2FQoEr8fNmVEorH7ka7bRvx+Yvlngml2Xe1zpHAmQjaHLV611hW49QjVUu6nhnqulCAvLdSLV8G5nX0HblTRAZAIvDBKBYyI/AFDgVBUHUZY+9tD2tsW0obGVv3+je1q9Ad7POeCo+osGCkAALnn/vvv1w9/+EM1NDTo6KOP1n333acTTzzR6mEBWWW321TocarQM7y/MkWjpoxoVEbEVDgZhpkywrHALBQPzIzE9uTt9Puh+PMT20PhRAAXu230EbiFUp4Tjh/H6O3xcFRG1Ow1iAtHTYWjEXVmf3XpkLC6is5pz1xRl/Z6VNEB/UbghbxlmqZWbN6vNbv8ag2E1RowtHVfh9Y1tmrr/g6ZZt/PHVPk1iXHjdfEykIVe50aV+7TcRPKZRgj5Kc2AABD5Pe//70WLlyoJUuWaPbs2brnnnt0zjnnaO3ataqqqrJ6eMCIZ7fb5LE7NMw526CZpqlI1IyFbdF4EBbpGcr1qGgLRxWOxu/HnxML9Lqe0z1wyxTwGeFEUNhVRWf0UVGXb1V0ySCMKjogzQj5Ngr03+7WgP68eqceWbFNm/a297mfzab4ZbEdKitwq7LQrbHFHp19WLXOPaJGHqdjGEcNAMDIcPfdd+uqq67Sl7/8ZUnSkiVL9Le//U2/+MUvdNNNN6XtGwwGFQx2VU37/f5hHSuAoWez2eR02OR0SD6NjPPnRBVdeoVbV0VcoopuIBV13W8nK+oOsITVSAn1wonndquoo4qOKjoMDoEXcpZpmsm/9CR+ACX+ChSOb+8MRbWrpVO7WgJa29iqFZv2aeOerpCr0O3Qx6aPVXlh7IpDNSVeHVJTrOnVRRpb5OGbGAAAAxAKhbRq1SotWrQouc1ut2vu3Llavnx5j/0XL16sW2+9dTiHCAAHNFKr6MJRM1kdRxXd0MtUReey9x2+9ayw672iLtHjjiq6oTNC/hfHaLNhd6uu/NUb2rqvY8DPtdmko8aV6jMn1OvCY8apaKT8JAMAIMft3btXkUhE1dXVadurq6v14Ycf9th/0aJFWrhwYfK+3+9XfX39kI8TAPJJahWd10UV3VBU0RmRaI+WN1TRDayKzu1w6NgJZcPe9zCTIRsJzUwxWO3BsL72mzd7DbscdlsyaU/8D1pb6lVtqU/1FT6dMKlCsydXqrQgu5e8BgAAA+fxeOTxeKweBgBgmOVrFV1fFXWJKroeFXUpF5EwIlEF4ld+bQ+Gk1eBbU9c3TUUUWcorA4jkrHfdN/vwfoqupOnVuqRq06y7PW7G5KPH81MMVimaWrRn97Vht1tqi7x6P997WSNKfLE/qJhZx01AABWGjNmjBwOhxobG9O2NzY2qqamxqJRAQBGo7QWOOHUK4we4AqmKdVlqc9Lv0ppL8s5+1lF1usVTeMXdRhMkJVLHHZbrNLLHq8KS7sYgl3nHZFb5wJDEngNpJkpkGCaph5+ZYv+/PZOOew23f/541RfUWD1sAAAQJzb7dbxxx+vZcuW6aKLLpIkRaNRLVu2TNdcc421gwMAHJTEMkQjYsZCmtTbvYY9ZrxPWDzsie9nRFNux4OgcEqIFOq+tDC+bzhqxiujDlxdlbg90qX2CUv083J2C5G6X2GzxxLDeB+vRE+vvvqG9bmM0dl1PHdqf7GU/RLjHGl9wrIeeA20mWlfV+8xDEOGkf3FsoljDsWx88Fwz09bMKz3dvj1wro9evaD3dre1ClJunHedB09rjjn/jvx+cmM+cmM+cmM+cmM+cmMeRk+Cxcu1Pz58zVr1iydeOKJuueee9Te3p78QycAICYS7WUZXB/hTjg1ROq1Kqn3JXI9m84fuNF8b1VKRiS2pG+kS4ZCzviVGB1dgY7TntokvpdG8D36UsV7Ydm7bvfe26pbfyuuBJkzsh54DbSZaV9X73n22WdVUDB01T1Lly4dsmPng4Odn2BEajOkjrDUFrapPeV2mxF7bHenTQ2dkqmu/8ldNlOn1Jiqbv5Af//7Bwf7NoYMn5/MmJ/MmJ/MmJ/MmB9Y7bOf/az27NmjW265RQ0NDTrmmGP09NNP9zj3A4BsSvRY6t5wPC20SYQ7Kf2XMi9B6/1qg732Z+oeNvUaXKUvZRvp+VH3Jui9hj3OeDWRvet22lUI+wiFXM74sjhHalVSz+qmtObrfVwhMfE8WuCgO8tbyPV19Z558+appKQk669nGIaWLl2qs88+Wy4Xjc27O9D8mKapDbvbtbstmPyhkviG3h6K6P2dfq3a2qxNe9v7/Zp1pV7Nmliusw+r0mnTK1Xgtvxj2Sc+P5kxP5kxP5kxP5kxP5kZhqEnn3zS6mGMGtdccw1LGIERLtH/qPfeRT17HoWj6f2PDtjzqI+r6PV8ntnzueH4sriU/kehiHWNuLMl2f+oj3An7X4fVUI9r5bXtXQtbVlcWjVT9+qiXqqbnD33G2nL14Dusp4sDLSZaV9X73G5XEN6Qj/Uxx/pUucnYET07AeNWramUa9s2Ke9bcEDPDvG7bSrosCtsgKXygvcKi+Mfa0s8qiy0K3aUq+OqS9TVYl3KN/KkODzkxnzkxnzkxnzkxnzAwAjhz9gaH9bKK3nUdrStT56HhndKosSS9e69zxKf17/riKXbPKdB/2PEmHNQMOdZO+ilP5HLkfvPY+6VzOlPq+vnkc9+iTFm3zbCZCAYZX1wItmpiNf1JQ+aurUTn+Llq3ZrT+99ZGaO7p6o/hcDk2oKEiuY078cPE4HZpeXaTjJ5TruInlKi9wUVIKAACAUWnD7lZ94t6XFQqPnMqk1P5HsaVjqQ2tU6uJ+q5ScjkP3POoZ6+kvqqZ0iub0p/H8jUAmQ3J2jGameY+0zS1aW+7tuxt15Z9Hdq2r11b93do6952bdvvUOS1f6btX1fq1UXHjtNpM8bq2All8jgdFo0cAAAAyH2FHqdqSrza3tQh08JiqgK3Q6U+l0q8rthXn1MlPlfatlKfS8VepwrcTvncdvlcThW4HSpwO+RzO1TgdrK8DcCIMySBF81Mc9Oe1qBe3rBHL63bq3+u35thaWJsbXl9RYEOrS3Rp48fr9Omj+WHHAAAANBPtaU+vXTjmTJNUwEjqo5QWB2hiAJGRB2h2L9OI9x1O7ktos74voltHfFtnUa3fUORA/a2Shx/V0vgoN6P22mXz5UagjlU4HImb/tcXdt97pTAzOXIGKQVuB3yOO1UawHIuiHrDk4z0+EVMCL6qKlDLZ2GItHYJXCbOkLa1RLQR00dem3Tfq3Z5U97jtdl15QxRZpYWaCJlYWaWFmgcaVubXx7hT534XnyetwWvRsAAAAgP9hsNvniAU/lEBw/HImqw4goEOp/kJbc3keQFrsdVocRSVanhcKxHmAtnUbmAQ2Czab0MK2PIC0WnDlU4OoKzDIFad74MV0Oe9bHDCD35e7l8HBALZ2GfvzsWi1bs1s7Wzr7VSp9xLgSfWz6WJ02fayOm9hzaaJhGGr6UFRzAQAAACOA02FXicOuEm/2L2himqaC4WhaOHbAIM2I9F7NFoptTxyrIxRJ9jczza5KtKHgctiSlWaJUKyvIC1RoeZz2bu2dwvSUp/rdTpoRg/kKAKvEerp93bp5iff157WrmWJRR6nKovccthsstttKvE6VVvmU12pV0eMK9Up08ZoTFHPK2ICAAAAQHc2m01eV6xSqnwIjh+JmsmALL26LNJje2qQ1vV4V5DWkbItEA/VItFYRUDsqpVh+QPhIXgX6dVp6cs+ewvSui0LTVkCmhbCxavc3E6q04DBIvAaIX6/cpseemmTWgOxb/qtwdg36yljCvWdTx6qo8aXqbLQzdp3AAAAACOCw25TkcepIk/2fy01TVOhSPSAQVpHt55pmYK01OcHjK7eaZ1GbLvas/425LTbugVp3fujHSBIi4doffVdozoN+YzAK8eZpql7l23QT55bl7bdabfpq6dP0bUfny6viysmAgAAAECCzWaTx+mQx+lQWUH2jx9NVqcl+qGF5e8My99pqKXTkD8Q/9oZTt73dxryB2L7+DuNZBFDJuGoqdZAWK1DVJ02e3KFHr3qJIIv5CUCrxzUEv8G2GlE9MiKbfrlq1skSdecOU3nHVkjn8uhyiKPSn3ZX6cPAAAAAFYwTVORqKlwNFadFY6YMiLR+D9T4Ui02/bY13A0qlDYVDgaTd/efb9IVKH4cYxIVEbUlBGOprxe1769Pj9qKhSOxl8nfb9wtB8NlXNQoz+gqGnKLgIv5B8CLwuEI7GrmzR1GGr0B7SjuVMfNXVqzS6/3tvR0uslg//rgsP0pVMmWzBaAAAAACNJNGrKiPYW+MSDnWhURji2T3rg0zNg6m/g01cQlXydTPslbkej/boQ10jhctjkctjltNvkdtrltNvlctrksttj2+OPJ/dz2OV22OL72eWy2zLv57Cnb099neRzej4/9ThjizxychVL5CkCr2G0cU+b/uMPb+utbc0H3Nfrssvncqi80K3r587Qp46uG/oBAgAAAJAUqzYKR82UkKZbyJMMfw4QEnV/fjwECkf7qFY64Ov0sV8yzDKTzdrzgcNuiwU5acFNLLTpO/A5wH6O+PHSAihbPEjqPWBKC5AyvE7i+U67jf7KgMUIvIbJk6t36Nt/elftKZfaLfY6VVXs0bjyAo0r82paVbGOHFeqw+pKhqRxIwAAADDcotGuIMYIR+PVPj2DoZ7L0OL7RaN9Bj6hPgKm9OVqB9ovJYDq9jr5pPcgJ34/Hvw47QcIfFL3Swt8Dlx5FKtYOlCwlAi2ul6P3lIABotUZQi1Bgw9v3aP/vL2Ti39oFFSrCngXZ8+SuPKfJSOAgAAoF8S1UaZgqHUiqBYsHTgwCfUR/CTsUKpj8qjvvbLo2IjOey2lODnQIFPekVQ+hK1PoKhZOVRSkAUD5gyBkO9VR6lBEwOqo0AjEIEXlniDxh6+OUtWrqmQR2hiIJGVHtagwpFYpertdmka8+cpm+cNZ2gCwAAwCKRXpafDbQJdsZeR/1pgt1HENXnsaP5VW1ksyleVWSLh0Z9hETdq3/sdrkzhj/dlpf10gepe0jUW7DUZyUU1UYAMKIQeB2k7fs79MRbO/Szf26Sv5dLxU4ZW6hzDq/RJ4+q1eF1pRaMEAAAYHRYtXW/vvvn9+XvDHctd4s35TbiQVc+NcR29ghy+hv49NF3qJcqpO4BVK9h0AGDqPTjOAiNAADDgMBrEFo6Df3ombV6Yd1ubd/fmdw+rapIXzt9qurLfXI77aoodGtiZaGFIwUAABg9nly9U+/t8A/oOYlqo7QgJ1l51HdD6v4GPv29Aluvy9O6B0xpy9tYogYAQCYEXoPwpzc/0v+9tlVS7C9rx00o1xfmTNT5R9byFysAAACLROPlW5+ZNV5fOGliehDVLWCi2ggAgPxG4DUI6xrbJEmfnVWvWy44TIVcUREAACBn1Jb6dNT4MquHAQAALET39EHYtCcWeJ00tYKwCwAAAAAAIMcQeA3Cxj3tkqQpY4osHgkAAAAAAAC6I/AaoJZOQ3vbgpJiV2AEAAAAAABAbiHwGqDEcsaqYo+KvS6LRwMAAAAAAIDuCLwGaFN8OePUsSxnBAAAyEWm1QMAAACWI/AaoI3xCi+WMwIAAOQWh80mSYpGibwAABjtCLwGiAovAACA3GS3xwKviEngBQDAaEfgNUBUeAEAAOQmKrwAAEACgdcAhCNRbd3XIYkKLwAAgFzjSFR4EXgBADDqEXgNwEdNnQpFovI47RpX5rN6OAAAAEjBkkYAAJDgtHoAI8mmvbHljJPHFCZPqAAAAJAbEksa39/h1x9XfaSqEo+qir2qKvaorMAlm43zNwAARgsCrwHYuJuG9QAAALmq2Bs7tX19y369vmV/2mNuh11jiz0aW+xRVbEnLQxLvV1Z5EkujQQAACMXgdcAJCq8ptKwHgAAIOdcdsIEdRoRbdvfoT2tQe32B7W7NaCmDkOhSFQ7mju1o7kz4zHsNqmyKB6KJQMybzwU82hsSkjmcTqG6Z0BAICBIvAagESF1xQqvAAAAHJOaYFL18+d0WN7MBzR3raQdvsD2t0a1O7WoPak3N7dGtBuf1B724KKmtKe1qD2tAb1/oFez+fqUSE2ttijqhJvMjCrKvGqyMMpNwAAw42fvgOwcU+iwovACwAAYKTwOB0aV+Y74EWHIlFT+9pjlWF7UoKwZCjW2vVYKBJVS6ehlk5D63e3ZTxugdsRD8C8GlviSd5ODcvGFntUTp8xAACyhsCrn5o7QtrXHpIkTWFJIwAAQN5x2G3xIMqbcT/TNNXSaSQDsN2tgXhAFv/n77rfFgyrIxTRln0d2rKvI+NxXQ6bxhZ5NDa1QixlOWXidmWhW04HF1sHACATAq8MIlFTP3xmrd7e3pzs91BT4lUhZekAAACjls1mU1mBW2UFbs2oLs64b0conF4llnK7e58xI2JqZ0tAO1sCGY9pt0kVhanN93vvMza22COviz5jAIDRieQmg7++s1NLXtyYtu30GWMtGg0AAABGmgK3U5PGODVpTOYVAqFwVHvagv3uM7a3LXb7g12ZX7+vPmPnHVl7wCWeAACMZARefYhGTf3P87Gw63Mn1utTR49TdYlHkw9wsgIAAAAMlNtpH3CfsY+aOvXBzha9v9Ov93f61eDvWRnWV5+xl9bv1a+vODGr7wEAgFxC4NWHZR/u1trGVhV5nLrpvENV6nNZPSQAAACMAm3BcFqlV2pPsNRlkS2dRr+PabdJlUWeZG+wL548aejeAAAAOYDAqxemaeqnz2+QJP3rnImEXQAAADgopmmqqcPo1rur6+qPe/xdtztCkX4f1+2wa2x8mWL3pYuptytodA8AGGUIvHrxyoZ9ent7s7wuu648dbLVwwEAAECOCkei2tceSjafT71yY2ovrj1tQRkRs9/HLXQ7VFXi7QqyerlaY1WxR6U+l2w22xC+QwAARiYCr7gPdvr16sa92tMa1HNrGiVJl50wQWOKPBaPDAAAAMMtYESSywj39BZkxSu09rUHZfY/x1J5gSsZWI2Nh1ddoZZHVSWxiiyuCg4AwMEZ9T9Jg+GI/vu59Vry4kZFU05W3A67rj5tinUDAwAAQFaZphnrj5USXu1J6ZOV2jPLHwj3+7gOu01jitzJACsRXo0t8aYFWWOK3PI4HUP4DgEAQMKoC7zW727TmmabzHcb1BaK6jevbdWHDa2SpNNmjNX0qiKNLfboxMkVquNSzQAAACPG9v0d2ry3Pa25+57W9KqsTmMA/bGc9q7AqtuSwrEptysK3XLYWVYIAEAuGXWB1w+fXafn1zqkNe8kt1UUuvX9i47QeUfWWjgyAAAAHIz/XrZe/2/VRwfcr9jjTAusUhu8p/bMKvE56Y8FAMAINeoCr0mVhRpXsEfjqytU6nNrQkWBvnbGVHp1AQAAjHCTxxRqZk1x19LCXpq8VxV75XOzrBAAgHw36gKvb593iI4xN+oTnzhBLpfL6uEAAAAgSxacOU0Lzpxm9TAAAEAOsFs9AAAAAAAAACCbCLwAAAAAAACQVwi8AAAAAAAAkFcIvAAAAAAAAJBXCLwAAAAAAACQVwi8AAAAAAAAkFcIvAAAAAAAAJBXCLwAAAAAAACQVwi8AAAAAAAAkFcIvAAAAAAAAJBXCLwAAAAAAACQVwi8AAAAAAAAkFcIvAAAAEa5LVu26Morr9TkyZPl8/k0depUffe731UoFErb75133tHHPvYxeb1e1dfX66677rJoxAAAAJk5rR4AAAAArPXhhx8qGo3qwQcf1LRp0/Tee+/pqquuUnt7u370ox9Jkvx+v+bNm6e5c+dqyZIlevfdd3XFFVeorKxMV199tcXvAAAAIF3OBV6maUqKnVQNBcMw1NHRIb/fL5fLNSSvMZIxP5kxP5kxP5kxP5kxP5kxP5kl5kfqOpdA/5177rk699xzk/enTJmitWvX6oEHHkgGXr/97W8VCoX0i1/8Qm63W4cffrhWr16tu+++u8/AKxgMKhgMJu+3tLRIGrrzPAAAkJ8S5w4DOc/LucCrtbVVklRfX2/xSAAAwEjU2tqq0tJSq4cx4rW0tKiioiJ5f/ny5TrttNPkdruT28455xzdeeedampqUnl5eY9jLF68WLfeemuP7ZznAQCAwRjIeV7OBV51dXXavn27iouLZbPZsn58v9+v+vp6bd++XSUlJVk//kjH/GTG/GTG/GTG/GTG/GTG/GSWmJ8PPvhAdXV1Vg9nxNuwYYPuu+++ZHWXJDU0NGjy5Mlp+1VXVycf6y3wWrRokRYuXJi8H41GtX//flVWVnKel+OYy+xiPrOHucwe5jJ7mMvs6WsuTdNUa2vrgM7zci7wstvtGj9+/JC/TklJCR/EDJifzJifzJifzJifzJifzJifzMaNGye7nWvyJNx000268847M+6zZs0azZw5M3l/x44dOvfcc3XppZfqqquuOqjX93g88ng8advKysoO6pj9wf8n2cNcZhfzmT3MZfYwl9nDXGZPb3M50Ar+nAu8AAAAkB3//u//ri996UsZ95kyZUry9s6dO3XmmWfq5JNP1kMPPZS2X01NjRobG9O2Je7X1NRkZ8AAAABZQuAFAACQp8aOHauxY8f2a98dO3bozDPP1PHHH6+HH364R6XcnDlz9J//+Z8yDCN54YSlS5fqkEMO6XU5IwAAgJVGXc2/x+PRd7/73R7l9YhhfjJjfjJjfjJjfjJjfjJjfjJjfg7Ojh07dMYZZ2jChAn60Y9+pD179qihoUENDQ3JfT7/+c/L7Xbryiuv1Pvvv6/f//73+u///u+0Hl1W43OQPcxldjGf2cNcZg9zmT3MZfZkcy5tJtfuBgAAGNV++ctf6stf/nKvj6WeKr7zzjtasGCBVq5cqTFjxujaa6/Vt771reEaJgAAQL8ReAEAAAAAACCvjLoljQAAAAAAAMhvBF4AAAAAAADIKwReAAAAAAAAyCsEXgAAAAAAAMgroy7wuv/++zVp0iR5vV7Nnj1br7/+utVDssTixYt1wgknqLi4WFVVVbrooou0du3atH3OOOMM2Wy2tH9f+9rXLBrx8Pqv//qvHu995syZyccDgYAWLFigyspKFRUV6ZJLLlFjY6OFIx4+kyZN6jE3NptNCxYskDT6PjcvvfSSLrjgAtXV1clms+mJJ55Ie9w0Td1yyy2qra2Vz+fT3LlztX79+rR99u/fr8svv1wlJSUqKyvTlVdeqba2tmF8F0Mn0/wYhqFvfetbOvLII1VYWKi6ujp98Ytf1M6dO9OO0dtn7gc/+MEwv5OhcaDPz5e+9KUe7/3cc89N22e0fn4k9fq9yGaz6Yc//GFyn3z+/KAnzvMOXn/OETE4P/jBD2Sz2XT99ddbPZQRaceOHfrCF76gyspK+Xw+HXnkkXrjjTesHtaIE4lEdPPNN2vy5Mny+XyaOnWqvve974lr2R1YNs77EZON3xH6Y1QFXr///e+1cOFCffe739Wbb76po48+Wuecc452795t9dCG3YsvvqgFCxbotdde09KlS2UYhubNm6f29va0/a666irt2rUr+e+uu+6yaMTD7/DDD0977y+//HLysRtuuEF/+ctf9Ic//EEvvviidu7cqYsvvtjC0Q6flStXps3L0qVLJUmXXnppcp/R9Llpb2/X0Ucfrfvvv7/Xx++66y7de++9WrJkiVasWKHCwkKdc845CgQCyX0uv/xyvf/++1q6dKn++te/6qWXXtLVV189XG9hSGWan46ODr355pu6+eab9eabb+pPf/qT1q5dq0996lM99r3tttvSPlPXXnvtcAx/yB3o8yNJ5557btp7f/TRR9MeH62fH0lp87Jr1y794he/kM1m0yWXXJK2X75+fpCO87zs6O85IgZm5cqVevDBB3XUUUdZPZQRqampSaeccopcLpeeeuopffDBB/rxj3+s8vJyq4c24tx555164IEH9NOf/lRr1qzRnXfeqbvuukv33Xef1UPLedk470dMtn5HOCBzFDnxxBPNBQsWJO9HIhGzrq7OXLx4sYWjyg27d+82JZkvvvhictvpp59uXnfdddYNykLf/e53zaOPPrrXx5qbm02Xy2X+4Q9/SG5bs2aNKclcvnz5MI0wd1x33XXm1KlTzWg0aprm6P7cSDIff/zx5P1oNGrW1NSYP/zhD5PbmpubTY/HYz766KOmaZrmBx98YEoyV65cmdznqaeeMm02m7ljx45hG/tw6D4/vXn99ddNSebWrVuT2yZOnGj+5Cc/GdrB5YDe5mf+/PnmhRde2Odz+Pyku/DCC82Pf/zjadtGy+cHnOcNld7OETEwra2t5vTp082lS5eO6vOkg/Gtb33LPPXUU60eRl44//zzzSuuuCJt28UXX2xefvnlFo1oZBrMeT96N9jfEfpj1FR4hUIhrVq1SnPnzk1us9vtmjt3rpYvX27hyHJDS0uLJKmioiJt+29/+1uNGTNGRxxxhBYtWqSOjg4rhmeJ9evXq66uTlOmTNHll1+ubdu2SZJWrVolwzDSPkszZ87UhAkTRt1nKRQK6Te/+Y2uuOIK2Wy25PbR/LlJtXnzZjU0NKR9VkpLSzV79uzkZ2X58uUqKyvTrFmzkvvMnTtXdrtdK1asGPYxW62lpUU2m01lZWVp23/wgx+osrJSxx57rH74wx8qHA5bM0ALvPDCC6qqqtIhhxyif/u3f9O+ffuSj/H56dLY2Ki//e1vuvLKK3s8Npo/P6MF53lDp69zRPTfggULdP7556d9PjEwf/7znzVr1ixdeumlqqqq0rHHHquf/exnVg9rRDr55JO1bNkyrVu3TpL09ttv6+WXX9Z5551n8chGtv6c92Pw+vod4UCcQzOc3LN3715FIhFVV1enba+urtaHH35o0ahyQzQa1fXXX69TTjlFRxxxRHL75z//eU2cOFF1dXV655139K1vfUtr167Vn/70JwtHOzxmz56tX/7ylzrkkEO0a9cu3XrrrfrYxz6m9957Tw0NDXK73T3+Z6uurlZDQ4M1A7bIE088oebmZn3pS19KbhvNn5vuEp+H3r7vJB5raGhQVVVV2uNOp1MVFRWj7vMUCAT0rW99S5/73OdUUlKS3P6Nb3xDxx13nCoqKvTqq69q0aJF2rVrl+6++24LRzs8zj33XF188cWaPHmyNm7cqG9/+9s677zztHz5cjkcDj4/KX71q1+puLi4x/Ly0fz5GU04zxsafZ0jov9+97vf6c0339TKlSutHsqItmnTJj3wwANauHChvv3tb2vlypX6xje+Ibfbrfnz51s9vBHlpptukt/v18yZM+VwOBSJRPT9739fl19+udVDG9H6c96Pwenrd4T+GDWBF/q2YMECvffee2k9qiSl9YA58sgjVVtbq7POOksbN27U1KlTh3uYwyr1LxxHHXWUZs+erYkTJ+qxxx6Tz+ezcGS55ec//7nOO+881dXVJbeN5s8NBs8wDH3mM5+RaZp64IEH0h5buHBh8vZRRx0lt9utr371q1q8eLE8Hs9wD3VYXXbZZcnbRx55pI466ihNnTpVL7zwgs466ywLR5Z7fvGLX+jyyy+X1+tN2z6aPz/AwerrHBH9s337dl133XVaunRpj+9NGJhoNKpZs2bpjjvukCQde+yxeu+997RkyRICrwF67LHH9Nvf/laPPPKIDj/8cK1evVrXX3+96urqmEvknEy/I/THqFnSOGbMGDkcjh5X0mtsbFRNTY1Fo7LeNddco7/+9a96/vnnNX78+Iz7zp49W5K0YcOG4RhaTikrK9OMGTO0YcMG1dTUKBQKqbm5OW2f0fZZ2rp1q5577jl95StfybjfaP7cJD4Pmb7v1NTU9GioHA6HtX///lHzeUr8INu6dauWLl16wL/czJ49W+FwWFu2bBmeAeaQKVOmaMyYMcn/n/j8xPzzn//U2rVrD/j9SBrdn598xnle9g3kHBG9W7VqlXbv3q3jjjtOTqdTTqdTL774ou699145nU5FIhGrhzhi1NbW6rDDDkvbduihhyZbjqD/vvnNb+qmm27SZZddpiOPPFL/+q//qhtuuEGLFy+2emgjWn/O+zEwA/0doTejJvByu906/vjjtWzZsuS2aDSqZcuWac6cORaOzBqmaeqaa67R448/rn/84x+aPHnyAZ+zevVqSbEfOKNNW1ubNm7cqNraWh1//PFyuVxpn6W1a9dq27Zto+qz9PDDD6uqqkrnn39+xv1G8+dm8uTJqqmpSfus+P1+rVixIvlZmTNnjpqbm7Vq1arkPv/4xz8UjUaTYWE+S/wgW79+vZ577jlVVlYe8DmrV6+W3W7vsZRvNPjoo4+0b9++5P9Po/3zk/Dzn/9cxx9/vI4++ugD7juaPz/5jPO87BnMOSJ6d9ZZZ+ndd9/V6tWrk/9mzZqlyy+/XKtXr5bD4bB6iCPGKaecorVr16ZtW7dunSZOnGjRiEaujo4O2e3pMYDD4VA0GrVoRPmhP+f96L/B/I7Qm1G1pHHhwoWaP3++Zs2apRNPPFH33HOP2tvb9eUvf9nqoQ27BQsW6JFHHtGTTz6p4uLi5Lri0tJS+Xw+bdy4UY888og+8YlPqLKyUu+8845uuOEGnXbaaaPicsr/8R//oQsuuEATJ07Uzp079d3vflcOh0Of+9znVFpaqiuvvFILFy5URUWFSkpKdO2112rOnDk66aSTrB76sIhGo3r44Yc1f/58OZ1d30ZG4+emra0trXpt8+bNWr16tSoqKjRhwgRdf/31uv322zV9+nRNnjxZN998s+rq6nTRRRdJiv118txzz9VVV12lJUuWyDAMXXPNNbrsssvSloqOVJnmp7a2Vp/+9Kf15ptv6q9//asikUjye1FFRYXcbreWL1+uFStW6Mwzz1RxcbGWL1+uG264QV/4whfy4lLkmeanoqJCt956qy655BLV1NRo48aNuvHGGzVt2jSdc845kkb352fChAmSYieTf/jDH/TjH/+4x/Pz/fODdJznZceBzhHRf8XFxT16nxUWFqqyspKeaAN0ww036OSTT9Ydd9yhz3zmM3r99df10EMP6aGHHrJ6aCPOBRdcoO9///uaMGGCDj/8cL311lu6++67dcUVV1g9tJx3sOf96HKwvyP02yCvHDli3XfffeaECRNMt9ttnnjiieZrr71m9ZAsIanXfw8//LBpmqa5bds287TTTjMrKipMj8djTps2zfzmN79ptrS0WDvwYfLZz37WrK2tNd1utzlu3Djzs5/9rLlhw4bk452dnebXv/51s7y83CwoKDD/5V/+xdy1a5eFIx5ezzzzjCnJXLt2bdr20fi5ef7553v9f2n+/PmmacYuUXzzzTeb1dXVpsfjMc8666we87Zv3z7zc5/7nFlUVGSWlJSYX/7yl83W1lYL3k32ZZqfzZs39/m96PnnnzdN0zRXrVplzp492ywtLTW9Xq956KGHmnfccYcZCASsfWNZkml+Ojo6zHnz5pljx441XS6XOXHiRPOqq64yGxoa0o4xWj8/CQ8++KDp8/nM5ubmHs/P988PeuI87+Ad6BwRB+f00083r7vuOquHMSL95S9/MY844gjT4/GYM2fONB966CGrhzQi+f1+87rrrjMnTJhger1ec8qUKeZ//ud/msFg0Oqh5bxsnPcj5mB/R+gvm2maZv/jMQAAAAAAACC3jZoeXgAAAAAAABgdCLwAAAAAAACQVwi8AAAAAAAAkFcIvAAAAAAAAJBXCLwAAAAAAACQVwi8AAAAAAAAkFcIvAAAAAAAAJBXCLwAAAAAAACQVwi8AAAAAAAAkFcIvAAAAAAAAJBXCLwAAAAAAACQV/4/EQqA7c66zqwAAAAASUVORK5CYII=",
"text/plain": [
"<Figure size 1500x500 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"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": [
"<Figure size 1500x500 with 3 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"plt.figure(figsize=(15, 5))\n",
"plt.subplot(121)\n",
"plt.plot(np.sort(dh.eigh()))\n",
"plt.grid()\n",
"plt.subplot(122)\n",
"DOS = sisl.physics.electron.DOS(np.linspace(-15, 85, 50), dh.eig())\n",
"plt.plot(DOS, np.linspace(-15, 85, 50))\n",
"\n",
"coords = dh.xyz[-3:]\n",
"\n",
"\n",
"plt.figure(figsize=(15, 5))\n",
"plt.subplot(131)\n",
"plt.scatter(coords[:, 0], coords[:, 2], color=[\"r\", \"g\", \"b\"])\n",
"plt.xlabel(\"x\")\n",
"plt.ylabel(\"z\")\n",
"plt.subplot(132)\n",
"plt.scatter(coords[:, 1], coords[:, 2], color=[\"r\", \"g\", \"b\"])\n",
"plt.xlabel(\"y\")\n",
"plt.ylabel(\"z\")\n",
"plt.subplot(133)\n",
"plt.scatter(coords[:, 0], coords[:, 1], color=[\"r\", \"g\", \"b\"])\n",
"plt.xlabel(\"x\")\n",
"plt.ylabel(\"y\")\n",
"print(\"xyz[-3:]: red, green, blue\")"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Hamiltonian and exchange field rotated. Elapsed time: 199.591100666 s\n",
"================================================================================================================================================================\n"
]
}
],
"source": [
"NO = dh.no # shorthand for number of orbitals in the unit cell\n",
"\n",
"# preprocessing Hamiltonian and overlap matrix elements\n",
"h11 = dh.tocsr(dh.M11r)\n",
"h11 += dh.tocsr(dh.M11i) * 1.0j\n",
"h11 = h11.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype(\"complex128\")\n",
"\n",
"h22 = dh.tocsr(dh.M22r)\n",
"h22 += dh.tocsr(dh.M22i) * 1.0j\n",
"h22 = h22.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype(\"complex128\")\n",
"\n",
"h12 = dh.tocsr(dh.M12r)\n",
"h12 += dh.tocsr(dh.M12i) * 1.0j\n",
"h12 = h12.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype(\"complex128\")\n",
"\n",
"h21 = dh.tocsr(dh.M21r)\n",
"h21 += dh.tocsr(dh.M21i) * 1.0j\n",
"h21 = h21.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype(\"complex128\")\n",
"\n",
"sov = (\n",
" dh.tocsr(dh.S_idx)\n",
" .toarray()\n",
" .reshape(NO, dh.n_s, NO)\n",
" .transpose(0, 2, 1)\n",
" .astype(\"complex128\")\n",
")\n",
"\n",
"\n",
"# Reorganization of Hamiltonian and overlap matrix elements to SPIN BOX representation\n",
"U = np.vstack(\n",
" [np.kron(np.eye(NO, dtype=int), [1, 0]), np.kron(np.eye(NO, dtype=int), [0, 1])]\n",
")\n",
"# This is the permutation that transforms ud1ud2 to u12d12\n",
"# That is this transforms FROM SPIN BOX to ORBITAL BOX => U\n",
"# the inverse transformation is U.T u12d12 to ud1ud2\n",
"# That is FROM ORBITAL BOX to SPIN BOX => U.T\n",
"\n",
"# From now on everything is in SPIN BOX!!\n",
"hh, ss = np.array(\n",
" [\n",
" U.T @ np.block([[h11[:, :, i], h12[:, :, i]], [h21[:, :, i], h22[:, :, i]]]) @ U\n",
" for i in range(dh.lattice.nsc.prod())\n",
" ]\n",
"), np.array(\n",
" [\n",
" U.T\n",
" @ np.block([[sov[:, :, i], sov[:, :, i] * 0], [sov[:, :, i] * 0, sov[:, :, i]]])\n",
" @ U\n",
" for i in range(dh.lattice.nsc.prod())\n",
" ]\n",
")\n",
"\n",
"\n",
"# symmetrizing Hamiltonian and overlap matrix to make them hermitian\n",
"for i in range(dh.lattice.sc_off.shape[0]):\n",
" j = dh.lattice.sc_index(-dh.lattice.sc_off[i])\n",
" h1, h1d = hh[i], hh[j]\n",
" hh[i], hh[j] = (h1 + h1d.T.conj()) / 2, (h1d + h1.T.conj()) / 2\n",
" s1, s1d = ss[i], ss[j]\n",
" ss[i], ss[j] = (s1 + s1d.T.conj()) / 2, (s1d + s1.T.conj()) / 2\n",
"\n",
"# identifying TRS and TRB parts of the Hamiltonian\n",
"TAUY = np.kron(np.eye(NO), tau_y)\n",
"hTR = np.array([TAUY @ hh[i].conj() @ TAUY for i in range(dh.lattice.nsc.prod())])\n",
"hTRS = (hh + hTR) / 2\n",
"hTRB = (hh - hTR) / 2\n",
"\n",
"# extracting the exchange field\n",
"traced = [spin_tracer(hTRB[i]) for i in range(dh.lattice.nsc.prod())] # equation 77\n",
"XCF = np.array(\n",
" [\n",
" np.array([f[\"x\"] for f in traced]),\n",
" np.array([f[\"y\"] for f in traced]),\n",
" np.array([f[\"z\"] for f in traced]),\n",
" ]\n",
") # equation 77\n",
"\n",
"# Check if exchange field has scalar part\n",
"max_xcfs = abs(np.array(np.array([f[\"c\"] for f in traced]))).max()\n",
"if max_xcfs > 1e-12:\n",
" warnings.warn(\n",
" f\"Exchange field has non negligible scalar part. Largest value is {max_xcfs}\"\n",
" )\n",
"\n",
"if rank == root_node:\n",
" times[\"H_and_XCF_time\"] = timer()\n",
" print(\n",
" f\"Hamiltonian and exchange field rotated. Elapsed time: {times['H_and_XCF_time']} s\"\n",
" )\n",
" print(\n",
" \"================================================================================================================================================================\"\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Site and pair dictionaries created. Elapsed time: 199.638437125 s\n",
"================================================================================================================================================================\n"
]
}
],
"source": [
"# for every site we have to store 3 Greens function (and the associated _tmp-s) in the 3 reference directions\n",
"for i, mag_ent in enumerate(magnetic_entities):\n",
" parsed = parse_magnetic_entity(dh, **mag_ent) # parse orbital indexes\n",
" 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",
" mag_ent[\"energies\"] = [] # we will store the second order energy derivations here\n",
"\n",
" mag_ent[\"Gii\"] = [] # Greens function\n",
" mag_ent[\"Gii_tmp\"] = [] # Greens function for parallelization\n",
" # These will be the perturbed potentials from eq. 100\n",
" mag_ent[\"Vu1\"] = [list([]) for _ in range(len(ref_xcf_orientations))]\n",
" mag_ent[\"Vu2\"] = [list([]) for _ in range(len(ref_xcf_orientations))]\n",
" for i in ref_xcf_orientations:\n",
" # Greens functions for every quantization axis\n",
" mag_ent[\"Gii\"].append(\n",
" np.zeros((eset, spin_box_shape, spin_box_shape), dtype=\"complex128\")\n",
" )\n",
" mag_ent[\"Gii_tmp\"].append(\n",
" np.zeros((eset, spin_box_shape, spin_box_shape), dtype=\"complex128\")\n",
" )\n",
"\n",
"# for every site we have to store 2x3 Greens function (and the associated _tmp-s)\n",
"# in the 3 reference directions, because G_ij and G_ji are both needed\n",
"for pair in pairs:\n",
" # 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",
" 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",
" pair[\"Gji\"] = []\n",
" 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",
" )\n",
" pair[\"Gji\"].append(\n",
" np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype=\"complex128\")\n",
" )\n",
" pair[\"Gji_tmp\"].append(\n",
" np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype=\"complex128\")\n",
" )\n",
"\n",
"if rank == root_node:\n",
" times[\"site_and_pair_dictionaries_time\"] = timer()\n",
" print(\n",
" f\"Site and pair dictionaries created. Elapsed time: {times['site_and_pair_dictionaries_time']} s\"\n",
" )\n",
" print(\n",
" \"================================================================================================================================================================\"\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"k loop: 0%| | 0/400 [00:00<?, ?it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"k set created. Elapsed time: 199.64839675 s\n",
"================================================================================================================================================================\n"
]
}
],
"source": [
"kset = make_kset(dirs=kdirs, NUMK=kset) # generate k space sampling\n",
"wkset = np.ones(len(kset)) / len(kset) # generate weights for k points\n",
"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\")\n",
"\n",
"if rank == root_node:\n",
" times[\"k_set_time\"] = timer()\n",
" print(f\"k set created. Elapsed time: {times['k_set_time']} s\")\n",
" print(\n",
" \"================================================================================================================================================================\"\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Rotations done perpendicular to quantization axis. Elapsed time: 199.979076 s\n",
"================================================================================================================================================================\n"
]
}
],
"source": [
"# this will contain the three hamiltonians in the reference directions needed to calculate the energy variations upon rotation\n",
"hamiltonians = []\n",
"\n",
"# iterate over the reference directions (quantization axes)\n",
"for i, orient in enumerate(ref_xcf_orientations):\n",
" # obtain rotated exchange field\n",
" R = RotMa2b(scf_xcf_orientation, orient[\"o\"])\n",
" rot_XCF = np.einsum(\"ij,jklm->iklm\", R, XCF)\n",
" rot_H_XCF = sum(\n",
" [np.kron(rot_XCF[i], tau) for i, tau in enumerate([tau_x, tau_y, tau_z])]\n",
" )\n",
" rot_H_XCF_uc = rot_H_XCF[uc_in_sc_idx]\n",
"\n",
" # obtain total Hamiltonian with the rotated exchange field\n",
" rot_H = (\n",
" hTRS + rot_H_XCF\n",
" ) # equation 76 #######################################################################################\n",
"\n",
" hamiltonians.append(\n",
" dict(orient=orient[\"o\"], H=rot_H)\n",
" ) # store orientation and rotated Hamiltonian\n",
"\n",
" # these are the rotations (for now) perpendicular to the quantization axis\n",
" for u in orient[\"vw\"]:\n",
" Tu = np.kron(np.eye(NO, dtype=int), tau_u(u)) # section 2.H\n",
"\n",
" Vu1 = 1j / 2 * commutator(rot_H_XCF_uc, Tu) # equation 100\n",
" 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",
" )\n",
" mag_ent[\"Vu2\"][i].append(\n",
" Vu2[:, mag_ent[\"spin_box_indeces\"]][mag_ent[\"spin_box_indeces\"], :]\n",
" )\n",
"\n",
"if rank == root_node:\n",
" times[\"reference_rotations_time\"] = timer()\n",
" print(\n",
" f\"Rotations done perpendicular to quantization axis. Elapsed time: {times['reference_rotations_time']} s\"\n",
" )\n",
" print(\n",
" \"================================================================================================================================================================\"\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Starting matrix inversions\n",
"Total number of k points: 400\n",
"Number of energy samples per k point: 600\n",
"Total number of directions: 3\n",
"Total number of matrix inversions: 720000\n",
"The shape of the Hamiltonian and the Greens function is 84x84=7056\n",
"Memory taken by a single Hamiltonian is: 0.015625 KB\n",
"Expected memory usage per matrix inversion: 0.5 KB\n",
"Expected memory usage per k point for parallel inversion: 900.0 KB\n",
"Expected memory usage on root node: 351.5625 MB\n",
"================================================================================================================================================================\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"k loop: 100%|██████████| 400/400 [2:12:33<00:00, 19.88s/it] "
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Calculated Greens functions. Elapsed time: 1375.51451525 s\n",
"================================================================================================================================================================\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"if rank == root_node:\n",
" print(\"Starting matrix inversions\")\n",
" print(f\"Total number of k points: {kset.shape[0]}\")\n",
" print(f\"Number of energy samples per k point: {eset}\")\n",
" print(f\"Total number of directions: {len(hamiltonians)}\")\n",
" print(\n",
" f\"Total number of matrix inversions: {kset.shape[0] * len(hamiltonians) * eset}\"\n",
" )\n",
" print(f\"The shape of the Hamiltonian and the Greens function is {NO}x{NO}={NO*NO}\")\n",
" # https://stackoverflow.com/questions/70746660/how-to-predict-memory-requirement-for-np-linalg-inv\n",
" # memory is O(64 n**2) for complex matrices\n",
" memory_size = getsizeof(hamiltonians[0][\"H\"].base) / 1024\n",
" print(\n",
" f\"Memory taken by a single Hamiltonian is: {getsizeof(hamiltonians[0]['H'].base) / 1024} KB\"\n",
" )\n",
" print(f\"Expected memory usage per matrix inversion: {memory_size * 32} KB\")\n",
" print(\n",
" f\"Expected memory usage per k point for parallel inversion: {memory_size * len(hamiltonians) * eset * 32} KB\"\n",
" )\n",
" print(\n",
" f\"Expected memory usage on root node: {len(np.array_split(kset, size)[0]) * memory_size * len(hamiltonians) * eset * 32 / 1024} MB\"\n",
" )\n",
" print(\n",
" \"================================================================================================================================================================\"\n",
" )\n",
"\n",
"comm.Barrier()\n",
"# ----------------------------------------------------------------------\n",
"\n",
"# make energy contour\n",
"# we are working in eV now !\n",
"# and sisil shifts E_F to 0 !\n",
"cont = make_contour(emin=ebot, enum=eset, p=esetp)\n",
"eran = cont.ze\n",
"\n",
"# ----------------------------------------------------------------------\n",
"# sampling the integrand on the contour and the BZ\n",
"for k in kpcs[rank]:\n",
" wk = wkset[rank] # weight of k point in BZ integral\n",
" # iterate over reference directions\n",
" for i, hamiltonian_orientation in enumerate(hamiltonians):\n",
" # calculate Greens function\n",
" H = hamiltonian_orientation[\"H\"]\n",
" HK, SK = hsk(H, ss, dh.sc_off, k)\n",
" # Gk = inv(SK * eran.reshape(eset, 1, 1) - HK)\n",
"\n",
" # solve Greens function sequentially for the energies, because of memory bound\n",
" Gk = np.zeros(shape=(eset, HK.shape[0], HK.shape[1]), dtype=\"complex128\")\n",
" for j in range(eset):\n",
" Gk[j] = inv(SK * eran[j] - HK)\n",
"\n",
" # store the Greens function slice of the magnetic entities (for now) based on the on-site projections\n",
" for mag_ent in magnetic_entities:\n",
" mag_ent[\"Gii_tmp\"][i] += (\n",
" Gk[:, mag_ent[\"spin_box_indeces\"]][..., mag_ent[\"spin_box_indeces\"]]\n",
" * wk\n",
" )\n",
"\n",
" for pair in pairs:\n",
" # add phase shift based on the cell difference\n",
" phase = np.exp(1j * 2 * np.pi * k @ pair[\"Ruc\"].T)\n",
"\n",
" # get the pair orbital sizes from the magnetic entities\n",
" ai = magnetic_entities[pair[\"ai\"]][\"spin_box_indeces\"]\n",
" aj = magnetic_entities[pair[\"aj\"]][\"spin_box_indeces\"]\n",
"\n",
" # store the Greens function slice of the magnetic entities (for now) based on the on-site projections\n",
" pair[\"Gij_tmp\"][i] += Gk[:, ai][..., aj] * phase * wk\n",
" pair[\"Gji_tmp\"][i] += Gk[:, aj][..., ai] * phase * wk\n",
"\n",
"# summ reduce partial results of mpi nodes\n",
"for i in range(len(hamiltonians)):\n",
" for mag_ent in magnetic_entities:\n",
" comm.Reduce(mag_ent[\"Gii_tmp\"][i], mag_ent[\"Gii\"][i], root=root_node)\n",
"\n",
" for pair in pairs:\n",
" 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",
"if rank == root_node:\n",
" times[\"green_function_inversion_time\"] = timer()\n",
" print(\n",
" f\"Calculated Greens functions. Elapsed time: {times['green_function_inversion_time']} s\"\n",
" )\n",
" print(\n",
" \"================================================================================================================================================================\"\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Magnetic entities integrated.\n",
"Pairs integrated.\n",
"Magnetic parameters calculated.\n",
"##################################################################### GROGU OUTPUT #############################################################################\n",
"================================================================================================================================================================\n",
"Input file: \n",
"Not yet specified.\n",
"Output file: \n",
"test.pickle\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: 20\n",
"k point directions: xy\n",
"Ebot: -13\n",
"Eset: 600\n",
"Esetp: 10000\n",
"================================================================================================================================================================\n",
"Atomic information: \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",
"[5]Fe(2) 1.8954667088117545 1.0943913231921656 10.285002698393109\n",
"\n",
"================================================================================================================================================================\n",
"Exchange [meV]\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Magnetic entity1 Magnetic entity2 [i j k] d [Ang]\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"[3]Fe(2) [4]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -61.575578560021754\n",
"DMI: [-9.32925576e-01 -6.69311832e-04 -1.27860983e-06]\n",
"Symmetric-anisotropy: [ 3.27134566e-01 -2.18826996e-06 1.11817934e-06 -2.18826996e-06\n",
" 2.93842716e+00 -2.84299672e-05 1.11817934e-06 -2.84299672e-05\n",
" -3.26556173e+00]\n",
"Energies for debugging: \n",
"array([[-6.21988435e-02, -9.32897146e-04, 9.32954006e-04,\n",
" -6.14460702e-02],\n",
" [-6.74834371e-02, 6.68193653e-07, -6.70430011e-07,\n",
" -6.66686592e-02],\n",
" [-5.58282326e-02, 9.09660134e-10, 3.46687979e-09,\n",
" -5.58282288e-02]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.06666866, -0.05582823, -0.06219884])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.06666865917296813 -0.055828228814894104\n",
"\n",
"[3]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -60.54611612367399\n",
"DMI: [ 3.78484434e+00 -6.14206310e+00 1.63244853e-04]\n",
"Symmetric-anisotropy: [ 0.20562843 0.07121269 -0.0931207 0.07121269 0.01213322 -0.04251871\n",
" -0.0931207 -0.04251871 -0.21776164]\n",
"Energies for debugging: \n",
"array([[-6.08734388e-02, 3.82736305e-03, -3.74232563e-03,\n",
" -6.11763023e-02],\n",
" [-6.06543168e-02, 6.23518381e-03, -6.04894240e-03,\n",
" -6.08718006e-02],\n",
" [-5.98916636e-02, -7.10494409e-05, -7.13759306e-05,\n",
" -5.98091748e-02]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.0608718 , -0.05989166, -0.06087344])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.06087180064052234 -0.05980917475210674\n",
"\n",
"[4]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -60.54211382677383\n",
"DMI: [-3.79969677e+00 6.15441158e+00 1.64058302e-04]\n",
"Symmetric-anisotropy: [ 0.20695438 0.07120989 0.08986466 0.07120989 0.00574537 0.03644448\n",
" 0.08986466 0.03644448 -0.21269975]\n",
"Energies for debugging: \n",
"array([[-6.08700078e-02, -3.83614125e-03, 3.76325229e-03,\n",
" -6.11810143e-02],\n",
" [-6.06396194e-02, -6.24427624e-03, 6.06454692e-03,\n",
" -6.08610806e-02],\n",
" [-5.98917226e-02, -7.10458351e-05, -7.13739518e-05,\n",
" -5.98092383e-02]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.06086108, -0.05989172, -0.06087001])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.06086108060807331 -0.059809238286325714\n",
"\n",
"================================================================================================================================================================\n",
"Runtime information: \n",
"Total runtime: 1177.247012041 s\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Initial setup: 0.12478720800001497 s\n",
"Hamiltonian conversion and XC field extraction: 0.923 s\n",
"Pair and site datastructure creatrions: 0.047 s\n",
"k set cration and distribution: 0.010 s\n",
"Rotating XC potential: 0.331 s\n",
"Greens function inversion: 1175.535 s\n",
"Calculate energies and magnetic components: 0.276 s\n"
]
}
],
"source": [
"if rank == root_node:\n",
" # iterate over the magnetic entities\n",
" for tracker, mag_ent in enumerate(magnetic_entities):\n",
" # iterate over the quantization axes\n",
" for i, Gii in enumerate(mag_ent[\"Gii\"]):\n",
" storage = []\n",
" # iterate over the first and second order local perturbations\n",
" for Vu1, Vu2 in zip(mag_ent[\"Vu1\"][i], mag_ent[\"Vu2\"][i]):\n",
" # The Szunyogh-Lichtenstein formula\n",
" traced = np.trace((Vu2 @ Gii + 0.5 * Gii @ Vu1 @ Gii), axis1=1, axis2=2)\n",
" # evaluation of the contour integral\n",
" storage.append(np.trapz(-1 / np.pi * np.imag(traced * cont.we)))\n",
"\n",
" # fill up the magnetic entities dictionary with the energies\n",
" magnetic_entities[tracker][\"energies\"].append(storage)\n",
" # convert to np array\n",
" magnetic_entities[tracker][\"energies\"] = np.array(\n",
" magnetic_entities[tracker][\"energies\"]\n",
" )\n",
" print(\"Magnetic entities integrated.\")\n",
"\n",
" # iterate over the pairs\n",
" for tracker, pair in enumerate(pairs):\n",
" # iterate over the quantization axes\n",
" for i, (Gij, Gji) in enumerate(zip(pair[\"Gij\"], pair[\"Gji\"])):\n",
" site_i = magnetic_entities[pair[\"ai\"]]\n",
" site_j = magnetic_entities[pair[\"aj\"]]\n",
"\n",
" storage = []\n",
" # iterate over the first order local perturbations in all possible orientations for the two sites\n",
" for Vui in site_i[\"Vu1\"][i]:\n",
" for Vuj in site_j[\"Vu1\"][i]:\n",
" # The Szunyogh-Lichtenstein formula\n",
" traced = np.trace((Vui @ Gij @ Vuj @ Gji), axis1=1, axis2=2)\n",
" # evaluation of the contour integral\n",
" storage.append(np.trapz(-1 / np.pi * np.imag(traced * cont.we)))\n",
" # fill up the pairs dictionary with the energies\n",
" pairs[tracker][\"energies\"].append(storage)\n",
" # convert to np array\n",
" pairs[tracker][\"energies\"] = np.array(pairs[tracker][\"energies\"])\n",
"\n",
" print(\"Pairs integrated.\")\n",
"\n",
" # calculate magnetic parameters\n",
" for pair in pairs:\n",
" J_iso, J_S, D = calculate_exchange_tensor(pair)\n",
" pair[\"J_iso\"] = J_iso * sisl.unit_convert(\"eV\", \"meV\")\n",
" pair[\"J_S\"] = J_S * sisl.unit_convert(\"eV\", \"meV\")\n",
" pair[\"D\"] = D * sisl.unit_convert(\"eV\", \"meV\")\n",
"\n",
" print(\"Magnetic parameters calculated.\")\n",
"\n",
" times[\"end_time\"] = timer()\n",
" print(\n",
" \"##################################################################### GROGU OUTPUT #############################################################################\"\n",
" )\n",
"\n",
" print_parameters(simulation_parameters)\n",
" print_atoms_and_pairs(magnetic_entities, pairs)\n",
" print_runtime_information(times)\n",
"\n",
" # remove clutter from magnetic entities and pair information\n",
" for pair in pairs:\n",
" del pair[\"Gij\"]\n",
" del pair[\"Gij_tmp\"]\n",
" del pair[\"Gji\"]\n",
" del pair[\"Gji_tmp\"]\n",
" for mag_ent in magnetic_entities:\n",
" del mag_ent[\"Gii\"]\n",
" del mag_ent[\"Gii_tmp\"]\n",
" del mag_ent[\"Vu1\"]\n",
" del mag_ent[\"Vu2\"]\n",
" # create output dictionary with all the relevant data\n",
" results = dict(\n",
" parameters=simulation_parameters,\n",
" magnetic_entities=magnetic_entities,\n",
" pairs=pairs,\n",
" runtime=times,\n",
" )\n",
" # save dictionary\n",
" with open(outfile, \"wb\") as output_file:\n",
" pickle.dump(results, output_file)"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"ename": "SyntaxError",
"evalue": "invalid syntax (1445642317.py, line 1)",
"output_type": "error",
"traceback": [
"\u001b[0;36m Cell \u001b[0;32mIn[35], line 1\u001b[0;36m\u001b[0m\n\u001b[0;31m Magnetic entities integrated.\u001b[0m\n\u001b[0m ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
]
}
],
"source": [
"Magnetic entities integrated.\n",
"Pairs integrated.\n",
"Magnetic parameters calculated.\n",
"##################################################################### GROGU OUTPUT #############################################################################\n",
"================================================================================================================================================================\n",
"Input file: \n",
"Not yet specified.\n",
"Output file: \n",
"test.pickle\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: 20\n",
"k point directions: x\n",
"Ebot: -13\n",
"Eset: 600\n",
"Esetp: 10000\n",
"================================================================================================================================================================\n",
"Atomic information: \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",
"[5]Fe(2) 1.8954667088117545 1.0943913231921656 10.285002698393109\n",
"\n",
"================================================================================================================================================================\n",
"Exchange [meV]\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Magnetic entity1 Magnetic entity2 [i j k] d [Ang]\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"[3]Fe(2) [4]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -37.73293438752643\n",
"DMI: [-1.95914479e+00 6.21772915e+00 -7.63210027e-06]\n",
"Symmetric-anisotropy: [ 1.76549627e+00 1.32018786e+00 8.70141772e-04 1.32018786e+00\n",
" 3.02084625e+00 -1.61370176e-02 8.70141772e-04 -1.61370176e-02\n",
" -4.78634252e+00]\n",
"Energies for debugging: \n",
"array([[-0.03971497, -0.00194301, 0.00197528, -0.0357954 ],\n",
" [-0.04532359, -0.0062186 , 0.00621686, -0.03983102],\n",
" [-0.03362878, -0.0013202 , -0.00132018, -0.03210385]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.03983102, -0.03362878, -0.03971497])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.0398310234439899 -0.03210385279422295\n",
"\n",
"[3]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -73.83015029521252\n",
"DMI: [ 7.76660921e+00 -1.22515694e+01 6.25989582e-04]\n",
"Symmetric-anisotropy: [-1.94961899 -0.25694518 1.17973546 -0.25694518 -0.9315989 0.68240007\n",
" 1.17973546 0.68240007 2.8812179 ]\n",
"Energies for debugging: \n",
"array([[-0.06985729, 0.00708421, -0.00844901, -0.06974285],\n",
" [-0.07204057, 0.01107183, -0.0134313 , -0.07148291],\n",
" [-0.07978065, 0.00025757, 0.00025632, -0.08007663]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.07148291, -0.07978065, -0.06985729])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.07148290532391263 -0.08007663325482241\n",
"\n",
"[4]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -73.78915203163939\n",
"DMI: [-7.75054625e+00 1.22185627e+01 6.29159323e-04]\n",
"Symmetric-anisotropy: [-1.96064992 -0.25693115 -1.19384894 -0.25693115 -0.94861315 -0.70360485\n",
" -1.19384894 -0.70360485 2.90926307]\n",
"Energies for debugging: \n",
"array([[-0.06978756, -0.00704694, 0.00845415, -0.06969528],\n",
" [-0.07197222, -0.01102471, 0.01341241, -0.07142338],\n",
" [-0.07978025, 0.00025756, 0.0002563 , -0.08007622]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.07142338, -0.07978025, -0.06978756])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.07142338168884173 -0.08007622222171976\n",
"\n",
"================================================================================================================================================================\n",
"Runtime information: \n",
"Total runtime: 60.559684874999995 s\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Initial setup: 0.13768908300000016 s\n",
"Hamiltonian conversion and XC field extraction: 0.822 s\n",
"Pair and site datastructure creatrions: 0.036 s\n",
"k set cration and distribution: 0.025 s\n",
"Rotating XC potential: 0.245 s\n",
"Greens function inversion: 59.066 s\n",
"Calculate energies and magnetic components: 0.228 s"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [
{
"ename": "SyntaxError",
"evalue": "invalid syntax (3012513526.py, line 1)",
"output_type": "error",
"traceback": [
"\u001b[0;36m Cell \u001b[0;32mIn[13], line 1\u001b[0;36m\u001b[0m\n\u001b[0;31m Magnetic entities integrated.\u001b[0m\n\u001b[0m ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
]
}
],
"source": [
"Magnetic entities integrated.\n",
"Pairs integrated.\n",
"Magnetic parameters calculated.\n",
"##################################################################### GROGU OUTPUT #############################################################################\n",
"================================================================================================================================================================\n",
"Input file: \n",
"Not yet specified.\n",
"Output file: \n",
"test.pickle\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: 1000\n",
"k point directions: x\n",
"Ebot: -13\n",
"Eset: 500\n",
"Esetp: 10000\n",
"================================================================================================================================================================\n",
"Atomic information: \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",
"[5]Fe(2) 1.8954667088117545 1.0943913231921656 10.285002698393109\n",
"\n",
"================================================================================================================================================================\n",
"Exchange [meV]\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Magnetic entity1 Magnetic entity2 [i j k] d [Ang]\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"[3]Fe(2) [4]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -29.367092812942317\n",
"DMI: [-8.89792781e-01 4.16580540e+00 -3.59948517e-06]\n",
"Symmetric-anisotropy: [-4.56670448 1.17667365 0.33235705 4.15578143 -0.03089926]\n",
"Energies for debugging: \n",
"array([[-0.02819042, -0.00085889, 0.00092069, -0.02597706],\n",
" [-0.03690207, -0.00415578, 0.00417583, -0.0339338 ],\n",
" [-0.03160427, -0.00033236, -0.00033235, -0.03122101]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.0339338 , -0.03160427, -0.02819042])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.03393379729617979 -0.03122101431812066\n",
"\n",
"[4]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -74.98944973858333\n",
"DMI: [-6.82212886e+00 1.12891691e+01 5.77778279e-04]\n",
"Symmetric-anisotropy: [-1.03411109 0.63716653 -0.16031601 10.67952816 -0.32649278]\n",
"Energies for debugging: \n",
"array([[-0.07435228, -0.00649564, 0.00714862, -0.07459251],\n",
" [-0.07593353, -0.01067953, 0.01189881, -0.07602356],\n",
" [-0.0795701 , 0.00016089, 0.00015974, -0.07975458]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.07602356, -0.0795701 , -0.07435228])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.07602356083253432 -0.07975458355080561\n",
"\n",
"================================================================================================================================================================\n",
"Runtime information: \n",
"Total runtime: 2369.9501148750005 s\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Initial setup: 0.17246508400057792 s\n",
"Hamiltonian conversion and XC field extraction: 0.904 s\n",
"Pair and site datastructure creatrions: 0.028 s\n",
"k set cration and distribution: 0.009 s\n",
"Rotating XC potential: 0.212 s\n",
"Greens function inversion: 2368.445 s\n",
"Calculate energies and magnetic components: 0.180 s"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"Magnetic entities integrated.\n",
"Pairs integrated.\n",
"Magnetic parameters calculated.\n",
"##################################################################### GROGU OUTPUT #############################################################################\n",
"================================================================================================================================================================\n",
"Input file: \n",
"Not yet specified.\n",
"Output file: \n",
"test.pickle\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: 1000\n",
"k point directions: x\n",
"Ebot: -13\n",
"Eset: 100\n",
"Esetp: 10000\n",
"================================================================================================================================================================\n",
"Atomic information: \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",
"[5]Fe(2) 1.8954667088117545 1.0943913231921656 10.285002698393109\n",
"\n",
"================================================================================================================================================================\n",
"Exchange [meV]\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Magnetic entity1 Magnetic entity2 [i j k] d [Ang]\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"[3]Fe(2) [4]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -29.36445059067917\n",
"DMI: [-8.91784089e-01 4.16579134e+00 -3.59901207e-06]\n",
"Symmetric-anisotropy: [-4.5685362 1.18301202 0.3323615 4.15576769 -0.02915828]\n",
"Energies for debugging: \n",
"array([[-0.02818144, -0.00086263, 0.00092094, -0.02597893],\n",
" [-0.03690122, -0.00415577, 0.00417581, -0.03393299],\n",
" [-0.03160383, -0.00033237, -0.00033236, -0.03122058]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.03393299, -0.03160383, -0.02818144])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.03393298678958339 -0.031220578858898954\n",
"\n",
"[4]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -74.99083333178646\n",
"DMI: [-6.82222254e+00 1.12891474e+01 5.95795483e-04]\n",
"Symmetric-anisotropy: [-1.03214008 0.64046775 -0.16031723 10.67951504 -0.32709547]\n",
"Energies for debugging: \n",
"array([[-0.07435037, -0.00649513, 0.00714932, -0.07459916],\n",
" [-0.07593293, -0.01067952, 0.01189878, -0.07602297],\n",
" [-0.0795695 , 0.00016091, 0.00015972, -0.07975399]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.07602297, -0.0795695 , -0.07435037])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.07602297341278706 -0.07975398680085131\n",
"\n",
"================================================================================================================================================================\n",
"Runtime information: \n",
"Total runtime: 501.2763705420002 s\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Initial setup: 0.12531970899999578 s\n",
"Hamiltonian conversion and XC field extraction: 0.888 s\n",
"Pair and site datastructure creatrions: 0.062 s\n",
"k set cration and distribution: 0.024 s\n",
"Rotating XC potential: 0.245 s\n",
"Greens function inversion: 499.875 s\n",
"Calculate energies and magnetic components: 0.056 s"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"Magnetic entities integrated.\n",
"Pairs integrated.\n",
"Magnetic parameters calculated.\n",
"##################################################################### GROGU OUTPUT #############################################################################\n",
"================================================================================================================================================================\n",
"Input file: \n",
"Not yet specified.\n",
"Output file: \n",
"test.pickle\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: 1000\n",
"k point directions: x\n",
"Ebot: -13\n",
"Eset: 50\n",
"Esetp: 10000\n",
"================================================================================================================================================================\n",
"Atomic information: \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",
"[5]Fe(2) 1.8954667088117545 1.0943913231921656 10.285002698393109\n",
"\n",
"================================================================================================================================================================\n",
"Exchange [meV]\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Magnetic entity1 Magnetic entity2 [i j k] d [Ang]\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"[3]Fe(2) [4]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -29.359042196157713\n",
"DMI: [-8.93813379e-01 4.16574774e+00 -3.59758916e-06]\n",
"Symmetric-anisotropy: [-4.57147523 1.18770348 0.33237646 4.15572505 -0.02850636]\n",
"Energies for debugging: \n",
"array([[-0.02817134, -0.00086531, 0.00092232, -0.02597527],\n",
" [-0.03689861, -0.00415573, 0.00417577, -0.03393052],\n",
" [-0.03160252, -0.00033238, -0.00033237, -0.03121927]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.03393052, -0.03160252, -0.02817134])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.03393051742760642 -0.031219274782673133\n",
"\n",
"[4]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -74.98947011735306\n",
"DMI: [-6.82126017e+00 1.12890803e+01 6.51055432e-04]\n",
"Symmetric-anisotropy: [-1.03170277 0.64211178 -0.16032097 10.67947464 -0.32714492]\n",
"Energies for debugging: \n",
"array([[-0.07434736, -0.00649412, 0.00714841, -0.07459988],\n",
" [-0.07593108, -0.01067947, 0.01189869, -0.07602117],\n",
" [-0.07956766, 0.00016097, 0.00015967, -0.07975216]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.07602117, -0.07956766, -0.07434736])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.07602117288317486 -0.07975215648287909\n",
"\n",
"[3]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -75.04936117965987\n",
"DMI: [ 6.85701206e+00 -1.13051160e+01 6.39732773e-04]\n",
"Symmetric-anisotropy: [ -1.05318057 0.68103129 -0.16032 -10.70364791 0.31624841]\n",
"Energies for debugging: \n",
"array([[-0.07436833, 0.00654076, -0.00717326, -0.07467721],\n",
" [-0.07601601, 0.01070365, -0.01190658, -0.07610254],\n",
" [-0.07956796, 0.00016096, 0.00015968, -0.07975245]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.07610254, -0.07956796, -0.07436833])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.07610254174486064 -0.07975245001820208\n",
"\n",
"================================================================================================================================================================\n",
"Runtime information: \n",
"Total runtime: 261.890385875 s\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Initial setup: 0.1216103749998183 s\n",
"Hamiltonian conversion and XC field extraction: 1.045 s\n",
"Pair and site datastructure creatrions: 0.011 s\n",
"k set cration and distribution: 0.006 s\n",
"Rotating XC potential: 0.224 s\n",
"Greens function inversion: 260.450 s\n",
"Calculate energies and magnetic components: 0.033 s"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"Magnetic entities integrated.\n",
"Pairs integrated.\n",
"Magnetic parameters calculated.\n",
"##################################################################### GROGU OUTPUT #############################################################################\n",
"================================================================================================================================================================\n",
"Input file: \n",
"Not yet specified.\n",
"Output file: \n",
"test.pickle\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: 3000\n",
"k point directions: x\n",
"Ebot: -13\n",
"Eset: 50\n",
"Esetp: 10000\n",
"================================================================================================================================================================\n",
"Atomic information: \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",
"[5]Fe(2) 1.8954667088117545 1.0943913231921656 10.285002698393109\n",
"\n",
"================================================================================================================================================================\n",
"Exchange [meV]\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Magnetic entity1 Magnetic entity2 [i j k] d [Ang]\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"[3]Fe(2) [4]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -29.437239263932703\n",
"DMI: [-8.91096571e-01 4.16386016e+00 -3.68076797e-06]\n",
"Symmetric-anisotropy: [-4.68204753 1.22980695 0.33079224 4.15200152 -0.03180686]\n",
"Energies for debugging: \n",
"array([[-0.02820743, -0.00085929, 0.0009229 , -0.025985 ],\n",
" [-0.03707978, -0.004152 , 0.00417572, -0.03411929],\n",
" [-0.03167345, -0.0003308 , -0.00033079, -0.03129215]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.03411929, -0.03167345, -0.02820743])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.03411928679618277 -0.03129215156448653\n",
"\n",
"[4]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -74.98086703681366\n",
"DMI: [-6.81158085e+00 1.13265508e+01 5.68487147e-04]\n",
"Symmetric-anisotropy: [-1.03582622 0.63764125 -0.15959679 10.72110977 -0.32326625]\n",
"Energies for debugging: \n",
"array([[-0.07434323, -0.00648831, 0.00713485, -0.07458268],\n",
" [-0.07592777, -0.01072111, 0.01193199, -0.07601669],\n",
" [-0.07959289, 0.00016017, 0.00015903, -0.07977653]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.07601669, -0.07959289, -0.07434323])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.07601669325607426 -0.07977653333847946\n",
"\n",
"================================================================================================================================================================\n",
"Runtime information: \n",
"Total runtime: 808.34989875 s\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Initial setup: 0.1594020829999998 s\n",
"Hamiltonian conversion and XC field extraction: 0.995 s\n",
"Pair and site datastructure creatrions: 0.021 s\n",
"k set cration and distribution: 0.021 s\n",
"Rotating XC potential: 0.232 s\n",
"Greens function inversion: 806.884 s\n",
"Calculate energies and magnetic components: 0.037 s"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"Magnetic entities integrated.\n",
"Pairs integrated.\n",
"Magnetic parameters calculated.\n",
"##################################################################### GROGU OUTPUT #############################################################################\n",
"================================================================================================================================================================\n",
"Input file: \n",
"Not yet specified.\n",
"Output file: \n",
"test.pickle\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: 2000\n",
"k point directions: x\n",
"Ebot: -13\n",
"Eset: 50\n",
"Esetp: 10000\n",
"================================================================================================================================================================\n",
"Atomic information: \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",
"[5]Fe(2) 1.8954667088117545 1.0943913231921656 10.285002698393109\n",
"\n",
"================================================================================================================================================================\n",
"Exchange [meV]\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Magnetic entity1 Magnetic entity2 [i j k] d [Ang]\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"[3]Fe(2) [4]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -29.313466612791018\n",
"DMI: [-8.85106686e-01 4.16843252e+00 -3.50602459e-06]\n",
"Symmetric-anisotropy: [-4.7530289 1.27630805 0.34469396 4.15760367 -0.03005581]\n",
"Energies for debugging: \n",
"array([[-0.02803716, -0.00085505, 0.00091516, -0.02583675],\n",
" [-0.03702114, -0.0041576 , 0.00417926, -0.0340665 ],\n",
" [-0.03178903, -0.0003447 , -0.00034469, -0.03139132]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.0340665 , -0.03178903, -0.02803716])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.03406649551707385 -0.031391324924774755\n",
"\n",
"[4]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -74.97640334609967\n",
"DMI: [-6.84069197e+00 1.12809486e+01 8.75098125e-04]\n",
"Symmetric-anisotropy: [-1.03102826 0.63307456 -0.16123618 10.67285449 -0.32447457]\n",
"Energies for debugging: \n",
"array([[-0.07434333, -0.00651622, 0.00716517, -0.07457845],\n",
" [-0.07591065, -0.01067285, 0.01188904, -0.07600743],\n",
" [-0.07959609, 0.00016211, 0.00016036, -0.07978159]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.07600743, -0.07959609, -0.07434333])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.07600743160328334 -0.07978158717380361\n",
"\n",
"================================================================================================================================================================\n",
"Runtime information: \n",
"Total runtime: 534.5927785000001 s\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Initial setup: 0.14542879100008577 s\n",
"Hamiltonian conversion and XC field extraction: 0.903 s\n",
"Pair and site datastructure creatrions: 0.016 s\n",
"k set cration and distribution: 0.014 s\n",
"Rotating XC potential: 0.242 s\n",
"Greens function inversion: 533.223 s\n",
"Calculate energies and magnetic components: 0.050 s"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"Magnetic entities integrated.\n",
"Pairs integrated.\n",
"Magnetic parameters calculated.\n",
"##################################################################### GROGU OUTPUT #############################################################################\n",
"================================================================================================================================================================\n",
"Input file: \n",
"Not yet specified.\n",
"Output file: \n",
"test.pickle\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: 1000\n",
"k point directions: x\n",
"Ebot: -13\n",
"Eset: 50\n",
"Esetp: 10000\n",
"================================================================================================================================================================\n",
"Atomic information: \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",
"[5]Fe(2) 1.8954667088117545 1.0943913231921656 10.285002698393109\n",
"\n",
"================================================================================================================================================================\n",
"Exchange [meV]\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Magnetic entity1 Magnetic entity2 [i j k] d [Ang]\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"[3]Fe(2) [4]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -29.359042196157713\n",
"DMI: [-8.93813379e-01 4.16574774e+00 -3.59758916e-06]\n",
"Symmetric-anisotropy: [-4.57147523 1.18770348 0.33237646 4.15572505 -0.02850636]\n",
"Energies for debugging: \n",
"array([[-0.02817134, -0.00086531, 0.00092232, -0.02597527],\n",
" [-0.03689861, -0.00415573, 0.00417577, -0.03393052],\n",
" [-0.03160252, -0.00033238, -0.00033237, -0.03121927]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.03393052, -0.03160252, -0.02817134])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.03393051742760642 -0.031219274782673133\n",
"\n",
"[4]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -74.98947011735306\n",
"DMI: [-6.82126017e+00 1.12890803e+01 6.51055432e-04]\n",
"Symmetric-anisotropy: [-1.03170277 0.64211178 -0.16032097 10.67947464 -0.32714492]\n",
"Energies for debugging: \n",
"array([[-0.07434736, -0.00649412, 0.00714841, -0.07459988],\n",
" [-0.07593108, -0.01067947, 0.01189869, -0.07602117],\n",
" [-0.07956766, 0.00016097, 0.00015967, -0.07975216]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.07602117, -0.07956766, -0.07434736])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.07602117288317486 -0.07975215648287909\n",
"\n",
"================================================================================================================================================================\n",
"Runtime information: \n",
"Total runtime: 272.760337666 s\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Initial setup: 0.13726933300000033 s\n",
"Hamiltonian conversion and XC field extraction: 0.970 s\n",
"Pair and site datastructure creatrions: 0.013 s\n",
"k set cration and distribution: 0.017 s\n",
"Rotating XC potential: 0.217 s\n",
"Greens function inversion: 271.373 s\n",
"Calculate energies and magnetic components: 0.034 s"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"##################################################################### GROGU OUTPUT #############################################################################\n",
"================================================================================================================================================================\n",
"Input file: \n",
"Not yet specified.\n",
"Output file: \n",
"test.pickle\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: 20\n",
"k point directions: xy\n",
"Ebot: -15\n",
"Eset: 100\n",
"Esetp: 1000\n",
"================================================================================================================================================================\n",
"Atomic information: \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",
"[5]Fe(2) 1.8954667088117545 1.0943913231921656 10.285002698393109\n",
"\n",
"================================================================================================================================================================\n",
"Exchange [meV]\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Magnetic entity1 Magnetic entity2 [i j k] d [Ang]\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"[3]Fe(2) [4]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -63.51146633376187\n",
"DMI: [-9.32966002e-01 -8.89800564e-04 -2.02862680e-06]\n",
"Symmetric-anisotropy: [-3.33721683e+00 1.29672703e+00 6.03687008e-04 -8.09695186e-04\n",
" -6.00702131e-06]\n",
"Energies for debugging: \n",
"array([[-6.22147393e-02, -9.32959995e-04, 9.32972009e-04,\n",
" -6.14709765e-02],\n",
" [-6.75755955e-02, 8.09695186e-07, -9.69905943e-07,\n",
" -6.68486832e-02],\n",
" [-5.36624493e-02, -6.05715634e-07, -6.01658381e-07,\n",
" -5.36632256e-02]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.06684868, -0.05366245, -0.06221474])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.06684868316242108 -0.053663225627735525\n",
"\n",
"[4]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -60.97776480541008\n",
"DMI: [-3.79946502e+00 6.15248803e+00 3.52644581e-03]\n",
"Symmetric-anisotropy: [0.09164177 0.11139034 0.07107136 6.24190815 0.03636875]\n",
"Energies for debugging: \n",
"array([[-6.08663745e-02, -3.83583377e-03, 3.76309627e-03,\n",
" -6.11807969e-02],\n",
" [-6.06347990e-02, -6.24190815e-03, 6.06306790e-03,\n",
" -6.08861230e-02],\n",
" [-5.98977554e-02, -6.75449092e-05, -7.45978008e-05,\n",
" -5.98154343e-02]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.06088612, -0.05989776, -0.06086637])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.060886123032493376 -0.05981543426889457\n",
"\n",
"[3]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n",
"Isotropic: -60.9832595805832\n",
"DMI: [ 3.78505688e+00 -6.13845903e+00 3.52688440e-03]\n",
"Symmetric-anisotropy: [ 0.08120798 0.11261323 0.07107254 -6.23157554 -0.04249828]\n",
"Energies for debugging: \n",
"array([[-6.08706463e-02, 3.82755516e-03, -3.74255860e-03,\n",
" -6.11770808e-02],\n",
" [-6.06458648e-02, 6.23157554e-03, -6.04534251e-03,\n",
" -6.09020516e-02],\n",
" [-5.98977110e-02, -6.75456589e-05, -7.45994277e-05,\n",
" -5.98153841e-02]])\n",
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
"array([-0.06090205, -0.05989771, -0.06087065])\n",
"Test J_xx = E(y,z) = E(z,y)\n",
"-0.060902051599344476 -0.05981538413498572\n",
"\n",
"================================================================================================================================================================\n",
"Runtime information: \n",
"Total runtime: 203.083672208 s\n",
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Initial setup: 0.1483790410000001 s\n",
"Hamiltonian conversion and XC field extraction: 0.861 s\n",
"Pair and site datastructure creatrions: 0.017 s\n",
"k set cration and distribution: 0.022 s\n",
"Rotating XC potential: 0.247 s\n",
"Greens function inversion: 201.715 s\n",
"Calculate energies and magnetic components: 0.073 s"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"========================================\n",
" \n",
"Atom Angstrom\n",
"# Label, x y z Sx Sy Sz #Q Lx Ly Lz Jx Jy Jz\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Te1 1.8955 1.0943 13.1698 -0.0000 0.0000 -0.1543 # 5.9345 -0.0000 0.0000 -0.0537 -0.0000 0.0000 -0.2080 \n",
"Te2 1.8955 1.0943 7.4002 0.0000 -0.0000 -0.1543 # 5.9345 0.0000 -0.0000 -0.0537 0.0000 -0.0000 -0.2080 \n",
"Ge3 -0.0000 2.1887 10.2850 0.0000 0.0000 -0.1605 # 3.1927 -0.0000 0.0000 0.0012 0.0000 0.0000 -0.1593 \n",
"Fe4 -0.0000 0.0000 11.6576 0.0001 -0.0001 2.0466 # 8.3044 0.0000 -0.0000 0.1606 0.0001 -0.0001 2.2072 \n",
"Fe5 -0.0000 0.0000 8.9124 -0.0001 0.0001 2.0466 # 8.3044 -0.0000 0.0000 0.1606 -0.0001 0.0001 2.2072 \n",
"Fe6 1.8955 1.0944 10.2850 0.0000 0.0000 1.5824 # 8.3296 -0.0000 -0.0000 0.0520 -0.0000 0.0000 1.6344 \n",
"==================================================================================================================================\n",
" \n",
"Exchange meV\n",
"--------------------------------------------------------------------------------\n",
"# at1 at2 i j k # d (Ang)\n",
"--------------------------------------------------------------------------------\n",
"Fe4 Fe5 0 0 0 # 2.7452\n",
"Isotropic -82.0854\n",
"DMI 0.12557 -0.00082199 6.9668e-08\n",
"Symmetric-anisotropy -0.60237 -0.83842 -0.00032278 -1.2166e-05 -3.3923e-05\n",
"--------------------------------------------------------------------------------\n",
"Fe4 Fe6 0 0 0 # 2.5835\n",
"Isotropic -41.9627\n",
"DMI 1.1205 -1.9532 0.0018386\n",
"Symmetric-anisotropy 0.26007 -0.00013243 0.12977 -0.069979 -0.042066\n",
"--------------------------------------------------------------------------------\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": ".venv",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}