VMPS++
Loading...
Searching...
No Matches
HeisenbergU1Complex.h
Go to the documentation of this file.
1#ifndef STRAWBERRY_HeisenbergU1Complex
2#define STRAWBERRY_HeisenbergU1Complex
3
4#include "symmetry/U1.h"
7#include "ParamReturner.h"
8#include "Geometry2D.h" // from TOOLS
9
10namespace VMPS
11{
12class HeisenbergU1Complex : public Mpo<Sym::U1<Sym::SpinU1>,complex<double> >, public HeisenbergObservables<Sym::U1<Sym::SpinU1>,complex<double> >, public ParamReturner
13{
14public:
15
18 typedef complex<double> Scalar;
19
20 static qarray<1> singlet (int N=0) {return qarray<1>{0};};
21 static constexpr MODEL_FAMILY FAMILY = HEISENBERG;
22
23public:
26
27public:
28
31
32 HeisenbergU1Complex(Mpo<Symmetry,complex<double> > &Mpo_input, const vector<Param> &params)
33 :Mpo<Symmetry,complex<double> >(Mpo_input),
36 {
37 ParamHandler P(params,HeisenbergU1Complex::defaults);
38 size_t Lcell = P.size();
39 N_phys = 0;
40 for (size_t l=0; l<N_sites; ++l) N_phys += P.get<size_t>("Ly",l%Lcell);
41 this->precalc_TwoSiteData();
42 this->HERMITIAN = true;
43 this->HAMILTONIAN = true;
44 };
45
46 HeisenbergU1Complex (const size_t &L, const BC &boundary=BC::OPEN, const DMRG::VERBOSITY::OPTION &VERB=DMRG::VERBOSITY::OPTION::ON_EXIT);
47
48 HeisenbergU1Complex (const size_t &L, const vector<Param> &params, const BC & boundary=BC::OPEN, const DMRG::VERBOSITY::OPTION &VERB=DMRG::VERBOSITY::OPTION::ON_EXIT);
49
50 template<typename Symmetry_>
51 static void set_operators (const std::vector<SpinBase<Symmetry_> > &B, const ParamHandler &P,
52 PushType<SiteOperator<Symmetry_,complex<double> >,complex<double> >& pushlist,
53 std::vector<std::vector<std::string>>& labellist,
54 const BC boundary=BC::OPEN);
56
61 bool validate (qarray<1> qnum) const;
62
63 static const std::map<string,std::any> defaults;
64 static const std::map<string,std::any> sweep_defaults;
65};
66
67const std::map<string,std::any> HeisenbergU1Complex::defaults =
68{
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.},
73 {"R",0.},
74 {"Dy",0.}, {"Dyprime",0.}, {"Dyrung",0.},
75 {"Bz",0.}, {"Kz",0.},
76 {"mu",0.}, {"nu",0.}, // couple to Sz_i-1/2 and Sz_i+1/2
77 {"D",2ul}, {"maxPower",2ul}, {"CYLINDER",false}, {"Ly",1ul}, {"mfactor",1}
78};
79
80const std::map<string,std::any> HeisenbergU1Complex::sweep_defaults =
81{
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},
88 {"savePeriod",0ul}, {"CALC_S_ON_EXIT", true}, {"CONVTEST", DMRG::CONVTEST::VAR_2SITE}
89};
90
92HeisenbergU1Complex (const size_t &L, const BC &boundary, const DMRG::VERBOSITY::OPTION &VERB)
93:Mpo<Symmetry,complex<double> > (L, qarray<Symmetry::Nq>({0}), "", PROP::HERMITIAN, PROP::NON_UNITARY, boundary,VERB),
96{}
97
99HeisenbergU1Complex (const size_t &L, const vector<Param> &params, const BC &boundary, const DMRG::VERBOSITY::OPTION &VERB)
100:Mpo<Symmetry,complex<double> > (L, qarray<Symmetry::Nq>({0}), "", PROP::HERMITIAN, PROP::NON_UNITARY, boundary, VERB),
103{
104 ParamHandler P(params,defaults);
105 size_t Lcell = P.size();
106
107 for (size_t l=0; l<N_sites; ++l)
108 {
109 N_phys += P.get<size_t>("Ly",l%Lcell);
110 setLocBasis(B[l].get_basis().qloc(),l);
111 }
112
113 if (P.HAS_ANY_OF({"Jxy", "Jxypara", "Jxyperp", "Jxyfull"}))
114 {
115 this->set_name("XXZcomplex");
116 }
117 else if (P.HAS_ANY_OF({"Jz", "Jzpara", "Jzperp", "Jzfull"}))
118 {
119 this->set_name("IsingComplex");
120 }
121 else
122 {
123 this->set_name("HeisenbergComplex");
124 }
125
126 PushType<SiteOperator<Symmetry,complex<double> >,complex<double> > pushlist;
127 std::vector<std::vector<std::string>> labellist;
128 HeisenbergU1Complex::set_operators(B, P, pushlist, labellist, boundary);
129
130 this->construct_from_pushlist(pushlist, labellist, Lcell);
131 this->finalize(PROP::COMPRESS, P.get<size_t>("maxPower"));
132
133 this->precalc_TwoSiteData();
134}
135
136template<typename Symmetry_>
138set_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)
139{
140 std::size_t Lcell = P.size();
141 std::size_t N_sites = B.size();
142
143 if (labellist.size() != N_sites) {labellist.resize(N_sites);}
144
145 for (std::size_t loc=0; loc<N_sites; ++loc)
146 {
147 size_t lp1 = (loc+1)%N_sites;
148 size_t lp2 = (loc+2)%N_sites;
149 size_t lp3 = (loc+3)%N_sites;
150 size_t lp4 = (loc+4)%N_sites;
151
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();
157
158 stringstream ss1, ss2;
159 ss1 << "S=" << print_frac_nice(frac(P.get<size_t>("D",loc%Lcell)-1,2));
160 ss2 << "Ly=" << P.get<size_t>("Ly",loc%Lcell);
161 labellist[loc].push_back(ss1.str());
162 labellist[loc].push_back(ss2.str());
163
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"));
167
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);
172
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);
177
178 labellist[loc].push_back(Jxyperp.label);
179 labellist[loc].push_back(Jzperp.label);
180 labellist[loc].push_back(Jperp.label);
181
182 auto sum_array = [] (const ArrayXXd& a1, const ArrayXXd& a2)
183 {
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)
187 {
188 res(i,j) = a1(i,j) + a2(i,j);
189 }
190 return res;
191 };
192
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));
196
197// // Nearest-neighbour terms: Jxy, Jz, J
198// param2d Jxypara = P.fill_array2d<complex<double> >("Jxy", "Jxypara", {orbitals, next_orbitals}, loc%Lcell);
199// param2d Jzpara = P.fill_array2d<complex<double> >("Jz", "Jzpara", {orbitals, next_orbitals}, loc%Lcell);
200// labellist[loc].push_back(Jxypara.label);
201// labellist[loc].push_back(Jzpara.label);
202//
203// if (loc < N_sites-1 or !static_cast<bool>(boundary))
204// {
205// for (std::size_t alfa=0; alfa < orbitals; ++alfa)
206// for (std::size_t beta=0; beta < next_orbitals; ++beta)
207// {
208// pushlist.push_back(std::make_tuple(loc, Mpo<Symmetry,complex<double> >::get_N_site_interaction(B[loc].Scomp(SP,alfa),
209// B[lp1].Scomp(SM,beta)),
210// 0.5*Jxypara(alfa,beta)));
211// pushlist.push_back(std::make_tuple(loc, Mpo<Symmetry,complex<double> >::get_N_site_interaction(B[loc].Scomp(SM,alfa),
212// B[lp1].Scomp(SP,beta)),
213// 0.5*Jxypara(alfa,beta)));
214//
215// pushlist.push_back(std::make_tuple(loc, Mpo<Symmetry,complex<double> >::get_N_site_interaction(B[loc].Scomp(SZ,alfa),
216// B[lp1].Scomp(SZ,beta)),
217// Jzpara(alfa,beta)));
218// }
219// }
220
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
227 {
228 ArrayXXcd Full = P.get<Eigen::ArrayXXcd>(xxxFull);
229 vector<vector<std::pair<size_t,complex<double> > > > R = Geometry2D::rangeFormat(Full);
230
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!");}
233
234 for (size_t j=0; j<first.size(); j++)
235 for (size_t h=0; h<R[loc].size(); ++h)
236 {
237 size_t range = R[loc][h].first;
238 complex<double> value = R[loc][h].second;
239
240 if (range != 0)
241 {
242 vector<SiteOperatorQ<Symmetry_,Eigen::MatrixXcd> > ops(range+1);
243 ops[0] = first[j];
244 for (size_t i=1; i<range; ++i)
245 {
246 ops[i] = B[(loc+i)%N_sites].Id().template cast<complex<double> >();
247 }
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));
252 }
253 }
254
255 stringstream ss;
256 ss << label << "(" << Geometry2D::hoppingInfo(Full) << ")";
257 labellist[loc].push_back(ss.str());
258 };
259
260 // Full J-matrices
261 if (P.HAS("Jzfull"))
262 {
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> >();}
266
267 vector<vector<SiteOperatorQ<Symmetry_,Eigen::MatrixXcd> > > last {Sz_ranges};
268 push_full("Jzfull", "Jzᵢⱼ", first, last, {1.0}, {false});
269 }
270 if (P.HAS("Jxyfull"))
271 {
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> >();}
278
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});
281 }
282 if (P.HAS("JzfullA"))
283 {
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> >();}
287
288 vector<vector<SiteOperatorQ<Symmetry_,Eigen::MatrixXcd> > > last {Sz_ranges};
289 push_full("JzfullA", "JzAᵢⱼ", first, last, {1.0}, {false});
290 }
291 if (P.HAS("JxyfullA"))
292 {
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> >();}
299
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});
302 }
303
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);
306
307 labellist[loc].push_back(Jxy3site.label);
308 labellist[loc].push_back(Jxy4site.label);
309
310 if (loc < N_sites-3 or !static_cast<bool>(boundary))
311 {
312 assert(orbitals == next3_orbitals);
313 for (std::size_t alfa=0; alfa<orbitals; ++alfa)
314 {
315// if (std::abs(complex<double>(+0.25*Jxy3site(alfa,alfa),0.)) > 1e-10)
316// {
317// cout << "loc=" << loc << ", alfa=" << alfa << ", 0.25*val3=" << complex<double>(+0.25*Jxy3site(alfa,alfa),0.) << endl;
319// cout << endl;
320// }
321 // +-+-, sign:+ i,A-i,B--i+1,A-i+1,B
322 pushlist.push_back(std::make_tuple(loc,Mpo<Symmetry,complex<double> >::get_N_site_interaction(B[loc].Sp(alfa).template cast<complex<double> >(),
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.)
327 ));
328 // -+-+, sign:+ i,A-i,B--i+1,A-i+1,B
329 pushlist.push_back(std::make_tuple(loc,Mpo<Symmetry,complex<double> >::get_N_site_interaction(B[loc].Sm(alfa).template cast<complex<double> >(),
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.)
334 ));
335 // +--+, sign:- i,A-i,B--i+1,A-i+1,B
336 pushlist.push_back(std::make_tuple(loc,Mpo<Symmetry,complex<double> >::get_N_site_interaction(B[loc].Sp(alfa).template cast<complex<double> >(),
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.)
341 ));
342 // -++-, sign:- i,A-i,B--i+1,A-i+1,B
343 pushlist.push_back(std::make_tuple(loc,Mpo<Symmetry,complex<double> >::get_N_site_interaction(B[loc].Sm(alfa).template cast<complex<double> >(),
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.)
348 ));
349 }
350 }
351
352 if (loc < N_sites-4 or !static_cast<bool>(boundary))
353 {
354 assert(orbitals == next4_orbitals);
355 for (std::size_t alfa=0; alfa<orbitals; ++alfa)
356 {
357 //cout << "loc=" << loc << ", alfa=" << alfa << ", 0.25*val4=" << complex<double>(+0.25*Jxy4site(alfa,alfa)) << endl;
358 pushlist.push_back(std::make_tuple(lp1,Mpo<Symmetry,complex<double> >::get_N_site_interaction(B[lp1].Sp(alfa).template cast<complex<double> >(),
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.)
363 ));
364 pushlist.push_back(std::make_tuple(lp1,Mpo<Symmetry,complex<double> >::get_N_site_interaction(B[lp1].Sm(alfa).template cast<complex<double> >(),
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.)
369 ));
370 pushlist.push_back(std::make_tuple(lp1,Mpo<Symmetry,complex<double> >::get_N_site_interaction(B[lp1].Sp(alfa).template cast<complex<double> >(),
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.)
375 ));
376 pushlist.push_back(std::make_tuple(lp1,Mpo<Symmetry,complex<double> >::get_N_site_interaction(B[lp1].Sm(alfa).template cast<complex<double> >(),
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.)
381 ));
382 }
383 }
384 }
385}
386
387} //end namespace VMPS
388
389#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
BC
Definition: DmrgTypedefs.h:161
@ B
Definition: DmrgTypedefs.h:130
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
std::size_t N_phys
Definition: MpoTerms.h:400
std::size_t N_sites
Definition: MpoTerms.h:395
DMRG::VERBOSITY::OPTION VERB
Definition: MpoTerms.h:102
MpoTerms< Symmetry, OtherScalar > cast()
Definition: MpoTerms.h:3076
std::string label
Definition: MpoTerms.h:615
Definition: Mpo.h:40
static std::vector< T > get_N_site_interaction(T const &Op0, Operator const &... Ops)
Definition: Mpo.h:117
Definition: U1.h:25
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 > &params)
static const std::map< string, std::any > sweep_defaults
#define MAKE_TYPEDEFS(MODEL)
Definition: macros.h:4
const bool COMPRESS
Definition: DmrgTypedefs.h:499
const bool NON_UNITARY
Definition: DmrgTypedefs.h:495
const bool HERMITIAN
Definition: DmrgTypedefs.h:492
void finalize(bool PRINT_STATS=false)
Definition: functions.h:127
Definition: qarray.h:26