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 \(|\Psi_{\mathrm{VQE}}\rangle = U(\theta^{\mathrm{opt}})|\Psi_{\mathrm{HF}}\rangle\). [16] It satisfies the vacuum annihilation condition which ensures that the ground state cannot be de-excited. The excitation energies are obtained with classical diagonalization of the matrix \(M\) whose elements are evaluated with a quantum computer as:
where \(\hat{H}\) is the Hamiltonian matrix, \(\hat{S} =U(\theta^{\mathrm{opt}}) \hat{G} U^\dagger(\theta^{\mathrm{opt}})\) are the self-consistent excitation operators, and \(E_{\mathrm{gr}}\) is the VQE ground state energy. \(\hat{G}\) are the excitation operators that construct the SCEOM subspace.
Eq. (21) can be further simplified for the off-diagonal elements as:
where \(M_{i+j, i+j} = \langle \Psi_{\mathrm{HF}}|\frac{1}{\sqrt{2}} \left(\hat{G}_i + \hat{G}_j \right)^\dagger U^\dagger(\theta^{\mathrm{opt}}) \hat{H} U(\theta^{\mathrm{opt}}) \frac{1}{\sqrt{2}} \left(\hat{G}_i + \hat{G}_j \right) |\Psi_{\mathrm{HF}}\rangle\)
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. \(\langle \hat{S}_z \rangle\)).
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 \(M\) matrix in Eq. (21), using 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 2024-09-23 17:04:36.870603
# TIMER Measure and calculate M matrix elements: ENDS - DURATION (s): 0.1020898 [0:00:00.102090]
# TIMER Find eigenvalues of M matrix: BEGINS AT 2024-09-23 17:04:36.972742
# TIMER Find eigenvalues of M matrix: ENDS - DURATION (s): 0.0004657 [0:00:00.000466]
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:
Coefficients Basis states State index Probabilities
0 0.994256+0.000000j 0011 3 0.988545
1 0.107028+0.000000j 1100 12 0.011455
1th excited state:
Coefficients Basis states State index Probabilities
0 0.707107+0.000000j 0110 6 0.5
1 -0.707107+0.000000j 1001 9 0.5
2th excited state:
Coefficients Basis states State index Probabilities
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 \(\langle \hat{S}_z \rangle\), \(\langle \hat{S}^2 \rangle\), and the overlap with the VQE ground state for each excited state.
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 \(C_1\) group so
that we include the complete set of expansion_operators. In addition, we have also implemented a symmetry filter that prevents the calculation of redundant off-diagonal elements of the \(M\) matrix
based on comparison of the irreducible representation of the correlated excited states \(i,j\). This feature is optional and it can be used by passing an additional string argument 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 \(M\) matrix. This results to a higher number of measurements compared to QSE (see AlgorithmQSE). Please note that currently
only statevector calculations are supported with AlgorithmSCEOM. Shot-based calculations can still be performed with the SCEOMMatrixComputable class.