From 9ce291d640e0d86865dbc20031fecece94ea437e Mon Sep 17 00:00:00 2001 From: Olivier Parcollet Date: Tue, 5 Nov 2013 08:56:42 +0100 Subject: [PATCH] gf: update gf imfreq for g( om - nu) case DRAFT : to be tested further... - update gf - write a specific mesh for matsubara frequencies - now the cast series is : mesh_pt --> matsubara_freq --> complex - 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 | 20 +++++ triqs/gfs/domains/matsubara.hpp | 104 +++++++++++++++------ triqs/gfs/imfreq.hpp | 117 +++++++++++++++++------- triqs/gfs/local/tail.hpp | 14 +++ triqs/gfs/meshes/matsubara_freq.hpp | 135 ++++++++++++++++++++++++++++ 5 files changed, 330 insertions(+), 60 deletions(-) create mode 100644 triqs/gfs/meshes/matsubara_freq.hpp diff --git a/test/triqs/gfs/gfv2.cpp b/test/triqs/gfs/gfv2.cpp index e55e68ca..c3ad15a9 100644 --- a/test/triqs/gfs/gfv2.cpp +++ b/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{{beta, Fermion, 100}}; + auto D0w = gf{{beta, Boson, 100}}; + auto Sigma_w = gf{{beta, Fermion, 100}}; + G0w(om_) << 1 / (om_ - 3); + + for (auto const& nu : D0w.mesh()) Sigma_w(om_) << 2 * G0w(om_ - nu); + } + + // + { + auto G0w = gf{{beta, Fermion, 100}, {1,1}}; + auto D0w = gf{{beta, Boson, 100}, {1,1}}; + auto Sigma_w = gf{{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; diff --git a/triqs/gfs/domains/matsubara.hpp b/triqs/gfs/domains/matsubara.hpp index 5ef02afe..4d6b2bf8 100644 --- a/triqs/gfs/domains/matsubara.hpp +++ b/triqs/gfs/domains/matsubara.hpp @@ -21,49 +21,101 @@ #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> { + 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; + 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 auto operator OP(matsubara_freq const &x, Y &&y)->decltype(std::complex{}) { \ + return std::complex(x) OP std::forward(y); \ + } \ + template \ + auto operator OP(Y &&y, matsubara_freq const &x)->decltype(std::forward(y) OP std::complex{}) { \ + return std::forward(y) OP std::complex(x); \ + } + IMPL_OP(+); + IMPL_OP(-); + IMPL_OP(*); + IMPL_OP(/ ); +#undef IMPL_OP +*/ + +/* + * template auto operator OP(matsubara_freq const &x, Y &&y)->decltype(std::complex{}) { \ + return std::complex(x) OP std::forward(y); \ + } \ + template \ + auto operator OP(Y &&y, matsubara_freq const &x) \ + ->TYPE_DISABLE_IF(decltype(std::forward(y) OP std::complex{}), \ + std::is_same::type>::type, matsubara_freq>) { \ + return std::forward(y) OP std::complex(x); \ + } +*/ + + //--------------------------------------------------------------------------------------------------------- /// The domain - template - struct matsubara_domain { - typedef typename mpl::if_c, double>::type point_t; + template struct matsubara_domain { + using point_t = typename std::conditional, double>::type; double beta; statistic_enum statistic; matsubara_domain() = default; - matsubara_domain (double Beta, statistic_enum s): beta(Beta), statistic(s){ - if(beta<0)TRIQS_RUNTIME_ERROR<<"Matsubara domain construction : beta <0 : beta ="<< beta <<"\n"; + matsubara_domain(double Beta, statistic_enum s) : beta(Beta), statistic(s) { + if (beta < 0) TRIQS_RUNTIME_ERROR << "Matsubara domain construction : beta <0 : beta =" << beta << "\n"; } matsubara_domain(matsubara_domain const &) = default; - matsubara_domain(matsubara_domain 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));} + matsubara_domain(matsubara_domain 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 - friend void h5_write (h5::group fg, std::string subgroup_name, matsubara_domain const & d) { - h5::group gr = fg.create_group(subgroup_name); - h5_write(gr,"beta",d.beta); - h5_write(gr,"statistic",(d.statistic==Fermion ? "F" : "B")); + friend void h5_write(h5::group fg, std::string subgroup_name, matsubara_domain const &d) { + h5::group gr = fg.create_group(subgroup_name); + h5_write(gr, "beta", d.beta); + h5_write(gr, "statistic", (d.statistic == Fermion ? "F" : "B")); } /// Read from HDF5 - friend void h5_read (h5::group fg, std::string subgroup_name, matsubara_domain & d){ + 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; - h5_read(gr,"beta",beta); - h5_read(gr,"statistic",statistic); - d = matsubara_domain(beta,(statistic=="F" ? Fermion : Boson)); + double beta; + std::string statistic; + h5_read(gr, "beta", beta); + h5_read(gr, "statistic", statistic); + d = matsubara_domain(beta, (statistic == "F" ? Fermion : Boson)); } - + // BOOST Serialization friend class boost::serialization::access; - template - void serialize(Archive & ar, const unsigned int version) { - ar & boost::serialization::make_nvp("beta",beta); - ar & boost::serialization::make_nvp("statistic",statistic); - } - + template void serialize(Archive &ar, const unsigned int version) { + ar &boost::serialization::make_nvp("beta", beta); + ar &boost::serialization::make_nvp("statistic", statistic); + } }; - -}} +} +} #endif diff --git a/triqs/gfs/imfreq.hpp b/triqs/gfs/imfreq.hpp index 8a16f731..e12f6bdd 100644 --- a/triqs/gfs/imfreq.hpp +++ b/triqs/gfs/imfreq.hpp @@ -24,56 +24,105 @@ #include "./gf.hpp" #include "./local/tail.hpp" #include "./local/no_tail.hpp" -#include "./domains/matsubara.hpp" -#include "./meshes/linear.hpp" +#include "./meshes/matsubara_freq.hpp" #include "./evaluators.hpp" -namespace triqs { namespace gfs { +namespace triqs { +namespace gfs { struct imfreq {}; - template struct gf_mesh : linear_mesh> { - typedef linear_mesh> B; - static double m1(double beta) { return std::acos(-1)/beta;} + template struct gf_mesh : matsubara_freq_mesh { + 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) : - 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){} + gf_mesh(typename B::domain_t const &d, int Nmax = 1025) : B(d, Nmax, true) {} + gf_mesh(double beta, statistic_enum S, int Nmax = 1025) : gf_mesh({beta, S}, Nmax) {} }; - namespace gfs_implementation { + namespace gfs_implementation { - //singularity - template<> struct singularity { typedef local::tail type;}; - template<> struct singularity { typedef local::tail type;}; + // singularity + template <> struct singularity { + typedef local::tail type; + }; + template <> struct singularity { + typedef local::tail type; + }; - //h5 name - template struct h5_name { static std::string invoke(){ return "ImFreq";}}; + // h5 name + template struct h5_name { + static std::string invoke() { return "ImFreq"; } + }; /// --------------------------- evaluator --------------------------------- - template<> struct evaluator_fnt_on_mesh 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 { + long n; + evaluator_fnt_on_mesh() = default; + template evaluator_fnt_on_mesh(MeshType const &m, long p) { n = p; } + template evaluator_fnt_on_mesh(MeshType const &m, matsubara_freq_mesh::mesh_point_t const &p) { n = p.n; } + template auto operator()(F const &f) const DECL_AND_RETURN(f(n)); + }; - template - struct evaluator { // factorize and template on the long instead of double ? - public : - static constexpr int arity = 1; - template - auto operator()(G const * g, int n) const DECL_AND_RETURN((*g)[n]); + // ------------- evaluator ------------------- + // handle the case where the matsu. freq is out of grid... + template struct evaluator { + static constexpr int arity = 1; + template auto operator()(G const *g, int n) const DECL_AND_RETURN((*g)[n]); + template + auto operator()(G const *g, matsubara_freq_mesh::mesh_point_t const &p) const DECL_AND_RETURN((*g)[p.index()]); - template - auto operator() (G const * g, linear_mesh>::mesh_point_t const & p) - const DECL_AND_RETURN((*g)[p.index()]); - - template - typename G::singularity_t const & operator()(G const * g,freq_infty const &) const {return g->singularity();} - }; + // dispatch for 2x2 cases : matrix/scalar and tail/no_tail ( true means no_tail) + template + std::complex _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 + std::complex _call_impl(G const *g, matsubara_freq const &f, scalar_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 conj((*g)[-f.n]); + return 0; + } + + template + arrays::matrix_const_view> _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>{conj((*g)[-f.n]())}; + else + return g->singularity().evaluate(f); + } + + template + arrays::matrix_const_view> _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>{conj((*g)[-f.n]())}; + auto r = arrays::matrix>{get_target_shape(*g)}; + r() = 0; + return r; + } + + template + auto operator()(G const *g, matsubara_freq const &f) const + DECL_AND_RETURN(_call_impl(g, f, Target{}, std::integral_constant::value>{})); + + template typename G::singularity_t const &operator()(G const *g, freq_infty const &) const { + return g->singularity(); + } + }; /// --------------------------- data access --------------------------------- - template struct data_proxy : data_proxy_array,3> {}; - template struct data_proxy : data_proxy_array,1> {}; + template struct data_proxy : data_proxy_array, 3> {}; + template struct data_proxy : data_proxy_array, 1> {}; - } // gfs_implementation - -}} + } // gfs_implementation +} +} #endif diff --git a/triqs/gfs/local/tail.hpp b/triqs/gfs/local/tail.hpp index fa26a619..7bf728ae 100644 --- a/triqs/gfs/local/tail.hpp +++ b/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 evaluate(dcomplex const &omega) const { + auto r = arrays::matrix{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 {} diff --git a/triqs/gfs/meshes/matsubara_freq.hpp b/triqs/gfs/meshes/matsubara_freq.hpp new file mode 100644 index 00000000..2778bc8a --- /dev/null +++ b/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 . + * + ******************************************************************************/ +#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; + 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 { + 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; + 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 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; + }; + + } +} +