1#ifndef STRAWBERRY_HeisenbergU1Complex
2#define STRAWBERRY_HeisenbergU1Complex
38 size_t Lcell = P.size();
40 for (
size_t l=0; l<
N_sites; ++l)
N_phys += P.get<
size_t>(
"Ly",l%Lcell);
50 template<
typename Symmetry_>
53 std::vector<std::vector<std::string>>& labellist,
54 const BC boundary=BC::OPEN);
63 static const std::map<string,std::any>
defaults;
69 {
"J",0.}, {
"Jprime",0.}, {
"Jrung",0.},
70 {
"Jxy",0.}, {
"Jxyprime",0.}, {
"Jxyrung",0.},
71 {
"Jxy3site",0.}, {
"Jxy4site",0.},
72 {
"Jz",0.}, {
"Jzprime",0.}, {
"Jzrung",0.},
74 {
"Dy",0.}, {
"Dyprime",0.}, {
"Dyrung",0.},
77 {
"D",2ul}, {
"maxPower",2ul}, {
"CYLINDER",
false}, {
"Ly",1ul}, {
"mfactor",1}
82 {
"max_alpha",100.}, {
"min_alpha",1.e-11}, {
"lim_alpha",10ul}, {
"eps_svd",1.e-7},
83 {
"Mincr_abs", 50ul}, {
"Mincr_per", 2ul}, {
"Mincr_rel", 1.1},
84 {
"min_Nsv",0ul}, {
"max_Nrich",-1},
85 {
"max_halfsweeps",20ul}, {
"min_halfsweeps",4ul},
86 {
"Minit",1ul}, {
"Qinit",1ul}, {
"Mlimit",1000ul},
87 {
"tol_eigval",1e-7}, {
"tol_state",1e-6},
104 ParamHandler P(params,defaults);
105 size_t Lcell = P.size();
107 for (
size_t l=0; l<N_sites; ++l)
109 N_phys += P.get<
size_t>(
"Ly",l%Lcell);
110 setLocBasis(
B[l].get_basis().qloc(),l);
113 if (P.HAS_ANY_OF({
"Jxy",
"Jxypara",
"Jxyperp",
"Jxyfull"}))
115 this->set_name(
"XXZcomplex");
117 else if (P.HAS_ANY_OF({
"Jz",
"Jzpara",
"Jzperp",
"Jzfull"}))
119 this->set_name(
"IsingComplex");
123 this->set_name(
"HeisenbergComplex");
127 std::vector<std::vector<std::string>> labellist;
130 this->construct_from_pushlist(pushlist, labellist, Lcell);
133 this->precalc_TwoSiteData();
136template<
typename Symmetry_>
140 std::size_t Lcell = P.size();
145 for (std::size_t loc=0; loc<
N_sites; ++loc)
152 std::size_t orbitals =
B[loc].orbitals();
153 std::size_t next_orbitals =
B[lp1].orbitals();
154 std::size_t nextn_orbitals =
B[lp2].orbitals();
155 std::size_t next3_orbitals =
B[lp3].orbitals();
156 std::size_t next4_orbitals =
B[lp4].orbitals();
158 stringstream ss1, ss2;
160 ss2 <<
"Ly=" << P.get<
size_t>(
"Ly",loc%Lcell);
161 labellist[loc].push_back(ss1.str());
162 labellist[loc].push_back(ss2.str());
164 param2d Jxyperp = P.fill_array2d<
double>(
"Jxyrung",
"Jxy",
"Jxyperp", orbitals, loc%Lcell, P.get<
bool>(
"CYLINDER"));
165 param2d Jzperp = P.fill_array2d<
double>(
"Jzrung",
"Jz",
"Jzperp", orbitals, loc%Lcell, P.get<
bool>(
"CYLINDER"));
166 param2d Jperp = P.fill_array2d<
double>(
"Jrung",
"J",
"Jperp", orbitals, loc%Lcell, P.get<
bool>(
"CYLINDER"));
168 param1d Bz = P.fill_array1d<
double>(
"Bz",
"Bzorb", orbitals, loc%Lcell);
169 param1d mu = P.fill_array1d<
double>(
"mu",
"muorb", orbitals, loc%Lcell);
170 param1d nu = P.fill_array1d<
double>(
"nu",
"nuorb", orbitals, loc%Lcell);
171 param1d Kz = P.fill_array1d<
double>(
"Kz",
"Kzorb", orbitals, loc%Lcell);
173 labellist[loc].push_back(Bz.label);
174 labellist[loc].push_back(mu.label);
175 labellist[loc].push_back(nu.label);
176 labellist[loc].push_back(Kz.label);
178 labellist[loc].push_back(Jxyperp.label);
179 labellist[loc].push_back(Jzperp.label);
180 labellist[loc].push_back(Jperp.label);
182 auto sum_array = [] (
const ArrayXXd& a1,
const ArrayXXd& a2)
184 ArrayXXd res(a1.rows(), a1.cols());
185 for (
int i=0; i<a1.rows(); ++i)
186 for (
int j=0; j<a1.rows(); ++j)
188 res(i,j) = a1(i,j) + a2(i,j);
194 B[loc].HeisenbergHamiltonian(Jxyperp.a+Jperp.a, Jzperp.a+Jperp.a, Bz.a, mu.a, nu.a, Kz.a).template
cast<complex<double> >());
195 pushlist.push_back(std::make_tuple(loc, Hloc, 1.+0.i));
221 auto push_full = [&
N_sites, &loc, &
B, &P, &pushlist, &labellist, &boundary]
222 (
string xxxFull,
string label,
223 const vector<SiteOperatorQ<Symmetry_,Eigen::MatrixXcd> > &first,
224 const vector<vector<SiteOperatorQ<Symmetry_,Eigen::MatrixXcd> > > &last,
225 vector<complex<double> > factor,
226 vector<bool> CONJ) ->
void
228 ArrayXXcd Full = P.get<Eigen::ArrayXXcd>(xxxFull);
229 vector<vector<std::pair<size_t,complex<double> > > > R = Geometry2D::rangeFormat(Full);
231 if (
static_cast<bool>(boundary)) {assert(R.size() ==
N_sites and
"Use an (N_sites)x(N_sites) hopping matrix for open BC!");}
232 else {assert(R.size() >= 2*
N_sites and
"Use at least a (2*N_sites)x(N_sites) hopping matrix for infinite BC!");}
234 for (
size_t j=0; j<first.size(); j++)
235 for (
size_t h=0; h<R[loc].size(); ++h)
237 size_t range = R[loc][h].first;
238 complex<double> value = R[loc][h].second;
242 vector<SiteOperatorQ<Symmetry_,Eigen::MatrixXcd> > ops(range+1);
244 for (
size_t i=1; i<range; ++i)
246 ops[i] =
B[(loc+i)%
N_sites].Id().template cast<complex<double> >();
248 ops[range] = last[j][(loc+range)%
N_sites];
249 complex<double> total_value = factor[j] * value;
250 if (CONJ[j]) total_value = conj(total_value);
251 pushlist.push_back(std::make_tuple(loc, ops, total_value));
256 ss <<
label <<
"(" << Geometry2D::hoppingInfo(Full) <<
")";
257 labellist[loc].push_back(ss.str());
263 vector<SiteOperatorQ<Symmetry_,Eigen::MatrixXcd> > first {
B[loc].Sz(0).template cast<complex<double> >()};
264 vector<SiteOperatorQ<Symmetry_,Eigen::MatrixXcd> > Sz_ranges(
N_sites);
265 for (
size_t i=0; i<
N_sites; i++) {Sz_ranges[i] =
B[i].Sz(0).template cast<complex<double> >();}
267 vector<vector<SiteOperatorQ<Symmetry_,Eigen::MatrixXcd> > > last {Sz_ranges};
268 push_full(
"Jzfull",
"Jzᵢⱼ", first, last, {1.0}, {
false});
270 if (P.HAS(
"Jxyfull"))
272 vector<SiteOperatorQ<Symmetry_,Eigen::MatrixXcd> > first {
B[loc].Sp(0).template cast<complex<double> >(),
273 B[loc].
Sm(0).template cast<complex<double> >()};
274 vector<SiteOperatorQ<Symmetry_,Eigen::MatrixXcd> > Sp_ranges(
N_sites);
275 vector<SiteOperatorQ<Symmetry_,Eigen::MatrixXcd> > Sm_ranges(
N_sites);
276 for (
size_t i=0; i<
N_sites; i++) {Sp_ranges[i] =
B[i].Sp(0).template cast<complex<double> >();
277 Sm_ranges[i] =
B[i].Sm(0).template cast<complex<double> >();}
279 vector<vector<SiteOperatorQ<Symmetry_,Eigen::MatrixXcd> > > last {Sm_ranges, Sp_ranges};
280 push_full(
"Jxyfull",
"Jxyᵢⱼ", first, last, {0.5,0.5}, {
false,
true});
282 if (P.HAS(
"JzfullA"))
284 vector<SiteOperatorQ<Symmetry_,Eigen::MatrixXcd> > first {
B[loc].Sz(1).template cast<complex<double> >()};
285 vector<SiteOperatorQ<Symmetry_,Eigen::MatrixXcd> > Sz_ranges(
N_sites);
286 for (
size_t i=0; i<
N_sites; i++) {Sz_ranges[i] =
B[i].Sz(1).template cast<complex<double> >();}
288 vector<vector<SiteOperatorQ<Symmetry_,Eigen::MatrixXcd> > > last {Sz_ranges};
289 push_full(
"JzfullA",
"JzAᵢⱼ", first, last, {1.0}, {
false});
291 if (P.HAS(
"JxyfullA"))
293 vector<SiteOperatorQ<Symmetry_,Eigen::MatrixXcd> > first {
B[loc].Sp(1).template cast<complex<double> >(),
294 B[loc].
Sm(1).template cast<complex<double> >()};
295 vector<SiteOperatorQ<Symmetry_,Eigen::MatrixXcd> > Sp_ranges(
N_sites);
296 vector<SiteOperatorQ<Symmetry_,Eigen::MatrixXcd> > Sm_ranges(
N_sites);
297 for (
size_t i=0; i<
N_sites; i++) {Sp_ranges[i] =
B[i].Sp(1).template cast<complex<double> >();
298 Sm_ranges[i] =
B[i].Sm(1).template cast<complex<double> >();}
300 vector<vector<SiteOperatorQ<Symmetry_,Eigen::MatrixXcd> > > last {Sm_ranges, Sp_ranges};
301 push_full(
"JxyfullA",
"JxyAᵢⱼ", first, last, {0.5,0.5}, {
false,
true});
304 param2d Jxy3site = P.fill_array2d<
double>(
"Jxy3site",
"Jxy3site_array", {orbitals, next3_orbitals}, loc%Lcell);
305 param2d Jxy4site = P.fill_array2d<
double>(
"Jxy4site",
"Jxy4site_array", {orbitals, next4_orbitals}, loc%Lcell);
307 labellist[loc].push_back(Jxy3site.label);
308 labellist[loc].push_back(Jxy4site.label);
310 if (loc <
N_sites-3 or !
static_cast<bool>(boundary))
312 assert(orbitals == next3_orbitals);
313 for (std::size_t alfa=0; alfa<orbitals; ++alfa)
323 B[lp1].
Sm(alfa).
template cast<complex<double> >(),
324 B[lp2].
Sp(alfa).
template cast<complex<double> >(),
325 B[lp3].
Sm(alfa).
template cast<complex<double> >()),
326 complex<double>(+0.25*Jxy3site(alfa,alfa),0.)
330 B[lp1].
Sp(alfa).
template cast<complex<double> >(),
331 B[lp2].
Sm(alfa).
template cast<complex<double> >(),
332 B[lp3].
Sp(alfa).
template cast<complex<double> >()),
333 complex<double>(+0.25*Jxy3site(alfa,alfa),0.)
337 B[lp1].
Sm(alfa).
template cast<complex<double> >(),
338 B[lp2].
Sm(alfa).
template cast<complex<double> >(),
339 B[lp3].
Sp(alfa).
template cast<complex<double> >()),
340 complex<double>(-0.25*Jxy3site(alfa,alfa),0.)
344 B[lp1].
Sp(alfa).
template cast<complex<double> >(),
345 B[lp2].
Sp(alfa).
template cast<complex<double> >(),
346 B[lp3].
Sm(alfa).
template cast<complex<double> >()),
347 complex<double>(-0.25*Jxy3site(alfa,alfa),0.)
352 if (loc <
N_sites-4 or !
static_cast<bool>(boundary))
354 assert(orbitals == next4_orbitals);
355 for (std::size_t alfa=0; alfa<orbitals; ++alfa)
359 B[lp2].
Sm(alfa).
template cast<complex<double> >(),
360 B[lp3].
Sp(alfa).
template cast<complex<double> >(),
361 B[lp4].
Sm(alfa).
template cast<complex<double> >()),
362 complex<double>(+0.25*Jxy4site(alfa,alfa),0.)
365 B[lp2].
Sp(alfa).
template cast<complex<double> >(),
366 B[lp3].
Sm(alfa).
template cast<complex<double> >(),
367 B[lp4].
Sp(alfa).
template cast<complex<double> >()),
368 complex<double>(+0.25*Jxy4site(alfa,alfa),0.)
371 B[lp2].
Sm(alfa).
template cast<complex<double> >(),
372 B[lp3].
Sm(alfa).
template cast<complex<double> >(),
373 B[lp4].
Sp(alfa).
template cast<complex<double> >()),
374 complex<double>(-0.25*Jxy4site(alfa,alfa),0.)
377 B[lp2].
Sp(alfa).
template cast<complex<double> >(),
378 B[lp3].
Sp(alfa).
template cast<complex<double> >(),
379 B[lp4].
Sm(alfa).
template cast<complex<double> >()),
380 complex<double>(-0.25*Jxy4site(alfa,alfa),0.)
boost::rational< int > frac
std::string print_frac_nice(frac r)
vector< SpinBase< Sym::U1< Sym::SpinU1 > > > B
std::enable_if<!Dummy::IS_SPIN_SU2(), Mpo< Sym::U1< Sym::SpinU1 >, complex< double > > >::type Sm(size_t locx, size_t locy=0) const
std::enable_if<!Dummy::IS_SPIN_SU2(), Mpo< Sym::U1< Sym::SpinU1 >, complex< double > > >::type Sp(size_t locx, size_t locy=0) const
DMRG::VERBOSITY::OPTION VERB
MpoTerms< Symmetry, OtherScalar > cast()
static std::vector< T > get_N_site_interaction(T const &Op0, Operator const &... Ops)
void precalc_TwoSiteData(bool FORCE=false)
static const std::map< string, std::any > defaults
static constexpr MODEL_FAMILY FAMILY
static qarray< 1 > singlet(int N=0)
SiteOperator< Symmetry, SparseMatrix< complex< double > > > OperatorType
bool validate(qarray< 1 > qnum) const
static void set_operators(const std::vector< SpinBase< Symmetry_ > > &B, const ParamHandler &P, PushType< SiteOperator< Symmetry_, complex< double > >, complex< double > > &pushlist, std::vector< std::vector< std::string > > &labellist, const BC boundary=BC::OPEN)
Sym::U1< Sym::SpinU1 > Symmetry
HeisenbergU1Complex(Mpo< Symmetry, complex< double > > &Mpo_input, const vector< Param > ¶ms)
static const std::map< string, std::any > sweep_defaults
#define MAKE_TYPEDEFS(MODEL)
void finalize(bool PRINT_STATS=false)