imfreq.hpp

Go to the documentation of this file.

// Copyright (c) 2013-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
// Copyright (c) 2013-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: Thomas Ayral, Philipp Dumitrescu, Olivier Parcollet, Nils Wentzell

 
#pragma once
 
#include "./matsubara_freq.hpp"
#include "./mesh_iterator.hpp"
#include "./tail_fitter.hpp"
#include "./utils.hpp"
#include "../utility/macros.hpp"
 
#include <fmt/format.h>
#include <h5/h5.hpp>
#include <nda/nda.hpp>
 
#include <complex>
#include <cstdint>
#include <iostream>
#include <numbers>
#include <string>
#include <tuple>
 
namespace triqs::mesh {

  struct energy_t {
    double value = 0;

    explicit operator double() const { return value; }
  };


  class C2PY_RENAME(MeshImFreq) imfreq : public tail_fitter_handle {
    public:
    using value_t = matsubara_freq;

    using index_t = long;

    using data_index_t = long;

    class C2PY_IGNORE mesh_point_t : public matsubara_freq {
      public:
      using mesh_t = imfreq;

      mesh_point_t() = default;

      mesh_point_t(double b, statistic_enum stat, index_t n_idx, long d, uint64_t mhash)
         : matsubara_freq(n_idx, b, stat), data_index_(d), mesh_hash_(mhash) {}

      [[nodiscard]] long index() const { return n; }

      [[nodiscard]] long data_index() const { return data_index_; }

      [[nodiscard]] matsubara_freq const &value() const { return *this; }

      [[nodiscard]] uint64_t mesh_hash() const noexcept { return mesh_hash_; }
 
      private:
      long data_index_    = 0;
      uint64_t mesh_hash_ = 0;
    };

    enum class option { all_frequencies, positive_frequencies_only };

    imfreq() = default;

    C2PY_DEPRECATED_PARAMETER_NAME(S : statistic, n_max : n_iw)
    imfreq(double beta, statistic_enum statistic, long n_iw = 1025, option opt = option::all_frequencies)
       : beta_(beta), stat_(statistic), N_iw_(n_iw), opt_(opt), mesh_hash_(hash(beta, statistic, n_iw, opt)) {
      EXPECTS(beta_ > 0);
      EXPECTS(N_iw_ >= 0);
      if (opt == option::positive_frequencies_only) first_index_ = 0;
    }

    C2PY_IGNORE imfreq(double beta, statistic_enum statistic, energy_t w_max, option opt = option::all_frequencies)
       : imfreq(beta, statistic, static_cast<long>((w_max.value * beta / std::numbers::pi - (statistic == Fermion ? 1 : 0)) / 2) + 1, opt) {}

    bool operator==(imfreq const &m) const { return (std::tie(beta_, stat_, N_iw_, opt_) == std::tie(m.beta_, m.stat_, m.N_iw_, m.opt_)); }

    [[nodiscard]] bool is_index_valid(index_t n) const { return first_index() <= n and n <= last_index(); }

    [[nodiscard]] data_index_t to_data_index(index_t n) const noexcept {
      EXPECTS(is_index_valid(n));
      return n - first_index();
    }

    [[nodiscard]] data_index_t to_data_index(matsubara_freq const &iw) const noexcept {
      EXPECTS(beta_ == iw.beta and stat_ == iw.statistic);
      return to_data_index(iw.n);
    }

    [[nodiscard]] C2PY_IGNORE data_index_t to_data_index(closest_mesh_point_t<value_t> const &cmp) const {
      EXPECTS(beta_ == cmp.value.beta and stat_ == cmp.value.statistic);
      return to_data_index(to_index(cmp));
    }

    [[nodiscard]] index_t to_index(data_index_t d) const {
      EXPECTS(0 <= d and d < size());
      return d + first_index();
    }

    [[nodiscard]] C2PY_IGNORE index_t to_index(closest_mesh_point_t<value_t> const &cmp) const {
      EXPECTS(is_index_valid(cmp.value.n));
      return cmp.value.n;
    }

    [[nodiscard]] mesh_point_t operator[](long d) const { return {beta_, stat_, to_index(d), d, mesh_hash_}; }

    [[nodiscard]] C2PY_IGNORE mesh_point_t operator[](closest_mesh_point_t<value_t> const &cmp) const { return (*this)[to_data_index(cmp)]; }

    [[nodiscard]] mesh_point_t operator()(long n) const { return {beta_, stat_, n, to_data_index(n), mesh_hash_}; }

    [[nodiscard]] matsubara_freq to_value(index_t n) const {
      EXPECTS(is_index_valid(n));
      return {n, beta_, stat_};
    }

    imfreq get_positive_freq() const { return {beta_, stat_, N_iw_, option::positive_frequencies_only}; }

    std::complex<double> w_max() const { return to_value(last_index_); }

    [[nodiscard]] C2PY_PROPERTY_GET(beta) double beta() const noexcept { return beta_; }

    [[nodiscard]] C2PY_PROPERTY_GET(statistic) statistic_enum statistic() const noexcept { return stat_; }

    [[nodiscard]] C2PY_PROPERTY_GET(n_iw) long n_iw() const noexcept { return N_iw_; }

    [[nodiscard]] C2PY_PROPERTY_GET(mesh_hash) uint64_t mesh_hash() const noexcept { return mesh_hash_; }

    [[nodiscard]] long size() const noexcept { return last_index_ - first_index_ + 1; }

    [[nodiscard]] long first_index() const { return first_index_; }

    [[nodiscard]] long last_index() const { return last_index_; }

    [[nodiscard]] bool positive_only() const { return opt_ == option::positive_frequencies_only; }

    [[nodiscard]] auto begin() const { return mesh_iterator<imfreq>{.mesh_ptr = this, .data_index = 0}; }

    [[nodiscard]] auto cbegin() const { return begin(); }

    [[nodiscard]] auto end() const { return mesh_iterator<imfreq>{.mesh_ptr = this, .data_index = size()}; }

    [[nodiscard]] auto cend() const { return end(); }

    friend std::ostream &operator<<(std::ostream &sout, imfreq const &m) {
      auto stat_cstr = (m.stat_ == Boson ? "Boson" : "Fermion");
      return sout << fmt::format("Imaginary frequency mesh with beta = {}, statistics = {}, N_iw = {}, positive_only = {}", m.beta_, stat_cstr,
                                 m.N_iw_, m.positive_only());
    }

    void serialize(auto &ar) const { ar & beta_ & stat_ & N_iw_ & opt_ & last_index_ & first_index_ & mesh_hash_; }

    void deserialize(auto &ar) { ar & beta_ & stat_ & N_iw_ & opt_ & last_index_ & first_index_ & mesh_hash_; }

    [[nodiscard]] static std::string hdf5_format() { return "MeshImFreq"; }

    friend void h5_write(h5::group g, std::string name, imfreq const &m) {
      h5::group gr = g.create_group(name);
      h5::write_hdf5_format(gr, m); // NOLINT (downcasting to base class)
      h5::write(gr, "beta", m.beta_);
      h5::write(gr, "statistic", (m.stat_ == Fermion ? "F" : "B"));
      h5::write(gr, "size", m.size());
      h5::write(gr, "positive_freq_only", (m.positive_only() ? 1 : 0));
    }

    friend void h5_read(h5::group g, std::string name, imfreq &m) {
      h5::group gr = g.open_group(name);
      h5::assert_hdf5_format(gr, m, true); // NOLINT (downcasting to base class)
 
      // enum option: positive_frequencies_only = 1, all_frequencies = 0
      int pos_freq = 0;
      if (gr.has_key("positive_freq_only")) h5::read(gr, "positive_freq_only", pos_freq);
      if (gr.has_key("start_at_0")) h5::read(gr, "start_at_0", pos_freq); // backward compatibility
      auto opt = (pos_freq == 1 ? option::positive_frequencies_only : option::all_frequencies);
 
      // size N and the number of positive frequencies N_iw
      long N    = h5::read<long>(gr, "size");
      long N_iw = (pos_freq ? N : (N + 1) / 2);
 
      // Matsubara frequency domain for backward compatibility
      if (gr.has_key("domain")) { gr = gr.open_group("domain"); }
 
      // beta and statistics
      auto beta = h5::read<double>(gr, "beta");
      auto stat = (h5::read<std::string>(gr, "statistic") == "F" ? Fermion : Boson);
 
      m = imfreq{beta, stat, N_iw, opt};
    }

    C2PY_IGNORE bool eval_to_zero(index_t n) const { return !is_index_valid(n); }

    C2PY_IGNORE bool eval_to_zero(matsubara_freq iw) const { return eval_to_zero(iw.n); }

    C2PY_IGNORE bool eval_to_zero(mesh_point_t mp) const { return eval_to_zero(mp.value()); }
 
    private:
    double beta_         = 1.0;
    statistic_enum stat_ = Fermion;
    long N_iw_           = 0;
    option opt_          = option::all_frequencies;
    long last_index_     = N_iw_ - 1;
    long first_index_    = -(last_index_ + ((stat_ == Fermion) ? 1 : 0));
    uint64_t mesh_hash_  = 0;
  };

  auto evaluate([[maybe_unused]] imfreq const &m, auto const &f, matsubara_freq const &iw) {
    EXPECTS(m.beta() == iw.beta and m.statistic() == iw.statistic);
    return f(iw.n);
  }

 
} // namespace triqs::mesh