diff --git a/electronic_structure_algorithms/__init__.py b/electronic_structure_algorithms/__init__.py new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/electronic_structure_algorithms/__init__.py @@ -0,0 +1 @@ + diff --git a/electronic_structure_algorithms/excited_states_eigensolvers/__init__.py b/electronic_structure_algorithms/excited_states_eigensolvers/__init__.py new file mode 100644 index 0000000..9711027 --- /dev/null +++ b/electronic_structure_algorithms/excited_states_eigensolvers/__init__.py @@ -0,0 +1,3 @@ +from .ssvqe import SSVQE, SSVQEResult +from .vqd import VQD, VQDResult +from .mcvqe import MCVQE, MCVQEResult \ No newline at end of file diff --git a/electronic_structure_algorithms/excited_states_eigensolvers/mcvqe.py b/electronic_structure_algorithms/excited_states_eigensolvers/mcvqe.py new file mode 100644 index 0000000..f9be54f --- /dev/null +++ b/electronic_structure_algorithms/excited_states_eigensolvers/mcvqe.py @@ -0,0 +1,524 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2018, 2022. +# +# 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 multi-configuration Variational Quantum Eigensolver algorithm. + +See https://arxiv.org/abs/1810.09434. This implementation is derived +from Qiskit's implemtation of VQE. +""" + +from __future__ import annotations + +import logging +import warnings +from time import time +from collections.abc import Callable, Sequence + +import numpy as np + +from qiskit_algorithms.gradients import BaseEstimatorGradient +from qiskit.circuit import QuantumCircuit +from qiskit.opflow import PauliSumOp +from qiskit.primitives import BaseEstimator +from qiskit.quantum_info.operators.base_operator import BaseOperator +from qiskit.utils import algorithm_globals +from qiskit.quantum_info import Statevector + +from qiskit_algorithms.exceptions import AlgorithmError +from qiskit_algorithms.list_or_dict import ListOrDict +from qiskit_algorithms.optimizers import Optimizer, Minimizer, OptimizerResult +from qiskit_algorithms.eigensolvers.eigensolver import Eigensolver, EigensolverResult + +from qiskit_algorithms.observables_evaluator import estimate_observables +from .ssvqe import SSVQE, SSVQEResult +from ..initializations.configuration_interaction_states import get_CIS_states, get_CISD_states + +logger = logging.getLogger(__name__) + +class MCVQE(SSVQE): + + r"""The Multi-Configuration Variational Quantum Eigensolver algorithm. + `MCVQE `__ is a hybrid quantum-classical + algorithm that uses a variational technique to find the low-lying eigenvalues + of the Hamiltonian :math:`H` of a given system. MCVQE can be seen as a + special case of the Subspace Search Variational Quantum Eigensolver (SSVQE). + The primary difference is that MCVQE chooses all of the weights in the + objective function to be equal. Thus, its global minimia consist of + states which span the low-lying eigenspace of the Hamiltonian in question + rather than the eigenvectors themselves. Because of this, an extra + post-processing step involving the diagonalization of the contracted + Hamiltonian in this subspace is required to obtain the eigenstates themselves. + If the initial states are chosen to be a set of mutually orthogonal states + whose configuration interaction (CI) coefficients are known, then + the off-diagonal elements of this contracted Hamiltonian can be computed without + (potentially expensive) Hamard test-style circuits for computing inner product + terms. It is worth notating that mutually orthogonal sets of Slater + determinants fulfill this criteria. + + An instance of MCVQE requires defining four algorithmic sub-components: + + An :attr:`estimator` to compute the expectation values of operators, an integer ``k`` denoting + the number of eigenstates that the algorithm will attempt to find, an ansatz which is a + :class:`QuantumCircuit`, and one of the classical :mod:`~qiskit.algorithms.optimizers`. + + The ansatz is varied, via its set of parameters, by the optimizer, such that it works towards + a set of mutually orthogonal states which span the low-lying eigenspace of + the Hamiltonian :math:`H`. An optional array of parameter values, via the + ``initial_point``, may be provided as the starting point for the search of the + low-lying eigenvalues. This feature is particularly useful such as + when there are reasons to believe that the solution point is close to a particular + point. The length of the ``initial_point`` list value must match the number of the + parameters expected by the ansatz being used. If the ``initial_point`` is left at the + default of ``None``, then MCVQE will look to the ansatz for a preferred value, based + on its given initial state. If the ansatz returns ``None``, then a random point will + be generated within the parameter bounds set, as per above. If the ansatz provides + ``None`` as the lower bound, then MCVQE will default it to :math:`-2\pi`; similarly, + if the ansatz returns ``None`` as the upper bound, the default value will be :math:`2\pi`. + + There are two ways to provide initial states to MCVQE. The first is to provide values + for ``excitations``, ``one_body_integrals``, and ``two_body_integrals``, and + ``num_particles``. The one and two body integrals must be the same array + instance as those used to calculate the Hamiltonian and the two body integrals must be + in the Physicist notation. These integrals are used to compute the matrix elements of a + contracted configuration interaction Hamiltonian to be diagonalized, yielding CI + statevectors. The value assigned to ``excitations`` determins the excitation level + of this CI problem. Setting ``excitations`` to ``s`` will use configuration interaction + singles (CIS) states and setting it to ``sd`` will use configuration interaction singles + and doubles (CISD). Similarly, ``num_spin_orbitals`` and ``num_particles`` are also used + in this CI statevector construction and must match the number of spin-orbitals and + particles of the problem in question. + + An optional list of initial states, via the ``initial_states``, may also be provided + if the method described above is not used. Choosing these states appropriately + is a critical part of the algorithm. They must be mutually orthogonal because this + is how the algorithm enforces the mutual orthogonality of the solution states. + Additionally, if one wishes to find the low-lying excited states of a molecular + Hamiltonian, then we expect the output states to belong to a particular + particle-number subspace. If an ansatz that preserves particle number such as + :class:`UCCSD` is used, then states belonging to the incorrect particle + number subspace will be returned if the ``initial_states`` are not in the correct + particle number subspace. A similar statement can often be made for the + spin-magnetization quantum number. + + The following attributes can be set via the initializer but can also be read and + updated once the MCVQE object has been constructed. + + Attributes: + estimator (BaseEstimator): The primitive instance used to perform the expectation + estimation of observables. + k (int): The number of eigenstates that SSVQE will attempt to find. + ansatz (QuantumCircuit): A parameterized circuit used as an ansatz for the + wave function. + optimizer (Optimizer): A classical optimizer, which can be either a Qiskit optimizer + or a callable that takes an array as input and returns a Qiskit or SciPy optimization + result. + gradient (BaseEstimatorGradient | None): An optional estimator gradient to be used with the + optimizer. + initial_states (Sequence[QuantumCircuit]): An optional list of mutually orthogonal + initial states. If ``None``, then MCVQE will set these to be a list of mutually + orthogonal computational basis states. + callback (Callable[[int, np.ndarray, Sequence[float], dict[str, Any]], None] | None): A + function that can access the intermediate data at each optimization step. These data are + the evaluation count, the optimizer parameters for the ansatz, the evaluated mean + energies, and the metadata dictionary. + check_input_states_orthogonality: A boolean that sets whether or not to check + that the value of initial_states passed consists of a mutually orthogonal + set of states. If ``True``, then MCVQE will check that these states are mutually + orthogonal and return an :class:`AlgorithmError` if they are not. + This is set to ``True`` by default, but setting this to ``False`` may be desirable + for larger numbers of qubits to avoid exponentially large computational overhead. + """ + + def __init__( + self, + estimator: BaseEstimator, + num_particles: tuple | None = None, + one_body_integrals: np.ndarray | None = None, + two_body_integrals: np.ndarray | None = None, + k: int | None = 2, + ansatz: QuantumCircuit | None = None, + optimizer: Optimizer | Minimizer | None = None, + initial_point: Sequence[float] | None = None, + gradient: BaseEstimatorGradient | None = None, + callback: Callable[[int, np.ndarray, Sequence[float], float], None] | None = None, + check_input_states_orthogonality: bool = True, + excitations: str = None, + initial_states: list[QuantumCircuit] = None, + ) -> None: + """ + Args: + estimator: The estimator primitive. + num_particles: A tuple whose entries denote the number of alpha and beta electrons + in the problem Hamiltonian. + one_body_integrals: The one body integrals of the problem Hamiltonian. + two_body_integrals: The two body integrals of the problem Hamiltonian in the + Physicist's notation. + k: The number of eigenstates that the algorithm will attempt to find. + ansatz: A parameterized circuit used as Ansatz for the wave function. + optimizer: A classical optimizer. Can either be a Qiskit optimizer or a callable + that takes an array as input and returns a Qiskit or SciPy optimization result. + initial_point: An optional initial point (i.e. initial parameter values) + for the optimizer. If ``None`` then VQE will look to the ansatz for a preferred + point and if not will simply compute a random one. + gradient: An optional gradient function or operator for optimizer. + callback: A callback that can access the intermediate data at each optimization step. + These data are: the evaluation count, the optimizer ansatz parameters, + the evaluated mean energies, and the metadata dictionary. + check_input_states_orthogonality: A boolean that sets whether or not to check + that the value of ``initial_states`` passed consists of a mutually orthogonal + set of states. If ``True``, then SSVQE will check that these states are mutually + orthogonal and return an error if they are not. This is set to ``True`` by default, + but setting this to ``False`` may be desirable for larger numbers of qubits to avoid + exponentially large computational overhead before the simulation even starts. + excitations: One of `s` or `sd` to determine whether CIS or CISD states will be used for + initialization. + initial_states: An optional list of mutually orthogonal initial states. + """ + + super().__init__(estimator=estimator, + k=k, + ansatz=ansatz, + optimizer=optimizer, + initial_point=initial_point, + gradient=gradient, + callback=callback, + check_input_states_orthogonality=check_input_states_orthogonality) + + self.k = k + self.weight_vector = [1 for n in range(k)] + self.ansatz = ansatz + self.optimizer = optimizer + self.initial_point = initial_point + self.gradient = gradient + self.callback = callback + self.estimator = estimator + self.check_initial_states_orthogonal = check_input_states_orthogonality + self._excitations = excitations + self._initial_states = initial_states + self._one_body_integrals = one_body_integrals + self._two_body_integrals = two_body_integrals + self._off_diagonal_plus_states = None + self._off_diagonal_minus_states = None + self.num_particles = num_particles + + @property + def initial_states(self) -> Sequence[QuantumCircuit] | np.ndarray | None: + + return self._initial_states + + @initial_states.setter + def initial_states(self, states) -> None: + + self._initial_states = states + + @property + def excitations(self) -> str | None: + + return self._excitations + + @excitations.setter + def excitations(self, excitations_string) -> None: + + self._excitations = excitations_string + + @property + def one_body_integrals(self) -> np.ndarray | None: + + return self._one_body_integrals + + @one_body_integrals.setter + def one_body_integrals(self, integrals) -> None: + + self._one_body_integrals = integrals + + @property + def two_body_integrals(self) -> np.ndarray | None: + + return self._two_body_integrals + + @two_body_integrals.setter + def two_body_integrals(self, integrals) -> None: + + self._two_body_integrals = integrals + + def initialize_mcvqe(self) -> None: + + if self.excitations is None and self.initial_states is not None: + + for n, state in enumerate(initial_states): + + if isinstance(state, QuantumCircuit): + + initial_states[n] = np.asarray(Statevector(state)) + + self.initial_states = initial_states + + if self.excitations == 's': + + num_spin_orbitals = self.one_body_integrals.shape[0] + initial_states = get_CIS_states(one_body_integrals=self.one_body_integrals, + two_body_integrals=self.two_body_integrals, + num_particles=self.num_particles, + state_representation='dense')[:self.k] + + if self.excitations == 'sd': + + num_spin_orbitals = self.one_body_integrals.shape[0] + initial_states = get_CISD_states(one_body_integrals=self.one_body_integrals, + two_body_integrals=self.two_body_integrals, + num_particles=self.num_particles, + num_spin_orbitals=num_spin_orbitals, + state_representation='dense')[:self.k] + + self.initial_states = [QuantumCircuit(self.ansatz.num_qubits) for n in range(self.k)] + for n in range(self.k): + self.initial_states[n].initialize(initial_states[n], self.initial_states[n].qubits) + + self._off_diagonal_plus_states = np.full(fill_value=None, shape=(self.k,self.k)) + self._off_diagonal_minus_states = np.full(fill_value=None, shape=(self.k,self.k)) + + for i in range(self.k): + for j in range(i): + if i == j: + pass + else: + plus_state = QuantumCircuit(self.ansatz.num_qubits) + minus_state = QuantumCircuit(self.ansatz.num_qubits) + + plus_statevector = (initial_states[i] + initial_states[j])/np.sqrt(2) + plus_statevector = plus_statevector/np.sqrt(sum(amplitude**2 for amplitude in plus_statevector)) + + minus_statevector = (initial_states[i] - initial_states[j])/np.sqrt(2) + minus_statevector = minus_statevector/np.sqrt(sum(amplitude**2 for amplitude in minus_statevector)) + + self._off_diagonal_plus_states[i,j] = plus_state + self._off_diagonal_plus_states[i,j].initialize(plus_statevector, plus_state.qubits) + + self._off_diagonal_minus_states[i,j] = minus_state + self._off_diagonal_minus_states[i,j].initialize(minus_statevector, minus_state.qubits) + + + def compute_eigenvalues( + self, + operator: BaseOperator | PauliSumOp, + aux_operators: ListOrDict[BaseOperator | PauliSumOp] | None = None, + ) -> EigensolverResult: + + ansatz = self._check_operator_ansatz(operator) + + initial_point = _validate_initial_point(self.initial_point, ansatz) + + self.initialize_mcvqe() + + initial_states = self._check_operator_initial_states(self.initial_states, operator) + + bounds = _validate_bounds(ansatz) + + initialized_ansatz_list = [initial_states[n].compose(ansatz) for n in range(self.k)] + + self.weight_vector = self._check_weight_vector(self.weight_vector) + + evaluate_weighted_energy_sum = self._get_evaluate_weighted_energy_sum( + initialized_ansatz_list, operator + ) + + if self.gradient is not None: # need to implement _get_evaluate_gradient + evaluate_gradient = self._get_evalute_gradient(initialized_ansatz_list, operator) + else: + evaluate_gradient = None + + if aux_operators: + zero_op = PauliSumOp.from_list([("I" * self.ansatz.num_qubits, 0)]) + + # Convert the None and zero values when aux_operators is a list. + # Drop None and convert zero values when aux_operators is a dict. + if isinstance(aux_operators, list): + key_op_iterator = enumerate(aux_operators) + converted = [zero_op] * len(aux_operators) + else: + key_op_iterator = aux_operators.items() + converted = {} + for key, op in key_op_iterator: + if op is not None: + converted[key] = zero_op if op == 0 else op + + aux_operators = converted + + else: + aux_operators = None + + start_time = time() + + if callable(self.optimizer): + optimizer_result = self.optimizer( + fun=evaluate_weighted_energy_sum, + x0=initial_point, + jac=evaluate_gradient, + bounds=bounds, + ) + else: + optimizer_result = self.optimizer.minimize( + fun=evaluate_weighted_energy_sum, + x0=initial_point, + jac=evaluate_gradient, + bounds=bounds, + ) + + optimizer_time = time() - start_time + + logger.info( + "Optimization complete in %s seconds.\nFound opt_params %s", + optimizer_time, + optimizer_result.x, + ) + + if aux_operators is not None: + bound_ansatz_list = [ + initialized_ansatz_list[n].bind_parameters(optimizer_result.x) + for n in range(self.k) + ] + + aux_values_list = [ + estimate_observables( + self.estimator, + bound_ansatz_list[n], + aux_operators, + ) + for n in range(self.k) + ] + else: + aux_values_list = None + + return self._build_mcvqe_result( + optimizer_result, aux_values_list, optimizer_time, operator, initialized_ansatz_list + ) + + def _build_mcvqe_result( + self, + optimizer_result: OptimizerResult, + aux_operators_evaluated: ListOrDict[tuple[complex, tuple[complex, int]]], + optimizer_time: float, + operator: BaseOperator | PauliSumOp, + initialized_ansatz_list: list[QuantumCircuit], + ) -> MCVQEResult: + result = MCVQEResult() + + for i in range(self.k): + for j in range(i): + if i == j: + pass + else: + + self._off_diagonal_plus_states[i,j] = self._off_diagonal_plus_states[i,j].compose(self.ansatz) + self._off_diagonal_minus_states[i,j] = self._off_diagonal_minus_states[i,j].compose(self.ansatz) + + contracted_Hamiltonian = np.empty(shape=(self.k, self.k)) + for i in range(self.k): + for j in range(i): + + try: + + matrix_element_result = self.estimator.run([self._off_diagonal_plus_states[i,j], self._off_diagonal_minus_states[i,j]], + [operator]*2, + [optimizer_result.x]*2 + ).result().values + + contracted_Hamiltonian[i,j] = 0.5*(matrix_element_result[0] - matrix_element_result[1]) + if np.isclose(contracted_Hamiltonian[i,j].imag, 0.0): + contracted_Hamiltonian[j,i] = contracted_Hamiltonian[i,j] + else: + contracted_Hamiltonian[j,i] = np.conjugate(contracted_Hamiltonian[i,j]) + + except Exception as exc: + raise AlgorithmError("The primitive job to evaluate the eigenvalues failed!") from exc + try: + diagonal_elements = ( + self.estimator.run( + initialized_ansatz_list, [operator] * self.k, [optimizer_result.x] * self.k + ) + .result() + .values + ) + + for i in range(self.k): + + contracted_Hamiltonian[i,i] = diagonal_elements[i] + + except Exception as exc: + raise AlgorithmError("The primitive job to evaluate the eigenvalues failed!") from exc + + result.eigenvalues = np.linalg.eigh(np.asarray(contracted_Hamiltonian))[0] + + result.cost_function_evals = optimizer_result.nfev + result.optimal_point = optimizer_result.x + result.optimal_parameters = dict(zip(self.ansatz.parameters, optimizer_result.x)) + result.optimal_value = optimizer_result.fun + result.optimizer_time = optimizer_time + result.aux_operators_evaluated = aux_operators_evaluated + result.optimizer_result = optimizer_result + + return result + +class MCVQEResult(SSVQEResult): + """MCVQE Result.""" + + def __init__(self) -> None: + super().__init__() + + +def _validate_initial_point(point, ansatz): + expected_size = ansatz.num_parameters + + # try getting the initial point from the ansatz + if point is None and hasattr(ansatz, "preferred_init_points"): + point = ansatz.preferred_init_points + # if the point is None choose a random initial point + + if point is None: + # get bounds if ansatz has them set, otherwise use [-2pi, 2pi] for each parameter + bounds = getattr(ansatz, "parameter_bounds", None) + if bounds is None: + bounds = [(-2 * np.pi, 2 * np.pi)] * expected_size + + # replace all Nones by [-2pi, 2pi] + lower_bounds = [] + upper_bounds = [] + for lower, upper in bounds: + lower_bounds.append(lower if lower is not None else -2 * np.pi) + upper_bounds.append(upper if upper is not None else 2 * np.pi) + + # sample from within bounds + point = algorithm_globals.random.uniform(lower_bounds, upper_bounds) + + elif len(point) != expected_size: + raise ValueError( + f"The dimension of the initial point ({len(point)}) does not match the " + f"number of parameters in the circuit ({expected_size})." + ) + + return point + + +def _validate_bounds(ansatz): + if hasattr(ansatz, "parameter_bounds") and ansatz.parameter_bounds is not None: + bounds = ansatz.parameter_bounds + if len(bounds) != ansatz.num_parameters: + raise ValueError( + f"The number of bounds ({len(bounds)}) does not match the number of " + f"parameters in the circuit ({ansatz.num_parameters})." + ) + else: + bounds = [(None, None)] * ansatz.num_parameters + + return bounds + + \ No newline at end of file diff --git a/electronic_structure_algorithms/excited_states_eigensolvers/ssvqe.py b/electronic_structure_algorithms/excited_states_eigensolvers/ssvqe.py new file mode 100644 index 0000000..844454d --- /dev/null +++ b/electronic_structure_algorithms/excited_states_eigensolvers/ssvqe.py @@ -0,0 +1,627 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2018, 2022. +# +# 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 Subspace Search Variational Quantum Eigensolver algorithm. + +See https://arxiv.org/abs/1810.09434. This implementation is derived +from Qiskit's implementation of VQE. +""" + +from __future__ import annotations + +import logging +import warnings +from time import time +from collections.abc import Callable, Sequence + +import numpy as np + +from qiskit_algorithms.gradients import BaseEstimatorGradient +from qiskit.circuit import QuantumCircuit +from qiskit.circuit.library import RealAmplitudes +from qiskit.opflow import PauliSumOp +from qiskit.primitives import BaseEstimator +from qiskit.quantum_info.operators.base_operator import BaseOperator +from qiskit.utils import algorithm_globals +from qiskit.quantum_info import Statevector + +from qiskit_algorithms.exceptions import AlgorithmError +from qiskit_algorithms.list_or_dict import ListOrDict +from qiskit_algorithms.optimizers import Optimizer, Minimizer, OptimizerResult +from qiskit_algorithms.variational_algorithm import VariationalAlgorithm, VariationalResult +from qiskit_algorithms.eigensolvers.eigensolver import Eigensolver, EigensolverResult + +from qiskit_algorithms.observables_evaluator import estimate_observables + +logger = logging.getLogger(__name__) + + +class SSVQE(VariationalAlgorithm, Eigensolver): + r"""The Subspace Search Variational Quantum Eigensolver algorithm. + `SSVQE `__ is a hybrid quantum-classical + algorithm that uses a variational technique to find the low-lying eigenvalues + of the Hamiltonian :math:`H` of a given system. SSVQE can be seen as + a natural generalization of VQE. Whereas VQE minimizes the expectation + value of :math:`H` with respect to the ansatz state, SSVQE takes a set + of mutually orthogonal input states :math:`\{| \psi_{i}} \rangle\}_{i=0}^{k-1}`, + applies the same parameterized ansatz circuit :math:`U(\vec\theta)` to all of them, + then minimizes a weighted sum of the expectation values of :math:`H` with respect to these states: + + .. math:: + + \min_{\vec\theta}\sum_{i=0}^{k-1} w_{i}\langle\psi_{i}(\vec\theta)|H|\psi{i}(\vec\theta)\rangle + + where :math:`|\psi{i} (\vec\theta)\rangle` is shorthand for :math:`U (\vec\theta)| \psi_{i} \rangle` + and :math:`\{ w_{i} \}_{i=0}^{k-1}` are the components of the ``weight_vector``. + + An instance of SSVQE requires defining four algorithmic sub-components: + + An :attr:`estimator` to compute the expectation values of operators, an integer ``k`` denoting + the number of eigenstates that the algorithm will attempt to find, an ansatz which is a + :class:`QuantumCircuit`, and one of the classical :mod:`~qiskit.algorithms.optimizers`. + + The ansatz is varied, via its set of parameters, by the optimizer, such that it works towards + a set of mutually orthogonal states, as determined by the parameters applied to the ansatz, + that will result in the minimum weighted sum of expectation values being measured + of the input operator (Hamiltonian) with respect to these states. The weights given + to this list of expectation values is given by ``weight_vector``. An optional array of + parameter values, via the ``initial_point``, may be provided as the starting point for + the search of the low-lying eigenvalues. This feature is particularly useful such as + when there are reasons to believe that the solution point is close to a particular + point. The length of the ``initial_point`` list value must match the number of the + parameters expected by the ansatz being used. If the ``initial_point`` is left at the + default of ``None``, then SSVQE will look to the ansatz for a preferred value, based + on its given initial state. If the ansatz returns ``None``, then a random point will + be generated within the parameter bounds set, as per above. If the ansatz provides + ``None`` as the lower bound, then SSVQE will default it to :math:`-2\pi`; similarly, + if the ansatz returns ``None`` as the upper bound, the default value will be :math:`2\pi`. + + An optional list of initial states, via the ``initial_states``, may also be provided. + Choosing these states appropriately is a critical part of the algorithm. They must + be mutually orthogonal because this is how the algorithm enforces the mutual + orthogonality of the solution states. If the ``initial_states`` is left as ``None``, + then SSVQE will automatically generate a list of computational basis states and use + these as the initial states. For many physically-motivated problems, it is advised + to not rely on these default values as doing so can easily result in an unphysical + solution being returned. For example, if one wishes to find the low-lying excited + states of a molecular Hamiltonian, then we expect the output states to belong to a + particular particle-number subspace. If an ansatz that preserves particle number + such as :class:`UCCSD` is used, then states belonging to the incorrect particle + number subspace will be returned if the ``initial_states`` are not in the correct + particle number subspace. A similar statement can often be made for the + spin-magnetization quantum number. + + A minimal example of how one may initialize an instance of ``SSVQE`` and use it + to compute the low-lying eigenvalues of an operator: + + .. code-block:: python + + from qiskit import QuantumCircuit + from qiskit.quantum_info import Pauli + from qiskit.primitives import Estimator + from qiskit.circuit.library import RealAmplitudes + from qiskit.algorithms.optimizers import SPSA + from qiskit.algorithms.eigensolvers import SSVQE + + operator = Pauli("ZZ") + input_states = [QuantumCircuit(2), QuantumCircuit(2)] + input_states[0].x(0) + input_states[1].x(1) + + ssvqe_instance = SSVQE(k=2, + estimator=Estimator(), + optimizer=SPSA(), + ansatz=RealAmplitudes(2), + initial_states=input_states) + + result = ssvqe_instance.compute_eigenvalues(operator) + + The following attributes can be set via the initializer but can also be read and + updated once the SSVQE object has been constructed. + + Attributes: + estimator (BaseEstimator): The primitive instance used to perform the expectation + estimation of observables. + k (int): The number of eigenstates that SSVQE will attempt to find. + ansatz (QuantumCircuit): A parameterized circuit used as an ansatz for the + wave function. + optimizer (Optimizer): A classical optimizer, which can be either a Qiskit optimizer + or a callable that takes an array as input and returns a Qiskit or SciPy optimization + result. + gradient (BaseEstimatorGradient | None): An optional estimator gradient to be used with the + optimizer. + initial_states (Sequence[QuantumCircuit]): An optional list of mutually orthogonal + initial states. If ``None``, then SSVQE will set these to be a list of mutually + orthogonal computational basis states. + weight_vector (Sequence[float]): A 1D array of real positive-valued numbers to assign + as weights to each of the expectation values. If ``None``, then SSVQE will default + to [k, k-1, ..., 1]. + callback (Callable[[int, np.ndarray, Sequence[float], dict[str, Any]], None] | None): A + function that can access the intermediate data at each optimization step. These data are + the evaluation count, the optimizer parameters for the ansatz, the evaluated mean + energies, and the metadata dictionary. + check_input_states_orthogonality: A boolean that sets whether or not to check + that the value of initial_states passed consists of a mutually orthogonal + set of states. If ``True``, then SSVQE will check that these states are mutually + orthogonal and return an :class:`AlgorithmError` if they are not. + This is set to ``True`` by default, but setting this to ``False`` may be desirable + for larger numbers of qubits to avoid exponentially large computational overhead. + """ + + def __init__( + self, + estimator: BaseEstimator, + k: int | None = 2, + ansatz: QuantumCircuit | None = None, + optimizer: Optimizer | Minimizer | None = None, + initial_point: Sequence[float] | None = None, + initial_states: list[QuantumCircuit] | None = None, + weight_vector: Sequence[float] | Sequence[int] | None = None, + gradient: BaseEstimatorGradient | None = None, + callback: Callable[[int, np.ndarray, Sequence[float], float], None] | None = None, + check_input_states_orthogonality: bool = True, + ) -> None: + """ + Args: + estimator: The estimator primitive. + k: The number of eigenstates that the algorithm will attempt to find. + ansatz: A parameterized circuit used as Ansatz for the wave function. + optimizer: A classical optimizer. Can either be a Qiskit optimizer or a callable + that takes an array as input and returns a Qiskit or SciPy optimization result. + initial_point: An optional initial point (i.e. initial parameter values) + for the optimizer. If ``None`` then VQE will look to the ansatz for a preferred + point and if not will simply compute a random one. + initial_states: An optional list of mutually orthogonal initial states. + If ``None``, then SSVQE will set these to be a list of mutually orthogonal + computational basis states. + weight_vector: An optional list or array of real positive numbers with length + equal to the value of ``num_states`` to be used in the weighted energy summation + objective function. This fixes the ordering of the returned eigenstate/eigenvalue + pairs. If ``None``, then SSVQE will default to [n, n-1, ..., 1] for `k` = n. + gradient: An optional gradient function or operator for optimizer. + callback: A callback that can access the intermediate data at each optimization step. + These data are: the evaluation count, the optimizer ansatz parameters, + the evaluated mean energies, and the metadata dictionary. + check_input_states_orthogonality: A boolean that sets whether or not to check + that the value of ``initial_states`` passed consists of a mutually orthogonal + set of states. If ``True``, then SSVQE will check that these states are mutually + orthogonal and return an error if they are not. This is set to ``True`` by default, + but setting this to ``False`` may be desirable for larger numbers of qubits to avoid + exponentially large computational overhead before the simulation even starts. + """ + + super().__init__() + + self.k = k + self.initial_states = initial_states + self.weight_vector = weight_vector + self.ansatz = ansatz + self.optimizer = optimizer + self.initial_point = initial_point + self.gradient = gradient + self.callback = callback + self.estimator = estimator + self.check_initial_states_orthogonal = check_input_states_orthogonality + + @property + def initial_point(self) -> Sequence[float] | None: + """Returns initial point""" + return self._initial_point + + @initial_point.setter + def initial_point(self, initial_point: Sequence[float] | None): + """Sets initial point""" + self._initial_point = initial_point + + @classmethod + def supports_aux_operators(cls) -> bool: + return True + + def compute_eigenvalues( + self, + operator: BaseOperator | PauliSumOp, + aux_operators: ListOrDict[BaseOperator | PauliSumOp] | None = None, + ) -> EigensolverResult: + + ansatz = self._check_operator_ansatz(operator) + + initial_point = _validate_initial_point(self.initial_point, ansatz) + + initial_states = self._check_operator_initial_states(self.initial_states, operator) + + bounds = _validate_bounds(ansatz) + + initialized_ansatz_list = [initial_states[n].compose(ansatz) for n in range(self.k)] + + self.weight_vector = self._check_weight_vector(self.weight_vector) + + evaluate_weighted_energy_sum = self._get_evaluate_weighted_energy_sum( + initialized_ansatz_list, operator + ) + + if self.gradient is not None: # need to implement _get_evaluate_gradient + evaluate_gradient = self._get_evalute_gradient(initialized_ansatz_list, operator) + else: + evaluate_gradient = None + + if aux_operators: + zero_op = PauliSumOp.from_list([("I" * self.ansatz.num_qubits, 0)]) + + # Convert the None and zero values when aux_operators is a list. + # Drop None and convert zero values when aux_operators is a dict. + if isinstance(aux_operators, list): + key_op_iterator = enumerate(aux_operators) + converted = [zero_op] * len(aux_operators) + else: + key_op_iterator = aux_operators.items() + converted = {} + for key, op in key_op_iterator: + if op is not None: + converted[key] = zero_op if op == 0 else op + + aux_operators = converted + + else: + aux_operators = None + + start_time = time() + + if callable(self.optimizer): + optimizer_result = self.optimizer( + fun=evaluate_weighted_energy_sum, + x0=initial_point, + jac=evaluate_gradient, + bounds=bounds, + ) + else: + optimizer_result = self.optimizer.minimize( + fun=evaluate_weighted_energy_sum, + x0=initial_point, + jac=evaluate_gradient, + bounds=bounds, + ) + + optimizer_time = time() - start_time + + logger.info( + "Optimization complete in %s seconds.\nFound opt_params %s", + optimizer_time, + optimizer_result.x, + ) + + if aux_operators is not None: + bound_ansatz_list = [ + initialized_ansatz_list[n].bind_parameters(optimizer_result.x) + for n in range(self.k) + ] + + aux_values_list = [ + estimate_observables( + self.estimator, + bound_ansatz_list[n], + aux_operators, + ) + for n in range(self.k) + ] + else: + aux_values_list = None + + return self._build_ssvqe_result( + optimizer_result, aux_values_list, optimizer_time, operator, initialized_ansatz_list + ) + + def _get_evaluate_weighted_energy_sum( + self, + initialized_ansatz_list: list[QuantumCircuit], + operator: BaseOperator | PauliSumOp, + ) -> tuple[Callable[[np.ndarray], float | list[float]], dict]: + """Returns a function handle to evaluate the weighted energy sum at given parameters + for the ansatz. This is the objective function to be passed to the optimizer + that is used for evaluation. + Args: + initialized_anastz_list: A list consisting of copies of the ansatz initialized + in the initial states. + operator: The operator whose expectation value with respect to each of the + states in ``initialzed_ansatz_list`` is being measured. + Returns: + Weighted expectation value sum of the operator for each parameter. + Raises: + AlgorithmError: If the primitive job to evaluate the weighted energy + sum fails. + """ + num_parameters = initialized_ansatz_list[0].num_parameters + + eval_count = 0 + + def evaluate_weighted_energy_sum(parameters): + nonlocal eval_count + parameters = np.reshape(parameters, (-1, num_parameters)).tolist() + batchsize = len(parameters) + + try: + job = self.estimator.run( + [initialized_ansatz_list[m] for n in range(batchsize) for m in range(self.k)], + [operator] * self.k * batchsize, + [parameters[n] for n in range(batchsize) for m in range(self.k)], + ) + result = job.result() + values = result.values + + energies = np.reshape(values, (batchsize, self.k)) + weighted_energy_sums = np.dot(energies, self.weight_vector).tolist() + energies = energies.tolist() + + except Exception as exc: + raise AlgorithmError("The primitive job to evaluate the energy failed!") from exc + + if self.callback is not None: + metadata = result.metadata + for params, energies_value, metadata in zip(parameters, energies, metadata): + eval_count += 1 + self.callback(eval_count, params, energies_value, metadata) + + return ( + weighted_energy_sums[0] if len(weighted_energy_sums) == 1 else weighted_energy_sums + ) + + return evaluate_weighted_energy_sum + + def _get_evalute_gradient( + self, + initialized_ansatz_list: list[QuantumCircuit], + operator: BaseOperator | PauliSumOp, + ) -> tuple[Callable[[np.ndarray], np.ndarray]]: + """Get a function handle to evaluate the gradient at given parameters for the ansatz. + Args: + initialized_ansatz_list: The list of initialized ansatz preparing the quantum states. + operator: The operator whose energy to evaluate. + Returns: + A function handle to evaluate the gradient at given parameters for the initialized + ansatz list. + Raises: + AlgorithmError: If the primitive job to evaluate the gradient fails. + """ + + def evaluate_gradient(parameters): + # broadcasting not required for the estimator gradients + try: + job = self.gradient.run( + initialized_ansatz_list, + [operator] * self.k, + [parameters for n in range(self.k)], + ) + energy_gradients = job.result().gradients + weighted_energy_sum_gradient = sum( + [self.weight_vector[n] * energy_gradients[n] for n in range(self.k)] + ) + except Exception as exc: + raise AlgorithmError("The primitive job to evaluate the gradient failed!") from exc + + return weighted_energy_sum_gradient + + return evaluate_gradient + + def _check_circuit_num_qubits( + self, operator: BaseOperator | PauliSumOp, circuit: QuantumCircuit, circuit_type: str + ) -> QuantumCircuit: + """Check that the number of qubits for the circuit passed matches + the number of qubits of the operator. + """ + if operator.num_qubits != circuit.num_qubits: + try: + logger.info( + "Trying to resize %s to match operator on %s qubits.", + circuit_type, + operator.num_qubits, + ) + circuit.num_qubits = operator.num_qubits + except AttributeError as error: + raise AlgorithmError( + f"The number of qubits of the {circuit_type} does not match the ", + f"operator, and the {circuit_type} does not allow setting the " + "number of qubits using `num_qubits`.", + ) from error + return circuit + + def _check_operator_ansatz(self, operator: BaseOperator | PauliSumOp) -> QuantumCircuit: + """Check that the number of qubits of operator and ansatz match and that the ansatz is + parameterized. + """ + # set defaults + if self.ansatz is None: + ansatz = RealAmplitudes(num_qubits=operator.num_qubits, reps=6) + else: + ansatz = self.ansatz + + ansatz = self._check_circuit_num_qubits( + operator=operator, circuit=ansatz, circuit_type="ansatz" + ) + + if ansatz.num_parameters == 0: + raise AlgorithmError("The ansatz must be parameterized, but has no free parameters.") + + return ansatz + + def _check_operator_initial_states( + self, list_of_states: list[QuantumCircuit] | None, operator: BaseOperator | PauliSumOp + ) -> QuantumCircuit: + + """Check that the number of qubits of operator and all the initial states match.""" + + if list_of_states is None: + initial_states = [QuantumCircuit(operator.num_qubits) for n in range(self.k)] + for n in range(self.k): + initial_states[n].initialize(Statevector.from_int(n, 2**operator.num_qubits)) + + warnings.warn( + "No initial states have been provided to SSVQE, so they have been set to " + "a subset of the computational basis states. This may result in unphysical " + "results for some chemistry and physics problems." + ) + + else: + initial_states = list_of_states + if self.check_initial_states_orthogonal is True: + stacked_states_array = np.column_stack( + [np.asarray(Statevector(state)) for state in initial_states] + ) + if not np.allclose( + stacked_states_array.transpose() @ stacked_states_array, np.eye(self.k) + ): + print(stacked_states_array.transpose() @ stacked_states_array) + raise AlgorithmError( + "The set of initial states provided is not mutually orthogonal." + ) + + for initial_state in initial_states: + initial_state = self._check_circuit_num_qubits( + operator=operator, circuit=initial_state, circuit_type="initial state" + ) + + return initial_states + + def _check_weight_vector(self, weight_vector: Sequence[float]) -> Sequence[float]: + """Check that the number of weights matches the number of states.""" + if weight_vector is None: + weight_vector = [self.k - n for n in range(self.k)] + elif len(weight_vector) != self.k: + raise AlgorithmError( + "The number of weights provided does not match the number of states." + ) + + return weight_vector + + def _eval_aux_ops( + self, + ansatz: QuantumCircuit, + aux_operators: ListOrDict[BaseOperator | PauliSumOp], + ) -> ListOrDict[tuple(complex, complex)]: + """Compute auxiliary operator eigenvalues.""" + + if isinstance(aux_operators, dict): + aux_ops = list(aux_operators.values()) + else: + aux_ops = aux_operators + + num_aux_ops = len(aux_ops) + + try: + aux_job = self.estimator.run([ansatz] * num_aux_ops, aux_ops) + aux_values = aux_job.result().values + aux_values = list(zip(aux_values, [0] * len(aux_values))) + + if isinstance(aux_operators, dict): + aux_values = dict(zip(aux_operators.keys(), aux_values)) + + except Exception as exc: + raise AlgorithmError( + "The primitive job to evaluate the aux operator values failed!" + ) from exc + + return aux_values + + def _build_ssvqe_result( + self, + optimizer_result: OptimizerResult, + aux_operators_evaluated: ListOrDict[tuple[complex, tuple[complex, int]]], + optimizer_time: float, + operator: BaseOperator | PauliSumOp, + initialized_ansatz_list: list[QuantumCircuit], + ) -> SSVQEResult: + result = SSVQEResult() + + try: + result.eigenvalues = ( + self.estimator.run( + initialized_ansatz_list, [operator] * self.k, [optimizer_result.x] * self.k + ) + .result() + .values + ) + + except Exception as exc: + raise AlgorithmError("The primitive job to evaluate the eigenvalues failed!") from exc + + result.cost_function_evals = optimizer_result.nfev + result.optimal_point = optimizer_result.x + result.optimal_parameters = dict(zip(self.ansatz.parameters, optimizer_result.x)) + result.optimal_value = optimizer_result.fun + result.optimizer_time = optimizer_time + result.aux_operators_evaluated = aux_operators_evaluated + result.optimizer_result = optimizer_result + + return result + + +class SSVQEResult(VariationalResult, EigensolverResult): + """SSVQE Result.""" + + def __init__(self) -> None: + super().__init__() + self._cost_function_evals = None + + @property + def cost_function_evals(self) -> int: + """Returns number of cost optimizer evaluations""" + return self._cost_function_evals + + @cost_function_evals.setter + def cost_function_evals(self, value: int) -> None: + """Sets number of cost function evaluations""" + self._cost_function_evals = value + + +def _validate_initial_point(point, ansatz): + expected_size = ansatz.num_parameters + + # try getting the initial point from the ansatz + if point is None and hasattr(ansatz, "preferred_init_points"): + point = ansatz.preferred_init_points + # if the point is None choose a random initial point + + if point is None: + # get bounds if ansatz has them set, otherwise use [-2pi, 2pi] for each parameter + bounds = getattr(ansatz, "parameter_bounds", None) + if bounds is None: + bounds = [(-2 * np.pi, 2 * np.pi)] * expected_size + + # replace all Nones by [-2pi, 2pi] + lower_bounds = [] + upper_bounds = [] + for lower, upper in bounds: + lower_bounds.append(lower if lower is not None else -2 * np.pi) + upper_bounds.append(upper if upper is not None else 2 * np.pi) + + # sample from within bounds + point = algorithm_globals.random.uniform(lower_bounds, upper_bounds) + + elif len(point) != expected_size: + raise ValueError( + f"The dimension of the initial point ({len(point)}) does not match the " + f"number of parameters in the circuit ({expected_size})." + ) + + return point + + +def _validate_bounds(ansatz): + if hasattr(ansatz, "parameter_bounds") and ansatz.parameter_bounds is not None: + bounds = ansatz.parameter_bounds + if len(bounds) != ansatz.num_parameters: + raise ValueError( + f"The number of bounds ({len(bounds)}) does not match the number of " + f"parameters in the circuit ({ansatz.num_parameters})." + ) + else: + bounds = [(None, None)] * ansatz.num_parameters + + return bounds diff --git a/electronic_structure_algorithms/excited_states_eigensolvers/vqd.py b/electronic_structure_algorithms/excited_states_eigensolvers/vqd.py new file mode 100644 index 0000000..f7e293a --- /dev/null +++ b/electronic_structure_algorithms/excited_states_eigensolvers/vqd.py @@ -0,0 +1,537 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2022. +# +# 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 Variational Quantum Deflation Algorithm for computing higher energy states. + +See https://arxiv.org/abs/1805.08138. This implementation is derived from Qiskit's +implementation of VQD and has been modified to be able to handle the use +of different ansatzes for each excited state. +""" + +from __future__ import annotations + +from collections.abc import Callable, Sequence +from typing import Optional, Union, List, Any +import logging +from time import time + +import numpy as np + +from qiskit.circuit import QuantumCircuit +from qiskit.quantum_info.operators.base_operator import BaseOperator +from qiskit.opflow import PauliSumOp +from qiskit.primitives import BaseEstimator +from qiskit.algorithms.state_fidelities import BaseStateFidelity + +from qiskit.algorithms.list_or_dict import ListOrDict +from qiskit.algorithms.optimizers import Optimizer, Minimizer, OptimizerResult +from qiskit.algorithms.variational_algorithm import VariationalAlgorithm +from qiskit.algorithms.eigensolvers.eigensolver import Eigensolver, EigensolverResult +from qiskit.algorithms.utils import validate_bounds, validate_initial_point +from qiskit.algorithms.exceptions import AlgorithmError +from qiskit.algorithms.observables_evaluator import estimate_observables +from qiskit.algorithms.utils.set_batching import _set_default_batchsize + +logger = logging.getLogger(__name__) + + +class VQD(VariationalAlgorithm, Eigensolver): + r"""The Variational Quantum Deflation algorithm. Implementation using primitives. + + `VQD `__ is a quantum algorithm that uses a + variational technique to find + the k eigenvalues of the Hamiltonian :math:`H` of a given system. + + The algorithm computes excited state energies of generalised hamiltonians + by optimising over a modified cost function where each succesive eigenvalue + is calculated iteratively by introducing an overlap term with all + the previously computed eigenstates that must be minimised, thus ensuring + higher energy eigenstates are found. + + An instance of VQD requires defining three algorithmic sub-components: + an integer k denoting the number of eigenstates to calculate, a trial + state (a.k.a. ansatz) which is a :class:`QuantumCircuit` or a list of + :class:`QuantumCircuit`, and the classical :mod:`~qiskit.algorithms.optimizers`. + The optimizer varies the circuit parameters + The trial state :math:`|\psi(\vec\theta)\rangle` is varied by the optimizer, + which modifies the set of ansatz parameters :math:`\vec\theta` + such that the expectation value of the operator on the corresponding + state approaches a minimum. The algorithm does this by iteratively refining + each excited state to be orthogonal to all the previous excited states. + + An optional array of parameter values, via the *initial_point*, may be provided + as the starting point for the search of the minimum eigenvalue. This feature is + particularly useful when there are reasons to believe that the solution point + is close to a particular point. + + The length of the *initial_point* list value must match the number of the parameters + expected by the ansatz. If the *initial_point* is left at the default + of ``None``, then VQD will look to the ansatz for a preferred value, based on its + given initial state. If the ansatz returns ``None``, + then a random point will be generated within the parameter bounds set, as per above. + It is also possible to give a list of initial points, one for every kth eigenvalue. + If the ansatz provides ``None`` as the lower bound, then VQD + will default it to :math:`-2\pi`; similarly, if the ansatz returns ``None`` + as the upper bound, the default value will be :math:`2\pi`. + + The following attributes can be set via the initializer but can also be read and + updated once the VQD object has been constructed. + + Attributes: + estimator (BaseEstimator): The primitive instance used to perform the expectation + estimation as indicated in the VQD paper. + fidelity (BaseStateFidelity): The fidelity class instance used to compute the + overlap estimation as indicated in the VQD paper. + ansatz (QuantumCircuit): A parameterized circuit used as ansatz for the wave function. + optimizer(Optimizer | list[Optimizer]): A classical optimizer or a list of optimizers, + one for every k-th eigenvalue. Can either be a Qiskit optimizer or a callable + that takes an array as input and returns a Qiskit or SciPy optimization result. + k (int): the number of eigenvalues to return. Returns the lowest k eigenvalues. + betas (list[float]): Beta parameters in the VQD paper. + Should have length k - 1, with k the number of excited states. + These hyper-parameters balance the contribution of each overlap term to the cost + function and have a default value computed as the mean square sum of the + coefficients of the observable. + initial point (list[float] | list[list[float]]): An optional initial point (i.e. + initial parameter values) for the optimizer or a list of initial points, one + for every k-th eigenvalue. If ``None`` then VQD will look to the ansatz for a + preferred point and if not will simply compute a random one. + callback (Callable[[int, np.ndarray, float, dict[str, Any]], None] | None): + A callback that can access the intermediate data + during the optimization. Four parameter values are passed to the callback as + follows during each evaluation by the optimizer: the evaluation count, + the optimizer parameters for the ansatz, the estimated value, the estimation + metadata, and the current step. + """ + + def __init__( + self, + estimator: BaseEstimator, + fidelity: BaseStateFidelity, + ansatz: QuantumCircuit | list[QuantumCircuit], + optimizer: Optimizer | Minimizer | Sequence[Optimizer | Minimizer], + *, + k: int = 2, + betas: Sequence[float] | None = None, + initial_point: Sequence[float] | Sequence[Sequence[float]] | None = None, + callback: Callable[[int, np.ndarray, float, dict[str, Any]], None] | None = None, + ) -> None: + """ + + Args: + estimator: The estimator primitive. + fidelity: The fidelity class using primitives. + ansatz: A parameterized circuit used as ansatz for the wave function. + optimizer: A classical optimizer or a list of optimizers, one for every k-th eigenvalue. + Can either be a Qiskit optimizer or a callable + that takes an array as input and returns a Qiskit or SciPy optimization result. + k: The number of eigenvalues to return. Returns the lowest k eigenvalues. + betas: Beta parameters in the VQD paper. + Should have length k - 1, with k the number of excited states. + These hyperparameters balance the contribution of each overlap term to the cost + function and have a default value computed as the mean square sum of the + coefficients of the observable. + initial_point: An optional initial point (i.e. initial parameter values) + for the optimizer or a list of initial points, one for every k-th eigenvalue. + If ``None`` then VQD will look to the ansatz for a preferred + point and if not will simply compute a random one. + callback: A callback that can access the intermediate data + during the optimization. Four parameter values are passed to the callback as + follows during each evaluation by the optimizer: the evaluation count, + the optimizer parameters for the ansatz, the estimated value, + the estimation metadata, and the current step. + """ + super().__init__() + + self.estimator = estimator + self.fidelity = fidelity + self.ansatz = ansatz + self.optimizer = optimizer + self.k = k + self.betas = betas + self.initial_point = initial_point + self.callback = callback + + self._eval_count = 0 + + @property + def initial_point(self) -> Optional[Union[np.ndarray, List[np.ndarray]]]: + """Returns initial point.""" + return self._initial_point + + @initial_point.setter + def initial_point(self, initial_point: Optional[Union[np.ndarray, List[np.ndarray]]]): + """Sets initial point""" + self._initial_point = initial_point + + def _check_operator_ansatz(self, operator: BaseOperator | PauliSumOp): + """Check that the number of qubits of operator and ansatz match.""" + if operator is not None and self.ansatz is not None: + for ansatz in self.ansatz: + if operator.num_qubits != ansatz.num_qubits: + # try to set the number of qubits on the ansatz, if possible + try: + ansatz.num_qubits = operator.num_qubits + except AttributeError as ex: + raise AlgorithmError( + "The number of qubits of the ansatz does not match the " + "operator, and the ansatz does not allow setting the " + "number of qubits using `num_qubits`." + ) from ex + + @classmethod + def supports_aux_operators(cls) -> bool: + return True + + def compute_eigenvalues( + self, + operator: BaseOperator | PauliSumOp, + aux_operators: ListOrDict[BaseOperator | PauliSumOp] | None = None, + ) -> VQDResult: + + super().compute_eigenvalues(operator, aux_operators) + + if isinstance(self.ansatz, QuantumCircuit): + self.ansatz = [self.ansatz.copy()]*self.k + elif len(self.ansatz) != self.k: + raise AlgorithmError("The number of ansatz circuits provided does not match k, the number of states.") + # this sets the size of the ansatz, so it must be called before the initial point + # validation + self._check_operator_ansatz(operator) + + bounds = [validate_bounds(ansatz) for ansatz in self.ansatz] + + # We need to handle the array entries being zero or Optional i.e. having value None + if aux_operators: + zero_op = PauliSumOp.from_list([("I" * self.ansatz.num_qubits, 0)]) + + # Convert the None and zero values when aux_operators is a list. + # Drop None and convert zero values when aux_operators is a dict. + if isinstance(aux_operators, list): + key_op_iterator = enumerate(aux_operators) + converted = [zero_op] * len(aux_operators) + else: + key_op_iterator = aux_operators.items() + converted = {} + for key, op in key_op_iterator: + if op is not None: + converted[key] = zero_op if op == 0 else op + + aux_operators = converted + + else: + aux_operators = None + + if self.betas is None: + if isinstance(operator, PauliSumOp): + upper_bound = abs(operator.coeff) * sum( + abs(operation.coeff) for operation in operator + ) + betas = [upper_bound * 10] * (self.k) + logger.info("beta autoevaluated to %s", betas[0]) + else: + raise NotImplementedError( + r"Beta autoevaluation is only supported for operators" + f"of type PauliSumOp, found {type(operator)}." + ) + else: + betas = self.betas + + result = self._build_vqd_result() + + if aux_operators is not None: + aux_values = [] + + # We keep a list of the bound circuits with optimal parameters, to avoid re-binding + # the same parameters to the ansatz if we do multiple steps + prev_states = [] + + num_initial_points = 0 + if self.initial_point is not None: + if isinstance(self.initial_point[0], np.ndarray) or isinstance(self.initial_point[0], list): + initial_points = self.initial_point + else: + initial_points = [self.initial_point] + + num_initial_points = len(initial_points) + + # 0 just means the initial point is ``None`` and ``validate_initial_point`` + # will select a random point + #TODO: determine how to validate initial point for case where self.ansatz + # is a list, but num_initial_points <= 1. + if num_initial_points <= 1: + initial_point = validate_initial_point(self.initial_point, self.ansatz) + + for step in range(1, self.k + 1): + if num_initial_points > 1: + initial_point = validate_initial_point(initial_points[step - 1], self.ansatz[step - 1]) + + if step > 1: + prev_states.append(self.ansatz[step - 2].bind_parameters(result.optimal_points[-1])) + + self._eval_count = 0 + energy_evaluation = self._get_evaluate_energy( + step, operator, betas, prev_states=prev_states + ) + + start_time = time() + + # TODO: add gradient support after FidelityGradients are implemented + if isinstance(self.optimizer, list): + optimizer = self.optimizer[step - 1] + else: + optimizer = self.optimizer # fall back to single optimizer if not list + + if callable(optimizer): + opt_result = optimizer( # pylint: disable=not-callable + fun=energy_evaluation, x0=initial_point, bounds=bounds[step - 1] + ) + else: + # we always want to submit as many estimations per job as possible for minimal + # overhead on the hardware + was_updated = _set_default_batchsize(optimizer) + + opt_result = optimizer.minimize( + fun=energy_evaluation, x0=initial_point, bounds=bounds[step - 1] + ) + + # reset to original value + if was_updated: + optimizer.set_max_evals_grouped(None) + + eval_time = time() - start_time + + self._update_vqd_result(result, opt_result, eval_time, self.ansatz[step - 1].copy()) + + if aux_operators is not None: + aux_value = estimate_observables( + self.estimator, self.ansatz[step - 1], aux_operators, result.optimal_points[-1] + ) + aux_values.append(aux_value) + + if step == 1: + logger.info( + "Ground state optimization complete in %s seconds.\n" + "Found opt_params %s in %s evals", + eval_time, + result.optimal_points, + self._eval_count, + ) + else: + logger.info( + ( + "%s excited state optimization complete in %s s.\n" + "Found opt_params %s in %s evals" + ), + str(step - 1), + eval_time, + result.optimal_points, + self._eval_count, + ) + + # To match the signature of EigensolverResult + result.eigenvalues = np.array(result.eigenvalues) + result.optimal_points = np.array(result.optimal_points) + result.optimal_values = np.array(result.optimal_values) + result.cost_function_evals = np.array(result.cost_function_evals) + result.optimizer_times = np.array(result.optimizer_times) + + if aux_operators is not None: + result.aux_operators_evaluated = aux_values + + return result + + def _get_evaluate_energy( + self, + step: int, + operator: BaseOperator | PauliSumOp, + betas: Sequence[float], + prev_states: list[QuantumCircuit] | None = None, + ) -> Callable[[np.ndarray], float | list[float]]: + """Returns a function handle to evaluate the ansatz's energy for any given parameters. + This is the objective function to be passed to the optimizer that is used for evaluation. + + Args: + step: level of energy being calculated. 0 for ground, 1 for first excited state... + operator: The operator whose energy to evaluate. + prev_states: List of optimal circuits from previous rounds of optimization. + + Returns: + A callable that computes and returns the energy of the hamiltonian + of each parameter. + + Raises: + AlgorithmError: If the circuit is not parameterized (i.e. has 0 free parameters). + AlgorithmError: If operator was not provided. + RuntimeError: If the previous states array is of the wrong size. + """ + + num_parameters = self.ansatz[step - 1].num_parameters + if num_parameters == 0: + raise AlgorithmError("The ansatz must be parameterized, but has no free parameters.") + + if step > 1 and (len(prev_states) + 1) != step: + raise RuntimeError( + f"Passed previous states of the wrong size." + f"Passed array has length {str(len(prev_states))}" + ) + + self._check_operator_ansatz(operator) + + def evaluate_energy(parameters: np.ndarray) -> np.ndarray | float: + + try: + estimator_job = self.estimator.run( + circuits=[self.ansatz[step - 1]], observables=[operator], parameter_values=[parameters] + ) + estimator_result = estimator_job.result() + values = estimator_result.values + + except Exception as exc: + raise AlgorithmError("The primitive job to evaluate the energy failed!") from exc + + costs = None + if step > 1: + # Compute overlap cost + fidelity_job = self.fidelity.run( + [self.ansatz[step - 1]] * (step - 1), + prev_states, + [parameters] * (step - 1), + ) + + costs = fidelity_job.result().fidelities + + for (state, cost) in zip(range(step - 1), costs): + values += np.real(betas[state] * cost) + + if self.callback is not None: + metadata = estimator_result.metadata + for params, value, meta in zip([parameters], values, metadata): + self._eval_count += 1 + self.callback(self._eval_count, params, value, meta, step, costs) + else: + self._eval_count += len(values) + + return values if len(values) > 1 else values[0] + + return evaluate_energy + + @staticmethod + def _build_vqd_result() -> VQDResult: + + result = VQDResult() + result.optimal_points = [] + result.optimal_parameters = [] + result.optimal_values = [] + result.cost_function_evals = [] + result.optimizer_times = [] + result.eigenvalues = [] + result.optimizer_results = [] + result.optimal_circuits = [] + return result + + @staticmethod + def _update_vqd_result(result, opt_result, eval_time, ansatz) -> VQDResult: + + result.optimal_points.append(opt_result.x) + result.optimal_parameters.append(dict(zip(ansatz.parameters, opt_result.x))) + result.optimal_values.append(opt_result.fun) + result.cost_function_evals.append(opt_result.nfev) + result.optimizer_times.append(eval_time) + result.eigenvalues.append(opt_result.fun + 0j) + result.optimizer_results.append(opt_result) + result.optimal_circuits.append(ansatz) + return result + + +class VQDResult(EigensolverResult): + """VQD Result.""" + + def __init__(self) -> None: + super().__init__() + self._cost_function_evals = None + self._optimizer_times = None + self._optimal_values = None + self._optimal_points = None + self._optimal_parameters = None + self._optimizer_results = None + self._optimal_circuits = None + + @property + def cost_function_evals(self) -> Sequence[int] | None: + """Returns number of cost optimizer evaluations""" + return self._cost_function_evals + + @cost_function_evals.setter + def cost_function_evals(self, value: Sequence[int]) -> None: + """Sets number of cost function evaluations""" + self._cost_function_evals = value + + @property + def optimizer_times(self) -> Sequence[float] | None: + """Returns time taken for optimization for each step""" + return self._optimizer_times + + @optimizer_times.setter + def optimizer_times(self, value: Sequence[float]) -> None: + """Sets time taken for optimization for each step""" + self._optimizer_times = value + + @property + def optimal_values(self) -> Sequence[float] | None: + """Returns optimal value for each step""" + return self._optimal_values + + @optimal_values.setter + def optimal_values(self, value: Sequence[float]) -> None: + """Sets optimal values""" + self._optimal_values = value + + @property + def optimal_points(self) -> Sequence[np.ndarray] | None: + """Returns optimal point for each step""" + return self._optimal_points + + @optimal_points.setter + def optimal_points(self, value: Sequence[np.ndarray]) -> None: + """Sets optimal points""" + self._optimal_points = value + + @property + def optimal_parameters(self) -> Sequence[dict] | None: + """Returns the optimal parameters for each step""" + return self._optimal_parameters + + @optimal_parameters.setter + def optimal_parameters(self, value: Sequence[dict]) -> None: + """Sets optimal parameters""" + self._optimal_parameters = value + + @property + def optimizer_results(self) -> Sequence[OptimizerResult] | None: + """Returns the optimizer results for each step""" + return self._optimizer_results + + @optimizer_results.setter + def optimizer_results(self, value: Sequence[OptimizerResult]) -> None: + """Sets optimizer results""" + self._optimizer_results = value + + @property + def optimal_circuits(self) -> list[QuantumCircuit]: + """The optimal circuits. Along with the optimal parameters, + these can be used to retrieve the different eigenstates.""" + return self._optimal_circuits + + @optimal_circuits.setter + def optimal_circuits(self, optimal_circuits: list[QuantumCircuit]) -> None: + self._optimal_circuits = optimal_circuits \ No newline at end of file diff --git a/electronic_structure_algorithms/initializations/HF_permutation_matrix.py b/electronic_structure_algorithms/initializations/HF_permutation_matrix.py new file mode 100644 index 0000000..a1a77f8 --- /dev/null +++ b/electronic_structure_algorithms/initializations/HF_permutation_matrix.py @@ -0,0 +1,14 @@ +import torch + +def get_HF_permutation_matrix(num_original_spin_orbitals: int, + num_spin_orbitals: int): + + num_original_molecular_orbitals = int(num_original_spin_orbitals/2) + num_molecular_orbitals = int(num_spin_orbitals/2) + + initial_partial_unitary_guess = torch.zeros(size=(num_original_molecular_orbitals, num_molecular_orbitals), dtype=torch.float64) + for n in range(int(num_molecular_orbitals)): + initial_partial_unitary_guess[n,n] = 1.0 + + return initial_partial_unitary_guess + diff --git a/electronic_structure_algorithms/initializations/__init__.py b/electronic_structure_algorithms/initializations/__init__.py new file mode 100644 index 0000000..08256a5 --- /dev/null +++ b/electronic_structure_algorithms/initializations/__init__.py @@ -0,0 +1,2 @@ +from .configuration_interaction_states import get_CIS_states, get_CISD_states +from .HF_permutation_matrix import get_HF_permutation_matrix diff --git a/electronic_structure_algorithms/initializations/configuration_interaction_states.py b/electronic_structure_algorithms/initializations/configuration_interaction_states.py new file mode 100644 index 0000000..b1759ca --- /dev/null +++ b/electronic_structure_algorithms/initializations/configuration_interaction_states.py @@ -0,0 +1,530 @@ +import numpy as np +from typing import Optional +from qiskit.quantum_info import Statevector + +def count_mismatches(bitstring1: str, + bitstring2: str) -> int: + + """Determines the number of pairs of orbital occupation mismatches of two + Slater determinants represented by bitstrings. + + Args: + + bitstring1, bitstring2: the two Slater determinant bitstrings. + + Returns: + + The number of pairs of orbital occupation mismatches between + bitstring1 and bitstring2. + + """ + + return int(sum(index1 != index2 for index1, index2 in zip(bitstring1, bitstring2))/2) + +def get_occupation_indices(bitstring: str): + + """Determines the indices corresponding to occupied orbitals + in a Slater determinant represented by a bitstring. + + Args: + + bitstring: A bitstring representing a Slater determinant. + + Returns: + + A list of the indices of the orbitals which are occupied + in the Slater determinant. + + """ + + indices = [] + for n in range(len(bitstring)): + + if bitstring[-(n+1)] == '1': + + indices.append(n) + + return indices + +def get_one_mismatched_orbital_pair(bitstring1: str, + bitstring2: str): + + """Given two bitstrings bitstring1 and bitstring2 which + represent Slater determinants and are known to have two + pairs of mismatched orbital occupancies, return a tuple + where each entry denotes the one orbital index for each + bitstring involved in one of the mismatched pairs. + + Args: + + bitstring1, bitstring2: The two Slater determinant bitstrings. + + Returns: + + A tuple storing the occupied index involved in a mismatched pair + for each bitstring. + + """ + + for n in range(len(bitstring1)): + + if bitstring1[-(n+1)] != bitstring2[-(n+1)]: + + if bitstring1[-(n+1)] == '1': + + first_orbital_occupation = n + + else: + + second_orbital_occuptation = n + + return first_orbital_occupation, second_orbital_occuptation + +def get_two_mismatched_orbital_pairs(bitstring1: str, + bitstring2: str): + + """Given two bitstrings bitstring1 and bitstring2 which represent + Slater determinants and are known to have two mismatched pairs + of mismatched orbital occupancies, return a 4-tuple of the + occupied indices involved in each pair. + + Args: + + bitstring1, bitstring2: The two Slater determinant bitstrings. + + Returns: + + A 4-tuple of the occupied indices involved in each pair. + + """ + + state1_occupied_mismatches = [] + state2_occupied_mismatches = [] + + for n in range(len(bitstring1)): + + if bitstring1[-(n+1)] != bitstring2[-(n+1)]: + + if bitstring1[-(n+1)] == '1': + + state1_occupied_mismatches.append(n) + + else: + + state2_occupied_mismatches.append(n) + + return min(state1_occupied_mismatches), max(state1_occupied_mismatches), min(state2_occupied_mismatches), max(state2_occupied_mismatches) + +def gamma(bitstring: str, + index: int): + + """Given an index and a Slater determinant represented by a bitstring, + calculated the gamma factor representing the parity of the occupation + numbers of orbitals 0 through index - 1. For a parity p, the gamma + factor is given by (-1)^p. + + Args: + + bitstring: The slater determinant bitstring. + index: The orbital index. + + Returns: + + The gamma factor for the bitstring and orbital index. + + """ + + if index == 0: + + bitstring = '' + + else: + + bitstring = bitstring[-index:] + + gamma_factor = 1 + + for _ in bitstring: + + if int(_) == 1: + gamma_factor *= -1 + else: + pass + + return gamma_factor + +def get_CIS_states(one_body_integrals: np.ndarray, + two_body_integrals: np.ndarray, + num_particles: tuple[int, int], + state_representation: Optional[str] = 'sparse', # must be either 'sparse' or 'dense' + truncation_threshold: Optional[float] = 10**-10): + + """Calculates the configuration interaction singles states in either a sparse or dense representation. + + Args: + + one_body_integrals: The one body integrals in the spin-orbital representation. + two_body_integrals: The two body integrals in the spin-orbital representation in + the physicist's notation. + num_particles: A tuple representation the number of alpha and beta electrons. + state_representation: Either 'sparse' or 'dense'. Setting equal to 'sparse' + returns a representation mapping Slater determinant bitstrings to coefficients. + Setting equal to 'dense' returns a dense statevector. + truncation_threshold: The threshold for discarding CI coefficients. + + Returns: + + The CIS states. + + """ + num_spin_orbitals = one_body_integrals.shape[0] + num_molecular_orbitals = int(num_spin_orbitals/2) + num_alpha = num_particles[0] + num_beta = num_particles[1] + + alpha_string = '0'*(num_molecular_orbitals - num_alpha) + '1'*num_alpha + beta_string = '0'*(num_molecular_orbitals - num_beta) + '1'*num_beta + + HF_string = alpha_string + beta_string + HF_alpha_occupied_list = [] + HF_alpha_unoccupied_list = [] + HF_beta_occupied_list = [] + HF_beta_unoccupied_list = [] + + possible_alpha_excitations = [] + possible_beta_excitations = [] + excited_bitstrings = [HF_string] + + for n in range(num_molecular_orbitals): + + if alpha_string[-(n+1)] == '1': + + HF_alpha_occupied_list.append(n + num_molecular_orbitals) + + elif alpha_string[-(n+1)] == '0': + + HF_alpha_unoccupied_list.append(n + num_molecular_orbitals) + + + if beta_string[-(n+1)] == '1': + + HF_beta_occupied_list.append(n) + + elif beta_string[-(n+1)] == '0': + + HF_beta_unoccupied_list.append(n) + + + for occupied_index in HF_alpha_occupied_list: + + for unoccupied_index in HF_alpha_unoccupied_list: + + possible_alpha_excitations.append((unoccupied_index, occupied_index)) + + for occupied_index in HF_beta_occupied_list: + + for unoccupied_index in HF_beta_unoccupied_list: + + possible_beta_excitations.append((unoccupied_index, occupied_index)) + + + for possible_excitation in possible_alpha_excitations + possible_beta_excitations: + + excited_bitstring = '' + for n in range(len(HF_string)): + + if n == possible_excitation[0]: + excited_bitstring = '1' + excited_bitstring + elif n == possible_excitation[1]: + excited_bitstring = '0' + excited_bitstring + else: + excited_bitstring = HF_string[-(n+1)] + excited_bitstring + + excited_bitstrings.append(excited_bitstring) + + def get_matrix_element(i: int, + j: int): + + """Given two indices i and j, compute the Hamiltonian matrix element , + where i and j index Slater determinants. + + Args: + + i, j: The slater determinant indices. + + Returns: + + The Hamiltonian matrix eleent . + """ + + ith_state_occupations = get_occupation_indices(excited_bitstrings[i]) + jth_state_occupations = get_occupation_indices(excited_bitstrings[j]) + + num_mismatches = count_mismatches(excited_bitstrings[i], excited_bitstrings[j]) + + if num_mismatches == 0: + + matrix_element = 0 + + for index in ith_state_occupations: + + matrix_element += one_body_integrals[index, index] + + for index1 in ith_state_occupations: + for index2 in jth_state_occupations: + + matrix_element += two_body_integrals[index1, index2, index1, index2] + matrix_element -= two_body_integrals[index1, index2, index2, index1] + + elif num_mismatches == 1: + + m, p = get_one_mismatched_orbital_pair(excited_bitstrings[i], excited_bitstrings[j]) + matrix_element = one_body_integrals[m,p] + + for n in ith_state_occupations: + + matrix_element += 2*two_body_integrals[m, n, p, n] + matrix_element -= 2*two_body_integrals[m, n, n, p] + + matrix_element *= gamma(excited_bitstrings[i], m)*gamma(excited_bitstrings[j], p) + + elif num_mismatches == 2: + + m, n, p, q = get_two_mismatched_orbital_pairs(excited_bitstrings[i], excited_bitstrings[j]) + + matrix_element = two_body_integrals[m,n,p,q] - two_body_integrals[m,n,q,p] + matrix_element *= 2*gamma(excited_bitstrings[i], m)*gamma(excited_bitstrings[i], n)*gamma(excited_bitstrings[j], p)*gamma(excited_bitstrings[j], q) + + elif num_mismatches > 2: + + matrix_element = 0 + + + return matrix_element + + + CIS_matrix = np.empty(shape=(len(excited_bitstrings), len(excited_bitstrings))) + for x in range(len(excited_bitstrings)): + for y in range(len(excited_bitstrings)): + + CIS_matrix[x,y] = get_matrix_element(x,y) + + CIS_results = np.linalg.eigh(CIS_matrix) + + sparse_results = [dict(zip(excited_bitstrings, CIS_results[1].transpose()[n])) for n in range(len(excited_bitstrings))] + for state in sparse_results: + for key in state: + if np.abs(state[key]) <= truncation_threshold: + state[key] = 0.0 + + if state_representation == 'sparse': + + return sparse_results + + if state_representation =='dense': + + dense_results = [] + + for n, state in enumerate(sparse_results): + + key_counter = 0 + for key in state: + + if key_counter == 0: + + dense_results.append(state[key]*Statevector.from_label(key)) + + else: + + dense_results[n] += state[key]*Statevector.from_label(key) + + key_counter += 1 + + for n, state in enumerate(dense_results): + + if not state.is_valid(): + + norm = np.sqrt(np.sum([elem**2 for elem in state])) + dense_results[n] = Statevector(np.array([elem/norm for elem in state])) + + + return dense_results + + +def get_CISD_states(one_body_integrals: np.ndarray, + two_body_integrals: np.ndarray, + num_particles: tuple[int, int], + num_spin_orbitals: int, + state_representation: Optional[str] = 'sparse', # must be either 'sparse' or 'dense' + truncation_threshold: Optional[float] = 10 ** -10): + + """Calculates the configuration interaction singles and doubles (CISD) states in either a sparse + or dense representation. + + Args: + + one_body_integrals: The one body integrals in the spin-orbital representation. + two_body_integrals: The two body integrals in the spin-orbital representation in + the physicist's notation. + num_particles: A tuple representation the number of alpha and beta electrons. + state_representation: Either 'sparse' or 'dense'. Setting equal to 'sparse' + returns a representation mapping Slater determinant bitstrings to coefficients. + Setting equal to 'dense' returns a dense statevector. + truncation_threshold: The threshold for discarding CI coefficients. + + Returns: + + The CIS states. + + """ + + num_molecular_orbitals = int(num_spin_orbitals / 2) + num_alpha = num_particles[0] + num_beta = num_particles[1] + + alpha_string = '0' * (num_molecular_orbitals - num_alpha) + '1' * num_alpha + beta_string = '0' * (num_molecular_orbitals - num_beta) + '1' * num_beta + + HF_string = alpha_string + beta_string + HF_alpha_occupied_list = [] + HF_alpha_unoccupied_list = [] + HF_beta_occupied_list = [] + HF_beta_unoccupied_list = [] + + excited_bitstrings = [] + + for n in range(num_molecular_orbitals): + + if alpha_string[-(n + 1)] == '1': + + HF_alpha_occupied_list.append(n + num_molecular_orbitals) + + elif alpha_string[-(n + 1)] == '0': + + HF_alpha_unoccupied_list.append(n + num_molecular_orbitals) + + if beta_string[-(n + 1)] == '1': + + HF_beta_occupied_list.append(n) + + elif beta_string[-(n + 1)] == '0': + + HF_beta_unoccupied_list.append(n) + + for n in range(2**num_spin_orbitals): + + bitstr = '0' * (num_spin_orbitals - len(bin(n)[2:])) + bin(n)[2:] + alpha_bitstr, beta_bitstr = bitstr[:num_molecular_orbitals], bitstr[num_molecular_orbitals:] + + if (alpha_bitstr.count('1') == num_alpha and beta_bitstr.count('1') == num_beta): + + num_alpha_excitations, num_beta_excitations = alpha_bitstr[:-num_alpha].count('1'), beta_bitstr[:-num_beta].count('1') + if num_alpha_excitations + num_beta_excitations <= 2: + + excited_bitstrings.append(bitstr) + + + def get_matrix_element(i: int, + j: int): + + """Given two indices i and j, compute the Hamiltonian matrix element , + where i and j index Slater determinants. + + Args: + + i, j: The slater determinant indices. + + Returns: + + The Hamiltonian matrix eleent . + """ + + ith_state_occupations = get_occupation_indices(excited_bitstrings[i]) + jth_state_occupations = get_occupation_indices(excited_bitstrings[j]) + + num_mismatches = count_mismatches(excited_bitstrings[i], excited_bitstrings[j]) + + if num_mismatches == 0: + + matrix_element = 0 + + for index in ith_state_occupations: + matrix_element += one_body_integrals[index, index] + + for index1 in ith_state_occupations: + for index2 in jth_state_occupations: + matrix_element += two_body_integrals[index1, index2, index1, index2] + matrix_element -= two_body_integrals[index1, index2, index2, index1] + + elif num_mismatches == 1: + + m, p = get_one_mismatched_orbital_pair(excited_bitstrings[i], excited_bitstrings[j]) + matrix_element = one_body_integrals[m, p] + + for n in ith_state_occupations: + matrix_element += 2 * two_body_integrals[m, n, p, n] + matrix_element -= 2 * two_body_integrals[m, n, n, p] + + matrix_element *= gamma(excited_bitstrings[i], m) * gamma(excited_bitstrings[j], p) + + elif num_mismatches == 2: + + m, n, p, q = get_two_mismatched_orbital_pairs(excited_bitstrings[i], excited_bitstrings[j]) + + matrix_element = two_body_integrals[m, n, p, q] - two_body_integrals[m, n, q, p] + matrix_element *= 2 * gamma(excited_bitstrings[i], m) * gamma(excited_bitstrings[i], n) * gamma( + excited_bitstrings[j], p) * gamma(excited_bitstrings[j], q) + + elif num_mismatches > 2: + + matrix_element = 0 + + return matrix_element + + CISD_matrix = np.empty(shape=(len(excited_bitstrings), len(excited_bitstrings))) + for x in range(len(excited_bitstrings)): + for y in range(len(excited_bitstrings)): + CISD_matrix[x, y] = get_matrix_element(x, y) + + CISD_results = np.linalg.eigh(CISD_matrix) + + sparse_results = [dict(zip(excited_bitstrings, CISD_results[1].transpose()[n])) for n in + range(len(excited_bitstrings))] + for state in sparse_results: + for key in state: + if np.abs(state[key]) <= truncation_threshold: + state[key] = 0.0 + + if state_representation == 'sparse': + + return sparse_results + + elif state_representation == 'dense': + + dense_results = [] + + for n, state in enumerate(sparse_results): + + key_counter = 0 + for key in state: + + if key_counter == 0: + + dense_results.append(state[key] * Statevector.from_label(key)) + + else: + + dense_results[n] += state[key] * Statevector.from_label(key) + + key_counter += 1 + + for n, state in enumerate(dense_results): + + if not state.is_valid(): + + norm = np.sqrt(np.sum([elem**2 for elem in dense_results[n]])) + dense_results[n] = Statevector(np.array([elem/norm for elem in dense_results[n]])) + + return dense_results + + diff --git a/electronic_structure_algorithms/orbital_optimization/__init__.py b/electronic_structure_algorithms/orbital_optimization/__init__.py new file mode 100644 index 0000000..1861b4d --- /dev/null +++ b/electronic_structure_algorithms/orbital_optimization/__init__.py @@ -0,0 +1,10 @@ +from .opt_orb_vqe import OptOrbVQE, OptOrbVQEResult +from .opt_orb_adapt_vqe import OptOrbAdaptVQE, OptOrbAdaptVQEResult +from .opt_orb_ssvqe import OptOrbSSVQE, OptOrbSSVQEResult +from .opt_orb_mcvqe import OptOrbMCVQE, OptOrbMCVQEResult +from .opt_orb_vqd import OptOrbVQD, OptOrbVQDResult +from .partial_unitary_projection_optimizer import PartialUnitaryProjectionOptimizer + +from .base_opt_orb_solver import BaseOptOrbSolver, BaseOptOrbResult +from .opt_orb_minimum_eigensolver import OptOrbMinimumEigensolver, OptOrbMinimumEigensolverResult +from .opt_orb_eigensolver import OptOrbEigensolver, OptOrbEigensolverResult diff --git a/electronic_structure_algorithms/orbital_optimization/base_opt_orb_solver.py b/electronic_structure_algorithms/orbital_optimization/base_opt_orb_solver.py new file mode 100644 index 0000000..a77628f --- /dev/null +++ b/electronic_structure_algorithms/orbital_optimization/base_opt_orb_solver.py @@ -0,0 +1,657 @@ +from typing import Callable, Optional, Union, Dict, List +from qiskit.primitives import BaseEstimator +from .partial_unitary_projection_optimizer import PartialUnitaryProjectionOptimizer +from qiskit_algorithms.variational_algorithm import VariationalResult +from qiskit.opflow import PauliSumOp, OperatorBase +from qiskit.quantum_info import SparsePauliOp + +from qiskit_nature.second_q.mappers import QubitMapper +from qiskit_nature.second_q.operators import FermionicOp +from qiskit_nature.second_q.hamiltonians import ElectronicEnergy +from qiskit_nature.second_q.problems import ElectronicStructureProblem +from qiskit_nature.second_q.operators.tensor_ordering import to_physicist_ordering + +import torch +import copy +import numpy as np + +class BaseOptOrbSolver(): + + def __init__(self, + problem: Optional[ElectronicStructureProblem], + integral_tensors: Optional[Union[tuple[torch.Tensor, torch.Tensor], tuple[np.ndarray, np.ndarray]]], + num_spin_orbitals: int, + mapper: QubitMapper, + estimator: BaseEstimator, + partial_unitary_optimizer: PartialUnitaryProjectionOptimizer, + initial_partial_unitary: Optional[Union[torch.Tensor, np.ndarray]] = None, + maxiter: int = 10, + stopping_tolerance: float = 10**-5, + spin_conserving: bool = False, + wavefuntion_real: bool = False, + outer_loop_callback: Optional[Callable] = None, + partial_unitary_random_perturbation: Optional[float] = None, + RDM_ops_batchsize: Optional[int] = 100): + + """ + + Args: + problem: An ElectronicStructureProblem from which molecule information such as one and two body integrals is obtained. + integral_tensors: A tuple storing the one and two body integrals of the full orbital space. The first + entry stores the one body integrals and the second entry stores the two body integrals in the Physicist's notation in + the dense spin-orbital representation, represented as either :class:`torch.Tensor` or :class:`np.ndarray`. + num_spin_orbitals: The number of spin-orbitals to use for the active space. + mapper: A QubitMapper to use for the RDM calculations. + partial_unitary_optimizer: An instance of :class:`PartialUnitaryProjectionOptimizer` to use for + the basis optimization subroutine. + initial_partial_unitary: The initial guess for the orbital rotation matrix. If ``None``, then a permutation matrix + selecting the spatial orbitals with the lowest energy will be generated. + maxiter: The maximum number of outerloop iterations. (The number of times the wavefunction optimization is run.) + stopping_tolerance: The stopping tolerance used to determine if the algorithm should be stopped. + spin_conserving: A boolean flag that indicates whether or not we are assuming that spin is conserved + in the system being studied. Setting to ``True`` will skip the calculation of RDM element operators + which do not conserve spin, automatically setting these elements to zero. This should only be + set to ``True`` if it is known that the circuit used to construct the wavefunction in the eigensolver + subroutine conserves the spin-magnetization quantum number (e.g. UCCSD). + wavefunction_real: A boolean flag that indicates whether or not we are assuming that the wavefunction is real. + This allows for the exploitation of symmetries + in the RDM calculation to reduce runtime, but should only be set to ``True`` + if it is known that the circuit used in the eigensolver subroutine produces wavefunction coefficients + which are real-valued. + outer_loop_callback: A callback function that tracks the outerloop progress. + It takes the outerloop iteration, latest eigensolver results, and latest outer loop results as arguments. + This can be used to save the most recent results to a file after each outer loop iteration. + partial_unitary_random_perturbation: A float representing the standard deviation of a normal distribution from which + the matrix elements of a perturbation matrix will be sampled. This perturbation matrix will be added to the initial + partial unitary (and re-orthonormalized) at the beginning of each orbital optimization run. + RDM_ops_batchsize: The number of Pauli ops to store in an :class:`Estimator` at any given time before deleting + the Estimator and replacing it with a blank copy. Increasing this number will lead to increased memory consumption. Setting + this number to be too low will hinder runtime performance. + + """ + + self.mapper = mapper + + # generate copies of the PartialUnitaryProjectionOptimizer instance to use for every outerloop iteration. + self._partial_unitary_optimizer_list = [copy.deepcopy(partial_unitary_optimizer) for n in range(int(maxiter+1))] + + if integral_tensors is not None: + + if isinstance(integral_tensors[0], torch.Tensor) and isinstance(integral_tensors[1], torch.Tensor): + + self.one_body_integrals, self.two_body_integrals = integral_tensors + + elif isinstance(integral_tensors[0], np.ndarray) and isinstance(integral_tensors[1], np.ndarray): + + self.one_body_integrals, self.two_body_integrals = torch.from_numpy(integral_tensors[0]), torch.from_numpy(integral_tensors[1]) + + elif problem is not None: + + self.one_body_integrals = torch.from_numpy(np.asarray(problem.hamiltonian.electronic_integrals.second_q_coeffs()["+-"].to_dense())) + self.two_body_integrals = torch.from_numpy(np.asarray(-1*to_physicist_ordering(problem.hamiltonian.electronic_integrals.second_q_coeffs()["++--"].to_dense()))) + del problem + + if initial_partial_unitary is None: + + num_original_spin_orbitals = self.one_body_integrals.size()[0] + num_original_molecular_orbitals = int(num_original_spin_orbitals/2) + num_molecular_orbitals = int(num_spin_orbitals/2) + + initial_partial_unitary_guess = torch.zeros(size=(num_original_molecular_orbitals, num_molecular_orbitals), dtype=torch.float64) + for n in range(int(num_molecular_orbitals)): + initial_partial_unitary_guess[n,n] = 1.0 + + self.initial_partial_unitary = initial_partial_unitary_guess + + else: + + if isinstance(initial_partial_unitary, torch.Tensor): + + self.initial_partial_unitary = initial_partial_unitary + + elif isinstance(initial_partial_unitary, np.ndarray): + + self.initial_partial_unitary = torch.from_numpy(initial_partial_unitary) + + self.estimator = estimator + self.estimator_list = [copy.deepcopy(estimator) for n in range(maxiter)] + self.maxiter = maxiter + self.spin_conserving = spin_conserving + self.wavefunction_real = wavefuntion_real + self.outer_loop_callback = outer_loop_callback + + self._hamiltonian = None + self.stopping_tolerance = stopping_tolerance + self._current_partial_unitary = self.initial_partial_unitary + self._pauli_op_dict = None + self.num_spin_orbitals = num_spin_orbitals + self.partial_unitary_random_perturbation = partial_unitary_random_perturbation + self.RDM_ops_batchsize = RDM_ops_batchsize + + @property + def partial_unitary_optimizer_list(self) -> List[PartialUnitaryProjectionOptimizer]: + """Returns the list of partial unitary optimizers used for each outer loop iteration.""" + return self._partial_unitary_optimizer_list + + @partial_unitary_optimizer_list.setter + def partial_unitary_optimizer_list(self, optimizer_list: List[PartialUnitaryProjectionOptimizer]) -> None: + """Sets the list of partial unitary optimizers used for each outer loop iteration.""" + self._partial_unitary_optimizer_list = optimizer_list + + @property + def current_partial_unitary(self) -> torch.Tensor: + """Returns the current basis set rotation partial unitary matrix.""" + return self._current_partial_unitary + + @current_partial_unitary.setter + def current_partial_unitary(self, unitary: torch.Tensor) -> None: + """Sets the current basis set rotation partial unitary matrix.""" + self._current_partial_unitary = unitary + + @property + def hamiltonian(self) -> OperatorBase: + """Returns the Hamiltonian in the current basis.""" + return self._hamiltonian + + @hamiltonian.setter + def hamiltonian(self, op: OperatorBase) -> None: + """Sets the Hamiltonian in the current basis.""" + self._hamiltonian = op + + @property + def pauli_op_dict(self) -> Dict: + """Returns the dictonary containing all of the Pauli string operators + necessary for calculating the RDMs.""" + return self._pauli_op_dict + + @pauli_op_dict.setter + def pauli_op_dict(self, some_dict: Dict) -> None: + """Sets the dictionary containing all of the Pauli string operators + necessary for calculating the RDMs.""" + self._pauli_op_dict = some_dict + + def is_2body_op_spin_conserving(self, P,Q,R,S) -> bool: + + """Determines whether or not the two body fermionic excitation operator + involved in the 2-RDM element indexed by (P,Q,R,S) conserves spin or not. + + Args: + + P,Q,R,S: the RDM element index. + + Returns: + True if fermionic operator conserves spin, False if it does not conserve spin. + + """ + + N = self.num_spin_orbitals + + spin_change = 0 + if 0 <= P <= N/2 - 1: + spin_change += 1 + else: + spin_change -= 1 + + if 0 <= Q <= N/2 - 1: + spin_change += 1 + else: + spin_change -= 1 + + if 0 <= R <= N/2 - 1: + spin_change -= 1 + else: + spin_change += 1 + + if 0 <= S <= N/2 - 1: + spin_change -= 1 + else: + spin_change += 1 + + if spin_change == 0: + return True + else: + return False + + def is_1body_op_spin_conserving(self, P,Q) -> bool: + + """Determines whether or not the one body fermionic excitation operator + involved in the 1-RDM element indexed by (P,Q) conserves spin or not. + + Args: + + P,Q: The index of the 1-RDM. + + Returns: + True if the fermionic operator conserves spin, False if it does not conserve + spin. + + """ + + N = self.num_spin_orbitals + + spin_change = 0 + if 0 <= P <= N/2 - 1: + spin_change += 1 + else: + spin_change -= 1 + + if 0 <= Q <= N/2 - 1: + spin_change -= 1 + else: + spin_change += 1 + + if spin_change == 0: + return True + else: + return False + + def construct_pauli_op_dict(self, mapper: QubitMapper) -> Dict: + + """Constructs a dictionary of all the Pauli string operators necessary for computing the RDMs. + The dictionary key/value pairs are of the form str(pauli_op): pauli_op. The uniqueness of + python dictionary keys ensures that no redundant operator evaluations are done. + + Args: + state: The state with respect to which the RDMs are being calculated. + qubit_converter: The QubitConverter used to map the fermionic excitation + operators to qubit operators. + + Returns: + + The dictionary consisting of all the Pauli string operators necessary for calculating + The RDMs. + + """ + + N = self.num_spin_orbitals + pauli_op_dict = {} + + def oneRDM_add_pauli_ops_to_set(p,q): + + op = mapper.map(FermionicOp(data={f'+_{p} -_{q}': 1.0}, + num_spin_orbitals=N)) + + if op.equiv(op.adjoint()): + + pauli_op_dict[str(op)] = op + + else: + + pauli_string_list = op.to_list() + for op_tuple in pauli_string_list: + pauli_op_dict[str(op_tuple[0])] = PauliSumOp(SparsePauliOp(op_tuple[0])) + + return None + + def twoRDM_add_pauli_ops_to_set(p,q,r,s): + + op = mapper.map(FermionicOp(data={f'+_{p} +_{q} -_{s} -_{r}': 1.0}, + num_spin_orbitals=N)) + + if op.equiv(op.adjoint()): + + pauli_op_dict[str(op)] = op + + else: + + pauli_string_list = op.to_list() + for op_tuple in pauli_string_list: + pauli_op_dict[str(op_tuple[0])] = PauliSumOp(SparsePauliOp(op_tuple[0])) + + return None + + + + # if element is True, then we still need to generate the pauli operators relevant to this fermionic operator. + # if element is False, then we do not need to generate any additional pauli operators for this fermionic operator. + # needs_evaluation_array keeps track of this to avoid redundancy. + oneRDM_needs_evaluation_array = np.full((N,N), fill_value=True, dtype=bool) + twoRDM_needs_evaluation_array = np.full((N,N,N,N),fill_value=True, dtype=bool) + + for p in range(N): + for q in range(N): + for r in range(N): + for s in range(N): + + if twoRDM_needs_evaluation_array[p,q,r,s] == True: + + if p == q or r == s or (self.spin_conserving == True and self.is_2body_op_spin_conserving(p,q,r,s) == False): + + # we do not need to evaluate these operators as these entries will be zero in the 2-RDM. + twoRDM_needs_evaluation_array[p,q,r,s] = False + + else: + + # we add the relevant pauli ops to the set, then we set this entry to False + # so we can set other elements to False to take advantage of symmetries in the 2RDM + # to avoid redundant evaluations. + twoRDM_add_pauli_ops_to_set(p,q,r,s) + twoRDM_needs_evaluation_array[p,q,r,s] = False + + twoRDM_needs_evaluation_array[q,p,r,s] = twoRDM_needs_evaluation_array[p,q,r,s] + twoRDM_needs_evaluation_array[p,q,s,r] = twoRDM_needs_evaluation_array[p,q,r,s] + twoRDM_needs_evaluation_array[q,p,s,r] = twoRDM_needs_evaluation_array[p,q,r,s] + + + + twoRDM_needs_evaluation_array[r,s,p,q] = twoRDM_needs_evaluation_array[p,q,r,s] + twoRDM_needs_evaluation_array[r,s,q,p] = twoRDM_needs_evaluation_array[q,p,r,s] + twoRDM_needs_evaluation_array[s,r,p,q] = twoRDM_needs_evaluation_array[p,q,s,r] + twoRDM_needs_evaluation_array[s,r,q,p] = twoRDM_needs_evaluation_array[q,p,s,r] + + for p in range(N): + for q in range(N): + + if oneRDM_needs_evaluation_array[p,q] == True: + + if self.spin_conserving == True and self.is_1body_op_spin_conserving(p,q) == False: + + # we don't need the pauli operators from these terms because + # they are just zero in the 1-RDM + oneRDM_needs_evaluation_array[p,q] = False + + else: + + oneRDM_add_pauli_ops_to_set(p,q) + oneRDM_needs_evaluation_array[p,q] = False + oneRDM_needs_evaluation_array[q,p] = oneRDM_needs_evaluation_array[p,q] + + return pauli_op_dict + + def get_two_RDM_tensor(self, expectval_dict: dict, mapper: QubitMapper) -> torch.Tensor: + + """Constructs and returns the 2-RDM tensor. The class attribute pauli_ops_expectation_values_dict stores the expectation values + of all the Pauli operators necessary for this calculation. get_two_RDM_tensor simply retrieves these expectation values + and constructs the 2-RDM. + + Args: + + state: The state with respect to which the 2-RDM is being calculated. + qubit_converter: The QubitConverter used for mapping fermionic operators to qubit operators. + + Returns: + The 2-RDM with respect to the given state. + + """ + + #N = state.num_qubits # Associating this number with the number of qubits could fail if symmetry reductions for some mappings is used. + N = self.num_spin_orbitals # This should change if this is something we want to do eventually. + global two_RDM_found_complex_value_flag + two_RDM_found_complex_value_flag = False + def get_two_RDM_element(p: int,q: int,r: int,s: int, mapper: QubitMapper): + + global two_RDM_found_complex_value_flag + + op = mapper.map(FermionicOp(data={f'+_{p} +_{q} -_{s} -_{r}': 1.0}, + num_spin_orbitals=N)) + + if op.equiv(op.adjoint()): + + mean = expectval_dict[str(op)] + else: + pauli_string_list = op.to_list() + mean = 0 + for op_tuple in pauli_string_list: + + mean += op_tuple[1]*expectval_dict[str(op_tuple[0])] + + if not np.isclose(np.imag(mean), 0.0, rtol=10e-12, atol=10e-12): + + two_RDM_found_complex_value_flag = True + + if self.wavefunction_real == True: + return np.real(mean) + else: + return mean + + if self.wavefunction_real == True: + tensor = np.full(shape=(N,N,N,N), fill_value=None, dtype=np.float64) + else: + tensor = np.full(shape=(N,N,N,N), fill_value=None, dtype=np.complex128) + + for p in range(N): + for q in range(N): + for r in range(N): + for s in range(N): + + if np.isnan(tensor[p,q,r,s]): + + if p == q or r == s or (self.spin_conserving == True and self.is_2body_op_spin_conserving(p,q,r,s) == False): + + tensor[p,q,r,s] = 0 + + else: + + tensor[p,q,r,s] = get_two_RDM_element(p=p,q=q,r=r,s=s, mapper=mapper) + + tensor[q,p,r,s] = -1*tensor[p,q,r,s] + tensor[p,q,s,r] = -1*tensor[p,q,r,s] + tensor[q,p,s,r] = tensor[p,q,r,s] + + if self.wavefunction_real == True: + + tensor[r,s,p,q] = tensor[p,q,r,s] + tensor[r,s,q,p] = tensor[q,p,r,s] + tensor[s,r,p,q] = tensor[p,q,s,r] + tensor[s,r,q,p] = tensor[q,p,s,r] + + elif self.wavefunction_real == False: + + tensor[r,s,p,q] = np.conj(tensor[p,q,r,s]) + tensor[r,s,q,p] = np.conj(tensor[q,p,r,s]) + tensor[s,r,p,q] = np.conj(tensor[p,q,s,r]) + tensor[s,r,q,p] = np.conj(tensor[q,p,s,r]) + + if two_RDM_found_complex_value_flag == False: + + tensor = np.real(tensor) + + tensor = torch.from_numpy(tensor) + tensor.requires_grad = False + + return tensor + + def get_one_RDM_tensor(self, expectval_dict: dict, mapper: QubitMapper) -> torch.Tensor: + + """Constructs and returns the 1-RDM tensor. The class attribute pauli_ops_expectation_values_dict stores the expectation values + of all the Pauli operators necessary for this calculation. get_one_RDM_tensor simply retrieves these expectation values + as needed and constructs the 1-RDM. + + Args: + state: The state with respect to which the 1-RDM is being calculated. + qubit_converter: The QubitConverter used to map fermionic operators to qubit operators. + + Returns: + The 1-RDM tensor. + + """ + + global one_RDM_found_complex_value_flag + one_RDM_found_complex_value_flag = False + N = self.num_spin_orbitals + def get_one_RDM_element(p: int,q: int, mapper: QubitMapper) -> torch.Tensor: + + global one_RDM_found_complex_value_flag + + op = mapper.map(FermionicOp(data={f'+_{p} -_{q}': 1.0}, + num_spin_orbitals=N)) + + if op.equiv(op.adjoint()): + + mean = expectval_dict[str(op)] + + else: + pauli_string_list = op.to_list() + mean = 0 + for op_tuple in pauli_string_list: + + mean += op_tuple[1]*expectval_dict[str(op_tuple[0])] + + if not np.isclose(np.imag(mean), 0.0, rtol=10e-12, atol=10e-12): + + one_RDM_found_complex_value_flag = True + + if self.wavefunction_real == True: + return np.real(mean) + else: + return mean + + if self.wavefunction_real == True: + tensor = np.full(shape=(N,N), fill_value=None, dtype=np.float64) + else: + tensor = np.full(shape=(N,N), fill_value=None, dtype=np.complex128) + + for p in range(N): + for q in range(N): + + if np.isnan(tensor[p,q]): + + if self.spin_conserving == True and self.is_1body_op_spin_conserving(p,q) == False: + + tensor[p,q] = 0 + + else: + + tensor[p,q] = get_one_RDM_element(p=p,q=q, mapper=mapper) + + if self.wavefunction_real == True: + + tensor[q,p] = tensor[p,q] + + else: + + tensor[q,p] = np.conj(tensor[p,q]) + + if one_RDM_found_complex_value_flag == False: + tensor = np.real(tensor) + + tensor = torch.from_numpy(tensor) + tensor.requires_grad = False + + return tensor + + def compute_rotated_energy(self, partial_unitary: torch.Tensor, + oneRDM: torch.Tensor, + twoRDM: torch.Tensor, + one_body_integrals: torch.Tensor, + two_body_integrals: torch.Tensor) -> float: + """ + Calculates the energy functional with varied U, but fixed wavefunction. + + Args: + partial_unitary: The partial unitary matrix U. + + Returns: + P(U), the energy functional for a given rotation U. + """ + + partial_unitary = torch.block_diag(partial_unitary, partial_unitary) + + + if self.wavefunction_real == True or (oneRDM.dtype !=torch.complex128 and twoRDM.dtype != torch.complex128): + + energy = torch.einsum('pq,pi,qj,ij', one_body_integrals, + partial_unitary, + partial_unitary, + oneRDM) + energy += torch.einsum('pqrs,pi,qj,rk,sl,ijkl', two_body_integrals, + partial_unitary, + partial_unitary, + partial_unitary, + partial_unitary, + twoRDM) + + else: + + partial_unitary = partial_unitary.cdouble() + temp_one_body_integrals = one_body_integrals.cdouble() + temp_two_body_integrals = two_body_integrals.cdouble() + + energy = torch.einsum('pq,pi,qj,ij', temp_one_body_integrals, + partial_unitary, + partial_unitary, + oneRDM) + energy -= torch.einsum('pqrs,pi,qj,rk,sl,ijkl', temp_two_body_integrals, + partial_unitary, + partial_unitary, + partial_unitary, + partial_unitary, + twoRDM) + + return np.real(energy) + + def get_rotated_hamiltonian(self, partial_unitary: torch.Tensor) -> OperatorBase: + + """Transforms the one and two body integrals from the initial larger basis and transforms them according to + a partial unitary matrix U. The transformed Hamiltonian is then constructed from these new integrals. + + Args: + partial_unitary: The partial unitary transformation U. + + Returns: + The transformed Hamiltonian. + + """ + + partial_unitary = torch.block_diag(partial_unitary, partial_unitary) + + rotated_one_body_integrals = torch.einsum('pq,pi,qj->ij', + self.one_body_integrals, + partial_unitary, partial_unitary) + rotated_two_body_integrals = torch.einsum('pqrs,pi,qj,rk,sl->ijkl', + self.two_body_integrals, + partial_unitary, partial_unitary, partial_unitary, partial_unitary) + + num_MO = int(self.num_spin_orbitals/2) + electronic_energy_from_ints = ElectronicEnergy.from_raw_integrals(h1_a=rotated_one_body_integrals.detach().numpy()[0:num_MO, 0:num_MO], + h2_aa=-2*rotated_two_body_integrals.detach().numpy()[0:num_MO, 0:num_MO, 0:num_MO, 0:num_MO]) + + fermionic_op = electronic_energy_from_ints.second_q_op().normal_order() + + return self.mapper.map(fermionic_op) + + def orth(self, V: torch.Tensor) -> torch.Tensor: + """ + Generate the orthonormal projection of the matrix V. + + Args: + V: The matrix to be orthonormalized. + + Returns: + orth(V), the orthogonal projection of the matrix V. + """ + L, Q = torch.linalg.eigh(torch.t(V) @ V) + result = V @ Q @ (torch.float_power(torch.inverse(torch.diag(L)), 0.5)) @ torch.t(Q).double() + return result + +class BaseOptOrbResult(VariationalResult): + + def __init__(self) -> None: + super().__init__() + self._num_vqe_evaluations = 0 + self._optimal_partial_unitary = None + + @property + def num_vqe_evaluations(self) -> int: + """Returns the number of times VQE was run in OptOrbVQE.""" + return self._num_vqe_evaluations + + @num_vqe_evaluations.setter + def num_vqe_evaluations(self, some_int: int) -> None: + """Sets the number of times VQE was run in OptOrbVQE.""" + self._num_vqe_evaluations = some_int + + @property + def optimal_partial_unitary(self) -> torch.Tensor: + """Returns the optimal partial unitary basis transformation.""" + return self._optimal_partial_unitary + + @optimal_partial_unitary.setter + def optimal_partial_unitary(self, some_tensor: torch.Tensor) -> None: + """Sets the optimal partial unitary basis transformation.""" + self._optimal_partial_unitary = some_tensor + + + + + \ No newline at end of file diff --git a/electronic_structure_algorithms/orbital_optimization/opt_orb_adapt_vqe.py b/electronic_structure_algorithms/orbital_optimization/opt_orb_adapt_vqe.py new file mode 100644 index 0000000..67bc327 --- /dev/null +++ b/electronic_structure_algorithms/orbital_optimization/opt_orb_adapt_vqe.py @@ -0,0 +1,95 @@ +from typing import Callable, Optional, Union +from qiskit.primitives import BaseEstimator +from qiskit_algorithms.minimum_eigensolvers import MinimumEigensolver, MinimumEigensolverResult +from qiskit_algorithms.exceptions import AlgorithmError +from qiskit_nature.second_q.mappers import QubitMapper +from qiskit_nature.second_q.problems import ElectronicStructureProblem +import torch +import numpy as np + +from .opt_orb_minimum_eigensolver import OptOrbMinimumEigensolver, OptOrbMinimumEigensolverResult +from .partial_unitary_projection_optimizer import PartialUnitaryProjectionOptimizer + +class OptOrbAdaptVQE(OptOrbMinimumEigensolver): + + def __init__(self, + problem: Optional[ElectronicStructureProblem], + integral_tensors: Optional[Union[tuple[torch.Tensor, torch.Tensor], tuple[np.ndarray, np.ndarray]]], + num_spin_orbitals: int, + ground_state_solver: MinimumEigensolver, + mapper: QubitMapper, + estimator: BaseEstimator, + partial_unitary_optimizer: PartialUnitaryProjectionOptimizer, + initial_partial_unitary: Optional[Union[torch.Tensor, np.ndarray]] = None, + maxiter: int = 10, + stopping_tolerance: float = 10**-5, + spin_conserving: bool = False, + wavefuntion_real: bool = False, + outer_loop_callback: Optional[Callable] = None, + partial_unitary_random_perturbation: Optional[float] = None, + RDM_ops_batchsize: Optional[int] = None): + + """ + + Args: + problem: The ElectronicStructureProblem from which molecule information such as one and two body integrals is obtained. + integral_tensors: A tuple storing the one and two body integrals of the full orbital space. The first + entry stores the one body integrals and the second entry stores the two body integrals in the Physicist's notation in + the dense spin-orbital representation, represented as either :class:`torch.Tensor` or :class:`np.ndarray`. + num_spin_orbitals: The number of spin-orbitals to use for the active space. + ground_state_solver: An instance of VQE to use for the wavefunction optimization. + mapper: A QubitMapper to use for the RDM calculations. + partial_unitary_optimizer: An instance of PartialUnitaryProjectionOptimizer to use for the basis optimization. + initial_partial_unitary: The initial guess for the orbital rotation matrix. If ``None``, then a permutation matrix + selecting the spatial orbitals with the lowest energy will be generated. + maxiter: The maximum number of outerloop iterations. (The number of times the wavefunction optimization is run.) + stopping tolerance: The stopping tolerance used to determine if the algorithm should be stopped. + spin_conserving: A boolean flag that indicates whether or not we are assuming that spin is conserved + in the system being studied. Setting to True will skip the calculation of RDM element operators + which do not conserve spin, automatically setting these elements to zero. + wavefunction_real: A boolean flag that indicates whether or not we are assuming that the wavefunction is real. + outer_loop_callback: A callback function that tracks the outerloop progress. + It takes the outerloop iteration, latest VQE results, and latest outer loop results as arguments. + This can be used to save the most recent results to a file after each outer loop iteration. + partial_unitary_random_perturbation: A float representing the standard deviation of a normal distribution from which + the matrix elements of a perturbation matrix will be sampled. This perturbation matrix will be added to the initial + partial unitary (and re-orthonormalized) at the beginning of each orbital optimization run. + minimum_eigensolver_random_perturbation: A float representation the standard deviation of a normal distribution + from which the elements of a random set of parameters will be sampled. This perturbation will be added + to the initial parameters of the eigensolver at each outer loop iteration. + RDM_ops_batchsize: The number of Pauli ops to store in an :class:`Estimator` at any given time before deleting + the Estimator and replacing it with a blank copy. Increasing this number will lead to increased memory consumption. Setting + this number to be too low will hinder runtime performance. + + """ + + super().__init__(problem=problem, + integral_tensors=integral_tensors, + num_spin_orbitals=num_spin_orbitals, + ground_state_solver=ground_state_solver, + mapper=mapper, + estimator=estimator, + partial_unitary_optimizer=partial_unitary_optimizer, + initial_partial_unitary=initial_partial_unitary, + maxiter=maxiter, + stopping_tolerance=stopping_tolerance, + spin_conserving=spin_conserving, + wavefuntion_real=wavefuntion_real, + outer_loop_callback=outer_loop_callback, + partial_unitary_random_perturbation=partial_unitary_random_perturbation, + RDM_ops_batchsize=RDM_ops_batchsize) + + if ground_state_solver.__class__.__name__ != 'AdaptVQE': + + raise AlgorithmError(f"The ground state solver needs to be of type AdaptVQE, not {ground_state_solver.__class__.__name__}") + + def parameter_update_rule(self, result: MinimumEigensolverResult, + iteration: int): + + pass + +class OptOrbAdaptVQEResult(OptOrbMinimumEigensolverResult): + + def __init__(self) -> None: + super().__init__() + diff --git a/electronic_structure_algorithms/orbital_optimization/opt_orb_eigensolver.py b/electronic_structure_algorithms/orbital_optimization/opt_orb_eigensolver.py new file mode 100644 index 0000000..0432687 --- /dev/null +++ b/electronic_structure_algorithms/orbital_optimization/opt_orb_eigensolver.py @@ -0,0 +1,282 @@ +from typing import Callable, Optional, Union, Dict, List +from qiskit.primitives import BaseEstimator +from .partial_unitary_projection_optimizer import PartialUnitaryProjectionOptimizer +from qiskit_algorithms.eigensolvers import Eigensolver, EigensolverResult +from qiskit_nature.second_q.mappers import QubitMapper +from qiskit_nature.second_q.problems import ElectronicStructureProblem +import torch +import copy +import numpy as np +from electronic_structure_algorithms.excited_states_eigensolvers import MCVQE, SSVQE + +from .base_opt_orb_solver import BaseOptOrbSolver, BaseOptOrbResult + +class OptOrbEigensolver(BaseOptOrbSolver): + + def __init__(self, + problem: Optional[ElectronicStructureProblem], + integral_tensors: Optional[Union[tuple[torch.Tensor, torch.Tensor], tuple[np.ndarray, np.ndarray]]], + num_spin_orbitals: int, + excited_states_solver: Eigensolver, + mapper: QubitMapper, + estimator: BaseEstimator, + partial_unitary_optimizer: PartialUnitaryProjectionOptimizer, + initial_partial_unitary: Optional[Union[torch.Tensor, np.ndarray]] = None, + maxiter: int = 10, + stopping_tolerance: float = 10**-5, + spin_conserving: bool = False, + wavefuntion_real: bool = False, + outer_loop_callback: Optional[Callable] = None, + partial_unitary_random_perturbation: Optional[float] = None, + RDM_ops_batchsize: Optional[int] = None, + weight_vector: Optional[Union[list, np.ndarray]] = None): + + """ + + Args: + problem: The ElectronicStructureProblem from which molecule information such as one and two body integrals is obtained. + integral_tensors: A tuple storing the one and two body integrals of the full orbital space. The first + entry stores the one body integrals and the second entry stores the two body integrals in the Physicist's notation in + the dense spin-orbital representation, represented as either :class:`torch.Tensor` or :class:`np.ndarray`. + num_spin_orbitals: The number of spin-orbitals to use for the active space. + excited_states_solver: An instance of :class:`Eigensolver` for the excited states solver subroutine. + mapper: A QubitMapper to use for the RDM calculations. + partial_unitary_optimizer: An instance of PartialUnitaryProjectionOptimizer to use for the basis optimization. + initial_partial_unitary: The initial guess for the orbital rotation matrix. If ``None``, then a permutation matrix + selecting the spatial orbitals with the lowest energy will be generated. + maxiter: The maximum number of outerloop iterations. (The number of times the wavefunction optimization is run.) + stopping tolerance: The stopping tolerance used to determine if the algorithm should be stopped. + spin_conserving: A boolean flag that indicates whether or not we are assuming that spin is conserved + in the system being studied. Setting to True will skip the calculation of RDM element operators + which do not conserve spin, automatically setting these elements to zero. + wavefunction_real: A boolean flag that indicates whether or not we are assuming that the wavefunction is real. + outer_loop_callback: A callback function that tracks the outerloop progress. + It takes the outerloop iteration, latest eigensolver results, and latest outer loop results as arguments. + This can be used to save the most recent results to a file after each outer loop iteration. + partial_unitary_random_perturbation: A float representing the standard deviation of a normal distribution from which + the matrix elements of a perturbation matrix will be sampled. This perturbation matrix will be added to the initial + partial unitary (and re-orthonormalized) at the beginning of each orbital optimization run. + RDM_ops_batchsize: The number of Pauli ops to store in an :class:`Estimator` at any given time before deleting + the Estimator and replacing it with a blank copy. Increasing this number will lead to increased memory consumption. Setting + this number to be too low will hinder runtime performance. + weight_vector: A list or 1-D array of positive-real valued weights for the weighted sum orbital optimization + objective function whose length is equal to the number of excited states. + + """ + super().__init__(problem=problem, + integral_tensors=integral_tensors, + num_spin_orbitals=num_spin_orbitals, + mapper=mapper, + estimator=estimator, + partial_unitary_optimizer=partial_unitary_optimizer, + initial_partial_unitary=initial_partial_unitary, + maxiter=maxiter, + stopping_tolerance=stopping_tolerance, + spin_conserving=spin_conserving, + wavefuntion_real=wavefuntion_real, + outer_loop_callback=outer_loop_callback, + partial_unitary_random_perturbation=partial_unitary_random_perturbation, + RDM_ops_batchsize=RDM_ops_batchsize) + + # generate copies of the eigensolver instance to use for every outerloop iteration. + self._excited_states_solver_list = [copy.deepcopy(excited_states_solver) for n in range(int(maxiter+1))] + + self.num_states = excited_states_solver.k + if weight_vector is not None: + self.weight_vector = weight_vector + elif hasattr(excited_states_solver, 'weight_vector'): + self.weight_vector = excited_states_solver.weight_vector + else: + self.weight_vector = [self.num_states - n for n in range(self.num_states)] + self._energy_convergence_list = [] + self._pauli_ops_expectation_values_dict_list = None + + @property + def energy_convergence_list(self) -> List[float]: + """Returns the list of outerloop iteration energy values.""" + return self._energy_convergence_list + + @energy_convergence_list.setter + def energy_convergence_list(self, energy_list: List[float]) -> None: + """Sets the list of outerloop iteration energy values.""" + self._energy_convergence_list = energy_list + + @property + def pauli_ops_expectation_values_dict_list(self) -> list[Dict]: + """Returns the dictionary containing all of the expectation values + of all the Pauli string operators necessary for calculating + the RDMs with respect to a particular wavefunction.""" + return self._pauli_ops_expectation_values_dict + + @pauli_ops_expectation_values_dict_list.setter + def pauli_ops_expectation_values_dict_list(self, some_dict_list: list[Dict]) -> None: + """Sets the dictionary containing all of the expectation values + of all the Pauli string operators necessary for calculating + the RDMs with respect to a particular wavefunction.""" + self._pauli_ops_expectation_values_dict_list = some_dict_list + + def stopping_condition(self, iteration) -> bool: + + """Evaluates whether or not the stopping condition is True. + Returns True if the algorithm should be stopped, otherwise returns False. + """ + + if len(self._energy_convergence_list) >= 2: + if iteration == self.maxiter or np.abs(self._energy_convergence_list[-1] - self._energy_convergence_list[-2]) < self.stopping_tolerance: + return True + else: + return False + + else: + return False + + def compute_rotated_weighted_energy_sum(self, + partial_unitary: torch.Tensor, + oneRDM: list[torch.Tensor], + twoRDM: list[torch.Tensor], + one_body_integrals, + two_body_integrals): + + weighted_energy_sum = 0 + state_index = 0 + + for one_RDM, two_RDM in zip(oneRDM, twoRDM): + + weighted_energy_sum += torch.tensor(self.weight_vector[state_index], + dtype=torch.float64)*self.compute_rotated_energy(partial_unitary=partial_unitary, + oneRDM=one_RDM, + twoRDM=two_RDM, + one_body_integrals=one_body_integrals, + two_body_integrals=two_body_integrals) + state_index += 1 + + return weighted_energy_sum + + def compute_energies(self) -> EigensolverResult: + + self.outer_loop_iteration = 0 + + optorb_result = OptOrbEigensolverResult() + self._hamiltonian = self.get_rotated_hamiltonian(self._current_partial_unitary) + + self._pauli_op_dict = self.construct_pauli_op_dict(mapper=self.mapper) + + while self.stopping_condition(self.outer_loop_iteration) == False: + + result = self._excited_states_solver_list[self.outer_loop_iteration].compute_eigenvalues(operator=self._hamiltonian) + energies = np.real(result.eigenvalues) + opt_params = result.optimal_parameters + + # update the optorb result to hold the most recent VQE value. + optorb_result.eigenvalues = energies + + # update the optorb result to hold the most recent VQE parameters. + optorb_result.optimal_parameters = opt_params + + if isinstance(opt_params, dict): + opt_params = [copy.deepcopy(opt_params)] * self.num_states + + # update the optorb result to hold the most recent partial unitary basis transformation. + optorb_result.optimal_partial_unitary = self._current_partial_unitary + optorb_result.num_vqe_evaluations += 1 + + if self.outer_loop_callback is not None: + self.outer_loop_callback(self.outer_loop_iteration, result, optorb_result) + + self._energy_convergence_list.append(np.dot(self.weight_vector, energies)) + if isinstance(self._excited_states_solver_list[self.outer_loop_iteration], SSVQE) or isinstance(self._excited_states_solver_list[self.outer_loop_iteration], MCVQE): + states = [self._excited_states_solver_list[self.outer_loop_iteration].initial_states[n].compose(self._excited_states_solver_list[self.outer_loop_iteration].ansatz).bind_parameters(opt_params[n]) for n in range(self.num_states)] + else: + states = [self._excited_states_solver_list[self.outer_loop_iteration].ansatz[n].bind_parameters(opt_params[n]) for n in range(self.num_states)] + + if self.stopping_condition(self.outer_loop_iteration) == True: + break + + string_op_tuple_list = [(key, self._pauli_op_dict[key]) for key in self._pauli_op_dict] + + results = [[] for n in range(self.num_states)] + for n in range(self.num_states): + ops_counter = 1 + num_ops = len(string_op_tuple_list) + for op_tuple in string_op_tuple_list: + + results[n].append(self.estimator_list[self.outer_loop_iteration].run([states[n]], [op_tuple[1]]).result().values[0]) + + if self.RDM_ops_batchsize is not None: + if ops_counter % self.RDM_ops_batchsize == 0: + self.estimator_list[self.outer_loop_iteration] = copy.deepcopy(self.estimator) + + ops_counter += 1 + + self._pauli_ops_expectation_values_dict_list = [dict(zip([op_tuple[0] for op_tuple in string_op_tuple_list], results[n])) for n in range(self.num_states)] + + self.estimator_list[self.outer_loop_iteration] = None + + oneRDM_list = [self.get_one_RDM_tensor(mapper=self.mapper, expectval_dict=self._pauli_ops_expectation_values_dict_list[n]) for n in range(self.num_states)] + twoRDM_list = [self.get_two_RDM_tensor(mapper=self.mapper, expectval_dict=self._pauli_ops_expectation_values_dict_list[n]) for n in range(self.num_states)] + + if self.partial_unitary_random_perturbation is not None: + + initial_partial_unitary = self.orth(self._current_partial_unitary + torch.Tensor(np.random.normal(loc=0.0, + scale=self.partial_unitary_random_perturbation, size=(self._current_partial_unitary.size()[0], + self._current_partial_unitary.size()[1])))) + else: + + initial_partial_unitary = self._current_partial_unitary + + oneRDM_list = [oneRDM.to(self._partial_unitary_optimizer_list[self.outer_loop_iteration].device) for oneRDM in oneRDM_list] + twoRDM_list = [twoRDM.to(self._partial_unitary_optimizer_list[self.outer_loop_iteration].device) for twoRDM in twoRDM_list] + self.one_body_integrals = self.one_body_integrals.to(self._partial_unitary_optimizer_list[self.outer_loop_iteration].device) + self.two_body_integrals = self.two_body_integrals.to(self._partial_unitary_optimizer_list[self.outer_loop_iteration].device) + self._current_partial_unitary = self._partial_unitary_optimizer_list[self.outer_loop_iteration].compute_optimal_rotation(fun=self.compute_rotated_weighted_energy_sum, + oneRDM=oneRDM_list, + twoRDM=twoRDM_list, + one_body_integrals=self.one_body_integrals, + two_body_integrals=self.two_body_integrals, + initial_partial_unitary=initial_partial_unitary)[0] + oneRDM_list = [oneRDM.to('cpu') for oneRDM in oneRDM_list] + twoRDM_list = [twoRDM.to('cpu') for twoRDM in twoRDM_list] + del oneRDM_list + del twoRDM_list + del string_op_tuple_list + self.one_body_integrals = self.one_body_integrals.to('cpu') + self.two_body_integrals = self.two_body_integrals.to('cpu') + + self._hamiltonian = self.get_rotated_hamiltonian(self._current_partial_unitary) + self.outer_loop_iteration += 1 + + self.parameter_update_rule(result, self.outer_loop_iteration) + + self._partial_unitary_optimizer_list[self.outer_loop_iteration - 1] = None + self._excited_states_solver_list[self.outer_loop_iteration - 1] = None + + return optorb_result + +class OptOrbEigensolverResult(BaseOptOrbResult, EigensolverResult): + + def __init__(self) -> None: + super().__init__() + self._num_vqe_evaluations = 0 + self._optimal_partial_unitary = None + + @property + def num_vqe_evaluations(self) -> int: + """Returns the number of times the eigensolver was run.""" + return self._num_vqe_evaluations + + @num_vqe_evaluations.setter + def num_vqe_evaluations(self, some_int: int) -> None: + """Sets the number of times the eigensolver was run.""" + self._num_vqe_evaluations = some_int + + @property + def optimal_partial_unitary(self) -> torch.Tensor: + """Returns the optimal partial unitary basis transformation.""" + return self._optimal_partial_unitary + + @optimal_partial_unitary.setter + def optimal_partial_unitary(self, some_tensor: torch.Tensor) -> None: + """Sets the optimal partial unitary basis transformation.""" + self._optimal_partial_unitary = some_tensor + + diff --git a/electronic_structure_algorithms/orbital_optimization/opt_orb_mcvqe.py b/electronic_structure_algorithms/orbital_optimization/opt_orb_mcvqe.py new file mode 100644 index 0000000..6fe2df8 --- /dev/null +++ b/electronic_structure_algorithms/orbital_optimization/opt_orb_mcvqe.py @@ -0,0 +1,123 @@ +from typing import Callable, Optional, Union +from qiskit.primitives import BaseEstimator +from .partial_unitary_projection_optimizer import PartialUnitaryProjectionOptimizer +from qiskit.algorithms.eigensolvers import Eigensolver +from qiskit.algorithms.exceptions import AlgorithmError +from qiskit_nature.second_q.mappers import QubitMapper +from qiskit_nature.second_q.problems import ElectronicStructureProblem +import torch +import numpy as np + +from .opt_orb_eigensolver import OptOrbEigensolver, OptOrbEigensolverResult + +class OptOrbMCVQE(OptOrbEigensolver): + + def __init__(self, + problem: ElectronicStructureProblem, + integral_tensors: Optional[Union[tuple[torch.Tensor, torch.Tensor], tuple[np.ndarray, np.ndarray]]], + num_spin_orbitals: int, + excited_states_solver: Eigensolver, + mapper: QubitMapper, + estimator: BaseEstimator, + partial_unitary_optimizer: PartialUnitaryProjectionOptimizer, + initial_partial_unitary: Optional[Union[torch.Tensor, np.ndarray]] = None, + maxiter: int = 10, + stopping_tolerance: float = 10**-5, + spin_conserving: bool = False, + wavefuntion_real: bool = False, + outer_loop_callback: Optional[Callable] = None, + partial_unitary_random_perturbation: Optional[float] = None, + eigensolver_random_perturbation: Optional[float] = None, + RDM_ops_batchsize: Optional[int] = None): + + """ + + Args: + problem: The ElectronicStructureProblem from which molecule information such as one and two body integrals is obtained. + integral_tensors: A tuple storing the one and two body integrals of the full orbital space. The first + entry stores the one body integrals and the second entry stores the two body integrals in the Physicist's notation in + the dense spin-orbital representation, represented as either :class:`torch.Tensor` or :class:`np.ndarray`. + num_spin_orbitals: The number of spin-orbitals to use for the active space. + excited_states_solver: An instance of VQE to use for the wavefunction optimization. + mapper: A QubitMapper to use for the RDM calculations. + partial_unitary_optimizer: An instance of PartialUnitaryProjectionOptimizer to use for the basis optimization. + initial_partial_unitary: The initial guess for the orbital rotation matrix. If ``None``, then a permutation matrix + selecting the spatial orbitals with the lowest energy will be generated. + maxiter: The maximum number of outerloop iterations. (The number of times the wavefunction optimization is run.) + stopping tolerance: The stopping tolerance used to determine if the algorithm should be stopped. + spin_conserving: A boolean flag that indicates whether or not we are assuming that spin is conserved + in the system being studied. Setting to True will skip the calculation of RDM element operators + which do not conserve spin, automatically setting these elements to zero. + wavefunction_real: A boolean flag that indicates whether or not we are assuming that the wavefunction is real. + outer_loop_callback: A callback function that tracks the outerloop progress. + It takes the outerloop iteration, latest VQE results, and latest outer loop results as arguments. + This can be used to save the most recent results to a file after each outer loop iteration. + partial_unitary_random_perturbation: A float representing the standard deviation of a normal distribution from which + the matrix elements of a perturbation matrix will be sampled. This perturbation matrix will be added to the initial + partial unitary (and re-orthonormalized) at the beginning of each orbital optimization run. + eigensolver_random_perturbation: A float representation the standard deviation of a normal distribution from which + the elements of a set of perturbation parameters will be sampled. These perturbation elements will be added + to the initial parameters of the eigensolver at each outer loop iteration. + RDM_ops_batchsize: The number of Pauli ops to store in an :class:`Estimator` at any given time before deleting + the Estimator and replacing it with a blank copy. Increasing this number will lead to increased memory consumption. Setting + this number to be too low will hinder runtime performance. + + """ + + super().__init__(problem=problem, + integral_tensors=integral_tensors, + num_spin_orbitals=num_spin_orbitals, + excited_states_solver=excited_states_solver, + mapper=mapper, + estimator=estimator, + partial_unitary_optimizer=partial_unitary_optimizer, + initial_partial_unitary=initial_partial_unitary, + maxiter=maxiter, + stopping_tolerance=stopping_tolerance, + spin_conserving=spin_conserving, + wavefuntion_real=wavefuntion_real, + outer_loop_callback=outer_loop_callback, + partial_unitary_random_perturbation=partial_unitary_random_perturbation, + RDM_ops_batchsize=RDM_ops_batchsize) + + if excited_states_solver.__class__.__name__ != 'MCVQE': + + raise AlgorithmError(f"The ground state solver needs to be of type MCVQE, not {excited_states_solver.__class__.__name__}") + + + self.eigensolver_random_perturbation = eigensolver_random_perturbation + print(self.initial_partial_unitary) + print() + print(self.one_body_integrals) + print() + init_partial_unitary = torch.block_diag(self.initial_partial_unitary, self.initial_partial_unitary) + + rotated_one_body_integrals = torch.einsum('pq,pi,qj->ij', + self.one_body_integrals, + init_partial_unitary, init_partial_unitary).detach().numpy() + rotated_two_body_integrals = torch.einsum('pqrs,pi,qj,rk,sl->ijkl', + self.two_body_integrals, + init_partial_unitary, init_partial_unitary, init_partial_unitary, + init_partial_unitary).detach().numpy() + + for solver in self._excited_states_solver_list: + + solver.one_body_integrals = rotated_one_body_integrals + solver.two_body_integrals = rotated_two_body_integrals + print(f'set rotated integrals') + + def parameter_update_rule(self, result: OptOrbEigensolverResult, + iteration: int): + + if self.eigensolver_random_perturbation == None or self.eigensolver_random_perturbation == 0.0: + self._excited_states_solver_list[iteration].initial_point = result.optimal_point + else: + + self._excited_states_solver_list[iteration].initial_point = result.optimal_point + np.random.normal(loc=0.0, + scale=self.eigensolver_random_perturbation, size=result.optimal_point.size) + +class OptOrbMCVQEResult(OptOrbEigensolverResult): + + def __init__(self) -> None: + super().__init__() + diff --git a/electronic_structure_algorithms/orbital_optimization/opt_orb_minimum_eigensolver.py b/electronic_structure_algorithms/orbital_optimization/opt_orb_minimum_eigensolver.py new file mode 100644 index 0000000..bb324fe --- /dev/null +++ b/electronic_structure_algorithms/orbital_optimization/opt_orb_minimum_eigensolver.py @@ -0,0 +1,253 @@ +from typing import Callable, Optional, Union, Dict, List +from abc import abstractmethod +from qiskit.primitives import BaseEstimator +from .partial_unitary_projection_optimizer import PartialUnitaryProjectionOptimizer +from qiskit.algorithms.minimum_eigensolvers import MinimumEigensolverResult, MinimumEigensolver +from qiskit.algorithms.minimum_eigensolvers import MinimumEigensolverResult +from qiskit_nature.second_q.mappers import QubitMapper +from qiskit_nature.second_q.problems import ElectronicStructureProblem +import torch +import copy +import numpy as np + +from .base_opt_orb_solver import BaseOptOrbSolver, BaseOptOrbResult + +class OptOrbMinimumEigensolver(BaseOptOrbSolver): + + def __init__(self, + problem: ElectronicStructureProblem, + integral_tensors: Optional[Union[tuple[torch.Tensor, torch.Tensor], tuple[np.ndarray, np.ndarray]]], + num_spin_orbitals: int, + ground_state_solver: MinimumEigensolver, + mapper: QubitMapper, + estimator: BaseEstimator, + partial_unitary_optimizer: PartialUnitaryProjectionOptimizer, + initial_partial_unitary: Optional[Union[torch.Tensor, np.ndarray]] = None, + maxiter: int = 10, + stopping_tolerance: float = 10**-5, + spin_conserving: bool = False, + wavefuntion_real: bool = False, + outer_loop_callback: Optional[Callable] = None, + partial_unitary_random_perturbation: Optional[float] = None, + RDM_ops_batchsize: Optional[int] = None): + + """ + + Args: + problem: An ElectronicStructureProblem from which molecule information such as one and two body integrals is obtained. + integral_tensors: A tuple storing the one and two body integrals of the full orbital space. The first + entry stores the one body integrals and the second entry stores the two body integrals in the Physicist's notation in + the dense spin-orbital representation, represented as either :class:`torch.Tensor` or :class:`np.ndarray`. + num_spin_orbitals: The number of spin-orbitals to use for the active space. + ground_state_solver: A :class:`MinimumEigensolver` used for the ground state solver subroutine. + mapper: A QubitMapper to use for the RDM calculations. + partial_unitary_optimizer: An instance of :class:`PartialUnitaryProjectionOptimizer` to use for + the basis optimization subroutine. + initial_partial_unitary: The initial guess for the orbital rotation matrix. If ``None``, then a permutation matrix + selecting the spatial orbitals with the lowest energy will be generated. + maxiter: The maximum number of outerloop iterations. (The number of times the wavefunction optimization is run.) + stopping tolerance: The stopping tolerance used to determine if the algorithm should be stopped. + spin_conserving: A boolean flag that indicates whether or not we are assuming that spin is conserved + in the system being studied. Setting to ``True`` will skip the calculation of RDM element operators + which do not conserve spin, automatically setting these elements to zero. This should only be + set to ``True`` if it is known that the circuit used to construct the wavefunction in the eigensolver + subroutine conserves the spin-magnetization quantum number (e.g. UCCSD). + wavefunction_real: A boolean flag that indicates whether or not we are assuming that the wavefunction is real. + This allows for the exploitation of symmetries + in the RDM calculation to reduce runtime, but should only be set to ``True`` + if it is known that the circuit used in the eigensolver subroutine produces wavefunction coefficients + which are real-valued. + outer_loop_callback: A callback function that tracks the outerloop progress. + It takes the outerloop iteration, latest eigensolver results, and latest outer loop results as arguments. + This can be used to save the most recent results to a file after each outer loop iteration. + partial_unitary_random_perturbation: A float representing the standard deviation of a normal distribution from which + the matrix elements of a perturbation matrix will be sampled. This perturbation matrix will be added to the initial + partial unitary (and re-orthonormalized) at the beginning of each orbital optimization run. + RDM_ops_batchsize: The number of Pauli ops to store in an :class:`Estimator` at any given time before deleting + the Estimator and replacing it with a blank copy. Increasing this number will lead to increased memory consumption. Setting + this number to be too low will hinder runtime performance. + + """ + + super().__init__(problem=problem, + integral_tensors=integral_tensors, + num_spin_orbitals=num_spin_orbitals, + mapper=mapper, + estimator=estimator, + partial_unitary_optimizer=partial_unitary_optimizer, + initial_partial_unitary=initial_partial_unitary, + maxiter=maxiter, + stopping_tolerance=stopping_tolerance, + spin_conserving=spin_conserving, + wavefuntion_real=wavefuntion_real, + outer_loop_callback=outer_loop_callback, + partial_unitary_random_perturbation=partial_unitary_random_perturbation, + RDM_ops_batchsize=RDM_ops_batchsize) + + self._ground_state_solver_list = [copy.deepcopy(ground_state_solver) for n in range(int(maxiter+1))] + self._energy_convergence_list = [] + self._pauli_ops_expectation_values_dict = None + + @property + def ground_state_solver_list(self) -> List[MinimumEigensolver]: + """Returns the list of ground state solver instances.""" + return self._ground_state_solver_list + + @ground_state_solver_list.setter + def ground_state_solver_list(self, instance_list: List[MinimumEigensolver]) -> None: + """Sets the list of ground state solver instances.""" + self._ground_state_solver_list = instance_list + + @property + def energy_convergence_list(self) -> List[float]: + """Returns the list of outerloop iteration energy values.""" + return self._energy_convergence_list + + @energy_convergence_list.setter + def energy_convergence_list(self, energy_list: List[float]) -> None: + """Sets the list of outerloop iteration energy values.""" + self._energy_convergence_list = energy_list + + @property + def pauli_ops_expectation_values_dict(self) -> Dict: + """Returns the dictionary containing all of the expectation values + of all the Pauli string operators necessary for calculating + the RDMs with respect to a particular wavefunction.""" + return self._pauli_ops_expectation_values_dict + + @pauli_ops_expectation_values_dict.setter + def pauli_ops_expectation_values_dict(self, some_dict: Dict) -> None: + """Sets the dictionary containing all of the expectation values + of all the Pauli string operators necessary for calculating + the RDMs with respect to a particular wavefunction.""" + self._pauli_ops_expectation_values_dict = some_dict + + def stopping_condition(self, iteration) -> bool: + + """Evaluates whether or not the stopping condition is True. + Returns True if the algorithm should be stopped, otherwise returns False. + """ + + if len(self._energy_convergence_list) >= 2: + if iteration == self.maxiter or np.abs(self._energy_convergence_list[-1] - self._energy_convergence_list[-2]) < self.stopping_tolerance: + return True + else: + return False + + else: + return False + + @abstractmethod + def parameter_update_rule(self, result: MinimumEigensolverResult, iteration: int): + + raise NotImplementedError("Minimum eigensolver needs to implement a way to update parameters after each orbital optimization.") + + @abstractmethod + def return_RDM_circuit(self, result: MinimumEigensolverResult, iteration: int): + + NotImplementedError("Minimum eigensolver needs to implement a way to return the circuit used to calculate the one and two RDM.") + + def compute_minimum_energy(self) -> MinimumEigensolverResult: + + outer_loop_iteration = 0 + + optorb_result = OptOrbMinimumEigensolverResult() + self._hamiltonian = self.get_rotated_hamiltonian(self._current_partial_unitary) + + self._pauli_op_dict = self.construct_pauli_op_dict(mapper=self.mapper) + + while self.stopping_condition(outer_loop_iteration) == False: + + result = self._ground_state_solver_list[outer_loop_iteration].compute_minimum_eigenvalue(operator=self._hamiltonian) + energy = np.real(result.eigenvalue) + opt_params = result.optimal_parameters + + # update the optorb result to hold the most recent ground state energy. + optorb_result.eigenvalue = energy + + # update the optorb result to hold the most recent ground state parameters. + optorb_result.optimal_parameters = opt_params + + # update the optorb result to hold the most recent partial unitary basis transformation. + optorb_result.optimal_partial_unitary = self._current_partial_unitary + optorb_result.num_vqe_evaluations += 1 + optorb_result.optimal_circuit = result.optimal_circuit + optorb_result.optimal_point = result.optimal_point + optorb_result.optimal_value = result.optimal_value + + if self.outer_loop_callback is not None: + self.outer_loop_callback(outer_loop_iteration, result, optorb_result) + + self._energy_convergence_list.append(energy) + state = copy.deepcopy(result.optimal_circuit).bind_parameters(opt_params) + + if self.stopping_condition(outer_loop_iteration) == True: + break + + string_op_tuple_list = [(key, self._pauli_op_dict[key]) for key in self._pauli_op_dict] + + results = [] + ops_counter = 1 + num_ops = len(string_op_tuple_list) + for op_tuple in string_op_tuple_list: + + results.append(self.estimator_list[outer_loop_iteration].run([state], [op_tuple[1]]).result().values[0]) + + if self.RDM_ops_batchsize is not None: + if ops_counter % self.RDM_ops_batchsize == 0: + self.estimator_list[outer_loop_iteration] = copy.deepcopy(self.estimator) + + ops_counter += 1 + + self._pauli_ops_expectation_values_dict = dict(zip([op_tuple[0] for op_tuple in string_op_tuple_list], results)) + + + self.estimator_list[outer_loop_iteration] = None + + oneRDM = self.get_one_RDM_tensor(mapper=self.mapper, expectval_dict=self._pauli_ops_expectation_values_dict) + twoRDM = self.get_two_RDM_tensor(mapper=self.mapper, expectval_dict=self._pauli_ops_expectation_values_dict) + + if self.partial_unitary_random_perturbation is not None: + + initial_partial_unitary = self.orth(self._current_partial_unitary + torch.Tensor(np.random.normal(loc=0.0, + scale=self.partial_unitary_random_perturbation, size=(self._current_partial_unitary.size()[0], + self._current_partial_unitary.size()[1])))) + else: + + initial_partial_unitary = self._current_partial_unitary + + oneRDM = oneRDM.to(self._partial_unitary_optimizer_list[outer_loop_iteration].device) + twoRDM = twoRDM.to(self._partial_unitary_optimizer_list[outer_loop_iteration].device) + self.one_body_integrals = self.one_body_integrals.to(self._partial_unitary_optimizer_list[outer_loop_iteration].device) + self.two_body_integrals = self.two_body_integrals.to(self._partial_unitary_optimizer_list[outer_loop_iteration].device) + self._current_partial_unitary = self._partial_unitary_optimizer_list[outer_loop_iteration].compute_optimal_rotation(fun=self.compute_rotated_energy, + oneRDM=oneRDM, + twoRDM=twoRDM, + one_body_integrals=self.one_body_integrals, + two_body_integrals=self.two_body_integrals, + initial_partial_unitary=initial_partial_unitary)[0] + oneRDM = oneRDM.to('cpu') + twoRDM = twoRDM.to('cpu') + del oneRDM + del twoRDM + del string_op_tuple_list + self.one_body_integrals = self.one_body_integrals.to('cpu') + self.two_body_integrals = self.two_body_integrals.to('cpu') + + self._hamiltonian = self.get_rotated_hamiltonian(self._current_partial_unitary) + outer_loop_iteration += 1 + + self.parameter_update_rule(result, outer_loop_iteration) + self._ground_state_solver_list[outer_loop_iteration].initial_point = result.optimal_point + + self._partial_unitary_optimizer_list[outer_loop_iteration - 1] = None + self._ground_state_solver_list[outer_loop_iteration - 1] = None + + return optorb_result + +class OptOrbMinimumEigensolverResult(BaseOptOrbResult, MinimumEigensolverResult): + + def __init__(self) -> None: + super().__init__() + + diff --git a/electronic_structure_algorithms/orbital_optimization/opt_orb_ssvqe.py b/electronic_structure_algorithms/orbital_optimization/opt_orb_ssvqe.py new file mode 100644 index 0000000..d152550 --- /dev/null +++ b/electronic_structure_algorithms/orbital_optimization/opt_orb_ssvqe.py @@ -0,0 +1,106 @@ +from typing import Callable, Optional, Union +from qiskit.primitives import BaseEstimator +from .partial_unitary_projection_optimizer import PartialUnitaryProjectionOptimizer +from qiskit.algorithms.eigensolvers import Eigensolver +from qiskit.algorithms.exceptions import AlgorithmError +from qiskit_nature.second_q.mappers import QubitMapper +from qiskit_nature.second_q.problems import ElectronicStructureProblem +import torch +import numpy as np + +from .opt_orb_eigensolver import OptOrbEigensolver, OptOrbEigensolverResult + +class OptOrbSSVQE(OptOrbEigensolver): + + def __init__(self, + problem: ElectronicStructureProblem, + integral_tensors: Optional[Union[tuple[torch.Tensor, torch.Tensor], tuple[np.ndarray, np.ndarray]]], + num_spin_orbitals: int, + excited_states_solver: Eigensolver, + mapper: QubitMapper, + estimator: BaseEstimator, + partial_unitary_optimizer: PartialUnitaryProjectionOptimizer, + initial_partial_unitary: Optional[Union[torch.Tensor, np.ndarray]] = None, + maxiter: int = 10, + stopping_tolerance: float = 10**-5, + spin_conserving: bool = False, + wavefuntion_real: bool = False, + outer_loop_callback: Optional[Callable] = None, + partial_unitary_random_perturbation: Optional[float] = None, + eigensolver_random_perturbation: Optional[float] = None, + RDM_ops_batchsize: Optional[int] = None, + weight_vector: Optional[Union[list, np.ndarray]] = None): + + """ + + Args: + problem: The ElectronicStructureProblem from which molecule information such as one and two body integrals is obtained. + integral_tensors: A tuple storing the one and two body integrals of the full orbital space. The first + entry stores the one body integrals and the second entry stores the two body integrals in the Physicist's notation in + the dense spin-orbital representation, represented as either :class:`torch.Tensor` or :class:`np.ndarray`. + num_spin_orbitals: The number of spin-orbitals to use for the active space. + excited_states_solver: An instance of VQE to use for the wavefunction optimization. + mapper: A QubitMapper to use for the RDM calculations. + partial_unitary_optimizer: An instance of PartialUnitaryProjectionOptimizer to use for the basis optimization. + initial_partial_unitary: The initial guess for the orbital rotation matrix. If ``None``, then a permutation matrix + selecting the spatial orbitals with the lowest energy will be generated. + maxiter: The maximum number of outerloop iterations. (The number of times the wavefunction optimization is run.) + stopping tolerance: The stopping tolerance used to determine if the algorithm should be stopped. + spin_conserving: A boolean flag that indicates whether or not we are assuming that spin is conserved + in the system being studied. Setting to True will skip the calculation of RDM element operators + which do not conserve spin, automatically setting these elements to zero. + wavefunction_real: A boolean flag that indicates whether or not we are assuming that the wavefunction is real. + outer_loop_callback: A callback function that tracks the outerloop progress. + It takes the outerloop iteration, latest VQE results, and latest outer loop results as arguments. + This can be used to save the most recent results to a file after each outer loop iteration. + partial_unitary_random_perturbation: A float representing the standard deviation of a normal distribution from which + the matrix elements of a perturbation matrix will be sampled. This perturbation matrix will be added to the initial + partial unitary (and re-orthonormalized) at the beginning of each orbital optimization run. + eigensolver_random_perturbation: A float representation the standard deviation of a normal distribution from which + the elements of a set of perturbation parameters will be sampled. These perturbation elements will be added + to the initial parameters of the eigensolver at each outer loop iteration. + RDM_ops_batchsize: The number of Pauli ops to store in an :class:`Estimator` at any given time before deleting + the Estimator and replacing it with a blank copy. Increasing this number will lead to increased memory consumption. Setting + this number to be too low will hinder runtime performance. + weight_vector: A list or 1-D array of positive-real valued weights for the weighted sum orbital optimization + objective function whose length is equal to the number of excited states. + """ + + super().__init__(problem=problem, + integral_tensors=integral_tensors, + num_spin_orbitals=num_spin_orbitals, + excited_states_solver=excited_states_solver, + mapper=mapper, + estimator=estimator, + partial_unitary_optimizer=partial_unitary_optimizer, + initial_partial_unitary=initial_partial_unitary, + maxiter=maxiter, + stopping_tolerance=stopping_tolerance, + spin_conserving=spin_conserving, + wavefuntion_real=wavefuntion_real, + outer_loop_callback=outer_loop_callback, + partial_unitary_random_perturbation=partial_unitary_random_perturbation, + RDM_ops_batchsize=RDM_ops_batchsize, + weight_vector=weight_vector) + + if excited_states_solver.__class__.__name__ != 'SSVQE': + + raise AlgorithmError(f"The ground state solver needs to be of type SSVQE, not {excited_states_solver.__class__.__name__}") + + self.eigensolver_random_perturbation = eigensolver_random_perturbation + + def parameter_update_rule(self, result: OptOrbEigensolverResult, + iteration: int): + + if self.eigensolver_random_perturbation == None or self.eigensolver_random_perturbation == 0.0: + self._excited_states_solver_list[iteration].initial_point = result.optimal_point + else: + print('Perturb SSVQE parameters!') + self._excited_states_solver_list[iteration].initial_point = result.optimal_point + np.random.normal(loc=0.0, + scale=self.eigensolver_random_perturbation, size=result.optimal_point.size) + +class OptOrbSSVQEResult(OptOrbEigensolverResult): + + def __init__(self) -> None: + super().__init__() + diff --git a/electronic_structure_algorithms/orbital_optimization/opt_orb_vqd.py b/electronic_structure_algorithms/orbital_optimization/opt_orb_vqd.py new file mode 100644 index 0000000..3285909 --- /dev/null +++ b/electronic_structure_algorithms/orbital_optimization/opt_orb_vqd.py @@ -0,0 +1,110 @@ +from typing import Callable, Optional, Union +from qiskit.primitives import BaseEstimator +from .partial_unitary_projection_optimizer import PartialUnitaryProjectionOptimizer +from qiskit.algorithms.eigensolvers import Eigensolver +from qiskit.algorithms.exceptions import AlgorithmError +from qiskit_nature.second_q.mappers import QubitMapper +from qiskit_nature.second_q.problems import ElectronicStructureProblem +import torch +import numpy as np + +from .opt_orb_eigensolver import OptOrbEigensolver, OptOrbEigensolverResult + +class OptOrbVQD(OptOrbEigensolver): + + def __init__(self, + problem: ElectronicStructureProblem, + integral_tensors: Optional[Union[tuple[torch.Tensor, torch.Tensor], tuple[np.ndarray, np.ndarray]]], + num_spin_orbitals: int, + excited_states_solver: Eigensolver, + mapper: QubitMapper, + estimator: BaseEstimator, + partial_unitary_optimizer: PartialUnitaryProjectionOptimizer, + initial_partial_unitary: Optional[Union[torch.Tensor, np.ndarray]] = None, + maxiter: int = 10, + stopping_tolerance: float = 10**-5, + spin_conserving: bool = False, + wavefuntion_real: bool = False, + outer_loop_callback: Optional[Callable] = None, + partial_unitary_random_perturbation: Optional[float] = None, + eigensolver_random_perturbation: Optional[float] = None, + RDM_ops_batchsize: Optional[int] = None, + weight_vector: Optional[Union[list, np.ndarray]] = None): + + """ + + Args: + problem: The ElectronicStructureProblem from which molecule information such as one and two body integrals is obtained. + integral_tensors: A tuple storing the one and two body integrals of the full orbital space. The first + entry stores the one body integrals and the second entry stores the two body integrals in the Physicist's notation in + the dense spin-orbital representation, represented as either :class:`torch.Tensor` or :class:`np.ndarray`. + num_spin_orbitals: The number of spin-orbitals to use for the active space. + excited_states_solver: An instance of VQE to use for the wavefunction optimization. + mapper: A QubitMapper to use for the RDM calculations. + partial_unitary_optimizer: An instance of PartialUnitaryProjectionOptimizer to use for the basis optimization. + initial_partial_unitary: The initial guess for the orbital rotation matrix. If ``None``, then a permutation matrix + selecting the spatial orbitals with the lowest energy will be generated. + maxiter: The maximum number of outerloop iterations. (The number of times the wavefunction optimization is run.) + stopping tolerance: The stopping tolerance used to determine if the algorithm should be stopped. + spin_conserving: A boolean flag that indicates whether or not we are assuming that spin is conserved + in the system being studied. Setting to True will skip the calculation of RDM element operators + which do not conserve spin, automatically setting these elements to zero. + wavefunction_real: A boolean flag that indicates whether or not we are assuming that the wavefunction is real. + outer_loop_callback: A callback function that tracks the outerloop progress. + It takes the outerloop iteration, latest VQE results, and latest outer loop results as arguments. + This can be used to save the most recent results to a file after each outer loop iteration. + partial_unitary_random_perturbation: A float representing the standard deviation of a normal distribution from which + the matrix elements of a perturbation matrix will be sampled. This perturbation matrix will be added to the initial + partial unitary (and re-orthonormalized) at the beginning of each orbital optimization run. + eigensolver_random_perturbation: A float representation the standard deviation of a normal distribution from which + the elements of a set of perturbation parameters will be sampled. These perturbation elements will be added + to the initial parameters of the eigensolver at each outer loop iteration. + RDM_ops_batchsize: The number of Pauli ops to store in an :class:`Estimator` at any given time before deleting + the Estimator and replacing it with a blank copy. Increasing this number will lead to increased memory consumption. Setting + this number to be too low will hinder runtime performance. + weight_vector: A list or 1-D array of positive-real valued weights for the weighted sum orbital optimization + objective function whose length is equal to the number of excited states. + + """ + + super().__init__(problem=problem, + integral_tensors=integral_tensors, + num_spin_orbitals=num_spin_orbitals, + excited_states_solver=excited_states_solver, + mapper=mapper, + estimator=estimator, + partial_unitary_optimizer=partial_unitary_optimizer, + initial_partial_unitary=initial_partial_unitary, + maxiter=maxiter, + stopping_tolerance=stopping_tolerance, + spin_conserving=spin_conserving, + wavefuntion_real=wavefuntion_real, + outer_loop_callback=outer_loop_callback, + partial_unitary_random_perturbation=partial_unitary_random_perturbation, + RDM_ops_batchsize=RDM_ops_batchsize, + weight_vector=weight_vector) + + if excited_states_solver.__class__.__name__ != 'VQD': + + raise AlgorithmError(f"The ground state solver needs to be of type VQD, not {excited_states_solver.__class__.__name__}") + + self.eigensolver_random_perturbation = eigensolver_random_perturbation + + def parameter_update_rule(self, result: OptOrbEigensolverResult, + iteration: int): + + if self.eigensolver_random_perturbation == None or self.eigensolver_random_perturbation == 0.0: + + for n in range(self._excited_states_solver_list[iteration].k): + self._excited_states_solver_list[iteration].initial_point[n] = result.optimal_points[n] + else: + + for n in range(self._excited_states_solver_list[iteration].k): + self._excited_states_solver_list[iteration].initial_point[n] = result.optimal_points[n] + np.random.normal(loc=0.0, + scale=self.eigensolver_random_perturbation, size=result.optimal_points[n].size) + +class OptOrbVQDResult(OptOrbEigensolverResult): + + def __init__(self) -> None: + super().__init__() + diff --git a/electronic_structure_algorithms/orbital_optimization/opt_orb_vqe.py b/electronic_structure_algorithms/orbital_optimization/opt_orb_vqe.py new file mode 100644 index 0000000..8064069 --- /dev/null +++ b/electronic_structure_algorithms/orbital_optimization/opt_orb_vqe.py @@ -0,0 +1,103 @@ +from typing import Callable, Optional, Union +from qiskit.primitives import BaseEstimator +from .partial_unitary_projection_optimizer import PartialUnitaryProjectionOptimizer +from qiskit.algorithms.minimum_eigensolvers import MinimumEigensolver +from qiskit.algorithms.exceptions import AlgorithmError +from qiskit_nature.second_q.mappers import QubitMapper +from qiskit_nature.second_q.problems import ElectronicStructureProblem +import torch +import numpy as np + +from .opt_orb_minimum_eigensolver import OptOrbMinimumEigensolver, OptOrbMinimumEigensolverResult + +class OptOrbVQE(OptOrbMinimumEigensolver): + + def __init__(self, + problem: ElectronicStructureProblem, + integral_tensors: Optional[Union[tuple[torch.Tensor, torch.Tensor], tuple[np.ndarray, np.ndarray]]], + num_spin_orbitals: int, + ground_state_solver: MinimumEigensolver, + mapper: QubitMapper, + estimator: BaseEstimator, + partial_unitary_optimizer: PartialUnitaryProjectionOptimizer, + initial_partial_unitary: Optional[Union[torch.Tensor, np.ndarray]] = None, + maxiter: int = 10, + stopping_tolerance: float = 10**-5, + spin_conserving: bool = False, + wavefuntion_real: bool = False, + outer_loop_callback: Optional[Callable] = None, + partial_unitary_random_perturbation: Optional[float] = None, + minimum_eigensolver_random_perturbation: Optional[float] = None, + RDM_ops_batchsize: Optional[int] = None): + + """ + + Args: + problem: The ElectronicStructureProblem from which molecule information such as one and two body integrals is obtained. + integral_tensors: A tuple storing the one and two body integrals of the full orbital space. The first + entry stores the one body integrals and the second entry stores the two body integrals in the Physicist's notation in + the dense spin-orbital representation, represented as either :class:`torch.Tensor` or :class:`np.ndarray`. + num_spin_orbitals: The number of spin-orbitals to use for the active space. + ground_state_solver: An instance of VQE to use for the wavefunction optimization. + mapper: A QubitMapper to use for the RDM calculations. + partial_unitary_optimizer: An instance of PartialUnitaryProjectionOptimizer to use for the basis optimization. + initial_partial_unitary: The initial guess for the orbital rotation matrix. If ``None``, then a permutation matrix + selecting the spatial orbitals with the lowest energy will be generated. + maxiter: The maximum number of outerloop iterations. (The number of times the wavefunction optimization is run.) + stopping tolerance: The stopping tolerance used to determine if the algorithm should be stopped. + spin_conserving: A boolean flag that indicates whether or not we are assuming that spin is conserved + in the system being studied. Setting to True will skip the calculation of RDM element operators + which do not conserve spin, automatically setting these elements to zero. + wavefunction_real: A boolean flag that indicates whether or not we are assuming that the wavefunction is real. + outer_loop_callback: A callback function that tracks the outerloop progress. + It takes the outerloop iteration, latest VQE results, and latest outer loop results as arguments. + This can be used to save the most recent results to a file after each outer loop iteration. + partial_unitary_random_perturbation: A float representing the standard deviation of a normal distribution from which + the matrix elements of a perturbation matrix will be sampled. This perturbation matrix will be added to the initial + partial unitary (and re-orthonormalized) at the beginning of each orbital optimization run. + minimum_eigensolver_random_perturbation: A float representation the standard deviation of a normal distribution + from which the elements of a random set of parameters will be sampled. This perturbation will be added + to the initial parameters of the eigensolver at each outer loop iteration. + RDM_ops_batchsize: The number of Pauli ops to store in an :class:`Estimator` at any given time before deleting + the Estimator and replacing it with a blank copy. Increasing this number will lead to increased memory consumption. Setting + this number to be too low will hinder runtime performance. + """ + + super().__init__(problem=problem, + integral_tensors=integral_tensors, + num_spin_orbitals=num_spin_orbitals, + ground_state_solver=ground_state_solver, + mapper=mapper, + estimator=estimator, + partial_unitary_optimizer=partial_unitary_optimizer, + initial_partial_unitary=initial_partial_unitary, + maxiter=maxiter, + stopping_tolerance=stopping_tolerance, + spin_conserving=spin_conserving, + wavefuntion_real=wavefuntion_real, + outer_loop_callback=outer_loop_callback, + partial_unitary_random_perturbation=partial_unitary_random_perturbation, + RDM_ops_batchsize=RDM_ops_batchsize) + + if ground_state_solver.__class__.__name__ != 'VQE': + + raise AlgorithmError(f"The ground state solver needs to be of type VQE, not {ground_state_solver.__class__.__name__}") + + self.minimum_eigensolver_random_perturbation = minimum_eigensolver_random_perturbation + + def parameter_update_rule(self, result: OptOrbMinimumEigensolverResult, + iteration: int): + + if self.minimum_eigensolver_random_perturbation == None or self.minimum_eigensolver_random_perturbation == 0.0: + + self._ground_state_solver_list[iteration].initial_point = result.optimal_point + else: + + self._ground_state_solver_list[iteration].initial_point = result.optimal_point + np.random.normal(loc=0.0, + scale=self.minimum_eigensolver_random_perturbation, size=result.optimal_point.size) + +class OptOrbVQEResult(OptOrbMinimumEigensolverResult): + + def __init__(self) -> None: + super().__init__() + diff --git a/electronic_structure_algorithms/orbital_optimization/partial_unitary_projection_optimizer.py b/electronic_structure_algorithms/orbital_optimization/partial_unitary_projection_optimizer.py new file mode 100644 index 0000000..3e9e95d --- /dev/null +++ b/electronic_structure_algorithms/orbital_optimization/partial_unitary_projection_optimizer.py @@ -0,0 +1,351 @@ +import numpy as np +from time import perf_counter +from typing import Optional, Callable, Tuple +import torch +from functools import partial + +class PartialUnitaryProjectionOptimizer(): + + r"""An implementation a gradient projection optimizer with alternating Barzilai-Borwein stepsize. + See https://epubs.siam.org/doi/10.1137/16M1098759 for reference. This algorithm optimizes an objective + function over the space of MxN partial unitary matrices. In this instance we are interested in minimizing + functions involving quantities measured on a quantum computer relevant to the electronic structure problem. + e.g. the expectation values of observables.""" + + def __init__(self, + initial_BBstepsize: float, + stopping_tolerance: float, + maxiter: int, + callback: Optional[Callable] = None, + decay_factor: int = 0.8, + gradient_method: Optional[str] = 'autograd', + device: Optional[str] = 'cpu' + ) -> None: + """ + Args: + initial_BBstepsize: The initial stepsize to be used in the optimization. + stopping_tolerance: The stopping tolerance that determines when to end the optimization. + maxiter: The maximum number of optimization iterations. + callback: A callback function to track the progress of the optimization. This callback takes the iteration number and + the value of the objective function as arguments. + decay_factor: The decay factor used in constructing the stopping condition: + S_t = (1 - decay_factor)*abs(delta E_t) + decay_factor*S_t-1 + wavefunction_real: A boolean that flags whether or not we can assume the RDM elements to be real. + gradient_method: One of "finite_difference", "autograd". Setting to "finite_difference" tells the optimization + to use a finite difference approximation to calculate the gradient. Setting to "autograd" uses + automatic differentiation functionality from Pytorch, which is much faster but can possibly use + more memory. + device: Which device type to run the optimization on. e.g. setting to 'cpu' runs the optimization + on the CPU. Setting to 'cuda:n' runs one the nth GPU in a GPU array. + See https://pytorch.org/docs/stable/tensor_attributes.html#torch.torch.device for complete documentation. + """ + self._callback = callback + self.stopping_tolerance = stopping_tolerance + self.maxiter = maxiter + self._BBstepsize = initial_BBstepsize + self.decay_factor = decay_factor + self.device = device + self.gradient_method = gradient_method + + @property + def callback(self) -> Callable: + """Returns the callback function.""" + return self._callback + + @callback.setter + def callback(self, func: Callable) -> None: + """Sets the callback function.""" + self._callback = func + + @property + def BBstepsize(self) -> float: + """Returns the current BB step size.""" + return self._BBstepsize + + @BBstepsize.setter + def BBstepsize(self, stepsize: float) -> None: + """Sets the current BB step size.""" + self._BBstepsize = stepsize + + def orth(self, V: torch.Tensor) -> torch.Tensor: + """ + Generate the orthonormal projection of the matrix V. + + Args: + V: The matrix to be orthonormalized. + + Returns: + orth(V), the orthogonal projection of the matrix V. + """ + L, Q = torch.linalg.eigh(torch.t(V) @ V) + result = V @ Q @ (torch.float_power(torch.inverse(torch.diag(L)), 0.5)) @ torch.t(Q).double() + + return result + + def compute_rotated_energy_automatic_gradient(self, partial_unitary: torch.Tensor, + func: Callable) -> torch.Tensor: + + """ + Calculates the gradient of the energy function for a given partial unitary matrix U using Pytorch automatic differentiation. + + Args: + partial_unitary: The partial unitary matrix U at which the gradient is being evaluated. + + Returns: The gradient at the point U. + """ + + partial_unitary.requires_grad = True + + rotated_gradient = torch.autograd.grad([func(partial_unitary)], inputs=[partial_unitary])[0] + + partial_unitary.requires_grad = False + + return rotated_gradient + + def compute_rotated_energy_gradient(self, partial_unitary: torch.Tensor, + func: Callable) -> torch.Tensor: + + """ + Calculates the gradient of the energy functional for a given partial unitary matrix U using a finite difference approximation. + + Args: + partial_unitary: The partial unitary matrix U at which the gradient is being evaluated. + + Returns: + The gradient at the point U. + """ + + point = partial_unitary + point_size = point.size() + rotated_gradient = torch.empty(point_size, dtype=torch.float64, device=self.device) + for i in range(point_size[0]): + for j in range(point_size[1]): + temp_tensor = torch.zeros(point_size[0], point_size[1], dtype=torch.float64, device=self.device) + temp_tensor[i,j] = 10**-8 + rotated_gradient[i,j] = (func(point + temp_tensor) - func(point - temp_tensor))/(2*10**-8) + + return rotated_gradient + + def compute_updated_partial_unitary(self, iteration_number: int, + current_partial_unitary: torch.Tensor, + previous_partial_unitary: torch.Tensor, + current_rotated_energy_gradient: torch.Tensor, + previous_rotated_energy_gradient: torch.Tensor) -> torch.Tensor: + + """ + Computes the updated partial unitary using the stored values of the current and previous unitaries and gradients. + + Returns: The updated partial unitary matrix. + """ + + if iteration_number == 0: + pass + if iteration_number % 2 != 0: # iteration number is odd + + delta_U = current_partial_unitary - previous_partial_unitary + delta_G = current_rotated_energy_gradient - previous_rotated_energy_gradient + + self._BBstepsize = (torch.trace(torch.t(delta_U) @ delta_U))/torch.abs(torch.trace(torch.t(delta_U) @ delta_G)) + + if iteration_number % 2 == 0 and iteration_number != 0: # iteration number is even, but not zero. + + delta_U = current_partial_unitary - previous_partial_unitary + delta_G = current_rotated_energy_gradient - previous_rotated_energy_gradient + + self._BBstepsize = torch.abs(torch.trace(torch.t(delta_U) @ delta_G))/(torch.trace(torch.t(delta_G) @ delta_G)) + + updated_partial_unitary = self.orth(current_partial_unitary - torch.mul(current_rotated_energy_gradient, self._BBstepsize)) + + return updated_partial_unitary + + def compute_optimal_rotation(self, fun: Callable, + initial_partial_unitary: torch.Tensor, + oneRDM: torch.Tensor, + twoRDM: torch.Tensor, + one_body_integrals: torch.Tensor, + two_body_integrals: torch.Tensor) -> Tuple[torch.Tensor, float]: + + """ + Performs the optimization with fixed wavefunction and varied partial unitary U in order + to find the optimal rotation and energy. The first few iterations are written out explicitly + in order to first fill the arrays Pt_array and St_array with values that are not of type ``None``. + Then the rest of the optimization is handled in a while loop. + + Returns: A tuple of the form optimal_partial_unitary, optimal_energy. + """ + objfunc = partial(fun, oneRDM=oneRDM, twoRDM=twoRDM, one_body_integrals=one_body_integrals, two_body_integrals=two_body_integrals) + P4_array = np.array([None, None, None]) + St_array = np.array([None, 1.5*self.stopping_tolerance]) + current_partial_unitary = initial_partial_unitary.to(self.device) + previous_partial_unitary = None + current_rotated_energy_gradient = None + previous_rotated_energy_gradient = None + + iteration_number = 0 + + if self.gradient_method == 'finite_difference': + gradient_function = self.compute_rotated_energy_gradient + if self.gradient_method == 'autograd': + gradient_function = self.compute_rotated_energy_automatic_gradient + + # calculates f(U_0) and stores it in self._P4_array[2] + P4_array[2] = objfunc(partial_unitary=current_partial_unitary) + if self._callback is not None: + self._callback(iteration_number, P4_array[2].item()) + + # computes gradf(U_0) and stores it in self._current_rotated_energy_gradient + current_rotated_energy_gradient = gradient_function(current_partial_unitary, objfunc) + + # calculates U_1 (from U_0 and gradf(U_0)) and stores it in temporary variable + new_partial_unitary = self.compute_updated_partial_unitary(iteration_number=iteration_number, + current_partial_unitary=current_partial_unitary, + previous_partial_unitary=previous_partial_unitary, + current_rotated_energy_gradient=current_rotated_energy_gradient, + previous_rotated_energy_gradient=previous_rotated_energy_gradient) + + # calculates gradf(U_1) + new_rotated_energy_gradient = gradient_function(new_partial_unitary, objfunc) + + # copies U_0 into self._previous_partial_unitary for k=1 iteration + previous_partial_unitary = current_partial_unitary + + # copies gradf(U_0) into self._previous_rotated_energy_gradient for k=1 iteration + previous_rotated_energy_gradient = current_rotated_energy_gradient + + # sets the value of self._current_partial_unitary to U_1 for k=1 iteration + current_partial_unitary = new_partial_unitary + + # sets the value of self._current_rotated_energy_gradient to gradf(U_1) for k=1 iteration + current_rotated_energy_gradient = new_rotated_energy_gradient + + # updates iteration number to 1 + iteration_number += 1 + + # calculates f(U_1) and stores it in self._P4_array[1]. + P4_array[1] = objfunc(partial_unitary=current_partial_unitary) + if self._callback is not None: + self._callback(iteration_number, P4_array[1].item()) + + # calculate S_1 and store it in self._St_array[0] + St_array[0] = (1 - self.decay_factor)*torch.abs(P4_array[1] - P4_array[2]) + self.decay_factor*St_array[1] + + # calculates U_2 (from U_1 and gradf(U_1)) and stores it in temporary variable + new_partial_unitary = self.compute_updated_partial_unitary(iteration_number=iteration_number, + current_partial_unitary=current_partial_unitary, + previous_partial_unitary=previous_partial_unitary, + current_rotated_energy_gradient=current_rotated_energy_gradient, + previous_rotated_energy_gradient=previous_rotated_energy_gradient) + + # calculates gradf(U_1) and stores it in temporary variable + new_rotated_energy_gradient = gradient_function(new_partial_unitary, objfunc) + + # copies U_1 into self._previous_partial_unitary for k=2 iteration + previous_partial_unitary = current_partial_unitary + + # copies gradf(U_1) into self._previous_rotated_energy_gradient for k=2 iteration + previous_rotated_energy_gradient = current_rotated_energy_gradient + + # override value of self._current_partial_unitary to make it U_2 + current_partial_unitary = new_partial_unitary + + # override value of self._current_rotated_energy_gradient to make it gradf(U_2) + current_rotated_energy_gradient = new_rotated_energy_gradient + + # update iteration number to make it 2 + iteration_number += 1 + + + # calculates and stores f(U_2) in P4_array[0]. + P4_array[0] = objfunc(partial_unitary=current_partial_unitary) + if self._callback is not None: + self._callback(iteration_number, P4_array[0].item()) + + + # roll self._St_array so that S_1 is stored in self._St_array[1] + St_array = np.roll(St_array, shift = 1) + + # calculate S_2 and store it in self._St_array[0] + St_array[0] = (1 - self.decay_factor)*torch.abs(P4_array[0] - P4_array[1]) + self.decay_factor*St_array[1] + + # calculates U_3 (from U_2 and gradf(U_2)) and stores it in temporary variable + new_partial_unitary = self.compute_updated_partial_unitary(iteration_number=iteration_number, + current_partial_unitary=current_partial_unitary, + previous_partial_unitary=previous_partial_unitary, + current_rotated_energy_gradient=current_rotated_energy_gradient, + previous_rotated_energy_gradient=previous_rotated_energy_gradient) + + # calculates gradf(U_2) and stores it in temporary variable + new_rotated_energy_gradient = gradient_function(new_partial_unitary, objfunc) + + # copies U_2 into self._previous_partial_unitary for k=3 iteration + previous_partial_unitary = current_partial_unitary + + # copies gradf(U_2) into self._previous_rotated_energy_gradient for k=3 iteration + previous_rotated_energy_gradient = current_rotated_energy_gradient + + # override value of self._current_partial_unitary to make it U_3 + current_partial_unitary = new_partial_unitary + + # override value of self._current_rotated_energy_gradient to make it gradf(U_3) + current_rotated_energy_gradient = new_rotated_energy_gradient + + # update iteration number to make it 3 + iteration_number += 1 + + # if everything is done properly: + # self._St_array = [S_2, S_1] + # self._P4_array = [f(U_2), f(U_1), f(U_0)] + # self._current_partial_unitary = U_3 + # self._previous_partial_unitary = U_2 + # self._current_rotated_energy_gradient = gradf(U_3) + # self._previous_rotated_energy_gradient = gradf(U_2) + # everything should now be properly initialized to proceed with general iteration k loop + + while St_array[0] > self.stopping_tolerance and iteration_number <= self.maxiter: + + # roll self._P4_array so that f(U_k-1) is stored in self._P4_array[1] and f(U_k-2) is stored in self._P4_array[2] + P4_array = np.roll(P4_array, shift = 1) + + # set self._P4_array[0] to f(U_k) + P4_array[0] = objfunc(partial_unitary=current_partial_unitary) + + if self._callback is not None: + self._callback(iteration_number, P4_array[1].item()) + + # self._P4_array should now be [f(U_k), f(U_k-1), f(U_k-2)] + # roll self._St_array so that S_k-1 is stored in self._St_array[0] + St_array = np.roll(St_array, shift = 1) + + # set self._St_array[0] to S_k + St_array[0] = (1 - self.decay_factor)*torch.abs(P4_array[1] - P4_array[2]) + self.decay_factor*St_array[1] + + # self._St_array should now be [S_k, S_k-1] + # calculates U_k+1 (from U_k and gradf(U_k)) and stores it in temporary variable + new_partial_unitary = self.compute_updated_partial_unitary(iteration_number=iteration_number, + current_partial_unitary=current_partial_unitary, + previous_partial_unitary=previous_partial_unitary, + current_rotated_energy_gradient=current_rotated_energy_gradient, + previous_rotated_energy_gradient=previous_rotated_energy_gradient) + + # calculates gradf(U_k+1) and stores it in temporary variable + new_rotated_energy_gradient = gradient_function(new_partial_unitary, objfunc) + + # copies U_k into self._previous_partial_unitary for k=k+1 iteration + previous_partial_unitary = current_partial_unitary + + # copies gradf(U_k) into self._previous_rotated_energy_gradient for k=k+1 iteration + previous_rotated_energy_gradient = current_rotated_energy_gradient + + # override value of self._current_partial_unitary to make it U_K+1 + current_partial_unitary = new_partial_unitary + + # override value of self._current_rotated_energy_gradient to make it gradf(U_K+1) + current_rotated_energy_gradient = new_rotated_energy_gradient + + # update iteration number to make it k+1 + iteration_number += 1 + + current_partial_unitary = current_partial_unitary.to(torch.device('cpu')) + + return current_partial_unitary, P4_array[0] + \ No newline at end of file diff --git a/electronic_structure_algorithms/partial_unitary_initializations.py/__init__.py b/electronic_structure_algorithms/partial_unitary_initializations.py/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/examples/H2_OptOrbAdaptVQE.py b/examples/H2_OptOrbAdaptVQE.py new file mode 100644 index 0000000..ec4871b --- /dev/null +++ b/examples/H2_OptOrbAdaptVQE.py @@ -0,0 +1,95 @@ +import numpy as np +from qiskit_nature.second_q.drivers import PySCFDriver +from qiskit_nature.units import DistanceUnit +from qiskit_nature.second_q.mappers import JordanWignerMapper +from qiskit_algorithms.optimizers import L_BFGS_B, COBYLA +from qiskit_aer.primitives import Estimator +from qiskit_nature.second_q.circuit.library import HartreeFock, UCCSD + +from qiskit_algorithms import AdaptVQE, VQE + +from electronic_structure_algorithms.orbital_optimization import PartialUnitaryProjectionOptimizer, OptOrbAdaptVQE + +from time import perf_counter + +estimator = Estimator(approximation=True, + run_options={'shots': 10**30}) +mapper=JordanWignerMapper() + +driver = PySCFDriver(atom=f'H 0 0 0; H 0 0 {0.735}', + charge=0, + spin=0, + unit=DistanceUnit.ANGSTROM, + basis='cc-pVTZ') + +q_molecule = driver.run() +num_particles = q_molecule.num_particles + +l_bfgs_b = L_BFGS_B(maxfun=10**6, maxiter=10**6) +cobyla = COBYLA(maxiter=10**6) + +num_reduced_qubits = 4 + +HF_state = HartreeFock(qubit_mapper=mapper, + num_spatial_orbitals=int(num_reduced_qubits/2), + num_particles=num_particles) + +ansatz = UCCSD(qubit_mapper=mapper, + num_spatial_orbitals=int(num_reduced_qubits/2), + num_particles=num_particles, + initial_state=HF_state) + + +outer_iteration = 0 +vqe_start_time = perf_counter() +def vqe_callback(eval_count, parameters, mean, std): + global vqe_start_time + print(f'Outer loop iteration: {outer_iteration}, function evaluation: {eval_count}, energy: {mean}, time = {perf_counter() - vqe_start_time}') + + vqe_start_time = perf_counter() + + +orbital_rotation_start_time = perf_counter() +def orbital_rotation_callback(orbital_rotation_iteration, energy): + global orbital_rotation_start_time + print(f'Outer loop iteration: {outer_iteration}, Iteration: {orbital_rotation_iteration}, energy: {energy}, time: {perf_counter() - orbital_rotation_start_time}') + orbital_rotation_start_time = perf_counter() + + +def outer_loop_callback(optorb_iteration, vqe_result, optorb_result): + global outer_iteration + outer_iteration += 1 + + +partial_unitary_optimizer = PartialUnitaryProjectionOptimizer(initial_BBstepsize=10**-3, + stopping_tolerance=10**-5, + maxiter=10000, + gradient_method='autograd', + callback=orbital_rotation_callback) + +vqe_instance = VQE(ansatz=ansatz, + initial_point=np.zeros(ansatz.num_parameters), + optimizer=l_bfgs_b, + estimator=estimator, + callback=vqe_callback) + +adapt_vqe_instance = AdaptVQE(solver=vqe_instance) +adapt_vqe_instance.supports_aux_operators = lambda: True + +optorbvqe_instance = OptOrbAdaptVQE(problem=q_molecule, + integral_tensors = None, + num_spin_orbitals=num_reduced_qubits, + ground_state_solver=adapt_vqe_instance, + mapper=mapper, + estimator=estimator, + partial_unitary_optimizer=partial_unitary_optimizer, + maxiter=20, + wavefuntion_real=True, + spin_conserving=True, + stopping_tolerance=10**-5, + outer_loop_callback=outer_loop_callback, + partial_unitary_random_perturbation=0.01) + +ground_state_energy_result = optorbvqe_instance.compute_minimum_energy() +print(ground_state_energy_result) + diff --git a/examples/H2_OptOrbMCVQE.py b/examples/H2_OptOrbMCVQE.py new file mode 100644 index 0000000..15a59fd --- /dev/null +++ b/examples/H2_OptOrbMCVQE.py @@ -0,0 +1,89 @@ +import numpy as np +from qiskit_nature.second_q.drivers import PySCFDriver +from qiskit_nature.units import DistanceUnit +from qiskit_nature.second_q.mappers import JordanWignerMapper +from qiskit_algorithms.optimizers import L_BFGS_B, COBYLA +from qiskit_aer.primitives import Estimator +from qiskit_nature.second_q.circuit.library import UCCSD + +from electronic_structure_algorithms.orbital_optimization import PartialUnitaryProjectionOptimizer, OptOrbMCVQE +from electronic_structure_algorithms.excited_states_eigensolvers import MCVQE + +from time import perf_counter + +estimator = Estimator(approximation=True) +mapper=JordanWignerMapper() + +driver = PySCFDriver(atom=f'H 0 0 0; H 0 0 {0.735}', + charge=0, + spin=0, + unit=DistanceUnit.ANGSTROM, + basis='cc-pVTZ') + +q_molecule = driver.run() +num_particles = q_molecule.num_particles + +l_bfgs_b = L_BFGS_B(maxfun=10**6, maxiter=10**6) + +num_reduced_qubits = 4 + +ansatz = UCCSD(qubit_mapper=mapper, + num_spatial_orbitals=int(num_reduced_qubits/2), + num_particles=num_particles, + reps=2) + + +outer_iteration = 0 +vqe_start_time = perf_counter() +def mcvqe_callback(eval_count, parameters, mean, std): + global vqe_start_time + print(f'Outer loop iteration: {outer_iteration}, function evaluation: {eval_count}, energy: {mean}, time = {perf_counter() - vqe_start_time}') + + vqe_start_time = perf_counter() + + +orbital_rotation_start_time = perf_counter() +def orbital_rotation_callback(orbital_rotation_iteration, energy): + global orbital_rotation_start_time + print(f'Outer loop iteration: {outer_iteration}, Iteration: {orbital_rotation_iteration}, energy: {energy}, time: {perf_counter() - orbital_rotation_start_time}') + orbital_rotation_start_time = perf_counter() + + +def outer_loop_callback(optorb_iteration, vqe_result, optorb_result): + global outer_iteration + outer_iteration += 1 + + +partial_unitary_optimizer = PartialUnitaryProjectionOptimizer(initial_BBstepsize=10**-3, + stopping_tolerance=10**-5, + maxiter=10000, + gradient_method='autograd', + callback=orbital_rotation_callback) + +mcvqe_instance = MCVQE(k=2, + excitations='s', + num_particles=q_molecule.num_particles, + ansatz=ansatz, + initial_point=np.zeros(ansatz.num_parameters), + optimizer=l_bfgs_b, + estimator=estimator, + callback=mcvqe_callback) + +optorbvqe_instance = OptOrbMCVQE(problem=q_molecule, + integral_tensors = None, + num_spin_orbitals=num_reduced_qubits, + excited_states_solver=mcvqe_instance, + mapper=mapper, + estimator=estimator, + partial_unitary_optimizer=partial_unitary_optimizer, + maxiter=20, + wavefuntion_real=True, + spin_conserving=True, + stopping_tolerance=10**-5, + outer_loop_callback=outer_loop_callback, + partial_unitary_random_perturbation=0.01, + eigensolver_random_perturbation=0.0) + +ground_state_energy_result = optorbvqe_instance.compute_energies() +print(ground_state_energy_result) + diff --git a/examples/H2_OptOrbSSVQE.py b/examples/H2_OptOrbSSVQE.py new file mode 100644 index 0000000..a6dac0f --- /dev/null +++ b/examples/H2_OptOrbSSVQE.py @@ -0,0 +1,123 @@ +import numpy as np +import torch +from qiskit import QuantumCircuit +from qiskit_nature.second_q.drivers import PySCFDriver +from qiskit_nature.units import DistanceUnit +from qiskit_nature.second_q.mappers import JordanWignerMapper +from qiskit_algorithms.optimizers import L_BFGS_B, COBYLA +from qiskit_aer.primitives import Estimator +from qiskit_nature.second_q.circuit.library import UCCSD +from qiskit_nature.second_q.operators.tensor_ordering import to_physicist_ordering + +from electronic_structure_algorithms.orbital_optimization import PartialUnitaryProjectionOptimizer, OptOrbSSVQE +from electronic_structure_algorithms.excited_states_eigensolvers import SSVQE +from electronic_structure_algorithms.initializations import get_HF_permutation_matrix, get_CIS_states + +from time import perf_counter + +estimator = Estimator(approximation=True) +mapper=JordanWignerMapper() + +driver = PySCFDriver(atom=f'H 0 0 0; H 0 0 {0.735}', + charge=0, + spin=0, + unit=DistanceUnit.ANGSTROM, + basis='cc-pVTZ') + +q_molecule = driver.run() +num_particles = q_molecule.num_particles +num_original_spin_orbitals = q_molecule.num_spin_orbitals + +l_bfgs_b = L_BFGS_B(maxfun=10**6, maxiter=10**6) +cobyla = COBYLA(maxiter=10**6) + +num_reduced_qubits = 4 +num_states = 2 + +ansatz = UCCSD(qubit_mapper=mapper, + num_spatial_orbitals=int(num_reduced_qubits/2), + num_particles=num_particles, + reps=2) + +one_body_integrals = torch.from_numpy(np.asarray(q_molecule.hamiltonian.electronic_integrals.second_q_coeffs()["+-"].to_dense())) +two_body_integrals = torch.from_numpy(np.asarray(-1*to_physicist_ordering(q_molecule.hamiltonian.electronic_integrals.second_q_coeffs()["++--"].to_dense()))) + +initial_partial_unitary = get_HF_permutation_matrix(num_original_spin_orbitals=num_original_spin_orbitals, + num_spin_orbitals=num_reduced_qubits) + +transformed_one_body_integrals = torch.einsum('pq,pi,qj->ij', + one_body_integrals, + torch.block_diag(initial_partial_unitary, initial_partial_unitary), + torch.block_diag(initial_partial_unitary, initial_partial_unitary)).detach().numpy() + +transformed_two_body_integrals = torch.einsum('pqrs,pi,qj,rk,sl->ijkl', + two_body_integrals, + torch.block_diag(initial_partial_unitary, initial_partial_unitary), + torch.block_diag(initial_partial_unitary, initial_partial_unitary), + torch.block_diag(initial_partial_unitary, initial_partial_unitary), + torch.block_diag(initial_partial_unitary, initial_partial_unitary)).detach().numpy() + +cis_states = get_CIS_states(one_body_integrals=transformed_one_body_integrals, + two_body_integrals=transformed_two_body_integrals, + num_particles=q_molecule.num_particles, + state_representation='dense')[:num_states] + +initial_states = [QuantumCircuit(ansatz.num_qubits) for n in range(num_states)] +for n in range(num_states): + initial_states[n].initialize(cis_states[n], initial_states[n].qubits) + +outer_iteration = 0 +vqe_start_time = perf_counter() +def ssvqe_callback(eval_count, parameters, mean, std): + global vqe_start_time + print(f'Outer loop iteration: {outer_iteration}, function evaluation: {eval_count}, energy: {mean}, time = {perf_counter() - vqe_start_time}') + + vqe_start_time = perf_counter() + + +orbital_rotation_start_time = perf_counter() +def orbital_rotation_callback(orbital_rotation_iteration, energy): + global orbital_rotation_start_time + print(f'Outer loop iteration: {outer_iteration}, Iteration: {orbital_rotation_iteration}, energy sum: {energy}, time: {perf_counter() - orbital_rotation_start_time}') + orbital_rotation_start_time = perf_counter() + + +def outer_loop_callback(optorb_iteration, vqe_result, optorb_result): + global outer_iteration + outer_iteration += 1 + + +partial_unitary_optimizer = PartialUnitaryProjectionOptimizer(initial_BBstepsize=10**-3, + stopping_tolerance=10**-5, + maxiter=10000, + gradient_method='autograd', + callback=orbital_rotation_callback) + +ssvqe_instance = SSVQE(k=num_states, + initial_states=initial_states, + ansatz=ansatz, + initial_point=np.zeros(ansatz.num_parameters), + optimizer=l_bfgs_b, + estimator=estimator, + callback=ssvqe_callback, + weight_vector=[2,1]) + +optorbssvqe_instance = OptOrbSSVQE(problem=q_molecule, + integral_tensors=(one_body_integrals, two_body_integrals), + num_spin_orbitals=num_reduced_qubits, + excited_states_solver=ssvqe_instance, + mapper=mapper, + estimator=estimator, + partial_unitary_optimizer=partial_unitary_optimizer, + initial_partial_unitary=initial_partial_unitary, + maxiter=20, + wavefuntion_real=True, + spin_conserving=True, + stopping_tolerance=10**-5, + outer_loop_callback=outer_loop_callback, + partial_unitary_random_perturbation=0.01, + eigensolver_random_perturbation=0.0) + +ground_state_energy_result = optorbssvqe_instance.compute_energies() +print(ground_state_energy_result) + diff --git a/examples/H2_OptOrbVQD.py b/examples/H2_OptOrbVQD.py new file mode 100644 index 0000000..f83254f --- /dev/null +++ b/examples/H2_OptOrbVQD.py @@ -0,0 +1,136 @@ +import numpy as np +import torch +from qiskit import QuantumCircuit, transpile +from qiskit.algorithms.state_fidelities import ComputeUncompute +from qiskit_nature.second_q.drivers import PySCFDriver +from qiskit_nature.units import DistanceUnit +from qiskit_nature.second_q.mappers import JordanWignerMapper +from qiskit_algorithms.optimizers import L_BFGS_B, COBYLA +from qiskit_aer.primitives import Estimator, Sampler +from qiskit_nature.second_q.circuit.library import UCCSD, HartreeFock +from qiskit_nature.second_q.operators.tensor_ordering import to_physicist_ordering + +from electronic_structure_algorithms.orbital_optimization import PartialUnitaryProjectionOptimizer, OptOrbVQD +from electronic_structure_algorithms.excited_states_eigensolvers import VQD +from electronic_structure_algorithms.initializations import get_HF_permutation_matrix, get_CIS_states + +from time import perf_counter + +estimator = Estimator(approximation=True) +sampler = Sampler(run_options={'shots': None}) +fidelity = ComputeUncompute(sampler=sampler) +mapper=JordanWignerMapper() + +driver = PySCFDriver(atom=f'H 0 0 0; H 0 0 {0.735}', + charge=0, + spin=0, + unit=DistanceUnit.ANGSTROM, + basis='cc-pVTZ') + +q_molecule = driver.run() +num_particles = q_molecule.num_particles +num_original_spin_orbitals = q_molecule.num_spin_orbitals + +l_bfgs_b = L_BFGS_B(maxfun=10**6, maxiter=10**6) +cobyla = COBYLA(maxiter=10**6) + +num_reduced_qubits = 4 +num_states = 2 + +ansatz = UCCSD(qubit_mapper=mapper, + num_spatial_orbitals=int(num_reduced_qubits/2), + num_particles=num_particles, + reps=2) + +one_body_integrals = torch.from_numpy(np.asarray(q_molecule.hamiltonian.electronic_integrals.second_q_coeffs()["+-"].to_dense())) +two_body_integrals = torch.from_numpy(np.asarray(-1*to_physicist_ordering(q_molecule.hamiltonian.electronic_integrals.second_q_coeffs()["++--"].to_dense()))) + +initial_partial_unitary = get_HF_permutation_matrix(num_original_spin_orbitals=num_original_spin_orbitals, + num_spin_orbitals=num_reduced_qubits) + +transformed_one_body_integrals = torch.einsum('pq,pi,qj->ij', + one_body_integrals, + torch.block_diag(initial_partial_unitary, initial_partial_unitary), + torch.block_diag(initial_partial_unitary, initial_partial_unitary)).detach().numpy() + +transformed_two_body_integrals = torch.einsum('pqrs,pi,qj,rk,sl->ijkl', + two_body_integrals, + torch.block_diag(initial_partial_unitary, initial_partial_unitary), + torch.block_diag(initial_partial_unitary, initial_partial_unitary), + torch.block_diag(initial_partial_unitary, initial_partial_unitary), + torch.block_diag(initial_partial_unitary, initial_partial_unitary)).detach().numpy() + +cis_states = get_CIS_states(one_body_integrals=transformed_one_body_integrals, + two_body_integrals=transformed_two_body_integrals, + num_particles=q_molecule.num_particles, + state_representation='dense')[:num_states] + +initial_states = [QuantumCircuit(num_reduced_qubits) for n in range(num_states)] +for n in range(num_states): + initial_states[n].initialize(cis_states[n], initial_states[n].qubits) + initial_states[n] = transpile(initial_states[n], optimization_level=3).decompose() + +for state in initial_states: + + for n, instruction in enumerate(state.data): + + if instruction.operation.name == 'reset': + del state.data[n] + +ansatz_list = [UCCSD(qubit_mapper=mapper, + num_spatial_orbitals=int(num_reduced_qubits/2), + num_particles=num_particles, + reps=2, + initial_state=initial_states[n]) for n in range(num_states)] + +outer_iteration = 0 +vqd_start_time = perf_counter() +def vqd_callback(eval_count, parameters, mean, std, step, costs): + global vqd_start_time + print(f'OptOrb iteration: {outer_iteration}, function evaluation: {eval_count}, energy level {step - 1}: {mean}, time = {perf_counter() - vqd_start_time}') + vqd_start_time = perf_counter() + +orbital_rotation_start_time = perf_counter() +def orbital_rotation_callback(orbital_rotation_iteration, energy): + global orbital_rotation_start_time + print(f'Outer loop iteration: {outer_iteration}, Iteration: {orbital_rotation_iteration}, energy sum: {energy}, time: {perf_counter() - orbital_rotation_start_time}') + orbital_rotation_start_time = perf_counter() + + +def outer_loop_callback(optorb_iteration, vqe_result, optorb_result): + global outer_iteration + outer_iteration += 1 + + +partial_unitary_optimizer = PartialUnitaryProjectionOptimizer(initial_BBstepsize=10**-3, + stopping_tolerance=10**-5, + maxiter=10000, + gradient_method='autograd', + callback=orbital_rotation_callback) + +vqd_instance = VQD(k=num_states, + ansatz=ansatz_list, + initial_point=[np.zeros(ansatz_list[n].num_parameters) for n in range(num_states)], + optimizer=l_bfgs_b, + estimator=estimator, + callback=vqd_callback, + fidelity=fidelity, + betas=[2,2]) + +optorbssvqe_instance = OptOrbVQD(problem=q_molecule, + integral_tensors=(one_body_integrals, two_body_integrals), + num_spin_orbitals=num_reduced_qubits, + excited_states_solver=vqd_instance, + mapper=mapper, + estimator=estimator, + partial_unitary_optimizer=partial_unitary_optimizer, + maxiter=20, + wavefuntion_real=True, + spin_conserving=True, + stopping_tolerance=10**-5, + outer_loop_callback=outer_loop_callback, + partial_unitary_random_perturbation=0.01, + eigensolver_random_perturbation=0.0) + +ground_state_energy_result = optorbssvqe_instance.compute_energies() +print(ground_state_energy_result) \ No newline at end of file diff --git a/examples/H2_OptOrbVQE.py b/examples/H2_OptOrbVQE.py new file mode 100644 index 0000000..25f9568 --- /dev/null +++ b/examples/H2_OptOrbVQE.py @@ -0,0 +1,91 @@ +import numpy as np +from qiskit_nature.second_q.drivers import PySCFDriver +from qiskit_nature.units import DistanceUnit +from qiskit_nature.second_q.mappers import JordanWignerMapper +from qiskit_algorithms.minimum_eigensolvers import VQE +from qiskit_algorithms.optimizers import L_BFGS_B, COBYLA +from qiskit_aer.primitives import Estimator +from qiskit_nature.second_q.circuit.library import HartreeFock, UCCSD + +from electronic_structure_algorithms.orbital_optimization import PartialUnitaryProjectionOptimizer, OptOrbVQE + +from time import perf_counter + +estimator = Estimator(approximation=True) +mapper=JordanWignerMapper() + +driver = PySCFDriver(atom=f'H 0 0 0; H 0 0 {0.735}', + charge=0, + spin=0, + unit=DistanceUnit.ANGSTROM, + basis='6-31G') + +q_molecule = driver.run() +num_particles = q_molecule.num_particles + +l_bfgs_b = L_BFGS_B(maxfun=10**6, maxiter=10**6) +cobyla = COBYLA(maxiter=10**6) + +num_reduced_qubits = 4 + +HF_state = HartreeFock(qubit_mapper=mapper, + num_spatial_orbitals=int(num_reduced_qubits/2), + num_particles=num_particles) + +ansatz = UCCSD(qubit_mapper=mapper, + num_spatial_orbitals=int(num_reduced_qubits/2), + num_particles=num_particles, + initial_state=HF_state) + + +outer_iteration = 0 +vqe_start_time = perf_counter() +def vqe_callback(eval_count, parameters, mean, std): + global vqe_start_time + print(f'Outer loop iteration: {outer_iteration}, function evaluation: {eval_count}, energy: {mean}, time = {perf_counter() - vqe_start_time}') + + vqe_start_time = perf_counter() + + +orbital_rotation_start_time = perf_counter() +def orbital_rotation_callback(orbital_rotation_iteration, energy): + global orbital_rotation_start_time + print(f'Outer loop iteration: {outer_iteration}, Iteration: {orbital_rotation_iteration}, energy: {energy}, time: {perf_counter() - orbital_rotation_start_time}') + orbital_rotation_start_time = perf_counter() + + +def outer_loop_callback(optorb_iteration, vqe_result, optorb_result): + global outer_iteration + outer_iteration += 1 + + +partial_unitary_optimizer = PartialUnitaryProjectionOptimizer(initial_BBstepsize=10**-3, + stopping_tolerance=10**-5, + maxiter=10000, + gradient_method='autograd', + callback=orbital_rotation_callback) + +vqe_instance = VQE(ansatz=ansatz, + initial_point=np.zeros(ansatz.num_parameters), + optimizer=l_bfgs_b, + estimator=estimator, + callback=vqe_callback) + +optorbvqe_instance = OptOrbVQE(problem=q_molecule, + integral_tensors = None, + num_spin_orbitals=num_reduced_qubits, + ground_state_solver=vqe_instance, + mapper=mapper, + estimator=estimator, + partial_unitary_optimizer=partial_unitary_optimizer, + maxiter=20, + wavefuntion_real=True, + spin_conserving=True, + stopping_tolerance=10**-5, + outer_loop_callback=outer_loop_callback, + partial_unitary_random_perturbation=0.01, + minimum_eigensolver_random_perturbation=0.0) + +ground_state_energy_result = optorbvqe_instance.compute_minimum_energy() +print(ground_state_energy_result) + diff --git a/examples/H4_OptOrbAdaptVQE.py b/examples/H4_OptOrbAdaptVQE.py new file mode 100644 index 0000000..b24df13 --- /dev/null +++ b/examples/H4_OptOrbAdaptVQE.py @@ -0,0 +1,95 @@ +import numpy as np +from qiskit_nature.second_q.drivers import PySCFDriver +from qiskit_nature.units import DistanceUnit +from qiskit_nature.second_q.mappers import JordanWignerMapper +from qiskit_algorithms.optimizers import L_BFGS_B, COBYLA +from qiskit_aer.primitives import Estimator +from qiskit_nature.second_q.circuit.library import HartreeFock, UCCSD + +from qiskit_algorithms import AdaptVQE, VQE + +from electronic_structure_algorithms.orbital_optimization import PartialUnitaryProjectionOptimizer, OptOrbAdaptVQE + +from time import perf_counter + +estimator = Estimator(approximation=True, + run_options={'shots': 10**30}) +mapper=JordanWignerMapper() + +interatomic_distance = 1.23 +driver = PySCFDriver(atom=f'H 0 0 0; H 0 {interatomic_distance} 0; H {interatomic_distance} 0 0; H {interatomic_distance} {interatomic_distance} 0', + charge=0, + spin=0, + unit=DistanceUnit.ANGSTROM, basis='cc-pVDZ') + +q_molecule = driver.run() +num_particles = q_molecule.num_particles + +l_bfgs_b = L_BFGS_B(maxfun=10**6, maxiter=10**6) +cobyla = COBYLA(maxiter=10**6) + +num_reduced_qubits = 8 + +HF_state = HartreeFock(qubit_mapper=mapper, + num_spatial_orbitals=int(num_reduced_qubits/2), + num_particles=num_particles) + +ansatz = UCCSD(qubit_mapper=mapper, + num_spatial_orbitals=int(num_reduced_qubits/2), + num_particles=num_particles, + initial_state=HF_state) + + +outer_iteration = 0 +vqe_start_time = perf_counter() +def vqe_callback(eval_count, parameters, mean, std): + global vqe_start_time + print(f'Outer loop iteration: {outer_iteration}, function evaluation: {eval_count}, energy: {mean}, time = {perf_counter() - vqe_start_time}') + + vqe_start_time = perf_counter() + + +orbital_rotation_start_time = perf_counter() +def orbital_rotation_callback(orbital_rotation_iteration, energy): + global orbital_rotation_start_time + print(f'Outer loop iteration: {outer_iteration}, Iteration: {orbital_rotation_iteration}, energy: {energy}, time: {perf_counter() - orbital_rotation_start_time}') + orbital_rotation_start_time = perf_counter() + + +def outer_loop_callback(optorb_iteration, vqe_result, optorb_result): + global outer_iteration + outer_iteration += 1 + + +partial_unitary_optimizer = PartialUnitaryProjectionOptimizer(initial_BBstepsize=10**-3, + stopping_tolerance=10**-5, + maxiter=10000, + gradient_method='autograd', + callback=orbital_rotation_callback) + +vqe_instance = VQE(ansatz=ansatz, + initial_point=np.zeros(ansatz.num_parameters), + optimizer=l_bfgs_b, + estimator=estimator, + callback=vqe_callback) + +adapt_vqe_instance = AdaptVQE(solver=vqe_instance) +adapt_vqe_instance.supports_aux_operators = lambda: True + +optorbvqe_instance = OptOrbAdaptVQE(problem=q_molecule, + integral_tensors = None, + num_spin_orbitals=num_reduced_qubits, + ground_state_solver=adapt_vqe_instance, + mapper=mapper, + estimator=estimator, + partial_unitary_optimizer=partial_unitary_optimizer, + maxiter=20, + wavefuntion_real=True, + spin_conserving=True, + stopping_tolerance=10**-5, + outer_loop_callback=outer_loop_callback, + partial_unitary_random_perturbation=0.01) + +ground_state_energy_result = optorbvqe_instance.compute_minimum_energy() +print(ground_state_energy_result) + diff --git a/examples/H4_OptOrbVQE.py b/examples/H4_OptOrbVQE.py new file mode 100644 index 0000000..f463c65 --- /dev/null +++ b/examples/H4_OptOrbVQE.py @@ -0,0 +1,92 @@ +import numpy as np +from qiskit_nature.second_q.drivers import PySCFDriver +from qiskit_nature.units import DistanceUnit +from qiskit_nature.second_q.mappers import JordanWignerMapper +from qiskit_algorithms.minimum_eigensolvers import VQE +from qiskit_algorithms.optimizers import L_BFGS_B, COBYLA +from qiskit_aer.primitives import Estimator +from qiskit_nature.second_q.circuit.library import HartreeFock, UCCSD + +from electronic_structure_algorithms.orbital_optimization import PartialUnitaryProjectionOptimizer, OptOrbVQE + +from time import perf_counter + +estimator = Estimator(approximation=True) +mapper=JordanWignerMapper() + +interatomic_distance = 1.23 +driver = PySCFDriver(atom=f'H 0 0 0; H 0 {interatomic_distance} 0; H {interatomic_distance} 0 0; H {interatomic_distance} {interatomic_distance} 0', + charge=0, + spin=0, + unit=DistanceUnit.ANGSTROM, basis='cc-pVDZ') + +q_molecule = driver.run() +num_particles = q_molecule.num_particles + +l_bfgs_b = L_BFGS_B(maxfun=10**6, maxiter=10**6) +cobyla = COBYLA(maxiter=10**6) + +num_reduced_qubits = 8 + +HF_state = HartreeFock(qubit_mapper=mapper, + num_spatial_orbitals=int(num_reduced_qubits/2), + num_particles=num_particles) + +ansatz = UCCSD(qubit_mapper=mapper, + num_spatial_orbitals=int(num_reduced_qubits/2), + num_particles=num_particles, + initial_state=HF_state, + reps=2) + + +outer_iteration = 0 +vqe_start_time = perf_counter() +def vqe_callback(eval_count, parameters, mean, std): + global vqe_start_time + print(f'Outer loop iteration: {outer_iteration}, function evaluation: {eval_count}, energy: {mean}, time = {perf_counter() - vqe_start_time}') + + vqe_start_time = perf_counter() + + +orbital_rotation_start_time = perf_counter() +def orbital_rotation_callback(orbital_rotation_iteration, energy): + global orbital_rotation_start_time + print(f'Outer loop iteration: {outer_iteration}, Iteration: {orbital_rotation_iteration}, energy: {energy}, time: {perf_counter() - orbital_rotation_start_time}') + orbital_rotation_start_time = perf_counter() + + +def outer_loop_callback(optorb_iteration, vqe_result, optorb_result): + global outer_iteration + outer_iteration += 1 + + +partial_unitary_optimizer = PartialUnitaryProjectionOptimizer(initial_BBstepsize=10**-3, + stopping_tolerance=10**-5, + maxiter=10000, + gradient_method='autograd', + callback=orbital_rotation_callback) + +vqe_instance = VQE(ansatz=ansatz, + initial_point=np.zeros(ansatz.num_parameters), + optimizer=l_bfgs_b, + estimator=estimator, + callback=vqe_callback) + +optorbvqe_instance = OptOrbVQE(problem=q_molecule, + integral_tensors = None, + num_spin_orbitals=num_reduced_qubits, + ground_state_solver=vqe_instance, + mapper=mapper, + estimator=estimator, + partial_unitary_optimizer=partial_unitary_optimizer, + maxiter=20, + wavefuntion_real=True, + spin_conserving=True, + stopping_tolerance=10**-5, + outer_loop_callback=outer_loop_callback, + partial_unitary_random_perturbation=0.01, + minimum_eigensolver_random_perturbation=0.0) + +ground_state_energy_result = optorbvqe_instance.compute_minimum_energy() +print(ground_state_energy_result) + diff --git a/requirements-dev.txt b/requirements-dev.txt index 291e411..0fb4706 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -2,3 +2,5 @@ coverage==5.5 pylint==2.9.5 tox==3.24.0 black==22.3.0 +ddt +qiskit_algorithms \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index e581560..b0e81c0 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,7 @@ pytest==6.2.4 qiskit -qiskit-aer +qiskit-aer>=0.12.0 +qiskit-nature>=0.7.1 +qiskit-terra>=0.45.0 +torch>=1.12.0 +qiskit_algorithms diff --git a/setup.py b/setup.py index 2aa41fe..e980952 100644 --- a/setup.py +++ b/setup.py @@ -9,10 +9,11 @@ install_requires = fp.read() setuptools.setup( - name="my_super_project", - description="My Super Project", + name="orbital-optimization", + description="A collection of quantum algorithms to solve the electronic structure problem using orbital optimization.", long_description=long_description, packages=setuptools.find_packages(), install_requires=install_requires, - python_requires='>=3.6' + python_requires='>=3.8', + version='0.1.0' ) diff --git a/tests/__init__.py b/tests/__init__.py index e69de29..201c768 100644 --- a/tests/__init__.py +++ b/tests/__init__.py @@ -0,0 +1,17 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2018, 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. + +"""Algorithms test module""" + +from .algorithms_test_case import QiskitAlgorithmsTestCase + +__all__ = ["QiskitAlgorithmsTestCase"] \ No newline at end of file diff --git a/tests/algorithms_test_case.py b/tests/algorithms_test_case.py new file mode 100644 index 0000000..0008a76 --- /dev/null +++ b/tests/algorithms_test_case.py @@ -0,0 +1,21 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2018, 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. + +"""Algorithms Test Case""" + +from qiskit.test import QiskitTestCase + + +class QiskitAlgorithmsTestCase(QiskitTestCase): + """Algorithms test Case""" + + pass \ No newline at end of file diff --git a/tests/test_mcvqe.py b/tests/test_mcvqe.py new file mode 100644 index 0000000..bd7a376 --- /dev/null +++ b/tests/test_mcvqe.py @@ -0,0 +1,490 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2022 +# +# 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 SSVQE """ + +import unittest +from .algorithms_test_case import QiskitAlgorithmsTestCase + +from functools import partial +import numpy as np +from ddt import data, ddt + +from qiskit import QuantumCircuit +from qiskit_algorithms import AlgorithmError +from qiskit_algorithms.optimizers import ( + COBYLA, + L_BFGS_B, + SLSQP, + CG, + TNC, + P_BFGS, + GradientDescent, + OptimizerResult, +) + +from qiskit_algorithms.gradients import ParamShiftEstimatorGradient +from qiskit.circuit.library import TwoLocal, RealAmplitudes +from qiskit.quantum_info import SparsePauliOp +from qiskit.primitives import Sampler, Estimator +from qiskit_algorithms.state_fidelities import ComputeUncompute +from qiskit_algorithms.utils import algorithm_globals +from qiskit.quantum_info.operators import Operator + +from electronic_structure_algorithms.excited_states_eigensolvers import MCVQE + +# pylint: disable=invalid-name +def _mock_optimizer(fun, x0, jac=None, bounds=None, inputs=None) -> OptimizerResult: + """A mock of a callable that can be used as minimizer in SSVQE.""" + result = OptimizerResult() + result.x = np.zeros_like(x0) + result.fun = fun(result.x) + result.nit = 0 + + if inputs is not None: + inputs.update({"fun": fun, "x0": x0, "jac": jac, "bounds": bounds}) + return result + + +I = SparsePauliOp.from_list([("I", 1)]) # pylint: disable=invalid-name +X = SparsePauliOp.from_list([("X", 1)]) # pylint: disable=invalid-name +Z = SparsePauliOp.from_list([("Z", 1)]) # pylint: disable=invalid-name + +H2_PAULI = ( + -1.052373245772859 * (I ^ I) + + 0.39793742484318045 * (I ^ Z) + - 0.39793742484318045 * (Z ^ I) + - 0.01128010425623538 * (Z ^ Z) + + 0.18093119978423156 * (X ^ X) +) + +H2_OP = Operator(H2_PAULI.to_matrix()) + + +@ddt +class TestMCVQE(QiskitAlgorithmsTestCase): + """Test MCVQE""" + + def setUp(self): + super().setUp() + self.seed = 50 + algorithm_globals.random_seed = self.seed + + self.h2_energy = -1.85727503 + self.h2_energy_excited = [-1.85727503, -1.24458455] + + self.ryrz_wavefunction = TwoLocal( + rotation_blocks=["ry", "rz"], entanglement_blocks="cz", reps=1 + ) + self.ry_wavefunction = TwoLocal(rotation_blocks="ry", entanglement_blocks="cz") + + self.estimator = Estimator() + self.estimator_shots = Estimator(options={"shots": 2048, "seed": self.seed}) + self.fidelity = ComputeUncompute(Sampler()) + self.betas = [50, 50] + + @data(H2_PAULI, H2_OP) + def test_basic_operator(self, op): + """Test SSVQE without aux_operators.""" + wavefunction = self.ryrz_wavefunction + ssvqe = MCVQE( + estimator=self.estimator, + ansatz=wavefunction, + optimizer=COBYLA(), + ) + + result = ssvqe.compute_eigenvalues(operator=op) + + with self.subTest(msg="test eigenvalue"): + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=1 + ) + + with self.subTest(msg="test dimension of optimal point"): + self.assertEqual(len(result.optimal_point), 8) + + with self.subTest(msg="assert cost_function_evals is set"): + self.assertIsNotNone(result.cost_function_evals) + + with self.subTest(msg="assert optimizer_time is set"): + self.assertIsNotNone(result.optimizer_time) + + @data(H2_PAULI, H2_OP) + def test_mismatching_num_qubits(self, op): + """Ensuring circuit and operator mismatch is caught""" + wavefunction = QuantumCircuit(1) + optimizer = SLSQP(maxiter=50) + ssvqe = MCVQE( + estimator=self.estimator, k=1, ansatz=wavefunction, optimizer=optimizer + ) + with self.assertRaises(AlgorithmError): + _ = ssvqe.compute_eigenvalues(operator=op) + + @data(H2_PAULI, H2_OP) + def test_initial_states_mismatching_num_qubits(self, op): + """Ensuring initial states and operator mismatch is caught""" + initial_state = [QuantumCircuit(1)] + wavefunction = QuantumCircuit(2) + optimizer = SLSQP(maxiter=50) + ssvqe = MCVQE( + estimator=self.estimator, + k=1, + ansatz=wavefunction, + optimizer=optimizer, + initial_states=initial_state, + ) + with self.assertRaises(AlgorithmError): + _ = ssvqe.compute_eigenvalues(operator=op) + + @data(H2_PAULI, H2_OP) + def test_mismatching_weights(self, op): + """Ensuring incorrect number of weights is caught""" + optimizer = SLSQP(maxiter=50) + ssvqe = MCVQE( + estimator=self.estimator, k=2, weight_vector=[3, 2, 1], optimizer=optimizer + ) + with self.assertRaises(AlgorithmError): + _ = ssvqe.compute_eigenvalues(operator=op) + + @data(H2_PAULI, H2_OP) + def test_initial_states_nonorthogonal(self, op): + """Ensuring non-orthogonality of input initial states is caught""" + initial_states = [QuantumCircuit(2), QuantumCircuit(2)] + initial_states[0].x(0) + initial_states[1].x(0) + optimizer = SLSQP(maxiter=50) + ssvqe = MCVQE( + estimator=self.estimator, + k=1, + optimizer=optimizer, + initial_states=initial_states, + ) + with self.assertRaises(AlgorithmError): + _ = ssvqe.compute_eigenvalues(operator=op) + + @data(H2_PAULI, H2_OP) + def test_missing_varform_params(self, op): + """Test specifying a variational form with no parameters raises an error.""" + circuit = QuantumCircuit(op.num_qubits) + ssvqe = MCVQE(estimator=self.estimator, ansatz=circuit, optimizer=SLSQP(), k=1) + with self.assertRaises(AlgorithmError): + ssvqe.compute_eigenvalues(operator=op) + + @data(H2_PAULI, H2_OP) + def test_callback(self, op): + """Test the callback on SSVQE.""" + history = { + "eval_count": [], + "parameters": [], + "mean_energies": [], + "metadata": [], + } + + def store_intermediate_result(eval_count, parameters, mean, metadata): + history["eval_count"].append(eval_count) + history["parameters"].append(parameters) + history["mean_energies"].append(mean) + history["metadata"].append(metadata) + + optimizer = COBYLA(maxiter=3) + wavefunction = self.ry_wavefunction + + ssvqe = MCVQE( + estimator=self.estimator, + ansatz=wavefunction, + optimizer=optimizer, + callback=store_intermediate_result, + ) + + ssvqe.compute_eigenvalues(operator=op) + + self.assertTrue(all(isinstance(count, int) for count in history["eval_count"])) + # self.assertTrue(all(isinstance(mean, float) for mean in history["mean_energies"])) + self.assertTrue( + all( + isinstance(mean, float) + for mean_list in history["mean_energies"] + for mean in mean_list + ) + ) + self.assertTrue( + all(isinstance(metadata, dict) for metadata in history["metadata"]) + ) + for params in history["parameters"]: + self.assertTrue(all(isinstance(param, float) for param in params)) + + ref_eval_count = [1, 2, 3] + ref_mean = [[-1.07, -1.44], [-1.45, -1.06], [-1.37, -0.94]] + + np.testing.assert_array_equal(history["eval_count"], ref_eval_count) + np.testing.assert_array_almost_equal( + history["mean_energies"], ref_mean, decimal=2 + ) + + @data(H2_PAULI, H2_OP) + def test_ssvqe_optimizer(self, op): + """Test running same SSVQE twice to re-use optimizer, then switch optimizer""" + ssvqe = MCVQE( + estimator=self.estimator, + ansatz=RealAmplitudes(reps=6), + optimizer=SLSQP(), + k=2, + ) + + def run_check(): + result = ssvqe.compute_eigenvalues(operator=op) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=3 + ) + + run_check() + + with self.subTest("Optimizer re-use"): + run_check() + + with self.subTest("Optimizer replace"): + ssvqe.optimizer = L_BFGS_B() + run_check() + + @data(H2_PAULI, H2_OP) + def test_aux_operators_list(self, op): + """Test list-based aux_operators.""" + wavefunction = self.ry_wavefunction + ssvqe = MCVQE( + estimator=self.estimator, ansatz=wavefunction, optimizer=SLSQP(), k=2 + ) + + # Start with an empty list + result = ssvqe.compute_eigenvalues(op, aux_operators=[]) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=2 + ) + self.assertIsNone(result.aux_operators_evaluated) + + # Go again with two auxiliary operators + aux_op1 = SparsePauliOp.from_list([("II", 2.0)]) + aux_op2 = SparsePauliOp.from_list( + [("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)] + ) + aux_ops = [aux_op1, aux_op2] + result = ssvqe.compute_eigenvalues(op, aux_operators=aux_ops) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=2 + ) + self.assertEqual(len(result.aux_operators_evaluated), 2) + # expectation values + self.assertAlmostEqual(result.aux_operators_evaluated[0][0][0], 2, places=2) + self.assertAlmostEqual(result.aux_operators_evaluated[0][1][0], 0, places=2) + # metadata + self.assertIsInstance(result.aux_operators_evaluated[0][1][1], dict) + self.assertIsInstance(result.aux_operators_evaluated[0][1][1], dict) + + # Go again with additional None and zero operators + extra_ops = [*aux_ops, None, 0] + result = ssvqe.compute_eigenvalues(op, aux_operators=extra_ops) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=2 + ) + self.assertEqual(len(result.aux_operators_evaluated), 2) + # expectation values + self.assertAlmostEqual(result.aux_operators_evaluated[0][0][0], 2, places=2) + self.assertAlmostEqual(result.aux_operators_evaluated[0][1][0], 0, places=2) + self.assertEqual(result.aux_operators_evaluated[0][2][0], 0.0) + self.assertEqual(result.aux_operators_evaluated[0][3][0], 0.0) + # metadata + self.assertIsInstance(result.aux_operators_evaluated[0][0][1], dict) + self.assertIsInstance(result.aux_operators_evaluated[0][1][1], dict) + self.assertIsInstance(result.aux_operators_evaluated[0][3][1], dict) + + @data(H2_PAULI, H2_OP) + def test_aux_operators_dict(self, op): + """Test dictionary compatibility of aux_operators""" + wavefunction = self.ry_wavefunction + ssvqe = MCVQE(estimator=self.estimator, ansatz=wavefunction, optimizer=SLSQP()) + + # Start with an empty dictionary + result = ssvqe.compute_eigenvalues(op, aux_operators={}) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=2 + ) + self.assertIsNone(result.aux_operators_evaluated) + + # Go again with two auxiliary operators + aux_op1 = SparsePauliOp.from_list([("II", 2.0)]) + aux_op2 = SparsePauliOp.from_list( + [("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)] + ) + aux_ops = {"aux_op1": aux_op1, "aux_op2": aux_op2} + result = ssvqe.compute_eigenvalues(op, aux_operators=aux_ops) + self.assertEqual(len(result.eigenvalues), 2) + self.assertEqual(result.eigenvalues.dtype, np.float64) + self.assertAlmostEqual(result.eigenvalues[0], -1.85727503, 2) + self.assertEqual(len(result.aux_operators_evaluated), 2) + self.assertEqual(len(result.aux_operators_evaluated[0]), 2) + # expectation values + self.assertAlmostEqual( + result.aux_operators_evaluated[0]["aux_op1"][0], 2, places=6 + ) + self.assertAlmostEqual( + result.aux_operators_evaluated[0]["aux_op2"][0], 0, places=1 + ) + # metadata + self.assertIsInstance(result.aux_operators_evaluated[0]["aux_op1"][1], dict) + self.assertIsInstance(result.aux_operators_evaluated[0]["aux_op2"][1], dict) + + # Go again with additional None and zero operators + extra_ops = {**aux_ops, "None_operator": None, "zero_operator": 0} + result = ssvqe.compute_eigenvalues(op, aux_operators=extra_ops) + self.assertEqual(len(result.eigenvalues), 2) + self.assertEqual(result.eigenvalues.dtype, np.float64) + self.assertAlmostEqual(result.eigenvalues[0], -1.85727503, places=5) + self.assertEqual(len(result.aux_operators_evaluated), 2) + self.assertEqual(len(result.aux_operators_evaluated[0]), 3) + # expectation values + self.assertAlmostEqual( + result.aux_operators_evaluated[0]["aux_op1"][0], 2, places=6 + ) + self.assertAlmostEqual( + result.aux_operators_evaluated[0]["aux_op2"][0], 0.0, places=2 + ) + self.assertEqual(result.aux_operators_evaluated[0]["zero_operator"][0], 0.0) + self.assertTrue("None_operator" not in result.aux_operators_evaluated[0].keys()) + # metadata + self.assertIsInstance(result.aux_operators_evaluated[0]["aux_op1"][1], dict) + self.assertIsInstance(result.aux_operators_evaluated[0]["aux_op2"][1], dict) + self.assertIsInstance( + result.aux_operators_evaluated[0]["zero_operator"][1], dict + ) + + @data(H2_PAULI, H2_OP) + def test_aux_operator_std_dev(self, op): + """Test non-zero standard deviations of aux operators.""" + wavefunction = self.ry_wavefunction + ssvqe = MCVQE( + estimator=self.estimator, + ansatz=wavefunction, + initial_point=[ + 1.70256666, + -5.34843975, + -0.39542903, + 5.99477786, + -2.74374986, + -4.85284669, + 0.2442925, + -1.51638917, + ], + optimizer=COBYLA(maxiter=0), + ) + + # Go again with two auxiliary operators + aux_op1 = SparsePauliOp.from_list([("II", 2.0)]) + aux_op2 = SparsePauliOp.from_list( + [("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)] + ) + aux_ops = [aux_op1, aux_op2] + result = ssvqe.compute_eigenvalues(op, aux_operators=aux_ops) + self.assertEqual(len(result.aux_operators_evaluated), 2) + # expectation values + self.assertAlmostEqual(result.aux_operators_evaluated[0][0][0], 2.0, places=1) + self.assertAlmostEqual( + result.aux_operators_evaluated[0][1][0], 0.0019531249999999445, places=1 + ) + # metadata + self.assertIsInstance(result.aux_operators_evaluated[0][0][1], dict) + self.assertIsInstance(result.aux_operators_evaluated[0][1][1], dict) + + # Go again with additional None and zero operators + aux_ops = [*aux_ops, None, 0] + result = ssvqe.compute_eigenvalues(op, aux_operators=aux_ops) + self.assertEqual(len(result.aux_operators_evaluated[0]), 4) + # expectation values + self.assertAlmostEqual(result.aux_operators_evaluated[0][0][0], 2.0, places=1) + self.assertAlmostEqual( + result.aux_operators_evaluated[0][1][0], 0.0019531249999999445, places=1 + ) + self.assertEqual(result.aux_operators_evaluated[0][2][0], 0.0) + self.assertEqual(result.aux_operators_evaluated[0][3][0], 0.0) + # metadata + self.assertIsInstance(result.aux_operators_evaluated[0][0][1], dict) + self.assertIsInstance(result.aux_operators_evaluated[0][1][1], dict) + self.assertIsInstance(result.aux_operators_evaluated[0][2][1], dict) + self.assertIsInstance(result.aux_operators_evaluated[0][3][1], dict) + + @data( + CG(), + L_BFGS_B(), + P_BFGS(), + SLSQP(), + TNC(), + ) + def test_with_gradient(self, optimizer): + """Test SSVQE using gradient primitive.""" + estimator = Estimator() + ssvqe = MCVQE( + k=2, + estimator=estimator, + ansatz=self.ry_wavefunction, + optimizer=optimizer, + gradient=ParamShiftEstimatorGradient(estimator), + ) + result = ssvqe.compute_eigenvalues(operator=H2_PAULI) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=5 + ) + + def test_gradient_passed(self): + """Test the gradient is properly passed into the optimizer.""" + inputs = {} + estimator = Estimator() + ssvqe = MCVQE( + k=2, + estimator=estimator, + ansatz=RealAmplitudes(), + optimizer=partial(_mock_optimizer, inputs=inputs), + gradient=ParamShiftEstimatorGradient(estimator), + ) + _ = ssvqe.compute_eigenvalues(operator=H2_PAULI) + + self.assertIsNotNone(inputs["jac"]) + + def test_gradient_run(self): + """Test using the gradient to calculate the minimum.""" + estimator = Estimator() + ssvqe = MCVQE( + k=2, + estimator=estimator, + ansatz=RealAmplitudes(), + optimizer=GradientDescent(maxiter=200, learning_rate=0.1), + gradient=ParamShiftEstimatorGradient(estimator), + ) + result = ssvqe.compute_eigenvalues(operator=H2_PAULI) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=5 + ) + + def test_max_evals_grouped(self): + """Test with SLSQP with max_evals_grouped.""" + optimizer = SLSQP(maxiter=50, max_evals_grouped=5) + ssvqe = MCVQE( + k=2, + estimator=Estimator(), + ansatz=RealAmplitudes(reps=6), + optimizer=optimizer, + ) + result = ssvqe.compute_eigenvalues(operator=H2_PAULI) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=5 + ) + + +if __name__ == "__main__": + unittest.main() \ No newline at end of file diff --git a/tests/test_optorbvqe.py b/tests/test_optorbvqe.py new file mode 100644 index 0000000..e2af9c9 --- /dev/null +++ b/tests/test_optorbvqe.py @@ -0,0 +1,140 @@ +"""Test OptOrbVQE""" + +import unittest +from .algorithms_test_case import QiskitAlgorithmsTestCase + +from functools import partial +import numpy as np +import torch +from ddt import data, ddt + +from qiskit import QuantumCircuit +from qiskit_algorithms import AlgorithmError +from qiskit_algorithms.optimizers import L_BFGS_B + +import numpy as np +from qiskit_nature.second_q.drivers import PySCFDriver +from qiskit_nature.units import DistanceUnit +from qiskit_nature.second_q.mappers import JordanWignerMapper +from qiskit_algorithms.minimum_eigensolvers import VQE +from qiskit_aer.primitives import Estimator +from qiskit_nature.second_q.circuit.library import HartreeFock, UCCSD +from qiskit_nature.second_q.operators.tensor_ordering import to_physicist_ordering +from qiskit_algorithms.utils import algorithm_globals + +from electronic_structure_algorithms.orbital_optimization import PartialUnitaryProjectionOptimizer, OptOrbVQE + +estimator = Estimator(approximation=True) +mapper=JordanWignerMapper() + +driver = PySCFDriver(atom=f'H 0 0 0; H 0 0 {0.735}', + charge=0, + spin=0, + unit=DistanceUnit.ANGSTROM, + basis='6-31G') + +es_problem = driver.run() +num_particles = es_problem.num_particles +num_reduced_qubits = 4 + +one_body_integrals = torch.from_numpy(np.asarray(es_problem.hamiltonian.electronic_integrals.second_q_coeffs()["+-"].to_dense())) +two_body_integrals = torch.from_numpy(np.asarray(-1*to_physicist_ordering(es_problem.hamiltonian.electronic_integrals.second_q_coeffs()["++--"].to_dense()))) + +integrals = one_body_integrals, two_body_integrals + +HF_state = HartreeFock(qubit_mapper=mapper, + num_spatial_orbitals=int(num_reduced_qubits/2), + num_particles=num_particles) + +uccsd = UCCSD(qubit_mapper=mapper, + num_spatial_orbitals=int(num_reduced_qubits/2), + num_particles=num_particles, + initial_state=HF_state) + +@ddt +class TestOptOrbVQE(QiskitAlgorithmsTestCase): + """Test OptOrbVQE""" + + def setUp(self): + super().setUp() + self.seed = 50 + #algorithm_globals.random_seed = self.seed + + self.HF_state = HartreeFock(qubit_mapper=mapper, + num_spatial_orbitals=int(num_reduced_qubits/2), + num_particles=num_particles) + + self.uccsd = UCCSD(qubit_mapper=mapper, + num_spatial_orbitals=int(num_reduced_qubits/2), + num_particles=num_particles, + initial_state=HF_state) + + self.partial_unitary_optimizer = PartialUnitaryProjectionOptimizer(initial_BBstepsize=10**-3, + stopping_tolerance=10**-5, + maxiter=10000, + gradient_method='autograd') + + self.h2_energy = -1.8661038079694765 + self.estimator = Estimator(approximation=True) + + @data(es_problem) + def test_ground_state_problem(self, problem): + + vqe_instance = VQE(ansatz=self.uccsd, + initial_point=np.zeros(self.uccsd.num_parameters), + optimizer=L_BFGS_B(maxfun=10**6, maxiter=10**6), + estimator=estimator) + + optorbvqe_instance = OptOrbVQE(problem=problem, + integral_tensors = None, + num_spin_orbitals=num_reduced_qubits, + ground_state_solver=vqe_instance, + mapper=mapper, + estimator=estimator, + partial_unitary_optimizer=self.partial_unitary_optimizer, + maxiter=20, + wavefuntion_real=True, + spin_conserving=True, + stopping_tolerance=10**-5, + partial_unitary_random_perturbation=0.01, + minimum_eigensolver_random_perturbation=0.0) + + result = optorbvqe_instance.compute_minimum_energy() + + with self.subTest(msg="test eigenvalue, provided ElectronicStructureProblem"): + np.testing.almost_equal( + [result.eigenvalue.real], [self.h2_energy], decimal=3 + ) + + @data(integrals) + def test_ground_state_problem(self, integrals): + + vqe_instance = VQE(ansatz=self.uccsd, + initial_point=np.zeros(self.uccsd.num_parameters), + optimizer=L_BFGS_B(maxfun=10**6, maxiter=10**6), + estimator=estimator) + + optorbvqe_instance = OptOrbVQE(problem=None, + integral_tensors=integrals, + num_spin_orbitals=num_reduced_qubits, + ground_state_solver=vqe_instance, + mapper=mapper, + estimator=estimator, + partial_unitary_optimizer=self.partial_unitary_optimizer, + maxiter=20, + wavefuntion_real=True, + spin_conserving=True, + stopping_tolerance=10**-5, + partial_unitary_random_perturbation=0.01, + minimum_eigensolver_random_perturbation=0.0) + + result = optorbvqe_instance.compute_minimum_energy() + + with self.subTest(msg="test eigenvalue, provided integrals"): + np.testing.almost_equal( + [result.eigenvalue.real], [self.h2_energy], decimal=3 + ) + + +if __name__ == "__main__": + unittest.main() \ No newline at end of file diff --git a/tests/test_random.py b/tests/test_random.py deleted file mode 100644 index 275d488..0000000 --- a/tests/test_random.py +++ /dev/null @@ -1,14 +0,0 @@ -"""Docstring.""" -import unittest - -from src import Random - - -class TestRandom(unittest.TestCase): - """Tests Random class implementation.""" - - def test_run(self): - """Tests run method random.""" - random = Random() - - self.assertEqual(random.run(2), 4) diff --git a/tests/test_ssvqe.py b/tests/test_ssvqe.py new file mode 100644 index 0000000..ffc2bb7 --- /dev/null +++ b/tests/test_ssvqe.py @@ -0,0 +1,492 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2022 +# +# 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. + +# These tests are derived from Qiskit's VQD tests. + +""" Test SSVQE """ + +import unittest +from .algorithms_test_case import QiskitAlgorithmsTestCase + +from functools import partial +import numpy as np +from ddt import data, ddt + +from qiskit import QuantumCircuit +from qiskit_algorithms import AlgorithmError +from qiskit_algorithms.optimizers import ( + COBYLA, + L_BFGS_B, + SLSQP, + CG, + TNC, + P_BFGS, + GradientDescent, + OptimizerResult, +) + +from qiskit_algorithms.gradients import ParamShiftEstimatorGradient +from qiskit.circuit.library import TwoLocal, RealAmplitudes +from qiskit.quantum_info import SparsePauliOp +from qiskit.primitives import Sampler, Estimator +from qiskit_algorithms.state_fidelities import ComputeUncompute +from qiskit_algorithms.utils import algorithm_globals +from qiskit.quantum_info.operators import Operator + +from electronic_structure_algorithms.excited_states_eigensolvers import SSVQE + +# pylint: disable=invalid-name +def _mock_optimizer(fun, x0, jac=None, bounds=None, inputs=None) -> OptimizerResult: + """A mock of a callable that can be used as minimizer in SSVQE.""" + result = OptimizerResult() + result.x = np.zeros_like(x0) + result.fun = fun(result.x) + result.nit = 0 + + if inputs is not None: + inputs.update({"fun": fun, "x0": x0, "jac": jac, "bounds": bounds}) + return result + + +I = SparsePauliOp.from_list([("I", 1)]) # pylint: disable=invalid-name +X = SparsePauliOp.from_list([("X", 1)]) # pylint: disable=invalid-name +Z = SparsePauliOp.from_list([("Z", 1)]) # pylint: disable=invalid-name + +H2_PAULI = ( + -1.052373245772859 * (I ^ I) + + 0.39793742484318045 * (I ^ Z) + - 0.39793742484318045 * (Z ^ I) + - 0.01128010425623538 * (Z ^ Z) + + 0.18093119978423156 * (X ^ X) +) + +H2_OP = Operator(H2_PAULI.to_matrix()) + + +@ddt +class TestSSVQE(QiskitAlgorithmsTestCase): + """Test SSVQE""" + + def setUp(self): + super().setUp() + self.seed = 50 + algorithm_globals.random_seed = self.seed + + self.h2_energy = -1.85727503 + self.h2_energy_excited = [-1.85727503, -1.24458455] + + self.ryrz_wavefunction = TwoLocal( + rotation_blocks=["ry", "rz"], entanglement_blocks="cz", reps=1 + ) + self.ry_wavefunction = TwoLocal(rotation_blocks="ry", entanglement_blocks="cz") + + self.estimator = Estimator() + self.estimator_shots = Estimator(options={"shots": 2048, "seed": self.seed}) + self.fidelity = ComputeUncompute(Sampler()) + self.betas = [50, 50] + + @data(H2_PAULI, H2_OP) + def test_basic_operator(self, op): + """Test SSVQE without aux_operators.""" + wavefunction = self.ryrz_wavefunction + ssvqe = SSVQE( + estimator=self.estimator, + ansatz=wavefunction, + optimizer=COBYLA(), + ) + + result = ssvqe.compute_eigenvalues(operator=op) + + with self.subTest(msg="test eigenvalue"): + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=1 + ) + + with self.subTest(msg="test dimension of optimal point"): + self.assertEqual(len(result.optimal_point), 8) + + with self.subTest(msg="assert cost_function_evals is set"): + self.assertIsNotNone(result.cost_function_evals) + + with self.subTest(msg="assert optimizer_time is set"): + self.assertIsNotNone(result.optimizer_time) + + @data(H2_PAULI, H2_OP) + def test_mismatching_num_qubits(self, op): + """Ensuring circuit and operator mismatch is caught""" + wavefunction = QuantumCircuit(1) + optimizer = SLSQP(maxiter=50) + ssvqe = SSVQE( + estimator=self.estimator, k=1, ansatz=wavefunction, optimizer=optimizer + ) + with self.assertRaises(AlgorithmError): + _ = ssvqe.compute_eigenvalues(operator=op) + + @data(H2_PAULI, H2_OP) + def test_initial_states_mismatching_num_qubits(self, op): + """Ensuring initial states and operator mismatch is caught""" + initial_state = [QuantumCircuit(1)] + wavefunction = QuantumCircuit(2) + optimizer = SLSQP(maxiter=50) + ssvqe = SSVQE( + estimator=self.estimator, + k=1, + ansatz=wavefunction, + optimizer=optimizer, + initial_states=initial_state, + ) + with self.assertRaises(AlgorithmError): + _ = ssvqe.compute_eigenvalues(operator=op) + + @data(H2_PAULI, H2_OP) + def test_mismatching_weights(self, op): + """Ensuring incorrect number of weights is caught""" + optimizer = SLSQP(maxiter=50) + ssvqe = SSVQE( + estimator=self.estimator, k=2, weight_vector=[3, 2, 1], optimizer=optimizer + ) + with self.assertRaises(AlgorithmError): + _ = ssvqe.compute_eigenvalues(operator=op) + + @data(H2_PAULI, H2_OP) + def test_initial_states_nonorthogonal(self, op): + """Ensuring non-orthogonality of input initial states is caught""" + initial_states = [QuantumCircuit(2), QuantumCircuit(2)] + initial_states[0].x(0) + initial_states[1].x(0) + optimizer = SLSQP(maxiter=50) + ssvqe = SSVQE( + estimator=self.estimator, + k=1, + optimizer=optimizer, + initial_states=initial_states, + ) + with self.assertRaises(AlgorithmError): + _ = ssvqe.compute_eigenvalues(operator=op) + + @data(H2_PAULI, H2_OP) + def test_missing_varform_params(self, op): + """Test specifying a variational form with no parameters raises an error.""" + circuit = QuantumCircuit(op.num_qubits) + ssvqe = SSVQE(estimator=self.estimator, ansatz=circuit, optimizer=SLSQP(), k=1) + with self.assertRaises(AlgorithmError): + ssvqe.compute_eigenvalues(operator=op) + + @data(H2_PAULI, H2_OP) + def test_callback(self, op): + """Test the callback on SSVQE.""" + history = { + "eval_count": [], + "parameters": [], + "mean_energies": [], + "metadata": [], + } + + def store_intermediate_result(eval_count, parameters, mean, metadata): + history["eval_count"].append(eval_count) + history["parameters"].append(parameters) + history["mean_energies"].append(mean) + history["metadata"].append(metadata) + + optimizer = COBYLA(maxiter=3) + wavefunction = self.ry_wavefunction + + ssvqe = SSVQE( + estimator=self.estimator, + ansatz=wavefunction, + optimizer=optimizer, + callback=store_intermediate_result, + ) + + ssvqe.compute_eigenvalues(operator=op) + + self.assertTrue(all(isinstance(count, int) for count in history["eval_count"])) + # self.assertTrue(all(isinstance(mean, float) for mean in history["mean_energies"])) + self.assertTrue( + all( + isinstance(mean, float) + for mean_list in history["mean_energies"] + for mean in mean_list + ) + ) + self.assertTrue( + all(isinstance(metadata, dict) for metadata in history["metadata"]) + ) + for params in history["parameters"]: + self.assertTrue(all(isinstance(param, float) for param in params)) + + ref_eval_count = [1, 2, 3] + ref_mean = [[-1.07, -1.44], [-1.45, -1.06], [-1.37, -0.94]] + + np.testing.assert_array_equal(history["eval_count"], ref_eval_count) + np.testing.assert_array_almost_equal( + history["mean_energies"], ref_mean, decimal=2 + ) + + @data(H2_PAULI, H2_OP) + def test_ssvqe_optimizer(self, op): + """Test running same SSVQE twice to re-use optimizer, then switch optimizer""" + ssvqe = SSVQE( + estimator=self.estimator, + ansatz=RealAmplitudes(reps=6), + optimizer=SLSQP(), + k=2, + ) + + def run_check(): + result = ssvqe.compute_eigenvalues(operator=op) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=3 + ) + + run_check() + + with self.subTest("Optimizer re-use"): + run_check() + + with self.subTest("Optimizer replace"): + ssvqe.optimizer = L_BFGS_B() + run_check() + + @data(H2_PAULI, H2_OP) + def test_aux_operators_list(self, op): + """Test list-based aux_operators.""" + wavefunction = self.ry_wavefunction + ssvqe = SSVQE( + estimator=self.estimator, ansatz=wavefunction, optimizer=SLSQP(), k=2 + ) + + # Start with an empty list + result = ssvqe.compute_eigenvalues(op, aux_operators=[]) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=2 + ) + self.assertIsNone(result.aux_operators_evaluated) + + # Go again with two auxiliary operators + aux_op1 = SparsePauliOp.from_list([("II", 2.0)]) + aux_op2 = SparsePauliOp.from_list( + [("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)] + ) + aux_ops = [aux_op1, aux_op2] + result = ssvqe.compute_eigenvalues(op, aux_operators=aux_ops) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=2 + ) + self.assertEqual(len(result.aux_operators_evaluated), 2) + # expectation values + self.assertAlmostEqual(result.aux_operators_evaluated[0][0][0], 2, places=2) + self.assertAlmostEqual(result.aux_operators_evaluated[0][1][0], 0, places=2) + # metadata + self.assertIsInstance(result.aux_operators_evaluated[0][1][1], dict) + self.assertIsInstance(result.aux_operators_evaluated[0][1][1], dict) + + # Go again with additional None and zero operators + extra_ops = [*aux_ops, None, 0] + result = ssvqe.compute_eigenvalues(op, aux_operators=extra_ops) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=2 + ) + self.assertEqual(len(result.aux_operators_evaluated), 2) + # expectation values + self.assertAlmostEqual(result.aux_operators_evaluated[0][0][0], 2, places=2) + self.assertAlmostEqual(result.aux_operators_evaluated[0][1][0], 0, places=2) + self.assertEqual(result.aux_operators_evaluated[0][2][0], 0.0) + self.assertEqual(result.aux_operators_evaluated[0][3][0], 0.0) + # metadata + self.assertIsInstance(result.aux_operators_evaluated[0][0][1], dict) + self.assertIsInstance(result.aux_operators_evaluated[0][1][1], dict) + self.assertIsInstance(result.aux_operators_evaluated[0][3][1], dict) + + @data(H2_PAULI, H2_OP) + def test_aux_operators_dict(self, op): + """Test dictionary compatibility of aux_operators""" + wavefunction = self.ry_wavefunction + ssvqe = SSVQE(estimator=self.estimator, ansatz=wavefunction, optimizer=SLSQP()) + + # Start with an empty dictionary + result = ssvqe.compute_eigenvalues(op, aux_operators={}) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=2 + ) + self.assertIsNone(result.aux_operators_evaluated) + + # Go again with two auxiliary operators + aux_op1 = SparsePauliOp.from_list([("II", 2.0)]) + aux_op2 = SparsePauliOp.from_list( + [("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)] + ) + aux_ops = {"aux_op1": aux_op1, "aux_op2": aux_op2} + result = ssvqe.compute_eigenvalues(op, aux_operators=aux_ops) + self.assertEqual(len(result.eigenvalues), 2) + self.assertEqual(result.eigenvalues.dtype, np.float64) + self.assertAlmostEqual(result.eigenvalues[0], -1.85727503, 2) + self.assertEqual(len(result.aux_operators_evaluated), 2) + self.assertEqual(len(result.aux_operators_evaluated[0]), 2) + # expectation values + self.assertAlmostEqual( + result.aux_operators_evaluated[0]["aux_op1"][0], 2, places=6 + ) + self.assertAlmostEqual( + result.aux_operators_evaluated[0]["aux_op2"][0], 0, places=1 + ) + # metadata + self.assertIsInstance(result.aux_operators_evaluated[0]["aux_op1"][1], dict) + self.assertIsInstance(result.aux_operators_evaluated[0]["aux_op2"][1], dict) + + # Go again with additional None and zero operators + extra_ops = {**aux_ops, "None_operator": None, "zero_operator": 0} + result = ssvqe.compute_eigenvalues(op, aux_operators=extra_ops) + self.assertEqual(len(result.eigenvalues), 2) + self.assertEqual(result.eigenvalues.dtype, np.float64) + self.assertAlmostEqual(result.eigenvalues[0], -1.85727503, places=5) + self.assertEqual(len(result.aux_operators_evaluated), 2) + self.assertEqual(len(result.aux_operators_evaluated[0]), 3) + # expectation values + self.assertAlmostEqual( + result.aux_operators_evaluated[0]["aux_op1"][0], 2, places=6 + ) + self.assertAlmostEqual( + result.aux_operators_evaluated[0]["aux_op2"][0], 0.0, places=2 + ) + self.assertEqual(result.aux_operators_evaluated[0]["zero_operator"][0], 0.0) + self.assertTrue("None_operator" not in result.aux_operators_evaluated[0].keys()) + # metadata + self.assertIsInstance(result.aux_operators_evaluated[0]["aux_op1"][1], dict) + self.assertIsInstance(result.aux_operators_evaluated[0]["aux_op2"][1], dict) + self.assertIsInstance( + result.aux_operators_evaluated[0]["zero_operator"][1], dict + ) + + @data(H2_PAULI, H2_OP) + def test_aux_operator_std_dev(self, op): + """Test non-zero standard deviations of aux operators.""" + wavefunction = self.ry_wavefunction + ssvqe = SSVQE( + estimator=self.estimator, + ansatz=wavefunction, + initial_point=[ + 1.70256666, + -5.34843975, + -0.39542903, + 5.99477786, + -2.74374986, + -4.85284669, + 0.2442925, + -1.51638917, + ], + optimizer=COBYLA(maxiter=0), + ) + + # Go again with two auxiliary operators + aux_op1 = SparsePauliOp.from_list([("II", 2.0)]) + aux_op2 = SparsePauliOp.from_list( + [("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)] + ) + aux_ops = [aux_op1, aux_op2] + result = ssvqe.compute_eigenvalues(op, aux_operators=aux_ops) + self.assertEqual(len(result.aux_operators_evaluated), 2) + # expectation values + self.assertAlmostEqual(result.aux_operators_evaluated[0][0][0], 2.0, places=1) + self.assertAlmostEqual( + result.aux_operators_evaluated[0][1][0], 0.0019531249999999445, places=1 + ) + # metadata + self.assertIsInstance(result.aux_operators_evaluated[0][0][1], dict) + self.assertIsInstance(result.aux_operators_evaluated[0][1][1], dict) + + # Go again with additional None and zero operators + aux_ops = [*aux_ops, None, 0] + result = ssvqe.compute_eigenvalues(op, aux_operators=aux_ops) + self.assertEqual(len(result.aux_operators_evaluated[0]), 4) + # expectation values + self.assertAlmostEqual(result.aux_operators_evaluated[0][0][0], 2.0, places=1) + self.assertAlmostEqual( + result.aux_operators_evaluated[0][1][0], 0.0019531249999999445, places=1 + ) + self.assertEqual(result.aux_operators_evaluated[0][2][0], 0.0) + self.assertEqual(result.aux_operators_evaluated[0][3][0], 0.0) + # metadata + self.assertIsInstance(result.aux_operators_evaluated[0][0][1], dict) + self.assertIsInstance(result.aux_operators_evaluated[0][1][1], dict) + self.assertIsInstance(result.aux_operators_evaluated[0][2][1], dict) + self.assertIsInstance(result.aux_operators_evaluated[0][3][1], dict) + + @data( + CG(), + L_BFGS_B(), + P_BFGS(), + SLSQP(), + TNC(), + ) + def test_with_gradient(self, optimizer): + """Test SSVQE using gradient primitive.""" + estimator = Estimator() + ssvqe = SSVQE( + k=2, + estimator=estimator, + ansatz=self.ry_wavefunction, + optimizer=optimizer, + gradient=ParamShiftEstimatorGradient(estimator), + ) + result = ssvqe.compute_eigenvalues(operator=H2_PAULI) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=5 + ) + + def test_gradient_passed(self): + """Test the gradient is properly passed into the optimizer.""" + inputs = {} + estimator = Estimator() + ssvqe = SSVQE( + k=2, + estimator=estimator, + ansatz=RealAmplitudes(), + optimizer=partial(_mock_optimizer, inputs=inputs), + gradient=ParamShiftEstimatorGradient(estimator), + ) + _ = ssvqe.compute_eigenvalues(operator=H2_PAULI) + + self.assertIsNotNone(inputs["jac"]) + + def test_gradient_run(self): + """Test using the gradient to calculate the minimum.""" + estimator = Estimator() + ssvqe = SSVQE( + k=2, + estimator=estimator, + ansatz=RealAmplitudes(), + optimizer=GradientDescent(maxiter=200, learning_rate=0.1), + gradient=ParamShiftEstimatorGradient(estimator), + ) + result = ssvqe.compute_eigenvalues(operator=H2_PAULI) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=5 + ) + + def test_max_evals_grouped(self): + """Test with SLSQP with max_evals_grouped.""" + optimizer = SLSQP(maxiter=50, max_evals_grouped=5) + ssvqe = SSVQE( + k=2, + estimator=Estimator(), + ansatz=RealAmplitudes(reps=6), + optimizer=optimizer, + ) + result = ssvqe.compute_eigenvalues(operator=H2_PAULI) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=5 + ) + + +if __name__ == "__main__": + unittest.main() \ No newline at end of file diff --git a/tox.ini b/tox.ini index 944dbcc..4cb50da 100644 --- a/tox.ini +++ b/tox.ini @@ -30,7 +30,7 @@ commands = [testenv:black] envdir = .tox/lint -commands = black {posargs} src tests --check +commands = black {posargs} electronic_structure_algorithms tests --check [testenv:json] allowlist_externals = /bin/bash