Demo VQE in different spin sectors doesn't seem to work for O2

Demos
1

I have successfully reproduced the demo ‘VQE in different spin sectors’ here

However I have run into difficulties attempting to replace the input H2 by O2.

At first the results for O2 triplet seem to make sense

import numpy as np
from pennylane import qchem
active_electrons = 2 #12 or 8
active_orbitals = 3  #8 or 6
symbols = ["O", "O"]
coordinates = np.array([0.0, 0.0, 0.0000, 0.0, 0.0, 2.5000])
molecule = qchem.Molecule(symbols, coordinates, mult=3)

import pennylane as qml

H, qubits = qml.qchem.molecular_hamiltonian(molecule, active_electrons=active_electrons, active_orbitals=active_orbitals, method='openfermion')
print("Number of qubits = ", qubits)
#print("The Hamiltonian is ", H)
H_sparse = qml.SparseHamiltonian(H.sparse_matrix(), wires=H.wires)
eigenvals = qml.eigvals(H_sparse, k=20)
print ("Exact energy %s " % (eigenvals))


electrons = active_electrons
S2 = qml.qchem.spin2(electrons, qubits)
#print(S2)

hf = qml.qchem.hf_state(electrons, qubits)
# Ground state triplet looks like this in a (2e,3o) active space
mystate = np.array([1,0,1,0,0,0])
print(mystate)

# The below assumes a hf state that is the singlet to generate excitations so may not be appropriate
singles, doubles = qml.qchem.excitations(electrons, qubits, delta_sz=0)
print(singles)
print(doubles)

def circuit(params, wires):
    qml.AllSinglesDoubles(params, wires, mystate, singles, doubles)

dev = qml.device("lightning.qubit", wires=qubits)
@qml.qnode(dev, interface="jax")
def cost_fn(params):
    circuit(params, wires=range(qubits))
    drawer = qml.draw(circuit)
    return qml.expval(H)

@qml.qnode(dev, interface="jax")
def S2_exp_value(params):
    circuit(params, wires=range(qubits))
    return qml.expval(S2)

def total_spin(params):
    return -0.5 + np.sqrt(1 / 4 + S2_exp_value(params))

from jax import random
import optax

opt = optax.sgd(learning_rate=0.8) # sgd stands for StochasticGradientDescent
key = random.PRNGKey(0)
init_params = random.normal(key, shape=(len(singles) + len(doubles),)) * np.pi
print (init_params)

# VQE
import catalyst
max_iterations = 100

prev_energy = 0.0
@qml.qjit
def update_step(i, params, opt_state):
    """Perform a single gradient update step"""
    grads = qml.grad(cost_fn)(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(max_iterations):
    params, opt_state = update_step(i, params, opt_state)
    energy = cost_fn(params)
    if (i % 10 == 0):
        print(f"iteration {i} energy {energy:.8f} Ha") 

print("\n" f"Final value of the ground-state energy = {energy:.8f} Ha")
print("\n" f"Optimal value of the circuit parameters = {params}")
final_spin = total_spin(params)
print("\n" f"The ground state spin is  = {final_spin:.5f}")

The above code generates

Exact energy [-147.5092067 -147.38376755 -147.61107466 -147.61107466 -147.55836291
-147.55836291 -147.61107466 -147.17599662 -147.17599662 -147.24187354
-147.38376755 -147.38376755 -147.38376755 -147.24187354 -147.17599662
-147.24187354 -147.17599662 -147.24187354 -147.24187354 -147.24187354]
…
iteration 0 energy -147.51428140 Ha
iteration 10 energy -147.60551214 Ha
iteration 20 energy -147.61084304 Ha
iteration 30 energy -147.61106517 Ha
iteration 40 energy -147.61107427 Ha
iteration 50 energy -147.61107464 Ha
iteration 60 energy -147.61107466 Ha
iteration 70 energy -147.61107466 Ha
iteration 80 energy -147.61107466 Ha
iteration 90 energy -147.61107466 Ha

Final value of the ground-state energy = -147.61107466 Ha

Optimal value of the circuit parameters = [ 2.5405380e-01 -1.6501926e-07 -1.1801571e+00 5.2432365e+00
-4.0081062e+00 6.6576648e+00 -2.6961534e+00 3.5518634e+00]

The ground state spin is = 1.00000

Which suggests VQE found the desired lowest energy triplet state.

The code below was then run to find the energy of the lowest lying singlet state

# appropriate if ground state is a triplet?
singles, doubles = qml.qchem.excitations(electrons, qubits, delta_sz=-1)
print(singles)
print(doubles)
mystate = np.array([1,1,0,0,0,0])
print(mystate)


def circuit(params, wires):
    qml.AllSinglesDoubles(params, wires, mystate, singles, doubles)
@qml.qnode(dev, interface="jax")
def cost_fn(params):
    circuit(params, wires=range(qubits))
    return qml.expval(H)
@qml.qnode(dev, interface="jax")
def S2_exp_value(params):
    circuit(params, wires=range(qubits))
    return qml.expval(S2)

opt = optax.sgd(learning_rate=0.8) # sgd stands for StochasticGradientDescent
key = random.PRNGKey(0)
init_params = random.normal(key, shape=(len(singles) + len(doubles),)) * np.pi

max_iterations = 100

@qml.qjit
def update_step(i, params, opt_state):
    """Perform a single gradient update step"""
    grads = qml.grad(cost_fn)(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

print (init_params)
for i in range(max_iterations):
    params, opt_state = update_step(i, params, opt_state)
    energy = cost_fn(params)
    if (i % 10 == 0):
        print(f"iteration {i} energy {energy:.8f} Ha") 

print("\n" f"Final value of the energy = {energy:.8f} Ha")
print("\n" f"Optimal value of the circuit parameters = {params}")
final_spin = total_spin(params)
print("\n" f"The ground state spin is  = {final_spin:.5f}")

Which gives

iteration 0 energy -147.45650544 Ha
iteration 10 energy -147.56103809 Ha
iteration 20 energy -147.58332820 Ha
iteration 30 energy -147.58464770 Ha
iteration 40 energy -147.58471515 Ha
iteration 50 energy -147.58471859 Ha
iteration 60 energy -147.58471877 Ha
iteration 70 energy -147.58471878 Ha
iteration 80 energy -147.58471878 Ha
iteration 90 energy -147.58471878 Ha

Final value of the energy = -147.58471878 Ha

Optimal value of the circuit parameters = [ 8.99984130e-08 -3.14159332e+00 3.14159184e+00]

The ground state spin is = 0.61803

The energy doesn’t match any eigenvalues from the exact diagonalisation, and the spin is contaminated.

My questions)

  • Is the use of qml.qchem.excitations flawed given the electronic structure of the ground state of O2?
  • Shouldn’t the molecular hamiltonian be recomputed before running VQE on the singlet state?
  • How could I update the above code to give the correct behavior?
2

Hi @Julien ,

Let me check and get back to you on this.

3

Hi Julien,

Thank you for your interest in the demo “VQE in different spin sectors”. Please, see the answers to your questions below.

  • Is the use of qml.qchem.excitations flawed given the electronic structure of the ground state of O2?

No, you can still use it to generate excitations in different sectors of the total spin projection (S_z).

  • Shouldn’t the molecular hamiltonian be recomputed before running VQE on the singlet state?

No, this should not be necessary.

  • How could I update the above code to give the correct behavior?

The function qml.qchem.excitations assumes a closed-shell reference state. To generate excitations corresponding to singly- and doubly-excited configurations with S_z=1 you can specify \Delta S_z=1 :

singles, doubles = qml.qchem.excitations(electrons, qubits, delta_sz=1)

Then, you can define the VQE ansatz as

hf = qml.qchem.hf_state(electrons, qubits)
def circuit(params, wires):
    qml.AllSinglesDoubles(params, wires, hf, singles, doubles)

Note that the initial state of the qubit register is still defined as the closed-shell Hartree-Fock state for consistency. By running VQE using this ansatz you should be able to find the ground state with spin quantum numbers S=1, S_z=1.

Similarly, to find eigenstates in the spin sector with S_z=0 you can proceed as follows

singles, doubles = qml.qchem.excitations(electrons, qubits, delta_sz=0)

Thus, the VQE algorithm will find the ground state with S_z=0. A subtlety here is that the algorithm may end up finding the state with S=1, S_z=0 which is not the lowest-lying singlet state, but a spin sublevel of the degenerate (triplet (S=1)) ground state. In this case, you may need to modify the cost function to help the optimizer finding the lowest-energy singlet state (S=0). Please, see the demo “How to implement VQD with PennyLane”.

I hope this helps. Please, do not hesitate to get back to us if you have further questions.

Regards,

Alain.