Hi,
I have studied the ADAPT-VQE demo code here Adaptive circuits for quantum chemistry | PennyLane Demos and then put together a script to do a PES scan for LiH dissociation based on the demo here Modelling chemical reactions on a quantum computer | PennyLane Demos
I observe that the ADATP-VQE and VQE energies are in close agreement in the vicinity of the equilibrium bond distance (and shorter distances). However at longer separation I observe a steady drift between ADAPT-VQE and VQE.
This is more apparent when plotting the differences in energy w.r.t to exact FCI.
The increased VQE error may be due to the poorer performance of the optimiser that hits the maximum number of iterations (50) beyond distances of ca. 4 Bohr.
However the ADAPT-VQE error increases significantly more rapidly. I have repeated runs after adjusting the gradient tolerance which had no effect (most excitations have a zero gradient contribution). Changing the adaptive optimiser settings param_steps and step_size has a small effect on the accuracy of the ADAPT-VQE energies, but I am unable to reduce significantly further the gap in accuracy between VQE and ADAPT-VQE.
I am wondering if:
-
I made a mistake implementing ADAPT-VQE ?
-
Whether there is an intrinsic reason why it would be more difficult to minimise the ground state energies with ADAPT-VQE ? My understanding of the literature is that both methods should be able to give the same energies.
I have pasted below the complete code to reproduce these findings.
Hello! If applicable, put your complete code example down below. Make sure that your code:
import pennylane as qml
from pennylane import qchem
from jax import numpy as jnp
import numpy as np
# system setup
electrons = 4
orbitals = 12
active_electrons = 2
active_orbitals = 5
# atomic symbols defining the molecule
symbols = ['Li', 'H']
# set up a loop to change bond length
r_range = jnp.arange(1.25, 6.50, 0.50)
Full VQE
vqe_energies = []
# list to store number of VQE steps to reach convergence
vqe_steps = []
# keeps track of points in the potential energy surface
pes_point = 0
for r in r_range:
# Change only the z coordinate of one atom
coordinates = jnp.array([[0.0, 0.0, 0.0], [0.0, 0.0, r]])
# Construct the Molecule object
molecule = qchem.Molecule(symbols, coordinates)
# Obtain the qubit Hamiltonian
#H, qubits = qchem.molecular_hamiltonian(molecule, active_electrons=active_electrons, active_orbitals=active_orbitals, method='dhf')
H, qubits = qchem.molecular_hamiltonian(molecule, active_electrons=active_electrons, active_orbitals=active_orbitals, method='pyscf')
H_sparse = H.sparse_matrix()
singles, doubles = qchem.excitations(active_electrons, qubits)
hf = qchem.hf_state(active_electrons, qubits)
print(f"Total number of excitations = {len(singles) + len(doubles)}")
#sys.exit(-1)
# VQE define the device, optimizer and circuit
dev = qml.device("lightning.qubit", wires=qubits)
opt = optax.sgd(learning_rate=0.4) # sgd stands for StochasticGradientDescent
@qml.qnode(dev, interface='jax')
def circuit(parameters):
AllSinglesDoubles(parameters, range(qubits), hf, singles, doubles)
return qml.expval(H) # we are interested in minimizing this expectation value
# return qml.expval(qml.SparseHamiltonian(H_sparse, wires=range(qubits))) # Not supported with lightning.qubit
# initialize the gate parameters
# init_params = jnp.zeros(3)
init_params = jnp.zeros(len(singles) + len(doubles))
# initialize with converged parameters from previous point
if pes_point > 0:
init_params = params_old
prev_energy = 0.0
@qml.qjit
def update_step(i, params, opt_state):
"""Perform a single gradient update step"""
grads = qml.grad(circuit)(params)
updates, opt_state = opt.update(grads, opt_state)
params = optax.apply_updates(params, updates)
return (params, opt_state)
loss_history = []
opt_state = opt.init(init_params)
params = init_params
for i in range(50):
params, opt_state = update_step(i, params, opt_state)
energy = circuit(params)
if jnp.abs(energy - prev_energy) < 1e-6:
break
prev_energy = energy
vqe_steps.append(i)
# exact result
print ("Computed distance %s, VQE energy %s, Hamiltonian length %s, number of qubits %s, number of VQE steps %s " % (r,energy, len(H),qubits, i))
# store the converged parameters
params_old = params
pes_point = pes_point + 1
vqe_energies.append(energy)
Exact diagonalisation and HF energy
from pennylane import numpy as np
exact_energies = []
hf_energies = []
for r in r_range:
# Change only the z coordinate of one atom
coordinates = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, r]], requires_grad= False)
# Construct the Molecule object
molecule = qchem.Molecule(symbols, coordinates)
H, qubits = qchem.molecular_hamiltonian(molecule, active_electrons=active_electrons, active_orbitals=active_orbitals, method='pyscf')
H_sparse = qml.SparseHamiltonian(H.sparse_matrix(), wires=H.wires)
eigenvals = qml.eigvals(H_sparse)
args = []
hf_energy = qchem.hf_energy(molecule)(*args)
print ("Computed distance %s, exact energy %s HF energy %s" % (r,eigenvals[0], hf_energy))
exact_energies.append(eigenvals[0])
hf_energies.append(hf_energy)
Adapt VQE using Sparse Matrix representation
sparse_adapt_vqe_energies = []
for r in r_range:
# Change only the z coordinate of one atom
coordinates = jnp.array([[0.0, 0.0, 0.0], [0.0, 0.0, r]])
# Construct the Molecule object
molecule = qchem.Molecule(symbols, coordinates)
# Obtain the qubit Hamiltonian
H, qubits = qchem.molecular_hamiltonian(molecule, active_electrons=active_electrons, active_orbitals=active_orbitals, method='pyscf')
H_sparse = H.sparse_matrix()
singles, doubles = qchem.excitations(active_electrons, qubits)
hf_state = qchem.hf_state(active_electrons, qubits)
#print(f"Total number of excitations = {len(singles) + len(doubles)}")
#sys.exit(-1)
singles_excitations = [qml.SingleExcitation(0.0, x) for x in singles]
doubles_excitations = [qml.DoubleExcitation(0.0, x) for x in doubles]
operator_pool = doubles_excitations + singles_excitations
# VQE define the device, optimizer and circuit
dev = qml.device("lightning.qubit", wires=qubits)
# ADAPT VQE
@qml.qnode(dev, diff_method="parameter-shift")
#@qml.qnode(dev)
def circuit():
[qml.PauliX(i) for i in np.nonzero(hf_state)[0]]
return qml.expval(qml.SparseHamiltonian(H_sparse, wires=range(qubits)))
#return qml.expval(H)
opt = qml.optimize.AdaptiveOptimizer(param_steps=75, stepsize=0.6)#empiricall checks
for i in range(len(operator_pool)):
circuit, energy, gradient = opt.step_and_cost(circuit, operator_pool, drain_pool=True)
if i % 2 == 0:
print("n = {:}, E = {:.8f} H, Largest Gradient = {:.3f}".format(i, energy, gradient))
# print(qml.draw(circuit, decimals=None)())
# print()
if gradient < 1e-5:
break
print ("Computed distance %s, VQE energy %s, Hamiltonian length %s, number of qubits %s" % (r,energy, len(H),qubits))
sparse_adapt_vqe_energies.append(energy)
Plotting
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.plot(r_range, sparse_adapt_vqe_energies, label='adapt-vqe')
ax.plot(r_range, vqe_energies, label='vqe')
ax.plot(r_range, hf_energies, label='hf')
ax.plot(r_range, exact_energies, label='fci')
plt.legend(loc='upper right')
ax.set(
xlabel="Bond length (Bohr)",
ylabel="Total energy (Hartree)",
title="Potential energy surface for LiH dissociation",
)
ax.grid()
plt.show()
delta_vqe_energies = []
delta_sparse_adapt_vqe_energies = []
chem_accuracy = []
for i in range(0,len(exact_energies)):
delta1 = vqe_energies[i] - exact_energies[i]
delta_vqe_energies.append(delta1)
delta2 = sparse_adapt_vqe_energies[i] - exact_energies[i]
delta_sparse_adapt_vqe_energies.append(delta2)
chem_accuracy.append(0.0016)
fig, ax = plt.subplots()
ax.plot(r_range, delta_vqe_energies, label='vqe error')
ax.plot(r_range, delta_sparse_adapt_vqe_energies, label='adapt-vqe error')
ax.plot(r_range, chem_accuracy, label='chemical accuracy', linestyle='dashed')
plt.legend(loc='upper left')
ax.set(
xlabel="Bond length (Bohr)",
ylabel="Total energy (Hartree)",
title="VQE error for LiH dissociation",
)
ax.grid()
plt.show()
VQE convergence
# Plot the number of VQE steps to reach convergence per bond value
fig, ax = plt.subplots()
ax.plot(r_range, vqe_steps)
ax.set(
xlabel="Bond length (Bohr)",
ylabel="Number of VQE steps",
title="VQE convergence vs bond length for LiH dissociation",
)
ax.grid()
plt.show()
And, finally, make sure to include the versions of your packages. Specifically, show us the output of qml.about().
Name: PennyLane
Version: 0.40.0
Summary: PennyLane is a cross-platform Python library for quantum computing, quantum machine learning, and quantum chemistry. Train a quantum computer the same way as a neural network.
Home-page: GitHub - PennyLaneAI/pennylane: PennyLane is a cross-platform Python library for quantum computing, quantum machine learning, and quantum chemistry. Train a quantum computer the same way as a neural network.
Author:
Author-email:
License: Apache License 2.0
Location: /home/director/miniforge3/envs/pennylane/lib/python3.12/site-packages
Requires: appdirs, autograd, autoray, cachetools, diastatic-malt, networkx, numpy, packaging, pennylane-lightning, requests, rustworkx, scipy, tomlkit, typing-extensions
Required-by: PennyLane-Catalyst, PennyLane_Lightning, PennyLane_Lightning_GPU
Platform info: Linux-5.4.0-205-generic-x86_64-with-glibc2.31
Python version: 3.12.0
Numpy version: 1.26.4
Scipy version: 1.15.1
Installed devices:
- lightning.qubit (PennyLane_Lightning-0.40.0)
- default.clifford (PennyLane-0.40.0)
- default.gaussian (PennyLane-0.40.0)
- default.mixed (PennyLane-0.40.0)
- default.qubit (PennyLane-0.40.0)
- default.qutrit (PennyLane-0.40.0)
- default.qutrit.mixed (PennyLane-0.40.0)
- default.tensor (PennyLane-0.40.0)
- null.qubit (PennyLane-0.40.0)
- reference.qubit (PennyLane-0.40.0)
- lightning.gpu (PennyLane_Lightning_GPU-0.40.0)
- nvidia.custatevec (PennyLane-Catalyst-0.10.0)
- nvidia.cutensornet (PennyLane-Catalyst-0.10.0)
- oqc.cloud (PennyLane-Catalyst-0.10.0)
- softwareq.qpp (PennyLane-Catalyst-0.10.0)
