VMPS++
Loading...
Searching...
No Matches
Mpo.h
Go to the documentation of this file.
1#ifndef STRAWBERRY_Mpo_WITH_Q
2#define STRAWBERRY_Mpo_WITH_Q
3
5#include "boost/multi_array.hpp"
7
8#include "termcolor.hpp" //from https://github.com/ikalnytskyi/termcolor
9
11#include <Eigen/SparseCore>
13
14#ifndef EIGEN_DEFAULT_SPARSE_INDEX_TYPE
15#define EIGEN_DEFAULT_SPARSE_INDEX_TYPE int
16#endif
17typedef Eigen::SparseMatrix<double,Eigen::ColMajor,EIGEN_DEFAULT_SPARSE_INDEX_TYPE> SparseMatrixXd;
18using namespace Eigen;
19
21#include <unsupported/Eigen/KroneckerProduct>
23
24#include "util/macros.h"
25#include "MpoTerms.h"
26#include "tensors/Qbasis.h"
27#include "DmrgTypedefs.h"
28
29namespace VMPS{};
30
31template<typename Symmetry, typename Scalar> class Mps;
32template<typename Symmetry, typename Scalar> class Umps;
33template<typename Symmetry, typename MpHamiltonian, typename Scalar> class DmrgSolver;
34template<typename Symmetry, typename Scalar, typename MpoScalar> class MpsCompressor;
35template<typename Symmetry, typename MpHamiltonian, typename Scalar> class VumpsSolver;
36
37
38template<typename Symmetry, typename Scalar=double>
39class Mpo : public MpoTerms<Symmetry,Scalar>
40{
41public:
44 static constexpr size_t Nq = Symmetry::Nq;
45 typedef typename Symmetry::qType qType;
46 typename Symmetry::qType qVac = Symmetry::qvacuum();
47
48 template<typename Symmetry_, typename MpHamiltonian, typename Scalar_> friend class DmrgSolver;
49 template<typename Symmetry_, typename MpHamiltonian, typename Scalar_> friend class VumpsSolver;
50 template<typename Symmetry_, typename S1, typename S2> friend class MpsCompressor;
51 template<typename H, typename Symmetry_, typename S1, typename S2, typename V> friend class TDVPPropagator;
52 template<typename Symmetry_, typename S_> friend class Mpo;
53
54public:
55
56 typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
57 typedef Scalar Scalar_;
58
59 Mpo() : MpoTerms<Symmetry, Scalar>(){};
60
62 :MpoTerms<Symmetry,Scalar>(Terms)
63 {};
64
65 Mpo (size_t L_input);
66
67 Mpo (std::size_t L_input, qType Qtot_input, std::string label_input="Mpo", bool HERMITIAN_input=false, bool UNITARY_input=false, BC BC_input=BC::OPEN, DMRG::VERBOSITY::OPTION VERB_input=DMRG::VERBOSITY::OPTION::SILENT);
68
69 template<typename CouplScalar>
70 void construct_from_pushlist(const PushType<OperatorType,CouplScalar>& pushlist, const std::vector<std::vector<std::string>>& labellist, size_t Lcell);
71
72 template<typename CouplScalar>
73 void calc_reversedData_from_pushlist(const PushType<OperatorType,CouplScalar>& pushlist, double tolerance=::mynumeric_limits<double>::epsilon());
74
75 void push_qpath (const std::size_t loc, const std::vector<OperatorType> &opList, const std::vector<qType> &qList, const Scalar lambda = 1.0)
76 {
77 MpoTerms<Symmetry,Scalar>::push(loc,opList,qList,lambda);
78 }
79
80 void setLocal(std::size_t loc, const OperatorType& op);
81
82 void setLocal(std::size_t loc, const OperatorType& op, const OperatorType& signOp);
83
84 void setLocal(std::size_t loc, const OperatorType& op, const std::vector<OperatorType>& signOp);
85
86 void setLocal(const std::vector<std::size_t>& locs, const std::vector<OperatorType>& ops);
87
88 void setLocal(const std::vector<std::size_t>& locs, const std::vector<OperatorType>& ops, const OperatorType& signOp);
89
90 void setLocal(const std::vector<std::size_t>& locs, const std::vector<OperatorType>& ops, const std::vector<OperatorType>& signOps);
91
92 void setLocalStag(std::size_t loc, const OperatorType& op, const std::vector<OperatorType>& stagSignOps);
93
94 void setLocalSum(const OperatorType& op, Scalar (*f)(int)=localSumTrivial);
95
96 void setLocalSum(const std::vector<OperatorType>& op, std::vector<Scalar> coeffs);
97
98 void setProductSum(const OperatorType& op1, const OperatorType& op2);
99
100 // void scale(double factor=1., double offset=0.);
101
102 void precalc_TwoSiteData(bool FORCE=false);
103
104 std::string info (bool REDUCED=false) const;
105
106 int get_dAux_max (int power=1) const;
107
108 //double memory(MEMUNIT memunit=GB) const;
109
110 //double sparsity(bool USE_SQUARE=false, bool PER_MATRIX=true) const;
111
112 inline std::size_t length() const {return this->size();}
113
114 inline std::size_t volume() const {return this->N_phys;}
115
116 template<typename T, typename ... Operator>
117 static std::vector<T> get_N_site_interaction(T const & Op0, Operator const & ... Ops) {std::vector<T> out { {Op0, Ops ...} }; return out;};
118 //inline std::size_t auxrows(std::size_t loc) const {return this->get_qAux()[loc].fullM();}
119 //inline std::size_t auxcols(std::size_t loc) const {return this->get_qAux()[loc+1].fullM();}
120 //inline void setOpBasis (const vector<vector<qType> > &q) {qOp=q;}
121 //inline void setOpBasisSq (const vector<vector<qType> > &qOpSq_in) {qOpSq=qOpSq_in;}
122 //inline const unordered_map<tuple<size_t,size_t,size_t,qarray<Symmetry::Nq>,qarray<Symmetry::Nq> >,SparseMatrix<Scalar> > &Vsq_at (size_t loc) const {return Vsq[loc];};
123
124 inline int locality() const {return LocalSite;}
125 inline void set_locality(std::size_t LocalSite_input) {LocalSite = LocalSite_input;}
126 inline OperatorType localOperator() const {return LocalOp;}
127 inline void set_localOperator (OperatorType LocalOp_input) {LocalOp = LocalOp_input;}
128
129 static Mpo<Symmetry,Scalar> Identity(const std::vector<std::vector<qType>>& qPhys, const qType& Q=Symmetry::qvacuum()); // static Identity, needs basis
130 Mpo<Symmetry,Scalar> Identity() const {return Identity(this->locBasis());}; // Identity for basis of the given class
131
132 static Mpo<Symmetry,Scalar> Zero(const std::vector<std::vector<qType>>& qPhys);
133 Mpo<Symmetry,Scalar> Zero() const {return Zero(this->locBasis());};
134
135 inline bool IS_UNITARY() const {return UNITARY;};
136
137 inline bool IS_HERMITIAN() const {return HERMITIAN;};
138
139 inline bool IS_HAMILTONIAN() const {return HAMILTONIAN;};
140
141 inline bool HAS_TWO_SITE_DATA() const {return GOT_TWO_SITE_DATA;};
142
143 inline bool HAS_W() const {return this->GOT_W;};
144
145 boost::multi_array<Scalar,4> H2site(std::size_t loc, bool HALF_THE_LOCAL_TERM=false) const {assert(false and "Method H2site is deprecated.");}
146
154
155 bool GOT_TWO_SITE_DATA = false;
156 std::vector<std::vector<TwoSiteData<Symmetry,Scalar>>> TSD;
157
158 bool UNITARY = false;
159 bool HERMITIAN = false;
160 bool HAMILTONIAN = false;
161 bool GOT_SEMIOPEN_LEFT = false;
162 bool GOT_SEMIOPEN_RIGHT = false;
163
165 int LocalSite = -1;
166
167 //std::size_t N_phys = 0;
168
169 //void initialize();
170
171 void generate_label(std::size_t Lcell);
172
173
185 void push_width(const std::size_t n, const std::size_t loc, const Scalar lambda, const OperatorType& outOp, const std::vector<OperatorType>& trans, const OperatorType& inOp);
186
195 void push_local(const std::size_t loc, const Scalar lambda, const OperatorType& op);
196
206 void push_tight(const std::size_t loc, const Scalar lambda, const OperatorType& op1, const OperatorType& op2);
207
208
219 void push_nextn(const std::size_t loc, const Scalar lambda, const OperatorType& op1, const OperatorType& trans, const OperatorType& op2);
220
222};
223
224template<typename Symmetry, typename Scalar>
226Mpo(std::size_t L_input)
227: MpoTerms<Symmetry, Scalar>(L_input){}
228
229template<typename Symmetry, typename Scalar>
231Mpo (std::size_t L_input, qType Qtot_input, string label_input, bool HERMITIAN_input, bool UNITARY_input, BC BC_input, DMRG::VERBOSITY::OPTION VERB_input)
232: MpoTerms<Symmetry, Scalar>(L_input, BC_input, Qtot_input, VERB_input), HERMITIAN(HERMITIAN_input), UNITARY(UNITARY_input)
233{
234 this->set_name(label_input);
235}
236
237template<typename Symmetry, typename Scalar> void Mpo<Symmetry,Scalar>::
238push_width(const std::size_t width, const std::size_t loc, const Scalar lambda, const OperatorType& outOp, const std::vector<OperatorType>& trans, const OperatorType& inOp)
239{
240 std::vector<OperatorType> oplist(0);
241 oplist.push_back(outOp);
242 for(std::size_t m=0; m<trans.size(); ++m)
243 {
244 oplist.push_back(trans[m]);
245 }
246 oplist.push_back(inOp);
247 this->push(loc, oplist, lambda);
248}
249
250template<typename Symmetry, typename Scalar> void Mpo<Symmetry,Scalar>::
251push_local(const std::size_t loc, const Scalar lambda, const OperatorType& op)
252{
253 this->push(loc, {op}, lambda);
254}
255
256template<typename Symmetry, typename Scalar> void Mpo<Symmetry,Scalar>::
257push_tight(const std::size_t loc, const Scalar lambda, const OperatorType& op1, const OperatorType& op2)
258{
259 this->push(loc, {op1, op2}, lambda);
260}
261
262template<typename Symmetry, typename Scalar> void Mpo<Symmetry,Scalar>::
263push_nextn(const std::size_t loc, const Scalar lambda, const OperatorType& op1, const OperatorType& trans, const OperatorType& op2)
264{
265 this->push(loc, {op1, trans, op2}, lambda);
266}
267
268
269template<typename Symmetry, typename Scalar>
271info (bool REDUCED) const
272{
273 std::stringstream ss;
274 ss << termcolor::colorize << termcolor::bold;
275 if (!REDUCED)
276 {
277 ss << this->get_name();
278 }
279 else
280 {
281 ss << "Mpo";
282 }
283 ss << termcolor::reset << "→ L=" << this->size();
284 if (this->N_phys > this->size()) ss << ",V=" << this->N_phys;
285 ss << ", " << Symmetry::name() << ", ";
286
287 ss << "UNITARY=" << boolalpha << UNITARY << ", ";
288 ss << "HERMITIAN=" << boolalpha << HERMITIAN << ", ";
289 ss << "maxPower=" << this->maxPower() << ", ";
290 ss << "BC=" << this->get_boundary_condition() << ", ";
291 ss << "2SITE_DATA=" << boolalpha << GOT_TWO_SITE_DATA << ", ";
292 if(LocalSite != -1)
293 {
294 ss << "locality=" << LocalSite << ", ";
295 }
296
297 auto print_qaux = [this,&ss] (int power) -> void
298 {
299 auto qAux = this->get_qAux_power(power);
300 std::vector<int> dAux(this->size()+1);
301// dAux[0] = this->auxBasis(0).fullM();
302 dAux[0] = qAux[0].fullM();
303 std::set<std::pair<int,int> > dAux_set;
304 for (std::size_t loc=0; loc<this->size(); ++loc)
305 {
306// dAux[loc+1] = this->auxBasis(loc+1).fullM();
307 dAux[loc+1] = qAux[loc+1].fullM();
308 dAux_set.insert(std::make_pair(dAux[loc],dAux[loc+1]));
309 }
310 ss << "dAux" << power << "=";
311 for (const auto& dAux_pair : dAux_set)
312 {
313 ss << dAux_pair.first << "x" << dAux_pair.second;
314 ss << ",";
315 }
316 ss << " ";
317 };
318
319 for (int power=1; power<=this->maxPower(); ++power) print_qaux(power);
320
321 ss << "mem=" << round(this->memory(GB),3) << "GB";
322 if (this->GOT_W and this->GOT_QAUX) ss << ", sparsity=" << this->sparsity();
323 // if(this->check_SQUARE())
324 // {
325 // ss << ", sparsity(sq)=" << sparsity(true);
326 // }
327 return ss.str();
328}
329
330template<typename Symmetry, typename Scalar>
332get_dAux_max (int power) const
333{
334 VectorXi res(this->size());
335 auto qAux = this->get_qAux_power(power);
336 for (std::size_t l=0; l<this->size(); ++l)
337 {
338 res(l) = qAux[l].fullM();
339 }
340 return res.maxCoeff();
341}
342
343/*template<typename Symmetry, typename Scalar>
344double Mpo<Symmetry,Scalar>::
345memory(MEMUNIT memunit) const
346{
347 double res = 0.;
348
349 for (std::size_t loc=0; loc<this->W.size(); ++loc)
350 {
351 std::size_t hd = this->get_hilbert_dimension(loc);
352 for (std::size_t n1=0; n1<this->W[loc].size(); ++n1)
353 {
354 for (std::size_t n2=0; n2<this->W[loc][n1].size(); ++n2)
355 {
356 for (std::size_t t=0; t<this->W[loc][n1][n2].size(); ++t)
357 {
358 res += this->W[loc][n1][n2][t].memory(memunit);
359 }
360 }
361 }
362 }
363 return res;
364}*/
365
366/*template<typename Symmetry, typename Scalar>
367double Mpo<Symmetry,Scalar>::
368sparsity (bool USE_SQUARE, bool PER_MATRIX) const
369{
370 if (USE_SQUARE) {assert(this->check_power(2ul));}
371 double N_nonZeros = 0.;
372 double N_elements = 0.;
373 double N_matrices = 0.;
374
375 for (std::size_t loc=0; loc<this->size(); ++loc)
376 {
377 std::size_t hd = this->get_hilbert_dimension(loc);
378 N_matrices += hd * hd * this->opBasis(loc).size();
379
380 for (size_t s1=0; s1<qloc[l].size(); ++s1)
381 for (size_t s2=0; s2<qloc[l].size(); ++s2)
382 for (size_t k=0; k<this->qOp[l].size(); ++k)
383 {
384 // if constexpr (Symmetry::NON_ABELIAN) {N_nonZeros += this->W[l][s1][s2][k].nonZeros();}
385 // if constexpr (Symmetry::NON_ABELIAN) {N_elements += this->W[l][s1][s2][k].rows() * this->W[l][s1][s2][k].cols();}
386 // else
387 // {
388 // N_nonZeros += (USE_SQUARE)? Wsq[l][s1][s2][k].nonZeros() : this->W[l][s1][s2][k].nonZeros();
389 // N_elements += (USE_SQUARE)? Wsq[l][s1][s2][k].rows() * Wsq[l][s1][s2][k].cols():
390 // this->W[l][s1][s2][k].rows() * this->W[l][s1][s2][k].cols();
391 // }
392 }
393 }
394
395 //return (PER_MATRIX)? N_nonZeros/N_matrices : N_nonZeros/N_elements;
396 return 0.;
397}*/
398
399template<typename Symmetry, typename Scalar> Mpo<Symmetry,Scalar> Mpo<Symmetry,Scalar>::
400Identity(const std::vector<std::vector<qType>>& qPhys, const typename Symmetry::qType& Q)
401{
402 Mpo<Symmetry,Scalar> out(qPhys.size(), Symmetry::qvacuum(), "Id", true, true, BC::OPEN, DMRG::VERBOSITY::OPTION::SILENT);
403 for(std::size_t loc=0; loc<qPhys.size(); ++loc)
404 {
405 out.set_qPhys(loc, qPhys[loc]);
406 }
407 out.set_Identity(Q);
408 return out;
409}
410
411template<typename Symmetry, typename Scalar> Mpo<Symmetry,Scalar> Mpo<Symmetry,Scalar>::
412Zero(const std::vector<std::vector<qType>>& qPhys)
413{
414 Mpo<Symmetry,Scalar> out(qPhys.size(), Symmetry::qvacuum(), "Zero", true, true, BC::OPEN, DMRG::VERBOSITY::OPTION::SILENT);
415 for(std::size_t loc=0; loc<qPhys.size(); ++loc)
416 {
417 out.set_qPhys(loc, qPhys[loc]);
418 }
419 out.set_Zero();
420 return out;
421}
422
423template<typename Symmetry, typename Scalar> void Mpo<Symmetry,Scalar>::
424generate_label(std::size_t Lcell)
425{
426 std::stringstream ss;
427 ss << this->get_name();
428 std::vector<std::string> info = this->get_info();
429
430 std::map<std::string,std::set<std::size_t> > cells;
431
432 for (std::size_t loc=0; loc<info.size(); ++loc)
433 {
434 cells[info[loc]].insert(loc%Lcell);
435 }
436
437 if (cells.size() == 1)
438 {
439 ss << "(" << info[0] << ")";
440 }
441 else
442 {
443 std::vector<std::pair<std::string,std::set<std::size_t> > > cells_resort(cells.begin(), cells.end());
444
445 // sort according to smallest l, not according to label
446 sort(cells_resort.begin(), cells_resort.end(),
447 [](const std::pair<std::string,std::set<std::size_t> > &a, const std::pair<std::string,std::set<std::size_t> > &b) -> bool
448 {
449 return *min_element(a.second.begin(),a.second.end()) < *min_element(b.second.begin(),b.second.end());
450 });
451
452 ss << ":";
453 for (auto c:cells_resort)
454 {
455 ss << std::endl << " •l=";
456 // for (auto s:c.second)
457 // {
458 // cout << s << ",";
459 // }
460 // cout << endl;
461 if (c.second.size() == 1)
462 {
463 ss << *c.second.begin(); // one site
464 }
465 else
466 {
467 // check mod 2
468 if (std::all_of(c.second.cbegin(), c.second.cend(), [](int i){ return i%2==0;}) and c.second.size() == this->size()/2)
469 {
470 ss << "even";
471 }
472 else if (std::all_of(c.second.cbegin(), c.second.cend(), [](int i){ return i%2==1;}) and c.second.size() == this->size()/2)
473 {
474 ss << "odd";
475 }
476 // check mod 4
477 else if (std::all_of(c.second.cbegin(), c.second.cend(), [](int i){ return i%4==0;}) and c.second.size() == this->size()/4)
478 {
479 ss << "0mod4";
480 }
481 else if (std::all_of(c.second.cbegin(), c.second.cend(), [](int i){ return i%4==1;}) and c.second.size() == this->size()/4)
482 {
483 ss << "1mod4";
484 }
485 else if (std::all_of(c.second.cbegin(), c.second.cend(), [](int i){ return i%4==2;}) and c.second.size() == this->size()/4)
486 {
487 ss << "2mod4";
488 }
489 else if (std::all_of(c.second.cbegin(), c.second.cend(), [](int i){ return i%4==3;}) and c.second.size() == this->size()/4)
490 {
491 ss << "3mod4";
492 }
493 else
494 {
495 if (c.second.size() == 2)
496 {
497 ss << *c.second.begin() << "," << *c.second.rbegin(); // two sites
498 }
499 else
500 {
501 bool CONSECUTIVE = true;
502 for (auto it=c.second.begin(); it!=c.second.end(); ++it)
503 {
504 if (next(it) != c.second.end() and *next(it)!=*it+1ul)
505 {
506 CONSECUTIVE = false;
507 }
508 }
509 if (CONSECUTIVE)
510 {
511 ss << *c.second.begin() << "-" << *c.second.rbegin(); // range of sites
512 }
513 else
514 {
515 for (auto it=c.second.begin(); it!=c.second.end(); ++it)
516 {
517 ss << *it << ","; // some unknown order
518 }
519 ss.seekp(-1,ios_base::end); // delete last comma
520 }
521 }
522 }
523 }
524 // ss.seekp(-1,ios_base::end); // delete last comma
525 ss << ": " << c.first;
526 }
527 }
528
529 this->set_name(ss.str());
530}
531
532template<typename Symmetry, typename Scalar>
533template<typename CouplScalar>
535construct_from_pushlist(const PushType<OperatorType,CouplScalar>& pushlist, const std::vector<std::vector<std::string>>& labellist, size_t Lcell)
536{
537 for(std::size_t i=0; i<pushlist.size(); ++i)
538 {
539 const auto& [loc, ops, coupling] = pushlist[i];
540 if ( std::abs(coupling) != 0. )
541 {
542 this->push(loc, ops, coupling);
543 }
544 }
545 for(std::size_t loc=0; loc<this->size(); ++loc)
546 {
547 for(std::size_t i=0; i<labellist[loc].size(); ++i)
548 {
549 this->save_label(loc, labellist[loc][i]);
550 }
551 }
552 generate_label(Lcell);
553}
554
555template<typename Symmetry, typename Scalar>
556template<typename CouplScalar>
559{
560 Mpo<Symmetry,Scalar> mpo_rev(this->N_sites, this->get_qTot(), "Reversed(temp)", false, false, this->get_boundary_condition(), DMRG::VERBOSITY::OPTION::SILENT);
561 for(std::size_t loc=0; loc<this->N_sites; ++loc)
562 {
563 mpo_rev.set_qPhys(loc, this->qPhys[this->N_sites-1-loc]);
564 }
565 for(std::size_t i=0; i<pushlist.size(); ++i)
566 {
567 const auto& [loc, ops, coupling] = pushlist[i];
568 std::size_t range = ops.size();
569 std::size_t loc_rev = this->N_sites-loc-range;
570 std::vector<OperatorType> ops_rev(range);
571 for(std::size_t j=0; j<range; ++j)
572 {
573 ops_rev[j] = ops[range-1-j];
574 }
575 if ( std::abs(coupling) != 0. )
576 {
577 mpo_rev.push(loc_rev, ops_rev, coupling);
578 }
579 }
580 mpo_rev.finalize(PROP::COMPRESS, 1, tolerance);
581 this->reversed.qAux.resize(this->N_sites+1);
582 this->reversed.W.resize(this->N_sites);
583 for(std::size_t loc=0; loc<this->N_sites; ++loc)
584 {
585 this->reversed.qAux[loc] = mpo_rev.qAux[this->N_sites-loc];
586 auto Wloc = mpo_rev.W[this->N_sites-loc-1];
587 decltype(Wloc) Wloc_rev;
588 Wloc_rev.resize(Wloc.size());
589 for(std::size_t m=0; m<Wloc.size(); ++m)
590 {
591 Wloc_rev[m].resize(Wloc[m].size());
592 for(std::size_t n=0; n<Wloc[m].size(); ++n)
593 {
594 Wloc_rev[m][n].resize(Wloc[m][n].size());
595 for(std::size_t t=0; t<Wloc[m][n].size(); ++t)
596 {
597 Wloc_rev[m][n][t] = Wloc[m][n][t].transpose();
598 }
599 }
600 }
601 this->reversed.W[loc] = Wloc_rev;
602 }
603 this->reversed.qAux[this->N_sites] = mpo_rev.qAux[0];
604 this->reversed.SET = true;
605}
606
607
608template<typename Symmetry, typename Scalar> void Mpo<Symmetry,Scalar>::
609setLocal(std::size_t loc, const OperatorType& op)
610{
611 assert(this->check_qPhys() and "Physical bases have to be set before");
612 LocalOp = op;
613 LocalSite = loc;
616 pushlist.push_back(std::make_tuple(loc, Ops, 1.));
617 std::vector<std::vector<std::string>> labellist(this->N_sites);
618 this->construct_from_pushlist(pushlist, labellist, 1);
619 this->finalize(PROP::COMPRESS, 1);
620 this->calc_reversedData_from_pushlist(pushlist);
621}
622
623template<typename Symmetry, typename Scalar> void Mpo<Symmetry,Scalar>::
624setLocal(std::size_t loc, const OperatorType& op, const OperatorType& signOp)
625{
626 std::vector<OperatorType> signOps(loc, signOp);
627 setLocal(loc, op, signOps);
628}
629
630template<typename Symmetry, typename Scalar> void Mpo<Symmetry,Scalar>::
631setLocal(std::size_t loc, const OperatorType& op, const std::vector<OperatorType>& signOps)
632{
633 assert(this->check_qPhys() and "Physical bases have to be set before");
634 assert(signOps.size() == loc and "Number of sign operators does not match the chosen lattice site");
635 LocalOp = op;
636 LocalSite = loc;
637 std::vector<OperatorType> ops = signOps;
638 ops.push_back(op);
639 this->push(0, ops);
640 this->finalize(PROP::COMPRESS, 1);
641}
642
643template<typename Symmetry, typename Scalar> void Mpo<Symmetry,Scalar>::
644setLocalStag(std::size_t loc, const OperatorType& op, const std::vector<OperatorType>& stagSignOps)
645{
646 assert(this->check_qPhys() and "Physical bases have to be set before");
647 assert(stagSignOps.size() == this->size() and "Number of staggered sign operators does not match the chosen lattice size");
648 LocalOp = op;
649 LocalSite = loc;
650
651 std::vector<OperatorType> ops(this->size());
652 ops[loc] = op;
653 for(std::size_t loc2=0; loc2<this->size(); ++loc2)
654 {
655 if(loc2 != loc)
656 {
657 ops[loc2] = stagSignOps[loc2];
658 }
659 }
660 this->push(0, ops);
661 this->finalize(PROP::COMPRESS, 1);
662}
663
664template<typename Symmetry, typename Scalar> void Mpo<Symmetry,Scalar>::
665setLocal(const std::vector<std::size_t>& locs, const std::vector<OperatorType>& ops)
666{
667 auto Id = ops[0];
668 Id.setIdentity();
669 #ifdef OPLABELS
670 Id.label = "id";
671 #endif
672 setLocal(locs, ops, Id);
673}
674
675template<typename Symmetry, typename Scalar> void Mpo<Symmetry,Scalar>::
676setLocal(const std::vector<std::size_t>& locs, const std::vector<OperatorType>& ops, const std::vector<OperatorType>& signOps)
677{
678 assert(this->check_qPhys() and "Physical bases have to be set before");
679 assert(locs.size() == 2 and "setLocal(...) only works for two operators!");
680 int left = 0;
681 int right = 1;
682 if(locs[left] > locs[right])
683 {
684 left = 1;
685 right = 0;
686 }
687 assert(locs[left] != locs[right] and "setLocal(...) needs to local operators at different sites!");
688 assert(signOps.size() == locs[right]-locs[left]-1);
689
690 std::vector<OperatorType> ops_with_signs;
691 ops_with_signs.push_back(ops[left]);
692 for(std::size_t pos=0; pos<signOps.size(); ++pos)
693 {
694 ops_with_signs.push_back(signOps[pos]);
695 }
696 ops_with_signs.push_back(ops[right]);
697 this->push(locs[left], ops_with_signs);
698 this->finalize(PROP::COMPRESS, 1);
699}
700
701// Warning: Needs homogenenous basis
702template<typename Symmetry, typename Scalar> void Mpo<Symmetry,Scalar>::
703setLocal(const std::vector<std::size_t>& locs, const std::vector<OperatorType>& ops, const OperatorType& signOp)
704{
705 assert(locs.size() == 2 and "setLocal(...) only works for two operators!");
706 assert(locs[0] != locs[1] and "setLocal(...) needs to local operators at different sites!");
707 int distance = locs[1]-locs[0];
708 if(distance<0)
709 {
710 distance *= -1;
711 }
712 std::vector<OperatorType> signOps(distance-1, signOp);
713 setLocal(locs, ops, signOps);
714}
715
716template<typename Symmetry, typename Scalar> void Mpo<Symmetry,Scalar>::
717setLocalSum(const OperatorType& op, Scalar (*f)(int))
718{
719 assert(this->check_qPhys() and "Physical bases have to be set before");
720 for (std::size_t loc=0; loc<this->size(); ++loc)
721 {
722 this->push(loc, {f(loc)*op});
723 }
724 this->finalize(PROP::COMPRESS, 1);
725}
726
727template<typename Symmetry, typename Scalar> void Mpo<Symmetry,Scalar>::
728setLocalSum(const std::vector<OperatorType>& ops, std::vector<Scalar> coeffs)
729{
730 assert(this->check_qPhys() and "Physical bases have to be set before");
732 for (std::size_t loc=0; loc<this->size(); ++loc)
733 {
734 auto Ops = Mpo<Symmetry,Scalar>::get_N_site_interaction(coeffs[loc]*ops[loc]);
735 pushlist.push_back(std::make_tuple(loc, Ops, 1.));
736 // this->push(loc, {coeffs[loc]*ops[loc]});
737 }
738 std::vector<std::vector<std::string>> labellist(this->N_sites);
739 this->construct_from_pushlist(pushlist, labellist, 1);
740 this->finalize(PROP::COMPRESS, 1);
741 this->calc_reversedData_from_pushlist(pushlist);
742}
743
744template<typename Symmetry, typename Scalar> void Mpo<Symmetry,Scalar>::
745setProductSum(const OperatorType& op1, const OperatorType& op2)
746{
747 assert(this->check_qPhys() and "Physical bases have to be set before");
748 for(std::size_t loc=0; loc<this->size()-1; ++loc)
749 {
750 this->push(loc, {op1, op2});
751 }
752 this->finalize(PROP::COMPRESS, 1);
753}
754
755// template<typename Symmetry, typename Scalar> void Mpo<Symmetry,Scalar>::
756// scale(double factor, double offset)
757// {
758// if (LocalSite != -1)
759// {
760// auto Id = LocalOp;
761// Id.setIdentity();
762// #ifdef OPLABELS
763// Id.label = "id";
764// #endif
765// if(std::abs(factor - 1.0) > ::mynumeric_limits<Scalar>::epsilon())
766// {
767// LocalOp = factor * LocalOp;
768// }
769// if(std::abs(offset) > ::mynumeric_limits<Scalar>::epsilon())
770// {
771// LocalOp += offset * Id;
772// }
773// setLocal(LocalSite, LocalOp);
774// }
775// else
776// {
777// this->scale(factor, offset);
778// this->finalize(PROP::COMPRESS, 1);
779// }
780// }
781
782template<typename Symmetry, typename Scalar>
784precalc_TwoSiteData (bool FORCE)
785{
786 if(!GOT_TWO_SITE_DATA or FORCE)
787 {
788 TSD.clear();
789 TSD = this->calc_TwoSiteData();
790 GOT_TWO_SITE_DATA = true;
791 }
792}
793
794template<typename Symmetry, typename Scalar>
795std::ostream &operator<<(std::ostream& os, const Mpo<Symmetry,Scalar>& O)
796{
797 os << setfill('-') << setw(30) << "-" << setfill(' ');
798 os << "Mpo: L=" << O.length();
799 os << setfill('-') << setw(30) << "-" << endl << setfill(' ');
800
801 for(std::size_t loc=0; loc<O.length(); ++loc)
802 {
803 for(std::size_t s1=0; s1<O.locBasis(loc).size(); ++s1)
804 {
805 for(std::size_t s2=0; s2<O.locBasis(loc).size(); ++s2)
806 {
807 for(std::size_t k=0; k<O.opBasis(loc).size(); ++k)
808 {
809 // TODO: Angepasster Code für W=Biped
810 /*
811 if(O.W_at(loc)[s1][s2][k].nonZeros() > 0)
812 {
813 std::array<typename Symmetry::qType,3> qCheck = {O.locBasis(loc)[s2],O.opBasis(l)[k],O.locBasis(loc)[s1]};
814 if(!Symmetry::validate(qCheck))
815 {
816 continue;
817 }
818 os << "[l=" << l << "]\t|" << Sym::format<Symmetry>(O.locBasis(loc)[s1]) << "><" << Sym::format<Symmetry>(O.locBasis(loc)[s2]) << "|:" << std::endl;
819 os << Matrix<Scalar,Dynamic,Dynamic>(O.W_at(loc)[s1][s2][k]) << std::endl;
820 }
821 */
822
823 }
824 }
825 }
826 os << setfill('-') << setw(80) << "-" << setfill(' ');
827 if(loc != O.length()-1)
828 {
829 os << std::endl;
830 }
831 }
832 return os;
833}
834
835template<typename Symmetry, typename Scalar1, typename Scalar2>
837{
838 lout << setfill('-') << setw(30) << "-" << setfill(' ');
839 lout << "Mpo: L=" << O1.length();
840 lout << setfill('-') << setw(30) << "-" << endl << setfill(' ');
841
842 for (std::size_t loc=0; loc<O1.length(); ++loc)
843 {
844 for (std::size_t s1=0; s1<O1.locBasis(loc).size(); ++s1)
845 {
846 for (std::size_t s2=0; s2<O1.locBasis(loc).size(); ++s2)
847 {
848 for (std::size_t k=0; k<O1.opBasis(loc).size(); ++k)
849 {
850 // TODO: Angepasster Code für W=Biped
851 /*
852 lout << "[l=" << loc << "]\t|" << Sym::format<Symmetry>(O1.locBasis(loc)[s1]) << "><" << Sym::format<Symmetry>(O1.locBasis(loc)[s2]) << "|:" << std::endl;
853 auto M1 = Matrix<Scalar1,Dynamic,Dynamic>(O1.W_at(loc)[s1][s2][k]);
854 auto Mtmp = Matrix<Scalar2,Dynamic,Dynamic>(O2.W_at(loc)[s1][s2][k]);
855 auto M2 = Mtmp.template cast<Scalar1>();
856 lout << "norm(diff)=" << (M1-M2).norm() << endl;
857 if((M1-M2).norm() > ::mynumeric_limits<Scalar1>::epsilon() or (M1-M2).norm() > ::mynumeric_limits<Scalar2>::epsilon())
858 {
859 lout << "M1=" << endl << M1 << endl << endl;
860 lout << "M2=" << endl << Matrix<Scalar2,Dynamic,Dynamic>(O2.W_at(loc)[s1][s2][k]) << endl << std::endl;
861 }
862 */
863 }
864 }
865 }
866 lout << setfill('-') << setw(80) << "-" << setfill(' ');
867 if(loc != O1.length()-1)
868 {
869 lout << std::endl;
870 }
871 }
872 lout << std::endl;
873}
874
875template<typename Symmetry, typename Scalar>
878{
879 Mpo<Symmetry,Scalar> output(input.size(), input.get_qTot(), input.label, false, false, input.get_boundary_condition(), input.get_verbosity());
880 output.reconstruct(input.get_O(), input.get_qAux(), input.get_qPhys(), input.is_finalized(), input.get_boundary_condition(), input.get_qTot());
881 return output;
882}
883
884#endif
Scalar localSumTrivial(int i)
Definition: DmrgExternal.h:275
BC
Definition: DmrgTypedefs.h:161
@ OPEN
Definition: DmrgTypedefs.h:163
void compare(const Mpo< Symmetry, Scalar1 > &O1, const Mpo< Symmetry, Scalar2 > &O2)
Definition: Mpo.h:836
Eigen::SparseMatrix< double, Eigen::ColMajor, EIGEN_DEFAULT_SPARSE_INDEX_TYPE > SparseMatrixXd
Definition: Mpo.h:17
std::ostream & operator<<(std::ostream &os, const Mpo< Symmetry, Scalar > &O)
Definition: Mpo.h:795
const std::vector< Qbasis< Symmetry > > & get_qAux_power(std::size_t power) const
Definition: MpoTerms.h:2812
std::vector< std::vector< std::string > > info
Definition: MpoTerms.h:86
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
std::size_t N_sites
Definition: MpoTerms.h:395
std::vector< std::string > get_info() const
Definition: MpoTerms.h:2755
std::vector< std::map< std::array< qType, 2 >, std::vector< std::vector< std::map< qType, OperatorType > > > > > O
Definition: MpoTerms.h:47
std::vector< std::vector< qType > > qPhys
Definition: MpoTerms.h:366
bool check_qPhys() const
Definition: MpoTerms.h:2843
const std::vector< std::map< std::array< qType, 2 >, std::vector< std::vector< std::map< qType, OperatorType > > > > > & get_O() const
Definition: MpoTerms.h:526
virtual void push(const std::size_t loc, const std::vector< OperatorType > &opList, const std::vector< qType > &qList, const Scalar lambda=1.0)
Definition: MpoTerms.h:893
std::vector< Qbasis< Symmetry > > qAux
Definition: MpoTerms.h:360
double memory(MEMUNIT memunit=kB) const
Definition: MpoTerms.h:4217
std::size_t size() const
Definition: MpoTerms.h:610
bool GOT_W
Definition: MpoTerms.h:127
DMRG::VERBOSITY::OPTION get_verbosity() const
Definition: MpoTerms.h:693
std::vector< std::vector< TwoSiteData< Symmetry, Scalar > > > calc_TwoSiteData() const
Definition: MpoTerms.h:3196
const std::vector< std::vector< qType > > & locBasis() const
Definition: MpoTerms.h:702
const std::string get_name() const
Definition: MpoTerms.h:476
BC get_boundary_condition() const
Definition: MpoTerms.h:574
const qType & get_qTot() const
Definition: MpoTerms.h:584
bool is_finalized() const
Definition: MpoTerms.h:521
void set_name(const std::string &label_in)
Definition: MpoTerms.h:471
const std::vector< std::vector< qType > > & get_qPhys() const
Definition: MpoTerms.h:546
std::size_t maxPower() const
Definition: MpoTerms.h:605
bool GOT_QAUX
Definition: MpoTerms.h:122
const std::vector< Qbasis< Symmetry > > & get_qAux() const
Definition: MpoTerms.h:536
reversedData reversed
Definition: MpoTerms.h:409
void save_label(const std::size_t loc, const std::string &info_label)
Definition: MpoTerms.h:2745
const std::vector< std::vector< qType > > & opBasis() const
Definition: MpoTerms.h:710
std::string label
Definition: MpoTerms.h:615
double sparsity(const std::size_t power=1, const bool PER_MATRIX=false) const
Definition: MpoTerms.h:4342
Definition: Mpo.h:40
bool GOT_SEMIOPEN_LEFT
Definition: Mpo.h:161
void setLocal(std::size_t loc, const OperatorType &op)
int get_dAux_max(int power=1) const
Symmetry::qType qType
Definition: Mpo.h:45
void setLocal(const std::vector< std::size_t > &locs, const std::vector< OperatorType > &ops)
Mpo()
Definition: Mpo.h:59
static std::vector< T > get_N_site_interaction(T const &Op0, Operator const &... Ops)
Definition: Mpo.h:117
boost::multi_array< Scalar, 4 > H2site(std::size_t loc, bool HALF_THE_LOCAL_TERM=false) const
Definition: Mpo.h:145
Mpo< Symmetry, Scalar > Identity() const
Definition: Mpo.h:130
Mps< Symmetry, double > StateXd
Definition: Mpo.h:147
bool IS_HAMILTONIAN() const
Definition: Mpo.h:139
void push_nextn(const std::size_t loc, const Scalar lambda, const OperatorType &op1, const OperatorType &trans, const OperatorType &op2)
Mpo< Symmetry > Operator
Definition: Mpo.h:153
Mpo(size_t L_input)
void setLocal(std::size_t loc, const OperatorType &op, const std::vector< OperatorType > &signOp)
MpsCompressor< Symmetry, double, double > CompressorXd
Definition: Mpo.h:151
static constexpr size_t Nq
Definition: Mpo.h:44
OperatorType LocalOp
Definition: Mpo.h:164
static Mpo< Symmetry, Scalar > cast_Terms_to_Mpo(const MpoTerms< Symmetry, Scalar > &input)
void push_width(const std::size_t n, const std::size_t loc, const Scalar lambda, const OperatorType &outOp, const std::vector< OperatorType > &trans, const OperatorType &inOp)
static Mpo< Symmetry, Scalar > Identity(const std::vector< std::vector< qType > > &qPhys, const qType &Q=Symmetry::qvacuum())
Mpo(MpoTerms< Symmetry, Scalar > &Terms)
Definition: Mpo.h:61
bool GOT_TWO_SITE_DATA
Definition: Mpo.h:155
void set_locality(std::size_t LocalSite_input)
Definition: Mpo.h:125
SiteOperator< Symmetry, Scalar > OperatorType
Definition: Mpo.h:43
void setLocalSum(const OperatorType &op, Scalar(*f)(int)=localSumTrivial)
Mpo< Symmetry, Scalar > Zero() const
Definition: Mpo.h:133
void set_localOperator(OperatorType LocalOp_input)
Definition: Mpo.h:127
std::size_t volume() const
Definition: Mpo.h:114
std::size_t length() const
Definition: Mpo.h:112
void setLocalSum(const std::vector< OperatorType > &op, std::vector< Scalar > coeffs)
void push_local(const std::size_t loc, const Scalar lambda, const OperatorType &op)
Scalar Scalar_
Definition: Mpo.h:57
int locality() const
Definition: Mpo.h:124
void precalc_TwoSiteData(bool FORCE=false)
bool GOT_SEMIOPEN_RIGHT
Definition: Mpo.h:162
void generate_label(std::size_t Lcell)
void setLocal(const std::vector< std::size_t > &locs, const std::vector< OperatorType > &ops, const std::vector< OperatorType > &signOps)
static Mpo< Symmetry, Scalar > Zero(const std::vector< std::vector< qType > > &qPhys)
MpsCompressor< Symmetry, std::complex< double >, double > CompressorXcd
Definition: Mpo.h:152
SparseMatrixXd SparseMatrixType
Definition: Mpo.h:42
bool IS_UNITARY() const
Definition: Mpo.h:135
std::string info(bool REDUCED=false) const
friend class Mpo
Definition: Mpo.h:52
Mpo(std::size_t L_input, qType Qtot_input, std::string label_input="Mpo", bool HERMITIAN_input=false, bool UNITARY_input=false, BC BC_input=BC::OPEN, DMRG::VERBOSITY::OPTION VERB_input=DMRG::VERBOSITY::OPTION::SILENT)
void push_qpath(const std::size_t loc, const std::vector< OperatorType > &opList, const std::vector< qType > &qList, const Scalar lambda=1.0)
Definition: Mpo.h:75
void push_tight(const std::size_t loc, const Scalar lambda, const OperatorType &op1, const OperatorType &op2)
Mps< Symmetry, std::complex< double > > StateXcd
Definition: Mpo.h:149
bool HAMILTONIAN
Definition: Mpo.h:160
bool HAS_TWO_SITE_DATA() const
Definition: Mpo.h:141
Umps< Symmetry, double > StateUd
Definition: Mpo.h:148
Symmetry::qType qVac
Definition: Mpo.h:46
bool UNITARY
Definition: Mpo.h:158
void setLocal(const std::vector< std::size_t > &locs, const std::vector< OperatorType > &ops, const OperatorType &signOp)
void setProductSum(const OperatorType &op1, const OperatorType &op2)
OperatorType localOperator() const
Definition: Mpo.h:126
std::vector< std::vector< TwoSiteData< Symmetry, Scalar > > > TSD
Definition: Mpo.h:156
bool HERMITIAN
Definition: Mpo.h:159
bool IS_HERMITIAN() const
Definition: Mpo.h:137
Matrix< Scalar, Dynamic, Dynamic > MatrixType
Definition: Mpo.h:56
void construct_from_pushlist(const PushType< OperatorType, CouplScalar > &pushlist, const std::vector< std::vector< std::string > > &labellist, size_t Lcell)
void calc_reversedData_from_pushlist(const PushType< OperatorType, CouplScalar > &pushlist, double tolerance=::mynumeric_limits< double >::epsilon())
void setLocalStag(std::size_t loc, const OperatorType &op, const std::vector< OperatorType > &stagSignOps)
bool HAS_W() const
Definition: Mpo.h:143
int LocalSite
Definition: Mpo.h:165
Umps< Symmetry, std::complex< double > > StateUcd
Definition: Mpo.h:150
void setLocal(std::size_t loc, const OperatorType &op, const OperatorType &signOp)
Definition: Mps.h:39
Definition: Umps.h:42
const bool COMPRESS
Definition: DmrgTypedefs.h:499
std::vector< Qbasis< Symmetry > > qAux
Definition: MpoTerms.h:405
std::vector< std::vector< std::vector< std::vector< Biped< Symmetry, MatrixType > > > > > W
Definition: MpoTerms.h:404
std::size_t size() const
Definition: DmrgTypedefs.h:212
void push_back(const std::tuple< std::size_t, std::vector< OtherOperator >, Scalar > &elem)
Definition: DmrgTypedefs.h:193