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