-
Notifications
You must be signed in to change notification settings - Fork 1
/
SC_BathMPO.h
113 lines (89 loc) · 4.77 KB
/
SC_BathMPO.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
inline void Fill_SCBath_MPO(MPO& H, const double Eshift, const std::vector<double>& eps_,
const std::vector<double>& v_, double epseff, double Ueff, const params &p)
{
//QN objects are necessary to have abelian symmetries in MPS
QN qn0 ( {"Sz", 0},{"Nf", 0} ),
cupC ( {"Sz", +1},{"Nf",+1} ),
cdnC ( {"Sz", -1},{"Nf",+1} ),
cupA ( {"Sz", -1},{"Nf",-1} ),
cdnA ( {"Sz", +1},{"Nf",-1} );
std::vector<Index> links;
links.push_back( Index() );
//first we create the link indices which carry quantum number information
for(auto i : range1(length(H))){
links.push_back(Index( qn0, 2,
cupC, 1,
cdnC, 1,
cupA, 1,
cdnA, 1,
cupC+cdnC, 1,
cupA+cdnA, 1, Out, "Link" ));
}
//then we just fill the MPO tensors which can be viewed
//as matrices (vectors) of operators. if one multiplies
//all matrices togehter one obtains the hamiltonian.
//therefore the tensor on the first and last site must be column/ row vectors
//and all sites between matrices
//first site is a vector:
{
int i = 1;
ITensor& W = H.ref(i);
Index right = links.at(i);
W = ITensor(right, p.sites.si(i), p.sites.siP(i) );
W += p.sites.op("Id",i) * setElt(right(1));
W += p.sites.op("Ntot",i) * setElt(right(2)) * epseff;
W += p.sites.op("Nup",i) * setElt(right(2)) * p.qd->EZ()/2.0;
W += p.sites.op("Ndn",i) * setElt(right(2)) * (-1) * p.qd->EZ()/2.0;
W += p.sites.op("Nupdn",i) * setElt(right(2)) * Ueff;
W += p.sites.op("Id",i) * setElt(right(2)) * Eshift;
W += p.sites.op("Cup*F",i) * setElt(right(3)) * (-1);
W += p.sites.op("Cdn*F",i) * setElt(right(4)) * (-1);
W += p.sites.op("Cdagup*F",i) * setElt(right(5)) * (+1);
W += p.sites.op("Cdagdn*F",i) * setElt(right(6)) * (+1);
}
// sites 2 ... N-1 are matrices
for(auto i : range1(2,length(H)-1)){
ITensor& W = H.ref(i);
Index left = dag( links.at(i-1) );
Index right = links.at(i);
W = ITensor(left, right, p.sites.si(i), p.sites.siP(i) );
W += p.sites.op("Id",i) * setElt(left(1), right(1));
W += p.sites.op("Ntot",i) * setElt(left(1),right(2)) * eps_[i-1];
W += p.sites.op("Nup",i) * setElt(left(1),right(2)) * p.sc->EZ()/2.0;
W += p.sites.op("Ndn",i) * setElt(left(1),right(2)) * (-1) * p.sc->EZ()/2.0;
W += p.sites.op("Nupdn",i) * setElt(left(1),right(2)) * p.sc->g() * pow(p.sc->y(i-1),2);
W += p.sites.op("Cdn*Cup",i) * setElt(left(1),right(7)) * p.sc->g() * p.sc->y(i-1);
W += p.sites.op("Cdagup*Cdagdn",i) * setElt(left(1),right(8)) * p.sc->g() * p.sc->y(i-1);
W += p.sites.op("Id",i)*setElt(left(2),right(2));
W += p.sites.op("F" ,i)*setElt(left(3),right(3));
W += p.sites.op("F" ,i)*setElt(left(4),right(4));
W += p.sites.op("F" ,i)*setElt(left(5),right(5));
W += p.sites.op("F" ,i)*setElt(left(6),right(6));
W += p.sites.op("Id",i)*setElt(left(7),right(7));
W += p.sites.op("Id",i)*setElt(left(8),right(8));
W += p.sites.op("Cdagup",i)*setElt(left(3),right(2))* v_[i-1];
W += p.sites.op("Cdagdn",i)*setElt(left(4),right(2))* v_[i-1];
W += p.sites.op("Cup", i)*setElt(left(5),right(2))* v_[i-1];
W += p.sites.op("Cdn", i)*setElt(left(6),right(2))* v_[i-1];
W += p.sites.op("Cdagup*Cdagdn",i)*setElt(left(7),right(2)) * p.sc->y(i-1);
W += p.sites.op("Cdn*Cup",i) *setElt(left(8),right(2)) * p.sc->y(i-1);
}
//site N is a vector again
{
int i = length(H);
ITensor& W = H.ref(i);
Index left = dag( links.at(i-1) );
W = ITensor(left, p.sites.si(i), p.sites.siP(i) );
W += p.sites.op("Ntot", i) * setElt(left(1)) * eps_[i-1];
W += p.sites.op("Nup", i) * setElt(left(1)) * p.sc->EZ()/2.0;
W += p.sites.op("Ndn", i) * setElt(left(1)) * (-1) * p.sc->EZ()/2.0;
W += p.sites.op("Nupdn",i) * setElt(left(1)) * p.sc->g() * pow(p.sc->y(i-1), 2);
W += p.sites.op("Id", i) * setElt(left(2)) ;
W += p.sites.op("Cdagup",i) * setElt(left(3)) * v_[i-1];
W += p.sites.op("Cdagdn",i) * setElt(left(4)) * v_[i-1];
W += p.sites.op("Cup", i) * setElt(left(5)) * v_[i-1];
W += p.sites.op("Cdn", i) * setElt(left(6)) * v_[i-1];
W += p.sites.op("Cdagup*Cdagdn",i) * setElt(left(7)) * p.sc->y(i-1);
W += p.sites.op("Cdn*Cup", i) * setElt(left(8)) * p.sc->y(i-1);
}
}