From 508a57b17a11bc324baf5bb188e9da39d2f22bfa Mon Sep 17 00:00:00 2001 From: he-tianshen Date: Thu, 29 Sep 2022 00:02:00 +0800 Subject: [PATCH] test --- summer_ospp/2022/221cb0269/doc.ipynb | 776 ++++++++++++++++++ summer_ospp/2022/221cb0269/main.py | 207 +++++ ...\351\241\271\346\212\245\345\221\212.docx" | Bin 0 -> 17001 bytes 3 files changed, 983 insertions(+) create mode 100644 summer_ospp/2022/221cb0269/doc.ipynb create mode 100644 summer_ospp/2022/221cb0269/main.py create mode 100644 "summer_ospp/2022/221cb0269/\344\275\225\345\244\251\346\267\261-210610324-\347\273\223\351\241\271\346\212\245\345\221\212.docx" diff --git a/summer_ospp/2022/221cb0269/doc.ipynb b/summer_ospp/2022/221cb0269/doc.ipynb new file mode 100644 index 000000000..c3cc4462d --- /dev/null +++ b/summer_ospp/2022/221cb0269/doc.ipynb @@ -0,0 +1,776 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "a398f4dd-f015-4461-96b7-b607e02ce00e", + "metadata": {}, + "source": [ + "# 基于昇思MindQuantum，实现量子蒙特卡洛算法" + ] + }, + { + "cell_type": "markdown", + "id": "8124e84e-40b5-4f22-a910-cd85cad3ed99", + "metadata": {}, + "source": [ + "#### 本项目主要基于Quantum Computing Quantum Monte Carlo （arXiv:2206.10431v1)的思路实现。" + ] + }, + { + "cell_type": "markdown", + "id": "45e059b4-bb25-404d-8ce0-1257951d4956", + "metadata": {}, + "source": [ + "该文章基于Fermion Monte Carlo without fixed nodes: A game of life, death, and annihilation in Slater determinant space, J. Chem. Phys. 131, 054106 (2009), 将经典计算机上的量子蒙特卡洛方法迁移到了量子计算机上。" + ] + }, + { + "cell_type": "markdown", + "id": "6d805ff2-f59d-4e76-b532-7ca826339937", + "metadata": {}, + "source": [ + "安装所需要的包（ModelArts Notebook环境自带的包这里已经跳过，如有需要请自行先安装）" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fee30ab7-c28d-4e8a-88d8-dec73e778233", + "metadata": {}, + "outputs": [], + "source": [ + "!pip install --upgrade pip -i https://pypi.tuna.tsinghua.edu.cn/simple" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f22824e7-483b-4f82-9fc5-9cf42beb76f0", + "metadata": {}, + "outputs": [], + "source": [ + "!pip install https://hiq.huaweicloud.com/download/mindquantum/newest/linux/mindquantum-master-cp37-cp37m-linux_x86_64.whl -i https://pypi.tuna.tsinghua.edu.cn/simple" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "51b901e5-a1e0-4ada-9ab6-c64bba456956", + "metadata": {}, + "outputs": [], + "source": [ + "!pip install openfermionpyscf -i https://pypi.tuna.tsinghua.edu.cn/simple" + ] + }, + { + "cell_type": "markdown", + "id": "aeb7d88e-54fd-4f75-887f-563978d9913b", + "metadata": {}, + "source": [ + "导入所需要的包" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11e8aac4-de70-4cab-a34b-9602b547799d", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import math\n", + "import mindspore as ms\n", + "from mindspore.common.parameter import Parameter\n", + "import mindspore.context as context\n", + "\n", + "from openfermion.chem import MolecularData\n", + "from openfermionpyscf import run_pyscf\n", + "from mindquantum.core.gates import X\n", + "from mindquantum.core.circuit import Circuit, pauli_word_to_circuits, controlled, dagger\n", + "from mindquantum.core.operators import Hamiltonian\n", + "from mindquantum.simulator import Simulator\n", + "from mindquantum.framework import MQAnsatzOnlyLayer\n", + "from mindquantum.algorithm.nisq import generate_uccsd\n", + "\n", + "context.set_context(mode=context.PYNATIVE_MODE, device_target=\"CPU\")" + ] + }, + { + "cell_type": "markdown", + "id": "7c60bc18-f9cb-4c60-9a86-474f2d080b4b", + "metadata": {}, + "source": [ + "## 一、量子化学准备，VQE" + ] + }, + { + "cell_type": "markdown", + "id": "a7588b51-d2d2-4f9b-b098-87efe684937c", + "metadata": {}, + "source": [ + "定义所要求解的分子的几何构型，这里以直线型的H4分子为例。" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "2ba17b07-b4c5-466e-8cc0-cb07a562b4a9", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Geometry: \n", + " [['H', [0.0, 0.0, 2.0]], ['H', [0.0, 0.0, 4.0]], ['H', [0.0, 0.0, 6.0]], ['H', [0.0, 0.0, 8.0]]]\n" + ] + } + ], + "source": [ + "dist = 2.0\n", + "geometry = [\n", + " [\"H\", [0.0, 0.0, 1.0 * dist]],\n", + " [\"H\", [0.0, 0.0, 2.0 * dist]],\n", + " [\"H\", [0.0, 0.0, 3.0 * dist]],\n", + " [\"H\", [0.0, 0.0, 4.0 * dist]],\n", + "]\n", + "basis = \"sto3g\"\n", + "spin = 0\n", + "print(\"Geometry: \\n\", geometry)" + ] + }, + { + "cell_type": "markdown", + "id": "8842ce22-3b9b-4bc7-9dac-4b7c28909036", + "metadata": {}, + "source": [ + "利用openfermion(pyscf)自带的方法，计算出分子基态的Hartree-Fock, CCSD, FCI能量。其中FCI能量即为分子基态能量的精确值。" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "1e22c465-e294-44f7-aa36-ffa4979e6bea", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Hartree-Fock energy: -1.5756164767018666 Ha\n", + "CCSD energy: -1.9160861180298063 Ha\n", + "FCI energy: -1.8977806459898732 Ha\n" + ] + } + ], + "source": [ + "molecule_of = MolecularData(\n", + " geometry,\n", + " basis,\n", + " multiplicity=2 * spin + 1\n", + ")\n", + "molecule_of = run_pyscf(\n", + " molecule_of,\n", + " run_scf=1,\n", + " run_ccsd=1,\n", + " run_fci=1\n", + ")\n", + "\n", + "print(\"Hartree-Fock energy: %20.16f Ha\" % (molecule_of.hf_energy))\n", + "print(\"CCSD energy: %20.16f Ha\" % (molecule_of.ccsd_energy))\n", + "print(\"FCI energy: %20.16f Ha\" % (molecule_of.fci_energy))" + ] + }, + { + "cell_type": "markdown", + "id": "bc0aa49d-7bae-4f87-97d6-c9c9f8915dab", + "metadata": {}, + "source": [ + "保存计算结果，供后续调用" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "18b80b08-a5b2-4210-8d07-7b18e4699315", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "/home/ma-user/anaconda3/envs/MindSpore/lib/python3.7/site-packages/openfermion/testing/data/H4_sto3g_singlet\n" + ] + } + ], + "source": [ + "molecule_of.save()\n", + "molecule_file = molecule_of.filename\n", + "print(molecule_file)" + ] + }, + { + "cell_type": "markdown", + "id": "f673572e-e2f0-43ee-93a5-c150cdd3b36e", + "metadata": {}, + "source": [ + "构建VQE的初态，也即hartree-fock态，即在编号前molecule_of.n_electrons（也即总量子比特数的一半）个qubit上作用一个X门。" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "ab0a3510-b5de-44a9-a2cf-2ddb679c7e2d", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "q0: ──X──\n", + "\n", + "q1: ──X──\n", + "\n", + "q2: ──X──\n", + "\n", + "q3: ──X──\n" + ] + } + ], + "source": [ + "hartreefock_wfn_circuit = Circuit([X.on(i) for i in range(molecule_of.n_electrons)])\n", + "print(hartreefock_wfn_circuit)" + ] + }, + { + "cell_type": "markdown", + "id": "fad75005-65b2-43e8-aed7-074a82a71d94", + "metadata": {}, + "source": [ + "通过MindQuantum自带的方法构造基于uccsd的VQE线路ansatz" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "161f24bd-9def-44d6-99ed-116c9d2c2fe1", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ccsd:-1.9160861180298063.\n", + "fci:-1.8977806459898732.\n" + ] + } + ], + "source": [ + "ansatz_circuit, \\\n", + "init_amplitudes, \\\n", + "ansatz_parameter_names, \\\n", + "hamiltonian_QubitOp, \\\n", + "n_qubits, n_electrons = generate_uccsd(molecule_file)" + ] + }, + { + "cell_type": "markdown", + "id": "a96f45dd-28f0-4de5-a1c8-5ee05f59d37c", + "metadata": {}, + "source": [ + "由于制备hartree-fock态的线路不含参数，可以将hartreefock_wfn_circuit与ansatz_circuit合并，方便后续使用MQAnsatzOnlyLayer这一非常便捷的功能。" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "22e943c0-c76d-43ed-8594-20608a3ead88", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "==========================Circuit Summary==========================\n", + "|Total number of gates : 3332. |\n", + "|Parameter gates : 160. |\n", + "|with 14 parameters are : |\n", + "|p0, p4, p1, p5, p2, p6, p3, p7, p8, p9.. .|\n", + "|Number qubit of circuit: 8 |\n", + "===================================================================\n", + "Number of parameters: 14\n" + ] + } + ], + "source": [ + "total_circuit = hartreefock_wfn_circuit + ansatz_circuit\n", + "total_circuit.summary()\n", + "print(\"Number of parameters: %d\" % (len(ansatz_parameter_names)))" + ] + }, + { + "cell_type": "markdown", + "id": "2d69d9a5-b2af-41f9-be21-3f45646780b4", + "metadata": {}, + "source": [ + "定义模拟器和梯度更新算子。" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "648a8fdf-4743-42ec-9d2b-2571ecb56761", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Initial energy: -1.5756164789199829\n" + ] + } + ], + "source": [ + "sim = Simulator('projectq', total_circuit.n_qubits)\n", + "molecule_pqc = sim.get_expectation_with_grad(Hamiltonian(hamiltonian_QubitOp), total_circuit)\n", + "molecule_pqcnet = MQAnsatzOnlyLayer(molecule_pqc, 'Zeros')\n", + "initial_energy = molecule_pqcnet()\n", + "print(\"Initial energy: %20.16f\" % (initial_energy.asnumpy()))" + ] + }, + { + "cell_type": "markdown", + "id": "88b67676-0ac6-44f8-a7fc-5fe85eeb775b", + "metadata": {}, + "source": [ + "使用MindSpore的Adagrad算法进行梯度优化，标准的VQE流程。" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "b92765e2-9618-4f67-a866-384bb02721bb", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Step 0 energy -1.5756164789199829\n", + "Step 5 energy -1.7936779260635376\n", + "Step 10 energy -1.8394664525985718\n", + "Step 15 energy -1.8595349788665771\n", + "Step 20 energy -1.8704031705856323\n", + "Step 25 energy -1.8769905567169189\n", + "Step 30 energy -1.8812921047210693\n", + "Step 35 energy -1.8842602968215942\n", + "Step 40 energy -1.8863967657089233\n", + "Step 45 energy -1.8879827260971069\n", + "Step 50 energy -1.8891848325729370\n", + "Step 55 energy -1.8901070356369019\n", + "Step 60 energy -1.8908178806304932\n", + "Step 65 energy -1.8913655281066895\n", + "Step 70 energy -1.8917859792709351\n", + "Step 75 energy -1.8921068906784058\n", + "Step 80 energy -1.8923504352569580\n", + "Step 85 energy -1.8925340175628662\n", + "Step 90 energy -1.8926718235015869\n", + "Step 95 energy -1.8927748203277588\n", + "Step 100 energy -1.8928515911102295\n", + "Step 105 energy -1.8929085731506348\n", + "Step 110 energy -1.8929510116577148\n", + "Step 115 energy -1.8929824829101562\n", + "Step 120 energy -1.8930058479309082\n", + "Step 125 energy -1.8930232524871826\n", + "Step 130 energy -1.8930362462997437\n", + "Step 135 energy -1.8930459022521973\n", + "Step 140 energy -1.8930531740188599\n", + "Step 145 energy -1.8930586576461792\n", + "Step 150 energy -1.8930627107620239\n", + "Step 155 energy -1.8930658102035522\n", + "Step 160 energy -1.8930681943893433\n", + "Step 165 energy -1.8930699825286865\n", + "Step 170 energy -1.8930712938308716\n", + "Step 175 energy -1.8930723667144775\n", + "Step 180 energy -1.8930730819702148\n", + "Step 185 energy -1.8930737972259521\n", + "Optimization completed at step 186\n", + "Optimized energy: -1.8930737972259521\n", + "Optimized amplitudes: \n", + " [-7.2468901e-03 5.0181672e-02 -8.1651461e-17 2.1638997e-01\n", + " -2.9372861e-16 1.4097199e-01 -6.7032170e-03 4.8143927e-02\n", + " 5.7785130e-17 1.7130126e-17 -3.9021072e-01 -2.8981766e-01\n", + " 2.7006064e-17 -1.8819426e-17]\n" + ] + } + ], + "source": [ + "optimizer = ms.nn.Adagrad(molecule_pqcnet.trainable_params(), learning_rate=4e-2)\n", + "train_pqcnet = ms.nn.TrainOneStepCell(molecule_pqcnet, optimizer)\n", + "\n", + "eps = 1.e-8\n", + "energy_diff = eps * 1000\n", + "energy_last = initial_energy.asnumpy() + energy_diff\n", + "iter_idx = 0\n", + "while abs(energy_diff) > eps:\n", + " energy_i = train_pqcnet().asnumpy()\n", + " if iter_idx % 5 == 0:\n", + " print(\"Step %3d energy %20.16f\" % (iter_idx, float(energy_i)))\n", + " energy_diff = energy_last - energy_i\n", + " energy_last = energy_i\n", + " iter_idx += 1\n", + "\n", + "print(\"Optimization completed at step %3d\" % (iter_idx - 1))\n", + "print(\"Optimized energy: %20.16f\" % (energy_i))\n", + "print(\"Optimized amplitudes: \\n\", molecule_pqcnet.weight.asnumpy())" + ] + }, + { + "cell_type": "markdown", + "id": "99d5d711-fb3c-44b2-ab92-f7e42c2da22c", + "metadata": {}, + "source": [ + "可以看到上述的最优化能量-1.89307与精确值-1.89778仍有大于化学精度的误差。在VQE结果的基础上，我们进行量子蒙特卡洛算法Quantum Computing Quantum Monte Carlo (QQMC)的实现。" + ] + }, + { + "cell_type": "markdown", + "id": "ab666782-d631-46d5-9314-15d0167cf858", + "metadata": {}, + "source": [ + "## 二、量子蒙特卡洛，QQMC" + ] + }, + { + "cell_type": "markdown", + "id": "e44edfdb-ec39-4716-bb3d-5cb7a8fb09ff", + "metadata": {}, + "source": [ + "将上述计算得到的最优化参数导入VQE线路，作为U，需要注意这里的U不包含Hartree-Fock态的制备线路，只有ansatz" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "bfbe3da9-13fb-4b0d-8b47-0d6e937fc891", + "metadata": {}, + "outputs": [], + "source": [ + "U = ansatz_circuit.apply_value(dict(zip(ansatz_parameter_names, list(molecule_pqcnet.weight.asnumpy()))))\n", + "Udag = dagger(U)" + ] + }, + { + "cell_type": "markdown", + "id": "9e9d15b7-8943-49f4-b1a2-43973983bec9", + "metadata": {}, + "source": [ + "在论文中，作者使用一个量子线路来估计相似变换后的哈密顿量在qubit basis下的矩阵元。但在模拟器中，我们可以直接计算出这些矩阵元。将之前计算出的hamiltonian_QubitOp利用MindQuantum自带的split方法逐项拆分为一系列单个的Pauli串，将每一项先用pauli_word_to_circuits方法转化成Circuit类型，再用Circuit类型自带的matrix方法即可导入Pauli串对应的矩阵。将这些矩阵依照split时的权重重新加和，即可得到完整的哈密顿量。\n", + "\n", + "H4分子的哈密顿量还较小，我们可以做一个简单的验证，将其对角化与前文中的FCI能量做对比，发现是一致的。" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "e7d8fe9e-eba2-4ccb-b876-197b562b5dd7", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-1.8977806459900022\n" + ] + } + ], + "source": [ + "p = list(hamiltonian_QubitOp.split())\n", + "emptycirc = Circuit()\n", + "for i in range(n_qubits):\n", + " emptycirc += X.on(i)\n", + " emptycirc += X.on(i)\n", + "a = pauli_word_to_circuits(p[0][1])\n", + "circ = emptycirc + a\n", + "hamatrix = p[0][0].const * circ.matrix()\n", + "n_pau = len(p)\n", + "for i in range(1, n_pau):\n", + " a = pauli_word_to_circuits(p[i][1])\n", + " circ = emptycirc + a\n", + " hamatrix += p[i][0].const * circ.matrix()\n", + "hamatrix = np.dot(np.dot(Udag.matrix(), hamatrix),U.matrix())\n", + "hamatrix = hamatrix.real\n", + "\n", + "e, v = np.linalg.eig(hamatrix)\n", + "e.sort()\n", + "print(e[0])" + ] + }, + { + "cell_type": "markdown", + "id": "ac0d9405-7633-4bce-a3ff-78405f56783b", + "metadata": {}, + "source": [ + "对于较大的哈密顿量，这一计算可能耗时较长，可以启用以下的两个cell\n", + "中的四行代码存储矩阵文件，在后续使用时直接调用文件。" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "25592b52-969e-45c1-b7be-44e68d13b2e5", + "metadata": {}, + "outputs": [], + "source": [ + "# with open('H.npy', 'wb') as f:\n", + " # np.save(f, hamatrix)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "6f677822-7919-4db4-85fb-3fe2085fe195", + "metadata": {}, + "outputs": [], + "source": [ + "# with open('h4.npy', 'rb') as f:\n", + " # hamatrix = np.load(f)" + ] + }, + { + "cell_type": "markdown", + "id": "1639651a-dd66-44bf-8f25-b154aead1013", + "metadata": {}, + "source": [ + "可以将制备Hartree-Fock态的线路也看做一个酉变换。我们因此可以知道，如果定义一个二进制串00...0，将最低位的，数量为qubit数目一半的0翻转为1，得到的数会对应于Hartree-Fock态在哈密顿量的qubit basis表象下的量子态。在H4分子中，即对应第00001111个，也即第15个量子态。\n", + "\n", + "我们也可以做一个简单地验证，输出hamatrix[15,15]这个矩阵元，发现它与前文计算得到的Hartree-Fock能量是一致的。" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "id": "a8101859-c625-48b1-b41c-75b71e26f551", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-1.8930739379526567\n" + ] + } + ], + "source": [ + "hf = 2 ** molecule_of.n_electrons - 1\n", + "\n", + "print(hamatrix[hf,hf])" + ] + }, + { + "cell_type": "markdown", + "id": "d9f3879c-dbb2-4783-a125-04a17a793757", + "metadata": {}, + "source": [ + "接下来是QMC的主体，参考论文中的Algorithm 1 QC-FCIQMC，但需要注意的是，论文有多处错误：\n", + "\n", + "1）“Label the new walker $|\\phi_j \\rangle ...$”这一步，需要在两符号相乘的基础上再乘上一个负号，原因是式(6)的等式右端本身就带有负号\n", + "\n", + "2）没有给出S的更新方式。更新方式实际我参考了https://zhuanlan.zhihu.com/p/369346554 中的在每$A$个循环后，$S(\\tau)=S(\\tau-A \\delta \\tau)-\\frac{\\zeta}{A \\delta \\tau} \\ln \\frac{N_{\\mathrm{w}}(\\tau)}{N_{\\mathrm{w}}(\\tau-A \\delta \\tau)}$，其中损耗因子$\\zeta$我设置为了1\n", + "\n", + "3）没有给出更新S的条件。更新条件我参考了论文中的控制popolation至10000，设置为了每10步判断总population是否超过10000，在超过是更新S\n", + "\n", + "4）同种反号walker相互湮灭这一步骤，将for $i$ in $\\mathcal{D}$ do错误地写为了in $S$\n", + "\n", + "5）这一步骤循环的位置也写错了，应该写在与上一个for $i$ in $\\mathcal{D}$同级的地方，仅在n的循环的内层。" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "id": "01222bb7-3a05-4cd2-a961-5192279916d9", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "step 50: mixed energy = -1.8930737972, total pupulation = 1\n", + "step 100: mixed energy = -1.8930737972, total pupulation = 1\n", + "step 150: mixed energy = -1.8930737972, total pupulation = 3\n", + "step 200: mixed energy = -1.8930737972, total pupulation = 14\n", + "step 250: mixed energy = -1.8970708654, total pupulation = 72\n", + "step 300: mixed energy = -1.8967936314, total pupulation = 191\n", + "step 350: mixed energy = -1.8966750774, total pupulation = 521\n", + "step 400: mixed energy = -1.8966766404, total pupulation = 1286\n", + "step 450: mixed energy = -1.8965401927, total pupulation = 3283\n", + "step 500: mixed energy = -1.8966937984, total pupulation = 8485\n", + "step 550: mixed energy = -1.8969660751, total pupulation = 10304\n", + "step 600: mixed energy = -1.8972964796, total pupulation = 10272\n", + "step 650: mixed energy = -1.8976165738, total pupulation = 10254\n", + "step 700: mixed energy = -1.8979787387, total pupulation = 10279\n", + "VQE: -1.8930737972\n", + "QQMC: -1.8979787387\n", + "FCI: -1.8977806460\n" + ] + } + ], + "source": [ + "T = 700\n", + "dt = 0.01\n", + "zeta = 1.0\n", + "A = 10\n", + "maxpop = 10000\n", + "\n", + "D = set([hf])\n", + "Nmax = 2 ** n_qubits\n", + "pospop = [0 for i in range(Nmax)]\n", + "pospop[hf] = 1\n", + "negpop = [0 for i in range(Nmax)]\n", + "Slast = 0\n", + "S = 0\n", + "Nlast = 1\n", + "N = 1\n", + "for n in range(T):\n", + " toadd = []\n", + " for i in D:\n", + " if max(pospop[i], negpop[i]) == 0:\n", + " continue\n", + " isign = 1\n", + " if negpop[i] > 0:\n", + " isign = -1\n", + " for w in range(Nmax):\n", + " if w == i:\n", + " hii = hamatrix[i,i]\n", + " p_i = dt * (hii - S)\n", + " ran = np.random.uniform(0, 1, max(pospop[i], negpop[i]))\n", + " n_new = np.sum(ran < np.abs(p_i))\n", + " if n_new > 0:\n", + " if p_i < 0:\n", + " if isign > 0:\n", + " pospop[i] += n_new\n", + " else:\n", + " negpop[i] += n_new\n", + " N += n_new\n", + " else:\n", + " if isign > 0:\n", + " pospop[i] -= n_new\n", + " else:\n", + " negpop[i] -= n_new\n", + " N -= n_new\n", + " continue\n", + " Hjisign = 1\n", + " Hji = hamatrix[w,i]\n", + " if Hji < 0:\n", + " Hjisign = -1\n", + " if np.abs(Hji) > 1e-8:\n", + " ran = np.random.uniform(0, 1, max(pospop[i], negpop[i]))\n", + " prob = dt * np.abs(Hji)\n", + " n_new = np.sum(ran < prob)\n", + " if n_new > 0:\n", + " toadd.append(w)\n", + " if isign * Hjisign > 0:\n", + " negpop[w] += n_new\n", + " else:\n", + " pospop[w] += n_new\n", + " N += n_new\n", + " for i in D:\n", + " if pospop[i] > negpop[i]:\n", + " pospop[i] -= negpop[i]\n", + " N -= 2 * negpop[i]\n", + " negpop[i] = 0\n", + " else:\n", + " negpop[i] -= pospop[i]\n", + " N -= 2 * pospop[i]\n", + " pospop[i] = 0\n", + " \n", + " if n % A == A - 1:\n", + " if N > maxpop:\n", + " S = Slast - zeta / A / dt * np.log(N / Nlast)\n", + " Slast = S\n", + " Nlast = N\n", + " \n", + " for i in toadd:\n", + " D.add(i)\n", + "\n", + " if n % 50 == 49:\n", + " temp = sum(pospop) + sum(negpop)\n", + " mixed_energy = energy_i[0]\n", + " for i in range(Nmax):\n", + " if i != 15:\n", + " if pospop[i] > 0:\n", + " mixed_energy += hamatrix[i,hf] * pospop[i] / pospop[hf]\n", + " else:\n", + " mixed_energy -= hamatrix[i,hf] * negpop[i] / pospop[hf]\n", + " print('step %d: mixed energy = %.10f, total pupulation = %d' % (n + 1, mixed_energy, temp))\n", + " \n", + "mixed_energy = energy_i[0]\n", + "\n", + "for i in range(Nmax):\n", + " if i != 15:\n", + " if pospop[i] > 0:\n", + " mixed_energy += hamatrix[i,hf] * pospop[i] / pospop[hf]\n", + " else:\n", + " mixed_energy -= hamatrix[i,hf] * negpop[i] / pospop[hf]\n", + " \n", + "print('VQE: %.10f' % energy_i[0])\n", + "print('QQMC: %.10f' % mixed_energy)\n", + "print('FCI: %.10f' % molecule_of.fci_energy)" + ] + }, + { + "cell_type": "markdown", + "id": "19e1508e-70f6-4ccb-81ad-4664e8feb1e8", + "metadata": {}, + "source": [ + "可以发现，QQMC的结果比VQE的结果更接近FCI的精确值，验证了算法的有效性。\n", + "\n", + "此外，若将U改为空线路，则实现了经典计算机上的量子蒙特卡洛FCIQMC" + ] + }, + { + "cell_type": "markdown", + "id": "fb62ac99-4604-493f-8955-b2bbf37ae624", + "metadata": {}, + "source": [ + "## 三、可调参数" + ] + }, + { + "cell_type": "markdown", + "id": "f824d555-120b-45bf-9401-39069e76538a", + "metadata": {}, + "source": [ + "1. 原子/分子构型\n", + "2. 演化步数T\n", + "3. 演化步长dt\n", + "4. 损耗因子zeta\n", + "5. S的更新频率A\n", + "6. 最大总walker规模maxpop\n", + "\n", + "改变任何参数都只需要在变量的定义处作修改，不需要改变其他代码，会自动相应调整。" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "MindSpore", + "language": "python", + "name": "mindspore" + }, + "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.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/summer_ospp/2022/221cb0269/main.py b/summer_ospp/2022/221cb0269/main.py new file mode 100644 index 000000000..892f7662c --- /dev/null +++ b/summer_ospp/2022/221cb0269/main.py @@ -0,0 +1,207 @@ +import numpy as np +import math +import mindspore as ms +from mindspore.common.parameter import Parameter +import mindspore.context as context + +from openfermion.chem import MolecularData +from openfermionpyscf import run_pyscf +from mindquantum.core.gates import X +from mindquantum.core.circuit import Circuit, pauli_word_to_circuits, controlled, dagger +from mindquantum.core.operators import Hamiltonian +from mindquantum.simulator import Simulator +from mindquantum.framework import MQAnsatzOnlyLayer +from mindquantum.algorithm.nisq import generate_uccsd + +context.set_context(mode=context.PYNATIVE_MODE, device_target="CPU") + +dist = 2.0 +geometry = [ + ["H", [0.0, 0.0, 1.0 * dist]], + ["H", [0.0, 0.0, 2.0 * dist]], + ["H", [0.0, 0.0, 3.0 * dist]], + ["H", [0.0, 0.0, 4.0 * dist]], +] +basis = "sto3g" +spin = 0 + +molecule_of = MolecularData( + geometry, + basis, + multiplicity=2 * spin + 1 +) +molecule_of = run_pyscf( + molecule_of, + run_scf=1, + run_ccsd=1, + run_fci=1 +) +print("Hartree-Fock energy: %20.16f Ha" % (molecule_of.hf_energy)) +print("CCSD energy: %20.16f Ha" % (molecule_of.ccsd_energy)) +print("FCI energy: %20.16f Ha" % (molecule_of.fci_energy)) + +molecule_of.save() +molecule_file = molecule_of.filename + +hartreefock_wfn_circuit = Circuit([X.on(i) for i in range(molecule_of.n_electrons)]) +print(hartreefock_wfn_circuit) + +ansatz_circuit, \ +init_amplitudes, \ +ansatz_parameter_names, \ +hamiltonian_QubitOp, \ +n_qubits, n_electrons = generate_uccsd(molecule_file) + +total_circuit = hartreefock_wfn_circuit + ansatz_circuit +total_circuit.summary() + +sim = Simulator('projectq', total_circuit.n_qubits) +molecule_pqc = sim.get_expectation_with_grad(Hamiltonian(hamiltonian_QubitOp), total_circuit) +molecule_pqcnet = MQAnsatzOnlyLayer(molecule_pqc, 'Zeros') +initial_energy = molecule_pqcnet() + +optimizer = ms.nn.Adagrad(molecule_pqcnet.trainable_params(), learning_rate=4e-2) +train_pqcnet = ms.nn.TrainOneStepCell(molecule_pqcnet, optimizer) + +eps = 1.e-8 +energy_diff = eps * 1000 +energy_last = initial_energy.asnumpy() + energy_diff +iter_idx = 0 +while abs(energy_diff) > eps: + energy_i = train_pqcnet().asnumpy() + if iter_idx % 5 == 0: + print("Step %3d energy %20.16f" % (iter_idx, float(energy_i))) + energy_diff = energy_last - energy_i + energy_last = energy_i + iter_idx += 1 + +print("Optimization completed at step %3d" % (iter_idx - 1)) +print("Optimized energy: %20.16f" % (energy_i)) +print("Optimized amplitudes: \n", molecule_pqcnet.weight.asnumpy()) + +U = ansatz_circuit.apply_value(dict(zip(ansatz_parameter_names, list(molecule_pqcnet.weight.asnumpy())))) +Udag = dagger(U) + +p = list(hamiltonian_QubitOp.split()) +emptycirc = Circuit() +for i in range(n_qubits): + emptycirc += X.on(i) + emptycirc += X.on(i) +a = pauli_word_to_circuits(p[0][1]) +circ = emptycirc + a +hamatrix = p[0][0].const * circ.matrix() +n_pau = len(p) +for i in range(1, n_pau): + a = pauli_word_to_circuits(p[i][1]) + circ = emptycirc + a + hamatrix += p[i][0].const * circ.matrix() +hamatrix = np.dot(np.dot(Udag.matrix(), hamatrix),U.matrix()) +hamatrix = hamatrix.real + +# with open('H.npy', 'wb') as f: + # np.save(f, hamatrix) +# with open('h4.npy', 'rb') as f: + # hamatrix = np.load(f) + +hf = 2 ** molecule_of.n_electrons - 1 + +T = 700 +dt = 0.01 +zeta = 1.0 +A = 10 +maxpop = 10000 + +D = set([hf]) +Nmax = 2 ** n_qubits +pospop = [0 for i in range(Nmax)] +pospop[hf] = 1 +negpop = [0 for i in range(Nmax)] +Slast = 0 +S = 0 +Nlast = 1 +N = 1 +for n in range(T): + toadd = [] + for i in D: + if max(pospop[i], negpop[i]) == 0: + continue + isign = 1 + if negpop[i] > 0: + isign = -1 + for w in range(Nmax): + if w == i: + hii = hamatrix[i,i] + p_i = dt * (hii - S) + ran = np.random.uniform(0, 1, max(pospop[i], negpop[i])) + n_new = np.sum(ran < np.abs(p_i)) + if n_new > 0: + if p_i < 0: + if isign > 0: + pospop[i] += n_new + else: + negpop[i] += n_new + N += n_new + else: + if isign > 0: + pospop[i] -= n_new + else: + negpop[i] -= n_new + N -= n_new + continue + Hjisign = 1 + Hji = hamatrix[w,i] + if Hji < 0: + Hjisign = -1 + if np.abs(Hji) > 1e-8: + ran = np.random.uniform(0, 1, max(pospop[i], negpop[i])) + prob = dt * np.abs(Hji) + n_new = np.sum(ran < prob) + if n_new > 0: + toadd.append(w) + if isign * Hjisign > 0: + negpop[w] += n_new + else: + pospop[w] += n_new + N += n_new + for i in D: + if pospop[i] > negpop[i]: + pospop[i] -= negpop[i] + N -= 2 * negpop[i] + negpop[i] = 0 + else: + negpop[i] -= pospop[i] + N -= 2 * pospop[i] + pospop[i] = 0 + + if n % A == A - 1: + if N > maxpop: + S = Slast - zeta / A / dt * np.log(N / Nlast) + Slast = S + Nlast = N + + for i in toadd: + D.add(i) + + if n % 50 == 49: + temp = sum(pospop) + sum(negpop) + mixed_energy = energy_i[0] + for i in range(Nmax): + if i != 15: + if pospop[i] > 0: + mixed_energy += hamatrix[i,hf] * pospop[i] / pospop[hf] + else: + mixed_energy -= hamatrix[i,hf] * negpop[i] / pospop[hf] + print('step %d: mixed energy = %.10f, total pupulation = %d' % (n + 1, mixed_energy, temp)) + +mixed_energy = energy_i[0] + +for i in range(Nmax): + if i != 15: + if pospop[i] > 0: + mixed_energy += hamatrix[i,hf] * pospop[i] / pospop[hf] + else: + mixed_energy -= hamatrix[i,hf] * negpop[i] / pospop[hf] + +print('VQE: %.10f' % energy_i[0]) +print('QQMC: %.10f' % mixed_energy) +print('FCI: %.10f' % molecule_of.fci_energy) \ No newline at end of file diff --git "a/summer_ospp/2022/221cb0269/\344\275\225\345\244\251\346\267\261-210610324-\347\273\223\351\241\271\346\212\245\345\221\212.docx" "b/summer_ospp/2022/221cb0269/\344\275\225\345\244\251\346\267\261-210610324-\347\273\223\351\241\271\346\212\245\345\221\212.docx" new file mode 100644 index 0000000000000000000000000000000000000000..bf6a90ff461c61e680d24f25ab7bffeab4825438 GIT binary patch literal 17001 zcmeHu1DhsG(stXnZClfrwr!i!wr!i!p6;G@PusR_+cv*GyJz=c_xk?8yY*a^^;Blw zk&#uI6&V>3a+1IxC;(sp5C8xGgaAx&HmW9o000G`0077U5I~xOHr9^D){eSLZnnk_ z+H|f~mIQeqKomIuK%e#hcl;lI0=0=F7Cj6I!Vf|30aNvh63r9^y+e61#@J@hfZ;4M zLU#iBt?%tS9E=4dj4VU_1QV+s6A+96Q_Ja=kcf3Iq$lV+$$kks=Iqo<3%j&FMwnz< zi45!fLu|g5EUj&rQfz<>#C1LKV~S9+-}Ao*psGCrL{lZjsE}~=W5Nq~M_^vQ=v??o zULZ9l7A*(Bg44@|1&J@6qoYR%8Wd`&=@Oq~UdjmY3oy=3*^9 zrBXI&J>iqsLwkxsL`e@nt>Ge|w~Wy_g3Qvsb`^^RSpit(a4w|Qb$erh(lPplpK3TW zwAgRNuy#HZA#5IiW#fGbM>9h)Zk98em!Lf9G3M8UN50Adi>+_!?UCvNS^VA74S!oa z%1_`W&laN)zR!{M@c|4V_m}C37l+kw@mV7CIs2eLr>CyHv84k&-EZmtrse!zZx7 z;vPd&$2yVssD$qKilEJjC?{qONow-5?35QCQoWgT(irFFEyOfEz)N_9+k6vfo75za|XeryTpN{eSx?OBnm@J05ru)E=r3$o7=R4vm1Ds2~rF!#cC|#?%1+zmD&<8OTwi`mlW;F56A#nk;%N2-- zpAalm-+^%Mw7k%2)=!6-exr*tC10TxL8|L>H+S4oPX% zrAwpi7<;kI%ow%%RF;ohv0VWTJ&N&-$Vu-5hO_&4oK^Zf6?RMwxC;V?BTB`B^swK9 zw@L5o@I$a^(6lzBK*he!0i5lrWu2i_xG(tzot+o9B^G4wd7@<`RZ(sD#1zH*-_#-WcYAM_!C)lea zsB#A47$f>??aA3Ei&DneR$Nmo*;CG?hz?~*4r${R;BbK21CE1Hf~h#yX!)XSk}eoR z#;*3rRly;_#RIJgznJR-OxHHEFp-RL(*zJNQi&7lJ!GZ|w2^i`Q*^8q$?Io~5=NUV ztt^nP5YJ=VdE5zp>q;Wxv$ofEf(ze8i@xq@*TZa2fg_{kmZM`;TPj67q9zm4?W|l| zz$cMlNf9_d;yNEj=~7r!ec{_2_>NvArY1QFK@?{x4+WhzG0Oq@6{QeQOJ(7Q6P z@D#+NJw#)xV_59G+2(rSb!W1jKLT zeBmwRNd@^7YslB5)TR|luv|$9S8SU*1{zW4CqhvKK{5#7ggn4oe2c!o$Z=2DLMZJr z8vrjb>9;6`xyvr)VjPAC3ehf=%4+H7iJFI$z5=8Ube@g*dq|!mi*T` zSoh{F{SSOY87SHd-Opa3c*55L%(4lfHka2=r8H1Aj~yYZylrLKdwL1H(;G(h)IN{( zAt(c8dn#JJABn9~xk;`FE#}5(${3Oo^$sWqDtFT`r$$lTZsDsN6k1<&){@Zkx2H;2 zFa%M7D*E5zfDut4G~~fUpjIGE9UL4GSAVkNnzI-2Q@W~RTq4X}Wm_<47 zI!%0azPWi90K1)OCV+3hTcy_Hr0Xv#STIz0Ja7(u8Alj`$n89M6GwOx&RXZj%K*eT zL7v#(yIMnPma=JMuA{P|9iCF19zZvKmS~t;+7w4cdq}EqFX}O7_cT6xlu&Nt>S?>R zSJj9TV`dJ^r+*%1Q;lem&bGdCwBM-7_s0*F3Y(FPD?3^^pe;B)olqXFL4%6=h2^$8 zkUq*Ey-isW!=%gXw)mBSr1MvTIO_5&+*7U;Ui?i>N-pqi^j@U|f*f^sMN3_lS*c4B zK^<-7tb}>`Y<*2i+-`f(3zijI=?Rs-(_$EJ9+5^H!?l-J#kLaQk0ws1E*NIX=M%P! zqTG1Ba{pUfxwy~@0dG#5N6-s(UTWP2ZY!=>^F>ozy5-9cgbeZ9vdp=T1mI1j(?Jp|~{tREGh zSUz+HZ@+(Iuv01v`92m{G^W17b8G4fF`?aQ4M3dT+m^!N9M&?T=d7Ty#0v*Mkb_s9(jsDjY2P0#khosIbATfT|XzZ;4^7QC&JzcU_C>rJ}s+f%L zXbyUj3JIMcSdQw`7TAh_JHCN;orh+U8yWf{vk6=f-2iqVBH< z8i5ya!z8VMvV=XF#}K7oK4@ua4)Qs4nhfHCC@gc?H@u@?URqu^;(#q7!~^bTrTOw^ z7jAzdxi^s1Z{p__KL&rX?0AOYOY<&Ja#!=A==?Efio$R@h88d0taR_^oMglRS&xl# zCaCMNN&Z7$!B)7j{5aw=(ja5yd{oiWWo4Ce7>gfJT>-GeQ})W7c1&ea>B`*SgdB%u zf@KU@WzR3P+BpRMW1&LHP@Y5P-gX-4IsMI1cG|HGg;?&ph#cwF-evA}F=`f_=g)FR zkK&Xsl=I``C#Vm0sBS3k2rl2rPkfaSwRBm#`!zmqL~sm+Z?JcUhFRf7W2R|Si|GpnCg z_;JjnqLSx-Vd|7J*$imo zU<@#Yv4;7zY!aj=1lG=_;D@jh=A;mQs!+9Txq zSAcg=U~gmn^LZ~vMzFK~oVMrRsEuLNF8CH0?=pebz#2|^t5-+W^k>t7pJkq9b?_Wo zvaaDAZkuJFdkv-p3xTeMGHS4|hHi<#u>1G%Y@{VUyGH>Un)KN!nCi8`HnXdJ{!dje zZXZ8}5iIRWOB!i-ZQD;$88$JX@FYjB_oqx^^m0^JD5+s%a{C~~$Q<~xx3iXAkYc#K zPTq4JiHY8GJnNfq0-bW!sN5N-)7f`DaK5Ct+ zyn&_X3qByU+R9UCbefbsw)BZ{y{H_QYaVPm}PD+Ko)IM1zQ`9L3d9J%Gm6FUe);eWPe`7NzRF9?*6Oj1=L8yISbVH3j87Oe)?F3x z`_9NajE9MMj`d_Q>chNGMTs|INuywy|_8Rh%yk~1}gIgY;*$BARWK&9Ey78 zsfaVjTG0vx%p|}+q_Kk9>aNCUg&H1AGstsBm*$)@k%7|>|MDr2ft#=OhJwwIZ)d8e z#p52h5ePfO;J^X9;zT)N5AF>>AK0=x8tRAuwPxS4=L$svdH_Aons^yOy%#Yb@WWm| z^qeckK@?%+4bS()HG-e(mfsaDP}ctfAbBxMwI^n7PDytzMX=L}Jr;P7qb^oysAt8m zu@t<=FhFgBjLHYo4#h@D8@2NQ9`*!!G)Hdk16O!0cZ;+;A8zZecs0)9Fm=JK%@Q`= z%6$W_$ZUglxsyCO1XAV z8#ieQSBoOhX3?gCq#Xik8pGzX=~v(R?#sRAP9X0Mn_WCGoo~RExGwGNPV1k`>l;j7ZN2gE>{Gb#)XD&My`t#Nad!l9N@jT{UxR9ErC5@eN z2X{=GiqJ3m6jxUhHsN8VAt^#;@AzTZ0~!Ms$p#K#Jp=YEG*++@l|_V1lh-i|6IXI7 zS?fZ4Z{L^1oRLCDO{e+yV&VD`H0bgyYlD+lxOoA-jX;=wv)OZz{ptxk%AG{E(4~e$ zYtRwl4GEhiT}W#dvWWbk8xWQsuaf?gYASqM$JBTpqq3$1mRqylcuwpC9VV%7VNJ*8 z*D}u?D;`^H*EmVE336$|bno-kYU?%{++H2WSSMSMOC)y!-p}ZCX_A3BVBLG6l12?g znP{jEtBKY%TUq6ph`6EwZ2B$;sc1ALQA0fkvMXmb_XeuBD-uKx!nVc4!aNI&A0TW9 z+6eH4a&QG0VD2OBwUqE@-Julc<&s-GSDS9Ju)(iCzSh9dfX2yFvIXT~M<(OX#~dsu z2`J(hTnQp9jph`PkTUcaiuneZ4n&y2gI(#NFj1dN5A|gCW>T7)3Tfzej?5qgP`4;} zqSd23bZ=XOfXQRzOc^=uniJL2wxx>3gmDYI$VB@sdJ}O&I1>$o+}OMFKkuQo2O?tb z1IL=XIjUn}sWNxNV^3e~Fhu7kRYCWJjj>O^qG2onDMD{9d=Q9>c!Sp@Oy4q_yW!pd zXTv40Z`qT~^?ES;qB!T~w+0Ej#2y#bmjG~p8FV1R5*4M&huP#=R|B$YV+nCpNj~o< z=n;Rvf_+5XmPaX*&Q&6555=cF=k+V5^^Sc`>`}F#zbKdJVHV39w2VV7=iWo*x^}B= zC96=gn$Qm9s?DdK(o+VLXuE6~v2rnl0P|QD7n0WtINZbyQO*rYPN2%@q)O2)BdGFI zV%4BeFF7j^!w1biai%JkW*>n%tj^_UOwWqV(JBcIx%Wq?Y-R}=3ybh7djJ#7?<;YQ zA>$j9YrUp;EC@hrjWz$0OPkuamg9M5-~~L4E=GuwgR!;z@n82~e3tviCcyvz1_%HE zQ2(VoJDM3=8PoqGW&CXi=2TrW0$UWJ9sC+UxM{?F?C_Tw+q77<;l*-)i}U3#Rg1CP z5>pMw24~tx!o?yjAVBh61VL-gKmVt z%tS>v`r!u24ZxjN#fFz;nrK8AB6~sEdLT!oF=FBZ>EK5Uqh-F@NR8A<6v#&-iffVc zGyZJjOp8TyhD(%qn3j0*?d}kQVr9tl)Gh804(yLEQobWaBg$Fj6sg_m<}!4)_q~u@ z$Q2VyJ#0jMJ4mM2rl}12NfSFZ1berkU`@%8HD8MZ55M4_b>Ip-Mya=V`c_oWbxnXwLCP>sU*5v zS5LA#-X4xr$N4%wZl8x-oLdKXqRQ2`rgQGzYE&n(+g`UaI(*y={9HJ>Mtk7JB7+=& zvpE{nVX$jyeBqy3#Bqc|cS)PYDdhYUBSCp$$nf)b}F^gjqK9H~?Qx z;qbf#1pD0kW~m+6d=QVo!erb^6o~nTo9+-qKjWG%CJJR`uNMoq&CM+Eo~|DyGl#K% zJd@tim5tdlE3O@lU!^JBs=g_BR89E}1ezjA$;!}<22o`2d?yz5a!J#YBnwo+pf`eL zBQ>3nnPb-$!|F;xRFOCdsRFJ%vsW4s`-OVeomxOdt+>5v;KU#M-6Mk;j7ju@FknL+ z4CzI4WgkG^(2qtB4#CdZ>Ih8;=E0Qcn8W;6ld`n|-ts8yUQ{|c*-!>lPSo_khid^m z3I@H8r$twmbdw>tw^*)Xq9of$2GOdCJTe0>S+=+IjOfiSTJX1CARc=2amn-#FL6c3 ztO@9}8-e}S*k2KzfdXWu+PAxMgywHn$yb0B~u zHFgxODm!uoLyommqj($Tb0<1K@oOp|D3)9Z<1|M}e6l)b8VLK5a>1IzXpWN=J!XqJ z&|f9gPBtfXY+tzj2QLKcxQZ0aN-j%-^7ZZnn8rZd9h0=rB~8c02fXC2T6E-Y0u&1N z;3ag{N)t7_8+4AJz1z4hOOw1S?-3%2@*9tfzQJz|@M_a+ET3GJvp-yw%(gQa*rgqD zHI@b4WxQCFKaRBtb~3?Z)#7fv(T@1V+UhLM{{dahu_3st?Gu&3I~ig-{yM=!9+^Gv*iJJGHe)fCGN6SYRgKEC3RRu zwULtGBDY5Z_0`vdCP*{|5Z8n=y(%V(ZZC>^Bw%m==W z!T|9Oa#RD-yzg6XAECw{V<&tcu1=}qWhg&~@Z=acOT*oKij;V;;f#OvpvI%O_QQ0- zqN6HYQr7oCg|!sop|qC75pm`U;fDAf^t~kqPrO9w)g&9p#R=FBu}fr(KAH$Zt-d>; zZIgPJ4~3&k2c;GnsS(6alh@hQmxiCo)J>D&^=_UHY zbJl5r{)I#iP`=?D*_M#QNu!w-qOlGiTjggzis?7g=O7(%hF(ha>`E&g1mUtyOk4R@ zZYqjlvW8OnA4*@a^T(H2MLb%C`Z9lrouf5pVpiog03jy8Q!^}|ykMp4)DvUhK~hu5 z1Lt06iArmA86Qn6NStK{u@Kbyg;a9Vz)73|i?-=Y5M+~c8oCLdV5(WwF+#suMV|wMYYC!S z->7L{;OgX7<`R4I;Syh%8gk$SM(?lC+(&-M5>ukX1UExb8)tDF7)hqjWiI(fi4;6# zFjXIEDs~EfqD;<&Nt!Q=jzs%PVpaK}Xmt8g`g6u)TX^ON^KlrysD8$i=^zCj{g4SwYM|cgvYm+}Z4x6}Vx15o6U51x z?q}?ZvS&k>&qigAV^Yiy%!?#_ut+n-u;|Ijj<s#+QkbJWThXX^4Ma^RgGO+iBY-Uol3Su{rwuxQLa@B;Lq734E zJ&%)=b){MFE~}t*NBCJ^K2-duzo>x(=1^d7zqFfHV{kdm5ukVJvk09~r05RE3|e5V z@yg02?KBFO&8)qPAzl{?GWG)B!hgNK-;qN!gA25g(_9}#w()$!JC{GUU;@=H;rc5h z#uK|%O|PocU>_LOy(VXkPNuqZX#7@0hK*}TfE-%3dT=N&fyTL7ZN!;D1BCPx&OJ-h z;we~h7{0NUdut=q5@{E1&eb?u)^?{6wf2H?X=UP)hkM3W9OV_EL|i@BgM*Ep0d>kk zk@+tB&_arOK^9Zuhuq3uh-~Coxer=|_dF_ywV89-^$hV*H<(OVAf$YchlJf$RR_#| z$Ym&HMNXrh@h*S0PE28#m87@Py8bW>BFBTQYmE?63hllW@waCuAmer%(3vid%N_pI z{w_U4l6AdW36XUJs6npRhH}k^jLSQpPkrG(&2jtp$qyQz9uMvq001ce(kNVv4HW%OgjkOFrnS&8kA8 zJiR-tAh!DqDYcM^anZ1l&U)O#R(fXCwp8Kdl?M8_wqXJ zG>a%EzE0BB&GO<$l+Snsa;GmDXN@FaYVUuGR`m}1X5V^9hnLMJpi=2CK-urB=lhZa z^7K*p$}ol~$rwO0$R~1-JV*c(>qqc`aa}|LP9O*Y`ZUn^RFzmgd!ge0oKM1Mkce>n z(9ZwPxH|!%M*-WA<@58sof9X70cJkIwhBfM#Lg8VMjpm38lVrwl>hD1WDyAf1qo5d zAVMqGt0$t4=PaB**X>c2#r++^KX#XMrK2E(UX|ErSe%?ShL@?(ms*RXZ?tg6{g zAdK9*Xsr8i-hM{6JXWD$qr+?c1^>;E;nVr%zl=r06-W|nm`~?8M*x6N3-Lcq!w!ya zmd5`O?5&w=&MS?G0~a1L8-AjyXFWMwOi=HJx?}XfAgA%5LI08~g;s{3h#0Z$*4okH!CjqN1}{8rLbQu-tVYI7 z8~Mhu^NVPhB|Ab`XC~RpOfpS4LMoVhLb)iRahp`?2lePWApj}Ss1coMANCyP3Q>T< zzDLV`C5aAjw5_if00*3KpdpTCznVH52j_r(w1w3B_sb&r*ubI%_~Ow;78Md!7QbU{ zz0T*4XT1w$TOp>`LDRCx2rZw{bEcOmZ>A5e;o8A~J<1V5!f}xYBKa%Fk(4ggGzQk! z5ykQAZEo~a%i9f~w;jHZ_o8(+KG~c5ww_iw9)&T9K#MJK#IPwf{KGiY;^$}bhzXIG zxx&|nx!^YUksxcyy0LTAn%ms<{+UdhgfX@^-i@Zlxs1YNIaA`~GpX^qtwrZw&ka*o zeaDJ*3lZ`Pgk?bP3P%!NBc|ELkAwN*x$PX-XUBrYQ5zK(OoYq@nKe@_F3vlFon;j{w`|z<#Z_ zpqvC?-*hZ6e2FfVG&CVy9U@QhlCTH$!!){%IODM=9vR-SRMLSU_V`RTw@Q$GXP(#T$(>SP(sz6Oau!B1o~6MNSg65tk2N^>~|D)9zKs#>8KZ4s#Vt(*S_Uwm#*Y*^D+mB{Yn1%b6f)FwF4k$2$m! z%nI-b=d=>e(^$wrXp(C~Gw1aD%5J8#t#K8)Tepn}Y1xNO_Y1pDY96|lI^EqwK_)_F z3QK{0N=8jf1$ueQ!m#fWLKAvZ!H6>rX89KgqGxGJ&2lD6n$6nq7OZAX@W`>T3XxTo zFj6g3@)of)R=c~-rBsq}On`U_Q}kZEsAI2q+z{1#di(}E-4W@N?2y5@y}m^Otk6es zdwU_?YE|DR+*+U14;Js2+5`qfa7XIU(wB^ht~RYhCPjD{*vkMKsI*wdv$Re=uY(!B z8J0^e)UG3Dd7Fta8FMQRyxcY8*5;^!O`%>$ zNwm4VfUlZV^mNMagle#u>lXEJdfkwh)~cFMqyBVlV}Z+DeRHzQ>1!j$(p|>_GzO+Xwi` zKXrDB;WE6!WrA5gT5-WHUV1oDVN4*3I=FKP;W>`a+ z`Kr*T1uV`!^#fEEjmdC{8Sz?i6d$7tayr~OWhl+NoxgBSdKE;>E1HW;%TkqtirnWw zvAWP4&A0_(JhYl4IsR~U*}a=FIpd+iin(QBLKDWiTYFT1H^r`(Su}gKzCgtEkNeROKAS&+VZu8Nt@g-cXL38ncnGkD-XyqD|q=t?4ZQt>;(-KtD$s$ zUVD5&Sxb!KGdMO*4S)ap@IoisS@BjkX0sQJUonfN1! zkOF+IswCwovR_2KmX?&;#dQ{TpFt|)f^Vs|nL&~Y2+ z*IsJYhYM+>7G;Oof^ae34<^LftngfA;D&Ov_cr&Mb@(2w`j7%7pvu7pL2fOeswSX2 zcZ30zUpccH*uL(7{MIZm3EBP@a~xjHk2!tgA-DzrDmOx->am=Cd2_6Vit# zBJ6;BvgBYFYV1&F`s~9KbBz|oW54&MUs+N{b*H4@-wA0$sok7H#y^wad!$u2dt~u@ zvg^!Rp<&_hM{?H*xfJ%<3nGS|P8)r!wnt^S()yZwK1 zf^AF$p##V&uBHOwcoETJ@YJf;qbu`BJLXEalg$k^sf|m|?p-5+-}6fKeYrypcdnAw zhT9NDxzeicYt2p_5#vpbVl$)PTQTj+H^suS`|FY?5BWFiMFk=I;_RkzBdv_Z#%IKN z?855Q6*8=vUpz#Xh8(I+=kAk>*;k9nU*lQ7C%mK7eHwdbhn)iO=u;16Q!6_sP%?6^z-a za}A#hM_Lg2_G}gUS%+m)vKubIow0r;{FZ=R_rCT|Pe9h$7g|B9b7@Sj;d5sz2=ldg z;7O!}r}>zB<(@&5_lV@yr9L=3%YDDL+KpMQWazUK38NNt@{kl!XZlhReFObM-HvDV zAFl$%lneBkEbkgtD$IlJQ-{tx){d++4yb@+ZLN+0s5z3UYy$PH9mHu>)v?b@_U~5- zHyJA_8$x_a5go^kWkS{O)#pt&5AX|ob{PAmvcXu#n3ij>Z@PEFHz}}Zh5hu{%URa1 z)wxUakA4CtfzOq%Zy*0vr(HXxDzy8wpB;YcwTPdd%Aek)a`ra14)lgL_Qt=5+5FEv z`Sa05#%s%VGoS=tNWa2MUWhjMg^yr_VbMjBkr}Y8&)owgm490XHvCw(HUv(Dz8D{H zSkG#=Vp`ivYXCJg&Ltzk0qF(qgC?A_7g*nC_L1U_kb+D`g(CX_Bw2X0C*18c3JpwO zIH@P+fA>AUnsEWST40v3Ai0vj)OizPFFm4r7e?D6S&6ybL^4MnL2Cvl7BV<;-oO1E zbx23yQb)l?Y90bYhA;pJMFR^?z+Def%6bBhC6kybMa`Te!g(xfU41FQSRHOyKcM+d z@cHWXc3EpKrQ|*wSRA>_Ffo-seNym{sO9MyEZV9`R%;ZNfkx^XqF;%+7IqT#Tg!r6 zQV`;NRF)AOvR+dh)*5krH#w-N>F9Qi$6OGIohokTRfFmdSbg<7qPgoFfmhb)n-x0n z7N_p^d|$;4Z>vx0gcug6{4RKdH%nHCJkDA#|E2ncGm{U(=mYkvq8&<%H>ERBIZu%f z^V@&5=BmVP-9>(C!X9`407(DRcdebQ42%sXg0LV!ZIL$ZB8cqz_x9dQ~1!Ymz$P@F?sGeV+VfA#p#EHd}fDsW)9~xDOG{qWm2%9R0WJIZeVb&t# zvTxoI-*ox{JitnADbD)jtD3rQiDYSrmQhS{C0Zjgrh!(fE@0I=2@dC5joQov4H8`T zQ;r;^!n1M*mS^H-25TD2syTbl54qXkY1dJla-Q=FY-3NdqvK=nU_;b7d3%YC`-x`k zYvWR38XjjRLAjCNBln%rKBmDMx0-bbbypqJc1dWb`3u-oe_t_+2h23jcVe7N#w@}~ z(+DYbYmIe5o>zaoGe@7o4l?^S72hKh%-dq`{cqjXz(g?OQ%e*P0GUuUBjN zn4mLvnW#*6;r?bKJ&?p!^ENUf&X@Y#^rh0mB!J&$TibOT~p8^HqA(thYu6(%vR&-{v zS|dPh3E(PHPSjB&R(>Wk4MtP1{12=`Wxd(&afUN`jW>GTK|=WmVmbZ3YnYTAPqJhu zd!{|9P8oaJ6{3Ro|9rTW(WTrV(Z1lt{XaaMx&X2t<0^5IUGh11gTj3FC0wmLfLH7X8BNq8!kS;*j*xlc83WQ(*iN zNwMQQ8ad!k*sc-^ZL?zSn@aYoMZvlNJ;YNE6-U_+9_OW2({k|wO$;1r=~nJNG}v|$ z`9qm}$I^QDba$*_#-$hn&;V-#e_bgpY zAN9DI_nEiS}(LPYG=yPO^lL&?uBfclD7yIkC+nv4_RrKtQiiZ0ArZ`r*Te*HC zju^3}?Nks7xVt2znK4r&s|gu6UG?oc_K9&4l~%V>B%3$ocY9XPpY68q!JMr55^p!h ztC~&`gTcq+e6746w{xyy3^>Owz>MBrT}G|{5#K2K{fF?tP%daY|^*2{oktT^O602d{-Ee`3&l6 zCBA}}KQB+XmgCnXo`tg8lrIJxPi_;Z!C9$>&GVTEPDJR9`tEGny0a9vb)3(<^dmYOOG^JIB^1IU7XtEnxLX`xDsa?_qPxYNs_o#h(W`o`c!L#^s;oJ7@CXWV8prr*if!^{K7hPl6r!WCn_L> z$i{;s+q*D}Aim!&qLPW^TVHiXHuFjq?;rO#Pe6U8V*t*62|l4YF5dtx)&_mC8tOt? zl#8t+NGI_0xMX*(O6cB9N7}4oK=gY9sB+FkPu^Wri(CKZd7ZR!tB)h%9gsk)>uBMk zBUUX(H#~FHjtu1-n3o4C0#6LkaI5KmM4zL+lJm@alX8gciMLKpKXCx+r5opESo`>B z`$+!z?Eg7p0fA^g4ORa=rv6`V=wH|W5@9bV`Bwsejr#o;@V9IB=gRC)F~ENZ{xux( z&%lmP2fhDaaOB@9{XHV?4_aW5|0O=|@9@9-1pWaR{Iu!)gOR^`2>zYK-+e0nAR+Q8 z{r*AX&mI