Source code for qiskit_nature.second_q.algorithms.excited_states_solvers.qeom

# This code is part of a Qiskit project.
#
# (C) Copyright IBM 2020, 2024.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"""The calculation of excited states via the qEOM algorithm."""

from __future__ import annotations

from typing import Any, Callable, Sequence
from enum import Enum
import itertools
import logging


import numpy as np
from scipy import linalg

from qiskit_algorithms import EigensolverResult, MinimumEigensolver
from qiskit_algorithms.list_or_dict import ListOrDict as ListOrDictType
from qiskit_algorithms.observables_evaluator import estimate_observables
from qiskit.circuit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit.primitives import BaseEstimator

from qiskit_nature.second_q.algorithms.ground_state_solvers import GroundStateSolver
from qiskit_nature.second_q.algorithms.excited_states_solvers.excited_states_solver import (
    ExcitedStatesSolver,
)
from qiskit_nature.second_q.mappers import QubitMapper, TaperedQubitMapper
from qiskit_nature.second_q.operators import SparseLabelOp
from qiskit_nature.second_q.problems import (
    BaseProblem,
    ElectronicStructureProblem,
    VibrationalStructureProblem,
    EigenstateResult,
    ElectronicStructureResult,
)
from qiskit_nature.utils import _parallel_map

from .qeom_electronic_ops_builder import build_electronic_ops
from .qeom_vibrational_ops_builder import build_vibrational_ops


logger = logging.getLogger(__name__)


[docs]class EvaluationRule(Enum): """An enumeration of the available evaluation rules for the excited states solvers. This ``Enum`` simply names the available evaluation rules. """ ALL = "all" DIAG = "diag"
def _commutator(op_a: SparsePauliOp, op_b: SparsePauliOp) -> SparsePauliOp: r"""Compute commutator of `op_a` and `op_b`. .. math:: AB - BA. Args: op_a: Operator A. op_b: Operator B. Returns: The computed commutator. """ return (op_a @ op_b - op_b @ op_a).simplify(atol=0) def _double_commutator( op_a: SparsePauliOp, op_b: SparsePauliOp, op_c: SparsePauliOp, sign: bool = False, ) -> SparsePauliOp: r"""Compute symmetric double commutator of `op_a`, `op_b` and `op_c`. See also Equation (13.6.18) in [1]. If `sign` is `False`, it returns .. math:: [[A, B], C]/2 + [A, [B, C]]/2 = (2ABC + 2CBA - BAC - CAB - ACB - BCA)/2. If `sign` is `True`, it returns .. math:: \lbrace[A, B], C\rbrace/2 + \lbrace A, [B, C]\rbrace/2 = (2ABC - 2CBA - BAC + CAB - ACB + BCA)/2. Args: op_a: Operator A. op_b: Operator B. op_c: Operator C. sign: False anti-commutes, True commutes. Returns: The computed double commutator. References: [1]: R. McWeeny. Methods of Molecular Quantum Mechanics. 2nd Edition, Academic Press, 1992. ISBN 0-12-486552-6. """ sign_num = 1 if sign else -1 op_ab = op_a @ op_b op_ba = op_b @ op_a op_ac = op_a @ op_c op_ca = op_c @ op_a op_abc = op_ab @ op_c op_cba = op_c @ op_ba op_bac = op_ba @ op_c op_cab = op_c @ op_ab op_acb = op_ac @ op_b op_bca = op_b @ op_ca res = ( op_abc - sign_num * op_cba + 0.5 * (-op_bac + sign_num * op_cab - op_acb + sign_num * op_bca) ) return res.simplify(atol=0)
[docs]class QEOM(ExcitedStatesSolver): """The calculation of excited states via the qEOM algorithm. This algorithm approximates the excited-state properties of a problem using additional measurements on the ground state provided by a ``GroundStateSolver`` object. The precision of the ``GroundStateSolver.solve`` method for the ground state approximate directly affects the precision of the qEOM algorithm for the same problem. The ``excitations`` are used to build a linear subspace in which an eigenvalue problem for the projected Hamiltonian will be solved. This method typically works well for calculating the lowest-lying excited states of a problem. The excited-state energies are calculated by default in this algorithm for all excited states. Auxiliary observables can be specified to the ``solve`` method along with auxiliary evaluation rules of the ``QEOM`` object. For more details, please refer to https://arxiv.org/abs/1910.12890. The following attributes can be read and updated once the ``QEOM`` object has been constructed. Attributes: excitations: The excitations to be included in the eom pseudo-eigenvalue problem. aux_eval_rules: The rules determining how observables should be evaluated on excited states. tol: The tolerance threshold for the qEOM eigenvalues. """ def __init__( self, ground_state_solver: GroundStateSolver, estimator: BaseEstimator, excitations: str | int | list[int] | Callable[ [int, tuple[int, int]], list[tuple[tuple[int, ...], tuple[int, ...]]], ] = "sd", aux_eval_rules: EvaluationRule | dict[str, list[tuple[int, int]]] | None = None, *, tol: float = 1e-6, ) -> None: """ Args: ground_state_solver: A ``GroundStateSolver`` object. The qEOM algorithm will use this ground state to compute the EOM matrix elements. estimator: The ``BaseEstimator`` to use for the evaluation of the qubit operators at the ground state ansatz. If the internal solver provided to the ``GroundStateSolver`` also uses a ``BaseEstimator`` primitive, you can provide the same estimator instance here. excitations: The excitations to be included in the eom pseudo-eigenvalue problem. :`str`: which contains the types of excitations. Allowed characters are + `s` for singles + `d` for doubles + `t` for triples + `q` for quadruples :`int`: a single, positive integer which denotes the number of excitations (1 == `s`, etc.) :`list[int]`: a list of positive integers generalizing the above :`Callable`: a function which is used to generate the excitations. The callable must take the __keyword__ arguments `num_spin_orbitals` and `num_particles` (with identical types to those explained above) and must return a `list[tuple[tuple[int, ...], tuple[int, ...]]]`. For more information on how to write such a callable refer to the default method :meth:`~qiskit_nature.circuit.library.ansatzes.utils.generate_fermionic_excitations`. aux_eval_rules: The rules determining how observables should be evaluated on excited states. By default, none of the auxiliary operators are evaluated on none of the excited states. :`Enum`: specific predefined rules. Allowed rules are: + ALL to compute all expectation values and all transition amplitudes. + DIAG to only compute expectation values. :`dict[str, list[tuple[int, int]]]`: Dictionary mapping valid auxiliary operator's name to lists of tuple (i, j) specifying the indices of the excited states to be evaluated on. tol: Tolerance threshold for the qEOM eigenvalues. This plays a role when one excited state approaches the ground state, in which case it is best to avoid manipulating very small absolute values. """ self._gsc = ground_state_solver self._estimator = estimator self.excitations = excitations self.aux_eval_rules = aux_eval_rules self.tol = tol self._problem_generated_aux_op_names: set[str] = set() @property def qubit_mapper(self) -> QubitMapper: """Returns the qubit_mapper object defined in the ground state solver.""" return self._gsc.qubit_mapper @property def solver(self) -> MinimumEigensolver: """Returns the solver object defined in the ground state solver.""" return self._gsc.solver def _map_operators( self, operators: SparseLabelOp | ListOrDictType[SparseLabelOp] ) -> SparsePauliOp | ListOrDictType[SparsePauliOp]: if isinstance(self.qubit_mapper, TaperedQubitMapper): mapped_ops = self.qubit_mapper.map_clifford(operators) else: mapped_ops = self.qubit_mapper.map(operators) return mapped_ops def _taper_operators( self, operators: SparsePauliOp | ListOrDictType[SparsePauliOp] ) -> SparsePauliOp | ListOrDictType[SparsePauliOp]: if isinstance(self.qubit_mapper, TaperedQubitMapper): tapered_ops = self.qubit_mapper.taper_clifford(operators, suppress_none=True) else: tapered_ops = operators return tapered_ops
[docs] def get_qubit_operators( self, problem: BaseProblem, aux_operators: dict[str, SparseLabelOp | SparsePauliOp] | None = None, ) -> tuple[SparsePauliOp, dict[str, SparsePauliOp] | None]: """ Gets the operator and auxiliary operators, and transforms the provided auxiliary operators. If the user-provided ``aux_operators`` contain a name which clashes with an internally constructed auxiliary operator, then the corresponding internal operator will be overridden by the user-provided operator. Note that this methods performs a specific treatment of the symmetries required by the qEOM calculation. Args: problem: A class encoding a problem to be solved. aux_operators: Additional auxiliary operators to evaluate. Returns: Tuple of the form (Qubit operator, Auxiliary operators). """ main_operator, aux_second_q_operators = problem.second_q_ops() self._problem_generated_aux_op_names = set(aux_second_q_operators.keys()) # 1. Convert the main operator (hamiltonian) to a Qubit Operator and apply two qubit reduction if isinstance(self.qubit_mapper, TaperedQubitMapper): main_op = self.qubit_mapper.map_clifford(main_operator) else: main_op = self.qubit_mapper.map(main_operator) # 3. Convert the auxiliary operators. # aux_ops set to None if the solver does not support auxiliary operators. aux_ops = None if self.solver.supports_aux_operators(): aux_ops = self._map_operators(aux_second_q_operators) if aux_operators is not None: for name, op in aux_operators.items(): if isinstance(op, (SparseLabelOp)): converted_aux_op = self._map_operators(op) else: converted_aux_op = op if name in aux_ops.keys(): logger.warning( "The key '%s' was already taken by an internally constructed auxiliary" "operator!" "The internal operator was overridden by the one provided manually." "If this was not the intended behavior, please consider renaming" "this operator.", name, ) # The custom op overrides the default op if the key is already taken. aux_ops[name] = converted_aux_op untap_main_op = main_op untap_aux_ops = aux_ops return untap_main_op, untap_aux_ops
[docs] def solve( self, problem: BaseProblem, aux_operators: dict[str, SparseLabelOp | SparsePauliOp] | None = None, ) -> EigenstateResult: """Run the excited-states calculation. Construct and solve the EOM pseudo-eigenvalue problem to obtain the excitation energies and the excitation operators expansion coefficients. Args: problem: A class encoding a problem to be solved. aux_operators: Additional auxiliary operators to evaluate. Returns: An interpreted :class:`~.EigenstateResult`. For more information see also :meth:`~.BaseProblem.interpret`. The :class:`~.EigenstateResult` is constructed from a :class:`~qiskit_nature.second_q.algorithms.excited_states_solvers.qeom.QEOMResult` instance which holds additional information specific to the qEOM problem. """ # 1. Prepare all operators and set the particle number in the qubit mapper ( untap_main_op, # Hamiltonian untap_aux_ops, # Auxiliary observables ) = self.get_qubit_operators(problem, aux_operators) # before we taper our operators we filter the ones which come from the problem internally as # to not trigger a bunch of warnings being raised about overwritten auxiliary operators filtered_aux_ops = { k: v for k, v in untap_aux_ops.items() if k not in self._problem_generated_aux_op_names and k not in aux_operators.keys() } # 2. Run ground state calculation with fully tapered custom auxiliary operators # Note that the solve() method includes the `second_q' auxiliary operators tap_aux_operators = self._taper_operators(filtered_aux_ops) groundstate_result = self._gsc.solve(problem, tap_aux_operators) ground_state = groundstate_result.eigenstates[0] # 3. Prepare the expansion operators for the excited state calculation expansion_basis_data = self._prepare_expansion_basis(problem) # 4. Obtain the representation of the Hamiltonian in the linear subspace h_mat, s_mat, h_mat_std, s_mat_std = self._build_qeom_pseudoeigenvalue_problem( untap_main_op, expansion_basis_data, ground_state ) # 5. Solve the pseudo-eigenvalue problem energy_gaps, expansion_coefs, commutator_metric = self._compute_excitation_energies( h_mat, s_mat ) gammas_square: np.ndarray = np.abs(np.diagonal(commutator_metric)) logger.info("Gamma square = %s", gammas_square) scaling_matrix: np.ndarray = np.diag( np.divide(np.ones_like(gammas_square), np.sqrt(gammas_square)) ) expansion_coefs_rescaled: np.ndarray = expansion_coefs @ scaling_matrix # 6. Evaluate auxiliary operators on the excited states ( aux_operators_eigenvalues, transition_amplitudes, ) = self._evaluate_observables_excited_states( untap_aux_ops, expansion_basis_data, ground_state, expansion_coefs_rescaled, ) result = self._build_qeom_result( problem, groundstate_result, expansion_coefs, energy_gaps, h_mat, s_mat, h_mat_std, s_mat_std, aux_operators_eigenvalues, transition_amplitudes, gammas_square, ) return result
def _build_hopping_ops( self, problem: BaseProblem ) -> tuple[ dict[str, SparsePauliOp], dict[str, list[bool]], dict[str, tuple[tuple[int, ...], tuple[int, ...]]], ]: """Builds the product of raising and lowering operators for a given problem. Args: problem: The problem for which to build out the operators. Raises: NotImplementedError: For an unsupported problem type. Returns: Dict of hopping operators, dict of commutativity types and dict of excitation indices. """ if isinstance(problem, ElectronicStructureProblem): return build_electronic_ops( problem.num_spatial_orbitals, (problem.num_alpha, problem.num_beta), self.excitations, self.qubit_mapper, ) elif isinstance(problem, VibrationalStructureProblem): return build_vibrational_ops( problem.num_modals, self.excitations, self.qubit_mapper, ) else: raise NotImplementedError( "The building of qEOM hopping operators is not yet implemented for a problem of " f"type {type(problem)}" ) def _build_all_eom_operators( self, untap_operator: SparsePauliOp, expansion_basis_data: tuple[dict[str, SparsePauliOp], dict[str, list[bool]], int], ) -> dict[str, SparsePauliOp]: """Building all commutators for Q, W, M, V matrices. Args: untap_operator: Not yet tapered Hamiltonian operator expansion_basis_data: all hopping operators based on excitations_list, key is the string of single/double excitation; value is corresponding operator. Returns: A dictionary that contains the operators for each matrix element. """ untap_hopping_ops, type_of_commutativities, size = expansion_basis_data to_be_computed_list = [] mus, nus = np.triu_indices(size) def _build_one_sector(available_hopping_ops): for idx, m_u in enumerate(mus): n_u = nus[idx] left_op_1 = available_hopping_ops.get(f"E_{m_u}") right_op_1 = available_hopping_ops.get(f"E_{n_u}") right_op_2 = available_hopping_ops.get(f"Edag_{n_u}") to_be_computed_list.append((m_u, n_u, left_op_1, right_op_1, right_op_2)) if ( isinstance(self.qubit_mapper, TaperedQubitMapper) and not self.qubit_mapper.z2symmetries.is_empty() ): combinations = itertools.product( [1, -1], repeat=len(self.qubit_mapper.z2symmetries.symmetries) ) for targeted_tapering_values in combinations: logger.info( "In sector: (%s)", ",".join([str(x) for x in targeted_tapering_values]), ) # remove the excited operators which are not suitable for the sector available_hopping_ops = {} targeted_sector = np.asarray(targeted_tapering_values) == 1 for key, value in type_of_commutativities.items(): if np.all(np.asarray(value) == targeted_sector): available_hopping_ops[key] = untap_hopping_ops[key] _build_one_sector(available_hopping_ops) else: _build_one_sector(untap_hopping_ops) if logger.isEnabledFor(logging.INFO): logger.info("Building all commutators:") results = _parallel_map( self._build_commutator_routine, to_be_computed_list, task_args=(untap_operator,), ) all_matrix_operators = {} for result in results: m_u, n_u, eom_operators = result for index_op, op in eom_operators.items(): if op is not None: all_matrix_operators[f"{index_op}_{m_u}_{n_u}"] = op return all_matrix_operators @staticmethod def _build_commutator_routine( params: list, operator: SparsePauliOp ) -> tuple[int, int, dict[str, SparsePauliOp]]: """Numerically computes the commutator / double commutator between operators. Args: params: list containing the indices of matrix element and the corresponding excitation operators. operator: The hamiltonian. Returns: The indices of the matrix element and the corresponding qubit operator for each of the EOM matrices. """ m_u, n_u, left_op_1, right_op_1, right_op_2 = params if left_op_1 is None or right_op_1 is None and right_op_2 is None: q_mat_op = None w_mat_op = None m_mat_op = None v_mat_op = None else: if right_op_1 is not None: # The sign which we use in the case of the double commutator is arbitrary. In # theory, one would choose this according to the nature of the problem (i.e. # whether it is fermionic or bosonic), but in practice, always choosing the # anti-commutator has proven to be more robust. q_mat_op = -_double_commutator(left_op_1, operator, right_op_1, sign=False) # In the case of the single commutator, we are always interested in the energy # difference of two states. Thus, regardless of the problem's nature, we will # always use the commutator. w_mat_op = -_commutator(left_op_1, right_op_1) q_mat_op = None if len(q_mat_op) == 0 else q_mat_op w_mat_op = None if len(w_mat_op) == 0 else w_mat_op else: q_mat_op = None w_mat_op = None if right_op_2 is not None: # For explanations on the choice of commutation relation, please refer to the # comments above. m_mat_op = _double_commutator(left_op_1, operator, right_op_2, sign=False) v_mat_op = _commutator(left_op_1, right_op_2) m_mat_op = None if len(m_mat_op) == 0 else m_mat_op v_mat_op = None if len(v_mat_op) == 0 else v_mat_op else: m_mat_op = None v_mat_op = None eom_operators = {"q": q_mat_op, "w": w_mat_op, "m": m_mat_op, "v": v_mat_op} return m_u, n_u, eom_operators def _build_eom_matrices( self, gs_results: dict[str, tuple[complex, dict[str, Any]]], size: int ) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """Constructs the H and S matrices from the results on the ground state. Args: gs_results: A ground state result object. size: Size of eigenvalue problem. Returns: H and S matrices and their standard deviation """ mus, nus = np.triu_indices(size) m_mat = np.zeros((size, size), dtype=complex) v_mat = np.zeros((size, size), dtype=complex) q_mat = np.zeros((size, size), dtype=complex) w_mat = np.zeros((size, size), dtype=complex) m_mat_std, v_mat_std, q_mat_std, w_mat_std = 0.0, 0.0, 0.0, 0.0 # evaluate results for m_u, n_u in zip(mus, nus): if gs_results.get(f"q_{m_u}_{n_u}") is not None: q_mat[m_u][n_u] = gs_results[f"q_{m_u}_{n_u}"][0] q_mat_std += gs_results[f"q_{m_u}_{n_u}"][1].get("variance", 0) if gs_results.get(f"w_{m_u}_{n_u}") is not None: w_mat[m_u][n_u] = gs_results[f"w_{m_u}_{n_u}"][0] w_mat_std += gs_results[f"w_{m_u}_{n_u}"][1].get("variance", 0) if gs_results.get(f"m_{m_u}_{n_u}") is not None: m_mat[m_u][n_u] = gs_results[f"m_{m_u}_{n_u}"][0] m_mat_std += gs_results[f"m_{m_u}_{n_u}"][1].get("variance", 0) if gs_results.get(f"v_{m_u}_{n_u}") is not None: v_mat[m_u][n_u] = gs_results[f"v_{m_u}_{n_u}"][0] v_mat_std += gs_results[f"v_{m_u}_{n_u}"][1].get("variance", 0) # these matrices are numpy arrays and therefore have the ``shape`` attribute # Matrix building rules # M.adjoint() = M, V.adjoint() = V # Q.T = Q, W.T = -W q_mat = q_mat + q_mat.T - np.identity(q_mat.shape[0]) * q_mat w_mat = w_mat - w_mat.T - np.identity(w_mat.shape[0]) * w_mat m_mat = m_mat + m_mat.T.conj() - np.identity(m_mat.shape[0]) * m_mat v_mat = v_mat + v_mat.T.conj() - np.identity(v_mat.shape[0]) * v_mat q_mat = np.real(q_mat) w_mat = np.real(w_mat) m_mat = np.real(m_mat) v_mat = np.real(v_mat) q_mat_std = q_mat_std / float(size**2) w_mat_std = w_mat_std / float(size**2) m_mat_std = m_mat_std / float(size**2) v_mat_std = v_mat_std / float(size**2) logger.debug("\nQ:=========================\n%s", q_mat) logger.debug("\nW:=========================\n%s", w_mat) logger.debug("\nM:=========================\n%s", m_mat) logger.debug("\nV:=========================\n%s", v_mat) # Matrix building rules # h_mat = [[M, Q], [P, N]] and s_mat = [[V, W], [T, U]] # N = M.conj() = M.T # P = Q.adjoint() = Q.conj() because Q = Q.T # U = -V.conj() = -V.T # T = W.adjoint() h_mat: np.ndarray = np.block([[m_mat, q_mat], [q_mat.T.conj(), m_mat.T]]) s_mat: np.ndarray = np.block([[v_mat, w_mat], [w_mat.T.conj(), -v_mat.T]]) h_mat_std: np.ndarray = np.array([[m_mat_std, q_mat_std], [q_mat_std, m_mat_std]]) s_mat_std: np.ndarray = np.array([[v_mat_std, w_mat_std], [w_mat_std, v_mat_std]]) return h_mat, s_mat, h_mat_std, s_mat_std def _prepare_expansion_basis( self, problem: BaseProblem ) -> tuple[dict[str, SparsePauliOp], dict[str, list[bool]], int]: """Prepares the basis expansion operators by calling the builder for second quantized operator and applying transformations (Mapping, Reduction, First step of the tapering). Args: problem: the problem for which to build out the operators. Returns: Dict of transformed hopping operators, dict of commutativity types, size of the qEOM problem """ logger.debug("Building expansion basis data...") data = self._build_hopping_ops(problem) hopping_operators, type_of_commutativities, excitation_indices = data size = int(len(list(excitation_indices.keys())) // 2) return hopping_operators, type_of_commutativities, size def _build_qeom_pseudoeigenvalue_problem( self, untap_operator: SparsePauliOp, expansion_basis_data: tuple[dict[str, SparsePauliOp], dict[str, list[bool]], int], reference_state: tuple[QuantumCircuit, Sequence[float]], ) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """Builds the matrices for the qEOM pseudo-eigenvalue problem Args: untap_operator: Not yet tapered hamiltonian expansion_basis_data: Dict of transformed hopping operators, dict of commutativity types, size of the qEOM problem reference_state: Reference state (often the VQE ground state) to be used for the evaluation of EOM operators. Returns: Matrices of the Pseudo-eigenvalue problem H @ X = S @ X @ E with the associated standard deviation errors. """ logger.debug("Build qEOM pseudoeigenvalue problem...") # 1. Build all EOM operators to evaluate on the ground state untap_eom_matrix_ops = self._build_all_eom_operators( untap_operator, expansion_basis_data, ) tap_eom_matrix_ops = self._taper_operators(untap_eom_matrix_ops) # 2. Evaluate all EOM operators on the ground state measurement_results = estimate_observables( self._estimator, reference_state[0], tap_eom_matrix_ops, reference_state[1], ) # 4. Post-process the measurement results to construct eom matrices _, _, size = expansion_basis_data h_mat, s_mat, h_mat_std, s_mat_std = self._build_eom_matrices(measurement_results, size) return h_mat, s_mat, h_mat_std, s_mat_std def _compute_excitation_energies( self, h_mat: np.ndarray, s_mat: np.ndarray ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """Diagonalizing H, S matrices for excitation energies. Args: h_mat : H matrix s_mat : S matrix Returns: 1-D vector stores all energy gap to reference state. 2-D array storing the X and Y expansion coefficients. """ logger.debug("Diagonalizing qeom matrices for excited states...") res = linalg.eig(h_mat, s_mat) # convert nan value into 0 res[0][np.where(np.isnan(res[0]))] = 0.0 # Only the positive eigenvalues are physical. We need to take care # though of very small values # should an excited state approach ground state. Here the small values # may be both negative or # positive. We should take just one of these pairs as zero. # Since we may now have # small values (positive or negative) take the absolute and then threshold zero. logger.debug("... %s", res[0]) # We keep only the negative half of the eigenvalues in decreasing order. Negative eigenvalues # correspond to the excitations whereas positive eigenvalues correspond to de-excitations order = np.argsort(np.real(res[0]))[::-1][len(res[0]) // 2 : :] w = np.real(res[0])[order] logger.debug("Sorted real parts %s", w) w = np.abs(w) w[w < self.tol] = 0 excitation_energies_gap = w expansion_coefs = res[1][:, order] commutator_metric = expansion_coefs.T.conj() @ s_mat @ expansion_coefs return excitation_energies_gap, expansion_coefs, commutator_metric def _build_excitation_operators( self, expansion_basis_data: tuple[ dict[str, SparsePauliOp], dict[str, list[bool]], int, ], reference_state: tuple[QuantumCircuit, Sequence[float]], expansion_coefs_rescaled: np.ndarray, ) -> list[SparsePauliOp]: """Build the excitation operators O_k such that O_k applied on the reference ground state gives the k-th excited state. Args: expansion_basis_data: Dict of transformed hopping operators, dict of commutativity types, size of the qEOM problem. reference_state : Reference ground state expansion_coefs_rescaled: Expansion coefficient matrix X such that H @ X = S @ X @ E and X^dag @ S @ X is the identity Returns: list of excitation operators [Identity, O_1, O_2, ...] """ untap_hopping_ops, _, size = expansion_basis_data tap_hopping_ops = self._taper_operators(untap_hopping_ops) additionnal_measurements = estimate_observables( self._estimator, reference_state[0], tap_hopping_ops, reference_state[1] ) num_qubits = list(untap_hopping_ops.values())[0].num_qubits identity_op = SparsePauliOp(["I" * num_qubits], [1.0]) ordered_keys = [f"E_{k}" for k in range(size)] + [f"Edag_{k}" for k in range(size)] ordered_signs = [1 for k in range(size)] + [-1 for k in range(size)] translated_hopping_ops = {} for key, sign in zip(ordered_keys, ordered_signs): tap_hopping_ops_eval = ( additionnal_measurements.get(key)[0] if additionnal_measurements.get(key) is not None else 0.0 ) translated_hopping_ops[key] = sign * ( untap_hopping_ops[key] - identity_op * tap_hopping_ops_eval ) # From the matrix of coefficients and the vector of basis operators, we create the vector of # excitation operators. hopping_ops_vector = list(translated_hopping_ops.values()) excitations_ops = [ SparsePauliOp.sum( [ expansion_coefs_rescaled[k, i] * hopping_ops_vector[k] for k in range(expansion_coefs_rescaled.shape[0]) ] ).simplify() for i in range(expansion_coefs_rescaled.shape[1]) ] excitations_ops_reduced = [identity_op] + excitations_ops return excitations_ops_reduced def _prepare_excited_states_observables( self, untap_aux_ops: dict[str, SparsePauliOp], operators_reduced: list[SparsePauliOp], size: int, ) -> dict[tuple[str, int, int], SparsePauliOp]: """Prepare the operators O_k^dag @ Aux @ O_l associated to properties of the excited states k,l defined in the aux_eval_rules. By default, the expectation value of all observables on all excited states are evaluated while no transition amplitudes are computed. Args: untap_aux_ops: Dict of auxiliary operators for which properties will be computed. expansion_basis_data: Dict of transformed hopping operators, dict of commutativity types, size of the qEOM problem. operators_reduced: list of excitation operators [Identity, O_1, O_2, ...] size: size of the qEOM problem. Raises: ValueError: For when the aux_eval_rules do not correspond to any previously defined observable and excited state. Returns: Dict of operators of the form O_k^dag @ Aux @ O_l as specified in the constraints. """ indices = np.diag_indices(size + 1) eval_rules: dict[str, list[tuple[Any, Any]]] if self.aux_eval_rules is None: eval_rules = {} elif isinstance(self.aux_eval_rules, dict): eval_rules = self.aux_eval_rules elif self.aux_eval_rules == EvaluationRule.ALL: indices = np.triu_indices(size + 1) aux_names = untap_aux_ops.keys() indices_list = list(zip(indices[0], indices[1])) eval_rules = {aux_name: indices_list for aux_name in aux_names} elif self.aux_eval_rules == EvaluationRule.DIAG: indices = np.diag_indices(size + 1) aux_names = untap_aux_ops.keys() indices_list = list(zip(indices[0], indices[1])) eval_rules = {aux_name: indices_list for aux_name in aux_names} else: raise ValueError("Aux evaluation rules are ill-defined") op_aux_op_dict: dict[tuple[str, int, int], SparsePauliOp] = {} for op_name, indices_constraint in eval_rules.items(): if op_name not in untap_aux_ops.keys(): raise ValueError("Evaluation constrains cannot be satisfied") aux_op = untap_aux_ops[op_name] for i, j in indices_constraint: if i >= len(operators_reduced) or j >= len(operators_reduced): raise ValueError("Evaluation constrains cannot be satisfied") opi, opj = operators_reduced[i], operators_reduced[j] op_aux_op_dict[(op_name, i, j)] = (opi.adjoint() @ aux_op @ opj).simplify() return op_aux_op_dict def _evaluate_observables_excited_states( self, untap_aux_ops: dict[str, SparsePauliOp], expansion_basis_data: tuple[dict[str, SparsePauliOp], dict[str, list[bool]], int], reference_state: tuple[QuantumCircuit, Sequence[float]], expansion_coefs_rescaled: np.ndarray, ) -> tuple[dict[tuple[int, int], dict[str, Any]], dict[tuple[int, int], dict[str, Any]]]: """Evaluate the expectation values and transition amplitudes of the auxiliary operators on the excited states. Custom rules can be used to define which expectation values and transition amplitudes to compute. A typical rule is specified in the form of a dictionary {'hamiltonian':[(1,1)]} Args: untap_aux_ops: Dict of auxiliary operators for which properties will be computed. expansion_basis_data: Dict of transformed hopping operators, dict of commutativity types, size of the qEOM problem. reference_state: Reference ground state. expansion_coefs_rescaled: Expansion coefficient matrix X such that H @ X = S @ X @ E and X^dag @ S @ X is the identity. Returns: Auxiliary operators eigenvalues and transition amplitudes, following the evaluation rules defined as attributes of the qEOM class. """ aux_operators_eigenvalues: dict[tuple[int, int], dict[str, Any]] = {} transition_amplitudes: dict[tuple[int, int], dict[str, Any]] = {} _, _, size = expansion_basis_data if untap_aux_ops is not None: # 1. Build excitation operators O_l such that O_l |0> = |l> excitations_ops_reduced = self._build_excitation_operators( expansion_basis_data, reference_state, expansion_coefs_rescaled ) # 2. Prepare observables O_k^\dag @ Aux @ O_l op_aux_op_dict = self._prepare_excited_states_observables( untap_aux_ops, excitations_ops_reduced, size ) # 3. Measure observables tap_op_aux_op_dict = self._taper_operators(op_aux_op_dict) aux_measurements = estimate_observables( self._estimator, reference_state[0], tap_op_aux_op_dict, reference_state[1] ) # 4. Format aux_operators_eigenvalues indices_diag = np.diag_indices(size + 1) indices_diag_as_list = list(zip(indices_diag[0], indices_diag[1])) for indice in indices_diag_as_list: aux_operators_eigenvalues[indice] = {} for aux_name in untap_aux_ops.keys(): aux_operators_eigenvalues[indice][aux_name] = aux_measurements.get( (aux_name, indice[0], indice[1]), (0.0, {}) ) # 5. Format transition_amplitudes indices_offdiag = np.triu_indices(size + 1, k=1) indices_offdiag_as_list = list(zip(indices_offdiag[0], indices_offdiag[1])) for indice in indices_offdiag_as_list: transition_amplitudes[indice] = {} for aux_name in untap_aux_ops.keys(): transition_amplitudes[indice][aux_name] = aux_measurements.get( (aux_name, indice[0], indice[1]), (0.0, {}) ) return aux_operators_eigenvalues, transition_amplitudes def _build_qeom_result( self, problem, groundstate_result, expansion_coefs, energy_gaps, h_mat, s_mat, h_mat_std, s_mat_std, aux_operators_eigenvalues, transition_amplitudes, gammas_square, ) -> ElectronicStructureResult: qeom_result = QEOMResult() qeom_result.ground_state_raw_result = groundstate_result.raw_result qeom_result.expansion_coefficients = expansion_coefs qeom_result.excitation_energies = energy_gaps qeom_result.h_matrix = h_mat qeom_result.s_matrix = s_mat qeom_result.h_matrix_std = h_mat_std qeom_result.s_matrix_std = s_mat_std qeom_result.gamma_square = gammas_square qeom_result.aux_operators_evaluated = list(aux_operators_eigenvalues.values()) qeom_result.transition_amplitudes = transition_amplitudes groundstate_energy_reference = groundstate_result.eigenvalues[0] excited_eigenenergies = energy_gaps + groundstate_energy_reference qeom_result.eigenvalues = np.append( groundstate_result.eigenvalues[0], excited_eigenenergies ) eigenstate_result = EigenstateResult.from_result(qeom_result) result = problem.interpret(eigenstate_result) return result
[docs]class QEOMResult(EigensolverResult): """The results class for the qEOM algorithm. For more details about the definitions, please refer to https://arxiv.org/abs/1910.12890. Attributes: ground_state_raw_result (EigenstateResult): The raw results of the ground state eigensolver. excitation_energies (np.ndarray): The excitation energies approximated by the qEOM algorithm. expansion_coefficients (np.ndarray): The expansion coefficients matrix of the excitation operators onto the set of basis operators spanning the linear qEOM subspace. h_matrix (np.ndarray): Matrix representing the Hamiltonian in the qEOM subspace. Because of our choice for the expansion basis, the two square sub-matrices on the diagonal are related by a transposition and the two submatrices on the anti diagonal are hermitian conjugates. s_matrix (np.ndarray): Matrix representing the geometry of the qEOM subspace. Because of our choice for the expansion basis, the two square submatrices on the diagonal are related by a transposition (with a sign) and the two submatrices on the anti diagonal are hermitian conjugates. h_matrix_std (np.ndarray): 2 by 2 matrix representing the sums of standard deviations in the four square submatrices of H. s_matrix_std (np.ndarray): 2 by 2 matrix representing the sums of standard deviations in the four square submatrices of S. transition_amplitudes (list[ListOrDictType[tuple[complex, dict[str, Any]]]): Transition amplitudes of the auxiliary operators computed following the evaluation rules specified when the qEOM class was created. """ def __init__(self) -> None: super().__init__() self.ground_state_raw_result = None self.excitation_energies: np.ndarray | None = None self.expansion_coefficients: np.ndarray | None = None self.eigenvalues: np.ndarray | None = None self.h_matrix: np.ndarray | None = None self.s_matrix: np.ndarray | None = None self.h_matrix_std: np.ndarray = np.zeros((2, 2)) self.s_matrix_std: np.ndarray = np.zeros((2, 2)) self.transition_amplitudes: list[ ListOrDictType[tuple[complex, dict[str, Any]]] ] | None = None self.gamma_square: np.ndarray = None @property def m_matrix(self) -> np.ndarray | None: """returns the M matrix""" if self.h_matrix is None: return None return self.h_matrix[: len(self.h_matrix) // 2, : len(self.h_matrix) // 2] @property def v_matrix(self) -> np.ndarray | None: """returns the V matrix""" if self.s_matrix is None: return None return self.s_matrix[: len(self.s_matrix) // 2, : len(self.s_matrix) // 2] @property def q_matrix(self) -> np.ndarray | None: """returns the Q matrix""" if self.h_matrix is None: return None return self.h_matrix[len(self.h_matrix) // 2 :, : len(self.h_matrix) // 2] @property def w_matrix(self) -> np.ndarray | None: """returns the W matrix""" if self.s_matrix is None: return None return self.s_matrix[len(self.s_matrix) // 2 :, : len(self.s_matrix) // 2] @property def m_matrix_std(self) -> float: """returns the M matrix standard deviation""" if np.isclose(self.h_matrix_std[0, 0], 0.0): return 0.0 return self.h_matrix_std[0, 0] @property def v_matrix_std(self) -> float: """returns the V matrix standard deviation""" if np.isclose(self.s_matrix_std[0, 0], 0.0): return 0.0 return self.s_matrix_std[0, 0] @property def q_matrix_std(self) -> float: """returns the Q matrix standard deviation""" if np.isclose(self.h_matrix_std[0, 1], 0.0): return 0.0 return self.h_matrix_std[0, 0] @property def w_matrix_std(self) -> float: """returns the W matrix standard deviation""" if np.isclose(self.s_matrix_std[0, 1], 0.0): return 0.0 return self.s_matrix_std[0, 0]