How to run PES simulation using ADAPT-vqe?

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(

        # Pennylane Device = qml.device(name="default.qubit", wires=self.qubits), 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)
                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,, 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,, optimize=True)
    circuit_gradient = qml.grad(vqe.cost_fn, argnum=0)
    params = [0.0] * len(

    grads = circuit_gradient(params,, gates_select=double_select, params_select=params_doubles)
    singles_select = [[i] for i in range(len( if grads[i] > 1.0e-5]

    vqe.cost_fn = qml.ExpvalCost(vqe.circuit_1, vqe.hamiltonian,, 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
    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,, 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)
    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.

Hi @Hatedfate, it’s great that you’re trying out the adapt-vqe demo.

In order for it to be easier to help you, could you please point out the specific functions that you have changed with respect to the adapt-vqe demo?

There can be various reasons why it doesn’t seem to be converging so this can help us get a sense of the possible problem.

Hello, I actually did not change anything as shown in the demo. Instead, I tried to use the circuit made by adapt-vqe to map out a potential energy surface, but this did not converge as I thought it would. I assume the problem is in the pes() function because I tried something similar in Jupyter Notebook and it did in fact converge. Highly appreciate your response!

Also, as a side note, I am learning teh Quantum Natural Gradient Descent algorithm and I am curious if there is a mathematical derivation of its complexity (i.e. convergence rate/big o notation). l.m.k. if i am misunderstanding anything because I am just a beginner at the moment. Thanks.

Hi @Hatedfate, you need to define the variable active_electrons.

In the Adapt VQE demo you will notice that active_electrons=2.

This should help. I also encourage you to increase the number of iterations.

Please let me know if this helps!

1 Like