Source code for qiskit_nature.second_q.algorithms.initial_points.mp2_initial_point

# This code is part of a Qiskit project.
# (C) Copyright IBM 2019, 2023.
# 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
# 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.

"""Second-order Møller-Plesset perturbation theory (MP2) computation."""

from __future__ import annotations

import numpy as np

from qiskit_nature.exceptions import QiskitNatureError
from qiskit_nature.second_q.circuit.library import UCC
from qiskit_nature.second_q.operators import ElectronicIntegrals
from qiskit_nature.second_q.operators.tensor import Tensor
from qiskit_nature.second_q.operators.symmetric_two_body import SymmetricTwoBodyIntegrals, unfold
from qiskit_nature.second_q.problems import BaseProblem, ElectronicStructureProblem
from qiskit_nature.utils import get_einsum

from .initial_point import InitialPoint

def _compute_mp2(
    num_occ: int, integral_matrix: Tensor, orbital_energies: np.ndarray
) -> tuple[np.ndarray, float]:
    """Compute the T2 amplitudes and MP2 energy correction.

        num_occ: The number of occupied molecular orbitals.
        integral_matrix: The two-body molecular orbitals tensor.
        orbital_energies: The orbital energies.

        A tuple consisting of the:
        - T amplitudes t2[i, j, a, b] (i, j in occupied, a, b in virtual).
        - The MP2 energy correction.

    if isinstance(integral_matrix, SymmetricTwoBodyIntegrals):
        integral_matrix = unfold(integral_matrix)

    # We use NumPy broadcasting to compute the matrix of occupied - virtual energy deltas with
    # shape (num_occ, num_vir), such that
    # energy_deltas[i, a] = orbital_energy[i] - orbital_energy[a].
    # NOTE: In the unrestricted-spin calculation, the orbital energies will be a 2D array, and this
    # logic will need to be revisited.
    energy_deltas = orbital_energies[:num_occ, np.newaxis] - orbital_energies[num_occ:]

    # We now want to compute a 4D tensor of (occupied, occupied) - (virtual, virtual)
    # energy deltas with shape (num_occ, num_vir, num_occ, num_vir), such that
    # double_deltas[i, a, j, b] = orbital_energies[i] + orbital_energies[j]
    #                             - orbital_energies[a] - orbital_energies[b].
    # Again we can use NumPy broadcasting to speed this up.
    double_deltas = energy_deltas[:, :, np.newaxis, np.newaxis] + energy_deltas

    # Now we transpose this matrix into (num_occ, num_occ, num_vir, num_vir) which is the expected
    # ordering of T2 amplitudes following the convention set out by PySCF
    double_deltas = double_deltas.transpose(0, 2, 1, 3)

    # We must swap the last two axes in order to match the ordering of double_deltas as
    # explained above. We use the _reverse_label_template routine here to handle non-default label
    # templates in the integral_matrix Tensor, too.
    integral_matrix = integral_matrix.transpose(
        integral_matrix._reverse_label_template((0, 1, 3, 2))
    # Create integral matrix that uses occupied and virtual indices rather than MO indices.
    integral_matrix_oovv = integral_matrix[:num_occ, :num_occ, num_occ:, num_occ:]

    # Compute T2 amplitudes and transpose to num_occ, num_occ, num_vir, num_vir.
    t2_amplitudes = integral_matrix_oovv / double_deltas

    # Compute MP2 energy correction.
    einsum_func, _ = get_einsum()
    energy_correction = einsum_func("ijab,ijab", t2_amplitudes, integral_matrix_oovv) * 2
    energy_correction -= einsum_func("ijab,ijba", t2_amplitudes, integral_matrix_oovv)

    return t2_amplitudes, energy_correction

[docs]class MP2InitialPoint(InitialPoint): """Compute the second-order Møller-Plesset perturbation theory (MP2) initial point. The computed MP2 correction coefficients are intended for use as an initial point for the VQE parameters in combination with a :class:`~qiskit_nature.second_q.circuit.library.ansatzes.ucc.UCC` ansatz. The initial point parameters are computed using the :meth:`compute` method, which requires the :attr:`problem` and :attr:`ansatz` to be passed as arguments or the :attr:`problem` and :attr:`ansatz` attributes to be set already. The :attr:`problem` is required to be an :class:`~qiskit_nature.second_q.problems.ElectronicStructureProblem` which contains an :class:`~qiskit_nature.second_q.hamiltonians.ElectronicEnergy` Hamiltonian. It also must have its ``num_particles`` and ``orbital_energies`` attributes specified. If its ``reference_energy`` attribute is provided, this will be used to compute the :attr:`total_energy`. :class:`~qiskit_nature.second_q.hamiltonians.ElectronicEnergy` must contain the two-body, molecular-orbital ``electronic_integrals``. Setting the :attr:`problem` will compute the :attr:`t2_amplitudes` and :attr:`energy_correction`. Following computation, one can obtain the initial point array via the :meth:`to_numpy_array` method. The initial point parameters that correspond to double excitations in the ``excitation_list`` will equal the appropriate T2 amplitude, while those below :attr:`threshold` or that correspond to single, triple, or higher excitations will be zero. """ def __init__(self, *, threshold: float = 1e-12) -> None: super().__init__() self.threshold: float = threshold self._ansatz: UCC | None = None self._t2_amplitudes: np.ndarray | None = None self._parameters: np.ndarray | None = None self._energy_correction: float = 0.0 self._total_energy: float = 0.0 @property def ansatz(self) -> UCC: """The UCC ansatz. The ``excitation_list`` and ``reps`` used by the :class:`~qiskit_nature.circuit.library.ansatzes.ucc.UCC` ansatz is obtained to ensure that the shape of the initial point is appropriate. """ return self._ansatz @ansatz.setter def ansatz(self, ansatz: UCC) -> None: self._invalidate() self._ansatz = ansatz @property def problem(self) -> BaseProblem | None: """The problem instance. See :class:`~qiskit_nature.second_q.algorithms.initial_points.MP2InitialPoint` for more information on the required properties. Raises: QiskitNatureError: If :class:`~qiskit_nature.second_q.hamiltonians.ElectronicEnergy` is missing or the two-body molecular orbitals matrix. QiskitNatureError: If :attr:`~qiskit_nature.second_q.problems.ElectronicStructureProblem.num_particles` or :attr:`~qiskit_nature.second_q.problems.ElectronicStructureProblem.orbital_energies` is ``None``. NotImplementedError: If alpha and beta spin molecular orbitals are not identical. """ return self._problem @problem.setter def problem(self, problem: BaseProblem) -> None: if not isinstance(problem, ElectronicStructureProblem): raise QiskitNatureError( "Only an `ElectronicStructureProblem` is compatible with the MP2InitialPoint, not a" f" problem of type, {type(problem)}." ) electronic_energy = problem.hamiltonian if electronic_energy is None: raise QiskitNatureError("The `ElectronicEnergy` cannot be obtained from the `problem`.") two_body_mo_integral: ElectronicIntegrals = electronic_energy.electronic_integrals.two_body if two_body_mo_integral.alpha.is_empty(): raise QiskitNatureError( "The alpha-alpha spin two-body MO `electronic_integrals` cannot be empty." ) orbital_energies: np.ndarray | None = problem.orbital_energies if orbital_energies is None: raise QiskitNatureError("The `orbital_energies` cannot be obtained from the `problem`.") if two_body_mo_integral.beta.get("++--", None) is not None: raise NotImplementedError( "`MP2InitialPoint` only supports restricted-spin setups. " "Alpha and beta spin orbitals must be identical. " "See" ) # Get number of occupied molecular orbitals as the number of alpha particles. # Only valid for restricted-spin setups. num_occ = problem.num_alpha if num_occ is None: raise QiskitNatureError( "The `num_particles` attribute of the `ElectronicStructureProblem` is required by " "the MP2InitialPoint." ) integral_matrix = two_body_mo_integral.alpha.get("++--") reference_energy = problem.reference_energy if not None else 0.0 self._invalidate() t2_amplitudes, energy_correction = _compute_mp2(num_occ, integral_matrix, orbital_energies) # Save state. self._problem = problem self._t2_amplitudes = t2_amplitudes self._energy_correction = energy_correction self._total_energy = reference_energy + energy_correction @property def t2_amplitudes(self) -> np.ndarray: """T2 amplitudes. Given ``t2[i, j, a, b]`` ``i, j`` carry virtual indices, while ``a, b`` carry occupied indices. """ return self._t2_amplitudes @property def energy_correction(self) -> float: """The MP2 energy correction.""" return self._energy_correction @property def total_energy(self) -> float: """The total energy including the Hartree-Fock energy. If the reference energy was not obtained from :class:`~qiskit_nature.second_q.hamiltonians.ElectronicEnergy` this will be equal to :attr:`energy_correction`. """ return self._total_energy @property def threshold(self) -> float: """Amplitudes below this vanish in the initial point array.""" return self._threshold @threshold.setter def threshold(self, threshold: float) -> None: try: threshold = abs(float(threshold)) except TypeError: threshold = 0.0 self._invalidate() self._threshold = threshold
[docs] def compute( self, ansatz: UCC | None = None, problem: BaseProblem | None = None, ) -> None: """Compute the initial point parameter for each excitation. See class documentation for more information. Args: problem: The :attr:`problem`. ansatz: The :attr:`ansatz`. Raises: QiskitNatureError: If :attr:`ansatz` or :attr:`problem` is not set. """ if ansatz is not None: self.ansatz = ansatz if self._ansatz is None: raise QiskitNatureError( "The ansatz property has not been set. " "Not enough information has been provided to compute the initial point. " "Set the ansatz or call compute with it as an argument. " ) if problem is not None: self.problem = problem if self._problem is None: raise QiskitNatureError("The problem has not been set.") self._compute()
def _compute(self) -> None: """Compute the MP2 amplitudes given an excitation list. Non-double excitations will have zero coefficient. Returns: The MP2 T2 amplitudes for each excitation. """ # Operators must be built to compute the excitation list. _ = self._ansatz.operators num_occ = self._t2_amplitudes.shape[0] amplitudes = np.zeros(len(self._ansatz.excitation_list), dtype=float) for index, excitation in enumerate(self._ansatz.excitation_list): if len(excitation[0]) == 2: # Get the amplitude of the double excitation. [[i, j], [a, b]] = np.asarray(excitation) % num_occ amplitude = self._t2_amplitudes[i, j, a - num_occ, b - num_occ] amplitudes[index] = amplitude if abs(amplitude) > self._threshold else 0.0 self._parameters = np.tile(amplitudes, self._ansatz.reps)
[docs] def to_numpy_array(self) -> np.ndarray: """The initial point as a NumPy array.""" if self._parameters is None: self.compute() return self._parameters
def _invalidate(self): """Invalidate any previous computation.""" self._parameters = None