mirror of
https://github.com/triqs/dft_tools
synced 2024-12-26 06:14:14 +01:00
88f8e4cce4
I also added a test to make sure the time mesh is twice as long as the frequency mesh. Obviously now some tests don't pass... I will fix them in the next commit.
245 lines
9.3 KiB
C++
245 lines
9.3 KiB
C++
/*******************************************************************************
|
|
*
|
|
* TRIQS: a Toolbox for Research in Interacting Quantum Systems
|
|
*
|
|
* Copyright (C) 2011 by M. Ferrero, 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 <http://www.gnu.org/licenses/>.
|
|
*
|
|
******************************************************************************/
|
|
#include "fourier_base.hpp"
|
|
#include "fourier_matsubara.hpp"
|
|
#include <fftw3.h>
|
|
|
|
namespace triqs {
|
|
namespace gfs {
|
|
|
|
template <typename GfElementType> GfElementType convert_green(dcomplex const& x) { return x; }
|
|
template <> double convert_green<double>(dcomplex const& x) { return real(x); }
|
|
|
|
//--------------------------------------------------------------------------------------
|
|
|
|
struct impl_worker {
|
|
|
|
arrays::vector<dcomplex> g_in, g_out;
|
|
|
|
dcomplex oneFermion(dcomplex a, double b, double tau, double beta) {
|
|
return -a * (b >= 0 ? exp(-b * tau) / (1 + exp(-beta * b)) : exp(b * (beta - tau)) / (1 + exp(beta * b)));
|
|
}
|
|
|
|
dcomplex oneBoson(dcomplex a, double b, double tau, double beta) {
|
|
return a * (b >= 0 ? exp(-b * tau) / (exp(-beta * b) - 1) : exp(b * (beta - tau)) / (1 - exp(b * beta)));
|
|
}
|
|
|
|
//-------------------------------------
|
|
|
|
void direct(gf_view<imfreq, scalar_valued> gw, gf_const_view<imtime, scalar_valued> gt) {
|
|
auto ta = gt(freq_infty());
|
|
direct_impl(make_gf_view_without_tail(gw), make_gf_view_without_tail(gt), ta);
|
|
gw.singularity() = gt.singularity(); // set tail
|
|
}
|
|
|
|
void direct(gf_view<imfreq, scalar_valued, no_tail> gw, gf_const_view<imtime, scalar_valued, no_tail> gt) {
|
|
auto ta = local::tail{1,1};
|
|
direct_impl(gw, gt, ta);
|
|
}
|
|
|
|
//-------------------------------------
|
|
|
|
private:
|
|
void direct_impl(gf_view<imfreq, scalar_valued, no_tail> gw, gf_const_view<imtime, scalar_valued, no_tail> gt,
|
|
local::tail const& ta) {
|
|
// TO BE MODIFIED AFTER SCALAR IMPLEMENTATION TODO
|
|
dcomplex d = ta(1)(0, 0), A = ta.get_or_zero(2)(0, 0), B = ta.get_or_zero(3)(0, 0);
|
|
double b1 = 0, b2 = 0, b3 = 0;
|
|
dcomplex a1, a2, a3;
|
|
double beta = gt.mesh().domain().beta;
|
|
auto L = (gt.mesh().kind() == full_bins ? gt.mesh().size() - 1 : gt.mesh().size());
|
|
if (L < 2*gw.mesh().size()) TRIQS_RUNTIME_ERROR << "The time mesh mush be at least twice as long as the freq mesh";
|
|
double fact = beta / L;
|
|
dcomplex iomega = dcomplex(0.0, 1.0) * std::acos(-1) / beta;
|
|
dcomplex iomega2 = iomega * 2 * gt.mesh().delta() * (gt.mesh().kind() == half_bins ? 0.5 : 0.0);
|
|
g_in.resize(L);
|
|
g_out.resize(gw.mesh().size());
|
|
if (gw.domain().statistic == Fermion) {
|
|
b1 = 0;
|
|
b2 = 1;
|
|
b3 = -1;
|
|
a1 = d - B;
|
|
a2 = (A + B) / 2;
|
|
a3 = (B - A) / 2;
|
|
} else {
|
|
b1 = -0.5;
|
|
b2 = -1;
|
|
b3 = 1;
|
|
a1 = 4 * (d - B) / 3;
|
|
a2 = B - (d + A) / 2;
|
|
a3 = d / 6 + A / 2 + B / 3;
|
|
}
|
|
if (gw.domain().statistic == Fermion) {
|
|
for (auto& t : gt.mesh())
|
|
if(t.index() < L) {
|
|
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())
|
|
if(t.index() < L) {
|
|
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);
|
|
}
|
|
}
|
|
|
|
public:
|
|
//-------------------------------------
|
|
|
|
void inverse(gf_view<imtime, scalar_valued> gt, gf_const_view<imfreq, scalar_valued> gw) {
|
|
static bool Green_Function_Are_Complex_in_time = false;
|
|
// If the Green function are NOT complex, then one use the symmetry property
|
|
// fold the sum and get a factor 2
|
|
auto ta = gw(freq_infty());
|
|
// TO BE MODIFIED AFTER SCALAR IMPLEMENTATION TODO
|
|
dcomplex d = ta(1)(0, 0), A = ta.get_or_zero(2)(0, 0), B = ta.get_or_zero(3)(0, 0);
|
|
double b1, b2, b3;
|
|
dcomplex a1, a2, a3;
|
|
|
|
double beta = gw.domain().beta;
|
|
size_t L = gt.mesh().size() - (gt.mesh().kind() == full_bins ? 1 : 0); // L can be different from gt.mesh().size() (depending
|
|
// on the mesh kind) and is given to the FFT algorithm
|
|
if (L < 2*gw.mesh().size()) TRIQS_RUNTIME_ERROR << "The time mesh mush be at least twice as long as the freq mesh";
|
|
dcomplex iomega = dcomplex(0.0, 1.0) * std::acos(-1) / beta;
|
|
dcomplex iomega2 = -iomega * 2 * gt.mesh().delta() * (gt.mesh().kind() == half_bins ? 0.5 : 0.0);
|
|
double fact = (Green_Function_Are_Complex_in_time ? 1 : 2) / beta;
|
|
g_in.resize(gw.mesh().size());
|
|
g_out.resize(L);
|
|
|
|
if (gw.domain().statistic == Fermion) {
|
|
b1 = 0;
|
|
b2 = 1;
|
|
b3 = -1;
|
|
a1 = d - B;
|
|
a2 = (A + B) / 2;
|
|
a3 = (B - A) / 2;
|
|
} else {
|
|
b1 = -0.5;
|
|
b2 = -1;
|
|
b3 = 1;
|
|
a1 = 4 * (d - B) / 3;
|
|
a2 = B - (d + A) / 2;
|
|
a3 = d / 6 + A / 2 + B / 3;
|
|
}
|
|
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)));
|
|
}
|
|
// 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;
|
|
|
|
details::fourier_base(g_in, g_out, L, false);
|
|
|
|
// CORRECT FOR COMPLEX G(tau) !!!
|
|
typedef double gt_result_type;
|
|
// typedef typename gf<imtime>::mesh_type::gf_result_type gt_result_type;
|
|
if (gw.domain().statistic == Fermion) {
|
|
for (auto& t : gt.mesh()) {
|
|
if (t.index() < L) {
|
|
gt[t] =
|
|
convert_green<gt_result_type>(g_out(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())
|
|
if (t.index() < L) {
|
|
gt[t] = convert_green<gt_result_type>(g_out(t.index()) + oneBoson(a1, b1, t, beta) +
|
|
oneBoson(a2, b2, t, beta) + oneBoson(a3, b3, t, beta));
|
|
}
|
|
}
|
|
double pm = (gw.domain().statistic == Fermion ? -1.0 : 1.0);
|
|
if (gt.mesh().kind() == full_bins) gt.on_mesh(L) = pm * (gt.on_mesh(0) + convert_green<gt_result_type>(ta(1)(0, 0)));
|
|
// set tail
|
|
gt.singularity() = gw.singularity();
|
|
}
|
|
|
|
}; // class worker
|
|
|
|
//--------------------------------------------
|
|
|
|
template <typename Opt>
|
|
void fourier_impl(gf_view<imfreq, scalar_valued, Opt> gw, gf_const_view<imtime, scalar_valued, Opt> gt) {
|
|
impl_worker w;
|
|
w.direct(gw, gt);
|
|
}
|
|
|
|
template <typename Opt>
|
|
void fourier_impl(gf_view<imfreq, matrix_valued, Opt> gw, gf_const_view<imtime, matrix_valued, Opt> gt) {
|
|
impl_worker w;
|
|
for (size_t n1 = 0; n1 < gt.data().shape()[1]; n1++)
|
|
for (size_t n2 = 0; n2 < gt.data().shape()[2]; n2++) {
|
|
auto gw_sl = slice_target_to_scalar(gw, n1, n2);
|
|
auto gt_sl = slice_target_to_scalar(gt, n1, n2);
|
|
w.direct(gw_sl, gt_sl);
|
|
}
|
|
}
|
|
|
|
//---------------------------------------------------------------------------
|
|
|
|
void inverse_fourier_impl(gf_view<imtime, scalar_valued> gt, gf_const_view<imfreq, scalar_valued> gw) {
|
|
impl_worker w;
|
|
w.inverse(gt, gw);
|
|
}
|
|
|
|
void inverse_fourier_impl(gf_view<imtime, matrix_valued> gt, gf_const_view<imfreq, matrix_valued> gw) {
|
|
impl_worker w;
|
|
for (size_t n1 = 0; n1 < gw.data().shape()[1]; n1++)
|
|
for (size_t n2 = 0; n2 < gw.data().shape()[2]; n2++) {
|
|
auto gt_sl = slice_target_to_scalar(gt, n1, n2);
|
|
auto gw_sl = slice_target_to_scalar(gw, n1, n2);
|
|
w.inverse(gt_sl, gw_sl);
|
|
}
|
|
}
|
|
|
|
//---------------------------------------------------------------------------
|
|
void triqs_gf_view_assign_delegation(gf_view<imfreq, scalar_valued> g,
|
|
gf_keeper<tags::fourier, imtime, scalar_valued> const& L) {
|
|
fourier_impl(g, L.g);
|
|
}
|
|
void triqs_gf_view_assign_delegation(gf_view<imfreq, matrix_valued> g,
|
|
gf_keeper<tags::fourier, imtime, matrix_valued> const& L) {
|
|
fourier_impl(g, L.g);
|
|
}
|
|
void triqs_gf_view_assign_delegation(gf_view<imtime, scalar_valued> g,
|
|
gf_keeper<tags::fourier, imfreq, scalar_valued> const& L) {
|
|
inverse_fourier_impl(g, L.g);
|
|
}
|
|
void triqs_gf_view_assign_delegation(gf_view<imtime, matrix_valued> g,
|
|
gf_keeper<tags::fourier, imfreq, matrix_valued> const& L) {
|
|
inverse_fourier_impl(g, L.g);
|
|
}
|
|
void triqs_gf_view_assign_delegation(gf_view<imfreq, scalar_valued, no_tail> g,
|
|
gf_keeper<tags::fourier, imtime, scalar_valued, no_tail> const& L) {
|
|
fourier_impl(g, L.g);
|
|
}
|
|
void triqs_gf_view_assign_delegation(gf_view<imfreq, matrix_valued, no_tail> g,
|
|
gf_keeper<tags::fourier, imtime, matrix_valued, no_tail> const& L) {
|
|
fourier_impl(g, L.g);
|
|
}
|
|
}
|
|
}
|
|
|