1 #ifndef CNTR_EQUILIBRIUM_IMPL_H 2 #define CNTR_EQUILIBRIUM_IMPL_H 46 if (fabs(arg) > EXPMAX) {
47 return (arg > 0.0 ? 0.0 : 1.0);
49 return 1.0 / (1.0 + exp(arg));
71 dvector
fermi(T beta,dvector &omega){
72 int size=omega.size();
74 for(
int i=0;i<size;i++){
75 tmp(i)=
fermi(beta,omega(i));
81 T fermi_exp(T beta, T tau, T omega) {
83 return exp(omega * tau) *
86 return exp((tau - beta) * omega) *
92 dvector fermi_exp(T beta,T tau,dvector &omega){
93 int size=omega.size();
96 for(
int i=0;i<size;i++){
97 tmp(i)=fermi_exp(beta,tau,omega(i));
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));
130 template <
typename T>
132 T arg = omega * beta;
134 return (-1.0 - bose<T>(beta, -omega));
135 if (fabs(arg) > EXPMAX) {
137 }
else if (arg < 1e-10) {
140 return 1.0 / (exp(arg) - 1.0);
163 dvector
bose(T beta,dvector &omega){
164 int size=omega.size();
166 for(
int i=0;i<size;i++){
167 tmp(i)=
bose(beta,omega(i));
175 template <
typename T>
176 T bose_exp(T beta, T tau, T omega) {
178 return exp(tau * omega) * bose<T>(beta, omega);
180 return -exp((tau - beta) * omega) * bose<T>(beta, -omega);
184 dvector bose_exp(T beta,T tau,dvector &omega){
185 int size=omega.size();
188 for(
int i=0;i<size;i++){
189 tmp(i)=bose_exp(beta,tau,omega(i));
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);
205 return 1.0 / (1.0 + exp(arg));
210 double distribution_eq(T beta,T omega,
int sign)
213 return green_cntr_full_equilibrium_fermi(beta,omega);
215 return bose(beta,omega);
220 double distribution_exp_eq(T beta,T tau,T omega,
int sign)
223 return fermi_exp(beta,tau,omega);
225 return bose_exp(beta,tau,omega);
233 enum flavor { ret, adv, mat, tv, vt, les, gtr };
251 double a=
A_(omega),fm;
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);
338 ohmic(
double omegac) {
343 double operator()(
double x) {
344 return x*x*std::exp(-x/omegac_)/(2.0*omegac_*omegac_*omegac_);
380 lo_ = A_ - 20.0 / nu_;
381 hi_ =
B_ + 20.0 / nu_;
388 if (nu_ > 100 || nu_ <= 0.1 || A_ >=
B_) {
389 std::cerr <<
"nu=" << nu_ <<
" in smooth_box " << std::endl;
392 lo_ = A_ - 20.0 / nu_;
393 hi_ =
B_ + 20.0 / nu_;
396 typedef std::complex<double>
cplx;
400 adft.
sample(0.0,
lo_, hi_, *
this, 20, 100);
401 adft.
dft(0.0, res, err);
407 double arg1 = fermi<double>(nu_, x -
B_);
408 double arg2 = fermi<double>(nu_, A_ - x);
409 return arg1 * arg2 /
N_;
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)
418 typedef std::complex<double>
cplx;
419 int nt=G.nt(),i,l,size1=G.size1();
424 dos_wrapper<dos_function> dos1(dos,sign,mu);
427 adft.
sample(0.0,dos.lo_,dos.hi_,dos1,nn,limit);
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));
437 template <
typename T,
class dos_function>
440 typedef std::complex<double>
cplx;
451 for(m=0;m<=ntau;m++){
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));
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)
462 typedef std::complex<double>
cplx;
463 int ntau=G.ntau(),nt=G.nt(),n,m,size1=G.size1();
468 dos_wrapper<dos_function> dos1(dos,sign,mu);
473 for(m=0;m<=ntau;m++){
475 adft.
sample(0.0,dos.lo_,dos.hi_,dos1,nn,limit);
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));
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)
487 typedef std::complex<double>
cplx;
488 int nt=G.nt(),i,l,size1=G.size1();
493 dos_wrapper<dos_function> dos1(dos,sign,mu);
497 adft.
sample(0.0,dos.lo_-0.0,dos.hi_-0.0,dos1,nn,limit);
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);
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)
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);
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)
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);
575 template <
typename T>
608 template <
typename T>
610 int limit,
int nn,
double mu) {
612 green_equilibrium(G, dos, beta, h, limit, nn,mu);
650 cdmatrix tmp1(size,size),tmp2(size,size);
651 int ktextrap=(tstp<=kt ? tstp : kt);
652 int n1=(tstp<=kt ? kt : tstp);
654 if(!fixHam && tstp>kt){
700 cdmatrix tmp(size,size);
701 int ktextrap=(tstp<=kt ? tstp : kt);
702 int n1=(tstp<=kt ? kt : tstp);
704 if(!fixHam && tstp>kt){
748 cdmatrix prop(size,size);
751 assert(order==2 || order==4);
754 if(tstp==-1 || tstp==0){
759 cdmatrix arg(size,size);
762 arg=arg*std::complex<double>(0.0,-1.0)*dt;
767 cdmatrix H1(size,size),H2(size,size);
768 cdmatrix arg1(size,size),arg2(size,size);
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);
776 prop=arg1.exp()*arg2.exp()*prop;
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();
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;
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();
802 for(
int m=0;m<=ntau;m++){
805 value=(-1.0)*evec0*fermi_exp(beta,tau,eval0).asDiagonal()*evec0.adjoint();
807 value=evec0*bose_exp(beta,tau,eval0).asDiagonal()*evec0.adjoint();
813 cdmatrix idm(size,size);
814 idm = MatrixXcd::Identity(size,size);
815 Ut.set_value(-1,idm);
817 for(
int tstp=1;tstp<=nt;tstp++){
820 for(
int tstp=0;tstp<=nt;tstp++){
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);
828 cdmatrix expp(size,size);
829 for(
int m=0;m<=ntau;m++){
831 for(
int n=0;n<=nt;n++){
832 Ut.get_value(n,expp);
834 value=iu*expp*evec0*fermi_exp(beta,tau,eval0m).asDiagonal()*evec0.adjoint();
836 value=-iu*expp*evec0*bose_exp(beta,tau,eval0m).asDiagonal()*evec0.adjoint();
842 value=evec0*
fermi(beta,eval0m).asDiagonal()*evec0.adjoint();
844 value=-1.0*evec0*
bose(beta,eval0m).asDiagonal()*evec0.adjoint();
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();
855 tmp=iu*exppt2*value*exppt1.adjoint();
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();
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);
877 idm = MatrixXcd::Identity(size,size);
878 Hmu = -H0 + mu * idm;
880 Eigen::SelfAdjointEigenSolver<cdmatrix> eigensolver(Hmu);
881 evec0=eigensolver.eigenvectors();
882 eval0=eigensolver.eigenvalues();
885 for(
int m=0;m<=ntau;m++){
888 value=(-1.0)*evec0*fermi_exp(beta,tau,eval0).asDiagonal()*evec0.adjoint();
890 value=(1.0)*evec0*bose_exp(beta,tau,eval0).asDiagonal()*evec0.adjoint();
901 cdmatrix Un(size,size);
902 Ut.set_value(-1,idm);
904 for(
int n=1;n<=nt;n++){
905 Ut.get_value(n-1,Un);
910 cdmatrix expp(size,size);
911 for(
int m=0;m<=ntau;m++){
913 for(
int n=0;n<=nt;n++){
914 Ut.get_value(n,expp);
916 value=iu*expp*evec0*fermi_exp(beta,tau,eval0m).asDiagonal()*evec0.adjoint();
918 value=-iu*expp*evec0*bose_exp(beta,tau,eval0m).asDiagonal()*evec0.adjoint();
925 value=evec0*
fermi(beta,eval0m).asDiagonal()*evec0.adjoint();
927 value=-1.0*evec0*
bose(beta,eval0m).asDiagonal()*evec0.adjoint();
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();
938 tmp = iu*exppt2*value*exppt1.adjoint();
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();
952 cdmatrix H0(size,size),H1(size,size);
953 eps.get_value(-1,H0);
954 H1=mu*cdmatrix::Identity(size,size)-H0;
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();
964 for(
int m=0;m<=ntau;m++){
965 double tau,t,dtau=beta/ntau;
968 value=(-1.0)*evec0*fermi_exp(beta,tau,eval0).asDiagonal()*evec0.adjoint();
970 value=evec0*bose_exp(beta,tau,eval0).asDiagonal()*evec0.adjoint();
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);
982 for(
int t=1;t<=tstp;t++){
986 for(
int n=0; n<=tstp;n++){
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);
995 for(
int m=0;m<=ntau;m++){
996 double tau,dtau=beta/ntau;
998 Ut.get_value(tstp,expp);
1000 value=iu*expp*evec0*fermi_exp(beta,tau,eval0m).asDiagonal()*evec0.adjoint();
1002 value=-iu*expp*evec0*bose_exp(beta,tau,eval0m).asDiagonal()*evec0.adjoint();
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();
1016 value=evec0*
fermi(beta,eval0m).asDiagonal()*evec0.adjoint();
1018 value=-1.0*evec0*
bose(beta,eval0m).asDiagonal()*evec0.adjoint();
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();
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();
1037 cdmatrix evec0(size,size),value(size,size);
1038 dvector eval0(size),eval0m(size);
1040 H1 = mu*cdmatrix::Identity(size,size)-H0;
1042 Eigen::SelfAdjointEigenSolver<cdmatrix> eigensolver(H1);
1043 evec0=eigensolver.eigenvectors();
1044 eval0=eigensolver.eigenvalues();
1045 eval0m=(-1.0)*eval0;
1048 for(
int m=0;m<=ntau;m++){
1049 double tau,t,dtau=beta/ntau;
1052 value=(-1.0)*evec0*fermi_exp(beta,tau,eval0).asDiagonal()*evec0.adjoint();
1054 value=evec0*bose_exp(beta,tau,eval0).asDiagonal()*evec0.adjoint();
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);
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);
1081 for(
int m=0;m<=ntau;m++){
1082 double tau,dtau=beta/ntau;
1084 Ut.get_value(tstp,expp);
1086 value=iu*expp*evec0*fermi_exp(beta,tau,eval0m).asDiagonal()*evec0.adjoint();
1088 value=-iu*expp*evec0*bose_exp(beta,tau,eval0m).asDiagonal()*evec0.adjoint();
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());
1102 value=evec0*
fermi(beta,eval0m).asDiagonal()*evec0.adjoint();
1104 value=-1.0*evec0*
bose(beta,eval0m).asDiagonal()*evec0.adjoint();
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());
1141 template<
typename T>
1143 assert(G.
size1()==eps.rows());
1144 assert(eps.rows()==eps.cols());
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);
1152 template<
typename T>
1154 assert(G.size1()==eps.
size2_);
1156 assert(SolveOrder <= MAX_SOLVE_ORDER);
1157 assert(cf_order == 2 || cf_order == 4);
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);
1195 template<
typename T>
1199 assert(SolveOrder <= MAX_SOLVE_ORDER);
1200 assert(cf_order == 2 || cf_order == 4);
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);
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());
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);
1248 template<
typename T>
1250 assert(tstp == G.
tstp());
1251 assert(G.
size1()==eps.rows());
1252 assert(eps.rows()==eps.cols());
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);
1287 template<
typename T>
1289 assert(G.
size1()==eps.rows());
1290 assert(eps.rows()==eps.cols());
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);
1302 template<
typename T>
1304 assert(G.size1()==eps.
size2_);
1306 assert(SolveOrder <= MAX_SOLVE_ORDER);
1307 assert(cf_order == 2 || cf_order == 4);
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);
1347 template<
typename T>
1349 assert(tstp == G.
tstp());
1352 assert(SolveOrder <= MAX_SOLVE_ORDER);
1353 assert(cf_order == 2 || cf_order == 4);
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);
1393 template<
typename T>
1395 assert(tstp <= G.
nt());
1398 assert(SolveOrder <= MAX_SOLVE_ORDER);
1399 assert(cf_order == 2 || cf_order == 4);
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);
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) {
1454 double bw =
bose(beta,w);
1455 double c,s,ewv,emwv;
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;
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 );
1466 d_gtr = -i * ( (c-i*s) + 2.0 * c * bw );
1467 d_ret = d_gtr + conj(d_les);
1468 les[tp*s2]=0.5*d_les;
1469 ret[tp*s2]=0.5*d_ret;
1471 for(
int v=0;v<=ntau;v++) {
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 );
1481 template<
typename T>
1482 void green_single_pole_XX_mat(
int ntau,
int size1,std::complex<T> *mat,T w,T beta) {
1485 std::complex<T> d_mat;
1486 std::complex<T> i = std::complex<T>(0.0, 1.0);
1487 double mat_axisunit = beta/ntau;
1490 for(
int v=0;v<=ntau;v++) {
1491 ewv=bose_exp(beta,v*mat_axisunit,w);
1492 emwv=-bose_exp(beta,v*mat_axisunit,-w);
1493 d_mat = -0.5*(ewv+emwv);
1498 template<
typename T>
1501 if(tstp==-1) green_single_pole_XX_mat(D0.ntau(),D0.size1(),D0.matptr(0),w,beta);
1505 template<
typename T>
1507 assert(tstp == D0.
tstp());
1508 if(tstp==-1) green_single_pole_XX_mat(D0.
ntau(),D0.
size1(),D0.
matptr(0),w,beta);
1512 template<
typename T>
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);
1540 template<
typename T>
1541 void green_single_pole_XX(herm_matrix<T> &D0,T w,T beta,T h) {
1543 if(nt>=-1) green_single_pole_XX_mat(D0.ntau(),D0.size1(),D0.matptr(0),w,beta);
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);
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...
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
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.
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 '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.
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