diff --git a/.pylintdict b/.pylintdict index 119b1949cb..b4866d9588 100644 --- a/.pylintdict +++ b/.pylintdict @@ -218,6 +218,7 @@ fcompiler fd fermionic fermionicoperator +fermionictransformation fermionicdriver fermions fi diff --git a/qiskit/aqua/algorithms/eigen_solvers/numpy_eigen_solver.py b/qiskit/aqua/algorithms/eigen_solvers/numpy_eigen_solver.py index 71cbac3186..9a4db287d8 100755 --- a/qiskit/aqua/algorithms/eigen_solvers/numpy_eigen_solver.py +++ b/qiskit/aqua/algorithms/eigen_solvers/numpy_eigen_solver.py @@ -285,6 +285,8 @@ def _run(self): self._ret['eigvecs'] = np.array(eigvecs) self._ret['eigvals'] = np.array(eigvals) self._ret['energies'] = np.array(energies) + # conversion to np.array breaks in case of aux_ops + self._ret['aux_ops'] = aux_ops self._k = k_orig diff --git a/qiskit/aqua/operators/legacy/common.py b/qiskit/aqua/operators/legacy/common.py index b661199ace..4c31099fda 100644 --- a/qiskit/aqua/operators/legacy/common.py +++ b/qiskit/aqua/operators/legacy/common.py @@ -384,7 +384,7 @@ def evolution_instruction(pauli_list, evo_time, num_time_slices, return qc.to_instruction() -def commutator(op_a, op_b, op_c=None, threshold=1e-12): +def commutator(op_a, op_b, op_c=None, sign=-1, threshold=1e-12): r""" Compute commutator of `op_a` and `op_b` or the symmetric double commutator of `op_a`, `op_b` and `op_c`. @@ -401,6 +401,7 @@ def commutator(op_a, op_b, op_c=None, threshold=1e-12): op_a (WeightedPauliOperator): operator a op_b (WeightedPauliOperator): operator b op_c (Optional(WeightedPauliOperator)): operator c + sign (int): -1 is anti-commute, 1 is commute threshold (float): the truncation threshold Returns: @@ -413,7 +414,7 @@ def commutator(op_a, op_b, op_c=None, threshold=1e-12): op_ba = op_b * op_a if op_c is None: - res = op_ab - op_ba + res = op_ab + sign * op_ba else: op_ac = op_a * op_c op_ca = op_c * op_a @@ -427,7 +428,7 @@ def commutator(op_a, op_b, op_c=None, threshold=1e-12): tmp = (op_bac + op_cab + op_acb + op_bca) tmp = 0.5 * tmp - res = op_abc + op_cba - tmp + res = op_abc + op_cba + sign * tmp res.simplify() res.chop(threshold) diff --git a/qiskit/chemistry/algorithms/excited_states_solvers/__init__.py b/qiskit/chemistry/algorithms/excited_states_solvers/__init__.py new file mode 100644 index 0000000000..e248adb377 --- /dev/null +++ b/qiskit/chemistry/algorithms/excited_states_solvers/__init__.py @@ -0,0 +1,25 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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. + +"""Excited states calculation algorithms.""" + +from .excited_states_solver import ExcitedStatesSolver +from .qeom import QEOM +from .eigensolver_factories import EigensolverFactory, NumPyEigensolverFactory +from .excited_states_eigensolver import ExcitedStatesEigensolver + +__all__ = ['ExcitedStatesSolver', + 'ExcitedStatesEigensolver', + 'EigensolverFactory', + 'NumPyEigensolverFactory', + 'QEOM' + ] diff --git a/qiskit/chemistry/algorithms/excited_states_solvers/eigensolver_factories/__init__.py b/qiskit/chemistry/algorithms/excited_states_solvers/eigensolver_factories/__init__.py new file mode 100644 index 0000000000..693d826bf5 --- /dev/null +++ b/qiskit/chemistry/algorithms/excited_states_solvers/eigensolver_factories/__init__.py @@ -0,0 +1,21 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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. + +"""Factories that create an eigensolver based on a qubit transformation.""" + + +from .eigensolver_factory import EigensolverFactory +from .numpy_eigensolver_factory import NumPyEigensolverFactory + +__all__ = ['EigensolverFactory', + 'NumPyEigensolverFactory' + ] diff --git a/qiskit/chemistry/algorithms/excited_states_solvers/eigensolver_factories/eigensolver_factory.py b/qiskit/chemistry/algorithms/excited_states_solvers/eigensolver_factories/eigensolver_factory.py new file mode 100644 index 0000000000..20f07dc4f9 --- /dev/null +++ b/qiskit/chemistry/algorithms/excited_states_solvers/eigensolver_factories/eigensolver_factory.py @@ -0,0 +1,34 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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 eigensolver factory for excited states calculation algorithms.""" + +from abc import ABC, abstractmethod +from qiskit.aqua.algorithms import Eigensolver +from qiskit.chemistry.transformations import Transformation + + +class EigensolverFactory(ABC): + """A factory to construct a eigensolver based on a qubit operator transformation.""" + + @abstractmethod + def get_solver(self, transformation: Transformation) -> Eigensolver: + """Returns a eigensolver, based on the qubit operator transformation. + + Args: + transformation: The qubit operator transformation. + + Returns: + An eigensolver suitable to compute the excited states of the molecule transformed + by ``transformation``. + """ + raise NotImplementedError diff --git a/qiskit/chemistry/algorithms/excited_states_solvers/eigensolver_factories/numpy_eigensolver_factory.py b/qiskit/chemistry/algorithms/excited_states_solvers/eigensolver_factories/numpy_eigensolver_factory.py new file mode 100644 index 0000000000..a31efd7295 --- /dev/null +++ b/qiskit/chemistry/algorithms/excited_states_solvers/eigensolver_factories/numpy_eigensolver_factory.py @@ -0,0 +1,97 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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 numpy eigensolver factory for ground+excited states calculation algorithms.""" + +from typing import Optional, Union, List, Callable +import numpy as np + +from qiskit.aqua.algorithms import Eigensolver, NumPyEigensolver +from qiskit.chemistry.transformations import FermionicTransformation +from qiskit.aqua.utils.validation import validate_min + +from .eigensolver_factory import EigensolverFactory + + +class NumPyEigensolverFactory(EigensolverFactory): + """A factory to construct a NumPyEigensolver.""" + + def __init__(self, + filter_criterion: Callable[[Union[List, np.ndarray], float, Optional[List[float]]], + bool] = None, + k: int = 100, + use_default_filter_criterion: bool = False) -> None: + """ + Args: + filter_criterion: callable that allows to filter eigenvalues/eigenstates. The minimum + eigensolver is only searching over feasible states and returns an eigenstate that + has the smallest eigenvalue among feasible states. The callable has the signature + `filter(eigenstate, eigenvalue, aux_values)` and must return a boolean to indicate + whether to consider this value or not. If there is no + feasible element, the result can even be empty. + k: How many eigenvalues are to be computed, has a min. value of 1. + use_default_filter_criterion: whether to use the transformation's default filter + criterion if ``filter_criterion`` is ``None``. + """ + self._filter_criterion = filter_criterion + self._k = k # pylint:disable=invalid-name + self._use_default_filter_criterion = use_default_filter_criterion + + @property + def filter_criterion(self) -> Callable[[Union[List, np.ndarray], float, Optional[List[float]]], + bool]: + """ returns filter criterion """ + return self._filter_criterion + + @filter_criterion.setter + def filter_criterion(self, value: Callable[[Union[List, np.ndarray], float, + Optional[List[float]]], bool]) -> None: + """ sets filter criterion """ + self._filter_criterion = value + + @property + def k(self) -> int: + """ returns k (number of eigenvalues requested) """ + return self._k + + @k.setter + def k(self, k: int) -> None: + """ set k (number of eigenvalues requested) """ + validate_min('k', k, 1) + self._k = k + + @property + def use_default_filter_criterion(self) -> bool: + """ returns whether to use the default filter criterion """ + return self._use_default_filter_criterion + + @use_default_filter_criterion.setter + def use_default_filter_criterion(self, value: bool) -> None: + """ sets whether to use the default filter criterion """ + self._use_default_filter_criterion = value + + def get_solver(self, transformation: FermionicTransformation) -> Eigensolver: + """Returns a NumPyEigensolver with the desired filter + + Args: + transformation: a fermionic qubit operator transformation. + + Returns: + A NumPyEigensolver suitable to compute the ground state of the molecule + transformed by ``transformation``. + """ + filter_criterion = self._filter_criterion + if not filter_criterion and self._use_default_filter_criterion: + filter_criterion = transformation.get_default_filter_criterion() + + npe = NumPyEigensolver(filter_criterion=filter_criterion, k=self.k) + return npe diff --git a/qiskit/chemistry/algorithms/excited_states_solvers/excited_states_eigensolver.py b/qiskit/chemistry/algorithms/excited_states_solvers/excited_states_eigensolver.py new file mode 100644 index 0000000000..d74c9c6264 --- /dev/null +++ b/qiskit/chemistry/algorithms/excited_states_solvers/excited_states_eigensolver.py @@ -0,0 +1,115 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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 an Eigensolver algorithm""" + +import logging + +from typing import List, Union, Optional, Any + +from qiskit.aqua.algorithms import Eigensolver +from qiskit.aqua.operators import WeightedPauliOperator +from qiskit.chemistry import FermionicOperator +from qiskit.chemistry.drivers import BaseDriver +from qiskit.chemistry.results import EigenstateResult +from qiskit.chemistry.transformations import Transformation +from qiskit.chemistry.results import ElectronicStructureResult, VibronicStructureResult + +from .excited_states_solver import ExcitedStatesSolver +from .eigensolver_factories import EigensolverFactory + +logger = logging.getLogger(__name__) + + +class ExcitedStatesEigensolver(ExcitedStatesSolver): + """The calculation of excited states via an Eigensolver algorithm""" + + def __init__(self, transformation: Transformation, + solver: Union[Eigensolver, EigensolverFactory]) -> None: + """ + + Args: + transformation: Qubit Operator Transformation + solver: Minimum Eigensolver or MESFactory object. + """ + self._transformation = transformation + self._solver = solver + + @property + def solver(self) -> Union[Eigensolver, EigensolverFactory]: + """Returns the minimum eigensolver or factory.""" + return self._solver + + @solver.setter + def solver(self, solver: Union[Eigensolver, EigensolverFactory]) -> None: + """Sets the minimum eigensolver or factory.""" + self._solver = solver + + @property + def transformation(self) -> Transformation: + """Returns the transformation used to obtain a qubit operator from the molecule.""" + return self._transformation + + @transformation.setter + def transformation(self, transformation: Transformation) -> None: + """Sets the transformation used to obtain a qubit operator from the molecule.""" + self._transformation = transformation + + def solve(self, driver: BaseDriver, + aux_operators: Optional[List[Any]] = None + ) -> Union[ElectronicStructureResult, VibronicStructureResult]: + """Compute Ground and Excited States properties. + + Args: + driver: a chemistry driver object which defines the chemical problem that is to be + solved by this calculation. + aux_operators: Additional auxiliary operators to evaluate. Must be of type + ``FermionicOperator`` if the qubit transformation is fermionic and of type + ``BosonicOperator`` it is bosonic. + + Raises: + NotImplementedError: If an operator in ``aux_operators`` is not of type + ``FermionicOperator``. + + Returns: + An eigenstate result. Depending on the transformation this can be an electronic + structure or bosonic result. + """ + if aux_operators is not None: + if any(not isinstance(op, (WeightedPauliOperator, FermionicOperator)) + for op in aux_operators): + raise NotImplementedError('Currently only fermionic problems are supported.') + + # get the operator and auxiliary operators, and transform the provided auxiliary operators + # note that ``aux_operators`` contains not only the transformed ``aux_operators`` passed + # by the user but also additional ones from the transformation + operator, aux_operators = self.transformation.transform(driver, aux_operators) + + if isinstance(self._solver, EigensolverFactory): + # this must be called after transformation.transform + solver = self._solver.get_solver(self.transformation) + else: + solver = self._solver + + # if the eigensolver does not support auxiliary operators, reset them + if not solver.supports_aux_operators(): + aux_operators = None + + raw_es_result = solver.compute_eigenvalues(operator, aux_operators) + + eigenstate_result = EigenstateResult() + eigenstate_result.raw_result = raw_es_result + eigenstate_result.eigenenergies = raw_es_result.eigenvalues + eigenstate_result.eigenstates = raw_es_result.eigenstates + eigenstate_result.aux_operator_eigenvalues = raw_es_result.aux_operator_eigenvalues + result = self.transformation.interpret(eigenstate_result) + return result diff --git a/qiskit/chemistry/algorithms/excited_states_solvers/excited_states_solver.py b/qiskit/chemistry/algorithms/excited_states_solvers/excited_states_solver.py new file mode 100644 index 0000000000..0e4a08270f --- /dev/null +++ b/qiskit/chemistry/algorithms/excited_states_solvers/excited_states_solver.py @@ -0,0 +1,41 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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 excited states calculation interface """ + +from abc import ABC, abstractmethod +from typing import List, Optional, Union + +from qiskit.chemistry.drivers import BaseDriver +from qiskit.chemistry import FermionicOperator, BosonicOperator +from qiskit.chemistry.results import ElectronicStructureResult, VibronicStructureResult + + +class ExcitedStatesSolver(ABC): + """The excited states calculation interface""" + + @abstractmethod + def solve(self, driver: BaseDriver, + aux_operators: Optional[Union[List[FermionicOperator], + List[BosonicOperator]]] = None + ) -> Union[ElectronicStructureResult, VibronicStructureResult]: + """Compute the excited states energies of the molecule that was supplied via the driver. + Args: + driver: a chemistry driver object which defines the chemical problem that is to be + solved by this calculation. + aux_operators: Additional auxiliary operators to evaluate. Must be of type + ``FermionicOperator`` if the qubit transformation is fermionic and of type + ``BosonicOperator`` it is bosonic. + Returns: + an eigenstate result + """ + raise NotImplementedError() diff --git a/qiskit/chemistry/algorithms/excited_states_solvers/qeom.py b/qiskit/chemistry/algorithms/excited_states_solvers/qeom.py new file mode 100644 index 0000000000..c7f35914e1 --- /dev/null +++ b/qiskit/chemistry/algorithms/excited_states_solvers/qeom.py @@ -0,0 +1,523 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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 typing import List, Union, Optional, Tuple, Dict, cast +import itertools +import logging +import sys +import numpy as np +from scipy import linalg + +from qiskit.tools import parallel_map +from qiskit.tools.events import TextProgressBar +from qiskit.aqua import aqua_globals +from qiskit.aqua.algorithms import AlgorithmResult +from qiskit.aqua.operators import Z2Symmetries, commutator, WeightedPauliOperator +from qiskit.chemistry import FermionicOperator, BosonicOperator +from qiskit.chemistry.drivers import BaseDriver +from qiskit.chemistry.results import (ElectronicStructureResult, VibronicStructureResult, + EigenstateResult) + +from .excited_states_solver import ExcitedStatesSolver +from ..ground_state_solvers import GroundStateSolver + +logger = logging.getLogger(__name__) + + +class QEOM(ExcitedStatesSolver): + """The calculation of excited states via the qEOM algorithm""" + + def __init__(self, ground_state_solver: GroundStateSolver, + excitations: Union[str, List[List[int]]] = 'sd') -> None: + """ + Args: + ground_state_solver: a GroundStateSolver object. The qEOM algorithm + will use this ground state to compute the EOM matrix elements + excitations: The excitations to be included in the eom pseudo-eigenvalue problem. + If a string ('s', 'd' or 'sd') then all excitations of the given type will be used. + Otherwise a list of custom excitations can directly be provided. + """ + self._gsc = ground_state_solver + self.excitations = excitations + + @property + def excitations(self) -> Union[str, List[List[int]]]: + """Returns the excitations to be included in the eom pseudo-eigenvalue problem.""" + return self._excitations + + @excitations.setter + def excitations(self, excitations: Union[str, List[List[int]]]) -> None: + """The excitations to be included in the eom pseudo-eigenvalue problem. If a string then + all excitations of given type will be used. Otherwise a list of custom excitations can + directly be provided.""" + if isinstance(excitations, str): + if excitations not in ['s', 'd', 'sd']: + raise ValueError('Excitation type must be s (singles), d (doubles) or sd ' + '(singles and doubles)') + self._excitations = excitations + + def solve(self, driver: BaseDriver, + aux_operators: Optional[Union[List[FermionicOperator], + List[BosonicOperator]]] = None + ) -> Union[ElectronicStructureResult, VibronicStructureResult]: + """Run the excited-states calculation. + + Construct and solves the EOM pseudo-eigenvalue problem to obtain the excitation energies + and the excitation operators expansion coefficients. + + Args: + driver: a chemistry driver object which defines the chemical problem that is to be + solved by this calculation. + aux_operators: Additional auxiliary operators to evaluate. Must be of type + ``FermionicOperator`` if the qubit transformation is fermionic and of type + ``BosonicOperator`` it is bosonic. + + Returns: + The excited states result. In case of a fermionic problem a + ``ElectronicStructureResult`` is returned and in the bosonic case a + ``VibronicStructureResult``. + """ + + if aux_operators is not None: + logger.warning("With qEOM the auxiliary operators can currently only be " + "evaluated on the ground state.") + + # 1. Run ground state calculation + groundstate_result = self._gsc.solve(driver, aux_operators) + + # 2. Prepare the excitation operators + matrix_operators_dict, size = self._prepare_matrix_operators() + + # 3. Evaluate eom operators + measurement_results = self._gsc.evaluate_operators( + groundstate_result.raw_result['eigenstate'], + matrix_operators_dict) + measurement_results = cast(Dict[str, List[float]], measurement_results) + + # 4. Post-process ground_state_result to construct eom matrices + m_mat, v_mat, q_mat, w_mat, m_mat_std, v_mat_std, q_mat_std, w_mat_std = \ + self._build_eom_matrices(measurement_results, size) + + # 5. solve pseudo-eigenvalue problem + energy_gaps, expansion_coefs = self._compute_excitation_energies(m_mat, v_mat, q_mat, w_mat) + + 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.m_matrix = m_mat + qeom_result.v_matrix = v_mat + qeom_result.q_matrix = q_mat + qeom_result.w_matrix = w_mat + qeom_result.m_matrix_std = m_mat_std + qeom_result.v_matrix_std = v_mat_std + qeom_result.q_matrix_std = q_mat_std + qeom_result.w_matrix_std = w_mat_std + + eigenstate_result = EigenstateResult() + eigenstate_result.eigenstates = groundstate_result.eigenstates + eigenstate_result.aux_operator_eigenvalues = groundstate_result.aux_operator_eigenvalues + eigenstate_result.raw_result = qeom_result + + eigenstate_result.eigenenergies = np.append(groundstate_result.eigenenergies, + np.asarray([groundstate_result.eigenenergies[0] + + gap for gap in energy_gaps])) + + result = self._gsc.transformation.interpret(eigenstate_result) + + return result + + def _prepare_matrix_operators(self) -> Tuple[dict, int]: + """Construct the excitation operators for each matrix element. + + Returns: + a dictionary of all matrix elements operators and the number of excitations + (or the size of the qEOM pseudo-eigenvalue problem) + """ + data = self._gsc.transformation.build_hopping_operators(self._excitations) + hopping_operators, type_of_commutativities, excitation_indices = data + + size = int(len(list(excitation_indices.keys()))/2) + + eom_matrix_operators = self._build_all_commutators( + hopping_operators, type_of_commutativities, size) + + return eom_matrix_operators, size + + def _build_all_commutators(self, hopping_operators: dict, type_of_commutativities: dict, + size: int) -> dict: + """Building all commutators for Q, W, M, V matrices. + + Args: + hopping_operators: all hopping operators based on excitations_list, + key is the string of single/double excitation; + value is corresponding operator. + type_of_commutativities: if tapering is used, it records the commutativities of + hopping operators with the + Z2 symmetries found in the original operator. + size: the number of excitations (size of the qEOM pseudo-eigenvalue problem) + + Returns: + a dictionary that contains the operators for each matrix element + """ + + all_matrix_operators = {} + + mus, nus = np.triu_indices(size) + + def _build_one_sector(available_hopping_ops, untapered_op, z2_symmetries, sign): + + to_be_computed_list = [] + for idx, _ in enumerate(mus): + m_u = mus[idx] + n_u = nus[idx] + left_op = available_hopping_ops.get('E_{}'.format(m_u)) + right_op_1 = available_hopping_ops.get('E_{}'.format(n_u)) + right_op_2 = available_hopping_ops.get('Edag_{}'.format(n_u)) + to_be_computed_list.append((m_u, n_u, left_op, right_op_1, right_op_2)) + + if logger.isEnabledFor(logging.INFO): + logger.info("Building all commutators:") + TextProgressBar(sys.stderr) + results = parallel_map(self._build_commutator_routine, + to_be_computed_list, + task_args=(untapered_op, z2_symmetries, sign), + num_processes=aqua_globals.num_processes) + for result in results: + m_u, n_u, q_mat_op, w_mat_op, m_mat_op, v_mat_op = result + + if q_mat_op is not None: + all_matrix_operators['q_{}_{}'.format(m_u, n_u)] = q_mat_op + if w_mat_op is not None: + all_matrix_operators['w_{}_{}'.format(m_u, n_u)] = w_mat_op + if m_mat_op is not None: + all_matrix_operators['m_{}_{}'.format(m_u, n_u)] = m_mat_op + if v_mat_op is not None: + all_matrix_operators['v_{}_{}'.format(m_u, n_u)] = v_mat_op + + try: + # The next step only works in the case of the FermionicTransformation. Thus, it is done + # in a try-except block. However, mypy doesn't detect this and thus we ignore it. + z2_symmetries = self._gsc.transformation.molecule_info['z2_symmetries'] # type: ignore + except AttributeError: + z2_symmetries = Z2Symmetries([], [], []) + + if not z2_symmetries.is_empty(): + combinations = itertools.product([1, -1], repeat=len(z2_symmetries.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(): + value = np.asarray(value) + if np.all(value == targeted_sector): + available_hopping_ops[key] = hopping_operators[key] + # untapered_qubit_op is a WeightedPauliOperator and should not be exposed. + _build_one_sector(available_hopping_ops, + self._gsc.transformation.untapered_qubit_op, # type: ignore + z2_symmetries, + self._gsc.transformation.commutation_rule) + + else: + # untapered_qubit_op is a WeightedPauliOperator and should not be exposed. + _build_one_sector(hopping_operators, + self._gsc.transformation.untapered_qubit_op, # type: ignore + z2_symmetries, + self._gsc.transformation.commutation_rule) + + return all_matrix_operators + + @staticmethod + def _build_commutator_routine(params: List, operator: WeightedPauliOperator, + z2_symmetries: Z2Symmetries, sign: int + ) -> Tuple[int, int, WeightedPauliOperator, WeightedPauliOperator, + WeightedPauliOperator, WeightedPauliOperator]: + """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 + z2_symmetries: z2_symmetries in case of tapering + sign: commute or anticommute + + Returns: + The indices of the matrix element and the corresponding qubit + operator for each of the EOM matrices + """ + m_u, n_u, left_op, right_op_1, right_op_2 = params + if left_op is None: + q_mat_op = None + w_mat_op = None + m_mat_op = None + v_mat_op = None + else: + if 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: + q_mat_op = commutator(left_op, operator, right_op_1, sign=sign) + w_mat_op = commutator(left_op, right_op_1, sign=sign) + q_mat_op = None if q_mat_op.is_empty() else q_mat_op + w_mat_op = None if w_mat_op.is_empty() else w_mat_op + else: + q_mat_op = None + w_mat_op = None + + if right_op_2 is not None: + m_mat_op = commutator(left_op, operator, right_op_2, sign=sign) + v_mat_op = commutator(left_op, right_op_2, sign=sign) + m_mat_op = None if m_mat_op.is_empty() else m_mat_op + v_mat_op = None if v_mat_op.is_empty() else v_mat_op + else: + m_mat_op = None + v_mat_op = None + + if not z2_symmetries.is_empty(): + if q_mat_op is not None and not q_mat_op.is_empty(): + q_mat_op = z2_symmetries.taper(q_mat_op) + if w_mat_op is not None and not w_mat_op.is_empty(): + w_mat_op = z2_symmetries.taper(w_mat_op) + if m_mat_op is not None and not m_mat_op.is_empty(): + m_mat_op = z2_symmetries.taper(m_mat_op) + if v_mat_op is not None and not v_mat_op.is_empty(): + v_mat_op = z2_symmetries.taper(v_mat_op) + + return m_u, n_u, q_mat_op, w_mat_op, m_mat_op, v_mat_op + + def _build_eom_matrices(self, gs_results: Dict[str, List[float]], size: int + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, + float, float, float, float]: + """Constructs the M, V, Q and W matrices from the results on the ground state + + Args: + gs_results: a ground state result object + size: size of eigenvalue problem + + Returns: + the 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. + + # evaluate results + for idx, _ in enumerate(mus): + m_u = mus[idx] + n_u = nus[idx] + + q_mat[m_u][n_u] = gs_results['q_{}_{}'.format(m_u, n_u)][0] if gs_results.get( + 'q_{}_{}'.format(m_u, n_u)) is not None else q_mat[m_u][n_u] + w_mat[m_u][n_u] = gs_results['w_{}_{}'.format(m_u, n_u)][0] if gs_results.get( + 'w_{}_{}'.format(m_u, n_u)) is not None else w_mat[m_u][n_u] + m_mat[m_u][n_u] = gs_results['m_{}_{}'.format(m_u, n_u)][0] if gs_results.get( + 'm_{}_{}'.format(m_u, n_u)) is not None else m_mat[m_u][n_u] + v_mat[m_u][n_u] = gs_results['v_{}_{}'.format(m_u, n_u)][0] if gs_results.get( + 'v_{}_{}'.format(m_u, n_u)) is not None else v_mat[m_u][n_u] + + q_mat_std += gs_results['q_{}_{}_std'.format(m_u, n_u)][0] if gs_results.get( + 'q_{}_{}_std'.format(m_u, n_u)) is not None else 0 + w_mat_std += gs_results['w_{}_{}_std'.format(m_u, n_u)][0] if gs_results.get( + 'w_{}_{}_std'.format(m_u, n_u)) is not None else 0 + m_mat_std += gs_results['m_{}_{}_std'.format(m_u, n_u)][0] if gs_results.get( + 'm_{}_{}_std'.format(m_u, n_u)) is not None else 0 + v_mat_std += gs_results['v_{}_{}_std'.format(m_u, n_u)][0] if gs_results.get( + 'v_{}_{}_std'.format(m_u, n_u)) is not None else 0 + + # these matrices are numpy arrays and therefore have the ``shape`` attribute + # pylint: disable=unsubscriptable-object + 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 - np.identity(m_mat.shape[0]) * m_mat + v_mat = v_mat + v_mat.T - 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) + + return m_mat, v_mat, q_mat, w_mat, m_mat_std, v_mat_std, q_mat_std, w_mat_std + + @staticmethod + def _compute_excitation_energies(m_mat: np.ndarray, v_mat: np.ndarray, q_mat: np.ndarray, + w_mat: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: + """Diagonalizing M, V, Q, W matrices for excitation energies. + + Args: + m_mat : M matrices + v_mat : V matrices + q_mat : Q matrices + w_mat : W matrices + + 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...') + a_mat = np.bmat([[m_mat, q_mat], [q_mat.T.conj(), m_mat.T.conj()]]) + b_mat = np.bmat([[v_mat, w_mat], [-w_mat.T.conj(), -v_mat.T.conj()]]) + # pylint: disable=too-many-function-args + res = linalg.eig(a_mat, b_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. So to get the values we want we + # sort the real parts and then take the upper half of the sorted values. + # Since we may now have + # small values (positive or negative) take the absolute and then threshold zero. + logger.debug('... %s', res[0]) + w = np.sort(np.real(res[0])) + logger.debug('Sorted real parts %s', w) + w = np.abs(w[len(w) // 2:]) + w[w < 1e-06] = 0 + excitation_energies_gap = w + + return excitation_energies_gap, res[1] + + +class QEOMResult(AlgorithmResult): + """The results class for the QEOM algorithm.""" + + @property + def ground_state_raw_result(self): + """ returns ground state raw result """ + return self.get('ground_state_raw_result') + + @ground_state_raw_result.setter + def ground_state_raw_result(self, value) -> None: + """ sets ground state raw result """ + self.data['ground_state_raw_result'] = value + + @property + def excitation_energies(self) -> np.ndarray: + """ returns the excitation energies (energy gaps) """ + return self.get('excitation_energies') + + @excitation_energies.setter + def excitation_energies(self, value: np.ndarray) -> None: + """ sets the excitation energies (energy gaps) """ + self.data['excitation_energies'] = value + + @property + def expansion_coefficients(self) -> np.ndarray: + """ returns the X and Y expansion coefficients """ + return self.get('expansion_coefficients') + + @expansion_coefficients.setter + def expansion_coefficients(self, value: np.ndarray) -> None: + """ sets the X and Y expansion coefficients """ + self.data['expansion_coefficients'] = value + + @property + def m_matrix(self) -> np.ndarray: + """ returns the M matrix """ + return self.get('m_matrix') + + @m_matrix.setter + def m_matrix(self, value: np.ndarray) -> None: + """ sets the M matrix """ + self.data['m_matrix'] = value + + @property + def v_matrix(self) -> np.ndarray: + """ returns the V matrix """ + return self.get('v_matrix') + + @v_matrix.setter + def v_matrix(self, value: np.ndarray) -> None: + """ sets the V matrix """ + self.data['v_matrix'] = value + + @property + def q_matrix(self) -> np.ndarray: + """ returns the Q matrix """ + return self.get('q_matrix') + + @q_matrix.setter + def q_matrix(self, value: np.ndarray) -> None: + """ sets the Q matrix """ + self.data['q_matrix'] = value + + @property + def w_matrix(self) -> np.ndarray: + """ returns the W matrix """ + return self.get('w_matrix') + + @w_matrix.setter + def w_matrix(self, value: np.ndarray) -> None: + """ sets the W matrix """ + self.data['w_matrix'] = value + + @property + def m_matrix_std(self) -> float: + """ returns the M matrix standard deviation """ + return self.get('m_matrix_std') + + @m_matrix_std.setter + def m_matrix_std(self, value: float) -> None: + """ sets the M matrix standard deviation """ + self.data['m_matrix_std'] = value + + @property + def v_matrix_std(self) -> float: + """ returns the V matrix standard deviation """ + return self.get('v_matrix_std') + + @v_matrix_std.setter + def v_matrix_std(self, value: float) -> None: + """ sets the V matrix standard deviation """ + self.data['v_matrix_std'] = value + + @property + def q_matrix_std(self) -> float: + """ returns the Q matrix standard deviation """ + return self.get('q_matrix_std') + + @q_matrix_std.setter + def q_matrix_std(self, value: float) -> None: + """ sets the Q matrix standard deviation """ + self.data['q_matrix_std'] = value + + @property + def w_matrix_std(self) -> float: + """ returns the W matrix standard deviation """ + return self.get('w_matrix_std') + + @w_matrix_std.setter + def w_matrix_std(self, value: float) -> None: + """ sets the W matrix standard deviation """ + self.data['w_matrix_std'] = value diff --git a/qiskit/chemistry/algorithms/ground_state_solvers/adapt_vqe.py b/qiskit/chemistry/algorithms/ground_state_solvers/adapt_vqe.py index cc1de106e6..9414a5ba4d 100644 --- a/qiskit/chemistry/algorithms/ground_state_solvers/adapt_vqe.py +++ b/qiskit/chemistry/algorithms/ground_state_solvers/adapt_vqe.py @@ -240,7 +240,7 @@ def solve(self, result.final_max_gradient = max_grad[0] result.finishing_criterion = finishing_criterion - logger.info('The final energy is: %s', str(result.computed_electronic_energy)) + logger.info('The final energy is: %s', str(result.computed_energies[0])) return result diff --git a/qiskit/chemistry/algorithms/ground_state_solvers/ground_state_eigensolver.py b/qiskit/chemistry/algorithms/ground_state_solvers/ground_state_eigensolver.py index 8a8aac131a..b489ca4d18 100644 --- a/qiskit/chemistry/algorithms/ground_state_solvers/ground_state_eigensolver.py +++ b/qiskit/chemistry/algorithms/ground_state_solvers/ground_state_eigensolver.py @@ -115,7 +115,7 @@ def evaluate_operators(self, QuantumCircuit, Instruction, OperatorBase], operators: Union[WeightedPauliOperator, OperatorBase, list, dict] - ) -> Union[float, List[float], Dict[str, float]]: + ) -> Union[float, List[float], Dict[str, List[float]]]: """Evaluates additional operators at the given state. Args: diff --git a/qiskit/chemistry/algorithms/ground_state_solvers/ground_state_solver.py b/qiskit/chemistry/algorithms/ground_state_solvers/ground_state_solver.py index 612878567d..a67f6664b8 100644 --- a/qiskit/chemistry/algorithms/ground_state_solvers/ground_state_solver.py +++ b/qiskit/chemistry/algorithms/ground_state_solvers/ground_state_solver.py @@ -88,7 +88,7 @@ def evaluate_operators(self, QuantumCircuit, Instruction, OperatorBase], operators: Union[WeightedPauliOperator, OperatorBase, list, dict] - ) -> Union[float, List[float], Dict[str, float]]: + ) -> Union[float, List[float], Dict[str, List[float]]]: """Evaluates additional operators at the given state. Args: diff --git a/qiskit/chemistry/algorithms/pes_samplers/bopes_sampler.py b/qiskit/chemistry/algorithms/pes_samplers/bopes_sampler.py index 9ce821b4f1..229430999d 100644 --- a/qiskit/chemistry/algorithms/pes_samplers/bopes_sampler.py +++ b/qiskit/chemistry/algorithms/pes_samplers/bopes_sampler.py @@ -112,7 +112,7 @@ def sample(self, driver: BaseDriver, points: List[float]) -> BOPESSamplerResult: self._points = list(self._raw_results.keys()) self._energies = [] for key in self._raw_results: - energy = self._raw_results[key].computed_electronic_energy + \ + energy = self._raw_results[key].computed_energies[0] + \ self._raw_results[key].nuclear_repulsion_energy self._energies.append(energy) diff --git a/qiskit/chemistry/results/electronic_structure_result.py b/qiskit/chemistry/results/electronic_structure_result.py index 06b92fd748..6d01188052 100644 --- a/qiskit/chemistry/results/electronic_structure_result.py +++ b/qiskit/chemistry/results/electronic_structure_result.py @@ -12,7 +12,7 @@ """The electronic structure result.""" -from typing import List, Optional, Tuple, cast +from typing import List, Optional, Tuple, cast, Union import logging import numpy as np @@ -67,31 +67,31 @@ def nuclear_dipole_moment(self, value: DipoleTuple) -> None: # construct the circuit of the GS from here (if the algorithm supports this) @property - def energy(self) -> Optional[float]: + def total_energies(self) -> np.ndarray: """ Returns ground state energy if nuclear_repulsion_energy is available from driver """ - # TODO the fact that this property is computed on the fly breaks the `.combine()` - # functionality - nre = self.nuclear_repulsion_energy - return self.electronic_energy + nre if nre is not None else None + nre = self.nuclear_repulsion_energy if self.nuclear_repulsion_energy is not None else 0 + # Adding float to np.ndarray adds it to each entry + return self.electronic_energies + nre @property - def electronic_energy(self) -> float: + def electronic_energies(self) -> np.ndarray: """ Returns electronic part of ground state energy """ # TODO the fact that this property is computed on the fly breaks the `.combine()` # functionality - return (self.computed_electronic_energy + # Adding float to np.ndarray adds it to each entry + return (self.computed_energies + self.ph_extracted_energy + self.frozen_extracted_energy) @property - def computed_electronic_energy(self) -> float: + def computed_energies(self) -> np.ndarray: """ Returns computed electronic part of ground state energy """ - return self.get('computed_electronic_energy') + return self.get('computed_energies') - @computed_electronic_energy.setter - def computed_electronic_energy(self, value: float) -> None: + @computed_energies.setter + def computed_energies(self, value: np.ndarray) -> None: """ Sets computed electronic part of ground state energy """ - self.data['computed_electronic_energy'] = value + self.data['computed_energies'] = value @property def ph_extracted_energy(self) -> float: @@ -131,73 +131,84 @@ def reverse_dipole_sign(self, value: bool) -> None: self.data['reverse_dipole_sign'] = value @property - def total_dipole_moment(self) -> Optional[float]: + def total_dipole_moment(self) -> Optional[List[float]]: """ Returns total dipole of moment """ if self.dipole_moment is None: return None # No dipole at all - if np.any(np.equal(list(self.dipole_moment), None)): - return None # One or more components in the dipole is None - return np.sqrt(np.sum(np.power(list(self.dipole_moment), 2))) + tdm: List[float] = [] + for dip in self.dipole_moment: + if np.any(np.equal(list(dip), None)): + tdm.append(None) # One or more components in the dipole is None + else: + tdm.append(np.sqrt(np.sum(np.power(list(dip), 2)))) + return tdm @property - def total_dipole_moment_in_debye(self) -> Optional[float]: + def total_dipole_moment_in_debye(self) -> Optional[List[float]]: """ Returns total dipole of moment in Debye """ tdm = self.total_dipole_moment - return tdm / QMolecule.DEBYE if tdm is not None else None + if tdm is None: + return None + return [dip / QMolecule.DEBYE for dip in tdm] @property - def dipole_moment(self) -> Optional[DipoleTuple]: + def dipole_moment(self) -> Optional[List[DipoleTuple]]: """ Returns dipole moment """ edm = self.electronic_dipole_moment if self.reverse_dipole_sign: - edm = cast(DipoleTuple, tuple(-1 * x if x is not None else None for x in edm)) - return _dipole_tuple_add(edm, self.nuclear_dipole_moment) + edm = [cast(DipoleTuple, tuple(-1 * x if x is not None else None for x in dip)) + for dip in edm] + return [_dipole_tuple_add(dip, self.nuclear_dipole_moment) for dip in edm] @property - def dipole_moment_in_debye(self) -> Optional[DipoleTuple]: + def dipole_moment_in_debye(self) -> Optional[List[DipoleTuple]]: """ Returns dipole moment in Debye """ dipm = self.dipole_moment if dipm is None: return None - dipmd0 = dipm[0]/QMolecule.DEBYE if dipm[0] is not None else None - dipmd1 = dipm[1]/QMolecule.DEBYE if dipm[1] is not None else None - dipmd2 = dipm[2]/QMolecule.DEBYE if dipm[2] is not None else None - return dipmd0, dipmd1, dipmd2 + dipmd = [] + for dip in dipm: + dipmd0 = dip[0]/QMolecule.DEBYE if dip[0] is not None else None + dipmd1 = dip[1]/QMolecule.DEBYE if dip[1] is not None else None + dipmd2 = dip[2]/QMolecule.DEBYE if dip[2] is not None else None + dipmd += [(dipmd0, dipmd1, dipmd2)] + return dipmd @property - def electronic_dipole_moment(self) -> Optional[DipoleTuple]: + def electronic_dipole_moment(self) -> Optional[List[DipoleTuple]]: """ Returns electronic dipole moment """ - return _dipole_tuple_add(self.computed_dipole_moment, - _dipole_tuple_add(self.ph_extracted_dipole_moment, - self.frozen_extracted_dipole_moment)) + return [_dipole_tuple_add(comp_dip, _dipole_tuple_add(ph_dip, frozen_dip)) for + comp_dip, ph_dip, frozen_dip in zip(self.computed_dipole_moment, + self.ph_extracted_dipole_moment, + self.frozen_extracted_dipole_moment)] @property - def computed_dipole_moment(self) -> Optional[DipoleTuple]: + def computed_dipole_moment(self) -> Optional[List[DipoleTuple]]: """ Returns computed electronic part of dipole moment """ return self.get('computed_dipole_moment') @computed_dipole_moment.setter - def computed_dipole_moment(self, value: DipoleTuple) -> None: + def computed_dipole_moment(self, value: List[DipoleTuple]) -> None: """ Sets computed electronic part of dipole moment """ self.data['computed_dipole_moment'] = value @property - def ph_extracted_dipole_moment(self) -> Optional[DipoleTuple]: + def ph_extracted_dipole_moment(self) -> Optional[List[DipoleTuple]]: """ Returns particle hole extracted part of dipole moment """ return self.get('ph_extracted_dipole_moment') @ph_extracted_dipole_moment.setter - def ph_extracted_dipole_moment(self, value: DipoleTuple) -> None: + def ph_extracted_dipole_moment(self, value: List[DipoleTuple]) -> None: """ Sets particle hole extracted part of dipole moment """ self.data['ph_extracted_dipole_moment'] = value @property - def frozen_extracted_dipole_moment(self) -> Optional[DipoleTuple]: + def frozen_extracted_dipole_moment(self) -> Optional[List[DipoleTuple]]: """ Returns frozen extracted part of dipole moment """ return self.get('frozen_extracted_dipole_moment') @frozen_extracted_dipole_moment.setter - def frozen_extracted_dipole_moment(self, value: DipoleTuple) -> None: + def frozen_extracted_dipole_moment(self, value: List[DipoleTuple]) -> None: """ Sets frozen extracted part of dipole moment """ self.data['frozen_extracted_dipole_moment'] = value @@ -211,39 +222,42 @@ def has_observables(self): or self.magnetization is not None @property - def total_angular_momentum(self) -> Optional[float]: + def total_angular_momentum(self) -> Optional[List[float]]: """ Returns total angular momentum (S^2) """ return self.get('total_angular_momentum') @total_angular_momentum.setter - def total_angular_momentum(self, value: float) -> None: + def total_angular_momentum(self, value: List[float]) -> None: """ Sets total angular momentum """ self.data['total_angular_momentum'] = value @property - def spin(self) -> Optional[float]: + def spin(self) -> Optional[List[float]]: """ Returns computed spin """ if self.total_angular_momentum is None: return None - return (-1.0 + np.sqrt(1 + 4 * self.total_angular_momentum)) / 2 + spin = [] + for total_angular_momentum in self.total_angular_momentum: + spin.append((-1.0 + np.sqrt(1 + 4 * total_angular_momentum)) / 2) + return spin @property - def num_particles(self) -> Optional[float]: + def num_particles(self) -> Optional[List[float]]: """ Returns measured number of particles """ return self.get('num_particles') @num_particles.setter - def num_particles(self, value: float) -> None: + def num_particles(self, value: List[float]) -> None: """ Sets measured number of particles """ self.data['num_particles'] = value @property - def magnetization(self) -> Optional[float]: + def magnetization(self) -> Optional[List[float]]: """ Returns measured magnetization """ return self.get('magnetization') @magnetization.setter - def magnetization(self, value: float) -> None: + def magnetization(self, value: List[float]) -> None: """ Sets measured magnetization """ self.data['magnetization'] = value @@ -258,9 +272,9 @@ def formatted(self) -> List[str]: lines.append('=== GROUND STATE ENERGY ===') lines.append(' ') lines.append('* Electronic ground state energy (Hartree): {}'. - format(round(self.electronic_energy, 12))) + format(round(self.electronic_energies[0], 12))) lines.append(' - computed part: {}'. - format(round(self.computed_electronic_energy, 12))) + format(round(self.computed_energies[0], 12))) lines.append(' - frozen energy part: {}'. format(round(self.frozen_extracted_energy, 12))) lines.append(' - particle hole part: {}' @@ -269,40 +283,68 @@ def formatted(self) -> List[str]: lines.append('~ Nuclear repulsion energy (Hartree): {}'. format(round(self.nuclear_repulsion_energy, 12))) lines.append('> Total ground state energy (Hartree): {}'. - format(round(self.energy, 12))) + format(round(self.total_energies[0], 12))) + + if len(self.computed_energies) > 1: + lines.append(' ') + lines.append('=== EXCITED STATE ENERGIES ===') + lines.append(' ') + for idx, (elec_energy, total_energy) in enumerate(zip(self.electronic_energies[1:], + self.total_energies[1:])): + lines.append('{: 3d}: '.format(idx+1)) + lines.append('* Electronic excited state energy (Hartree): {}'. + format(round(elec_energy, 12))) + lines.append('> Total excited state energy (Hartree): {}'. + format(round(total_energy, 12))) + if self.has_observables(): - line = ' Measured::' - if self.num_particles is not None: - line += ' # Particles: {:.3f}'.format(self.num_particles) - if self.spin is not None: - line += ' S: {:.3f}'.format(self.spin) - if self.total_angular_momentum is not None: - line += ' S^2: {:.3f}'.format(self.total_angular_momentum) - if self.magnetization is not None: - line += ' M: {:.5f}'.format(self.magnetization) - lines.append(line) + lines.append(' ') + lines.append('=== MEASURED OBSERVABLES ===') + lines.append(' ') + for idx, (num_particles, spin, total_angular_momentum, magnetization) in enumerate(zip( + self.num_particles, self.spin, self.total_angular_momentum, + self.magnetization)): + line = '{: 3d}: '.format(idx) + if num_particles is not None: + line += ' # Particles: {:.3f}'.format(num_particles) + if spin is not None: + line += ' S: {:.3f}'.format(spin) + if total_angular_momentum is not None: + line += ' S^2: {:.3f}'.format(total_angular_momentum) + if magnetization is not None: + line += ' M: {:.3f}'.format(magnetization) + lines.append(line) if self.has_dipole(): lines.append(' ') - lines.append('=== DIPOLE MOMENT ===') + lines.append('=== DIPOLE MOMENTS ===') lines.append(' ') - lines.append('* Electronic dipole moment (a.u.): {}' - .format(_dipole_to_string(self.electronic_dipole_moment))) - lines.append(' - computed part: {}' - .format(_dipole_to_string(self.computed_dipole_moment))) - lines.append(' - frozen energy part: {}' - .format(_dipole_to_string(self.frozen_extracted_dipole_moment))) - lines.append(' - particle hole part: {}' - .format(_dipole_to_string(self.ph_extracted_dipole_moment))) if self.nuclear_dipole_moment is not None: lines.append('~ Nuclear dipole moment (a.u.): {}' .format(_dipole_to_string(self.nuclear_dipole_moment))) - lines.append('> Dipole moment (a.u.): {} Total: {}' - .format(_dipole_to_string(self.dipole_moment), - _float_to_string(self.total_dipole_moment))) - lines.append(' (debye): {} Total: {}' - .format(_dipole_to_string(self.dipole_moment_in_debye), - _float_to_string(self.total_dipole_moment_in_debye))) + lines.append(' ') + for idx, (elec_dip, comp_dip, frozen_dip, ph_dip, dip, tot_dip, dip_db, tot_dip_db) in \ + enumerate(zip( + self.electronic_dipole_moment, self.computed_dipole_moment, + self.frozen_extracted_dipole_moment, self.ph_extracted_dipole_moment, + self.dipole_moment, self.total_dipole_moment, + self.dipole_moment_in_debye, self.total_dipole_moment_in_debye)): + lines.append('{: 3d}: '.format(idx)) + lines.append(' * Electronic dipole moment (a.u.): {}' + .format(_dipole_to_string(elec_dip))) + lines.append(' - computed part: {}' + .format(_dipole_to_string(comp_dip))) + lines.append(' - frozen energy part: {}' + .format(_dipole_to_string(frozen_dip))) + lines.append(' - particle hole part: {}' + .format(_dipole_to_string(ph_dip))) + if self.nuclear_dipole_moment is not None: + lines.append(' > Dipole moment (a.u.): {} Total: {}' + .format(_dipole_to_string(dip), _float_to_string(tot_dip))) + lines.append(' (debye): {} Total: {}' + .format(_dipole_to_string(dip_db), _float_to_string(tot_dip_db))) + lines.append(' ') + return lines diff --git a/qiskit/chemistry/transformations/fermionic_transformation.py b/qiskit/chemistry/transformations/fermionic_transformation.py index 794f7b3c09..1d22cd39bd 100644 --- a/qiskit/chemistry/transformations/fermionic_transformation.py +++ b/qiskit/chemistry/transformations/fermionic_transformation.py @@ -21,12 +21,15 @@ from enum import Enum import numpy as np -from qiskit.aqua.algorithms import EigensolverResult, MinimumEigensolverResult +from qiskit.tools import parallel_map +from qiskit.aqua import AquaError, aqua_globals from qiskit.aqua.operators import Z2Symmetries, WeightedPauliOperator, OperatorBase +from qiskit.aqua.algorithms import EigensolverResult, MinimumEigensolverResult from qiskit.chemistry import QiskitChemistryError, QMolecule from qiskit.chemistry.fermionic_operator import FermionicOperator from qiskit.chemistry.drivers import BaseDriver from qiskit.chemistry.results import DipoleTuple, EigenstateResult, ElectronicStructureResult +from qiskit.chemistry.components.variational_forms import UCCSD from .transformation import Transformation from ..components.initial_states import HartreeFock @@ -96,6 +99,7 @@ def __init__(self, raise QiskitChemistryError('Invalid z2symmetry_reduction value') self._z2symmetry_reduction = z2symmetry_reduction self._has_dipole_moments = False + self._untapered_qubit_op = None # Store values that are computed by the classical logic in order # that later they may be combined with the quantum result @@ -116,6 +120,11 @@ def __init__(self, self._molecule_info: Dict[str, Any] = {} + @property + def commutation_rule(self) -> int: + """Getter of the commutation rule""" + return -1 + @property def molecule_info(self) -> Dict[str, Any]: """Getter of the molecule information.""" @@ -359,6 +368,7 @@ def _dipole_op(dipole_integrals: np.ndarray, axis: str) \ self._molecule_info['num_orbitals'] = new_nspinorbs reduction = self._two_qubit_reduction if self._qubit_mapping == 'parity' else False self._molecule_info['two_qubit_reduction'] = reduction + self._untapered_qubit_op = qubit_op z2symmetries = Z2Symmetries([], [], [], None) if self._z2symmetry_reduction is not None: @@ -369,6 +379,11 @@ def _dipole_op(dipole_integrals: np.ndarray, axis: str) \ logger.debug('Processing complete ready to run algorithm') return qubit_op, aux_ops + @property + def untapered_qubit_op(self): + """Getter for the untapered qubit operator""" + return self._untapered_qubit_op + def _process_z2symmetry_reduction(self, qubit_op: WeightedPauliOperator, aux_ops: List[WeightedPauliOperator]) -> Tuple: @@ -527,10 +542,10 @@ def interpret(self, raw_result: Union[EigenstateResult, EigensolverResult, eigenstate_result.raw_result = raw_result eigenstate_result.eigenenergies = np.asarray([raw_result.eigenvalue]) eigenstate_result.eigenstates = [raw_result.eigenstate] - eigenstate_result.aux_operator_eigenvalues = raw_result.aux_operator_eigenvalues + eigenstate_result.aux_operator_eigenvalues = [raw_result.aux_operator_eigenvalues] result = ElectronicStructureResult(eigenstate_result.data) - result.computed_electronic_energy = eigenstate_result.groundenergy + result.computed_energies = np.asarray([e.real for e in eigenstate_result.eigenenergies]) result.hartree_fock_energy = self._hf_energy result.nuclear_repulsion_energy = self._nuclear_repulsion_energy if self._nuclear_dipole_moment is not None: @@ -540,40 +555,48 @@ def interpret(self, raw_result: Union[EigenstateResult, EigensolverResult, if result.aux_operator_eigenvalues is not None: # the first three values are hardcoded to number of particles, angular momentum # and magnetization in this order - if result.aux_operator_eigenvalues[0] is not None: - result.num_particles = result.aux_operator_eigenvalues[0][0].real + result.num_particles = [] + result.total_angular_momentum = [] + result.magnetization = [] + result.computed_dipole_moment = [] + result.ph_extracted_dipole_moment = [] + result.frozen_extracted_dipole_moment = [] + if not isinstance(result.aux_operator_eigenvalues, list): + aux_operator_eigenvalues = [result.aux_operator_eigenvalues] else: - result.num_particles = None - - if result.aux_operator_eigenvalues[1] is not None: - result.total_angular_momentum = result.aux_operator_eigenvalues[1][0].real - else: - result.total_angular_momentum = None - - if result.aux_operator_eigenvalues[2] is not None: - result.magnetization = result.aux_operator_eigenvalues[2][0].real - else: - result.magnetization = None - - # the next three are hardcoded to Dipole moments, if they are set - if len(result.aux_operator_eigenvalues) >= 6 and self._has_dipole_moments: - # check if the names match - # extract dipole moment in each axis - dipole_moment = [] - for moment in result.aux_operator_eigenvalues[3:6]: - if moment is not None: - dipole_moment += [moment[0].real] - else: - dipole_moment += [None] - - result.reverse_dipole_sign = self._reverse_dipole_sign - result.computed_dipole_moment = cast(DipoleTuple, tuple(dipole_moment)) - result.ph_extracted_dipole_moment = (self._ph_x_dipole_shift, - self._ph_y_dipole_shift, - self._ph_z_dipole_shift) - result.frozen_extracted_dipole_moment = (self._x_dipole_shift, - self._y_dipole_shift, - self._z_dipole_shift) + aux_operator_eigenvalues = result.aux_operator_eigenvalues + for aux_op_eigenvalues in aux_operator_eigenvalues: + if aux_op_eigenvalues is None: + continue + if aux_op_eigenvalues[0] is not None: + result.num_particles.append(aux_op_eigenvalues[0][0].real) + + if aux_op_eigenvalues[1] is not None: + result.total_angular_momentum.append(aux_op_eigenvalues[1][0].real) + + if aux_op_eigenvalues[2] is not None: + result.magnetization.append(aux_op_eigenvalues[2][0].real) + + # the next three are hardcoded to Dipole moments, if they are set + if len(aux_op_eigenvalues) >= 6 and self._has_dipole_moments: + # check if the names match + # extract dipole moment in each axis + dipole_moment = [] + for moment in aux_op_eigenvalues[3:6]: + if moment is not None: + dipole_moment += [moment[0].real] + else: + dipole_moment += [None] + + result.reverse_dipole_sign = self._reverse_dipole_sign + result.computed_dipole_moment.append(cast(DipoleTuple, + tuple(dipole_moment))) + result.ph_extracted_dipole_moment.append( + (self._ph_x_dipole_shift, self._ph_y_dipole_shift, + self._ph_z_dipole_shift)) + result.frozen_extracted_dipole_moment.append( + (self._x_dipole_shift, self._y_dipole_shift, + self._z_dipole_shift)) return result @@ -623,3 +646,98 @@ def _map_fermionic_operator_to_qubit(fer_op: FermionicOperator, if qubit_mapping == 'parity' and two_qubit_reduction: qubit_op = Z2Symmetries.two_qubit_reduction(qubit_op, num_particles) return qubit_op + + @staticmethod + def _build_single_hopping_operator(index, num_particles, num_orbitals, qubit_mapping, + two_qubit_reduction, z2_symmetries): + + h_1 = np.zeros((num_orbitals, num_orbitals), dtype=complex) + h_2 = np.zeros((num_orbitals, num_orbitals, num_orbitals, num_orbitals), dtype=complex) + if len(index) == 2: + i, j = index + h_1[i, j] = 4.0 + elif len(index) == 4: + i, j, k, m = index + h_2[i, j, k, m] = 16.0 + fer_op = FermionicOperator(h_1, h_2) + qubit_op = fer_op.mapping(qubit_mapping) + if qubit_mapping == 'parity' and two_qubit_reduction: + qubit_op = Z2Symmetries.two_qubit_reduction(qubit_op, num_particles) + + commutativities = [] + if not z2_symmetries.is_empty(): + for symmetry in z2_symmetries.symmetries: + symmetry_op = WeightedPauliOperator(paulis=[[1.0, symmetry]]) + commuting = qubit_op.commute_with(symmetry_op) + anticommuting = qubit_op.anticommute_with(symmetry_op) + + if commuting != anticommuting: # only one of them is True + if commuting: + commutativities.append(True) + elif anticommuting: + commutativities.append(False) + else: + raise AquaError("Symmetry {} is nor commute neither anti-commute " + "to exciting operator.".format(symmetry.to_label())) + + return qubit_op, commutativities + + def build_hopping_operators(self, excitations: Union[str, List[List[int]]] = 'sd' + ) -> Tuple[Dict[str, WeightedPauliOperator], + Dict[str, List[bool]], + Dict[str, List[Any]]]: + """Builds the product of raising and lowering operators (basic excitation operators) + + Args: + excitations: The excitations to be included in the eom pseudo-eigenvalue problem. + If a string ('s', 'd' or 'sd') then all excitations of the given type will be used. + Otherwise a list of custom excitations can directly be provided. + + Returns: + A tuple containing the hopping operators, the types of commutativities and the + excitation indices. + """ + + num_alpha, num_beta = self._molecule_info['num_particles'] + num_orbitals = self._molecule_info['num_orbitals'] + + if isinstance(excitations, str): + se_list, de_list = UCCSD.compute_excitation_lists([num_alpha, num_beta], num_orbitals, + excitation_type=excitations) + excitations_list = se_list+de_list + else: + excitations_list = excitations + + size = len(excitations_list) + + # # get all to-be-processed index + # mus, nus = np.triu_indices(size) + + # build all hopping operators + hopping_operators: Dict[str, WeightedPauliOperator] = {} + type_of_commutativities: Dict[str, List[bool]] = {} + excitation_indices = {} + to_be_executed_list = [] + for idx in range(size): + to_be_executed_list += [excitations_list[idx], list(reversed(excitations_list[idx]))] + hopping_operators['E_{}'.format(idx)] = None + hopping_operators['Edag_{}'.format(idx)] = None + type_of_commutativities['E_{}'.format(idx)] = None + type_of_commutativities['Edag_{}'.format(idx)] = None + excitation_indices['E_{}'.format(idx)] = excitations_list[idx] + excitation_indices['Edag_{}'.format(idx)] = list(reversed(excitations_list[idx])) + + result = parallel_map(self._build_single_hopping_operator, + to_be_executed_list, + task_args=(num_alpha + num_beta, + num_orbitals, + self._qubit_mapping, + self._two_qubit_reduction, + self._molecule_info['z2_symmetries']), + num_processes=aqua_globals.num_processes) + + for key, res in zip(hopping_operators.keys(), result): + hopping_operators[key] = res[0] + type_of_commutativities[key] = res[1] + + return hopping_operators, type_of_commutativities, excitation_indices diff --git a/qiskit/chemistry/transformations/transformation.py b/qiskit/chemistry/transformations/transformation.py index 6348257606..305b410dfc 100644 --- a/qiskit/chemistry/transformations/transformation.py +++ b/qiskit/chemistry/transformations/transformation.py @@ -13,12 +13,12 @@ """Base class for transformation to qubit operators for chemistry problems""" from abc import ABC, abstractmethod -from typing import Tuple, List, Optional, Union, Callable +from typing import Tuple, List, Optional, Union, Callable, Dict, Any import numpy as np from qiskit.aqua.algorithms import EigensolverResult, MinimumEigensolverResult -from qiskit.aqua.operators import OperatorBase +from qiskit.aqua.operators import OperatorBase, WeightedPauliOperator from qiskit.chemistry import FermionicOperator, BosonicOperator from qiskit.chemistry.drivers import BaseDriver from qiskit.chemistry.results import EigenstateResult @@ -65,3 +65,26 @@ def interpret(self, raw_result: Union[EigenstateResult, EigensolverResult, An "interpreted" eigenstate result. """ raise NotImplementedError + + @property + @abstractmethod + def commutation_rule(self) -> int: + """Getter of the commutation rule""" + raise NotImplementedError + + @abstractmethod + def build_hopping_operators(self, excitations: Union[str, List[List[int]]] = 'sd' + ) -> Tuple[Dict[str, WeightedPauliOperator], + Dict[str, List[bool]], + Dict[str, List[Any]]]: + """Builds the product of raising and lowering operators (basic excitation operators) + + Args: + excitations: The excitations to be included in the eom pseudo-eigenvalue problem. + If a string ('s', 'd' or 'sd') then all excitations of the given type will be used. + Otherwise a list of custom excitations can directly be provided. + + Returns: + + """ + raise NotImplementedError diff --git a/releasenotes/notes/excited-states-interface-a7ac8ae621b35daa.yaml b/releasenotes/notes/excited-states-interface-a7ac8ae621b35daa.yaml new file mode 100644 index 0000000000..c38d36d441 --- /dev/null +++ b/releasenotes/notes/excited-states-interface-a7ac8ae621b35daa.yaml @@ -0,0 +1,20 @@ +--- +prelude: > + This release introduces an interface for excited states calculations. + It is now easier for the user to create a general excited states calculation. + This calculation is based on a Driver which provides the relevant information + about the molecule, a Transformation which provides the information about the + mapping of the problem into a qubit Hamiltonian, and finally a Solver. + The Solver is the specific way which the excited states calculation is done + (the algorithm). This structure follows the one of the ground state calculations. + The results are modified to take lists of expectation values instead of a single one. + The QEOM and NumpyEigensolver are adapted to the new structure. + A factory is introduced to run a numpy eigensolver with a specific filter + (to target states of specific symmetries). +features: + - ExcitedStateSolver + - ExcitedStateEigensolver + - QEOM + - EigensolverFactory + - NumPyEigenSolverFactory + - ElectronicStructureResult diff --git a/releasenotes/notes/ground_state_interface-42576cb6658a46e0.yaml b/releasenotes/notes/ground_state_interface-42576cb6658a46e0.yaml index d4862f9cc9..41f08d63bc 100644 --- a/releasenotes/notes/ground_state_interface-42576cb6658a46e0.yaml +++ b/releasenotes/notes/ground_state_interface-42576cb6658a46e0.yaml @@ -10,12 +10,12 @@ features: Introduces ``chemistry/results`` where the eigenstate_result and the electronic_structure_result are also used for the algorithms. Introduces Minimum Eigensolver factories ``minimum_eigensolver_factories`` where chemistry specific - minimum eigensolvers can be initialised + minimum eigensolvers can be initialized Introduces orbital optimization vqe ``oovqe`` as a ground state solver for chemistry applications deprecations: - | ``Core Hamiltonian`` class is deprecated in favor of the ``FermionicTransformation`` - ``Chemistry Operator`` class is deprecated in facor of the ``tranformations`` - ``minimum_eigen_solvers/vqe_adapt`` is also deprecate and moved as an implementation + ``Chemistry Operator`` class is deprecated in favor of the ``tranformations`` + ``minimum_eigen_solvers/vqe_adapt`` is also deprecated and moved as an implementation of the ground_state_solver interface - ``applications/molecular_ground_state_energy`` is deprecated in favor of ``ground_state_solver`` \ No newline at end of file + ``applications/molecular_ground_state_energy`` is deprecated in favor of ``ground_state_solver`` diff --git a/test/chemistry/test_adapt_vqe.py b/test/chemistry/test_adapt_vqe.py index 573dd2bdea..eee14a0a64 100644 --- a/test/chemistry/test_adapt_vqe.py +++ b/test/chemistry/test_adapt_vqe.py @@ -49,7 +49,7 @@ def test_default(self): solver = VQEUCCSDFactory(QuantumInstance(BasicAer.get_backend('statevector_simulator'))) calc = AdaptVQE(self.transformation, solver) res = calc.solve(self.driver) - self.assertAlmostEqual(res.electronic_energy, self.expected, places=6) + self.assertAlmostEqual(res.electronic_energies[0], self.expected, places=6) def test_custom_minimum_eigensolver(self): """ Test custom MES """ @@ -80,7 +80,7 @@ def get_solver(self, transformation): calc = AdaptVQE(self.transformation, solver) res = calc.solve(self.driver) - self.assertAlmostEqual(res.electronic_energy, self.expected, places=6) + self.assertAlmostEqual(res.electronic_energies[0], self.expected, places=6) def test_custom_excitation_pool(self): """ Test custom excitation pool """ @@ -101,7 +101,7 @@ def get_solver(self, transformation): solver = CustomFactory(QuantumInstance(BasicAer.get_backend('statevector_simulator'))) calc = AdaptVQE(self.transformation, solver) res = calc.solve(self.driver) - self.assertAlmostEqual(res.electronic_energy, self.expected, places=6) + self.assertAlmostEqual(res.electronic_energies[0], self.expected, places=6) def test_vqe_adapt_check_cyclicity(self): """ VQEAdapt index cycle detection """ diff --git a/test/chemistry/test_bopes_sampler.py b/test/chemistry/test_bopes_sampler.py index dfec35f7c2..97a7ec183e 100644 --- a/test/chemistry/test_bopes_sampler.py +++ b/test/chemistry/test_bopes_sampler.py @@ -18,6 +18,7 @@ import numpy as np from qiskit import BasicAer from qiskit.aqua import QuantumInstance +from qiskit.aqua import aqua_globals from qiskit.aqua.algorithms import VQE, NumPyMinimumEigensolver from qiskit.aqua.components.optimizers import AQGD from qiskit.aqua.operators import PauliExpectation @@ -35,7 +36,8 @@ class TestBOPES(unittest.TestCase): def test_h2_bopes_sampler(self): """Test BOPES Sampler on H2""" - np.random.seed(100) + seed = 50 + aqua_globals.random_seed = seed # Molecule dof = partial(Molecule.absolute_distance, atom_pair=(1, 0)) @@ -52,8 +54,8 @@ def test_h2_bopes_sampler(self): shots = 1 backend = 'statevector_simulator' quantum_instance = QuantumInstance(BasicAer.get_backend(backend), shots=shots) - quantum_instance.run_config.seed_simulator = 50 - quantum_instance.compile_config['seed_transpiler'] = 50 + quantum_instance.run_config.seed_simulator = seed + quantum_instance.compile_config['seed_transpiler'] = seed # Variational form i_state = HartreeFock(num_orbitals=f_t._molecule_info['num_orbitals'], @@ -100,7 +102,8 @@ def test_h2_bopes_sampler(self): def test_potential_interface(self): """Tests potential interface.""" - np.random.seed(100) + seed = 50 + aqua_globals.random_seed = seed stretch = partial(Molecule.absolute_distance, atom_pair=(1, 0)) # H-H molecule near equilibrium geometry diff --git a/test/chemistry/test_driver_methods_gsc.py b/test/chemistry/test_driver_methods_gsc.py index f10dfd9657..7f18805cf6 100644 --- a/test/chemistry/test_driver_methods_gsc.py +++ b/test/chemistry/test_driver_methods_gsc.py @@ -57,11 +57,11 @@ def _run_driver(driver, transformation=TransformationType.FULL, return result def _assert_energy(self, result, mol): - self.assertAlmostEqual(self.ref_energies[mol], result.energy, places=3) + self.assertAlmostEqual(self.ref_energies[mol], result.total_energies[0], places=3) def _assert_energy_and_dipole(self, result, mol): self._assert_energy(result, mol) - self.assertAlmostEqual(self.ref_dipoles[mol], result.total_dipole_moment, places=3) + self.assertAlmostEqual(self.ref_dipoles[mol], result.total_dipole_moment[0], places=3) if __name__ == '__main__': diff --git a/test/chemistry/test_excited_states_solvers.py b/test/chemistry/test_excited_states_solvers.py new file mode 100644 index 0000000000..9b323481b1 --- /dev/null +++ b/test/chemistry/test_excited_states_solvers.py @@ -0,0 +1,102 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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. + +""" Test NumericalqEOM excited states calculation """ + +import unittest +from test.chemistry import QiskitChemistryTestCase +import numpy as np + +from qiskit import BasicAer +from qiskit.aqua import aqua_globals, QuantumInstance + +from qiskit.aqua.algorithms import NumPyMinimumEigensolver, NumPyEigensolver +from qiskit.chemistry import QiskitChemistryError +from qiskit.chemistry.drivers import PySCFDriver, UnitsType +from qiskit.chemistry.transformations import FermionicTransformation +from qiskit.chemistry.transformations.fermionic_transformation import QubitMappingType +from qiskit.chemistry.algorithms.ground_state_solvers import (GroundStateEigensolver, + VQEUCCSDFactory) +from qiskit.chemistry.algorithms.excited_states_solvers import ( + NumPyEigensolverFactory, ExcitedStatesEigensolver, QEOM +) + + +class TestNumericalQEOMESCCalculation(QiskitChemistryTestCase): + """ Test NumericalqEOM excited states calculation """ + + def setUp(self): + super().setUp() + aqua_globals.random_seed = 8 + try: + self.driver = PySCFDriver(atom='H .0 .0 .0; H .0 .0 0.75', + unit=UnitsType.ANGSTROM, + charge=0, + spin=0, + basis='sto3g') + except QiskitChemistryError: + self.skipTest('PYSCF driver does not appear to be installed') + + self.reference_energies = [-1.8427016, -1.8427016 + 0.5943372, -1.8427016 + 0.95788352, + -1.8427016 + 1.5969296] + self.transformation = FermionicTransformation(qubit_mapping=QubitMappingType.JORDAN_WIGNER) + solver = NumPyEigensolver() + self.ref = solver + self.quantum_instance = QuantumInstance(BasicAer.get_backend('statevector_simulator'), + seed_transpiler=90, seed_simulator=12) + + def test_numpy_mes(self): + """ Test NumPyMinimumEigenSolver with QEOM """ + solver = NumPyMinimumEigensolver() + gsc = GroundStateEigensolver(self.transformation, solver) + esc = QEOM(gsc, 'sd') + results = esc.solve(self.driver) + + for idx in range(len(self.reference_energies)): + self.assertAlmostEqual(results.computed_energies[idx], self.reference_energies[idx], + places=4) + + def test_vqe_mes(self): + """ Test VQEUCCSDFactory with QEOM """ + solver = VQEUCCSDFactory(self.quantum_instance) + gsc = GroundStateEigensolver(self.transformation, solver) + esc = QEOM(gsc, 'sd') + results = esc.solve(self.driver) + + for idx in range(len(self.reference_energies)): + self.assertAlmostEqual(results.computed_energies[idx], self.reference_energies[idx], + places=4) + + def test_numpy_factory(self): + """ Test NumPyEigenSolverFactory with ExcitedStatesEigensolver """ + + # pylint: disable=unused-argument + def filter_criterion(eigenstate, eigenvalue, aux_values): + return np.isclose(aux_values[0][0], 2.) + + solver = NumPyEigensolverFactory(filter_criterion=filter_criterion) + esc = ExcitedStatesEigensolver(self.transformation, solver) + results = esc.solve(self.driver) + + # filter duplicates from list + computed_energies = [results.computed_energies[0]] + for comp_energy in results.computed_energies[1:]: + if not np.isclose(comp_energy, computed_energies[-1]): + computed_energies.append(comp_energy) + + for idx in range(len(self.reference_energies)): + self.assertAlmostEqual(computed_energies[idx], self.reference_energies[idx], + places=4) + + +if __name__ == '__main__': + unittest.main() diff --git a/test/chemistry/test_groundstate_eigensolver.py b/test/chemistry/test_groundstate_eigensolver.py index 2627d85a9f..7f9244b838 100644 --- a/test/chemistry/test_groundstate_eigensolver.py +++ b/test/chemistry/test_groundstate_eigensolver.py @@ -50,21 +50,21 @@ def test_npme(self): solver = NumPyMinimumEigensolverFactory() calc = GroundStateEigensolver(self.transformation, solver) res = calc.solve(self.driver) - self.assertAlmostEqual(res.energy, self.reference_energy, places=6) + self.assertAlmostEqual(res.total_energies[0], self.reference_energy, places=6) def test_npme_with_default_filter(self): """ Test NumPyMinimumEigensolver with default filter """ solver = NumPyMinimumEigensolverFactory(use_default_filter_criterion=True) calc = GroundStateEigensolver(self.transformation, solver) res = calc.solve(self.driver) - self.assertAlmostEqual(res.energy, self.reference_energy, places=6) + self.assertAlmostEqual(res.total_energies[0], self.reference_energy, places=6) def test_vqe_uccsd(self): """ Test VQE UCCSD case """ solver = VQEUCCSDFactory(QuantumInstance(BasicAer.get_backend('statevector_simulator'))) calc = GroundStateEigensolver(self.transformation, solver) res = calc.solve(self.driver) - self.assertAlmostEqual(res.energy, self.reference_energy, places=6) + self.assertAlmostEqual(res.total_energies[0], self.reference_energy, places=6) def _setup_evaluation_operators(self): # first we run a ground state calculation diff --git a/test/chemistry/test_swaprz.py b/test/chemistry/test_swaprz.py index 4cf2419a0b..76febe3a42 100644 --- a/test/chemistry/test_swaprz.py +++ b/test/chemistry/test_swaprz.py @@ -67,8 +67,7 @@ def test_excitation_preserving(self): gsc = GroundStateEigensolver(fermionic_transformation, solver) result = gsc.solve(driver) - - self.assertAlmostEqual(result.energy, self.reference_energy, places=4) + self.assertAlmostEqual(result.total_energies[0], self.reference_energy, places=4) if __name__ == '__main__': diff --git a/test/chemistry/test_symmetries.py b/test/chemistry/test_symmetries.py index 38a5a62aa9..930907aaf0 100644 --- a/test/chemistry/test_symmetries.py +++ b/test/chemistry/test_symmetries.py @@ -103,8 +103,7 @@ def test_tapered_op(self): gsc = GroundStateEigensolver(self.fermionic_transformation, solver) result = gsc.solve(self.driver) - - self.assertAlmostEqual(result.energy, self.reference_energy, places=6) + self.assertAlmostEqual(result.total_energies[0], self.reference_energy, places=6) if __name__ == '__main__': diff --git a/test/chemistry/test_uccsd_advanced.py b/test/chemistry/test_uccsd_advanced.py index 6e13d0a8d8..e85699c9c9 100644 --- a/test/chemistry/test_uccsd_advanced.py +++ b/test/chemistry/test_uccsd_advanced.py @@ -27,7 +27,6 @@ from qiskit.chemistry.core import TransformationType, QubitMappingType from qiskit.chemistry.transformations import FermionicTransformation - # pylint: disable=invalid-name @@ -100,7 +99,7 @@ def test_uccsd_hf_qpUCCD(self): result = gsc.solve(self.driver) - self.assertAlmostEqual(result.energy, self.reference_energy_pUCCD, places=6) + self.assertAlmostEqual(result.total_energies[0], self.reference_energy_pUCCD, places=6) def test_uccsd_hf_qUCCD0(self): """ singlet uccd test """ @@ -133,7 +132,7 @@ def test_uccsd_hf_qUCCD0(self): result = gsc.solve(self.driver) - self.assertAlmostEqual(result.energy, self.reference_energy_UCCD0, places=6) + self.assertAlmostEqual(result.total_energies[0], self.reference_energy_UCCD0, places=6) def test_uccsd_hf_qUCCD0full(self): """ singlet full uccd test """ @@ -166,8 +165,7 @@ def test_uccsd_hf_qUCCD0full(self): gsc = GroundStateEigensolver(self.fermionic_transformation, solver) result = gsc.solve(self.driver) - - self.assertAlmostEqual(result.energy, self.reference_energy_UCCD0full, places=6) + self.assertAlmostEqual(result.total_energies[0], self.reference_energy_UCCD0full, places=6) def test_uccsd_hf_qUCCSD(self): """ uccsd tapering test using all double excitations """ @@ -215,7 +213,7 @@ def test_uccsd_hf_qUCCSD(self): raw_result = solver.compute_minimum_eigenvalue(qubit_op, None) result = fermionic_transformation.interpret(raw_result) - self.assertAlmostEqual(result.energy, self.reference_energy_UCCSD, places=6) + self.assertAlmostEqual(result.total_energies[0], self.reference_energy_UCCSD, places=6) def test_uccsd_hf_excitations(self): """ uccsd tapering test using all double excitations """ diff --git a/test/chemistry/test_uccsd_hartree_fock.py b/test/chemistry/test_uccsd_hartree_fock.py index 5a94152f9e..1104bec3bd 100644 --- a/test/chemistry/test_uccsd_hartree_fock.py +++ b/test/chemistry/test_uccsd_hartree_fock.py @@ -68,7 +68,7 @@ def test_uccsd_hf(self): result = gsc.solve(self.driver) - self.assertAlmostEqual(result.energy, self.reference_energy, places=6) + self.assertAlmostEqual(result.total_energies[0], self.reference_energy, places=6) def test_uccsd_hf_qasm(self): """ uccsd hf test with qasm_simulator. """ @@ -83,8 +83,7 @@ def test_uccsd_hf_qasm(self): gsc = GroundStateEigensolver(self.fermionic_transformation, solver) result = gsc.solve(self.driver) - - self.assertAlmostEqual(result.energy, -1.138, places=2) + self.assertAlmostEqual(result.total_energies[0], -1.138, places=2) def test_uccsd_hf_aer_statevector(self): """ uccsd hf test with Aer statevector """ @@ -101,8 +100,7 @@ def test_uccsd_hf_aer_statevector(self): gsc = GroundStateEigensolver(self.fermionic_transformation, solver) result = gsc.solve(self.driver) - - self.assertAlmostEqual(result.energy, self.reference_energy, places=6) + self.assertAlmostEqual(result.total_energies[0], self.reference_energy, places=6) def test_uccsd_hf_aer_qasm(self): """ uccsd hf test with Aer qasm_simulator. """ @@ -123,8 +121,7 @@ def test_uccsd_hf_aer_qasm(self): gsc = GroundStateEigensolver(self.fermionic_transformation, solver) result = gsc.solve(self.driver) - - self.assertAlmostEqual(result.energy, -1.138, places=2) + self.assertAlmostEqual(result.total_energies[0], -1.138, places=2) def test_uccsd_hf_aer_qasm_snapshot(self): """ uccsd hf test with Aer qasm_simulator snapshot. """ @@ -143,7 +140,7 @@ def test_uccsd_hf_aer_qasm_snapshot(self): gsc = GroundStateEigensolver(self.fermionic_transformation, solver) result = gsc.solve(self.driver) - self.assertAlmostEqual(result.energy, self.reference_energy, places=3) + self.assertAlmostEqual(result.total_energies[0], self.reference_energy, places=3) EXCITATION_RESULTS = \ [[[[0, 1], [0, 2], [3, 4], [3, 5]],