LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-10-31 11:13:46 +01:00
Code Issues Releases Wiki Activity

[API BREAK] gf : mesh_pt into []

Browse Source 
- All access to grid is now using the [] operator.
() always returns const, and is an evaluation over the domain.

- Now :
 * () is always a call to evaluator (or a lazy expression).
 * mesh_pt accepted by [] , both const and non const,
and not by ().

- Ported the libs (fourier et al.) but this break API.
This commit is contained in:
Olivier Parcollet 2013-07-28 13:47:42 +02:00
parent 40f8cb5c18
commit 7847b71552
20 changed files with 106 additions and 82 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
doc/reference/c++/gf/concepts.rst
View File

@ -168,7 +168,8 @@ MeshPoint concept
+------------------------------------------------+-----------------------------------------------------------------------------+
| void reset() | Reset the mesh point to the first point |
+------------------------------------------------+-----------------------------------------------------------------------------+
| operator mesh_t::domain_t::point_t() const | cast to the corresponding domain point |
| cast_t | type of the corresponding domain point |
| operator cast_t() const | cast to the corresponding domain point |
+------------------------------------------------+-----------------------------------------------------------------------------+
| Implements the basic operations on the domain | Only for non composite mesh |
| by using the cast operation | |

2
test/triqs/gf/gf_2times.cpp
View File

@ -34,7 +34,7 @@ int main() {
//G2(t_,tp_) << t_ + 3*tp_;
TEST( G(1,1) );
TEST( G(G.mesh()(1,1)) );
TEST( G[G.mesh()(1,1) ]);
TEST( G.on_mesh(1,1));
//G2(t_,tp_) << G(tp_,t_);

2
test/triqs/gf/gf_2times_b.cpp
View File

@ -20,7 +20,7 @@ int main(){
int nt=tmax/dt;
auto R= make_gf<two_real_times> (tmax,nt,make_shape(1,1));//results
for(auto point:R.mesh()) R(point)=0;
for(auto point:R.mesh()) R[point]=0;
const auto rg = on_mesh(R);
R.on_mesh(1,1) = 10;

4
test/triqs/gf/gf_re_im_freq_time.cpp
View File

@ -47,7 +47,7 @@ int main() {
G_w.singularity()(2)=triqs::arrays::matrix<double>{{0}};
auto G_w2 = slice_mesh_imtime(G_w_tau, index);
for(auto& w:G_w.mesh())
if ( std::abs(G_w(w)-G_w2(w)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_slice error : w="<< w <<" ,G_w="<< G_w(w)<<" ,G_w2="<< G_w2(w) <<"\n";
if ( std::abs(G_w[w]-G_w2[w]) > precision) TRIQS_RUNTIME_ERROR<<" fourier_slice error : w="<< w <<" ,G_w="<< G_w[w]<<" ,G_w2="<< G_w2[w] <<"\n";
//test of the interpolation
std::cout << G_t_tau(0.789,0.123) << std::endl;
@ -70,6 +70,6 @@ int main() {
auto G_t = inverse_fourier(G_w);
auto G_t2 = inverse_fourier(slice_mesh_imtime(G_w_tau, index) );
for(auto& t:G_t.mesh())
if ( std::abs(G_t(t)-G_t2(t)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_slice_re_time error : t="<< t <<" ,G_t="<< G_t(t) <<" ,G_t2="<< G_t2(t) <<"\n";
if ( std::abs(G_t[t]-G_t2[t]) > precision) TRIQS_RUNTIME_ERROR<<" fourier_slice_re_time error : t="<< t <<" ,G_t="<< G_t[t] <<" ,G_t2="<< G_t2[t] <<"\n";
}

4
test/triqs/gf/gf_retw.cpp
View File

@ -26,8 +26,8 @@ int main() {
auto Gt2 = make_gf<retime,scalar_valued> (0, tmax, N);
int i =0;
for (auto & t : Gt.mesh()) Gt(t) = 1.0*t;
for (auto & w : Gw.mesh()) Gw(w) = 1.0*w;
for (auto & t : Gt.mesh()) Gt[t] = 1.0*t;
for (auto & w : Gw.mesh()) Gw[w] = 1.0*w;
triqs::clef::placeholder<0> t_;
triqs::clef::placeholder<1> w_;

2
test/triqs/gf/pb_affichage.cpp
View File

@ -20,7 +20,7 @@ int main(){
int nt=tmax/dt;
auto R= make_gf<two_real_times> (tmax,nt,make_shape(1,1));//results
for(auto point:R.mesh()) R(point)=0;
for(auto point:R.mesh()) R[point]=0;
const auto rg = on_mesh(R);
R.on_mesh(1,1) = 10;

8
test/triqs/gf/test_fourier_matsubara.cpp
View File

@ -27,7 +27,7 @@ int main() {
auto Gw1 = make_gf<imfreq> (beta, Fermion, make_shape(1,1), N);
Gw1(om_) << 1/(om_-E);
// for(auto const& w:Gw1.mesh()){
// std::cout<<"w="<<std::complex<double>(w)<<", Gw1=" << Gw1(w)(0,0)<<std::endl;
// std::cout<<"w="<<std::complex<double>(w)<<", Gw1=" << Gw1[w](0,0)<<std::endl;
// }
h5_write(file, "Gw1", Gw1); // the original lorentzian
@ -44,14 +44,14 @@ int main() {
for(auto const& w:Gw1.mesh()){
// std::cout<<"w="<<std::complex<double>(w)<<",Gw1b=" << Gw1b(w)(0,0)<<std::endl;
// std::cout<<"w="<<std::complex<double>(w)<<",Delta Gw1b=" << Gw1b(w)(0,0)-Gw1(w)(0,0)<<std::endl;
if ( std::abs(Gw1b(w)(0,0)-Gw1(w)(0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_matsubara error : w="<<std::complex<double>(w)<<" ,Gw1b="<<std::abs(Gw1b(w)(0,0))<<"\n";
if ( std::abs(Gw1b[w](0,0)-Gw1[w](0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_matsubara error : w="<<std::complex<double>(w)<<" ,Gw1b="<<std::abs(Gw1b[w](0,0))<<"\n";
}
h5_write(file,"Gw1b",Gw1b); // must be 0
///verification that the TF is OK
for(auto const & t:Gt1.mesh()){
Gt1(t)-= exp(-E*t) * ( (t>0?-1:0)+1/(1+exp(E*beta)) );
if ( std::abs(Gt1(t)(0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_matsubara error : t="<<t<<" ,G1="<<std::abs(Gt1(t)(0,0))<<"\n";
Gt1[t]-= exp(-E*t) * ( (t>0?-1:0)+1/(1+exp(E*beta)) );
if ( std::abs(Gt1[t](0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_matsubara error : t="<<t<<" ,G1="<<std::abs(Gt1[t](0,0))<<"\n";
}
h5_write(file,"Gt1b",Gt1); // must be 0

24
test/triqs/gf/test_fourier_real_time.cpp
View File

@ -35,7 +35,7 @@ int main() {
auto Gw1 = make_gf<refreq> (-wmax, wmax, Nw, make_shape(1,1),triqs::gfs::full_bins);
double a = Gw1.mesh().delta() * sqrt( Gw1.mesh().size() );
for(auto const & w:Gw1.mesh()) Gw1(w)=lorentzian(w,a);
for(auto const & w:Gw1.mesh()) Gw1[w]=lorentzian(w,a);
Gw1.singularity()(2)=triqs::arrays::matrix<double>{{2.0*a}};
h5_write(file,"Gw1",Gw1); // the original lorentzian
@ -45,15 +45,15 @@ int main() {
// verification that TF(TF^-1)=Id
auto Gw1b = fourier(Gt1);
for(auto const & w:Gw1b.mesh()){
Gw1b(w)-=Gw1(w);
if ( std::abs(Gw1b(w)(0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_real_time error : w="<<w<<" ,G1="<<std::abs(Gw1b(w)(0,0))<<"\n";
Gw1b[w]-=Gw1[w];
if ( std::abs(Gw1b[w](0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_real_time error : w="<<w<<" ,G1="<<std::abs(Gw1b[w](0,0))<<"\n";
}
h5_write(file,"Gw1b",Gw1b); // must be 0
// verification that TF is the lorentzian_inverse function
for(auto const & t:Gt1.mesh()){
Gt1(t)-=lorentzian_inverse(t,a);
if ( std::abs(Gt1(t)(0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_real_time error : t="<<t<<" ,G1="<<std::abs(Gt1(t)(0,0))<<"\n";
Gt1[t]-=lorentzian_inverse(t,a);
if ( std::abs(Gt1[t](0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_real_time error : t="<<t<<" ,G1="<<std::abs(Gt1[t](0,0))<<"\n";
}
h5_write(file,"Gt1b",Gt1); // must be 0
@ -66,8 +66,8 @@ int main() {
auto Gt2 = make_gf<retime> (-tmax, tmax, Nt, make_shape(1,1));
a = 2*acos(-1.) / ( Gt2.mesh().delta() *sqrt( Gt2.mesh().size() ) );
for(auto const & t:Gt2.mesh()) Gt2(t) = 0.5 *I * ( lorentzian_inverse(-t,a)*theta(-t)-lorentzian_inverse(t,a)*theta(t) );
//for(auto const & t:Gt2.mesh()) Gt2(t) = 0.5_j * ( lorentzian_inverse(-t,a)*theta(-t)-lorentzian_inverse(t,a)*theta(t) );
for(auto const & t:Gt2.mesh()) Gt2[t] = 0.5 *I * ( lorentzian_inverse(-t,a)*theta(-t)-lorentzian_inverse(t,a)*theta(t) );
//for(auto const & t:Gt2.mesh()) Gt2[t] = 0.5_j * ( lorentzian_inverse(-t,a)*theta(-t)-lorentzian_inverse(t,a)*theta(t) );
Gt2.singularity()(1)=triqs::arrays::matrix<double>{{1.0}};
h5_write(file,"Gt2",Gt2);
@ -75,9 +75,9 @@ int main() {
h5_write(file,"Gw2",Gw2);
for(auto const & w:Gw2.mesh()){
Gw2(w)-= 0.5/(w+a*I)+0.5/(w-a*I);
//Gw2(w)-= 0.5/(w+a*1_j)+0.5/(w-a*1_j);
if ( std::abs(Gw2(w)(0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_real_time error : w="<<w<<" ,G2="<<std::abs(Gw2(w)(0,0))<<"\n";
Gw2[w]-= 0.5/(w+a*I)+0.5/(w-a*I);
//Gw2[w]-= 0.5/(w+a*1_j)+0.5/(w-a*1_j);
if ( std::abs(Gw2[w](0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_real_time error : w="<<w<<" ,G2="<<std::abs(Gw2[w](0,0))<<"\n";
}
h5_write(file,"Gw2b",Gw2);
@ -88,8 +88,8 @@ int main() {
tmax=4*acos(-1.);
auto Gt3 = make_gf<retime> (-tmax, tmax, Nt, make_shape(1,1));
for(auto const & t:Gt3.mesh()) Gt3(t) = 1.0 * std::cos(10*t) + 0.25*std::sin(4*t) + 0.5 * I*std::sin(8*t+0.3*acos(-1.)) ;
//for(auto const & t:Gt3.mesh()) Gt3(t) = 1.0 * std::cos(10*t) + 0.25*std::sin(4*t) + 0.5_j*std::sin(8*t+0.3*acos(-1.)) ;
for(auto const & t:Gt3.mesh()) Gt3[t] = 1.0 * std::cos(10*t) + 0.25*std::sin(4*t) + 0.5 * I*std::sin(8*t+0.3*acos(-1.)) ;
//for(auto const & t:Gt3.mesh()) Gt3[t] = 1.0 * std::cos(10*t) + 0.25*std::sin(4*t) + 0.5_j*std::sin(8*t+0.3*acos(-1.)) ;
h5_write(file,"Gt3",Gt3);
auto Gw3 = fourier(Gt3);

14
triqs/clef/clef.hpp
View File

@ -382,6 +382,20 @@ namespace triqs { namespace clef {
make_expr_call(Obj&& obj, Args &&... args)
{ return typename result_of::make_expr_call<Obj,Args...>::type (tags::function(),std::forward<Obj>(obj), std::forward<Args>(args)...);}
/* --------------------------------------------------------------------------------------------------
* Create a [] call (subscript) node of an object
* The object can be kept as a : a ref, a copy, a view
* --------------------------------------------------------------------------------------------------- */
namespace result_of {
template< typename Obj, typename Arg> struct make_expr_subscript :
std::enable_if< is_any_lazy<Arg>::value, expr<tags::subscript,typename remove_cv_ref<Obj>::type, typename remove_cv_ref<Arg>::type> > {};
}
template< typename Obj, typename Arg>
typename result_of::make_expr_subscript<Obj,Arg>::type
make_expr_subscript(Obj&& obj, Arg && arg)
{ return typename result_of::make_expr_subscript<Obj,Arg>::type (tags::subscript(),std::forward<Obj>(obj), std::forward<Arg>(arg));}
/* --------------------------------------------------------------------------------------------------
* function class : stores any expression polymorphically
* f(x_,y_ ) = an expression associates this expression dynamically to f, which

2
triqs/gf/block.hpp
View File

@ -142,7 +142,7 @@ namespace triqs { namespace gfs {
big_gf_t big_gf;
mesh_iterator_t mesh_it;
typename Target::view_type const & dereference() const { return big_gf(*mesh_it);}
typename Target::view_type const & dereference() const { return big_gf[*mesh_it];}
bool equal(block_gf_iterator const & other) const { return ((mesh_it == other.mesh_it));}
public:
block_gf_iterator(gf_view<block_index,Target,Opt> bgf, bool at_end = false): big_gf(std::move(bgf)), mesh_it(&big_gf.mesh(),at_end) {}

71
triqs/gf/gf.hpp
View File

@ -172,7 +172,6 @@ namespace triqs { namespace gfs {
typename std::add_const<
typename boost::lazy_disable_if< // disable the template if one the following conditions it true
boost::mpl::or_< // starting condition [OR]
std::is_base_of< tag::mesh_point, typename std::remove_reference<Arg0>::type>, // Arg0 is (a & or a &&) to a mesh_point_t
clef::is_any_lazy<Arg0, Args...> // One of Args is a lazy expression
, boost::mpl::bool_<(sizeof...(Args)!= evaluator_t::arity -1 ) && (evaluator_t::arity !=-1)> // if -1 : no check
>, // end of OR
@ -188,18 +187,8 @@ namespace triqs { namespace gfs {
operator()(Arg0 arg0, Args... args) const {
return clef::make_expr_call(view_type(*this),arg0, args...);
}
typedef typename std::result_of<data_proxy_t(data_t &,size_t)>::type r_type;
typedef typename std::result_of<data_proxy_t(data_t const &,size_t)>::type cr_type;
r_type operator() (mesh_point_t const & x) { return _data_proxy(_data, x.linear_index());}
cr_type operator() (mesh_point_t const & x) const { return _data_proxy(_data, x.linear_index());}
template<typename ... U>
r_type operator() (closest_pt_wrap<U...> const & p) { return _data_proxy(_data, _mesh.index_to_linear( gfs_implementation::get_closest_point<Variable,Target,Opt>::invoke(this,p)));}
template<typename ... U>
cr_type operator() (closest_pt_wrap<U...> const & p) const { return _data_proxy(_data, _mesh.index_to_linear( gfs_implementation::get_closest_point<Variable,Target,Opt>::invoke(this,p)));}
/*
// on mesh component for composite meshes
// enable iif the first arg is a mesh_point_t for the first component of the mesh_t
template<typename Arg0, typename ... Args, bool MeshIsComposite = std::is_base_of<tag::composite, mesh_t>::value >
@ -211,6 +200,30 @@ namespace triqs { namespace gfs {
typename std::enable_if< MeshIsComposite && std::is_base_of< tag::mesh_point, Arg0>::value, cr_type>::type
operator() (Arg0 const & arg0, Args const & ... args) const
{ return _data_proxy(_data, _mesh.mesh_pt_components_to_linear(arg0, args...));}
*/
//// [] and access to the grid point
typedef typename std::result_of<data_proxy_t(data_t &,size_t)>::type r_type;
typedef typename std::result_of<data_proxy_t(data_t const &,size_t)>::type cr_type;
r_type operator[] (mesh_index_t const & arg) { return _data_proxy(_data,_mesh.index_to_linear(arg));}
cr_type operator[] (mesh_index_t const & arg) const { return _data_proxy(_data,_mesh.index_to_linear(arg));}
r_type operator[] (mesh_point_t const & x) { return _data_proxy(_data, x.linear_index());}
cr_type operator[] (mesh_point_t const & x) const { return _data_proxy(_data, x.linear_index());}
template<typename ... U>
r_type operator[] (closest_pt_wrap<U...> const & p) { return _data_proxy(_data, _mesh.index_to_linear( gfs_implementation::get_closest_point<Variable,Target,Opt>::invoke(this,p)));}
template<typename ... U>
cr_type operator[] (closest_pt_wrap<U...> const & p) const { return _data_proxy(_data, _mesh.index_to_linear( gfs_implementation::get_closest_point<Variable,Target,Opt>::invoke(this,p)));}
// Interaction with the CLEF library : calling the gf with any clef expression as argument build a new clef expression
template<typename Arg>
typename boost::lazy_enable_if< // enable the template if
clef::is_any_lazy<Arg>, // One of Args is a lazy expression
clef::result_of::make_expr_subscript<view_type,Arg>
>::type // end of lazy_enable_if
operator[](Arg && arg) const { return clef::make_expr_subscript(view_type(*this),std::forward<Arg>(arg));}
/// A direct access to the grid point
template<typename... Args>
@ -231,18 +244,7 @@ namespace triqs { namespace gfs {
};
_on_mesh_wrapper_const friend on_mesh(gf_impl const & f) { return f;}
_on_mesh_wrapper friend on_mesh(gf_impl & f) { return f;}
public:
r_type operator[] (mesh_index_t const & arg) { return _data_proxy(_data,_mesh.index_to_linear(arg));}
cr_type operator[] (mesh_index_t const & arg) const { return _data_proxy(_data,_mesh.index_to_linear(arg));}
// Interaction with the CLEF library : calling the gf with any clef expression as argument build a new clef expression
template<typename Arg>
typename boost::lazy_enable_if< // enable the template if
clef::is_any_lazy<Arg>, // One of Args is a lazy expression
clef::result_of::make_expr_call<view_type,Arg>
>::type // end of lazy_enable_if
operator[](Arg && arg) const { return clef::make_expr_call(view_type(*this),std::forward<Arg>(arg));}
//----------------------------- HDF5 -----------------------------
@ -362,10 +364,11 @@ namespace triqs { namespace gfs {
private:
template<typename RHS> void triqs_clef_auto_assign_impl (RHS const & rhs, std::integral_constant<bool,false>) {
for (auto const & w: this->mesh()) (*this)(w) = rhs(w);
for (auto const & w: this->mesh()) (*this)[w] = rhs(w);
//for (auto const & w: this->mesh()) (*this)[w] = rhs(typename B::mesh_t::mesh_point_t::cast_t(w));
}
template<typename RHS> void triqs_clef_auto_assign_impl (RHS const & rhs, std::integral_constant<bool,true>) {
for (auto const & w: this->mesh()) (*this)(w) = triqs::tuple::apply(rhs,w.components_tuple());
for (auto const & w: this->mesh()) (*this)[w] = triqs::tuple::apply(rhs,w.components_tuple());
//for (auto w: this->mesh()) triqs::tuple::apply(*this,w.components_tuple()) = triqs::tuple::apply(rhs,w.components_tuple());
}
@ -389,14 +392,14 @@ namespace triqs { namespace gfs {
template<typename Tag, typename D, typename Target = matrix_valued> struct gf_keeper{ gf_view<D,Target> g; gf_keeper (gf_view<D,Target> const & g_) : g(g_) {} };
// ---------------------------------- slicing ------------------------------------
//slice
template<typename Variable, typename Target, typename Opt, bool V, typename... Args>
gf_view<Variable,matrix_valued,Opt> slice_target (gf_impl<Variable,Target,Opt,V> const & g, Args... args) {
static_assert(std::is_same<Target,matrix_valued>::value, "slice_target only for matrix_valued GF's");
using arrays::range;
//auto sg=slice_target (g.singularity(),range(args,args+1)...);
return gf_view<Variable,matrix_valued,Opt>(g.mesh(), g.data()(tqa::range(), args... ), slice_target (g.singularity(),args...) , g.symmetry());
return gf_view<Variable,matrix_valued,Opt>(g.mesh(), g.data()(range(), args... ), slice_target (g.singularity(),args...) , g.symmetry());
}
template<typename Variable, typename Target, typename Opt, bool V, typename... Args>
@ -404,14 +407,14 @@ namespace triqs { namespace gfs {
static_assert(std::is_same<Target,matrix_valued>::value, "slice_target only for matrix_valued GF's");
using arrays::range;
auto sg=slice_target (g.singularity(),range(args,args+1)...);
return gf_view<Variable,scalar_valued,Opt>(g.mesh(), g.data()(tqa::range(), args... ), sg , g.symmetry());
return gf_view<Variable,scalar_valued,Opt>(g.mesh(), g.data()(range(), args... ), sg , g.symmetry());
}
/*
template<typename Variable1,typename Variable2, typename Target, typename Opt, bool V, typename... Args>
gf_view<Variable2,Target,Opt> slice_mesh (gf_impl<Variable1,Target,Opt,V> const & g, Args... args) {
return gf_view<Variable2,Target,Opt>(g.mesh().slice(args...), g.data()(g.mesh().slice_get_range(args...),arrays::ellipsis()), g.singularity(), g.symmetry());
}*/
/*
template<typename Variable1,typename Variable2, typename Target, typename Opt, bool V, typename... Args>
gf_view<Variable2,Target,Opt> slice_mesh (gf_impl<Variable1,Target,Opt,V> const & g, Args... args) {
return gf_view<Variable2,Target,Opt>(g.mesh().slice(args...), g.data()(g.mesh().slice_get_range(args...),arrays::ellipsis()), g.singularity(), g.symmetry());
}*/
}}

3
triqs/gf/imfreq.hpp
View File

@ -57,6 +57,9 @@ namespace triqs { namespace gfs {
typedef typename std::conditional < std::is_same<Target, matrix_valued>::value, arrays::matrix_view<std::complex<double>>, std::complex<double>>::type rtype;
template<typename G>
rtype operator() (G const * g, long n) const {return g->data()(n, arrays::ellipsis()); }
// crucial because the mesh_point is cast in a complex, not an int !
template<typename G>
rtype operator() (G const * g, linear_mesh<matsubara_domain<true>>::mesh_point_t const & p) const { return (*this)(g,p.index());}
template<typename G>
local::tail_view operator()(G const * g, freq_infty const &) const {return g->singularity();}
};

12
triqs/gf/local/fourier_matsubara.cpp
View File

@ -68,15 +68,15 @@ namespace triqs { namespace gfs {
}
if (gw.domain().statistic == Fermion){
for (auto & t : gt.mesh())
g_in(t.index()) = fact * exp(iomega*t) * ( gt(t) - ( oneFermion(a1,b1,t,beta) + oneFermion(a2,b2,t,beta)+ oneFermion(a3,b3,t,beta) ) );
g_in[t.index()] = fact * exp(iomega*t) * ( gt[t] - ( oneFermion(a1,b1,t,beta) + oneFermion(a2,b2,t,beta)+ oneFermion(a3,b3,t,beta) ) );
}
else {
for (auto & t : gt.mesh())
g_in(t.index()) = fact * ( gt(t) - ( oneBoson(a1,b1,t,beta) + oneBoson(a2,b2,t,beta) + oneBoson(a3,b3,t,beta) ) );
g_in[t.index()] = fact * ( gt[t] - ( oneBoson(a1,b1,t,beta) + oneBoson(a2,b2,t,beta) + oneBoson(a3,b3,t,beta) ) );
}
details::fourier_base(g_in, g_out, L, true);
for (auto & w : gw.mesh()) {
gw(w) = g_out( w.index() ) * exp(iomega2*w.index() ) + a1/(w-b1) + a2/(w-b2) + a3/(w-b3);
gw[w] = g_out( w.index() ) * exp(iomega2*w.index() ) + a1/(w-b1) + a2/(w-b2) + a3/(w-b3);
}
gw.singularity() = gt.singularity();// set tail
}
@ -110,7 +110,7 @@ namespace triqs { namespace gfs {
}
g_in() = 0;
for (auto & w: gw.mesh()) {
g_in( w.index() ) = fact * exp(w.index()*iomega2) * ( gw(w) - (a1/(w-b1) + a2/(w-b2) + a3/(w-b3)) );
g_in[ w.index() ] = fact * exp(w.index()*iomega2) * ( gw[w] - (a1/(w-b1) + a2/(w-b2) + a3/(w-b3)) );
}
// for bosons GF(w=0) is divided by 2 to avoid counting it twice
if (gw.domain().statistic == Boson && !Green_Function_Are_Complex_in_time ) g_in(0) *= 0.5;
@ -122,13 +122,13 @@ namespace triqs { namespace gfs {
//typedef typename gf<imtime>::mesh_type::gf_result_type gt_result_type;
if (gw.domain().statistic == Fermion){
for (auto & t : gt.mesh()){
gt(t) = convert_green<gt_result_type> ( g_out( t.index() == L ? 0 : t.index() ) * exp(-iomega*t)
gt[t] = convert_green<gt_result_type> ( g_out( t.index() == L ? 0 : t.index() ) * exp(-iomega*t)
+ oneFermion(a1,b1,t,beta) + oneFermion(a2,b2,t,beta)+ oneFermion(a3,b3,t,beta) );
}
}
else {
for (auto & t : gt.mesh())
gt(t) = convert_green<gt_result_type> ( g_out( t.index() == L ? 0 : t.index() )
gt[t] = convert_green<gt_result_type> ( g_out( t.index() == L ? 0 : t.index() )
+ oneBoson(a1,b1,t,beta) + oneBoson(a2,b2,t,beta) + oneBoson(a3,b3,t,beta) );
}
if (gt.mesh().kind() == full_bins)

10
triqs/gf/local/fourier_real.cpp
View File

@ -58,12 +58,12 @@ namespace triqs { namespace gfs {
dcomplex a1 = (t1 + I * t2/a )/2., a2 = (t1 - I * t2/a )/2.;
for (auto const & t : gt.mesh())
g_in(t.index()) = (gt(t) - (a1*th_expo(t,a) + a2*th_expo_neg(t,a))) * std::exp(I*t*wmin);
g_in[t.index()] = (gt[t] - (a1*th_expo(t,a) + a2*th_expo_neg(t,a))) * std::exp(I*t*wmin);
details::fourier_base(g_in, g_out, L, true);
for (auto const & w : gw.mesh())
gw(w) = gt.mesh().delta() * std::exp(I*(w-wmin)*tmin) * g_out(w.index())
gw[w] = gt.mesh().delta() * std::exp(I*(w-wmin)*tmin) * g_out(w.index())
+ a1*th_expo_inv(w,a) + a2*th_expo_neg_inv(w,a);
gw.singularity() = gt.singularity();// set tail
@ -90,14 +90,14 @@ namespace triqs { namespace gfs {
g_in() = 0;
for (auto const & w: gw.mesh())
g_in(w.index()) = (gw(w) - a1*th_expo_inv(w,a) - a2*th_expo_neg_inv(w,a) ) * std::exp(-I*w*tmin);
g_in(w.index()) = (gw[w] - a1*th_expo_inv(w,a) - a2*th_expo_neg_inv(w,a) ) * std::exp(-I*w*tmin);
details::fourier_base(g_in, g_out, L, false);
const double corr = 1.0/(gt.mesh().delta()*L);
for (auto const & t : gt.mesh())
gt(t) = corr * std::exp(I*wmin*(tmin-t)) *
g_out( t.index() ) + a1 * th_expo(t,a) + a2 * th_expo_neg(t,a) ;
gt[t] = corr * std::exp(I*wmin*(tmin-t)) *
g_out[ t.index() ] + a1 * th_expo(t,a) + a2 * th_expo_neg(t,a) ;
// set tail
gt.singularity() = gw.singularity();

8
triqs/gf/local/functions.cpp
View File

@ -44,7 +44,7 @@ namespace triqs { namespace gfs {
dcomplex d= t(1)(n1,n2) , A=t(2)(n1,n2),B = t(3)(n1,n2) ;
double b1 = 0,b2 =1, b3 =-1;
dcomplex a1 = d-B, a2 = (A+B)/2, a3 = (B-A)/2;
for (auto & w : G.mesh()) dens_part(n1,n2)+= G(w)(n1,n2) - (a1/(w - b1) + a2 / (w-b2) + a3/(w-b3));
for (auto & w : G.mesh()) dens_part(n1,n2)+= G[w](n1,n2) - (a1/(w - b1) + a2 / (w-b2) + a3/(w-b3));
dens_part(n1,n2) = dens_part(n1,n2)/Beta;
dens_tail(n1,n2) = d + F(a1,b1,Beta) + F(a2,b2,Beta)+ F(a3,b3,Beta);
// If the Green function are NOT complex, then one use the symmetry property
@ -77,7 +77,7 @@ namespace triqs { namespace gfs {
res() = 0.0;
for (auto l : gl.mesh()) {
res -= sqrt(2*l.index()+1) * gl(l);
res -= sqrt(2*l.index()+1) * gl[l];
}
res /= gl.domain().beta;
@ -95,7 +95,7 @@ namespace triqs { namespace gfs {
for (int p=1; p<=t.order_max(); p++)
for (auto l : gl.mesh())
t(p) += (triqs::utility::legendre_t(l.index(),p)/pow(gl.domain().beta,p)) * gl(l);
t(p) += (triqs::utility::legendre_t(l.index(),p)/pow(gl.domain().beta,p)) * gl[l];
return t;
@ -113,7 +113,7 @@ namespace triqs { namespace gfs {
tqa::array<double,2> corr(disc.shape()); corr() = 0;
for (auto l : gl.mesh()) {
corr += t(l.index()) * gl(l);
corr += t(l.index()) * gl[l];
}
tqa::range R;

6
triqs/gf/local/legendre_matsubara.cpp
View File

@ -37,7 +37,7 @@ void legendre_matsubara_direct(gf_view<imfreq> & gw, gf_view<legendre> const & g
// Use the transformation matrix
for (auto om: gw.mesh()) {
for (auto l: gl.mesh()) {
gw(om) += legendre_T(om.index(),l.index()) * gl(l);
gw[om] += legendre_T(om.index(),l.index()) * gl[l];
}
}
@ -72,7 +72,7 @@ void legendre_matsubara_direct (gf_view<imtime> & gt, gf_view<legendre> const &
for (auto t : gt.mesh()) {
L.reset( 2*t/gt.domain().beta-1 );
for (auto l : gl.mesh()) {
gt(t) += sqrt(2*l.index()+1) / gt.domain().beta * gl(l) * L.next();
gt[t] += sqrt(2*l.index()+1) / gt.domain().beta * gl[l] * L.next();
}
}
@ -89,7 +89,7 @@ void legendre_matsubara_inverse (gf_view<legendre> & gl, gf_view<imtime> const &
for (auto t : gt.mesh()) {
L.reset( 2*t/gt.domain().beta-1 );
for (auto l : gl.mesh()) {
gl(l) += sqrt(2*l.index()+1) * L.next() * gt(t);
gl[l] += sqrt(2*l.index()+1) * L.next() * gt[t];
}
}
gl.data() *= gt.mesh().delta();

2
triqs/gf/local/pade.cpp
View File

@ -53,7 +53,7 @@ namespace triqs { namespace gfs {
gr() = 0.0;
for (auto om : gr.mesh()) {
dcomplex e = om + dcomplex(0.0,1.0)*freq_offset;
gr(om)(n1,n2) = PA(e);
gr[om](n1,n2) = PA(e);
}
}

3
triqs/gf/meshes/discrete.hpp
View File

@ -49,7 +49,8 @@ namespace triqs { namespace gfs {
public:
mesh_point_t( discrete_mesh const & mesh, index_t const & index_): m(&mesh), _index(index_) {}
void advance() { ++_index;}
operator size_t () const { return m->index_to_point(_index);}
typedef size_t cast_t;
operator cast_t() const { return m->index_to_point(_index);}
size_t linear_index() const { return _index;}
size_t index() const { return _index;}
bool at_end() const { return (_index == m->size());}

5
triqs/gf/meshes/linear.hpp
View File

@ -78,11 +78,12 @@ namespace triqs { namespace gfs {
/// The wrapper for the mesh point
class mesh_point_t : tag::mesh_point, public arith_ops_by_cast<mesh_point_t, domain_pt_t > {
linear_mesh const * m;
index_t _index;
index_t _index;
public:
mesh_point_t( linear_mesh const & mesh, index_t const & index_): m(&mesh), _index(index_) {}
void advance() { ++_index;}
operator domain_pt_t () const { return m->index_to_point(_index);}
typedef domain_pt_t cast_t;
operator cast_t () const { return m->index_to_point(_index);}
size_t linear_index() const { return _index;}
size_t index() const { return _index;}
bool at_end() const { return (_index == m->size());}

3
triqs/gfs.hpp
View File

@ -22,11 +22,12 @@
#define TRIQS_GFS_ALL_H
// The basic classes
#include <triqs/gf/block.hpp>
#include <triqs/gf/imtime.hpp>
#include <triqs/gf/imfreq.hpp>
#include <triqs/gf/retime.hpp>
#include <triqs/gf/refreq.hpp>
#include <triqs/gf/block.hpp>
#include <triqs/gf/legendre.hpp>
//#include <triqs/gf/product.hpp>