diff --git a/.github/workflows/config/spelling_allowlist.txt b/.github/workflows/config/spelling_allowlist.txt index b300464829..ecff28b403 100644 --- a/.github/workflows/config/spelling_allowlist.txt +++ b/.github/workflows/config/spelling_allowlist.txt @@ -1,4 +1,5 @@ ABI +AFQMC API APIs AST diff --git a/.github/workflows/docker_images.yml b/.github/workflows/docker_images.yml index 797ddadc50..30cb6d25d2 100644 --- a/.github/workflows/docker_images.yml +++ b/.github/workflows/docker_images.yml @@ -778,7 +778,7 @@ jobs: docker cp docs/notebook_validation.py cuda-quantum:"/home/cudaq/notebook_validation.py" # In containers without GPU support, UCX does not work properly since it is configured to work with GPU-support. # Hence, don't enforce UCX when running these tests. - docker exec cuda-quantum bash -c "python3 -m pip install pandas scipy pandas seaborn" + docker exec cuda-quantum bash -c "python3 -m pip install pandas scipy pandas seaborn 'h5py<3.11'" (docker exec cuda-quantum bash -c "unset OMPI_MCA_pml && set -o pipefail && bash validate_container.sh | tee /tmp/validation.out") && passed=true || passed=false docker cp cuda-quantum:"/tmp/validation.out" /tmp/validation.out docker stop cuda-quantum diff --git a/.github/workflows/python_wheels.yml b/.github/workflows/python_wheels.yml index f3f95a457a..c60ba3baaf 100644 --- a/.github/workflows/python_wheels.yml +++ b/.github/workflows/python_wheels.yml @@ -316,7 +316,7 @@ jobs: docker run --rm -dit --name wheel-validation-examples wheel_validation:local status_sum=0 - for ex in `find docs/sphinx/examples/python -name '*.py' -not -path '*/providers/*' -not -path '*/divisive_clustering_src/*'`; do + for ex in `find docs/sphinx/examples/python -name '*.py' -not -path '*/providers/*' -not -path '*/divisive_clustering_src/*' -not -path '*/utils_ipie.py' -not -path '*/vqe_cudaq_qnp.py'`; do file="${ex#docs/sphinx/examples/python/}" echo "__Example ${file}:__" >> /tmp/validation.out (docker exec wheel-validation-examples bash -c "python${{ inputs.python_version }} /tmp/examples/$file" >> /tmp/validation.out) && success=true || success=false diff --git a/docker/test/installer/linux.Dockerfile b/docker/test/installer/linux.Dockerfile index 8b9c196e8e..288974fa0e 100644 --- a/docker/test/installer/linux.Dockerfile +++ b/docker/test/installer/linux.Dockerfile @@ -10,7 +10,7 @@ ARG base_image=redhat/ubi8:8.0 ARG base_image_mpibuild=amd64/almalinux:8 # [OpenMPI Installation] -FROM ${base_image_mpibuild} as mpibuild +FROM ${base_image_mpibuild} AS mpibuild ARG base_image_mpibuild SHELL ["/bin/bash", "-c"] ARG DEBIAN_FRONTEND=noninteractive diff --git a/docs/sphinx/examples/python/tutorials/afqmc.ipynb b/docs/sphinx/examples/python/tutorials/afqmc.ipynb new file mode 100644 index 0000000000..5b020c74e5 --- /dev/null +++ b/docs/sphinx/examples/python/tutorials/afqmc.ipynb @@ -0,0 +1,531 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Quantum Enhanced Auxiliary Field Quantum Monte Carlo\n", + "\n", + "This work was done in collaboration with the Next Generation Computing team at [BASF](https://www.basf.com/global/en.html).\n", + "\n", + "In this tutorial we implement a quantum-classical hybrid workflow for computing the ground state energies of a strongly interacting molecular system. The algorithm consists of two parts:\n", + "\n", + "\n", + "1. A variational quantum eigensolver that uses the quantum-number-preserving ansatz proposed by [Anselmetti et al. (2021)](https://doi.org/10.1088/1367-2630/ac2cb3) to generate a quantum trial wave function $|\\Psi_T\\rangle$ using CUDA Quantum.\n", + "\n", + "2. An Auxiliary-Field Quantum Monte Carlo simulation that realizes a classical imaginary time evolution and collects the ground state energy estimates.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Package installs\n", + "!pip install pyscf==2.6.2\n", + "!pip install openfermion==1.6.1\n", + "!pip install ipie==0.7.1" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# Relevant imports\n", + "\n", + "import cudaq\n", + "\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from pyscf import gto, scf, ao2mo, mcscf\n", + "\n", + "from src.vqe_cudaq_qnp import VQE, get_cudaq_hamiltonian\n", + "from src.utils_ipie import get_coeff_wf, gen_ipie_input_from_pyscf_chk\n", + "\n", + "from ipie.hamiltonians.generic import Generic as HamGeneric\n", + "from ipie.qmc.afqmc import AFQMC\n", + "from ipie.systems.generic import Generic\n", + "from ipie.trial_wavefunction.particle_hole import ParticleHole\n", + "from ipie.analysis.extraction import extract_observable\n", + "\n", + "from ipie.config import config\n", + "\n", + "cudaq.set_target(\"nvidia\")\n", + "\n", + "# Ipie has recently added GPU support however this remains a bit tricky to use as it requires manual installation of several packages.\n", + "# Once this is streamlined, we can set the GPU option to True in the tutorial.\n", + "config.update_option(\"use_gpu\", False)\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We start by defining the structure of the molecule, the basis set, and its spin. We build the molecule object with PySCF and run a preliminary Hartree-Fock computation. Here we choose as an example a [chelating agent](https://doi.org/10.1021/acs.jctc.3c01375) representing a relevant class of substances industrially produced at large scales. Their use ranges, among the others, from water softeners in cleaning applications, modulators of redox behaviour in oxidative bleaching, scale suppressants, soil remediation and ligands for catalysts. In particular we focus here in a Fe(III)-NTA complex whose structure is given in the file imported below.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# Define the molecular structure and the basis set for the Fenta molecule.\n", + "\n", + "atom = \"src/geo_fenta.xyz\"\n", + "basis = \"cc-pVTZ\"\n", + "spin = 1\n", + "num_active_orbitals = 5\n", + "num_active_electrons = 5\n", + "\n", + "# You can swap to O3 which is a smaller system and takes less computational resources and time to run.\n", + "\n", + "atom = \"src/geo_o3.xyz\"\n", + "basis = \"cc-pVTZ\"\n", + "spin = 0\n", + "num_active_orbitals = 9\n", + "num_active_electrons = 12\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "-224.34048064812245" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# PYSCF helps to build the molecule and run Hartree-Fock.\n", + "\n", + "# Define the molecule.\n", + "molecule = gto.M(atom=atom, spin=spin, basis=basis, verbose=0)\n", + "\n", + "# Restriced open shell HF.\n", + "hartee_fock = scf.ROHF(molecule)\n", + "hartee_fock.chkfile = \"src/output.chk\"\n", + "\n", + "# Run Hartree-Fock.\n", + "hartee_fock.kernel()\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Hamiltonian preparation for VQE\n", + "\n", + "Since this molecule contains of around 600 orbitals (which would correspond to 1200 qubits) and 143 total electrons, it is impossible to perform a full VQE with full statevector simulation. Therefore, we need to identify an active space with fewer orbitals and electrons that contribute to the strongly interacting part of the whole molecule. We then run a post Hartree-Fock computation with the PySCF's built-in CASCI method in order to obtain the one-body ($t_{pq}$) and two-body \n", + "($V_{prqs}$) integrals that define the molecular Hamiltonian in the active space:\n", + "\n", + "$$ H= \\sum_{pq}t_{pq}\\hat{a}_{p}^\\dagger \\hat {a}_{q}+\\sum_{pqrs} V_{prqs}\\hat a_{p}^\\dagger \\hat a_{q}^\\dagger \\hat a_{s}\\hat a_{r} \\tag{1}$$\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "from openfermion.transforms import jordan_wigner\n", + "from openfermion import generate_hamiltonian\n", + "\n", + "# Run a CASCI simulation for computing the Hamiltonian in the active space.\n", + "casci = mcscf.CASCI(hartee_fock, num_active_orbitals, num_active_electrons)\n", + "casci.fix_spin_(ss=(molecule.spin / 2 * (molecule.spin / 2 + 1)))\n", + "\n", + "\n", + "# Executes the kernel to compute the hamiltonian in the active space.\n", + "casci.kernel()\n", + "\n", + "# Compute the one-body (h1) and two-body integrals (tbi) as shown in equation 1.\n", + "h1, energy_core = casci.get_h1eff()\n", + "\n", + "h2 = casci.get_h2eff()\n", + "h2_no_symmetry = ao2mo.restore('1', h2, num_active_orbitals)\n", + "\n", + "# V_pqrs terms in H.\n", + "tbi = np.asarray(h2_no_symmetry.transpose(0, 2, 3, 1), order='C')\n", + "\n", + "# Compute the hamiltonian and convert it to a CUDA-Q operator.\n", + "mol_ham = generate_hamiltonian(h1, tbi, energy_core.item())\n", + "jw_hamiltonian = jordan_wigner(mol_ham)\n", + "hamiltonian, constant_term = get_cudaq_hamiltonian(jw_hamiltonian)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Run VQE with CUDA-Q\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can now execute the VQE algorithm using the quantum number preserving ansatz. At the end of the VQE, we store the final statevector that will be used in the classical AFQMC computation as an initial guess.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "# Using cudaq optimizer\n", + "# Num Params: 16\n", + "# Qubits: 18\n", + "# N_layers: 1\n", + "# Energy after the VQE: -224.3926998830352\n" + ] + } + ], + "source": [ + "# Define some options for the VQE.\n", + "options = {\n", + " 'n_vqe_layers': 1,\n", + " 'maxiter': 100,\n", + " 'energy_core': constant_term,\n", + " 'return_final_state_vec': True\n", + "}\n", + "\n", + "n_qubits = 2 * num_active_orbitals\n", + "\n", + "vqe = VQE(n_qubits=n_qubits,\n", + " num_active_electrons=num_active_electrons,\n", + " spin=spin,\n", + " options=options\n", + ")\n", + "\n", + "results = vqe.execute(hamiltonian)\n", + "\n", + "# Best energy from VQE.\n", + "optimized_energy = results['energy_optimized']\n", + "\n", + "# Final state vector.\n", + "final_state_vector = results[\"state_vec\"]\n", + "\n", + "# Energies during the optimization loop.\n", + "vqe_energies = results[\"callback_energies\"]\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Auxiliary Field Quantum Monte Carlo (AFQMC)\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "AFQMC is a numerical method for computing relevant properties of strongly interacting molecules. AFQMC is a type of Quantum Monte Carlo method that combines the use of random walks with an auxiliary field to simulate the imaginary-time evolution of a quantum system and drive it to the lowest energy state. This method can provide accurate results for ground-state properties of a wide range of physical systems, including atoms, molecules, and solids. Here we summarize the main features of AFQMC while a detailed introduction to can be found [here](https://www.cond-mat.de/events/correl13/manuscripts/zhang.pdf).\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We consider the electronic Hamiltonian in the second quantization\n", + "\\begin{equation}\n", + "H = {H}_1 + {H}_2 \n", + "=\\sum_{pq} h_{pq} {a}_{p}^{\\dagger} {a}_{q} + \\frac{1}{2} \\sum_{pqrs} v_{pqrs}{a}_{p}^{\\dagger} {a}_r {a}^{\\dagger}_{q} {a}_s \\tag{2}\n", + "\\end{equation}\n", + "where ${a}_{p}^{\\dagger}$ and ${a}_{q}$ are fermionic creation and annihilation operators of orbitals $p$ and $q$, respectively. The terms $h_{pq} $ and \n", + "$v_{pqrs}$ are the matrix elements of the one-body, $H_1$, and two-body, $H_2$, interactions of $H$, respectively. Here, we omit the spin indices for simplicity.\n", + "\n", + "AFQMC realizes an imaginary time propagation of an initial state (chosen as a Slater determinant) $\\ket{\\Psi_{I}}$ towards the ground state $\\ket{\\Psi_0}$ of a given hamiltonian, $H$, with\n", + "\\begin{equation}\n", + "\\ket{\\Psi_0} \\sim\\lim_{n \\to \\infty} \\left[ e^{-\\Delta\\tau H } \\right]^{n} \\ket{\\Psi_{I}}\n", + "\\tag{3}\n", + "\\end{equation} \n", + "where $\\Delta\\tau$ is the imaginary time step.\n", + "\n", + "AFQMC relies on decomposing the two-body interactions $H_2$ in terms of sum of squares of one-body operators ${v}_\\gamma$ such that the Hamiltonian ${H}$ becomes\n", + "\\begin{equation}\n", + "H = v_0 - \\frac{1}{2}\\sum_{\\gamma=1}^{N_\\gamma} {v}_\\gamma^2\n", + "\\tag{4}\n", + "\\end{equation}\n", + "with ${v}_0 = {H}_1 $ and $\n", + "{v}_\\gamma = i \\sum_{pq} L^{\\gamma}_{pq} {a}_{p}^{\\dagger}{a}_{q}.\n", + "$\n", + "The $N_\\gamma$ matrices $L^{\\gamma}_{pq}$ are called Cholesky vectors as they are obtained via a Cholesky decomposition of the two-body matrix elements \n", + "$v_{pqrs}$ via $v_{pqrs} = \\sum_{\\gamma=1}^{N_\\gamma} L^{\\gamma}_{pr} L^{\\gamma}_{qs}$.\n", + "\n", + "The imaginary time propagation evolves an ensemble of walkers $\\{\\phi^{(n)}\\}$ (that are Slater determinants) and allows one to access observables of the system. For example, the local energy\n", + "\\begin{equation}\n", + "\\mathcal{E}_{\\text{loc}}(\\phi^{(n)}) = \\frac{\\bra{\\Psi_\\mathrm{T}}H\\ket{\\phi^{(n)}}}{\\braket{\\Psi_\\mathrm{T}| \\phi^{(n)}}}\n", + "\\tag{5}\n", + "\\end{equation}\n", + "defined as the mixed expectation value of the Hamiltonian with the trial wave function $\\ket{\\Psi_\\mathrm{T}}$.\n", + "\n", + "\n", + "The trial wavefunction can be in general a single or a multi-Slater determinant coming from VQE for example. This might help in achieving more accurate ground state energy estimates.\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "The implementation of AFQMC we use here is from [ipie](https://github.com/JoonhoLee-Group/ipie) that supports both CPUs and GPUs and requires the following steps:\n", + "\n", + "\n", + "1. Preparation of the molecular Hamiltonian by performing the Cholesky decomposition\n", + "\n", + "2. Preparation of the trial state from the VQE wavefunction\n", + "\n", + "3. Executing AFQMC\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Preparation of the molecular Hamiltonian\n" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "# Number of electrons in simulation: (12, 12)\n" + ] + } + ], + "source": [ + "# AFQMC.\n", + "\n", + "# Generate the input Hamiltonian for ipie from the checkpoint file from pyscf.\n", + "ipie_hamiltonian = gen_ipie_input_from_pyscf_chk(hartee_fock.chkfile,\n", + " mcscf=True,\n", + " chol_cut=1e-5)\n", + "\n", + "\n", + "h1e, cholesky_vectors, e0 = ipie_hamiltonian\n", + "\n", + "\n", + "num_basis = cholesky_vectors.shape[1]\n", + "num_chol = cholesky_vectors.shape[0]\n", + "\n", + "\n", + "system = Generic(nelec=molecule.nelec)\n", + "\n", + "\n", + "afqmc_hamiltonian = HamGeneric(\n", + " np.array([h1e, h1e]),\n", + " cholesky_vectors.transpose((1, 2, 0)).reshape((num_basis * num_basis, num_chol)),\n", + " e0\n", + ")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Preparation of the trial wave function\n" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "# Build the trial wavefunction from the state vector computed via VQE.\n", + "wavefunction = get_coeff_wf(final_state_vector,\n", + " n_active_elec=num_active_electrons,\n", + " spin=spin\n", + ")\n", + "\n", + "\n", + "trial = ParticleHole(wavefunction,\n", + " molecule.nelec,\n", + " num_basis,\n", + " num_dets_for_props=len(wavefunction[0]),\n", + " verbose=False\n", + ")\n", + "\n", + "\n", + "trial.compute_trial_energy = True\n", + "trial.build()\n", + "trial.half_rotate(afqmc_hamiltonian)\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Setup of the AFQMC parameters\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here we can choose the input options like the timestep $\\Delta\\tau$, the total number of walkers `num_walkers` and the total number of AFQMC iterations `num_blocks`.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "# random seed is 96264512\n", + " Block Weight WeightFactor HybridEnergy ENumer EDenom ETotal E1Body E2Body\n", + " 0 1.0000000000000000e+02 1.0000000000000000e+02 0.0000000000000000e+00 -2.2440474833475844e+04 1.0000000000000000e+02 -2.2440474833475844e+02 -3.7639365216971441e+02 1.5198890383495583e+02\n", + " 1 4.2243136724295840e+02 1.4111352538015647e+03 -1.1706587453647199e+02 -2.2475815411995434e+04 1.0000000000000000e+02 -2.2475815411995436e+02 -3.7647553637385266e+02 1.5171738225389828e+02\n", + " 2 1.0031170675347687e+02 3.8278928957673060e+02 -1.1737451166143732e+02 -2.2491281828803054e+04 1.0000000000000000e+02 -2.2491281828803054e+02 -3.7652044141962244e+02 1.5160762313159196e+02\n", + " 3 9.9897251153181060e+01 1.0007109151776878e+02 -1.1731110994412086e+02 -2.2497799156719888e+04 1.0000000000000000e+02 -2.2497799156719887e+02 -3.7661015848515297e+02 1.5163216691795418e+02\n", + " 4 1.0008903318642547e+02 1.0004229364465529e+02 -1.1743127719028180e+02 -2.2497508617456722e+04 1.0000000000000000e+02 -2.2497508617456722e+02 -3.7677595649331067e+02 1.5180087031874356e+02\n", + " 5 9.9996605525745039e+01 1.0010272949474228e+02 -1.1746990296069042e+02 -2.2504582061507692e+04 1.0000000000000001e+02 -2.2504582061507688e+02 -3.7665153235346526e+02 1.5160571173838829e+02\n", + " 6 1.0011254258887075e+02 1.0017483684345871e+02 -1.1763876429459457e+02 -2.2514581314082043e+04 1.0000000000000000e+02 -2.2514581314082042e+02 -3.7663426702613344e+02 1.5148845388531296e+02\n", + " 7 9.9928457034724943e+01 9.9908499343007676e+01 -1.1758362687039269e+02 -2.2517447174632791e+04 1.0000000000000000e+02 -2.2517447174632792e+02 -3.7662538302479038e+02 1.5145091127846254e+02\n", + " 8 9.9903126164668137e+01 9.9911410200683974e+01 -1.1754847066192167e+02 -2.2518432403785569e+04 1.0000000000000000e+02 -2.2518432403785570e+02 -3.7677549496249651e+02 1.5159117092464081e+02\n", + " 9 1.0011461067810892e+02 1.0011099130959934e+02 -1.1773137871282010e+02 -2.2513737823028405e+04 1.0000000000000000e+02 -2.2513737823028404e+02 -3.7679200521685516e+02 1.5165462698657115e+02\n", + " 10 9.9639714467961483e+01 9.9243154673312134e+01 -1.1743325794078407e+02 -2.2518897510231673e+04 1.0000000000000000e+02 -2.2518897510231673e+02 -3.7689860950318899e+02 1.5170963440087232e+02\n" + ] + } + ], + "source": [ + "# Setup the AFQMC parameters.\n", + "afqmc_msd = AFQMC.build(molecule.nelec,\n", + " afqmc_hamiltonian,\n", + " trial,\n", + " num_walkers=100,\n", + " num_steps_per_block=25,\n", + " num_blocks=10,\n", + " timestep=0.005,\n", + " stabilize_freq=5,\n", + " seed=96264512,\n", + " pop_control_freq=5,\n", + " verbose=False\n", + ")\n", + "\n", + "\n", + "# Run the AFQMC.\n", + "afqmc_msd.run(estimator_filename = 'src/estimates.0.h5')\n", + "afqmc_msd.finalise(verbose=False)\n", + "\n", + "\n", + "# Extract the energies.\n", + "qmc_data = extract_observable(afqmc_msd.estimators.filename, \"energy\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlQAAAGwCAYAAABvpfsgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB3XUlEQVR4nO3deXgT1f4G8HeSNume7hstpZRCWQXBpSgCiheUIvzkKm7si9zLVVD0CopwxaWigAsqIrJ5RbmgKIoLIFIEQTZB1hYKlKU7pW3SvU3m90c606RJ2nQjKbyf58nTZmYyc2ZA8nrOd84IoiiKICIiIqJGUzi6AUREREStHQMVERERURMxUBERERE1EQMVERERURMxUBERERE1EQMVERERURMxUBERERE1kYujG3AjMBgMyMjIgLe3NwRBcHRziIiIyA6iKEKn0yE8PBwKRd19UAxU10BGRgYiIyMd3QwiIiJqhEuXLiEiIqLObRiorgFvb28Axj8QHx8fB7eGiIiI7KHVahEZGSl/j9eFgeoakIb5fHx8GKiIiIhaGXvKdViUTkRERNREDFRERERETcRARURERNRErKEiIiJyEL1ej8rKSkc344amUqnqnRLBHgxURERE15goisjKykJBQYGjm3LDUygUiI6OhkqlatJ+GKiIiIiuMSlMBQcHw8PDg5M+O4g08XZmZibatm3bpD8HBioiIqJrSK/Xy2EqICDA0c254QUFBSEjIwNVVVVwdXVt9H5YlE5ERHQNSTVTHh4eDm4JAZCH+vR6fZP2w0BFRETkABzmcw7N9efAQEVERETURAxURERERE3EQEVERETURAxU5JRKK/QQRdHRzSAiomrDhg3DkCFDrK7btWsXBEHA0aNHAQBr1qzBLbfcAg8PD3h7e6N///7YvHmz2WeSkpIgCILVV1ZWVoufT3NjoCKnk5pThJ7zt2LuphOObgoREVWbOHEitm3bhsuXL1usW7VqFfr06YMePXrgueeew5NPPolRo0bh6NGj2L9/P+68804MHz4cH3zwgcVnU1JSkJmZafYKDg6+FqfUrDgPFTmdExmFKK8y4JdT2Xh1RDdHN4eIqMWJoojSyqbdtt9Y7q5Ku+50S0hIQFBQEFavXo05c+bIy4uKirBhwwa8/fbb+OOPP7Bo0SK8//77eOqpp+RtXn/9dZSVleHZZ5/F8OHDERkZKa8LDg6Gr69vs56TIzBQkdMpLjf+o5JZWIbCkkpoPBo/0RoRUWtQWqlHl7lbHHLsk/MHw0NVfxxwcXHBmDFjsHr1arz00ktyCNuwYQP0ej0effRRzJ07F15eXnjyySctPj9z5kwsXrwYX3/9NWbMmNHcp+FwHPIjp1NcXiX/npKtc2BLiIjI1IQJE3D27Fns3LlTXrZq1SqMHDkSGo0Gp0+fRkxMjNXn4oWHh8PHxwenT582Wx4REQEvLy/51bVr1xY/j5bAHipyOkWmgSpLi1uj/R3YGiKilufuqsTJ+YMddmx7xcXFoW/fvli5ciUGDBiA1NRU7Nq1C/Pnz5e3qe+Gotpha9euXfD29pbfN+XxL47EQEVOx7SHKjmLPVREdP0TBMGuYTdnMHHiRDz11FP48MMPsWrVKsTExKB///4AgNjYWOzevRsVFRUWwSkjIwNarRYdO3Y0Wx4dHX1d1FBxyI+cTnGFaQ8VAxURkTN5+OGHoVAo8MUXX+Czzz7DhAkT5HqqRx99FEVFRVi2bJnF5xYuXAg3NzeMGjXqWjf5mmgdcZhuKFJROmAMVKIo8plXREROwsvLC6NGjcLs2bOh1Woxbtw4eV18fDymT5+O559/HhUVFRgxYgQqKyvx+eef4/3338fq1asREBBgtr+cnByUlZWZLQsICGh1Q38MVOR0TIf8dOVVSC8oRYQfn8pOROQsJk6ciBUrVuD+++9HeHi42bp3330XPXr0wEcffYQ5c+agrKwMKpUKv/76K+666y6LfXXq1Mli2d69e3H77be3WPtbAof8yOmYFqUDHPYjInI28fHxEEURP/zwg9X1EyZMwMGDB1FaWorz588jNDQUH330EfT6mhGIAQMGQBRFq6/WFqYABipyQlINlcbd2N3LwnQiotarXbt2SEpKQlxcHI4cOeLo5rQYDvmR05FqqHpH+eHX5Bz2UBERtXLR0dH4z3/+4+hmtCj2UJHTkYb8ekf5AeCQHxEROT8GKnI6UlH6zW2NgepsbhEqqgyObBIREVGdGKjIqRgMIkoqjEN+sSFe8Fa7oMog4tyVIge3jIiIyDYGKnIqppN6eqld0DHU+DgCDvsREZEzY6AipyIVpCsVAtQuCnSqDlS804+IiJwZAxU5Fakg3VOlhCAIiGMPFRERtQIMVORUSqqH/LzUxhk9OoUwUBERkfNjoCKnIvdQVQequFAfAEB6QSm0ZZUOaxcREVFdGKjIqUg1VFKg0ni4IkzjBgA4zV4qIiKnsHfvXiiVSgwdOtRseVpaGgRBsHg98cQTZtutWbMGt9xyCzw8PODt7Y3+/ftj8+bNZtskJSVBEAT4+flZPDz5wIED8r5NiaKITz75BLfddhu8vLzg6+uLPn364N1330VJSUkzXgFLDFTkVKQ5qKQhPwByYfopBqoWUaXnHF9E1DArVqzAU089hd9++w0ZGRkW63/55RdkZmbKrw8//FBe99xzz+HJJ5/EqFGjcPToUezfvx933nknhg8fjg8++MBiX97e3vjmm28sjt+2bVuLbUePHo0ZM2Zg+PDh2LFjB44cOYKXX34ZmzZtwtatW5vhzG3jo2fIqdQM+SnlZZ1CvZGUkouULK2jmnXduphXgqFLdmFUn0jMSeji6OYQ3bhEEahs2R4Um1w9gFo9PXUpKirC//73Pxw8eBBZWVlYvXo1XnzxRbNtAgICEBoaavHZP/74A4sWLcL777+Pp556Sl7++uuvo6ysDM8++yyGDx+OyMhIed3YsWOxcuVKPProowCA0tJSrFu3Dk8//TReffVVebv169dj7dq1+PbbbzF8+HB5ebt27fDAAw9Aq23Z7xAGKnIqxbVqqADwTr8WdPhSPnRlVfjpeBYDFZEjVZYAb4Q75tgvZgAqT7s3X79+PeLi4tCpUyc88cQTmDFjBmbPnm0x/GbNl19+CS8vLzz55JMW62bOnInFixfj66+/xowZM+Tlo0ePxttvv42LFy+ibdu2+Prrr9GuXTvcfPPNZp9fu3YtOnXqZBamJIIgQKPR2H2OjdEqhvzS0tIwceJEREdHw93dHTExMZg3bx4qKirkbZKSkjB8+HCEhYXB09MTPXv2xNq1a23uc926dRAEASNGjKj3+ElJSbj55puhVqvRoUMHrF69uhnOiqyRA5WqJlDFBhsD1fkrxQ5p0/WsoMRY6J9eUCrfYUlEVJcVK1bINVFDhgxBYWEhdu7cabZN37594eXlJb8OHz4MADh9+jRiYmKgUqks9hseHg4fHx+cPn3abHlwcDDuu+8++bt35cqVmDBhgsXnz5w5g06dOjXHKTZKq+ihSk5OhsFgwLJly9ChQwccP34ckydPRnFxMRYuXAgA2LNnD3r06IEXXngBISEh2Lx5M8aMGQONRoOEhASz/aWlpeG5555Dv3796j32+fPnMXToUEydOhVr167F9u3bMWnSJISFhWHw4MEtcr43sqJaRekAEOilBmD88hdF0a7/CyL7SIEKAM7lFqNbm5b9PzgissHVw9hT5Khj2yklJQX79++Xa5pcXFwwatQorFixAgMGDJC3+9///ofOnTvL702H8ERRrPMY1sLWhAkTMH36dDzxxBPYu3cvNmzYgF27dpltU99+W1qrCFRDhgzBkCFD5Pft27dHSkoKli5dKgeq2uO306dPx9atW7Fx40azQKXX6/H444/jlVdewa5du1BQUFDnsT/++GNER0dj0aJFAIDOnTtj9+7deOedd2wGqvLycpSXl8vvW3rc9npSU5ReU0Pl6+EKAKgyiCiu0JsVrFPT5JfU9PKm5hQxUBE5iiA0aNjNUVasWIGqqiqEh9cMT4qiCLVabVZQHhkZiQ4dOlh8PjY2Frt370ZFRYVFcMrIyIBWq0XHjh0tPnffffdhypQpmDhxIoYNG4aAgACLbTp27Ijk5OSmnF6TtIohP2sKCwvh7+/f4G3mz5+P4OBgTJw40a7j7N27F4MGDTJbNnjwYOzdu9fmZxITE6HRaOSXaTKnuhVVWNZQubkqoXYx/lXNL66w+jlqnMLSmh6q1Bw+gJqIbKuqqsJnn32GRYsW4ciRI/Lrr7/+Qnh4OL788st69/Hoo4+iqKgIy5Yts1i3cOFCuLm5YdSoURbrXFxcMGbMGCQlJVkd7gOAxx57DKdPn8amTZss1omiiMLCQjvOsvFaZaBKTU3FkiVLrBa1SdavX48DBw5g/Pjx8rLdu3djxYoVWL58ud3HysrKQkhIiNmykJAQaLValJaWWv3M7NmzUVhYKL8uXbpk9/FudCVWitKBml4q0wBATVdQq4eKiMiWzZs3Iz8/HxMnTkS3bt3MXiNHjsSKFSvq3Ud8fDymT5+O559/HosWLcLZs2eRnJyMOXPm4P3338fy5cut9j4BwKuvvorc3Fybo0MPP/wwRo0ahUcffRRvvPEGDh48iAsXLmDz5s0YNGgQduzY0aTzr49Dx05mzZqFBQsW1LnNqVOnEBcXJ79PT0/HkCFD8NBDD2Hy5MlWP7Njxw6MHz8ey5cvR9euXQEAOp0Oo0ePxvLlyxEYGNh8J2GFWq2GWq1u0WNcr6SJPWsP6/m6q5CtLTer+aGmyze5nqm5DFREZNuKFSswaNAgq3fLjRw5Em+99ZZdJS7vvvsuevTogY8++ghz5sxBWVkZVCoVfv31V9x11102P6dSqer8/hYEAV988QU++eQTrFy5Eq+//jpcXFwQGxuLMWPGtHjds0MD1cyZMzFu3Lg6t2nfvr38e0ZGBgYOHIi+ffvik08+sbr9zp07MWzYMLzzzjsYM2aMvPzs2bNIS0vDsGHD5GUGg3FCQxcXF6SkpCAmJsZif6GhocjOzjZblp2dDR8fH7i7u9d7jtQwtR89I9FU91AVlHLIrzmZ9vilXSlGpd4AV2Wr7Lgmohb2/fff21x36623ykXh9hSHT5gwQR66S0tLQ//+/fHRRx/hjjvugFJprKEdMGBAnfsaMWKExXqFQoGpU6di6tSp9bahuTk0UAUFBSEoKMiubdPT0zFw4ED07t0bq1atgkJh+Y9+UlISEhISsGDBAkyZMsVsXVxcHI4dO2a2bM6cOdDpdHjvvfds1jnFx8fjxx9/NFu2bds2xMfH29VuapjiCsuidADwda8OVOyhalamQ35VBhEX8krQIdjLgS0iohtNu3btkJSUhDVr1uDIkSPo3bu3o5vUKK3idqn09HQMGDAAUVFRWLhwIXJzc+V10kysO3bsQEJCAqZPn46RI0ciKysLgLGL0N/fH25ubujWrZvZfn19fQHAbPns2bORnp6Ozz77DAAwdepUfPDBB/j3v/+NCRMm4Ndff8X69evxww8/tOQp37CsTewJsIaqJRgMonw9Q3zUyNaWIzWniIGKiK656Oho/Oc//3F0M5qkVfTtb9u2Dampqdi+fTsiIiIQFhYmvyRr1qxBSUkJEhMTzdY/+OCDDTpWZmYmLl68KL+Pjo7GDz/8gG3btuGmm27CokWL8Omnn3IOqhZSZGViTwDw9TDeXmvao0JNoyurgqG6t7xPlPFu2LOsoyIiapRW0UM1bty4emutVq9e3eAZzK1tb23ZgAED5FleqeVU6Q0oqzTWtVnUUHHIr9lJ9WgeKiW6hPvgh2OZvNOP6Bpy9ESUZNRcfw6tooeKbgzFFXr5d8/aNVRyUToDVXORwqmvuytigozDfAxURC3P1dX471lJiYMehkxmpMfYScXwjdUqeqjoxiDVT7kqBahdahelG4f8CtlD1WykWdI1Hiq5bupsbhEMBhEKBR/vQ9RSlEolfH19kZOTAwDw8PDgI7UcxGAwIDc3Fx4eHnBxaVokYqAip2GrIB0w7aFiDVVzkQrS/TxcERXgAReFgJIKPTK1ZWjjyylBiFqSdEOVFKrIcRQKBdq2bdvkUMtARU7DVkE6wBqqliAP+Xm4wlWpQLtAT6TmFCE1p4iBiqiFCYKAsLAwBAcHo7KS/645kkqlsjoVU0MxUJHTKKmwPks6APh5Vt/lV1oJURTZPd4M5CG/6uHUDkFecqDq39G++eGIqGmUSmWTa3fIObAonZxGzSzplv+4SBN7VlTV3AlITSP1UPlVD6dKdVQsTCciajgGKnIaddVQeaiUcFUae6XyORdVs5BqqHxrBaqzDFRERA3GQEVOQwpU1ob8BEGQh6ZYR9U8pGAq3UEp91Bxck8iogZjoCKnUVRurKHysFKUDvBOv+ZmWpQOAO2DPAEAV4srcLWY15iIqCEYqMhp1PRQWS/QlOqoOBdV86gZ8jP2UHmoXOS7+1hHRUTUMAxU5DSK6qihAjhbenOTh/yqryvAwnQiosZioCKnUVdROgDWUDVCjq4Mn+46Z9GrZzCIFkXpAAMVEVFjMVCR0yiusF2UDrCGqjGW7TyH1344hc/2ppkt15VVQXoeqDRpKoCaZ/qxMJ2IqEEYqMhpSEXpNof8WEPVYBfyigEA568Umy2XQqmHSmn23EROnUBE1DgMVOQ0SuorSvfg42caKktbBgC4nF9qtjxfntRTZbY8KsBD/lyVnhOoEhHZi4GKnEZ9RekaD+nxMxzys1dWYTkA4HJ+idnyAvmxM65mywO91FAqBOgNIq4U8ToTEdmLgYqchlRDVd+QH3uo7FOpNyCv2BiosrRlqKiq6XGyVpAOAEqFgBBvNQAgs9C8V4uIiGxjoCKnUSzVUNUzsWchp02wS46uXC48N4hAVmGZvC6/euLO2kN+ABCqcQNgvj0REdWNgYqcRl0PRwZqHpHCHir71A5EpsN+0lxemlo9VAAQpjFO7pnJQEVEZDcGKnIKlXqDPCRla9oE6cu/tFKPskr9NWtba5WtrRWoCmqG8OTHzrhbBiq5h0rLQEVEZC8GKnIK0qSegO0aKh83FygVAgAO+9nDsofKNFDZHvILqw5U7KEiIrIfAxU5BWm4T+WigKvS+l9LQRDku9I47Fc/qYdKCqH2DvnV1FCxKJ2IyF4MVOQUpIJ0W8N9kpo7/XhLf32kIbsuYT4AavdQ2R7yYw8VEVHDMVCRU6ivIF2i4QOS7SYN+fWO8gMApFsb8vO0dpefsSg9W1sGg0Fs6WYSEV0XGKjIKcgPRrYxZYKEj5+xn9RDJQWqzMJSVFbPfi4FUms9VMHeaggCUKkXkVfMnkAiInswUJFTKKnnwcgSX86WbhdRFOUequ5tNFC5KOS5qAwGUS7qt1ZD5apUIMjLOLkn56IiIrIPAxU5hfoejCxhUbp9CksrUV49DUWoxg0RvsZhvMv5pdCWVcoTfkpze9VWU0fFwnQiInswUJFTKLazhsqXNVR2kYb7fD1c4eaqRBs/KVCVyGHUU6WEysX6PwGci4qIqGEYqMgpFLGGqllJQ3WhPsZgFOFX00Ml109ZmYNKwtnSiYgahoGKnEJNDxVrqJqDNAeV1NMU4ecBoDpQVd/hV/vByKb4PD8iooZhoCKnIAWq+orS5WkT2ENVp6zCcgDWeqhqhvzqClSsoSIiahgGKnIK9hal+7Io3S5S7VOItSE/qYfKRkE6UBPE2ENFRGQfBipyCjU9VPUVpRtDAJ/lVzdbQ35Z2jJ5bqm6e6hqaqhEkZN7EhHVh4GKnEJxhZ01VNU9VEXlVfIklWSpdlF6kJcaKqUCeoOIU5k6AHUHqmAf4zxU5VUG9gYSEdmBgYqcQpGdRek+JjN784vetuxaQ34KhYBwX+PvJzIKAdQ95OfmqkRA9WNpeKcfEVH9GKjIKZTY+XBkpUKAj5txm0Le6WdVeZVeHtaThvyAmmE/KSDV1UNl+tksLQvTiYjqw0BFTkHqofJQ1V1DBdQ80Jc9VNblaI13+KlcFPAzCU1SYbqkrnmoANM7/dhDRURUHwYqcgrFdj7LD+CdfvWpGe5TQxAEeblloLKzh4qBioioXgxU5BTsndgTADTy5J4MVNZIUyZIBekSachP4ldPoOJs6URE9mOgIocrr9KjUm+8Nd+eQFXTQ8UaqhxdGY6nF5otk3qUQiwClXkPlaaOonSAc1ERETUEAxU5XHF1QTpgfGBvfaShqhthLqrUnCJkFFgvCjcYRDzx6T488MFuHLlUIC+XAlCYpu4eKo17fT1UnC2diMherSJQpaWlYeLEiYiOjoa7uztiYmIwb948VFTU9FAkJSVh+PDhCAsLg6enJ3r27Im1a9fa3Oe6desgCAJGjBhR57E3btyIe++9F0FBQfDx8UF8fDy2bNnSXKdGqBnuc3NVwEVZ/1/JG6WG6lxuEe5/bxeGf/g7SqprzEz9fvYKTmcXwSAC6/ZflJfXniVdEuythqvSWFPlpXaByqXuax1qUpTOyT2JiOrWKgJVcnIyDAYDli1bhhMnTuCdd97Bxx9/jBdffFHeZs+ePejRowe+/vprHD16FOPHj8eYMWOwefNmi/2lpaXhueeeQ79+/eo99m+//YZ7770XP/74Iw4dOoSBAwdi2LBhOHz4cLOe442syM7n+ElulBqq5bvOoUJvQK6uHBsOXrZY//kfF+TfNx/NlENX7VnSJca5qIzDfvX1Tpl+vqRCD125ZaAjIqIa9n2DOdiQIUMwZMgQ+X379u2RkpKCpUuXYuHChQBgFq4AYPr06di6dSs2btyIhIQEebler8fjjz+OV155Bbt27UJBQUGdx3733XfN3r/xxhvYtGkTvv/+e/Tq1atpJ0YAGlaQDrTeGiptWSV2nb6CO2MD6w00OboyfH0oXX7/6e5zePy2tnIPXmZhKbadzAYA+HuqcLW4Aj8fz8KDN0fYLEoHjHVUF/JK6r3DDwA8VC7QuLuisLQSWYVl8HGr/zNERDeqVtFDZU1hYSH8/f0bvM38+fMRHByMiRMnNuq4BoMBOp2uzmOXl5dDq9Wavcg2eZZ0lZ2BqhXWUG0/lY17F+/EtC/+xN/e2YlfqsOQLat/T0OF3oAeERr4ebji0tVS/HwiS17/5f5LMIjAbdH+GNe3HQBgw8HLEEUR2dXzUNUe8gOACF9jHZVfPXNQSTgXFRGRfVploEpNTcWSJUvw5JNP2txm/fr1OHDgAMaPHy8v2717N1asWIHly5c3+tgLFy5EUVERHn74YZvbJCYmQqPRyK/IyMhGH+9GUFJh3yzpEilQtYYaqvziCjzzvyOYuOYgsrXlcFEIyNaWY9JnBzF93WFcLbbsZSsqr8J/q4fzpg3sgNHx7QAAy387B1EUUak3yDVTT9wehQdvbgNBAPaey8Ox9EJUVBmfcWg1UFXf6WdPDxVgOhcVC9OJiOri0CG/WbNmYcGCBXVuc+rUKcTFxcnv09PTMWTIEDz00EOYPHmy1c/s2LED48ePx/Lly9G1a1cAgE6nw+jRo7F8+XIEBgY2qr1ffPEFXnnlFWzatAnBwcE2t5s9ezaeffZZ+b1Wq2WoqoM8S7q6/jv8gJrb/WsP+RWUVOB0dhFOZ+twJluH/JJKRAd6omOINzqGeCFU44YLeSXG9TlFuHS1BOG+7ugQ7IWOId5oH+SJK7pynM4uwplsHVJzi+CpdkHH6vWxId4wiCJOZ+uQkqXDmewiVOoNiKle3ynEGxp3V6TmGj9/OrsI3/2VjitFFVAIwOR+7fGPATFYuvMslv92DpuOZGD3mStYMLIHBnUJkc9j3f6L0JVVoX2QJ+7tHII+UX5YtvMs/rpciH3nryK/uAI5unIEeqkxuGsoVC4K9I0JwO+pefhwRyoAIMBTZbXofHjPNjh0MR+P3dbWrmvNHioiIvs4NFDNnDkT48aNq3Ob9u3by79nZGRg4MCB6Nu3Lz755BOr2+/cuRPDhg3DO++8gzFjxsjLz549i7S0NAwbNkxeZjAY/0/excUFKSkpiImJsdmOdevWYdKkSdiwYQMGDRpUZ5vVajXUanWd2zSHHF0Zbn19e7Psy2RC7WtOuoHM7hqq6t4VbVkVomf/YLEfZ9MxxAtv/f0m9Iz0BQDMvq8z7u8Whue/+guns4sw6bODmDEoFk/fHYsqg4hPd50HADx5V3soFAICvNT4e+8IrN13EZ/8dg5llcYevVG3RMih6aHekfg9NQ9bThiHEq31TgFA2wAPrB5/q91tD/Ux9mhxLioioro5NFAFBQUhKCjIrm3T09MxcOBA9O7dG6tWrYJCYfl/30lJSUhISMCCBQswZcoUs3VxcXE4duyY2bI5c+ZAp9Phvffeq7MH6csvv8SECROwbt06DB061K72tjbOEEZuj667Jk7i76FCXKg3krN0Fu1u4+uOjiHGHiM/TxXO5hThdE4RUrN1KK7Qw8/DFbHVPVZt/T2QUVCG09W9SVeKyuHmqkBssDdig73QIcQLRWVVxh6rHB0uXi2BQhAQFeCBjsHGfahcFHKv2LncYlToDQj1cUNsdRu6hvtgaI8wqF3Me99uivTF90/didd/OIXP9l7Au7+cwYkMLe6ICUCWtgxB3mqM6NVG3n5yv/b4Yv9F/JqcA8AYgB+9taaXaXDXUHirXeS78Wrf4ddY7KEiIrJPq7jLLz09HQMGDEBUVBQWLlyI3NxceV1oaCgA4zBfQkICpk+fjpEjRyIry1jAq1Kp4O/vDzc3N3Tr1s1sv76+vgBgtnz27NlIT0/HZ599BsA4zDd27Fi89957uO222+T9uru7Q6PRtNg52yPAU42Dcyx7y5oajkSIENC0LisRDWuEWqmExs66HoVCwA9P97OoP/JQKW32chkMIooqquCtdjF7vp0pXVklPFUuUCisry+r1EMQYBGOJFV6A8qqDHbXgqldlJg/vBu6tdFgzjfHse1ktnzn3oQ7os2O0y7QE0O6huKn48a/f3d3CjabqNNdpUTCTeH4srq2ylYPVUPxeX5ERPZpFYFq27ZtSE1NRWpqKiIiIszWSRMOrlmzBiUlJUhMTERiYqK8vn///khKSrL7WJmZmbh4sWaSxE8++QRVVVWYNm0apk2bJi8fO3YsVq9e3bgTaiZKhYBAr5YfWnRGSoWAIG/7z12hEOq97d+7nvVurnXXeLkoFfCyY2LS2h7uE4mOId6Y+t9DyNKWwUvtgsdvt6xxmnJXezlQPREfZbH+oT4RcqCyNmVCY3C2dCIi+wgip0BucVqtFhqNBoWFhfDx8XF0c8hJ5ejKsDTpLPrFBuLuuBCr2yzamoLC0kr8Z1hXi540URRx7zu/ITWnCIseugkje0dY3UdD6Moq0f0/WwEAf839m929iERE14OGfH+3ih4qohtBsLcb5g3rWuc2M//WyeY6QRDwzsM98d1f6RjaI6xZ2uTt5op2AR5IyyvBoYtXbQY9IqIbXauch4qIrOseocFLQ7vUOzzZELdW3yyw/3x+s+2TiOh6w0BFRHW6NToAALD/fJ6DW0JE5LwYqIioTrdV91AdvVyI0upZ7YmIyBwDFRHVKcLPHWEaN1QZRBy+yGE/IiJrGKiIqE6CIMh1VPvOX3Vwa8jRjl0uxEvfHENeUbmjm0LkVBioiKheNYXpDFQ3uo+SUrF230V8czjd0U0hcioMVERUL6mO6s+L+aioMji4NeRIl/NLzX4SkREDFRHVKybIC/6eKpRXGXAsvcDRzSEHkmbN5+z5ROYYqIioXoIg4NZ2rKO60ZVV6nGlyPgMTT4wm8gcAxUR2UWqozrAQHXDMn1IdkaB43qocnRlWLbzrMUD0iV6g4iVu8/jeHrhNW4Z3cgYqIjILlKgOpiWD72BjwC9EWWYDPNdKapAeVXj5iUzGEQ05TGyy387h8SfkrF6T5rV9btTr2D+5pOYvfFYo49B1FAMVERkl85hPvBWu0BXXoVTmVpHN4ccILPAfJgvqxHDfuVVetz7zk48sWJfo9txNrfY+DOnyPr66uXJWVreREHXDAMVEdlFqRDQu50fAE6fcKOqPcyX3ohhvzPZRTibW4zfU/OgLatsVDsu5BkD1cWrJVbXS8sr9SLO5loPXUTNjYGKiOzG+ahubBm1eqRq91jZ40JeTQi6mGc9ENXFYBBxqXrKhvoCFQD2ptI14+LoBhBR63GbPGN6HracyJKXB3mrcXNbP0c1i64RaaoEV6WASr3YqKkT0qp7l6Tfu7XRNOjzObpyeRivsLQShSWV0Hi4mm3DQEWOwEBFRHbr3sYXbq4K5JdU4sn/HjJbt/7JeLkHi65PUo9U9zYa/HmxwKLHyh4XTALVBRs9VLqySuTqytE+yMtiXe1eqYtXS9DdoyaUGQwiLpkFKl2D20jUGBzyIyK7qVwUeOn+zugd5Se/Ar3UAIyzqNP1Taqh6h3lZ/a+IdJMQlTalWKr20z9/BDuWbwTqVaKzq0FKlM5unKUmxSin8rUNumOQiJ7sYeKiBpkdHw7jI5vJ79fsv0MFm07jTPZLP69nunKKqErrwKA6uHd842qoTKtm7LWQ6U3iDiYlg9RBP68kI8Owea9VLUD1IWrxVbXB3urcaWoHHnFFcjVlSPYx63BbSVqCPZQEVGTxIYYv/DO5NgeWimpqKqzl6CsUo8qve3b2yv1hjrnPDIYRJRW1D0nUmmFHoY65s8qr9Kjso42VOkNKKu0fQxRFFFSUVVnG8oq9XXO4VVRZajzNn+DQWz03E9NJc2MrnF3lf/MMxpYQ1VaoUeWtiaEmdZTSS5eLZF7mE5nW/6dkobzPFVKs/emnweAjiHeiA70BACcZB0VXQMMVETUJB2CvQEAqTlFVgPL9lPZ6DJ3C9buu2j18zm6Mtz2xnZM++JPm8d44tN9uOutHTZvs1+68yy6zPvZ5t2Hf17MR9d5P+O97Wesri8ur8LAt5Mwatlem8Fvxv+O4JbXfrFZiL3h4GV0mbsFPx3LtLr+XG4Res7fipc3Hbe6Xm8QMWzJbgx59zeboer1H0+h+3+2Wg0aAHA8vRB3L0rC9lPZVtc3hTRFQpjGDWEadwCArqwKugZMfSCFHbWL8asnR1duEUJNz+20lSE/qQbr9vYBZvusfYxIfw90DvMBwDoqujYYqIioSdoFeMBVKaCkQm+1x+KH6oDxo42gsSc1D4Wlldh+KsdqD1BeUTn2nb+KbG05Dl2wXqf1w9FMiCLw8/Esq+u3nMiCQbTdhiOXjAXWf14sQLa23GJ9pd6ArSezoSuvwq4zV6y3QTpPG234NTkHZZUG/Hgs02poO5tbhJRsHc5dKUZKlvUA8MPRTFRUGbDtpPXA9L8Dl3Autxif/3HB6vqmkIb3wn3d4al2gY+bsWKkIc/0k3qk4kK94Vt9Z17tYb8zJoHqjJXgePGq8e/YnbGBVj9/sfoYUQGmgYo9VNTyGKiIqElclAq0D5SG/Sx7FE6kG7/MjqcXWg0S0vPWqgyi1SBxPKPmy/CElWezlVfp5V6N4xnWn90mteFsbpHVYTnTZ75Ze/7bmewiudfIWhtEUZQ/Z209AJyoPo+CkkpczrcMnmZtsHIeubpyebjshI3zPFq9j2PpzV+ILfXMhWmMtUjhvsZeqoYUpku9S20DPBEV4Gm2THLapBYvs7DMrAespKIKV4qMgffODoHy8U2HaqUeqrb+HujCQEXXEAMVETVZh+qamtRahemlFXq5tkpbVmU1SByrJ0iYhx3LL8aULB2qqocaT2ZoLYYdRVGUj2EQrX+5mrbhmJVAdLye9VnaMuRVP6j33JViq8Ngpp+zFojqbUNG3esr9Qb53K4UlVvtaWuKDJMeKtOfDemhknqT2gV4oF2ABwDzu/4Ay7op05B+qbp3SuPuig7BXnBzVcAgmoc6qQerrcmQ37krxXXWvxE1BwYqImqy2GDrhemnsrQwzTe1g4DBIOKkSQ+Utd6h+sKMacgqKq+yKHS+nF+KwtJKq9tLTpj2glkLdSbLTmZqLQrLa+/T9JwAY8+K6SNQ6gtt1nq5jl+uWXbpaikKS8xDm2kvmq1jNIUUWsJ9jT1UUk9VZoN6qIzhKcpGD1WV3oBzuTVDdoD5sJ9p75MgCGjr72G23+Lymh6stgEeCPFRw8/DFXqDyLtQqcUxUBFRk8VWF6bXHvKrHZBqv79wtUS+Fd+4vu7eo/SCUuRX9wRZW2/tfe2AVHu9rqwS503mQ6ov7JRVGiyeD1dfG05maCGaBUvz8zQYRLNQdypLZ3HHYe191u7Nq31tmztQ1Qz5mfdQpTdg6gQp7Jr1UF0xmUbhagkq9Aa4uyoxsFMwAJgFIdNAZfpTWn4p3/jT18MVPm6uEASBdVR0zTBQEVGTxZoM+ZnW7khf8sHexsk/LUJBrfUpWTqzXpaCkgp5mDCoepvaQUIKTNI+TtTqHTpW6xi1g4e0faCXCoIAZGvLkaOrCQlVeoN8273NfdRzDOl9iI9a3t70Op27UoySCj3cXZXwdnNBRZXBYuhLaqetaym9l4rFrfX2NZYoivLQXht5yK+6h8rOqRPKq/RyL5etHiqpN6pDsBc6hRpD+mmzIb+aO/hMf0rLpTmupKAFQA5UnDqBWhoDFRE1WbsATygVAnTlVWbzDEk9MQ/1iQBgDAXWAte9XUKgcXdFhd48SEg9VlEBHvJjbUx7sSqqDEiuviVeOsaxy7XDjHkbzuQUmdXTSG24ua0fYqofdXLC5BjG+hsDPFVK3N89rPq8rPcWPdwnsvp97VBnfP9/vSKgVAjIK64wqz2S2tAl3AfdwjUWbbhaXCFPW/D33hFmn6k5hvH9gzdHWG1jU1wtrkB5lQGCAIT4SEN+DauhupxfCoMIeKiUCPRSyT1UGYVl8p+HVJAeG+JVM4xsY8jP9Kc05Fd7PQD2UNE1w0BFRE2mclHIX5DSEE1ZpV7+Mvx770i4KARcLa4we/6bFES6t9GgWxvjF5+1u926hWvQvfohuqbrz+ToUKE3wMfNBfd1C5M/I4U207vvBnUOQYCnCnqDiOQs09Bm0oZwyzZIAa2rSRtMw06OrgzZ2nIIQk2gOptbhGKToUypF613lJ8cFMyK8dOl8/SRr4O19dGBnoiPCbBoY5VJQfpDfSKgEIx3BWZrGz6TuTVSQXqglxqq6jmkwjU1d/nZc0fhBXk6A08IggB/TxW81cbeNKmHSRoy7hjijdgQYw9VZmGZPP+YfJdgdWCS6qykIGU9UBn3w0fQUEtjoCKiZlG7jkq6+87f09gbIX1BSkHAGHaMIaBbG43cM2PtbjZb64+brO8Y4g2VUgFdWZX8xSrdfadUGGtpulUHIvM7C03aYHV9daBq44PuEdWBKqNQvptQClcxQV5yIbRocjdhWaVevibd21gPhmbn2cb2dehq0oOVllciB40zOUUorzLAS+2CzqE+8uNaavfWNZY0v1i4pubxLSEa49BjeZUBV2vVtVljeocfAAiCgKhA8zv9pADeMcQLGndXeYhUmjT2Un7NHXymPy9dLYEoilYDVYdgL7goBGjLqhr1MGciezFQEVGzkOuocsznhOoa7gNBENC9Vg+UdPedSqlAxxBvkzBjOe9UtzY1PTcX8krku/ZMA5nKRYG4MG+z5VKgiA32gpurUt6HtN/i8pq770zDzAkrdx52b6NB+0BPuLkqUFyhx/nq3pJjJutNf0rLT1XfFRjopUKIj1oOZdJ+TQvSu0fUtOFUplZ+HM8Jk548P0+VXMckhTnTwKVQCFaDYVNkynf4ucvL1C5Kua7NnmE/0zv8JKZ1VKZ3+EnhvGN1CD+TrUOOrhwVVQYoFYJcvxXhZwxOuvIqFJRUWq2hUrso5YB5KoPDftRyGKiIqFnEyl9+xoByvFbQ6FarZ0b6su8U6g2Vi8IsSFTqDdCWVco9F93CNfD1UCHCrzpIZJjvQ/ps13DzIGHa+2TaFtOwI4rGYvEgbzW6Vg/5pReU4mpxhVnY6dZGAxelQp4sUjqP4yZhxrQNx00mNJU+LwiCSRuNQ1AXrpagqLwKahcFOgR5ITrAE54qZfXdhHWHttptsLW+qaSeHaluSiL1WNkzuafpHX6SmrmoipGWZ7zDz0OllAOjFIROZxfJvU9tfN3hojR+dbm5KhFaXdN1Pq9YvoGhrckxANZR0bXBQEVEzaJmLqoii+E8059SkDhu0vsEAFH+HvBWG+9wS80pkntf2vi6w89TBcA8KJjWDUm1T3KNU4b1oCGFmdPZOpRX6S2Cirebq/xA3ePphfLdd26uCrlgvXYwrC/MyNeh+thdwnygEGom35Ta0DnMBy5KBRQKwSwYFpZUyhNaSsvlerNawVLq/aodHJuq9hxUkjCN/bOlSz1UpmGnpoeqRB7uiw32gkIhADDpocopsjqcZ/r+wPmrqNAb4KIQLIKfXEeVxUBFLYeBioiaRXSgJxQCUFhaiYzCMvkxMtKXe+fQmiCRoyu36F1SKAR0MSkKrx1UTLc9nq5Fam5N3VC76i9m04Jua6Etws8dvh6uqNSLOJ1VJIcdKaiYHuNYeqEczLqE+UBZ/SVvuj6vqFzuvZHaLq0/k6NDaYW+prC+erm7SikPaR1PL5SHH62fZ6H8+bb+HtBUP//OtA2mwVLaR5dw47XO0ZUjpxkK0zNt9FCFyVMn1H2MKr1BLjxvZzLkJ/2ellcs3+EnPWwbMNZSAcYhv4u1pkyQSO93pxqfsRjh5y7/WUmkHirp7wVRS2CgIqJm4eaqlHscfjiagQq9ARp3V3mYzjRIHLtcWDOUZhJmanqYtCZBxEdeb1qwLYWhLtV1Q4Bx+NBFIaCgpBKHLxUgR1cOhVDzhSoIgllxu2ltknyM6mB0IsNGqDOZ1kAKhe0DPeHtZgw7IT5qBHqpYRCBvy4XyNNAmJ5HV5Pgdyzd8jy7R1gGS2vX4fyVYvx1uQBllebB0kPlIveoNUcvVaaNHippaK6+Yu+MgjJUGUSoXBTyEB1QM+SXnl+Kk5nGdkohCqgJV5mFZThpEixNSXf67T9/1bjeJLBJekT4ws1VgUtXS20+3Lo+l66WYNHWFJsF+GWVeny4IxVHLxfY3MeGg5fw7eF0m+sPXcjHhztSUV5l/TE5V4rKsXjbaZtzf4miiG8OX7b5rEcA2HcuD1tOWH+ANwBka8vw1aHLFhPLUv0YqIio2Ug1Lxv/NH5pdGtjLEiXSEFi28lsXC2ugItCkCdwNG5f0/NSuwcLqAk7568U449zedXLatarXZTyMNH/9l8CYLz7zkPlYnGMA2lX5bvvulnpHTJtQ1eT9bEhXlC5KKArr8KPxzIt1guCIIefjX9eRqVehJ+Hqxw+APMhueNWz7MmWP5V/QVtuj7QS40wjRtEEfjfAeN5djUJlqbHsDb7fENU6Q3y3GKmRemAyVxU9Qz5XbhaPWWCv4dZG4O81XB3VcIgArurg4705wfA7E4/qQfK1pBfefWEsG39zdso7eexW6MAAB/sSK2zrdYYDCL+ufZPLPk1FS99c8zqNu9vP4O3t6Rg0pqDZlNmSPakXsHzXx3FjP8dwcG0qxbrC0sqMfmzg3h7SwqWJp21eowXvjqK97efwVNfHLba07b+4CU887+/8MSn+1BQYhn80q4UY/SK/Xjyv4fk621KbxAxcc0BPLfhL7z7y2mrbSDbGKiIqNlIdVTSPE+mYQeo+ZL/7q8M4/Yh3nBzVcrrTYe6pMfBmA7HBXipEV4dJH44agwzUm+OrWOY9i4Zj2Hc/qdjWdV336nlL23TNl+6Woq/Lln2ULkqFehcHQK//6u6DW2st0FaLxWk117/e+oVaMuq5DsdJe2DvODuqkRppR47knOtnod0XWraUPs8m6eOKkdXDoMIuCgEBHqpzdZJQ3711VClWbnDD6ieOqG6h6m4wtgrE2vSQwXUBKyySikwWR/yk9ReL5lyV3uolArsP39V7s0ydTa3CDuSc6wGla/+vCxfx5+OZ2HPWfMwcjGvBJ/uOg/AeL0+SjIPbVV6A+ZvPim//8/3Jywe4v3u9tNy79fHO8/KE7lKdqTkYHtyDgDg4IV8fHvEvKeroKQCC35OAQDkl1Ri0VbLQDR/80lUVPc8zd98Qr6LVLL+4CU5gC//7bzZI5mofgxURNRsan8ZdrXxJV9aPTO21OMkiQ70hIdKifIqA0QRCPVxk2/Nr73Pmn1YD0zSeos2hNf6fK1eNI2HKyKrezlKK/VQuSjknjfb52E97NQcw3x95zAfCEJNr0pcmDdclTX/HCtN6slsHaN7rTZIBeny+lrTMzSWNLwUqnGzqE2Set2ydeUWD4w2deGK+QOPTZnWVHma3OEniTWpqQJsD/nVrLcc8pPaP7J6lvnavVQpWTo8sGQ3xq8+gBW7z5ut05VV4q3qoCINX8///qRZGHntB2NQkdq+fNd5eQoHAPjywCUkZ+mgcXeFt5sLjqdrseHQJXn9mWwdPtt7AQAQ6e+OskoDEn88Ja+vqDLg1e9Pmp1/4o/JKDLpCVu01RjIpEcTrd13wWzob/upbPyanANXpQCNuytOZxfhi/0X5fWFpZV4e4vxPP09VajQG/Cf706w5qwBGKiIqNnU/vKr3WvSpTpIyOtrhQClQpCnHwAsg0jtfbq7KtE+yHrYsdWGqAAPeLu52Fxfe1nnMB+zsGPtGLVDW+3zqh2GPNU1NU7W9le7DaZ3OtYcwzyM1t6HdK2ztGXI1ZXLy6v0hgZ9SUqzpIdrLIfSAr3UcFEI0BtEs+cf1pZWa1JPU9LkngDQIcTbLNwC5jVVGndXuTBfEuCpgoeqppfTVg8VAPyjfwyUCgG/nc6Va53yiysw6bMDcg/ZGz+ewo6UHPkzH/yaiitF5YgO9MTGf/SFxt0VyVk6rKseat195gq2nsyGUiFg1fhbcGeHQFRUGfBGdSAqLKnE4q3GoDLzbx0x/Z5YAMDbW1KgLauEKIp45fuT0BtE/K1LCJY90QcKAdh8NBP7qoe11+xJw7krxQj0UuGbf/ZFuwAP5OjKsWT7GQDG0Lx2nzGQvfdILyT0CINBhByIyir1cg/ZhDui8dzgTgCMIUx62Ph7v5zB1eIKdAj2wv+m3A6VUoGdp3Ox9WS2zetJ5hioiKjZxAR5yYHJS+2CqFpfbp5qF7QPrOlB6BpuGSTM77jzsVhvuqxLuI/VO7qkZYJQc/edxDgXlEmBeH1tCLdsg2nYaevvAY27+Zd8uMYNfiZf/PWFNmvrzYOlletg0kYvtQuiaw2nmYa24+mFSC8oxfzvT+KmV7ZiyLu7rNbQWCMN54XVKkgHjAFYerafFLysMX3sTG2mPVQda/UEAua9ntbCkiAIZssjrdRQyZ8P8MDwm8IBAB/uSEWl3oB/rD2ES1dLEenvjhE9w2EQgae/OIzUHB3OXynGyt+NPVYvJ3RGsI8bnr23IwBg0dYU5BWVY/7mEwCA0bdHoWOIN15O6AKlQsDPJ7KwJ/UK3t1+GvkllegY4oXHbm2LsX3bISbIE1eKKvD+L2ew9WQ2dqdegcpFgTlDu6BLuA8evbUtAOA/359EtrYM71UHp38PjkOAlxrzhnUFAKzYfR6pOTrM3XQcBhF44KZwxMcE4KWhneHuqsSBtHxsOpKBFbvP40JeCYK91Xjqnlg8dmtbxIV6o7C0Eu/8chqpOTp8tjcNADA3oQtiQ7wx+a5oAMbeuNIK60XyZI6BioiajbtKKQ+L1C6SlkjhQSFAniTT2vrav0u61bPezVUp13JFB3rCS+1isY3ZMSIaHnZiQ7zgqhRsrjcWphuX+7i5WP2SNwtMVkKdabusHSPYx00e3ulSz7Wev/kk7nprB1b+fh7FFXqkZOvwxIp9mPzZQTns2GJrygRJuEkdVaXegL8uFWDd/ov47XQudGWVMBhqHgnTzkqgMh2yM60jk5hOo2Cr90mqo/L3VMl3W9ryz4ExEARgy4ls/OPzQ/jj3FV4qpRYMfYWvPX3m3BrO3/oyqswac1BzN10HJV6Ef07BmFgp2AAwOO3tUXHEC/kl1TioY/34nR2Efw8XDFjkLHnqVOoN564zRiIZm08Jg/lzU3oChelAq5KBeZWB6LVe9Iwb5MxkE3uFy3P0TXzb53g4+aCU5laPPTxXhSVV6FHhEZ+MPbAuGDcExeMKoOIx5bvw58XC+CpUuKloZ0BGP+s/nV3BwDA6z+ewge/Goc4X7y/M7zULlAqBMwd1gUA8PkfFzB93RFUGUQM6hyCuzoGAQCmDeyAcI0b0gtKsTSp4YX8N6JWEajS0tIwceJEREdHw93dHTExMZg3bx4qKmruYkhKSsLw4cMRFhYGT09P9OzZE2vXrrW5z3Xr1kEQBIwYMcLudvz+++9wcXFBz549m3A2RNe3jtVfgNaGsUyXdwj2grvJUE3t9bb2EexdEyS6Wuk9Mi6vngDTSlAx3a+fh6vZ8+nsbYPp3YRdrfQemX6udkG6RAo7rkoBHUMte2Y6BHlBXf0g4vqupbXAZbr+/JVi6A0i+sYEYPmYPhh/RzsoFQK2nczGvYt/w9T/HsLzG/7C3E3H8caPp/D2lmQs3JKChVtS8Hv13XVtrPRQATV3/r3x4yl0m7cFwz/8HbM2HsOYlfvR45WtGPzubyivMk64WXvaBcA8ZNWuwQOMw3zSVAu1C9AlUbWe7VeXDsHeuK9bKADgl1M5EATjMFnHEOOM/UufuBltfN2RlleCXWeuwEUh4OWEzvKfoYtSgbkJxkB0rro27Nm/dYKvR82Q7DP3doSvhysuXi2B3iDi3i4huDM2UF7fv2MQBnU2BqIsbRlCfNT454AO8np/TxWeqe4Jk8LovGFdzULzywldoFIqkFM9nDt9UKzcWwgAk/pFo12AB3J15Sit1OOWdn4Y3jNcXt83JhD3dQuFQTTeTapSKvByQmd5vYfKBS8nGEPXx7+dqzd4E2D5v25OKDk5GQaDAcuWLUOHDh1w/PhxTJ48GcXFxVi4cCEAYM+ePejRowdeeOEFhISEYPPmzRgzZgw0Gg0SEhLM9peWlobnnnsO/fr1s7sNBQUFGDNmDO655x5kZ3NMmciW0fFRyNGV49FbI62uf+CmcPx4LBOP3RZldX1ssBeG9wyHh8rF7AvC1D8GxOCn41n4W5dQm204lanFmHjrx7g7Lhjx7QNwd1yw1bDj76nC2PgoXCmukOewqm3KXe2xek8a/q9XG6vrH7klEnvP5mHindFW1/eO8sNDvSPQIdgLahfLYOmiVODpe2Jx5FIBbm8fYHUfk/u1R1FZFR6v7hGpbdhNYfjhaAYi/Dww5a72csC6t0sIHru1LeZvPoldZ67g5zrmJZJYm98JqAlEUk+Wxt0VXcN9cCm/BJeulspTU0QHesqPjDEV6uMGjbsrisqrrPZYAsaZzrO0ZWbDxabiqj8nzYhen38O6IAfjxnP+bm/dcKgLiHyugAvNT4d2wcjl+5BSYUeY+LbmfWSAcCdsYG4t0sItp3MRlyoNx69xfzvuq+HCs/e2xFzN52ASqnAS/d3Rm0vDe2CnadzUakXMeu+OHjW6kl94vYofLHvIs7kFOH/erVB7yg/s/XtAj0x5a72+GBHKjoEe2H8HeZ/z9QuSswd1gUTVh+EQgD+80BXi7/rL97fGduTc1BRZcDEftEWQ7JDuoXizg6B2J16BW/8eArLRvep67Le8ATRjurEm2++uWE7FQR89913aNPG+j80zeHtt9/G0qVLce7cOZvbDB06FCEhIVi5cqW8TK/X46677sKECROwa9cuFBQU4Ntvv633eI888ghiY2OhVCrx7bff4siRI3a3VavVQqPRoLCwED4+1v/BICK61kRRxN6zeUjJ1qG0Uo+ySgPKKvWoqDK/nT5M44ZJ/dpb1KsBxqLutfsuIMTHDTdH+SE6wFPuScnRluHQhXyczNRiQKcg9I7yt9qOI5cKUFhaif7Vw021nc7W4efjWZjUL9psTjFJpd6A7adycFu0v0Xxvi1f7LuI4vIqTOoXbTVU7z9/FduTs/H03bEWYQcAcnXl+HjnWTx6a1uLu0ABY/H/xzvPokOwN4Z0sx78f03OxuX8Uoy+PcpqG87mFuHbw+mY1K+9RZ0eYDzvrw5dxp0dAm323m04eAnebq422/D9XxnYf/4qZt8fZ/XapuYU4c2fkjFnaGe0sxFor2cN+f62K1ApFArMnDkTXl6Wf2lqE0URb775Jk6ePIn27dvb3+oGmjNnDn7++WccPHjQ5jZ33nknbr/9drkXCwDmzZuHo0eP4ptvvsG4cePsClSrVq3C0qVLsWfPHrz22mv1Bqry8nKUl9fcVaPVahEZGclARURErUvBJUCbDniFAP7We1uvZw0JVHYP+T3//PMIDg62a9tFixbZu9tGSU1NxZIlS8yCUm3r16/HgQMHsGzZMnnZ7t27sWLFigb1Lp05cwazZs3Crl274OJi3+VKTEzEK6+8YvcxiIiInNKeJcD+ZUC/mcA9cx3dGqdmV1H6+fPnERRkvSvWmpMnTyIqynrtgqlZs2ZBEIQ6X8nJyWafSU9Px5AhQ/DQQw9h8uTJVve7Y8cOjB8/HsuXL0fXrsbiQZ1Oh9GjR2P58uUIDAy0+rna9Ho9HnvsMbzyyivo2LGjXZ8BgNmzZ6OwsFB+Xbp0qf4PERERORtV9TBfRUnd25F9Q34tJTc3F3l5eXVu0759e6hUxjHxjIwMDBgwALfffjtWr14NhcIyD+7cuRNDhw7F4sWLMWXKFHn5kSNH0KtXLyiVNcWfBoOxTkChUCAlJQUxMTFm+yooKICfn5/FZ0RRhFKpxNatW3H33XfXe56soSIiolbpt7eBX18Deo0Ghn/g6NZccy0y5FdbSUkJLl68aDZ1AQD06NHD7n0EBQXZ3fOVnp6OgQMHonfv3li1apXVMJWUlISEhAQsWLDALEwBQFxcHI4dM3+o5Zw5c6DT6fDee+8hMtLyjiQfHx+Lz3z00Uf49ddf8dVXXyE6+sYbTyYiohuIqrp2upI9VPVpcKDKzc3F+PHj8dNPP1ldr9c3/4yq6enpGDBgAKKiorBw4ULk5ubK60JDjXcu7NixAwkJCZg+fTpGjhyJrCzjLbEqlQr+/v5wc3NDt27dzPbr6+sLAGbLZ8+ejfT0dHz22WdQKBQWnwkODra6LyIiouuOa/XdgxWch6o+DZ7Yc8aMGSgoKMC+ffvg7u6On3/+GWvWrEFsbCy+++67lmgjtm3bhtTUVGzfvh0REREICwuTX5I1a9agpKQEiYmJZusffPDBBh0rMzMTFy9erH9DIiKi651cQ8VAVZ8G11CFhYVh06ZNuPXWW+Hj44ODBw+iY8eO+O677/DWW29h9+7dLdXWVos1VERE1Cql/AR8+QgQfjMwZYejW3PNNeT7u8E9VMXFxfL0CX5+fvLwW/fu3fHnn382orlERETklKQeKtZQ1avBgapTp05ISUkBANx0001YtmwZ0tPT8fHHH5sNwREREVEr58ohP3s1uCh9+vTpyMzMBGCcdXzIkCFYu3YtVCoVVq9e3dztIyIiIkdhDZXdGhyonnjiCfn33r1748KFC0hOTkbbtm3tnjCTiIiIWgEV7/KzV6PnoZJ4eHg0+OHJRERE1ApI81DpywF9FaBscmy4btl9ZZ599lm7tlu8eHGjG0NERERORJqHCgAqiwGlxnFtcXJ2B6rDhw+bvd+9ezd69+4Nd3d3eZkgCM3XMiIiInIsFzUgKAFRb3yenxsDlS12B6odO8znn/D29sYXX3yB9u3bN3ujiIiIyAkIgrEwvVzLOqp6NHjaBCIiIrqByHNRMVDVhYGKiIiIbOPz/OzCQEVERES2yXNRcbb0uthdQ3X06FGz96IoIjk5GUVFRWbLe/To0TwtIyIiIseTA1VR3dvd4OwOVD179oQgCDB9lnJCQgIAyMsFQYBer2/+VhIREZFj8Hl+drE7UJ0/f74l20FERETOiDVUdrE7UEVFRbVkO4iIiMgZSbOlM1DVya6i9KNHj8JgMNi90xMnTqCqqqrRjSIiIiInwef52cWuQNWrVy/k5eXZvdP4+HhcvHix0Y0iIiIiJ8EaKrvYNeQniiJefvlleHh41L8xgIqKiiY1ioiIiJyEK+/ys4ddgequu+5CSkqK3TuNj483e8YfERERtVKch8oudgWqpKSkFm4GEREROSXWUNmFM6UTERGRbdJdfnyWX50YqIiIiMg2zkNlFwYqIiIiso01VHZhoCIiIiLb+Cw/uzQ4UBUXs8uPiIjohsF5qOzS4EAVEhKCCRMmYPfu3S3RHiIiInImrKGyS4MD1eeff46rV6/i7rvvRseOHfHmm28iIyOjJdpGREREjibf5VcCNOAxdDeaBgeqESNG4Ntvv0V6ejqmTp2KL774AlFRUUhISMDGjRv5DD8iIqLricrkKSkc9rOp0UXpQUFBePbZZ3H06FEsXrwYv/zyC/7+978jPDwcc+fORUkJLzoREVGr5+IOQDD+zkBlk10zpVuTnZ2NNWvWYPXq1bhw4QL+/ve/Y+LEibh8+TIWLFiAP/74A1u3bm3OthIREdG1plAY66gqi6vv9At2dIucUoMD1caNG7Fq1Sps2bIFXbp0wT//+U888cQT8PX1lbfp27cvOnfu3JztJCIiIkdReVYHKvZQ2dLgQDV+/Hg88sgj+P3333HLLbdY3SY8PBwvvfRSkxtHRERETkDlARSDd/rVocGBKjMzEx4eHnVu4+7ujnnz5jW6UURERORE+Dy/ejU4UFVVVUGr1VosFwQBarUaKpWqWRpGREREToJzUdWrwYHK19cXgiDYXB8REYFx48Zh3rx5UCj4ZBsiIqJWj8/zq1eDA9Xq1avx0ksvYdy4cbj11lsBAPv378eaNWswZ84c5ObmYuHChVCr1XjxxRebvcFERER0jfF5fvVqcKBas2YNFi1ahIcfflheNmzYMHTv3h3Lli3D9u3b0bZtW7z++usMVERERNcDPs+vXg0ek9uzZw969eplsbxXr17Yu3cvAODOO+/ExYsXm946IiIicjy5h4o1VLY0OFBFRkZixYoVFstXrFiByMhIAEBeXh78/Pya3joiIiJyPBal16vBQ34LFy7EQw89hJ9++kmeh+rgwYNITk7GV199BQA4cOAARo0a1bwtJSIiIseQpk1goLKpwYHqgQceQEpKCpYtW4aUlBQAwH333Ydvv/0W7dq1AwD84x//aNZGEhERkQNJD0hmDZVNDRryq6ysxD333IPKykokJiZi48aN2LhxIxITE+Uw1RLS0tIwceJEREdHw93dHTExMZg3bx4qKirkbZKSkjB8+HCEhYXB09MTPXv2xNq1a23uc926dRAEASNGjKj3+OXl5XjppZcQFRUFtVqNdu3aYeXKlc1xakRERM6Pd/nVq0E9VK6urjh69GhLtcWm5ORkGAwGLFu2DB06dMDx48cxefJkFBcXY+HChQCMxfI9evTACy+8gJCQEGzevBljxoyBRqNBQkKC2f7S0tLw3HPPoV+/fnYd/+GHH0Z2djZWrFiBDh06IDMzEwaDodnPk4iIyCm5ch6q+giiKIoN+cAzzzwDtVqNN998s6XaZJe3334bS5cuxblz52xuM3ToUISEhJj1Jun1etx1112YMGECdu3ahYKCAnz77bc29/Hzzz/jkUcewblz5+Dv79+otmq1Wmg0GhQWFsLHx6dR+yAiInKYk98B60cDkbcDE7c4ujXXTEO+vxv16JmVK1fil19+Qe/eveHp6Wm2fvHixQ3dZaMUFhbWG3AKCwvRuXNns2Xz589HcHAwJk6ciF27dtV7nO+++w59+vTBW2+9hf/+97/w9PTEAw88gFdffRXu7u5WP1NeXo7y8nL5vbVH9RAREbUacg0Vi9JtaXCgOn78OG6++WYAwOnTp83W1fVImuaUmpqKJUuWyMN91qxfvx4HDhzAsmXL5GW7d+/GihUrcOTIEbuPde7cOezevRtubm745ptvcOXKFfzzn/9EXl4eVq1aZfUziYmJeOWVV+w+BhERkVPjXX71anCg2rFjR7MdfNasWViwYEGd25w6dQpxcXHy+/T0dAwZMgQPPfQQJk+ebLON48ePx/Lly9G1a1cAgE6nw+jRo7F8+XIEBgba3UaDwQBBELB27VpoNBoAxl64v//97/joo4+s9lLNnj0bzz77rPxeq9XKc3QRERG1OvI8VKyhsqXBgUqSmpqKs2fP4q677oK7uztEUWxwD9XMmTMxbty4Ordp3769/HtGRgYGDhyIvn374pNPPrG6/c6dOzFs2DC88847GDNmjLz87NmzSEtLw7Bhw+RlUmG5i4sLUlJSEBMTY7G/sLAwtGnTRg5TANC5c2eIoojLly8jNjbW4jNqtRpqtbrO8yIiImo1OFN6vRocqPLy8vDwww9jx44dEAQBZ86cQfv27TFx4kT4+flh0aJFdu8rKCgIQUFBdm2bnp6OgQMHonfv3li1ahUUCssZH5KSkpCQkIAFCxZgypQpZuvi4uJw7Ngxs2Vz5syBTqfDe++9Z7MH6Y477sCGDRtQVFQELy9jl+fp06ehUCgQERFhV9uJiIhaNflZfsWAKALXqMSnNWnwo2eeeeYZuLq64uLFi/Dw8JCXjxo1Cj///HOzNk6Snp6OAQMGoG3btli4cCFyc3ORlZWFrKwseZsdO3Zg6NChePrppzFy5Eh5/dWrVwEAbm5u6Natm9nL19cX3t7e6NatG1QqFQDjcJ1pz9Zjjz2GgIAAjB8/HidPnsRvv/2G559/HhMmTLBZlE5ERHRdkQKVaACqyhzbFifV4B6qrVu3YsuWLRa9M7Gxsbhw4UKzNczUtm3bkJqaitTUVIvjSrM+rFmzBiUlJUhMTERiYqK8vn///khKSrL7WJmZmWYPdvby8sK2bdvw1FNPoU+fPggICMDDDz+M1157rWknRURE1Fq41nSgoKIYcGWHQm0NnofK29sbf/75J2JjY+Ht7Y2//voL7du3x8GDBzF48GDk5eW1VFtbLc5DRURErd5roUBVKTD9KOAX5ejWXBMN+f5u8JBfv3798Nlnn8nvBUGAwWDAW2+9hYEDBza8tUREROT8pLmoWJhuVYOH/N566y3cc889OHjwICoqKvDvf/8bJ06cwNWrV/H777+3RBuJiIjI0VSeQEkeH5BsQ4N7qLp164bTp0/jzjvvxPDhw1FcXIwHH3wQhw8ftjrtABEREV0HXPmA5Lo0ah4qjUaDl156qbnbQkRERM5KxQck16VRgaqgoAD79+9HTk6OPDmmxHTKASIiIrpOsIaqTg0OVN9//z0ef/xxFBUVwcfHx2x2dEEQGKiIiIiuR9Lz/PiAZKsaXEM1c+ZMTJgwAUVFRSgoKEB+fr78kibRJCIiouuMK3uo6tLgQJWeno6nn37abJZ0IiIius6xhqpODQ5UgwcPxsGDB1uiLUREROSsVLzLry4NrqEaOnQonn/+eZw8eRLdu3eHq6ur2foHHnig2RpHRERETkJ+QDJ7qKxpcKCaPHkyAGD+/PkW6wRBgF6vb3qriIiIyLmwhqpODQ5UtadJICIiohuAdJcfA5VVDa6hIiIiohsQ56Gqk92B6v7770dhYaH8/s0330RBQYH8Pi8vD126dGnWxhEREZGTYA1VnewOVFu2bEF5ebn8/o033jCbd6qqqgopKSnN2zoiIiJyDnyWX53sDlSiKNb5noiIiK5jnIeqTqyhIiIiovqxhqpOdgcqQRDMntsnLSMiIqIbAJ/lVye7p00QRRHjxo2DWq0GAJSVlWHq1Knw9DR2AZrWVxEREdF1hvNQ1cnuQDV27Fiz90888YTFNmPGjGl6i4iIiMj5SDVUhiqgqgJwUTm2PU7G7kC1atWqlmwHEREROTMpUAHGO/1c/B3XFifEonQiIiKqn9IVUFb3SnEuKgsMVERERGQf1lHZxEBFRERE9uHz/GxioCIiIiL7cC4qmxioiIiIyD58np9NDFRERERkHz7PzyYGKiIiIrIPn+dnEwMVERER2Yc1VDYxUBEREZF95BoqBqraGKiIiIjIPnINFQNVbQxUREREZB/WUNnEQEVERET2kWuoeJdfbQxUREREZB/OlG4TAxURERHZxyPQ+LM417HtcEIMVERERGQf71DjT12mY9vhhBioiIiIyD4+4cafuizHtsMJMVARERGRfbxCjD8rioBynWPb4mQYqIiIiMg+ai9A7WP8XcthP1MMVERERGQ/1lFZxUBFRERE9vMOM/5kHZWZVhGo0tLSMHHiRERHR8Pd3R0xMTGYN28eKioq5G2SkpIwfPhwhIWFwdPTEz179sTatWtt7nPdunUQBAEjRoyo9/hr167FTTfdBA8PD4SFhWHChAnIy8trjlMjIiJqXeRAleHYdjiZVhGokpOTYTAYsGzZMpw4cQLvvPMOPv74Y7z44ovyNnv27EGPHj3w9ddf4+jRoxg/fjzGjBmDzZs3W+wvLS0Nzz33HPr161fvsX///XeMGTMGEydOxIkTJ7Bhwwbs378fkydPbtZzJCIiahXkIT/2UJkSRFEUHd2Ixnj77bexdOlSnDt3zuY2Q4cORUhICFauXCkv0+v1uOuuuzBhwgTs2rULBQUF+Pbbb23uY+HChVi6dCnOnj0rL1uyZAkWLFiAy5cv29VWrVYLjUaDwsJC+Pj42PUZIiIip7RvGfDTv4Euw4GHP3N0a1pUQ76/W0UPlTWFhYXw9/dv8Dbz589HcHAwJk6caNdx4uPjcenSJfz4448QRRHZ2dn46quvcP/999v8THl5ObRardmLiIjousAeKqtaZaBKTU3FkiVL8OSTT9rcZv369Thw4ADGjx8vL9u9ezdWrFiB5cuX232sO+64A2vXrsWoUaOgUqkQGhoKjUaDDz/80OZnEhMTodFo5FdkZKTdxyMiInJqUg0Vp00w49BANWvWLAiCUOcrOTnZ7DPp6ekYMmQIHnroIZt1TDt27MD48eOxfPlydO3aFQCg0+kwevRoLF++HIGBgXa38eTJk5g+fTrmzp2LQ4cO4eeff0ZaWhqmTp1q8zOzZ89GYWGh/Lp06ZLdxyMiInJqclF6JtA6q4ZahENrqHJzc+u9W659+/ZQqVQAgIyMDAwYMAC33347Vq9eDYXCMg/u3LkTQ4cOxeLFizFlyhR5+ZEjR9CrVy8olUp5mcFgAAAoFAqkpKQgJibGYn+jR49GWVkZNmzYIC/bvXs3+vXrh4yMDISFhdV7nqyhIiKi60ZVBfBakPH3588BngGObU8Lasj3t8s1apNVQUFBCAoKsmvb9PR0DBw4EL1798aqVaushqmkpCQkJCRgwYIFZmEKAOLi4nDs2DGzZXPmzIFOp8N7771nc1iupKQELi7ml0kKZa20np+IiKjxXFSARyBQcsU4dcJ1HKgawqGByl7p6ekYMGAAoqKisHDhQuTm5srrQkONxXE7duxAQkICpk+fjpEjRyIry1gsp1Kp4O/vDzc3N3Tr1s1sv76+vgBgtnz27NlIT0/HZ58Z71wYNmwYJk+ejKVLl2Lw4MHIzMzEjBkzcOuttyI8PLwlT5uIiMg5+YRVB6osILS7o1vjFFpFoNq2bRtSU1ORmpqKiIgIs3VSL9GaNWtQUlKCxMREJCYmyuv79++PpKQku4+VmZmJixcvyu/HjRsHnU6HDz74ADNnzoSvry/uvvtuLFiwoGknRURE1Fp5hwFZx/j4GROtdh6q1oQ1VEREdF357ingz8+AAS8CA15wdGtazA0xDxURERE5iHd1yQt7qGQMVERERNQwnNzTAgMVERERNQwfkGyBgYqIiIgahj1UFhioiIiIqGF8qmuoinIAfZVj2+IkGKiIiIioYTwCAUEJQASKsh3dGqfAQEVEREQNo1Bw2K8WBioiIiJqONOHJBMDFRERETWC3EPFQAUwUBEREVFjsIfKDAMVERERNZyPFKhYQwUwUBEREVFjSD1UWk7uCTBQERERUWPwLj8zDFRERETUcHxAshkGKiIiImo4qYeqrACoLHVoU5wBAxURERE1nJsGcHE3/s5eKgYqIiIiagRB4J1+JhioiIiIqHF4p5+MgYqIiIgah3f6yRioiIiIqHE4W7qMgYqIiIgah4FKxkBFREREjSMN+WkZqBioiIiIqHECYow/s08ABr1j2+JgDFRERETUOCHdAZU3UF5oDFU3MAYqIiIiahylC9D2NuPvF353bFscjIGKiIiIGi+qr/EnAxURERFRI0Xdafx5YQ8gio5tiwMxUBEREVHjhfcyPtOvJA/ITXF0axyGgYqIiIgaz0UFRN5i/P3Cbse2xYEYqIiIiKhpou4w/rywx7HtcCAGKiIiImoa00B1g9ZRMVARERFR00T0ARSuxkfQXD3n6NY4BAMVERERNY2rO9Cmt/H3G3TYj4GKiIiImq6dNOx3Y85HxUBFRERETXeDT/DJQEVERERNF3kbICiBgotAwSVHt+aaY6AiIiKiplN7A2E3GX+/AeuoXBzdACIiIrpOtLsDyPgTSN4MKJRA9nEg+yQQGAv87TVAEBzdwhbDQEVERETNI+oOYM8S4NR3xpfkzBbg5jFAUCfHta2FcciPiIiImke7foB/DODqCUTcAvQeDwTEGtdd2u/YtrUw9lARERFR81B7AU8dMs6Wrqjus9k2D/j9XeDSPuDm0Q5tXktqFT1UaWlpmDhxIqKjo+Hu7o6YmBjMmzcPFRUV8jZJSUkYPnw4wsLC4OnpiZ49e2Lt2rVm+1m9ejUEQTB7ubm51Xv8pKQk3HzzzVCr1ejQoQNWr17d3KdIRER0fRCEmjAFGO/+A4DLBxzTnmukVfRQJScnw2AwYNmyZejQoQOOHz+OyZMno7i4GAsXLgQA7NmzBz169MALL7yAkJAQbN68GWPGjIFGo0FCQoK8Lx8fH6SkpMjvhXoK5M6fP4+hQ4di6tSpWLt2LbZv345JkyYhLCwMgwcPbpkTJiIiul5E3GL8mZsMlOYD7n6ObU8LEUSxdT7F8O2338bSpUtx7pztZwYNHToUISEhWLlyJQBjD9WMGTNQUFBg93FeeOEF/PDDDzh+/Li87JFHHkFBQQF+/vlnu/ah1Wqh0WhQWFgIHx8fu49NRER0XXi/l/EZf49/DcQOcnRr7NaQ7+9WMeRnTWFhIfz9/Ru8TVFREaKiohAZGYnhw4fjxIkTde5j7969GDTI/A9/8ODB2Lt3r83PlJeXQ6vVmr2IiIhuWBG3Gn9e2ufYdrSgVhmoUlNTsWTJEjz55JM2t1m/fj0OHDiA8ePHy8s6deqElStXYtOmTfj8889hMBjQt29fXL582eZ+srKyEBISYrYsJCQEWq0WpaWlVj+TmJgIjUYjvyIjIxt4hkRERNeRyOpAdfn6vdPPoYFq1qxZFkXitV/Jyclmn0lPT8eQIUPw0EMPYfLkyVb3u2PHDowfPx7Lly9H165d5eXx8fEYM2YMevbsif79+2Pjxo0ICgrCsmXLmvW8Zs+ejcLCQvl16dKNNwU/ERGRTA5UBwGD3rFtaSEOLUqfOXMmxo0bV+c27du3l3/PyMjAwIED0bdvX3zyySdWt9+5cyeGDRuGd955B2PGjKlz366urujVqxdSU1NtbhMaGors7GyzZdnZ2fDx8YG7u7vVz6jVaqjV6jqPTUREdMMI7gKovICKIiDnFBDazdEtanYODVRBQUEICgqya9v09HQMHDgQvXv3xqpVq6BQWHauJSUlISEhAQsWLMCUKVPq3ader8exY8dw//3329wmPj4eP/74o9mybdu2IT4+3q52ExER3fAUSqBNb+D8TmMd1XUYqFpFDVV6ejoGDBiAtm3bYuHChcjNzUVWVhaysrLkbXbs2IGhQ4fi6aefxsiRI+X1V69elbeZP38+tm7dinPnzuHPP//EE088gQsXLmDSpEnyNrNnzzbr2Zo6dSrOnTuHf//730hOTsZHH32E9evX45lnnrk2J09ERHQ9uM7no2oV81Bt27YNqampSE1NRUREhNk6adaHNWvWoKSkBImJiUhMTJTX9+/fH0lJSQCA/Px8TJ48GVlZWfDz80Pv3r2xZ88edOnSRd4+MzMTFy9elN9HR0fjhx9+wDPPPIP33nsPERER+PTTTzkHFRERUUNEXt93+rXaeahaE85DRUREN7zSfGBBO+Pvz58FPAMd2hx73BDzUBEREVEr4u4HBMUZf78OH5TMQEVERETXhvQYmutwPioGKiIiIro2pMJ09lARERERNZJUmJ7+J6CvdGxbmhkDFREREV0bAbGAmy9QVQpk1/0s3daGgYqIiIiuDYUCCKl+JFxuimPb0swYqIiIiOjaCexo/HmFgYqIiIiocYI6GX9eOe3YdjQzBioiIiK6dgJjjT9zGaiIiIiIGiewuofq6tnr6k4/BioiIiK6dnzaAK6egKEKyE9zdGuaDQMVERERXTsKBRDYwfj7dXSnHwMVERERXVvSsN91dKcfAxURERFdW0HS1AlnHNuOZsRARURERNeW1EPFIT8iIiKiRgo06aESRce2pZkwUBEREdG15d8eEJRAhQ7QZTq6Nc2CgYqIiIiuLReVMVQB182wHwMVERERXXvysN/1MWM6AxURERFde9KdfuyhIiIiImqkwOvrIckMVERERHTtcciPiIiIqIkCY40/i7KB0gKHNqU5MFARERHRtefmA3iHG3+/DmZMZ6AiIiIix5AfQdP6C9MZqIiIiMgxAq+fO/0YqIiIiMgxAq+fhyQzUBEREZFjBElTJzSgh8qgB/YvB3KSW6ZNjcRARURERI4h9VDlpwGVZfZ95q91wI/PARsnt1izGoOBioiIiBzDKwRQawDRAFw9a99njn9l/Jl1FMi/0HJtayAGKiIiInIMQai50y/t9/q3L74CnNtZ8/70zy3TrkZgoCIiIiLH6fZ348/f3gbKi+re9tT3gKiveZ/8Q8u1q4EYqIiIiMhx+kwA/KKB4hxgz/t1b3tio/HnzWONPy/87jSzrDNQERERkeO4qIB7XzH+/vv7gDbD+nZFOUDabuPv/WYCQXGAoQo4s+3atLMeDFRERETkWJ0fACJvA6pKgR2vW9/m5CZj8Xqb3oBfFNDpfuPyFOcY9mOgIiIiIscSBOBv1UHq8Fog67jlNie+Mf7s+qDxZ9xQ488zvwBV5S3fxnowUBEREZHjRd4CdP0/ACKw7WXzddpM4MIe4+9dRxh/ht8MeIUCFTogbde1bKlVDFRERETkHO6ZByhcgbO/Ar+/Z5wVHTAO90E0DgtqIozLFAqg0xDj78k/OqS5phioiIiIyDn4RwN9/2X8fdtcYPndQMbhmrv7pOE+iVxH9RMgiteunVYwUBEREZHzuHsuMHSxcQb1zCPGUHVpHwAB6DLcfNvo/oCrJ6DLMAYvB2KgIiIiIuehUAC3TAT+dQDo/pDxzj4AiOoL+ISZb+vqBnS42/h7imOH/VpFoEpLS8PEiRMRHR0Nd3d3xMTEYN68eaioqJC3SUpKwvDhwxEWFgZPT0/07NkTa9euNdvP6tWrIQiC2cvNza3OY2/cuBH33nsvgoKC4OPjg/j4eGzZsqVFzpOIiIiqeYcAIz8FRn9jDFaD37C+Xafqu/1Sfrp2bbPCxaFHt1NycjIMBgOWLVuGDh064Pjx45g8eTKKi4uxcOFCAMCePXvQo0cPvPDCCwgJCcHmzZsxZswYaDQaJCQkyPvy8fFBSkqK/F4QhDqP/dtvv+Hee+/FG2+8AV9fX6xatQrDhg3Dvn370KtXr5Y5YSIiIjKKudv4sqXjYEBQAhCMj65Re12zppkSRNHBVVyN9Pbbb2Pp0qU4d+6czW2GDh2KkJAQrFy5EoCxh2rGjBkoKCho0rG7du2KUaNGYe7cuXZtr9VqodFoUFhYCB8fnyYdm4iIiGopygG8gpt9tw35/m4VQ37WFBYWwt/fv8HbFBUVISoqCpGRkRg+fDhOnDjRoOMaDAbodLo6j11eXg6tVmv2IiIiohbSAmGqoVploEpNTcWSJUvw5JNP2txm/fr1OHDgAMaPHy8v69SpE1auXIlNmzbh888/h8FgQN++fXH58mW7j71w4UIUFRXh4YcftrlNYmIiNBqN/IqMjLR7/0RERNT6OHTIb9asWViwYEGd25w6dQpxcXHy+/T0dPTv3x8DBgzAp59+avUzO3bsQEJCApYuXYoxY8bY3HdlZSU6d+6MRx99FK+++mq97f3iiy8wefJkbNq0CYMGDbK5XXl5OcrLa6bB12q1iIyM5JAfERFRK9KQIT+HFqXPnDkT48aNq3Ob9u3by79nZGRg4MCB6Nu3Lz755BOr2+/cuRPDhg3DO++8U2eYAgBXV1f06tULqamp9bZ13bp1mDRpEjZs2FBnmAIAtVoNtVpd7z6JiIjo+uDQQBUUFISgoCC7tk1PT8fAgQPRu3dvrFq1CgqF5WhlUlISEhISsGDBAkyZMqXefer1ehw7dgz3339/ndt9+eWXmDBhAtatW4ehQ4fa1V4iIiK6cbSKaRPS09MxYMAAREVFYeHChcjNzZXXhYaGAqgZ5ps+fTpGjhyJrKwsAIBKpZILyOfPn4/bb78dHTp0QEFBAd5++21cuHABkyZNkvc3e/ZspKen47PPPgNgHOYbO3Ys3nvvPdx2223yft3d3aHRaK7J+RMREZFzaxVF6du2bUNqaiq2b9+OiIgIhIWFyS/JmjVrUFJSgsTERLP1Dz5Y89yf/Px8TJ48GZ07d8b9998PrVaLPXv2oEuXLvI2mZmZuHjxovz+k08+QVVVFaZNm2a23+nTp1+bkyciIiKn12rnoWpNOA8VERFR63NDzENFRERE5CwYqIiIiIiaiIGKiIiIqIkYqIiIiIiaiIGKiIiIqIkYqIiIiIiaqFVM7NnaSTNTaLVaB7eEiIiI7CV9b9szwxQD1TWg0+kAAJGRkQ5uCRERETWUTqer9+konNjzGjAYDMjIyIC3tzcEQWjWfWu1WkRGRuLSpUucNLSJeC2bB69j8+G1bD68ls3nRrqWoihCp9MhPDzc6jOETbGH6hpQKBSIiIho0WP4+Phc93+xrxVey+bB69h8eC2bD69l87lRrqW9z+1lUToRERFREzFQERERETURA1Urp1arMW/ePKjVakc3pdXjtWwevI7Nh9ey+fBaNh9eS+tYlE5ERETUROyhIiIiImoiBioiIiKiJmKgIiIiImoiBioiIiKiJmKgasU+/PBDtGvXDm5ubrjtttuwf/9+RzfJ6SUmJuKWW26Bt7c3goODMWLECKSkpJhtU1ZWhmnTpiEgIABeXl4YOXIksrOzHdTi1uHNN9+EIAiYMWOGvIzX0X7p6el44oknEBAQAHd3d3Tv3h0HDx6U14uiiLlz5yIsLAzu7u4YNGgQzpw548AWOye9Xo+XX34Z0dHRcHd3R0xMDF599VWz57DxWlr322+/YdiwYQgPD4cgCPj222/N1ttz3a5evYrHH38cPj4+8PX1xcSJE1FUVHQNz8KxGKhaqf/973949tlnMW/ePPz555+46aabMHjwYOTk5Di6aU5t586dmDZtGv744w9s27YNlZWV+Nvf/obi4mJ5m2eeeQbff/89NmzYgJ07dyIjIwMPPvigA1vt3A4cOIBly5ahR48eZst5He2Tn5+PO+64A66urvjpp59w8uRJLFq0CH5+fvI2b731Ft5//318/PHH2LdvHzw9PTF48GCUlZU5sOXOZ8GCBVi6dCk++OADnDp1CgsWLMBbb72FJUuWyNvwWlpXXFyMm266CR9++KHV9fZct8cffxwnTpzAtm3bsHnzZvz222+YMmXKtToFxxOpVbr11lvFadOmye/1er0YHh4uJiYmOrBVrU9OTo4IQNy5c6coiqJYUFAgurq6ihs2bJC3OXXqlAhA3Lt3r6Oa6bR0Op0YGxsrbtu2Tezfv784ffp0URR5HRvihRdeEO+8806b6w0GgxgaGiq+/fbb8rKCggJRrVaLX3755bVoYqsxdOhQccKECWbLHnzwQfHxxx8XRZHX0l4AxG+++UZ+b891O3nypAhAPHDggLzNTz/9JAqCIKanp1+ztjsSe6haoYqKChw6dAiDBg2SlykUCgwaNAh79+51YMtan8LCQgCAv78/AODQoUOorKw0u7ZxcXFo27Ytr60V06ZNw9ChQ82uF8Dr2BDfffcd+vTpg4ceegjBwcHo1asXli9fLq8/f/48srKyzK6lRqPBbbfdxmtZS9++fbF9+3acPn0aAPDXX39h9+7duO+++wDwWjaWPddt79698PX1RZ8+feRtBg0aBIVCgX379l3zNjsCH47cCl25cgV6vR4hISFmy0NCQpCcnOygVrU+BoMBM2bMwB133IFu3boBALKysqBSqeDr62u2bUhICLKyshzQSue1bt06/Pnnnzhw4IDFOl5H+507dw5Lly7Fs88+ixdffBEHDhzA008/DZVKhbFjx8rXy9p/77yW5mbNmgWtVou4uDgolUro9Xq8/vrrePzxxwGA17KR7LluWVlZCA4ONlvv4uICf3//G+baMlDRDWvatGk4fvw4du/e7eimtDqXLl3C9OnTsW3bNri5uTm6Oa2awWBAnz598MYbbwAAevXqhePHj+Pjjz/G2LFjHdy61mX9+vVYu3YtvvjiC3Tt2hVHjhzBjBkzEB4ezmtJLY5Dfq1QYGAglEqlxR1T2dnZCA0NdVCrWpd//etf2Lx5M3bs2IGIiAh5eWhoKCoqKlBQUGC2Pa+tuUOHDiEnJwc333wzXFxc4OLigp07d+L999+Hi4sLQkJCeB3tFBYWhi5dupgt69y5My5evAgA8vXif+/1e/755zFr1iw88sgj6N69O0aPHo1nnnkGiYmJAHgtG8ue6xYaGmpxU1RVVRWuXr16w1xbBqpWSKVSoXfv3ti+fbu8zGAwYPv27YiPj3dgy5yfKIr417/+hW+++Qa//voroqOjzdb37t0brq6uZtc2JSUFFy9e5LU1cc899+DYsWM4cuSI/OrTpw8ef/xx+XdeR/vccccdFlN3nD59GlFRUQCA6OhohIaGml1LrVaLffv28VrWUlJSAoXC/GtNqVTCYDAA4LVsLHuuW3x8PAoKCnDo0CF5m19//RUGgwG33XbbNW+zQzi6Kp4aZ926daJarRZXr14tnjx5UpwyZYro6+srZmVlObppTu0f//iHqNFoxKSkJDEzM1N+lZSUyNtMnTpVbNu2rfjrr7+KBw8eFOPj48X4+HgHtrp1ML3LTxR5He21f/9+0cXFRXz99dfFM2fOiGvXrhU9PDzEzz//XN7mzTffFH19fcVNmzaJR48eFYcPHy5GR0eLpaWlDmy58xk7dqzYpk0bcfPmzeL58+fFjRs3ioGBgeK///1veRteS+t0Op14+PBh8fDhwyIAcfHixeLhw4fFCxcuiKJo33UbMmSI2KtXL3Hfvn3i7t27xdjYWPHRRx911CldcwxUrdiSJUvEtm3biiqVSrz11lvFP/74w9FNcnoArL5WrVolb1NaWir+85//FP38/EQPDw/x//7v/8TMzEzHNbqVqB2oeB3t9/3334vdunUT1Wq1GBcXJ37yySdm6w0Gg/jyyy+LISEholqtFu+55x4xJSXFQa11XlqtVpw+fbrYtm1b0c3NTWzfvr340ksvieXl5fI2vJbW7dixw+q/jWPHjhVF0b7rlpeXJz766KOil5eX6OPjI44fP17U6XQOOBvHEETRZApZIiIiImow1lARERERNREDFREREVETMVARERERNREDFREREVETMVARERERNREDFREREVETMVARERERNREDFREREVETMVARUavwn//8Bz179mzSPtLS0iAIAo4cOdIsbbJlwIABmDFjRoseg4icCwMVETWLS5cuYcKECQgPD4dKpUJUVBSmT5+OvLy8Bu9LEAR8++23Zsuee+45s4ezNkZkZCQyMzPRrVu3Ju1HkpSUBEEQUFBQYLZ848aNePXVV5vlGI1xrYIjEdVgoCKiJjt37hz69OmDM2fO4Msvv0Rqaio+/vhjbN++HfHx8bh69WqTj+Hl5YWAgIAm7UOpVCI0NBQuLi5Nbk9d/P394e3t3aLHICLnwkBFRE02bdo0qFQqbN26Ff3790fbtm1x33334ZdffkF6ejpeeukledt27drh1VdfxaOPPgpPT0+0adMGH374odl6APi///s/CIIgv6895Ddu3DiMGDECb7zxBkJCQuDr64v58+ejqqoKzz//PPz9/REREYFVq1bJn6ndczNu3DgIgmDxSkpKAgD897//RZ8+feDt7Y3Q0FA89thjyMnJkfc1cOBAAICfnx8EQcC4ceMAWA755efnY8yYMfDz84OHhwfuu+8+nDlzRl6/evVq+Pr6YsuWLejcuTO8vLwwZMgQZGZm2rzm+fn5ePzxxxEUFAR3d3fExsbK5xodHQ0A6NWrFwRBwIABA+TPffrpp+jcuTPc3NwQFxeHjz76yOL6rFu3Dn379oWbmxu6deuGnTt32nVcohuao5/OTEStW15enigIgvjGG29YXT958mTRz89PNBgMoiiKYlRUlOjt7S0mJiaKKSkp4vvvvy8qlUpx69atoiiKYk5OjghAXLVqlZiZmSnm5OSIoiiK8+bNE2+66SZ5v2PHjhW9vb3FadOmicnJyeKKFStEAOLgwYPF119/XTx9+rT46quviq6uruKlS5dEURTF8+fPiwDEw4cPi6IoigUFBWJmZqb8mj59uhgcHCxmZmaKoiiKK1asEH/88Ufx7Nmz4t69e8X4+HjxvvvuE0VRFKuqqsSvv/5aBCCmpKSImZmZYkFBgSiKoti/f39x+vTpclsfeOABsXPnzuJvv/0mHjlyRBw8eLDYoUMHsaKiQhRFUVy1apXo6uoqDho0SDxw4IB46NAhsXPnzuJjjz1m87pPmzZN7Nmzp3jgwAHx/Pnz4rZt28TvvvtOFEVR3L9/vwhA/OWXX8TMzEwxLy9PFEVR/Pzzz8WwsDDx66+/Fs+dOyd+/fXXor+/v7h69Wqz6xMRESF+9dVX4smTJ8VJkyaJ3t7e4pUrV+o9LtGNjIGKiJrkjz/+EAGI33zzjdX1ixcvFgGI2dnZoigaA9WQIUPMthk1apQcVERRtLo/a4EqKipK1Ov18rJOnTqJ/fr1k99XVVWJnp6e4pdffimKomWgMvX111+Lbm5u4u7du22e64EDB0QAok6nE0VRFHfs2CECEPPz8822Mw1Up0+fFgGIv//+u7z+ypUroru7u7h+/XpRFI2BCoCYmpoqb/Phhx+KISEhNtsybNgwcfz48VbX2TrPmJgY8YsvvjBb9uqrr4rx8fFmn3vzzTfl9ZWVlWJERIS4YMGCeo9LdCPjkB8RNQtRFO3eNj4+3uL9qVOnGnzMrl27QqGo+WcsJCQE3bt3l98rlUoEBATIw3S2HD58GKNHj8YHH3yAO+64Q15+6NAhDBs2DG3btoW3tzf69+8PALh48aLdbTx16hRcXFxw2223ycsCAgLQqVMns3P28PBATEyM/D4sLKzOdv/jH//AunXr0LNnT/z73//Gnj176mxHcXExzp49i4kTJ8LLy0t+vfbaazh79qzZtqZ/Pi4uLujTp4/c1oYel+hGwUBFRE3SoUMHCIJgMxCdOnUKfn5+CAoKavZju7q6mr0XBMHqMoPBYHMfWVlZeOCBBzBp0iRMnDhRXl5cXIzBgwfDx8cHa9euxYEDB/DNN98AACoqKprxLIystbuukHrffffhwoULeOaZZ5CRkYF77rkHzz33nM3ti4qKAADLly/HkSNH5Nfx48fxxx9/2N3Ohh6X6EbBQEVETRIQEIB7770XH330EUpLS83WZWVlYe3atRg1ahQEQZCX1/4C/+OPP9C5c2f5vaurK/R6fcs2HEBZWRmGDx+OuLg4LF682GxdcnIy8vLy8Oabb6Jfv36Ii4uz6DFSqVQAUGdbO3fujKqqKuzbt09elpeXh5SUFHTp0qVJ7Q8KCsLYsWPx+eef491338Unn3xis10hISEIDw/HuXPn0KFDB7OXVMQuMf3zqaqqwqFDh8z+fGwdl+hG1rL3DhPRDeGDDz5A3759MXjwYLz22muIjo7GiRMn8Pzzz6NNmzZ4/fXXzbb//fff8dZbb2HEiBHYtm0bNmzYgB9++EFe365dO2zfvh133HEH1Go1/Pz8WqTdTz75JC5duoTt27cjNzdXXu7v74+2bdtCpVJhyZIlmDp1Ko4fP24xt1RUVBQEQcDmzZtx//33w93dHV5eXmbbxMbGYvjw4Zg8eTKWLVsGb29vzJo1C23atMHw4cMb3fa5c+eid+/e6Nq1K8rLy7F582Y59AQHB8Pd3R0///wzIiIi4ObmBo1Gg1deeQVPP/00NBoNhgwZgvLychw8eBD5+fl49tln5X1/+OGHiI2NRefOnfHOO+8gPz8fEyZMqPe4RDcy9lARUZPFxsbi4MGDaN++PR5++GHExMRgypQpGDhwIPbu3Qt/f3+z7WfOnImDBw+iV69eeO2117B48WIMHjxYXr9o0SJs27YNkZGR6NWrV4u1e+fOncjMzESXLl0QFhYmv/bs2YOgoCCsXr0aGzZsQJcuXfDmm29i4cKFZp9v06YNXnnlFcyaNQshISH417/+ZfU4q1atQu/evZGQkID4+HiIoogff/zRYpivIVQqFWbPno0ePXrgrrvuglKpxLp16wAY657ef/99LFu2DOHh4XJwmzRpEj799FOsWrUK3bt3R//+/bF69WqLHqo333wTb775Jm666Sbs3r0b3333HQIDA+s9LtGNTBAbUklKRNRE7dq1w4wZM/hoFieUlpaG6OhoHD58uMmP+SG60bCHioiIiKiJGKiIiIiImohDfkRERERNxB4qIiIioiZioCIiIiJqIgYqIiIioiZioCIiIiJqIgYqIiIioiZioCIiIiJqIgYqIiIioiZioCIiIiJqov8HsFAj6RBpY2AAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Plot the energies.\n", + "\n", + "vqe_y = vqe_energies\n", + "vqe_x = list(range(len(vqe_y)))\n", + "plt.plot(vqe_x, vqe_y, label=\"VQE\")\n", + "\n", + "afqmc_y = list(qmc_data[\"ETotal\"])\n", + "afqmc_x = [i + vqe_x[-1] for i in list(range(len(afqmc_y)))]\n", + "plt.plot(afqmc_x, afqmc_y, label=\"AFQMC\")\n", + "\n", + "plt.xlabel(\"Optimization steps\")\n", + "plt.ylabel(\"Energy [Ha]\")\n", + "plt.legend()\n" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CUDA-Q Version (https://github.com/NVIDIA/cuda-quantum 9c4acf80526257f03f567668a3b42822451cfb7e)\n" + ] + } + ], + "source": [ + "print(cudaq.__version__)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/sphinx/examples/python/tutorials/images/qaoa-circuit-layers.png b/docs/sphinx/examples/python/tutorials/images/qaoa-circuit-layers.png index 3fe20e0cee..52671468c6 100644 Binary files a/docs/sphinx/examples/python/tutorials/images/qaoa-circuit-layers.png and b/docs/sphinx/examples/python/tutorials/images/qaoa-circuit-layers.png differ diff --git a/docs/sphinx/examples/python/tutorials/images/qaoa-problem-kernel.png b/docs/sphinx/examples/python/tutorials/images/qaoa-problem-kernel.png index 00406a6007..c6ccd197df 100644 Binary files a/docs/sphinx/examples/python/tutorials/images/qaoa-problem-kernel.png and b/docs/sphinx/examples/python/tutorials/images/qaoa-problem-kernel.png differ diff --git a/docs/sphinx/examples/python/tutorials/src/geo_fenta.xyz b/docs/sphinx/examples/python/tutorials/src/geo_fenta.xyz new file mode 100644 index 0000000000..a8ef211dfc --- /dev/null +++ b/docs/sphinx/examples/python/tutorials/src/geo_fenta.xyz @@ -0,0 +1,28 @@ +26 +Energy = -2155.935897674 +C -1.0290441 -0.9471117 -1.9865887 +C -2.2505436 -0.7922495 0.1582250 +C -0.2593185 -2.2544373 0.0162178 +O -1.3454085 0.8783278 -3.5342318 +O -1.1347371 1.1483476 1.0151618 +O 0.8453748 -2.7999130 2.0922098 +O -0.2780587 1.3040005 -1.6070564 +O -3.2498211 0.8000765 1.6534429 +C -0.9188318 0.5115628 -2.4417824 +C -2.2565356 0.4618854 1.0159078 +C 0.5342286 -1.9274094 1.2848894 +N -0.9221422 -1.0235290 -0.5013883 +Fe 0.2526900 0.4497192 0.0039380 +O 1.8221685 -0.2815654 -1.0936915 +H -1.9555949 -1.4065777 -2.3473606 +H -0.1824638 -1.4867161 -2.4264081 +H -2.5130292 -1.6526686 0.7833244 +H -3.0320380 -0.6959329 -0.6040366 +H 0.4617231 -2.6044401 -0.7316186 +H -0.9840969 -3.0551564 0.1969690 +H 2.1227503 0.3853924 -1.7373941 +H 2.5891208 -0.4447721 -0.5150575 +O 0.9052095 -0.6681693 1.3941789 +O 1.4573655 1.9796179 0.5273073 +H 1.0868407 2.8424014 0.2685815 +H 1.5840897 2.0352695 1.4913642 \ No newline at end of file diff --git a/docs/sphinx/examples/python/tutorials/src/geo_o3.xyz b/docs/sphinx/examples/python/tutorials/src/geo_o3.xyz new file mode 100644 index 0000000000..a0c771460c --- /dev/null +++ b/docs/sphinx/examples/python/tutorials/src/geo_o3.xyz @@ -0,0 +1,5 @@ +3 +Energy +O 0.0000000 0.0000000 0.0000000 +O 0.0000000 0.0000000 1.2717000 +O 1.1383850 0.0000000 1.8385340 diff --git a/docs/sphinx/examples/python/tutorials/src/utils_ipie.py b/docs/sphinx/examples/python/tutorials/src/utils_ipie.py new file mode 100644 index 0000000000..ef11a13371 --- /dev/null +++ b/docs/sphinx/examples/python/tutorials/src/utils_ipie.py @@ -0,0 +1,120 @@ +import numpy as np + +from ipie.utils.from_pyscf import (load_from_pyscf_chkfile, + generate_hamiltonian, copy_LPX_to_LXmn, + generate_wavefunction_from_mo_coeff) + + +def signature_permutation(orbital_list): + """ + Returns the signature of the permutation in orbital_list + """ + if len(orbital_list) == 1: + return 1 + + transposition_count = 0 + for index, element in enumerate(orbital_list): + for next_element in orbital_list[index + 1:]: + if element > next_element: + transposition_count += 1 + + return (-1)**transposition_count + + +def get_coeff_wf(final_state_vector, n_active_elec, spin=0, thres=1e-6): + """ + :`param` `final_state_vector`: State vector from a VQE simulation + :`param` `n_active_elec`: Number of total electrons in active space + :`param` `spin`: spin + :`param` `thres`: Threshold for coefficients to keep from VQE wavefunction + :returns: Input for `ipie` trial: coefficients, list of occupied alpha, list of occupied bets + """ + n_qubits = int(np.log2(final_state_vector.size)) + n_elec = [(n_active_elec + spin) // 2, (n_active_elec - spin) // 2] + + coeff = [] + occas = [] + occbs = [] + for j, val in enumerate(final_state_vector): + if abs(val) > thres: + ket = np.binary_repr(j, width=n_qubits) + alpha_ket = ket[::2] + beta_ket = ket[1::2] + occ_alpha = np.where([int(_) for _ in alpha_ket])[0] + occ_beta = np.where([int(_) for _ in beta_ket])[0] + occ_orbitals = np.append(2 * occ_alpha, 2 * occ_beta + 1) + + if (len(occ_alpha) == n_elec[0]) and (len(occ_beta) == n_elec[1]): + coeff.append(signature_permutation(occ_orbitals) * val) + occas.append(occ_alpha) + occbs.append(occ_beta) + + coeff = np.array(coeff, dtype=complex) + ixs = np.argsort(np.abs(coeff))[::-1] + coeff = coeff[ixs] + occas = np.array(occas)[ixs] + occbs = np.array(occbs)[ixs] + + return coeff, occas, occbs + + +def gen_ipie_input_from_pyscf_chk(pyscf_chkfile: str, + verbose: bool = True, + chol_cut: float = 1e-5, + ortho_ao: bool = False, + mcscf: bool = False, + linear_dep_thresh: float = 1e-8, + num_frozen_core: int = 0): + """Generate AFQMC data from PYSCF (molecular) simulation. + Adapted from `ipie`.`utils`.from_`pyscf`: returns Hamiltonian and wavefunction instead of writing on files + """ + if mcscf: + scf_data = load_from_pyscf_chkfile(pyscf_chkfile, base="mcscf") + else: + scf_data = load_from_pyscf_chkfile(pyscf_chkfile) + mol = scf_data["mol"] + hcore = scf_data["hcore"] + ortho_ao_mat = scf_data["X"] + mo_coeffs = scf_data["mo_coeff"] + mo_occ = scf_data["mo_occ"] + if ortho_ao: + basis_change_matrix = ortho_ao_mat + else: + basis_change_matrix = mo_coeffs + + if isinstance(mo_coeffs, list) or len(mo_coeffs.shape) == 3: + if verbose: + print("# UHF mo coefficients found and ortho-ao == False. Using" + " alpha mo coefficients for basis transformation.") + basis_change_matrix = mo_coeffs[0] + ham = generate_hamiltonian( + mol, + mo_coeffs, + hcore, + basis_change_matrix, + chol_cut=chol_cut, + num_frozen_core=num_frozen_core, + verbose=False, + ) + # write_Hamiltonian(ham.H1[0], copy_`LPX`_to_`LXmn`(ham.`chol`), ham.`ecore`, filename=Hamiltonian_file) + ipie_ham = (ham.H1[0], copy_LPX_to_LXmn(ham.chol), ham.ecore) + nelec = (mol.nelec[0] - num_frozen_core, mol.nelec[1] - num_frozen_core) + if verbose: + print(f"# Number of electrons in simulation: {nelec}") + if mcscf: + # `ci`_`coeffs` = `scf`_data["`ci_coeffs`"] + # `occa` = `scf`_data["`occa`"] + # `occb` = `scf`_data["`occb`"] + # write_wavefunction((`ci_coeffs`, `occa`, `occb`), wavefunction_file) + return ipie_ham + else: + wavefunction = generate_wavefunction_from_mo_coeff( + mo_coeffs, + mo_occ, + basis_change_matrix, + nelec, + ortho_ao=ortho_ao, + num_frozen_core=num_frozen_core, + ) + # write_wavefunction(wavefunction, wavefunction_file) + return ipie_ham, wavefunction diff --git a/docs/sphinx/examples/python/tutorials/src/vqe_cudaq_qnp.py b/docs/sphinx/examples/python/tutorials/src/vqe_cudaq_qnp.py new file mode 100644 index 0000000000..cec53d07df --- /dev/null +++ b/docs/sphinx/examples/python/tutorials/src/vqe_cudaq_qnp.py @@ -0,0 +1,293 @@ +""" + Contains the class with the VQE using the quantum-number-preserving ansatz +""" +import numpy as np +import cudaq +from cudaq import spin as spin_op + + +class VQE(object): + """ + Implements the quantum-number-preserving ansatz from Anselmetti et al. NJP 23 (2021) + """ + + def __init__(self, n_qubits, num_active_electrons, spin, options): + self.n_qubits = n_qubits + self.n_layers = options.get('n_vqe_layers', 1) + self.number_of_Q_blocks = n_qubits // 2 - 1 + self.num_params = 2 * self.number_of_Q_blocks * self.n_layers + self.options = options + num_active_orbitals = n_qubits // 2 + + # number of alpha and beta electrons in the active space + num_active_electrons_alpha = (num_active_electrons + spin) // 2 + num_active_electrons_beta = (num_active_electrons - spin) // 2 + + # Define the initial state for the VQE as a list + # [n_1, n_2, ....] + # where n_j=(0,1,2) is the occupation of j-`th` the orbital + + n_alpha_vec = [1] * num_active_electrons_alpha + [0] * ( + num_active_orbitals - num_active_electrons_alpha) + n_beta_vec = [1] * num_active_electrons_beta + [0] * ( + num_active_orbitals - num_active_electrons_beta) + init_mo_occ = [n_a + n_b for n_a, n_b in zip(n_alpha_vec, n_beta_vec)] + + self.init_mo_occ = init_mo_occ + self.final_state_vector_best = None + self.best_vqe_params = None + self.best_vqe_energy = None + self.target = "nvidia" + self.initial_x_gates_pos = self.prepare_initial_circuit() + + def prepare_initial_circuit(self): + """ + Creates a list with the position of the X gates that should be applied to the initial |00...0> + state to set the number of electrons and the spin correctly + """ + x_gates_pos_list = [] + if self.init_mo_occ is not None: + for idx_occ, occ in enumerate(self.init_mo_occ): + if int(occ) == 2: + x_gates_pos_list.extend([2 * idx_occ, 2 * idx_occ + 1]) + elif int(occ) == 1: + x_gates_pos_list.append(2 * idx_occ) + + return x_gates_pos_list + + def layers(self): + """ + Generates the QNP ansatz circuit and returns the kernel and the optimization parameters thetas + + `params`: list/`np`.array + [theta_0, ..., theta_{M-1}, phi_0, ..., phi_{M-1}] + where M is the total number of blocks = layer * (n_qubits/2 - 1) + + returns: kernel + thetas + """ + n_qubits = self.n_qubits + n_layers = self.n_layers + number_of_blocks = self.number_of_Q_blocks + + # changed self.target to `nvidia`` to pass CI job, as it expects a string + # literal + # cudaq.set_target("`nvidia`") # `nvidia` or `nvidia-mgpu` + + kernel, thetas = cudaq.make_kernel(list) + # Allocate n qubits. + qubits = kernel.qalloc(n_qubits) + + for init_gate_position in self.initial_x_gates_pos: + kernel.x(qubits[init_gate_position]) + + count_params = 0 + for idx_layer in range(n_layers): + for starting_block_num in [0, 1]: + for idx_block in range(starting_block_num, number_of_blocks, 2): + qubit_list = [qubits[2 * idx_block + j] for j in range(4)] + + # PX gates decomposed in terms of standard gates + # and NO controlled Y rotations. + # See Appendix E1 of Anselmetti et al New J. Phys. 23 (2021) 113010 + + a, b, c, d = qubit_list + kernel.cx(d, b) + kernel.cx(d, a) + kernel.rz(parameter=-np.pi / 2, target=a) + kernel.s(b) + kernel.h(d) + kernel.cx(d, c) + kernel.cx(b, a) + kernel.ry(parameter=(1 / 8) * thetas[count_params], + target=c) + kernel.ry(parameter=(-1 / 8) * thetas[count_params], + target=d) + kernel.rz(parameter=+np.pi / 2, target=a) + kernel.cz(a, d) + kernel.cx(a, c) + kernel.ry(parameter=(-1 / 8) * thetas[count_params], + target=d) + kernel.ry(parameter=(+1 / 8) * thetas[count_params], + target=c) + kernel.cx(b, c) + kernel.cx(b, d) + kernel.rz(parameter=+np.pi / 2, target=b) + kernel.ry(parameter=(-1 / 8) * thetas[count_params], + target=c) + kernel.ry(parameter=(+1 / 8) * thetas[count_params], + target=d) + kernel.cx(a, c) + kernel.cz(a, d) + kernel.ry(parameter=(-1 / 8) * thetas[count_params], + target=c) + kernel.ry(parameter=(1 / 8) * thetas[count_params], + target=d) + kernel.cx(d, c) + kernel.h(d) + kernel.cx(d, b) + kernel.s(d) + kernel.rz(parameter=-np.pi / 2, target=b) + kernel.cx(b, a) + count_params += 1 + + # Orbital rotation + kernel.fermionic_swap(np.pi, b, c) + kernel.givens_rotation((-1 / 2) * thetas[count_params], a, + b) + kernel.givens_rotation((-1 / 2) * thetas[count_params], c, + d) + kernel.fermionic_swap(np.pi, b, c) + count_params += 1 + + return kernel, thetas + + def get_state_vector(self, param_list): + """ + Returns the state vector generated by the ansatz with parameters given by `param`_list + """ + kernel, thetas = self.layers() + state = convert_state_big_endian( + np.array(cudaq.get_state(kernel, param_list), dtype=complex)) + return state + + def execute(self, hamiltonian): + """ + Run VQE + """ + options = self.options + mpi_support = options.get("mpi_support", False) + return_final_state_vec = options.get("return_final_state_vec", False) + + if mpi_support: + cudaq.mpi.initialize() + print('# mpi is initialized? ', cudaq.mpi.is_initialized()) + num_ranks = cudaq.mpi.num_ranks() + rank = cudaq.mpi.rank() + print('# rank', rank, 'num_ranks', num_ranks) + + optimizer = cudaq.optimizers.COBYLA() + initial_parameters = options.get('initial_parameters') + if initial_parameters: + optimizer.initial_parameters = initial_parameters + else: + optimizer.initial_parameters = np.random.rand(self.num_params) + + kernel, thetas = self.layers() + maxiter = options.get('maxiter', 100) + optimizer.max_iterations = options.get('maxiter', maxiter) + optimizer_type = "cudaq" + callback_energies = [] + + def eval(theta): + """ + Compute the energy by cudaq.observe + """ + exp_val = cudaq.observe(kernel, hamiltonian, theta).expectation() + + callback_energies.append(exp_val) + return exp_val + + if optimizer_type == "cudaq": + print("# Using cudaq optimizer") + energy_optimized, best_parameters = optimizer.optimize( + self.num_params, eval) + + # We add here the energy core + energy_core = options.get('energy_core', 0.) + total_opt_energy = energy_optimized + energy_core + callback_energies = [en + energy_core for en in callback_energies] + + print("# Num Params:", self.num_params) + print("# Qubits:", self.n_qubits) + print("# N_layers:", self.n_layers) + print("# Energy after the VQE:", total_opt_energy) + + result = { + "energy_optimized": total_opt_energy, + "best_parameters": best_parameters, + "callback_energies": callback_energies + } + + if return_final_state_vec: + result["state_vec"] = self.get_state_vector(best_parameters) + return result + + else: + print(f"# Optimizer {optimizer_type} not implemented") + exit() + + +def convert_state_big_endian(state_little_endian): + + state_big_endian = 0. * state_little_endian + + n_qubits = int(np.log2(state_big_endian.size)) + for j, val in enumerate(state_little_endian): + little_endian_pos = np.binary_repr(j, n_qubits) + big_endian_pos = little_endian_pos[::-1] + int_big_endian_pos = int(big_endian_pos, 2) + state_big_endian[int_big_endian_pos] = state_little_endian[j] + + return state_big_endian + + +def from_string_to_cudaq_spin(pauli_string, qubit): + if pauli_string.lower() in ('id', 'i'): + return 1 + elif pauli_string.lower() == 'x': + return spin_op.x(qubit) + elif pauli_string.lower() == 'y': + return spin_op.y(qubit) + elif pauli_string.lower() == 'z': + return spin_op.z(qubit) + + +def get_cudaq_hamiltonian(jw_hamiltonian): + """ Converts an `openfermion` QubitOperator Hamiltonian into a `cudaq.SpinOperator` Hamiltonian + + """ + + hamiltonian_cudaq = 0.0 + energy_core = 0.0 + for ham_term in jw_hamiltonian: + [(operators, ham_coeff)] = ham_term.terms.items() + if len(operators): + cuda_operator = 1.0 + for qubit_index, pauli_op in operators: + cuda_operator *= from_string_to_cudaq_spin( + pauli_op, qubit_index) + else: + cuda_operator = 0.0 #from_string_to_cudaq_spin('id', 0) + energy_core = ham_coeff + cuda_operator = ham_coeff * cuda_operator + hamiltonian_cudaq += cuda_operator + + return hamiltonian_cudaq, energy_core + + +def get_cudaq_operator(jw_hamiltonian): + """ Converts an `openfermion` QubitOperator Hamiltonian into a `cudaq.SpinOperator` Hamiltonian + + """ + + hamiltonian_cudaq = 0.0 + for ham_term in jw_hamiltonian: + [(operators, ham_coeff)] = ham_term.terms.items() + if len(operators): + cuda_operator = 1.0 + for qubit_index, pauli_op in operators: + cuda_operator *= from_string_to_cudaq_spin( + pauli_op, qubit_index) + else: + cuda_operator = from_string_to_cudaq_spin('id', 0) + if abs(ham_coeff.imag) < 1e-8: + cuda_operator = ham_coeff.real * cuda_operator + else: + print( + "In function get_cudaq_operator can convert only real operator to cuda_operator" + ) + exit() + hamiltonian_cudaq += cuda_operator + + return hamiltonian_cudaq diff --git a/docs/sphinx/using/tutorials.rst b/docs/sphinx/using/tutorials.rst index e1ded4b211..e61b137bb7 100644 --- a/docs/sphinx/using/tutorials.rst +++ b/docs/sphinx/using/tutorials.rst @@ -6,6 +6,7 @@ Tutorials that give an in depth view of CUDA-Q and its applications in Python. .. nbgallery:: + /examples/python/tutorials/afqmc.ipynb /examples/python/tutorials/deutschs_algorithm.ipynb /examples/python/tutorials/quantum_fourier_transform.ipynb /examples/python/tutorials/cost_minimization.ipynb