Excited states methods

Beyond the ground state, the excited states of molecules are often of interest in applications such as light harvesting and bioimaging as well. In this example, we demonstrate a relatively simple excited state approach, the folded spectrum Hamiltonian (REF!!)

As usual, we start by first defining and building the molecule of interest:

from qibochem.driver import Molecule

h2 = Molecule(xyz_file="../tests/data/h2.xyz")
h2.run_pyscf()

hamiltonian = h2.hamiltonian()
eigenvalues = hamiltonian.eigenvalues()
print("Eigenvalues of the H2 Hamiltonian:")
print(eigenvalues)

Eigenvalues of the H2 Hamiltonian:
[-1.13618945 -0.52188556 -0.52188556 -0.47845306 -0.47845306 -0.47845306
 -0.40318375 -0.40318375 -0.1204519   0.30766775  0.30766775  0.44908566
  0.44908566  0.5833141   0.75596744  1.01608717]

The folded spectrum Hamiltonian requires an additional hyperparameter: \(\omega\), the number to `fold’ the original Hamiltonian about. Ideally, \(\omega\) should be as close to the value of the target excited state as possible. Here, we’ll use the exact ground state, which we will obtain by running VQE normally, as a reference starting point for \(\omega\).

import numpy as np

from qibo.models import VQE

from qibochem.ansatz import ucc_ansatz

# We'll start by trying out the UCC circuit ansatz first
ucc_circuit = ucc_ansatz(h2)

# Run VQE to find the ground state for the original molecular Hamiltonian
vqe = VQE(ucc_circuit, hamiltonian)
ucc_parameters = np.random.rand(len(ucc_circuit.get_parameters()))
vqe_ground, ucc_parameters, _result = vqe.minimize(ucc_parameters)
print(f"VQE result: {vqe_ground:.9f}")

VQE result: -1.136189454

Does the VQE result match the lowest eigenvalue from the above cell? For the small \(H_2\) molecule, VQE shouldn’t have any problem finding the exact ground state.

Now, let’s try using the (ground state) VQE result in the folded spectrum Hamiltonian. To try and find the lowest excited state, we will set the hyperparameter \(\omega\) to be slightly less negative than the VQE ground state, so that we can target the first excited state. Then, we carry out VQE with the folded spectrum Hamiltonian to obtain the circuit parameters that yield the ‘ground state’.

# Define the folded spectrum Hamiltonian
omega = vqe_ground + 0.4  # Hyperparameter
fs_hamiltonian = h2.fs_hamiltonian(omega)

# Run VQE to find the ground state for the folded spectrum Hamiltonian
fs_vqe = VQE(ucc_circuit, fs_hamiltonian)
fs_parameters = np.random.rand(len(ucc_circuit.get_parameters()))
vqe_fs, fs_parameters, _result = fs_vqe.minimize(fs_parameters)
print(f"VQE folded spectrum result: {vqe_fs:.9f}")

VQE folded spectrum result: 0.066428051

Lastly, the expectation value of the original molecular Hamiltonian is calculated using the circuit parameters that were obtained after running VQE with the folded spectrum Hamiltonian. If everything works, this expectation value should give us an excited state of the molecular Hamiltonian.

from qibochem.measurement import expectation

ucc_circuit.set_parameters(fs_parameters)
print(f"Expectation value: {expectation(ucc_circuit, hamiltonian):.8f} (w.r.t. the original Hamiltonian)")

Expectation value: -0.47845306 (w.r.t. the original Hamiltonian)

Comparing this result to the initial set of eigenvalues, we can see that the folded spectum Hamiltonian approach successfully found a higher eigenvalue of the molecular Hamiltonian (\(-0.478...\) vs \(-1.136...\)).

However, one might ask why wasn’t the obtained solution the second (degenerate) eigenvalue (\(-0.521...\))? The answer lies in the circuit ansatz that was used here. More specifically, the first two \(X\) gates in the UCCSD circuit ansatz define the quantum simulation to have exactly two electrons, restricting VQE to searching this particular subspace of the entire Hilbert space. In contrast, the second (degenerate) eigenvalue refers to the case in which there is only one electron in the quantum simulation, which is why the default UCCSD circuit ansatz is unable to locate this solution.

To show this, let’s make some slight modifications to the UCCSD circuit ansatz, and run a normal ground state VQE with it. Essentially, we’ll just set the initial Hartree-Fock reference circuit (and by extension, the VQE search space) to have only one electron.

from qibochem.ansatz import hf_circuit

# Changing the HF reference circuit to have only one electron, and make sure it isn't included in the ucc_ansatz function
ucc_circuit2 = hf_circuit(h2.nso, n_electrons=1) + ucc_ansatz(h2, include_hf=False)

# Run VQE to find the ground state w.r.t. the modified UCCSD circuit ansatz
vqe2 = VQE(ucc_circuit2, hamiltonian)
ucc_parameters2 = np.random.rand(len(ucc_circuit2.get_parameters()))
vqe_ground2, ucc_parameters2, _result = vqe2.minimize(ucc_parameters2)
print(f"VQE result: {vqe_ground2:.9f} (with a one-electron reference HF circuit)")

VQE result: -0.521885562 (with a one-electron reference HF circuit)

So that’s an example of how the circuit ansatz can affect the obtained solution in VQE.

Lastly, what happens if we use a circuit ansatz that doesn’t conserve the number of electrons in the quantum simulation; e.g. the hardware-efficient ansatz?

from qibochem.ansatz import he_circuit

hardware_efficient_circuit = he_circuit(h2.nso, n_layers=3, coupling_gates="CNOT")
hardware_efficient_circuit.draw()

0: ─RY─RZ─o─────X─RY─RZ─o─────X─RY─RZ─o─────X─
1: ─RY─RZ─X─o───|─RY─RZ─X─o───|─RY─RZ─X─o───|─
2: ─RY─RZ───X─o─|─RY─RZ───X─o─|─RY─RZ───X─o─|─
3: ─RY─RZ─────X─o─RY─RZ─────X─o─RY─RZ─────X─o─

# Run VQE with the hardware-efficient circuit ansatz and the folded spectrum Hamiltonian
fs_he_vqe = VQE(hardware_efficient_circuit, fs_hamiltonian)
fs_he_parameters = np.random.rand(len(hardware_efficient_circuit.get_parameters()))
vqe_he_fs, fs_he_parameters, _result = fs_he_vqe.minimize(fs_he_parameters)

# Expectation value of the molecular Hamiltonian
hardware_efficient_circuit.set_parameters(fs_he_parameters)
print(f"Expectation value of molecular Hamiltonian: {expectation(hardware_efficient_circuit, hamiltonian)}")

Expectation value of molecular Hamiltonian: -0.5218851719525988

This time, we should be able to find the second lowest eigenvalue for the molecular Hamiltonian (\(-0.521...\)). (If not, try re-running VQE in the cell immediately above)

Finding a good value of \(\omega\)

As explained above, \(\omega\) is a hyperparameter in the folded spectrum Hamiltonian approach. Unfortunately, there is no easy way of pre-determining a ‘good’ value of \(\omega\) for the excited state of interest. At best, one can use chemical intuition to make an educated guess for \(\omega\); but if one has completely no pre-exising knowledge or experience regarding the system of interest, then a brute force approach of scanning through a range of values of \(\omega\).

An example of how to carry out this scan is given below:

# Hyperparameters to perform the scanning of \omega
delta_step = 0.05
delta_max_multiplier = 1.0
max_iterations = 10

delta_multiplier = 0.1  # Fraction of ground state energy for determining the folding

# Perform the scan
for _i in range(max_iterations):
    delta = delta_multiplier * abs(vqe_ground)
    omega = vqe_ground + delta
    fs_ham = h2.fs_hamiltonian(omega, hamiltonian)

    # Construct and run VQE
    scan_vqe = VQE(ucc_circuit, fs_ham)
    scan_vqe_parameters = np.random.rand(len(ucc_circuit.get_parameters()))
    scan_vqe_min, scan_vqe_parameters, _result = scan_vqe.minimize(scan_vqe_parameters)

    # Print some results
    ucc_circuit.set_parameters(scan_vqe_parameters)
    energy = expectation(ucc_circuit, hamiltonian)
    print(
        f"[FS-VQE] Iteration {_i+1}: delta_multiplier={delta_multiplier:.4f}, delta={delta:.6f}, folded_energy={scan_vqe_min:.8f}, "
        f"Electronic energy: {energy:.7f}"
    )

    if scan_vqe_min < delta**2:
        print(f"\nExcited state found: E1 = {energy:.7f}")
        break
    delta_multiplier += delta_step
    if delta_multiplier > delta_max_multiplier:
        print(f"Failed to find an excited state!")
        break

[FS-VQE] Iteration 1: delta_multiplier=0.1000, delta=0.113619, folded_energy=0.01290926, Electronic energy: -1.1361895
[FS-VQE] Iteration 2: delta_multiplier=0.1500, delta=0.170428, folded_energy=0.02904585, Electronic energy: -1.1361895
[FS-VQE] Iteration 3: delta_multiplier=0.2000, delta=0.227238, folded_energy=0.05163706, Electronic energy: -1.1361895
[FS-VQE] Iteration 4: delta_multiplier=0.2500, delta=0.284047, folded_energy=0.08068301, Electronic energy: -1.1361889
[FS-VQE] Iteration 5: delta_multiplier=0.3000, delta=0.340857, folded_energy=0.10041266, Electronic energy: -0.4784531

Excited state found: E1 = -0.4784531

Longer simulations

Lastly, let’s try to perform a similar scan with the hardware efficient ansatz, and see if we can find all the eigenvalues of the molecular Hamiltonian.

In other words, we’re trying to find all possible excited states, without any considerations regarding the number of electrons in the system.

# Define a function to find the eigenvalue w.r.t. \omega

def find_excited_state(hamiltonian, omega, circuit):
    """
    For a given Hamiltonian, construct the folded spectum Hamiltonian w.r.t. omega, and run VQE to find the excited state
    for the original Hamiltonian

    Args:
        hamiltonian: Original Hamiltonian of interest
        omega: Value to fold the original Hamiltonian about
        circuit: Circuit ansatz for running VQE

    Returns:
        float: Result of folded spectrum Hamiltonian approach
    """
    dummy = Molecule()  # Dummy Molecule for calling the fs_hamiltonian() method
    fs_ham = dummy.fs_hamiltonian(omega, hamiltonian)

    fs_ham_min = None
    optimal_parameters = None
    # Run VQE a number of times to reduce randomness
    for _ in range(3):
        fs_vqe = VQE(circuit, fs_ham)
        fs_parameters = np.random.rand(len(circuit.get_parameters()))
        vqe_fs, fs_parameters, _result = fs_vqe.minimize(fs_parameters)
        if fs_ham_min is None or vqe_fs < fs_ham_min:
            fs_ham_min = vqe_fs
            optimal_parameters = fs_parameters.copy()

    # Expectation value of the original Hamiltonian with the optimal circuit parameters
    circuit.set_parameters(fs_parameters)
    return expectation(circuit, hamiltonian)

Warning: The following cell might take a few minutes to run!

for shift in np.arange(0.0, 2.1, 0.2):
    omega = vqe_ground + shift
    print(f"Omega: {omega:.3f} - Result: {find_excited_state(hamiltonian, omega, hardware_efficient_circuit)}")

Omega: -1.136 - Result: -1.1361894540659196
Omega: -0.936 - Result: -1.1361048278357495
Omega: -0.736 - Result: -0.5216290656815866
Omega: -0.536 - Result: -0.4784525261646579
Omega: -0.336 - Result: -0.40319651894602726
Omega: -0.136 - Result: -0.12045190380706042
Omega: 0.064 - Result: -0.12045165045058437
Omega: 0.264 - Result: 0.4490873502350595
Omega: 0.464 - Result: 0.44908639833377056
Omega: 0.664 - Result: 0.6693673680293584
Omega: 0.864 - Result: 0.7559673673329969

Compare the results to the eigenvalues obtained by diagonalizing the Hamiltonian directly. Did we succeed in finding all of them?