Line data    Source code

       1             : # SPDX-FileCopyrightText: 2025 PairInteraction Developers
       2             : # SPDX-License-Identifier: LGPL-3.0-or-later
       3           1 : from __future__ import annotations
       4             : 
       5           1 : import logging
       6           1 : from typing import TYPE_CHECKING
       7             : 
       8           1 : import numpy as np
       9           1 : from scipy import sparse
      10           1 : from scipy.sparse import csr_matrix
      11             : 
      12             : if TYPE_CHECKING:
      13             :     from scipy.sparse import csr_matrix
      14             : 
      15             :     from pairinteraction.units import NDArray
      16             : 
      17             : 
      18           1 : logger = logging.getLogger(__name__)
      19             : 
      20             : 
      21           1 : def calculate_perturbative_hamiltonian(
      22             :     hamiltonian: csr_matrix,
      23             :     model_inds: list[int],
      24             :     perturbation_order: int,
      25             : ) -> tuple[dict[int, NDArray], csr_matrix]:
      26             :     r"""Calculate the perturbative Hamiltonian up to a given order.
      27             : 
      28             :     This function calculates both the effective Hamiltonian, spanned up by the states of the model space,
      29             :     as well as the effective (perturbed) basisvectors due to interactions with the exterior space
      30             :     up to the desired order of perturbation theory.
      31             : 
      32             :     Args:
      33             :         hamiltonian: Hermitian matrix.
      34             :         model_inds: List of indices corresponding to the states that span up the model space.
      35             :         perturbation_order: Order up to which the perturbation theory is expanded.
      36             :             Support up to third order. Default is second order.
      37             : 
      38             :     Returns:
      39             :         Effective hamiltonians as dictionaries mapping perturbation orders to :math:`m \times m` matrices.
      40             :         Effective basisvectors in perturbation theory, returned as a sparse matrix. Each row represents one basisvector.
      41             : 
      42             :     """
      43           1 :     m_inds = np.asarray(model_inds, dtype=int)
      44           1 :     o_inds = np.setdiff1d(np.arange(hamiltonian.shape[0]), m_inds)
      45           1 :     eff_h_dict, eff_vecs = _calculate_unsorted_perturbative_hamiltonian(hamiltonian, m_inds, o_inds, perturbation_order)
      46             : 
      47             :     # resort eigvec to original order
      48           1 :     all_inds = np.append(m_inds, o_inds)
      49           1 :     all_inds_positions = np.argsort(all_inds)
      50           1 :     eff_vecs = eff_vecs[:, all_inds_positions]
      51             : 
      52             :     # include the hermitian conjugate part of the effective Hamiltonian
      53           1 :     for order, h_eff in eff_h_dict.items():
      54           1 :         eff_h_dict[order] = 0.5 * (h_eff + h_eff.conj().T)
      55             : 
      56           1 :     return eff_h_dict, eff_vecs
      57             : 
      58             : 
      59           1 : def _calculate_unsorted_perturbative_hamiltonian(
      60             :     hamiltonian: csr_matrix,
      61             :     m_inds: NDArray,
      62             :     o_inds: NDArray,
      63             :     perturbation_order: int,
      64             : ) -> tuple[dict[int, NDArray], csr_matrix]:
      65             :     # This function is outsourced from calculate_perturbative_hamiltonian to allow for better type checking
      66           1 :     energies = np.real_if_close(hamiltonian.diagonal())
      67           1 :     if any(np.iscomplex(energies)):
      68           0 :         logger.error("The Hamiltonian has complex entries on the diagonal, this might lead to unexpected results.")
      69             : 
      70           1 :     energies_m = energies[m_inds]
      71             : 
      72           1 :     eff_h_dict: dict[int, NDArray] = {}  # perturbation_order -> h_eff
      73             : 
      74           1 :     eff_h_dict[0] = np.diag(energies_m)
      75           1 :     eff_vecs = sparse.csr_matrix(
      76             :         sparse.hstack(
      77             :             [sparse.eye(len(m_inds), len(m_inds), format="csr"), sparse.csr_matrix((len(m_inds), len(o_inds)))]
      78             :         )
      79             :     )
      80             : 
      81           1 :     if perturbation_order == 0:
      82           1 :         return eff_h_dict, eff_vecs
      83             : 
      84           1 :     v_offdiag = hamiltonian - sparse.diags(energies, dtype=float)
      85           1 :     v_mm = v_offdiag[np.ix_(m_inds, m_inds)]  # type: ignore [index]
      86           1 :     eff_h_dict[1] = v_mm.toarray()
      87             : 
      88           1 :     if perturbation_order == 1:
      89           1 :         return eff_h_dict, eff_vecs
      90             : 
      91           1 :     energies_diff = energies_m[np.newaxis, :] - energies[o_inds, np.newaxis]
      92           1 :     with np.errstate(divide="ignore"):
      93           1 :         delta_e_em = 1 / energies_diff
      94           1 :     v_me = v_offdiag[np.ix_(m_inds, o_inds)]  # type: ignore [index]
      95           1 :     eff_h_dict[2] = (v_me @ ((v_me.conj().T).multiply(delta_e_em))).toarray()
      96             : 
      97           1 :     addition_mm = sparse.csr_matrix((len(m_inds), len(m_inds)))
      98           1 :     addition_me = sparse.csr_matrix(((v_me.conj().T).multiply(delta_e_em)).T)
      99           1 :     if np.isinf(np.isinf(addition_me.data)).any():
     100           0 :         logger.critical(
     101             :             "Detected 'inf' entries in the effective basisvectors. "
     102             :             "This might happen, if you forgot to include a degenerate state in the model space. "
     103             :         )
     104             : 
     105           1 :     nan_idx = np.where(np.isnan(addition_me.data))[0]
     106           1 :     nonzero_inds = addition_me.nonzero()
     107           1 :     for i in nan_idx:
     108           0 :         value = addition_me.data[i]
     109           0 :         if np.isinf(value.real):
     110           0 :             row, col = nonzero_inds[0][i], nonzero_inds[1][i]
     111           0 :             addition_me[row, col] = value.real
     112             :         else:
     113           0 :             logger.error("Detected unexpected 'nan' entries in the effective basisvectors.")
     114             : 
     115           1 :     eff_vecs = eff_vecs + sparse.hstack([addition_mm, addition_me])
     116             : 
     117           1 :     if perturbation_order == 2:
     118           1 :         return eff_h_dict, eff_vecs
     119             : 
     120           1 :     diff = energies_m[np.newaxis, :] - energies_m[:, np.newaxis]
     121           1 :     diff = np.where(diff == 0, np.inf, diff)
     122           1 :     delta_e_mm = 1 / diff
     123           1 :     v_ee = v_offdiag[np.ix_(o_inds, o_inds)]  # type: ignore [index]
     124           1 :     if len(m_inds) > 1:
     125           1 :         logger.debug(
     126             :             "At third order, the effective basisvectors are currently only valid, "
     127             :             "when only one state is in the model space. "
     128             :             "Take care with interpretation of the effective basisvectors."
     129             :         )
     130           1 :     _intermediate_term = v_ee @ ((v_me.conj().T).multiply(delta_e_em)) - ((v_me.conj().T).multiply(delta_e_em)) @ v_mm
     131           1 :     eff_h_dict[3] = (v_me @ (_intermediate_term.multiply(delta_e_em))).toarray()
     132             : 
     133           1 :     addition_mm_diag = -0.5 * sparse.csr_matrix(
     134             :         sparse.diags((v_me @ ((v_me.conj().T).multiply(np.square(delta_e_em)))).diagonal())
     135             :     )
     136           1 :     addition_mm_offdiag = sparse.csr_matrix(((v_me @ (v_me.conj().T).multiply(delta_e_em)).multiply(delta_e_mm)).T)
     137           1 :     addition_me = sparse.csr_matrix(((v_ee @ ((v_me.conj().T).multiply(delta_e_em))).multiply(delta_e_em)).T)
     138           1 :     addition_me_2 = sparse.csr_matrix(((v_me.conj().T @ ((v_mm.conj().T).multiply(delta_e_mm))).multiply(delta_e_em)).T)
     139           1 :     eff_vecs = eff_vecs + sparse.hstack([addition_mm_diag + addition_mm_offdiag, addition_me + addition_me_2])
     140             : 
     141           1 :     if perturbation_order == 3:
     142           1 :         return eff_h_dict, eff_vecs
     143             : 
     144           0 :     raise ValueError("Perturbation theory is only implemented for orders [0, 1, 2, 3].")