{ "cells": [ { "cell_type": "markdown", "id": "38bf69d0", "metadata": {}, "source": [ "# Spontaneous emission of TLS: TCP socket\n", "\n", "Here, we introduce the two-level system (TLS) spontaneous emission tutorial using `maxwelllink.Molecule` and a TCP socket." ] }, { "cell_type": "markdown", "id": "f71d059c", "metadata": {}, "source": [ "## 1. Setting up the socket communication layer\n", "\n", "Using the TCP socket requires setting the `hostname` and `port number`. On a local machine, we can use the helper function `get_available_host_port()` from **MaxwellLink** to obtain these two pieces of information. Then, we initialize a `SocketHub` instance to provide the socket communication in **MaxwellLink**." ] }, { "cell_type": "code", "execution_count": 1, "id": "2c8428dc", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SocketHub listening on 127.0.0.1:59062\n" ] } ], "source": [ "import numpy as np\n", "import maxwelllink as mxl\n", "from maxwelllink import sockets as mxs\n", "\n", "try:\n", " import meep as mp\n", "except ImportError as exc: \n", " raise RuntimeError(\n", " \"Meep is required for this tutorial.\"\n", " \"Install via conda: conda install -c conda-forge pymeep=*=mpi_mpich_*\"\n", " ) from exc\n", "\n", "host, port = mxs.get_available_host_port()\n", "hub = mxl.SocketHub(host=host, port=port, timeout=10.0, latency=1e-5)\n", "\n", "print(f\"SocketHub listening on {host}:{port}\")\n" ] }, { "cell_type": "markdown", "id": "e0f5c662", "metadata": {}, "source": [ "Before the simulation, users can understand the units system in MEEP with the following built-in function:" ] }, { "cell_type": "code", "execution_count": 2, "id": "d3400074", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Now let's take a moment to understand the Meep unit system based on your provided units:\n", "Given that you have specified:\n", "- Time unit: 0.1 fs\n", "By further assuming the MEEP simulation resolution is 10 pixels per length unit, we can derive the following Meep units:\n", "\n", "\n", " ######### MaxwellLink Units Helper #########\n", " MEEP uses its own units system, which is based on the speed of light in vacuum (c=1), \n", " the permittivity of free space (epsilon_0=1), and the permeability of free space (mu_0=1). \n", " To couple MEEP with molecular dynamics, we set [c] = [epsilon_0] = [mu_0] = [hbar] = 1. \n", " By further defining the time unit or length unit, we can fix the units system of MEEP (mu).\n", "\n", " - Time [t]: 1 mu = 1.0000E-01 fs = 4.1341E+00 a.u.\n", " - Length [x]: 1 mu = 2.9979E+01 nm\n", " - EM wavelength of 1 mu, angular frequency omega = 2pi mu = 4.1375E+01 eV = 3.3371E+05 cm-1 = 1.5205E+00 a.u.\n", " - Note that sources and dielectrics defined in MEEP use rotational frequency (f=omega/2pi), \n", " - so probabably we need covert 1 eV photon energy to rotational frequency f = 2.4169E-02 mu\n", " - Electric field [E]: 1 mu = 6.6486E+07 V/m = 1.2930E-04 a.u.\n", "\n", " Given the simulation resolution = 10,\n", " - FDTD dt = 5.0000E-02 mu (0.5/resolution) = 5.0000E-03 fs\n", " - FDTD dx = 1.0000E-01 mu (1.0/resolution) = 2.9979E+00 nm\n", " Hope this helps!\n", " ############################################\n", "\n", "\n" ] } ], "source": [ "# or use: mxl.meep_units_helper(length_units_nm=2.9979E+01)\n", "mxl.meep_units_helper(time_units_fs=0.1)" ] }, { "cell_type": "markdown", "id": "3f43f478", "metadata": {}, "source": [ "## 2. Bind Molecule and EM solver to the SocketHub\n", "\n", "Then, we create a `Molecule` instance to define the information of this molecule in the EM simulation environment, including the `center`, `size`, `sigma` (width of the molecular polarization distribution), and `dimensions`. \n", "\n", "We also need to setup the EM solver (**MEEP**) using `mxl.MeepSimulation`. This class is a wrapper of the `meep.Simulation` object with extended parameters for **MaxwellLink**." ] }, { "cell_type": "code", "execution_count": 3, "id": "a5f84aa1", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[Init Molecule] Under socket mode, registered molecule with ID 0\n", "\n", "\n", " ######### MaxwellLink Units Helper #########\n", " MEEP uses its own units system, which is based on the speed of light in vacuum (c=1), \n", " the permittivity of free space (epsilon_0=1), and the permeability of free space (mu_0=1). \n", " To couple MEEP with molecular dynamics, we set [c] = [epsilon_0] = [mu_0] = [hbar] = 1. \n", " By further defining the time unit or length unit, we can fix the units system of MEEP (mu).\n", "\n", " - Time [t]: 1 mu = 1.0000E-01 fs = 4.1341E+00 a.u.\n", " - Length [x]: 1 mu = 2.9979E+01 nm\n", " - EM wavelength of 1 mu, angular frequency omega = 2pi mu = 4.1375E+01 eV = 3.3371E+05 cm-1 = 1.5205E+00 a.u.\n", " - Note that sources and dielectrics defined in MEEP use rotational frequency (f=omega/2pi), \n", " - so probabably we need covert 1 eV photon energy to rotational frequency f = 2.4169E-02 mu\n", " - Electric field [E]: 1 mu = 6.6486E+07 V/m = 1.2930E-04 a.u.\n", "\n", " Given the simulation resolution = 10,\n", " - FDTD dt = 5.0000E-02 mu (0.5/resolution) = 5.0000E-03 fs\n", " - FDTD dx = 1.0000E-01 mu (1.0/resolution) = 2.9979E+00 nm\n", " Hope this helps!\n", " ############################################\n", "\n", "\n" ] } ], "source": [ "molecule = mxl.Molecule(\n", " hub=hub,\n", " center=mp.Vector3(0, 0, 0),\n", " size=mp.Vector3(1, 1, 1),\n", " sigma=0.1,\n", " dimensions=2,\n", ")\n", "\n", "sim = mxl.MeepSimulation(\n", " hub=hub,\n", " molecules=[molecule],\n", " cell_size=mp.Vector3(8, 8, 0),\n", " boundary_layers=[mp.PML(3.0)],\n", " resolution=10,\n", " # fix a units system \n", " time_units_fs=0.1,\n", ")\n" ] }, { "cell_type": "markdown", "id": "20cf3a66", "metadata": {}, "source": [ "## 3. Python way to launch `mxl_driver` on a separate terminal\n", "\n", "Generally, using the Socket Interface requires launching the EM simulation in one terminal and then starting the molecular driver simulation in a separate terminal. To avoid opening a second terminal, below we introduce a Python helper function `launch_tls_driver(...)`, which will launch `mxl_driver` from Python (so we can stay within this notebook to finish this tutorial).\n", "\n", "Here, we set the TLS starting at the initial excited-state population of 1e-4.\n", "\n", "Immediately after launching this driver in the background, we run the simulation using `sim.run(...)`. This function is a wrapper of the `meep.Simulation.run(...)` function, which can accept user-defined step functions." ] }, { "cell_type": "code", "execution_count": 4, "id": "d2c4f2ad", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Launching TLS driver via subprocess...\n", "If you prefer to run it manually, execute:\n", " /Users/taoli/miniforge3/envs/mxl/bin/mxl_driver --model tls --address 127.0.0.1 --port 59062 --param \"omega=0.242, mu12=187, orientation=2, pe_initial=1e-4\"\n", "-----------\n", "Initializing structure...\n", "time for choose_chunkdivision = 0.000616074 s\n", "Working in 2D dimensions.\n", "Computational cell is 8 x 8 x 0 with resolution 10\n", "time for set_epsilon = 0.00299597 s\n", "-----------\n", "[SocketHub][initialization] Time step in atomic units: CONNECTED: mol 0 <- 127.0.0.1:59065\n", "0.20670686667500004\n", "[initialization] Assigned a molecular ID: 0\n", "init TLSModel with dt = 0.206707 a.u., molecule ID = 0\n", "[initialization] Finished initialization for molecular ID: 0\n", "Meep progress: 379.15000000000003/400.0 = 94.8% done in 4.0s, 0.2s to go\n", "on time step 7592 (time=379.6), 0.000526919 s/step\n", "run 0 finished at t = 400.0 (8000 timesteps)\n", "Received STOP, exiting\n", "[SocketHub] DISCONNECTED: mol 0 from 127.0.0.1:59065\n" ] } ], "source": [ "import shlex\n", "import shutil\n", "import subprocess\n", "import time\n", "\n", "\n", "def launch_tls_driver(host: str, port: int, sleep_time: float = 0.5):\n", " executable = shutil.which('mxl_driver')\n", " if executable is None:\n", " raise RuntimeError('mxl_driver executable not found in PATH.')\n", " cmd = (\n", " f\"{executable} --model tls --address {host} --port {port}\"\n", " f' --param \"omega=0.242, mu12=187, orientation=2, pe_initial=1e-4\"'\n", " )\n", " print('Launching TLS driver via subprocess...')\n", " print('If you prefer to run it manually, execute:')\n", " print(' ' + cmd)\n", " argv = shlex.split(cmd)\n", " proc = subprocess.Popen(argv)\n", " time.sleep(sleep_time)\n", " return proc\n", "\n", "launch_tls_driver(host, port)\n", "\n", "sim.run(until=400)" ] }, { "cell_type": "markdown", "id": "cdabe7d1", "metadata": {}, "source": [ "## 4. Retrieve molecular simulation data\n", "\n", "After the simulation, we can retrieve molecular simulation data from `molecule.extra`, a Python dictionary which stores the molecular information sent from the driver code at each step of the simulation." ] }, { "cell_type": "code", "execution_count": 5, "id": "9120e20d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Collected 8001 samples.\n" ] } ], "source": [ "# users can also use molecule.additional_data_history to access the time-resolved data recorded during the simulation, \n", "# population = np.array([entry[\"Pe\"] for entry in molecule.additional_data_history])\n", "# time_au = np.array([entry[\"time_au\"] for entry in molecule.additional_data_history])\n", "# but here we demonstrate the use of molecule.extra which is more convenient for post-processing and plotting. \n", "\n", "population = molecule.extra[\"Pe\"]\n", "time_au = molecule.extra[\"time_au\"]\n", "\n", "print(f\"Collected {population.size} samples.\")" ] }, { "cell_type": "markdown", "id": "ae529503", "metadata": {}, "source": [ "## 5. Compare with the Analytical Golden-Rule Decay\n", "\n", "Finally, we can compare this numerical simulation with analytical golden-rule rate calculations, with the 2D spontaneus emission in vacuum as:\n", "\n", "$$\n", "k_{\\rm {FGR}} = \\frac{|\\mu_{12}|^2 \\omega^2}{2\\hbar\\epsilon_0 c^2}\n", "$$\n", "\n", "The corresponding TLS excited-state population decay dynamics obey:\n", "\n", "$$\n", "P_{\\rm e}^{\\rm {QM}}(t) = P_{\\rm e}(0) e^{-k_{\\rm{FGR}} t}\n", "$$\n", "\n", "When the EM field is described entirely classically, more than half a century ago, Jaynes and collaborators (https://ieeexplore.ieee.org/document/1443594) calculated the semiclassical spontaneous emission rate. More recently, we have also reproduced this semiclassical excited-state population decay (https://doi.org/10.1103/PhysRevA.97.032105):\n", "\n", "$$\n", "P_{\\rm e}^{\\rm {sc}}(t) = \\frac{e^{-k_{\\rm{FGR}} t}}{e^{-k_{\\rm{FGR}} t} + \\frac{1 - P_{\\rm e}(0) }{P_{\\rm e}(0) }}\n", "$$\n", "\n", "When $P_{\\rm e}(0)\\rightarrow 0$, the semiclassical decay dynamics exactly agree with the quantum correspondance $P_{\\rm e}^{\\rm {QM}}(t)$.\n", "\n", "As shown below, using $P_{\\rm e}(0)= 10^{-4}$, our semiclassical simulation exactly reproduces the quantum golden-rule decay." ] }, { "cell_type": "code", "execution_count": 6, "id": "ebe9c2ed", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "std_dev=2.012e-03, max_abs_diff=7.928e-03\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAGGCAYAAACNCg6xAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAapNJREFUeJzt3Qd4FEUbB/B/7nLpIaQQemih916k9yZNBBERRBSwgQ0QUcBCUwHp0hSwgFQRBAGFEHrvTSB0QiAhCZCe2+95J1/OJCR4BwmXS/6/51nubm9ud3Y35N7MzL5jp2maBiIiIiL6T7r/LkJEREREDJyIiIiILMAWJyIiIiIzMXAiIiIiMhMDJyIiIiIzMXAiIiIiMhMDJyIiIiIzMXAiIiIiMpO9uQUpaxiNRty4cQPu7u6ws7PjaSYiInrKJBf4vXv3UKhQIeh0j25TYuBkZRI0FS1a1NrVICIiyvWuXr2KIkWKPPI8MHCyMmlpSr5YefLksXZ1iIiIcp3IyEjViJH8nfwoDJysLLl7ToImBk5ERETWY86QGQ4OJyIiIjITAyciIiIiMzFwIiIiIjITxzgREVGul5iYiPj4+Fx/HnIqg8EAvV6fKdti4ERERLk6f09wcDDCw8OtXRXKYnnz5kWBAgWeOGciAyciIsq1koMmX19fuLi4MBFxDg2Oo6KiEBISol4XLFjwibbHwImIiHJt91xy0OTt7W3t6lAWcnZ2Vo8SPMn1fpJuO6sPDp81axZKlCgBJycn1KxZE4GBgY8sHxAQoMpJ+ZIlS2LOnDkPlVm5ciUqVKgAR0dH9bh69WqL97tq1Sq0adMGPj4+6i+QI0eOPLSN2NhYvP3226qMq6srOnXqhGvXrj3WeSAioqcreUyTtDRRzufy/+v8pGPZrBo4LVu2DEOHDsXHH3+Mw4cPo1GjRmjXrh2uXLmSbvmgoCC0b99elZPyI0eOxDvvvKMCpWS7d+9Gz5490adPHxw9elQ99ujRA3v37rVovw8ePMAzzzyDCRMmZFh/2YYEZUuXLsWOHTtw//59dOzYUf0VQ0REtoHzhOYOdpk1H6xmRXXq1NEGDRqUal25cuW0ESNGpFt+2LBh6v2UBg4cqNWrV8/0ukePHlrbtm1TlWnTpo32wgsvPNZ+g4KCNDlNhw8fTrU+PDxcMxgM2tKlS03rrl+/rul0Om3jxo2auSIiItT25ZGIiJ6e6Oho7dSpU+qRcvf1jrDgu9hqLU5xcXE4ePAgWrdunWq9vN61a1e6n5HWpLTlpTvtwIEDpqa3jMokb/Nx9pse2YbsM+V2ZFblSpUqWbSdrPLX6VuIiOattUREubV1Zc2aNVm+n+LFi2Pq1KnZZjtPg9UCpzt37qgurfz586daL6/lLof0yPr0yickJKjtPapM8jYfZ78Z1cXBwQGenp4WbUfGRclkgimXzHY1LAoDF+3B9Ekf4ULw3UzfPhERWZcMch44cCD8/PzUeF65zV4aCaTxQNy8eVMNQclufvjhB5UWIK39+/fj9ddfhy2wz259jnLb4KP6IdMrn3a9Odu0dL/m+q/tjB8/HmPHjkVWComMwc9Ok1BHO4GNa3QoNWhSlu6PiIierueee071eixatEjdKHXr1i389ddfCAsLU+9LIGVL8uXLB1thtRYnuRNNbgdM2zojUXTa1qBk8oOQXnl7e3vTraQZlUne5uPsN6O6SLff3bt3LdrORx99hIiICNNy9epVZLaaxb3g2aCvet4seCHibp3J9H0QEZF1SAoFuSFp4sSJaNasGYoVK4Y6deqo75cOHTo81FV36dIl9frXX39VN0PJrfm1a9fGuXPnVEtPrVq14ObmhrZt2+L27dum/TRt2lTdBJVSly5d0K9fvwzrNnnyZFSuXFndaV60aFG88cYb6sYpsW3bNrzyyivqu0/qI8uYMWPS7aqTm7U6d+6s6pUnTx51k5cEh8nkc9WqVcOSJUvUZz08PPDCCy/g3r17yLGBk3RzSRqAzZs3p1ovrxs0aJDuZ+rXr/9Q+U2bNqmLLunUH1UmeZuPs9/0yDZknym3I02jJ06ceOR2pElVfghSLlmhZIsB2GVXHY6Ix4PlgxERFZsl+yEiynHJEuMSnvqS3HtiDgkmZJHASIZ/mGv06NEYNWoUDh06pBocevXqhWHDhuHbb79VKXkuXLiATz/9FE9Cp9Nh2rRp6rtQWsP+/vtvtQ8h340SHMn3nnxfyvLBBx88tA05FxKgSeuZpCCS71mpm9wxn5Ksk3Owbt06tUjZR90JnyO66t577z2VLkACHwl45s6dq6LMQYMGqfcler5+/ToWL16sXsv6GTNmqM+99tprqi93wYIF+OWXX0zbHDJkCBo3bqwicYlWf/vtN2zZskVF5+buV8gFk3U3btxQr8+ePWtqaZJFottXX30V77//vmrt8vLyUj8AEmm3bNkS1qbX67C34ieocvxFeN45hDHjhqFy1w/xXM0i1q4aEVG2FR2fiAqf/vnU93vqszZwcTDvK1mCHhkrJN+DksuwRo0aaNKkiWpxqVKlSoafk+8oGQeV/F0pgZN070nqHSHfabLdJzE0RQuV5Er8/PPPMXjwYJU7URou5LtTWpoe1ZUo39nHjh1TKYik1UpIy1LFihVVC5m0lgmj0ajq6+7url7L97ocz5dffomsZNU8ThI9SvT52WefqSa37du3448//lDNjkKi0ZS5leQiyPvS3Cfl5YJIZCt9vckkopW8St9//736AZKTKnmb6tata/Z+xdq1a1G9enVTs6f8QMrrlAk3p0yZoqJiaUKUHzxJrvX7779n2kSCT6pZ3ZqYkNBLPf9Q/wt+/DMQiUbz/6ohIqLsSb735A97+a6SYEi+FyWAelTgkzKoSh5SIn/sp1yXPC3J49q6dStatWqFwoULq4Dm5ZdfRmhoqMqNaK7Tp0+rgCk5aBKSzFoGlct7yaSLLjloSp5K5UnrbxODw6X/U5b0pPcDIFG1NDM+Svfu3dXyuPsV0of7qH5cIVnHp0+frpbsqFrRvLjQ9T0EbT2EEg+O4o3ouTh9sz0qFfawdtWIiLIlZ4Netf5YY7+Wku8gCVJkkS62AQMGqO64jL67koe0iOSbmNKuk1aclN1uabsQ4x+Rdfvy5csqSbX03kjDhvTESG+PtGRZkq07o5us0q5PWff06p9VrD7lCmWt52r5ocQrC3HSqSbGJfTGzvNJaRuIiOhh8uUrXWZPe8mMu7qlVcaSlh1z7nSTnp9kksrnxIkTGZaXnIqSHuibb75BvXr1UKZMGdNwl2TSXfdfs2vIcUhvU8qbp06dOqUGlZcvXx7WxsApN/Dxx+5n5iNIK4idF0KtXRsiInoC0vXVvHlz/Pjjj6axQMuXL8ekSZPU2N7MIvtYv369Ws6cOaN6acLDwzMsX6pUKRU4SS/MxYsX1biktPPJSvea3GUnY5Ekr2JUVNRD25FxwtKt2Lt3b9XDtG/fPtXlJz1OMjbZ2hg45RJNyiTlyNh+7jY+nL0ci3ZdsnaViIjoMcgddTJuV8bZys1QMmPFJ598ogaLyw1UmaV///7o27evKWiRccbNmjXLsLyMGZZ0BHJzltTpp59+UrkLU5JxyNKVJ2ONpUVLgr20klMpSIJpOT4JpCRXlYxXzg7sZN4Va1ciN5PM4XKXgTRBZlVqgmRtp25Hlzvf4XX9erwRPwQDBw1Fdb/Umc+JiHKLmJgY1VojAYGMF6Lce70jLfguZotTLvJZ50pwcXaGzk7D54bv8ef+f+9OICIiov/GwCkXqVPCCy8Pn4n7eUohn10Eqp2cYFHSNSIiotyOgVNuY+8IQ7dZSNTs0NYYgGu7l1u7RkRERDaDgVMu5Fi8HjblTUpd77b5Q0z5bTfCHsRZu1pERETZHgOnXOpB/WE4aywCTy0cpQ+MwWe/n7R2lYiIiLI9Bk65VIcaxbHA50OEaW7YmlgNfxy/qSaaJCIioowxcMqlnB30mPROP3iOPIs9edogLlHD3qAwa1eLiIgoW2PglMvZObqhcRkf9XzP6SsyGZC1q0RERJRtMXAiNC6dD010R9H/cHf8vXwm9l9iyxMRESVNkTJ16tQnOhXbtm1T2cAfNV2LJS5duqS2d+TIEatcIgZOhKZlfVHbEIT8dndR4+Q4vL9gA0IiY3hmiIiyuV27dkGv16Nt27bIDpo2bYqhQ4c+NM2KTBYsmblzAgZOpMY7ebQehhNaSeS1e4DRmIc1h6/xzBARZXMLFy7E22+/jR07duDKlSvIjhwcHFCgQAHVSpQTMHAipc8zpVF24BIk2hnQQn8Y2pGfeWaIiLKxBw8e4Ndff8XgwYPRsWNH/PDDDw91j/3111+oVasWXFxcVMvP2bNnTWUuXLiAzp07I3/+/Gri4Nq1a2PLli2PnPS3Y8eOqdYlJCSooEgCuH79+iEgIADffvut2rcs0q2WXlfdzp071cTBUi+ZzLdNmza4e/euem/jxo1o2LAh8ubNC29vb7VPqWt2wcCJTAyFKiG87gfq+YthsxB9+xLPDhFRNrVs2TKULVtWLS+99BK+//77h6bR+vjjj/HNN9/gwIEDsLe3V8FPsvv376N9+/YqWDp8+LAKXp599tkMW64GDBigghrpdkv2xx9/qO306NFDBUz169fHa6+9psrIUrRo0Ye2I2OTWrRogYoVK2L37t2qtUz2m5iYaAoI33vvPezfv18FfjqdDl27doXRaER2YG/tClD24tXqfZzYtxKVjOdwdmFfbK41F6828lfdeUREuUbcg4zfs9MDBiczy+oAg/Ojyzq4PlYVFyxYoAImIWOcJICRQKNly5amMl9++aVq2REjRoxAhw4dEBMTAycnJ1StWlUtyb744gusXr0aa9euxVtvvfXQ/qTFqmzZsliyZAmGDRum1kmw9vzzz6sWK3UoDg6qFUlaoTIyadIk1Qo2a9Ys0zoJopI999xzDx2nr68vTp06hUqVKsHaGDhRKnZ6A/ZWHYeSh3rj1H1XzNh8EpGxRoxsX55niohyj3GFMn6vdGugd4p5Pr/yB+Kj0i9brCHwyvp/X0+tDESFpi4zJsLi6kmX2759+7Bq1Sr1WlqTevbsqbrMUgZOVapUMT0vWLCgegwJCYGfn59q2Rk7dizWrVuHGzduqG636OjoR46VGjBgAObOnasCJ9nO+vXrVbBmCWlxkmArI9It98knn2DPnj24c+eOqaVJ6sXAibKlzi0aY8ilWdh8M+mvpOUHrmJ423LQ63LGwD4iIlsnrTAS6BQuXNi0TrrpDAaDaayQkNfJkgdnJwciH374If788098/fXX8Pf3h7OzM7p37464uIznLn355ZdVy5V0scki6QoaNWpkUd1lP48i3XbSxTdv3jwUKlRI1VcCpkfV62liixM9xMfNEfOGdEd8ohE1Pt+Mu1FxOH41FNWKJSXKJCLK8UbeeHRXXUofnn9E2TRDiYcef8KKJQ3IXrx4sRq71Lp161TvSTfXTz/9ZFbLTGBgoBrQLeOHhHT1yWDuR/H29kaXLl1UF50ETq+88kqq96WrLnmsUkakFUxaqaS1K63Q0FCcPn0a3333nSkgkzFQ2QkDJ8qQQa9Du+J6tLowEQmbKgKvTefZIqLcwZJxR1lVNgPStSatSq+++upDuZGkxUhao6ZMmfKf25FWJunqkxYeaY2S7jFzBmAPGDBA3ekmAVLfvn1TvSctUHv37lUBmIx78vLyeujzH330ESpXrow33ngDgwYNUsHW1q1bVfedlJfgTLoDpWtRuuekhSs74V119EjP+txEK/1B1Li2BHMXL8bWMyE8Y0REViSBkYxjSi+hpLQ4yRiiQ4cO/ed2JLiSVAAy6FuCJ7mrrkaNGv/5uZYtW6qgRspLV1pKH3zwgUrIWaFCBeTLly/d8VJlypTBpk2bcPToUdSpU0fdiffbb7+pcVpyB93SpUtx8OBB1Wr27rvv4quvvkJ2YqelvXeRnqrIyEj1wx8REYE8efJku7MfGROPP8c9j+d1W3FN80HnxEn4c0RH1Z1HRGTL5O6yoKAglChRQt1lRuaJiopSAZMMRO/WrVuOuN6WfBezxYkeKY+TAfbtJ+C6XX4UsbuDkXbfY/2xf3N4EBFR7iDdeDdu3FBdehJkdOrUCbkRAyf6T13rlUPhV5bACB2e0wfi/qEVPGtERLmMdLsVLlxYZSuX1ibpWsuNGDiRefzq4m6NpIRoL96eguhQzmVHRJSbyMBvTdNw9epVlfk7t2LgRGbzaj8KZ+1K4r7mjH4z1uOl+XsRci+GZ5CIiHINBk5kNjt7R+yu9S3axY3H3ugi2HH+DqZu+YdnkIiIcg0GTmSRnq2eQa9GFVHQI+mOhA1HryHRyBszich28eby3EHLpCQCDJzIIjLZ78cdKiDww6Z4zelv/GQcjqMXON6JiGxP8nQkcns95XxR/7/OKaeheRy5c0g8PTH7hAd4w34NPBPu4OimjxBdbJEKqoiIbIUkasybN6+arFa4uLiY5nOjnNXSFBUVpa6zXG+57k+CgRM9Hqc8ONvgG9QJ6Ieqt3/Hm2PGIE+tHhjfrTLPKBHZjAIFCqjH5OCJcq68efOarveTYOBEj61G405YvPM59EtcgfGG+Wi/vxTONywOf193nlUisgnSwiTTh/j6+iI+Pt7a1aEsIt1zT9rSlIyBEz02B3sdWr8xBSGLz8E34himGmZi04l68G9ejmeViGyKfKlm1hcr5WwcHE5PpJB3Hvj2W4I4ezfU0p2D94FveUaJiCjHYuBET86zOB60+goJmg5XIhPw25HrCPznNs8sERHlOOyqo0zhWfdFDDnggN+uOgNLj6h181+uhZYV8vMMExFRjsEWJ8o0L7RtBidD0o+UI+KwcMdFnl0iIspR2OJEmaZ+KW8c+qQVrv1zFInL+uLXK80REV0LHs5PlmyMiIgou2CLE2UqFwd7lLl/EOV1V/GR/kcc2buVZ5iIiHIMBk6U+eq8hnOeTeFgl4gSW99Av9mb8feZWzzTRERk8xg4Ueazs4NT99m4quWDn91t9LwxEe8tO4KouASebSIismkMnChL+BUuhAtNZiAB9min34/OcesRcJYpCoiIyLYxcKIs07R5W9i3/VI9H2n/E04fCuDZJiIim8bAibJW3YEI82uNY1pJbAjScO7WPVwNi+JZJyIim2SnaZpm7UrkZpGRkfDw8EBERATy5MmDnCgxOhJNp+zG1cikMU6O9jqsfashyhbgZMBERGRb38VscaIsp3fOg0Et/p341yfhFuYHMjkmERHZHgZO9FT0rlsMW4Y2wKbyG/G34we4fmoXEhKNPPtERGRTGDjRU+OfPy/8DXfgaBePiYnf4PDZIJ59IiKyKQyc6Omxs4Ou6yyEGgqhqO42opf1R885O7HnYiivAhER2QQGTvR0OXviTof5iNEMaGx3BHWvLsT7vx5ltx0REdkEBk701JWt9gwOV/lUPR9qvxL+kXuwLyiMV4KIiLI9Bk5kFfWfeweo2Q86Ow1TDDMRcPQfXgkiIsr27K1dAcrF2k5E2NUzGHWtDvaeuo/qZYNR0MMJVYvmtXbNiIiI0sUWJ7IegxPcXvsD+12bIPRBHAb9eBDPzd6lsosTERFlRwycyKocDHp82KYs7HV28MVdNMQhfL/zEq8KERFlS+yqI6vrUasoOhSNg/33raFFR6D/SV8kdqkEvc7O2lUjIiLKXi1Os2bNQokSJeDk5ISaNWsiMDDwkeUDAgJUOSlfsmRJzJkz56EyK1euRIUKFeDo6KgeV69ebfF+ZQq/MWPGoFChQnB2dkbTpk1x8uTJVGWCg4PRp08fFChQAK6urqhRowZWrFjx2OciN3PNVwKGItXhZBeP8fFf4eCZIMTEJ1q7WkRERJkTOMXFxeHatWu4cuVKqsUSy5Ytw9ChQ/Hxxx/j8OHDaNSoEdq1a5fhdoKCgtC+fXtVTsqPHDkS77zzjgqUku3evRs9e/ZUAc3Ro0fVY48ePbB3716L9jtp0iRMnjwZM2bMwP79+1Vw1KpVK9y79+/4G9n22bNnsXbtWhw/fhzdunVT+5ZtkoV0Oui6zUWYoSCK6ULw4Jd+qDx6A5YfuMpTSURE2YdmoXPnzmkNGzbUdDpdqsXOzk49WqJOnTraoEGDUq0rV66cNmLEiHTLDxs2TL2f0sCBA7V69eqZXvfo0UNr27ZtqjJt2rTRXnjhBbP3azQatQIFCmgTJkwwvR8TE6N5eHhoc+bMMa1zdXXVFi9enGo7Xl5e2vz58zVzRUREaHIZ5JE0Lej4Li3qUx9NG51Hm/5xH63m55u0mPgEnhoiIsoylnwXW9zi1K9fP+h0Oqxbtw4HDx7EoUOH1CKtLPJoSYuVfL5169ap1svrXbt2pfsZaU1KW75NmzY4cOAA4uPjH1kmeZvm7FdatqQbLmUZ6fZr0qRJqro1bNhQtV6FhYXBaDRi6dKliI2NVd169HiKV6qP0BZfq+dv2f+GOlHbsenkLZ5OIiKyzcHhR44cUYFHuXLlnmjHd+7cQWJiIvLnz59qvbyWoCU9sj698gkJCWp7BQsWzLBM8jbN2W/yY3plLl++bHotQZN0zXl7e8Pe3h4uLi5qPFWpUqUyPG4JrGRJFhkZmWHZ3KpI475AzD/ArukYYP8HpuzvgGerFrJ2tYiIiCwf4ySDrSX4yCx2dnYPDcpOu+6/yqddb842M6PMqFGjcPfuXWzZskW1er333nt4/vnn1XinjIwfPx4eHh6mpWjRohmWzdVajEF4/Y/QJ34kAs+HovGkrRi/4bTpehMREdlE4DRx4kQMGzYM27ZtQ2hoqGoxSbmYy8fHB3q9/qHWpZCQkIdaepLJAO30yktrj7T6PKpM8jbN2a9sQzyqzIULF9TA8YULF6JFixaoWrUqRo8ejVq1amHmzJkZHvdHH32EiIgI03L1Kgc/p0tvj7xtRqCGfxH18kpYFL4LuIht525neG6JiIiyXeDUsmVL7NmzRwULvr6+8PT0VEvevHnVo7kcHBxUGoDNmzenWi+vGzRokO5n6tev/1D5TZs2qWDFYDA8skzyNs3Zr6QpkOApZRkZGyWpEJLLREVFqUcZ75WSBGUy3ikjMlYqT548qRbK2LiuldGtWiG87fA7Rtj/jKX7LLtzk4iIyKpjnLZu3ZppO5euLbmlXwIfCXjmzp2rUgIMGjTI1Dpz/fp1LF68WL2W9dLKI5977bXX1EDwBQsW4JdffjFtc8iQIWjcuLFqGevcuTN+++031ZW2Y8cOs/cr3XGSrmDcuHEoXbq0WuS5jGF68cUXVRkZ4+Xv74+BAwfi66+/Vi1ea9asUcGWDJynzFHUywWTG8QBZ35RYf77Z4shPKoK8ro48BQTEdHTp1nZzJkztWLFimkODg5ajRo1tICAANN7ffv21Zo0aZKq/LZt27Tq1aur8sWLF9dmz5790DaXL1+ulS1bVjMYDCrNwMqVKy3ab3JKgtGjR6u0BI6Ojlrjxo2148ePP5SaoVu3bpqvr6/m4uKiValS5aH0BP+F6QjMtGWsSlEQ/am3NnDifO2LdSe1uw9iLTrXRERET/pdbCf/WBpshYeHq5ae06dPq9YZGTDev39/NdiZLCPjwuS8yXgndts9gjERV2Z2hl9oIK5r3ugU+wUaVC2P6b2q80eOiIie2nexxWOc5O4xud1+ypQpKn+R3GEnGbZlnSV5nIgsotPDp+8i3DIUQWG7UMxxmIK/jl9ByL0YnkgiInpqLG5xkulJZGzPvHnz1N1sQvIoDRgwABcvXsT27duzqq45ElucLHT7HDC/JRAbgRWJjRHeaioGNM44bxYREZHVW5yGDx9uCpqEPJcUBfIeUZbKVwZ4/nsk6Bxx0Fgak7f8g7nbL+DE9QieeCIiynIWB04SiaU3Ca/kI3J3d8+sehFlzL8FHgw6iJVohai4RIz74wx6frcbIZHstiMiomwWOMkUI6+++qqabkSCpWvXrqk52qSrrlevXllTS6I0PHyLYljbsnAy6JAX91Ag/gqW7mcyUSIiymZ5nCRnkdxJ9/LLL6uxTUKSTw4ePBgTJkzIijoSpWtAo5LoX15D1MLOCH8QjTcP5MPbzf0fOWUPERHRk3isdATJmbNl2hH5uAwWl+SQZDkODn9C0XdhnNcCurALOGAsg611F6CSXz60rVSAARQREVl/cHgyCZQqV66MKlWqMGgi63H2hO7FZYjSuaGW7hxK7hmJwT8dxHfbL/KqEBGRdbrqunXrhh9++EFFYfL8UVatWpVZdSMyj09pRHScB4e1vfGcPhD/GAtjfqATXm1YAgb9Y/9tQERE9HiBkzRfJY8bkeCJY0gouylYoz0S4ycCGz7EMMMyBEUVxN9nKqFNxQLWrhoREeUgjz3GiTIHxzhlsvUfAPvn4bTRDx96TcfLz5RCw9I+KJTXObP3REREOUSWjnFq3ry5mqsuvZ3Ke0RW1XYCwmq8hRfjRuJE8AMMW3kM3WfvQnRcIi8MERE9MYsDp23btiEuLu6h9TExMQgMDHzyGhE9Cb09vDp9iY71KptW3YyIwm9HrvO8EhHR08vjdOzYMdPzU6dOITg42PQ6MTERGzduROHChZ+8RkSZ4PMulTCiXTkcWDUZdqfWYtq+L/BCHT+eWyIiejqBU7Vq1dSgcFnS65JzdnbG9OnTn6w2RJnINS4UjYK+hU5/H7dvfYXxf/iiYuG8eLZKQd7gQEREWRs4BQUFqWSXJUuWxL59+5AvXz7Tew4ODvD19YVer3+8WhBlBff80PVYhMQfn8dz+h2YuusrvJPQHRFRcehTvzjPORERZV3gVKxYMfVoNBot3wuRtfi3xO3G41Fg+3AMtV+Fq0ZffLfdGb3rFoNOx6lZiIgoi+eqSznO6cqVKw8NFO/UqdPjbpIoSxRoPgha4i3Y7ZyMCQ7z8HKEN3ZeqIxGpf9tNSUiIsqSwOnixYvo2rUrjh8/rsaJJKeBSk6KKQPFibIbuxafAOGXYDi5Ct8ZpqD7L8XQqEZlPF+rCMoVeHTODiIiosdORzBkyBCUKFECt27dUnPUnTx5Etu3b0etWrVUqgKibEmnA7rMRlSB2piT8CzORbliwY4gvDR/H6LiEqxdOyIiyqmB0+7du/HZZ5+pweE6nU4tDRs2xPjx4/HOO+9kTS2JMoPBCS6vbYBLiw/h5+WqVt25H4u1R27w/BIRUdYETtIV5+bmpp77+Pjgxo0bpsHjZ8+etXRzRE+X3oC3mpfG9mHN8GmrohhmvxQr9p7nVSAioqwZ41SpUiWVDFPSEtStWxeTJk1S6Qjmzp2r1hHZBE3DS5dHwsE+ECVCbqLDVB2q+Hlh9LMV4WRgWg0iIsqkFqdRo0aZUhJ88cUXuHz5Mho1aoQ//vgD06ZNs3RzRNZhZweH5iOQYGdAO/1+9LwzA7/su4JZ2y7wihARUYbstOTb4p5AWFgYPD09mY05i2dkpswXeXAF3H4fAB00fBXfA8uce2DXiBZwsLf4bwoiIsoF38WZ8u3g5eXFoIlsUp6a3aFrN0k9/9DwK5pFb8KmU//Ow0hERGTxGKdu3brBXKtWrTK7LFG2UPd14H4wEPgNxtvPx+tL8+DD5bUxrG1ZvPJMCWvXjoiIbC1wkuYrohyt+SeICr2GByf/RLDmhej4RIz74zQ6VC4I3zxO1q4dERHlpDFO9Pg4xikbSYzH+aCLCAh2wMIdQbgeHo2hLUtjaMsy1q4ZERFlk+/ix56rjijH0Rvg718W/v6Aj5sDlixbisC94Xi2aiGVoqBwXmdr15CIiKzM4sBJpltJnpcuo7nsiGxdO7d/0MZxAkJiPdD9Gw0R9t5Y8+YzKF+Qdz4SEeVmFgdOQ4cOTfU6Pj4ehw8fxsaNG/Hhhx9mZt2IrMYhf3lEuuSHX9RV/OgwDj3iPsXULefwXZ9avCpERLmY/eNM8puemTNn4sCBA5lRJyLrc8+PPK+tR+L81ijz4Dp+cJiEl06NxNWwCijq5WLt2hERkZVkWpa/du3aYeXKlZm1OSLr8ywGfd/fAGcvVNNdwHf2k9F37nZ89vspXAmNsnbtiIjIlgOnFStWqESYRDmKbzngpRVIsHfFM/qTGP7gKyzaeR6vLzkAo5E3pBIR5TYWd9VVr1491eBwyWYQHByM27dvY9asWZldPyLrK1wT+hd/QcKS7nBzdoQuXsOZ4HvYcf4OGpfJZ+3aERFRdg6cunTpkuq1TqdDvnz50LRpU5QrVy4z60aUbdiVbAL71zbjmQKV0XvdGfyw6xIW777MwImIKJdhAkwrYwJM23M+5D5aTd6KlrpD2GVfB6V83fHTgLpwdzJYu2pERJQdE2AmJiZi9erVOH36tOq2K1++PDp37gx7e+bTpJzPP58rfvT9Gc9E/oEp8c/h22vPYcGOIGYYJyLKBSyOdE6cOKGCJBnXVLZsWbXu3Llzqrtu7dq1qFy5clbUkyj7sLNDnbrPAJv/wLuGlYiDAYt2PYdBTUqpDONERJRzWdxVV69ePfj6+mLRokXw9PRU6+7evYt+/fohJCQEu3fvzqq65kjsqrNhgZOBv8aqp5/F98G1sv1Qu7gXetX1g5sjW1+JiHLid7HFgZOzs7NKdFmxYsWHWqJq166N6Ojox6t1LsXAycZtHQ8ETFBPR8W/gh8TW6FtxQKY06emtWtGRERZ8F1scR4n6Z67devWQ+ultclfZkclyk2ajkB8/XfU0y8M3+N5/TZsPBmMi7fvW7tmRESUBSwOnMaNG4d33nlHJby8du2aWuS5zGE3ceJEFbUlL0Q5np0dDK0/A+oOBnQGFCtUSK1esueytWtGRERZwOKuOsnbZPrw/xNhJm8i5Wt5Lnff0aOxqy6HkP8Dt05ie2R+vLxwn1rl6WJAlSJ5Me/lWnCwz7Qk/UREZEvpCLZu3fokdSPKmeSPhgKV0Ci/hprFPHHrylmUibmGv8/VwLL9V9CnfnFr15CIiDKBxYFTkyZNMmO/RDmStLQuer4o7Bb0h1N0CN6IexvfbXdGrzp+sNez1YmIyNY91j3T4eHhWLBggSkBZoUKFdC/f3/VzEWU27l5FQL8GwHHf8VMh2l4KwIYsNgNlQt74PXGJZlhnIjIhln8J7CkIihVqhSmTJmCsLAw3LlzB5MnT1brDh06lDW1JLIlOj3QdQ5QpSfsYcQMwzQ4/7MO0/8+j5GrT1i7dkRE9DQHhzdq1EilHZg3b55pipWEhAQMGDAAFy9exPbt25+kPrkOB4fnYMZExK8aBMOJX5EIHd6OewsbtHrY/mEzFPVysXbtiIjoaSXAPHz4MMqVK5dq/alTp1CrVi1ERUVZsrlcj4FTDmdMBH57Czj6swqeBsS9j5INuuHVhiXg7eYAR3tO0UJElKMTYMoGr1y58tD6q1evwt3d3dLNEeX8brvOM4BqvRHtURpHjaXUhMANJvyNLjN3ISaeKTuIiGyJxYFTz5498eqrr2LZsmUqWJIEmEuXLlVddb169cqaWhLZevDUaTpcB/2JKmVLmVafvhmJn/c+/EcIERHloLvqvv76a3Un3csvv6zGNgmDwYDBgwdjwoSkObuIKA2dHnbOnljYtzbOBN/D7YDvsPp4GBbscEaf+sVgYKoCIiKbYPEYp2QylunChQsqS7gMFndx4WDXx8ExTrnQpZ3AD+1hhB0+ju+PPZ6dkM/NEV89XwXFvF2tXTsiolwnMivHOCWTQClv3rzw8vJi0ERkCb/6QJ3XoYOG8YYFaHZ3BfZdCsNHq47zPBIRZXMWB07SPffJJ5+oyKx48eIoVqyYej5q1CjEx8dnTS2JchKZ77HdJBgbDFEvPzUswZv6Ndh1IRQHL9+1du2IiCgzxzi99dZbWL16NSZNmoT69eurdbt378aYMWNUMsw5c+ZYukmi3MfODrpWYwEHV2DbOHxo+BXOdrF4b5kzWlUogHaVC6o574iIyMZbnH755Rf88MMPGDhwIKpUqaIWeb5w4UL1nqVmzZqFEiVKwMnJCTVr1kRgYOAjywcEBKhyUr5kyZLpBmorV65U08A4OjqqRwn0LN2vjN2SYLBQoUIqd1XTpk1x8uTJh7YjQWPz5s3h6uqqui6lXHR0tMXngXLpxMBNhwOtv1Av37L/DX7hezF/RxBeXrAXwREx1q4hERE9aeAkgYZ00aUl6xwcHCzalqQ0GDp0KD7++GOVVFOykrdr1y7dPFEiKCgI7du3V+Wk/MiRI/HOO++oQCllICMpE/r06YOjR4+qxx49emDv3r0W7Vda1GQqmRkzZmD//v0oUKAAWrVqhXv37qXaV9u2bdG6dWvs27dPlZMWOZ10xRCZq8HbQIdvcKLkq7hXqJFa9SAuEdP+/ofnkIgou9EsNHbsWK1Xr15aTEyMaZ087927tzZmzBiLtlWnTh1t0KBBqdaVK1dOGzFiRLrlhw0bpt5PaeDAgVq9evVMr3v06KG1bds2VZk2bdpoL7zwgtn7NRqNWoECBbQJEyakOkYPDw9tzpw5pnV169bVRo0apT2JiIgIuatRPRKJvRdDtXLDV2glh/+mLdoVpF4TEVHWseS72OKmEWmhWbduHYoUKYKWLVuqRZ7//vvvqoWnW7dupuVR4uLicPDgQdVak5K83rVrV7qfkRaetOXbtGmjJh5OHpieUZnkbZqzX2nZCg4OTlVGuv2aNGliKhMSEqJasXx9fdGgQQPkz59fvb9jx45HHndsbKy67THlQpRSnSLOWOnxLWYZvsWXvx1Gj+92Y+XBazxJRES2ODhcxvE899xzqdYVLVrU4h3LQPLExEQVcKQkryVoSY+sT6+83Okn2ytYsGCGZZK3ac5+kx/TK3P58mX1XCY0FjIOSpKCVqtWDYsXL0aLFi1w4sQJlC5dOt1jGD9+PMaOHWvGGaJc6+YxlIs/jQr6OHyPSXg9/j18veksnq1aCA727AYmIrKpwOn777/P1ApIFvK0g7LTrvuv8mnXm7PNJy1jNBrVowyMf+WVV9Tz6tWr46+//lID5SVASs9HH32E9957z/RaWpweJ/CkHMyvLnR9VkL75UU0iDuF5fbj0TviQ6w5fB31S3mjcF5n6HQZ/x8hIqKsY7U/X318fKDX6x9qXZIusLQtPclkgHZ65e3t7eHt7f3IMsnbNGe/sg3xqDLSuiXkrr2Uypcvn+Hg9uQuP8lKmnIhekiJxrDr9zvg4o3y2gUsdxiLqSv/RqNJWzHwx4OmPxiIiCiXBE5yB56kAdi8eXOq9fJaxgylR/JGpS2/adMm1KpVS82X96gyyds0Z7+SpkCCp5RlZGyUpEJILiN3EUqqgrNnz6bazrlz51RSUKInVqg60P9PGPMUQSndTax0HAN/u2vYfOoW/jodwhNMRGQNmhUtXbpUMxgM2oIFC7RTp05pQ4cO1VxdXbVLly6p9+Uutz59+pjKX7x4UXNxcdHeffddVV4+J59fsWKFqczOnTs1vV6v7og7ffq0erS3t9f27Nlj9n6FfE7uolu1apV2/PhxdSdhwYIFtcjISFOZKVOmaHny5NGWL1+u/fPPP+oOOycnJ+38+fNmnwPeVUf/KfyqZpxeW0sYV1SbtWytVmz4Oq3DtO1afEKiFpeQyBNIRPSELPkutmrgJGbOnKkVK1ZMc3Bw0GrUqKEFBASY3uvbt6/WpEmTVOW3bdumVa9eXZUvXry4Nnv27Ie2KYFM2bJlVXAkaQZWrlxp0X6TUxKMHj1apSVwdHTUGjdurAKotMaPH68VKVJEBXT169fXAgMDLTp+Bk5klgehmnb9kHbnXoxWbtQGFTzJUm/cFi3o9n2eRCKiJ2DJd7Gd/PO4rVUxMTEqISY9nRmZicSUzeew5++18LKLxAZjXTQpkw+L+tfhySEiegrfxRaPcZK7yT7//HMULlwYbm5uptvyZeLfBQsWPG6dichMQ2vo8ZPbFMxymIa++k0IOHcbey+G8vwRET0FFgdOX3zxhZqrTqYkSTnFSuXKlTF//vzMrh8RpWHnWRz2VXvADhrGGn7AcPtf8MLcXaj22SYs2nWJ54uIKDsFTpLkce7cuejdu7e6rT+ZTPZ75syZzK4fEaWl0wMdJgPNR6mXg+1/x2TDbDyIisZn607hcugDnjMiouwSOF2/fh3+/v7pduElT3tCRFlMErE2/hDoPAuazh5d9Tvxq9s3cDE+wHfbk7rPiYgoGwROFStWRGBg4EPrly9frjJnE9FTVL037F5cBji4oXrCUQywX4+f915Bo0l/44W5uxERxT9miIisOuXK6NGj0adPH9XyJK1Mq1atUkkgpQtPJv8loqfMvyXQbz2weyZO3+8PnA7D1bBotUzefBZjO1fiJSEiyiSPlY7gzz//xLhx43Dw4EEVPNWoUQOffvopWrdunVn1yjWYjoAyU0x8ItYfu4kT18Kwa89OXNQVw9/vN0VRLxeeaCKiTPgufqI8TvTkGDhRppP/0htHIH7vAgyNG4w/UR95nA0Y/WwFdK5WmCeciOhp5nEqWbIkQkMfzhkTHh6u3iMiKzMmABHXYEA8ZjpMw+t2axD2IBYjVh5HcESMtWtHRGTTLA6cLl26hMTExIfWx8bGqnFPRGRlegPQYzFQd7B6OcywDDNcFyAhPlaNeWIjMxHRUxgcvnbt2lRjnKRJK5kEUn/99ReKFy/+BFUhokzN9dRuAuBdCtgwDB0T/4a3IRiDDgzFrweuoWYxTzVNi5ujxfeHEBHlamaPcdLpkhqn7OzsHvqL1WAwqKDpm2++QceOHbOmpjkUxzhRlvtnC7C8HxB3DyeMxfFs3BfQoMPrjUtiZPvyvABElOtFZsUYJ7l7ThY/Pz+EhISYXssi3XSSkoBBE1E2VLol8Oqf0PL6IbHhBxjWtoJa/eOeywi5F6PuxCMiIvPwrjorY4sTPTXxMYDBSbUYd5y+A1dv3EAk3KCzA8Z0qoiX67OrnYhyp0gLWpwea4DDgwcPEBAQgCtXriAuLi7Ve++8887jbJKIsprBydTdPrGlJ3yXvYRfEpthSkJ3fLH+NJqV9WW+JyKi/2Bx4HT48GG0b98eUVFRKoDy8vLCnTt34OLiAl9fXwZORDag0r1dgF04htivRh23UPQP74exv5/E641Lwd/XDV6uDtauIhFRzkhH8O677+LZZ59FWFgYnJ2dsWfPHly+fBk1a9bE119/nTW1JKLMVfd1oNMMQGdA/ZjtWOEwFqdOn0KP73ajzdTtuBXJfE9ERJkSOB05cgTvv/8+9Hq9WmRgeNGiRTFp0iSMHDnS0s0RkbXU6AP0XQu4+KCi7jLWOo5CLbszuH0vFlO3nON1ISLKjMBJUg/IGAmRP39+Nc5JyKCq5OdEZCOKNQBe3woUqAwfu0j86jwez+iOY9n+q5i59Tx+3X8VcQlGa9eSiMh2xzhVr14dBw4cQJkyZdCsWTM1ua+McVqyZAkqV66cNbUkoqyT1w/o/yew5g3owi7Axb4+jOfv46s/z6q3D1wOw6TuVXkFiIgeJx2BBE337t1TQdPt27fRt29f7NixA/7+/li4cCGqVavGE2sBpiOgbEN+FcSE41a8M95ddgRXwx4g/G4Y7sEFG4c2QrkCj75Fl4goN3wXM4+TlTFwomxr+1e4HbgQL94filjPMqhSxAMty+dHl+qFrV0zIqLsnzk8WfPmzREeHp7uTuU9IsoB4h4Ah39EvvgbWOPwKcqHB2DdsZsYuuwIdl8ItXbtiIisxuLAadu2bQ8lvRQxMTEIDAzMrHoRkTU5uAID/gaKN4KrXQy+c5iCCR6roIMRo9Yc54BxIsq1zB4cfuzYMdPzU6dOITg42PQ6MTERGzduROHCbMInyjFcvYE+q4HNo4E9M/FC7AoUdz6HwbffQJlRG1A4rzNmv1QDVYrktXZNiYieGrPHOOl0OlMagvQ+Iskwp0+fjv79+2d+LXMwjnEim3B8BbD2bSA+CheMBdE2biLiYY+S+VyxaWhj2OstbrwmIsrZc9UFBQWpgKlkyZLYt28f8uXLZ3rPwcFBTbciCTGJKAeq3B3wrQAs6w19udcwrWAdDFtxDBdvP8CcgAuo7uepBo+7OxmsXVMioizFu+qsjC1OZFPiogAHF/V03vaLWLxhG25pXoiDQc1xt+bNZ+Dm+FhzhxMR5cy76hYtWoT169ebXg8bNgx58+ZFgwYN1Jx1RJSD/T9oEq9Uz4PVrhOxzOFzFEAozofcx+RNnKqFiHI2iwOncePGqfFMYvfu3ZgxY4aap87Hx0dNAExEuYN9+EX42Eejuu48tnuMVlO1LNwZhN7z92D8htOIjku0dhWJiKzfVefi4oIzZ87Az88Pw4cPx82bN7F48WKcPHkSTZs2VdnEyXzsqiObFhYELOsD3DoOI+wwLaErpiV0gxE69KpTFOO7VbF2DYmIrNtV5+bmhtDQpAR4mzZtQsuWLdVzJycnREdHW7o5IrJlXiWAAZuBGi9DBw1D7VfhD8/J8EGEmij4bPA9RMUlWLuWRESZxuLAqVWrVhgwYIBazp07hw4dOqj10uJUvHjxzKsZEdkGgzPQaTrQdS5gcEG56EP4Nt9vMGpAm6nbUeHTP1XSTAsbt4mIckbgNHPmTNSvX191ya1cuRLe3t5q/cGDB9GrV6+sqCMR2YKqPYHXtgL+rVCy9xT4ujua3vpxzxVsOR1i1eoREWUGpiOwMo5xopwqIjoe52/dw+0tU/HR+fJw8vDFi3X84Oftgk5VC5kS6hIR5cgEmOmpXLky/vjjDxQtWvRJNkNEOZCHswE1Q38Hrn2Lak7eeDPyLXyzOUa9F/YgDq88U8LaVSQistgTzZNw6dIlxMfHP8kmiCgnK1IL8C6t8jz96vg5PvX4Q00UPHnzORU8ERHZGk4wRURZJ39F4PWtQOXnoYcR/WN/xGq3SXCJCUGNzzej3CcbMD/wIq8AEeWOwKlRo0amZJhEROlydAe6zQM6zwIMrqiacAwbHT9CE91RxMQb8cX60zh6NZwnj4hsAgeHWxkHh1Oucuc8sLI/tJvHcLT5Ynx7oQC2nr2NUvlc0aVaYZQp4I42FQtYu5ZElMtEWjA43KzAae3atWbvvFOnTmaXJQZOlAslxAJB24HSrXAjPBqtp2xHQuwDxCApfcGUnlXRtXoRa9eSiHKRyMwOnHS61D16chtxyo+lvK04MZHzU2XVxSLKiQ4fPoCSv3fD9w4vYmp4Q/i4OWLLe03gYK+Di8MT3fhLRGSdKVeMRqNpkWlWqlWrhg0bNiA8PFztRFIS1KhRAxs3bjSvhkRE/1c9ZDU8jBEYGjMbS1ynI+F+KKp9thmVRv+Judsv8DwRkW2PcapUqRLmzJmDhg0bplofGBiI119/HadPn87sOuZobHGiXM9oBPbMAraMAYzxCNE88WH86wgwVlWnZnH/OmhcJl+uP01EZKOT/F64cEFtPC1ZJ3mdiIgsIkMBGrwFDNgC+JSBr91dLHKYiF8KL4cTYtHv+31455fD+HHPZc53R0RWZ3HgVLt2bQwdOhQ3b940rQsODsb777+POnXqZHb9iCi3KFQNGLgdqDtIvawfuhrv5t2uJgtee/QGRq05gZ/3XbF2LYkol7N45OXChQvRtWtXFCtWDH5+fmrdlStXUKZMGaxZsyYr6khEuYXBGWg3ESjTBtg3D22aj8HFwCsIuvMA+y6F4cv1p+Ht6ghvNwfUKubJ+e6IyDbyOMlHNm/ejDNnzqjnFSpUQMuWLflL7DFwjBPRfzPGx+Lvb/rgy4jWCNIKqnW96vhhfLfKPH1ElP3SEWQkJiYGjo6ODJieAAMnIjNsmwBsG49oOGJC4ktYFN9cfn1hWq/q6FS1EE8hEWXfweGSkuDzzz9H4cKF4ebmhqCgILX+k08+wYIFCx6/1kREGan+ElCiCZwRi7H6BdhWYDoK4Y4aNN7s6234cLlM38IcckSU9SwOnL744gv88MMPmDRpEhwcHEzrK1eujPnz52d2/YiIAI8iQJ81QNsJgL0TiofvwRbn4eip34qgO/ex/OA1jFl7kmeKiLJf4LR48WLMnTsXvXv3hl6vN62vUqWKGvNERJRlaQvqDQYG7QCK1oWLFo2JhnlY4b8ZMnnB0v1X0X32LhVAhUTG8CIQUfYInK5fvw5/f/90u/Di4+Mzq15EROnzKQ28sgFo/SXg7IVa3YZgSIvS6q0Dl+/ih12XMGDxAcQnGnkGicj66QgqVqyosoRLOoKUli9fjurVq2dm3YiI0qfTJyXNrPUK4OCKt5tr0NvZodCFZZh2uRiOXYMa91SvpDcalPKBn7cLzyQRWSdwGj16NPr06aNanqSVadWqVTh79qzqwlu3bl3m1IqIyBwOrupBr7PD2/63gcCv0cnZBZ9E98LSI82w5sgNuDvZY/UbDeDv685zSkRPv6vu2WefxbJly9TEvnZ2dvj000/V/HS///47WrVq9eQ1IiJ6HG75gaJ1YEh4gAmG+Vjh+hUqOIXhXkwCBv94CNfDoxEZw+EERPRkniiPEz055nEiykTGRGDvHOCvz4CEGGj2zphmfB7TolohEXrVMvX181XQtXoRnnYiejp5nEqWLInQ0NCH1oeHh6v3LDVr1iyUKFECTk5OqFmzpho/9SgBAQGqnJSX/c2ZM+ehMitXrlTZzCU5pzyuXr3a4v1KPDlmzBgUKlQIzs7OaNq0KU6eTP92Zynbrl071QLHaWeIrDz2qf6bwOBdQPFGsEuIxhDjYvzo9JX8T0WiUcPHq0/g0p0HvExE9FgsDpwuXbqExMSHE83FxsaqcU+WkC4/mTD4448/xuHDh9GoUSMVgMjcd+mRZJvt27dX5aT8yJEj8c4776hAKdnu3bvRs2dPNQ7r6NGj6rFHjx7Yu3evRfuVPFWTJ0/GjBkzsH//fhQoUEB1Rd67d++hek2dOpXZ04myE+9SQN/fgU4zACcP1Hx2II6PaYO6JbwQFZeIlpMDUO6TDfho1XEYZRZhIqLM7qpbu3ateuzSpQsWLVqkmrSSSSD1119/qfnrZKC4uerWrYsaNWpg9uzZpnXly5dX+xg/fvxD5YcPH67qIWOqkg0aNEgFSBIwCQmapMltw4YNpjJt27aFp6cnfvnlF7P2K6dEWpokuJJ9JgeG+fPnx8SJEzFw4EDT52TfHTt2VMFVwYIFVeuWbMdc7KojymIPQgEXL0iyJxnnNH7aDNyNScROY9I8d8PblsPgpqV4GYhysUgLuurMvqsuORiQ7qi+ffumes9gMKB48eL45ptvzK5kXFwcDh48iBEjRqRa37p1a+zatSvdz0hwJO+n1KZNGzXVi+SQknpImXffffehMtIqZO5+pWUrODg41b6k269JkyaqTHLgFBUVhV69eqlWKWmRIqJsyNXb9LSwYyymOc+DzngLp3w74MUrnTBx4xmsPHQNfl4uatLg/HmcrFpdIsohXXWSekAWPz8/hISEmF7LIq0x0tIkLS/munPnjmqpklaclOS1BC3pkfXplU9ISFDbe1SZ5G2as9/kx/+qmwRoDRo0QOfOnc0+bjlXEtmmXIjoKdHpoaso/1/tUCFkPQJdh6OrLhDnQ+7h7zMheG3xAcQmcM47IsrEMU7SGuPj44PMIi1YKUk3Wdp1/1U+7XpztvmkZaTL8O+//za1ZJlLugKlOTB5KVq0qEWfJ6In4OgOtP8KeHUTkK883BPDMcVhNvYUnILqTsE4di0CnabvxIBFBxBw7jZPNRE9XlfdtGnT8Prrr6s70OT5o8hgbXNI8CVz3aVtXZLWrLQtPcmkOyy98vb29vD29n5kmeRtmrPf5G43KSPjltIrI0HThQsXkDdv3lTbee6559Rg823btqV7DB999BHee+8902tpcWLwRPSUFa0DDNwO7J4OBHyFAncPYKXuKBrbfYWzt4Czt+5h+7nbWDqwHmr4efLyEJFlgdOUKVPUpL4SOMnzjEhrjLmBk4ODg0oDIAPKu3btalovrzPq+qpfv75KtJnSpk2bUKtWLTW+KbmMbCPlOCcpI11q5u5X0hRI8CTrkqeRkbFRkgpBBocLGSM1YMCAVHWpXLmyOj+SJDQjMlZKFiKyMnsHoNH7QKXuwIbh0Dm64bOK7XHg0l3sCwpT8971/2E/WlfIj5L53PBqwxIw6C1upCeinEazoqVLl2oGg0FbsGCBdurUKW3o0KGaq6urdunSJfX+iBEjtD59+pjKX7x4UXNxcdHeffddVV4+J59fsWKFqczOnTs1vV6vTZgwQTt9+rR6tLe31/bs2WP2foV8zsPDQ1u1apV2/PhxrVevXlrBggW1yMjIDI9HTufq1astOgcRERHqc/JIRFYUH2N6eu9WkPbXZ+20Z0Z8rxUbvk4tI1cd4+UhyqEs+S62eK66zCSpAySZ5meffYabN2+iUqVKaiqX5AmEZV3K3ErSEiTvS2vSzJkzVcoA6TqU7rFk0rK0dOlSjBo1Cp988glKlSql8jZJCgJz9yuGDRuG6OhovPHGG7h79676vLRcubtzviuiHMn+35Zgt4AxaJ64Ew2dDmCTz8t492oj/LT3ipq+pWBeJ7xYxw/FvJPmySOi3MXiKVe6d++uusbS3s7/1VdfYd++fVi+fHlm1zFHYx4nomzo9llg/fvApaQZBUKdi2NoxAsINFZRr/O5O+L3txqigAdTFxDltu9iiwOnfPnyqYHRMp4npePHj6Nly5a4devW49U6l2LgRJRNya/GY78Cmz4GHiTdYXfMtQE+j38J+yPzolQ+V7SpWADlC+ZBxyoFOXsAkQ3L0rnq7t+/rwZYpyWDs5mTiIhyDEk9UrUn8NYBoN6bgM4eVR7swvwqZ+HuZI8Ltx9g1rYLePuXw5gXeNHatSWip8TiwEnGA8mYobRkXJFMqEtElKM45wXajkuaOLhKT3i0GoYfXqmj7rZrU0KGiWqYsOEMRq05jimbzyEkMsbaNSaiLGRxV50kfpTB2C+++CKaN2+u1sk8dTIPnIxvsmSeNmJXHZHNMhqhzW+BS+HxeOtuL5zUiqvVxb1dsObNZ5DX5eGWeSLKhWOcxPr16zFu3DgcOXIEzs7OqFKlCkaPHq3mcqOsu1hElI0EnwAWtALio6DBDjvytMeY+91wIcoZ5Qq4o34pb5U889mqhaxdUyKyduBEmYeBE5ENi7gObP4UOLFCvUx0yIPJMc9iflwrxCKpxWlsp4ro2yCpNYqIcuHgcMmNJJPkpiU769Wrl6WbIyKyXR6Fge4LgFc2AAUqQx8XiQ91P2GX+wh0L2VURcb+fhJv/nQIn6w5gevh0dauMRE9IYsDp8WLF+OZZ55R87Qlk3nZJD3BpUuXnrQ+RES2p1gD4PUAoPMswL0gvPMVxFevtke36oVh1ID1x29iyZ7LeGHuboTej7V2bYnoaQZOx44dQ/HixVGtWjXMmzcPH374IVq3bo1+/fphx44dT1IXIiLbpdMD1XsDbx8Eun8PO50e47pVxsiWRbGy8FLUzBOOq2HReHb6DvRZsBfztl+UKa+sXWsistBjj3H6+OOPMX78eNjb22PDhg1o0aLF42wm1+MYJ6Icbut4IGACjDoHLE5sjcmxnRAJN/XWkBal8W6rMtauIVGuF5nVg8OnT5+O4cOHo2vXrjh48CD0ej1+/vlnVK1aNdeffEsxcCLK4W6dBP78GLi4Vb2MtvfAtvz9MORCDcTBgCKezvBxc8RnnSuiSpG81q4tUa4UmZWDw9u1a4exY8eqsU4//fQTDh8+jMaNG6NevXqYNGnSk9SbiCjnyV8R6LMa6L0CyFcOzgkRaHf9W+zLMwJddDtw7W40jlwNR58F+3Dx9n0kyqAoIsq2LG5xatWqFRYtWoRChQo9lNtpwIABuHnzZmbXMUdjixNRLpKYABxeAmybANwPxq1iHbGvxldqypZj1yJMxdpWLICpL1SDk0Fv1eoS5RaR1srjdOfOHfj4+GTW5nIFBk5EuVBcFLB3DlCxC+BVErfvxeLtWasRH34DB7WyqkirCvkxrE1Z5HE2IH8eJ2vXmChHi8zKrjoRGBiIl156CfXr18f169fVuiVLluDMmTOPV2MiotzEwQVo9J4KmkQ+d0f8WHIzVjqOxb4S81BRfx2bT91CqynbUX/8X/hp72Vr15iIHjdwWrlyJdq0aaOmWpHxTbGxSTlJ7t27p6ZhISIiCxmNsHdyB+z08L25FesMwzDVcQ4K47bKAyXJM6f99Q+W7b+CW5xEmMiqLO6qq169Ot599128/PLLcHd3x9GjR1GyZEk1b13btm0RHBycdbXNgdhVR0Qmd/4B/v4cOPWbeqnpDNjj2RFDr7fALXipdd6uDlg+qD5K5ktKaUBE2byr7uzZs+ouurRkR+Hh4ZZujoiIkvmUBnosBl77GyjRBHbGeNQPXY3xJY+iYqE88HV3ROiDOPScuwfv/3oUkzedRVRcAs8f0VNkb+kHChYsiPPnz6vs4SlJ1nBpeSIioidUuCbQdy0QFAjsmYXm3UajuaM77tyPxXszluJkuAErDyUNkzh8NRwL+taGg/1jDVkloqwOnAYOHIghQ4Zg4cKFsLOzw40bN7B792588MEH+PTTTy3dHBERZaREo6Tl/3xcHTA/z3wg7h/s8+2OD683QeA/d1Bn3Ba4GPTo26A4BjYpxfNJlJ0Cp2HDhqk+wGbNmiEmJkZ12zk6OqrA6a233sqaWhIRERAVCgd7PWCMRsPgJQh0XIU5aIO5UW1xA24Yv+EMouIS0a5yAeR3d4KnqwPPGlEme+w8TlFRUTh16hSMRiMqVKgANzcOVHwcHBxORBaRX9nnNgJbvwSCj6tV8QZ37Pd9HoMv1EPE/+fBczLoMLt3TTQr58sTTJRdE2CS5Rg4EdFjMRqBM+uAbeOBkFNq1bpyX+Hd40WQYNRUfOVor0P/hiVUN95zNYugUF5nnmyidDBwsiEMnIjoiQOo02uTlm7zkaAB8tfwt/PmYfklF1MaA7kj79eB9VHcx5UnnCgNBk42hIETEWW6uChoUysjMToCu9zbYnpcR+wPd4ergx4FPJxQvmAefNm1MjycDTz5RHgKU64QEVE29uA27PKVhb0Wj8aRv+PX2Dcwx20+fOOv4cLtB1h37CZeXrAXh6/cxembkeCIDSLzcYyTlbHFiYiyzKWdwPavgItb1UsNOlwv3BZv3GiLY9H/TsjeoUpBTO1ZDQY9/5am3CmSg8NtBwMnIspy1w4A278Gzm1QL889twWvb7yPkHuxKn2BqF3cEyV93FCzuCeer1lE5ekjyi0iGTjZDgZORPTU3DwGBAUADd42rbq44lN8e0TDuoTaSIRerRvYuCSGtS2nuvDs2QpFuUAkAyfbwcCJiKwm/ArwbTVAS0SYQyEEePfEiKCqiEVS4kxngx4jO5RHn3rFeJEoR4tk4GQ7GDgRkdVE3wX2fpe0RIclrTJ4YnZ0CyxKaG1KptmvQXHkz+OE+qW8Ua1oXl4wynEYONkQBk5EZHVxUcDhH4Hd05NaoSQ9lMEFS0tNwsgjSXmghM4OmNyjGrpUL2zFyhJlPgZONoSBExFlG4kJwKk1wM6pQOhFaO+dxJIjEfjrdAi0qFBsvyYDye3g7eoAdyd7jOlUEU3LckoXsn0MnGwIAyciynZkvpbQC4CPv+m1Nq85bobdw9eRLbDOWB9xMMCgt8NbzUojr4sBjcvkQwlmJScbxcDJhjBwIqJs7+4lYGY9ICFavYx3zoeNzh0x+kZdhCEpy7LMi7egb200LP1vfigiW8HAyYYwcCIimxAVBhz8Adg3D7h3Q61KsHNAgFMz/KTvir/vJAVQXq4OahnXtTLqlPh3fBRRdsbAyYYwcCIim5IYD5xcA+yZCdw4rFbFP78Ebx8qjI0ng03FpAVqYJNSKqVBm4r5UTJf0h16RNkRAycbwsCJiGx2HNTVvcCxX4H2X0Gz0+HE9Uh4nPgBB84EYVxwHdyBhyoqwdO8l2uhdglPGHQ66OT2PKJshIGTDWHgREQ5qjVqSiXgfjAS7eyx36UxftW1w6rbhdTdeMLX3VHNi9fAn2OhKPtg4GRDGDgRUY4KnE6sTBoHdf2AafVVx9KYfr8p1iY2QAwc4WCvQ+eqhdR0Ls/XKoIafp5WrTZRJDOH2w4GTkSUI8n4p33zgRMrgIQYtep+1f54/35v/HnylqmYpDT4oksllMnvjiKeLsjn7mjFSlNuFcnAyXYwcCKiHH833uElwP4FwAs/IdG3Epbuv4L7l48i4dYpfHu9nMoJJRz0OkzsXhldqxexdq0pl4lk4GQ7GDgRUa5gNAI63b+v17wJHPkRUfYe+DW+EVbZtcSxmKQs5P6+bjDodXi1YQl0r8kgirLXd7H9U6gPERHldimDJuFdCnAvBJd7N9DPbh36YR0u+1THN2EN8GdIbcTCAR8sP4rj18JRqbCH6saTSYaJrM1O0+SeUrIWtjgRUa6eG+/85qTEmv9sAjSjWh3hWREz/OdjXmBQquK96/pheLtyqkvPyaC3UqUpJ2JXnQ1h4EREJNHSdeDwj8ChxUCtfkDjD7Hh+E2sPRSEiqGbMTukIh7AWZ0qOzugW/Ui+LJrJQZQlCkYONkQBk5ERCkYE4HEOMCQFCSpLOXL+yJB74yNxjr4MbYR9hrLQYMOfl4uKOrlDP98bni/TVnkcUoaZE5kKY5xIiIi26TTA7r/B03JvErBPuwCOiIAHR0CEO1aFN8/qI+f7z6DnWH5sPN8KALP30Gv2n5wctCjfaUC8HZjWgPKGhzjZGVscSIiMmd6l33AkZ+AE6uAuHtJq2GHFQ3X4+t90bgVGWsqXiCPE8Y/V1m1SBXycIazA8dD0aOxq86GMHAiIrJAXBRwZl3SeKj4aGDAZoRExmBe4EWUurIcW0M9sel+CdWVJzycDfj2hWpoWjYp1QFRehg42RAGTkREjykhFrD/f5dc9F3g6zJqfFSYoQDWJDTAzzH1cF5LygPl4+YIN0c9BjUphRfq+PGUUyoMnGwIAycioky6K2/rl8CptaauPHHduQy+j6yt5skLQdKceHWKe8HH3UHNkffKMyWg1yVNQEy5VyQzh9sOBk5ERJnclXduA3BseVKOKGOCWh1c71MsRgfMDrighkwlq1PCCy3K+SKPswHPVi0EN0fmhc6NIhk42Q4GTkREWeRBKHBqNXB8BfD8D4B7AZy+GYmwPT/C++pfmHa7OjbHV0b8/yfRKJnPFcMkrYGzARUK5kFeFwdemlwikoGT7WDgRET0lC3qBAQFqKfRencccWuERZE1sDm6LBKhNw0qH9+tMuqW8IKroz0TbeZwkQycbAcDJyKip+zGEeD48qTl/i3T6vt6D2zVP4NRcf0QEZPUxSdkipchLUtjcJNS0HE8VI7EwMmGMHAiIrJilvLLO4GTq4FTvwFRoUCJJojtvRpf/3kWP++9gmLxF3BG84MROjgb9ConVPNyvvj02QrMVJ6DMHCyIQyciIiyyYTDlwIBvQNQ/JmkdfeCoX1TDjGOPlgZUxOr4+rikFZa5YjycnVAUS8X+Lo7YnjbsvD3dbf2EdBT+i5OyhBmRbNmzUKJEiXg5OSEmjVrIjAw8JHlAwICVDkpX7JkScyZM+ehMitXrkSFChXg6OioHlevXm3xfjVNw5gxY1CoUCE4OzujadOmOHnypOn9sLAwvP322yhbtixcXFzg5+eHd955R510IiKyMXp7oFSzf4MmEXIado554Bx7Gy/ZbcRKx7E47fkBvnL9CWWij+DE1VBsPnULHafvwMAlBzBk6WH8eTLYmkdBT4FVA6dly5Zh6NCh+Pjjj3H48GE0atQI7dq1w5UrV9ItHxQUhPbt26tyUn7kyJEqWJFAKdnu3bvRs2dP9OnTB0ePHlWPPXr0wN69ey3a76RJkzB58mTMmDED+/fvR4ECBdCqVSvcu5eUH+TGjRtq+frrr3H8+HH88MMP2LhxI1599dUsPWdERPSUSCD14T9Ar2VAlRcAB3c4RQfj+cT1WOrwBVY2uokGpbwRE2/Enydv4bcjNzBwyUEMWLQfn687hYU7ghCbkMjLlcNYda66unXrokaNGpg9e7ZpXfny5dGlSxeMHz/+ofLDhw/H2rVrcfr0adO6QYMGqQBJAiYhQZM0uW3YsMFUpm3btvD09MQvv/xi1n7llEhLkwRXsk8RGxuL/PnzY+LEiRg4cGC6x7N8+XK89NJLePDgAeztzcsFwq46IiIbER8DXPg7acqXfzYDb+2D0TEvAv65DafD38P55j7Mu1MBWxOrIQpO6iNl8ruhc7XCcLTXoU3FAqp7j7Ifm+iqi4uLw8GDB9G6detU6+X1rl270v2MBEdpy7dp0wYHDhxAfHz8I8skb9Oc/UrLVnBwcKoy0u3XpEmTDOsmkk+4uUETERHZEIMTUK490GUW8P5ZwNlT3WXXrKwv6t/fgmoRWzDTMA3HnAdhY/5ZeNl5J27dCsZXf57FF+tPo+3U7Zi65Rx+2XcFey+GWvto6DFZ7Rv+zp07SExMVK04KclrCVrSI+vTK5+QkKC2V7BgwQzLJG/TnP0mP6ZX5vLly+nWLTQ0FJ9//nmGrVHJpOVKlpRRLhER2RhdmnaHNuOA02uB07/D/m4QykXswGfYgTFOepxyr4cPdMNx5tZ9TN3yj+kjHSoXRLNyvjDo7dCqQn64OPCPbltg9atkZ5d6jiDpJku77r/Kp11vzjYzq0xy8NOhQwc1EH306NF4FOkKHDt27CPLEBGRjSlaO2lp9RkQckoFUDi9Drpbx1GpsCfW92iMpfuvYMc/d9Ao5CesuVMYG44bsf74TfXxwnmdMbhpKTXlS5UiHiiZz83aR0TZLXDy8fGBXq9/qHUpJCTkoZaeZDJAO73y0jXm7e39yDLJ2zRnv7INIWWkFetRdZPB4jKGys3NTd29ZzAYHnncH330Ed57771UQVfRokUf+RkiIrIR8sd1/opJS9MRQNhFICFWTSTcu24x9PZPAKYvwIsOknAzDw471sG62GpYH14Oo9ZEq01Ijk0pW7+Ut8pa3tDfhxMRZyNWG+Pk4OCg0gBs3rw51Xp53aBBg3Q/U79+/YfKb9q0CbVq1TIFLBmVSd6mOfuVNAUSPKUsI2OjJBVCyrpJ0CPjoGSbMmhdUhv8FxkrJeOgUi5ERJRDeZUEfMunTrpZqTvg5AG3xEg0itqCiYlf46jzIPyW52v09r0IowYs2XMZb/x0CH0X7kOHaYH49cBVrD16A7ciY6x5NGTtrjppeZF0ARL4SMAzd+5clRJA7pRLbp25fv06Fi9erF7LekkPIJ977bXX1EDwBQsWmO6WE0OGDEHjxo3V3W+dO3fGb7/9hi1btmDHjh1m71e64+SOunHjxqF06dJqkeeSr+nFF180tTRJ0BQVFYUff/xRBVHJ45Xy5cunWrWIiIhSyVcG6L4ASIwHruwBzm0Ezm6APuwCqsYdQtWu76O1vg7mbb8I96iriAu7jO3B/hi24p5p+pcX6/qhalEP5Hd3Qr2S3pwGJjcFTpI6QAZVf/bZZ7h58yYqVaqEP/74A8WKFVPvy7qUuZWkJUjef/fddzFz5kyVMmDatGl47rnnTGWkRWjp0qUYNWoUPvnkE5QqVUrlbZIUBObuVwwbNgzR0dF44403cPfuXfV5ablyd0/KDit35iXnhvL39091XHJXXvHixbPwzBERkU3TG4ASjZKWNl8Cd/5RAZTkjmri4IomZfIBf28Gtn+FOBdnHLavir26GlgWXhY/7DKaNlOpcB70rFVUDSxvWNoH+fP8d88H2XAeJ2IeJyIiykDAV8C+ucCDkFSrr9sXxWHHWvjsfheExP47rtbBXqeCKMkdVdjTGU3K+HJslJk4V50NYQJMIiLKkNEIBB8Dzm9JWq7uA7REwNkLoYNPYvb2SzhxIwLFIw9hV6gzrmj/3sBUNr87OlYpqCYmlnQHxbxdeaIzwMDJhjBwIiIis0WHAxe3AdF3gVqvJK3TNGhTKsIu8jpu2xfEKefqWHuvLP6KKYdwJA0vkbv6JIiStAfFvV3RqVohOBk4FjcZAycbwsCJiIieiARRS18Cru4BjAmm1RrscNWxNP5ybIGxIY1SfaSghxMalfZRwdOzVQuhdnGvXH0RIi2YcoVjnKyMgRMREWWK2HvA5V1JLVKySCJOUe8N7C3zgUq2aZcQA99T32NDVDmc1IpD+39WoqpF88LH1QGl87uj/zPF4ZvLBplHMnCyHQyciIgoS9wLBi4GJOWRKlglaZ0EVIs7q6fR9h4461IDK8JKYmdiBQRpkvzZTg0y98/npsZGda1eGM/XKgJH+5zdrRfJwMl2MHAiIqKn5speYOdUICgQiEvKDZUsyjEfpji/hXnBpVOtd7TXqS49SX0wsHEpVPfLC4M+aV1OwcDJhjBwIiKip04ScF4/BFzcClzakXS3XmIstIHbcSjOD+FRcdBOrkHiybXYFlcWe4z/tkgJmZhYxkZJ+gNHgx6lfd3U9DC2ioGTDWHgREREVhcfA1zbDxR7BtD9fza21YOAo//OzHHP4INtsWWwK7H8Q4GUu6O9ymheJr87PJwNaFwmn+rysxUMnGwIAyciIsqWru4Hzm9OapGSoCoxLtXbQ/xWYfdNDfGJRuij7iAM7jD+f7C53LXXpmIBFTzVK+mFpmV8s/XUMAycbAgDJyIiyvbio5OCp0s7kpb4KOD1beoto1FD+KzWcA49gX8cymNXfGkExJTCYaM/YuCoyhT3dkERTxd4uzmgVx0/1C3hpeaFzS4YONkQBk5ERGSTGc11un+ff1P2oalhEqHHDecyWBtVGV/Fdkn1nnTtSQuU5JLqU6+YmmPP09VBdfNZAwMnG8LAiYiIbJ4xMSlv1JU9wJXdwOXdwL0b6q14/zb4o9IUJCRqOHjlLvwPj8fpxMKqReqCVsiUS0oGnHeqWhgNSnnDXm+HBqV8kM89qcUqqzFwsiEMnIiIKMfRNCDialIg5eIN+LdIWh9+FZhayVQsWueKw4mlcEQrjX0JJXHYWBoRcFPvOeh1aFHeV7VCFfVyUfmkfN2zJjEnAycbwsCJiIhyjcgbwN7vgGsHgBuHksZKpRDg2R1zXV9D6P04BAWHwt/uBs5oRVW3n73ODgXzOiHsfhwOfdoqU5NyWvJdbLtJF4iIiMi25CkEtBqb9DwxIal7Twad/39p0qwDmlSqp94+u3cDym4YiXidE87p/REYXQKHw/2x3VgZIZGxqhXKGhg4ERER0dOnt0+aCkaW2q/+28X3f2VdHgCOHjDERqCi8QQq2p9Q66/03afSHVgLAyciIiLKHuxSpCio3B2o2A24c+7fVqmwi/ArXiZ1uaeMgRMRERFlTzod4FsuaanRB9mB7eRDJyIiIrIyBk5EREREZmLgRERERGQmBk5EREREZmLgRERERGQmBk5EREREZmLgRERERGQmBk5EREREZmLgRERERMTAiYiIiChzscWJiIiIyEwMnIiIiIjMxEl+rUzTNPUYGRlp7aoQERHlSpH//w5O/k5+FAZOVnbv3j31WLRoUWtXhYiICLn9O9nDw+ORZew0c8IryjJGoxE3btyAu7s77OzsMjV6lmDs6tWryJMnD3ISHpvt4TWzTbxutonXzXISCknQVKhQIeh0jx7FxBYnK5MLVKRIkSzbvgRNOS1wSsZjsz28ZraJ18028bpZ5r9ampJxcDgRERGRmRg4EREREZmJgVMO5ejoiNGjR6vHnIbHZnt4zWwTr5tt4nXLWhwcTkRERGQmtjgRERERmYmBExEREZGZGDgRERERmYmBUw40a9YslChRAk5OTqhZsyYCAwNha8aMGaMSgqZcChQokCpZmZSRZGXOzs5o2rQpTp48iexo+/btePbZZ1Vd5TjWrFmT6n1zjiU2NhZvv/02fHx84Orqik6dOuHatWvI7sfWr1+/h65jvXr1sv2xjR8/HrVr11aJaX19fdGlSxecPXs2R1w3c47NVq/b7NmzUaVKFVP+ovr162PDhg02f83MOTZbvWbp/XxK3YcOHZptrxsDpxxm2bJl6gfu448/xuHDh9GoUSO0a9cOV65cga2pWLEibt68aVqOHz9uem/SpEmYPHkyZsyYgf3796ugqlWrVqYpbLKTBw8eoGrVqqqu6THnWOSarl69GkuXLsWOHTtw//59dOzYEYmJicjOxybatm2b6jr+8ccfqd7PjscWEBCAN998E3v27MHmzZuRkJCA1q1bq+O19etmzrHZ6nWTZMITJkzAgQMH1NK8eXN07tzZ9CVrq9fMnGOz1WuWklyTuXPnqgAxpWx33WTKFco56tSpow0aNCjVunLlymkjRozQbMno0aO1qlWrpvue0WjUChQooE2YMMG0LiYmRvPw8NDmzJmjZWfyX2716tUWHUt4eLhmMBi0pUuXmspcv35d0+l02saNG7Xsemyib9++WufOnTP8jK0cW0hIiDq+gICAHHfd0h5bTrpuwtPTU5s/f36OumZpjy0nXLN79+5ppUuX1jZv3qw1adJEGzJkiFqfHa8bW5xykLi4OBw8eFD99ZiSvN61axdszT///KOaZqXb8YUXXsDFixfV+qCgIAQHB6c6Tslb0qRJE5s7TnOORa5pfHx8qjJyXipVqmQTx7tt2zbVJVSmTBm89tprCAkJMb1nK8cWERGhHr28vHLcdUt7bDnluklLg7Q+SEuadGvlpGuW9thywjV788030aFDB7Rs2TLV+ux43ThXXQ5y584d9R8qf/78qdbLa/nBsyV169bF4sWL1S+AW7du4YsvvkCDBg1Us3TysaR3nJcvX4YtMedYpIyDgwM8PT1t7rpKN/Hzzz+PYsWKqV+An3zyiepikF908svPFo5NGtPee+89NGzYUP0izknXLb1js/XrJl36EkzExMTAzc1Ndd9UqFDB9AVqy9cso2Oz9Wu2dOlSHDp0SHXDpZUd/68xcMqBZGBd2l+Oaddld/JLIFnlypXVL4tSpUph0aJFpgGPOeE4kz3OsdjC8fbs2dP0XL6Ya9WqpX6xr1+/Ht26dbOJY3vrrbdw7NgxNW4ip123jI7Nlq9b2bJlceTIEYSHh2PlypXo27evGteVE65ZRscmwZOtXrOrV69iyJAh2LRpk7qhKSPZ6bqxqy4HkbsJ9Hr9QxG2NNemjdZtjdwlIQGUdN8l312XE47TnGORMtINe/fu3QzL2IqCBQuqX+ZyHW3h2OQunbVr12Lr1q1qcG5Oum4ZHZutXzdpefD391eBg9yhJTcvfPvttznimmV0bLZ8zQ4ePKjqIHeA29vbq0WCwWnTpqnnyXXLTteNgVMOIv+p5IdP7pRJSV5LN5ctk1tNT58+rX4ZyJgn+Y+S8jjlP438Z7O14zTnWOSaGgyGVGXkjpkTJ07Y3PGGhoaqvzDlOmbnY5O/VKU1ZtWqVfj777/Vdcop1+2/js2Wr1tGxyu/P2z5mv3XsdnyNWvRooXqgpSWtORFAsPevXur5yVLlsx+1y3Th5uTVcldBXJ3wYIFC7RTp05pQ4cO1VxdXbVLly7Z1JV5//33tW3btmkXL17U9uzZo3Xs2FFzd3c3HYfcYSF3VaxatUo7fvy41qtXL61gwYJaZGSklt3I3SKHDx9Wi/yXmzx5snp++fJls49F7pQsUqSItmXLFu3QoUNa8+bN1V2HCQkJ2fbY5D25jrt27dKCgoK0rVu3avXr19cKFy6c7Y9t8ODB6prIz+DNmzdNS1RUlKmMrV63/zo2W75uH330kbZ9+3ZV72PHjmkjR45Ud1Zt2rTJpq/Zfx2bLV+z9KS8qy47XjcGTjnQzJkztWLFimkODg5ajRo1Ut1mbCt69uyp/mNIEFioUCGtW7du2smTJ03vyy2qkrJAblN1dHTUGjdurP5DZUfyS0yCirSL3D5s7rFER0drb731lubl5aU5OzurQPLKlStadj42+SJu3bq1li9fPnUd/fz81Pq09c6Ox5beMcny/fffm8rY6nX7r2Oz5evWv39/0+8+qX+LFi1MQZMtX7P/OjZbvmbmBE7Z7brZyT+Z345FRERElPNwjBMRERGRmRg4EREREZmJgRMRERGRmRg4EREREZmJgRMRERGRmRg4EREREZmJgRMRERGRmRg4EREREZmJgRMRPRXbtm1TM5XLzO7WIPOylStXDkajEblR06ZNMXTo0Mf+vMwnJpMBP3jwIFPrRWRrGDgR0VP5kpbJNmXiTQ8PD6uc8WHDhuHjjz+GTpf0a++HH35QgVz58uUfKvvrr7+q94oXL46cQib1/fzzz02v5dimTp1q9ucrV66MOnXqYMqUKVlUQyLbwMCJiJ4KBwcHNcu5BCRP265du/DPP//g+eefT7Xe1dUVISEh2L17d6r1CxcuhJ+fH3ISLy8vuLu7P9E2XnnlFcyePRuJiYmZVi8iW8PAiYgyVb9+/RAQEIBvv/1WBUmyXLp06aGuOmnxyZs3L9atW4eyZcvCxcUF3bt3V11BixYtUi0inp6eePvtt1N9UcfFxanWo8KFC6vAp27dumrbj7J06VK0bt0aTk5Oqdbb29vjxRdfVIFSsmvXrqntyfq0fv/9d9SsWVNtp2TJkhg7diwSEhJM70dEROD111+Hr68v8uTJg+bNm+Po0aOm98eMGYNq1arhu+++Q9GiRdUxSzD3qO7L5POU0po1a1IFoMnbXbJkiTpv0qr3wgsv4N69e+m2Asrzy5cv49133zVdIyHrnn32WXXe5dxWrFgRf/zxh2kbbdq0QWhoqLq+RLkVAyciylQSMNWvXx+vvfaa6pqTRYKE9ERFRWHatGkqsNm4caMKWLp166a+rGWRQGDu3LlYsWJFqlaPnTt3qs8cO3ZMBR5t27ZVLUoZ2b59O2rVqpXue6+++iqWLVum6pIcqMj28ufPn6rcn3/+iZdeegnvvPMOTp06pYIfKfvll1+q92W+9A4dOiA4OFjV/eDBg6hRowZatGiBsLAw03bOnz+vugIlCJNjPnLkCN588008qQsXLqiASgJRWSS4mTBhQobddjJe6bPPPjNdIyH1iI2NVedLxjRNnDgRbm5uqVoNq1atisDAwCeuL5GtYuBERJlKWjvkC1ZaU6RrTha9Xp9u2fj4eNX1U716dTRu3Fi1OO3YsQMLFixAhQoV0LFjRzRr1gxbt241BQe//PILli9fjkaNGqFUqVL44IMP0LBhQ3z//fcZ1klavAoVKpTue9JSI9uR4EyCHwmG+vfv/1A5CZBGjBiBvn37qtamVq1aqTFDEkAJqaMEG1I3CdJKly6Nr7/+WrUWpQz8YmJiVIua7FeOefr06SoIlIDrScigd6l7pUqV1Lnp06cP/vrrrwy77eSaSNdd8jUSV65cwTPPPKPGM8kxyvmXOqYkLX1yPolyK3trV4CIci8JriRoSSatPNLVlLKVQ9bJOCRx6NAhFdyUKVMm1XaklcTb2zvD/URHRz/UTZeSBEoSeMm4pvv376N9+/aYMWNGqjLSgrR//35TC5OQLkQJhKS1St6Xz6ath+xbAr5ksg9p7UkmrXMS9Jw9e9YUwDwOOW8pxzAVLFjQdN7MJa1pgwcPxqZNm9CyZUs899xzqFKlSqoyzs7OptY5otyIgRMRWY3BYEj1WsbapLcuOYWAPEpLiQQpaVuxUgZbafn4+ODu3bsZvt+7d281bkrGCr388stq7FNasm8Z0yRdiWlJUCbvS7CS3nirtGOU0h5fyse05C5ACRbTttSl9ajzZq4BAwaocUzr169XwdP48ePxzTffqHFmyaTbMWWwS5TbMHAiokwnXXVZceeVdOnJdqUlRbqjLPmcjEvKiHRdderUSY09mjNnTrplZLyStAr5+/tn+L50t0nQ9ag0BtIdduPGDVPXodzRJ8FR2la0ZPny5VODvGXQvAzYFjIuKquukYxHGzRokFo++ugjzJs3L1XgdOLECdWlSpRbcYwTEWU6CRz27t2rxsLcuXMn05JOSnAhrUPSKiQDnIOCglT3mQxiTnn3V1rSiiJjpx5FxgdJXSVJZno+/fRTLF68WLVKnTx5EqdPn1aDykeNGqXel64t6Xbr0qWLGkguxy5pEOT9AwcOpGqdknFScredDLKW7rEePXpk2E0ndw1Kl+bIkSPVwPKff/5Z1TUzrpEMAr9+/bo6biF33Und5bxKt6gkDU2Z50qOScrLsRLlVgyciCjTyYBt6UqTAd7SYiKtLJlFxiJJ4PT++++rNAbSUiRBWkZ37gm5G05anKTFKCMydudR46Qk+JK71TZv3ozatWujXr16mDx5MooVK2bqGpPgTQZTy5gpCfIkJYAEGynv0JMWK+nuk3FUkiJBBnPPmjXrka1hP/74o9q2DNqWwfESvD0puaNO6ibdbnKNhLRAyZ11EizJnYVyflPWTfYtdU4+ZqLcyE5L23lORJQDyRgmybOUfBecNUjAIykDMqOr7WmTAfhyp6AET3LnHVFuxRYnIsoVZLoVaSlh1uvHI8kx5RwyaKLcjoPDiSjX5JeScUL0eKTrMaMB7ES5CbvqiIiIiMzErjoiIiIiMzFwIiIiIjITAyciIiIiMzFwIiIiIjITAyciIiIiMzFwIiIiIjITAyciIiIiMzFwIiIiIjITAyciIiIimOd/FO7P5GY0L+cAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "time_fs = time_au * 0.02418884254\n", "time_meep = time_fs / 0.1\n", "initial = population[0]\n", "dipole_moment = 0.1 # meep units of mu12\n", "frequency = 1.0 # meep units of omega\n", "\n", "# analytical golden-rule decay rate\n", "gamma = dipole_moment**2 * frequency**2 / 2.0\n", "\n", "# simple exponential decay reference\n", "reference = initial * np.exp(-time_meep * gamma)\n", "\n", "# below is a more accurate reference which works under any initial population\n", "#reference = np.exp(-time_meep * gamma) / (\n", "# np.exp(-time_meep * gamma) + (1.0 - initial) / initial\n", "#)\n", "\n", "std_rel = np.std(population - reference) / initial\n", "max_rel = np.max(np.abs(population - reference)) / initial\n", "print(f\"std_dev={std_rel:.3e}, max_abs_diff={max_rel:.3e}\")\n", "\n", "import matplotlib.pyplot as plt\n", "\n", "plt.figure(figsize=(6, 4))\n", "plt.plot(time_meep, population, label=\"Simulation\")\n", "plt.plot(time_meep, reference, label=\"Analytical\", linestyle=\"--\")\n", "plt.xlabel(\"time (Meep units)\")\n", "plt.ylabel(\"excited-state population\")\n", "plt.legend()\n", "plt.tight_layout()\n", "plt.show()\n" ] } ], "metadata": { "kernelspec": { "display_name": "mxl", "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.13.9" }, "title": "Socket TLS Workflow" }, "nbformat": 4, "nbformat_minor": 5 }