LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-12-25 13:53:40 +01:00
Code Issues Releases Wiki Activity

gf. Clean Fourier

Browse Source 
- lazy_fourier and co --> fourier
- ex fourier --> make_gf_from_fourier to make a new gf
- = fourier (g) works only iif lhs is a view, like scalar.
- updated python (commented fourier method).
This commit is contained in:
Olivier Parcollet 2013-10-23 15:51:51 +02:00
parent 1d5ea99d4f
commit 579368f24b
29 changed files with 728 additions and 492 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

48
doc/reference/c++/gf/clef.rst Normal file
View File

@ -0,0 +1,48 @@
.. highlight:: c
.. _gf_clef:
Interaction with CLEF expressions
============================================
* The gf containers and their view classes can be used with the :doc:`../clef/contents` library :
* They can be called with CLEF expressions.
* :doc:`Automatic assignment<../clef/assign>` has been set up.
* Using the CLEF library offers a quick and efficient way to fill an array with multiple advantages :
* It is simpler and more readeable than a series of for loops.
* It is usually more optimal since the for loops are automatically written in the TraversalOrder of the array.
* **Example** :
.. compileblock::
#include <triqs/gfs.hpp>
using namespace triqs::gfs; using triqs::clef::placeholder;
int main(){
// Cf gf<imfreq> specialisation page for the constructor
double beta=10; int Nfreq =100;
auto g = gf<imfreq> { {beta,Fermion,Nfreq}, {1,1} };
// Filling the gf with something...
placeholder<0> wn_;
g(wn_) << 1/ (wn_ + 2);
g(wn_) << 1/ (wn_ + 2 + g(wn_) );
}
.. note::
The syntax uses a <<, not = since the array is not assigned to an expression
but filled by the evaluation thereof.
The LHS uses () and not brackets, even though it is on the mesh, because of the strange C++ limitation
that [] can not be overloaded for several variables...

175
doc/reference/c++/gf/fourier.rst
View File

@ -3,151 +3,78 @@
Fourier transforms Fourier transforms
################### ###################
The Fourier transforms from real and imaginary frequencies to times, and inverse, are currently implemented,
along with the analogous transformation from the Legendre expansion to imaginary time and frequencies.
Synopsis and example
======================
**Synopsis** ::
lazy_object fourier (gf<imfreq,Target,Opt> const &)
lazy_object fourier (gf_const_view<imfreq,Target,Opt> const &)
lazy_object inverse_fourier (gf<imfreq,Target,Opt> const &)
lazy_object inverse_fourier (gf_const_view<imfreq,Target,Opt> const &)
The fourier/inverse_fourier functions do **not** perform the Fourier transformation,
but returns a small lazy object (basically saying "Fourier Transform of XXX"),
which is then used in an assignment of a *view* of a gf.
Example
.. compileblock::
#include <triqs/gfs.hpp>
using namespace triqs::gfs;
int main() {
double beta =1, a=1;
int N=10000;
auto gw = gf<imfreq> {{beta, Fermion, N}, {1,1}};
auto gt = gf<imtime> {{beta, Fermion, N}, {1,1}};
triqs::clef::placeholder<0> om_;
gw (om_) << 1/(om_-a);
gt() = inverse_fourier(gw); // fills the *View* with the contents of the FFT.
// NB : the mesh must have the same size.
// make a new fresh gf, with the same size mesh, from the FFT of gt
auto gw2 = make_gf_from_fourier(gt);
}
Note that :
* the LHS of the = must be a view, since the RHS can not compute the domain of the function
(how many points to use ?).
* In the make_gf_from_fourier function, choice is explicitly made to generate a new gf with the same number of points in the mesh.
Convention Convention
========== ===========
For real time/frequency: For real time/frequency:
:label: _TF_R
.. math:: \tilde G(\omega)=\int_{-\infty}^\infty dt G(t)e^{i\omega t} .. math:: \tilde G(\omega)=\int_{-\infty}^\infty dt G(t)e^{i\omega t}
:label: _inv_TF_R
.. math:: G(t)=\int_{-\infty}^\infty \frac{d\omega}{2\pi} \tilde G(\omega)e^{-i\omega t} .. math:: G(t)=\int_{-\infty}^\infty \frac{d\omega}{2\pi} \tilde G(\omega)e^{-i\omega t}
For Matsubara (imaginary) time/frequency: For Matsubara (imaginary) time/frequency:
:label: _TF_I
.. math:: \tilde G(i\omega_n)=\int_{0}^\beta d\tau G(t)e^{i\omega_n \tau} .. math:: \tilde G(i\omega_n)=\int_{0}^\beta d\tau G(t)e^{i\omega_n \tau}
:label: _inv_TF_I
.. math:: G(\tau)=\sum_{n=-\infty}^\infty \frac{1}{\beta} \tilde G(i\omega_n)e^{-i\omega_n \tau} .. math:: G(\tau)=\sum_{n=-\infty}^\infty \frac{1}{\beta} \tilde G(i\omega_n)e^{-i\omega_n \tau}
The :math:`\omega_n`'s are :math:`\frac{(2n+1)\pi}{\beta}` for fermions, :math:`\frac{2n\pi}{\beta}` for bosons (as :math:`G(\tau+\beta)=-G(\tau)` for fermions, :math:`G(\tau)` for bosons). The :math:`\omega_n`'s are :math:`\frac{(2n+1)\pi}{\beta}` for fermions, :math:`\frac{2n\pi}{\beta}` for bosons (as :math:`G(\tau+\beta)=-G(\tau)` for fermions, :math:`G(\tau)` for bosons).
*
.. toctree::
The FFTW library :maxdepth: 1
================
Documentation on FFTW is on https://www.fftw.org.
FFTW is a C subroutine library for computing the discrete Fourier transform (DFT) in one or more dimensions, of arbitrary input size, and of both real and complex data.
It will be used to calculate the (inverse) Fourier transform, in real/imaginary time/frequency, using the fact that the GF values are stored for a finite amount of regularly spaced values.
The DFT transforms of a sequence of :math:`N` complex numbers :math:`f_0...,f_{N-1}` into a sequence of :math:`N` complex numbers :math:`\tilde f_0...,\tilde f_{N-1}` according to the formula:
:label: _DFT
.. math:: \tilde f_k = \sum_{n=0}^{N-1} f_n e^{-i 2 \pi k n / N}.
The inverse DFT formula is
:label: _inv_DFT
.. math:: f_n = \frac{1}{N} \sum_{k=0}^{N-1} \tilde f_k e^{i 2 \pi k n / N}.
Implementation in real time/frequency using FFTW
================================================
The real time mesh parameters are :math:`t_{min}`, :math:`\delta t` and :math:`N_t`.
For the real frequency mesh, they are :math:`\omega_{min}`, :math:`\delta \omega` and :math:`N_\omega`.
The Fourier transform requires :math:`N_\omega=N_t` and :math:`\delta t \delta \omega= \frac{2\pi}{N_t}`.
The times are :math:`t_k=t_{min}+k\delta t` and the frequencies :math:`\omega_m=\omega_{min}+m\delta \omega`.
By approximating Eq. :ref:`TF_R` by
.. math:: \tilde G(\omega_m) = \delta t \sum_{k=0}^{N_t} G(t_k) e^{i\omega_m t_k},
we recognize an inverse DFT (Eq. :ref:`inv_DFT`). To calculate it using FFTW, we first need to prepare the input :math:`\tilde f_k`, then to do the DFT and finally to modify the output to obtain :math:`\tilde G(\omega_m)` using the two formulas:
.. math:: \tilde f_k = G(t_k) e^{i \omega_{min}t_k},
.. math:: \tilde G(\omega_m) = \delta t f_m e^{i t_{min}(\omega_m-\omega_{min})}.
Similarly, the inverse transformation is obtained by approximating Eq. :ref:`eq_inv_TF_R` by
.. math:: G(t_k)=\frac{\delta\omega}{2\pi}\sum_{m=0}^{N_\omega} \tilde G(\omega_m)e^{-i\omega_m t_k},
we recognize a DFT (Eq. :ref:`DFT`). To calculate it using FFTW, we first need to prepare the input :math:`f_m`, then to do the inverse DFT and finally to modify the output to obtain :math:`G(t_k)`:
.. math:: f_m = \tilde G(\omega_m) e^{-i t_{min}\omega_m},
.. math:: G(t_k) = \frac{1}{N_t \delta t}\tilde f_k e^{-i \omega_{min}(t_k-t_{min})}.
Implementation in imaginary time/frequency using FFTW
=====================================================
The imaginary time mesh parameters are :math:`\beta` and :math:`N_\tau`, plus a tag ``half_bins``, ``full_bins`` or ``without_last``.
In the ``full_bins`` case, one point of the time GF has to be removed for the fourier transform.
From these parameters, we deduce :math:`\delta\tau=\beta/N_\tau`.
For the imaginary frequency mesh, the mesh parameters are :math:`\beta`, :math:`n_{min}` and :math:`N_{\omega_n}`.
From them, we deduce :math:`\delta\omega=\frac{2\pi}{\beta}`.
The Fourier transform requires :math:`N_\omega=N_\tau`.
The times are :math:`\tau_k=\tau_{min}+k\delta\tau` and the frequencies :math:`\omega_n=\omega_{min}+n\delta \omega`.
:math:`\tau_{min}` is either 0 or :math:`\delta\tau/2` depending on the mesh kind.
:math:`\omega_{min}` is either :math:`\frac{2\pi(n_{min}+1)}{\beta}` or :math:`\frac{2\pi n_{min}}{\beta}` depending on the statistics.
We approximate the TF and its inverse by
.. math:: \tilde G(i\omega_n) = \delta\tau \sum_{k=0}^{N_\tau} G(\tau_k)e^{i\omega_n \tau_k}
.. math:: G(\tau_k) = \sum_{n=n_{min}}^{N_\tau} \frac{1}{\beta} \tilde G(i\omega_n)e^{-i\omega_n \tau_k}
We use for the TF:
.. math:: \tilde f_k = G(\tau_k) e^{i \omega_{min}\tau_k},
.. math:: \tilde G(i\omega_n) = \frac{\beta}{N_\tau} f_n e^{i \tau_{min}(\omega_n-\omega_{min})}.
and for the inverse TF:
.. math:: f_m = \frac{1}{\beta}\tilde G(i\omega_n) e^{-i t_{min}\omega_n},
.. math:: G(t_k) = \tilde f_k e^{-i \omega_{min}(\tau_k-\tau_{min})},
Special case of real functions in time for fermions
----------------------------------------------------
In this case, :math:`G(i\omega_n)=conj(G(i\omega_n))` and we only store the values of :math:`G(i\omega_n)` for :math:`\omega_n > 0`.
The Eq. :ref:`inv_DFT_I` becomes:
:label: _inv_TF_I_real_fermion
.. math:: G(\tau)=\sum_{n=0}^\infty \frac{2}{\beta} \tilde G(i\omega_n)\cos(\omega_n \tau)
The inverse TF formulas are in this case
.. math:: f_m = \frac{2}{\beta}\tilde G(i\omega_n) e^{-i t_{min}\omega_n},
.. math:: G(t_k) = \tilde f_k e^{-i \omega_{min}(\tau_k-\tau_{min})},
Special case of real functions in time for bosons
--------------------------------------------------
In this case, :math:`G(i\omega_n)=conj(G(i\omega_n))` and we only store the values of :math:`G(i\omega_n)` for :math:`\omega_n \ge 0`.
The Eq. :ref:`inv_DFT_I` becomes:
:label: _inv_TF_I_real_bosons
.. math:: G(\tau)=\frac{1}{\beta} \tilde G(0)+\sum_{n=1}^\infty \frac{2}{\beta} \tilde G(i\omega_n)\cos(\omega_n \tau)
The inverse TF formulas are in this case
.. math:: f_0 = \frac{1}{\beta}\tilde G(0),
.. math:: f_m = \frac{2}{\beta}\tilde G(i\omega_n) \cos(t_{min}\omega_n),
.. math:: G(t_k) = \tilde f_k e^{-i \omega_{min}(\tau_k-\tau_{min})},
Effect of a TF on the tail
===========================
The tail is unchanged during a TF, but its value is used to limit the errors.
The components 1 and 2 of the tail (:math:`t_1` and :math:`t_2`) are used to improve the computation of the GF in the following way:
in the large :math:`\omega` limit,
.. math:: G(\omega)\simeq \frac{t_1}{\omega}+\frac{t_2}{\omega^2}\simeq \frac{a_1}{\omega+i}+\frac{a_2}{\omega-i}
with :math:`a_1=\frac{t_1+it_2}{2}` and :math:`a_2=\frac{t_1-it_2}{2}`.
As these large w terms are badly taken into account if we naively Fourier transform the function described by its values on the mesh in w, we substract them, do the Fourier transform and add their Fourier transform to the result.
We use the following Fourier tranforms:
.. math:: \frac{2a}{\omega^2+a^2} \leftrightarrow e^{-a|x|}
.. math:: \frac{1}{\omega+i} \leftrightarrow -i e^{-x} \theta(x)
.. math:: \frac{1}{\omega-i} \leftrightarrow i e^{x} \theta(-x)
For the inverse Fourier transform, the inverse procedure is used.
In the library, :math:`a` is optimized according to the mesh properties (its size :math:`L=G.mesh().size()` and its precision :math:`\delta = G.mesh().delta()`).
The requirements are :math:`a \gg \delta\omega` and :math:`a \ll L\delta\omega`, or equivalently :math:`a \gg \delta t` and :math:`a \ll L\delta t`.
Thus, we chose :math:`a=\sqrt{L}\delta\omega`
fourier_impl_notes

126
doc/reference/c++/gf/fourier_impl_notes.rst Normal file
View File

@ -0,0 +1,126 @@
.. highlight:: c
Fourier : implementation notes
#######################################
The FFTW library
---------------------
Documentation on FFTW is on https://www.fftw.org.
FFTW is a C subroutine library for computing the discrete Fourier transform (DFT) in one or more dimensions, of arbitrary input size, and of both real and complex data.
It will be used to calculate the (inverse) Fourier transform, in real/imaginary time/frequency, using the fact that the GF values are stored for a finite amount of regularly spaced values.
The DFT transforms of a sequence of :math:`N` complex numbers :math:`f_0...,f_{N-1}` into a sequence of :math:`N` complex numbers :math:`\tilde f_0...,\tilde f_{N-1}` according to the formula:
.. math:: \tilde f_k = \sum_{n=0}^{N-1} f_n e^{-i 2 \pi k n / N}.
The inverse DFT formula is
.. math:: f_n = \frac{1}{N} \sum_{k=0}^{N-1} \tilde f_k e^{i 2 \pi k n / N}.
Implementation in real time/frequency using FFTW
---------------------------------------------------------------
The real time mesh parameters are :math:`t_{min}`, :math:`\delta t` and :math:`N_t`.
For the real frequency mesh, they are :math:`\omega_{min}`, :math:`\delta \omega` and :math:`N_\omega`.
The Fourier transform requires :math:`N_\omega=N_t` and :math:`\delta t \delta \omega= \frac{2\pi}{N_t}`.
The times are :math:`t_k=t_{min}+k\delta t` and the frequencies :math:`\omega_m=\omega_{min}+m\delta \omega`.
By approximating Eq. :ref:`TF_R` by
.. math:: \tilde G(\omega_m) = \delta t \sum_{k=0}^{N_t} G(t_k) e^{i\omega_m t_k},
we recognize an inverse DFT (Eq. :ref:`inv_DFT`). To calculate it using FFTW, we first need to prepare the input :math:`\tilde f_k`, then to do the DFT and finally to modify the output to obtain :math:`\tilde G(\omega_m)` using the two formulas:
.. math:: \tilde f_k = G(t_k) e^{i \omega_{min}t_k},
.. math:: \tilde G(\omega_m) = \delta t f_m e^{i t_{min}(\omega_m-\omega_{min})}.
Similarly, the inverse transformation is obtained by approximating Eq. :ref:`eq_inv_TF_R` by
.. math:: G(t_k)=\frac{\delta\omega}{2\pi}\sum_{m=0}^{N_\omega} \tilde G(\omega_m)e^{-i\omega_m t_k},
we recognize a DFT (Eq. :ref:`DFT`). To calculate it using FFTW, we first need to prepare the input :math:`f_m`, then to do the inverse DFT and finally to modify the output to obtain :math:`G(t_k)`:
.. math:: f_m = \tilde G(\omega_m) e^{-i t_{min}\omega_m},
.. math:: G(t_k) = \frac{1}{N_t \delta t}\tilde f_k e^{-i \omega_{min}(t_k-t_{min})}.
Implementation in imaginary time/frequency using FFTW
---------------------------------------------------------------
The imaginary time mesh parameters are :math:`\beta` and :math:`N_\tau`, plus a tag ``half_bins``, ``full_bins`` or ``without_last``.
In the ``full_bins`` case, one point of the time GF has to be removed for the fourier transform.
From these parameters, we deduce :math:`\delta\tau=\beta/N_\tau`.
For the imaginary frequency mesh, the mesh parameters are :math:`\beta`, :math:`n_{min}` and :math:`N_{\omega_n}`.
From them, we deduce :math:`\delta\omega=\frac{2\pi}{\beta}`.
The Fourier transform requires :math:`N_\omega=N_\tau`.
The times are :math:`\tau_k=\tau_{min}+k\delta\tau` and the frequencies :math:`\omega_n=\omega_{min}+n\delta \omega`.
:math:`\tau_{min}` is either 0 or :math:`\delta\tau/2` depending on the mesh kind.
:math:`\omega_{min}` is either :math:`\frac{2\pi(n_{min}+1)}{\beta}` or :math:`\frac{2\pi n_{min}}{\beta}` depending on the statistics.
We approximate the TF and its inverse by
.. math:: \tilde G(i\omega_n) = \delta\tau \sum_{k=0}^{N_\tau} G(\tau_k)e^{i\omega_n \tau_k}
.. math:: G(\tau_k) = \sum_{n=n_{min}}^{N_\tau} \frac{1}{\beta} \tilde G(i\omega_n)e^{-i\omega_n \tau_k}
We use for the TF:
.. math:: \tilde f_k = G(\tau_k) e^{i \omega_{min}\tau_k},
.. math:: \tilde G(i\omega_n) = \frac{\beta}{N_\tau} f_n e^{i \tau_{min}(\omega_n-\omega_{min})}.
and for the inverse TF:
.. math:: f_m = \frac{1}{\beta}\tilde G(i\omega_n) e^{-i t_{min}\omega_n},
.. math:: G(t_k) = \tilde f_k e^{-i \omega_{min}(\tau_k-\tau_{min})},
Special case of real functions in time for fermions
...........................................................
In this case, :math:`G(i\omega_n)=conj(G(i\omega_n))` and we only store the values of :math:`G(i\omega_n)` for :math:`\omega_n > 0`.
The Eq. :ref:`inv_DFT_I` becomes:
.. math:: G(\tau)=\sum_{n=0}^\infty \frac{2}{\beta} \tilde G(i\omega_n)\cos(\omega_n \tau)
The inverse TF formulas are in this case
.. math:: f_m = \frac{2}{\beta}\tilde G(i\omega_n) e^{-i t_{min}\omega_n},
.. math:: G(t_k) = \tilde f_k e^{-i \omega_{min}(\tau_k-\tau_{min})},
Special case of real functions in time for bosons
...........................................................
In this case, :math:`G(i\omega_n)=conj(G(i\omega_n))` and we only store the values of :math:`G(i\omega_n)` for :math:`\omega_n \ge 0`.
The Eq. :ref:`inv_DFT_I` becomes:
.. math:: G(\tau)=\frac{1}{\beta} \tilde G(0)+\sum_{n=1}^\infty \frac{2}{\beta} \tilde G(i\omega_n)\cos(\omega_n \tau)
The inverse TF formulas are in this case
.. math:: f_0 = \frac{1}{\beta}\tilde G(0),
.. math:: f_m = \frac{2}{\beta}\tilde G(i\omega_n) \cos(t_{min}\omega_n),
.. math:: G(t_k) = \tilde f_k e^{-i \omega_{min}(\tau_k-\tau_{min})},
Effect of a TF on the tail
------------------------------------------
The tail is unchanged during a TF, but its value is used to limit the errors.
The components 1 and 2 of the tail (:math:`t_1` and :math:`t_2`) are used to improve the computation of the GF in the following way:
in the large :math:`\omega` limit,
.. math:: G(\omega)\simeq \frac{t_1}{\omega}+\frac{t_2}{\omega^2}\simeq \frac{a_1}{\omega+i}+\frac{a_2}{\omega-i}
with :math:`a_1=\frac{t_1+it_2}{2}` and :math:`a_2=\frac{t_1-it_2}{2}`.
As these large w terms are badly taken into account if we naively Fourier transform the function described by its values on the mesh in w, we substract them, do the Fourier transform and add their Fourier transform to the result.
We use the following Fourier tranforms:
.. math:: \frac{2a}{\omega^2+a^2} \leftrightarrow e^{-a|x|}
.. math:: \frac{1}{\omega+i} \leftrightarrow -i e^{-x} \theta(x)
.. math:: \frac{1}{\omega-i} \leftrightarrow i e^{x} \theta(-x)
For the inverse Fourier transform, the inverse procedure is used.
In the library, :math:`a` is optimized according to the mesh properties (its size :math:`L=G.mesh().size()` and its precision :math:`\delta = G.mesh().delta()`).
The requirements are :math:`a \gg \delta\omega` and :math:`a \ll L\delta\omega`, or equivalently :math:`a \gg \delta t` and :math:`a \ll L\delta t`.
Thus, we chose :math:`a=\sqrt{L}\delta\omega`

20
pytriqs/gf/local/functions.pxd
View File

@ -8,22 +8,22 @@ cdef extern from "triqs/gfs/local/functions.hpp":
############### Fourier ######################### ############### Fourier #########################
cdef extern from "triqs/gfs/local/fourier_matsubara.hpp" : cdef extern from "triqs/gfs/local/fourier_matsubara.hpp" :
gf_imfreq lazy_fourier (gf_imtime & ) gf_imfreq fourier (gf_imtime & )
gf_imtime lazy_inverse_fourier (gf_imfreq & ) gf_imtime inverse_fourier (gf_imfreq & )
cdef extern from "triqs/gfs/local/fourier_real.hpp" : cdef extern from "triqs/gfs/local/fourier_real.hpp" :
gf_refreq lazy_fourier (gf_retime & ) gf_refreq fourier (gf_retime & )
gf_retime lazy_inverse_fourier (gf_refreq & ) gf_retime inverse_fourier (gf_refreq & )
gf_refreq fourier (gf_retime & ) except + #gf_refreq fourier (gf_retime & ) except +
gf_retime inverse_fourier (gf_refreq & ) except + #gf_retime inverse_fourier (gf_refreq & ) except +
############### Legendre ######################### ############### Legendre #########################
cdef extern from "triqs/gfs/local/legendre_matsubara.hpp" : cdef extern from "triqs/gfs/local/legendre_matsubara.hpp" :
gf_imfreq lazy_legendre_imfreq (gf_legendre &) gf_imfreq legendre_to_imfreq (gf_legendre &)
gf_imtime lazy_legendre_imtime (gf_legendre &) gf_imtime legendre_to_imtime (gf_legendre &)
gf_legendre lazy_imfreq_legendre (gf_imfreq &) gf_legendre imfreq_to_legendre (gf_imfreq &)
gf_legendre lazy_imtime_legendre (gf_imtime &) gf_legendre imtime_to_legendre (gf_imtime &)
############### Pade ######################### ############### Pade #########################

4
pytriqs/gf/local/imfreq.pyx
View File

@ -12,11 +12,11 @@ cdef class GfImFreq_cython:
def set_from_fourier(self,GfImTime_cython gt) : def set_from_fourier(self,GfImTime_cython gt) :
"""Fills self with the Fourier transform of gt""" """Fills self with the Fourier transform of gt"""
self._c << lazy_fourier( gt._c ) self._c << fourier( gt._c )
def set_from_legendre(self, GfLegendre_cython gl) : def set_from_legendre(self, GfLegendre_cython gl) :
"""Fills self with the Legendre transform of gl""" """Fills self with the Legendre transform of gl"""
self._c << lazy_legendre_imfreq(gl._c) self._c << legendre_to_imfreq(gl._c)
def density(self): def density(self):
return density(self._c).to_python() return density(self._c).to_python()

4
pytriqs/gf/local/imtime.pyx
View File

@ -12,11 +12,11 @@ cdef class GfImTime_cython:
def set_from_inverse_fourier(self,GfImFreq_cython gw) : def set_from_inverse_fourier(self,GfImFreq_cython gw) :
"""Fills self with the Inverse Fourier transform of gw""" """Fills self with the Inverse Fourier transform of gw"""
self._c << lazy_inverse_fourier( gw._c) self._c << inverse_fourier( gw._c)
def set_from_legendre(self, GfLegendre_cython gl) : def set_from_legendre(self, GfLegendre_cython gl) :
"""Fills self with the Legendre transform of gl""" """Fills self with the Legendre transform of gl"""
self._c << lazy_legendre_imtime(gl._c) self._c << legendre_to_imtime(gl._c)
def __dealloc__ (self): def __dealloc__ (self):
pass pass

4
pytriqs/gf/local/legendre.pyx
View File

@ -10,11 +10,11 @@ cdef class GfLegendre_cython:
def set_from_imtime(self, GfImTime_cython gt) : def set_from_imtime(self, GfImTime_cython gt) :
"""Fills self with the Legendre transform of gt""" """Fills self with the Legendre transform of gt"""
self._c << lazy_imtime_legendre(gt._c) self._c << imtime_to_legendre(gt._c)
def set_from_imfreq(self, GfImFreq_cython gw) : def set_from_imfreq(self, GfImFreq_cython gw) :
"""Fills self with the Legendre transform of gw""" """Fills self with the Legendre transform of gw"""
self._c << lazy_imfreq_legendre(gw._c) self._c << imfreq_to_legendre(gw._c)
def density(self): def density(self):
return density(self._c).to_python() return density(self._c).to_python()

7
pytriqs/gf/local/refreq.pyx
View File

@ -10,10 +10,11 @@ cdef class GfReFreq_cython:
def set_from_fourier(self, GfReTime_cython gt) : def set_from_fourier(self, GfReTime_cython gt) :
"""Fills self with the Fourier transform of gt""" """Fills self with the Fourier transform of gt"""
self._c << lazy_fourier( gt._c ) self._c << fourier( gt._c )
def inverse_fourier(self): # put if back with make_gf_from_fourier when approved
return make_GfReTime(inverse_fourier(self._c)) #def inverse_fourier(self):
# return make_GfReTime(inverse_fourier(self._c))
def set_from_pade(self, GfImFreq_cython gw, n_points = 100, freq_offset = 0.0) : def set_from_pade(self, GfImFreq_cython gw, n_points = 100, freq_offset = 0.0) :
pade(self._c, gw._c, n_points, freq_offset) pade(self._c, gw._c, n_points, freq_offset)

7
pytriqs/gf/local/retime.pyx
View File

@ -8,12 +8,13 @@ cdef class GfReTime_cython:
def __write_hdf5_cython__ (self, gr , char * key) : def __write_hdf5_cython__ (self, gr , char * key) :
h5_write (make_h5_group(gr), key, self._c) h5_write (make_h5_group(gr), key, self._c)
def fourier(self): # Put it back ?
return make_GfReFreq(fourier(self._c)) #def fourier(self):
# return make_GfReFreq(fourier(self._c))
def set_from_inverse_fourier(self, GfReFreq_cython gw) : def set_from_inverse_fourier(self, GfReFreq_cython gw) :
"""Fills self with the Inverse Fourier transform of gw""" """Fills self with the Inverse Fourier transform of gw"""
self._c << lazy_inverse_fourier( gw._c) self._c << inverse_fourier( gw._c)
def __dealloc__ (self): def __dealloc__ (self):
pass pass

6
test/triqs/gfs/curry_and_fourier.cpp
View File

@ -35,11 +35,11 @@ try {
auto G_w_wn_curry0 = curry<0>(G_w_wn); auto G_w_wn_curry0 = curry<0>(G_w_wn);
auto G_w_tau_curry0 = curry<0>(G_w_tau); auto G_w_tau_curry0 = curry<0>(G_w_tau);
for (auto const & w : G_w_wn_curry0.mesh()) G_w_wn_curry0[w] = lazy_fourier(G_w_tau_curry0[w]); for (auto const & w : G_w_wn_curry0.mesh()) G_w_wn_curry0[w] = fourier(G_w_tau_curry0[w]);
G_w_wn_curry0[w_] << lazy_fourier(G_w_tau_curry0[w_]); G_w_wn_curry0[w_] << fourier(G_w_tau_curry0[w_]);
curry<0>(G_w_wn) [w_] << lazy_fourier(curry<0>(G_w_tau)[w_]); curry<0>(G_w_wn) [w_] << fourier(curry<0>(G_w_tau)[w_]);
} }
// temp fix : THE TEST DOES NOT RUN !! // temp fix : THE TEST DOES NOT RUN !!
//TRIQS_CATCH_AND_ABORT; //TRIQS_CATCH_AND_ABORT;

10
test/triqs/gfs/fourier1.cpp
View File

@ -12,13 +12,13 @@ int main() {
auto G = gf<imfreq> {{beta, Fermion}, {2,2}}; auto G = gf<imfreq> {{beta, Fermion}, {2,2}};
auto Gc = G; auto Gc = G;
auto G3 = G; auto G3 = G;
auto Gt = gf<imtime> {{beta, Fermion}, {2,2}}; auto Gt = gf<imtime> {{beta, Fermion,100}, {2,2}};
auto gt = inverse_fourier(G); auto gt = make_gf_from_inverse_fourier(G);
auto gw = fourier(gt); auto gw = make_gf_from_fourier(gt);
//gw() = lazy_fourier(gt); //gw() = fourier(gt);
G() = lazy_fourier(Gt); G() = fourier(Gt);
} }

38
test/triqs/gfs/gf_notail.cpp Normal file
View File

@ -0,0 +1,38 @@
#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
#include <triqs/gfs.hpp>
using namespace triqs::gfs;
using namespace triqs::arrays;
#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> " << (X) << std::endl << std::endl;
#include <triqs/gfs/local/functions.hpp>
int main() {
try {
double beta = 1, a = 1;
// construct
auto gw_n = gf<imfreq, matrix_valued, no_tail>{{beta, Fermion}, {2, 2}};
auto gt_n = gf<imtime, matrix_valued, no_tail>{{beta, Fermion, 100}, {2, 2}};
// make a view from a g with a tail
auto gw = gf<imfreq, matrix_valued>{{beta, Fermion}, {2, 2}};
auto gt = gf<imtime, matrix_valued>{gf_mesh<imtime>{beta, Fermion, 100}, {2, 2}};
auto vw = make_gf_view_without_tail(gw);
auto vt = make_gf_view_without_tail(gt);
triqs::clef::placeholder<0> tau_;
// should not and does not compile : improve error message ???
// gt(tau_) << exp(-a * tau_) / (1 + exp(-beta * a));
vt(tau_) << exp(-a * tau_) / (1 + exp(-beta * a));
// test hdf5
H5::H5File file("ess_g_notail.h5", H5F_ACC_TRUNC);
h5_write(file, "g", vt);
// rebuilding a new gf...
auto g3 = make_gf_from_g_and_tail(vw, gw.singularity());
// need to test all this....
}
TRIQS_CATCH_AND_ABORT;
}

2
test/triqs/gfs/gfv2.cpp
View File

@ -20,7 +20,7 @@ int main() {
auto G = gf<imfreq>{ {beta, Fermion}, {2,2} }; auto G = gf<imfreq>{ {beta, Fermion}, {2,2} };
auto Gc = gf<imfreq>{ {beta, Fermion}, {2,2} }; auto Gc = gf<imfreq>{ {beta, Fermion}, {2,2} };
auto G3 = gf<imfreq>{ {beta, Fermion}, {2,2} }; auto G3 = gf<imfreq>{ {beta, Fermion}, {2,2} };
auto Gt = gf<imtime>{ {beta, Fermion}, {2,2} }; auto Gt = gf<imtime>{ {beta, Fermion,100}, {2,2} };
auto Gv = G(); auto Gv = G();
TEST( G( 0) ) ; TEST( G( 0) ) ;

12
test/triqs/gfs/test_fourier_matsubara.cpp
View File

@ -5,6 +5,14 @@ using namespace triqs::arrays;
#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> "<< (X) <<std::endl<<std::endl; #define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> "<< (X) <<std::endl<<std::endl;
#include <triqs/gfs/local/fourier_matsubara.hpp> #include <triqs/gfs/local/fourier_matsubara.hpp>
namespace triqs { namespace gfs {
// defined in the cpp file
void fourier_impl (gf_view<imfreq,scalar_valued> gw , gf_const_view<imtime,scalar_valued> gt, scalar_valued);
void fourier_impl (gf_view<imfreq,matrix_valued> gw , gf_const_view<imtime,matrix_valued> gt, matrix_valued);
void inverse_fourier_impl (gf_view<imtime,scalar_valued> gt, gf_const_view<imfreq,scalar_valued> gw, scalar_valued);
void inverse_fourier_impl (gf_view<imtime,matrix_valued> gt, gf_const_view<imfreq,matrix_valued> gw, matrix_valued);
}}
int main() { int main() {
double precision=10e-9; double precision=10e-9;
@ -45,9 +53,9 @@ int main() {
} }
h5_write(file,"Gt1b",Gt1); // must be 0 h5_write(file,"Gt1b",Gt1); // must be 0
///to verify that lazy_fourier computes ///to verify that fourier computes
auto Gw2 = gf<imfreq> {{beta, Fermion}, {1,1}}; auto Gw2 = gf<imfreq> {{beta, Fermion}, {1,1}};
Gw2() = lazy_fourier(Gt1); Gw2() = fourier(Gt1);
} }

8
test/triqs/gfs/test_fourier_real_time.cpp
View File

@ -33,11 +33,11 @@ int main() {
Gw1.singularity()(2)=triqs::arrays::matrix<double>{{2.0*a}}; Gw1.singularity()(2)=triqs::arrays::matrix<double>{{2.0*a}};
h5_write(file,"Gw1",Gw1); // the original lorentzian h5_write(file,"Gw1",Gw1); // the original lorentzian
auto Gt1 = inverse_fourier(Gw1); auto Gt1 = make_gf_from_inverse_fourier(Gw1);
h5_write(file,"Gt1",Gt1); // the lorentzian TF : lorentzian_inverse h5_write(file,"Gt1",Gt1); // the lorentzian TF : lorentzian_inverse
// verification that TF(TF^-1)=Id // verification that TF(TF^-1)=Id
auto Gw1b = fourier(Gt1); auto Gw1b = make_gf_from_fourier(Gt1);
for(auto const & w:Gw1b.mesh()){ for(auto const & w:Gw1b.mesh()){
Gw1b[w]-=Gw1[w]; 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"; 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";
@ -65,7 +65,7 @@ int main() {
Gt2.singularity()(1)=triqs::arrays::matrix<double>{{1.0}}; Gt2.singularity()(1)=triqs::arrays::matrix<double>{{1.0}};
h5_write(file,"Gt2",Gt2); h5_write(file,"Gt2",Gt2);
auto Gw2 = fourier(Gt2); auto Gw2 = make_gf_from_fourier(Gt2);
h5_write(file,"Gw2",Gw2); h5_write(file,"Gw2",Gw2);
for(auto const & w:Gw2.mesh()){ for(auto const & w:Gw2.mesh()){
@ -86,7 +86,7 @@ int main() {
//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_j*std::sin(8*t+0.3*acos(-1.)) ;
h5_write(file,"Gt3",Gt3); h5_write(file,"Gt3",Gt3);
auto Gw3 = fourier(Gt3); auto Gw3 = make_gf_from_fourier(Gt3);
h5_write(file,"Gw3",Gw3); h5_write(file,"Gw3",Gw3);
} }

4
test/triqs/gfs/test_gf_triqs.cpp
View File

@ -33,7 +33,7 @@ void test_0(){
/* ---------- Fourier transform ---------------------*/ /* ---------- Fourier transform ---------------------*/
auto Gt = gf<imtime> {{beta, Fermion, Ntau, full_bins}, {1,1}}; auto Gt = gf<imtime> {{beta, Fermion, Ntau, full_bins}, {1,1}};
Gt() = lazy_inverse_fourier(G); Gt() = inverse_fourier(G);
TEST(Gt(0.0)); TEST(Gt(0.0));
@ -52,7 +52,7 @@ void test_1(){
auto Gw = gf<imfreq> {{beta, Fermion, 100}, {1,1}}; auto Gw = gf<imfreq> {{beta, Fermion, 100}, {1,1}};
Gw.singularity()(1) = 1; Gw.singularity()(1) = 1;
Gt() = lazy_inverse_fourier(Gw); Gt() = inverse_fourier(Gw);
} }
int main() { int main() {

3
triqs/gfs.hpp
View File

@ -30,5 +30,8 @@
#include <triqs/gfs/product.hpp> #include <triqs/gfs/product.hpp>
#include <triqs/gfs/curry.hpp> #include <triqs/gfs/curry.hpp>
#include <triqs/gfs/local/fourier_matsubara.hpp>
#include <triqs/gfs/local/fourier_real.hpp>
#include <triqs/gfs/local/legendre_matsubara.hpp>
#endif #endif

3
triqs/gfs/domains/matsubara.hpp
View File

@ -30,7 +30,8 @@ namespace triqs { namespace gfs {
typedef typename mpl::if_c<IsComplex, std::complex<double>, double>::type point_t; typedef typename mpl::if_c<IsComplex, std::complex<double>, double>::type point_t;
double beta; double beta;
statistic_enum statistic; statistic_enum statistic;
matsubara_domain (double Beta=1, statistic_enum s = Fermion): beta(Beta), statistic(s){ matsubara_domain() = default;
matsubara_domain (double Beta, statistic_enum s): beta(Beta), statistic(s){
if(beta<0)TRIQS_RUNTIME_ERROR<<"Matsubara domain construction : beta <0 : beta ="<< beta <<"\n"; if(beta<0)TRIQS_RUNTIME_ERROR<<"Matsubara domain construction : beta <0 : beta ="<< beta <<"\n";
} }
matsubara_domain(matsubara_domain const &) = default; matsubara_domain(matsubara_domain const &) = default;

16
triqs/gfs/gf.hpp
View File

@ -94,6 +94,9 @@ namespace triqs { namespace gfs {
template <typename G> auto get_gf_data_shape(G const &g) DECL_AND_RETURN(g.get_data_shape()); template <typename G> auto get_gf_data_shape(G const &g) DECL_AND_RETURN(g.get_data_shape());
template <typename Variable, typename Target, typename Opt, bool IsView, bool IsConst>
auto get_target_shape(gf_impl<Variable, Target, Opt, IsView, IsConst> const &g) DECL_AND_RETURN(g.data().shape().front_pop());
// ---------------------- implementation -------------------------------- // ---------------------- implementation --------------------------------
// overload get_shape for a vector to simplify code below in gf block case. // overload get_shape for a vector to simplify code below in gf block case.
@ -334,8 +337,10 @@ namespace triqs { namespace gfs {
template<typename RHS> friend void triqs_clef_auto_assign_subscript (gf_impl & g, RHS rhs) { triqs_clef_auto_assign(g,rhs);} template<typename RHS> friend void triqs_clef_auto_assign_subscript (gf_impl & g, RHS rhs) { triqs_clef_auto_assign(g,rhs);}
private: private:
template<typename RHS> void triqs_clef_auto_assign_impl (RHS const & rhs, std::integral_constant<bool,false>) { 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)); //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>) { template <typename RHS> void triqs_clef_auto_assign_impl(RHS const &rhs, std::integral_constant<bool, true>) {
@ -582,7 +587,11 @@ namespace triqs { namespace gfs {
template <typename Var, typename Opt> struct factories<Var, scalar_valued, Opt> { template <typename Var, typename Opt> struct factories<Var, scalar_valued, Opt> {
typedef gf<Var, scalar_valued, Opt> gf_t; typedef gf<Var, scalar_valued, Opt> gf_t;
struct target_shape_t {}; struct target_shape_t {
target_shape_t front_pop() const { // this make the get_target_shape function works in this case...
return {};
}
};
typedef typename gf_t::mesh_t mesh_t; typedef typename gf_t::mesh_t mesh_t;
static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t shape) { static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t shape) {
@ -620,6 +629,7 @@ namespace triqs { namespace gfs {
}; };
} // gfs_implementation } // gfs_implementation
}} }}
// same as for arrays : views can not be swapped by the std::swap. Delete it // same as for arrays : views can not be swapped by the std::swap. Delete it

6
triqs/gfs/imfreq.hpp
View File

@ -23,6 +23,7 @@
#include "./tools.hpp" #include "./tools.hpp"
#include "./gf.hpp" #include "./gf.hpp"
#include "./local/tail.hpp" #include "./local/tail.hpp"
#include "./local/no_tail.hpp"
#include "./domains/matsubara.hpp" #include "./domains/matsubara.hpp"
#include "./meshes/linear.hpp" #include "./meshes/linear.hpp"
#include "./evaluators.hpp" #include "./evaluators.hpp"
@ -34,6 +35,7 @@ namespace triqs { namespace gfs {
typedef linear_mesh<matsubara_domain<true>> B; typedef linear_mesh<matsubara_domain<true>> B;
static double m1(double beta) { return std::acos(-1)/beta;} static double m1(double beta) { return std::acos(-1)/beta;}
gf_mesh() = default; gf_mesh() = default;
gf_mesh(B const &x) : B(x) {} // enables also construction from another Opt
gf_mesh (typename B::domain_t const & d, int Nmax = 1025) : gf_mesh (typename B::domain_t const & d, int Nmax = 1025) :
B(d, d.statistic==Fermion?m1(d.beta):0, d.statistic==Fermion?(2*Nmax+1)*m1(d.beta): 2*Nmax*m1(d.beta), Nmax, without_last){} B(d, d.statistic==Fermion?m1(d.beta):0, d.statistic==Fermion?(2*Nmax+1)*m1(d.beta): 2*Nmax*m1(d.beta), Nmax, without_last){}
gf_mesh (double beta, statistic_enum S, int Nmax = 1025) : gf_mesh({beta,S}, Nmax){} gf_mesh (double beta, statistic_enum S, int Nmax = 1025) : gf_mesh({beta,S}, Nmax){}
@ -42,8 +44,8 @@ namespace triqs { namespace gfs {
namespace gfs_implementation { namespace gfs_implementation {
//singularity //singularity
template<typename Opt> struct singularity<imfreq,matrix_valued,Opt> { typedef local::tail type;}; template<> struct singularity<imfreq,matrix_valued,void> { typedef local::tail type;};
template<typename Opt> struct singularity<imfreq,scalar_valued,Opt> { typedef local::tail type;}; template<> struct singularity<imfreq,scalar_valued,void> { typedef local::tail type;};
//h5 name //h5 name
template<typename Opt> struct h5_name<imfreq,matrix_valued,Opt> { static std::string invoke(){ return "ImFreq";}}; template<typename Opt> struct h5_name<imfreq,matrix_valued,Opt> { static std::string invoke(){ return "ImFreq";}};

9
triqs/gfs/imtime.hpp
View File

@ -23,6 +23,7 @@
#include "./tools.hpp" #include "./tools.hpp"
#include "./gf.hpp" #include "./gf.hpp"
#include "./local/tail.hpp" #include "./local/tail.hpp"
#include "./local/no_tail.hpp"
#include "./domains/matsubara.hpp" #include "./domains/matsubara.hpp"
#include "./meshes/linear.hpp" #include "./meshes/linear.hpp"
#include "./evaluators.hpp" #include "./evaluators.hpp"
@ -35,15 +36,16 @@ namespace triqs { namespace gfs {
template <typename Opt> struct gf_mesh<imtime, Opt> : linear_mesh<matsubara_domain<false>> { template <typename Opt> struct gf_mesh<imtime, Opt> : linear_mesh<matsubara_domain<false>> {
typedef linear_mesh<matsubara_domain<false>> B; typedef linear_mesh<matsubara_domain<false>> B;
gf_mesh() = default; gf_mesh() = default;
gf_mesh(B const &x) : B(x) {} // enables also construction from another Opt
gf_mesh(typename B::domain_t d, int n_time_slices, mesh_kind mk = half_bins) : B(d, 0, d.beta, n_time_slices, mk) {} gf_mesh(typename B::domain_t d, int n_time_slices, mesh_kind mk = half_bins) : B(d, 0, d.beta, n_time_slices, mk) {}
gf_mesh(double beta, statistic_enum S, int n_time_slices, mesh_kind mk = half_bins) : gf_mesh({beta, S}, n_time_slices, mk) {} gf_mesh(double beta, statistic_enum S, int n_time_slices, mesh_kind mk = half_bins) : gf_mesh({beta, S}, n_time_slices, mk) {}
}; };
namespace gfs_implementation { namespace gfs_implementation {
// singularity // singularity. If no_tail is given, then it is the default (nothing)
template<typename Opt> struct singularity<imtime,matrix_valued,Opt> { typedef local::tail type;}; template<> struct singularity<imtime,matrix_valued,void> { typedef local::tail type;};
template<typename Opt> struct singularity<imtime,scalar_valued,Opt> { typedef local::tail type;}; template<> struct singularity<imtime,scalar_valued,void> { typedef local::tail type;};
// h5 name // h5 name
template<typename Opt> struct h5_name<imtime,matrix_valued,Opt> { static std::string invoke(){ return "ImTime";}}; template<typename Opt> struct h5_name<imtime,matrix_valued,Opt> { static std::string invoke(){ return "ImTime";}};
@ -92,6 +94,7 @@ namespace triqs { namespace gfs {
template<typename Opt, typename Target> struct evaluator<imtime,Target,Opt> : evaluator_one_var<imtime>{}; template<typename Opt, typename Target> struct evaluator<imtime,Target,Opt> : evaluator_one_var<imtime>{};
} // gfs_implementation. } // gfs_implementation.
}} }}
#endif #endif

251
triqs/gfs/local/fourier_matsubara.cpp
View File

@ -18,25 +18,15 @@
* TRIQS. If not, see <http://www.gnu.org/licenses/>. * TRIQS. If not, see <http://www.gnu.org/licenses/>.
* *
******************************************************************************/ ******************************************************************************/
#include "fourier_base.hpp"
#include "fourier_matsubara.hpp" #include "fourier_matsubara.hpp"
#include <fftw3.h> #include <fftw3.h>
namespace triqs { namespace gfs { namespace triqs {
namespace gfs {
namespace impl_local_matsubara { template <typename GfElementType> GfElementType convert_green(dcomplex const& x) { return x; }
template <> double convert_green<double>(dcomplex const& x) { return real(x); }
inline 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)) );
}
inline 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)) );
}
}
template<typename GfElementType> GfElementType convert_green ( dcomplex const& x) { return x;}
template<> double convert_green<double> ( dcomplex const& x) { return real(x);}
//-------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------
@ -44,142 +34,187 @@ namespace triqs { namespace gfs {
tqa::vector<dcomplex> g_in, g_out; tqa::vector<dcomplex> g_in, g_out;
void direct (gf_view<imfreq,scalar_valued> gw, gf_const_view<imtime,scalar_valued> gt) { dcomplex oneFermion(dcomplex a, double b, double tau, double beta) {
using namespace impl_local_matsubara; return -a * (b >= 0 ? exp(-b * tau) / (1 + exp(-beta * b)) : exp(b * (beta - tau)) / (1 + exp(beta * b)));
auto ta = gt(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=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());
double fact= beta/ gt.mesh().size();
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(gt.mesh().size());
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())
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) ) );
}
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.singularity() = gt.singularity();// set tail
} }
void inverse(gf_view<imtime,scalar_valued> gt, gf_const_view<imfreq,scalar_valued> gw){ dcomplex oneBoson(dcomplex a, double b, double tau, double beta) {
using namespace impl_local_matsubara; 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) =delete;
//-------------------------------------
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());
double fact = beta / gt.mesh().size();
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(gt.mesh().size());
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())
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)));
}
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; static bool Green_Function_Are_Complex_in_time = false;
// If the Green function are NOT complex, then one use the symmetry property // If the Green function are NOT complex, then one use the symmetry property
// fold the sum and get a factor 2 // fold the sum and get a factor 2
auto ta = gw(freq_infty()); auto ta = gw(freq_infty());
//TO BE MODIFIED AFTER SCALAR IMPLEMENTATION TODO // 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); 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; double b1, b2, b3;
dcomplex a1, a2, a3; dcomplex a1, a2, a3;
double beta=gw.domain().beta; 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 size_t L = gt.mesh().size() - (gt.mesh().kind() == full_bins ? 1 : 0); // L can be different from gt.mesh().size() (depending
dcomplex iomega = dcomplex(0.0,1.0) * std::acos(-1) / beta; // on the mesh kind) and is given to the FFT algorithm
dcomplex iomega2 = -iomega * 2 * gt.mesh().delta() * (gt.mesh().kind() == half_bins ? 0.5 : 0.0) ; dcomplex iomega = dcomplex(0.0, 1.0) * std::acos(-1) / beta;
double fact = (Green_Function_Are_Complex_in_time ? 1 : 2)/beta; dcomplex iomega2 = -iomega * 2 * gt.mesh().delta() * (gt.mesh().kind() == half_bins ? 0.5 : 0.0);
g_in.resize( gw.mesh().size()); double fact = (Green_Function_Are_Complex_in_time ? 1 : 2) / beta;
g_in.resize(gw.mesh().size());
g_out.resize(gt.mesh().size()); g_out.resize(gt.mesh().size());
if (gw.domain().statistic == Fermion){ if (gw.domain().statistic == Fermion) {
b1 = 0; b2 =1; b3 =-1; b1 = 0;
a1 = d-B; a2 = (A+B)/2; a3 = (B-A)/2; b2 = 1;
} b3 = -1;
else { a1 = d - B;
b1 = -0.5; b2 =-1; b3 =1; a2 = (A + B) / 2;
a1=4*(d-B)/3; a2=B-(d+A)/2; a3=d/6+A/2+B/3; 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; g_in() = 0;
for (auto & w: gw.mesh()) { 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 // 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; 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); details::fourier_base(g_in, g_out, L, false);
// CORRECT FOR COMPLEX G(tau) !!! // CORRECT FOR COMPLEX G(tau) !!!
typedef double gt_result_type; typedef double gt_result_type;
//typedef typename gf<imtime>::mesh_type::gf_result_type gt_result_type; // typedef typename gf<imtime>::mesh_type::gf_result_type gt_result_type;
if (gw.domain().statistic == Fermion){ if (gw.domain().statistic == Fermion) {
for (auto & t : gt.mesh()){ 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] =
+ oneFermion(a1,b1,t,beta) + oneFermion(a2,b2,t,beta)+ oneFermion(a3,b3,t,beta) ); 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()) + oneBoson(a1, b1, t, beta) +
oneBoson(a2, b2, t, beta) + oneBoson(a3, b3, t, beta));
} }
else { if (gt.mesh().kind() == full_bins) gt.on_mesh(L) = -gt.on_mesh(0) - convert_green<gt_result_type>(ta(1)(0, 0));
for (auto & t : gt.mesh())
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)
gt.on_mesh(L) = -gt.on_mesh(0)-convert_green<gt_result_type>(ta(1)(0,0));
// set tail // set tail
gt.singularity() = gw.singularity(); gt.singularity() = gw.singularity();
} }
}; }; // class worker
void fourier_impl (gf_view<imfreq,scalar_valued> gw , gf_const_view<imtime,scalar_valued> gt, scalar_valued){ //--------------------------------------------
void fourier_impl(gf_view<imfreq, scalar_valued> gw, gf_const_view<imtime, scalar_valued> gt, scalar_valued) {
impl_worker w; impl_worker w;
w.direct(gw,gt); w.direct(gw, gt);
} }
void fourier_impl (gf_view<imfreq,matrix_valued> gw , gf_const_view<imtime,matrix_valued> gt, matrix_valued){ void fourier_impl(gf_view<imfreq, matrix_valued> gw, gf_const_view<imtime, matrix_valued> gt, matrix_valued) {
impl_worker w; impl_worker w;
for (size_t n1=0; n1<gt.data().shape()[1];n1++) for (size_t n1 = 0; n1 < gt.data().shape()[1]; n1++)
for (size_t n2=0; n2<gt.data().shape()[2];n2++){ for (size_t n2 = 0; n2 < gt.data().shape()[2]; n2++) {
auto gw_sl=slice_target_to_scalar(gw,n1,n2); auto gw_sl = slice_target_to_scalar(gw, n1, n2);
auto gt_sl=slice_target_to_scalar(gt,n1,n2); auto gt_sl = slice_target_to_scalar(gt, n1, n2);
w.direct(gw_sl, gt_sl); w.direct(gw_sl, gt_sl);
} }
} }
//--------------------------------------------------------------------------- //---------------------------------------------------------------------------
void inverse_fourier_impl (gf_view<imtime,scalar_valued> gt , gf_const_view<imfreq,scalar_valued> gw, scalar_valued){ void inverse_fourier_impl(gf_view<imtime, scalar_valued> gt, gf_const_view<imfreq, scalar_valued> gw, scalar_valued) {
impl_worker w; impl_worker w;
w.inverse(gt,gw); w.inverse(gt, gw);
} }
void inverse_fourier_impl (gf_view<imtime,matrix_valued> gt , gf_const_view<imfreq,matrix_valued> gw, matrix_valued){ void inverse_fourier_impl(gf_view<imtime, matrix_valued> gt, gf_const_view<imfreq, matrix_valued> gw, matrix_valued) {
impl_worker w; impl_worker w;
for (size_t n1=0; n1<gw.data().shape()[1];n1++) for (size_t n1 = 0; n1 < gw.data().shape()[1]; n1++)
for (size_t n2=0; n2<gw.data().shape()[2];n2++){ for (size_t n2 = 0; n2 < gw.data().shape()[2]; n2++) {
auto gt_sl=slice_target_to_scalar(gt, n1, n2); auto gt_sl = slice_target_to_scalar(gt, n1, n2);
auto gw_sl=slice_target_to_scalar(gw, n1, n2); auto gw_sl = slice_target_to_scalar(gw, n1, n2);
w.inverse( gt_sl, gw_sl); 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, scalar_valued() ); } void triqs_gf_view_assign_delegation(gf_view<imfreq, scalar_valued> 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, matrix_valued() ); } gf_keeper<tags::fourier, imtime, scalar_valued> const& L) {
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, scalar_valued() ); } fourier_impl(g, L.g, scalar_valued());
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, matrix_valued() ); } }
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, matrix_valued());
}
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, scalar_valued());
}
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, matrix_valued());
}
}
}

89
triqs/gfs/local/fourier_matsubara.hpp
View File

@ -25,60 +25,55 @@
#include <triqs/gfs/imfreq.hpp> #include <triqs/gfs/imfreq.hpp>
#include <triqs/gfs/imtime.hpp> #include <triqs/gfs/imtime.hpp>
namespace triqs { namespace gfs { namespace triqs {
namespace gfs {
// First the implementation of the fourier transform template <typename Target, typename Opt, bool V, bool C>
void fourier_impl (gf_view<imfreq,scalar_valued> gw , gf_const_view<imtime,scalar_valued> gt, scalar_valued); gf_keeper<tags::fourier, imtime, Target> fourier(gf_impl<imtime, Target, Opt, V, C> const& g) {
void fourier_impl (gf_view<imfreq,matrix_valued> gw , gf_const_view<imtime,matrix_valued> gt, matrix_valued); return {g};
void inverse_fourier_impl (gf_view<imtime,scalar_valued> gt, gf_const_view<imfreq,scalar_valued> gw, scalar_valued);
void inverse_fourier_impl (gf_view<imtime,matrix_valued> gt, gf_const_view<imfreq,matrix_valued> gw, matrix_valued);
inline gf_view<imfreq,matrix_valued> fourier (gf_const_view<imtime,matrix_valued> gt) {
int L = (gt.mesh().kind() == full_bins ? gt.mesh().size()-1 : gt.mesh().size() );
auto gw = gf<imfreq,matrix_valued>{ {gt.domain(),L}, gt.data().shape().front_pop() };
auto V = gw();
fourier_impl(V, gt, matrix_valued());
return gw;
} }
inline gf_view<imfreq,scalar_valued> fourier (gf_const_view<imtime,scalar_valued> gt) { template <typename Target, typename Opt, bool V, bool C>
int L = (gt.mesh().kind() == full_bins ? gt.mesh().size()-1 : gt.mesh().size() ); gf_keeper<tags::fourier, imfreq, Target> inverse_fourier(gf_impl<imfreq, Target, Opt, V, C> const& g) {
auto gw = gf<imfreq,scalar_valued>{ {gt.domain(),L} }; return {g};
auto V = gw(); }
fourier_impl(V, gt, scalar_valued());
void triqs_gf_view_assign_delegation(gf_view<imfreq, scalar_valued> g, gf_keeper<tags::fourier, imtime, scalar_valued> const& L);
void triqs_gf_view_assign_delegation(gf_view<imfreq, matrix_valued> g, gf_keeper<tags::fourier, imtime, matrix_valued> const& L);
void triqs_gf_view_assign_delegation(gf_view<imtime, scalar_valued> g, gf_keeper<tags::fourier, imfreq, scalar_valued> const& L);
void triqs_gf_view_assign_delegation(gf_view<imtime, matrix_valued> g, gf_keeper<tags::fourier, imfreq, matrix_valued> const& L);
template <typename Opt> gf_mesh<imfreq, Opt> make_mesh_fourier_compatible(gf_mesh<imtime, Opt> const& m) {
int L = m.size() - (m.kind() == full_bins ? 1 : 0);
return {m.domain(), L};
}
template <typename Opt>
gf_mesh<imtime, Opt> make_mesh_fourier_compatible(gf_mesh<imfreq, Opt> const& m, mesh_kind mk = full_bins) {
int L = m.size() + (mk == full_bins ? 1 : 0);
return {m.domain(), L};
}
template <typename Target, typename Opt, bool V, bool C>
gf_view<imfreq, Target> make_gf_from_fourier(gf_impl<imtime, Target, Opt, V, C> const& gt) {
auto gw = gf<imfreq, Target>{make_mesh_fourier_compatible(gt.mesh()), get_target_shape(gt)};
gw() = fourier(gt);
return gw; return gw;
} }
inline gf_view<imtime, matrix_valued> inverse_fourier (gf_const_view<imfreq, matrix_valued> gw, mesh_kind mk = half_bins) { template <typename Target, typename Opt, bool V, bool C>
double pi = std::acos(-1); gf_view<imtime, Target> make_gf_from_inverse_fourier(gf_impl<imfreq, Target, Opt, V, C> const& gw, mesh_kind mk = full_bins) {
int L = (mk == full_bins ? gw.mesh().size()+1 : gw.mesh().size() ); auto gt = gf<imtime, Target>{make_mesh_fourier_compatible(gw.mesh(), mk), get_target_shape(gw)};
auto gt = gf<imtime,matrix_valued>{ {gw.domain(),L}, gw.data().shape().front_pop()}; gt() = inverse_fourier(gw);
auto V = gt();
inverse_fourier_impl(V, gw, matrix_valued());
return gt;
}
inline gf_view<imtime,scalar_valued> inverse_fourier (gf_const_view<imfreq,scalar_valued> gw, mesh_kind mk = half_bins) {
double pi = std::acos(-1);
int L = (mk == full_bins ? gw.mesh().size()+1 : gw.mesh().size() );
auto gt = gf<imtime,scalar_valued>{ {gw.domain(),L} };
auto V = gt();
inverse_fourier_impl(V, gw,scalar_valued());
return gt; return gt;
} }
}
}
inline gf_keeper<tags::fourier,imtime,scalar_valued> lazy_fourier (gf_const_view<imtime,scalar_valued> g) { return {g};} namespace triqs {
inline gf_keeper<tags::fourier,imfreq,scalar_valued> lazy_inverse_fourier (gf_const_view<imfreq,scalar_valued> g) { return {g};} namespace clef {
inline gf_keeper<tags::fourier,imtime,matrix_valued> lazy_fourier (gf_const_view<imtime,matrix_valued> g) { return {g};} TRIQS_CLEF_MAKE_FNT_LAZY(fourier);
inline gf_keeper<tags::fourier,imfreq,matrix_valued> lazy_inverse_fourier (gf_const_view<imfreq,matrix_valued> g) { return {g};} TRIQS_CLEF_MAKE_FNT_LAZY(inverse_fourier);
}
void triqs_gf_view_assign_delegation( gf_view<imfreq,scalar_valued> g, gf_keeper<tags::fourier,imtime,scalar_valued> const & L); }
void triqs_gf_view_assign_delegation( gf_view<imfreq,matrix_valued> g, gf_keeper<tags::fourier,imtime,matrix_valued> const & L);
void triqs_gf_view_assign_delegation( gf_view<imtime,scalar_valued> g, gf_keeper<tags::fourier,imfreq,scalar_valued> const & L);
void triqs_gf_view_assign_delegation( gf_view<imtime,matrix_valued> g, gf_keeper<tags::fourier,imfreq,matrix_valued> const & L);
}}
namespace triqs { namespace clef {
TRIQS_CLEF_MAKE_FNT_LAZY (lazy_fourier);
TRIQS_CLEF_MAKE_FNT_LAZY (lazy_inverse_fourier);
}}
#endif #endif

89
triqs/gfs/local/fourier_real.hpp
View File

@ -2,7 +2,7 @@
* *
* TRIQS: a Toolbox for Research in Interacting Quantum Systems * TRIQS: a Toolbox for Research in Interacting Quantum Systems
* *
* Copyright (C) 2011 by M. Ferrero, O. Parcollet * Copyright (C) 2011-2014 by M. Ferrero, O. Parcollet
* *
* TRIQS is free software: you can redistribute it and/or modify it under the * 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 * terms of the GNU General Public License as published by the Free Software
@ -27,63 +27,50 @@
namespace triqs { namespace gfs { namespace triqs { namespace gfs {
// First the implementation of the fourier transform template <typename Target, typename Opt, bool V, bool C>
void fourier_impl (gf_view<refreq,scalar_valued> gw , gf_const_view<retime,scalar_valued> gt, scalar_valued); gf_keeper<tags::fourier, retime, Target> fourier(gf_impl<retime, Target, Opt, V, C> const& g) {
void fourier_impl (gf_view<refreq,matrix_valued> gw , gf_const_view<retime,matrix_valued> gt, matrix_valued); return {g};
void inverse_fourier_impl (gf_view<retime,scalar_valued> gt, gf_const_view<refreq,scalar_valued> gw, scalar_valued);
void inverse_fourier_impl (gf_view<retime,matrix_valued> gt, gf_const_view<refreq,matrix_valued> gw, matrix_valued);
inline gf_view<refreq,matrix_valued> fourier (gf_const_view<retime, matrix_valued> gt) {
double pi = std::acos(-1);
int L = gt.mesh().size();
double wmin = -pi * (L-1) / (L*gt.mesh().delta());
double wmax = pi * (L-1) / (L*gt.mesh().delta());
auto gw = gf<refreq,matrix_valued>{ {wmin, wmax, L}, gt.data().shape().front_pop()};
auto V = gw();
fourier_impl(V, gt, matrix_valued());
return gw;
} }
inline gf_view<refreq,scalar_valued> fourier (gf_const_view<retime, scalar_valued> gt) { template <typename Target, typename Opt, bool V, bool C>
gf_keeper<tags::fourier, refreq, Target> inverse_fourier(gf_impl<refreq, Target, Opt, V, C> const& g) {
return {g};
}
void triqs_gf_view_assign_delegation(gf_view<refreq, scalar_valued> g, gf_keeper<tags::fourier, retime, scalar_valued> const& L);
void triqs_gf_view_assign_delegation(gf_view<refreq, matrix_valued> g, gf_keeper<tags::fourier, retime, matrix_valued> const& L);
void triqs_gf_view_assign_delegation(gf_view<retime, scalar_valued> g, gf_keeper<tags::fourier, refreq, scalar_valued> const& L);
void triqs_gf_view_assign_delegation(gf_view<retime, matrix_valued> g, gf_keeper<tags::fourier, refreq, matrix_valued> const& L);
template <typename Opt> gf_mesh<refreq, Opt> make_mesh_fourier_compatible(gf_mesh<retime, Opt> const& m) {
int L = m.size();
double pi = std::acos(-1); double pi = std::acos(-1);
int L = gt.mesh().size(); double wmin = -pi * (L - 1) / (L * m.delta());
double wmin = -pi * (L-1) / (L*gt.mesh().delta()); double wmax = pi * (L - 1) / (L * m.delta());
double wmax = pi * (L-1) / (L*gt.mesh().delta()); return {wmin, wmax, L};
auto gw = gf<refreq,scalar_valued>{ {wmin, wmax, L} }; }
auto V = gw();
fourier_impl(V, gt, scalar_valued()); template <typename Opt>
gf_mesh<retime, Opt> make_mesh_fourier_compatible(gf_mesh<refreq, Opt> const& m, mesh_kind mk = full_bins) {
double pi = std::acos(-1);
int L = m.size();
double tmin = -pi * (L-1) / (L*m.delta());
double tmax = pi * (L-1) / (L*m.delta());
return {tmin, tmax, L};
}
template <typename Target, typename Opt, bool V, bool C>
gf_view<refreq, Target> make_gf_from_fourier(gf_impl<retime, Target, Opt, V, C> const& gt) {
auto gw = gf<refreq, Target>{make_mesh_fourier_compatible(gt.mesh()), get_target_shape(gt)};
gw() = fourier(gt);
return gw; return gw;
} }
inline gf_view<retime,matrix_valued> inverse_fourier (gf_const_view<refreq,matrix_valued> gw) { template <typename Target, typename Opt, bool V, bool C>
double pi = std::acos(-1); gf_view<retime, Target> make_gf_from_inverse_fourier(gf_impl<refreq, Target, Opt, V, C> const& gw) {
int L = gw.mesh().size(); auto gt = gf<retime, Target>{make_mesh_fourier_compatible(gw.mesh()), get_target_shape(gw)};
double tmin = -pi * (L-1) / (L*gw.mesh().delta()); gt() = inverse_fourier(gw);
double tmax = pi * (L-1) / (L*gw.mesh().delta());
auto gt = gf<retime,matrix_valued>{{ tmin, tmax, L} , gw.data().shape().front_pop()};
auto V = gt();
inverse_fourier_impl(V, gw, matrix_valued());
return gt; return gt;
} }
inline gf_view<retime,scalar_valued> inverse_fourier (gf_const_view<refreq,scalar_valued> gw) {
double pi = std::acos(-1);
int L = gw.mesh().size();
double tmin = -pi * (L-1) / (L*gw.mesh().delta());
double tmax = pi * (L-1) / (L*gw.mesh().delta());
auto gt = gf<retime,scalar_valued>{ {tmin, tmax, L} };
auto V = gt();
inverse_fourier_impl(V, gw, scalar_valued());
return gt;
}
inline gf_keeper<tags::fourier, retime, scalar_valued> lazy_fourier(gf_const_view<retime, scalar_valued> g) { return {g}; }
inline gf_keeper<tags::fourier, refreq, scalar_valued> lazy_inverse_fourier(gf_const_view<refreq, scalar_valued> g) { return {g}; }
inline gf_keeper<tags::fourier, retime, matrix_valued> lazy_fourier(gf_const_view<retime, matrix_valued> g) { return {g}; }
inline gf_keeper<tags::fourier, refreq, matrix_valued> lazy_inverse_fourier(gf_const_view<refreq, matrix_valued> g) { return {g}; }
void triqs_gf_view_assign_delegation( gf_view<refreq,scalar_valued> g, gf_keeper<tags::fourier,retime,scalar_valued> const & L);
void triqs_gf_view_assign_delegation( gf_view<refreq,matrix_valued> g, gf_keeper<tags::fourier,retime,matrix_valued> const & L);
void triqs_gf_view_assign_delegation( gf_view<retime,scalar_valued> g, gf_keeper<tags::fourier,refreq,scalar_valued> const & L);
void triqs_gf_view_assign_delegation( gf_view<retime,matrix_valued> g, gf_keeper<tags::fourier,refreq,matrix_valued> const & L);
}} }}
#endif #endif

116
triqs/gfs/local/legendre_matsubara.cpp
View File

@ -1,9 +1,8 @@
/******************************************************************************* /*******************************************************************************
* *
* TRIQS: a Toolbox for Research in Interacting Quantum Systems * TRIQS: a Toolbox for Research in Interacting Quantum Systems
* *
* Copyright (C) 2011 by L. Boehnke, M. Ferrero, O. Parcollet * Copyright (C) 2011-2014 by L. Boehnke, M. Ferrero, O. Parcollet
* *
* TRIQS is free software: you can redistribute it and/or modify it under the * 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 * terms of the GNU General Public License as published by the Free Software
@ -27,25 +26,63 @@
using namespace triqs::utility; using namespace triqs::utility;
namespace triqs { namespace gfs { namespace triqs {
namespace gfs {
void legendre_matsubara_direct(gf_view<imfreq> & gw, gf_const_view<legendre> gl) { // ----------------------------
void legendre_matsubara_direct(gf_view<imfreq> gw, gf_const_view<legendre> gl) {
gw() = 0.0; gw() = 0.0;
triqs::arrays::range R; triqs::arrays::range R;
// Use the transformation matrix // Use the transformation matrix
for (auto om: gw.mesh()) { for (auto om : gw.mesh()) {
for (auto l: gl.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];
} }
} }
gw.singularity() = get_tail(gl, gw.singularity().size(), gw.singularity().order_min()); gw.singularity() = get_tail(gl, gw.singularity().size(), gw.singularity().order_min());
}
} // ----------------------------
void legendre_matsubara_inverse (gf_view<legendre> & gl, gf_const_view<imfreq> gw) { void legendre_matsubara_direct(gf_view<imtime> gt, gf_const_view<legendre> gl) {
gt() = 0.0;
legendre_generator L;
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.singularity() = get_tail(gl, gt.singularity().size(), gt.singularity().order_min());
}
// ----------------------------
void legendre_matsubara_inverse(gf_view<legendre> gl, gf_const_view<imtime> gt) {
gl() = 0.0;
legendre_generator L;
// Do the integral over imaginary time
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.data() *= gt.mesh().delta();
}
// ----------------------------
void legendre_matsubara_inverse(gf_view<legendre> gl, gf_const_view<imfreq> gw) {
gl() = 0.0; gl() = 0.0;
@ -53,64 +90,25 @@ void legendre_matsubara_inverse (gf_view<legendre> & gl, gf_const_view<imfreq>
// I set Nt time bins. This is ugly, one day we must code the direct // I set Nt time bins. This is ugly, one day we must code the direct
// transformation without going through imaginary time // transformation without going through imaginary time
int Nt = 50000; int Nt = 50000;
auto gt = gf<imtime>{ {gw.domain(),Nt, half_bins}, gw.data().shape().front_pop() }; auto gt = gf<imtime>{{gw.domain(), Nt, half_bins}, gw.data().shape().front_pop()};
// We first transform to imaginary time because it's been coded with the knowledge of the tails // We first transform to imaginary time because it's been coded with the knowledge of the tails
gt() = lazy_inverse_fourier(gw); gt() = inverse_fourier(gw);
legendre_matsubara_inverse(gl, gt()); legendre_matsubara_inverse(gl, gt());
}
void legendre_matsubara_direct (gf_view<imtime> & gt, gf_const_view<legendre> gl) {
gt() = 0.0;
legendre_generator L;
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.singularity() = get_tail(gl, gt.singularity().size(), gt.singularity().order_min());
} void triqs_gf_view_assign_delegation(gf_view<imfreq> gw, gf_keeper<tags::legendre, legendre> const& L) {
void legendre_matsubara_inverse (gf_view<legendre> & gl, gf_const_view<imtime> gt) {
gl() = 0.0;
legendre_generator L;
// Do the integral over imaginary time
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.data() *= gt.mesh().delta();
}
gf_keeper<tags::legendre,legendre> lazy_legendre_imfreq (gf_const_view<legendre> gl) { return {gl}; }
gf_keeper<tags::legendre,legendre> lazy_legendre_imtime (gf_const_view<legendre> gl) { return {gl}; }
gf_keeper<tags::legendre,imfreq> lazy_imfreq_legendre (gf_const_view<imfreq> gw) { return {gw}; }
gf_keeper<tags::legendre,imtime> lazy_imtime_legendre (gf_const_view<imtime> gt) { return {gt}; }
void triqs_gf_view_assign_delegation( gf_view<imfreq> &gw, gf_keeper<tags::legendre,legendre> const & L) {
legendre_matsubara_direct(gw, L.g); legendre_matsubara_direct(gw, L.g);
} }
void triqs_gf_view_assign_delegation( gf_view<imtime> &gt, gf_keeper<tags::legendre,legendre> const & L) { void triqs_gf_view_assign_delegation(gf_view<imtime> gt, gf_keeper<tags::legendre, legendre> const& L) {
legendre_matsubara_direct(gt, L.g); legendre_matsubara_direct(gt, L.g);
} }
void triqs_gf_view_assign_delegation( gf_view<legendre> &gl, gf_keeper<tags::legendre,imfreq> const & L) { void triqs_gf_view_assign_delegation(gf_view<legendre> gl, gf_keeper<tags::legendre, imfreq> const& L) {
legendre_matsubara_inverse(gl, L.g); legendre_matsubara_inverse(gl, L.g);
} }
void triqs_gf_view_assign_delegation( gf_view<legendre> &gl, gf_keeper<tags::legendre,imtime> const & L) { void triqs_gf_view_assign_delegation(gf_view<legendre> gl, gf_keeper<tags::legendre, imtime> const& L) {
legendre_matsubara_inverse(gl, L.g); legendre_matsubara_inverse(gl, L.g);
}
}
} }
}}

45
triqs/gfs/local/legendre_matsubara.hpp
View File

@ -1,9 +1,8 @@
/******************************************************************************* /*******************************************************************************
* *
* TRIQS: a Toolbox for Research in Interacting Quantum Systems * TRIQS: a Toolbox for Research in Interacting Quantum Systems
* *
* Copyright (C) 2011 by L. Boehnke, M. Ferrero, O. Parcollet * Copyright (C) 2011-2014 by L. Boehnke, M. Ferrero, O. Parcollet
* *
* TRIQS is free software: you can redistribute it and/or modify it under the * 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 * terms of the GNU General Public License as published by the Free Software
@ -26,26 +25,30 @@
#include <triqs/gfs/imtime.hpp> #include <triqs/gfs/imtime.hpp>
#include <triqs/gfs/legendre.hpp> #include <triqs/gfs/legendre.hpp>
namespace triqs { namespace gfs { namespace triqs {
namespace gfs {
void legendre_matsubara_direct (gf_view<imfreq> &gw, gf_const_view<legendre> gl); namespace tags {
void legendre_matsubara_inverse (gf_view<legendre> &gl, gf_const_view<imfreq> gw); struct legendre {};
}
void legendre_matsubara_direct (gf_view<imtime> &gt, gf_const_view<legendre> gl); inline gf_keeper<tags::legendre, legendre> legendre_to_imfreq(gf_const_view<legendre> gl) {
void legendre_matsubara_inverse (gf_view<legendre> &gl, gf_const_view<imtime> gt); return {gl};
}
inline gf_keeper<tags::legendre, legendre> legendre_to_imtime(gf_const_view<legendre> gl) {
return {gl};
}
inline gf_keeper<tags::legendre, imfreq> imfreq_to_legendre(gf_const_view<imfreq> gw) {
return {gw};
}
inline gf_keeper<tags::legendre, imtime> imtime_to_legendre(gf_const_view<imtime> gt) {
return {gt};
}
namespace tags { struct legendre{}; } void triqs_gf_view_assign_delegation(gf_view<imfreq> gw, gf_keeper<tags::legendre, legendre> const &L);
void triqs_gf_view_assign_delegation(gf_view<imtime> gt, gf_keeper<tags::legendre, legendre> const &L);
gf_keeper<tags::legendre,legendre> lazy_legendre_imfreq (gf_const_view<legendre> gl); void triqs_gf_view_assign_delegation(gf_view<legendre> gl, gf_keeper<tags::legendre, imfreq> const &L);
gf_keeper<tags::legendre,legendre> lazy_legendre_imtime (gf_const_view<legendre> gl); void triqs_gf_view_assign_delegation(gf_view<legendre> gl, gf_keeper<tags::legendre, imtime> const &L);
gf_keeper<tags::legendre,imfreq> lazy_imfreq_legendre (gf_const_view<imfreq> gw); }
gf_keeper<tags::legendre,imtime> lazy_imtime_legendre (gf_const_view<imtime> gt); }
void triqs_gf_view_assign_delegation( gf_view<imfreq> &gw, gf_keeper<tags::legendre,legendre> const & L);
void triqs_gf_view_assign_delegation( gf_view<imtime> &gt, gf_keeper<tags::legendre,legendre> const & L);
void triqs_gf_view_assign_delegation( gf_view<legendre> &gl, gf_keeper<tags::legendre,imfreq> const & L);
void triqs_gf_view_assign_delegation( gf_view<legendre> &gl, gf_keeper<tags::legendre,imtime> const & L);
}}
#endif #endif

48
triqs/gfs/local/no_tail.hpp Normal file
View File

@ -0,0 +1,48 @@
/*******************************************************************************
*
* TRIQS: a Toolbox for Research in Interacting Quantum Systems
*
* Copyright (C) 2013-2014 by 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/>.
*
******************************************************************************/
#pragma once
#include "../tools.hpp"
#include "../gf.hpp"
#include "./tail.hpp"
namespace triqs {
namespace gfs {
struct no_tail {};
template <typename Variable, typename Target, bool V, bool C>
gf_view<Variable, Target, no_tail, C> make_gf_view_without_tail(gf_impl<Variable, Target, void, V, C> const &g) {
return {g.mesh(), g.data(), {}, g.symmetry()};
}
template <typename Variable, typename Target, bool V, bool C>
gf_view<Variable, Target> make_gf_from_g_and_tail(gf_impl<Variable, Target, no_tail, V, C> const &g, local::tail t) {
auto g2 = gf<Variable, Target, no_tail>{g}; // copy the function without tail
return {std::move(g2.mesh()), std::move(g2.data()), std::move(t), g2.symmetry()};
}
template <typename Variable, typename Target, bool V, bool C>
gf_view<Variable, Target, void, C> make_gf_view_from_g_and_tail(gf_impl<Variable, Target, no_tail, V, C> const &g,
local::tail_view t) {
return {g.mesh(), g.data(), t, g.symmetry()};
}
}
}

2
triqs/gfs/meshes/linear.hpp
View File

@ -43,7 +43,7 @@ namespace gfs {
linear_mesh() : _dom(), L(0), a_pt(0), b_pt(0), xmin(0), xmax(0), del(0), meshk(half_bins) {} linear_mesh() : _dom(), L(0), a_pt(0), b_pt(0), xmin(0), xmax(0), del(0), meshk(half_bins) {}
linear_mesh(domain_t dom, double a, double b, size_t n_pts, mesh_kind mk) explicit linear_mesh(domain_t dom, double a, double b, size_t n_pts, mesh_kind mk)
: _dom(std::move(dom)), L(n_pts), a_pt(a), b_pt(b), meshk(mk) { : _dom(std::move(dom)), L(n_pts), a_pt(a), b_pt(b), meshk(mk) {
switch (mk) { switch (mk) {
case half_bins: case half_bins:

2
triqs/gfs/tools.hpp
View File

@ -50,6 +50,8 @@ namespace triqs { namespace gfs {
struct freq_infty{}; // the point at infinity struct freq_infty{}; // the point at infinity
//------------------------------------------------------
inline std::vector<std::string> split(const std::string &s, char delim){ inline std::vector<std::string> split(const std::string &s, char delim){
std::vector<std::string> elems; std::vector<std::string> elems;
std::stringstream ss(s); std::stringstream ss(s);