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)