mirror of
https://github.com/triqs/dft_tools
synced 2024-10-31 19:23:45 +01:00
101 lines
3.3 KiB
C++
101 lines
3.3 KiB
C++
#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
|
|
|
|
#include <triqs/gfs/refreq.hpp>
|
|
#include <triqs/gfs/retime.hpp>
|
|
#include <triqs/gfs/local/fourier_real.hpp>
|
|
#include <triqs/arrays.hpp>
|
|
|
|
using triqs::arrays::make_shape;
|
|
using triqs::gfs::refreq;
|
|
using triqs::gfs::retime;
|
|
using triqs::gfs::make_gf;
|
|
|
|
double lorentzian(double w, double a){
|
|
return 2*a / (w*w + a*a) ;
|
|
};
|
|
std::complex<double> lorentzian_inverse(double t, double a){
|
|
return std::exp(-a*std::abs(t)) ;
|
|
};
|
|
double theta(double x){
|
|
return x>0 ? 1.0 : ( x<0 ? 0.0 : 0.5 ) ;
|
|
};
|
|
|
|
int main() {
|
|
|
|
|
|
double precision=10e-10;
|
|
H5::H5File file("fourier_real_time.h5",H5F_ACC_TRUNC);
|
|
|
|
std::complex<double> I(0,1);
|
|
|
|
//Test on the tail: GF in frequency that is a lorentzian, with its singularity, TF and TF^-1.
|
|
|
|
double wmax=10;
|
|
int Nw=1001;
|
|
|
|
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);
|
|
Gw1.singularity()(2)=triqs::arrays::matrix<double>{{2.0*a}};
|
|
h5_write(file,"Gw1",Gw1); // the original lorentzian
|
|
|
|
auto Gt1 = inverse_fourier(Gw1);
|
|
h5_write(file,"Gt1",Gt1); // the lorentzian TF : lorentzian_inverse
|
|
|
|
// 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";
|
|
}
|
|
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";
|
|
}
|
|
h5_write(file,"Gt1b",Gt1); // must be 0
|
|
|
|
|
|
|
|
//Test on the tail: GF in time that is a decreasing exponential
|
|
|
|
double tmax=10.;
|
|
int Nt=501;
|
|
|
|
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) );
|
|
Gt2.singularity()(1)=triqs::arrays::matrix<double>{{1.0}};
|
|
h5_write(file,"Gt2",Gt2);
|
|
|
|
auto Gw2 = fourier(Gt2);
|
|
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";
|
|
}
|
|
h5_write(file,"Gw2b",Gw2);
|
|
|
|
|
|
|
|
//Test : GF in time is a simple trigonometric function, the result is a sum of Dirac functions
|
|
|
|
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.)) ;
|
|
h5_write(file,"Gt3",Gt3);
|
|
|
|
auto Gw3 = fourier(Gt3);
|
|
h5_write(file,"Gw3",Gw3);
|
|
|
|
}
|
|
|
|
|