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