{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# This cell is added by sphinx-gallery\n# It can be customized to whatever you like\n%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Differentiable Hartree-Fock\n===========================\n\n::: {.meta}\n:property=\\\"og:description\\\": Learn how to use the differentiable\nHartree-Fock solver :property=\\\"og:image\\\":\n\n:::\n\n::: {.related}\ntutorial\\_quantum\\_chemistry Building molecular Hamiltonians\ntutorial\\_vqe A brief overview of VQE tutorial\\_givens\\_rotations Givens\nrotations for quantum chemistry tutorial\\_adaptive\\_circuits Adaptive\ncircuits for quantum chemistry\n:::\n\n*Author: Soran Jahangiri --- Posted: 09 May 2022. Last updated: 09 May\n2022.*\n\nIn this tutorial, you will learn how to use PennyLane\\'s differentiable\nHartree-Fock solver . The quantum chemistry module in PennyLane,\n`qml.qchem `{.interpreted-text role=\"mod\"}, provides\nbuilt-in methods for constructing atomic and molecular orbitals,\nbuilding Fock matrices and solving the self-consistent field equations\nto obtain optimized orbitals which can be used to construct\nfully-differentiable molecular Hamiltonians. PennyLane allows users to\nnatively compute derivatives of all these objects with respect to the\nunderlying parameters using the methods of [automatic\ndifferentiation](https://pennylane.ai/qml/glossary/automatic_differentiation.html).\nWe introduce a workflow to jointly optimize circuit parameters, nuclear\ncoordinates and basis set parameters in a variational quantum\neigensolver algorithm. You will also learn how to visualize the atomic\nand molecular orbitals which can be used to create an animation like\nthis:\n\n![The bonding molecular orbital of hydrogen visualized during a full\ngeometry, circuit and basis set\noptimization.](/demonstrations/differentiable_HF/h2.gif){.align-center\nwidth=\"60.0%\"}\n\nLet\\'s get started!\n\nDifferentiable Hamiltonians\n---------------------------\n\nVariational quantum algorithms aim to calculate the energy of a molecule\nby constructing a parameterized quantum circuit and finding a set of\nparameters that minimize the expectation value of the electronic\n[molecular\nHamiltonian](https://en.wikipedia.org/wiki/Molecular_Hamiltonian). The\noptimization can be carried out by computing the gradients of the\nexpectation value with respect to these parameters and iteratively\nupdating them until convergence is achieved. In principle, the\noptimization process is not limited to the circuit parameters and can be\nextended to include the parameters of the Hamiltonian that can be\noptimized concurrently with the circuit parameters. The aim is now to\nobtain the set of parameters that minimize the following expectation\nvalue\n\n$$\\left \\langle \\Psi(\\theta) | H(\\beta) | \\Psi(\\theta) \\right \\rangle,$$\n\nwhere $\\theta$ and $\\beta$ represent the circuit and Hamiltonian\nparameters, respectively.\n\nComputing the gradient of a molecular Hamiltonian is challenging because\nthe dependency of the Hamiltonian on the molecular parameters is\ntypically not very straightforward. This makes symbolic differentiation\nmethods, which obtain derivatives of an input function by direct\nmathematical manipulation, of limited scope. Furthermore, numerical\ndifferentiation methods based on [finite\ndifferences](https://en.wikipedia.org/wiki/Finite_difference_method) are\nnot always reliable due to their intrinsic instability, especially when\nthe number of differentiable parameters is large. These limitations can\nbe alleviated by using automatic differentiation methods which can be\nused to compute exact gradients of a function, implemented with computer\ncode, using resources comparable to those required to evaluate the\nfunction itself.\n\nEfficient optimization of the molecular Hamiltonian parameters in a\nvariational quantum algorithm is essential for tackling problems such as\n[geometry\noptimization](https://pennylane.ai/qml/demos/tutorial_mol_geo_opt.html)\nand vibrational frequency calculations. These problems require computing\nthe first- and second-order derivatives of the molecular energy with\nrespect to nuclear coordinates which can be efficiently obtained if the\nvariational workflow is automatically differentiable. Another important\nexample is the simultaneous optimization of the parameters of the basis\nset used to construct the atomic orbitals which can in principle\nincrease the accuracy of the computed energy without increasing the\nnumber of qubits in a quantum simulation. The joint optimization of the\ncircuit and Hamiltonian parameters can also be used when the chemical\nproblem involves optimizing the parameters of external potentials.\n\nThe Hartree-Fock method\n-----------------------\n\nThe main goal of the Hartree-Fock method is to obtain molecular orbitals\nthat minimize the energy of a system where electrons are treated as\nindependent particles that experience a mean field generated by the\nother electrons. These optimized molecular orbitals are then used to\nconstruct one- and two-body electron integrals in the basis of molecular\norbitals, $\\phi$,\n\n$$\\begin{aligned}\nh_{pq} =\\int dx \\,\\phi_p^*(x)\\left(-\\frac{\\nabla^2}{2}-\\sum_{i=1}^N\\frac{Z_i}{|r-R_i|}\\right)\\phi_q(x),\\\\\\\\\n\\end{aligned}$$\n\n$$h_{pqrs} = \\int dx_1 dx_2\\, \\frac{\\phi_p^*(x_1)\\phi_q^*(x_2)\\phi_r(x_2)\\phi_s(x_1)}{|r_1-r_2|}.$$\n\nThese integrals are used to generate a differentiable\n[second-quantized](https://en.wikipedia.org/wiki/Second_quantization)\nmolecular Hamiltonian as\n\n$$H=\\sum_{pq} h_{pq}a_p^\\dagger a_q +\\frac{1}{2}\\sum_{pqrs}h_{pqrs}a_p^\\dagger a_q^\\dagger a_r a_s,$$\n\nwhere $a^\\dagger$ and $a$ are the fermionic creation and annihilation\noperators, respectively. This Hamiltonian is then transformed to the\nqubit basis. Let\\'s see how this can be done in PennyLane.\n\nTo get started, we need to define the atomic symbols and the nuclear\ncoordinates of the molecule. For the hydrogen molecule we have\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from autograd import grad\nimport pennylane as qml\nfrom pennylane import numpy as np\nimport matplotlib.pyplot as plt\nnp.set_printoptions(precision=5)\n\nsymbols = [\"H\", \"H\"]\n# optimized geometry at the Hartree-Fock level\ngeometry = np.array([[-0.672943567415407, 0.0, 0.0],\n [ 0.672943567415407, 0.0, 0.0]], requires_grad=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The use of `requires_grad=True` specifies that the nuclear coordinates\nare differentiable parameters. We can now compute the Hartree-Fock\nenergy and its gradient with respect to the nuclear coordinates. To do\nthat, we create a molecule object that stores all the molecular\nparameters needed to perform a Hartree-Fock calculation.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"mol = qml.qchem.Molecule(symbols, geometry)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Hartree-Fock energy can now be computed with the\n`~.pennylane.qchem.hf_energy`{.interpreted-text role=\"func\"} function\nwhich is a function transform\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"qml.qchem.hf_energy(mol)(geometry)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We now compute the gradient of the energy with respect to the nuclear\ncoordinates\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"grad(qml.qchem.hf_energy(mol))(geometry)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The obtained gradients are equal or very close to zero because the\ngeometry we used here has been previously optimized at the Hartree-Fock\nlevel. You can use a different geometry and verify that the newly\ncomputed gradients are not all zero.\n\nWe can also compute the values and gradients of several other quantities\nthat are obtained during the Hartree-Fock procedure. These include all\nintegrals over basis functions, matrices formed from these integrals and\nthe one- and two-body integrals over molecular orbitals. Let\\'s look at\na few examples.\n\nWe first compute the overlap integral between the two S-type atomic\norbitals of the hydrogen atoms. Here we are using the\n[STO-3G](https://en.wikipedia.org/wiki/STO-nG_basis_sets) basis set in\nwhich each of the atomic orbitals is represented by one basis function\ncomposed of three primitive Gaussian functions. These basis functions\ncan be accessed from the molecule object as\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"S1 = mol.basis_set[0]\nS2 = mol.basis_set[1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can check the parameters of the basis functions as\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"for param in S1.params:\n print(param)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This returns the exponents, contraction coefficients and the centres of\nthe three Gaussian functions of the STO-3G basis set. These parameters\ncan be also obtained individually by using `S1.alpha`, `S1.coeff` and\n`S1.r`, respectively. You can verify that both of the `S1` an `S2`\natomic orbitals have the same exponents and contraction coefficients but\nare centred on different hydrogen atoms. You can also verify that the\norbitals are S-type by printing the angular momentum quantum numbers\nwith\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"S1.l"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This gives us a tuple of three integers, representing the exponents of\nthe $x$, $y$ and $z$ components in the Gaussian functions.\n\nWe can now compute the overlap integral,\n\n$$S_{\\mu \\nu} = \\int \\chi_\\mu^* (r) \\chi_\\nu (r) dr$$\n\nbetween the atomic orbitals $\\chi$, by passing the orbitals and the\ninitial values of their centres to the\n`~.pennylane.qchem.overlap_integral`{.interpreted-text role=\"func\"}\nfunction. The centres of the orbitals are those of the hydrogen atoms by\ndefault and are therefore treated as differentiable parameters by\nPennyLane.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"qml.qchem.overlap_integral(S1, S2)([geometry[0], geometry[1]])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can verify that the overlap integral between two identical atomic\norbitals is equal to one. We can now compute the gradient of the overlap\nintegral with respect to the orbital centres\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"grad(qml.qchem.overlap_integral(S1, S2))([geometry[0], geometry[1]])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Can you explain why some of the computed gradients are zero?\n\nLet\\'s now plot the atomic orbitals and their overlap. We can do it by\nusing the\n:py`~.pennylane.qchem.Molecule.atomic_orbital`{.interpreted-text\nrole=\"meth\"} function, which evaluates the atomic orbital at a given\ncoordinate. For instance, the value of the S orbital on the first\nhydrogen atom can be computed at the origin as\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"V1 = mol.atomic_orbital(0)\nV1(0.0, 0.0, 0.0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can evaluate this orbital at different points along the $x$ axis and\nplot it.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x = np.linspace(-5, 5, 1000)\nplt.plot(x, V1(x, 0.0, 0.0), color='teal')\nplt.xlabel('X [Bohr]')\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can also plot the second S orbital and visualize the overlap between\nthem\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"V2 = mol.atomic_orbital(1)\nplt.plot(x, V1(x, 0.0, 0.0), color='teal')\nplt.plot(x, V2(x, 0.0, 0.0), color='teal')\nplt.fill_between(\n x, np.minimum(V1(x, 0.0, 0.0), V2(x, 0.0, 0.0)), color = 'red', alpha = 0.5, hatch = '||')\nplt.xlabel('X [Bohr]')\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"By looking at the orbitals, can you guess at what distance the value of\nthe overlap becomes negligible? Can you verify your guess by computing\nthe overlap at that distance?\n\nSimilarly, we can plot the molecular orbitals of the hydrogen molecule\nobtained from the Hartree-Fock calculations. We plot the cross-section\nof the bonding orbital on the $x-y$ plane.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"n = 30 # number of grid points along each axis\n\nmol.mo_coefficients = mol.mo_coefficients.T\nmo = mol.molecular_orbital(0)\nx, y = np.meshgrid(np.linspace(-2, 2, n),\n np.linspace(-2, 2, n))\nval = np.vectorize(mo)(x, y, 0)\nval = np.array([val[i][j]._value for i in range(n) for j in range(n)]).reshape(n, n)\n\nfig, ax = plt.subplots()\nco = ax.contour(x, y, val, 10, cmap='summer_r', zorder=0)\nax.clabel(co, inline=2, fontsize=10)\nplt.scatter(mol.coordinates[:,0], mol.coordinates[:,1], s = 80, color='black')\nax.set_xlabel('X [Bohr]')\nax.set_ylabel('Y [Bohr]')\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"VQE simulations\n===============\n\nBy performing the Hartree-Fock calculations, we obtain a set of one- and\ntwo-body integrals over molecular orbitals that can be used to construct\nthe molecular Hamiltonian with the\n`~.pennylane.qchem.molecular_hamiltonian`{.interpreted-text role=\"func\"}\nfunction.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"hamiltonian, qubits = qml.qchem.molecular_hamiltonian(mol.symbols,\n mol.coordinates,\n args=[mol.coordinates])\nprint(hamiltonian)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Hamiltonian contains 15 terms and, importantly, the coefficients of\nthe Hamiltonian are all differentiable. We can construct a circuit and\nperform a VQE simulation in which both of the circuit parameters and the\nnuclear coordinates are optimized simultaneously by using the computed\ngradients. We will have two sets of differentiable parameters: the first\nset is the rotation angles of the excitation gates which are applied to\nthe reference Hartree-Fock state to construct the exact ground state.\nThe second set contains the nuclear coordinates of the hydrogen atoms.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"dev = qml.device(\"default.qubit\", wires=4)\ndef energy(mol):\n @qml.qnode(dev)\n def circuit(*args):\n qml.BasisState(np.array([1, 1, 0, 0]), wires=range(4))\n qml.DoubleExcitation(*args[0][0], wires=[0, 1, 2, 3])\n H = qml.qchem.molecular_hamiltonian(mol.symbols, mol.coordinates, alpha=mol.alpha,\n coeff=mol.coeff, args=args[1:])[0]\n return qml.expval(H)\n return circuit"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that we only use the\n`~.pennylane.DoubleExcitation`{.interpreted-text role=\"class\"} gate as\nthe `~.pennylane.SingleExcitation`{.interpreted-text role=\"class\"} ones\ncan be neglected in this particular example . We now compute the\ngradients of the energy with respect to the circuit parameter and the\nnuclear coordinates and update the parameters iteratively. Note that the\nnuclear coordinate gradients are simply the forces on the atomic nuclei.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# initial value of the circuit parameter\ncircuit_param = [np.array([0.0], requires_grad=True)]\n\nfor n in range(36):\n\n args = [circuit_param, geometry]\n mol = qml.qchem.Molecule(symbols, geometry)\n\n # gradient for circuit parameters\n g_param = grad(energy(mol), argnum = 0)(*args)\n circuit_param = circuit_param - 0.25 * g_param[0]\n\n # gradient for nuclear coordinates\n forces = -grad(energy(mol), argnum = 1)(*args)\n geometry = geometry + 0.5 * forces\n\n if n % 5 == 0:\n print(f'n: {n}, E: {energy(mol)(*args):.8f}, Force-max: {abs(forces).max():.8f}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"After 35 steps of optimization, the forces on the atomic nuclei and the\ngradient of the circuit parameter are both approaching zero, and the\nenergy of the molecule is that of the optimized geometry at the\n[full-CI](https://en.wikipedia.org/wiki/Full_configuration_interaction)\nlevel: $-1.1373060483$ Ha. You can print the optimized geometry and\nverify that the final bond length of hydrogen is identical to the one\ncomputed with full-CI which is $1.389$\n[Bohr](https://en.wikipedia.org/wiki/Bohr_radius).\n\nWe are now ready to perform a full optimization where the circuit\nparameters, the nuclear coordinates and the basis set parameters are all\ndifferentiable parameters that can be optimized simultaneously.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# initial values of the nuclear coordinates\ngeometry = np.array([[0.0, 0.0, -0.672943567415407],\n [0.0, 0.0, 0.672943567415407]], requires_grad=True)\n\n# initial values of the basis set exponents\nalpha = np.array([[3.42525091, 0.62391373, 0.1688554],\n [3.42525091, 0.62391373, 0.1688554]], requires_grad=True)\n\n# initial values of the basis set contraction coefficients\ncoeff = np.array([[0.1543289673, 0.5353281423, 0.4446345422],\n [0.1543289673, 0.5353281423, 0.4446345422]], requires_grad=True)\n\n# initial value of the circuit parameter\ncircuit_param = [np.array([0.0], requires_grad=True)]\n\nfor n in range(36):\n\n args = [circuit_param, geometry, alpha, coeff]\n mol = qml.qchem.Molecule(symbols, geometry, alpha=alpha, coeff=coeff)\n\n # gradient for circuit parameters\n g_param = grad(energy(mol), argnum=0)(*args)\n circuit_param = circuit_param - 0.25 * g_param[0]\n\n # gradient for nuclear coordinates\n forces = -grad(energy(mol), argnum=1)(*args)\n geometry = geometry + 0.5 * forces\n\n # gradient for basis set exponents\n g_alpha = grad(energy(mol), argnum=2)(*args)\n alpha = alpha - 0.25 * g_alpha\n\n # gradient for basis set contraction coefficients\n g_coeff = grad(energy(mol), argnum=3)(*args)\n coeff = coeff - 0.25 * g_coeff\n\n if n % 5 == 0:\n print(f'n: {n}, E: {energy(mol)(*args):.8f}, Force-max: {abs(forces).max():.8f}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can also print the gradients of the circuit and basis set parameters\nand confirm that they are approaching zero. The computed energy of\n$-1.14040160$ Ha is lower than the full-CI energy, $-1.1373060483$ Ha\n(obtained with the STO-3G basis set for the hydrogen molecule) because\nwe have optimized the basis set parameters in our example. This means\nthat we can reach a lower energy for hydrogen without increasing the\nbasis set size, which would otherwise lead to a larger number of qubits.\n\nConclusions\n===========\n\nThis tutorial introduces an important feature of PennyLane that allows\nperforming fully-differentiable Hartree-Fock and subsequently VQE\nsimulations. This feature provides two major benefits: i) All gradient\ncomputations needed for parameter optimization can be carried out with\nthe elegant methods of automatic differentiation which facilitates\nsimultaneous optimizations of circuit and Hamiltonian parameters in\napplications such as VQE molecular geometry optimizations. ii) By\noptimizing the molecular parameters such as the exponent and contraction\ncoefficients of Gaussian functions of the basis set, one can reach a\nlower energy without increasing the number of basis functions. Can you\nthink of other interesting molecular parameters that can be optimized\nalong with the nuclear coordinates and basis set parameters that we\noptimized in this tutorial?\n\nReferences\n==========\n\nAbout the author\n================\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.15"
}
},
"nbformat": 4,
"nbformat_minor": 0
}