{ "cells": [ { "cell_type": "markdown", "id": "78ef79ac", "metadata": {}, "source": [ "# 在量子化学计算中应用量子变分求解器\n", "\n", "[![下载Notebook](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/website-images/r2.5.0/resource/_static/logo_notebook.svg)](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/notebook/r2.5.0/mindquantum/zh_cn/case_library/mindspore_vqe_for_quantum_chemistry.ipynb) [![下载样例代码](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/website-images/r2.5.0/resource/_static/logo_download_code.svg)](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/notebook/r2.5.0/mindquantum/zh_cn/case_library/mindspore_vqe_for_quantum_chemistry.py) [![查看源文件](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/website-images/r2.5.0/resource/_static/logo_source.svg)](https://gitee.com/mindspore/docs/blob/r2.5.0/docs/mindquantum/docs/source_zh_cn/case_library/vqe_for_quantum_chemistry.ipynb)\n", "\n", "## 概述\n", "\n", "量子化学，指的是运用量子力学的基本理论及方法，求解含时或定态薛定谔方程的数值解。在高性能计算机上进行量子化学模拟已成为研究材料的物理、化学性质的重要手段。然而，精确求解薛定谔方程具有指数级的复杂度，可模拟的化学体系规模严重受制于此。近年量子计算的发展为解决这个问题提供了一条可行的路，有望在量子计算机上实现多项式复杂度下对薛定谔方程的高精度求解。\n", "\n", "[Peruzzo等人](https://doi.org/10.1038/ncomms5213)在2014年首次将量子变分求解器(Variational quantum eigensolver, VQE)结合[幺正耦合簇理论](https://linkinghub.elsevier.com/retrieve/pii/S0009261489873725)用于量子化学的模拟中，实现了He-H+基态能量的求解。量子变分求解器是一个量子——经典混合算法，在基于量子算法的化学模拟中应用广泛，本教程将介绍使用量子变分求解器求解分子体系基态能量的方法。\n", "\n", "本教程的主要内容包括如下几个部分：\n", "\n", "1. 量子化学原理简介。\n", "2. 量子变分求解器的应用。\n", "3. 使用MindSpore Quantum实现高效自动求导的VQE模拟。\n", "\n", "> 本文档适用于CPU环境。\n", "\n", "## 环境准备\n", "\n", "本教程需要安装以下环境：\n", "\n", "- NumPy\n", "- SciPy\n", "- [mindquantum](https://gitee.com/mindspore/mindquantum)\n", "- [mindspore](https://gitee.com/mindspore/mindspore)\n", "- PySCF\n", "- openfermion\n", "- openfermionpyscf\n", "\n", "> 以上依赖都可通过`pip`命令来安装。\n", "\n", "PySCF、openfermion这两个化学模拟的包只支持在Linux系统或Mac系统上安装，Windows 10系统则需要在Ubuntu子系统中才能安装。\n", "\n", "## 导入依赖\n", "\n", "导入本教程所依赖模块" ] }, { "cell_type": "code", "execution_count": 1, "id": "0fa38838", "metadata": { "scrolled": true }, "outputs": [], "source": [ "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\n", "from mindquantum.core.operators import Hamiltonian\n", "from mindquantum.simulator import Simulator\n", "from mindquantum.algorithm.nisq import generate_uccsd\n", "import mindspore as ms\n", "\n", "ms.set_context(mode=ms.PYNATIVE_MODE, device_target=\"CPU\")" ] }, { "cell_type": "markdown", "id": "3ca57bba", "metadata": {}, "source": [ "## 量子化学计算方法\n", "\n", "量子化学的核心问题在于求解薛定谔方程（Schrödinger Equation）。一般来说，求解含时薛定谔方程（Time-dependent Schrödinger Equation）较为复杂，故引入玻恩-奥本海默近似（Born-Oppenheimer approximation, BO approximation）。BO近似认为，原子核质量远大于电子、运动速度远低于电子，故可以将两者进行分离变量，单独讨论原子核或电子的运动，于是可得到如下不含时的电子运动方程，也称为定态薛定谔方程：\n", "\n", "$$\n", "\\hat{H} |\\Psi\\rangle = E |\\Psi\\rangle\n", "$$\n", "\n", "其中$\\hat{H}$包含以下三项：\n", "\n", "$$\n", "\\hat{H} = \\hat{K} _{e} + \\hat{V} _{ee} + \\hat{V} _{Ne}\n", "$$\n", "\n", "分别为电子动能、电子-电子势能和电子-核势能。\n", "\n", "有多种数值算法可以求解定态薛定谔方程。本教程将介绍其中的一类：波函数方法。波函数方法直接求解给定分子哈密顿量的本征波函数和本征能量，目前有大量的开源软件包可实现，如[PySCF](http://pyscf.org/)等。此处从一个简单的例子：氢化锂分子开始，使用openfermion结合openfermionpyscf插件进行。首先定义分子结构：" ] }, { "cell_type": "code", "execution_count": 2, "id": "4d322d8a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Geometry: \n", " [['Li', [0.0, 0.0, 0.0]], ['H', [0.0, 0.0, 1.5]]]\n" ] } ], "source": [ "dist = 1.5\n", "geometry = [\n", " [\"Li\", [0.0, 0.0, 0.0 * dist]],\n", " [\"H\", [0.0, 0.0, 1.0 * dist]],\n", "]\n", "basis = \"sto3g\"\n", "spin = 0\n", "print(\"Geometry: \\n\", geometry)" ] }, { "cell_type": "markdown", "id": "cda1a0b2", "metadata": {}, "source": [ "上面的代码定义了一个Li-H键长为1.5Å分子。使用STO-3G基组进行计算。接下来使用openfermionpyscf，调用PySCF进行HF、CCSD和FCI计算。这三种方法属于波函数方法，开始计算之前，先对这些方法作一个简单的介绍。\n", "\n", "### 波函数方法\n", "\n", "求解定态薛定谔方程的方法之一是[Hartree-Fock（HF）](https://doi.org/10.1098/rspa.1935.0085)方法，该方法在二十世纪三十年代左右由Hartree等人提出，是量子化学计算中的基本方法。HF方法引入了单行列式近似，即 $N$ 电子体系的波函数由一个行列式形式的波函数表示：\n", "\n", "$$\n", "| \\Psi \\rangle = | \\psi_{1} \\psi_{2} \\psi_{3} \\dots \\psi_{N} \\rangle\n", "$$\n", "\n", "其中 $| \\psi_{1} \\psi_{2} \\psi_{3} \\dots \\rangle$ 代表由一组自旋轨道波函数 $\\{ \\pi_{i} \\}$ 构成的N阶行列式。\n", "自旋轨道波函数 $\\psi_{i}$ 可进一步用一组形式已知的基函数展开：\n", "\n", "$$\\psi_{i} = \\phi_{i} \\eta_{i}$$\n", "\n", "$$\\phi_{i} = \\sum_{\\mu}{C_{\\mu i} \\chi_{\\mu}}$$\n", "\n", "其中 $\\{\\chi_{\\mu}\\}$ 被称为基函数，可以是高斯函数等。\n", "该近似考虑了电子间的交换作用，但是忽略了电子间的关联作用，故无法正确计算如解离能等性质。\n", "\n", "HF方法的改进可以从波函数展开定理出发。波函数展开定理可以表述为，若 $\\{ \\psi_{i} \\}$ 是一组完备的自旋轨道波函数，则 $N$ 电子体系波函数可以由 $\\{ \\psi_{i} \\}$ 构成的行列式波函数精确展开：\n", "\n", "$$\n", "| \\Psi \\rangle = \\sum^{\\infty} _ {i_{1} < i_{2} < \\dots < i_{N}} {C_{i_{1} i_{2} \\dots i_{N}} | \\psi_{i_{1}} \\psi_{i_{2}} \\dots \\psi_{i_{N}} \\rangle}\n", "$$\n", "\n", "由此可得到Configuration Interaction（CI）方法：\n", "\n", "$$\n", "| \\Psi_{CI} \\rangle = C_{0} | \\Psi_{HF} \\rangle + \\sum^{a\\rightarrow\\infty} _{i\\in occ, a\\not\\in occ}{C^{a} _{i} | \\Psi^{a} _{i} \\rangle } + \\sum^{ab\\rightarrow\\infty} _{ij\\in occ, ab\\not\\in occ}{C^{ab} _{ij} | \\Psi^{ab} _{ij} \\rangle }\n", "$$\n", "\n", "上式中的 $| \\Psi^{a}_{i} \\rangle + \\dots$ 代表电子由轨道 $i$ 激发到轨道 $a$ 的单激发波函数，以此类推。只考虑单激发和双激发的CI被称为CISD，即Configuration Interaction with singles and doubles。将基态HF波函数一直到N激发波函数全部考虑在内的Configuration Interaction被称为Full Configuration Interaction（FCI），FCI波函数是定态薛定谔方程在给定基函数下的精确解。\n", "\n", "### 二次量子化\n", "\n", "在二次量子化表述下，体系的哈密顿量具有如下形式：\n", "\n", "$$\n", "\\hat{H} = \\sum_{p, q}{h^p_q E^p_q} + \\sum_{p, q, r, s}{\\frac{1}{2} g^{pq} _ {rs} E^{pq}_{rs}}\n", "$$\n", "\n", "其中 $E^p_q$ 和 $E^{pq}_{rs}$ 分别为：\n", "\n", "$$\n", "\\begin{align*}\n", "E^{pq} _{rs} &= a^{\\dagger} _{p} a^{\\dagger} _{q} a _ {r} a _ {s}\\\\\n", "E^p_q &= a^{\\dagger}_pa_q\n", "\\end{align*}\n", "$$\n", "\n", "$a^{\\dagger}_p$ 和 $a_q$ 分别为产生算符（Creation Operator）和湮灭算符（Annihilation Operator）。\n", "\n", "使用二次量子化的表述方法，可以非常方便地表示激发态波函数：\n", "\n", "$$\n", "| \\Psi^{abc\\dots} _ {ijk\\dots} \\rangle = a^{\\dagger} _ {a} a^{\\dagger} _ {b} a^{\\dagger} _ {c} \\dots a _ {i} a_{j} a_{k} \\dots | \\Psi \\rangle\n", "$$\n", "\n", "CI方法的一个改进是耦合簇理论（Coupled-Cluster theory, CC）。CC引入指数化算符：\n", "\n", "$$\n", "| \\Psi_{CC} \\rangle = \\exp{(\\hat{T})} | \\Psi_{HF} \\rangle\n", "$$\n", "\n", "其中耦合簇算符 $\\hat{T}$ 为对激发算符的求和：\n", "\n", "$$\n", "\\hat{T} = \\sum_{p\\not\\in occ, q\\in occ}{\\theta^{p} _ {q} E^{p} _ {q}} + \\sum_{pq\\not\\in occ, rs\\in occ}{\\theta^{pq} _ {rs} E^{pq} _ {rs}} + \\dots\n", "$$\n", "\n", "其中 $\\theta$ 和CI方法中的 $C$ 类似，是待求解的参数。由指数的泰勒展开易知，即使耦合簇算符 $\\hat{T}$ 中只包含低阶激发项，$\\exp{(\\hat{T})}$ 也可以隐含部分高阶激发，这也使得CC方法向FCI波函数收敛的速度要远快于CI，同样截断到K激发，如K=2，CCSD的精度会超过CISD。\n", "\n", "\n", "\n", "电子关联作用的效果是使得总能量降低，故HF得到的基态能量会略高于CCSD和FCI。另外，从上述理论不难发现，FCI的计算量远大于CCSD和HF。我们使用openfermion封装的`MolecularData`和openfermionpyscf封装的`run_pyscf`函数来进行演示：" ] }, { "cell_type": "code", "execution_count": 3, "id": "5e6b2fab", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Hartree-Fock energy: -7.8633576215351200 Ha\n", "CCSD energy: -7.8823529091526812 Ha\n", "FCI energy: -7.8823622867987284 Ha\n" ] } ], "source": [ "molecule_of = MolecularData(\n", " geometry,\n", " basis,\n", " multiplicity=2 * spin + 1,\n", " data_directory='./'\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": "a3e28ac4", "metadata": {}, "source": [ "在上面的例子中，我们运行了Hartree-Fock（HF）、CCSD、FCI进行总能量的计算。若对运行时间进行统计，会发现$T_{HF}E_{CCSD}>E_{FCI}$。计算完成后，我们将结果保存到`molecule_file`文件（即`molecule_of.filename`）中：" ] }, { "cell_type": "code", "execution_count": 4, "id": "d983c107", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "H1-Li1_sto3g_singlet\n" ] } ], "source": [ "molecule_of.save()\n", "molecule_file = molecule_of.filename\n", "print(molecule_file.split('/')[-1])" ] }, { "cell_type": "markdown", "id": "2faf9510", "metadata": {}, "source": [ "量子化学计算的一大阻碍是计算量。随着体系大小（电子数、原子数）的增加，求解FCI波函数和基态能量的时间消耗大约以$2^{N}$增长，即使是较小的分子如乙烯分子等，进行FCI计算也并不容易。量子计算机的出现为此提供了一条可能的解决途径，已有的研究表明，量子计算机可以多项式的时间复杂度模拟哈密顿量的含时演化，在量子处理器上进行化学模拟相较于经典计算机有指数级的加速。本教程将介绍其中一类量子算法：量子变分求解器。\n", "\n", "## 量子变分求解器\n", "\n", "量子变分求解器（Variational Quantum Eigensolver, VQE）是一类量子-经典混合（Hybrid quantum-classical）算法，应用变分原理实现对基态波函数的求解。其中，变分参数的优化步在经典计算机上进行。\n", "\n", "### 变分原理\n", "\n", "变分原理可使用如下形式表述：\n", "\n", "$$\n", "E_{0} \\le \\frac{\\langle \\Psi_{t} | \\hat{H} | \\Psi_{t} \\rangle}{\\langle \\Psi_{t} | \\Psi_{t} \\rangle}\n", "$$\n", "\n", "上式中的$| \\Psi_{t} \\rangle$代表试探波函数。变分原理表明，在满足一定的条件下，任意试探波函数得到的基态能量总是大于等于真实的基态能量。变分原理为求解分子基态薛定谔方程提供了一种方法：使用一个参数化的函数$f(\\theta)$作为精确基态波函数的近似，通过优化参数$\\theta$来逼近精确的基态能量。\n", "\n", "### 初态制备\n", "\n", "在二次量子化表述下，$N$-电子HF波函数也具有非常简洁的形式：\n", "\n", "$$\n", "| \\Psi_{HF} \\rangle = \\prod^{i=0} _{N-1}{a^{\\dagger} _{i}| 0 \\rangle}\n", "$$\n", "\n", "上式搭建了一个由量子化学波函数到量子计算的桥梁：用$|0\\rangle$代表非占据轨道，用$|1\\rangle$代表电子占据的轨道，由此可以将$N$-电子HF波函数映射为由一串$M+N$个量子比特$| 00\\dots 11\\dots \\rangle$，$M$代表非占据轨道的数量。\n", "\n", "以下代码构造了对应于LiH分子的HF初态波函数。在Jordan-Wigner变换下，相当于将$N$个$\\text{X}$门作用于$|000\\dots\\rangle$上。" ] }, { "cell_type": "code", "execution_count": 5, "id": "32a3cafd", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " ┏━━━┓ \n", "q0: ──┨╺╋╸┠───\n", " ┗━━━┛ \n", " ┏━━━┓ \n", "q1: ──┨╺╋╸┠───\n", " ┗━━━┛ \n", " ┏━━━┓ \n", "q2: ──┨╺╋╸┠───\n", " ┗━━━┛ \n", " ┏━━━┓ \n", "q3: ──┨╺╋╸┠───\n", " ┗━━━┛ \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": "eab886b7", "metadata": {}, "source": [ "基于此，我们可以构造如下形式的试探波函数：\n", "\n", "$$\n", "| \\Psi_{t} \\rangle = U(\\theta) | \\Psi_{HF} \\rangle\n", "$$\n", "\n", "其中$U(\\theta)$代表一个可通过量子线路模拟的幺正变换，$| \\Psi_{HF} \\rangle$作为初态，可通过多个单比特$\\text{X}$门来方便地制备。$U(\\theta) | \\Psi_{HF} \\rangle$的具体形式也被称为波函数拟设。\n", "\n", "### 波函数拟设\n", "\n", "前文提到的耦合簇理论是一个非常高效的波函数拟设。在量子计算机上使用，需要作一些修改：\n", "\n", "$$\n", "| \\Psi_{UCC} \\rangle = \\exp{(\\hat{T} - \\hat{T}^{\\dagger})} | \\Psi_{HF} \\rangle\n", "$$\n", "\n", "UCC即幺正耦合簇(Unitary Coupled-Cluster theory)，$\\hat{T}^{\\dagger}$代表$\\hat{T}$的厄米共轭。如此，$\\exp{(\\hat{T} - \\hat{T}^{\\dagger})}$即为幺正算符。[Peruzzo等人](https://doi.org/10.1038/ncomms5213)在2014年首次使用VQE结合UCCSD(Unitary coupled-cluster with singles and doubles)拟设进行了量子计算机上的化学模拟实验。值得注意的是幺正耦合簇默认了耦合簇算符中的参数$\\{\\theta\\}$是实数。在分子体系中该假设不会有问题；在周期性体系中，[刘杰等人](https://doi.org/10.1021/acs.jctc.0c00881)的研究表明幺正耦合簇会因为忽略复数部分而造成误差。本教程暂时不讨论幺正耦合簇在周期性体系中的应用。" ] }, { "cell_type": "markdown", "id": "d5da4c30", "metadata": {}, "source": [ "使用mindquantum的circuit模块中的`generate_uccsd`函数可读取先前保存在`molecule_file`的计算结果，“一键”构造UCCSD波函数拟设，以及其对应的量子线路：" ] }, { "cell_type": "code", "execution_count": 6, "id": "002f8dd8", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "ccsd:-7.882352909152681.\n", "fci:-7.882362286798728.\n" ] } ], "source": [ "ansatz_circuit, \\\n", "init_amplitudes, \\\n", "ansatz_parameter_names, \\\n", "hamiltonian_QubitOp, \\\n", "n_qubits, n_electrons = generate_uccsd(molecule_file, threshold=-1)" ] }, { "cell_type": "markdown", "id": "6082cca2", "metadata": {}, "source": [ "[generate_uccsd](https://www.mindspore.cn/mindquantum/docs/zh-CN/r0.10/algorithm/nisq/mindquantum.algorithm.nisq.generate_uccsd.html) 将幺正耦合簇相关的函数打包了起来，包括导出分子哈密度量、构造幺正耦合簇拟设算符、提取CCSD计算的耦合簇系数等多个步骤。该函数通过输入分子的文件路径来读取该分子，参数`th`是表示量子线路中哪些参数需要更新梯度的阈值。在[分步构造幺正耦合簇拟设](#分步构造幺正耦合簇拟设)章节，我们会演示如何使用mindquantum的相关接口分步完成其中包含的步骤。完整的量子线路包含HF初态+UCCSD拟设，如下代码所示：" ] }, { "cell_type": "code", "execution_count": 7, "id": "d2210b4f", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
                            Circuit Summary                            \n",
       "╭──────────────────────┬──────────────────────────────────────────────╮\n",
       "│ Info                 │ value                                        │\n",
       "├──────────────────────┼──────────────────────────────────────────────┤\n",
       "│ Number of qubit      │ 12                                           │\n",
       "├──────────────────────┼──────────────────────────────────────────────┤\n",
       "│ Total number of gate │ 12612                                        │\n",
       "│ Barrier              │ 2560                                         │\n",
       "│ Noise Channel        │ 0                                            │\n",
       "│ Measurement          │ 0                                            │\n",
       "├──────────────────────┼──────────────────────────────────────────────┤\n",
       "│ Parameter gate       │ 640                                          │\n",
       "│ 44 ansatz parameters │ p0, p8, p1, p9, p2, p10, p3, p11, p4, p12... │\n",
       "╰──────────────────────┴──────────────────────────────────────────────╯\n",
       "
\n" ], "text/plain": [ "\u001b[1;38;2;255;0;0m Circuit Summary \u001b[0m\n", "╭──────────────────────┬──────────────────────────────────────────────╮\n", "│\u001b[1m \u001b[0m\u001b[1;38;2;59;59;149mInfo\u001b[0m\u001b[1m \u001b[0m\u001b[1m \u001b[0m│\u001b[1m \u001b[0m\u001b[1;38;2;59;59;149mvalue\u001b[0m\u001b[1m \u001b[0m\u001b[1m \u001b[0m│\n", "├──────────────────────┼──────────────────────────────────────────────┤\n", "│ \u001b[1mNumber of qubit\u001b[0m │ 12 │\n", "├──────────────────────┼──────────────────────────────────────────────┤\n", "│ \u001b[1mTotal number of gate\u001b[0m │ 12612 │\n", "│ Barrier │ 2560 │\n", "│ Noise Channel │ 0 │\n", "│ Measurement │ 0 │\n", "├──────────────────────┼──────────────────────────────────────────────┤\n", "│ \u001b[1mParameter gate\u001b[0m │ 640 │\n", "│ 44 ansatz parameters │ \u001b[38;2;72;201;176mp0, p8, p1, p9, p2, p10, p3, p11, p4, p12...\u001b[0m │\n", "╰──────────────────────┴──────────────────────────────────────────────╯\n" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Number of parameters: 44\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": "b1eab60c", "metadata": {}, "source": [ "对于LiH分子而言，其UCCSD波函数拟设中包含44个变分参数。该线路总共的量子比特门数量为12612，总共需要12个量子比特进行模拟。" ] }, { "cell_type": "markdown", "id": "06c43609", "metadata": {}, "source": [ "### VQE的一般流程\n", "\n", "使用VQE进行分子基态求解的一般流程如下：\n", "\n", "1. 制备HF初态：$| 00\\dots11\\dots \\rangle$；\n", "2. 定义波函数拟设，如UCCSD等；\n", "3. 将波函数拟设转化为参数化的量子线路；\n", "4. 初始化变分参数，如全设为0等；\n", "5. 在量子计算机上多次测量得到分子哈密顿量在该套变分参数下的能量$E(\\theta)$以及能量关于参数的导数$\\{ {\\partial E} / {\\partial \\theta_{i}} \\}$\n", "6. 在经典计算机上使用优化算法，如梯度下降、BFGS等更新变分参数；\n", "7. 将新的变分参数传入量子线路中进行更新；\n", "8. 重复步骤(5)到(7)，直到满足收敛标准；\n", "9. 结束\n", "\n", "在第5步中，求取能量关于参数的导数$\\{ {\\partial E} / {\\partial \\theta_{i}} \\}$在量子计算机上可通过parameter-shift rule来进行，在模拟器中也可通过模拟parameter-shift rule或者有限差分法来计算，是个较为耗时的过程。mindquantum基于mindspore框架，提供了类似于机器学习的自动求导功能，可以在模拟中可以高效计算变分量子线路的导数。以下使用mindquantum构造带自动求导功能的参数化UCCSD量子线路：" ] }, { "cell_type": "code", "execution_count": 8, "id": "16b5d1c3", "metadata": {}, "outputs": [], "source": [ "sim = Simulator('mqvector', total_circuit.n_qubits)\n", "molecule_pqc = sim.get_expectation_with_grad(Hamiltonian(hamiltonian_QubitOp), total_circuit)" ] }, { "attachments": {}, "cell_type": "markdown", "id": "6c47ddee", "metadata": {}, "source": [ "通过将参数的具体数值传入`molecule_pqc`，即可得到对应于此变分参数的能量$E(\\theta)=\\langle \\Psi_{UCC}(\\theta) | \\hat{H} | \\Psi_{UCC}(\\theta) \\rangle$以及关于每个变分参数的导数。\n", "\n", "例如，我们可以利用下面的代码计算当变分量子线路中所有参数都为零时，哈密顿量的期望值和期望值关于梯度的导数。" ] }, { "cell_type": "code", "execution_count": 9, "id": "46d4958c", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Energy: [[-7.86335762+0.j]] \n", "shape: (1, 1) \n", "\n", "Gradient: [[[-5.76513698e-11+0.j -8.60518134e-02+0.j 1.19273071e-08+0.j\n", " -4.85545124e-02+0.j -6.07652177e-15+0.j -3.92769093e-02+0.j\n", " -1.00082540e-14+0.j -9.59481745e-02+0.j 3.21896998e-15+0.j\n", " -3.92769093e-02+0.j -2.06907796e-15+0.j -9.59481745e-02+0.j\n", " -8.59782560e-11+0.j -2.89649669e-02+0.j 2.97637789e-10+0.j\n", " -4.91813233e-01+0.j -9.35293242e-04+0.j -1.69582252e-16+0.j\n", " -3.24387386e-17+0.j 1.60229175e-17+0.j 1.54862635e-16+0.j\n", " 2.66601400e-17+0.j -2.14277745e-17+0.j 5.06813117e-03+0.j\n", " 1.08542346e-02+0.j -1.28614265e-02+0.j -5.51985604e-17+0.j\n", " 9.81307914e-18+0.j 1.33973787e-01+0.j -3.03063680e-02+0.j\n", " 1.83995210e-18+0.j -1.83995210e-18+0.j -2.21362064e-29+0.j\n", " 5.61185390e-17+0.j -3.80815167e-17+0.j 3.33490478e-17+0.j\n", " 1.14905747e-28+0.j 1.15391965e-16+0.j -3.03063680e-02+0.j\n", " -5.51985630e-17+0.j 4.93146033e-17+0.j -4.22508287e-17+0.j\n", " -2.44510199e-17+0.j -1.68035039e-03+0.j]]] \n", "shape: (1, 1, 44)\n" ] } ], "source": [ "import numpy as np\n", "\n", "n_params = len(total_circuit.params_name)\n", "p0 = np.zeros(n_params)\n", "f, g = molecule_pqc(p0)\n", "print(\"Energy: \", f, \"\\nshape: \", f.shape, '\\n')\n", "print(\"Gradient: \", g, \"\\nshape: \", g.shape)" ] }, { "attachments": {}, "cell_type": "markdown", "id": "f2c4c562", "metadata": {}, "source": [ "通过上述计算，我们获得了能量值和梯度值，用户可以根据实际应用需求来提取这些值。\n", "接下来需要进行VQE优化的(5)~(7)步，即对变分量子线路进行优化。我们可以借助scipy中的优化器来对线路参数进行优化。首先我们构造适用与scipy优化器的优化函数：" ] }, { "cell_type": "code", "execution_count": 10, "id": "2bd09b88", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(-7.86335762153696,\n", " array([-5.76513698e-11, -8.60518134e-02, 1.19273071e-08, -4.85545124e-02,\n", " -6.07652177e-15, -3.92769093e-02, -1.00082540e-14, -9.59481745e-02,\n", " 3.21896998e-15, -3.92769093e-02, -2.06907796e-15, -9.59481745e-02,\n", " -8.59782560e-11, -2.89649669e-02, 2.97637789e-10, -4.91813233e-01,\n", " -9.35293242e-04, -1.69582252e-16, -3.24387386e-17, 1.60229175e-17,\n", " 1.54862635e-16, 2.66601400e-17, -2.14277745e-17, 5.06813117e-03,\n", " 1.08542346e-02, -1.28614265e-02, -5.51985604e-17, 9.81307914e-18,\n", " 1.33973787e-01, -3.03063680e-02, 1.83995210e-18, -1.83995210e-18,\n", " -2.21362064e-29, 5.61185390e-17, -3.80815167e-17, 3.33490478e-17,\n", " 1.14905747e-28, 1.15391965e-16, -3.03063680e-02, -5.51985630e-17,\n", " 4.93146033e-17, -4.22508287e-17, -2.44510199e-17, -1.68035039e-03]))" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def fun(p0, molecule_pqc, energy_list=None):\n", " f, g = molecule_pqc(p0)\n", " f = np.real(f)[0, 0]\n", " g = np.real(g)[0, 0]\n", " if energy_list is not None:\n", " energy_list.append(f)\n", " if len(energy_list) % 5 == 0:\n", " print(f\"Step: {len(energy_list)},\\tenergy: {f}\")\n", " return f, g\n", "\n", "fun(p0, molecule_pqc)" ] }, { "attachments": {}, "cell_type": "markdown", "id": "1d64afd0", "metadata": {}, "source": [ "此时，我们定义出来的 `fun` 函数能够正确返回符合要求的数据类型：一个实数格式的能量值，一个跟参数长度一直的梯度数组。接下来，我们利用 scipy 中的 `bfgs` 二阶优化器来进行优化。" ] }, { "cell_type": "code", "execution_count": 11, "id": "a2b6dea1", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Step: 5,\tenergy: -7.880227726111784\n", "Step: 10,\tenergy: -7.88181712406493\n", "Step: 15,\tenergy: -7.882213242904179\n", "Step: 20,\tenergy: -7.882345336993505\n", "Step: 25,\tenergy: -7.882352494990301\n", "Step: 30,\tenergy: -7.882352691272338\n", "Step: 35,\tenergy: -7.882352706023384\n" ] } ], "source": [ "from scipy.optimize import minimize\n", "\n", "energy_list = []\n", "res = minimize(fun, p0, args=(molecule_pqc, energy_list), method='bfgs', jac=True)" ] }, { "attachments": {}, "cell_type": "markdown", "id": "b69ea772", "metadata": {}, "source": [ "如上，我们便完成了变分量子线路的梯度优化。在这里 `energy_list` 用于存储收敛过程中的能量。这里我们对 `minimize` 函数做一个简单介绍：\n", "\n", "- `fun`: 第一个参数表示想要优化的函数\n", "- `p0`：第二个参数表示变量的初始值\n", "- `args`：`fun` 函数中除了第一个参数以外的其他参数，根据我们的 `fun` 函数的定义，这里选择 `args=(molecule_pqc, energy_list)`\n", "- `method`：所使用的优化算法，我们这里选择二阶优化算法 `bfgs`。更多优化算法请参考 [scipy 官方文档](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html)\n", "- `jac`：`fun` 函数是否返回导数值，这里我们选择 `True`，因为 MindSpore Quantum 框架能够算出变分量子线路中参数的精确梯度值。如选 `False`，则 `minimize` 框架将在内部利用差分方法计算近似梯度值。\n", "\n", "`res` 为 `scipy` 中优化算法得到的优化结果，包括优化得到的参数，最优化函数值和迭代次数等。" ] }, { "cell_type": "code", "execution_count": 12, "id": "82e40e12", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Ground state: \n", "-7.8823527083497105\n", "\n", "FCI: \n", "-7.882362286798721\n", "\n", "Optimized amplitudes: \n", "[ 2.38544548e-04 1.89071786e-03 3.52375725e-02 1.60369224e-02\n", " -5.38925036e-09 9.09472571e-04 -1.56807151e-10 1.41642144e-02\n", " 1.51399639e-08 9.08742805e-04 -1.74968590e-10 1.41698985e-02\n", " -5.47895177e-04 4.26807722e-04 2.87128385e-03 5.38110483e-02\n", " 2.34823599e-04 1.08613207e-08 -1.31353895e-07 1.29567542e-07\n", " -3.56002770e-08 3.17328351e-07 -3.10806922e-07 1.32372340e-05\n", " -1.04187318e-04 7.98992314e-04 4.80948977e-10 -1.09620509e-09\n", " -5.50006189e-02 3.09112335e-03 2.82662610e-10 1.68729874e-09\n", " -1.70596920e-09 -1.22679172e-08 -1.60711122e-08 1.29352307e-08\n", " 8.29267936e-13 -5.33349675e-11 3.09111737e-03 4.61835836e-08\n", " -5.41032067e-08 5.57952993e-08 -9.64149879e-11 3.72726345e-04]\n" ] } ], "source": [ "print(f\"Ground state: \\n{res.fun}\\n\")\n", "print(f\"FCI: \\n-7.882362286798721\\n\")\n", "print(f\"Optimized amplitudes: \\n{res.x}\")" ] }, { "attachments": {}, "cell_type": "markdown", "id": "9d8fcd7f", "metadata": {}, "source": [ "可以看到，幺正耦合簇给出的计算结果和FCI非常接近，具有良好的精度。" ] }, { "cell_type": "markdown", "id": "d435c8c6", "metadata": {}, "source": [ "## 分步构造幺正耦合簇拟设\n", "\n", "" ] }, { "cell_type": "markdown", "id": "61acab6d", "metadata": {}, "source": [ "在上文中，我们使用了 [generate_uccsd](https://www.mindspore.cn/mindquantum/docs/zh-CN/r0.10/algorithm/nisq/mindquantum.algorithm.nisq.generate_uccsd.html) 一步构造出了幺正耦合簇拟设所需要的所有内容，此处我们将步骤拆分，分别得到我们需要的耦合簇算符、对应的量子线路以及取自于经典CCSD计算结果的变分参数初猜值。\n", "首先，导入部分额外依赖，主要包含mindquantum中hiqfermion模块的相关函数：" ] }, { "cell_type": "code", "execution_count": 13, "id": "69da8b15", "metadata": {}, "outputs": [], "source": [ "from mindquantum.algorithm.nisq import Transform\n", "from mindquantum.algorithm.nisq import get_qubit_hamiltonian\n", "from mindquantum.algorithm.nisq import uccsd_singlet_generator, uccsd_singlet_get_packed_amplitudes\n", "from mindquantum.core.operators import TimeEvolution\n" ] }, { "cell_type": "markdown", "id": "b758f50f", "metadata": {}, "source": [ "分子哈密顿量使用 [get_qubit_hamiltonian](https://www.mindspore.cn/mindquantum/docs/zh-CN/r0.10/algorithm/nisq/mindquantum.algorithm.nisq.get_qubit_hamiltonian.html) ，读取之前的计算结果得到：" ] }, { "cell_type": "code", "execution_count": 14, "id": "fd9d413f", "metadata": {}, "outputs": [], "source": [ "hamiltonian_QubitOp = get_qubit_hamiltonian(molecule_of)" ] }, { "cell_type": "markdown", "id": "aecb0ba1", "metadata": {}, "source": [ "对于幺正耦合簇算符 $\\hat{T} - \\hat{T}^{\\dagger}$ ，可以使用 [uccsd_singlet_generator](https://www.mindspore.cn/mindquantum/docs/zh-CN/r0.10/algorithm/nisq/mindquantum.algorithm.nisq.uccsd_singlet_generator.html) 进行构造。提供总量子比特数（总自旋轨道数）和总电子数，并设置参数`anti_hermitian=True`：" ] }, { "cell_type": "code", "execution_count": 15, "id": "aa2b81fa", "metadata": {}, "outputs": [], "source": [ "ucc_fermion_ops = uccsd_singlet_generator(\n", " molecule_of.n_qubits, molecule_of.n_electrons, anti_hermitian=True)" ] }, { "cell_type": "markdown", "id": "ca2213b0", "metadata": {}, "source": [ "上一步构造的`ucc_fermion_ops`是参数化的。使用Jordan-Wigner变换将费米子激发算符映射为Pauli算符：" ] }, { "cell_type": "code", "execution_count": 16, "id": "17b6c897", "metadata": {}, "outputs": [], "source": [ "ucc_qubit_ops = Transform(ucc_fermion_ops).jordan_wigner()" ] }, { "cell_type": "markdown", "id": "657c9f22", "metadata": {}, "source": [ "接下来，我们需要得到幺正算符 $\\exp \\left(\\hat{T}-\\hat{T}^{\\dagger}\\right)$ 所对应的量子线路。[TimeEvolution](https://www.mindspore.cn/mindquantum/docs/zh-CN/r0.10/core/operators/mindquantum.core.operators.TimeEvolution.html) 可生成 $\\exp (-i \\hat{H} t)$ 所对应的线路，其中 $\\hat{H}$ 是一个厄米算符， $t$ 是实数。需要注意的是，使用[TimeEvolution](https://www.mindspore.cn/mindquantum/docs/zh-CN/r0.10/core/operators/mindquantum.core.operators.TimeEvolution.html) 时，`ucc_qubit_ops`中已经包含了复数因子 $i$ ，所以我们需要将`ucc_qubit_ops`除以 $i$ ，或者提取其虚部：" ] }, { "cell_type": "code", "execution_count": 17, "id": "55aa2ceb", "metadata": {}, "outputs": [], "source": [ "ansatz_circuit = TimeEvolution(ucc_qubit_ops.imag, 1.0).circuit\n", "ansatz_parameter_names = ansatz_circuit.params_name" ] }, { "cell_type": "markdown", "id": "6c216bb3", "metadata": {}, "source": [ "我们使用`ansatz_parameter_names`记录该线路中的参数名。到目前为止，我们已经得到了VQE量子线路所需要内容，包括哈密顿量`hamiltonian_QubitOp`、参数化的波函数拟设线路`ansatz_circuit`，故可仿照前文，得到完整的态制备线路。其中Hartree-Fock参考态复用之前的`hartreefock_wfn_circuit`：" ] }, { "cell_type": "code", "execution_count": 18, "id": "61cf22a6", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
                                  Circuit Summary                                  \n",
       "╭──────────────────────┬──────────────────────────────────────────────────────────╮\n",
       "│ Info                 │ value                                                    │\n",
       "├──────────────────────┼──────────────────────────────────────────────────────────┤\n",
       "│ Number of qubit      │ 12                                                       │\n",
       "├──────────────────────┼──────────────────────────────────────────────────────────┤\n",
       "│ Total number of gate │ 12612                                                    │\n",
       "│ Barrier              │ 2560                                                     │\n",
       "│ Noise Channel        │ 0                                                        │\n",
       "│ Measurement          │ 0                                                        │\n",
       "├──────────────────────┼──────────────────────────────────────────────────────────┤\n",
       "│ Parameter gate       │ 640                                                      │\n",
       "│ 44 ansatz parameters │ s_0, d1_0, s_1, d1_1, s_2, d1_2, s_3, d1_3, s_4, d1_4... │\n",
       "╰──────────────────────┴──────────────────────────────────────────────────────────╯\n",
       "
\n" ], "text/plain": [ "\u001b[1;38;2;255;0;0m Circuit Summary \u001b[0m\n", "╭──────────────────────┬──────────────────────────────────────────────────────────╮\n", "│\u001b[1m \u001b[0m\u001b[1;38;2;59;59;149mInfo\u001b[0m\u001b[1m \u001b[0m\u001b[1m \u001b[0m│\u001b[1m \u001b[0m\u001b[1;38;2;59;59;149mvalue\u001b[0m\u001b[1m \u001b[0m\u001b[1m \u001b[0m│\n", "├──────────────────────┼──────────────────────────────────────────────────────────┤\n", "│ \u001b[1mNumber of qubit\u001b[0m │ 12 │\n", "├──────────────────────┼──────────────────────────────────────────────────────────┤\n", "│ \u001b[1mTotal number of gate\u001b[0m │ 12612 │\n", "│ Barrier │ 2560 │\n", "│ Noise Channel │ 0 │\n", "│ Measurement │ 0 │\n", "├──────────────────────┼──────────────────────────────────────────────────────────┤\n", "│ \u001b[1mParameter gate\u001b[0m │ 640 │\n", "│ 44 ansatz parameters │ \u001b[38;2;72;201;176ms_0, d1_0, s_1, d1_1, s_2, d1_2, s_3, d1_3, s_4, d1_4...\u001b[0m │\n", "╰──────────────────────┴──────────────────────────────────────────────────────────╯\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "total_circuit = hartreefock_wfn_circuit + ansatz_circuit\n", "total_circuit.summary()" ] }, { "attachments": {}, "cell_type": "markdown", "id": "f2988e2f", "metadata": {}, "source": [ "下一步，需要为变分参数提供一个合理的初始值。前文构造的优化器默认使用0作为初猜，在大多数情况下是可行的。不过，使用CCSD的计算数据作为UCC的出发点，可能会有更好的结果。使用 [uccsd_singlet_get_packed_amplitudes](https://www.mindspore.cn/mindquantum/docs/zh-CN/r0.10/algorithm/nisq/mindquantum.algorithm.nisq.uccsd_singlet_get_packed_amplitudes.html) 函数从`molecule_of`提取CCSD的参数：" ] }, { "cell_type": "code", "execution_count": 19, "id": "0a6c69ca", "metadata": {}, "outputs": [], "source": [ "init_amplitudes_ccsd = uccsd_singlet_get_packed_amplitudes(\n", " molecule_of.ccsd_single_amps, molecule_of.ccsd_double_amps, molecule_of.n_qubits, molecule_of.n_electrons)\n", "init_amplitudes_ccsd = [init_amplitudes_ccsd[param_i] for param_i in ansatz_parameter_names]" ] }, { "attachments": {}, "cell_type": "markdown", "id": "92507b54", "metadata": {}, "source": [ "根据之前类似的方式，我们可以利用 MindSpore Quantum 来得到梯度算子，并利用 scipy 来进行优化。" ] }, { "cell_type": "code", "execution_count": 20, "id": "87bd142e", "metadata": {}, "outputs": [], "source": [ "grad_ops = Simulator('mqvector', total_circuit.n_qubits).get_expectation_with_grad(\n", " Hamiltonian(hamiltonian_QubitOp.real),\n", " total_circuit)" ] }, { "cell_type": "markdown", "id": "b75268c1", "metadata": {}, "source": [ "使用`init_amplitudes_ccsd`（即CCSD计算的耦合簇系数）作为初始变分参数：" ] }, { "cell_type": "code", "execution_count": 21, "id": "0cb42dc9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Step: 5,\tenergy: -7.87822328311087\n", "Step: 10,\tenergy: -7.880288481802379\n", "Step: 15,\tenergy: -7.882035668318941\n", "Step: 20,\tenergy: -7.882302370884022\n", "Step: 25,\tenergy: -7.882349803535785\n", "Step: 30,\tenergy: -7.882352702053598\n", "Step: 35,\tenergy: -7.882352707882877\n", "Step: 40,\tenergy: -7.882352708334483\n" ] } ], "source": [ "energy_list = []\n", "res = minimize(fun, init_amplitudes_ccsd, args=(grad_ops, energy_list), method='bfgs', jac=True)" ] }, { "attachments": {}, "cell_type": "markdown", "id": "33ec7112", "metadata": {}, "source": [ "最后，我们可以得到优化的结果如下：" ] }, { "cell_type": "code", "execution_count": 22, "id": "5c3764a9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Ground state: \n", "-7.882352708353205\n", "\n", "FCI: \n", "-7.882362286798721\n", "\n", "Optimized amplitudes: \n", "[-2.38523911e-04 1.89069822e-03 -3.52373476e-02 1.60367877e-02\n", " 1.01252245e-09 9.09441478e-04 -3.44110570e-10 1.41641951e-02\n", " -3.55437460e-09 9.08694499e-04 1.51899574e-10 1.41696901e-02\n", " 5.47582129e-04 4.26818549e-04 -2.87197044e-03 5.38109082e-02\n", " 2.34573340e-04 9.61091660e-11 3.26047168e-08 -3.28565811e-08\n", " -3.01353861e-08 2.57229375e-08 -2.35406397e-08 1.32955992e-05\n", " 7.98917451e-04 -1.04228796e-04 -1.41632464e-10 -1.72636257e-10\n", " -5.50006376e-02 3.09132360e-03 -4.33208163e-10 5.56668226e-09\n", " -5.53725165e-09 -9.18193176e-08 3.36714462e-08 -2.47305566e-08\n", " -1.26595719e-12 -1.96392557e-10 3.09104547e-03 8.41723416e-08\n", " -3.28071662e-09 -2.20592853e-09 2.57550150e-10 3.72870642e-04]\n" ] } ], "source": [ "print(f\"Ground state: \\n{res.fun}\\n\")\n", "print(f\"FCI: \\n-7.882362286798721\\n\")\n", "print(f\"Optimized amplitudes: \\n{res.x}\")" ] }, { "attachments": {}, "cell_type": "markdown", "id": "d4589eb4", "metadata": {}, "source": [ "## 总结\n", "\n", "在本案例中，我们通过两种方法，利用 `scipy` 中的优化器，得到了LiH分子的基态能量。在第一种方法中，我们利用MindSpore Quantum打包好的 [generate_uccsd](https://www.mindspore.cn/mindquantum/docs/zh-CN/r0.10/algorithm/nisq/mindquantum.algorithm.nisq.generate_uccsd.html) 函数生成了能够解决该问题的量子神经网络，而在第二种方法中，我们一步一步的构造出了类似的梯度算子。最终得到的结果是一致的。" ] }, { "cell_type": "code", "execution_count": 23, "id": "1fe3c3cd", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", " \n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", " \n", " \n", "\n", "\n", "
SoftwareVersion
mindquantum0.9.11
scipy1.9.3
numpy1.23.5
SystemInfo
Python3.8.17
OSLinux x86_64
Memory16.62 GB
CPU Max Thread16
DateTue Jan 2 17:03:00 2024
\n" ], "text/plain": [ "" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from mindquantum.utils.show_info import InfoTable\n", "\n", "InfoTable('mindquantum', 'scipy', 'numpy')" ] } ], "metadata": { "kernelspec": { "display_name": "MindSpore", "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.8.13" } }, "nbformat": 4, "nbformat_minor": 5 }