Computing energy of O2 in ground triplet state or first excited singlet state

Hi,

I am trying to compute the energy of O2 in two different electronc states, the ground triplet state T0 and the first singlet excited state S1. The code below uses exact diagonalisation, and ADAPT-VQE. The code below doesn’t give the expected results.

-For the exact diagonalisation section the lowest eigenvalue of the two Hamiltonians is the same. Is it necessary to project the Hamiltonian into a singlet subspace to find the energy of the lowest singlet state? If so is this possible with pennylane ?

-For the ADAPT-VQE section, both calculations give the same energies. It seems that even though the multiplicity was varied and the Hamiltonian used in the calculation is different, the same optimisation was done. Is it necessary to replace the array hf_state by another array to perform a calculation of the energy of the O2 singlet state?

Output for ADAPT-VQE run on T0

Output for ADAPT-VQE run on S1

Code to reproduce these results

import pennylane as qml
from pennylane import numpy as np
from pennylane import qchem
active_electrons = 8 #12 or 8
active_orbitals = 6  #8 or 6
# atomic symbols defining the molecule
symbols = ['O', 'O']
#O2 near equilibrium bond length in bohr
coordinates = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.50]], requires_grad= False)
# Triplet T0 ground state
molecule = qchem.Molecule(symbols, coordinates, mult=3)
H, qubits = qchem.molecular_hamiltonian(molecule, active_electrons=active_electrons, active_orbitals=active_orbitals, method='openfermion')
print ('H size %s qubits %s' % (len(H),qubits))
H_sparse = qml.SparseHamiltonian(H.sparse_matrix(), wires=H.wires)
eigenvals = qml.eigvals(H_sparse)
print ("Exact energy %s " % (eigenvals))
# Singlet S1 state
molecule = qchem.Molecule(symbols, coordinates, mult=1)
H, qubits = qchem.molecular_hamiltonian(molecule, active_electrons=active_electrons, active_orbitals=active_orbitals, method='openfermion')
print ('H size %s qubits %s' % (len(H),qubits))
H_sparse = qml.SparseHamiltonian(H.sparse_matrix(), wires=H.wires)
eigenvals = qml.eigvals(H_sparse, k=10)
print ("Exact energy %s " % (eigenvals))
# Question why exact diagonalisation gives the same energy for a singlet state. Should the Hamiltonian be restricted to a subspace. How?
#ADAPT VQE on T0
# Construct the Molecule object
molecule = qchem.Molecule(symbols, coordinates, mult=3)

# Obtain the qubit Hamiltonian
H, qubits = qchem.molecular_hamiltonian(molecule, active_electrons=active_electrons, active_orbitals=active_orbitals, method='openfermion')
H_sparse = H.sparse_matrix()
singles, doubles = qchem.excitations(active_electrons, qubits)
hf_state = qchem.hf_state(active_electrons, qubits)
print (hf_state)
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
print ('size of operator pool %s' % len(operator_pool))
# 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=False)
    if i % 2 == 0:
        print("n = {:},  E = {:.8f} H, Largest Gradient = {:.3f}".format(i, energy, gradient))
    if gradient < 3e-3:
        break 
    
print ("VQE energy %s, Hamiltonian length %s, number of qubits %s" % (energy, len(H),qubits))
#ADAPT VQE on S1
# Construct the Molecule object
molecule = qchem.Molecule(symbols, coordinates, mult=1)

# Obtain the qubit Hamiltonian
H, qubits = qchem.molecular_hamiltonian(molecule, active_electrons=active_electrons, active_orbitals=active_orbitals, method='openfermion')
H_sparse = H.sparse_matrix()
singles, doubles = qchem.excitations(active_electrons, qubits)
hf_state = qchem.hf_state(active_electrons, qubits)
print (hf_state)
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
print ('size of operator pool %s' % len(operator_pool))
# 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=False)
    if i % 2 == 0:
        print("n = {:},  E = {:.8f} H, Largest Gradient = {:.3f}".format(i, energy, gradient))
    if gradient < 3e-3:
        break 
    
print ("VQE energy %s, Hamiltonian length %s, number of qubits %s" % (energy, len(H),qubits))
#Why getting the same results as for T0. Should ADAPT-VQE be initialised with a different array than hf_state?

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: https://github.com/PennyLaneAI/pennylane
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-216-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)

Hi @Julien ,

Thank you for reporting this.
Let me check and get back to you. In the meantime you can try using PennyLane v0.41.1 just in case.

Hi @Julien ,

As you noticed changing the spin multiplicity isn’t enough to get the first excited state.

Instead you might want to use VQD (for example as shown in the VQD demo) to get that state.

Let me know if you have questions about VQD.

I hope this helps!

1 Like

Hi @CatalinaAlbornoz

Thanks, I will take a look at VQD.
Could you comment why ADAPT-VQE on the S1 openfermion generated Hamiltonian gives the same energy as ADAPT-VQE on the T0 openfermion generated Hamiltonian even though the qubit Hamiltonians have different number of terms (247 vs 383).
Is it because it is possible to use a RHF formalism to compute the MO coefficients of the S1 Hamiltonian ?

Hi @Julien. The backend you are using implements OpenFermion-PySCF functions. Could you please try to generate the Hamiltonian following the steps here to see if you get the same Hamiltonian or not. I suspect that the size of the basis set and the selection of the active space affects the energy gap you compute for the singlet and triplet states. Please let me know if you still get unexpected Hamiltonians.

1 Like