I am very new to quantum computing and quantum chemistry, and i am following the tutorial here to design an ADAPT-vqe for LiH molecule. Our of curiosity, I want to try to run this circuit to map out the potential energy surface, but I am not sure how to begin. For every iteration (every time the distance increases), do I have to design a new ADAPT-vqe ansatz, or can i just use the ansatz created by ADAPT-vqe for the ground state energy?

Here is one of my attempt to map PES using ADAPT-vqe but I know I made several mistakes because it looks like the energy is not converging. Here is the snippet of code:

```
from pennylane import qchem
import pennylane as qml
import pennylane.numpy as np
class AdaptVQE:
def __init__(self, symbols: list, geometry: np.array):
"""
:param symbols: Molecule(s)
:param geometry: Nuclear Coordinates
:param active_electrons: Total Active Electrons
:param active_orbitals: Total Active Orbitals
"""
# Getting Hamiltonian and Qubits
self.hamiltonian, self.qubits = qchem.molecular_hamiltonian(
symbols,
geometry,
)
# Pennylane Device
self.dev = qml.device(name="default.qubit", wires=self.qubits)
self.singles, self.doubles = qchem.excitations(active_electrons, self.qubits)
# Hartree-Fock Initialization
self.hf_state = qchem.hf_state(active_electrons, self.qubits)
# Cost Function
self.cost_fn = None
# Selected Gates
self.gate_select = None
# Quantum Gradient Optimizer
self.opt = qml.QNGOptimizer(stepsize=0.05, approx="block-diag")
def circuit_1(self, params: list, wires, excitations) -> None:
qml.BasisState(self.hf_state, wires=wires)
for i, excitation in enumerate(excitations):
if len(excitation) == 4:
qml.DoubleExcitation(params[i], wires=excitation)
else:
qml.SingleExcitation(params[i], wires=excitation)
def circuit_2(self, params: list, wires, excitations, gates_select, params_select) -> None:
qml.BasisState(self.hf_state, wires=wires)
for i, gate in enumerate(gates_select):
if len(gate) == 4:
qml.DoubleExcitation(params_select[i], wires=gate)
elif len(gate) == 2:
qml.SingleExcitation(params_select[i], wires=gate)
for i, gate in enumerate(excitations):
if len(gate) == 4:
qml.DoubleExcitation(params[i], wires=gate)
elif len(gate) == 2:
qml.SingleExcitation(params[i], wires=gate)
def calculate(vqe: AdaptVQE):
vqe = vqe
vqe.cost_fn = qml.ExpvalCost(vqe.circuit_1, vqe.hamiltonian, vqe.dev, optimize=True)
circuit_gradient = qml.grad(vqe.cost_fn, argnum=0)
params = [0.0] * len(vqe.doubles)
grads = circuit_gradient(params, excitations=vqe.doubles)
double_select = [vqe.doubles[i] for i in range(len(vqe.doubles)) if abs(grads[i]) > 1.0e-5]
params_doubles = np.zeros(len(double_select), requires_grad=True)
for i in range(5):
params_doubles = vqe.opt.step(vqe.cost_fn, params_doubles, excitations=double_select)
vqe.cost_fn = qml.ExpvalCost(vqe.circuit_2, vqe.hamiltonian, vqe.dev, optimize=True)
circuit_gradient = qml.grad(vqe.cost_fn, argnum=0)
params = [0.0] * len(vqe.singles)
grads = circuit_gradient(params, excitations=vqe.singles, gates_select=double_select, params_select=params_doubles)
singles_select = [vqe.singles[i] for i in range(len(vqe.singles)) if grads[i] > 1.0e-5]
vqe.cost_fn = qml.ExpvalCost(vqe.circuit_1, vqe.hamiltonian, vqe.dev, optimize=True)
params = np.zeros(len(double_select + singles_select), requires_grad=True)
gate_select = double_select + singles_select
vqe.gate_select = gate_select
max_iteration = 5
print(gate_select)
for i in range(max_iteration):
params, energy = vqe.opt.step_and_cost(vqe.cost_fn, params, excitations=gate_select)
print("Energy: {}".format(energy))
def pes(vqe: AdaptVQE, symbols, coord):
vqe.hamiltonian, vqe.qubits = qchem.molecular_hamiltonian(symbols, coord)
for i in range(5):
vqe.cost_fn = qml.ExpvalCost(vqe.circuit_1, vqe.hamiltonian, vqe.dev, optimize=True)
params, energy = vqe.opt.step_and_cost(vqe.cost_fn, np.zeros(len(vqe.gate_select)), excitations=vqe.gate_select)
print(params, energy)
if __name__ == '__main__':
m = ["H", "H"]
c = np.array([0.0, 0.0, -0.6614, 0.0, 0.0, 0.6614])
v = AdaptVQE(symbols=m, geometry=c)
calculate(v)
r = np.arange(0.5, 5, 0.25)
for j in r:
c = np.array([0.0, 0.0, 0.0, 0.0, 0.0, j])
pes(v, m, c)
```

In the `pes`

function, I attempted to create a new electronic hamiltonian after the coordinates changed and feed it into the ADAPT-vqe that was created for simulating the ground sate energy.