evaluator.hpp

Go to the documentation of this file.

// Copyright (c) 2016-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
// Copyright (c) 2016-2018 Centre national de la recherche scientifique (CNRS)
// Copyright (c) 2018-2023 Simons Foundation
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You may obtain a copy of the License at
//     https://www.gnu.org/licenses/gpl-3.0.txt
//
// Authors: Michel Ferrero, Olivier Parcollet, Nils Wentzell

 
#pragma once
 
#include "./functions/functions2.hpp"
#include "./gf/defs.hpp"
#include "./gf/gf.hpp"
 
#include "../mesh/imfreq.hpp"
 
namespace triqs::gfs {
 
  // evaluator by default forwards everything to evaluate
  // specialize e.g. for tail where gf data is used


  template <Mesh M> struct gf_evaluator {

    template <typename G, typename... XS>
      requires(is_gf_v<G>)
    auto operator()(G const &g, XS &&...xs) const {
      auto l    = [&g](auto &&...ys) -> decltype(auto) { return g.operator[](ys...); };
      using r_t = std::decay_t<decltype(make_regular(evaluate(g.mesh(), l, std::forward<XS>(xs)...)))>;
      if constexpr (nda::Array<r_t> or nda::is_scalar_v<r_t>) {
        // Return zero if any mesh evaluates to zero
        if (detail::eval_to_zero(g.mesh(), xs...)) { return r_t{nda::zeros<typename G::target_t::scalar_t>(g.target_shape())}; }
      }
      return make_regular(evaluate(g.mesh(), l, std::forward<XS>(xs)...));
    }
  };
 
  // Evaluator specialization for Matsubara-frequency Green's functions (internal; the primary gf_evaluator is the
  // documented API). Off the stored grid it uses the analytic structure: on a positive-only mesh negative
  // frequencies are reconstructed by conjugation, and outside the grid the high-frequency tail is fitted and summed.
  template <> struct gf_evaluator<mesh::imfreq> {
 
    // Evaluate at a Matsubara frequency: stored value on the grid, conjugate on a positive-only mesh, else tail sum.
    template <typename G> auto operator()(G const &g, matsubara_freq const &f) const {
 
      using r_t = std::decay_t<decltype(make_regular(g[0]))>;
 
      if (g.mesh().is_index_valid(f.n)) return r_t{g[f.n]};
      if (g.mesh().positive_only()) {
        int sh = (g.mesh().statistic() == Fermion ? 1 : 0);
        if (g.mesh().is_index_valid(-f.n - sh)) return r_t{conj(g[-f.n - sh])};
        TRIQS_RUNTIME_ERROR << " ERROR: Cannot evaluate Green function with positive only mesh outside grid ";
      }
 
      auto [tail, err] = fit_tail_no_normalize(g);
      dcomplex x       = std::abs(g.mesh().w_max()) / f;
      auto res         = r_t{nda::zeros<dcomplex>(g.target_shape())}; // a new array
 
      dcomplex z = 1.0;
      for (int n : range(tail.extent(0))) {
        res += tail(n, ellipsis()) * z;
        z = z * x;
      }
 
      return res;
    }
 
    // Integer overload: convert n to a matsubara_freq and delegate to the matsubara_freq overload.
    template <typename G> decltype(auto) operator()(G const &g, int n) const { return g(matsubara_freq(n, g.mesh().beta(), g.mesh().statistic())); }
  };

 
} // namespace triqs::gfs