From 1ecec0b9339d8b468a7575ee4ccfba173eb719dd Mon Sep 17 00:00:00 2001 From: Olivier Parcollet Date: Sun, 5 Jan 2014 20:10:33 +0100 Subject: [PATCH] gfs: Fix bug in imfreq with <0 freq. - there was a confusion in gf imfreq, in the new case where freq can be <0 (non real gf, or for product gf). - index: is the matsubara n, as in the struct matsubara_freq index can be >0 or <0 - linear_index : is the shift from the 0. It is always >0. Fixed function to compute it. - Also changed the construction of mesh_point in the generic iterator. Before, was constructed with a mesh point of index 0 Now, added a new constructor on mesh_point_t, just taking the mesh which construct the *first* mesh_point. Fixed linear, discrete, product accordingly. Added to the documentation of the concepts of gf. --- doc/reference/c++/gf/concepts.rst | 50 +++++++++++++++-------------- triqs/gfs/meshes/discrete.hpp | 1 + triqs/gfs/meshes/linear.hpp | 1 + triqs/gfs/meshes/matsubara_freq.hpp | 5 +-- triqs/gfs/meshes/mesh_tools.hpp | 2 +- triqs/gfs/meshes/product.hpp | 4 +-- 6 files changed, 34 insertions(+), 29 deletions(-) diff --git a/doc/reference/c++/gf/concepts.rst b/doc/reference/c++/gf/concepts.rst index af13b50e..6ffbb2d3 100644 --- a/doc/reference/c++/gf/concepts.rst +++ b/doc/reference/c++/gf/concepts.rst @@ -146,30 +146,32 @@ MeshPoint * **Definition** : -+------------------------------------------------+-----------------------------------------------------------------------------+ -| Elements | Comment | -+================================================+=============================================================================+ -| mesh_t | Type of the mesh | -+------------------------------------------------+-----------------------------------------------------------------------------+ -| mesh_t const * m | A pointer to the mesh to which the point belongs. | -+------------------------------------------------+-----------------------------------------------------------------------------+ -| mesh_t::index_t index | The index of the point | -+------------------------------------------------+-----------------------------------------------------------------------------+ -| mesh_point_t( mesh_t const &, index_t const &) | Constructor | -+------------------------------------------------+-----------------------------------------------------------------------------+ -| mesh_t::index_t [const &|] index() const | The index corresponding to the point | -+------------------------------------------------+-----------------------------------------------------------------------------+ -| size_t linear_index() const | The linear index of the point (same as m->index_to_linear(index()) | -+------------------------------------------------+-----------------------------------------------------------------------------+ -| void advance() | Advance to the next point on the mesh (used by iterators). | -+------------------------------------------------+-----------------------------------------------------------------------------+ -| void at_end() | Is the point at the end of the grid | -+------------------------------------------------+-----------------------------------------------------------------------------+ -| void reset() | Reset the mesh point to the first point | -+------------------------------------------------+-----------------------------------------------------------------------------+ -| cast_t | == mesh_t::domain_t::point_t | -| operator cast_t() const | *implicit* cast to the corresponding domain point | -+------------------------------------------------+-----------------------------------------------------------------------------+ ++------------------------------------------------+--------------------------------------------------------------------+ +| Elements | Comment | ++================================================+====================================================================+ +| mesh_t | Type of the mesh | ++------------------------------------------------+--------------------------------------------------------------------+ +| mesh_t const * m | A pointer to the mesh to which the point belongs. | ++------------------------------------------------+--------------------------------------------------------------------+ +| mesh_t::index_t index | The index of the point | ++------------------------------------------------+--------------------------------------------------------------------+ +| mesh_point_t( mesh_t const &, index_t const &) | Constructor : a mesh point at the given index | ++------------------------------------------------+--------------------------------------------------------------------+ +| mesh_point_t( mesh_t const &) | Constructor : the first mesh point | ++------------------------------------------------+--------------------------------------------------------------------+ +| mesh_t::index_t [const &,] index() const | The index corresponding to the point | ++------------------------------------------------+--------------------------------------------------------------------+ +| size_t linear_index() const | The linear index of the point (same as m->index_to_linear(index()) | ++------------------------------------------------+--------------------------------------------------------------------+ +| void advance() | Advance to the next point on the mesh (used by iterators). | ++------------------------------------------------+--------------------------------------------------------------------+ +| void at_end() | Is the point at the end of the grid | ++------------------------------------------------+--------------------------------------------------------------------+ +| void reset() | Reset the mesh point to the first point | ++------------------------------------------------+--------------------------------------------------------------------+ +| cast_t operator cast_t() const | == mesh_t::domain_t::point_t *implicit* cast to the corresponding | +| | domain point | ++------------------------------------------------+--------------------------------------------------------------------+ For one dimensional mesh, we also require that the MeshPoint implement the basic arithmetic operations using the cast. diff --git a/triqs/gfs/meshes/discrete.hpp b/triqs/gfs/meshes/discrete.hpp index 4bf9a5a6..6e15738f 100644 --- a/triqs/gfs/meshes/discrete.hpp +++ b/triqs/gfs/meshes/discrete.hpp @@ -47,6 +47,7 @@ namespace gfs { public: mesh_point_t() = default; mesh_point_t(discrete_mesh const &mesh, index_t const &index_) : m(&mesh), _index(index_) {} + mesh_point_t(discrete_mesh const &mesh) : mesh_point_t(mesh, 0){} void advance() { ++_index; } using cast_t = long; operator cast_t() const { return m->index_to_point(_index); } diff --git a/triqs/gfs/meshes/linear.hpp b/triqs/gfs/meshes/linear.hpp index 89835b78..32465594 100644 --- a/triqs/gfs/meshes/linear.hpp +++ b/triqs/gfs/meshes/linear.hpp @@ -84,6 +84,7 @@ namespace gfs { public: mesh_point_t() : m(nullptr) {} mesh_point_t(linear_mesh const &mesh, index_t const &index_) : m(&mesh), _index(index_) {} + mesh_point_t(linear_mesh const &mesh) : mesh_point_t(mesh,0) {} void advance() { ++_index; } using cast_t = domain_pt_t; operator cast_t() const { return m->index_to_point(_index); } diff --git a/triqs/gfs/meshes/matsubara_freq.hpp b/triqs/gfs/meshes/matsubara_freq.hpp index 16b7ed60..617e6685 100644 --- a/triqs/gfs/meshes/matsubara_freq.hpp +++ b/triqs/gfs/meshes/matsubara_freq.hpp @@ -67,7 +67,7 @@ namespace gfs { domain_pt_t index_to_point(index_t ind) const { return 1_j * M_PI * (2 * ind + (_dom.statistic == Fermion)) / _dom.beta; } /// Flatten the index in the positive linear index for memory storage (almost trivial here). - long index_to_linear(index_t ind) const { return ind + index_start(); } + long index_to_linear(index_t ind) const { return ind - index_start(); } /** * The mesh point @@ -80,8 +80,9 @@ namespace gfs { : matsubara_freq(index_, mesh.domain().beta, mesh.domain().statistic), index_start(mesh.index_start()), index_stop(mesh.index_start() + mesh.size() - 1) {} + mesh_point_t(matsubara_freq_mesh const &mesh) : mesh_point_t(mesh, mesh.index_start()) {} void advance() { ++n; } - long linear_index() const { return n; } + long linear_index() const { return n - index_start; } long index() const { return n; } bool at_end() const { return (n == index_stop); } void reset() { n = index_start; } diff --git a/triqs/gfs/meshes/mesh_tools.hpp b/triqs/gfs/meshes/mesh_tools.hpp index d7b1c75c..8288c766 100644 --- a/triqs/gfs/meshes/mesh_tools.hpp +++ b/triqs/gfs/meshes/mesh_tools.hpp @@ -40,7 +40,7 @@ namespace gfs { //bool equal(mesh_pt_generator const & other) const { return ((mesh == other.mesh) && (other.u==u) );} public: mesh_pt_generator(): mesh(nullptr), u(0) {} - mesh_pt_generator( MeshType const * m, bool atEnd = false): mesh(m), u(atEnd ? m->size(): 0), pt((*m)[typename MeshType::index_t()]) {} + mesh_pt_generator( MeshType const * m, bool atEnd = false): mesh(m), u(atEnd ? m->size(): 0), pt(*m) {} void increment() { ++u; pt.advance(); } bool at_end() const { return (u>=mesh->size());} typename MeshType::domain_t::point_t to_point() const { return pt;} diff --git a/triqs/gfs/meshes/product.hpp b/triqs/gfs/meshes/product.hpp index ac8f4a42..96cfbb86 100644 --- a/triqs/gfs/meshes/product.hpp +++ b/triqs/gfs/meshes/product.hpp @@ -128,14 +128,14 @@ namespace gfs { template typename M::mesh_point_t operator()(M const &m, typename M::index_t const &i) const { return m[i]; } }; struct F1 { - template typename M::mesh_point_t operator()(M const &m) const { return m[typename M::index_t()]; } + template typename M::mesh_point_t operator()(M const &m) const { return {m}; } }; public: mesh_point_t() = default; mesh_point_t(mesh_product const &m_, index_t index_) : m(&m_), _c(triqs::tuple::apply_on_zip(F2(), m_.m_tuple, index_)), _atend(false) {} - mesh_point_t(mesh_product const &m_) : m(&m_), _c(triqs::tuple::apply(F1(), m_.m_tuple)), _atend(false) {} + mesh_point_t(mesh_product const &m_) : m(&m_), _c(triqs::tuple::apply_on_tuple(F1(), m_.m_tuple)), _atend(false) {} m_pt_tuple_t const &components_tuple() const { return _c; } size_t linear_index() const { return m->mp_to_linear(_c); } const mesh_product *mesh() const { return m; }