r"""
Optimization of molecular geometries
====================================
.. meta::
:property="og:description": Find the equilibrium geometry of a molecule
:property="og:image": https://pennylane.ai/qml/_images/fig_pes.png
.. related::
tutorial_quantum_chemistry Building molecular Hamiltonians
tutorial_vqe A brief overview of VQE
tutorial_givens_rotations Givens rotations for quantum chemistry
*Author: Alain Delgado — Posted: 30 June 2021. Last updated: 25 June 2022.*
Predicting the most stable arrangement of atoms in a molecule is one of the most important tasks
in quantum chemistry. Essentially, this is an optimization problem where the total energy of the
molecule is minimized with respect to the positions of the atomic nuclei. The molecular
geometry obtained from this calculation is in fact the starting point for many simulations of
molecular properties. If the geometry is inaccurate, then any calculations that rely on it may
also be inaccurate.
Since the nuclei are much heavier than the electrons, we can treat them as point particles
clamped to their positions. Under this assumption, the total energy of the molecule :math:`E(x)`
depends on the nuclear coordinates :math:`x`, which define the potential energy surface.
Solving the stationary problem :math:`\nabla_x E(x) = 0` corresponds to molecular geometry
optimization and the optimized nuclear coordinates determine the equilibrium geometry of the
molecule. The figure below illustrates these concepts for the
`trihydrogen cation `_. Its equilibrium
geometry in the electronic ground state corresponds to the minimum energy of the potential
energy surface. At this minimum, the three hydrogen atoms are located at the vertices of an
equilateral triangle whose side length is the optimized bond length :math:`d`.

.. figure:: /demonstrations/mol_geo_opt/fig_pes.png
:width: 50%
:align: center

In this tutorial, you will learn how to recast the problem of finding the equilibrium
geometry of a molecule in terms of a general variational quantum algorithm. The
central idea is to consider explicitly that the target electronic Hamiltonian :math:`H(x)`
is a **parametrized** observable that depends on the nuclear coordinates :math:`x`. This
implies that the objective function, defined by the expectation value of the Hamiltonian
computed in the trial state prepared by a quantum computer, depends on both the quantum
circuit and the Hamiltonian parameters.
The quantum algorithm in a nutshell

The goal of the variational algorithm is to find the global minimum of
the cost function :math:`g(\theta, x) = \langle \Psi(\theta) \vert H(x) \vert \Psi(\theta) \rangle`
with respect to the circuit parameters :math:`\theta` and the
nuclear coordinates :math:`x` entering the electronic Hamiltonian of the molecule. To that end,
we use a gradientdescent method and follow a **joint** optimization scheme where the gradients of
the cost function with respect to circuit and Hamiltonian parameters are simultaneously computed
at each step. This approach does not require nested optimization of the state
parameters for each set of nuclear coordinates, as occurs in classical algorithms for
molecular geometry optimization, where the energy minimum is searched for along the potential energy
surface of the electronic state [#jensenbook]_.
In this tutorial we demonstrate how to use PennyLane to implement
quantum optimization of molecular geometries. The algorithm consists of the following steps:
#. Build the parametrized electronic Hamiltonian :math:`H(x)` of the molecule.
#. Design the variational quantum circuit to prepare the electronic trial state of the
molecule, :math:`\vert \Psi(\theta) \rangle`.
#. Define the cost function :math:`g(\theta, x) = \langle \Psi(\theta) \vert H(x) \vert
\Psi(\theta) \rangle`.
#. Initialize the variational parameters :math:`\theta` and :math:`x`. Perform a joint
optimization of the circuit and Hamiltonian parameters to minimize the cost function
:math:`g(\theta, x)`. The gradient with respect to the circuit parameters can be obtained
using a variety of established methods, which are natively supported in PennyLane. The
gradients with respect to the nuclear coordinates can be computed using the formula
.. math::
\nabla_x g(\theta, x) = \langle \Psi(\theta) \vert \nabla_x H(x) \vert \Psi(\theta) \rangle.
Once the optimization is finalized, the circuit parameters determine the energy of the
electronic state, and the nuclear coordinates determine the equilibrium geometry of the
molecule in this state.
Let's get started! ⚛️
Building the parametrized electronic Hamiltonian

In this example, we want to optimize the geometry of the trihydrogen cation
:math:`\mathrm{H}_3^+`, described in a minimal basis set, where two electrons are shared
between three hydrogen atoms (see figure above). The molecule is specified by providing a list
with the symbols of the atomic species and a onedimensional array with the initial
set of nuclear coordinates in `atomic units
`_ .
"""
from pennylane import numpy as np
symbols = ["H", "H", "H"]
x = np.array([0.028, 0.054, 0.0, 0.986, 1.610, 0.0, 1.855, 0.002, 0.0], requires_grad=True)
##############################################################################
# Next, we need to build the parametrized electronic Hamiltonian :math:`H(x)`.
# We use the JordanWigner transformation [#seeley2012]_ to represent the fermionic
# Hamiltonian as a linear combination of Pauli operators,
#
# .. math::
#
# H(x) = \sum_j h_j(x) \prod_i^{N} \sigma_i^{(j)}.
#
# The expansion coefficients :math:`h_j(x)` carry the dependence on the coordinates :math:`x`,
# the operators :math:`\sigma_i` represent the Pauli group :math:`\{I, X, Y, Z\}`, and
# :math:`N` is the number of qubits required to represent the electronic wave function.
#
# We define the function ``H(x)`` to build the parametrized Hamiltonian
# of the trihydrogen cation using the
# :func:`~.pennylane.qchem.molecular_hamiltonian` function.
import pennylane as qml
def H(x):
return qml.qchem.molecular_hamiltonian(symbols, x, charge=1)[0]
##############################################################################
# The variational quantum circuit
# 
#
# Here, we describe the second step of the quantum algorithm: define the quantum circuit
# to prepare the electronic groundstate :math:`\vert \Psi(\theta)\rangle` of the
# :math:`\mathrm{H}_3^+` molecule.
#
# Six qubits are required to encode the occupation number of the molecular spinorbitals.
# To capture the effects of electronic correlations [#kohanoff2006]_, we need to prepare
# the :math:`N`qubit system in a superposition of the HartreeFock state
# :math:`\vert 110000 \rangle` with other states that differ by a double or singleexcitation.
# For example, the state :math:`\vert 000011 \rangle` is obtained by exciting two particles
# from qubits 0, 1 to 4, 5. Similarly, the state :math:`\vert 011000 \rangle` corresponds to a
# single excitation from qubit 0 to 2. This can be done using the singleexcitation and
# doubleexcitation gates :math:`G` and :math:`G^{(2)}` [#qchemcircuits]_ implemented
# in the form of Givens rotations in PennyLane. For more details see the tutorial
# :doc:`tutorial_givens_rotations`.
#
# In addition, we use an adaptive algorithm [#geo_opt_paper]_ to select the excitation
# operations included in the variational quantum circuit. The algorithm proceeds as follows:
#
# #. Generate the indices of the qubits involved in all single and
# doubleexcitations.
# For example, the indices of the singlyexcited state :math:`\vert 011000 \rangle`
# are given by the list ``[0, 2]``. Similarly, the indices of the doublyexcited
# state :math:`\vert 000011 \rangle` are ``[0, 1, 4, 5]``.
#
# #. Construct the circuit using all doubleexcitation gates. Compute the gradient
# of the cost function :math:`g(\theta, x)` with respect to each doubleexcitation
# gate and retain only those with nonzero gradient.
#
# #. Include the selected doubleexcitation gates and repeat the process for the
# singleexcitation gates.
#
# #. Build the final variational quantum circuit by including the selected gates.
#
# For the :math:`\mathrm{H}_3^+` molecule in a minimal basis set we have a total of eight
# excitations of the reference state. After applying the adaptive algorithm the final
# quantum circuit contains only two doubleexcitation operations that act on the qubits
# ``[0, 1, 2, 3]`` and ``[0, 1, 4, 5]``. The circuit is illustrated in the figure below.
#
# 
#
# .. figure:: /demonstrations/mol_geo_opt/fig_circuit.png
# :width: 60%
# :align: center
#
# 
#
# To implement this quantum circuit, we use the
# :func:`~.pennylane.qchem.hf_state` function to generate the
# occupationnumber vector representing the HartreeFock state
hf = qml.qchem.hf_state(electrons=2, orbitals=6)
print(hf)
##############################################################################
# The ``hf`` array is used by the :class:`~.pennylane.BasisState` operation to initialize
# the qubit register. Then, the :class:`~.pennylane.DoubleExcitation` operations are applied
# First, we define the quantum device used to compute the expectation value.
# In this example, we use the ``default.qubit`` simulator:
num_wires = 6
dev = qml.device("default.qubit", wires=num_wires)
@qml.qnode(dev)
def circuit(params, obs, wires):
qml.BasisState(hf, wires=wires)
qml.DoubleExcitation(params[0], wires=[0, 1, 2, 3])
qml.DoubleExcitation(params[1], wires=[0, 1, 4, 5])
return qml.expval(obs)
##############################################################################
# This circuit prepares the trial state
#
# .. math::
#
# \vert\Psi(\theta_1, \theta_2)\rangle =
# \mathrm{cos}(\theta_1)\mathrm{cos}(\theta_2)\vert110000\rangle 
# \mathrm{cos}(\theta_1)\mathrm{sin}(\theta_2)\vert000011\rangle 
# \mathrm{sin}(\theta_1)\vert001100\rangle,
#
# where :math:`\theta_1` and :math:`\theta_2` are the circuit parameters that need to be
# optimized to find the groundstate energy of the molecule.
#
# The cost function and the nuclear gradients
# 
#
# The third step of the algorithm is to define the cost function
# :math:`g(\theta, x) = \langle \Psi(\theta) \vert H(x) \vert\Psi(\theta) \rangle`. It
# evaluates the expectation value of the parametrized Hamiltonian :math:`H(x)` in the
# trial state :math:`\vert\Psi(\theta)\rangle`.
##############################################################################
# Next, we define the ``cost`` function :math:`g(\theta, x)` which depends on
# both the circuit and the Hamiltonian parameters. Specifically we consider the
# expectation values of the Hamiltonian.
def cost(params, x):
hamiltonian = H(x)
return circuit(params, obs=hamiltonian, wires=range(num_wires))
##############################################################################
#
# We minimize the cost function :math:`g(\theta, x)` using a gradientbased
# method, and compute the gradients with respect to both the
# circuit parameters :math:`\theta` and the nuclear coordinates :math:`x`.
# The circuit gradients are computed analytically using the automatic differentiation
# techniques available in PennyLane. The nuclear gradients are evaluated
# by taking the expectation value of the gradient of the electronic Hamiltonian,
#
# .. math::
#
# \nabla_x g(\theta, x) = \langle \Psi(\theta) \vert \nabla_x H(x) \vert \Psi(\theta) \rangle.
#
# We use the :func:`finite_diff` function to compute the gradient of
# the Hamiltonian using a centraldifference approximation. Then, we evaluate the expectation
# value of the gradient components :math:`\frac{\partial H(x)}{\partial x_i}`. This is implemented by
# the function ``grad_x``:
def finite_diff(f, x, delta=0.01):
"""Compute the centraldifference finite difference of a function"""
gradient = []
for i in range(len(x)):
shift = np.zeros_like(x)
shift[i] += 0.5 * delta
res = (f(x + shift)  f(x  shift)) * delta**1
gradient.append(res)
return gradient
def grad_x(params, x):
grad_h = finite_diff(H, x)
grad = [circuit(params, obs=obs, wires=range(num_wires)) for obs in grad_h]
return np.array(grad)
##############################################################################
# Optimization of the molecular geometry
# 
#
# Finally, we proceed to minimize our cost function to find the ground state equilibrium
# geometry of the :math:`\mathrm{H}_3^+` molecule. As a reminder,
# the circuit parameters and the nuclear coordinates will be jointly optimized at
# each optimization step. This approach does not require nested VQE
# optimization of the circuit parameters for each set of nuclear coordinates.
#
# We start by defining the classical optimizers:
opt_theta = qml.GradientDescentOptimizer(stepsize=0.4)
opt_x = qml.GradientDescentOptimizer(stepsize=0.8)
##############################################################################
# Next, we initialize the circuit parameters :math:`\theta`. The angles
# :math:`\theta_1` and :math:`\theta_2` are set to zero so that the
# initial state :math:`\vert\Psi(\theta_1, \theta_2)\rangle`
# is the HartreeFock state.
theta = np.array([0.0, 0.0], requires_grad=True)
##############################################################################
# The initial set of nuclear coordinates :math:`x`, defined at
# the beginning of the tutorial, was computed classically within the HartreeFock
# approximation using the GAMESS program [#ref_gamess]_. This is a natural choice
# for the starting geometry that we are aiming to improve due to the electronic
# correlation effects included in the trial state :math:`\vert\Psi(\theta)\rangle`.
#
# We carry out the optimization over a maximum of 100 steps.
# The circuit parameters and the nuclear coordinates are optimized until the
# maximum component of the nuclear gradient :math:`\nabla_x g(\theta,x)` is
# less than or equal to :math:`10^{5}` Hartree/Bohr. Typically, this is the
# convergence criterion used for optimizing molecular geometries in
# quantum chemistry simulations.
from functools import partial
# store the values of the cost function
energy = []
# store the values of the bond length
bond_length = []
# Factor to convert from Bohrs to Angstroms
bohr_angs = 0.529177210903
for n in range(100):
# Optimize the circuit parameters
theta.requires_grad = True
x.requires_grad = False
theta, _ = opt_theta.step(cost, theta, x)
# Optimize the nuclear coordinates
x.requires_grad = True
theta.requires_grad = False
_, x = opt_x.step(cost, theta, x, grad_fn=grad_x)
energy.append(cost(theta, x))
bond_length.append(np.linalg.norm(x[0:3]  x[3:6]) * bohr_angs)
if n % 4 == 0:
print(f"Step = {n}, E = {energy[1]:.8f} Ha, bond length = {bond_length[1]:.5f} A")
# Check maximum component of the nuclear gradient
if np.max(grad_x(theta, x)) <= 1e05:
break
print("\n" f"Final value of the groundstate energy = {energy[1]:.8f} Ha")
print("\n" "Groundstate equilibrium geometry")
print("%s %4s %8s %8s" % ("symbol", "x", "y", "z"))
for i, atom in enumerate(symbols):
print(f" {atom} {x[3 * i]:.4f} {x[3 * i + 1]:.4f} {x[3 * i + 2]:.4f}")
##############################################################################
# Next, we plot the values of the ground state energy of the molecule
# and the bond length as a function of the optimization step.
import matplotlib.pyplot as plt
fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(12)
# Add energy plot on column 1
E_fci = 1.27443765658
E_vqe = np.array(energy)
ax1 = fig.add_subplot(121)
ax1.plot(range(n + 1), E_vqe  E_fci, "go", ls="dashed")
ax1.plot(range(n + 1), np.full(n + 1, 0.001), color="red")
ax1.set_xlabel("Optimization step", fontsize=13)
ax1.set_ylabel("$E_{VQE}  E_{FCI}$ (Hartree)", fontsize=13)
ax1.text(5, 0.0013, r"Chemical accuracy", fontsize=13)
plt.yscale("log")
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
# Add bond length plot on column 2
d_fci = 0.986
ax2 = fig.add_subplot(122)
ax2.plot(range(n + 1), bond_length, "go", ls="dashed")
ax2.plot(range(n + 1), np.full(n + 1, d_fci), color="red")
ax2.set_ylim([0.965, 0.99])
ax2.set_xlabel("Optimization step", fontsize=13)
ax2.set_ylabel("bond length ($\AA$)", fontsize=13)
ax2.text(5, 0.9865, r"Equilibrium bond length", fontsize=13)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.subplots_adjust(wspace=0.3)
plt.show()
##############################################################################
# 
# Notice that despite the fact that the groundstate energy is already converged
# within chemical accuracy (:math:`0.0016` Ha) after the fourth iteration, more
# optimization steps are required to find the equilibrium bond length of the
# molecule.
#
# The figure below animates snapshots of the atomic structure of the
# trihydrogen cation as the quantum algorithm was searching for the equilibrium
# geometry. For visualization purposes, the initial nuclear coordinates were
# generated by perturbing the HF geometry. The quantum algorithm
# is able to find the correct equilibrium geometry of the :math:`\mathrm{H}_3^+`
# molecule where the three H atoms are located at the vertices of an equilateral triangle.
#
# 
#
# .. figure:: /demonstrations/mol_geo_opt/fig_movie.gif
# :width: 50%
# :align: center
#
# 
#
# To summarize, we have shown how the scope of variational quantum algorithms can be
# extended to perform quantum simulations of molecules involving both the electronic and
# the nuclear degrees of freedom. The joint optimization scheme described here
# is a generalization of the usual VQE algorithm where only the electronic
# state is parametrized. Extending the applicability of the variational quantum algorithms to
# target parametrized Hamiltonians could be also relevant to simulate the optical properties of
# molecules where the fermionic observables depend also on the electric field of the
# incoming radiation [#pulay]_.
#
# References
# 
#
# .. [#jensenbook]
#
# F. Jensen. "Introduction to computational chemistry".
# (John Wiley & Sons, 2016).
#
# .. [#seeley2012]
#
# Jacob T. Seeley, Martin J. Richard, Peter J. Love. "The BravyiKitaev transformation for
# quantum computation of electronic structure". `Journal of Chemical Physics 137, 224109
# (2012).
# `__
#
# .. [#kohanoff2006]
#
# Jorge Kohanoff. "Electronic structure calculations for solids and molecules: theory and
# computational methods". (Cambridge University Press, 2006).
#
# .. [#qchemcircuits]
#
# J.M. Arrazola, O. Di Matteo, N. Quesada, S. Jahangiri, A. Delgado, N. Killoran.
# "Universal quantum circuits for quantum chemistry". arXiv:2106.13839, (2021)
#
# .. [#ref_gamess]
#
# M.W. Schmidt, K.K. Baldridge, J.A. Boatz, S.T. Elbert, M.S. Gordon, J.H. Jensen,
# S. Koseki, N. Matsunaga, K.A. Nguyen, S.Su, *et al.* "General atomic and molecular
# electronic structure system". `Journal of Computational Chemistry 14, 1347 (1993)
# `__
#
# .. [#geo_opt_paper]
#
# A. Delgado, J.M. Arrazola, S. Jahangiri, Z. Niu, J. Izaac, C. Roberts, N. Killoran.
# "Variational quantum algorithm for molecular geometry optimization".
# arXiv:2106.13840, (2021)
#
# .. [#pulay]
#
# P. Pulay.
# "Analytical derivative methods in quantum chemistry".
# `Advances in Chemical Sciences (1987)
# `__
#
# About the author
# 
# .. include:: ../_static/authors/alain_delgado.txt