From 8c1f7945e458fdc1173827b614ca97dcf032eb07 Mon Sep 17 00:00:00 2001 From: tayral Date: Wed, 19 Feb 2014 11:59:47 +0000 Subject: [PATCH] Various small fixes in imfreq Gfs. - documentation of matsubara_freq_mesh - fixing doc of fit_tail - fixing right conventions for matsubara_freq_mesh --- doc/reference/c++/gf/gf_imfreq.rst | 4 +- doc/reference/c++/gf/matsubara_freq_mesh.rst | 89 ++++++++++++++++++++ doc/reference/c++/gf/set_tail_from_fit.rst | 4 +- test/triqs/gfs/test_fit_tail.output | 6 +- triqs/gfs/imfreq.hpp | 18 ++-- triqs/gfs/local/fit_tail.hpp | 16 ++-- triqs/gfs/meshes/matsubara_freq.hpp | 24 +++--- 7 files changed, 127 insertions(+), 34 deletions(-) create mode 100644 doc/reference/c++/gf/matsubara_freq_mesh.rst diff --git a/doc/reference/c++/gf/gf_imfreq.rst b/doc/reference/c++/gf/gf_imfreq.rst index 1fe8b2ca..e786b5d7 100644 --- a/doc/reference/c++/gf/gf_imfreq.rst +++ b/doc/reference/c++/gf/gf_imfreq.rst @@ -2,7 +2,7 @@ .. _gf_imfreq: -Matsubara frequencies +Green function on Matsubara frequencies ========================================================== This is a specialisation of :ref:`gf` for imaginary Matsubara frequencies. @@ -31,7 +31,7 @@ The domain is :doxy:`matsubara_freq_domain`. The Matsubara frequencies are :math:`\omega_n=\frac{(2n+1)\pi}{\beta}` for fermions and :math:`\omega_n=\frac{2n\pi}{\beta}` for bosons. -The mesh is :doxy:`matsubara_freq_mesh`. +The mesh is :doc:`matsubara_freq_mesh`. Singularity diff --git a/doc/reference/c++/gf/matsubara_freq_mesh.rst b/doc/reference/c++/gf/matsubara_freq_mesh.rst new file mode 100644 index 00000000..95f32e88 --- /dev/null +++ b/doc/reference/c++/gf/matsubara_freq_mesh.rst @@ -0,0 +1,89 @@ +Mesh of Matsubara frequencies +============================== + + +Synopsis +------------ +The class `matsubara_freq_mesh` is the implementation of a mesh of fermionic (:math:`\frac{(2n+1)\pi}{\beta}`) or bosonic (:math:`\frac{2n\pi}{\beta}`) Matsubara frequencies. + +The mesh can span either only positive frequencies or both positive and negative frequencies (which is necessary for complex functions :math:`G(\tau)`). + +Public member functions +------------------------ + ++------------------+---------------------------------------------------------------------------------------------+------------------------------------------------------------------------------------------+ +| type | name | description | ++==================+=============================================================================================+==========================================================================================+ +| | matsubara_freq_mesh() | constructor | ++------------------+---------------------------------------------------------------------------------------------+------------------------------------------------------------------------------------------+ +| | matsubara_freq_mesh(domain_t dom, int n_pts=1025, bool positive_only=true) | constructor | ++------------------+---------------------------------------------------------------------------------------------+------------------------------------------------------------------------------------------+ +| | matsubara_freq_mesh(double beta, statistic_enum S, int n_pts=1025, bool positive_only=true) | constructor | ++------------------+---------------------------------------------------------------------------------------------+------------------------------------------------------------------------------------------+ +| domain_t const & | domain() const | corresponding :doxy:`domain` | ++------------------+---------------------------------------------------------------------------------------------+------------------------------------------------------------------------------------------+ +| int | first_index() const | first Matsubara index | ++------------------+---------------------------------------------------------------------------------------------+------------------------------------------------------------------------------------------+ +| int | last_index() const | last Matsubara index | ++------------------+---------------------------------------------------------------------------------------------+------------------------------------------------------------------------------------------+ +| long | size() const | Size (linear) of the mesh | ++------------------+---------------------------------------------------------------------------------------------+------------------------------------------------------------------------------------------+ +| domain_pt_t | index_to_point(index_t ind) const | From an index of a point in the mesh, returns the corresponding point in the domain. | ++------------------+---------------------------------------------------------------------------------------------+------------------------------------------------------------------------------------------+ +| long | index_to_linear(index_t ind) const | Flatten the index in the positive linear index for memory storage (almost trivial here). | ++------------------+---------------------------------------------------------------------------------------------+------------------------------------------------------------------------------------------+ +| bool | positive_only() const | Is the mesh only for positive omega_n (G(tau) real)) | ++------------------+---------------------------------------------------------------------------------------------+------------------------------------------------------------------------------------------+ +| mesh_point_t | operator[](index_t i) const | Accessing a point of the mesh from its index. | ++------------------+---------------------------------------------------------------------------------------------+------------------------------------------------------------------------------------------+ +| const_iterator | begin() const | | ++------------------+---------------------------------------------------------------------------------------------+------------------------------------------------------------------------------------------+ +| const_iterator | end() const | | ++------------------+---------------------------------------------------------------------------------------------+------------------------------------------------------------------------------------------+ +| const_iterator | cbegin() const | | ++------------------+---------------------------------------------------------------------------------------------+------------------------------------------------------------------------------------------+ +| const_iterator | cend() const | | ++------------------+---------------------------------------------------------------------------------------------+------------------------------------------------------------------------------------------+ +| bool | operator== (matsubara_freq_mesh const &M) const | | ++------------------+---------------------------------------------------------------------------------------------+------------------------------------------------------------------------------------------+ +| bool | operator!= (matsubara_freq_mesh const &M) const | | ++------------------+---------------------------------------------------------------------------------------------+------------------------------------------------------------------------------------------+ + +Public member types +------------------------ + ++--------------+-------------------------+ +| type | description | ++==============+=========================+ +| mesh_point_t | type of the mesh point | ++--------------+-------------------------+ + + +Example +---------- + +.. compileblock:: + + #include + using namespace triqs::gfs; + int main(){ + double beta=1; + int n_pts=4; + matsubara_freq_mesh m(beta,Fermion, n_pts); + std::cout << "Mesh m with only positive Matsubara frequencies : " << std::endl; + for(auto & w:m) std::cout << w << std::endl; + std::cout << "m.first_index() = " << m.first_index() << std::endl; + std::cout << "m.last_index() = " << m.last_index() << std::endl; + std::cout << "m.size() = " << m.size() << std::endl; + matsubara_freq_mesh m2(beta,Fermion, n_pts, false); + std::cout << "Mesh m2 with positive and negative Matsubara frequencies : " << std::endl; + for(auto & w:m2) std::cout << w << std::endl; + std::cout << "m2.first_index() = " << m2.first_index() << std::endl; + std::cout << "m2.last_index() = " << m2.last_index() << std::endl; + std::cout << "m2.size() = " << m2.size() << std::endl; + return 0; + } + + + + diff --git a/doc/reference/c++/gf/set_tail_from_fit.rst b/doc/reference/c++/gf/set_tail_from_fit.rst index e10826dc..992b3db2 100644 --- a/doc/reference/c++/gf/set_tail_from_fit.rst +++ b/doc/reference/c++/gf/set_tail_from_fit.rst @@ -22,9 +22,9 @@ where +------------+----------------+----------------------------------------------------------------------+---------+ | int | n_moments | number of moments in the final tail (including known ones) | no | +------------+----------------+----------------------------------------------------------------------+---------+ -| size_t | n_min | linear index on mesh to start the fit | no | +| size_t | n_min | Matsubara index on mesh to start the fit | no | +------------+----------------+----------------------------------------------------------------------+---------+ -| size_t | n_max | final linear index for fit (included) | no | +| size_t | n_max | final Matsubara index for fit (included) | no | +------------+----------------+----------------------------------------------------------------------+---------+ | bool | replace_by_fit | if true, replace the gf data in the fitting range by the tail values | true | +------------+----------------+----------------------------------------------------------------------+---------+ diff --git a/test/triqs/gfs/test_fit_tail.output b/test/triqs/gfs/test_fit_tail.output index 81f5231e..15b38302 100644 --- a/test/triqs/gfs/test_fit_tail.output +++ b/test/triqs/gfs/test_fit_tail.output @@ -66,7 +66,7 @@ ... Order 2 = [[(1,0)]] ... Order 3 = -[[(0.999247,0)]] +[[(0.999288,0)]] ... Order 4 = [[(0.998655,0)]] @@ -80,7 +80,7 @@ ... Order 2 = [[(1,0)]] ... Order 3 = -[[(0.999186,0)]] +[[(0.999209,0)]] ... Order 4 = [[(0.998631,0)]] @@ -94,7 +94,7 @@ ... Order 2 = [[(1,0)]] ... Order 3 = -[[(0.999214,0)]] +[[(0.999251,0)]] ... Order 4 = [[(0.998655,0)]] diff --git a/triqs/gfs/imfreq.hpp b/triqs/gfs/imfreq.hpp index f1f7e558..885a08e6 100644 --- a/triqs/gfs/imfreq.hpp +++ b/triqs/gfs/imfreq.hpp @@ -59,7 +59,7 @@ namespace gfs { evaluator_fnt_on_mesh() = default; template evaluator_fnt_on_mesh(MeshType const &m, long p) { n = p; w=1; } template evaluator_fnt_on_mesh(MeshType const &m, matsubara_freq const &p) { - if ((p.n >= m.index_start()) && (p.n < m.size()+m.index_start())) {w=1; n =p.n;} + if ((p.n >= m.first_index()) && (p.n < m.size()+m.first_index())) {w=1; n =p.n;} else {w=0; n=0;} } template auto operator()(F const &f) const DECL_AND_RETURN(w*f(n)); @@ -76,24 +76,24 @@ namespace gfs { // 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 (g->mesh().index_start()==0){//only positive Matsubara frequencies + if (g->mesh().positive_only()){//only positive Matsubara frequencies 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]); } else{ - if ((f.n >= g->mesh().index_start()) && (f.n < g->mesh().size()+g->mesh().index_start())) return (*g)[f.n]; + if ((f.n >= g->mesh().first_index()) && (f.n < g->mesh().size()+g->mesh().first_index())) return (*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 (g->mesh().index_start()==0){//only positive Matsubara frequencies + if (g->mesh().positive_only()){//only positive Matsubara frequencies 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]); } else{ - if ((f.n >= g->mesh().index_start()) && (f.n < g->mesh().size()+g->mesh().index_start())) return (*g)[f.n]; + if ((f.n >= g->mesh().first_index()) && (f.n < g->mesh().size()+g->mesh().first_index())) return (*g)[f.n]; } return 0; } @@ -101,13 +101,13 @@ namespace gfs { template arrays::matrix_const_view> _call_impl(G const *g, matsubara_freq const &f, matrix_valued, std::false_type) const { - if (g->mesh().index_start()==0){//only positive Matsubara frequencies + if (g->mesh().positive_only()){//only positive Matsubara frequencies 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{ - if ((f.n >= g->mesh().index_start()) && (f.n < g->mesh().size()+g->mesh().index_start())) return (*g)[f.n]; + if ((f.n >= g->mesh().first_index()) && (f.n < g->mesh().size()+g->mesh().first_index())) return (*g)[f.n]; } return g->singularity().evaluate(f); } @@ -115,13 +115,13 @@ namespace gfs { template arrays::matrix_const_view> _call_impl(G const *g, matsubara_freq const &f, matrix_valued, std::true_type) const { - if (g->mesh().index_start()==0){//only positive Matsubara frequencies + if (g->mesh().positive_only()){//only positive Matsubara frequencies 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{ - if ((f.n >= g->mesh().index_start()) && (f.n < g->mesh().size()+g->mesh().index_start())) return (*g)[f.n]; + if ((f.n >= g->mesh().first_index()) && (f.n < g->mesh().size()+g->mesh().first_index())) return (*g)[f.n]; } auto r = arrays::matrix>{get_target_shape(*g)}; r() = 0; diff --git a/triqs/gfs/local/fit_tail.hpp b/triqs/gfs/local/fit_tail.hpp index 86e07fb1..f6c0ffcd 100644 --- a/triqs/gfs/local/fit_tail.hpp +++ b/triqs/gfs/local/fit_tail.hpp @@ -78,8 +78,8 @@ namespace triqs { namespace gfs { namespace local { const double rcond = 0.0; int rank; - for (size_t i = 0; i < get_target_shape(gf)[0]; i++) { - for (size_t j = 0; j < get_target_shape(gf)[1]; j++) { + for (int i = 0; i < get_target_shape(gf)[0]; i++) { + for (int j = 0; j < get_target_shape(gf)[1]; j++) { // fit the odd moments // S.resize(size_odd); @@ -132,12 +132,12 @@ namespace triqs { namespace gfs { namespace local { return res; // return tail } - void set_tail_from_fit(gf_view gf, tail_view known_moments, int n_moments, size_t n_min, size_t n_max, + void set_tail_from_fit(gf_view gf, tail_view known_moments, int n_moments, int n_min, int n_max, bool replace_by_fit = false) { if (get_target_shape(gf) != known_moments.shape()) TRIQS_RUNTIME_ERROR << "shape of tail does not match shape of gf"; gf.singularity() = fit_tail_impl(gf, known_moments, n_moments, n_min, n_max); if (replace_by_fit) { // replace data in the fitting range by the values from the fitted tail - size_t i = 0; + int i = 0; for (auto iw : gf.mesh()) { // (arrays::range(n_min,n_max+1)) { if ((i >= n_min) && (i <= n_max)) gf[iw] = gf.singularity().evaluate(iw); i++; @@ -145,13 +145,13 @@ namespace triqs { namespace gfs { namespace local { } } - void set_tail_from_fit(gf_view> block_gf, tail_view known_moments, int n_moments, size_t n_min, - size_t n_max, bool replace_by_fit = false) { + void set_tail_from_fit(gf_view> block_gf, tail_view known_moments, int n_moments, int n_min, + int n_max, bool replace_by_fit = false) { // for(auto &gf : block_gf) set_tail_from_fit(gf, known_moments, n_moments, n_min, n_max, replace_by_fit); - for (size_t i = 0; i < block_gf.mesh().size(); i++) + for (int i = 0; i < block_gf.mesh().size(); i++) set_tail_from_fit(block_gf[i], known_moments, n_moments, n_min, n_max, replace_by_fit); } - void set_tail_from_fit(gf_view gf, tail_view known_moments, int n_moments, size_t n_min, size_t n_max, bool replace_by_fit = false) { + void set_tail_from_fit(gf_view gf, tail_view known_moments, int n_moments, int n_min, int n_max, bool replace_by_fit = false) { set_tail_from_fit(reinterpret_scalar_valued_gf_as_matrix_valued(gf), known_moments, n_moments, n_min, n_max, replace_by_fit ); } diff --git a/triqs/gfs/meshes/matsubara_freq.hpp b/triqs/gfs/meshes/matsubara_freq.hpp index 6db0eccf..62ded5c7 100644 --- a/triqs/gfs/meshes/matsubara_freq.hpp +++ b/triqs/gfs/meshes/matsubara_freq.hpp @@ -57,8 +57,13 @@ namespace gfs { * else : * For fermions : -Nmax - 1 * For Bosons : -Nmax - **/ - int index_start() const { return -(_positive_only ? 0 : n_max() + (_dom.statistic == Fermion)); } + **/ + + /// last Matsubara index + int last_index() const { return (_positive_only ? _n_pts : (_n_pts - (_dom.statistic == Boson ? 1 : 2))/2);} + + /// first Matsubara index + int first_index() const { return -(_positive_only ? 0 : last_index() + (_dom.statistic == Fermion)); } /// Size (linear) of the mesh long size() const { return _n_pts;} @@ -67,7 +72,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 - first_index(); } /// Is the mesh only for positive omega_n (G(tau) real)) bool positive_only() const { return _positive_only;} @@ -81,17 +86,17 @@ namespace gfs { mesh_point_t() = default; mesh_point_t(matsubara_freq_mesh const &mesh, index_t const &index_) : 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()) {} + first_index(mesh.first_index()), + index_stop(mesh.first_index() + mesh.size() - 1) {} + mesh_point_t(matsubara_freq_mesh const &mesh) : mesh_point_t(mesh, mesh.first_index()) {} void advance() { ++n; } - long linear_index() const { return n - index_start; } + long linear_index() const { return n - first_index; } long index() const { return n; } bool at_end() const { return (n == index_stop); } - void reset() { n = index_start; } + void reset() { n = first_index; } private: - index_t index_start, index_stop; + index_t first_index, index_stop; }; /// Accessing a point of the mesh from its index @@ -155,7 +160,6 @@ namespace gfs { private: domain_t _dom; int _n_pts; - long n_max() const { return (_positive_only ? _n_pts : (_n_pts - (_dom.statistic == Boson ? 1 : 2))/2);} bool _positive_only; };