From 87ef3190a8d628675672a040a95a4a1403230382 Mon Sep 17 00:00:00 2001 From: Guen Prawiroatmodjo Date: Mon, 14 Sep 2020 14:49:45 -0700 Subject: [PATCH] Add session 2 Chemistry demo --- .devcontainer/devcontainer.json | 5 +- Dockerfile | 5 + session2/Chemistry/data/hydrogen_0.2.yaml | 57 ++ session2/Chemistry/vqe.ipynb | 629 ++++++++++++++++++++++ 4 files changed, 695 insertions(+), 1 deletion(-) create mode 100644 session2/Chemistry/data/hydrogen_0.2.yaml create mode 100644 session2/Chemistry/vqe.ipynb diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json index 45ffd48..257ebdf 100644 --- a/.devcontainer/devcontainer.json +++ b/.devcontainer/devcontainer.json @@ -4,6 +4,9 @@ "IQSHARP_HOSTING_ENV": "ieee-quantum-week-2020-devcontainer", }, "extensions": [ - "quantum.quantum-devkit-vscode" + "quantum.quantum-devkit-vscode", + "editorconfig.editorconfig", + "ms-python.python", + "robbowen.synthwave-vscode" ] } diff --git a/Dockerfile b/Dockerfile index 023b022..50a59f5 100644 --- a/Dockerfile +++ b/Dockerfile @@ -9,6 +9,11 @@ RUN chown -R ${USER} ${HOME} # Install RISE. RUN pip install rise +RUN pip install ipykernel +RUN pip install scipy +RUN pip install matplotlib +RUN pip install seaborn +RUN pip install pyscf # Finish by dropping back to the notebook user. USER ${USER} diff --git a/session2/Chemistry/data/hydrogen_0.2.yaml b/session2/Chemistry/data/hydrogen_0.2.yaml new file mode 100644 index 0000000..0e479ce --- /dev/null +++ b/session2/Chemistry/data/hydrogen_0.2.yaml @@ -0,0 +1,57 @@ + +"$schema": https://raw.githubusercontent.com/Microsoft/Quantum/master/Chemistry/Schema/broombridge-0.2.schema.json + +bibliography: +- {url: 'https://www.nwchem-sw.org'} +format: {version: '0.2'} +generator: {source: nwchem, version: '6.8'} +problem_description: +- basis_set: {name: sto-3g, type: gaussian} + coulomb_repulsion: {units: hartree, value: 0.713776188} + energy_offset: {units: hartree, value: 0.0} + fci_energy: {lower: 0.0, units: hartree, upper: 0.0, value: 0.0} + geometry: + atoms: + - coords: [0.0, 0.0, 0.0] + name: H + - coords: [0.0, 0.0, 1.624] + name: H + coordinate_system: cartesian + symmetry: c1 + units: angstrom + hamiltonian: + one_electron_integrals: + format: sparse + units: hartree + values: + - [1, 1, -1.252477495] + - [2, 2, -0.475934275] + two_electron_integrals: + format: sparse + index_convention: mulliken + units: hartree + values: + - [1,1,1,1,0.674493166] + - [1,2,2,1,0.181287518] + - [1,1,2,2,0.663472101] + - [2,2,2,2,0.697398010] + initial_state_suggestions: + - energy: {units: hartree, value: -1.13} + method: sparse_multi_configurational + label: '|G>' + superposition: + - [1.0, (1a)+, (1b)+, '|vacuum>'] + - label: "UCCSD |G>" + method: unitary_coupled_cluster + cluster_operator: # Initial state that cluster operator is applied to. + reference_state: [1.0, (1a)+, (1b)+, '|vacuum>'] # A one-body cluster term is t^{q}_{p} a^\dag_q a_p # A one-body unitary cluster term is t^{q}_{p}(a^\dag_q a_p- a^\dag_p a_q) + one_body_amplitudes: # t^{q}_{p} p q + - [0.001, "(2a)+", "(1a)"] + - [-0.001, "(2b)+", "(1b)"] + two_body_amplitudes: # t^{pq}_{rs} p q r s # If this is a PQQR term, the middle two indices must coincide. + - [-0.001, "(2a)+", "(2b)+", "(1a)", "(1b)"] + metadata: {molecule_name: H2} + n_electrons: 2 + n_orbitals: 2 + scf_energy: {units: hartree, value: -7.86099161} + scf_energy_offset: {units: hartree, value: 0.0} diff --git a/session2/Chemistry/vqe.ipynb b/session2/Chemistry/vqe.ipynb new file mode 100644 index 0000000..ab7ff47 --- /dev/null +++ b/session2/Chemistry/vqe.ipynb @@ -0,0 +1,629 @@ +{ + "metadata": { + "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.7.9-final" + }, + "orig_nbformat": 2, + "kernelspec": { + "name": "python_defaultSpec_1600118182129", + "display_name": "Python 3.7.8 64-bit" + } + }, + "nbformat": 4, + "nbformat_minor": 2, + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Quantum Chemistry: Variational Quantum Eigensolver\n", + "## Simulate the dissociation curve (energy vs. bond length) for a $H_2$ molecule\n", + "Based on http://openqemist.1qbit.com/docs/vqe_microsoft_qsharp.html#the-variational-quantum-eigensolver-vqe" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": "3.63 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n" + } + ], + "source": [ + "%%timeit -n 1 -r 1\n", + "import qsharp" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "output_type": "display_data", + "data": { + "text/plain": "Adding package microsoft.quantum.chemistry.jupyter.", + "application/json": "{\"LastUpdated\":\"2020-09-14T21:22:52.968931+00:00\",\"IsCompleted\":false,\"Description\":\"Adding package microsoft.quantum.chemistry.jupyter\",\"Subtask\":null}" + }, + "metadata": {} + }, + { + "output_type": "stream", + "name": "stdout", + "text": "1.89 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n" + } + ], + "source": [ + "%%timeit -n 1 -r 1\n", + "from qsharp.chemistry import load_broombridge, load_fermion_hamiltonian, load_input_state, encode" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Load computed $H_2$ from file\n", + "To generate:\n", + "\n", + "1. Go to: https://arrows.emsl.pnnl.gov/api/qsharp_chem\n", + "\n", + "2. Enter `HH theory{qsharp_chem}`\n", + "\n", + "3. Under Datafiles, click \"download\"" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import qsharp\n", + "from qsharp.chemistry import load_broombridge, load_fermion_hamiltonian, load_input_state, encode\n", + "\n", + "filename = 'data/hydrogen_0.2.yaml'\n", + "broombridge_data = load_broombridge(filename)\n", + "problem_description = broombridge_data.problem_description[0]\n", + "ferm_hamiltonian = problem_description.load_fermion_hamiltonian()\n", + "input_state = load_input_state(filename, \"UCCSD |G>\")\n", + "num_qubits, hamiltonian_term_list, input_state_terms, energy_offset = encode(ferm_hamiltonian, input_state)" + ] + }, + { + "source": [ + "### Hamiltonian terms\n", + "## $H = \\sum_{pq}{h_{pq}} a_p^{\\dagger}a_q^{\\dagger} + \\frac{1}{2} \\sum_{pqrs}{h_{pqrs}} a_p^{\\dagger}a_q^{\\dagger}a_{r}a_{s}$" + ], + "cell_type": "markdown", + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": "([([0], [0.17120128499999998]),\n ([1], [0.17120128499999998]),\n ([2], [-0.222796536]),\n ([3], [-0.222796536])],\n [([0, 1], [0.1686232915]),\n ([0, 2], [0.12054614575]),\n ([0, 3], [0.16586802525]),\n ([1, 2], [0.16586802525]),\n ([1, 3], [0.12054614575]),\n ([2, 3], [0.1743495025])],\n [],\n [([0, 1, 2, 3], [0.0, -0.0453218795, 0.0, 0.0453218795])])\n" + } + ], + "source": [ + "import pprint as pp\n", + "pp.pprint(hamiltonian_term_list)" + ] + }, + { + "source": [ + "### Input trial state terms\n", + "## $T(\\vec{\\theta}) = \\sum_{ij} \\theta_{ij} a_{i}^{\\dagger} a_{j} + \\theta_{ijkl} a_{i}^{\\dagger} a_{j}^{\\dagger} a_{k} a_{l}$" + ], + "cell_type": "markdown", + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": "(3,\n [((0.001, 0.0), [2, 0]),\n ((-0.001, 0.0), [3, 1]),\n ((0.001, 0.0), [2, 3, 1, 0]),\n ((1.0, 0.0), [0, 1])])\n" + } + ], + "source": [ + "pp.pprint(input_state_terms)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": "[0.001, -0.001, 0.001]" + }, + "metadata": {}, + "execution_count": 6 + } + ], + "source": [ + "def get_var_params(input_state_terms):\n", + " \"\"\"Get variational parameters from the input state\"\"\"\n", + " _, var_params = input_state_terms\n", + " return [param for ((param, _), _) in var_params[:-1]]\n", + "\n", + "def new_input_state(var_params, input_state_terms):\n", + " \"\"\"Returns a new input state with updated variational parameters\"\"\"\n", + " val, terms = input_state_terms\n", + " new_amps = [((param, val), indices) for (((_, val), indices), param) in zip(terms, list(var_params) + [1.0])]\n", + " return (val, new_amps)\n", + "\n", + "var_params = get_var_params(input_state_terms)\n", + "var_params" + ] + }, + { + "source": [ + "## 2. Import Q# library operation\n", + "To be able to access the namespace, we need to add it to the qsharp packages list and reload.\n", + "\n", + "We can then import the wrapper for the `EstimateEnergy` operation directly into our Python namespace using `QSharpCallable`." + ], + "cell_type": "markdown", + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "output_type": "display_data", + "data": { + "text/plain": "Adding package Microsoft.Quantum.Chemistry.", + "application/json": "{\"LastUpdated\":\"2020-09-14T21:22:55.4999385+00:00\",\"IsCompleted\":false,\"Description\":\"Adding package Microsoft.Quantum.Chemistry\",\"Subtask\":null}" + }, + "metadata": {} + }, + { + "output_type": "display_data", + "data": { + "text/plain": "Reloading workspace.", + "application/json": "{\"LastUpdated\":\"2020-09-14T21:22:56.2890218+00:00\",\"IsCompleted\":false,\"Description\":\"Reloading workspace\",\"Subtask\":null}" + }, + "metadata": {} + } + ], + "source": [ + "qsharp.packages.add(\"Microsoft.Quantum.Chemistry\")\n", + "qsharp.reload()\n", + "estimate_energy = qsharp.QSharpCallable(\"Microsoft.Quantum.Chemistry.JordanWigner.VQE.EstimateEnergy\", \"\")" + ] + }, + { + "source": [ + "### Resource estimation\n", + "\n", + "Estimate quantum resources required to run this algorithm." + ], + "cell_type": "markdown", + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": "{'CNOT': 896,\n 'QubitClifford': 1148,\n 'R': 168,\n 'Measure': 70,\n 'T': 0,\n 'Depth': 0,\n 'Width': 4,\n 'BorrowedWidth': 0}" + }, + "metadata": {}, + "execution_count": 8 + } + ], + "source": [ + "estimate_energy.estimate_resources(jwHamiltonian=(num_qubits, hamiltonian_term_list, input_state_terms, energy_offset), nSamples=1)" + ] + }, + { + "source": [ + "## 3. Quantum evaluation\n", + "Here we create a function that evaluates the Hamiltonian on a given number of qubits using the VQE algorithm and returns the energy result." + ], + "cell_type": "markdown", + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "def evaluate_H(var_params, num_qubits, hamiltonian_term_list, input_state_terms, energy_offset, n_samples, verbose=False):\n", + " \"\"\"Runs EstimateEnergy with the input state modified by new variational parameters\"\"\"\n", + " input_state_terms = new_input_state(var_params, input_state_terms)\n", + " energy = estimate_energy.simulate(jwHamiltonian=(num_qubits, hamiltonian_term_list, input_state_terms, energy_offset), nSamples=1e18)\n", + " print(f\"Parameters: {var_params}, calculated energy: {energy}\") if verbose else _\n", + " return energy" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": "Parameters: [0.001, -0.001, 0.001], calculated energy: -1.1163206832928165\n" + }, + { + "output_type": "execute_result", + "data": { + "text/plain": "-1.1163206832928165" + }, + "metadata": {}, + "execution_count": 10 + } + ], + "source": [ + "evaluate_H(var_params, num_qubits, hamiltonian_term_list, input_state_terms, energy_offset, n_samples=1e18, verbose=True)" + ] + }, + { + "source": [ + "### Minimize the energy\n", + "The goal of VQE is to estimate the ground state energy of a given Hamiltonian variationally. Here, we run the VQE circuit multiple times to get the expectation value of the energy with varying input parameters, until the global minimum is found." + ], + "cell_type": "markdown", + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.optimize import minimize\n", + "\n", + "def VQE(initial_var_params, num_qubits, hamiltonian_term_list, input_state_terms, energy_offset, n_samples, verbose):\n", + " opt_result = minimize(evaluate_H,\n", + " initial_var_params,\n", + " args=(num_qubits, hamiltonian_term_list, input_state_terms, energy_offset, n_samples, verbose),\n", + " method=\"COBYLA\",\n", + " tol=0.000001,\n", + " options={'disp': True, 'maxiter': 200,'rhobeg' : 0.05})\n", + " return opt_result" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": "Parameters: [ 0.001 -0.001 0.001], calculated energy: -1.1163206834000161\nParameters: [ 0.051 -0.001 0.001], calculated energy: -1.1143693480147085\nParameters: [0.001 0.049 0.001], calculated energy: -1.1144494156823095\nParameters: [ 0.001 -0.001 0.051], calculated energy: -1.0941301538175583\nParameters: [-0.0033645 -0.00518542 -0.04863299], calculated energy: -1.1305266765693185\nParameters: [-0.01021579 -0.01175558 -0.09772365], calculated energy: -1.1366498431269147\nParameters: [-0.04496727 -0.04625581 -0.10782658], calculated energy: -1.133169622130144\nParameters: [ 0.00741028 -0.02948449 -0.09781082], calculated energy: -1.1362581504555616\nParameters: [ 0.00679366 0.01160861 -0.13852552], calculated energy: -1.1360580830871494\nParameters: [-1.71106807e-03 -7.34886915e-05 -1.18124585e-01], calculated energy: -1.1372266413690135\nParameters: [ 0.01101093 0.02025182 -0.11105121], calculated energy: -1.1367639233856788\nParameters: [-0.01406918 0.00127264 -0.13981479], calculated energy: -1.135961700306503\nParameters: [ 0.00830554 -0.00728485 -0.12010293], calculated energy: -1.1371149026083591\nParameters: [ 0.00150415 0.0047795 -0.11585007], calculated energy: -1.1372355787290442\nParameters: [ 0.0017774 0.00595035 -0.11873452], calculated energy: -1.1371844502573945\nParameters: [-0.00125662 0.00459899 -0.11024578], calculated energy: -1.137241855571082\nParameters: [ 0.00140629 0.00312955 -0.10952794], calculated energy: -1.1372393005840267\nParameters: [-0.00492482 0.00126103 -0.10644251], calculated energy: -1.137181594775063\nParameters: [-0.00032518 0.00719314 -0.10877323], calculated energy: -1.1372007749607476\nParameters: [-0.0010556 0.00376601 -0.11325101], calculated energy: -1.1372597707651262\nParameters: [-0.00254376 0.00111338 -0.1139683 ], calculated energy: -1.1372640475482676\nParameters: [-0.0024217 0.00103689 -0.11552415], calculated energy: -1.1372560571424941\nParameters: [-0.00321375 0.00150872 -0.1140403 ], calculated energy: -1.1372607055177557\nParameters: [-0.00212284 0.00020201 -0.11277096], calculated energy: -1.1372668576685356\nParameters: [-0.00105628 0.00115993 -0.11214948], calculated energy: -1.1372675764115958\nParameters: [ 0.00017566 0.0002499 -0.11245859], calculated energy: -1.1372697332775559\nParameters: [ 0.00108608 -0.00070472 -0.11162119], calculated energy: -1.1372660135422867\nParameters: [ 0.00045258 0.00050012 -0.11314493], calculated energy: -1.1372699641057173\nParameters: [ 0.00021951 0.0002478 -0.11333095], calculated energy: -1.1372701921393849\nParameters: [ 0.00059892 -0.00038644 -0.1135842 ], calculated energy: -1.1372696598328536\nParameters: [-0.00041218 0.00053469 -0.11369014], calculated energy: -1.1372695004400217\nParameters: [ 5.60496883e-05 9.09562614e-05 -1.13012727e-01], calculated energy: -1.1372703992925959\nParameters: [-0.00033171 0.00012435 -0.11297928], calculated energy: -1.1372703188995201\nParameters: [ 9.79181262e-05 -2.93785092e-04 -1.12959727e-01], calculated energy: -1.1372703320693396\nParameters: [ 0.00011357 0.00015267 -0.11283657], calculated energy: -1.1372702962526056\nParameters: [-4.13604342e-05 9.09699775e-05 -1.13019656e-01], calculated energy: -1.1372704050329556\nParameters: [-4.01815908e-05 1.36917258e-04 -1.13036137e-01], calculated energy: -1.1372703986045076\nParameters: [-4.58184023e-05 3.74773718e-05 -1.13101236e-01], calculated energy: -1.1372704096470567\nParameters: [-8.70963830e-05 -4.86184787e-05 -1.13080733e-01], calculated energy: -1.137270404691916\nParameters: [ 4.19913960e-07 2.18723292e-05 -1.13099597e-01], calculated energy: -1.1372704120926564\nParameters: [ 3.22473501e-05 2.86754004e-05 -1.13135997e-01], calculated energy: -1.1372704037877162\nParameters: [-6.39551299e-06 6.35119960e-07 -1.13109526e-01], calculated energy: -1.1372704105543345\nParameters: [ 5.79237035e-06 1.21039465e-05 -1.13052059e-01], calculated energy: -1.1372704136692455\nParameters: [ 4.44582445e-05 3.40478120e-05 -1.13031870e-01], calculated energy: -1.137270410354743\nParameters: [-1.44437113e-05 2.29595809e-05 -1.13043770e-01], calculated energy: -1.1372704123860577\nParameters: [ 4.17335425e-06 -1.22375090e-05 -1.13051100e-01], calculated energy: -1.1372704140423169\nParameters: [ 2.63195354e-05 -1.88264596e-05 -1.13043215e-01], calculated energy: -1.137270413145167\nParameters: [ 4.07756103e-06 -1.60379657e-05 -1.13075217e-01], calculated energy: -1.1372704141568484\nParameters: [-7.62804061e-06 -1.94286876e-05 -1.13075919e-01], calculated energy: -1.1372704132368008\nParameters: [ 5.70637047e-06 -2.18494168e-05 -1.13074307e-01], calculated energy: -1.1372704138186829\nParameters: [ 1.16691464e-05 -6.71527857e-06 -1.13077330e-01], calculated energy: -1.1372704136422263\nParameters: [ 7.89737958e-06 -2.07630441e-05 -1.13074637e-01], calculated energy: -1.137270413435715\nParameters: [ 3.76731817e-06 -1.47338096e-05 -1.13069262e-01], calculated energy: -1.137270414012915\nParameters: [ 1.28980049e-06 -1.51401832e-05 -1.13076074e-01], calculated energy: -1.1372704133957372\nParameters: [ 6.10522230e-06 -1.38094276e-05 -1.13075702e-01], calculated energy: -1.1372704130345634\nParameters: [ 3.86696384e-06 -1.59061826e-05 -1.13073711e-01], calculated energy: -1.1372704140874856\nParameters: [ 4.29479349e-06 -1.75477273e-05 -1.13075175e-01], calculated energy: -1.1372704140517813\nParameters: [ 3.10102628e-06 -1.60643727e-05 -1.13075430e-01], calculated energy: -1.1372704141618564\nParameters: [ 3.15080028e-06 -1.52648313e-05 -1.13076029e-01], calculated energy: -1.1372704137160072\n fun: -1.1372704137160072\n maxcv: 0.0\n message: 'Optimization terminated successfully.'\n nfev: 59\n status: 1\n success: True\n x: array([ 3.15080028e-06, -1.52648313e-05, -1.13076029e-01])\n" + } + ], + "source": [ + "opt_result = VQE(var_params, num_qubits, hamiltonian_term_list, input_state_terms, energy_offset, n_samples=1e18, verbose=True)\n", + "print(opt_result)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": "Difference with exact FCI value: 8.376432880652374e-09\n" + } + ], + "source": [ + "# Print difference with exact FCI value known for this bond length\n", + "fci_value = -1.1372704220924401\n", + "print(\"Difference with exact FCI value: \", abs(opt_result.fun - fci_value))" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "def compute_energy_vqe(problem_description, verbose=False):\n", + " \"\"\"Parse the problem description and variationally find the ground state energy in Ha using VQE\"\"\"\n", + " ferm_hamiltonian = problem_description.load_fermion_hamiltonian()\n", + " input_state = load_input_state(filename, \"UCCSD |G>\")\n", + " num_qubits, hamiltonian_term_list, input_state_terms, energy_offset = encode(ferm_hamiltonian, input_state)\n", + " var_params = get_var_params(input_state_terms)\n", + " opt_result = VQE(var_params, num_qubits, hamiltonian_term_list, input_state_terms, energy_offset, n_samples=1e18, verbose=verbose)\n", + "\n", + " if opt_result.status == 1:\n", + " return opt_result.fun" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": "-1.1372704137575473" + }, + "metadata": {}, + "execution_count": 15 + } + ], + "source": [ + "problem_description = broombridge_data.problem_description[0]\n", + "compute_energy_vqe(problem_description)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 3. Calculate one and two-electron integrals using PySCF for a given bond length\n", + "\n", + "First, create a function that returns a molecule with a varying bond length." + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from pyscf import gto, scf, fci, ao2mo\n", + "\n", + "def create_H2_molecule(bond_length):\n", + " # Create molecule object with PySCF\n", + " H2 = [['H',[ 0, 0, 0]], ['H',[0,0, bond_length]]]\n", + " mol = gto.Mole()\n", + " mol.atom = H2\n", + " mol.basis = \"sto-3g\"\n", + " mol.charge = 0\n", + " mol.spin = 0\n", + " mol.build()\n", + " return mol" + ] + }, + { + "source": [ + "Here are some functions that wrangle the output from PySCF into Broombridge format that is compatible with the Q# operation `EstimateEnergy`." + ], + "cell_type": "markdown", + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [], + "source": [ + "from pyscf import gto, scf, fci, ao2mo\n", + "\n", + "def compute_integrals(problem_description, molecule, RHF):\n", + " \"\"\"\n", + " Compute one and two-electron integrals and output in Broombridge format\n", + " \"\"\"\n", + " one_electron_integrals = problem_description.hamiltonian['OneElectronIntegrals']['Values']\n", + " two_electron_integrals = problem_description.hamiltonian['TwoElectronIntegrals']['Values']\n", + " \n", + " # Get molecular orbital coefficients\n", + " h_pq = np.around(RHF.mo_coeff.T @ RHF.get_hcore() @ RHF.mo_coeff, decimals=6)\n", + " one_electron_integrals = [([j, i], h_pq[i-1, j-1]) for (j, i), _ in one_electron_integrals]\n", + "\n", + " # Get interaction terms\n", + " h_pqrs = np.around(ao2mo.restore(1, ao2mo.kernel(molecule, RHF.mo_coeff), 2), 6)\n", + " two_electron_integrals = [([l, k, j, i], h_pqrs[i-1, j-1, k-1, l-1]) for (l, k, j, i), _ in two_electron_integrals]\n", + "\n", + " return one_electron_integrals, two_electron_integrals\n", + "\n", + "def updated_problem_description(filename, one_electron_integrals, two_electron_integrals, nuclear_repulsion):\n", + " \"\"\"Update problem description with new one and two-electrom integral and nuclear repulsion terms\"\"\"\n", + " broombridge_data = load_broombridge(filename)\n", + " problem_description = broombridge_data.problem_description[0]\n", + " problem_description.hamiltonian['OneElectronIntegrals']['Values'] = one_electron_integrals\n", + " problem_description.hamiltonian['TwoElectronIntegrals']['Values'] = two_electron_integrals\n", + " problem_description.coulomb_repulsion['Value'] = nuclear_repulsion\n", + "\n", + " return problem_description" + ] + }, + { + "source": [ + "Run the PySCF Restricted Hartree-Fock simulation" + ], + "cell_type": "markdown", + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": "converged SCF energy = -1.10067441313401\n" + }, + { + "output_type": "execute_result", + "data": { + "text/plain": "([([1, 1], -1.184281), ([2, 2], -0.538244)],\n [([1, 1, 1, 1], 0.652123),\n ([1, 2, 2, 1], 0.188121),\n ([1, 1, 2, 2], 0.643609),\n ([2, 2, 2, 2], 0.676401)],\n 0.6157648096008939)" + }, + "metadata": {}, + "execution_count": 18 + } + ], + "source": [ + "molecule = create_H2_molecule(bond_length=1.624/1.88973)\n", + "RHF = scf.RHF(molecule)\n", + "RHF.scf()\n", + "energy_nuc = RHF.energy_nuc()\n", + "one_el, two_el = compute_integrals(problem_description, molecule, RHF)\n", + "one_el, two_el, energy_nuc" + ] + }, + { + "source": [ + "Use this as input to the VQE algorithm and evaluate quantumly" + ], + "cell_type": "markdown", + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": "-1.1270309926775335" + }, + "metadata": {}, + "execution_count": 19 + } + ], + "source": [ + "new_pd = updated_problem_description(filename, one_el, two_el, energy_nuc)\n", + "compute_energy_vqe(new_pd)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 4. Run VQE for varying bond length\n", + "\n", + "Now, we will take everything we've done and run it in a for-loop to calculate the energy for a list of bond lengths. " + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": "converged SCF energy = -0.904361394163539\nconverged SCF energy = -1.04299627454009\nconverged SCF energy = -1.1011282422677\nconverged SCF energy = -1.11734903499028\nconverged SCF energy = -1.11085039747659\nconverged SCF energy = -1.09191404102006\nconverged SCF energy = -1.06610864931794\nconverged SCF energy = -1.03653887502918\nconverged SCF energy = -1.00510670656849\nconverged SCF energy = -0.973110615777578\nconverged SCF energy = -0.941480654707798\nconverged SCF energy = -0.910873554594386\nconverged SCF energy = -0.881732449946057\n" + } + ], + "source": [ + "energy = []\n", + "bond_lengths = np.arange(0.4,1.7,0.1)\n", + "\n", + "for bond_length in bond_lengths:\n", + " molecule = create_H2_molecule(bond_length=bond_length)\n", + " RHF = scf.RHF(molecule)\n", + " RHF.scf()\n", + " energy_nuc = RHF.energy_nuc()\n", + " one_el, two_el = compute_integrals(problem_description, molecule, RHF)\n", + " new_pd = updated_problem_description(filename, one_el, two_el, energy_nuc)\n", + " _energy = compute_energy_vqe(new_pd)\n", + " energy.append(_energy)" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": "Text(0, 0.5, 'Energy (Ha)')" + }, + "metadata": {}, + "execution_count": 21 + }, + { + "output_type": "display_data", + "data": { + "text/plain": "
", + "image/svg+xml": "\n\n\n\n \n \n \n \n 2020-09-14T21:26:14.971675\n image/svg+xml\n \n \n Matplotlib v3.3.1, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEOCAYAAABbxmo1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA2T0lEQVR4nO3dd3TUVf7/8eeU9ExIIb0wIUgLUhOKKAiEEiEQWMGCqOjicVl3Xcv6VVkVWfdr2fWw2NafrAVk5Ytx6YYioCAtVAGBBAkkhHSS4KSXmfn9ETMSSBkgUzLzfpzDOcx8PjPzvklmXnM/9/O5V2E0Go0IIYQQFqC0dQFCCCEcl4SMEEIIi5GQEUIIYTESMkIIISxGQkYIIYTFSMgIIYSwGAkZIYQQFqO2dQH2pqysEoPhxi4dCgjwpqSkooMrsj5HaQdIW+yRo7QDpC0ASqUCPz+vVrdLyFzFYDDecMg0Pd4ROEo7QNpijxylHSBtaY8cLhNCCGExEjJCCCEsRkJGCCGExUjICCGEsBgZ+O8A+04WsHpnJqW6Wvx93JgxOoYRsSG2LksIIWxOQuYm7TtZwLJN6dQ1GAAo0dWybFM6gASNEMLpyeGym7R6Z6YpYJrUNRhYvTPTRhUJIYT9kJC5SSW62uu6XwghnImEzE0K8HG7rvuFEMKZSMjcpBmjY3BVN/8xuqqVzBgdY6OKhBDCfsjA/01qGtxfsfUM1bUN+GncuPtOObtMCCFAejIdYkRsCPOT+wHw6OQ+EjBCCPELCZkO0i1EA0BWQbmNKxFCCPshIdNBvD1cCAnwlJARQogrSMh0oJgIX7LydbYuQwgh7IaETAe6JcKXSz/XUFFdb+tShBDCLkjIdKAekb4AZMshMyGEACRkOlRMhC8AWQVyyEwIIUBCpkN5e7gQ5OdBVr70ZIQQAiRkOpw2RCM9GSGE+IWETAfThvhQoqtFV1Vn61KEEMLmJGQ6mPaXizJl8F8IISRkOpzpyn+5XkYIISRkOpqHm5oQf7nyXwghQELGIrShGgkZIYRAQsYitMEayspr+blCVscUQjg3CRkL0Ib6ADIjsxBCSMhYQFSwNwrkDDMhhJCQsQB3V7VM+y+EEEjIWIw2xIfzcuW/EMLJSchYiDZUw88VdZSVy+C/EMJ5SchYSHRI4+C/jMsIIZyZhIyFRAZ7o1DItP9CCOcmIWMhbi4qwrp6yeC/EMKpSchYkDZEQ1a+DqPRaOtShBDCJtS2LgCgurqaF154gZMnT6JSqfif//kfxowZc81+er2eN954g3379lFfX8/YsWN57rnnUCgUpKWl8dhjj6HVagFwdXUlJSXFyi1pThviw54TBZSV1+Lv427TWoQQwhbsImQ+/vhjvL29+eabb8jKymL27Nls3boVLy+vZvt99dVXnDt3jjVr1qBQKHj88cdJTU1l8uTJAMTExLB69WpbNKFFTdP+ZxWUS8gIIZySXRwu27RpE/fccw8AWq2Wfv36sWvXrmv2S09PZ8SIEbi4uKBWqxk5ciQbNmywdrlmiwzyRqlQyOC/EMJp2UXI5OXlER4ebrodGhpKQUHBNfvFxsayY8cOqqqqqKqqYtu2beTm5pq2Z2VlMX36dGbOnMmaNWusUntbXF1UhAfK4L8QwnlZ5XDZ9OnTycvLa3Hb3r17zX6eGTNmkJOTw3333Ye3tzf9+/dn//79QGMA7dy5E41GQ05ODnPnziU4OJjbbrvtumoNCPC+rv2vFhioaXa7t9aftJMFdO3qjUKhuKnntqar29GZSVvsj6O0A6Qt7bFKyLTXqwgLCyM3Nxd/f38A8vPzGTZs2DX7KZVKnnrqKZ566ikAli5dSkxMDADe3r+GQ2RkJAkJCRw5cuS6Q6akpAKD4cbOBgsM1FBc3LzXEuLrjq6yjvTMYrp28bih57W2ltrRWUlb7I+jtAOkLQBKpaLNL+d2cbhs0qRJrFq1Cmg85HXixAnuuOOOa/arra2lvLzxh5CXl8fKlSuZO3cuAEVFRaZThS9fvsyePXvo3bu3lVrQOtO0//mO8YcohBDXwy7OLnv00Ud5/vnnGT9+PEqlkkWLFpl6JkuWLCEoKIj77ruP8vJy5syZg1LZmI3PPvsssbGxAGzdupWVK1eiVqvR6/UkJyeTkJBgszY1iQj0RqVUkF1YTlzvIFuXI4QQVqUwypWCzXT04TKAhZ8eQOPhwjP3DrrZ8qxCDgHYJ0dpi6O0A6Qt0EkOlzk6bYgPWQXlcuW/EMLpSMhYgTZUQ2VNA5d+rrF1KUIIYVUSMlZw5ZX/QgjhTCRkrCC8qzdqlYKsfLnyXwjhXCRkrMBFrSQi0Ft6MkIIpyMhYyXaUBn8F0I4HwkZK9GGaKiubaDocrWtSxFCCKuRkLES0+C/XPkvhHAiEjJWEtbVC7VKKdP+CyGcioSMlahVSqKCvcmWwX8hhBORkLGibiEasgrKMcjgvxDCSUjIWJE2RENNnZ7C0ipblyKEEFYhIWNF0SGN0/7LITMhhLOQkLGi0K6euKqVclGmEMJpSMhYkUqpJDLYW6aXEUI4DQkZK9OG+JBdeONr1gghRGciIWNl2hANtfV6CmTwXwjhBCRkrOzXaf/lkJkQwvFJyFhZaIAXri5KmV5GCOEUJGSsTKlU0C1YI2eYCSGcgoSMDXQL0XChqBy9wWDrUoQQwqLUti7AGUWH+LDt0EXyS6qICPS2dTlCCCe272QBq3dmUqqrxd/HjRmjYxgRG9Jhzy89GRvQhsq0/0II29t3soBlm9Ip0dViBEp0tSzblM6+kwUd9hoSMjYQ7O+Jm6tKppcRQtjU6p2Z1DU0P2xf12Bg9c7MDnsNCRkbUCqaBv/lNGYhhO2U6Gqv6/4bISFjI9oQDReKKmjQy+C/EMK6Sn6u4bNN6a1uD/Bx67DXkoF/G9GGaqg/aCDvUiVRwRpblyOEcAKXK2r5el82O3/IBaCv1o+fLv5M/RWHzFzVSmaMjumw15SQsZErp/2XkBFCWFJ5VR2b9l9gx5GL6A1GRt4aStJtWgK6uFv87LJ2Q6akpITvv/+ejIwMdDodPj4+9OrVi5EjRxIYGNhhhTibQD8PPNxUZBWUc8cAW1cjhHBEVTX1bD6QwzeHcqir0zM8NoRpt2sJ8vM07TMiNoQRsSEEBmooLu74k5FaDZnMzEyWLFlCWloasbGxdO/ena5du1JZWcn69et5/fXXGTZsGH/84x/p0aNHhxfm6GTwXwhhKdW1DWw7fJEtaReoqm0grncQ026PJryrl9VraTVknn/+eR599FH+8Y9/4Orqes32uro6tm/fzoIFC1i1apVFi3RU2lAfth3KoUFvQK2SczCEEDenrl7PjiO5pO7PpqK6noE9upJ8R7RND8m3GjIpKSltPtDV1ZXExEQSExM7vChnoQ3R0KA3kltcSbcQGZcRQtyY+gYDu47lsXFfFj9X1BEb7c/0O7rTPczH1qXJwL8tXTntv4SMEOJ6NegN7P2xgA17zlOiq6VnRBcenxpLryg/W5dmYlbIVFRU8O6773Lw4EHKysowGn9d1fG7776zVG0OL9DXA083NVkF5Yy2dTFCiE7DYDCSdqqQdbvPU3S5muhQHx5O7ENfrR8KhcLW5TVjVsgsXLiQwsJC5s+fz5///Gf+/ve/8/HHHzNx4kRL1+fQFAoF2lCZ9l8IYR6D0ciRjGLWfH+O/JIqIoO8+ePd/RkQE2B34dLErJDZs2cPqamp+Pn5oVKpSEhI4NZbb+Xxxx/n4YcftnCJjq1biIatB3KobzDgopbBfyHEr5quYSnR1aLxdMFFpaS0vJbQAE/mJ/djcK9AlHYaLk3M+lQzGAxoNI1jBp6enpSXlxMYGEh2dnaHFLFu3TqSkpLo27cvK1asaHPfL7/8kvHjx5OQkMCiRYswXLEmS1vb7FV0iA96g5GLxRW2LkUIYUeunCEZoLyqntLyWsYMDOOvjw4jrneQ3QcMmBkyvXv35uDBgwDExcWxcOFCFi5ciFar7ZAi+vTpw+LFi5kyZUqb++Xk5PDee++xatUqtm7dSnZ2NuvXr293mz37dfBfDpkJIX6V8u3Za2ZIBjh+rgSl0v7DpYlZIfPaa68RHh4OwIIFC3B3d0en0/HWW291SBE9e/akR48eKJVtl7NlyxYSEhLw9/dHqVQyc+ZMUlNT291mzwK6uOPlriZbLsoUQtB4OvLGvVlcrqhrcXtHzpBsDWaNyURGRpr+HxAQwN/+9jeLFdSW/Px8wsLCTLfDwsLIz89vd9v1CAi4uZUqAwOv/1TknlF+XCyuuqHHWoo91XKzpC32x1HaAR3blmNnivnX6uPkFlfg5qKktv7ankygn4fFfn6WeN42Q2bfvn3tPsGIESPa3Wf69Onk5eW1uG3v3r2oVKp2n8NaSkoqMBiM7e/Yghud+ycswJPjZy+Rm3cZVxfb/ywsNYeRLUhb7I+jtAM6ri1l5bWs2vETB04XEeTrwZ9mDqCypp5lm9KbHTJzVStJvj3aIj+/G22LUqlo88t5myGzYMGCZrcLCgoICfl1dk6FQsH27dvbLWLNmjXt7mOO0NDQZmGVl5dHaGhou9vsnfaXwf+c4gpiwrrYuhwhhJXoDQa2H85l7ffnaNAbmXZ7NHcNj8JF/euXzaazywIsMEOyNbQZMjt27Gh2Oz4+/pr7rGnixInMnj2bJ554Al9fX1JSUkwnC7S1zd41Df5nF5RLyAjhJM5e/JnlWzK4WFxBv+7+zB7fk+ArZkeGX2dI7syua1oZS13ss3HjRt566y10Oh3bt2/no48+4pNPPqFHjx4sWbKEoKAg7rvvPiIjI5k/fz6zZs0CYOTIkUydOhWgzW32zt/HDY2nC1n5jnEIQQjRuvKqOlK+y2T38Xz8NG78fno/BvcMtNuLKW+WwnjlHDHtGDp0KAcOHLBkPTZnizEZgMVfHqOsvIZFjw67ocd3JDlmbp8cpS2O0g64vrYYjEZ2Hcvjv99lUlOnZ3x8JFNHanF3tY8pJG0yJiOsRxui4eT5Umrr9bjZweC/EKLjZBeU8/nWDM7l6egZ6cucCT0JD7y5M1k7izZD5v7772/WhausrGT27NnN9vnPf/5jmcqcjDZEg8FoJKeogh7hMi4jhCOoqqlnza7z7Dh6EY2HC7+d0ocRsSEOe2isJW2GzMyZM5vdvvvuuy1ajDPThjau+5CVr5OQEaKTMxqN7D9ZyKpvz1JeWceYweHMGNUdT3cXW5dmdW2GzPTp061Vh9Pz9Xali5erTC8jRCeXe6mSFVsyyMi5THSohj/N7I82xPaLh9lKq/O4mHP9y/XsJ9qmUCjoFqIhW0JGiE6ppq6BlG/PsvCTA1wsruDBib1YMCfOqQMG2ujJpKamsnjxYpKSkoiPjyc6OhovLy8qKyvJysri4MGDrF+/nt69ezNu3Dhr1uywtCEaTpwroaauwW7OOBFCXKtpCv5SXS3+Pm4M7hnI4TPFlOpquf3WUO4eE4OPp6uty7QLrX6Svf3222RkZLBq1Sqee+45Ll68aBqsioqKYtSoUSxevJhbbrnFasU6Om2oD0YjXCisoGekr63LEUK0oGkK/qbpXkp0tXxz6CJ+Gleenz1Y3rtXafPrcq9evXj55ZcBqK6uRqfT4ePjg4eHh1WKczZXTvsvf6hC2KfVOzNbnIJfoVDI+7YFZh+T8fDwkHCxMF9vN3y9XWXafyHsWGtT7Zd2sin4rUUO/NsZbYiPnGEmhB0yGo3sPNbybPIAAT5uVqym85CQsTPaUA3Hzl6iurYBDzf59QhhD36uqOXTTekczywhLMCT4p9rqL9qCv4Zo2NsWKH9kk8xO6MN0WAELhSW0yvKz9blCOH0DqUXsXxLBrX1eu4bdwvj4iJIO1XY7OyyzjgFv7WYFTLLli0jKSkJf39/S9fj9Lr9ck59VoGEjBC2VFXTwBfbzrD3xwK6BWv4bVJfwrt6Ab9Owe9Ik31ailkhs3//fv75z38ydOhQpk2bRkJCAq6ucg64JXTxcsXfx03GZYSwodPZZXzy9SlKy2uZcpuWqSO1qFWtXrsu2mBWyPzrX/+irKyM1NRUli1bxiuvvMKECRNITk4mPj7e0jU6HW2ID1n5coaZENZW36DnvzvPsfVgDkF+Hrz4wBBiZC7Bm2L2mIyfnx+zZ89m9uzZpKen89xzz7F69WpCQ0OZOXMmDz74IF5eXpas1Wl0C9Fw5EwxVTUNeLrLsJkQ1pBdUM7SjafIu1TJmEHhzBrTAzdXWXbjZl3XJ9i+fftYv34927dvp1+/fvz2t78lLCyM5cuXM2/ePL744gtL1elUopuWYy4sp083GZcRwpL0BgOp+y+wfvd5vD1deGrWAG7tHmDrshyGWSHz5ptv8vXXX6PRaJg2bRobNmwgODjYtH3AgAEMHTrUYkU6m26mK/91EjJCWFBhWRX/3niKzFwdcb2DeHBiL7w9nG86fksyK2Rqa2t577336N+/f4vbXVxc+Oqrrzq0MGem8XSlaxd3svJl8F8ISzAajez8IY9VO86iVCqYl9SX4X2DnWoxMWsxK2Sa5i9rS0yMXIjUkWTafyEs43JFLZ/9cmFln25+PDq5D/4+7rYuy2GZFTJXL8PcxNXVlZCQEMaPH8/YsWM7vDhnpg3RcDijmMqaeryccDU9ISyh2YWVCbcwbkgESum9WJRZJ34PHTqU3Nxc4uPjmTp1KvHx8eTl5dGvXz8CAgJ48cUXWbp0qaVrdSqm5ZilNyPETauqaWDphlN8sPZHArq488rD8YyPi5SAsQKzejJ79uzh448/bnZILCkpieeff56UlBQmTJjA008/zbx58yxWqLPpFvzLGWYF5cRqZaYFIW7U6axSPk49zeXyOqaO1DLlNrmw0prMCplz584RGRnZ7L7w8HDOnz8PQP/+/SkpKen46pyYt4cLgb7uclGmENepadXKEl0t7q4qaur0BPt58MKcwcSEyYWV1mZWnMfHx/PCCy+QnZ1NbW0t2dnZ/OUvf2HIkCEAZGRkEBgYaNFCnZFM+y/E9WlatbJpzZeaOj1KhYK7hneTgLERs0LmjTfewGAwMHnyZAYOHMjkyZMxGAy8/vrrQOMpzG+//bZFC3VG2lANl36uoaK63talCNEptLRqpcFoZP2e8zaqSLR7uEyv17Ns2TLeeOMN3n77bUpLS/H390ep/DWfunfvbtEinZU2+NeLMvtFyxXIQrSlqqah1VUrW7tfWF67PRmVSsUXX3yBi4sLSqWSrl27NgsYYTmmK//lokwh2nShsJxFyw62ul1WrbQds9IiOTmZlStXWroWcRVPdxeC/TxkXEaIVhiNRr77IZfXlh+mrl5P0m3dcFU3/1iTVStty6yzy44fP86KFSv4+OOPCQkJaXZh5n/+8x+LFScaezOZuT/bugwh7E5NXQPLt2Sw/2QhsdH+zJvSFx8vV0ICvExnlwXIqpU2Z1bIzJo1i1mzZlm6FtECbYgPB04Xoausw8dLFooTAiC3uIIP1v5IQWkVyXdEM2WEFqWy8ctv06qVwj6YFTLTp0+3dB2iFdGhTYP/5fSPkcF/IfacyOfzLRm4u6l59p6B9JGLle2aWSFjNBpJSUlh48aNlJWVsWHDBg4ePEhxcTF33XWXpWt0alHBGhQ0nmEmISOcWW29nv98c4bdx/PpHeXLY1Nj8fWWAX17Z9bA/5IlS/jqq6+45557yM/PByAkJIR///vfFi1OgIebmmB/T5mRWTi1/JJK/rb8ELuP5zPlNi3P3DtQAqaTMKsns2bNGtasWYO/vz8LFy4EICIigpycHEvWJn6hDdWQceGyrcsQwibSThXy2eZ0XFRKWbWyEzKrJ6PX6/Hy8gIwnVlWWVmJp6dnhxSxbt06kpKS6Nu3LytWrGhz3y+//JLx48eTkJDAokWLMBgar+5NS0tjwIABTJs2jWnTpjFz5swOqc0eaEN8KCuv5XKFXFAmnEd9g57Pt2Tw/9afJDLQm4Vz4yVgOiGzQmb06NG8/vrr1NXVAY1jNEuWLGHMmDEdUkSfPn1YvHgxU6ZMaXO/nJwc3nvvPVatWsXWrVvJzs5m/fr1pu0xMTGsW7eOdevWkZKS0iG12QNtyK+D/0I4g6KyKv72+WG+PZrLpGFRPHf/IFlYrJMyK2ReeOEFiouLGTJkCOXl5QwaNIi8vDyeffbZDimiZ8+e9OjRo92ZBLZs2UJCQoJpWpuZM2eSmpraITXYs6hgbxQg4zLCKRzOKOLVzw5y6XINf/jNrcwa00Om5u/EzBqT8fb25v3336ekpITc3FxCQ0NtMutyfn4+YWFhptthYWGmExEAsrKymD59Omq1mvvvv99hTr12d1UT2tVLpv0XDq1Bb+DLb8+y7dBFokM1/G5aP7r6eti6LHGTzAqZK/n5+VFTU2Ma9L96nZmWTJ8+nby8vBa37d27F5VKdb1lXCM2NpadO3ei0WjIyclh7ty5BAcHc9ttt13X8wQEeN9UHYGBmpt6fGt6a/05mlFksee/mrVexxqkLfbn6nYUlVbxjy+OcObCZZLu6M7cKbG4qDtH78VRfidgmbaYFTK7du1iwYIFFBcXN7tfoVBw+vTpdh+/Zs2aG6vuKqGhoc3CKi8vj9DQUKCxt9UkMjKShIQEjhw5ct0hU1JSgcFgvKH6AgM1FBdb5pBWdXUdZeW1JD2zzuJTZViyHdYmbbE/V7fjh7OX+HjjKfQGI/OT+xHXO4jLZZU2rNB8jvI7gRtvi1KpaPPLuVlfFRYtWsT8+fP54YcfSE9PN/0zJ2A60sSJE9m2bRulpaUYDAZSUlJITEwEoKioCKOxMRwuX77Mnj176N27t1Xrs5R9Jws4mF5kul2iq2XZpnT2nSywYVVC3JwGvYGUb8/yzlfHCfBx55W58cT1DrJ1WaKDmdWT0el03Hvvvc0mxuxIGzdu5K233kKn07F9+3Y++ugjPvnkE3r06MGSJUsICgrivvvuIzIykvnz55vmURs5ciRTp04FYOvWraxcuRK1Wo1eryc5OZmEhASL1Gttq3dm0qBv3ruqazCwememzNEkOo2mZZFLdbX4ervh6qKksKyaOweGcV/CLbiob/6wubA/CmPT1/82vPnmm8TExHD33XdboyabssfDZY+8saPVbZ88P7bDX08OAdinztyWpmWRr161csygMOZM7LxHHDrz7+RqljpcZlZP5tixY3z++ecsXbqUrl27NtsmU/1bXoCPW4sr+8lCTKKzaGlZZIDjmSU2qEZYk1khM3PmTIe6gr6zmTE65ppvgbIQk+hMZFlk5yVT/XcCTeMuTQsxAUwaFiXjMaJTuFhcgVKpaPEwtPTGHV+bZ5e99tprzW5fPVXLH/7wh46vSLRoRGwIf58/kg+eHoWXu5qcogpblyREuw6cLuS15YdwUytRq5qfOCS9cefQZsisXr262e2///3vzW7v2bOn4ysSbXJ3VTNuSARHf7pE3qXOcS2BcD56g4H/2/4TH647SVSwhr89Npy5d/UhwMcNBY09mIcSe0tv3Am0ebjs6hPPzDgRTVjBuCERbE67wOa0CzwyuY+tyxGiGV1lHR+u+5H0C5cZNziCe8Y1zj3WtCyyI52RJdrXZk/m6utiLHWdjLg+Gk9X7ugfxr6TBZTqamxdjhAmmXk/8+pnB8nM0/HbKX2YPaGnTG7p5Nrsyej1evbv32/qwTQ0NDS73bSWi7C+iUMj+fZoLlsP5nDvuFtsXY5wckajkZ3H8vjimzP4ervx4gND6BbiOHN6iRvXZsgEBATw4osvmm77+vo2u+3v72+5ykSbuvp6MLRvEDuP5THlNi3eHi62Lkk4qfoGPSu2nuH74/n0i/bnsamx8vcoTNoMmR07Wr/SXNhe4rBu7D9ZyLdHLpI0MtrW5QgnVPJzDe+vOUFWQTlTbtOSfHs0SqUcVhe/uu6p/oX9iAzypn9MANsOX2Ti0ChcXWTuJ2E9p7JK+XDdSfQGA3+YcSuDelp/jSlh/2RErpNLHBZFeVU9u0/kt7+zEB3AaDSyKS2bt1f9gI+XK395ME4CRrRKejKdXM9IX2LCfNicdoHRA8NQtbOEtRA3o7q2gU9TT3Moo5i4XoHMvasPHm7yMSJaJ59InZxCoeCu4d249HNNszVnhOho+SWVvLb8EIfPFDNrTA9+l9xPAka0S/5CHMCAW7oSGuDJpv0XGNYnWK5nEh3uyJli/r3xFGqVkmfvGUgfrZxZKswjPRkHoFQomDQsipyiCk6eL7V1OcKBGAxGVu/K5L3VJwgN8OSVh+MlYMR1kZBxECNiQ/DTuJG6P9vWpQgHUVFdzz9TjrFxbzajBoTy/OzBBHRxt3VZopORw2UOQq1SMiE+klU7znIuT0f3MB9blyQ6seyCct5fc4LLFbU8NKkXoweG27ok0UlJyDiQUQPC2LAni037s/n9jFttXY7oRPadLDCtV+Tloaa6poEu3m78z+zBxIR1sXV5ohOTw2UOxMNNzdgh4Rw5U0x+iSwDIMyz72QByzalmxbEq6xuwAhMHtFNAkbcNAkZB5MwJBK1WsnmtAu2LkV0Eqt3ZjZb2hvAaIRNMr4nOoCEjIPx8XLl9v6h7P2xgLJyWT9dtK+pB2Pu/UJcDwkZBzRxaBQGo5FvDubYuhRhxwxGIxv2ZrW6PcDHzXrFCIclIeOAgnw9GNonmO9+yKWqpt7W5Qg7VFlTzztfHWfNrnPEhPngqm7+UeCqVjJjdIyNqhOORELGQSUOi6KmTs+3R3NtXYqwM1kFOl799CAnz5cye3xPXpwzhIcSe5t6LgE+bjyU2JsRsSE2rlQ4AjmF2UFFBWvoF+3PNwdzGB8XKcsAiCtWr/wJHy8Xnn/g19OTR8SGSKgIi5CejAO7a3g3dFX17PmxwNalCBurrdfzydenWb45g15RvrzycLycniysQnoyDqxXlC/RoT5sSbvA6AFhsmKhkyosreL9NT+SW1zB1JFapo6U1SuF9UhPxoE1LgMQRdHlag5lyDIAzuhwRjGLlh2krLyGP80aQPId3SVghFVJT8bBDeoZSIi/J6n7s4nvHSTLADgJvcHAf787x+YDF4gO1fC75H507eJh67KEE5KejINrWgbgQmEFp7LKbF2OsILLFbX8feUPbD5wgTGDw3l+9hAJGGEzEjJOYERsCL7errIMgBPIuFDGwk8PklWgY15SX+ZM6IWLWt7mwnbkr88JuKiVjI+P5HR2GefzdbYuR1iA0WhkU1o2f1/5A55ual56ME5OSRZ2QULGSdw5MBwPN7VMeuiAqmrqeW/1CVK+zWRwr0BeeiiO8EBvW5clBCAD/07Dw03N2MHhpO7LprC0imB/T1uXJDrAhcJyPljzIyW6Gu4ddwvj4yLk5A5hV6Qn40QS4iJRqZRskmUAHMLu4/n87fPD1DXoee7+QUyIj5SAEXbHLkJm3bp1JCUl0bdvX1asWNHqfoWFhcyZM4chQ4YwY8aMa7Z/+eWXjB8/noSEBBYtWoTBYGjhWZxXF9MyAPlcrpBp3Dur+gY9n206zSepp+kR3oWFc4dyS4SvrcsSokV2cbisT58+LF68mI8++qjN/Tw9PXnyySepqKjgnXfeabYtJyeH9957j7Vr1+Lr68u8efNYv349ycnJFqy885k0NJKdP+TyzaEcZt7Zw9blCDM0LY1cqquli7crSqWCUl0tk0d0Y7pcXCnsnF30ZHr27EmPHj1QKtsuR6PREBcXh4fHtef8b9myhYSEBPz9/VEqlcycOZPU1FRLldxpBfl5EtcriO+O5lJV02DrckQ7rlwa2QhcrqijVFfLxKGR/GZ0jASMsHt20ZPpCPn5+YSFhZluh4WFkZ+ff93PExBwc2flBAZqburx1nB/Yh+eWryTgz9d4u6xt7S4T2doh7k6c1vW7t53zdLIAEd+usQT9wy2QUUdozP/Tq4mbWmbVUJm+vTp5OXltbht7969qFT2Mw19SUkFBoPxhh4bGKihuLi8gyvqeF3cVMRq/Vj73Vlu6xOIi7r5z7+ztMMcnb0txWXVrd7fWdvV2X8nV5K2gFKpaPPLuVVCZs2aNRZ/jdDQ0GZBlpeXR2hoqMVft7NKHN6Nf/zfD+z9sYDRA8NtXY64isFgZMvB1s8ClKWRRWdhF2MyHWHixIls27aN0tJSDAYDKSkpJCYm2rosu9Wnmx/aEA2b0i7ccM9NWEZhWRVvfHGElG8z6RaiuWZaGFkaWXQmdhEyGzduZNSoUWzevJklS5YwatQozp49C8CSJUtYuXIlAHq9nlGjRvHkk09y5swZRo0axbvvvgtAZGQk8+fPZ9asWUyYMIGIiAimTp1qszbZu8ZlALpRVFbNkTPFti5HAAajkR1HLvLKJwfILa5k3pS+vPxQHA//sjSyAlkaWXQ+CqPRKF9jr+AMYzJNDAYjC5bux8NNzUsPxZku5Ots7WhLZ2lLqa6GT1JPcyqrjH7R/jyc2Bt/H/dm+3SWtrTHUdoB0hawkzEZYZ+USgUTh0WxfHMGp7PL6Kv1t3VJTsdoNLL3xwK+2HYGgwEenNSL0QPC5Mp94TAkZJzcyH4hrPv+PJv2Z0vIWNnPlXUs35zO0Z8u0TOiC49M6UuQr6z7IhyLhIyTc1GrGB8fyVffZZJdUE63EMc559+eHUovYvmWDGrq9Nwztgfj4yNRSu9FOCC7GPgXttW4DICKTWmyDIClVVTX89H6k3yw9ke6dnFn4dx4Jg6NkoARDkt6MgJPdzV3Dgpnc9oFZoyqcqgrmO3J8cwSPtt0mvKqepLviOau4d1Qq+R7nnBs8hcuABgfF4lKqWDzgRxbl+Jwqmsb+GxTOv9MOYaXhwt/eTCOqSOjJWCEU5CejADA19uNmLAufHc0l++O5hLg48aM0TFyPcZNyrhQxsdfn6ZEV0Pi8CiSb+9+zcWVQjgyCRkBNM72ey5fZ7pdoqtl2aZ0AAmaG1BXr+e/O8/xzaEcgvw8eGH2EHpEdLF1WUJYnYSMAGD1zkzqr5rtt67BwOqdmRIy1+lcno5/bzxFQWkV4wZHcPedMbi52s8ksEJYk4SMABp7Ltdzv7hWg97A+j1ZpO7LxlfjyjP3DiRWrj0STk5CRgCNc2K1FChKReO4Qq8oPxtUZd+aVqws0dXSxcsVlVJBaXktt98ayr3jbsHTXd5eQsgIpABgxugYXK8akFarFHi5u/DWF0dJ+fbsNYfTnNmVK1ZC49X7peW1TIiP4JHJfSRghPiFhIwAGgf3H7pqtt+5d/Xhzd+N4I4BYWxKu8Bryw+RW1xh61Ltwn93Zra4YuXhDJnRWogrydctYTIiNoQRsSHXzMb6cGJvBvQI4LNN6bz62SHuvjOGhLgIp7xKvbq2ge9+yKVUxrCEMIuEjDDLoFsCiQnrwmeb0vm/7T9x7OwlHp3c55rp6B1VRXU92w7lsP3wRSprGnBRKanXX9uTkRUrhWhOQkaYzcfLlT/85lZ2Hcvj/7af5eWPD/DgpF4M7RNs69Ispqy8li0HLrDzhzxq6/UM7hnI5BHdKCitYtmm9GaHzGTFSiGuJSEjrotCoWD0wHB6d/Nj6YZTfLjuJD/8dIkHJvTE093F1uV1mKKyKjalXWDPiXwMBhjWN4i7hncjPLBxcaboUB8A09llMkOCEC2TkBE3JNjPkxceGMzX+7JZvzuLMxcv8+jkvvTp1rlPdb5YXEHqvmzSTheiUiq4vX8Yk4ZFtbjOS9MYlhCidRIy4oaplEqmjozm1u4BfLThFP9YeZQJQyOZMSqm083PlZn3M6n7sjn60yXcXFRMjI9iwtBIfL1ljEWImyEhI25adKgPCx+O58tvz7LlQA4nz5cyLymWyKDW1/22B0ajkfTsMjbuy+Z0dhle7mqm3R7NuCEReHs4zqE/IWxJQkZ0CDdXFXMm9mJAjwA+SU3nr8sOMmNUDBOG2t+KjwajkWNnL/H1vmzO5eno4uXKrDE9GD0wDA83eUsI0ZHkHSU6VP+Yrix6dCjLNqXz5bdnOZ55iUcn9yWgi+1PddYbDBw8XcTX+7PJLa6kaxd3HpzYi5G3huCilgkshbAECRnR4Xw8XXlixq3sPpHPF9t+4uVPDvDAhJ4M7xuMwkq9mqZ5xUp1tfj5uNFX60fGhcsUX64hrKsX86b0ZWjfIFTKzjV2JERnIyEjLEKhUHBH/zB6Rfnx742nWLrhFMfOXuKBCb04ca7Eoqf+Ns0r1nQNS6mult3HC+jq48YTM25l4C1d7e4QnhCOSkJGWFSQrwfP3z+Y1P3ZrNt9nh/PlVDXYKBBbwTaXxytQW+gqqaBiup6Kmvqqay+4v819VRUN1BZ3fT/xu0lupoWazECg3sGWqytQohrScgIi1MqFUy5TUu/7v68tvwwBoOx2fa6BgPLN2dw7OwlKqt/CY5fQqOmTt/q8yoU4OXugpeHC94eany93Qjv6s2+kwUt7i/ziglhfRIywmq0IT7XBEyT2no92YUVeLur6eLtSlhXL7w81Hj/EiLN/++Ct7sadzd1i4e9zuSUtRgoMq+YENYnISOsqrXF0QJ83Hj9seEd8hozRsfIvGJC2Ak5tUZYVUuLo3V0ALS0Ns5Dib1lChghbEB6MsKqmj7oLT2xZGtr4wghrEtCRlidTCwphPOQw2VCCCEsRkJGCCGExUjICCGEsBgJGSGEEBYjA/9XUSpvbk6rm328vXCUdoC0xR45SjtA2tLeYxRGo7HlS7CFEEKImySHy4QQQliMhIwQQgiLkZARQghhMRIyQgghLEZCRgghhMVIyAghhLAYCRkhhBAWIyEjhBDCYiRkhBBCWIyEzHU6f/4899xzDxMnTuSee+4hKyur1X3PnTvHgAEDePPNN61XoJnMbUdqaipJSUlMmTKFpKQkLl26ZN1CzWBOW0pKSnjsscdISkoiMTGRhQsX0tDQYP1i2/Dmm28yduxYevXqxZkzZ1rcR6/X8+qrr5KQkMD48eNJSUmxcpXmMact77//PpMnTyYpKYkZM2bw/fffW7lK85jTlib2/J43tx0d/p43iusyZ84c49q1a41Go9G4du1a45w5c1rcr6GhwfjAAw8Yn376aeMbb7xhzRLNYk47jh8/bkxMTDQWFRUZjUajUafTGWtqaqxapznMactrr71m+j3U1dUZ7777buPXX39t1Trbc/DgQWNeXp5xzJgxxoyMjBb3WbNmjfGRRx4x6vV6Y0lJifGOO+4w5uTkWLnS9pnTll27dhmrqqqMRqPRePr0aeOQIUOM1dXV1izTLOa0xWi0//e8Oe2wxHteejLXoaSkhFOnTjFlyhQApkyZwqlTpygtLb1m348++og777wTrVZr5SrbZ247PvvsMx555BECAwMB0Gg0uLm5Wb3etpjbFoVCQWVlJQaDgbq6Ourr6wkODrZFya2Ki4sjNDS0zX1SU1OZOXMmSqUSf39/EhIS2Lx5s5UqNJ85bbnjjjvw8PAAoFevXhiNRi5fvmyF6q6POW0B+37Pg3ntsMR7XkLmOuTn5xMcHIxKpQJApVIRFBREfn5+s/3S09PZvXs3Dz/8sA2qbJ+57cjMzCQnJ4fZs2czffp0PvjgA4x2Np+quW2ZP38+58+f5/bbbzf9GzJkiC1Kvin5+fmEhYWZboeGhlJQUGDDijrG2rVriYqKIiSkcy7Lbe/veXNZ4j0vIdPB6uvreemll3j11VdNH3ydlV6vJyMjg08//ZTPP/+cXbt2sW7dOluXdUM2b95Mr1692L17N7t27eLQoUN22QNwRgcOHGDJkiW8/fbbti7lhsh7vm2ynsx1CA0NpbCwEL1ej0qlQq/XU1RU1KwLWlxczIULF3jssccA0Ol0GI1GKioq+Otf/2qr0psxpx0AYWFhTJo0CVdXV1xdXRk3bhzHjx8nOTnZNoW3wNy2rFixgv/93/9FqVSi0WgYO3YsaWlpTJo0yUaV35jQ0FDy8vLo378/cG3PprM5evQof/7zn/nggw/o3r27rcu5IZ3hPW8uS7znpSdzHQICAujTpw8bN24EYOPGjfTp0wd/f3/TPmFhYaSlpbFjxw527NjBQw89xKxZs+zqj82cdkDj+Mbu3bsxGo3U19ezf/9+evfubYuSW2VuWyIiIti1axcAdXV17Nu3j1tuucXq9d6sSZMmkZKSgsFgoLS0lG3btjFx4kRbl3VDjh8/zlNPPcU777xDbGysrcu5YZ3hPW8uS7znJWSu08KFC1mxYgUTJ05kxYoVvPrqqwDMmzePEydO2Lg685nTjsmTJxMQEMBdd91FcnIyPXr04O6777Zl2S0ypy0vvvgihw8fJikpieTkZLRaLbNmzbJl2dd47bXXGDVqFAUFBcydO5fJkycDzdsxbdo0IiIimDBhArNmzeL3v/89kZGRtiy7Rea05dVXX6WmpoaXX36ZadOmMW3aNDIyMmxZdovMaUtnYE47LPGel5UxhRBCWIz0ZIQQQliMhIwQQgiLkZARQghhMRIyQgghLEZCRgghhMVIyAghbtgXX3zByJEjWbBgga1LEXZKQkYIccNWrlzJli1b0Ol0nD592tblCDskISOEBb377rs8++yzLW5LS0tj1KhRVq6oUVt1tebpp59m27Ztze7TarXU1dXh5ubWbHLLS5cukZiYSF1dXYfUKzovCRnhdMaOHUv//v0ZNGgQ8fHxPPbYY9fM2uxIOiLM0tPTSU9PZ9y4cc3uHzFiBCNGjKChoQE/Pz/T/V27dmXYsGGsWrXqpl5XdH4SMsIpffjhhxw9epTdu3cTEBDQKeeZsqZVq1aRlJSEQqFodv/mzZvx9fVl586d1NTUNNuWlJQkISMkZIRzc3NzY9KkSWRmZpruKy8v57nnnmP48OGMGTOGDz74AIPBAMDq1au57777ePPNN4mPj2fs2LHs3LnT9NicnBweeOABBg0axNy5cykrKzO7lsLCQv7whz8wfPhwxo4dy/Lly03b3n33XZ588kmee+45Bg0axOTJk5vNm3Xy5EmSk5MZNGgQf/zjH/nTn/7E4sWLqaqqYt68eRQVFTFo0CAGDRpEYWEh0DhFfWvPd7Vdu3YRHx9/Tb0HDx7kpZdeor6+nh07djTbPmDAAHJycsjNzTX7ZyAcj4SMcGrV1dWkpqYyYMAA031//etfKS8vZ9u2bXz++eesW7eO//73v6btx48fJzo6mv379/Pb3/6WBQsWmBZ2evbZZ4mNjSUtLY358+ezZs0as+owGAz87ne/o1evXuzatYtly5axbNmyZuve79ixg8mTJ3Po0CHGjh1r6n3V1dXxxBNPMH36dA4cOMCUKVNMYyeenp4sXbqUoKAgjh49ytGjR00rgrb2fFerqqri4sWL10zFv3HjRnx9fZk0aRKjR49m/fr1zbar1WqioqJIT08362cgHJOEjHBKv//974mLiyMuLo49e/bw6KOPAo2LNqWmpvLMM8/g7e1NREQEc+fObfYBGhYWxqxZs1CpVEyfPp3i4mIuXbpEXl4eJ06c4Mknn8TV1dXU0zHHiRMnKC0t5YknnsDV1ZXIyEhmzZpFamqqaZ8hQ4YwevRoVCoV06ZNM314Hzt2jIaGBh588EFcXFyYMGECt956a7uv2drzXa28vBwALy+vZvdv2LCBxMRE1Gq1aYr4q5dP9vLyMj1eOCdZtEw4pffff5/bbrsNvV7P9u3bmTNnDl9//TUKhYL6+vpmC4GFhYWZDjFB46B2k6Y16quqqigrK8PHxwdPT89mjzXnpILc3FyKioqIi4sz3afX65vdvvJ13d3dqa2tpaGhgaKiIoKDg5uNl5izJn1rz6dWN/9Y0Gg0AFRWVprWez979iynT5/mlVdeARpPpnBzc2Pz5s3ce++9psdWVlaaHi+ck4SMcGoqlYoJEybw8ssvc/jwYcaPH4+Liwt5eXn06NEDaFx9sukQU1sCAwPR6XRUVVWZgiYvL++awfKWhIaGEhERwdatW6+7DYGBgRQWFmI0Gk2vlZ+fb1pnxpzXb4unpydRUVGcP3/etBjchg0bAHjiiSdM+9XW1rJx40ZTyDQ0NHDhwgW7W+hOWJccLhNOzWg0sm3bNnQ6HTExMahUKiZNmsTixYupqKggNzeXTz/9lKlTp7b7XOHh4fTr1493332Xuro6Dh06xLfffmtWHf3798fLy4uPPvqImpoa9Ho9Z86c4fjx4+0+duDAgahUKlasWEFDQwPbtm1rNogfEBDA5cuXb+qw1ejRozl48CDQ+DPbsGEDTzzxBGvXrjX9W7p0KYcOHTL13I4fP054eDjh4eE3/Lqi85OQEU7p8ccfZ9CgQQwePJh//vOfvPHGG6blmF966SU8PDxISEjg/vvvZ8qUKfzmN78x63nffvttjh07xrBhw3j//ffNXhtdpVLx4Ycfmq5FGT58OH/5y1+oqKho97Gurq68++67fPXVV8THx7N+/XruvPNOXF1dAYiJiWHy5MkkJCQQFxfX7NCfuWbNmsWGDRswGo0cOXKE4uJiZs+eTWBgoOnfiBEj6Nevn2kp7A0bNjQ7dCack6yMKYQDmjlzJvfee6/Z4WiOZ555hsTERBISEtrdt6SkhAceeIC1a9eaxnGEc5KQEcIBHDhwgOjoaPz8/NiwYQOvvPIK27ZtIygoyNalCScnA/9COIDz58/zpz/9ierqaiIiInjnnXckYIRdkJ6MEEIIi5GBfyGEEBYjISOEEMJiJGSEEEJYjISMEEIIi5GQEUIIYTESMkIIISzm/wP79TvowyAolgAAAABJRU5ErkJggg==\n" + }, + "metadata": {} + } + ], + "source": [ + "import matplotlib as mpl\n", + "import seaborn as sns\n", + "import pylab as pl\n", + "%matplotlib inline\n", + "sns.set_theme(style=\"darkgrid\")\n", + "\n", + "pl.plot(bond_lengths, energy, \"o-\")\n", + "pl.xlabel(\"Bond length (Å)\")\n", + "pl.ylabel(\"Energy (Ha)\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ] +} \ No newline at end of file