rpa_uks.h

Go to the documentation of this file.

/*
 *            Copyright 2009-2021 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.
 *
 */
 
#pragma once
#ifndef VOTCA_XTP_RPA_UKS_H
#define VOTCA_XTP_RPA_UKS_H
 
#include <complex>
#include <vector>
 
#include "eigen.h"
#include "logger.h"
#include "out_of_window_shift.h"
#include "screening_kernels.h"
#include "votca/xtp/threecenter.h"
 
namespace votca {
namespace xtp {

class RPA_UKS {
 public:
  RPA_UKS(Logger& log, const TCMatrix_gwbse_spin& Mmn) : log_(log), Mmn_(Mmn) {}

  void setOutOfWindowShift(const OutOfWindowShift& shift) {
    out_of_window_ = shift;
  }
 
  void configure(Index homo_alpha, Index homo_beta, Index rpamin,
                 Index rpamax) {
    homo_alpha_ = homo_alpha;
    homo_beta_ = homo_beta;
    rpamin_ = rpamin;
    rpamax_ = rpamax;
    InvalidateH2pCache();
  }

  double getEta() const { return eta_; }

  void setTimings(DFTTimings* timings) { timings_ = timings; }

  Eigen::MatrixXd calculate_epsilon_i(double frequency) const {
    return calculate_epsilon<true>(frequency);
  }

  Eigen::MatrixXd calculate_epsilon_r(double frequency) const {
    return calculate_epsilon<false>(frequency);
  }

  Eigen::MatrixXd calculate_epsilon_r(std::complex<double> frequency) const;

  const Eigen::VectorXd& getRPAInputEnergiesAlpha() const {
    return energies_alpha_;
  }

  const Eigen::VectorXd& getRPAInputEnergiesBeta() const {
    return energies_beta_;
  }

  void setRPAInputEnergies(const Eigen::VectorXd& rpaenergies_alpha,
                           const Eigen::VectorXd& rpaenergies_beta) {
    energies_alpha_ = rpaenergies_alpha;
    energies_beta_ = rpaenergies_beta;
    InvalidateH2pCache();
  }

  void UpdateRPAInputEnergies(const Eigen::VectorXd& dftenergies_alpha,
                              const Eigen::VectorXd& dftenergies_beta,
                              const Eigen::VectorXd& gwaenergies_alpha,
                              const Eigen::VectorXd& gwaenergies_beta,
                              Index qpmin);


  struct rpa_eigensolution {
    Eigen::VectorXd omega;
    Eigen::MatrixXd XpY;
    double ERPA_correlation;
  };

  const rpa_eigensolution& Diagonalize_H2p() const;
 
  Index getHomoAlpha() const { return homo_alpha_; }
  Index getHomoBeta() const { return homo_beta_; }
 
  void GetCachedScreeningModes(
      const Eigen::VectorXd*& omegas,
      const std::vector<Eigen::VectorXd>*& modes) const;
 
 private:
  template <bool imag>
  Eigen::MatrixXd calculate_epsilon(double frequency) const;
 
  void InvalidateH2pCache();

  Eigen::VectorXd Calculate_H2p_AmB() const;

  Eigen::MatrixXd Calculate_H2p_ApB() const;

  SymmetricEigenSystem Diagonalize_H2p_C(const Eigen::MatrixXd& C) const;

  void ShiftUncorrectedEnergies(Eigen::VectorXd& energies,
                                const Eigen::VectorXd& dftenergies, Index homo,
                                Index qpmin, Index gwsize);
 
  OutOfWindowShift out_of_window_;
 
  // Spin-resolved HOMO indices in the full orbital numbering
  Index homo_alpha_;
  Index homo_beta_;
 
  // Shared orbital window used to build the RPA quantities
  Index rpamin_;
  Index rpamax_;
 
  // Small broadening parameter used in real-frequency expressions
  const double eta_ = 0.0001;
 
  // Current alpha/beta energies used in the response over [rpamin_, rpamax_]
  Eigen::VectorXd energies_alpha_;
  Eigen::VectorXd energies_beta_;
 
  Logger& log_;
 
  // Spin-resolved three-center integrals:
  // Mmn_.alpha for alpha MOs, Mmn_.beta for beta MOs
  const TCMatrix_gwbse_spin& Mmn_;
  DFTTimings* timings_ = nullptr;
 
  // weights(beta, v, w): weights of the virtual rows of hole slice v
  using SpinWeightsFn = std::function<void(bool, Index, Eigen::VectorXd&)>;
  // sum over both spins and hole slices v of M_v^T diag(w_v) M_v
  Eigen::MatrixXd ResponseSum(const SpinWeightsFn& weights) const;
 
  mutable bool h2p_cached_ = false;
  mutable rpa_eigensolution h2p_solution_cache_;
 
  // Cached screening modes (shared across alpha/beta exact GW)
  mutable bool screening_cached_ = false;
  mutable Eigen::VectorXd screening_omegas_;
  mutable std::vector<Eigen::VectorXd> screening_modes_;
  void BuildCachedScreeningModes() const;
};
 
}  // namespace xtp
}  // namespace votca
 
#endif