Quantum Self Consistent Equation of Motion AlgorithmSCEOM
¶
The Quantum Self Consistent Equation of Motion (QSCEOM) is an equation-of-motion method to compute excitation energies based on a VQE ground state
where
Eq. (21) can be further simplified for the off-diagonal elements as:
where
AlgorithmSCEOM
is a high level interface of InQuanto to run the SCEOM algorithm for given excitation operators and ansatz with optimal parameters to evaluate the eigenstates.
We take the minimal basis hydrogen molecule as an example to calculate the energy of single and double excited states from the VQE state with the UCCSD ansatz.
The first step is constructing the molecular Hamiltonian and the Fock state in fermion objects. For this, we pull from the inquanto.express
suite of chemical test systems.
from inquanto.express import get_system
hamiltonian, fermion_space, fermion_state = get_system("h2_sto3g.h5")
And then map them into the qubit objects. Here we use the simplified expression of the UCCSD ansatz for the hydrogen molecule with four qubits.
from inquanto.mappings import QubitMappingJordanWigner
from inquanto.ansatzes import TrotterAnsatz
from inquanto.operators import QubitOperator, QubitOperatorList
jw = QubitMappingJordanWigner()
qubit_hamiltonian = jw.operator_map(hamiltonian)
qubit_state = jw.state_map(fermion_state)
ansatz = TrotterAnsatz(
QubitOperatorList.from_list([QubitOperator("Y0 X1 X2 X3", 1j)]), qubit_state
)
Subsequently, we need to define the excitation operators to be applied to our ground state. Here we use single and double excitation operators (defined as expansion_operators
below) that conserve the particle number and the azimuthal spin (i.e.
expansion_operators = fermion_space.construct_single_excitation_operators(
fermion_state
)
expansion_operators += fermion_space.construct_double_excitation_operators(
fermion_state
)
We then define a new Computable object for the SCEOMMatrixComputable
class.
This computable is then passed to AlgorithmSCEOM
.
In general we first need to perform a VQE calculation to obtain the optimized value of the parameter of the ansatz. In this specific example, we already know it
from the previous VQE example (see in AlgorithmVQE
), so we provide the optimized value as an array.
from inquanto.computables.composite import SCEOMMatrixComputable
from inquanto.algorithms import AlgorithmSCEOM
vqe_parameters = ansatz.state_symbols.construct_from_array([-0.10723347230091537])
vqe_state = ansatz.to_CircuitAnsatz(vqe_parameters)
computable = SCEOMMatrixComputable(
space=fermion_space,
fermion_state=fermion_state,
ground_state=vqe_state,
mapping=jw,
hermitian_operator=qubit_hamiltonian,
expansion_operators=expansion_operators,
)
algorithm = AlgorithmSCEOM(
computable_sceom_matrix=computable
)
Here we perform a statevector calculation and we compare the generated energies of the excited states with the results obtained from the exact diagonalization of the hamiltonian.
from inquanto.protocols import SparseStatevectorProtocol
from pytket.extensions.qiskit import AerStateBackend
backend = AerStateBackend()
protocol = SparseStatevectorProtocol(backend)
algorithm.build(protocol=protocol).run()
report = algorithm.generate_report()
print(f'QSCEOM eigenvalues: {report["final_value"]}')
# Compare with the exact diagonalization results
qubit_state = jw.state_map(fermion_state)
e_exact = qubit_hamiltonian.eigenspectrum(qubit_state.single_term.hamming_weight)
print(f'Exact eigvals : {e_exact}')
# TIMER Measure and calculate M matrix elements: BEGINS AT 2025-06-02 14:21:42.425781
# TIMER Measure and calculate M matrix elements: ENDS - DURATION (s): 0.0867673 [0:00:00.086767]
# TIMER Find eigenvalues of M matrix: BEGINS AT 2025-06-02 14:21:42.512587
# TIMER Find eigenvalues of M matrix: ENDS - DURATION (s): 0.0003791 [0:00:00.000379]
QSCEOM eigenvalues: [ 0.552 -0.136 -0.495]
Exact eigvals : [-1.137 -0.495 -0.495 -0.495 -0.136 0.552]
The number of states generated by the QSCEOM method depends on the excitation operators that we use, while the full spectrum of excited states is obtained with the exact diagonalization of the hamiltonian. As discussed above, in this example, QSCEOM only generates excited states that preserve the particle number and the azimuthal spin of the ground state.
The eigenvalues and eigenvectors are stored as final_value()
and final_states()
, respectively.
Numerical instabilities in the matrix diagonalization are handled by removing near linear dependencies in the basis.
We can print the bitstrings and the corresponding coefficients for each excited state that is generated by QSCEOM as shown below:
print(algorithm.print_sceom_states())
0th excited state:
Coefficient Basis state State index Probability
0 0.994256+0.000000j 0011 3 0.988545
1 0.107028+0.000000j 1100 12 0.011455
1th excited state:
Coefficient Basis state State index Probability
0 -0.707107+0.000000j 0110 6 0.5
1 0.707107+0.000000j 1001 9 0.5
2th excited state:
Coefficient Basis state State index Probability
0 -0.707107+0.000000j 0110 6 0.5
1 -0.707107+0.000000j 1001 9 0.5
None
The excited states are sorted according to their energies.
We can also print a dataframe with information about the
print(algorithm.get_dataframe_sceom_analysis())
<Sz> <S2> <GS|ES>
0 -0.0 0.0 -0.0
1 -0.0 0.0 0.0
2 -0.0 2.0 -0.0
In this example, we neglected symmetry in the excited states calculation. If we want to utilize symmetries, we follow the same strategy as outlined above with the exception that
the space
from which we construct the set of expansion_operators
should correspond to the asymmetric expansion_operators
. In addition, we have also implemented a symmetry filter that prevents the calculation of redundant off-diagonal elements of the pointgroup
(which corresponds to the point group of the system) to AlgorithmSCEOM
.
Regarding shot-based measurements, please note that we need to use a ProtocolList
instead of a single protocol because the correlated excited state in the ExpectationValue
of Eq. (21) varies across the elements of the AlgorithmQSE
). Please note that currently
only statevector calculations are supported with AlgorithmSCEOM
. Shot-based calculations can still be performed with the SCEOMMatrixComputable
class.