VMPS++
Loading...
Searching...
No Matches
KondoNecklace.h
Go to the documentation of this file.
1#ifndef KONDONECKLACE_H_
2#define KONDONECKLACE_H_
3
4//#define OPLABELS
5
6#include<map>
7#include<string>
8
9#include "KondoNecklaceObservables.h"
10#include "symmetry/U0.h"
11#include "bases/SpinBase.h"
12#include "Mpo.h"
13#include "ParamReturner.h"
14#include "Geometry2D.h" // from TOOLS
15
16namespace VMPS
17{
18
19class KondoNecklace : public Mpo<Sym::U0,double>, public KondoNecklaceObservables<Sym::U0>, public ParamReturner
20{
21public:
22 typedef Sym::U0 Symmetry;
23 MAKE_TYPEDEFS(KondoNecklace)
24 static constexpr MODEL_FAMILY FAMILY = HEISENBERG;
25
26private:
27 typedef typename Symmetry::qType qType;
28 typedef Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> MatrixType;
29 typedef Eigen::SparseMatrix<double,Eigen::ColMajor,EIGEN_DEFAULT_SPARSE_INDEX_TYPE> SparseMatrixType;
30 typedef SiteOperatorQ<Symmetry,MatrixType> OperatorType;
31
32public:
36 KondoNecklace() : Mpo<Symmetry>(), ParamReturner(KondoNecklace::sweep_defaults) {};
37
43 KondoNecklace(const std::size_t L, const std::vector<Param>& params={}, const BC boundary=BC::OPEN, const DMRG::VERBOSITY::OPTION VERB=DMRG::VERBOSITY::ON_EXIT);
44
45
53 static void set_operators(const std::vector<SpinBase<Symmetry>>& Bsub, const std::vector<SpinBase<Symmetry>>& Bimp, const ParamHandler& P,
54 PushType<SiteOperator<Symmetry,double>,double>& pushlist, std::vector<std::vector<std::string>>& labellist, const BC boundary);
55
59 static const std::map<string,std::any> defaults;
60
64 static const std::map<string,std::any> sweep_defaults;
65
69 bool validate(qType qnum);
70};
71
72const std::map<string,std::any> KondoNecklace::defaults =
73{
74 {"Jlocxy",1.}, {"Jlocz",1.}, {"Jparaxy",1.}, {"Jparaz",1.}, {"Jperpxy",0.}, {"Jperpz",0.}, {"Jprimexy",1.}, {"Jprimez",1.},
75 {"Dimp",2ul}, {"Dsub",2ul}, {"Ly",1ul},
76 {"maxPower",2ul}, {"CYLINDER",false}
77};
78
79const std::map<string,std::any> KondoNecklace::sweep_defaults =
80{
81 {"max_alpha",100.}, {"min_alpha",1e-11}, {"lim_alpha",12ul}, {"eps_svd",1e-7},
82 {"Dincr_abs", 4ul}, {"Dincr_per", 2ul}, {"Dincr_rel", 1.1},
83 {"min_Nsv",1ul}, {"max_Nrich",-1},
84 {"max_halfsweeps",100ul}, {"min_halfsweeps",24ul},
85 {"Dinit",5ul}, {"Qinit",6ul}, {"Dlimit",250ul},
86 {"tol_eigval",1e-9}, {"tol_state",1e-8},
87 {"savePeriod",0ul}, {"CALC_S_ON_EXIT", true}, {"CONVTEST", DMRG::CONVTEST::VAR_HSQ}
88};
89
90KondoNecklace::KondoNecklace(const std::size_t L, const std::vector<Param>& params, const BC boundary, const DMRG::VERBOSITY::OPTION VERB)
91: Mpo<Symmetry>(L, Symmetry::qvacuum(), "KondoNecklace", PROP::HERMITIAN, PROP::NON_UNITARY, boundary, VERB),
92 KondoNecklaceObservables<Symmetry>(L,params,KondoNecklace::defaults),
93 ParamReturner(KondoNecklace::sweep_defaults)
94{
95 ParamHandler P(params,defaults);
96 this->set_verbosity(VERB);
97 std::size_t Lcell = P.size();
98 Bsub.resize(N_sites);
99 Bimp.resize(N_sites);
100 for (size_t loc=0; loc<N_sites; ++loc)
101 {
102 N_phys += P.get<size_t>("Ly",loc%Lcell);
103 setLocBasis((Bsub[loc].get_basis().combine(Bimp[loc].get_basis())).qloc(),loc);
104 }
105
106 PushType<SiteOperator<Symmetry,double>,double> pushlist;
107 std::vector<std::vector<std::string>> labellist(N_sites);
108 set_operators(Bsub, Bimp, P, pushlist, labellist, boundary);
109
110 this->construct_from_pushlist(pushlist, labellist, Lcell);
111 this->finalize(PROP::COMPRESS, P.get<std::size_t>("maxPower"));
112
113 this->precalc_TwoSiteData();
114}
115
116void KondoNecklace::set_operators(const std::vector<SpinBase<Symmetry>>& Bsub, const std::vector<SpinBase<Symmetry>>& Bimp, const ParamHandler &P,
117 PushType<SiteOperator<Symmetry,double>,double>& pushlist, std::vector<std::vector<std::string>>& labellist, const BC boundary)
118{
119 std::size_t Lcell = P.size();
120 std::size_t N_sites = Bsub.size();
121 for(std::size_t loc=0; loc<N_sites; ++loc)
122 {
123 std::size_t orbitals = Bsub[loc].orbitals();
124 std::size_t next_orbitals = Bsub[(loc+1)%N_sites].orbitals();
125 std::size_t nextn_orbitals = Bsub[(loc+2)%N_sites].orbitals();
126
127
128
129 std::stringstream ss1, ss2, ss3;
130 ss1 << "S_sub=" << print_frac_nice(frac(P.get<size_t>("Dsub",loc%Lcell)-1,2));
131 ss2 << "S_imp=" << print_frac_nice(frac(P.get<size_t>("Dimp",loc%Lcell)-1,2));
132 ss3 << "Ly=" << P.get<size_t>("Ly",loc%Lcell);
133 labellist[loc].push_back(ss1.str());
134 labellist[loc].push_back(ss2.str());
135 labellist[loc].push_back(ss3.str());
136
137
138 // Local Terms: J_Kondo
139 param1d Jlocxy = P.fill_array1d<double>("Jlocxy", "Jlocxy_array", orbitals, loc%Lcell);
140 param1d Jlocz = P.fill_array1d<double>("Jlocz", "Jlocz_array", orbitals, loc%Lcell);
141 if(Jlocxy.x != 0.)
142 {
143 labellist[loc].push_back(Jlocxy.label);
144 labellist[loc].push_back(Jlocz.label);
145 for(int alpha=0; alpha<orbitals; ++alpha)
146 {
147 std::vector<OperatorType> ops(1);
148 ops[0] = 0.5*Jlocxy(alpha) * kroneckerProduct(Bsub[loc].Scomp(SP,alpha), Bimp[loc].Scomp(SM,alpha));
149 ops[0] += 0.5*Jlocxy(alpha) * kroneckerProduct(Bsub[loc].Scomp(SM,alpha), Bimp[loc].Scomp(SP,alpha));
150 ops[0] += Jlocz(alpha) * kroneckerProduct(Bsub[loc].Scomp(SZ,alpha), Bimp[loc].Scomp(SZ,alpha));
151 pushlist.push_back(std::make_tuple(loc, ops, 1.));
152 }
153 }
154
155 auto push_full = [&N_sites, &loc, &Bimp, &Bsub, &P, &pushlist, &labellist, &boundary] (string xxxFull, string label,
156 const vector<SiteOperatorQ<Symmetry,Eigen::MatrixXd> > &first,
157 const vector<vector<SiteOperatorQ<Symmetry,Eigen::MatrixXd> > > &last,
158 vector<double> factor) -> void
159 {
160 ArrayXXd Full = P.get<Eigen::ArrayXXd>(xxxFull);
161 vector<vector<std::pair<size_t,double> > > R = Geometry2D::rangeFormat(Full);
162
163 if (static_cast<bool>(boundary)) {assert(R.size() == N_sites and "Use an (N_sites)x(N_sites) hopping matrix for open BC!");}
164 else {assert(R.size() >= 2*N_sites and "Use at least a (2*N_sites)x(N_sites) hopping matrix for infinite BC!");}
165
166 for (size_t j=0; j<first.size(); j++)
167 for (size_t h=0; h<R[loc].size(); ++h)
168 {
169 size_t range = R[loc][h].first;
170 double value = R[loc][h].second;
171
172 if (range != 0)
173 {
174 vector<SiteOperatorQ<Symmetry,Eigen::MatrixXd> > ops(range+1);
175 ops[0] = first[j];
176 for (size_t i=1; i<range; ++i)
177 {
178 ops[i] = kroneckerProduct(Bsub[(loc+i)%N_sites].Id(), Bimp[(loc+i)%N_sites].Id());
179 }
180 ops[range] = last[j][(loc+range)%N_sites];
181 pushlist.push_back(std::make_tuple(loc, ops, factor[j] * value));
182 }
183 }
184
185 stringstream ss;
186 ss << label << "(" << Geometry2D::hoppingInfo(Full) << ")";
187 labellist[loc].push_back(ss.str());
188 };
189
190 // Case where a full coupling matrix is providedf: Jᵢⱼ (all the code below this funtion will be skipped then.)
191 if (P.HAS("Jxyfull"))
192 {
193 vector<SiteOperatorQ<Symmetry,Eigen::MatrixXd> > first {kroneckerProduct(Bsub[loc].Scomp(SP,0),Bimp[loc].Id()),
194 kroneckerProduct(Bsub[loc].Scomp(SM,0),Bimp[loc].Id())};
195 vector<SiteOperatorQ<Symmetry,Eigen::MatrixXd> > Sm_ranges(N_sites); for (size_t i=0; i<N_sites; i++) {Sm_ranges[i] = kroneckerProduct(Bsub[i].Scomp(SM,0),Bimp[i].Id());}
196 vector<SiteOperatorQ<Symmetry,Eigen::MatrixXd> > Sp_ranges(N_sites); for (size_t i=0; i<N_sites; i++) {Sp_ranges[i] = kroneckerProduct(Bsub[i].Scomp(SP,0),Bimp[i].Id());}
197 vector<vector<SiteOperatorQ<Symmetry,Eigen::MatrixXd> > > last {Sm_ranges,Sp_ranges};
198 push_full("Jxyfull", "Jxyᵢⱼ", first, last, {0.5,0.5});
199 }
200 if (P.HAS("Jzfull"))
201 {
202 vector<SiteOperatorQ<Symmetry,Eigen::MatrixXd> > first {kroneckerProduct(Bsub[loc].Scomp(SZ,0),Bimp[loc].Id())};
203 vector<SiteOperatorQ<Symmetry,Eigen::MatrixXd> > Sz_ranges(N_sites); for (size_t i=0; i<N_sites; i++) {Sz_ranges[i] = kroneckerProduct(Bsub[i].Scomp(SZ,0),Bimp[i].Id());}
204 vector<vector<SiteOperatorQ<Symmetry,Eigen::MatrixXd> > > last {Sz_ranges};
205 push_full("Jzfull", "Jzᵢⱼ", first, last, {1.});
206 }
207 if (P.HAS("Ixyfull"))
208 {
209 vector<SiteOperatorQ<Symmetry,Eigen::MatrixXd> > first {kroneckerProduct(Bsub[loc].Id(),Bimp[loc].Scomp(SP,0)),
210 kroneckerProduct(Bsub[loc].Id(),Bimp[loc].Scomp(SM,0))};
211 vector<SiteOperatorQ<Symmetry,Eigen::MatrixXd> > Sm_ranges(N_sites); for (size_t i=0; i<N_sites; i++) {Sm_ranges[i] = kroneckerProduct(Bsub[i].Id(),Bimp[i].Scomp(SM,0));}
212 vector<SiteOperatorQ<Symmetry,Eigen::MatrixXd> > Sp_ranges(N_sites); for (size_t i=0; i<N_sites; i++) {Sp_ranges[i] = kroneckerProduct(Bsub[i].Id(),Bimp[i].Scomp(SP,0));}
213 vector<vector<SiteOperatorQ<Symmetry,Eigen::MatrixXd> > > last {Sm_ranges,Sp_ranges};
214 push_full("Ixyfull", "Ixyᵢⱼ", first, last, {0.5,0.5});
215 }
216 if (P.HAS("Izfull"))
217 {
218 vector<SiteOperatorQ<Symmetry,Eigen::MatrixXd> > first {kroneckerProduct(Bsub[loc].Id(),Bimp[loc].Scomp(SZ,0))};
219 vector<SiteOperatorQ<Symmetry,Eigen::MatrixXd> > Sz_ranges(N_sites); for (size_t i=0; i<N_sites; i++) {Sz_ranges[i] = kroneckerProduct(Bsub[i].Id(),Bimp[i].Scomp(SZ,0));}
220 vector<vector<SiteOperatorQ<Symmetry,Eigen::MatrixXd> > > last {Sz_ranges};
221 push_full("Izfull", "Izᵢⱼ", first, last, {1.});
222 }
223 if (P.HAS("Jzfull") or P.HAS("Jxyfull") or P.HAS("Izfull") or P.HAS("Ixyfull")) {continue;}
224
225
226 // Nearest-neighbour terms: J
227
228 // param2d Jpara = P.fill_array2d<double>("Jpara", "Jpara_array", {orbitals, next_orbitals}, loc%Lcell);
229 // if((Jpara.a != 0.).any())
230 // {
231 // labellist[loc].push_back(Jpara.label);
232 // if(loc < N_sites-1 or boundary == BC::INFINITE)
233 // {
234 // for(int alpha=0; alpha<orbitals; ++alpha)
235 // {
236 // for(int beta=0; beta<next_orbitals; ++beta)
237 // {
238 // double lambda = std::sqrt(3.)*Jpara.a(alpha,beta);
239 // std::vector<OperatorType> ops(2);
240 // ops[0] = OperatorType::outerprod(Bsub[loc].Sdag(alpha), Bimp[loc].Id(), {3});
241 // ops[1] = OperatorType::outerprod(Bsub[(loc+1)%N_sites].S(beta), Bimp[(loc+1)%N_sites].Id(), {3});
242 // pushlist.push_back(std::make_tuple(loc, ops, lambda));
243 // }
244 // }
245 // }
246 // }
247
248
249 // // Next-nearest-neighbour terms: J
250
251 // param2d Jprime = P.fill_array2d<double>("Jprime", "Jprime_array", {orbitals, nextn_orbitals}, loc%Lcell);
252 // if((Jprime.a != 0.).any())
253 // {
254 // labellist[loc].push_back(Jprime.label);
255 // if((N_sites > 1 and loc < N_sites-2) or boundary == BC::INFINITE)
256 // {
257 // for(int alpha=0; alpha<orbitals; ++alpha)
258 // {
259 // for(int beta=0; beta<nextn_orbitals; ++beta)
260 // {
261 // double lambda = std::sqrt(3.)*Jprime.a(alpha,beta);
262 // std::vector<OperatorType> ops(3);
263 // ops[0] = OperatorType::outerprod(Bsub[loc].Sdag(alpha), Bimp[loc].Id(), {3});
264 // ops[1] = OperatorType::outerprod(Bsub[(loc+1)%N_sites].Id(), Bimp[(loc+1)%N_sites].Id(), {1});
265 // ops[2] = OperatorType::outerprod(Bsub[(loc+2)%N_sites].S(beta), Bimp[(loc+2)%N_sites].Id(), {3});
266 // pushlist.push_back(std::make_tuple(loc, ops, lambda));
267 // }
268 // }
269 // }
270 // }
271 // }
272
273
274 // if(false) //(boundary == BC::PERIODIC)
275 // {
276 // std::size_t last_orbitals = Bsub[N_sites-1].orbitals();
277 // std::size_t first_orbitals = Bsub[0].orbitals();
278 // std::size_t previous_orbitals = Bsub[(2*N_sites-2)%N_sites].orbitals();
279 // std::size_t next_orbitals = Bsub[1%N_sites].orbitals();
280
281 // if(N_sites == 1)
282 // {
283 // param2d Jpara = P.fill_array2d<double>("Jpara", "Jpara_array", {first_orbitals,first_orbitals}, 0);
284 // if((Jpara.a != 0.).any())
285 // {
286 // labellist[0].push_back(Jpara.label);
287 // for(int alpha=0; alpha<last_orbitals; ++alpha)
288 // {
289 // for(int beta=0; beta<first_orbitals; ++beta)
290 // {
291 // double lambda = std::sqrt(3.)*Jpara.a(alpha,beta);
292 // std::vector<OperatorType> ops(1);
293 // ops[0] = OperatorType::outerprod(OperatorType::prod(Bsub[0].Sdag(alpha),Bsub[0].S(beta),{1}),Bimp[0].Id(),{1});
294 // pushlist.push_back(std::make_tuple(0ul, ops, lambda));
295 // }
296 // }
297 // }
298
299 // param2d Jprime = P.fill_array2d<double>("Jprime", "Jprime_array", {first_orbitals, first_orbitals}, 0%Lcell);
300 // if((Jprime.a != 0.).any())
301 // {
302 // labellist[0].push_back(Jprime.label);
303 // for(int alpha=0; alpha<first_orbitals; ++alpha)
304 // {
305 // for(int beta=0; beta<first_orbitals; ++beta)
306 // {
307 // double lambda = std::sqrt(3.)*Jprime.a(alpha,beta);
308 // std::vector<OperatorType> ops(1);
309 // ops[0] = OperatorType::outerprod(OperatorType::prod(Bsub[0].Sdag(alpha),Bsub[0].S(beta),{1}),Bimp[0].Id(),{1});
310 // pushlist.push_back(std::make_tuple(0ul, ops, lambda));
311 // }
312 // }
313 // }
314 // }
315 // else if(N_sites == 2)
316 // {
317 // param2d Jpara = P.fill_array2d<double>("Jpara", "Jpara_array", {last_orbitals,first_orbitals}, 1%Lcell);
318 // if((Jpara.a != 0.).any())
319 // {
320 // labellist[1].push_back(Jpara.label);
321 // for(int alpha=0; alpha<last_orbitals; ++alpha)
322 // {
323 // for(int beta=0; beta<first_orbitals; ++beta)
324 // {
325 // double lambda = std::sqrt(3.)*Jpara.a(alpha,beta);
326 // std::vector<OperatorType> ops(2);
327 // ops[0] = OperatorType::outerprod(Bsub[0].S(beta), Bimp[0].Id(), {3});
328 // ops[1] = OperatorType::outerprod(Bsub[1].Sdag(alpha), Bimp[1].Id(), {3});
329 // pushlist.push_back(std::make_tuple(0ul, ops, lambda));
330 // }
331 // }
332 // }
333
334 // param2d Jprime_prev_to_first = P.fill_array2d<double>("Jprime", "Jprime_array", {first_orbitals, first_orbitals}, 0%Lcell);
335 // if((Jprime_prev_to_first.a != 0.).any())
336 // {
337 // labellist[0].push_back(Jprime_prev_to_first.label);
338 // for(int alpha=0; alpha<first_orbitals; ++alpha)
339 // {
340 // for(int beta=0; beta<first_orbitals; ++beta)
341 // {
342 // double lambda = std::sqrt(3.)*Jprime_prev_to_first.a(alpha,beta);
343 // std::vector<OperatorType> ops(1);
344 // ops[0] = OperatorType::outerprod(OperatorType::prod(Bsub[0].Sdag(alpha),Bsub[0].S(beta),{1}),Bimp[0].Id(),{1});
345 // pushlist.push_back(std::make_tuple(0ul, ops, lambda));
346 // }
347 // }
348 // }
349
350 // param2d Jprime_last_to_next = P.fill_array2d<double>("Jprime", "Jprime_array", {last_orbitals, last_orbitals}, 1%Lcell);
351 // if((Jprime_last_to_next.a != 0.).any())
352 // {
353 // labellist[1].push_back(Jprime_last_to_next.label);
354 // for(int alpha=0; alpha<last_orbitals; ++alpha)
355 // {
356 // for(int beta=0; beta<last_orbitals; ++beta)
357 // {
358 // double lambda = std::sqrt(3.)*Jprime_last_to_next.a(alpha,beta);
359 // std::vector<OperatorType> ops(1);
360 // ops[0] = OperatorType::outerprod(OperatorType::prod(Bsub[1].Sdag(alpha),Bsub[1].S(beta),{1}),Bimp[1].Id(),{1});
361 // pushlist.push_back(std::make_tuple(1ul, ops, lambda));
362 // }
363 // }
364 // }
365 // }
366 // else
367 // {
368 // param2d Jpara = P.fill_array2d<double>("Jpara", "Jpara_array", {last_orbitals,first_orbitals}, (N_sites-1)%Lcell);
369 // if((Jpara.a != 0.).any())
370 // {
371 // labellist[N_sites-1].push_back(Jpara.label);
372 // for(int alpha=0; alpha<last_orbitals; ++alpha)
373 // {
374 // for(int beta=0; beta<first_orbitals; ++beta)
375 // {
376 // double lambda = std::sqrt(3.)*Jpara.a(alpha,beta);
377 // std::vector<OperatorType> ops(N_sites);
378 // ops[0] = OperatorType::outerprod(Bsub[0].S(beta), Bimp[0].Id(), {3});
379 // for(std::size_t loc=1; loc<N_sites-1; ++loc)
380 // {
381 // ops[loc] = OperatorType::outerprod(Bsub[loc].Id(), Bimp[loc].Id(), {1});
382 // }
383 // ops[N_sites-1] = OperatorType::outerprod(Bsub[N_sites-1].Sdag(alpha), Bimp[N_sites-1].Id(), {3});
384 // pushlist.push_back(std::make_tuple(0ul, ops, lambda));
385 // }
386 // }
387 // }
388
389 // param2d Jprime_prev_to_first = P.fill_array2d<double>("Jprime", "Jprime_array", {previous_orbitals, first_orbitals}, (N_sites-2)%Lcell);
390 // if((Jprime_prev_to_first.a != 0.).any())
391 // {
392 // labellist[N_sites-2].push_back(Jprime_prev_to_first.label);
393 // for(int alpha=0; alpha<previous_orbitals; ++alpha)
394 // {
395 // for(int beta=0; beta<first_orbitals; ++beta)
396 // {
397 // double lambda = std::sqrt(3.)*Jprime_prev_to_first.a(alpha,beta);
398 // std::vector<OperatorType> ops(N_sites-1);
399 // ops[0] = OperatorType::outerprod(Bsub[0].S(beta), Bimp[0].Id(), {3});
400 // for(std::size_t loc=1; loc<N_sites-2; ++loc)
401 // {
402 // ops[loc] = OperatorType::outerprod(Bsub[loc].Id(), Bimp[loc].Id(), {1});
403 // }
404 // ops[N_sites-2] = OperatorType::outerprod(Bsub[N_sites-2].Sdag(alpha), Bimp[N_sites-2].Id(), {3});
405 // pushlist.push_back(std::make_tuple(0ul, ops, lambda));
406 // }
407 // }
408 // }
409
410 // param2d Jprime_last_to_next = P.fill_array2d<double>("Jprime", "Jprime_array", {last_orbitals, next_orbitals}, (N_sites-1)%Lcell);
411 // if((Jprime_last_to_next.a != 0.).any())
412 // {
413 // labellist[N_sites-1].push_back(Jprime_last_to_next.label);
414 // for(int alpha=0; alpha<last_orbitals; ++alpha)
415 // {
416 // for(int beta=0; beta<next_orbitals; ++beta)
417 // {
418 // double lambda = std::sqrt(3.)*Jprime_last_to_next.a(alpha,beta);
419 // std::vector<OperatorType> ops(N_sites-1);
420 // ops[0] = OperatorType::outerprod(Bsub[1].S(beta), Bimp[1].Id(), {3});
421 // for(std::size_t loc=2; loc<N_sites-1; ++loc)
422 // {
423 // ops[loc-1] = OperatorType::outerprod(Bsub[loc].Id(), Bimp[loc].Id(), {1});
424 // }
425 // ops[N_sites-2] = OperatorType::outerprod(Bsub[N_sites-1].Sdag(alpha), Bimp[N_sites-1].Id(), {3});
426 // pushlist.push_back(std::make_tuple(1ul, ops, lambda));
427 // }
428 // }
429 // }
430 // }
431 }
432}
433
434// bool KondoNecklace::validate(qType qnum)
435// {
436// auto add = [](std::set<std::size_t>& left, std::set<std::size_t>& right) -> void
437// {
438// if(left.size() == 0)
439// {
440// left = right;
441// }
442// else
443// {
444// std::set<std::size_t> temp;
445// for(auto l : left) for(auto r : right)
446// {
447// std::size_t min = std::abs((static_cast<int>(l))-(static_cast<int>(r)))+1;
448// std::size_t max = l+r-1;
449// for(std::size_t i=min; i<=max; i+=2ul) temp.insert(i);
450// }
451// left = temp;
452// }
453// };
454// std::vector<std::set<std::size_t>> local(N_sites);
455// std::vector<std::set<std::size_t>> reachable(N_sites);
456// for(std::size_t loc=0; loc<N_sites; ++loc)
457// {
458// std::set<std::size_t> subspins;
459// if(Bsub[loc].orbitals()%2 == 0) subspins.insert(1ul);
460// for(std::size_t i = Bsub[loc].get_D(); i<=Bsub[loc].orbitals()*(Bsub[loc].get_D()-1)+1; i+=2ul) subspins.insert(i);
461// add(local[loc], subspins);
462// std::set<std::size_t> impspins;
463// if(Bimp[loc].orbitals()%2 == 0) impspins.insert(1ul);
464// for(std::size_t i = Bimp[loc].get_D(); i<=Bimp[loc].orbitals()*(Bimp[loc].get_D()-1)+1; i+=2ul) impspins.insert(i);
465// add(local[loc], impspins);
466// }
467// reachable[0] = local[0];
468// for(std::size_t loc=1; loc<N_sites; ++loc)
469// {
470// reachable[loc] = local[loc];
471// add(reachable[loc], reachable[loc-1]);
472// }
473
474// auto it = find(reachable[N_sites-1].begin(), reachable[N_sites-1].end(), qnum[0]);
475// return it!=reachable[N_sites-1].end();
476// }
477
478// Mpo<Sym::SU2<Sym::SpinSU2>> KondoNecklace::
479// Stot()
480// {
481// Mpo<Symmetry> Mout(this->N_sites, {3}, "S_tot", false, false, BC::OPEN, DMRG::VERBOSITY::OPTION::SILENT);
482// for(std::size_t loc=0; loc<this->N_sites; ++loc)
483// {
484// Mout.set_qPhys(loc, (Bsub[loc].get_basis().combine(Bimp[loc].get_basis())).qloc());
485// for(std::size_t alpha=0; alpha<Bsub[loc].orbitals(); ++alpha)
486// {
487// Mout.push(loc, {OperatorType::outerprod(Bsub[loc].S(alpha), Bimp[loc].Id(), {3}).plain<double>()});
488// Mout.push(loc, {OperatorType::outerprod(Bsub[loc].Id(), Bimp[loc].S(alpha), {3}).plain<double>()});
489// }
490// }
491// Mout.finalize();
492// return Mout;
493// }
494
495// Mpo<Sym::SU2<Sym::SpinSU2>> KondoNecklace::
496// Simp(std::size_t locx, std::size_t locy)
497// {
498// assert(locx<this->N_sites);
499// assert(locy<Bimp[locx].orbitals());
500// std::stringstream ss;
501// ss << "S_imp(" << locx;
502// if(Bimp[locx].orbitals() > 1)
503// {
504// ss << "," << locy;
505// }
506// ss << ")";
507// Mpo<Symmetry> Mout(N_sites, {3}, ss.str(), false, false, BC::OPEN, DMRG::VERBOSITY::OPTION::SILENT);
508// for(std::size_t loc=0; loc<this->N_sites; ++loc)
509// {
510// Mout.set_qPhys(loc, (Bsub[loc].get_basis().combine(Bimp[loc].get_basis())).qloc());
511// }
512// Mout.push(locx, {OperatorType::outerprod(Bsub[locx].Id(), Bimp[locx].S(locy), {3}).plain<double>()});
513// Mout.finalize();
514// return Mout;
515// }
516
517// Mpo<Sym::SU2<Sym::SpinSU2>> KondoNecklace::
518// Ssub(std::size_t locx, std::size_t locy)
519// {
520// assert(locx<this->N_sites);
521// assert(locy<Bsub[locx].orbitals());
522// std::stringstream ss;
523// ss << "S_sub(" << locx;
524// if(Bsub[locx].orbitals() > 1)
525// {
526// ss << "," << locy;
527// }
528// ss << ")";
529// Mpo<Symmetry> Mout(N_sites, {3}, ss.str(), false, false, BC::OPEN, DMRG::VERBOSITY::OPTION::SILENT);
530// for(std::size_t loc=0; loc<this->N_sites; ++loc)
531// {
532// Mout.set_qPhys(loc, (Bsub[loc].get_basis().combine(Bimp[loc].get_basis())).qloc());
533// }
534// Mout.push(locx, {OperatorType::outerprod(Bsub[locx].S(locy), Bimp[locx].Id(), {3}).plain<double>()});
535// Mout.finalize();
536// return Mout;
537// }
538
539// Mpo<Sym::SU2<Sym::SpinSU2>> KondoNecklace::
540// Simpdag(std::size_t locx, std::size_t locy)
541// {
542// assert(locx<this->N_sites);
543// assert(locy<Bimp[locx].orbitals());
544// std::stringstream ss;
545// ss << "S_imp†(" << locx;
546// if(Bimp[locx].orbitals() > 1)
547// {
548// ss << "," << locy;
549// }
550// ss << ")";
551// Mpo<Symmetry> Mout(N_sites, {3}, ss.str(), false, false, BC::OPEN, DMRG::VERBOSITY::OPTION::SILENT);
552// for(std::size_t loc=0; loc<this->N_sites; ++loc)
553// {
554// Mout.set_qPhys(loc, (Bsub[loc].get_basis().combine(Bimp[loc].get_basis())).qloc());
555// }
556// Mout.push(locx, {OperatorType::outerprod(Bsub[locx].Id(), Bimp[locx].Sdag(locy), {3}).plain<double>()});
557// Mout.finalize();
558// return Mout;
559// }
560
561// Mpo<Sym::SU2<Sym::SpinSU2>> KondoNecklace::
562// Ssubdag(std::size_t locx, std::size_t locy)
563// {
564// assert(locx<this->N_sites);
565// assert(locy<Bsub[locx].orbitals());
566// std::stringstream ss;
567// ss << "S_sub†(" << locx;
568// if(Bsub[locx].orbitals() > 1)
569// {
570// ss << "," << locy;
571// }
572// ss << ")";
573// Mpo<Symmetry> Mout(N_sites, {3}, ss.str(), false, false, BC::OPEN, DMRG::VERBOSITY::OPTION::SILENT);
574// for(std::size_t loc=0; loc<this->N_sites; ++loc)
575// {
576// Mout.set_qPhys(loc, (Bsub[loc].get_basis().combine(Bimp[loc].get_basis())).qloc());
577// }
578// Mout.push(locx, {OperatorType::outerprod(Bsub[locx].Sdag(locy), Bimp[locx].Id(), {3}).plain<double>()});
579// Mout.finalize();
580// return Mout;
581// }
582
583// Mpo<Sym::SU2<Sym::SpinSU2>> KondoNecklace::
584// SimpdagSimp(std::size_t locx1, std::size_t locx2, std::size_t locy1, std::size_t locy2)
585// {
586// assert(locx1<this->N_sites);
587// assert(locy1<Bimp[locx1].orbitals());
588// assert(locx2<this->N_sites);
589// assert(locy2<Bimp[locx2].orbitals());
590// std::stringstream ss;
591// ss << "S_imp†(" << locx1;
592// if(Bimp[locx1].orbitals() > 1)
593// {
594// ss << "," << locy1;
595// }
596// ss << ")*S_imp(" << locx2;
597// if(Bimp[locx2].orbitals() > 1)
598// {
599// ss << "," << locy2;
600// }
601// ss << ")";
602// Mpo<Symmetry> Mout(N_sites, {1}, ss.str(), true, false, BC::OPEN, DMRG::VERBOSITY::OPTION::SILENT);
603// for(std::size_t loc=0; loc<this->N_sites; ++loc)
604// {
605// Mout.set_qPhys(loc, (Bsub[loc].get_basis().combine(Bimp[loc].get_basis())).qloc());
606// }
607// std::vector<SiteOperator<Symmetry,double>> ops;
608// if(locx1 == locx2)
609// {
610// ops.push_back(OperatorType::outerprod(Bsub[locx1].Id(), OperatorType::prod(Bimp[locx1].Sdag(locy1), Bimp[locx1].S(locy2), {1}), {1}).plain<double>());
611// Mout.push(locx1, ops, {1});
612// }
613// else if(locx1 < locx2)
614// {
615// ops.push_back(OperatorType::outerprod(Bsub[locx1].Id(), Bimp[locx1].Sdag(locy1), {3}).plain<double>());
616// for(std::size_t loc=locx1+1; loc<locx2; ++loc)
617// {
618// ops.push_back(OperatorType::outerprod(Bsub[loc].Id(), Bimp[loc].Id(), {1}).plain<double>());
619// }
620// ops.push_back(OperatorType::outerprod(Bsub[locx2].Id(), Bimp[locx2].S(locy2), {3}).plain<double>());
621// Mout.push(locx1, ops, {1});
622// }
623// else
624// {
625// std::stringstream ss2;
626// ss2 << "S_imp(" << locx2;
627// if(Bimp[locx2].orbitals() > 1)
628// {
629// ss2 << "," << locy2;
630// }
631// ss2 << ")*S_imp†(" << locx1;
632// if(Bimp[locx1].orbitals() > 1)
633// {
634// ss2 << "," << locy1;
635// }
636// ss2 << ")";
637// Mout.set_name(ss2.str());
638// ops.push_back(OperatorType::outerprod(Bsub[locx2].Id(), Bimp[locx2].S(locy2), {3}).plain<double>());
639// for(std::size_t loc=locx2+1; loc<locx1; ++loc)
640// {
641// ops.push_back(OperatorType::outerprod(Bsub[loc].Id(), Bimp[loc].Id(), {1}).plain<double>());
642// }
643// ops.push_back(OperatorType::outerprod(Bsub[locx1].Id(), Bimp[locx1].Sdag(locy1), {3}).plain<double>());
644// Mout.push(locx2, ops, {1});
645// }
646// Mout.finalize();
647// return Mout;
648// }
649
650// Mpo<Sym::SU2<Sym::SpinSU2>> KondoNecklace::
651// SsubdagSsub(std::size_t locx1, std::size_t locx2, std::size_t locy1, std::size_t locy2)
652// {
653// assert(locx1<this->N_sites);
654// assert(locy1<Bsub[locx1].orbitals());
655// assert(locx2<this->N_sites);
656// assert(locy2<Bsub[locx2].orbitals());
657// std::stringstream ss;
658// ss << "S_sub†(" << locx1;
659// if(Bimp[locx1].orbitals() > 1)
660// {
661// ss << "," << locy1;
662// }
663// ss << ")*S_sub(" << locx2;
664// if(Bimp[locx2].orbitals() > 1)
665// {
666// ss << "," << locy2;
667// }
668// ss << ")";
669// Mpo<Symmetry> Mout(N_sites, {1}, ss.str(), true, false, BC::OPEN, DMRG::VERBOSITY::OPTION::SILENT);
670// for(std::size_t loc=0; loc<this->N_sites; ++loc)
671// {
672// Mout.set_qPhys(loc, (Bsub[loc].get_basis().combine(Bimp[loc].get_basis())).qloc());
673// }
674// std::vector<SiteOperator<Symmetry,double>> ops;
675// if(locx1 == locx2)
676// {
677// ops.push_back(OperatorType::outerprod(OperatorType::prod(Bsub[locx1].Sdag(locy1),Bsub[locx1].S(locy2),{1}), Bimp[locx1].Id(), {1}).plain<double>());
678// Mout.push(locx1, ops, {1});
679// }
680// else if(locx1 < locx2)
681// {
682// ops.push_back(OperatorType::outerprod(Bsub[locx1].Sdag(locy1), Bimp[locx1].Id(), {3}).plain<double>());
683// for(std::size_t loc=locx1+1; loc<locx2; ++loc)
684// {
685// ops.push_back(OperatorType::outerprod(Bsub[loc].Id(), Bimp[loc].Id(), {1}).plain<double>());
686// }
687// ops.push_back(OperatorType::outerprod(Bsub[locx2].S(locy2), Bimp[locx2].Id(), {3}).plain<double>());
688// Mout.push(locx1, ops, {1});
689// }
690// else
691// {
692// std::stringstream ss2;
693// ss2 << "S_sub(" << locx2;
694// if(Bimp[locx2].orbitals() > 1)
695// {
696// ss2 << "," << locy2;
697// }
698// ss2 << ")*S_sub†(" << locx1;
699// if(Bimp[locx1].orbitals() > 1)
700// {
701// ss2 << "," << locy1;
702// }
703// ss2 << ")";
704// Mout.set_name(ss2.str());
705// ops.push_back(OperatorType::outerprod(Bsub[locx2].S(locy2), Bimp[locx2].Id(), {3}).plain<double>());
706// for(std::size_t loc=locx2+1; loc<locx1; ++loc)
707// {
708// ops.push_back(OperatorType::outerprod(Bsub[loc].Id(), Bimp[loc].Id(), {1}).plain<double>());
709// }
710// ops.push_back(OperatorType::outerprod(Bsub[locx1].Sdag(locy1), Bimp[locx1].Id(), {3}).plain<double>());
711// Mout.push(locx2, ops, {1});
712// }
713// Mout.finalize();
714// return Mout;
715// }
716
717// Mpo<Sym::SU2<Sym::SpinSU2>> KondoNecklace::
718// SsubdagSimp(std::size_t locx1, std::size_t locx2, std::size_t locy1, std::size_t locy2)
719// {
720// assert(locx1<this->N_sites);
721// assert(locy1<Bsub[locx1].orbitals());
722// assert(locx2<this->N_sites);
723// assert(locy2<Bimp[locx2].orbitals());
724// std::stringstream ss;
725// ss << "S_sub†(" << locx1;
726// if(Bimp[locx1].orbitals() > 1)
727// {
728// ss << "," << locy1;
729// }
730// ss << ")*S_imp(" << locx2;
731// if(Bimp[locx2].orbitals() > 1)
732// {
733// ss << "," << locy2;
734// }
735// ss << ")";
736// Mpo<Symmetry> Mout(N_sites, {1}, ss.str(), false, false, BC::OPEN, DMRG::VERBOSITY::OPTION::SILENT);
737// for(std::size_t loc=0; loc<this->N_sites; ++loc)
738// {
739// Mout.set_qPhys(loc, (Bsub[loc].get_basis().combine(Bimp[loc].get_basis())).qloc());
740// }
741// std::vector<SiteOperator<Symmetry,double>> ops;
742// if(locx1 == locx2)
743// {
744// ops.push_back(OperatorType::outerprod(Bsub[locx1].Sdag(locy1), Bimp[locx1].S(locy2), {1}).plain<double>());
745// Mout.push(locx1, ops, {1});
746// }
747// else if(locx1 < locx2)
748// {
749// ops.push_back(OperatorType::outerprod(Bsub[locx1].Sdag(locy1), Bimp[locx1].Id(), {3}).plain<double>());
750// for(std::size_t loc=locx1+1; loc<locx2; ++loc)
751// {
752// ops.push_back(OperatorType::outerprod(Bsub[loc].Id(), Bimp[loc].Id(), {1}).plain<double>());
753// }
754// ops.push_back(OperatorType::outerprod(Bsub[locx2].Id(), Bimp[locx2].S(locy2), {3}).plain<double>());
755// Mout.push(locx1, ops, {1});
756// }
757// else
758// {
759// std::stringstream ss2;
760// ss2 << "S_imp(" << locx2;
761// if(Bimp[locx2].orbitals() > 1)
762// {
763// ss2 << "," << locy2;
764// }
765// ss2 << ")*S_sub†(" << locx1;
766// if(Bimp[locx1].orbitals() > 1)
767// {
768// ss2 << "," << locy1;
769// }
770// ss2 << ")";
771// Mout.set_name(ss2.str());
772// ops.push_back(OperatorType::outerprod(Bsub[locx2].Id(), Bimp[locx2].S(locy2), {3}).plain<double>());
773// for(std::size_t loc=locx2+1; loc<locx1; ++loc)
774// {
775// ops.push_back(OperatorType::outerprod(Bsub[loc].Id(), Bimp[loc].Id(), {1}).plain<double>());
776// }
777// ops.push_back(OperatorType::outerprod(Bsub[locx1].Sdag(locy1), Bimp[locx1].Id(), {3}).plain<double>());
778// Mout.push(locx2, ops, {1});
779// }
780// Mout.finalize();
781// return Mout;
782// }
783
784// Mpo<Sym::SU2<Sym::SpinSU2>> KondoNecklace::
785// SimpdagSsub(std::size_t locx1, std::size_t locx2, std::size_t locy1, std::size_t locy2)
786// {
787// assert(locx1<this->N_sites);
788// assert(locy1<Bimp[locx1].orbitals());
789// assert(locx2<this->N_sites);
790// assert(locy2<Bsub[locx2].orbitals());
791// std::stringstream ss;
792// ss << "S_imp†(" << locx1;
793// if(Bimp[locx1].orbitals() > 1)
794// {
795// ss << "," << locy1;
796// }
797// ss << ")*S_sub(" << locx2;
798// if(Bimp[locx2].orbitals() > 1)
799// {
800// ss << "," << locy2;
801// }
802// ss << ")";
803// Mpo<Symmetry> Mout(N_sites, {1}, ss.str(), false, false, BC::OPEN, DMRG::VERBOSITY::OPTION::SILENT);
804// for(std::size_t loc=0; loc<this->N_sites; ++loc)
805// {
806// Mout.set_qPhys(loc, (Bsub[loc].get_basis().combine(Bimp[loc].get_basis())).qloc());
807// }
808// std::vector<SiteOperator<Symmetry,double>> ops;
809// if(locx1 == locx2)
810// {
811// std::stringstream ss2;
812// ss2 << "S_sub(" << locx2;
813// if(Bimp[locx2].orbitals() > 1)
814// {
815// ss2 << "," << locy2;
816// }
817// ss2 << ")*S_imp†(" << locx1;
818// if(Bimp[locx1].orbitals() > 1)
819// {
820// ss2 << "," << locy1;
821// }
822// ss2 << ")";
823// Mout.set_name(ss2.str());
824// ops.push_back(OperatorType::outerprod(Bsub[locx1].S(locy2), Bimp[locx1].S(locy1), {1}).plain<double>());
825// Mout.push(locx1, ops, {1});
826// }
827// else if(locx1 < locx2)
828// {
829// ops.push_back(OperatorType::outerprod(Bsub[locx1].Id(), Bimp[locx1].Sdag(locy1), {3}).plain<double>());
830// for(std::size_t loc=locx1+1; loc<locx2; ++loc)
831// {
832// ops.push_back(OperatorType::outerprod(Bsub[loc].Id(), Bimp[loc].Id(), {1}).plain<double>());
833// }
834// ops.push_back(OperatorType::outerprod(Bsub[locx2].S(locy2), Bimp[locx2].Id(), {3}).plain<double>());
835// Mout.push(locx1, ops, {1});
836// }
837// else
838// {
839// std::stringstream ss2;
840// ss2 << "S_sub(" << locx2;
841// if(Bimp[locx2].orbitals() > 1)
842// {
843// ss2 << "," << locy2;
844// }
845// ss2 << ")*S_imp†(" << locx1;
846// if(Bimp[locx1].orbitals() > 1)
847// {
848// ss2 << "," << locy1;
849// }
850// ss2 << ")";
851// Mout.set_name(ss2.str());
852// ops.push_back(OperatorType::outerprod(Bsub[locx2].S(locy2), Bimp[locx2].Id(), {3}).plain<double>());
853// for(std::size_t loc=locx2+1; loc<locx1; ++loc)
854// {
855// ops.push_back(OperatorType::outerprod(Bsub[loc].Id(), Bimp[loc].Id(), {1}).plain<double>());
856// }
857// ops.push_back(OperatorType::outerprod(Bsub[locx1].Id(), Bimp[locx1].Sdag(locy1), {3}).plain<double>());
858// Mout.push(locx2, ops, {1});
859// }
860// Mout.finalize();
861// return Mout;
862// }
863
864} // end namespace VMPS
865
866
867
868#endif
frac
boost::rational< int > frac
Definition: DmrgExternal.h:11
print_frac_nice
std::string print_frac_nice(frac r)
Definition: DmrgExternal.h:32
MODEL_FAMILY
MODEL_FAMILY
Definition: DmrgTypedefs.h:96
HEISENBERG
@ HEISENBERG
Definition: DmrgTypedefs.h:96
SZ
@ SZ
Definition: DmrgTypedefs.h:60
SP
@ SP
Definition: DmrgTypedefs.h:60
SM
@ SM
Definition: DmrgTypedefs.h:60
BC
BC
Definition: DmrgTypedefs.h:161
EIGEN_DEFAULT_SPARSE_INDEX_TYPE
#define EIGEN_DEFAULT_SPARSE_INDEX_TYPE
Definition: DmrgTypedefs.h:181
kroneckerProduct
SiteOperator< Symmetry, Scalar_ > kroneckerProduct(const SiteOperator< Symmetry, Scalar_ > &O1, const SiteOperator< Symmetry, Scalar_ > &O2)
Definition: SiteOperator.h:164
KondoNecklaceObservables< Sym::U0 >::Scomp
std::enable_if<!Dummy::IS_SPIN_SU2(), Mpo< Sym::U0 > >::type Scomp(SPINOP_LABEL Sa, size_t locx, size_t locy=0, double factor=1.) const
Definition: KondoNecklaceObservables.h:41
MpoTerms::N_phys
std::size_t N_phys
Definition: MpoTerms.h:400
MpoTerms::finalize
void finalize(const bool COMPRESS=true, const std::size_t power=1, const double tolerance=::mynumeric_limits< double >::epsilon())
Definition: MpoTerms.h:1281
MpoTerms::set_verbosity
void set_verbosity(const DMRG::VERBOSITY::OPTION VERB_in)
Definition: MpoTerms.h:881
MpoTerms::N_sites
std::size_t N_sites
Definition: MpoTerms.h:395
MpoTerms::setLocBasis
void setLocBasis(const std::vector< std::vector< qType > > &q)
Definition: MpoTerms.h:715
MpoTerms::VERB
DMRG::VERBOSITY::OPTION VERB
Definition: MpoTerms.h:102
MpoTerms::label
std::string label
Definition: MpoTerms.h:615
Mpo
Definition: Mpo.h:40
Mpo< Sym::U0, double >::precalc_TwoSiteData
void precalc_TwoSiteData(bool FORCE=false)
Mpo< Sym::U0, double >::construct_from_pushlist
void construct_from_pushlist(const PushType< OperatorType, CouplScalar > &pushlist, const std::vector< std::vector< std::string > > &labellist, size_t Lcell)
Sym::U0
Definition: U0.h:28
VMPS::KondoNecklace::MatrixType
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixType
Definition: KondoNecklace.h:28
VMPS::KondoNecklace::defaults
static const std::map< string, std::any > defaults
Definition: KondoNecklace.h:59
VMPS::KondoNecklace::sweep_defaults
static const std::map< string, std::any > sweep_defaults
Definition: KondoNecklace.h:64
VMPS::KondoNecklace::FAMILY
static constexpr MODEL_FAMILY FAMILY
Definition: KondoNecklace.h:24
VMPS::KondoNecklace::SparseMatrixType
Eigen::SparseMatrix< double, Eigen::ColMajor, EIGEN_DEFAULT_SPARSE_INDEX_TYPE > SparseMatrixType
Definition: KondoNecklace.h:29
VMPS::KondoNecklace::validate
bool validate(qType qnum)
VMPS::KondoNecklace::set_operators
static void set_operators(const std::vector< SpinBase< Symmetry > > &Bsub, const std::vector< SpinBase< Symmetry > > &Bimp, const ParamHandler &P, PushType< SiteOperator< Symmetry, double >, double > &pushlist, std::vector< std::vector< std::string > > &labellist, const BC boundary)
Definition: KondoNecklace.h:116
MAKE_TYPEDEFS
#define MAKE_TYPEDEFS(MODEL)
Definition: macros.h:4
PROP::COMPRESS
const bool COMPRESS
Definition: DmrgTypedefs.h:499
qarray
Definition: qarray.h:26