{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# ADAPT-VQE algorithm\n", "\n", "\n", "This tutorial explains the implementation of Adaptive Derivative-Assembled Pseudo-Trotter VQE (ADAPT-VQE) algorithm introduced in this [paper](https://www.nature.com/articles/s41467-019-10988-2). \n", "\n", "In VQE (see this [tutorial](https://nvidia.github.io/cuda-quantum/latest/applications/python/vqe_advanced.html)), a parameterized wave-function using UCCSD ansatz is generated and variationally tuned to minimize the expectation value of the molecular electronic Hamiltonian. In VQE approach, we include all possible single and double excitations of electrons from the occupied spin molecular orbitals of a reference state (Hartree Fock) to the unoccupied spin molecular orbitals. The excessive depth of these quantum circuits make them ill-suited for applications in the NISQ regime. \n", "\n", "The VQE issue has led to the ADAPT-VQE proposal in which the ansatz wave-functions is constructed through the action of a selective subset of possible unitary operators , i.e., only those operators whose inclusion in the ansatz can potentially lead to the largest decrease in the expectation value of the molecular electronic Hamiltonian. In ADAPT-VQE, the ansatz is grown iteratively by appending a sequence of unitary operators to the reference Hartree-Fock state. At each iteration, the unitary operator to be applied is chosen according to a simple criterion based on the gradient of the expectation value of the Hamiltonian. Therefore, allowing us to build a compact quantum circuit which can lead to more efficient use of quantum resources.\n", "\n", "The ADAPT-VQE algorithm consists of 8 steps:\n", "\n", "1- On classical hardware, compute one- and two-electron integrals, and transform the fermionic Hamiltonian into a qubit representation using an appropriate transformation: Jordan–Wigner, Bravyi–Kitaev, etc. For this tutorial, we will use Jordan Wigner.\n", "\n", "2- Define an “Operator Pool”. This is simply a collection of operator definitions which will be used to construct the ansatz. For this tutorial, we will use UCCSD.\n", "\n", "3- Initialize qubits to an appropriate reference state. Here, we use HF state to initialize the qubits.\n", "\n", "4- Prepare a trial state with the current ansatz.\n", "\n", "5- Measure the commutator of the Hamiltonian with each operator in the pool to get the gradient. \n", "\n", "6- If the norm of the gradient vector is smaller than some threshold, ε, exit. otherwise, identify the operator with the largest gradient and add this single operator to the left end of the ansatz, with a new variational parameter.\n", "\n", "7- Perform a VQE experiment to re-optimize all parameters in the ansatz.\n", "\n", "8- go to step 4\n", "\n", "Below is a Schematic depiction of the ADAPT-VQE algorithm \n", "\n", "
\n", "\n", "
\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Defaulting to user installation because normal site-packages is not writeable\n", "Requirement already satisfied: pyscf in /home/cudaq/.local/lib/python3.10/site-packages (2.8.0)\n", "Requirement already satisfied: scipy>=1.6.0 in /usr/local/lib/python3.10/dist-packages (from pyscf) (1.12.0)\n", "Requirement already satisfied: h5py>=2.7 in /home/cudaq/.local/lib/python3.10/site-packages (from pyscf) (3.13.0)\n", "Requirement already satisfied: setuptools in /usr/lib/python3/dist-packages (from pyscf) (59.6.0)\n", "Requirement already satisfied: numpy!=1.16,!=1.17,>=1.13 in /usr/local/lib/python3.10/dist-packages (from pyscf) (1.26.4)\n", "Note: you may need to restart the kernel to use updated packages.\n" ] } ], "source": [ "# Requires pyscf to be installed\n", "%pip install pyscf" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import cudaq\n", "\n", "# When use mpi \n", "#cudaq.mpi.initialize()\n", "#print(f\"My rank {cudaq.mpi.rank()} of {cudaq.mpi.num_ranks()}\", flush=True)\n", "\n", "# Set the traget\n", "# Double precision is recommended for the best performance.\n", "cudaq.set_target(\"nvidia\", option = \"fp64\")\n", "\n", "#cudaq.set_target(\"nvidia\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Classical pre-processing\n", "\n", "Here, we compute one and two-electron intgrals using Hartree Fock molecular orbitals." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "converged SCF energy = -1.11675930739643\n", "E(CCSD) = -1.13728399861044 E_corr = -0.02052469121401448\n", "R-CCSD energy = -1.13728399861044\n" ] } ], "source": [ "from pyscf import gto, scf, cc, mcscf, ao2mo\n", "import numpy as np\n", "from functools import reduce\n", "\n", "mol = gto.Mole()\n", "mol.atom = ''' \n", "H 0.0 0.0 0.0\n", "H 0.0 0.0 0.74\n", "'''\n", "mol.basis = 'sto-3g'\n", "mol.charge = 0\n", "mol.spin = 0\n", "mol.build()\n", "mf = scf.RHF(mol)\n", "mf.kernel()\n", "\n", "mycc=cc.CCSD(mf)\n", "mycc.kernel()\n", "\n", "print('R-CCSD energy = ', mycc.e_tot)\n", " \n", "# Compute the electron integrals\n", "\n", "h1e_ao = mol.intor(\"int1e_kin\") + mol.intor(\"int1e_nuc\")\n", "h1e=reduce(np.dot, (mf.mo_coeff.T, h1e_ao, mf.mo_coeff))\n", "h2e_ao = mol.intor(\"int2e_sph\", aosym='1')\n", "h2e=ao2mo.incore.full(h2e_ao, mf.mo_coeff)\n", "h2e=h2e.transpose(0,2,3,1)\n", "nuclear_repulsion = mf.energy_nuc()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Jordan Wigner: \n", "\n", "Convert fermionic Hamiltonian to qubit Hamiltonian." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total number of pauli hamiltonian terms: 19\n", "[0.1659278503377034+0j] IZZI\n", "[0.12062523483390415+0j] IZIZ\n", "[0.17141282644776912+0j] IZII\n", "[0.17441287612261572+0j] IIZZ\n", "[0.04530261550379926+0j] YXXY\n", "[0.04530261550379926+0j] XYYX\n", "[0+0j] XYXY\n", "[0.12062523483390415+0j] ZIZI\n", "[-0.04530261550379926+0j] XXYY\n", "[0+0j] XXXX\n", "[-0.22343153690813433+0j] IIIZ\n", "[0+0j] YXYX\n", "[0.16868898170361205+0j] ZZII\n", "[-0.09706626816763109+0j] IIII\n", "[0.1659278503377034+0j] ZIIZ\n", "[-0.04530261550379926+0j] YYXX\n", "[0.17141282644776912+0j] ZIII\n", "[0+0j] YYYY\n", "[-0.22343153690813433+0j] IIZI\n", "\n" ] } ], "source": [ "from qchem.hamiltonian import jordan_wigner_fermion, generate_molecular_spin_ham_restricted\n", "\n", "# Get the spin molecular hamiltonian\n", "obi, tbi, nuclear_repulsion = generate_molecular_spin_ham_restricted(h1e,h2e, nuclear_repulsion)\n", "\n", "# Convert the fermionic Hamiltonian to a qubit Hamiltonian\n", "tolerance=1e-15\n", "spin_ham = jordan_wigner_fermion(obi, tbi, nuclear_repulsion,tolerance)\n", "print('Total number of pauli hamiltonian terms: ',spin_ham.get_term_count())\n", "\n", "print(spin_ham)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## UCCSD operator pool\n", "\n", "### Single excitation\n", "$$ T_{ij} = \\frac{i}{2} (X_i Y_j - Y_i X_j) \\prod_{p=i+1}^{j-1} Z_p$$\n", "\n", "### Double excitation\n", "$$ T_{ijkl} = \\frac{i}{8} (X_i Y_j X_k X_l + Y_i X_j X_k X_l + Y_i Y_j Y_k X_l + Y_i Y_j X_k Y_l − X_i X_j Y_k X_l − X_i X_j X_k Y_l − Y_i X_j Y_k Y_l − X_i Y_j Y_k Y_l) \\prod_{p=i+1}^{j-1} Zp \\prod_{r=k+1}^{l-1} Z_r $$\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of operator pool: 3\n", "[['XZY', 'YZX'], ['IXZY', 'IYZX'], ['XXYX', 'XXXY', 'XYYY', 'YYXY', 'XYXX', 'YXYY', 'YYYX', 'YXXX']]\n", "[[(-0.5-0j), (0.5+0j)], [(-0.5-0j), (0.5+0j)], [(0.125+0j), (0.125+0j), (0.125+0j), (-0.125-0j), (-0.125-0j), (0.125+0j), (-0.125-0j), (-0.125-0j)]]\n" ] } ], "source": [ "from qchem.operator_pool import get_uccsd_pool\n", "\n", "nelectrons = mol.nelectron\n", "n_qubits= mf.mo_coeff.shape[1] * 2\n", "\n", "pools = get_uccsd_pool(nelectrons, n_qubits)\n", "\n", "print('Number of operator pool: ', len(pools))\n", "\n", "sign_pool = []\n", "mod_pool = []\n", "for i in range(len(pools)):\n", " op_i = pools[i]\n", " temp_op = []\n", " temp_coef = []\n", " \n", " op_i.for_each_term(lambda term: temp_op.append(term.to_string(False)))\n", " op_i.for_each_term(lambda term: temp_coef.append(term.get_coefficient()))\n", " \n", " mod_pool.append(temp_op)\n", " sign_pool.append(temp_coef)\n", "print(mod_pool)\n", "print(sign_pool)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Commutator [$H$, $A_i$]" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of op for gradient: 3\n" ] } ], "source": [ "def commutator(pools, ham):\n", " com_op = []\n", " \n", " for i in range(len(pools)):\n", " # We add the imaginary number that we excluded when generating the operator pool.\n", " op = 1j * pools[i]\n", " \n", " com_op.append(ham * op - op * ham)\n", " \n", " return com_op\n", " \n", "grad_op = commutator(pools, spin_ham)\n", "print('Number of op for gradient: ', len(grad_op))\n", "\n", "#for op in grad_op:\n", "# print(op)\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reference State:\n", "\n", "Reference state here is Haretree Fock" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SV: [(0,0), (0,0), (0,0), (1,0), (0,0), (0,0), (0,0), (0,0), (0,0), (0,0), (0,0), (0,0), (0,0), (0,0), (0,0), (0,0)]\n", "\n" ] } ], "source": [ "# Get the initial state (reference state). \n", "\n", "@cudaq.kernel\n", "def initial_state(n_qubits:int, nelectrons:int):\n", " \n", " qubits = cudaq.qvector(n_qubits)\n", " \n", " for i in range(nelectrons):\n", " x(qubits[i])\n", "\n", "state = cudaq.get_state(initial_state, n_qubits, nelectrons)\n", "print(state)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Quantum kernels:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "###################################\n", "# Quantum kernels\n", "\n", "@cudaq.kernel\n", "def gradient(state:cudaq.State):\n", " q = cudaq.qvector(state)\n", "\n", "\n", "@cudaq.kernel\n", "def kernel(theta: list[float], qubits_num: int, nelectrons: int, pool_single: list[cudaq.pauli_word], \n", " coef_single: list[float], pool_double: list[cudaq.pauli_word], coef_double: list[float]):\n", " q = cudaq.qvector(qubits_num)\n", " \n", " for i in range(nelectrons):\n", " x(q[i])\n", " \n", " count=0\n", " for i in range(0, len(coef_single), 2):\n", " exp_pauli(coef_single[i] * theta[count], q, pool_single[i])\n", " exp_pauli(coef_single[i+1] * theta[count], q, pool_single[i+1])\n", " count+=1\n", "\n", " for i in range(0, len(coef_double), 8):\n", " exp_pauli(coef_double[i] * theta[count], q, pool_double[i])\n", " exp_pauli(coef_double[i+1] * theta[count], q, pool_double[i+1])\n", " exp_pauli(coef_double[i+2] * theta[count], q, pool_double[i+2])\n", " exp_pauli(coef_double[i+3] * theta[count], q, pool_double[i+3])\n", " exp_pauli(coef_double[i+4] * theta[count], q, pool_double[i+4])\n", " exp_pauli(coef_double[i+5] * theta[count], q, pool_double[i+5])\n", " exp_pauli(coef_double[i+6] * theta[count], q, pool_double[i+6])\n", " exp_pauli(coef_double[i+7] * theta[count], q, pool_double[i+7])\n", " count+=1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Beginning of ADAPT-VQE:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Beginning of ADAPT-VQE\n", "Step: 0\n", "Norm of the gradient: 0.36242092403039405\n", "max_grad: 0.36242092403039405\n", "Selected pool at current step: [['XXYX', 'XXXY', 'XYYY', 'YYXY', 'XYXX', 'YXYY', 'YYYX', 'YXXX']]\n", "pool single: []\n", "coef_single: []\n", "pool_double: ['XXYX', 'XXXY', 'XYYY', 'YYXY', 'XYXX', 'YXYY', 'YYYX', 'YXXX']\n", "coef_double: [0.125, 0.125, 0.125, -0.125, -0.125, 0.125, -0.125, -0.125]\n", "tot_single: 0\n", "tot_double: 1\n", "theta_single []\n", "theta_double: [0.0]\n", "theta [0.0]\n", "Optmized Energy: -1.1372838344885003\n", "Optimizer exited successfully: True\n", "CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH\n", "dE: -1.1372838344885003\n", "\n", "\n", "Step: 1\n", "Norm of the gradient: 1.1233558207257577e-07\n", "\n", " Final Result: \n", "\n", "Final parameters: [0.11278279920936399]\n", "Selected pools: [['XXYX', 'XXXY', 'XYYY', 'YYXY', 'XYXX', 'YXYY', 'YYYX', 'YXXX']]\n", "Number of pools: 1\n", "Final energy: -1.1372838344885003\n" ] } ], "source": [ "from scipy.optimize import minimize\n", "\n", "print('Beginning of ADAPT-VQE')\n", "\n", "threshold=1e-3\n", "E_prev=0.0\n", "e_stop=1e-5\n", "init_theta=0.0\n", "\n", "theta_single=[]\n", "theta_double=[]\n", "\n", "pool_single=[]\n", "pool_double=[]\n", "\n", "coef_single=[]\n", "coef_double=[]\n", "\n", "selected_pool=[]\n", "\n", "for i in range(10):\n", " \n", " print('Step: ', i)\n", " \n", " gradient_vec=[]\n", " \n", " for op in grad_op:\n", " grad=cudaq.observe(gradient, op, state).expectation()\n", " gradient_vec.append(grad)\n", " \n", " norm=np.linalg.norm(np.array(gradient_vec))\n", " print('Norm of the gradient: ', norm)\n", " \n", " \n", " # When using mpi to parallelize gradient calculation: uncomment the following lines\n", " \n", " #chunks=np.array_split(np.array(grad_op), cudaq.mpi.num_ranks())\n", " #my_rank_op=chunks[cudaq.mpi.rank()]\n", "\n", " #print('We have', len(grad_op), 'pool operators which we would like to split', flush=True)\n", " #print('We have', len(my_rank_op), 'pool operators on this rank', cudaq.mpi.rank(), flush=True)\n", "\n", " #gradient_vec_async=[]\n", " \n", " #for op in my_rank_op:\n", " #gradient_vec_async.append(cudaq.observe_async(gradient, op, state))\n", "\n", " #gradient_vec_rank=[]\n", " #for i in range(len(gradient_vec_async)):\n", " # get_result=gradient_vec_async[i].get()\n", " # get_expectation=get_result.expectation()\n", " # gradient_vec_rank.append(get_expectation)\n", " \n", " #print('My rank has', len(gradient_vec_rank), 'gradients', flush=True)\n", "\n", " #gradient_vec=cudaq.mpi.all_gather(len(gradient_vec_rank)*cudaq.mpi.num_ranks(), gradient_vec_rank)\n", " \n", " \n", " if norm <= threshold:\n", " print('\\n', 'Final Result: ', '\\n')\n", " print('Final parameters: ', theta)\n", " print('Selected pools: ', selected_pool)\n", " print('Number of pools: ', len(selected_pool))\n", " print('Final energy: ', result_vqe.fun)\n", " \n", " break\n", " \n", " else:\n", " \n", " max_grad=np.max(np.abs(gradient_vec))\n", " print('max_grad: ', max_grad)\n", " \n", " temp_pool = []\n", " temp_sign = []\n", " for i in range(len(mod_pool)):\n", " if np.abs(gradient_vec[i]) == max_grad:\n", " temp_pool.append(mod_pool[i])\n", " temp_sign.append(sign_pool[i])\n", " \n", " print('Selected pool at current step: ', temp_pool)\n", " \n", " selected_pool=selected_pool+temp_pool\n", " \n", " tot_single=0\n", " tot_double=0\n", " for p in temp_pool:\n", " if len(p) == 2:\n", " tot_single += 1\n", " for word in p:\n", " pool_single.append(word)\n", " else:\n", " tot_double += 1\n", " for word in p:\n", " pool_double.append(word)\n", " \n", " for coef in temp_sign:\n", " if len(coef) == 2:\n", " for value in coef:\n", " coef_single.append(value.real)\n", " else:\n", " for value in coef:\n", " coef_double.append(value.real)\n", " \n", " print('pool single: ', pool_single)\n", " print('coef_single: ', coef_single)\n", " print('pool_double: ', pool_double)\n", " print('coef_double: ', coef_double)\n", " print('tot_single: ', tot_single)\n", " print('tot_double: ', tot_double)\n", " \n", " init_theta_single = [init_theta] * tot_single\n", " init_theta_double = [init_theta] * tot_double\n", " \n", " theta_single = theta_single + init_theta_single\n", " theta_double = theta_double + init_theta_double\n", " print('theta_single', theta_single)\n", " print('theta_double: ', theta_double)\n", " \n", " theta = theta_single + theta_double\n", " print('theta', theta)\n", " \n", " def cost(theta):\n", " \n", " theta=theta.tolist()\n", " \n", " energy=cudaq.observe(kernel, spin_ham, theta, n_qubits, nelectrons, pool_single, \n", " coef_single, pool_double, coef_double).expectation()\n", " \n", " return energy\n", " \n", " result_vqe=minimize(cost, theta, method='L-BFGS-B', jac='3-point', tol=1e-7)\n", " \n", " theta=result_vqe.x.tolist()\n", " theta_single = theta[:tot_single]\n", " theta_double = theta[tot_single:]\n", " \n", " print('Optmized Energy: ', result_vqe.fun)\n", " print('Optimizer exited successfully: ',result_vqe.success, flush=True)\n", " print(result_vqe.message, flush=True)\n", " \n", " dE= result_vqe.fun-E_prev\n", " print('dE: ', dE)\n", " print('\\n')\n", " \n", " if np.abs(dE)<=e_stop:\n", " print('\\n', 'Final Result: ', '\\n')\n", " print('Final parameters: ', theta)\n", " print('Selected pools: ', selected_pool)\n", " print('Number of pools: ', len(selected_pool))\n", " print('Final energy: ', result_vqe.fun)\n", " \n", " break\n", " \n", " else:\n", " E_prev=result_vqe.fun\n", " \n", " # Prepare a trial state with the current ansatz.\n", " state=cudaq.get_state(kernel, theta, n_qubits, nelectrons, pool_single, \n", " coef_single, pool_double, coef_double)\n", " \n", "# When using mpi\n", "#cudaq.mpi.finalize()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We obtain the ground state energy of the H2 within chemical accuracy by having only 8 paulit string. This is less than the total number of pauli string (12) of H2. For larger molecules, building a compact quantum circuit can help to reduce the cost and improve perfomance." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.12" } }, "nbformat": 4, "nbformat_minor": 2 }