gw.h

Go to the documentation of this file.

/*
 *            Copyright 2009-2020 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_GW_H
#define VOTCA_XTP_GW_H
 
#include <unordered_map>
#include <unordered_set>
 
// Local VOTCA includes
#include "dfttimings.h"
#include "logger.h"
#include "orbitals.h"
#include "qp_solver_utils.h"
#include "rpa.h"
#include "sigma_base.h"
#include "threecenter.h"
 
namespace votca {
namespace xtp {
 
class GW {
 
  using EvalStage = qp_solver::EvalStage;
  using QPStats = qp_solver::Stats;
  using QPRootCandidate = qp_solver::RootCandidate;
  using QPWindowDiagnostics = qp_solver::WindowDiagnostics;
 
 public:
  GW(Logger& log, TCMatrix_gwbse& Mmn, const Eigen::MatrixXd& vxc,
     const Eigen::VectorXd& dft_energies)
      : log_(log),
        Mmn_(Mmn),
        vxc_(vxc),
        dft_energies_(dft_energies),
        rpa_(log, Mmn) {};
 
  struct options {
    Index homo;
    Index qpmin;
    Index qpmax;
    Index rpamin;
    Index rpamax;
    double eta;
    double g_sc_limit;
    Index g_sc_max_iterations;
    double gw_sc_limit;
    Index gw_sc_max_iterations;
    double shift = 0;
    double ScaHFX = 0.0;
    std::string sigma_integration;
    Index reset_3c;  // how often the 3c integrals in iterate should be
                     // rebuilt
    std::string qp_solver;
    double qp_solver_alpha = 0.75;
    // Legacy aliases from the original grid-search implementation.
    // They used to define both the total QP window and the scan resolution at
    // the same time. The new implementation keeps them only so that old XML
    // files and existing tests still work; configuration code maps them onto
    // the decoupled controls below.
    Index qp_grid_steps = 0;
    double qp_grid_spacing = 0.0;
 
    // Decoupled QP search controls:
    //  - qp_full_window_half_width defines the outer full search interval
    //  - qp_dense_spacing defines the dense scan used for robust sign changes
    //  - qp_adaptive_shell_width / count define the outward shell search
    // All values are normalized in configure(); negative values mean "unset"
    // and trigger either legacy mapping or robust defaults.
    double qp_full_window_half_width = -1.0;  // Ha; <=0 means "unset"
    double qp_dense_spacing = -1.0;           // Ha; <=0 means "unset"
    double qp_adaptive_shell_width = -1.0;    // Ha; <=0 means "unset"
    Index qp_adaptive_shell_count = 0;        // 0 means "use shell width"
 
    Index gw_mixing_order;          // mixing order
    double gw_mixing_alpha;         // mixing alpha, also linear mixing
    std::string quadrature_scheme;  // Kind of Gaussian-quadrature scheme to use
    Index order;   // only needed for complex integration sigma CDA
    double alpha;  // smooth tail in complex integration sigma CDA
    bool qp_restrict_search = true;
    // evGW: after the first iteration keep each level on the root branch of
    // the previous iteration (closest accepted root) instead of re-choosing
    // by weight
    bool qp_root_continuity = true;
    // evGW: when W is rebuilt (see ScreeningUpdate): "every" or "adaptive"
    std::string screening_update = "adaptive";
    // evGW: after the first iteration search each level's root only near
    // its previous one (falling back to the full search), verified by one
    // full search after convergence
    bool qp_root_tracking = false;
    double screening_update_ratio = 0.1;
    Index screening_update_max_inner = 10;
    double qp_zero_margin = 1e-6;
    double qp_virtual_min_energy = -0.1;
    std::string qp_root_finder = "brent";
    std::string qp_grid_search_mode = "adaptive_with_dense_fallback";
 
    // QSGW options
    // When true, run quasiparticle self-consistent GW instead of evGW.
    // The three-centre integrals are rotated at each iteration so that W and
    // Sigma are evaluated in the basis of the current QP wavefunctions.
    bool do_qsgw = false;
    Index qsgw_max_iterations = 50;
    double qsgw_sc_limit = 1e-5;  // Ha; convergence threshold on QP energies
 
    // Maximum allowed perturbative QP correction for virtual states to be
    // included in the QSGW self-consistency loop. Virtual states with
    // |e_QP - e_DFT| > qsgw_max_virt_correction are excluded from the QSGW
    // window and kept at their perturbative QP energies (DFT-MO wavefunctions).
    // This prevents pathological high-energy basis-set artefact states from
    // destabilising the QSGW rotation. Occupied states are never excluded.
    // Default: 0.5 Ha (13.6 eV, 1 Rydberg) -- corrections above this are
    // unlikely to be physically meaningful for typical molecular systems.
    // Set to a large value (e.g. 1e10) to disable the threshold entirely.
    double qsgw_max_virt_correction = 0.5;  // Ha
    bool qsgw_hartree_update = true;        // J follows the QP density
    // With a reaction field: the Hartree update is with u = v + v_reac, so
    // the environment's induced dipoles follow the QP density as well.
    // false: bare v (a control; inconsistent with the dressed W).
    bool qsgw_hartree_screened = true;
    // Static self-energy of QSGW: "eta" = the symmetrised
    // 1/2 Re[Sigma(e_p) + Sigma(e_q)] with denominators t/(t^2 + eta^2);
    // "srg" = the similarity-renormalization-group form of Marie and Loos,
    // Hermitian by construction, regularized by the flow parameter s (1/Ha^2)
    std::string qsgw_regularization = "eta";
    double qsgw_srg_s = 500.0;
    // RPA energies of levels outside the GW window (OutOfWindowShift)
    OutOfWindowShift out_of_window_shift;
  };
 
  void configure(const options& opt);

  Index ScreeningBuilds() const { return screening_builds_; }
  Index QPIterations() const { return qp_iterations_; }
  Index TrackedLocalSolves() const { return tracked_local_total_; }
  Eigen::VectorXd getGWAResults() const;

  void setReactionField(const Eigen::MatrixXd& R);
  void setCoreExchange(const Eigen::MatrixXd& X) { core_exchange_ = X; }
  const Eigen::MatrixXd& getSigmaReac() const { return Sigma_reac_; }
 
  // Calculates the diagonal elements up to self consistency
  void CalculateGWPerturbation();
 
  // Calculated offdiagonal elements as well
  void CalculateHQP();

  void CalculateQSGW();

  const Eigen::MatrixXd& getQSGWRotation() const { return qsgw_rotation_; }

  const Eigen::VectorXd& getQSGWSeedEnergies() const {
    return qsgw_seed_energies_;
  }

  void PrintQSGW_Energies(const std::string& seed_label,
                          const Eigen::VectorXd& seed_energies,
                          const Eigen::VectorXd& qsgw_energies) const;

  void PrintQSGW_Composition(double threshold = 0.01) const;
 
  Eigen::MatrixXd getHQP() const;
 
  // Diagonalize QP particle Hamiltonian
  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> DiagonalizeQPHamiltonian()
      const;

  void PlotSigma(std::string filename, Index steps, double spacing,
                 std::string states);

  struct sigma_dump_options {
    std::string filename = "sigma_dump.h5";
    double levels_emax_ev = 20.0;  // per-level spectra and imaginary axis
                                   // only up to this far above the DFT
                                   // LUMO (eV); <= 0: all QP levels
    Index imag_points = 24;        // Gauss-Legendre points on (0, inf)
    double imag_scale = 0.5;       // nu = s (1+x)/(1-x)
    double spectral_min = -3.0;    // pole histogram, absolute energies
    double spectral_max = 3.0;
    double spectral_binwidth = 0.02;
    Index plot_steps = 201;      // Sigma_c(omega) of all QP levels on a grid
    double plot_spacing = 0.01;  // around their DFT energies; 0 steps: off
  };

  void DumpSigma(const sigma_dump_options& dopt, const std::string& gw_mode,
                 const std::string& version);
  void setKeepScreeningData(bool keep) {
    keep_screening_data_ = keep;
    if (sigma_) {
      sigma_->setKeepScreeningData(keep);
    }
  }
 
  Eigen::VectorXd RPAInputEnergies() const {
    return rpa_.getRPAInputEnergies();
  }
 
 private:
  Index qptotal_;
  Eigen::MatrixXd SigmaCurves(const std::vector<Index>& levels, Index steps,
                              double spacing) const;
  bool keep_screening_data_ = false;
  // QP solution of the last solve and the energies in G it used (DumpSigma)
  Eigen::VectorXd qp_solution_;
  Eigen::VectorXd rpa_energies_solved_;
  // wall-clock accounting of the GW steps, reported at the end
  mutable DFTTimings timings_;
 
  Eigen::MatrixXd Sigma_x_;         // Fock exchange with the bare v, only
  Eigen::MatrixXd Sigma_c_;         // dynamic correlation
  Eigen::MatrixXd reaction_field_;  // R, fill-time aux frame; empty: off
  Eigen::MatrixXd dressing_;        // (1 + R)^(1/2), fill-time aux frame
  void DressForEnvironment();
  Eigen::MatrixXd Sigma_reac_;  // static correlation from v_reac (COH+SEX)
  void PrintReactionFieldShifts() const;
  Eigen::MatrixXd qsgw_rotation_;       // accumulated U: DFT MOs -> QSGW QP
                                        // wavefunctions
  Eigen::VectorXd qsgw_seed_energies_;  // evGW/G0W0 energies used as QSGW seed
  // Merged QSGW+seed energies for the full QP window. Set by CalculateQSGW()
  // when the virtual window is trimmed. getGWAResults() returns this when
  // non-empty, bypassing the sigma-matrix recomputation which would be wrong
  // for the excluded levels.
  Eigen::VectorXd qsgw_final_energies_;
 
  options opt_;
 
  std::unique_ptr<Sigma_base> sigma_ = nullptr;
  Logger& log_;
  TCMatrix_gwbse& Mmn_;
  const Eigen::MatrixXd& vxc_;
  const Eigen::VectorXd& dft_energies_;
 
  Eigen::MatrixXd core_exchange_;  // see setCoreExchange
  Index screening_builds_ = 0;
  Index qp_iterations_ = 0;
  Index gw_sc_iteration_;
 
  RPA rpa_;
  // small class which calculates f(w) with and df/dw(w)
  // f=Sigma_c(w)+offset-w
  // offset= e_dft+Sigma_x-Vxc
  class QPFunc {
   public:
    QPFunc(Index gw_level, const Sigma_base& sigma, double offset)
        : gw_level_(gw_level), offset_(offset), sigma_c_func_(sigma) {}
 
    std::pair<double, double> operator()(double frequency) const {
      std::pair<double, double> result;
      result.first = value(frequency, EvalStage::Other);
      result.second = deriv(frequency);
      return result;
    }
 
    double sigma(double frequency, EvalStage stage = EvalStage::Other) const {
      const std::uint64_t key = FrequencyKey(frequency);
 
      auto insert_result = seen_frequencies_.insert(key);
      if (!insert_result.second) {
        ++stats_.sigma_repeat_calls;
      } else {
        ++stats_.sigma_unique_frequencies;
      }
 
      CountSigmaStage(stage);
      auto cached = cache_.find(key);
      if (cached != cache_.end()) {
        return cached->second;
      }
      return sigma_c_func_.CalcCorrelationDiagElement(gw_level_, frequency);
    }

    void Prefetch(const std::vector<double>& frequencies) const {
      std::vector<double> todo;
      for (double f : frequencies) {
        if (cache_.find(FrequencyKey(f)) == cache_.end()) {
          todo.push_back(f);
        }
      }
      if (todo.empty()) {
        return;
      }
      const Eigen::VectorXd values = sigma_c_func_.CalcCorrelationDiagElements(
          gw_level_,
          Eigen::Map<const Eigen::VectorXd>(todo.data(), Index(todo.size())));
      for (std::size_t i = 0; i < todo.size(); ++i) {
        cache_[FrequencyKey(todo[i])] = values(Index(i));
      }
    }
 
    double value(double frequency, EvalStage stage = EvalStage::Other) const {
      return sigma(frequency, stage) + offset_ - frequency;
    }
 
    double deriv(double frequency) const {
      ++stats_.deriv_calls;
      return sigma_c_func_.CalcCorrelationDiagElementDerivative(gw_level_,
                                                                frequency) -
             1.0;
    }
 
    const QPStats& GetStats() const { return stats_; }
 
   private:
    static std::uint64_t FrequencyKey(double x) {
      std::uint64_t key = 0;
      static_assert(sizeof(double) == sizeof(std::uint64_t),
                    "Unexpected double size");
      std::memcpy(&key, &x, sizeof(double));
      return key;
    }
 
    void CountSigmaStage(EvalStage stage) const {
      switch (stage) {
        case EvalStage::Scan:
          ++stats_.sigma_scan_calls;
          break;
        case EvalStage::Refine:
          ++stats_.sigma_refine_calls;
          break;
        case EvalStage::Derivative:
          ++stats_.sigma_derivative_calls;
          break;
        case EvalStage::Other:
        default:
          ++stats_.sigma_other_calls;
          break;
      }
    }
 
    Index gw_level_;
    double offset_;
    const Sigma_base& sigma_c_func_;
 
    mutable std::unordered_set<std::uint64_t> seen_frequencies_;
    mutable std::unordered_map<std::uint64_t, double> cache_;
    mutable QPStats stats_;
  };
 
  double CalcHomoLumoShift(Eigen::VectorXd frequencies) const;
  Eigen::VectorXd ScissorShift_DFTlevel(
      const Eigen::VectorXd& dft_energies) const;
  void PrintQP_Energies(const Eigen::VectorXd& qp_diag_energies) const;
  void PrintGWA_Energies() const;
 
  Eigen::VectorXd SolveQP(const Eigen::VectorXd& frequencies) const;
  // HOMO-LUMO midpoint of the current QP guesses, set by SolveQP before
  // its parallel loop and only read inside it (SolveQP_Grid).
  mutable double qp_midgap_ = 0.0;
  // per level of the last SolveQP: whether two accepted roots carry
  // comparable weight (see qp_solver::CheckCompetingRoots)
  mutable std::vector<qp_solver::CompetingRoots> competing_roots_;
  // root tracking (opt_.qp_root_tracking): the roots of the last QP solve,
  // their weights and their last change; active from the second iteration
  // on until the verifying full search
  mutable Eigen::VectorXd tracked_root_;
  mutable Eigen::VectorXd tracked_Z_;
  mutable Eigen::VectorXd tracked_change_;
  mutable bool tracking_active_ = false;
  mutable Index tracked_local_ = 0;  // levels solved locally in the last solve
  mutable Index tracked_local_total_ = 0;
  boost::optional<qp_solver::RootCandidate> SolveQP_Tracked(
      double intercept0, double frequency0, Index gw_level,
      QPStats* stats) const;
  void RecordCompetingRoots(Index gw_level,
                            const std::vector<QPRootCandidate>& roots,
                            double chosen) const {
    if (gw_level < Index(competing_roots_.size())) {
      competing_roots_[std::size_t(gw_level)] =
          qp_solver::CheckCompetingRoots(roots, chosen);
    }
  }
  // logs the levels with competing roots of the last SolveQP
  void PrintCompetingRoots(Log::Level level) const;
  boost::optional<double> SolveQP_Grid(double intercept0, double frequency0,
                                       Index gw_level,
                                       QPStats* stats = nullptr) const;
 
  boost::optional<double> SolveQP_Grid_Windowed(
      double intercept0, double frequency0, Index gw_level, double left_limit,
      double right_limit, bool allow_rejected_return = true,
      QPStats* stats = nullptr) const;
 
  boost::optional<double> SolveQP_Grid_Windowed_Adaptive(
      double intercept0, double frequency0, Index gw_level, double left_limit,
      double right_limit, bool allow_rejected_return = true,
      QPStats* stats = nullptr) const;
 
  boost::optional<double> SolveQP_Grid_Windowed_Dense(
      double intercept0, double frequency0, Index gw_level, double left_limit,
      double right_limit, bool allow_rejected_return = true,
      QPStats* stats = nullptr) const;
 
  boost::optional<double> SolveQP_FixedPoint(double intercept0,
                                             double frequency0, Index gw_level,
                                             QPStats* stats = nullptr) const;
  boost::optional<double> SolveQP_Linearisation(double intercept0,
                                                double frequency0,
                                                Index gw_level,
                                                QPStats* stats = nullptr) const;
  bool Converged(const Eigen::VectorXd& e1, const Eigen::VectorXd& e2,
                 double epsilon) const;
 
  boost::optional<QPRootCandidate> RefineQPInterval(
      double lowerbound, double f_lowerbound, double upperbound,
      double f_upperbound, const QPFunc& f, double reference) const;
};
}  // namespace xtp
}  // namespace votca
 
#endif  // VOTCA_XTP_GW_H