gf: update gf imfreq for g( om - nu) case

DRAFT : to be tested further...

- update gf<imfreq>
- write a specific mesh for matsubara frequencies
- now the cast series is :
  mesh_pt --> matsubara_freq --> complex<double>
- matsubara_freq is just the matsubara frequency
- arithmetic of the mesh_pt casted to matsubara_freq
- arithmetic of matsubara_freq is casted to complex, except + and -,
  which are kept as matsubara_freq.
- evaluator now accept : int, mesh_pt, and matsubara_freq
  for matsubara_freq : for negative omega, use conjugation
  for omega outside windows, evaluate the tail on omega.
- as a result : g( om - nu) where om, nu are 2 meshes points,
  is the extrapolation outside the grid if necessary.
- updated tests
- added evaluation for tail.

test/triqs/gfs/gfv2.cpp

@@ -36,6 +36,7 @@ int main() {
   TEST( G( 0) ) ;
 
   triqs::clef::placeholder<0> om_;
+  triqs::clef::placeholder<1> nu_;
 
   TEST( G(om_) ) ;
   TEST( eval(G(om_), om_=0) ) ;
@@ -121,6 +122,25 @@ int main() {
   H5::H5File file("ess_gf.h5", H5F_ACC_TRUNC );
   h5_write(file, "g", G);
 
+  //
+  {
+   auto G0w = gf<imfreq, scalar_valued>{{beta, Fermion, 100}};
+   auto D0w = gf<imfreq, scalar_valued>{{beta, Boson, 100}};
+   auto Sigma_w = gf<imfreq, scalar_valued>{{beta, Fermion, 100}};
+   G0w(om_) << 1 / (om_ - 3);
+
+   for (auto const& nu : D0w.mesh()) Sigma_w(om_) << 2 * G0w(om_ - nu);
+  }
+
+  //
+  {
+   auto G0w = gf<imfreq, matrix_valued>{{beta, Fermion, 100}, {1,1}};
+   auto D0w = gf<imfreq, matrix_valued>{{beta, Boson, 100}, {1,1}};
+   auto Sigma_w = gf<imfreq, matrix_valued>{{beta, Fermion, 100}, {1,1}} ;
+   G0w(om_) << 1 / (om_ - 3);
+
+   for (auto const& nu : D0w.mesh()) Sigma_w(om_) << 2 * G0w(om_ - nu);
+  }
  }
  TRIQS_CATCH_AND_ABORT;
 

triqs/gfs/domains/matsubara.hpp

@@ -21,13 +21,66 @@
 #ifndef TRIQS_GF_DOMAIN_MATSUBARA_H
 #define TRIQS_GF_DOMAIN_MATSUBARA_H
 #include "../tools.hpp"
+#include "../meshes/mesh_tools.hpp"
 
-namespace triqs { namespace gfs { 
+namespace triqs {
+namespace gfs {
 
+ // --------------- One matsubara frequency, with its arithmetics -------------------------
+ // all operations are done by casting to complex, except addition and substraction of matsubara_freq
+
+ struct matsubara_freq :  public arith_ops_by_cast<matsubara_freq, std::complex<double>>  {
+  int n;
+  double beta;
+  statistic_enum statistic;
+  matsubara_freq(int const &n_, double beta_, statistic_enum stat_) : n(n_), beta(beta_), statistic(stat_) {}
+  using cast_t = std::complex<double>;
+  operator cast_t() const {
+   return {0, M_PI * (2 * n + (statistic == Fermion ? 1 : 0)) / beta};
+  }
+ };
+
+  inline matsubara_freq operator+(matsubara_freq const &x, matsubara_freq const &y) {
+  return {x.n + y.n + (x.statistic & y.statistic), x.beta, ((x.statistic ^ y.statistic) == 1 ? Fermion : Boson)};
+ }
+
+ inline matsubara_freq operator-(matsubara_freq const &x, matsubara_freq const &y) {
+  // std::cout  << x.n - y.n - (~x.statistic & y.statistic)<< std::endl;
+  return {x.n - y.n - (~x.statistic & y.statistic), x.beta, ((x.statistic ^ y.statistic) == 1 ? Fermion : Boson)};
+ }
+
+ /*
+#define IMPL_OP(OP)                                                                                                              \
+ template <typename Y> auto operator OP(matsubara_freq const &x, Y &&y)->decltype(std::complex<double>{}) {                      \
+  return std::complex<double>(x) OP std::forward<Y>(y);                                                                          \
+ }                                                                                                                               \
+ template <typename Y>                                                                                                           \
+ auto operator OP(Y &&y, matsubara_freq const &x)->decltype(std::forward<Y>(y) OP std::complex<double>{}) {                      \
+  return std::forward<Y>(y) OP std::complex<double>(x);                                                                          \
+ }
+ IMPL_OP(+);
+ IMPL_OP(-);
+ IMPL_OP(*);
+ IMPL_OP(/ );
+#undef IMPL_OP
+*/
+
+/*
+ * template <typename Y> auto operator OP(matsubara_freq const &x, Y &&y)->decltype(std::complex<double>{}) {                        \
+  return std::complex<double>(x) OP std::forward<Y>(y);                                                                          \
+ }                                                                                                                               \
+ template <typename Y>                                                                                                           \
+ auto operator OP(Y &&y, matsubara_freq const &x)                                                                                \
+     ->TYPE_DISABLE_IF(decltype(std::forward<Y>(y) OP std::complex<double>{}),                                                   \
+                       std::is_same<typename std::remove_cv<typename std::remove_reference<Y>::type>::type, matsubara_freq>) {   \
+  return std::forward<Y>(y) OP std::complex<double>(x);                                                                          \
+ }
+*/
+
+ //---------------------------------------------------------------------------------------------------------
  /// The domain
- template<bool IsComplex>
+ template <bool IsFreq> struct matsubara_domain {
- struct matsubara_domain {
+  using point_t = typename std::conditional<IsFreq, std::complex<double>, double>::type;
-  typedef typename mpl::if_c<IsComplex, std::complex<double>, double>::type point_t;
   double beta;
   statistic_enum statistic;
   matsubara_domain() = default;
@@ -35,7 +88,7 @@ namespace triqs { namespace gfs {
    if (beta < 0) TRIQS_RUNTIME_ERROR << "Matsubara domain construction :  beta <0 : beta =" << beta << "\n";
   }
   matsubara_domain(matsubara_domain const &) = default;
-  matsubara_domain(matsubara_domain<!IsComplex> const &x): beta(x.beta), statistic(x.statistic) {}
+  matsubara_domain(matsubara_domain<!IsFreq> const &x) : beta(x.beta), statistic(x.statistic) {}
   bool operator==(matsubara_domain const &D) const { return ((std::abs(beta - D.beta) < 1.e-15) && (statistic == D.statistic)); }
 
   /// Write into HDF5
@@ -48,7 +101,8 @@ namespace triqs { namespace gfs {
   /// Read from HDF5
   friend void h5_read(h5::group fg, std::string subgroup_name, matsubara_domain &d) {
    h5::group gr = fg.open_group(subgroup_name);
-   double beta; std::string statistic;
+   double beta;
+   std::string statistic;
    h5_read(gr, "beta", beta);
    h5_read(gr, "statistic", statistic);
    d = matsubara_domain(beta, (statistic == "F" ? Fermion : Boson));
@@ -56,14 +110,12 @@ namespace triqs { namespace gfs {
 
   //  BOOST Serialization
   friend class boost::serialization::access;
-  template<class Archive>
+  template <class Archive> void serialize(Archive &ar, const unsigned int version) {
-   void serialize(Archive & ar, const unsigned int version) {
    ar &boost::serialization::make_nvp("beta", beta);
    ar &boost::serialization::make_nvp("statistic", statistic);
   }
-
  };
-
+}
-}}
+}
 #endif
 

triqs/gfs/imfreq.hpp

@@ -24,49 +24,98 @@
 #include "./gf.hpp"
 #include "./local/tail.hpp"
 #include "./local/no_tail.hpp"
-#include "./domains/matsubara.hpp"
+#include "./meshes/matsubara_freq.hpp"
-#include "./meshes/linear.hpp"
 #include "./evaluators.hpp"
-namespace triqs { namespace gfs {
+namespace triqs {
+namespace gfs {
 
  struct imfreq {};
 
- template<typename Opt> struct gf_mesh<imfreq,Opt> : linear_mesh<matsubara_domain<true>> {
+ template <typename Opt> struct gf_mesh<imfreq, Opt> : matsubara_freq_mesh {
-  typedef  linear_mesh<matsubara_domain<true>> B;
+  using B = matsubara_freq_mesh;
   static double m1(double beta) { return std::acos(-1) / beta; }
   gf_mesh() = default;
   gf_mesh(B const &x) : B(x) {} // enables also construction from another Opt
-  gf_mesh (typename B::domain_t const & d, int Nmax = 1025) : 
+  gf_mesh(typename B::domain_t const &d, int Nmax = 1025) : B(d, Nmax, true) {}
-   B(d, d.statistic==Fermion?m1(d.beta):0, d.statistic==Fermion?(2*Nmax+1)*m1(d.beta): 2*Nmax*m1(d.beta), Nmax, without_last){}
   gf_mesh(double beta, statistic_enum S, int Nmax = 1025) : gf_mesh({beta, S}, Nmax) {}
  };
 
  namespace gfs_implementation {
 
   // singularity
-  template<> struct singularity<imfreq,matrix_valued,void>  { typedef local::tail type;};
+  template <> struct singularity<imfreq, matrix_valued, void> {
-  template<> struct singularity<imfreq,scalar_valued,void>  { typedef local::tail type;};
+   typedef local::tail type;
+  };
+  template <> struct singularity<imfreq, scalar_valued, void> {
+   typedef local::tail type;
+  };
 
   // h5 name
-  template<typename Opt> struct h5_name<imfreq,matrix_valued,Opt>      { static std::string invoke(){ return "ImFreq";}};
+  template <typename Opt> struct h5_name<imfreq, matrix_valued, Opt> {
+   static std::string invoke() { return "ImFreq"; }
+  };
 
   /// ---------------------------  evaluator ---------------------------------
 
- template<> struct evaluator_fnt_on_mesh<imfreq> TRIQS_INHERIT_AND_FORWARD_CONSTRUCTOR(evaluator_fnt_on_mesh, evaluator_grid_simple);
+  // simple evaluation : take the point on the grid...
+  template <> struct evaluator_fnt_on_mesh<imfreq> {
+   long n;
+   evaluator_fnt_on_mesh() = default;
+   template <typename MeshType> evaluator_fnt_on_mesh(MeshType const &m, long p) { n = p; }
+   template <typename MeshType> evaluator_fnt_on_mesh(MeshType const &m, matsubara_freq_mesh::mesh_point_t const &p) { n = p.n; }
+   template <typename F> auto operator()(F const &f) const DECL_AND_RETURN(f(n));
+  };
 
-  template<typename Opt, typename Target>
+  // ------------- evaluator  -------------------
-   struct evaluator<imfreq,Target,Opt> { // factorize and template on the long instead of double ?
+  // handle the case where the matsu. freq is out of grid...
-    public :
+  template <typename Target, typename Opt> struct evaluator<imfreq, Target, Opt> {
    static constexpr int arity = 1;
+   template <typename G> auto operator()(G const *g, int n) const DECL_AND_RETURN((*g)[n]);
    template <typename G>
-      auto operator()(G const * g, int n) const DECL_AND_RETURN((*g)[n]);
+   auto operator()(G const *g, matsubara_freq_mesh::mesh_point_t const &p) const DECL_AND_RETURN((*g)[p.index()]);
+
+   // dispatch for 2x2 cases : matrix/scalar and tail/no_tail ( true means no_tail)
+   template <typename G>
+   std::complex<double> _call_impl(G const *g, matsubara_freq const &f, scalar_valued, std::false_type) const {
+    if ((f.n >= 0) && (f.n < g->mesh().size())) return (*g)[f.n];
+    if ((f.n < 0) && (-f.n < g->mesh().size())) return conj((*g)[-f.n]);
+    return g->singularity().evaluate(f)(0, 0);
+   }
 
    template <typename G>
-      auto operator() (G const * g, linear_mesh<matsubara_domain<true>>::mesh_point_t const & p)  
+   std::complex<double> _call_impl(G const *g, matsubara_freq const &f, scalar_valued, std::true_type) const {
-      const DECL_AND_RETURN((*g)[p.index()]);
+    if ((f.n >= 0) && (f.n < g->mesh().size())) return (*g)[f.n];
+    if ((f.n < 0) && (-f.n < g->mesh().size())) return conj((*g)[-f.n]);
+    return 0;
+   }
 
    template <typename G>
-      typename G::singularity_t const & operator()(G const * g,freq_infty const &) const {return g->singularity();}
+   arrays::matrix_const_view<std::complex<double>> _call_impl(G const *g, matsubara_freq const &f, matrix_valued,
+                                                              std::false_type) const {
+    if ((f.n >= 0) && (f.n < g->mesh().size())) return (*g)[f.n]();
+    if ((f.n < 0) && (-f.n < g->mesh().size()))
+     return arrays::matrix<std::complex<double>>{conj((*g)[-f.n]())};
+    else
+     return g->singularity().evaluate(f);
+   }
+
+   template <typename G>
+   arrays::matrix_const_view<std::complex<double>> _call_impl(G const *g, matsubara_freq const &f, matrix_valued,
+                                                              std::true_type) const {
+    if ((f.n >= 0) && (f.n < g->mesh().size())) return (*g)[f.n]();
+    if ((f.n < 0) && (-f.n < g->mesh().size())) return arrays::matrix<std::complex<double>>{conj((*g)[-f.n]())};
+    auto r = arrays::matrix<std::complex<double>>{get_target_shape(*g)};
+    r() = 0;
+    return r;
+   }
+
+   template <typename G>
+   auto operator()(G const *g, matsubara_freq const &f) const
+       DECL_AND_RETURN(_call_impl(g, f, Target{}, std::integral_constant<bool, std::is_same<Opt, no_tail>::value>{}));
+
+   template <typename G> typename G::singularity_t const &operator()(G const *g, freq_infty const &) const {
+    return g->singularity();
+   }
   };
 
   /// ---------------------------  data access  ---------------------------------
@@ -74,6 +123,6 @@ namespace triqs { namespace gfs {
   template <typename Opt> struct data_proxy<imfreq, scalar_valued, Opt> : data_proxy_array<std::complex<double>, 1> {};
 
  } // gfs_implementation
-
+}
-}}
+}
 #endif

triqs/gfs/local/tail.hpp

@@ -126,6 +126,20 @@ namespace triqs { namespace gfs { namespace local {
 
       operator freq_infty() const { return freq_infty(); }
 
+      /// Evaluate the tail to  sum_{n=order_min}^ordermax M_n/omega^n 
+      arrays::matrix<dcomplex> evaluate(dcomplex const &omega) const {
+       auto r = arrays::matrix<dcomplex>{this->shape()};
+       r() = 0;
+       auto omin = this->order_min();
+       auto omax = this->order_max(); // precompute since long to do...
+       auto _ = arrays::range{};
+       for (int u = omax; u >= omin; --u)
+        r = r / omega +
+            const_mv_type{this->_data(u - omin, _, _)}; // need to make a matrix view because otherwise + is not defined
+       r /= pow(omega, omin);
+       return r;
+      }
+
       /// Save in txt file: doc the format  ? ---> prefer serialization or hdf5 !
       void save(std::string file, bool accumulate=false) const {}
 

triqs/gfs/meshes/matsubara_freq.hpp

@@ -0,0 +1,135 @@
+/*******************************************************************************
+ *
+ * TRIQS: a Toolbox for Research in Interacting Quantum Systems
+ *
+ * Copyright (C) 2013 by O. Parcollet
+ *
+ * TRIQS 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.
+ *
+ * TRIQS 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 should have received a copy of the GNU General Public License along with
+ * TRIQS. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#pragma once
+#include "./mesh_tools.hpp"
+#include "../domains/matsubara.hpp"
+
+namespace triqs {
+namespace gfs {
+
+ //----------------------------------- The mesh -----------------------------------------
+ struct matsubara_freq_mesh {
+
+  using domain_t = matsubara_domain<true>;
+  using index_t = int;
+  using domain_pt_t = typename domain_t::point_t;
+
+  matsubara_freq_mesh() : _dom(), n_max(0), _positive_only(true) {}
+
+  explicit matsubara_freq_mesh(domain_t dom, int n_pts, bool positive_only = true)
+     : _dom(std::move(dom)), n_max(n_pts), _positive_only(positive_only) {}
+
+  domain_t const &domain() const { return _dom; }
+  size_t size() const { return (_positive_only ? n_max : 2 * n_max + 1); }
+
+  /// Conversions point <-> index <-> linear_index
+  domain_pt_t index_to_point(index_t ind) const {
+   return 1_j * M_PI * (2 * ind + (_dom.statistic == Fermion ? 1 : 0)) / _dom.beta;
+  }
+
+  public:
+  size_t index_to_linear(index_t ind) const { return ind; }
+
+  /// The wrapper for the mesh point
+  struct mesh_point_t : tag::mesh_point, public arith_ops_by_cast<mesh_point_t, matsubara_freq> {
+   index_t n;
+   double beta;
+   statistic_enum statistic;
+   index_t n_max;
+
+   mesh_point_t()= default;
+   mesh_point_t(matsubara_freq_mesh const &mesh, index_t const &index_)
+      : n(index_), beta(mesh.domain().beta), statistic(mesh.domain().statistic), n_max(mesh.size()) {}
+   mesh_point_t(index_t const &index_, double beta_, statistic_enum stat_, index_t n_max_)
+      : n(index_), beta(beta_), statistic(stat_), n_max(n_max_) {}
+   void advance() { ++n; }
+   using cast_t = domain_pt_t;
+   operator cast_t() const { return 1_j * M_PI * (2 * n + (statistic == Fermion ? 1 : 0)) / beta; }
+   using cast1_t = matsubara_freq;
+   operator cast1_t() const { return {n,beta,statistic};}
+   size_t linear_index() const { return n; }
+   size_t index() const { return n; }
+   bool at_end() const { return (n == n_max); }
+   void reset() { n = 0; }
+  };
+
+  /// Accessing a point of the mesh
+  mesh_point_t operator[](index_t i) const { return mesh_point_t(*this, i); }
+
+  /// Iterating on all the points...
+  using const_iterator = mesh_pt_generator<matsubara_freq_mesh>;
+  const_iterator begin() const { return const_iterator(this); }
+  const_iterator end() const { return const_iterator(this, true); }
+  const_iterator cbegin() const { return const_iterator(this); }
+  const_iterator cend() const { return const_iterator(this, true); }
+
+  /// Mesh comparison
+  bool operator==(matsubara_freq_mesh const &M) const {
+   return (std::make_tuple(_dom, n_max, _positive_only, n_max) == std::make_tuple(M._dom, M.n_max, M._positive_only, n_max));
+  }
+  bool operator!=(matsubara_freq_mesh const &M) const { return !(operator==(M)); }
+
+  /// Write into HDF5
+  friend void h5_write(h5::group fg, std::string subgroup_name, matsubara_freq_mesh const &m) {
+   h5::group gr = fg.create_group(subgroup_name);
+   h5_write(gr, "domain", m.domain());
+   h5_write(gr, "size", m.size());
+   //h5_write(gr, "start_at_0", m._positive_only);
+   // kept for backward compatibility
+   auto beta = m.domain().beta;
+   h5_write(gr, "min", Fermion ? M_PI / beta : 0);
+   h5_write(gr, "max", ((Fermion ? 1 : 0) + 2 * m.size()) * M_PI / beta);
+   h5_write(gr, "kind", 2);
+  }
+
+  /// Read from HDF5
+  friend void h5_read(h5::group fg, std::string subgroup_name, matsubara_freq_mesh &m) {
+   h5::group gr = fg.open_group(subgroup_name);
+   typename matsubara_freq_mesh::domain_t dom;
+   int L;
+   bool s = true;
+   h5_read(gr, "domain", dom);
+   h5_read(gr, "size", L);
+   //h5_read(gr, "start_at_0", s);
+   m = matsubara_freq_mesh{std::move(dom), L, s};
+  }
+
+  //  BOOST Serialization
+  friend class boost::serialization::access;
+  template <class Archive> void serialize(Archive &ar, const unsigned int version) {
+   ar &boost::serialization::make_nvp("domain", _dom);
+   ar &boost::serialization::make_nvp("size", n_max);
+   ar &boost::serialization::make_nvp("kind", _positive_only);
+  }
+
+  friend std::ostream &operator<<(std::ostream &sout, matsubara_freq_mesh const &m) {
+   return sout << "Matsubara Freq Mesh of size " << m.size();
+  }
+
+  private:
+  domain_t _dom;
+  int n_max;
+  bool _positive_only;
+ };
+
+ }
+}
+