Source code for qiskit_nature.second_q.problems.vibrational_basis

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

"""The Vibrational basis base class."""

from __future__ import annotations

from abc import ABC, abstractmethod
from collections import Counter
from functools import lru_cache
from itertools import chain, cycle, permutations, product
from typing import Generator

import numpy as np

[docs]class VibrationalBasis(ABC): """The Vibrational basis base class. This class defines the interface which any vibrational basis must implement. A basis must be applied to the vibrational integrals in order to map them into a second-quantization form. The following attributes can be set via the initializer but can also be read and updated once the ``VibrationalBasis`` object has been constructed. Attributes: num_modals (list[int]): the number of modals into which each mode gets expanded in second-quantization. threshold (float): the threshold below which integral values will be dropped. """ def __init__( self, num_modals: list[int], *, threshold: float = 1e-6, ) -> None: """ Args: num_modals: the number of modals to be used for each mode. threshold: the threshold value below which an integral coefficient gets neglected. """ self.num_modals = num_modals self.threshold = threshold
[docs] @abstractmethod @lru_cache(maxsize=128) def eval_integral( self, mode: int, modal_1: int, modal_2: int, power: int, kinetic_term: bool = False, ) -> complex | None: """The integral evaluation method of this basis. Args: mode: the index of the mode. modal_1: the index of the first modal. modal_2: the index of the second modal. power: the exponent of the coordinate. kinetic_term: if this is True, the method should compute the integral of the kinetic term of the vibrational Hamiltonian, :math:``d^2/dQ^2``. Returns: The evaluated integral for the specified coordinate or ``None`` if this integral value falls below the threshold. Raises: ValueError: if an unsupported parameter is supplied. """
[docs] def map( self, coefficient: complex, modes: tuple[int, ...] ) -> Generator[tuple[complex, tuple[int, ...]], None, None]: """Maps the provided coefficient and mode index to this second-quantization basis. This applies the actual basis and expands each mode into the number of modals with which the basis instance was initialized. Args: coefficient: the initial coefficient associated with the mode indices. modes: the mode indices. If all of these are negative, the coefficient is treated as belonging to a kinetic term. Yields: Pairs of integral values and indices. The indices are now three times as long as the initially provided modes index. The reason for that is that each mode index gets expanded into three indices, denoting the ``(mode, modal_1, modal_2)`` indices. """ # negative indices may be treated specially by a basis kinetic_term = any(i < 0 for i in modes) # the number of times which an index occurs corresponds to the power of the operator powers = Counter(abs(i) for i in modes) # we generate the list of all possible modal permutations (lower triangular indices) for all # involved modes: each entry in this list is a list of tuples of the form: # (mode_index, modal_index_1, modal_index_2) index_list = list( zip(cycle([mode]), *np.tril_indices(self.num_modals[mode - 1])) for mode in powers ) # now we can iterate the product of all index lists (the cartesian product is equivalent to # nested for loops but has the benefit of being agnostic w.r.t. the number of body terms) for index in product(*index_list): index_permutations = [] coeff = coefficient for mode, m, n in index: integral = self.eval_integral( mode - 1, m, n, powers[mode], kinetic_term=kinetic_term ) if integral is None: break coeff *= integral # generate potentially symmetric permutations of the modal indices index_permutations.append( {(mode - 1, m_sub, n_sub) for (m_sub, n_sub) in permutations((m, n))} ) else: # update the matrix in all permuted locations for i in product(*index_permutations): yield (coeff, tuple(chain(*i)))