1#ifndef STRAWBERRY_Mpo_WITH_Q
2#define STRAWBERRY_Mpo_WITH_Q
5#include "boost/multi_array.hpp"
8#include "termcolor.hpp"
11#include <Eigen/SparseCore>
14#ifndef EIGEN_DEFAULT_SPARSE_INDEX_TYPE
15#define EIGEN_DEFAULT_SPARSE_INDEX_TYPE int
17typedef Eigen::SparseMatrix<double,Eigen::ColMajor,EIGEN_DEFAULT_SPARSE_INDEX_TYPE>
SparseMatrixXd;
21#include <unsupported/Eigen/KroneckerProduct>
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;
38template<
typename Symmetry,
typename Scalar=
double>
44 static constexpr size_t Nq = Symmetry::Nq;
45 typedef typename Symmetry::qType
qType;
46 typename Symmetry::qType
qVac = Symmetry::qvacuum();
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;
69 template<
typename CouplScalar>
72 template<
typename CouplScalar>
75 void push_qpath (
const std::size_t loc,
const std::vector<OperatorType> &opList,
const std::vector<qType> &qList,
const Scalar lambda = 1.0)
86 void setLocal(
const std::vector<std::size_t>& locs,
const std::vector<OperatorType>& ops);
88 void setLocal(
const std::vector<std::size_t>& locs,
const std::vector<OperatorType>& ops,
const OperatorType& signOp);
90 void setLocal(
const std::vector<std::size_t>& locs,
const std::vector<OperatorType>& ops,
const std::vector<OperatorType>& signOps);
96 void setLocalSum(
const std::vector<OperatorType>& op, std::vector<Scalar> coeffs);
104 std::string
info (
bool REDUCED=
false)
const;
116 template<
typename T,
typename ...
Operator>
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.");}
156 std::vector<std::vector<TwoSiteData<Symmetry,Scalar>>>
TSD;
224template<
typename Symmetry,
typename Scalar>
226Mpo(std::size_t L_input)
227:
MpoTerms<Symmetry, Scalar>(L_input){}
229template<
typename Symmetry,
typename Scalar>
232:
MpoTerms<Symmetry, Scalar>(L_input, BC_input, Qtot_input, VERB_input),
HERMITIAN(HERMITIAN_input),
UNITARY(UNITARY_input)
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)
240 std::vector<OperatorType> oplist(0);
241 oplist.push_back(outOp);
242 for(std::size_t m=0; m<trans.size(); ++m)
244 oplist.push_back(trans[m]);
246 oplist.push_back(inOp);
247 this->
push(loc, oplist, lambda);
253 this->
push(loc, {op}, lambda);
259 this->
push(loc, {op1, op2}, lambda);
265 this->
push(loc, {op1, trans, op2}, lambda);
269template<
typename Symmetry,
typename Scalar>
271info (
bool REDUCED)
const
273 std::stringstream ss;
274 ss << termcolor::colorize << termcolor::bold;
283 ss << termcolor::reset <<
"→ L=" << this->
size();
285 ss <<
", " << Symmetry::name() <<
", ";
287 ss <<
"UNITARY=" << boolalpha <<
UNITARY <<
", ";
288 ss <<
"HERMITIAN=" << boolalpha <<
HERMITIAN <<
", ";
289 ss <<
"maxPower=" << this->
maxPower() <<
", ";
297 auto print_qaux = [
this,&ss] (
int power) ->
void
300 std::vector<int> dAux(this->
size()+1);
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)
307 dAux[loc+1] =
qAux[loc+1].fullM();
308 dAux_set.insert(std::make_pair(dAux[loc],dAux[loc+1]));
310 ss <<
"dAux" << power <<
"=";
311 for (
const auto& dAux_pair : dAux_set)
313 ss << dAux_pair.first <<
"x" << dAux_pair.second;
319 for (
int power=1; power<=this->
maxPower(); ++power) print_qaux(power);
321 ss <<
"mem=" << round(this->
memory(GB),3) <<
"GB";
330template<
typename Symmetry,
typename Scalar>
334 VectorXi res(this->
size());
336 for (std::size_t l=0; l<this->
size(); ++l)
338 res(l) =
qAux[l].fullM();
340 return res.maxCoeff();
400Identity(
const std::vector<std::vector<qType>>&
qPhys,
const typename Symmetry::qType& Q)
403 for(std::size_t loc=0; loc<
qPhys.size(); ++loc)
405 out.set_qPhys(loc,
qPhys[loc]);
412Zero(
const std::vector<std::vector<qType>>&
qPhys)
415 for(std::size_t loc=0; loc<
qPhys.size(); ++loc)
417 out.set_qPhys(loc,
qPhys[loc]);
426 std::stringstream ss;
430 std::map<std::string,std::set<std::size_t> > cells;
432 for (std::size_t loc=0; loc<
info.size(); ++loc)
434 cells[
info[loc]].insert(loc%Lcell);
437 if (cells.size() == 1)
439 ss <<
"(" <<
info[0] <<
")";
443 std::vector<std::pair<std::string,std::set<std::size_t> > > cells_resort(cells.begin(), cells.end());
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
449 return *min_element(a.second.begin(),a.second.end()) < *min_element(b.second.begin(),b.second.end());
453 for (
auto c:cells_resort)
455 ss << std::endl <<
" •l=";
461 if (c.second.size() == 1)
463 ss << *c.second.begin();
468 if (std::all_of(c.second.cbegin(), c.second.cend(), [](
int i){ return i%2==0;}) and c.second.size() == this->size()/2)
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)
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)
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)
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)
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)
495 if (c.second.size() == 2)
497 ss << *c.second.begin() <<
"," << *c.second.rbegin();
501 bool CONSECUTIVE =
true;
502 for (
auto it=c.second.begin(); it!=c.second.end(); ++it)
504 if (next(it) != c.second.end() and *next(it)!=*it+1ul)
511 ss << *c.second.begin() <<
"-" << *c.second.rbegin();
515 for (
auto it=c.second.begin(); it!=c.second.end(); ++it)
519 ss.seekp(-1,ios_base::end);
525 ss <<
": " << c.first;
532template<
typename Symmetry,
typename Scalar>
533template<
typename CouplScalar>
537 for(std::size_t i=0; i<pushlist.
size(); ++i)
539 const auto& [loc, ops, coupling] = pushlist[i];
540 if ( std::abs(coupling) != 0. )
542 this->
push(loc, ops, coupling);
545 for(std::size_t loc=0; loc<this->
size(); ++loc)
547 for(std::size_t i=0; i<labellist[loc].size(); ++i)
555template<
typename Symmetry,
typename Scalar>
556template<
typename CouplScalar>
561 for(std::size_t loc=0; loc<this->
N_sites; ++loc)
563 mpo_rev.set_qPhys(loc, this->
qPhys[this->N_sites-1-loc]);
565 for(std::size_t i=0; i<pushlist.
size(); ++i)
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)
573 ops_rev[j] = ops[range-1-j];
575 if ( std::abs(coupling) != 0. )
577 mpo_rev.push(loc_rev, ops_rev, coupling);
583 for(std::size_t loc=0; loc<this->
N_sites; ++loc)
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)
591 Wloc_rev[m].resize(Wloc[m].
size());
592 for(std::size_t n=0; n<Wloc[m].size(); ++n)
594 Wloc_rev[m][n].resize(Wloc[m][n].
size());
595 for(std::size_t t=0; t<Wloc[m][n].size(); ++t)
597 Wloc_rev[m][n][t] = Wloc[m][n][t].transpose();
611 assert(this->
check_qPhys() and
"Physical bases have to be set before");
616 pushlist.
push_back(std::make_tuple(loc, Ops, 1.));
617 std::vector<std::vector<std::string>> labellist(this->N_sites);
626 std::vector<OperatorType> signOps(loc, signOp);
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");
637 std::vector<OperatorType> ops = signOps;
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");
651 std::vector<OperatorType> ops(this->
size());
653 for(std::size_t loc2=0; loc2<this->
size(); ++loc2)
657 ops[loc2] = stagSignOps[loc2];
665setLocal(
const std::vector<std::size_t>& locs,
const std::vector<OperatorType>& ops)
676setLocal(
const std::vector<std::size_t>& locs,
const std::vector<OperatorType>& ops,
const std::vector<OperatorType>& signOps)
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!");
682 if(locs[left] > locs[right])
687 assert(locs[left] != locs[right] and
"setLocal(...) needs to local operators at different sites!");
688 assert(signOps.size() == locs[right]-locs[left]-1);
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)
694 ops_with_signs.push_back(signOps[pos]);
696 ops_with_signs.push_back(ops[right]);
697 this->
push(locs[left], ops_with_signs);
703setLocal(
const std::vector<std::size_t>& locs,
const std::vector<OperatorType>& ops,
const OperatorType& signOp)
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];
712 std::vector<OperatorType> signOps(distance-1, signOp);
719 assert(this->
check_qPhys() and
"Physical bases have to be set before");
720 for (std::size_t loc=0; loc<this->
size(); ++loc)
722 this->
push(loc, {f(loc)*op});
728setLocalSum(
const std::vector<OperatorType>& ops, std::vector<Scalar> coeffs)
730 assert(this->
check_qPhys() and
"Physical bases have to be set before");
732 for (std::size_t loc=0; loc<this->
size(); ++loc)
735 pushlist.
push_back(std::make_tuple(loc, Ops, 1.));
738 std::vector<std::vector<std::string>> labellist(this->N_sites);
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)
750 this->
push(loc, {op1, op2});
782template<
typename Symmetry,
typename Scalar>
790 GOT_TWO_SITE_DATA =
true;
794template<
typename Symmetry,
typename Scalar>
797 os << setfill(
'-') << setw(30) <<
"-" << setfill(
' ');
798 os <<
"Mpo: L=" <<
O.length();
799 os << setfill(
'-') << setw(30) <<
"-" << endl << setfill(
' ');
801 for(std::size_t loc=0; loc<
O.length(); ++loc)
803 for(std::size_t s1=0; s1<
O.locBasis(loc).
size(); ++s1)
805 for(std::size_t s2=0; s2<
O.locBasis(loc).
size(); ++s2)
807 for(std::size_t k=0; k<
O.opBasis(loc).
size(); ++k)
826 os << setfill(
'-') << setw(80) <<
"-" << setfill(
' ');
827 if(loc !=
O.length()-1)
835template<
typename Symmetry,
typename Scalar1,
typename Scalar2>
838 lout << setfill(
'-') << setw(30) <<
"-" << setfill(
' ');
839 lout <<
"Mpo: L=" << O1.
length();
840 lout << setfill(
'-') << setw(30) <<
"-" << endl << setfill(
' ');
842 for (std::size_t loc=0; loc<O1.
length(); ++loc)
844 for (std::size_t s1=0; s1<O1.
locBasis(loc).
size(); ++s1)
846 for (std::size_t s2=0; s2<O1.
locBasis(loc).
size(); ++s2)
848 for (std::size_t k=0; k<O1.
opBasis(loc).
size(); ++k)
866 lout << setfill(
'-') << setw(80) <<
"-" << setfill(
' ');
875template<
typename Symmetry,
typename Scalar>
Scalar localSumTrivial(int i)
void compare(const Mpo< Symmetry, Scalar1 > &O1, const Mpo< Symmetry, Scalar2 > &O2)
Eigen::SparseMatrix< double, Eigen::ColMajor, EIGEN_DEFAULT_SPARSE_INDEX_TYPE > SparseMatrixXd
std::ostream & operator<<(std::ostream &os, const Mpo< Symmetry, Scalar > &O)
const std::vector< Qbasis< Symmetry > > & get_qAux_power(std::size_t power) const
std::vector< std::vector< std::string > > info
void finalize(const bool COMPRESS=true, const std::size_t power=1, const double tolerance=::mynumeric_limits< double >::epsilon())
std::vector< std::string > get_info() const
std::vector< std::map< std::array< qType, 2 >, std::vector< std::vector< std::map< qType, OperatorType > > > > > O
std::vector< std::vector< qType > > qPhys
const std::vector< std::map< std::array< qType, 2 >, std::vector< std::vector< std::map< qType, OperatorType > > > > > & get_O() const
virtual void push(const std::size_t loc, const std::vector< OperatorType > &opList, const std::vector< qType > &qList, const Scalar lambda=1.0)
std::vector< Qbasis< Symmetry > > qAux
double memory(MEMUNIT memunit=kB) const
DMRG::VERBOSITY::OPTION get_verbosity() const
std::vector< std::vector< TwoSiteData< Symmetry, Scalar > > > calc_TwoSiteData() const
const std::vector< std::vector< qType > > & locBasis() const
const std::string get_name() const
BC get_boundary_condition() const
const qType & get_qTot() const
bool is_finalized() const
void set_name(const std::string &label_in)
const std::vector< std::vector< qType > > & get_qPhys() const
std::size_t maxPower() const
const std::vector< Qbasis< Symmetry > > & get_qAux() const
void save_label(const std::size_t loc, const std::string &info_label)
const std::vector< std::vector< qType > > & opBasis() const
double sparsity(const std::size_t power=1, const bool PER_MATRIX=false) const
void setLocal(std::size_t loc, const OperatorType &op)
int get_dAux_max(int power=1) const
void setLocal(const std::vector< std::size_t > &locs, const std::vector< OperatorType > &ops)
static std::vector< T > get_N_site_interaction(T const &Op0, Operator const &... Ops)
boost::multi_array< Scalar, 4 > H2site(std::size_t loc, bool HALF_THE_LOCAL_TERM=false) const
Mpo< Symmetry, Scalar > Identity() const
Mps< Symmetry, double > StateXd
bool IS_HAMILTONIAN() const
void push_nextn(const std::size_t loc, const Scalar lambda, const OperatorType &op1, const OperatorType &trans, const OperatorType &op2)
void setLocal(std::size_t loc, const OperatorType &op, const std::vector< OperatorType > &signOp)
MpsCompressor< Symmetry, double, double > CompressorXd
static constexpr size_t Nq
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)
void set_locality(std::size_t LocalSite_input)
SiteOperator< Symmetry, Scalar > OperatorType
void setLocalSum(const OperatorType &op, Scalar(*f)(int)=localSumTrivial)
Mpo< Symmetry, Scalar > Zero() const
void set_localOperator(OperatorType LocalOp_input)
std::size_t volume() const
std::size_t length() const
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)
void precalc_TwoSiteData(bool FORCE=false)
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
SparseMatrixXd SparseMatrixType
std::string info(bool REDUCED=false) const
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)
void push_tight(const std::size_t loc, const Scalar lambda, const OperatorType &op1, const OperatorType &op2)
Mps< Symmetry, std::complex< double > > StateXcd
bool HAS_TWO_SITE_DATA() const
Umps< Symmetry, double > StateUd
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
std::vector< std::vector< TwoSiteData< Symmetry, Scalar > > > TSD
bool IS_HERMITIAN() const
Matrix< Scalar, Dynamic, Dynamic > MatrixType
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)
Umps< Symmetry, std::complex< double > > StateUcd
void setLocal(std::size_t loc, const OperatorType &op, const OperatorType &signOp)
std::vector< Qbasis< Symmetry > > qAux
std::vector< std::vector< std::vector< std::vector< Biped< Symmetry, MatrixType > > > > > W
void push_back(const std::tuple< std::size_t, std::vector< OtherOperator >, Scalar > &elem)