VMPS++
Loading...
Searching...
No Matches
HubbardKspace.h
Go to the documentation of this file.
1#ifndef HUBBARD_KSPACE
2#define HUBBARD_KSPACE
3
4template<typename Scalar>
5struct KspaceRawTerms
6{
7 // 2-site:
8 vector<tuple<int,int,Scalar> > spin_exchange;
9 vector<tuple<int,int,Scalar> > density_density; // only for U1xU1
10 vector<tuple<int,int,Scalar> > corr_hopping;
11 vector<tuple<int,int,Scalar> > corr_hoppingB;
12 vector<tuple<int,int,Scalar> > pair_hopping;
13
14 // 3-site:
15 vector<tuple<int,int,int,Scalar> > nonlocal_spin;
16 vector<tuple<int,int,int,Scalar> > nonlocal_spinB; // only for U1xU1
17 vector<tuple<int,int,int,Scalar> > corr_hopping3;
18 vector<tuple<int,int,int,Scalar> > corr_hopping3B;
19 vector<tuple<int,int,int,Scalar> > doublon_decay;
20
21 // 4-site
22 vector<tuple<int,int,int,int,Scalar> > foursite;
23
24 map<tuple<int,int,int,int>,Scalar> Umap;
25};
26
27template<typename MODEL>
28struct KspaceMpoTerms
29{
30 // 1-site:
31 vector<tuple<size_t,double> > HubbardU_kspace;
32
33 // 2-site:
34 vector<Mpo<typename MODEL::Symmetry,typename MODEL::Scalar_> > Hmpo_spin_exchange;
35 vector<Mpo<typename MODEL::Symmetry,typename MODEL::Scalar_> > Hmpo_density_density;
36 vector<Mpo<typename MODEL::Symmetry,typename MODEL::Scalar_> > Hmpo_pair_hopping;
37 vector<Mpo<typename MODEL::Symmetry,typename MODEL::Scalar_> > Hmpo_corr_hopping;
38
39 // 3-site:
40 vector<Mpo<typename MODEL::Symmetry,typename MODEL::Scalar_> > Hmpo_corr_hopping3;
41 vector<Mpo<typename MODEL::Symmetry,typename MODEL::Scalar_> > Hmpo_nonlocal_spin;
42 vector<Mpo<typename MODEL::Symmetry,typename MODEL::Scalar_> > Hmpo_doublon_decay;
43
44 // 4-site:
45 vector<Mpo<typename MODEL::Symmetry,typename MODEL::Scalar_> > Hmpo_foursite;
46};
47
48template<typename MODEL>
49struct KspaceHTerms
50{
51 // 1-site:
52 vector<tuple<size_t,double> > HubbardU_kspace;
53
54 // 2-site:
55 vector<MODEL> H2_spin_exchange;
56 vector<MODEL> H2_pair_hopping;
57 vector<MODEL> H2_corr_hopping;
58
59 // 3-site:
60 vector<MODEL> H3_corr_hopping3;
61 vector<MODEL> H3_nonlocal_spin;
62 vector<MODEL> H3_doublon_decay;
63
64 // 4-site:
65 vector<MODEL> H4;
66};
67
68template<typename MODEL>
69class HubbardKspace
70{
71public:
72
73 typedef Mpo<typename MODEL::Symmetry, typename MODEL::Scalar_> OPERATOR;
74
75 HubbardKspace(){};
76
77 // UU: unitary transformation of the hopping matrix
78 // U: Hubbard-U in real space
79 HubbardKspace (const MatrixXcd &UU_input, double U_input, DMRG::VERBOSITY::OPTION VERB_input=DMRG::VERBOSITY::SILENT, bool VUMPS_input=false, vector<int> x_input={}, vector<int> y_input={})
80 :UU(UU_input), U(U_input), VUMPS(VUMPS_input), x(x_input), y(y_input), VERB(VERB_input)
81 {
82 assert(UU.rows() == UU.cols());
83 L = static_cast<size_t>(UU.rows());
84
85 if (x.size() == 0)
86 {
87 x.resize(L);
88 y.resize(L);
89 for (int i=0; i<L; ++i)
90 {
91 x[i] = i;
92 y[i] = 0;
93 }
94 Lred = L;
95 }
96 else
97 {
98 Lred = *std::max_element(x.begin(), x.end())+1;
99 }
100
101 dummy_params.push_back({"maxPower",1ul});
102// typename MODEL::Scalar_ zero = 0.;
103// dummy.push_back({"t",zero});
104
105 Umap.clear();
106 compute_raw();
107 compute_MPO();
108 };
109
110 HubbardKspace (const MatrixXcd &UU_input, double U_input, const vector<Param> &params, DMRG::VERBOSITY::OPTION VERB_input=DMRG::VERBOSITY::SILENT, bool VUMPS_input=false, vector<int> x_input={}, vector<int> y_input={})
111 :UU(UU_input), U(U_input), dummy_params(params), VUMPS(VUMPS_input), x(x_input), y(y_input), VERB(VERB_input)
112 {
113 assert(UU.rows() == UU.cols());
114 L = static_cast<size_t>(UU.rows());
115
116 if (x.size() == 0)
117 {
118 x.resize(L);
119 y.resize(L);
120 for (int i=0; i<L; ++i)
121 {
122 x[i] = i;
123 y[i] = 0;
124 }
125 Lred = L;
126 }
127 else
128 {
129 Lred = *std::max_element(x.begin(), x.end())+1;
130 lout << "Lred=" << Lred << endl;
131 }
132
133 Umap.clear();
134 compute_raw();
135 compute_MPO();
136 };
137
138 string info() const;
139
140 template<class Dummy = typename MODEL::Symmetry>
141 typename std::enable_if<Dummy::IS_SPIN_SU2(),void>::type compute_raw();
142
143 template<class Dummy = typename MODEL::Symmetry>
144 typename std::enable_if<!Dummy::IS_SPIN_SU2(),void>::type compute_raw();
145
146 template<class Dummy = typename MODEL::Symmetry>
147 typename std::enable_if<Dummy::IS_SPIN_SU2(),void>::type compute_MPO();
148
149 template<class Dummy = typename MODEL::Symmetry>
150 typename std::enable_if<!Dummy::IS_SPIN_SU2(),void>::type compute_MPO();
151
152 KspaceHTerms<MODEL> get_Hterms() const {return Hterms;};
153
154 MODEL sum_all() const;
155 MODEL sum_all(const ArrayXXcd &hopping) const;
156 MODEL sum_2site() const;
157 MODEL sum_3site() const;
158 MODEL sum_4site() const;
159 Mpo<typename MODEL::Symmetry,typename MODEL::Scalar_> sum_all_mpo() const;
160 Mpo<typename MODEL::Symmetry,typename MODEL::Scalar_> sum_2site_mpo() const;
161
162private:
163
164 vector<Param> dummy_params;
165
166 DMRG::VERBOSITY::OPTION VERB;
167
168 size_t L;
169 size_t Lred;
170 MatrixXcd UU;
171 double U;
172 map<tuple<int,int,int,int>,typename MODEL::Scalar_> Umap;
173
174 bool VUMPS;
175
176 // If sites are blocked up, translate into x & y
177 vector<int> x;
178 vector<int> y;
179
180 KspaceRawTerms<typename MODEL::Scalar_> Raw;
181 KspaceMpoTerms<MODEL> Terms;
182 KspaceHTerms<MODEL> Hterms;
183};
184
185template<typename MODEL>
186string HubbardKspace<MODEL>::
187info() const
188{
189 stringstream ss;
190 ss << "HubbardKspace:" << endl;
191 ss << "#spin_exchange: " << Raw.spin_exchange.size()+Raw.density_density.size() << endl;
192 ss << "#corr_hopping: " << Raw.corr_hopping.size()+Raw.corr_hoppingB.size() << endl;
193 ss << "#pair_hopping: " << Raw.pair_hopping.size() << endl;
194 ss << "#nonlocal_spin: " << Raw.nonlocal_spin.size()+Raw.nonlocal_spinB.size() << endl;
195 ss << "#corr_hopping3: " << Raw.corr_hopping3.size()+Raw.corr_hopping3B.size() << endl;
196 ss << "#doublon_decay: " << Raw.doublon_decay.size() << endl;
197 ss << "#foursite: " << Raw.foursite.size() << endl;
198 ss << "#total: " << Raw.spin_exchange.size()+Raw.density_density.size()+Raw.corr_hopping.size()+Raw.corr_hoppingB.size()+Raw.pair_hopping.size()+Raw.nonlocal_spin.size()+Raw.nonlocal_spinB.size()+Raw.corr_hopping3.size()+Raw.corr_hopping3B.size()+Raw.doublon_decay.size()+Raw.foursite.size() <<endl;
199 return ss.str();
200}
201
202template<typename MODEL>
203template<typename Dummy>
204typename std::enable_if<Dummy::IS_SPIN_SU2(),void>::type HubbardKspace<MODEL>::
205compute_raw()
206{
207 vector<tuple<int,int,int,int> > terms_1site;
208 vector<tuple<int,int,int,int> > terms_2site;
209 vector<tuple<int,int,int,int> > terms_3site;
210 vector<tuple<int,int,int,int> > terms_4site;
211
212 // Carry out sum over local space, fill Umap and sort terms according to 1/2/3/4 site
213 for (int k=0; k<L; ++k)
214 for (int l=0; l<L; ++l)
215 for (int m=0; m<L; ++m)
216 for (int n=0; n<L; ++n)
217 {
218 std::set<int> s;
219 s.insert(k);
220 s.insert(l);
221 s.insert(m);
222 s.insert(n);
223
224 complex<double> Uelement_ = 0.;
225 for (int i=0; i<L; ++i)
226 {
227 Uelement_ += conj(UU(i,k))*UU(i,l)*conj(UU(i,m))*UU(i,n);
228 }
229 typename MODEL::Scalar_ Uelement;
230 if constexpr(is_same<typename MODEL::Scalar_,double>::value)
231 {
232 if (abs(Uelement_.imag()) > 1e-10)
233 {
234 lout << termcolor::red << "Warning: Non-zero imaginary part of transformed matrix element=" << Uelement_.imag() << termcolor::reset << endl;
235 }
236 Uelement = Uelement_.real();
237 }
238 else
239 {
240 Uelement = Uelement_;
241 }
242 if (abs(Uelement) > 1e-8)
243 {
244 auto Vklmn = make_tuple(k,l,m,n);
245 //if (k!=m and l!=n)
246 //{
247 // lout << "k=" << k << ", l=" << l << ", m=" << m << ", n=" << n << ", U=" << Uelement << endl;
248 //}
249 Umap[Vklmn] = U*Uelement;
250
251 if (s.size() == 1)
252 {
253 terms_1site.push_back(Vklmn);
254 }
255 else if (s.size() == 2)
256 {
257 terms_2site.push_back(Vklmn);
258 }
259 else if (s.size() == 3)
260 {
261 terms_3site.push_back(Vklmn);
262 }
263 else if (s.size() == 4)
264 {
265 terms_4site.push_back(Vklmn);
266 }
267 }
268 }
269
270 // Start filling Raw (not yet MPOs):
271 // Groups the 1/2/3/4-site contributions according to the various MPO terms
272
274 for (int t=0; t<terms_1site.size(); ++t)
275 {
276 size_t i = get<0>(terms_1site[t]);
277 Terms.HubbardU_kspace.push_back(make_tuple(i,real(Umap[make_tuple(i,i,i,i)])));
278 }
279 Hterms.HubbardU_kspace = Terms.HubbardU_kspace;
280
282 while (terms_2site.size() != 0)
283 {
284 for (int t=0; t<terms_2site.size(); ++t)
285 {
286 int k1 = get<0>(terms_2site[t]);
287 int k2 = get<1>(terms_2site[t]);
288 int k3 = get<2>(terms_2site[t]);
289 int k4 = get<3>(terms_2site[t]);
290
291 if (k1==k4 and k2==k3 and k1!=k2)
292 {
293 int i = k1;
294 int j = k2;
295
296 auto Vijji = make_tuple(i,j,j,i);
297 auto Vjiij = make_tuple(j,i,i,j);
298 auto Viijj = make_tuple(i,i,j,j);
299 auto Vjjii = make_tuple(j,j,i,i);
300
301 Raw.spin_exchange.push_back(make_tuple(i,j,Umap[Vijji]));
302 //lout << "i=" << i << ", j=" << j << ", spin_exchange: Vijji=" << Umap[Vijji] << ", Vjiij=" << Umap[Vjiij] << ", Viijj=" << Umap[Viijj] << ", Vjjii=" << Umap[Vjjii] << endl;
303
304 terms_2site.erase(std::remove(terms_2site.begin(), terms_2site.end(), Vijji), terms_2site.end());
305 terms_2site.erase(std::remove(terms_2site.begin(), terms_2site.end(), Vjiij), terms_2site.end());
306 terms_2site.erase(std::remove(terms_2site.begin(), terms_2site.end(), Viijj), terms_2site.end());
307 terms_2site.erase(std::remove(terms_2site.begin(), terms_2site.end(), Vjjii), terms_2site.end());
308 break;
309 }
310 else if (k1==k2 and k2==k3 and k1!=k4)
311 {
312 int i = k1;
313 int j = k4;
314
315 auto Viiij = make_tuple(i,i,i,j);
316 auto Vijii = make_tuple(i,j,i,i);
317
318 auto Viiji = make_tuple(i,i,j,i);
319 auto Vjiii = make_tuple(j,i,i,i);
320
321 Raw.corr_hopping.push_back(make_tuple(i,j,Umap[Viiij]));
322 //lout << "i=" << i << ", j=" << j << ", corr_hopping: Viiij=" << Umap[Viiij] << ", Vijii=" << Umap[Vijii] << ", Viiji" << Umap[Viiji] << ", Vjiii" << Umap[Vjiii] << endl;
323
324 terms_2site.erase(std::remove(terms_2site.begin(), terms_2site.end(), Vjiii), terms_2site.end());
325 terms_2site.erase(std::remove(terms_2site.begin(), terms_2site.end(), Vijii), terms_2site.end());
326 terms_2site.erase(std::remove(terms_2site.begin(), terms_2site.end(), Viiji), terms_2site.end());
327 terms_2site.erase(std::remove(terms_2site.begin(), terms_2site.end(), Viiij), terms_2site.end());
328 break;
329 }
330 else if (k1==k3 and k2==k4 and k1!=k2)
331 {
332 int i = k1;
333 int j = k2;
334
335 auto Vijij = make_tuple(i,j,i,j);
336 auto Vjiji = make_tuple(j,i,j,i);
337
338 Raw.pair_hopping.push_back(make_tuple(i,j,Umap[Vijij]));
339 //lout << "i=" << i << ", j=" << j << ", pair_hopping: Vijij=" << Umap[Vijij] << ", Vjiji=" << Umap[Vjiji] << endl;
340
341 terms_2site.erase(std::remove(terms_2site.begin(), terms_2site.end(), Vijij), terms_2site.end());
342 terms_2site.erase(std::remove(terms_2site.begin(), terms_2site.end(), Vjiji), terms_2site.end());
343 break;
344 }
345 }
346 }
347
349 while (terms_3site.size() != 0)
350 {
351 for (int t=0; t<terms_3site.size(); ++t)
352 {
353 int k1 = get<0>(terms_3site[t]);
354 int k2 = get<1>(terms_3site[t]);
355 int k3 = get<2>(terms_3site[t]);
356 int k4 = get<3>(terms_3site[t]);
357
358 if (k1==k2)
359 {
360 int i = k1;
361 int j = k3;
362 int k = k4;
363
364 auto Viijk = make_tuple(i,i,j,k);
365 auto Viikj = make_tuple(i,i,k,j);
366 auto Vjkii = make_tuple(j,k,i,i);
367 auto Vkjii = make_tuple(k,j,i,i);
368
369 Raw.corr_hopping3.push_back(make_tuple(i,j,k,Umap[Viijk]));
370 //lout << "i=" << i << ", j=" << j << ", k=" << k << ", corr_hopping3: Viijk=" << Umap[Viijk] << ", Viikj=" << Umap[Viikj] << ", Vjkii=" << Umap[Vjkii] << ", Vkjii=" << Umap[Vkjii] << endl;
371
372 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Viijk), terms_3site.end());
373 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Viikj), terms_3site.end());
374 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Vjkii), terms_3site.end());
375 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Vkjii), terms_3site.end());
376 break;
377 }
378 else if (k1==k3)
379 {
380 int i = k1;
381 int j = k2;
382 int k = k4;
383
384 auto Vijik = make_tuple(i,j,i,k);
385 auto Vikij = make_tuple(i,k,i,j);
386 auto Vjiki = make_tuple(j,i,k,i);
387 auto Vkiji = make_tuple(k,i,j,i);
388
389 Raw.doublon_decay.push_back(make_tuple(i,j,k,Umap[Vijik]));
390 //lout << "i=" << i << ", j=" << j << ", k=" << k << ", doublon_decay: Vijik=" << Umap[Vijik] << ", Vijik=" << Umap[Vikij] << ", Vjiki=" << Umap[Vjiki] << ", Vkiji=" << Umap[Vkiji] << endl;
391
392 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Vijik), terms_3site.end());
393 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Vikij), terms_3site.end());
394 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Vjiki), terms_3site.end());
395 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Vkiji), terms_3site.end());
396 break;
397 }
398 else if (k1==k4)
399 {
400 int i = k1;
401 int j = k2;
402 int k = k3;
403
404 auto Vijki = make_tuple(i,j,k,i);
405 auto Vkiij = make_tuple(k,i,i,j);
406 auto Vjiik = make_tuple(j,i,i,k);
407 auto Vikji = make_tuple(i,k,j,i);
408
409 Raw.nonlocal_spin.push_back(make_tuple(i,j,k,Umap[Vijki]));
410 //lout << "i=" << i << ", j=" << j << ", k=" << k << ", nonlocal_spin: Vijki=" << Umap[Vijki] << ", Vkiij=" << Umap[Vkiij] << ", Vjiik=" << Umap[Vjiik] << ", Vikji=" << Umap[Vikji] << endl;
411
412 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Vijki), terms_3site.end());
413 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Vkiij), terms_3site.end());
414 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Vjiik), terms_3site.end());
415 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Vikji), terms_3site.end());
416 break;
417 }
418 }
419 }
420
422 while (terms_4site.size() != 0)
423 {
424 for (int t=0; t<terms_4site.size(); ++t)
425 {
426 int k1 = get<0>(terms_4site[t]);
427 int k3 = get<1>(terms_4site[t]);
428 int k2 = get<2>(terms_4site[t]);
429 int k4 = get<3>(terms_4site[t]);
430
431 int i = k1;
432 int j = k2;
433 int k = k3;
434 int l = k4;
435
436 auto Vikjl = make_tuple(i,k,j,l);
437 auto Vjlik = make_tuple(j,l,i,k);
438 auto Vjkil = make_tuple(j,k,i,l);
439 auto Viljk = make_tuple(i,l,j,k);
440
441 auto Vkilj = make_tuple(k,i,l,j);
442 auto Vljki = make_tuple(l,j,k,i);
443 auto Vkjli = make_tuple(k,j,l,i);
444 auto Vlikj = make_tuple(l,i,k,j);
445
446 Raw.foursite.push_back(make_tuple(i,j,k,l,Umap[Vikjl]));
447 /*lout << "i=" << i << ", j=" << j << ", k=" << k << ", l=" << l
448 << ", 4s: "
449 << "Vikjl=" << Umap[Vikjl] << ", Vjlik=" << Umap[Vjlik] << ", Vjkil=" << Umap[Vjkil] << ", Viljk=" << Umap[Viljk]
450 << ", Vkilj=" << Umap[Vkilj] << ", Vljki=" << Umap[Vljki] << ", Vkjli=" << Umap[Vkjli] << ", Vkjli=" << Umap[Vkjli]
451 << endl;*/
452
453 terms_4site.erase(std::remove(terms_4site.begin(), terms_4site.end(), Vikjl), terms_4site.end());
454 terms_4site.erase(std::remove(terms_4site.begin(), terms_4site.end(), Vjlik), terms_4site.end());
455 terms_4site.erase(std::remove(terms_4site.begin(), terms_4site.end(), Vjkil), terms_4site.end());
456 terms_4site.erase(std::remove(terms_4site.begin(), terms_4site.end(), Viljk), terms_4site.end());
457
458 terms_4site.erase(std::remove(terms_4site.begin(), terms_4site.end(), Vkilj), terms_4site.end());
459 terms_4site.erase(std::remove(terms_4site.begin(), terms_4site.end(), Vljki), terms_4site.end());
460 terms_4site.erase(std::remove(terms_4site.begin(), terms_4site.end(), Vkjli), terms_4site.end());
461 terms_4site.erase(std::remove(terms_4site.begin(), terms_4site.end(), Vlikj), terms_4site.end());
462 break;
463 }
464 }
465}
466
467template<typename MODEL>
468template<typename Dummy>
469typename std::enable_if<Dummy::IS_SPIN_SU2(),void>::type HubbardKspace<MODEL>::
470compute_MPO()
471{
472 MODEL Hdummy(Lred,dummy_params,BC::OPEN,DMRG::VERBOSITY::SILENT);
473 //cout << Hdummy.info() << endl;
474
475 int s; // s controls how many terms are summed into bigger MPOs; default is: no additional summation
476
477 Terms.Hmpo_spin_exchange.resize(Raw.spin_exchange.size());
478 Terms.Hmpo_pair_hopping.resize(Raw.pair_hopping.size());
479 Terms.Hmpo_corr_hopping.resize(Raw.corr_hopping.size());
480
481 Hterms.H2_spin_exchange.resize(Raw.spin_exchange.size());
482 Hterms.H2_pair_hopping.resize(Raw.pair_hopping.size());
483 Hterms.H2_corr_hopping.resize(Raw.corr_hopping.size());
484
485 s = 0;
486 for (int t=0; t<Raw.spin_exchange.size(); ++t)
487 {
488 int i = get<0>(Raw.spin_exchange[t]);
489 int j = get<1>(Raw.spin_exchange[t]);
490 double lambda = real(get<2>(Raw.spin_exchange[t]));
491
492 if (VERB>0) lout << "spin exchange: " << i << ", " << j << ", λ=" << lambda << endl;
493
494 //auto SdagS = Hdummy.SdagS(i,j); SdagS.scale(-2.*lambda);
495 auto SdagS = Hdummy.SdagS(x[i],x[j],y[i],y[j]); SdagS.scale(-2.*lambda);
496 Terms.Hmpo_spin_exchange[s] = sum(Terms.Hmpo_spin_exchange[s],SdagS);
497
498 //auto nn = Hdummy.nn(i,j); nn.scale(0.5*lambda);
499 auto nn = Hdummy.nn(x[i],x[j],y[i],y[j]); nn.scale(0.5*lambda);
500 Terms.Hmpo_spin_exchange[s] = sum(Terms.Hmpo_spin_exchange[s],nn);
501 Hterms.H2_spin_exchange[s] = MODEL(Terms.Hmpo_spin_exchange[s],dummy_params);
502 s = (s+1)%Terms.Hmpo_spin_exchange.size();
503 }
504
505 s = 0;
506 for (int t=0; t<Raw.corr_hopping.size(); ++t)
507 {
508 int i = get<0>(Raw.corr_hopping[t]);
509 int j = get<1>(Raw.corr_hopping[t]);
510 typename MODEL::Scalar_ lambda = get<2>(Raw.corr_hopping[t]);
511
512 //auto T1 = Hdummy.cdagn_c(i,j); T1.scale(lambda);
513 //auto T2 = Hdummy.cdag_nc(j,i); T2.scale(conjIfcomplex(lambda));
514 auto T1 = Hdummy.cdagn_c(x[i],x[j],y[i],y[j]); T1.scale(lambda);
515 auto T2 = Hdummy.cdag_nc(x[j],x[i],y[j],y[i]); T2.scale(conjIfcomplex(lambda));
516 auto Term = sum(T1,T2);
517
518 if (VERB>0) lout << "corr_hopping: " << i << ", " << j << ", lambda=" << lambda << endl;
519
520 Terms.Hmpo_corr_hopping[s] = sum(Terms.Hmpo_corr_hopping[s],Term);
521 Hterms.H2_corr_hopping[s] = MODEL(Terms.Hmpo_corr_hopping[s],dummy_params);
522 s = (s+1)%Terms.Hmpo_corr_hopping.size();
523 }
524
525 s = 0;
526 for (int t=0; t<Raw.pair_hopping.size(); ++t)
527 {
528 int i = get<0>(Raw.pair_hopping[t]);
529 int j = get<1>(Raw.pair_hopping[t]);
530 typename MODEL::Scalar_ lambda = get<2>(Raw.pair_hopping[t]);
531
532 if (VERB>0)
533 {
534 lout << "pair hopping: " << i << ", " << j << ", λ=" << lambda << ", Q="
535 << Hdummy.cdagcdag(x[i],y[i]).Qtarget() << ", "
536 << Hdummy.cc(x[j],y[j]).Qtarget()
537 << endl;
538 }
539
540 //auto T1 = prod(Hdummy.cdagcdag(i),Hdummy.cc(j)); T1.scale(lambda);
541 //auto T2 = prod(Hdummy.cdagcdag(j),Hdummy.cc(i)); T2.scale(conjIfcomplex(lambda));
542 auto T1 = prod(Hdummy.cdagcdag(x[i],y[i]),Hdummy.cc(x[j],y[j])); T1.scale(lambda);
543 auto T2 = prod(Hdummy.cdagcdag(x[j],y[j]),Hdummy.cc(x[i],y[i])); T2.scale(conjIfcomplex(lambda));
544 auto Term = sum(T1,T2);
545
546 Terms.Hmpo_pair_hopping[s] = sum(Terms.Hmpo_pair_hopping[s],Term);
547 Hterms.H2_pair_hopping[s] = MODEL(Terms.Hmpo_pair_hopping[s],dummy_params);
548 s = (s+1)%Terms.Hmpo_pair_hopping.size();
549 }
550
551 Terms.Hmpo_corr_hopping3.resize(Raw.corr_hopping3.size());
552 Terms.Hmpo_nonlocal_spin.resize(Raw.nonlocal_spin.size());
553 Terms.Hmpo_doublon_decay.resize(Raw.doublon_decay.size());
554
555 Hterms.H3_corr_hopping3.resize(Raw.corr_hopping3.size());
556 Hterms.H3_nonlocal_spin.resize(Raw.nonlocal_spin.size());
557 Hterms.H3_doublon_decay.resize(Raw.doublon_decay.size());
558
559 s = 0;
560 for (int t=0; t<Raw.corr_hopping3.size(); ++t)
561 {
562 int i = get<0>(Raw.corr_hopping3[t]);
563 int j = get<1>(Raw.corr_hopping3[t]);
564 int k = get<2>(Raw.corr_hopping3[t]);
565 typename MODEL::Scalar_ lambda = get<3>(Raw.corr_hopping3[t]);
566
567 //auto T1 = prod(Hdummy.n(i), Hdummy.cdagc(j,k)); T1.scale(0.5*lambda);
568 //auto T2 = prod(Hdummy.cdagc(k,j), Hdummy.n(i)); T2.scale(0.5*conjIfcomplex(lambda));
569 auto T1 = prod(Hdummy.n(x[i],y[i]), Hdummy.cdagc(x[j],x[k],y[j],y[k])); T1.scale(0.5*lambda);
570 auto T2 = prod(Hdummy.cdagc(x[k],x[j],y[k],y[j]), Hdummy.n(x[i],y[i])); T2.scale(0.5*conjIfcomplex(lambda));
571 auto Term = sum(T1,T2);
572 Terms.Hmpo_corr_hopping3[s] = sum(Terms.Hmpo_corr_hopping3[s],Term);
573 Hterms.H3_corr_hopping3[s] = MODEL(Terms.Hmpo_corr_hopping3[s],dummy_params);
574 s = (s+1)%Terms.Hmpo_corr_hopping3.size();
575 }
576
577 s = 0;
578 for (int t=0; t<Raw.nonlocal_spin.size(); ++t)
579 {
580 int i = get<0>(Raw.nonlocal_spin[t]);
581 int j = get<1>(Raw.nonlocal_spin[t]);
582 int k = get<2>(Raw.nonlocal_spin[t]);
583 typename MODEL::Scalar_ lambda = get<3>(Raw.nonlocal_spin[t]);
584 {
585 double factor = sqrt(2)*sqrt(3)*sqrt(3);
586 //auto T1 = prod(Hdummy.S(i,0,factor), Hdummy.cdagc3(j,k)); T1.scale(lambda);
587 //auto T2 = prod(Hdummy.cdagc3(k,j), Hdummy.Sdag(i,0,factor)); T2.scale(conjIfcomplex(lambda));
588 auto T1 = prod(Hdummy.S(x[i],y[i],factor), Hdummy.cdagc3(x[j],x[k],y[j],y[k])); T1.scale(lambda);
589 auto T2 = prod(Hdummy.cdagc3(x[k],x[j],y[k],y[j]), Hdummy.Sdag(x[i],y[i],factor)); T2.scale(conjIfcomplex(lambda));
590 auto Term = sum(T1,T2);
591 Terms.Hmpo_nonlocal_spin[s] = sum(Terms.Hmpo_nonlocal_spin[s],Term);
592 Hterms.H3_nonlocal_spin[s] = MODEL(Terms.Hmpo_nonlocal_spin[s],dummy_params);
593 }
594 // This is compensated by the 0.5 factor in corr_hopping3:
595 /*{
596 auto T1 = prod(Hdummy.n(x[i],y[i]), Hdummy.cdagc(x[j],x[k],y[j],y[k]));
597 auto T2 = prod(Hdummy.cdagc(x[k],x[j],y[k],y[j]), Hdummy.n(x[i],y[i]));
598 auto Term = sum(T1,T2);
599 Term.scale(-0.5*lambda);
600 Terms.Hmpo_nonlocal_spin[s] = sum(Terms.Hmpo_nonlocal_spin[s],Term);
601 }*/
602 s = (s+1)%Terms.Hmpo_nonlocal_spin.size();
603 }
604
605 s = 0;
606 for (int t=0; t<Raw.doublon_decay.size(); ++t)
607 {
608 int i = get<0>(Raw.doublon_decay[t]);
609 int j = get<1>(Raw.doublon_decay[t]);
610 int k = get<2>(Raw.doublon_decay[t]);
611 typename MODEL::Scalar_ lambda = get<3>(Raw.doublon_decay[t]);
612
613 if (VERB>0)
614 {
615 lout << "doublon decay: cdagcdag_" << i << ", cc1_" << j << "_" << k << ", λ=" << lambda << ", Q="
616 << Hdummy.cdagcdag(x[i],y[i]).Qtarget() << ", "
617 << Hdummy.cc1(x[j],x[k],y[j],y[k]).Qtarget()
618 << endl;
619 }
620
621 //auto T1 = prod(Hdummy.cdagcdag(i), Hdummy.cc1(j,k)); T1.scale(lambda);
622 //auto T2 = prod(Hdummy.cdagcdag1(j,k), Hdummy.cc(i)); T2.scale(conjIfcomplex(lambda));
623 auto T1 = prod(Hdummy.cdagcdag(x[i],y[i]), Hdummy.cc1(x[j],x[k],y[j],y[k])); T1.scale(lambda);
624 auto T2 = prod(Hdummy.cdagcdag1(x[j],x[k],y[j],y[k]), Hdummy.cc(x[i],y[i])); T2.scale(conjIfcomplex(lambda));
625 auto Term = diff(T2,T1);
626 Terms.Hmpo_doublon_decay[s] = sum(Terms.Hmpo_doublon_decay[s],Term);
627 Hterms.H3_doublon_decay[s] = MODEL(Terms.Hmpo_doublon_decay[s],dummy_params);
628 s = (s+1)%Terms.Hmpo_doublon_decay.size();
629 }
630
631 Terms.Hmpo_foursite.resize(Raw.foursite.size());
632
633 Hterms.H4.resize(Raw.foursite.size());
634
635 s = 0;
636 for (int t=0; t<Raw.foursite.size(); ++t)
637 {
638 int i = get<0>(Raw.foursite[t]);
639 int j = get<1>(Raw.foursite[t]);
640 int k = get<2>(Raw.foursite[t]);
641 int l = get<3>(Raw.foursite[t]);
642 typename MODEL::Scalar_ lambda = get<4>(Raw.foursite[t]);
643
644 if (VERB>0)
645 {
646 lout << "foursite: " << i << ", " << j << ", " << k << ", " << l << ", λ=" << lambda << ", Q="
647 << Hdummy.cdagcdag1(x[i],x[j],y[i],y[j]).Qtarget() << ", "
648 << Hdummy.cc1(x[k],x[l],y[k],y[l]).Qtarget() << ", "
649 << Hdummy.cdagcdag1(x[l],x[k],y[l],y[k]).Qtarget() << ", "
650 << Hdummy.cc1(x[j],x[i],y[j],y[i]).Qtarget()
651 << endl;
652 }
653
654 //auto T1 = prod(Hdummy.cdagcdag1(i,j), Hdummy.cc1(k,l)); T1.scale(-lambda);
655 //auto T2 = prod(Hdummy.cdagcdag1(l,k), Hdummy.cc1(j,i)); T2.scale(-conjIfcomplex(lambda));
656 auto T1 = prod(Hdummy.cdagcdag1(x[i],x[j],y[i],y[j]), Hdummy.cc1(x[k],x[l],y[k],y[l])); T1.scale(-lambda);
657 auto T2 = prod(Hdummy.cdagcdag1(x[l],x[k],y[l],y[k]), Hdummy.cc1(x[j],x[i],y[j],y[i])); T2.scale(-conjIfcomplex(lambda));
658 auto Term = sum(T1,T2);
659 //Term.scale(-lambda);
660 Terms.Hmpo_foursite[s] = sum(Terms.Hmpo_foursite[s],Term);
661 Hterms.H4[s] = MODEL(Terms.Hmpo_foursite[s],dummy_params);
662 s = (s+1)%Terms.Hmpo_foursite.size();
663 }
664}
665
666template<typename MODEL>
667template<typename Dummy>
668typename std::enable_if<!Dummy::IS_SPIN_SU2(),void>::type HubbardKspace<MODEL>::
669compute_raw()
670{
671 vector<tuple<int,int,int,int> > terms_1site;
672 vector<tuple<int,int,int,int> > terms_2site;
673 vector<tuple<int,int,int,int> > terms_3site;
674 vector<tuple<int,int,int,int> > terms_4site;
675
676 // Carry out sum over local space, fill Umap and sort terms according to 1/2/3/4 site
677 for (int k=0; k<L; ++k)
678 for (int l=0; l<L; ++l)
679 for (int m=0; m<L; ++m)
680 for (int n=0; n<L; ++n)
681 {
682 std::set<int> s;
683 s.insert(k);
684 s.insert(l);
685 s.insert(m);
686 s.insert(n);
687
688 complex<double> Uelement_ = 0.;
689 for (int i=0; i<L; ++i)
690 {
691 Uelement_ += conj(UU(i,k))*UU(i,l)*conj(UU(i,m))*UU(i,n);
692 }
693 typename MODEL::Scalar_ Uelement;
694 if constexpr(is_same<typename MODEL::Scalar_,double>::value)
695 {
696 if (abs(Uelement_.imag()) > 1e-10)
697 {
698 lout << termcolor::red << "Warning: Non-zero imaginary part of transformed matrix element=" << Uelement_.imag() << termcolor::reset << endl;
699 }
700 Uelement = Uelement_.real();
701 }
702 else
703 {
704 Uelement = Uelement_;
705 }
706 if (abs(Uelement) > 1e-8)
707 {
708 auto Vklmn = make_tuple(k,l,m,n);
709 Umap[Vklmn] = U*Uelement;
710
711 if (s.size() == 1)
712 {
713 terms_1site.push_back(Vklmn);
714 }
715 else if (s.size() == 2)
716 {
717 terms_2site.push_back(Vklmn);
718 }
719 else if (s.size() == 3)
720 {
721 terms_3site.push_back(Vklmn);
722 }
723 else if (s.size() == 4)
724 {
725 terms_4site.push_back(Vklmn);
726 }
727 }
728 }
729
730 // Start filling Raw (not yet MPOs):
731 // Groups the 1/2/3/4-site contributions according to the various MPO terms
732
734 for (int t=0; t<terms_1site.size(); ++t)
735 {
736 size_t i = get<0>(terms_1site[t]);
737 Terms.HubbardU_kspace.push_back(make_tuple(i,real(Umap[make_tuple(i,i,i,i)])));
738 }
739 Hterms.HubbardU_kspace = Terms.HubbardU_kspace;
740
742 while (terms_2site.size() != 0)
743 {
744 for (int t=0; t<terms_2site.size(); ++t)
745 {
746 int k1 = get<0>(terms_2site[t]);
747 int k2 = get<1>(terms_2site[t]);
748 int k3 = get<2>(terms_2site[t]);
749 int k4 = get<3>(terms_2site[t]);
750
751 if (k1==k4 and k2==k3)
752 {
753 int i = k1;
754 int j = k2;
755
756 auto Vijji = make_tuple(i,j,j,i);
757 auto Vjiij = make_tuple(j,i,i,j);
758
759 Raw.spin_exchange.push_back(make_tuple(i,j,Umap[Vijji]));
760 //lout << "i=" << i << ", j=" << j << ", spin_exchange: Vijji=" << Umap[Vijji] << ", Vjiij=" << Umap[Vjiij] << endl;
761
762 terms_2site.erase(std::remove(terms_2site.begin(), terms_2site.end(), Vijji), terms_2site.end());
763 terms_2site.erase(std::remove(terms_2site.begin(), terms_2site.end(), Vjiij), terms_2site.end());
764 break;
765 }
766 else if (k1==k2 and k3==k4)
767 {
768 int i = k1;
769 int j = k3;
770
771 auto Viijj = make_tuple(i,i,j,j);
772
773 Raw.density_density.push_back(make_tuple(i,j,Umap[Viijj]));
774 //lout << "i=" << i << ", j=" << j << ", density_density: Viijj=" << Umap[Viijj] << endl;
775
776 terms_2site.erase(std::remove(terms_2site.begin(), terms_2site.end(), Viijj), terms_2site.end());
777 break;
778 }
779 else if (k1==k2 and k2==k3)
780 {
781 int i = k1;
782 int j = k4;
783
784 auto Viiij = make_tuple(i,i,i,j);
785 auto Viiji = make_tuple(i,i,j,i);
786
787 Raw.corr_hopping.push_back(make_tuple(i,j,Umap[Viiij]));
788 //lout << "i=" << i << ", j=" << j << ", corr_hopping: Viiij=" << Umap[Viiij] << ", Viiji" << Umap[Viiji] << endl;
789
790 terms_2site.erase(std::remove(terms_2site.begin(), terms_2site.end(), Viiij), terms_2site.end());
791 terms_2site.erase(std::remove(terms_2site.begin(), terms_2site.end(), Viiji), terms_2site.end());
792 break;
793 }
794 else if (k1==k3 and k3==k4)
795 {
796 int i = k1;
797 int j = k2;
798
799 auto Vijii = make_tuple(i,j,i,i);
800 auto Vjiii = make_tuple(j,i,i,i);
801
802 Raw.corr_hoppingB.push_back(make_tuple(i,j,Umap[Vijii]));
803 //lout << "i=" << i << ", j=" << j << ", corr_hoppingB: Vijii=" << Umap[Vijii] << ", Vjiii" << Umap[Vjiii] << endl;
804
805 terms_2site.erase(std::remove(terms_2site.begin(), terms_2site.end(), Vijii), terms_2site.end());
806 terms_2site.erase(std::remove(terms_2site.begin(), terms_2site.end(), Vjiii), terms_2site.end());
807 break;
808 }
809 else if (k1==k3 and k2==k4)
810 {
811 int i = k1;
812 int j = k2;
813
814 auto Vijij = make_tuple(i,j,i,j);
815 auto Vjiji = make_tuple(j,i,j,i);
816
817 Raw.pair_hopping.push_back(make_tuple(i,j,Umap[Vijij]));
818 //lout << "i=" << i << ", j=" << j << ", pair_hopping: Vijij=" << Umap[Vijij] << ", Vjiji=" << Umap[Vjiji] << endl;
819
820 terms_2site.erase(std::remove(terms_2site.begin(), terms_2site.end(), Vijij), terms_2site.end());
821 terms_2site.erase(std::remove(terms_2site.begin(), terms_2site.end(), Vjiji), terms_2site.end());
822 break;
823 }
824// else
825// {
826// lout << "2site other: k1=" << k1 << ", k2=" << k2 << ", k3=" << k3 << ", k4=" << k4 << endl;
827// }
828 }
829 }
830
832 while (terms_3site.size() != 0)
833 {
834 for (int t=0; t<terms_3site.size(); ++t)
835 {
836 int k1 = get<0>(terms_3site[t]);
837 int k2 = get<1>(terms_3site[t]);
838 int k3 = get<2>(terms_3site[t]);
839 int k4 = get<3>(terms_3site[t]);
840 //lout << "k1=" << k1 << ", k2=" << k2 << ", k3=" << k3 << ", k4=" << k4 << endl;
841
842 if (k1==k2)
843 {
844 int i = k1;
845 int j = k3;
846 int k = k4;
847
848 auto Viijk = make_tuple(i,i,j,k);
849 auto Viikj = make_tuple(i,i,k,j);
850
851 Raw.corr_hopping3.push_back(make_tuple(i,j,k,Umap[Viijk]));
852 //lout << "i=" << i << ", j=" << j << ", k=" << k << ", corr_hopping3: Viijk=" << Umap[Viijk] << ", Viikj=" << Umap[Viikj] << endl;
853
854 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Viijk), terms_3site.end());
855 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Viikj), terms_3site.end());
856 break;
857 }
858 else if (k3==k4)
859 {
860 int i = k1;
861 int j = k2;
862 int k = k3;
863
864 auto Vijkk = make_tuple(i,j,k,k);
865 auto Vjikk = make_tuple(j,i,k,k);
866
867 Raw.corr_hopping3B.push_back(make_tuple(i,j,k,Umap[Vijkk]));
868 //lout << "i=" << i << ", j=" << j << ", k=" << k << ", corr_hopping3B: Vjkii=" << Umap[Vjkii] << ", Vkjii=" << Umap[Vkjii] << endl;
869
870 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Vijkk), terms_3site.end());
871 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Vjikk), terms_3site.end());
872 break;
873 }
874 else if (k1==k4)
875 {
876 int i = k1;
877 int j = k2;
878 int k = k3;
879
880 auto Vijki = make_tuple(i,j,k,i);
881 auto Vjiik = make_tuple(j,i,i,k);
882
883 Raw.nonlocal_spin.push_back(make_tuple(i,j,k,Umap[Vijki]));
884 //lout << "i=" << i << ", j=" << j << ", k=" << k << ", nonlocal_spin: Vijki=" << Umap[Vijki] << ", Vjiik=" << Umap[Vjiik] << endl;
885
886 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Vijki), terms_3site.end());
887 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Vjiik), terms_3site.end());
888 break;
889 }
890 else if (k2==k3)
891 {
892 int i = k1;
893 int j = k2;
894 int k = k4;
895
896 auto Vijjk = make_tuple(i,j,j,k);
897 auto Vjikj = make_tuple(j,i,k,j);
898
899 Raw.nonlocal_spinB.push_back(make_tuple(i,j,k,Umap[Vijjk]));
900 //lout << "i=" << i << ", j=" << j << ", k=" << k << ", nonlocal_spinB: Vijjk=" << Umap[Vijjk] << ", Vjikj=" << Umap[Vjikj] << endl;
901
902 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Vijjk), terms_3site.end());
903 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Vjikj), terms_3site.end());
904 break;
905 }
906 else if (k1==k3)
907 {
908 int i = k1;
909 int j = k2;
910 int k = k4;
911
912 auto Vijik = make_tuple(i,j,i,k);
913 auto Vjiki = make_tuple(j,i,k,i);
914
915 Raw.doublon_decay.push_back(make_tuple(i,j,k,Umap[Vijik]));
916 //lout << "i=" << i << ", j=" << j << ", k=" << k << ", doublon_decay: Vijik=" << Umap[Vijik] << ", Vjiki=" << Umap[Vjiki] << endl;
917
918 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Vijik), terms_3site.end());
919 terms_3site.erase(std::remove(terms_3site.begin(), terms_3site.end(), Vjiki), terms_3site.end());
920 break;
921 }
922// else
923// {
924// lout << "3site other: k1=" << k1 << ", k2=" << k2 << ", k3=" << k3 << ", k4=" << k4 << endl;
925// }
926 }
927 }
928
930 while (terms_4site.size() != 0)
931 {
932 for (int t=0; t<terms_4site.size(); ++t)
933 {
934 int k1 = get<0>(terms_4site[t]);
935 int k2 = get<1>(terms_4site[t]);
936 int k3 = get<2>(terms_4site[t]);
937 int k4 = get<3>(terms_4site[t]);
938
939 int i = k1;
940 int j = k2;
941 int k = k3;
942 int l = k4;
943
944 auto Vijkl = make_tuple(i,j,k,l);
945 auto Vjilk = make_tuple(j,i,l,k);
946
947 Raw.foursite.push_back(make_tuple(i,j,k,l,Umap[Vijkl]));
948 lout << "i=" << i << ", j=" << j << ", k=" << k << ", l=" << l << ", 4s: " << "Vikjl=" << Umap[Vijkl] << ", Vjilk=" << Umap[Vjilk] << endl;
949
950 terms_4site.erase(std::remove(terms_4site.begin(), terms_4site.end(), Vijkl), terms_4site.end());
951 terms_4site.erase(std::remove(terms_4site.begin(), terms_4site.end(), Vjilk), terms_4site.end());
952 break;
953 }
954 }
955}
956
957template<typename MODEL>
958template<typename Dummy>
959typename std::enable_if<!Dummy::IS_SPIN_SU2(),void>::type HubbardKspace<MODEL>::
960compute_MPO()
961{
962 MODEL Hdummy(Lred,dummy_params,BC::OPEN,DMRG::VERBOSITY::SILENT);
963
965 for (int t=0; t<Raw.spin_exchange.size(); ++t)
966 {
967 int i = get<0>(Raw.spin_exchange[t]);
968 int j = get<1>(Raw.spin_exchange[t]);
969 typename MODEL::Scalar_ lambda = get<2>(Raw.spin_exchange[t]); // Vijji
970
971 OPERATOR hop1 = prod(Hdummy.template cdagc<UP,UP>(i,j),Hdummy.template cdagc<DN,DN>(j,i)); hop1.scale(lambda);
972 OPERATOR hop2 = prod(Hdummy.template cdagc<UP,UP>(j,i),Hdummy.template cdagc<DN,DN>(i,j)); hop2.scale(conj(lambda));
973 OPERATOR hop12 = sum(hop1,hop2);
974
975 Terms.Hmpo_spin_exchange.push_back(hop12);
976 Hterms.H2_spin_exchange.push_back(MODEL(hop12,dummy_params));
977 }
978 for (int t=0; t<Raw.density_density.size(); ++t)
979 {
980 int i = get<0>(Raw.density_density[t]);
981 int j = get<1>(Raw.density_density[t]);
982 double lambda = real(get<2>(Raw.density_density[t])); // Viijj
983
984 OPERATOR hop12 = prod(Hdummy.template n<UP>(i),Hdummy.template n<DN>(j)); hop12.scale(lambda);
985
986 Terms.Hmpo_spin_exchange.push_back(hop12);
987 Hterms.H2_spin_exchange.push_back(MODEL(hop12,dummy_params));
988 }
989 for (int t=0; t<Raw.corr_hopping.size(); ++t)
990 {
991 int i = get<0>(Raw.corr_hopping[t]);
992 int j = get<1>(Raw.corr_hopping[t]);
993 typename MODEL::Scalar_ lambda = get<2>(Raw.corr_hopping[t]); // Viiij
994
995 OPERATOR hop1 = prod(Hdummy.template n<UP>(i),Hdummy.template cdagc<DN,DN>(i,j)); hop1.scale(lambda);
996 OPERATOR hop2 = prod(Hdummy.template n<UP>(i),Hdummy.template cdagc<DN,DN>(j,i)); hop2.scale(conj(lambda));
997 OPERATOR hop12 = sum(hop1,hop2);
998
999 Terms.Hmpo_corr_hopping.push_back(hop12);
1000 Hterms.H2_corr_hopping.push_back(MODEL(hop12,dummy_params));
1001 }
1002 for (int t=0; t<Raw.corr_hoppingB.size(); ++t)
1003 {
1004 int i = get<0>(Raw.corr_hoppingB[t]);
1005 int j = get<1>(Raw.corr_hoppingB[t]);
1006 typename MODEL::Scalar_ lambda = get<2>(Raw.corr_hoppingB[t]); // Vijii
1007
1008 OPERATOR hop1 = prod(Hdummy.template cdagc<UP,UP>(i,j),Hdummy.template n<DN>(i)); hop1.scale(lambda);
1009 OPERATOR hop2 = prod(Hdummy.template cdagc<UP,UP>(j,i),Hdummy.template n<DN>(i)); hop2.scale(conj(lambda));
1010 OPERATOR hop12 = sum(hop1,hop2);
1011
1012 Terms.Hmpo_corr_hopping.push_back(hop12);
1013 Hterms.H2_corr_hopping.push_back(MODEL(hop12,dummy_params));
1014 }
1015 for (int t=0; t<Raw.pair_hopping.size(); ++t)
1016 {
1017 int i = get<0>(Raw.pair_hopping[t]);
1018 int j = get<1>(Raw.pair_hopping[t]);
1019 typename MODEL::Scalar_ lambda = get<2>(Raw.pair_hopping[t]); // Vijij
1020
1021 OPERATOR hop1 = prod(Hdummy.template cdagc<UP,UP>(i,j),Hdummy.template cdagc<DN,DN>(i,j)); hop1.scale(lambda);
1022 OPERATOR hop2 = prod(Hdummy.template cdagc<UP,UP>(j,i),Hdummy.template cdagc<DN,DN>(j,i)); hop2.scale(conj(lambda));
1023 OPERATOR hop12 = sum(hop1,hop2);
1024
1025 Terms.Hmpo_pair_hopping.push_back(hop12);
1026 Hterms.H2_pair_hopping.push_back(MODEL(hop12,dummy_params));
1027 }
1029 for (int t=0; t<Raw.corr_hopping3.size(); ++t)
1030 {
1031 int i = get<0>(Raw.corr_hopping3[t]);
1032 int j = get<1>(Raw.corr_hopping3[t]);
1033 int k = get<2>(Raw.corr_hopping3[t]);
1034 typename MODEL::Scalar_ lambda = get<3>(Raw.corr_hopping3[t]); // Viijk
1035
1036 OPERATOR hop1 = prod(Hdummy.template n<UP>(i),Hdummy.template cdagc<DN,DN>(j,k)); hop1.scale(lambda);
1037 OPERATOR hop2 = prod(Hdummy.template n<UP>(i),Hdummy.template cdagc<DN,DN>(k,j)); hop2.scale(conj(lambda));
1038 OPERATOR hop12 = sum(hop1,hop2);
1039
1040 Terms.Hmpo_corr_hopping3.push_back(hop12);
1041 Hterms.H3_corr_hopping3.push_back(MODEL(hop12,dummy_params));
1042 }
1043 for (int t=0; t<Raw.corr_hopping3B.size(); ++t)
1044 {
1045 int i = get<0>(Raw.corr_hopping3B[t]);
1046 int j = get<1>(Raw.corr_hopping3B[t]);
1047 int k = get<2>(Raw.corr_hopping3B[t]);
1048 typename MODEL::Scalar_ lambda = get<3>(Raw.corr_hopping3B[t]); // Vijkk
1049
1050 OPERATOR hop1 = prod(Hdummy.template cdagc<UP,UP>(i,j),Hdummy.template n<DN>(k)); hop1.scale(lambda);
1051 OPERATOR hop2 = prod(Hdummy.template cdagc<UP,UP>(j,i),Hdummy.template n<DN>(k)); hop2.scale(conj(lambda));
1052 OPERATOR hop12 = sum(hop1,hop2);
1053
1054 Terms.Hmpo_corr_hopping3.push_back(hop12);
1055 Hterms.H3_corr_hopping3.push_back(MODEL(hop12,dummy_params));
1056 }
1057 for (int t=0; t<Raw.nonlocal_spin.size(); ++t)
1058 {
1059 int i = get<0>(Raw.nonlocal_spin[t]);
1060 int j = get<1>(Raw.nonlocal_spin[t]);
1061 int k = get<2>(Raw.nonlocal_spin[t]);
1062 typename MODEL::Scalar_ lambda = get<3>(Raw.nonlocal_spin[t]); // Vijki
1063
1064 OPERATOR hop1 = prod(Hdummy.template cdagc<UP,UP>(i,j),Hdummy.template cdagc<DN,DN>(k,i)); hop1.scale(lambda);
1065 OPERATOR hop2 = prod(Hdummy.template cdagc<UP,UP>(j,i),Hdummy.template cdagc<DN,DN>(i,k)); hop2.scale(conj(lambda));
1066 OPERATOR hop12 = sum(hop1,hop2);
1067
1068 Terms.Hmpo_nonlocal_spin.push_back(hop12);
1069 Hterms.H3_nonlocal_spin.push_back(MODEL(hop12,dummy_params));
1070 }
1071 for (int t=0; t<Raw.nonlocal_spinB.size(); ++t)
1072 {
1073 int i = get<0>(Raw.nonlocal_spinB[t]);
1074 int j = get<1>(Raw.nonlocal_spinB[t]);
1075 int k = get<2>(Raw.nonlocal_spinB[t]);
1076 typename MODEL::Scalar_ lambda = get<3>(Raw.nonlocal_spinB[t]); // Vijjk
1077
1078 OPERATOR hop1 = prod(Hdummy.template cdagc<UP,UP>(i,j),Hdummy.template cdagc<DN,DN>(j,k)); hop1.scale(lambda);
1079 OPERATOR hop2 = prod(Hdummy.template cdagc<UP,UP>(j,i),Hdummy.template cdagc<DN,DN>(k,j)); hop2.scale(conj(lambda));
1080 OPERATOR hop12 = sum(hop1,hop2);
1081
1082 Terms.Hmpo_nonlocal_spin.push_back(hop12);
1083 Hterms.H3_nonlocal_spin.push_back(MODEL(hop12,dummy_params));
1084 }
1085 for (int t=0; t<Raw.doublon_decay.size(); ++t)
1086 {
1087 int i = get<0>(Raw.doublon_decay[t]);
1088 int j = get<1>(Raw.doublon_decay[t]);
1089 int k = get<2>(Raw.doublon_decay[t]);
1090 typename MODEL::Scalar_ lambda = get<3>(Raw.doublon_decay[t]); // Vijik
1091
1092 OPERATOR hop1 = prod(Hdummy.template cdagc<UP,UP>(i,j),Hdummy.template cdagc<DN,DN>(i,k)); hop1.scale(lambda);
1093 OPERATOR hop2 = prod(Hdummy.template cdagc<UP,UP>(j,i),Hdummy.template cdagc<DN,DN>(k,i)); hop2.scale(conj(lambda));
1094 OPERATOR hop12 = sum(hop1,hop2);
1095
1096 Terms.Hmpo_doublon_decay.push_back(hop12);
1097 Hterms.H3_doublon_decay.push_back(MODEL(hop12,dummy_params));
1098 }
1100 for (int t=0; t<Raw.foursite.size(); ++t)
1101 {
1102 int i = get<0>(Raw.foursite[t]);
1103 int j = get<1>(Raw.foursite[t]);
1104 int k = get<2>(Raw.foursite[t]);
1105 int l = get<3>(Raw.foursite[t]);
1106 typename MODEL::Scalar_ lambda = get<4>(Raw.foursite[t]); // Vijkl
1107
1108 OPERATOR hop1 = prod(Hdummy.template cdagc<UP,UP>(i,j),Hdummy.template cdagc<DN,DN>(k,l)); hop1.scale(lambda);
1109 OPERATOR hop2 = prod(Hdummy.template cdagc<UP,UP>(j,i),Hdummy.template cdagc<DN,DN>(l,k)); hop2.scale(conj(lambda));
1110 OPERATOR hop12 = sum(hop1,hop2);
1111
1112 Terms.Hmpo_foursite.push_back(hop12);
1113 Hterms.H4.push_back(MODEL(hop12,dummy_params));
1114 }
1115
1116// for (int k=0; k<L; ++k)
1117// for (int l=0; l<L; ++l)
1118// for (int m=0; m<L; ++m)
1119// for (int n=0; n<L; ++n)
1120// {
1121// std::set<int> s;
1122// s.insert(k);
1123// s.insert(l);
1124// s.insert(m);
1125// s.insert(n);
1126//
1127// typename MODEL::Scalar_ Uelement = 0.;
1128// for (int i=0; i<L; ++i)
1129// {
1130// Uelement += conj(UU(i,k))*UU(i,l)*conj(UU(i,m))*UU(i,n);
1131// }
1132// if (abs(Uelement) <= 1e-8) continue;
1133//
1134// auto Vklmn = make_tuple(k,l,m,n);
1135// Umap[Vklmn] = U*Uelement;
1136//
1137// bool LOCAL = false;
1138// bool CORRHOP = false;
1139// bool SPINEXC = false;
1140// bool PAIRHOP = false;
1141// bool NONLOCSPIN = false;
1142// bool CORRHOP3 = false;
1143// bool DOUBLON = false;
1144// bool FOURSITE = false;
1145//
1146// if (s.size() == 1)
1147// {
1148// //lout << termcolor::blue << "Viiii(local):\t" << k << "\t" << l << "\t" << m << "\t" << n << ":\t" << U*Uelement << ", s.size()=" << s.size() << termcolor::reset << endl;
1149// LOCAL = true;
1150// }
1151// else if (s.size() == 2)
1152// {
1153// if (l==m and m==n and l==n and k!=l)
1154// {
1155// //lout << termcolor::magenta << "Vjiii(corrhopA):\t" << k << "\t" << l << "\t" << m << "\t" << n << ":\t" << U*Uelement << ", s.size()=" << s.size() << termcolor::reset << endl;
1156// CORRHOP = true;
1157// }
1158// else if (k==m and m==n and k==n and k!=l)
1159// {
1160// //lout << termcolor::magenta << "Vijii(corrhopB):\t" << k << "\t" << l << "\t" << m << "\t" << n << ":\t" << U*Uelement << ", s.size()=" << s.size() << termcolor::reset << endl;
1161// CORRHOP = true;
1162// }
1163// else if (k==l and l==n and k==n and k!=m)
1164// {
1165// //lout << termcolor::magenta << "Viiji(corrhopC):\t" << k << "\t" << l << "\t" << m << "\t" << n << ":\t" << U*Uelement << ", s.size()=" << s.size() << termcolor::reset << endl;
1166// CORRHOP = true;
1167// }
1168// else if (k==l and l==m and k==m and k!=n)
1169// {
1170// //lout << termcolor::magenta << "Viiij(corrhopD):\t" << k << "\t" << l << "\t" << m << "\t" << n << ":\t" << U*Uelement << ", s.size()=" << s.size() << termcolor::reset << endl;
1171// CORRHOP = true;
1172// }
1173// else if (k==m and l==n and k!=l)
1174// {
1175// //lout << termcolor::red << "Viijj(pairhop):\t" << k << "\t" << l << "\t" << m << "\t" << n << ":\t" << U*Uelement << ", s.size()=" << s.size() << termcolor::reset << endl;
1176// PAIRHOP = true;
1177// }
1178// else if (k==n and l==m and k!=l)
1179// {
1180// //lout << termcolor::green << "Vijij(spinflip):\t" << k << "\t" << l << "\t" << m << "\t" << n << ":\t" << U*Uelement << ", s.size()=" << s.size() << termcolor::reset << endl;
1181// SPINEXC = true;
1182// }
1183// else if (k==l and m==n and k!=m)
1184// {
1185// //lout << termcolor::green << "Viijj(Ising):\t" << k << "\t" << l << "\t" << m << "\t" << n << ":\t" << U*Uelement << ", s.size()=" << s.size() << termcolor::reset << endl;
1186// SPINEXC = true;
1187// }
1188// }
1189// else if (s.size() == 3)
1190// {
1191// if (k==l or m==n)
1192// {
1193// CORRHOP3 = true;
1194// }
1195// else if (k==m or l==n)
1196// {
1197// DOUBLON = true;
1198// }
1199// else if (k==n or l==m)
1200// {
1201// NONLOCSPIN = true;
1202// }
1203// }
1204// else if (s.size() == 4)
1205// {
1206// FOURSITE = true;
1207// }
1208//
1209// if (s.size() == 1)
1210// {
1211// Terms.HubbardU_kspace.push_back(make_tuple(k,real(U*Uelement)));
1212// }
1213// else
1214// {
1215// OPERATOR hop12 = prod(Hdummy.template cdagc<UP,UP>(k,l),Hdummy.template cdagc<DN,DN>(m,n)); hop12.scale(U*Uelement);
1216//
1217// if (SPINEXC)
1218// {
1219// Terms.Hmpo_spin_exchange.push_back(hop12);
1220// Hterms.H2_spin_exchange.push_back(MODEL(hop12,dummy));
1221// }
1222// else if (CORRHOP)
1223// {
1224// Terms.Hmpo_corr_hopping.push_back(hop12);
1225// Hterms.H2_corr_hopping.push_back(MODEL(hop12,dummy));
1226// }
1227// else if (PAIRHOP)
1228// {
1229// Terms.Hmpo_pair_hopping.push_back(hop12);
1230// Hterms.H2_pair_hopping.push_back(MODEL(hop12,dummy));
1231// }
1232// else if (NONLOCSPIN)
1233// {
1234// Terms.Hmpo_nonlocal_spin.push_back(hop12);
1235// Hterms.H3_nonlocal_spin.push_back(MODEL(hop12,dummy));
1236// }
1237// else if (CORRHOP3)
1238// {
1239// Terms.Hmpo_corr_hopping3.push_back(hop12);
1240// Hterms.H3_corr_hopping3.push_back(MODEL(hop12,dummy));
1241// }
1242// else if (DOUBLON)
1243// {
1244// Terms.Hmpo_doublon_decay.push_back(hop12);
1245// Hterms.H3_doublon_decay.push_back(MODEL(hop12,dummy));
1246// }
1247// else if (FOURSITE)
1248// {
1249// Terms.Hmpo_foursite.push_back(hop12);
1250// Hterms.H4.push_back(MODEL(hop12,dummy));
1251// }
1252// }
1253// }
1254 Hterms.HubbardU_kspace = Terms.HubbardU_kspace;
1255}
1256
1257template<typename MODEL>
1258MODEL HubbardKspace<MODEL>::
1259sum_all() const
1260{
1261 auto res = Terms.Hmpo_spin_exchange[0];
1262
1263 for (int s=1; s<Terms.Hmpo_spin_exchange.size(); ++s) res = sum(res,Terms.Hmpo_spin_exchange[s]);
1264 for (int s=0; s<Terms.Hmpo_density_density.size(); ++s) res = sum(res,Terms.Hmpo_density_density[s]);
1265 for (int s=0; s<Terms.Hmpo_pair_hopping.size(); ++s) res = sum(res,Terms.Hmpo_pair_hopping[s]);
1266 for (int s=0; s<Terms.Hmpo_corr_hopping.size(); ++s) res = sum(res,Terms.Hmpo_corr_hopping[s]);
1267 for (int s=0; s<Terms.Hmpo_corr_hopping3.size(); ++s) res = sum(res,Terms.Hmpo_corr_hopping3[s]);
1268 for (int s=0; s<Terms.Hmpo_nonlocal_spin.size(); ++s) res = sum(res,Terms.Hmpo_nonlocal_spin[s]);
1269 for (int s=0; s<Terms.Hmpo_doublon_decay.size(); ++s) res = sum(res,Terms.Hmpo_doublon_decay[s]);
1270 for (int s=0; s<Terms.Hmpo_foursite.size(); ++s) res = sum(res,Terms.Hmpo_foursite[s]);
1271
1272 return MODEL(res,dummy_params);
1273}
1274
1275template<typename MODEL>
1276MODEL HubbardKspace<MODEL>::
1277sum_all (const ArrayXXcd &hoppingRealSpace) const
1278{
1279 ArrayXXcd hopping = (UU.adjoint() * hoppingRealSpace.matrix() * UU).array();
1280
1281 MODEL Hdummy(Lred,dummy_params,BC::OPEN,DMRG::VERBOSITY::SILENT);
1282 auto [i,Uval] = Terms.HubbardU_kspace[0];
1283 OPERATOR dtmp = Hdummy.d(0);
1284 dtmp.scale(Uval);
1285 auto res = dtmp;
1286
1287 for (int t=1; t<Terms.HubbardU_kspace.size(); ++t)
1288 {
1289 auto [i,Uval] = Terms.HubbardU_kspace[t];
1290 OPERATOR dtmp = Hdummy.d(t);
1291 dtmp.scale(Uval);
1292 res = sum(res,dtmp);
1293 }
1294
1295 for (int i=0; i<L; ++i)
1296 for (int j=0; j<L; ++j)
1297 {
1298 if (abs(hopping(i,j)) > 1e-10)
1299 {
1300 if (i==j)
1301 {
1302 OPERATOR ntmp = Hdummy.n(i);
1303 ntmp.scale(-hopping(i,i));
1304 res = sum(res,ntmp);
1305 //lout << "i=" << i << ", t0val=" << t0val << endl;
1306 }
1307 else
1308 {
1309 OPERATOR htmp;
1310 htmp = Hdummy.cdagc(x[i], x[j], y[i], y[j]); htmp.scale(-hopping(i,j));
1311 res = sum(res,htmp);
1312 // htmp = Hdummy.cdagc(x[j], x[i], y[j], y[i]); htmp.scale(conj(hopping(i,j)));
1313 // res = sum(res,htmp);
1314 //lout << "i=" << i << ", j=" << j << ", hopval=" << -hopping(i,j) << endl;
1315 }
1316 }
1317 }
1318
1319 for (int s=0; s<Terms.Hmpo_spin_exchange.size(); ++s) res = sum(res,Terms.Hmpo_spin_exchange[s]);
1320 for (int s=0; s<Terms.Hmpo_density_density.size(); ++s) res = sum(res,Terms.Hmpo_density_density[s]);
1321 for (int s=0; s<Terms.Hmpo_pair_hopping.size(); ++s) res = sum(res,Terms.Hmpo_pair_hopping[s]);
1322 for (int s=0; s<Terms.Hmpo_corr_hopping.size(); ++s) res = sum(res,Terms.Hmpo_corr_hopping[s]);
1323 for (int s=0; s<Terms.Hmpo_corr_hopping3.size(); ++s) res = sum(res,Terms.Hmpo_corr_hopping3[s]);
1324 for (int s=0; s<Terms.Hmpo_nonlocal_spin.size(); ++s) res = sum(res,Terms.Hmpo_nonlocal_spin[s]);
1325 for (int s=0; s<Terms.Hmpo_doublon_decay.size(); ++s) res = sum(res,Terms.Hmpo_doublon_decay[s]);
1326 for (int s=0; s<Terms.Hmpo_foursite.size(); ++s) res = sum(res,Terms.Hmpo_foursite[s]);
1327
1328 return MODEL(res,dummy_params);
1329}
1330
1331template<typename MODEL>
1332Mpo<typename MODEL::Symmetry,typename MODEL::Scalar_> HubbardKspace<MODEL>::
1333sum_all_mpo() const
1334{
1335 auto res = Terms.Hmpo_spin_exchange[0];
1336
1337 for (int s=1; s<Terms.Hmpo_spin_exchange.size(); ++s) res = sum(res,Terms.Hmpo_spin_exchange[s]);
1338 for (int s=0; s<Terms.Hmpo_density_density.size(); ++s) res = sum(res,Terms.Hmpo_density_density[s]);
1339 for (int s=0; s<Terms.Hmpo_pair_hopping.size(); ++s) res = sum(res,Terms.Hmpo_pair_hopping[s]);
1340 for (int s=0; s<Terms.Hmpo_corr_hopping.size(); ++s) res = sum(res,Terms.Hmpo_corr_hopping[s]);
1341
1342 for (int s=0; s<Terms.Hmpo_corr_hopping3.size(); ++s) res = sum(res,Terms.Hmpo_corr_hopping3[s]);
1343 for (int s=0; s<Terms.Hmpo_nonlocal_spin.size(); ++s) res = sum(res,Terms.Hmpo_nonlocal_spin[s]);
1344 for (int s=0; s<Terms.Hmpo_doublon_decay.size(); ++s) res = sum(res,Terms.Hmpo_doublon_decay[s]);
1345
1346 for (int s=0; s<Terms.Hmpo_foursite.size(); ++s) res = sum(res,Terms.Hmpo_foursite[s]);
1347
1348 return res;
1349}
1350
1351template<typename MODEL>
1352Mpo<typename MODEL::Symmetry,typename MODEL::Scalar_> HubbardKspace<MODEL>::
1353sum_2site_mpo() const
1354{
1355 auto res = Terms.Hmpo_spin_exchange[0];
1356
1357 for (int s=1; s<Terms.Hmpo_spin_exchange.size(); ++s) res = sum(res,Terms.Hmpo_spin_exchange[s]);
1358 for (int s=0; s<Terms.Hmpo_density_density.size(); ++s) res = sum(res,Terms.Hmpo_density_density[s]);
1359 for (int s=0; s<Terms.Hmpo_pair_hopping.size(); ++s) res = sum(res,Terms.Hmpo_pair_hopping[s]);
1360 for (int s=0; s<Terms.Hmpo_corr_hopping.size(); ++s) res = sum(res,Terms.Hmpo_corr_hopping[s]);
1361
1362 return res;
1363}
1364
1365template<typename MODEL>
1366MODEL HubbardKspace<MODEL>::
1367sum_2site() const
1368{
1369 auto res = Terms.Hmpo_spin_exchange[0];
1370
1371 for (int s=1; s<Terms.Hmpo_spin_exchange.size(); ++s) res = sum(res,Terms.Hmpo_spin_exchange[s]);
1372 for (int s=0; s<Terms.Hmpo_density_density.size(); ++s) res = sum(res,Terms.Hmpo_density_density[s]);
1373 for (int s=0; s<Terms.Hmpo_pair_hopping.size(); ++s) res = sum(res,Terms.Hmpo_pair_hopping[s]);
1374 for (int s=0; s<Terms.Hmpo_corr_hopping.size(); ++s) res = sum(res,Terms.Hmpo_corr_hopping[s]);
1375
1376 return MODEL(res,dummy_params);
1377}
1378
1379template<typename MODEL>
1380MODEL HubbardKspace<MODEL>::
1381sum_3site() const
1382{
1383 auto res = Terms.Hmpo_spin_exchange[0];
1384
1385 for (int s=0; s<Terms.Hmpo_corr_hopping3.size(); ++s) res = sum(res,Terms.Hmpo_corr_hopping3[s]);
1386 for (int s=0; s<Terms.Hmpo_nonlocal_spin.size(); ++s) res = sum(res,Terms.Hmpo_nonlocal_spin[s]);
1387 for (int s=0; s<Terms.Hmpo_doublon_decay.size(); ++s) res = sum(res,Terms.Hmpo_doublon_decay[s]);
1388
1389 return MODEL(res,dummy_params);
1390}
1391
1392template<typename MODEL>
1393MODEL HubbardKspace<MODEL>::
1394sum_4site() const
1395{
1396 auto res = Terms.Hmpo_spin_exchange[0];
1397
1398 for (int s=0; s<Terms.Hmpo_foursite.size(); ++s) res = sum(res,Terms.Hmpo_foursite[s]);
1399
1400 return MODEL(res,dummy_params);
1401}
1402
1403#endif
sum
Hamiltonian sum(const Hamiltonian &H1, const Hamiltonian &H2, DMRG::VERBOSITY::OPTION VERBOSITY=DMRG::VERBOSITY::SILENT)
Definition: DmrgLinearAlgebra.h:866
prod
Hamiltonian prod(const Hamiltonian &H1, const Hamiltonian &H2, const qarray< Hamiltonian::Symmetry::Nq > &Qtot=Hamiltonian::Symmetry::qvacuum(), DMRG::VERBOSITY::OPTION VERBOSITY=DMRG::VERBOSITY::SILENT)
Definition: DmrgLinearAlgebra.h:880
diff
Mpo< Symmetry, Scalar > diff(const Mpo< Symmetry, Scalar > &H1, const Mpo< Symmetry, Scalar > &H2, DMRG::VERBOSITY::OPTION VERBOSITY=DMRG::VERBOSITY::SILENT)
Definition: DmrgLinearAlgebra.h:914
conjIfcomplex
double conjIfcomplex(double x)
Definition: DmrgTypedefs.h:27
OPEN
@ OPEN
Definition: DmrgTypedefs.h:163
HubbardKspace::UU
MatrixXcd UU
Definition: HubbardKspace.h:170
HubbardKspace::sum_2site
MODEL sum_2site() const
Definition: HubbardKspace.h:1367
HubbardKspace::info
string info() const
Definition: HubbardKspace.h:187
HubbardKspace::HubbardKspace
HubbardKspace(const MatrixXcd &UU_input, double U_input, DMRG::VERBOSITY::OPTION VERB_input=DMRG::VERBOSITY::SILENT, bool VUMPS_input=false, vector< int > x_input={}, vector< int > y_input={})
Definition: HubbardKspace.h:79
HubbardKspace::sum_all_mpo
Mpo< typename MODEL::Symmetry, typename MODEL::Scalar_ > sum_all_mpo() const
Definition: HubbardKspace.h:1333
HubbardKspace::x
vector< int > x
Definition: HubbardKspace.h:177
HubbardKspace::get_Hterms
KspaceHTerms< MODEL > get_Hterms() const
Definition: HubbardKspace.h:152
HubbardKspace::sum_all
MODEL sum_all() const
Definition: HubbardKspace.h:1259
HubbardKspace::VERB
DMRG::VERBOSITY::OPTION VERB
Definition: HubbardKspace.h:166
HubbardKspace::dummy_params
vector< Param > dummy_params
Definition: HubbardKspace.h:164
HubbardKspace::HubbardKspace
HubbardKspace(const MatrixXcd &UU_input, double U_input, const vector< Param > &params, DMRG::VERBOSITY::OPTION VERB_input=DMRG::VERBOSITY::SILENT, bool VUMPS_input=false, vector< int > x_input={}, vector< int > y_input={})
Definition: HubbardKspace.h:110
HubbardKspace::compute_MPO
std::enable_if< Dummy::IS_SPIN_SU2(), void >::type compute_MPO()
Definition: HubbardKspace.h:470
HubbardKspace::Raw
KspaceRawTerms< typename MODEL::Scalar_ > Raw
Definition: HubbardKspace.h:180
HubbardKspace::y
vector< int > y
Definition: HubbardKspace.h:178
HubbardKspace::Umap
map< tuple< int, int, int, int >, typename MODEL::Scalar_ > Umap
Definition: HubbardKspace.h:172
HubbardKspace::sum_4site
MODEL sum_4site() const
Definition: HubbardKspace.h:1394
HubbardKspace::sum_3site
MODEL sum_3site() const
Definition: HubbardKspace.h:1381
HubbardKspace::compute_raw
std::enable_if< Dummy::IS_SPIN_SU2(), void >::type compute_raw()
Definition: HubbardKspace.h:205
HubbardKspace::OPERATOR
Mpo< typename MODEL::Symmetry, typename MODEL::Scalar_ > OPERATOR
Definition: HubbardKspace.h:73
HubbardKspace::Terms
KspaceMpoTerms< MODEL > Terms
Definition: HubbardKspace.h:181
HubbardKspace::sum_2site_mpo
Mpo< typename MODEL::Symmetry, typename MODEL::Scalar_ > sum_2site_mpo() const
Definition: HubbardKspace.h:1353
HubbardKspace::Hterms
KspaceHTerms< MODEL > Hterms
Definition: HubbardKspace.h:182
MpoTerms::scale
void scale(const Scalar factor, const Scalar offset=0., const std::size_t power=0ul, const double tolerance=1.e-14)
Definition: MpoTerms.h:2857
Mpo
Definition: Mpo.h:40
KspaceHTerms::H3_doublon_decay
vector< MODEL > H3_doublon_decay
Definition: HubbardKspace.h:62
KspaceHTerms::H2_corr_hopping
vector< MODEL > H2_corr_hopping
Definition: HubbardKspace.h:57
KspaceHTerms::H2_pair_hopping
vector< MODEL > H2_pair_hopping
Definition: HubbardKspace.h:56
KspaceHTerms::HubbardU_kspace
vector< tuple< size_t, double > > HubbardU_kspace
Definition: HubbardKspace.h:52
KspaceHTerms::H3_corr_hopping3
vector< MODEL > H3_corr_hopping3
Definition: HubbardKspace.h:60
KspaceHTerms::H4
vector< MODEL > H4
Definition: HubbardKspace.h:65
KspaceHTerms::H3_nonlocal_spin
vector< MODEL > H3_nonlocal_spin
Definition: HubbardKspace.h:61
KspaceHTerms::H2_spin_exchange
vector< MODEL > H2_spin_exchange
Definition: HubbardKspace.h:55
KspaceMpoTerms::Hmpo_foursite
vector< Mpo< typename MODEL::Symmetry, typename MODEL::Scalar_ > > Hmpo_foursite
Definition: HubbardKspace.h:45
KspaceMpoTerms::Hmpo_pair_hopping
vector< Mpo< typename MODEL::Symmetry, typename MODEL::Scalar_ > > Hmpo_pair_hopping
Definition: HubbardKspace.h:36
KspaceMpoTerms::Hmpo_spin_exchange
vector< Mpo< typename MODEL::Symmetry, typename MODEL::Scalar_ > > Hmpo_spin_exchange
Definition: HubbardKspace.h:34
KspaceMpoTerms::Hmpo_density_density
vector< Mpo< typename MODEL::Symmetry, typename MODEL::Scalar_ > > Hmpo_density_density
Definition: HubbardKspace.h:35
KspaceMpoTerms::HubbardU_kspace
vector< tuple< size_t, double > > HubbardU_kspace
Definition: HubbardKspace.h:31
KspaceMpoTerms::Hmpo_corr_hopping
vector< Mpo< typename MODEL::Symmetry, typename MODEL::Scalar_ > > Hmpo_corr_hopping
Definition: HubbardKspace.h:37
KspaceMpoTerms::Hmpo_corr_hopping3
vector< Mpo< typename MODEL::Symmetry, typename MODEL::Scalar_ > > Hmpo_corr_hopping3
Definition: HubbardKspace.h:40
KspaceMpoTerms::Hmpo_nonlocal_spin
vector< Mpo< typename MODEL::Symmetry, typename MODEL::Scalar_ > > Hmpo_nonlocal_spin
Definition: HubbardKspace.h:41
KspaceMpoTerms::Hmpo_doublon_decay
vector< Mpo< typename MODEL::Symmetry, typename MODEL::Scalar_ > > Hmpo_doublon_decay
Definition: HubbardKspace.h:42
KspaceRawTerms::Umap
map< tuple< int, int, int, int >, Scalar > Umap
Definition: HubbardKspace.h:24
KspaceRawTerms::doublon_decay
vector< tuple< int, int, int, Scalar > > doublon_decay
Definition: HubbardKspace.h:19
KspaceRawTerms::foursite
vector< tuple< int, int, int, int, Scalar > > foursite
Definition: HubbardKspace.h:22
KspaceRawTerms::nonlocal_spin
vector< tuple< int, int, int, Scalar > > nonlocal_spin
Definition: HubbardKspace.h:15
KspaceRawTerms::nonlocal_spinB
vector< tuple< int, int, int, Scalar > > nonlocal_spinB
Definition: HubbardKspace.h:16
KspaceRawTerms::spin_exchange
vector< tuple< int, int, Scalar > > spin_exchange
Definition: HubbardKspace.h:8
KspaceRawTerms::density_density
vector< tuple< int, int, Scalar > > density_density
Definition: HubbardKspace.h:9
KspaceRawTerms::corr_hopping3
vector< tuple< int, int, int, Scalar > > corr_hopping3
Definition: HubbardKspace.h:17
KspaceRawTerms::pair_hopping
vector< tuple< int, int, Scalar > > pair_hopping
Definition: HubbardKspace.h:12
KspaceRawTerms::corr_hopping
vector< tuple< int, int, Scalar > > corr_hopping
Definition: HubbardKspace.h:10
KspaceRawTerms::corr_hopping3B
vector< tuple< int, int, int, Scalar > > corr_hopping3B
Definition: HubbardKspace.h:18
KspaceRawTerms::corr_hoppingB
vector< tuple< int, int, Scalar > > corr_hoppingB
Definition: HubbardKspace.h:11