bse_uks.cc

Go to the documentation of this file.

/*
 *            Copyright 2009-2026 The VOTCA Development Team
 *                       (http://www.votca.org)
 *
 *      Licensed under the Apache License, Version 2.0 (the "License")
 *
 * You may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *              http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
 
#include <chrono>
#include <iostream>
 
#include <votca/tools/linalg.h>
 
#include "votca/xtp/bse_fullsolver.h"
#include "votca/xtp/bse_initialization.h"
#include "votca/xtp/bse_operator_uks.h"
#include "votca/xtp/bse_uks.h"
#include "votca/xtp/bseoperator_btda.h"
#include "votca/xtp/davidsonsolver.h"
#include "votca/xtp/rpa_uks.h"
 
using std::flush;
 
namespace {
 
template <class OP>
Eigen::VectorXd ExpValue(const Eigen::MatrixXd& state1, OP OPxstate2) {
  return state1.cwiseProduct(OPxstate2.eval()).colwise().sum().transpose();
}
 
Eigen::VectorXd ExpValue(const Eigen::MatrixXd& state1,
                         const Eigen::MatrixXd& OPxstate2) {
  return state1.cwiseProduct(OPxstate2).colwise().sum().transpose();
}
 
}  // namespace
 
namespace votca {
namespace xtp {
 
namespace {
double DavidsonToleranceValue(const std::string& tol) {
  if (tol == "loose") {
    return 1e-3;
  } else if (tol == "normal") {
    return 1e-4;
  } else if (tol == "strict") {
    return 1e-5;
  } else if (tol == "lapack") {
    return 1e-9;
  }
  throw std::runtime_error(tol + " is not a valid Davidson tolerance");
}
}  // namespace
 
void BSE_UKS::configure_with_precomputed_screening(
    const options& opt, Index homo_alpha, Index homo_beta,
    const Eigen::VectorXd& RPAInputEnergiesAlpha,
    const Eigen::VectorXd& RPAInputEnergiesBeta,
    const Eigen::MatrixXd& Hqp_alpha_in, const Eigen::MatrixXd& Hqp_beta_in,
    const Eigen::VectorXd& epsilon_0_inv,
    const Eigen::MatrixXd& epsilon_eigenvectors) {
 
  opt_ = opt;
  homo_alpha_ = homo_alpha;
  homo_beta_ = homo_beta;
 
  alpha_vtotal_ = homo_alpha_ - opt_.vmin + 1;
  alpha_ctotal_ = opt_.cmax - (homo_alpha_ + 1) + 1;
  alpha_size_ = alpha_vtotal_ * alpha_ctotal_;
 
  beta_vtotal_ = homo_beta_ - opt_.vmin + 1;
  beta_ctotal_ = opt_.cmax - (homo_beta_ + 1) + 1;
  beta_size_ = beta_vtotal_ * beta_ctotal_;
 
  if (opt_.use_Hqp_offdiag) {
    Hqp_alpha_ =
        AdjustHqpSize(Hqp_alpha_in, RPAInputEnergiesAlpha, homo_alpha_);
    Hqp_beta_ = AdjustHqpSize(Hqp_beta_in, RPAInputEnergiesBeta, homo_beta_);
  } else {
    Hqp_alpha_ = AdjustHqpSize(Hqp_alpha_in, RPAInputEnergiesAlpha, homo_alpha_)
                     .diagonal()
                     .asDiagonal();
    Hqp_beta_ = AdjustHqpSize(Hqp_beta_in, RPAInputEnergiesBeta, homo_beta_)
                    .diagonal()
                    .asDiagonal();
  }
 
  // Store the statically screened interaction in the working copy Mmn_.
  Mmn_ = Mmn_raw_;
  Mmn_.alpha.MultiplyRightWithAuxMatrix(epsilon_eigenvectors);
  Mmn_.beta.MultiplyRightWithAuxMatrix(epsilon_eigenvectors);
 
  epsilon_0_inv_ = epsilon_0_inv;
}
 
Eigen::MatrixXd BSE_UKS::AdjustHqpSize(const Eigen::MatrixXd& Hqp,
                                       const Eigen::VectorXd& RPAInputEnergies,
                                       Index homo) const {
  Index bse_vtotal = homo - opt_.vmin + 1;
  Index bse_ctotal = opt_.cmax - (homo + 1) + 1;
  Index hqp_size = bse_vtotal + bse_ctotal;
  Index gwsize = opt_.qpmax - opt_.qpmin + 1;
  Index RPAoffset = opt_.vmin - opt_.rpamin;
  Eigen::MatrixXd Hqp_BSE = Eigen::MatrixXd::Zero(hqp_size, hqp_size);
 
  if (opt_.vmin >= opt_.qpmin) {
    Index start = opt_.vmin - opt_.qpmin;
    if (opt_.cmax <= opt_.qpmax) {
      Hqp_BSE = Hqp.block(start, start, hqp_size, hqp_size);
    } else {
      Index virtoffset = gwsize - start;
      Hqp_BSE.topLeftCorner(virtoffset, virtoffset) =
          Hqp.block(start, start, virtoffset, virtoffset);
 
      Index virt_extra = opt_.cmax - opt_.qpmax;
      Hqp_BSE.diagonal().tail(virt_extra) =
          RPAInputEnergies.segment(RPAoffset + virtoffset, virt_extra);
    }
  }
 
  if (opt_.vmin < opt_.qpmin) {
    Index occ_extra = opt_.qpmin - opt_.vmin;
    Hqp_BSE.diagonal().head(occ_extra) =
        RPAInputEnergies.segment(RPAoffset, occ_extra);
 
    Hqp_BSE.block(occ_extra, occ_extra, gwsize, gwsize) = Hqp;
 
    if (opt_.cmax > opt_.qpmax) {
      Index virtoffset = occ_extra + gwsize;
      Index virt_extra = opt_.cmax - opt_.qpmax;
      Hqp_BSE.diagonal().tail(virt_extra) =
          RPAInputEnergies.segment(RPAoffset + virtoffset, virt_extra);
    }
  }
 
  return Hqp_BSE;
}
 
void BSE_UKS::SetupDirectInteractionOperator(
    const Eigen::VectorXd& RPAInputEnergiesAlpha,
    const Eigen::VectorXd& RPAInputEnergiesBeta, double energy) {
 
  RPA_UKS rpa(log_, Mmn_raw_);
  rpa.configure(homo_alpha_, homo_beta_, opt_.rpamin, opt_.rpamax);
  rpa.setRPAInputEnergies(RPAInputEnergiesAlpha, RPAInputEnergiesBeta);
 
  const SymmetricEigenSystem es =
      SymmetricEigen(rpa.calculate_epsilon_r(energy));
 
  Mmn_ = Mmn_raw_;
  Mmn_.alpha.MultiplyRightWithAuxMatrix(es.vectors);
  Mmn_.beta.MultiplyRightWithAuxMatrix(es.vectors);
 
  epsilon_0_inv_ = Eigen::VectorXd::Zero(es.values.size());
  for (Index i = 0; i < es.values.size(); ++i) {
    if (es.values(i) > 1e-8) {
      epsilon_0_inv_(i) = 1.0 / es.values(i);
    }
  }
}
 
template <typename BSE_OPERATOR>
void BSE_UKS::configureBSEOperator(BSE_OPERATOR& H) const {
  BSEOperatorUKS_Options opt;
  opt.homo_alpha = homo_alpha_;
  opt.homo_beta = homo_beta_;
  opt.rpamin = opt_.rpamin;
  opt.qpmin = opt_.qpmin;
  opt.vmin = opt_.vmin;
  opt.cmax = opt_.cmax;
  H.configure(opt);
}
 
template <typename BSE_OPERATOR>
BSE_UKS::ExpectationValues BSE_UKS::ExpectationValue_Operator(
    const Orbitals& orb, const BSE_OPERATOR& H) const {
 
  const tools::EigenSystem& BSECoefs = orb.BSEUKS();
 
  ExpectationValues expectation_values;
 
  const Eigen::MatrixXd temp = H * BSECoefs.eigenvectors();
 
  expectation_values.direct_term = ExpValue(BSECoefs.eigenvectors(), temp);
 
  if (!orb.getTDAApprox()) {
    expectation_values.direct_term +=
        ExpValue(BSECoefs.eigenvectors2(), H * BSECoefs.eigenvectors2());
    expectation_values.cross_term =
        2.0 * ExpValue(BSECoefs.eigenvectors2(), temp);
  } else {
    expectation_values.cross_term = Eigen::VectorXd::Zero(0);
  }
 
  return expectation_values;
}
 
template <typename BSE_OPERATOR>
BSE_UKS::ExpectationValues BSE_UKS::ExpectationValue_Operator_State(
    Index state, const Orbitals& orb, const BSE_OPERATOR& H) const {
 
  const tools::EigenSystem& BSECoefs = orb.BSEUKS();
 
  ExpectationValues expectation_values;
 
  const Eigen::MatrixXd Xstate = BSECoefs.eigenvectors().col(state);
  const Eigen::MatrixXd temp = H * Xstate;
 
  expectation_values.direct_term = ExpValue(Xstate, temp);
 
  if (!orb.getTDAApprox()) {
    const Eigen::MatrixXd Ystate = BSECoefs.eigenvectors2().col(state);
    expectation_values.direct_term += ExpValue(Ystate, H * Ystate);
    expectation_values.cross_term = 2.0 * ExpValue(Ystate, temp);
  } else {
    expectation_values.cross_term = Eigen::VectorXd::Zero(0);
  }
 
  return expectation_values;
}
 
tools::EigenSystem BSE_UKS::Solve_excitons_uks_TDA() const {
  ExcitonUKSOperator_TDA H(epsilon_0_inv_, Mmn_, Hqp_alpha_, Hqp_beta_);
  BSEOperatorUKS_Options opt;
  opt.homo_alpha = homo_alpha_;
  opt.homo_beta = homo_beta_;
  opt.rpamin = opt_.rpamin;
  opt.qpmin = opt_.qpmin;
  opt.vmin = opt_.vmin;
  opt.cmax = opt_.cmax;
  H.configure(opt);
 
  XTP_LOG(Log::error, log_)
      << TimeStamp() << " Setup combined UKS TDA Hamiltonian " << flush;
  return solve_hermitian(H);
}
 
/* tools::EigenSystem BSE_UKS::Solve_excitons_uks_BTDA() const {
  ExcitonUKSOperator_TDA A(epsilon_0_inv_, Mmn_, Hqp_alpha_, Hqp_beta_);
  ExcitonUKSOperator_BTDA_B B(epsilon_0_inv_, Mmn_, Hqp_alpha_, Hqp_beta_);
 
  BSEOperatorUKS_Options opt;
  opt.homo_alpha = homo_alpha_;
  opt.homo_beta = homo_beta_;
  opt.rpamin = opt_.rpamin;
  opt.qpmin = opt_.qpmin;
  opt.vmin = opt_.vmin;
  opt.cmax = opt_.cmax;
 
  A.configure(opt);
  B.configure(opt);
 
  XTP_LOG(Log::error, log_)
      << TimeStamp() << " Setup combined UKS full BSE exciton hamiltonian "
      << flush;
 
  // Cheap production initial guess:
  // rank excitations by local diagonal full-BSE estimate
  //   omega_i = sqrt(max(0, A_ii^2 - B_ii^2))
  // but keep the actual start vectors X-only for robustness with the
  // current HAM Davidson implementation.
  const Eigen::VectorXd adiag = A.diagonal();
  const Eigen::VectorXd bdiag = B.diagonal();
  Eigen::MatrixXd initial_guess =
      BuildFullBSEXRankedInitialGuess(adiag, bdiag, opt_.nmax);
 
  HamiltonianOperator<ExcitonUKSOperator_TDA, ExcitonUKSOperator_BTDA_B> Hop(A,
                                                                             B);
 
  DavidsonSolver DS(log_);
  DS.set_correction(opt_.davidson_correction);
  DS.set_tolerance(opt_.davidson_tolerance);
  DS.set_size_update(opt_.davidson_update);
  DS.set_iter_max(opt_.davidson_maxiter);
  DS.set_max_search_space(10 * opt_.nmax);
  DS.set_matrix_type("HAM");
  DS.solve(Hop, opt_.nmax, initial_guess);
 
  tools::EigenSystem result;
  result.eigenvalues() = DS.eigenvalues();
 
  Eigen::MatrixXd tmpX = DS.eigenvectors().topRows(A.rows());
  Eigen::MatrixXd tmpY = DS.eigenvectors().bottomRows(B.rows());
 
  Eigen::VectorXd normX = tmpX.colwise().squaredNorm();
  Eigen::VectorXd normY = tmpY.colwise().squaredNorm();
  Eigen::ArrayXd sqinvnorm = (normX - normY).array().inverse().cwiseSqrt();
 
  result.eigenvectors() = tmpX * sqinvnorm.matrix().asDiagonal();
  result.eigenvectors2() = tmpY * sqinvnorm.matrix().asDiagonal();
 
  return result;
}*/
 
tools::EigenSystem BSE_UKS::Solve_excitons_uks_BTDA() const {
  ExcitonUKSOperator_TDA A(epsilon_0_inv_, Mmn_, Hqp_alpha_, Hqp_beta_);
  ExcitonUKSOperator_BTDA_B B(epsilon_0_inv_, Mmn_, Hqp_alpha_, Hqp_beta_);
 
  BSEOperatorUKS_Options opt;
  opt.homo_alpha = homo_alpha_;
  opt.homo_beta = homo_beta_;
  opt.rpamin = opt_.rpamin;
  opt.qpmin = opt_.qpmin;
  opt.vmin = opt_.vmin;
  opt.cmax = opt_.cmax;
 
  A.configure(opt);
  B.configure(opt);
 
  XTP_LOG(Log::error, log_)
      << TimeStamp() << " Setup combined UKS full BSE exciton hamiltonian "
      << flush;
 
  // Dense fallback for small systems:
  // this avoids architecture-sensitive Davidson/HAM behavior in tiny
  // regression tests while keeping the iterative solver for production sizes.
  constexpr Index dense_threshold = 128;
 
  if (A.rows() <= dense_threshold) {
    XTP_LOG(Log::error, log_)
        << TimeStamp()
        << " Using dense full UKS-BSE solve for small system (dim=" << A.rows()
        << ")" << flush;
    std::chrono::time_point<std::chrono::system_clock> start =
        std::chrono::system_clock::now();
 
    const Eigen::MatrixXd Ad = A.dense_matrix();
    const Eigen::MatrixXd Bd = B.dense_matrix();
    const Index n = Ad.rows();
 
    Eigen::MatrixXd Hfull = Eigen::MatrixXd::Zero(2 * n, 2 * n);
    Hfull.topLeftCorner(n, n) = Ad;
    Hfull.topRightCorner(n, n) = Bd;
    Hfull.bottomLeftCorner(n, n) = -Bd;
    Hfull.bottomRightCorner(n, n) = -Ad;
 
    Eigen::EigenSolver<Eigen::MatrixXd> es(Hfull, true);
    if (es.info() != Eigen::Success) {
      throw std::runtime_error("Dense full UKS-BSE diagonalization failed.");
    }
 
    using RootInfo =
        std::tuple<double, Index, double>;  // (eval, col, imag_abs)
    std::vector<RootInfo> positive_roots;
    positive_roots.reserve(2 * n);
 
    for (Index i = 0; i < 2 * n; ++i) {
      const std::complex<double> ev = es.eigenvalues()(i);
      if (std::abs(ev.imag()) < 1e-8 && ev.real() > 0.0) {
        positive_roots.emplace_back(ev.real(), i, std::abs(ev.imag()));
      }
    }
 
    std::sort(positive_roots.begin(), positive_roots.end(),
              [](const RootInfo& a, const RootInfo& b) {
                return std::get<0>(a) < std::get<0>(b);
              });
 
    if (positive_roots.empty()) {
      throw std::runtime_error(
          "Dense full UKS-BSE diagonalization produced no positive real "
          "roots.");
    }
 
    const Index nroots =
        std::min<Index>(opt_.nmax, static_cast<Index>(positive_roots.size()));
 
    tools::EigenSystem result;
    result.eigenvalues() = Eigen::VectorXd::Zero(nroots);
    result.eigenvectors() = Eigen::MatrixXd::Zero(n, nroots);
    result.eigenvectors2() = Eigen::MatrixXd::Zero(n, nroots);
 
    for (Index iroot = 0; iroot < nroots; ++iroot) {
      const Index col = std::get<1>(positive_roots[iroot]);
      Eigen::VectorXd vec = es.eigenvectors().col(col).real();
 
      Eigen::VectorXd X = vec.head(n);
      Eigen::VectorXd Y = vec.tail(n);
 
      // Deterministic phase convention: make the largest |X_k| positive.
      Eigen::Index imax = 0;
      X.cwiseAbs().maxCoeff(&imax);
      if (X(imax) < 0.0) {
        X *= -1.0;
        Y *= -1.0;
      }
 
      const double norm = X.squaredNorm() - Y.squaredNorm();
      if (std::abs(norm) < 1e-12) {
        throw std::runtime_error(
            "Dense full UKS-BSE eigenvector has near-zero (X^2-Y^2) norm.");
      }
 
      const double invnorm = 1.0 / std::sqrt(std::abs(norm));
 
      result.eigenvalues()(iroot) = std::get<0>(positive_roots[iroot]);
      result.eigenvectors().col(iroot) = X * invnorm;
      result.eigenvectors2().col(iroot) = Y * invnorm;
    }
 
    std::chrono::time_point<std::chrono::system_clock> end =
        std::chrono::system_clock::now();
    std::chrono::duration<double> elapsed_time = end - start;
 
    XTP_LOG(Log::info, log_)
        << TimeStamp() << " Dense full UKS-BSE diagonalization done in "
        << elapsed_time.count() << " secs" << flush;
 
    return result;
  }
 
  // Hermitian (A-B)(A+B) form for larger systems; the non-Hermitian
  // Davidson below if A-B or A+B is not positive definite.
  try {
    FullBSEDavidson::Options fopt;
    fopt.tolerance = DavidsonToleranceValue(opt_.davidson_tolerance);
    fopt.max_iterations = opt_.davidson_maxiter;
    fopt.max_subspace = 10 * opt_.nmax;
    FullBSEDavidson solver(log_, fopt);
    return solver.Solve(A, B, opt_.nmax);
  } catch (const std::runtime_error& error) {
    XTP_LOG(Log::error, log_)
        << TimeStamp() << " " << error.what()
        << "; using the non-Hermitian Davidson solver" << flush;
  }
 
  const Eigen::VectorXd adiag = A.diagonal();
  const Eigen::VectorXd bdiag = B.diagonal();
  Eigen::MatrixXd initial_guess =
      BuildFullBSEXRankedInitialGuess(adiag, bdiag, opt_.nmax);
 
  HamiltonianOperator<ExcitonUKSOperator_TDA, ExcitonUKSOperator_BTDA_B> Hop(A,
                                                                             B);
 
  DavidsonSolver DS(log_);
  DS.set_correction(opt_.davidson_correction);
  DS.set_tolerance(opt_.davidson_tolerance);
  DS.set_size_update(opt_.davidson_update);
  DS.set_iter_max(opt_.davidson_maxiter);
  DS.set_max_search_space(10 * opt_.nmax);
  DS.set_matrix_type("HAM");
  DS.solve(Hop, opt_.nmax, initial_guess);
 
  tools::EigenSystem result;
  result.eigenvalues() = DS.eigenvalues();
 
  Eigen::MatrixXd tmpX = DS.eigenvectors().topRows(A.rows());
  Eigen::MatrixXd tmpY = DS.eigenvectors().bottomRows(B.rows());
 
  if (const Index bad = NormalizeExcitationVectors(tmpX, tmpY); bad > 0) {
    XTP_LOG(Log::error, log_)
        << TimeStamp() << " WARNING: " << bad
        << " BSE root(s) without positive X^2 - Y^2 norm (not converged?)"
        << flush;
  }
  result.eigenvectors() = tmpX;
  result.eigenvectors2() = tmpY;
 
  return result;
}
 
void BSE_UKS::Solve_excitons_uks(Orbitals& orb) const {
  orb.setTDAApprox(opt_.useTDA);
  if (opt_.useTDA) {
    orb.BSEUKS() = Solve_excitons_uks_TDA();
  } else {
    orb.BSEUKS() = Solve_excitons_uks_BTDA();
  }
  orb.CalcCoupledTransition_Dipoles(QMStateType(QMStateType::ExcitonUKS));
}
 
template <typename BSE_OPERATOR>
tools::EigenSystem BSE_UKS::solve_hermitian(BSE_OPERATOR& h) const {
  std::chrono::time_point<std::chrono::system_clock> start =
      std::chrono::system_clock::now();
 
  tools::EigenSystem result;
 
  DavidsonSolver DS(log_);
  DS.set_correction(opt_.davidson_correction);
  DS.set_tolerance(opt_.davidson_tolerance);
  DS.set_size_update(opt_.davidson_update);
  DS.set_iter_max(opt_.davidson_maxiter);
  DS.set_max_search_space(10 * opt_.nmax);
  DS.solve(h, opt_.nmax);
 
  result.eigenvalues() = DS.eigenvalues();
  result.eigenvectors() = DS.eigenvectors();
 
  std::chrono::time_point<std::chrono::system_clock> end =
      std::chrono::system_clock::now();
  std::chrono::duration<double> elapsed_time = end - start;
 
  XTP_LOG(Log::info, log_) << TimeStamp() << " Diagonalization done in "
                           << elapsed_time.count() << " secs" << flush;
 
  return result;
}
 
template <typename BSE_OPERATOR_A, typename BSE_OPERATOR_B>
tools::EigenSystem BSE_UKS::Solve_nonhermitian_Davidson(
    BSE_OPERATOR_A& Aop, BSE_OPERATOR_B& Bop) const {
  std::chrono::time_point<std::chrono::system_clock> start =
      std::chrono::system_clock::now();
 
  HamiltonianOperator<BSE_OPERATOR_A, BSE_OPERATOR_B> Hop(Aop, Bop);
 
  DavidsonSolver DS(log_);
  DS.set_correction(opt_.davidson_correction);
  DS.set_tolerance(opt_.davidson_tolerance);
  DS.set_size_update(opt_.davidson_update);
  DS.set_iter_max(opt_.davidson_maxiter);
  DS.set_max_search_space(10 * opt_.nmax);
  DS.set_matrix_type("HAM");
  DS.solve(Hop, opt_.nmax);
 
  tools::EigenSystem result;
  result.eigenvalues() = DS.eigenvalues();
 
  Eigen::MatrixXd tmpX = DS.eigenvectors().topRows(Aop.rows());
  Eigen::MatrixXd tmpY = DS.eigenvectors().bottomRows(Bop.rows());
 
  if (const Index bad = NormalizeExcitationVectors(tmpX, tmpY); bad > 0) {
    XTP_LOG(Log::error, log_)
        << TimeStamp() << " WARNING: " << bad
        << " BSE root(s) without positive X^2 - Y^2 norm (not converged?)"
        << flush;
  }
  result.eigenvectors() = tmpX;
  result.eigenvectors2() = tmpY;
 
  std::chrono::time_point<std::chrono::system_clock> end =
      std::chrono::system_clock::now();
  std::chrono::duration<double> elapsed_time = end - start;
 
  XTP_LOG(Log::info, log_) << TimeStamp() << " Diagonalization done in "
                           << elapsed_time.count() << " secs" << flush;
 
  return result;
}
 
void BSE_UKS::PrintWeightsUKS(const Eigen::VectorXd& weights) const {
  struct Contribution {
    double weight;
    bool is_alpha;
    Index v;
    Index c;
  };
 
  constexpr Index max_print = 8;
 
  double alpha_weight = weights.head(alpha_size_).sum();
  double beta_weight = weights.tail(beta_size_).sum();
 
  XTP_LOG(Log::error, log_)
      << boost::format(
             "           alpha-sector: %1$+6.2f%%   beta-sector: %2$+6.2f%%") %
             (100.0 * alpha_weight) % (100.0 * beta_weight)
      << flush;
 
  std::vector<Contribution> contributions;
  contributions.reserve(alpha_size_ + beta_size_);
 
  for (Index i = 0; i < alpha_size_; ++i) {
    if (std::abs(weights(i)) > opt_.min_print_weight) {
      contributions.push_back(
          {weights(i), true, i / alpha_ctotal_, i % alpha_ctotal_});
    }
  }
 
  for (Index i = 0; i < beta_size_; ++i) {
    const double w = weights(alpha_size_ + i);
    if (std::abs(w) > opt_.min_print_weight) {
      contributions.push_back({w, false, i / beta_ctotal_, i % beta_ctotal_});
    }
  }
 
  std::sort(contributions.begin(), contributions.end(),
            [](const Contribution& a, const Contribution& b) {
              return std::abs(a.weight) > std::abs(b.weight);
            });
 
  const Index nprint =
      std::min<Index>(max_print, static_cast<Index>(contributions.size()));
 
  for (Index i = 0; i < nprint; ++i) {
    const Contribution& c = contributions[i];
 
    if (c.is_alpha) {
      XTP_LOG(Log::error, log_)
          << boost::format(
                 "           [alpha] HOMO-%1$-3d -> LUMO+%2$-3d  : "
                 "%3$+6.2f%%") %
                 (homo_alpha_ - (opt_.vmin + c.v)) %
                 ((homo_alpha_ + 1 + c.c) - (homo_alpha_ + 1)) %
                 (100.0 * c.weight)
          << flush;
    } else {
      XTP_LOG(Log::error, log_)
          << boost::format(
                 "           [beta ] HOMO-%1$-3d -> LUMO+%2$-3d  : "
                 "%3$+6.2f%%") %
                 (homo_beta_ - (opt_.vmin + c.v)) %
                 ((homo_beta_ + 1 + c.c) - (homo_beta_ + 1)) %
                 (100.0 * c.weight)
          << flush;
    }
  }
 
  if (static_cast<Index>(contributions.size()) > nprint) {
    XTP_LOG(Log::error, log_)
        << boost::format(
               "           ... %1$d more contributions above threshold") %
               (static_cast<Index>(contributions.size()) - nprint)
        << flush;
  }
}
 
void BSE_UKS::Analyze_excitons_uks(
    std::vector<QMFragment<BSE_Population>> fragments,
    const Orbitals& orb) const {
  (void)fragments;
 
  const tools::EigenSystem& es = orb.BSEUKS();
  const Eigen::VectorXd oscs =
      orb.Oscillatorstrengths(QMStateType(QMStateType::ExcitonUKS));
 
  if (orb.getTDAApprox()) {
    XTP_LOG(Log::error, log_)
        << "  ====== combined UKS TDA exciton energies (eV) ====== " << flush;
  } else {
    XTP_LOG(Log::error, log_)
        << "  ====== combined UKS full-BSE exciton energies (eV) ====== "
        << flush;
  }
 
  for (Index i = 0; i < std::min<Index>(opt_.nmax, es.eigenvalues().size());
       ++i) {
    XTP_LOG(Log::error, log_) << boost::format("  XU%-4d %+1.6f") % (i + 1) %
                                     (es.eigenvalues()(i) * tools::conv::hrt2ev)
                              << flush;
 
    const Eigen::Vector3d& trdip = orb.TransitionDipoles()[i];
    const double osc = (i < oscs.size()) ? oscs(i) : 0.0;
 
    XTP_LOG(Log::error, log_)
        << boost::format(
               "           TrDipole length gauge[e*bohr]  dx = %1$+1.4f dy = "
               "%2$+1.4f dz = %3$+1.4f |d|^2 = %4$+1.4f f = %5$+1.4f") %
               trdip[0] % trdip[1] % trdip[2] % (trdip.squaredNorm()) % osc
        << flush;
 
    Eigen::VectorXd weights = es.eigenvectors().col(i).cwiseAbs2();
    if (!orb.getTDAApprox()) {
      weights -= es.eigenvectors2().col(i).cwiseAbs2();
    }
 
    PrintWeightsUKS(weights);
  }
}
 
void BSE_UKS::Perturbative_DynamicalScreening(Orbitals& orb) {
 
  const tools::EigenSystem& BSECoefs = orb.BSEUKS();
  const Eigen::VectorXd& BSEenergies = BSECoefs.eigenvalues();
 
  const Eigen::VectorXd& RPAInputEnergiesAlpha = orb.RPAInputEnergiesAlpha();
  const Eigen::VectorXd& RPAInputEnergiesBeta = orb.RPAInputEnergiesBeta();
 
  // Static reference contribution of the direct kernel.
  SetupDirectInteractionOperator(RPAInputEnergiesAlpha, RPAInputEnergiesBeta,
                                 0.0);
 
  HdUKSOperator Hd_static(epsilon_0_inv_, Mmn_, Hqp_alpha_, Hqp_beta_);
  configureBSEOperator(Hd_static);
 
  ExpectationValues expectation_values =
      ExpectationValue_Operator(orb, Hd_static);
  Eigen::VectorXd Hd_static_contribution = expectation_values.direct_term;
 
  if (!orb.getTDAApprox()) {
    Hd2UKSOperator Hd2_static(epsilon_0_inv_, Mmn_, Hqp_alpha_, Hqp_beta_);
    configureBSEOperator(Hd2_static);
    expectation_values = ExpectationValue_Operator(orb, Hd2_static);
    Hd_static_contribution += expectation_values.cross_term;
  }
 
  Eigen::VectorXd BSEenergies_dynamic = BSEenergies;
 
  for (Index i_exc = 0; i_exc < BSEenergies.size(); ++i_exc) {
    XTP_LOG(Log::info, log_)
        << "Dynamical Screening UKS BSE, Excitation " << i_exc << " static "
        << BSEenergies(i_exc) << flush;
 
    for (Index iter = 0; iter < opt_.max_dyn_iter; ++iter) {
      const double old_energy = BSEenergies_dynamic(i_exc);
 
      SetupDirectInteractionOperator(RPAInputEnergiesAlpha,
                                     RPAInputEnergiesBeta, old_energy);
 
      HdUKSOperator Hd_dyn(epsilon_0_inv_, Mmn_, Hqp_alpha_, Hqp_beta_);
      configureBSEOperator(Hd_dyn);
 
      expectation_values = ExpectationValue_Operator_State(i_exc, orb, Hd_dyn);
      Eigen::VectorXd Hd_dynamic_contribution = expectation_values.direct_term;
 
      if (!orb.getTDAApprox()) {
        Hd2UKSOperator Hd2_dyn(epsilon_0_inv_, Mmn_, Hqp_alpha_, Hqp_beta_);
        configureBSEOperator(Hd2_dyn);
        expectation_values =
            ExpectationValue_Operator_State(i_exc, orb, Hd2_dyn);
        Hd_dynamic_contribution += expectation_values.cross_term;
      }
 
      BSEenergies_dynamic(i_exc) = BSEenergies(i_exc) +
                                   Hd_static_contribution(i_exc) -
                                   Hd_dynamic_contribution(0);
 
      XTP_LOG(Log::info, log_) << "Dynamical Screening UKS BSE, excitation "
                               << i_exc << " iteration " << iter << " dynamic "
                               << BSEenergies_dynamic(i_exc) << flush;
 
      if (std::abs(BSEenergies_dynamic(i_exc) - old_energy) <
          opt_.dyn_tolerance) {
        break;
      }
    }
  }
 
  orb.BSEUKS_dynamic() = BSEenergies_dynamic;
 
  const double hrt2ev = tools::conv::hrt2ev;
  const bool has_dipoles = orb.hasTransitionDipoles();
 
  Eigen::VectorXd oscs = Eigen::VectorXd::Zero(0);
  if (has_dipoles) {
    oscs = orb.Oscillatorstrengths(QMStateType(QMStateType::ExcitonUKS));
  }
 
  if (orb.getTDAApprox()) {
    XTP_LOG(Log::error, log_)
        << "  ====== combined UKS TDA exciton energies with perturbative "
           "dynamical screening (eV) ====== "
        << flush;
  } else {
    XTP_LOG(Log::error, log_)
        << "  ====== combined UKS full-BSE exciton energies with perturbative "
           "dynamical screening (eV) ====== "
        << flush;
  }
 
  const Index nprint = std::min<Index>(opt_.nmax, BSEenergies_dynamic.size());
 
  for (Index i = 0; i < nprint; ++i) {
    double osc = 0.0;
    if (has_dipoles && i < oscs.size()) {
      osc = oscs(i) * BSEenergies_dynamic(i) / BSEenergies(i);
    }
 
    XTP_LOG(Log::error, log_)
        << boost::format(
               "  XU(dynamic) = %1$4d Omega = %2$+1.12f eV  lambda = %3$+3.2f "
               "nm"
               " f = %4$+1.4f") %
               (i + 1) % (hrt2ev * BSEenergies_dynamic(i)) %
               (1240.0 / (hrt2ev * BSEenergies_dynamic(i))) % osc
        << flush;
  }
}
 
template tools::EigenSystem BSE_UKS::solve_hermitian<ExcitonUKSOperator_TDA>(
    ExcitonUKSOperator_TDA&) const;
 
template tools::EigenSystem BSE_UKS::Solve_nonhermitian_Davidson<
    ExcitonUKSOperator_TDA, ExcitonUKSOperator_BTDA_B>(
    ExcitonUKSOperator_TDA&, ExcitonUKSOperator_BTDA_B&) const;
 
template void BSE_UKS::configureBSEOperator<ExcitonUKSOperator_TDA>(
    ExcitonUKSOperator_TDA&) const;
template void BSE_UKS::configureBSEOperator<ExcitonUKSOperator_BTDA_B>(
    ExcitonUKSOperator_BTDA_B&) const;
template void BSE_UKS::configureBSEOperator<HdUKSOperator>(
    HdUKSOperator&) const;
template void BSE_UKS::configureBSEOperator<Hd2UKSOperator>(
    Hd2UKSOperator&) const;
 
template BSE_UKS::ExpectationValues
    BSE_UKS::ExpectationValue_Operator<ExcitonUKSOperator_TDA>(
        const Orbitals&, const ExcitonUKSOperator_TDA&) const;
template BSE_UKS::ExpectationValues
    BSE_UKS::ExpectationValue_Operator<ExcitonUKSOperator_BTDA_B>(
        const Orbitals&, const ExcitonUKSOperator_BTDA_B&) const;
template BSE_UKS::ExpectationValues BSE_UKS::ExpectationValue_Operator<
    HdUKSOperator>(const Orbitals&, const HdUKSOperator&) const;
template BSE_UKS::ExpectationValues BSE_UKS::ExpectationValue_Operator<
    Hd2UKSOperator>(const Orbitals&, const Hd2UKSOperator&) const;
 
template BSE_UKS::ExpectationValues
    BSE_UKS::ExpectationValue_Operator_State<ExcitonUKSOperator_TDA>(
        Index, const Orbitals&, const ExcitonUKSOperator_TDA&) const;
template BSE_UKS::ExpectationValues
    BSE_UKS::ExpectationValue_Operator_State<ExcitonUKSOperator_BTDA_B>(
        Index, const Orbitals&, const ExcitonUKSOperator_BTDA_B&) const;
template BSE_UKS::ExpectationValues BSE_UKS::ExpectationValue_Operator_State<
    HdUKSOperator>(Index, const Orbitals&, const HdUKSOperator&) const;
template BSE_UKS::ExpectationValues BSE_UKS::ExpectationValue_Operator_State<
    Hd2UKSOperator>(Index, const Orbitals&, const Hd2UKSOperator&) const;
 
}  // namespace xtp
}  // namespace votca