NESSi  v1.0.2
The NonEquilibrium Systems Simulation Library
cntr_matsubara_impl.hpp
Go to the documentation of this file.
1 #ifndef CNTR_MATSUBARA_IMPL_H
2 #define CNTR_MATSUBARA_IMPL_H
3 
5 #include "fourier.hpp"
6 #include "cntr_elements.hpp"
7 
8 namespace cntr {
9 
10 // ----------------------------------------------------------------------
11 
12 // Helper functions for imaginary time and matsubara frequencies
14 template <typename F, typename I> F get_tau(const I tau_idx, const F beta, const I ntau) {
15  return tau_idx * (beta / ntau);
16 }
17 
19 template <typename F, typename I> F get_omega(const I m, const F beta, const I sig) {
20 
21  // matsub_one = 1 if sig = -1, 0 if sig = 1
22  // (shift is integer division by 2)
23  const I matsub_one = ((1 - sig) >> 1);
24  // const I matsub_one = (sig == -1 ? 1.0 : 0.0);
25 
26  // return fourier::pi/beta * (2*m + matsub_one);
27  return PI / beta * (2 * m + matsub_one);
28 }
29 
30 // ----------------------------------------------------------------------
32 
70 template <typename T, int SIZE1>
71 void set_first_order_tail(std::complex<T> *xmat, std::complex<T> *coeff, T beta, int sg,
72  int ntau, int sig, int size1) {
73 
74  // First order tail correction tau-polynomial for Bosons & Fermions
75  // Fermions : g_0(\tau) = -1/2
76  // Bosons : g_0(\tau) = -1/2 + \tau/\beta
77 
78  std::complex<T> *z1 = new std::complex<T>[sg];
79 
80  for (int r = 0; r <= ntau; r++) {
81  // -1/2 factor
82  element_set<T, SIZE1>(size1, xmat + r * sg, coeff);
83  element_smul<T, SIZE1>(size1, xmat + r * sg, -0.5);
84 
85  if (sig == 1) {
86  // tau/beta factor
87  double tau = r * (beta / ntau);
88  element_set<T, SIZE1>(size1, z1, coeff);
89  element_smul<T, SIZE1>(size1, z1, tau / beta);
90  element_incr<T, SIZE1>(size1, xmat + r * sg, z1);
91  }
92  }
93 
94  delete[] z1;
95 }
96 
97 // ----------------------------------------------------------------------
98 
99 /*-------------------------------------------------------
100 compute (WITHOUT FFT)
101 I=sum_{r=0}^m ftau[r] exp(i omega_n tau_r)
102 for n=0...m
103 tau_r=r*dtau,tau=beta/m
104 omega_n= (2n+1)pi/beta, (sgn=-1)
105 omega_n= 2n*pi/beta, (sgn=+1)
106 
107 I=sum_{r=0}^m ftau[r] exp(i pi (2n+1)r/m)
108 -------------------------------------------------------*/
109 
110 
133 template <typename T, class GG, int SIZE1>
134 void matsubara_dft(std::complex<T> *mdft, GG &G, int sig) {
135  typedef std::complex<T> cplx;
136  int ntau, r, m, sg, size1 = G.size1();
137  double arg, one;
138  cplx *z, *z1, expfac;
139  sg = G.element_size();
140  ntau = G.ntau();
141  z = new cplx[sg];
142  z1 = new cplx[sg];
143  one = (sig == -1 ? 1.0 : 0.0);
144  for (m = 0; m <= ntau; m++) {
145  element_set_zero<T, SIZE1>(size1, z);
146  for (r = 0; r <= ntau; r++) {
147  arg = ((2 * m + one) * r * PI) / ntau;
148  expfac = cplx(cos(arg), sin(arg));
149  element_set<T, SIZE1>(size1, z1, G.matptr(r));
150  element_smul<T, SIZE1>(size1, z1, expfac);
151  element_incr<T, SIZE1>(size1, z, z1);
152  }
153  element_set<T, SIZE1>(size1, mdft + m * sg, z);
154  }
155  delete[] z;
156  delete[] z1;
157 }
158 
190 template <typename T, class GG, int SIZE1>
191 void matsubara_ft(std::complex<T> *result, int m, GG &G, std::complex<T> *mdft, int sig,
192  T beta, int order) {
193  typedef std::complex<T> cplx;
194  double corfac, arg, dtau, one;
195  int ntau, m1, r, l, sg, size1 = G.size1();
196  cplx *z1, *z2, *dft, bcorr[4];
197 
198  one = (sig == -1 ? 1.0 : 0.0);
199  ntau = G.ntau();
200  sg = G.element_size();
201  dtau = beta / ntau;
202  z1 = new cplx[sg];
203  z2 = new cplx[sg];
204  dft = new cplx[sg];
205 
206  arg = (2 * m + one) * PI / ntau; /*arg=omn*dtau, omn=(2m+1)Pi/beta*/
207  if (order == 0) {
208  corfac = 1.0;
209  bcorr[0] = 0.0;
210  } else if (order == 1) {
211  fourier::get_dftcorr_linear(arg, &corfac, bcorr);
212  } else if (order == 3) {
213  fourier::get_dftcorr_cubic(arg, &corfac, bcorr);
214  } else {
215  std::cerr << "matsubara_ft: unknown order" << std::endl;
216  abort();
217  }
218  m1 = m; /*shift m1 to interval [0,ntau-1]*/
219  while (m1 < 0)
220  m1 += ntau;
221  while (m1 >= ntau)
222  m1 -= ntau;
223  for (l = 0; l < sg; l++)
224  dft[l] = mdft[m1 * sg + l] * corfac;
225  element_set_zero<T, SIZE1>(size1, z1);
226  for (r = 0; r <= order; r++) {
227  element_set<T, SIZE1>(size1, z2, G.matptr(r));
228  element_smul<T, SIZE1>(size1, z2, bcorr[r]);
229  element_incr<T, SIZE1>(size1, z1, z2);
230  element_set<T, SIZE1>(size1, z2, G.matptr(ntau - r));
231  // element_smul<T,SIZE1>(size1,z2,-cplx(bcorr[r].real(),-bcorr[r].imag()));
232  element_smul<T, SIZE1>(size1, z2,
233  (1. * sig) * cplx(bcorr[r].real(), -bcorr[r].imag()));
234  element_incr<T, SIZE1>(size1, z1, z2);
235  }
236  for (l = 0; l < sg; l++)
237  result[l] = dtau * (z1[l] + dft[l]);
238 
239  delete[] z1;
240  delete[] z2;
241  delete[] dft;
242 }
243 
244 // ----------------------------------------------------------------------
245 
246 } // namespace cntr
247 
248 #endif // CNTR_MATSUBARA_IMPL_H
std::complex< double > cplx
Definition: fourier.cpp:11
Integrator< T > & I(int k)
void get_dftcorr_linear(double th, double *corfac, cplx *endcor)
Returns correction factor and boundary correction term needed for linearly corrected DFT...
Definition: fourier.cpp:117
void matsubara_dft(std::complex< T > *mdft, GG &G, int sig)
Computes the Fourier series coefficients of a Matsubara function by plain DFT.
void matsubara_ft(std::complex< T > *result, int m, GG &G, std::complex< T > *mdft, int sig, T beta, int order)
Computes the Fourier series coefficients of a Matsubara function by cubically corrected DFT...
void get_dftcorr_cubic(double th, double *corfac, cplx *endcor)
Returns correction factor and boundary correction terms needed for cubically corrected DFT...
Definition: fourier.cpp:43