NESSi  v1.0.2
The NonEquilibrium Systems Simulation Library
cntr_equilibrium_impl.hpp
Go to the documentation of this file.
1 #ifndef CNTR_EQUILIBRIUM_IMPL_H
2 #define CNTR_EQUILIBRIUM_IMPL_H
3 
6 #include "fourier.hpp"
7 #include "integration.hpp"
8 #include "cntr_elements.hpp"
9 //#include "cntr_exception.hpp"
13 #include "cntr_function_decl.hpp"
14 #include "cntr_utilities_decl.hpp"
15 
16 namespace cntr {
17 
18 /*####################################################################################
19 #
20 # The computation of the equilibrium propagator from the density of states
21 #
22 ######################################################################################*/
23 
24 #define EXPMAX 100
25 
43 template <typename T>
44 T fermi(T beta, T omega) {
45  T arg = omega * beta;
46  if (fabs(arg) > EXPMAX) {
47  return (arg > 0.0 ? 0.0 : 1.0);
48  } else {
49  return 1.0 / (1.0 + exp(arg));
50  }
51 }
70 template<typename T>
71 dvector fermi(T beta,dvector &omega){
72  int size=omega.size();
73  dvector tmp(size);
74  for(int i=0;i<size;i++){
75  tmp(i)=fermi(beta,omega(i));
76  }
77  return tmp;
78 }
80 template <typename T>
81 T fermi_exp(T beta, T tau, T omega) {
82  if (omega < 0) {
83  return exp(omega * tau) *
84  fermi(beta, omega); // exp(w*t)/(1+exp(b*w)) always OK for w<0
85  } else {
86  return exp((tau - beta) * omega) *
87  fermi(beta, -omega); // exp((t-b)*w)/(1+exp(-w*b))
88  }
89 }
91 template<typename T>
92 dvector fermi_exp(T beta,T tau,dvector &omega){
93  int size=omega.size();
94  dvector tmp(size);
95 
96  for(int i=0;i<size;i++){
97  tmp(i)=fermi_exp(beta,tau,omega(i));
98  }
99  return tmp;
100 }
102 template<typename T>
103 cdmatrix diag_prop(T time,dvector &omega){
104  int size=omega.size();
105  cdmatrix tmp(size,size);
106  for(int i=0;i<size;i++){
107  tmp(i,i)=std::complex<T>(cos(omega(i)*time),sin(omega(i)*time));
108  }
109  return tmp;
110 }
111 
130 template <typename T>
131 T bose(T beta, T omega) {
132  T arg = omega * beta;
133  if (arg < 0)
134  return (-1.0 - bose<T>(beta, -omega));
135  if (fabs(arg) > EXPMAX) {
136  return 0.0;
137  } else if (arg < 1e-10) {
138  return 1.0 / arg;
139  } else {
140  return 1.0 / (exp(arg) - 1.0);
141  }
142 }
143 
162 template<typename T>
163 dvector bose(T beta,dvector &omega){
164  int size=omega.size();
165  dvector tmp(size);
166  for(int i=0;i<size;i++){
167  tmp(i)=bose(beta,omega(i));
168  }
169  return tmp;
170 }
171 
172 
173 // exp(tau*w)b(w) ... actually assuming positive tau
175 template <typename T>
176 T bose_exp(T beta, T tau, T omega) {
177  if (omega < 0)
178  return exp(tau * omega) * bose<T>(beta, omega);
179  else
180  return -exp((tau - beta) * omega) * bose<T>(beta, -omega);
181 }
183 template<typename T>
184 dvector bose_exp(T beta,T tau,dvector &omega){
185  int size=omega.size();
186  dvector tmp(size);
187 
188  for(int i=0;i<size;i++){
189  tmp(i)=bose_exp(beta,tau,omega(i));
190  }
191  return tmp;
192 }
193 
194 
195 #undef EXPMAX
196 
197 #define EXPMAX 100
198 template <typename T>
200 double green_cntr_full_equilibrium_fermi(T beta, T omega) {
201  T arg = omega * beta;
202  if (fabs(arg) > EXPMAX) {
203  return (arg > 0.0 ? 0.0 : 1.0);
204  } else {
205  return 1.0 / (1.0 + exp(arg));
206  }
207 }
209 template<typename T>
210 double distribution_eq(T beta,T omega,int sign)
211 {
212  if(sign==-1){
213  return green_cntr_full_equilibrium_fermi(beta,omega);
214  }else if(sign==1){
215  return bose(beta,omega);
216  }
217 }
219 template<typename T>
220 double distribution_exp_eq(T beta,T tau,T omega,int sign)
221 {
222  if(sign==-1){
223  return fermi_exp(beta,tau,omega);
224  }else if(sign==1){
225  return bose_exp(beta,tau,omega);
226  }
227 }
228 
229 
230 #undef EXPMAX
231 
233 enum flavor { ret, adv, mat, tv, vt, les, gtr };
234 
235 template < class dos > struct dos_wrapper{
236  flavor x_;
237  double beta_;
238  double tau_;
239  int sign_;
240  double mu_;
241  dos A_;
242  // init with a dos object:
243  dos_wrapper(const dos& f,int sign,double mu=0.0){
244  beta_=0;
245  tau_=0;
246  A_=f;
247  mu_=mu;
248  sign_=sign;
249  }
250  double operator()(double omega){
251  double a=A_(omega),fm;
252  if(x_==ret){
253  fm=1;
254  }else if(x_==les){
255  fm=sign_*distribution_eq(beta_,omega-mu_,sign_);
256  }else if (x_==tv){
257  fm=sign_*distribution_exp_eq(beta_,tau_,omega-mu_,sign_);
258  }else if(x_==mat){
259  if(omega-mu_>0){
260  fm=sign_*exp(-tau_*(omega-mu_))*distribution_eq(beta_,-(omega-mu_),sign_);
261  }else{
262  fm=(-1.0)*exp((beta_-tau_)*(omega-mu_))*distribution_eq(beta_,(omega-mu_),sign_);
263  }
264  }else{
265  fm=0;
266  }
267  return a*fm;
268  }
269 };
270 
271 
272 /*####################################################################################
273 #
274 # Here is the definition of some dos functions:
275 #
276 # The minimal functionality of the dos_function class is an
277 # operator()(double) which returns the values of the dos at some
278 # previously set parameters (such as the bandwidth etc.) and
279 # members lo_, hi_
280 #
281 ######################################################################################*/
282 
294 class bethedos {
295  public:
297  double hi_;
299  double lo_;
301  double V_;
303  V_ = 1;
304  lo_ = -2;
305  hi_ = 2;
306  }
307  double operator()(double x) {
308  double arg = 4.0 * V_ * V_ - x * x;
309  double num = V_ * V_ * 3.14159265358979323846 * 2;
310  return (arg < 0 ? 0.0 : sqrt(arg) / num);
311  }
312 };
314 
325 class ohmic{
326  public:
328  double hi_;
330  double lo_;
332  double omegac_;
333  ohmic() {
334  omegac_ = 1;
335  lo_ = 0.01;
336  hi_ = omegac_*20;
337  }
338  ohmic(double omegac) {
339  omegac_ = omegac;
340  lo_=0.01;
341  hi_=omegac*20; // value is 10^-6
342  }
343  double operator()(double x) {
344  return x*x*std::exp(-x/omegac_)/(2.0*omegac_*omegac_*omegac_);
345  }
346 };
347 
348 
361 class smooth_box {
362  public:
364  double hi_;
366  double lo_;
368  double A_;
370  double B_;
372  double nu_;
374  double N_;
376  A_ = -1.0;
377  B_ = 1.0;
378  nu_ = 10.0;
379  N_ = -1.0;
380  lo_ = A_ - 20.0 / nu_;
381  hi_ = B_ + 20.0 / nu_;
382  }
383  smooth_box(double a, double b, double nu) {
384  A_ = a;
385  B_ = b;
386  nu_ = nu;
387  N_ = -1.0;
388  if (nu_ > 100 || nu_ <= 0.1 || A_ >= B_) {
389  std::cerr << "nu=" << nu_ << " in smooth_box " << std::endl;
390  exit(0);
391  }
392  lo_ = A_ - 20.0 / nu_;
393  hi_ = B_ + 20.0 / nu_;
394  }
395  void set_norm(void) {
396  typedef std::complex<double> cplx;
397  cplx res, err;
398  N_ = 1.0;
399  fourier::adft_func adft;
400  adft.sample(0.0, lo_, hi_, *this, 20, 100);
401  adft.dft(0.0, res, err);
402  N_ = res.real();
403  }
404  double operator()(double x) {
405  if (N_ == -1)
406  set_norm();
407  double arg1 = fermi<double>(nu_, x - B_);
408  double arg2 = fermi<double>(nu_, A_ - x);
409  return arg1 * arg2 / N_;
410  }
411 };
412 
413 
415 template <typename T,class dos_function>
416 void green_equilibrium_ret(herm_matrix<T> &G,dos_function &dos,double h,int limit,int nn,double mu)
417 {
418  typedef std::complex<double> cplx;
419  int nt=G.nt(),i,l,size1=G.size1();
420  int sign=G.sig();
421  double t;
422  cplx res,err,cplx_i=cplx(0,1);
423  fourier::adft_func adft;
424  dos_wrapper<dos_function> dos1(dos,sign,mu);
425  dos1.mu_=mu;
426  dos1.x_=ret;
427  adft.sample(0.0,dos.lo_,dos.hi_,dos1,nn,limit);
428  for(l=0;l<=nt;l++){ // l = t-t'
429  t=h*l;
430  adft.dft(-t,res,err);
431  res *= std::complex<double>(0,-1.0);
432  for(i=l;i<=nt;i++) element_set<T,LARGESIZE>(size1,G.retptr(i,i-l),(std::complex<T>)(res));
433  }
434 }
435 
437 template <typename T,class dos_function>
438 void green_equilibrium_mat(herm_matrix<T> &G,dos_function &dos,double beta,int limit,int nn,double mu)
439 {
440  typedef std::complex<double> cplx;
441  int ntau=G.ntau(),m,size1=G.size1();
442  int sign=G.sig();
443  double dtau;
444  cplx res,err;
445  fourier::adft_func adft;
446  dos_wrapper<dos_function> dos1(dos,sign,mu);
447  dos1.beta_=beta;
448  dtau=beta/ntau;
449  dos1.mu_=mu;
450  dos1.x_=mat;
451  for(m=0;m<=ntau;m++){
452  dos1.tau_=m*dtau;
453  adft.sample(0.0,dos.lo_,dos.hi_,dos1,nn,limit);
454  adft.dft(0.0,res,err);
455  element_set<T,LARGESIZE>(size1,G.matptr(m),(std::complex<T>)(res));
456  }
457 }
459 template <typename T,class dos_function>
460 void green_equilibrium_tv(herm_matrix<T> &G,dos_function &dos,double beta,double h,int limit,int nn,double mu)
461 {
462  typedef std::complex<double> cplx;
463  int ntau=G.ntau(),nt=G.nt(),n,m,size1=G.size1();
464  double dtau;
465  int sign=G.sig();
466  cplx res,err,cplx_i=cplx(0,1);
467  fourier::adft_func adft;
468  dos_wrapper<dos_function> dos1(dos,sign,mu);
469  dos1.beta_=beta;
470  dtau=beta/ntau;
471  dos1.x_=tv;
472  dos1.mu_=mu;
473  for(m=0;m<=ntau;m++){
474  dos1.tau_=m*dtau;
475  adft.sample(0.0,dos.lo_,dos.hi_,dos1,nn,limit);
476  for(n=0;n<=nt;n++){
477  adft.dft(-n*h,res,err);
478  res *= std::complex<double>(0,-1.0);
479  element_set<T,LARGESIZE>(size1,G.tvptr(n,m),(std::complex<T>)(res));
480  }
481  }
482 }
484 template <typename T,class dos_function>
485 void green_equilibrium_les(herm_matrix<T> &G,dos_function &dos,double beta,double h,int limit,int nn,double mu)
486 {
487  typedef std::complex<double> cplx;
488  int nt=G.nt(),i,l,size1=G.size1();
489  double t;
490  int sign=G.sig();
491  cplx res,err,cplx_i=cplx(0,1);
492  fourier::adft_func adft;
493  dos_wrapper<dos_function> dos1(dos,sign,mu);
494  dos1.mu_=mu;
495  dos1.x_=les;
496  dos1.beta_=beta;
497  adft.sample(0.0,dos.lo_-0.0,dos.hi_-0.0,dos1,nn,limit);
498  for(l=0;l<=nt;l++){ // l = t'-t
499  t=h*l;
500  adft.dft(t,res,err);
501  res *= std::complex<double>(0,-1.0);
502  for(i=l;i<=nt;i++) element_set<T,LARGESIZE>(size1,G.lesptr(i-l,i),(std::complex<T>)res);
503  }
504 }
505 
507 template <typename T,class dos_function>
508 void green_equilibrium(herm_matrix<T> &G,dos_function &dos,double beta,double h,int limit,int nn,double mu)
509 {
510  green_equilibrium_mat(G,dos,beta,limit,nn,mu);
511  green_equilibrium_ret(G,dos,h,limit,nn,mu);
512  green_equilibrium_tv(G,dos,beta,h,limit,nn,mu);
513  green_equilibrium_les(G,dos,beta,h,limit,nn,mu);
514 }
515 
543 template <typename T,class dos_function>
544 void green_equilibrium(herm_matrix<T> &G,dos_function &dos,double beta,double h,double mu,int limit,int nn)
545 {
546  green_equilibrium_mat(G,dos,beta,limit,nn,mu);
547  green_equilibrium_ret(G,dos,h,limit,nn,mu);
548  green_equilibrium_tv(G,dos,beta,h,limit,nn,mu);
549  green_equilibrium_les(G,dos,beta,h,limit,nn,mu);
550 }
551 
552 
575 template <typename T>
576 void green_equilibrium_mat_bethe(herm_matrix<T> &G, double beta, int limit,
577  int nn,double mu) {
578  bethedos dos;
579  green_equilibrium_mat(G, dos, beta, limit, nn,mu);
580 }
581 
582 
608 template <typename T>
609 void green_equilibrium_bethe(herm_matrix<T> &G, double beta, double h,
610  int limit, int nn,double mu) {
611  bethedos dos;
612  green_equilibrium(G, dos, beta, h, limit, nn,mu);
613 }
614 
615 
616 
617 
618 
647 template<typename T>
648 void interpolate_CF2(int tstp,cntr::function<T> &H,cdmatrix &Hinter,int kt,bool fixHam=false){
649  int size=H.size1_;
650  cdmatrix tmp1(size,size),tmp2(size,size);
651  int ktextrap=(tstp<=kt ? tstp : kt);
652  int n1=(tstp<=kt ? kt : tstp); //Once you are at kt the next point needs to be included
653  // Extrapolate to get H(t + dt)
654  if(!fixHam && tstp>kt){
655  extrapolate_timestep(tstp-1,H,integration::I<double>(ktextrap));
656  }
657  // Interpolate to get H(t + dt/2)
658  H.get_value(tstp-1,tmp1);
659  H.get_value(tstp,tmp2);
660 
661  Hinter=interpolation(n1,(tstp-0.5),H,integration::I<double>(kt)); //ktextrap+1 since we added one more point
662 }
663 
664 
697 template<typename T>
698 void interpolate_CF4(int tstp,cntr::function<T> &H,cdmatrix &Hinte1,cdmatrix &Hinte2,int kt,bool fixHam=false){
699  int size=H.size1_;
700  cdmatrix tmp(size,size);
701  int ktextrap=(tstp<=kt ? tstp : kt);
702  int n1=(tstp<=kt ? kt : tstp);
703  // Extrapolate to get H(t)
704  if(!fixHam && tstp>kt){
705  extrapolate_timestep(tstp-1,H,integration::I<double>(ktextrap));
706  }
707  Hinte1=interpolation(n1,(tstp-0.5-sqrt(3)/6.0),H,integration::I<double>(kt)); //ktextrap+1 since we added one more point
708  Hinte2=interpolation(n1,(tstp-0.5+sqrt(3)/6.0),H,integration::I<double>(kt)); //ktextrap+1 since we added one more point
709 
710 }
711 
744 template<typename T>
745 void propagator_exp(int tstp,cntr::function<T> &U,cntr::function<T> &H,double dt,int order,int kt,bool fixHam=false){
746  int nt=U.nt_;
747  int size=U.size1_;
748  cdmatrix prop(size,size);
749 
750  assert(tstp<=nt);
751  assert(order==2 || order==4);
752  assert(size==H.size1_);
753 
754  if(tstp==-1 || tstp==0){
755  prop.setIdentity();
756  U.set_value(tstp,prop);
757  }else{
758  if(order==2){
759  cdmatrix arg(size,size);
760  // Get H(t+dt/2)-> Extrapolate and interpolate
761  interpolate_CF2(tstp,H,arg,kt,fixHam);
762  arg=arg*std::complex<double>(0.0,-1.0)*dt;
763  U.get_value(tstp-1,prop);
764  prop=arg.exp()*prop;
765  U.set_value(tstp,prop);
766  }else if(order==4){
767  cdmatrix H1(size,size),H2(size,size);
768  cdmatrix arg1(size,size),arg2(size,size);
769  // Get H(t+dt*c1) and H(t+dt*c2) -> Extrapolate and interpolate
770  interpolate_CF4(tstp,H,H1,H2,kt,fixHam);
771  double a1=(3.0-2.0*sqrt(3.0))/12.0;
772  double a2=(3.0+2.0*sqrt(3.0))/12.0;
773  arg1=std::complex<double>(0.0,-1.0)*dt*(a1*H1+a2*H2);
774  arg2=std::complex<double>(0.0,-1.0)*dt*(a2*H1+a1*H2);
775  U.get_value(tstp-1,prop);
776  prop=arg1.exp()*arg2.exp()*prop;
777  U.set_value(tstp,prop);
778  }
779  }
780 }
781 
783 template<typename T,int SIZE>
784 void green_from_H_dispatch(herm_matrix<T> &G,T mu,cntr::function<T> &eps,T beta,T h,int kt,int order,bool fixHam=false){
785  std::complex<T> iu = std::complex<T>(0.0, 1.0);
786  int nt=G.nt(),ntau=G.ntau();
787  int size=G.size1();
788  int sign=G.sig();
789  double tau,t,dtau=beta/ntau;
790  cdmatrix H0(size,size),H1(size,size);
791  eps.get_value(-1,H0);
792  H1=mu*cdmatrix::Identity(size,size)-H0;
793  cntr::function<T> Ut(nt,size);
794  cdmatrix evec0(size,size),value(size,size);
795  dvector eval0(size),eval0m(size);
796  Eigen::SelfAdjointEigenSolver<cdmatrix> eigensolver(H1);
797  evec0=eigensolver.eigenvectors();
798  eval0=eigensolver.eigenvalues();
799  eval0m=(-1.0)*eval0;
800 
801 
802  for(int m=0;m<=ntau;m++){
803  tau=m*dtau;
804  if(sign==-1){
805  value=(-1.0)*evec0*fermi_exp(beta,tau,eval0).asDiagonal()*evec0.adjoint();
806  }else if(sign==1){
807  value=evec0*bose_exp(beta,tau,eval0).asDiagonal()*evec0.adjoint();
808  }
809  G.set_mat(m,value);
810  }
811 
812  if(nt >= 0){
813  cdmatrix idm(size,size);
814  idm = MatrixXcd::Identity(size,size);
815  Ut.set_value(-1,idm);
816  Ut.set_value(0,idm);
817  for(int tstp=1;tstp<=nt;tstp++){
818  propagator_exp(tstp,Ut,eps,h,order,kt,fixHam);
819  }
820  for(int tstp=0;tstp<=nt;tstp++){
821  cdmatrix tmp;
822  Ut.get_value(tstp,tmp);
823  tmp = tmp * std::complex<T>(cos(mu * h * tstp),sin(mu * h * tstp));
824  Ut.set_value(tstp,tmp);
825  }
827 
828  cdmatrix expp(size,size);
829  for(int m=0;m<=ntau;m++){
830  tau=m*dtau;
831  for(int n=0;n<=nt;n++){
832  Ut.get_value(n,expp);
833  if(sign==-1){
834  value=iu*expp*evec0*fermi_exp(beta,tau,eval0m).asDiagonal()*evec0.adjoint();
835  }else if(sign==1){
836  value=-iu*expp*evec0*bose_exp(beta,tau,eval0m).asDiagonal()*evec0.adjoint();
837  }
838  G.set_tv(n,m,value);
839  }
840  }
841  if(sign==-1){
842  value=evec0*fermi(beta,eval0m).asDiagonal()*evec0.adjoint();
843  }else if(sign==1){
844  value=-1.0*evec0*bose(beta,eval0m).asDiagonal()*evec0.adjoint();
845  }
846  cdmatrix exppt1(size,size);
847  cdmatrix exppt2(size,size);
848  for(int m=0;m<=nt;m++){
849  for(int n=0;n<=m;n++){
850  cdmatrix tmp(size,size);
851  Ut.get_value(m,exppt1);
852  Ut.get_value(n,exppt2);
853  tmp = -iu*exppt1*exppt2.adjoint();
854  G.set_ret(m,n,tmp);
855  tmp=iu*exppt2*value*exppt1.adjoint();
856  G.set_les(n,m,tmp);
857  }
858  }
859  }
860 }
861 
863 template<typename T,int SIZE>
864 void green_from_H_const_dispatch(herm_matrix<T> &G,T mu,cdmatrix &H0,T beta,T h){
865  std::complex<T> iu = std::complex<T>(0.0, 1.0);
866  int nt=G.nt(),ntau=G.ntau();
867  int size=G.size1();
868  int sign=G.sig();
869  double tau,t,dtau=beta/ntau;
870  cdmatrix idm(size,size);
871  cdmatrix Udt(size,size);
872  cdmatrix IHdt(size,size);
873  cdmatrix Hmu(size,size);
874  cdmatrix evec0(size,size),value(size,size);
875  dvector eval0(size),eval0m(size);
876 
877  idm = MatrixXcd::Identity(size,size);
878  Hmu = -H0 + mu * idm;
879 
880  Eigen::SelfAdjointEigenSolver<cdmatrix> eigensolver(Hmu);
881  evec0=eigensolver.eigenvectors();
882  eval0=eigensolver.eigenvalues();
883  eval0m=(-1.0)*eval0;
884 
885  for(int m=0;m<=ntau;m++){
886  tau=m*dtau;
887  if(sign==-1){
888  value=(-1.0)*evec0*fermi_exp(beta,tau,eval0).asDiagonal()*evec0.adjoint();
889  }else if(sign==1){
890  value=(1.0)*evec0*bose_exp(beta,tau,eval0).asDiagonal()*evec0.adjoint();
891  }
892  G.set_mat(m,value);
893  }
894 
895  if(nt >=0 ){
896  // IHdt = -iu * h * H0;
897  IHdt = iu * h * Hmu;
898  Udt = IHdt.exp();
899 
900  cntr::function<T> Ut(nt,size);
901  cdmatrix Un(size,size);
902  Ut.set_value(-1,idm);
903  Ut.set_value(0,idm);
904  for(int n=1;n<=nt;n++){
905  Ut.get_value(n-1,Un);
906  Un = Un * Udt;
907  Ut.set_value(n,Un);
908  }
909 
910  cdmatrix expp(size,size);
911  for(int m=0;m<=ntau;m++){
912  tau=m*dtau;
913  for(int n=0;n<=nt;n++){
914  Ut.get_value(n,expp);
915  if(sign==-1){
916  value=iu*expp*evec0*fermi_exp(beta,tau,eval0m).asDiagonal()*evec0.adjoint();
917  }else if(sign==1){
918  value=-iu*expp*evec0*bose_exp(beta,tau,eval0m).asDiagonal()*evec0.adjoint();
919  }
920  G.set_tv(n,m,value);
921  }
922  }
923 
924  if(sign==-1){
925  value=evec0*fermi(beta,eval0m).asDiagonal()*evec0.adjoint();
926  }else if(sign==1){
927  value=-1.0*evec0*bose(beta,eval0m).asDiagonal()*evec0.adjoint();
928  }
929  cdmatrix exppt1(size,size);
930  cdmatrix exppt2(size,size);
931  for(int m=0;m<=nt;m++){
932  for(int n=0;n<=m;n++){
933  cdmatrix tmp(size,size);
934  Ut.get_value(m,exppt1);
935  Ut.get_value(n,exppt2);
936  tmp = -iu*exppt1*exppt2.adjoint();
937  G.set_ret(m,n,tmp);
938  tmp = iu*exppt2*value*exppt1.adjoint();
939  G.set_les(n,m,tmp);
940  }
941  }
942  }
943 }
944 
946 template<typename T,int SIZE>
947 void green_from_H_dispatch(herm_matrix_timestep<T> &G,T mu,cntr::function<T> &eps,T beta,T h,int kt,int order,bool fixHam=false){
948  std::complex<T> iu = std::complex<T>(0.0, 1.0);
949  int tstp=G.tstp(),ntau=G.ntau();
950  int size=G.size1();
951  int sign=G.sig();
952  cdmatrix H0(size,size),H1(size,size);
953  eps.get_value(-1,H0);
954  H1=mu*cdmatrix::Identity(size,size)-H0;
955  cntr::function<T> Ut(tstp,size);
956  cdmatrix evec0(size,size),value(size,size);
957  dvector eval0(size),eval0m(size);
958  Eigen::SelfAdjointEigenSolver<cdmatrix> eigensolver(H1);
959  evec0=eigensolver.eigenvectors();
960  eval0=eigensolver.eigenvalues();
961  eval0m=(-1.0)*eval0;
962 
963  if(tstp==-1){
964  for(int m=0;m<=ntau;m++){
965  double tau,t,dtau=beta/ntau;
966  tau=m*dtau;
967  if(sign==-1){
968  value=(-1.0)*evec0*fermi_exp(beta,tau,eval0).asDiagonal()*evec0.adjoint();
969  }else if(sign==1){
970  value=evec0*bose_exp(beta,tau,eval0).asDiagonal()*evec0.adjoint();
971  }
972  G.set_mat(m,value);
973  }
974  }else{
975  cdmatrix expp(size,size);
976  cdmatrix exppt1(size,size);
977  cdmatrix exppt2(size,size);
978  cdmatrix idm(size,size);
979  idm = MatrixXcd::Identity(size,size);
980  Ut.set_value(-1,idm);
981  Ut.set_value(0,idm);
982  for(int t=1;t<=tstp;t++){
983  propagator_exp(t,Ut,eps,h,order,kt,fixHam);
984  }
985 
986  for(int n=0; n<=tstp;n++){
987  cdmatrix tmp;
988  Ut.get_value(tstp,tmp);
989  tmp = tmp * std::complex<T>(cos(mu * h * n),sin(mu * h * n));
990  Ut.set_value(tstp,tmp);
991  }
992 
993 
994  // TV
995  for(int m=0;m<=ntau;m++){
996  double tau,dtau=beta/ntau;
997  tau=m*dtau;
998  Ut.get_value(tstp,expp);
999  if(sign==-1){
1000  value=iu*expp*evec0*fermi_exp(beta,tau,eval0m).asDiagonal()*evec0.adjoint();
1001  }else if(sign==1){
1002  value=-iu*expp*evec0*bose_exp(beta,tau,eval0m).asDiagonal()*evec0.adjoint();
1003  }
1004  G.set_tv(m,value);
1005  }
1006  // RET
1007  for(int n=0;n<=tstp;n++){
1008  cdmatrix tmp(size,size);
1009  Ut.get_value(tstp,exppt1);
1010  Ut.get_value(n,exppt2);
1011  tmp = -iu*exppt1*exppt2.adjoint();
1012  G.set_ret(n,tmp);
1013  }
1014  // LES
1015  if(sign==-1){
1016  value=evec0*fermi(beta,eval0m).asDiagonal()*evec0.adjoint();
1017  }else if(sign==1){
1018  value=-1.0*evec0*bose(beta,eval0m).asDiagonal()*evec0.adjoint();
1019  }
1020  for(int n=0;n<=tstp;n++){
1021  cdmatrix tmp(size,size);
1022  Ut.get_value(tstp,exppt1);
1023  Ut.get_value(n,exppt2);
1024  tmp = iu*exppt2*value*exppt1.adjoint();
1025  G.set_les(n,tmp);
1026  }
1027  }
1028 }
1029 
1031 template<typename T,int SIZE>
1032 void green_from_H_const_dispatch(herm_matrix_timestep<T> &G,T mu,cdmatrix &H0,T beta,T h){
1033  std::complex<T> iu = std::complex<T>(0.0, 1.0);
1034  int tstp=G.tstp(),ntau=G.ntau();
1035  int size=G.size1();
1036  int sign=G.sig();
1037  cdmatrix evec0(size,size),value(size,size);
1038  dvector eval0(size),eval0m(size);
1039  cdmatrix H1;
1040  H1 = mu*cdmatrix::Identity(size,size)-H0;
1041 
1042  Eigen::SelfAdjointEigenSolver<cdmatrix> eigensolver(H1);
1043  evec0=eigensolver.eigenvectors();
1044  eval0=eigensolver.eigenvalues();
1045  eval0m=(-1.0)*eval0;
1046 
1047  if(tstp==-1){
1048  for(int m=0;m<=ntau;m++){
1049  double tau,t,dtau=beta/ntau;
1050  tau=m*dtau;
1051  if(sign==-1){
1052  value=(-1.0)*evec0*fermi_exp(beta,tau,eval0).asDiagonal()*evec0.adjoint();
1053  }else if(sign==1){
1054  value=evec0*bose_exp(beta,tau,eval0).asDiagonal()*evec0.adjoint();
1055  }
1056  G.set_mat(m,value);
1057  }
1058  }else{
1059  cdmatrix Udt(size,size);
1060  cdmatrix IHdt(size,size);
1061  cdmatrix expp(size,size);
1062  cdmatrix exppt1(size,size);
1063  cdmatrix exppt2(size,size);
1064 
1065  // IHdt = -iu * h * H0;
1066  IHdt = iu * h * H1;
1067  Udt = IHdt.exp();
1068 
1069  cntr::function<T> Ut(tstp,size);
1070  cdmatrix idm(size,size);
1071  idm = MatrixXcd::Identity(size,size);
1072  Ut.set_value(-1,idm);
1073  Ut.set_value(0,idm);
1074  for(int n=1;n<=tstp;n++){
1075  cdmatrix Un(size,size);
1076  Ut.get_value(n-1,Un);
1077  Un = Un * Udt;
1078  Ut.set_value(n,Un);
1079  }
1080  // TV
1081  for(int m=0;m<=ntau;m++){
1082  double tau,dtau=beta/ntau;
1083  tau=m*dtau;
1084  Ut.get_value(tstp,expp);
1085  if(sign==-1){
1086  value=iu*expp*evec0*fermi_exp(beta,tau,eval0m).asDiagonal()*evec0.adjoint();
1087  }else if(sign==1){
1088  value=-iu*expp*evec0*bose_exp(beta,tau,eval0m).asDiagonal()*evec0.adjoint();
1089  }
1090  G.set_tv(m,value);
1091  }
1092  // RET
1093  for(int n=0;n<=tstp;n++){
1094  cdmatrix tmp(size,size);
1095  Ut.get_value(tstp,exppt1);
1096  Ut.get_value(n,exppt2);
1097  tmp = -iu * (exppt1 * exppt2.adjoint());
1098  G.set_ret(n,tmp);
1099  }
1100  // LES
1101  if(sign==-1){
1102  value=evec0*fermi(beta,eval0m).asDiagonal()*evec0.adjoint();
1103  }else if(sign==1){
1104  value=-1.0*evec0*bose(beta,eval0m).asDiagonal()*evec0.adjoint();
1105  }
1106  for(int n=0;n<=tstp;n++){
1107  cdmatrix tmp(size,size);
1108  Ut.get_value(tstp,exppt1);
1109  Ut.get_value(n,exppt2);
1110  tmp = iu * (exppt2 * value * exppt1.adjoint());
1111  G.set_les(n,tmp);
1112  }
1113  }
1114 }
1115 
1141 template<typename T>
1142 void green_from_H(herm_matrix<T> &G,T mu,cdmatrix &eps,T beta,T h){
1143  assert(G.size1()==eps.rows());
1144  assert(eps.rows()==eps.cols());
1145 
1146  int size=G.size1();
1147  if(size==1) green_from_H_const_dispatch<T,1>(G,mu,eps,beta,h);
1148  else green_from_H_const_dispatch<T,LARGESIZE>(G,mu,eps,beta,h);
1149 }
1150 
1152 template<typename T>
1153 void green_from_H(herm_matrix<T> &G,T mu,cntr::function<T> &eps,T beta,T h,bool fixHam,int SolveOrder,int cf_order){
1154  assert(G.size1()==eps.size2_);
1155  assert(eps.size1_==eps.size2_);
1156  assert(SolveOrder <= MAX_SOLVE_ORDER);
1157  assert(cf_order == 2 || cf_order == 4);
1158 
1159  int size=G.size1();
1160  if(size==1) green_from_H_dispatch<T,1>(G,mu,eps,beta,h,SolveOrder,cf_order,fixHam);
1161  else green_from_H_dispatch<T,LARGESIZE>(G,mu,eps,beta,h,SolveOrder,cf_order,fixHam);
1162 }
1163 
1164 
1195 template<typename T>
1196 void green_from_H(herm_matrix<T> &G,T mu,cntr::function<T> &eps,T beta,T h,int SolveOrder,int cf_order){
1197  assert(G.size1()==eps.size2_);
1198  assert(eps.size1_==eps.size2_);
1199  assert(SolveOrder <= MAX_SOLVE_ORDER);
1200  assert(cf_order == 2 || cf_order == 4);
1201 
1202  int size=G.size1();
1203  if(size==1) green_from_H_dispatch<T,1>(G,mu,eps,beta,h,SolveOrder,cf_order,true);
1204  else green_from_H_dispatch<T,LARGESIZE>(G,mu,eps,beta,h,SolveOrder,cf_order,true);
1205 }
1206 
1207 
1208 
1210 template<typename T>
1211 void green_from_H(herm_matrix_timestep<T> &G,T mu,cdmatrix &eps,T beta,T h){
1212  assert(G.size1()==eps.rows());
1213  assert(eps.rows()==eps.cols());
1214 
1215  int size=G.size1();
1216  if(size==1) green_from_H_const_dispatch<T,1>(G,mu,eps,beta,h);
1217  else green_from_H_const_dispatch<T,LARGESIZE>(G,mu,eps,beta,h);
1218 }
1219 
1220 
1248 template<typename T>
1249 void green_from_H(int tstp, herm_matrix_timestep<T> &G,T mu,cdmatrix &eps,T beta,T h){
1250  assert(tstp == G.tstp());
1251  assert(G.size1()==eps.rows());
1252  assert(eps.rows()==eps.cols());
1253 
1254  int size=G.size1();
1255  if(size==1) green_from_H_const_dispatch<T,1>(G,mu,eps,beta,h);
1256  else green_from_H_const_dispatch<T,LARGESIZE>(G,mu,eps,beta,h);
1257 }
1258 
1259 
1287 template<typename T>
1288 void green_from_H(int tstp, herm_matrix<T> &G,T mu,cdmatrix &eps,T beta,T h){
1289  assert(G.size1()==eps.rows());
1290  assert(eps.rows()==eps.cols());
1291  herm_matrix_timestep<T> Gstep(tstp,G.ntau(),G.size1(),G.sig());
1292 
1293  int size=G.size1();
1294  if(size==1) green_from_H_const_dispatch<T,1>(Gstep,mu,eps,beta,h);
1295  else green_from_H_const_dispatch<T,LARGESIZE>(Gstep,mu,eps,beta,h);
1296 
1297  G.set_timestep(tstp, Gstep);
1298 }
1299 
1300 
1302 template<typename T>
1303 void green_from_H(herm_matrix_timestep<T> &G,T mu,cntr::function<T> &eps,T beta,T h,bool fixHam,int SolveOrder,int cf_order){
1304  assert(G.size1()==eps.size2_);
1305  assert(eps.size1_==eps.size2_);
1306  assert(SolveOrder <= MAX_SOLVE_ORDER);
1307  assert(cf_order == 2 || cf_order == 4);
1308 
1309  int size=G.size1();
1310  if(size==1) green_from_H_dispatch<T,1>(G,mu,eps,beta,h,SolveOrder,cf_order,fixHam);
1311  else green_from_H_dispatch<T,LARGESIZE>(G,mu,eps,beta,h,SolveOrder,cf_order,fixHam);
1312 }
1313 
1347 template<typename T>
1348 void green_from_H(int tstp, herm_matrix_timestep<T> &G,T mu,cntr::function<T> &eps,T beta,T h,bool fixHam,int SolveOrder,int cf_order){
1349  assert(tstp == G.tstp());
1350  assert(G.size1()==eps.size2_);
1351  assert(eps.size1_==eps.size2_);
1352  assert(SolveOrder <= MAX_SOLVE_ORDER);
1353  assert(cf_order == 2 || cf_order == 4);
1354 
1355  int size=G.size1();
1356  if(size==1) green_from_H_dispatch<T,1>(G,mu,eps,beta,h,SolveOrder,cf_order,fixHam);
1357  else green_from_H_dispatch<T,LARGESIZE>(G,mu,eps,beta,h,SolveOrder,cf_order,fixHam);
1358 }
1359 
1393 template<typename T>
1394 void green_from_H(int tstp, herm_matrix<T> &G,T mu,cntr::function<T> &eps,T beta,T h,bool fixHam,int SolveOrder,int cf_order){
1395  assert(tstp <= G.nt());
1396  assert(G.size1()==eps.size2_);
1397  assert(eps.size1_==eps.size2_);
1398  assert(SolveOrder <= MAX_SOLVE_ORDER);
1399  assert(cf_order == 2 || cf_order == 4);
1400  herm_matrix_timestep<T> Gstep(tstp,G.ntau(),G.size1(),G.sig());
1401 
1402  int size=G.size1();
1403  if(size==1) green_from_H_dispatch<T,1>(Gstep,mu,eps,beta,h,SolveOrder,cf_order,fixHam);
1404  else green_from_H_dispatch<T,LARGESIZE>(Gstep,mu,eps,beta,h,SolveOrder,cf_order,fixHam);
1405 
1406  G.set_timestep(tstp, Gstep);
1407 }
1408 
1409 
1411 // BOSONIC GREENS FUNCTION:
1412 // G(t,t') = - ii * <TC b(t) bdag(t') >
1413 // with H = w * b^dag b
1414 /*template<typename T>
1415 void green_single_pole_bose(herm_matrix<T> &G,T *eps,T beta,T h){
1416  int nt=G.nt(),size1=G.size1();
1417  int ntau=G.ntau();
1418  int n,m,i0=0;
1419  double tau,fm,eps1;
1420  std::complex<T> x,ii=std::complex<T>(0,1.0);
1421  std::vector<std::complex<T> > expp;
1422  // expp[n] = -exp(ii*w*n)
1423  expp.resize(nt+1);
1424  G.clear();
1425  double dtau=beta/ntau;
1426  for(int i=0;i<size1;i++){
1427  i0=i*size1+i;
1428  eps1=-eps[i];
1429  for(m=0;m<=nt;m++){
1430  expp[m]=std::complex<T>(cos(eps1*m*h),sin(eps1*m*h));
1431  }
1432  fm=bose<T>(beta,-eps1);
1433  for(int m=0;m<=ntau;m++){
1434  tau=m*dtau;
1435  G.matptr(m)[i0]=bose_exp<T>(beta,tau,eps1);
1436  x=-ii*bose_exp<T>(beta,tau,-eps1);
1437  for(n=0;n<=nt;n++) G.tvptr(n,m)[i0]=x*expp[n];
1438  }
1439  x=bose<T>(beta,-eps1);
1440  for(m=0;m<=nt;m++){
1441  for(n=0;n<=m;n++){
1442  G.retptr(m,n)[i0]=-ii*expp[m-n];
1443  G.lesptr(n,m)[i0]=-ii*x*conj(expp[m-n]);
1444  }
1445  }
1446  }
1447 }
1448 */
1449 
1451 template<typename T>
1452 void green_single_pole_XX_timestep(int tstp,int ntau,int size1,std::complex<T> *ret,std::complex<T> *tv,std::complex<T> *les,T w,T beta,T h) {
1453 
1454  double bw = bose(beta,w);
1455  double c,s,ewv,emwv;
1456  int s2=size1*size1;
1457  std::complex<T> d_les, d_gtr, d_ret, d_tv;
1458  std::complex<T> i = std::complex<T>(0.0, 1.0);
1459  double mat_axisunit = beta/ntau;
1460  assert(w*beta>0);
1461 
1462  for(int tp=0;tp<=tstp;tp++){
1463  c = cos(w*(tstp-tp)*h);
1464  s = sin(w*(tstp-tp)*h);
1465  d_les = -i * ( (c-i*s) + 2.0 * c * bw ); // D_les(tp, t)
1466  d_gtr = -i * ( (c-i*s) + 2.0 * c * bw ); // D_gtr(t, tp)
1467  d_ret = d_gtr + conj(d_les); // D_ret(t, tp)
1468  les[tp*s2]=0.5*d_les;
1469  ret[tp*s2]=0.5*d_ret;
1470  }
1471  for(int v=0;v<=ntau;v++) {
1472  c = cos(w*tstp*h);
1473  s = sin(w*tstp*h);
1474  ewv = exp(w*v*mat_axisunit);
1475  emwv = exp(-w*v*mat_axisunit);
1476  d_tv = -i * ( (c+i*s)*emwv + ((c+i*s)*emwv + (c-i*s)*ewv)*bw ); // ok
1477  tv[v*s2]=0.5*d_tv;
1478  }
1479 }
1481 template<typename T>
1482 void green_single_pole_XX_mat(int ntau,int size1,std::complex<T> *mat,T w,T beta) {
1483  double ewv,emwv;
1484  int s2=size1*size1;
1485  std::complex<T> d_mat;
1486  std::complex<T> i = std::complex<T>(0.0, 1.0);
1487  double mat_axisunit = beta/ntau;
1488  assert(w*beta>0);
1489 
1490  for(int v=0;v<=ntau;v++) {
1491  ewv=bose_exp(beta,v*mat_axisunit,w); // exp(tau*w)b(w)
1492  emwv=-bose_exp(beta,v*mat_axisunit,-w); // -exp(-tau*w)b(-w)
1493  d_mat = -0.5*(ewv+emwv);
1494  mat[v*s2]=d_mat;
1495  }
1496 }
1498 template<typename T>
1499 void green_single_pole_XX_timestep(herm_matrix_timestep<T> &D0,T w,T beta,T h) {
1500  int tstp=D0.tstp_;
1501  if(tstp==-1) green_single_pole_XX_mat(D0.ntau(),D0.size1(),D0.matptr(0),w,beta);
1502  else green_single_pole_XX_timestep(tstp,D0.ntau(),D0.size1(),D0.retptr(0),D0.tvptr(0),D0.lesptr(0),w,beta,h);
1503 }
1505 template<typename T>
1506 void green_single_pole_XX_timestep(int tstp, herm_matrix_timestep<T> &D0,T w,T beta,T h) {
1507  assert(tstp == D0.tstp());
1508  if(tstp==-1) green_single_pole_XX_mat(D0.ntau(),D0.size1(),D0.matptr(0),w,beta);
1509  else green_single_pole_XX_timestep(tstp,D0.ntau(),D0.size1(),D0.retptr(0),D0.tvptr(0),D0.lesptr(0),w,beta,h);
1510 }
1512 template<typename T>
1513 void green_single_pole_XX_timestep(int tstp,herm_matrix<T> &D0,T w,T beta,T h) {
1514  if(tstp==-1) green_single_pole_XX_mat(D0.ntau(),D0.size1(),D0.matptr(0),w,beta);
1515  else green_single_pole_XX_timestep(tstp,D0.ntau(),D0.size1(),D0.retptr(tstp,0),D0.tvptr(tstp,0),D0.lesptr(0,tstp),w,beta,h);
1516 }
1517 
1518 
1540 template<typename T>
1541 void green_single_pole_XX(herm_matrix<T> &D0,T w,T beta,T h) {
1542  int nt=D0.nt(),t;
1543  if(nt>=-1) green_single_pole_XX_mat(D0.ntau(),D0.size1(),D0.matptr(0),w,beta);
1544  for(t=0;t<=nt;t++){
1545  green_single_pole_XX_timestep(t,D0.ntau(),D0.size1(),D0.retptr(t,0),D0.tvptr(t,0),D0.lesptr(0,t),w,beta,h);
1546  }
1547 }
1548 
1549 
1550 } // namespace cntr
1551 
1552 #endif // CNTR_EQUILIBRIUM_IMPL_H
void green_from_H(herm_matrix< T > &G, T mu, cdmatrix &eps, T beta, T h)
Propagator for time-independent free Hamiltonian
int size2_
Number of the rows in the Matrix form.
Class herm_matrix_timestep deals with contour objects at a particular timestep .
double operator()(double omega)
void green_equilibrium_mat(herm_matrix< T > &G, dos_function &dos, double beta, int limit=100, int nn=20, double mu=0.0)
Class adft_func contains functions for computing accurate Fourier transforms by adaptive splitting o...
Definition: fourier.hpp:42
double lo_
Internal value determining the higher edge of the Fourier transform region
double hi_
Higher edge of the density of states
double operator()(double x)
T bose(T beta, T omega)
Evaluates the Bose distribution function.
void green_single_pole_XX_timestep(int tstp, herm_matrix_timestep< T > &D0, T w, T beta, T h)
std::complex< double > cplx
Definition: fourier.cpp:11
void green_equilibrium_mat_bethe(herm_matrix< T > &G, double beta, int limit=100, int nn=20, double mu=0.0)
Equilibrium propagator for Bethe semicircular density of states. Matsubara only. ...
void green_equilibrium_bethe(herm_matrix< T > &G, double beta, double h, int limit=100, int nn=20, double mu=0.0)
Equilibrium propagator for Bethe semicircular density of states
double N_
Smoothness parameter of the box
void sample(double w, double a, double b, function &fz, int n, int limit)
Fourier samples the given function to determine the individual intervals.
Definition: fourier.hpp:366
int size1_
Number of the colums in the Matrix form.
Class function for objects with time on real axis.
void get_value(int tstp, EigenMatrix &M) const
Get matrix value of this function object at a specific time point
void interpolate_CF2(int tstp, cntr::function< T > &H, cdmatrix &Hinter, int kt, bool fixHam=false)
Interpolation for the second order propagator
Class &#39;smooth_box` represent a smooth box-like density of states
smooth_box(double a, double b, double nu)
Class herm_matrix for two-time contour objects with hermitian symmetry.
double operator()(double x)
void dft(double w, cplx &result, cplx &err)
Computes the cubically corrected DFT.
Definition: fourier.hpp:317
double B_
Lower border of the box
int nt_
Maximum number of the time steps.
dos_wrapper(const dos &f, int sign, double mu=0.0)
void set_timestep(int tstp, herm_matrix &g1)
Sets all components at time step tstp to the components of a given herm_matrix. ...
T fermi(T beta, T omega)
Evaluates the Fermi distribution function.
void interpolate_CF4(int tstp, cntr::function< T > &H, cdmatrix &Hinte1, cdmatrix &Hinte2, int kt, bool fixHam=false)
Interpolation for the forth order propagator
void set_value(int tstp, EigenMatrix &M)
Set matrix value at a specific time point
void propagator_exp(int tstp, cntr::function< T > &U, cntr::function< T > &H, double dt, int order, int kt, bool fixHam=false)
Propagator for time-dependent free Hamiltonian
double V_
Hopping integral and 4V corresponds to the bandwidth
template Integrator< double > & I< double >(int k)
double lo_
Lower edge of the density of states
Class bethe represent the semicircular density of states