VMPS++
Loading...
Searching...
No Matches
SpectralManager.h
Go to the documentation of this file.
1#ifndef SPECTRAL_MANAGER
2#define SPECTRAL_MANAGER
3
4#include "GreenPropagator.h"
5#include "DmrgLinearAlgebra.h"
6//#include "RootFinder.h" // from ALGS
7#include "VUMPS/VumpsSolver.h"
8
9template<typename Hamiltonian>
10class SpectralManager
11{
12public:
13
14 typedef typename Hamiltonian::Mpo::Scalar_ Scalar;
15 typedef typename Hamiltonian::Symmetry Symmetry;
16
17 SpectralManager(){};
18
19 // IBC
20 SpectralManager (const vector<string> &specs_input, const Hamiltonian &H, const vector<Param> &params, VUMPS::CONTROL::GLOB GlobSweepParams, qarray<Symmetry::Nq> Q,
21 int Ncells_input, const vector<Param> &params_hetero,
22 string gs_label="gs", bool LOAD_GS=false, bool SAVE_GS=false,
23 DMRG::VERBOSITY::OPTION VERB=DMRG::VERBOSITY::HALFSWEEPWISE,
24 double tol_OxV=2., int locyBra=0, int locyKet=0);
25
26 // OBC
27 SpectralManager (const vector<string> &specs_input, const Hamiltonian &H, const vector<Param> &params, DMRG::CONTROL::GLOB GlobSweepParams, qarray<Symmetry::Nq> Q,
28 string gs_label="gs", bool LOAD_GS=false, bool SAVE_GS=false,
29 DMRG::VERBOSITY::OPTION VERB=DMRG::VERBOSITY::HALFSWEEPWISE,
30 double tol_OxV=2., int locy=0);
31
32 SpectralManager (const vector<string> &specs_input, const Hamiltonian &H, DMRG::VERBOSITY::OPTION VERB=DMRG::VERBOSITY::HALFSWEEPWISE)
33 :specs(specs_input), Hwork(H), CHOSEN_VERB(VERB)
34 {};
35
36 template<typename HamiltonianThermal>
37 void beta_propagation (const Hamiltonian &Hprop, const HamiltonianThermal &Htherm, int Lcell, int dLphys,
38 double betamax_input, double dbeta_input, double tol_compr_beta_input, size_t Mlim, qarray<Symmetry::Nq> Q,
39 double s_betainit, double betaswitch,
40 string wd, string th_label, bool LOAD_BETA=false, bool SAVE_BETA=true,
41 DMRG::VERBOSITY::OPTION VERB=DMRG::VERBOSITY::HALFSWEEPWISE,
42 vector<double> stateSavePoints={}, vector<string> stateSaveLabels={},
43 int Ntaylor=0, bool CALC_C=true, bool CALC_CHI=true, bool USE_PHIT=false);
44
45 void continue_beta_propagation (const Hamiltonian &Hprop, int Lcell, int dLphys,
46 double s_betainit, double betainit, double betamax, double dbeta, double tol_compr_beta, size_t Mlim, qarray<Hamiltonian::Symmetry::Nq> Q,
47 double betaswitch,
48 string wd, string th_label, string LOAD_BETA, bool SAVE_BETA,
49 DMRG::VERBOSITY::OPTION VERB,
50 vector<double> stateSavePoints, vector<string> stateSaveLabels,
51 bool CALC_C, bool CALC_CHI);
52
53 void apply_operators_on_thermal_state (int Lcell, int dLphys, bool CHECK=true);
54
55 void compute (string wd, string label, int Ns, double tmax, double dt=0.2,
56 double wmin=-10., double wmax=10., int wpoints=501, Q_RANGE QR=ZERO_2PI, int qpoints=501, GREEN_INTEGRATION INT=OOURA,
57 size_t Mlim=500ul, double tol_DeltaS=1e-2, double tol_compr=1e-4);
58
59 void compute_finiteCell (int Lcell, int x0, string wd, string label, double tmax, double dt=0.1,
60 double wmin=-10., double wmax=10., int wpoints=501, Q_RANGE QR=ZERO_2PI, int qpoints=501, GREEN_INTEGRATION INT=OOURA,
61 size_t Mlim=500ul, double tol_DeltaS=1e-2, double tol_compr=1e-4);
62
63 void compute_finite (size_t j0, string wd, string label, int Ns, double tmax, double dt=0.1,
64 double wmin=-10., double wmax=10., int wpoints=501, GREEN_INTEGRATION INT=OOURA,
65 size_t Mlim=500ul, double tol_DeltaS=1e-2, double tol_compr=1e-4);
66
67 void compute_thermal (string wd, string label, int dLphys,
68 double tmax, double dt=0.1, double wmin=-10., double wmax=10., int wpoints=501, Q_RANGE QR=ZERO_2PI, int qpoints=501, GREEN_INTEGRATION INT=OOURA,
69 size_t Mlim=500ul, double tol_DeltaS=1e-2, double tol_compr=1e-4);
70
71 void reload (string wd, const vector<string> &specs_input, string label, int L, int Ncells, int Ns, double tmax,
72 double wmin=-10., double wmax=10., int wpoints=501, Q_RANGE QR=ZERO_2PI, int qpoints=501, GREEN_INTEGRATION INT=OOURA);
73
74 const Umps<Symmetry,Scalar> &ground() const {return g.state;};
75 const double &energy() const {return g.energy;};
76 const Mps<Symmetry,Scalar> &get_PhiT() const {return PhiT;};
77
78 void make_A1P (GreenPropagator<Hamiltonian,Symmetry,Scalar,complex<double> > &Gfull, string wd, string label, int Ns,
79 double tmax, double wmin=-10., double wmax=10., int wpoints=501, Q_RANGE QR=ZERO_2PI, int qpoints=501, GREEN_INTEGRATION INT=OOURA,
80 bool SAVE_N_MU=true);
81
82 void make_A1P_finite (GreenPropagator<Hamiltonian,Symmetry,Scalar,complex<double> > &Gfull, string wd, string label,
83 double tmax, double wmin=-10., double wmax=10., int wpoints=501, GREEN_INTEGRATION INT=OOURA);
84
85 Mpo<Symmetry,Scalar> get_Op (const Hamiltonian &H, size_t loc, std::string spec, double factor=1., size_t locy=0, int dLphys=1);
86
87 void set_measurement (int iz, string spec, double factor, int dLphys, qarray<Symmetry::Nq> Q, int Lcell, int measure_interval_input=10, string measure_name_input="M", string measure_subfolder_input=".", bool TRANSFORM=false);
88
89 void resize_Green (string wd, string label, int Ns, double tmax, double dt, double wmin, double wmax, int wpoints, Q_RANGE QR, int qpoints, GREEN_INTEGRATION INT);
90
91 void FTcell_xq();
92
93 static bool TIME_DIR (std::string spec)
94 {
95 // true=forwards in time
96 // false=backwards in time
97 return (spec=="PES" or spec=="PESUP" or spec=="PESDN" or spec=="AES" or spec=="IPSZ" or spec=="ICSF" or spec=="SDAGSF" or spec=="PDAGSF")? false:true;
98 }
99
100 static string DAG (std::string spec)
101 {
102 string res;
103 if (spec == "PES") res = "IPE";
104 if (spec == "PESUP") res = "IPEUP";
105 if (spec == "PESDN") res = "IPEDN";
106 else if (spec == "SSF") res = "SDAGSF";
107 else if (spec == "SPM") res = "SMP";
108 else if (spec == "SMP") res = "SPM";
109 else if (spec == "QSF") res = "QDAGSF";
110 else if (spec == "SSZ") res = "SSZ";
111 else if (spec == "IPE") res = "PES";
112 else if (spec == "IPEUP") res = "PESUP";
113 else if (spec == "IPEDN") res = "PESDN";
114 else if (spec == "AES") res = "APS";
115 else if (spec == "APS") res = "AES";
116 else if (spec == "CSF") res = "CSF";
117 else if (spec == "ICSF") res = "ICSF";
118 else if (spec == "PSZ") res = "PSZ";
119 else if (spec == "IPSZ") res = "IPSZ";
120 else if (spec == "PSF") res = "PDAGSF";
121 else if (spec == "PDAGSF")res = "PSF";
122 else if (spec == "HSF") res = "IHSF";
123 else if (spec == "IHSF") res = "HSF";
124 else if (spec == "IHS") res = "HSF";
125 else if (spec == "HTS") res = "IHTS";
126 else if (spec == "IHTS") res = "HTS";
127 return res;
128 }
129
130 static bool CHECK_SPEC (string spec)
131 {
132 std::array<string,26> possible_specs = {"PES","PESUP","PESDN", //3,3
133 "SSF","SSZ", //2,5
134 "IPE","IPEUP","IPEDN", //3,8
135 "AES","APS", //2,10
136 "CSF","ICSF","PSZ","IPSZ","PSF", //5,15
137 "HSF","IHSF", //2,17
138 "HTS","IHTS", // 2,19
139 "JCJC","JEJE", "JCJE", "JEJC", // 4,23
140 "QSF", // 1,24
141 "SPM","SMP" // 2,26
142 };
143 return find(possible_specs.begin(), possible_specs.end(), spec) != possible_specs.end();
144 }
145
146 vector<vector<MatrixXcd>> get_GwqCell (int z) const {return Green[z].get_GwqCell();};
147
148 void set_Ncells (int Ncells_input) {Ncells=Ncells_input;};
149
150 void set_PhiT (const Mps<Symmetry,Scalar> &PhiT_input)
151 {
152 PhiT = PhiT_input;
153 }
154
155private:
156
157 size_t L, Lhetero, Ncells;
158 size_t x0;
159 size_t Nspec;
160 vector<string> specs;
161
162 Hamiltonian Hwork;
163 Mps<Symmetry,Scalar> Phi;
164 double Eg;
165 vector<vector<Mps<Symmetry,complex<double>>>> OxPhiCellBra;
166 vector<vector<Mps<Symmetry,complex<double>>>> OxPhiCellKet;
167
168 Eigenstate<Umps<Symmetry,Scalar>> g;
169 Eigenstate<Mps<Symmetry,Scalar>> gfinite;
170 vector<GreenPropagator<Hamiltonian,typename Hamiltonian::Symmetry,Scalar,complex<double>>> Green;
171
172 double betamax, dbeta, tol_compr_beta;
173 Mps<Symmetry,Scalar> PhiT;
174 Mps<Symmetry,complex<double>> PhiTt;
175 vector<vector<Mps<Symmetry,complex<double>>>> OxPhiTt;
176 vector<vector<Mpo<typename Hamiltonian::Symmetry,Scalar>>> Odag;
177
178 DMRG::VERBOSITY::OPTION CHOSEN_VERB;
179 vector<vector<Scalar>> Oshift;
180};
181
182// non-thermal IBC
183template<typename Hamiltonian>
184SpectralManager<Hamiltonian>::
185SpectralManager (const vector<string> &specs_input, const Hamiltonian &H, const vector<Param> &params, VUMPS::CONTROL::GLOB GlobSweepParams, qarray<Symmetry::Nq> Q,
186 int Ncells_input, const vector<Param> &params_hetero,
187 string gs_label, bool LOAD_GS, bool SAVE_GS,
188 DMRG::VERBOSITY::OPTION VERB,
189 double tol_OxV, int locyBra, int locyKet)
190:specs(specs_input), Ncells(Ncells_input), CHOSEN_VERB(VERB)
191{
192 for (const auto &spec:specs)
193 {
194 assert(CHECK_SPEC(spec) and "Wrong spectral abbreviation!");
195 }
196 L = H.length();
197 Nspec = specs.size();
198 Lhetero = L*Ncells;
199
200 typename Hamiltonian::uSolver uDMRG(VERB);
201
202 if (LOAD_GS)
203 {
204 g.state.load(gs_label, g.energy);
205 lout << "loaded: " << g.state.info() << endl;
206 }
207 else
208 {
209 uDMRG.userSetGlobParam();
210 uDMRG.GlobParam = GlobSweepParams;
211 uDMRG.edgeState(H,g,Q);
212 if (SAVE_GS)
213 {
214 lout << "saving groundstate..." << endl;
215 g.state.save(gs_label, "groundstate", g.energy);
216 }
217 }
218
219 lout << setprecision(16) << "g.energy=" << g.energy << endl;
220
221 Hwork = Hamiltonian(Lhetero,params_hetero,BC::INFINITE);
222 lout << "H_hetero: " << Hwork.info() << endl;
223 Hwork.transform_base(Q,false,L); // PRINT=false
224 Hwork.precalc_TwoSiteData(true); // FORCE=true
225
226 // Phi
227 Phi = uDMRG.create_Mps(Ncells, g, H, Lhetero/2); // ground state as heterogenic MPS
228
229 // shift values O→O-<O>
230 vector<vector<Scalar>> OshiftKet(Nspec);
231 for (int z=0; z<Nspec; ++z)
232 {
233 OshiftKet[z].resize(L);
234 for (int l=0; l<L; ++l)
235 {
236 if (specs[z] == "HSF")
237 {
238 Hamiltonian Haux(2*L, {{"maxPower",1ul}}, BC::INFINITE, DMRG::VERBOSITY::SILENT);
239 OshiftKet[z][l] = avg(g.state, get_Op(Haux,l,specs[z],1.,locyKet), g.state);
240 }
241 else
242 {
243 OshiftKet[z][l] = avg(g.state, get_Op(H,l,specs[z],1.,locyKet), g.state);
244 }
245 lout << "spec=" << specs[z] << ", l=" << l << ", shift=" << OshiftKet[z][l] << endl;
246 }
247 }
248
249 // O and OxPhiCell for counterpropagate
250 vector<vector<Mpo<Symmetry,Scalar>>> Oket(Nspec);
251 for (int z=0; z<Nspec; ++z)
252 {
253 Oket[z].resize(L);
254 for (int l=0; l<L; ++l)
255 {
256 Oket[z][l] = get_Op(Hwork, Lhetero/2+l, specs[z], 1., locyKet);
257 if (abs(OshiftKet[z][l]) > 1e-6)
258 {
259 Oket[z][l].scale(1.,-OshiftKet[z][l]);
260 }
261 Oket[z][l].transform_base(Q,false,L); // PRINT=false
262 lout << Oket[z][l].info() << endl;
263 }
264 }
265
266 vector<vector<Mpo<Symmetry,Scalar>>> Obra(Nspec);
267 if (locyBra == locyKet)
268 {
269 Obra = Oket;
270 }
271 else
272 {
273 vector<vector<Scalar>> OshiftBra(Nspec);
274 for (int z=0; z<Nspec; ++z)
275 {
276 OshiftBra[z].resize(L);
277 for (int l=0; l<L; ++l)
278 {
279 if (specs[z] == "HSF")
280 {
281 Hamiltonian Haux(2*L, {{"maxPower",1ul}}, BC::INFINITE, DMRG::VERBOSITY::SILENT);
282 OshiftBra[z][l] = avg(g.state, get_Op(Haux,l,specs[z],1.,locyBra), g.state);
283 }
284 else
285 {
286 OshiftBra[z][l] = avg(g.state, get_Op(H,l,specs[z],1.,locyBra), g.state);
287 }
288 lout << "spec=" << specs[z] << ", l=" << l << ", shift=" << OshiftBra[z][l] << endl;
289 }
290 }
291
292 for (int z=0; z<Nspec; ++z)
293 {
294 Obra[z].resize(L);
295 for (int l=0; l<L; ++l)
296 {
297 Obra[z][l] = get_Op(Hwork, Lhetero/2+l, specs[z], 1., locyBra);
298 if (abs(OshiftBra[z][l]) > 1e-6)
299 {
300 Obra[z][l].scale(1.,-OshiftBra[z][l]);
301 }
302 Obra[z][l].transform_base(Q,false,L); // PRINT=false
303 lout << Obra[z][l].info() << endl;
304 }
305 }
306 }
307
308 OxPhiCellKet.resize(Nspec);
309 OxPhiCellBra.resize(Nspec);
310 for (int z=0; z<Nspec; ++z)
311 {
312 OxPhiCellKet[z].resize(L);
313 OxPhiCellBra[z].resize(L);
314
315 auto tmpKet = uDMRG.create_Mps(Ncells, g, H, Oket[z][0], Oket[z], tol_OxV); // O[z][0] for boundaries, O[z] is multiplied (vector in cell size)
316 auto tmpBra = uDMRG.create_Mps(Ncells, g, H, Obra[z][0], Obra[z], tol_OxV);
317
318 for (int l=0; l<L; ++l)
319 {
320 OxPhiCellKet[z][l] = tmpKet[l].template cast<complex<double>>();
321 OxPhiCellBra[z][l] = tmpBra[l].template cast<complex<double>>();
322 }
323 }
324
325 Eg = isReal(avg_hetero(Phi, Hwork, Phi, true)); // USE_BOUNDARY=true
326 lout << setprecision(16) << "Eg=" << Eg << ", eg=" << g.energy << endl;
327}
328
329// non-thermal OBC
330template<typename Hamiltonian>
331SpectralManager<Hamiltonian>::
332SpectralManager (const vector<string> &specs_input, const Hamiltonian &H, const vector<Param> &params, DMRG::CONTROL::GLOB GlobSweepParams, qarray<Symmetry::Nq> Q,
333 string gs_label, bool LOAD_GS, bool SAVE_GS, DMRG::VERBOSITY::OPTION VERB, double tol_OxV, int locy)
334:specs(specs_input), CHOSEN_VERB(VERB)
335{
336 for (const auto &spec:specs)
337 {
338 assert(CHECK_SPEC(spec) and "Wrong spectral abbreviation!");
339 }
340 L = H.length();
341 Ncells = 1;
342 Nspec = specs.size();
343 Lhetero = L*Ncells;
344
345 typename Hamiltonian::Solver fDMRG(VERB);
346
347 if (LOAD_GS)
348 {
349 gfinite.state.load(gs_label, gfinite.energy);
350 lout << "loaded: " << gfinite.state.info() << endl;
351 }
352 else
353 {
354 fDMRG.userSetGlobParam();
355 fDMRG.GlobParam = GlobSweepParams;
356 fDMRG.edgeState(H,gfinite,Q);
357 if (SAVE_GS)
358 {
359 lout << "saving groundstate..." << endl;
360 gfinite.state.save(gs_label, "groundstate", gfinite.energy);
361 }
362 }
363
364 lout << setprecision(16) << "gfinite.energy=" << gfinite.energy << endl;
365
366 Hwork = Hamiltonian(Lhetero,params,BC::OPEN);
367 lout << "H: " << Hwork.info() << endl;
368 Hwork.precalc_TwoSiteData(true); // FORCE=true
369
370 // Phi
371 Phi = gfinite.state;
372
373 // shift values O→O-<O>
374 Oshift.resize(Nspec);
375 for (int z=0; z<Nspec; ++z)
376 {
377 Oshift[z].resize(L);
378 for (int l=0; l<L; ++l)
379 {
380 if (specs[z] == "HSF")
381 {
382 Hamiltonian Haux(2*L, {{"maxPower",1ul}}, BC::OPEN, DMRG::VERBOSITY::SILENT);
383 Oshift[z][l] = avg(gfinite.state, get_Op(Haux,l,specs[z],1.,locy), gfinite.state);
384 }
385 else
386 {
387 Oshift[z][l] = avg(gfinite.state, get_Op(H,l,specs[z],1.,locy), gfinite.state);
388 }
389 lout << "spec=" << specs[z] << ", l=" << l << ", shift=" << Oshift[z][l] << endl;
390 }
391 }
392
393 // O and OxPhiCell
394 vector<vector<Mpo<Symmetry,Scalar>>> O(Nspec);
395 for (int z=0; z<Nspec; ++z)
396 {
397 O[z].resize(L);
398 for (int l=0; l<L; ++l)
399 {
400 O[z][l] = get_Op(Hwork, l, specs[z], 1., locy);
401 if (O[z][l].Qtarget() == Symmetry::qvacuum())
402 {
403 O[z][l].scale(1.,-Oshift[z][l]);
404 }
405 }
406 }
407
408 OxPhiCellKet.resize(Nspec);
409 for (int z=0; z<Nspec; ++z)
410 {
411 OxPhiCellKet[z].resize(L);
412 for (int l=0; l<L; ++l)
413 {
414 Mps<Symmetry,Scalar> tmp;
415 Mpo<Symmetry,Scalar> Otmp = O[z][l];
416 OxV_exact(Otmp, Phi, tmp, tol_OxV, DMRG::VERBOSITY::SILENT, 200, 1);
417 OxPhiCellKet[z][l] = tmp.template cast<complex<double>>();
418 }
419 }
420
421 Eg = gfinite.energy;
422}
423
424template<typename Hamiltonian>
425template<typename HamiltonianThermal>
426void SpectralManager<Hamiltonian>::
427beta_propagation (const Hamiltonian &Hprop, const HamiltonianThermal &Htherm, int Lcell, int dLphys,
428 double betamax, double dbeta, double tol_compr_beta, size_t Mlim, qarray<Hamiltonian::Symmetry::Nq> Q,
429 double s_betainit, double betaswitch,
430 string wd, string th_label, bool LOAD_BETA, bool SAVE_BETA,
431 DMRG::VERBOSITY::OPTION VERB,
432 vector<double> stateSavePoints, vector<string> stateSaveLabels,
433 int Ntaylor, bool CALC_C, bool CALC_CHI, bool USE_PHIT)
434{
435 for (const auto &spec:specs)
436 {
437 assert(CHECK_SPEC(spec) and "Wrong spectral abbreviation!");
438 }
439 L = Hprop.length()/dLphys;
440 Nspec = specs.size();
441
442 if (LOAD_BETA)
443 {
444 lout << "loading beta result from: " << make_string(wd,"/betaRes_",th_label) << endl;
445 PhiT.load(make_string(wd,"/betaRes_",th_label));
446 lout << "loaded: " << PhiT.info() << endl;
447 }
448 else
449 {
450 if (!USE_PHIT)
451 {
452 lout << "Computing T=infinity state" << endl;
453 typename HamiltonianThermal::Solver fDMRG(DMRG::VERBOSITY::ON_EXIT);
454 Eigenstate<Mps<Symmetry,typename HamiltonianThermal::Mpo::Scalar_> > th;
455
456 DMRG::CONTROL::GLOB GlobParam;
457 fDMRG.GlobParam.CALC_S_ON_EXIT = false;
458 DMRG::CONTROL::DYN DynParam;
459 if (dLphys == 2)
460 {
461 GlobParam.Minit = 10ul;
462 GlobParam.Qinit = 10ul;
463 GlobParam.INITDIR = DMRG::DIRECTION::RIGHT; // 1=left->right, 0=right->left
464
465 size_t start_2site = 0ul;
466 size_t end_2site = 10ul;
467 //int period_2site = 2ul;
468 DynParam.iteration = [start_2site,end_2site] (size_t i) {return (i>=start_2site and i<=end_2site)? DMRG::ITERATION::TWO_SITE : DMRG::ITERATION::ONE_SITE;};
469 DynParam.max_alpha_rsvd = [] (size_t i) {return (i<20ul)? 1e8:0.;};
470 fDMRG.DynParam = DynParam;
471 fDMRG.userSetDynParam();
472 }
473 fDMRG.GlobParam = GlobParam;
474 fDMRG.userSetGlobParam();
475 th.state.eps_truncWeight = 0.;
476 fDMRG.edgeState(Htherm, th, Q, LANCZOS::EDGE::GROUND, false);
477
478 //lout << Htherm.info() << endl;
479 th.state.entropy_skim();
480 lout << th.state.entropy().transpose() << endl;
481
482 bool ALL;
483
484 if (dLphys==2)
485 {
486 vector<bool> ENTROPY_CHECK1;
487 for (int l=0; l<dLphys*L-1; l+=dLphys) ENTROPY_CHECK1.push_back(abs(th.state.entropy()(l))>1e-10);
488
489 vector<bool> ENTROPY_CHECK2;
490 for (int l=1; l<dLphys*L-1; l+=dLphys) ENTROPY_CHECK2.push_back(abs(th.state.entropy()(l))<1e-10);
491
492 ALL = all_of(ENTROPY_CHECK1.begin(), ENTROPY_CHECK1.end(), [](const bool v){return v;}) and
493 all_of(ENTROPY_CHECK2.begin(), ENTROPY_CHECK2.end(), [](const bool v){return v;});
494
495 //ALL = true;
496 }
497 else
498 {
499 // vector<bool> ENTROPY_CHECK;
500 // for (int l=0; l<L-1; l+=1) ENTROPY_CHECK.push_back(abs(th.state.entropy()(l))<1e-10);
501 //
502 // ALL = all_of(ENTROPY_CHECK.begin(), ENTROPY_CHECK.end(), [](const bool v){return v;});
503
504 ALL = true;
505 }
506
507 while (ALL == false)
508 {
509 lout << termcolor::yellow << "restarting..." << termcolor::reset << endl;
510
511 typename HamiltonianThermal::Solver fDMRGnew(DMRG::VERBOSITY::ON_EXIT);
512
513 DMRG::CONTROL::GLOB GlobParam;
514 fDMRGnew.GlobParam.CALC_S_ON_EXIT = false;
515 if (GlobParam.INITDIR == DMRG::DIRECTION::RIGHT)
516 {
517 GlobParam.INITDIR = DMRG::DIRECTION::LEFT;
518 }
519 else
520 {
521 GlobParam.INITDIR = DMRG::DIRECTION::RIGHT;
522 }
523 fDMRGnew.GlobParam = GlobParam;
524 fDMRGnew.GlobParam.CALC_S_ON_EXIT = false;
525 fDMRGnew.GlobParam.min_halfsweeps = 30ul;
526 fDMRGnew.GlobParam.min_halfsweeps = 30ul;
527
528 DMRG::CONTROL::DYN DynParam;
529 if (dLphys == 2)
530 {
531 GlobParam.Minit = 10ul;
532 GlobParam.Qinit = 10ul;
533
534 size_t start_2site = 0ul;
535 size_t end_2site = 10ul;
536 //int period_2site = 2ul;
537 DynParam.iteration = [start_2site,end_2site] (size_t i) {return (i>=start_2site and i<=end_2site)? DMRG::ITERATION::TWO_SITE : DMRG::ITERATION::ONE_SITE;};
538 DynParam.max_alpha_rsvd = [] (size_t i) {return (i<30ul)? 1e8:0.;};
539 fDMRGnew.DynParam = DynParam;
540 }
541
542 fDMRGnew.userSetGlobParam();
543 fDMRGnew.userSetDynParam();
544
545 fDMRGnew.edgeState(Htherm, th, Q, LANCZOS::EDGE::GROUND, false);
546 th.state.entropy_skim();
547
548 if (dLphys==2)
549 {
550 lout << th.state.entropy().transpose() << endl;
551 vector<bool> ENTROPY_CHECK1, ENTROPY_CHECK2;
552
553 for (int l=0; l<2*L-1; l+=2) ENTROPY_CHECK1.push_back(abs(th.state.entropy()(l))>1e-10);
554 for (int l=1; l<2*L-1; l+=2) ENTROPY_CHECK2.push_back(abs(th.state.entropy()(l))<1e-10);
555
556 for (int l=1; l<2*L-1; l+=2)
557 {
558 bool TEST = abs(th.state.entropy()(l))<1e-10;
559 }
560 ALL = all_of(ENTROPY_CHECK1.begin(), ENTROPY_CHECK1.end(), [](const bool v){return v;}) and
561 all_of(ENTROPY_CHECK2.begin(), ENTROPY_CHECK2.end(), [](const bool v){return v;});
562 }
563 // else
564 // {
565 // vector<bool> ENTROPY_CHECK;
566 // for (int l=0; l<L-1; l+=1) ENTROPY_CHECK.push_back(abs(th.state.entropy()(l))<1e-10);
567 //
568 // ALL = all_of(ENTROPY_CHECK.begin(), ENTROPY_CHECK.end(), [](const bool v){return v;});
569 // }
570 }
571
572 for (int l=0; l<L-dLphys; l+=dLphys)
573 {
574 double val = 0.;
575 if (dLphys==1)
576 {
577 val = avg(th.state, Htherm.SdagS(l,l,0,1), th.state);
578 }
579 else
580 {
581 val = avg(th.state, Htherm.SdagS(l,l+1,0,0), th.state);
582 }
583 cout << termcolor::yellow
584 << "l=" << l
585 // << ", avg(th.state, Htherm.SdagS(NN), th.state)=" << avg(th.state, Htherm.SdagS(l,l+dLphys), th.state)
586 << ", avg_loc=" << val
587 << termcolor::reset
588 << endl;
589 }
590
591 PhiT = th.state.template cast<typename Hamiltonian::Mpo::Scalar_>();
592 }
593
594 // continue with PhiT
595 PhiT.eps_truncWeight = tol_compr_beta;
596 PhiT.min_Nsv = 0ul;
597 PhiT.max_Nsv = Mlim;
598 TDVPPropagator<Hamiltonian,
599 typename Hamiltonian::Symmetry,
600 typename Hamiltonian::Mpo::Scalar_,
601 typename Hamiltonian::Mpo::Scalar_,
602 Mps<Symmetry,typename Hamiltonian::Mpo::Scalar_>
603 > TDVPT(Hprop,PhiT);
604
605 vector<double> lnZvec;
606
607 vector<double> betavals;
608 vector<double> betasteps;
609 betavals.push_back(0.01);
610 betasteps.push_back(0.01);
611// double s_betainit = log((dLphys==2)?Hprop.locBasis(0).size():Hprop.locBasis(0).size()/2);
612
613 for (int i=1; i<22; ++i)
614 {
615 double beta_last = betavals[betavals.size()-1];
616 if (beta_last > betamax) break;
617 if (i>=20)
618 {
619 betasteps.push_back(0.05);
620 betavals.push_back(beta_last+0.05);
621 }
622 else
623 {
624 betasteps.push_back(0.01);
625 betavals.push_back(beta_last+0.01);
626 }
627 }
628
629 while (betavals[betavals.size()-1] < betamax)
630 {
631 betasteps.push_back(dbeta);
632 double beta_last = betavals[betavals.size()-1];
633 betavals.push_back(beta_last+dbeta);
634
635 }
636 if (betavals[betavals.size()-1] > betamax+0.005) // needs offset, otherwise random behaviour results from comparing floating point numbers
637 {
638 // lout << "popping last value " << betavals[betavals.size()-1] << "\t" << betamax << endl;
639 betavals.pop_back();
640 betasteps.pop_back();
641 }
642 // for (int i=0; i<betavals.size(); ++i)
643 // {
644 // lout << "betaval=" << betavals[i] << ", betastep=" << betasteps[i] << endl;
645 // }
646 lout << endl;
647
648 ofstream BetaFiler(make_string(wd,"/thermodyn_",th_label,".dat"));
649 BetaFiler << "#T\tβ\tc\te\tchi\tchiz\tchic\ts";
650 if constexpr (Hamiltonian::FAMILY == HUBBARD or Hamiltonian::FAMILY == KONDO) BetaFiler << "\tnphys";
651 if constexpr (Hamiltonian::FAMILY == HEISENBERG) BetaFiler << "\tSpSm\tSzSz";
652 BetaFiler << "truncWeightGlob\ttruncWeightLoc\tMmax";
653 BetaFiler << endl;
654
655 Mpo<Symmetry,typename Hamiltonian::Mpo::Scalar_> Hpropsq;
656 if (Ntaylor>0)
657 {
658 lout << "computing H^2..." << endl;
659 Hpropsq = prod(Hprop,Hprop);
660 lout << "H^2 done!" << endl;
661 }
662
663 Mpo<Symmetry,typename Hamiltonian::Mpo::Scalar_> U;
664
665 for (int i=0; i<betasteps.size(); ++i)
666 {
667 Stopwatch<> betaStepper;
668 double beta = betavals[i];
669
670 if (i<Ntaylor)
671 {
672 if (i==0 or betasteps[i]!=betasteps[i-1])
673 {
674 Mpo<Symmetry,typename Hamiltonian::Mpo::Scalar_> Id = Hprop.Identity(PhiT.locBasis());
675
676 Mpo<Symmetry,typename Hamiltonian::Mpo::Scalar_> Hmpo1 = Hprop;
677 Hmpo1.scale(-0.5*betasteps[i]); // -1/2*dβ*H
678
679 Mpo<Symmetry,typename Hamiltonian::Mpo::Scalar_> Hmpo2 = Hpropsq;
680 Hmpo2.scale(0.125*betasteps[i]*betasteps[i]); // +1/8*(dβ*H)^2
681
682 lout << "computing 1-1/2*dβ*H+1/8*(dβ*H)^2 for dβ=" << betasteps[i] << endl;
683 U = sum(Id,sum(Hmpo1,Hmpo2));
684 }
685
686 //Mps<Symmetry,Scalar> PhiTmp;
687 lout << "applying 1-1/2*dβ*H+1/8*(dβ*H)^2 and compressing for dβ=" << betasteps[i] << endl;
688 //OxV_exact(U,PhiT,PhiTmp,1e-9,DMRG::VERBOSITY::ON_EXIT,200,1,Mlim,DMRG::BROOM::QR);
689 //HxV(U,PhiT,PhiTmp,true,tol_compr_beta,50,Mlim);
690 //OxV(U,U,PhiT,true,tol_compr_beta,50,Mlim,16,6,false);
691 OxV(U,U,PhiT,true,tol_compr_beta,50,Mlim,16,6,false);
692 //PhiT = PhiTmp;
693
694 if (i==Ntaylor-1)
695 {
696 TDVPT = TDVPPropagator<Hamiltonian,
697 typename Hamiltonian::Symmetry,
698 typename Hamiltonian::Mpo::Scalar_,
699 typename Hamiltonian::Mpo::Scalar_,
700 Mps<Symmetry,typename Hamiltonian::Mpo::Scalar_>
701 >(Hprop,PhiT);
702 }
703 }
704 else
705 {
706 if (beta>betaswitch)
707 {
708 TDVPT.t_step0(Hprop, PhiT, -0.5*betasteps[i], 1);
709 }
710 else
711 {
712 //if (i==0)
713 //{
714 // PhiT.eps_truncWeight = 0.;
715 //}
716 //else
717 //{
718 PhiT.eps_truncWeight = tol_compr_beta;
719 //}
720 TDVPT.t_step(Hprop, PhiT, -0.5*betasteps[i], 1);
721 }
722 }
723 double norm = isReal(dot(PhiT,PhiT));
724 lnZvec.push_back(log(norm));
725 PhiT /= sqrt(norm);
726 lout << TDVPT.info() << endl;
727 size_t standard_precision = cout.precision();
728 lout << setprecision(16) << PhiT.info() << setprecision(standard_precision) << endl;
729
730 bool INTERMEDIATE_SAVEPOINT = false;
731 string saveLabel = "";
732 for (int is=0; is<stateSavePoints.size(); ++is)
733 {
734 if (abs(beta-stateSavePoints[is]) < 1e-8)
735 {
736 INTERMEDIATE_SAVEPOINT = true;
737 saveLabel = stateSaveLabels[is];
738 }
739 }
740 if (SAVE_BETA and INTERMEDIATE_SAVEPOINT)
741 {
742 lout << termcolor::green << "saving the intermediate beta result to: " << saveLabel << " in directory=" << wd << termcolor::reset << endl;
743 PhiT.save(make_string(wd,"/betaRes_",saveLabel));
744 }
745
746 double avg_H = isReal(avg(PhiT,Hprop,PhiT));
747
748 double e = avg_H/L;
749
750 double c = std::nan("0");
751 if (CALC_C)
752 {
753 double avg_Hsq = (Hprop.maxPower()==1)? isReal(avg(PhiT,Hprop,Hprop,PhiT)):isReal(avg(PhiT,Hprop,PhiT,2));
754 double avgH_sq = pow(avg_H,2);
755 c = isReal(beta*beta*(avg_Hsq-avgH_sq))/L;
756 }
757
758 double chi = std::nan("0");
759 double chiz = std::nan("0");
760 double chic = 0;
761 if (CALC_CHI)
762 {
763 #if defined(USING_SU2xU1) or defined(USING_SU2xSU2) or defined(USING_SU2)
764 chi = isReal(beta*avg(PhiT, Hprop.Sdagtot(0,sqrt(3.),dLphys), Hprop.Stot(0,1.,dLphys), PhiT))/L;
765 chiz = chi/3.;
766 for (int l=0; l<dLphys*L; l+=dLphys)
767 {
768 chic += isReal(beta*avg(PhiT, Hprop.SdagS(l,dLphys*L/2), PhiT))/3.;
769 }
770 #else
771 chi = 0.5*isReal(beta*avg(PhiT, Hprop.Scomptot(SP,0,1.,dLphys), Hprop.Scomptot(SM,0,1.,dLphys), PhiT))/L;
772 chi += 0.5*isReal(beta*avg(PhiT, Hprop.Scomptot(SM,0,1.,dLphys), Hprop.Scomptot(SP,0,1.,dLphys), PhiT))/L;
773 chiz = isReal(beta*avg(PhiT, Hprop.Scomptot(SZ,0,1.,dLphys), Hprop.Scomptot(SZ,0,1.,dLphys), PhiT))/L;
774 chi += chiz;
775 for (int l=0; l<dLphys*L; l+=dLphys)
776 {
777 chic += isReal(beta*avg(PhiT, Hprop.SzSz(l,dLphys*L/2), PhiT));
778 }
779 #endif
780 }
781
782// double chi_ = beta*(2.*avg(PhiT,Hchi,PhiT)/L+0.75);
783
784 VectorXd tmp = VectorXd::Map(lnZvec.data(), lnZvec.size());
785 double s = s_betainit + tmp.sum()/L + beta*e;
786// cout << "s_betainit=" << s_betainit << ", tmp.sum()/L=" << tmp.sum()/L << ", beta*e=" << beta*e << ", total=" << s << endl;
787// svec.push_back(s);
788
789 auto PhiTtmp = PhiT; PhiTtmp.entropy_skim();
790 cout << "aa" << endl;
791 lout << "S=" << PhiTtmp.entropy().transpose() << endl;
792
793 VectorXd nphys(Lcell); nphys.setZero();
794 //VectorXd SdagSphys(Lcell); for (int j=0; j<Lcell; ++j) SdagSphys(j) = 0.;
795 double nancl = 0.;
796 cout << "a" << endl;
797 #if not defined(USING_SU2xSU2)
798 if constexpr (Hamiltonian::FAMILY == HUBBARD or Hamiltonian::FAMILY == KONDO)
799 {
800 cout << "b" << endl;
801 for (int j=0, icell=0; j<dLphys*L; j+=dLphys, icell+=1)
802 {
803 nphys(icell%Lcell) += isReal(avg(PhiT, Hprop.n(j,0), PhiT));
804 }
805 for (int j=dLphys-1; j<dLphys*L; j+=dLphys)
806 {
807 nancl += isReal(avg(PhiT, Hprop.n(j,dLphys%2), PhiT));
808 }
809 }
810 #endif
811 #if not defined(USING_SU2)
812 if constexpr (Hamiltonian::FAMILY == HEISENBERG)
813 {
814 cout << "c" << endl;
815 for (int j=0, icell=0; j<dLphys*L; j+=dLphys, icell+=1)
816 {
817 nphys(icell%Lcell) += isReal(avg(PhiT, Hprop.Sz(j,0), PhiT));
818 }
819 for (int j=dLphys-1; j<dLphys*L; j+=dLphys)
820 {
821 nancl += isReal(avg(PhiT, Hprop.Sz(j,dLphys%2), PhiT));
822 }
823// for (int j=0, icell=0; j<dLphys*(L-1); j+=dLphys, icell+=1)
824// {
825// //cout << "j,j+dLphys=" << j << "," << j+dLphys << ", icell=" << icell << ", icellmodLcell=" << icell%Lcell << ", val=" << isReal(avg(PhiT, Hprop.SdagS(j,j+dLphys,0,0), PhiT)) << endl;
826// SdagSphys(icell%Lcell) += isReal(avg(PhiT, Hprop.SdagS(j,j+dLphys,0,0), PhiT));
827// }
828// }
829// double SzSz = 0.;
830// double SpSm = 0.;
831// if constexpr (Hamiltonian::FAMILY == HEISENBERG)
832// {
833// SpSm = isReal(avg(PhiT, Hprop.SpSm(L/4,3*L/4), PhiT));
834// SzSz = isReal(avg(PhiT, Hprop.SzSz(L/4,3*L/4), PhiT));
835 }
836 #endif
837 cout << "d" << endl;
838
839 nphys /= L;
840 double nphystot = nphys.sum();
841 nancl /= L;
842
843 lout << termcolor::bold << "β=" << beta << ", T=" << 1./beta << ", c=" << c << ", e=" << e << ", chi=" << chi;
844 #ifndef USING_SU2
845 lout << ", chiz=" << chiz << "(3x=" << 3.*chiz << ")";
846 #else
847 lout << ", chic=" << chic;
848 #endif
849 lout << ", s=" << s;
850 BetaFiler << setprecision(16) << 1./beta << "\t" << beta << "\t" << c << "\t" << e << "\t" << chi << "\t" << chiz << "\t" << chic << "\t" << s;
851 if constexpr (Hamiltonian::FAMILY == HUBBARD or Hamiltonian::FAMILY == KONDO)
852 {
853 lout << ", nphys=" << nphys.transpose() << ", sum=" << nphystot << ", nancl=" << nancl;
854 BetaFiler << "\t" << nphystot;
855 }
856 else if constexpr (Hamiltonian::FAMILY == HEISENBERG)
857 {
858 lout << ", Mphys=" << nphys.transpose() << ", sum=" << nphystot << ", Mancl=" << nancl;
859 BetaFiler << "\t" << nphystot;
860// lout << ", SpSm(" << L/4 << "," << 3*L/4 << ")=" << SpSm << ", SzSz(" << L/4 << "," << 3*L/4 << ")=" << SzSz;
861// BetaFiler << "\t" << SpSm << "\t" << SzSz;
862 }
863 BetaFiler << "\t" << PhiT.get_truncWeight().sum() << "\t" << PhiT.get_truncWeight().maxCoeff() << "\t" << PhiT.calc_Mmax();
864// if constexpr (Hamiltonian::FAMILY == HEISENBERG)
865// {
866// double Nb = L-1;
867// lout << ", SdagSphys/Nb=" << SdagSphys.transpose()/Nb << ", sum/Nb=" << SdagSphys.sum()/Nb;
868// BetaFiler << "\t" << SdagSphys.transpose()/Nb;
869// }
870 BetaFiler << endl;
871 lout << termcolor::reset << endl;
872 lout << betaStepper.info("βstep") << endl;
873 lout << endl;
874 }
875 BetaFiler.close();
876 }
877
878 if (SAVE_BETA and !LOAD_BETA)
879 {
880 lout << termcolor::green << "saving the final beta result to: " << th_label << " in directory=" << wd << termcolor::reset << endl;
881 PhiT.save(make_string(wd,"/betaRes_",th_label));
882 }
883}
884
885template<typename Hamiltonian>
886void SpectralManager<Hamiltonian>::
887continue_beta_propagation (const Hamiltonian &Hprop, int Lcell, int dLphys,
888 double s_betainit_input, double betainit, double betamax, double dbeta, double tol_compr_beta, size_t Mlim, qarray<Hamiltonian::Symmetry::Nq> Q, double betaswitch,
889 string wd, string th_label, string LOAD_BETA, bool SAVE_BETA,
890 DMRG::VERBOSITY::OPTION VERB,
891 vector<double> stateSavePoints, vector<string> stateSaveLabels,
892 bool CALC_C, bool CALC_CHI)
893{
894 for (const auto &spec:specs)
895 {
896 assert(CHECK_SPEC(spec) and "Wrong spectral abbreviation!");
897 }
898 L = Hprop.length()/dLphys;
899 Nspec = specs.size();
900
901 lout << "loading beta result from: " << LOAD_BETA << endl;
902 PhiT.load(LOAD_BETA);
903 lout << "loaded: " << PhiT.info() << endl;
904
905 PhiT.eps_truncWeight = tol_compr_beta;
906 PhiT.min_Nsv = 0ul;
907 PhiT.max_Nsv = Mlim;
908 TDVPPropagator<Hamiltonian,
909 typename Hamiltonian::Symmetry,
910 typename Hamiltonian::Mpo::Scalar_,
911 typename Hamiltonian::Mpo::Scalar_,
912 Mps<Symmetry,typename Hamiltonian::Mpo::Scalar_>
913 > TDVPT(Hprop,PhiT);
914
915 vector<double> lnZvec;
916 vector<double> betavals;
917 vector<double> betasteps;
918
919 assert(betainit>=0.2);
920 betavals.push_back(betainit+dbeta);
921 betasteps.push_back(dbeta);
922 lout << "betaval=" << betavals[betavals.size()-1] << ", betastep=" << betasteps[betasteps.size()-1] << endl;
923
924 double e = isReal(avg(PhiT,Hprop,PhiT)/static_cast<double>(L));
925 double s_betainit = s_betainit_input;
926 s_betainit -= betainit*e;
927
928 while (betavals[betavals.size()-1] < betamax)
929 {
930 betasteps.push_back(dbeta);
931 double beta_last = betavals[betavals.size()-1];
932 betavals.push_back(beta_last+dbeta);
933 lout << "betaval=" << betavals[betavals.size()-1] << ", betastep=" << betasteps[betasteps.size()-1] << endl;
934 }
935 if (betavals[betavals.size()-1] > betamax+0.005) // needs offset, otherwise random behaviour results from comparing floating point numbers
936 {
937 betavals.pop_back();
938 betasteps.pop_back();
939 }
940 lout << endl;
941
942 ofstream BetaFiler(make_string(wd,"/thermodyn_",th_label,".dat"));
943 BetaFiler << "#T\tβ\tc\te\tchi\tchiz\tchic\ts";
944 if constexpr (Hamiltonian::FAMILY == HUBBARD or Hamiltonian::FAMILY == KONDO) BetaFiler << "\tnphys";
945 if constexpr (Hamiltonian::FAMILY == HEISENBERG) BetaFiler << "\tSpSm\tSzSz";
946 BetaFiler << "truncWeightGlob\ttruncWeightLoc\tMmax";
947 BetaFiler << endl;
948
949 for (int i=0; i<betasteps.size(); ++i)
950 {
951 Stopwatch<> betaStepper;
952 double beta = betavals[i];
953
954 if (beta>betaswitch)
955 {
956 PhiT.eps_truncWeight = 0.;
957 TDVPT.t_step0(Hprop, PhiT, -0.5*betasteps[i], 1);
958 }
959 else
960 {
961 PhiT.eps_truncWeight = tol_compr_beta;
962 TDVPT.t_step(Hprop, PhiT, -0.5*betasteps[i], 1);
963 }
964 double norm = isReal(dot(PhiT,PhiT));
965 lnZvec.push_back(log(norm));
966 PhiT /= sqrt(norm);
967 lout << TDVPT.info() << endl;
968 size_t standard_precision = cout.precision();
969 lout << setprecision(16) << PhiT.info() << setprecision(standard_precision) << endl;
970
971 bool INTERMEDIATE_SAVEPOINT = false;
972 string saveLabel = "";
973 for (int is=0; is<stateSavePoints.size(); ++is)
974 {
975 if (abs(beta-stateSavePoints[is]) < 1e-8)
976 {
977 INTERMEDIATE_SAVEPOINT = true;
978 saveLabel = stateSaveLabels[is];
979 }
980 }
981 if (SAVE_BETA and INTERMEDIATE_SAVEPOINT)
982 {
983 lout << termcolor::green << "saving the intermediate beta result to: " << saveLabel << " in directory=" << wd << termcolor::reset << endl;
984 PhiT.save(make_string(wd,"/betaRes_",saveLabel));
985 }
986
987 double avg_H = isReal(avg(PhiT,Hprop,PhiT));
988
989 e = avg_H/L;
990
991 double c = std::nan("0");
992 if (CALC_C)
993 {
994 double avg_Hsq = (Hprop.maxPower()==1)? isReal(avg(PhiT,Hprop,Hprop,PhiT)):isReal(avg(PhiT,Hprop,PhiT,2));
995 double avgH_sq = pow(avg_H,2);
996 c = isReal(beta*beta*(avg_Hsq-avgH_sq))/L;
997 }
998
999 double chi = std::nan("0");
1000 double chiz = std::nan("0");
1001 double chic = 0;
1002 if (CALC_CHI)
1003 {
1004 #ifdef USING_SU2
1005 chi = isReal(beta*avg(PhiT, Hprop.Sdagtot(0,sqrt(3.),dLphys), Hprop.Stot(0,1.,dLphys), PhiT))/L;
1006 chiz = chi/3.;
1007 for (int l=0; l<dLphys*L; l+=dLphys)
1008 {
1009 chic += isReal(beta*avg(PhiT, Hprop.SdagS(l,dLphys*L/2), PhiT))/3.;
1010 }
1011 #else
1012 chi = 0.5*isReal(beta*avg(PhiT, Hprop.Scomptot(SP,0,1.,dLphys), Hprop.Scomptot(SM,0,1.,dLphys), PhiT))/L;
1013 chi += 0.5*isReal(beta*avg(PhiT, Hprop.Scomptot(SM,0,1.,dLphys), Hprop.Scomptot(SP,0,1.,dLphys), PhiT))/L;
1014 chiz = isReal(beta*avg(PhiT, Hprop.Scomptot(SZ,0,1.,dLphys), Hprop.Scomptot(SZ,0,1.,dLphys), PhiT))/L;
1015 chi += chiz;
1016 for (int l=0; l<dLphys*L; l+=dLphys)
1017 {
1018 chic += isReal(beta*avg(PhiT, Hprop.SzSz(l,dLphys*L/2), PhiT));
1019 }
1020 #endif
1021 }
1022
1023 VectorXd tmp = VectorXd::Map(lnZvec.data(), lnZvec.size());
1024 double s = s_betainit + tmp.sum()/L + beta*e;
1025
1026 auto PhiTtmp = PhiT; PhiTtmp.entropy_skim();
1027 lout << "S=" << PhiTtmp.entropy().transpose() << endl;
1028
1029 VectorXd nphys(Lcell); for (int j=0; j<Lcell; ++j) nphys(j) = 0.;
1030 double nancl = 0.;
1031 if constexpr (Hamiltonian::FAMILY == HUBBARD or Hamiltonian::FAMILY == KONDO)
1032 {
1033 for (int j=0, icell=0; j<dLphys*L; j+=dLphys, icell+=1)
1034 {
1035 nphys(icell%Lcell) += isReal(avg(PhiT, Hprop.n(j,0), PhiT));
1036 }
1037 for (int j=dLphys-1; j<dLphys*L; j+=dLphys)
1038 {
1039 nancl += isReal(avg(PhiT, Hprop.n(j,dLphys%2), PhiT));
1040 }
1041 }
1042 double SzSz = 0.;
1043 double SpSm = 0.;
1044 if constexpr (Hamiltonian::FAMILY == HEISENBERG)
1045 {
1046 SpSm = isReal(avg(PhiT, Hprop.SpSm(L/4,3*L/4), PhiT));
1047 SzSz = isReal(avg(PhiT, Hprop.SzSz(L/4,3*L/4), PhiT));
1048 }
1049
1050 nphys /= L;
1051 double nphystot = nphys.sum();
1052 nancl /= L;
1053
1054 lout << termcolor::bold << "β=" << beta << ", T=" << 1./beta << ", c=" << c << ", e=" << e << ", chi=" << chi;
1055 #ifndef USING_SU2
1056 lout << ", chiz=" << chiz << "(3x=" << 3.*chiz << ")";
1057 #else
1058 lout << ", chic=" << chic;
1059 #endif
1060 lout << ", s=" << s;
1061 BetaFiler << setprecision(16) << 1./beta << "\t" << beta << "\t" << c << "\t" << e << "\t" << chi << "\t" << chiz << "\t" << chic << "\t" << s;
1062 if constexpr (Hamiltonian::FAMILY == HUBBARD or Hamiltonian::FAMILY == KONDO)
1063 {
1064 lout << ", nphys=" << nphys.transpose() << ", sum=" << nphystot << ", nancl=" << nancl;
1065 BetaFiler << "\t" << nphystot;
1066 }
1067 if constexpr (Hamiltonian::FAMILY == HEISENBERG)
1068 {
1069 lout << ", SpSm(" << L/4 << "," << 3*L/4 << ")=" << SpSm << ", SzSz(" << L/4 << "," << 3*L/4 << ")=" << SzSz;
1070 BetaFiler << "\t" << SpSm << "\t" << SzSz;
1071 }
1072 BetaFiler << "\t" << PhiT.get_truncWeight().sum() << "\t" << PhiT.get_truncWeight().maxCoeff() << "\t" << PhiT.calc_Mmax();
1073 BetaFiler << endl;
1074 lout << termcolor::reset << endl;
1075 lout << betaStepper.info("βstep") << endl;
1076 lout << endl;
1077 }
1078
1079 BetaFiler.close();
1080
1081 if (SAVE_BETA)
1082 {
1083 lout << termcolor::green << "saving the final beta result to: " << th_label << " in directory=" << wd << termcolor::reset << endl;
1084 PhiT.save(make_string(wd,"/betaRes_",th_label));
1085 }
1086}
1087
1088template<typename Hamiltonian>
1089void SpectralManager<Hamiltonian>::
1090apply_operators_on_thermal_state (int Lcell, int dLphys, bool CHECK)
1091{
1092 PhiTt = PhiT.template cast<complex<double>>();
1093// PhiTt.eps_svd = tol_compr;
1094// PhiTt.min_Nsv = 0ul;
1095// PhiTt.max_Nsv = Mlim;
1096
1097 // OxV for time propagation
1098 vector<vector<Mpo<typename Hamiltonian::Symmetry,Scalar>>> O(Nspec);
1099 Odag.resize(Nspec);
1100 for (int z=0; z<Nspec; ++z) O[z].resize(L);
1101 for (int z=0; z<Nspec; ++z) Odag[z].resize(L);
1102
1103 for (int z=0; z<Nspec; ++z)
1104 for (int l=0; l<L; ++l)
1105 {
1106 O[z][l] = get_Op(Hwork, dLphys*l, specs[z], 1., 0, dLphys);
1107 double dagfactor;
1108 if (specs[z] == "SSF")
1109 {
1110 dagfactor = sqrt(3);
1111 }
1112 else if (specs[z] == "PES")
1113 {
1114 dagfactor = -sqrt(2);
1115 }
1116 else if (specs[z] == "IPE")
1117 {
1118 dagfactor = +sqrt(2);
1119 }
1120 else
1121 {
1122 dagfactor = 1.;
1123 }
1124 Odag[z][l] = get_Op(Hwork, dLphys*l, DAG(specs[z]), dagfactor, 0, dLphys);
1125 }
1126
1127 // shift values O→O-<O>
1128 vector<vector<Scalar>> Oshift(Nspec);
1129 vector<vector<Scalar>> Odagshift(Nspec);
1130 for (int z=0; z<Nspec; ++z)
1131 {
1132 Oshift[z].resize(L);
1133 Odagshift[z].resize(L);
1134 for (int l=0; l<L; ++l)
1135 {
1136 Oshift[z][l] = avg(PhiT, O[z][l], PhiT);
1137 Odagshift[z][l] = avg(PhiT, Odag[z][l], PhiT);
1138// lout << "spec=" << specs[z] << ", l=" << l << ", shift=" << Oshift[z][l] << endl;
1139 }
1140 }
1141
1142 for (int z=0; z<Nspec; ++z)
1143 for (int l=0; l<L; ++l)
1144 {
1145 O[z][l].scale (1.,-Oshift[z][l]);
1146 Odag[z][l].scale(1.,-Odagshift[z][l]);
1147 }
1148
1149 //---------check---------
1150 if (CHECK)
1151 {
1152 lout << endl;
1153 for (int z=0; z<Nspec; ++z)
1154 {
1155 lout << "check z=" << z
1156 << ", spec=" << specs[z] << ", dag=" << DAG(specs[z])
1157 << ", <O†[z][l]*O[z][l]>=" << avg(PhiTt, Odag[z][L/2], O[z][L/2], PhiTt)
1158 << ", <O†[z][l]>=" << avg(PhiTt, Odag[z][L/2], PhiTt)
1159 << ", <O[z][l]>=" << avg(PhiTt, O[z][L/2], PhiTt)
1160// << ", g.state: " << avg(g.state, Odag[z][L/2], O[z][L/2], g.state)
1161 << endl;
1162 }
1163 }
1164
1165 OxPhiTt.resize(Nspec);
1166 for (int z=0; z<Nspec; ++z)
1167 {
1168 OxPhiTt[z].resize(Lcell);
1169 for (int i=0; i<Lcell; ++i)
1170 {
1171 OxV_exact(O[z][L/2+i], PhiTt, OxPhiTt[z][i], 2., DMRG::VERBOSITY::ON_EXIT, 200, 1);
1172 }
1173 }
1174}
1175
1176template<typename Hamiltonian>
1177void SpectralManager<Hamiltonian>::
1178set_measurement (int iz, string spec, double factor, int dLphys, qarray<Symmetry::Nq> Q, int Lcell, int measure_interval, string measure_name, string measure_subfolder, bool TRANSFORM)
1179{
1180 assert(Green.size()>0);
1181
1182 vector<Mpo<Symmetry,Scalar>> Measure(Hwork.length()/dLphys);
1183 for (int l=0, iphys=0; l<Hwork.length(); l+=dLphys, iphys+=1)
1184 {
1185 Measure[iphys] = get_Op(Hwork, l, spec, factor, 0, dLphys);
1186 if (TRANSFORM) Measure[l].transform_base(Q,false,L); // PRINT=false
1187 }
1188 Green[iz].set_measurement(Measure, Lcell, measure_interval, measure_name, measure_subfolder);
1189}
1190
1191template<typename Hamiltonian>
1192void SpectralManager<Hamiltonian>::
1193resize_Green (string wd, string label, int Ns, double tmax, double dt, double wmin, double wmax, int wpoints, Q_RANGE QR, int qpoints, GREEN_INTEGRATION INT)
1194{
1195 int Nt = static_cast<int>(tmax/dt);
1196
1197 // GreenPropagator
1198 Green.resize(Nspec);
1199 for (int z=0; z<Nspec; ++z)
1200 {
1201 string spec = specs[z];
1202 Green[z] = GreenPropagator<Hamiltonian,Symmetry,Scalar,complex<double>>
1203 (wd+spec+"_"+label,tmax,Nt,Ns,wmin,wmax,wpoints,QR,qpoints,INT);
1204 Green[z].set_verbosity(DMRG::VERBOSITY::ON_EXIT);
1205 }
1206 Green[0].set_verbosity(CHOSEN_VERB);
1207}
1208
1209template<typename Hamiltonian>
1210void SpectralManager<Hamiltonian>::
1211compute (string wd, string label, int Ns, double tmax, double dt, double wmin, double wmax, int wpoints, Q_RANGE QR, int qpoints, GREEN_INTEGRATION INT, size_t Mlim, double tol_DeltaS, double tol_compr)
1212{
1213 if (Green.size()==0)
1214 {
1215 resize_Green(wd,label,Ns,tmax,dt,wmin,wmax,wpoints,QR,qpoints,INT);
1216 }
1217
1218 // Propagation
1219 #pragma omp parallel for
1220 for (int z=0; z<Nspec; ++z)
1221 {
1222 string spec = specs[z];
1223 Green[z].set_tol_DeltaS(tol_DeltaS);
1224 Green[z].set_lim_Nsv(Mlim);
1225 Green[z].set_tol_compr(tol_compr);
1226
1227 Green[z].compute_cell(Hwork, OxPhiCellBra[z], OxPhiCellKet[z], Eg, TIME_DIR(spec), true); // COUNTERPROPAGATE=true
1228 Green[z].save(false); // IGNORE_TX=false
1229 }
1230}
1231
1232template<typename Hamiltonian>
1233void SpectralManager<Hamiltonian>::
1234compute_finite (size_t j0, string wd, string label, int Ns, double tmax, double dt, double wmin, double wmax, int wpoints, GREEN_INTEGRATION INT, size_t Mlim, double tol_DeltaS, double tol_compr)
1235{
1236 if (Green.size()==0)
1237 {
1238 resize_Green(wd,label,Ns,tmax,dt,wmin,wmax,wpoints,ZERO_2PI,501,INT);
1239 }
1240
1241 // Propagation
1242 #pragma omp parallel for
1243 for (int z=0; z<Nspec; ++z)
1244 {
1245 string spec = specs[z];
1246 Green[z].set_tol_DeltaS(tol_DeltaS);
1247 Green[z].set_lim_Nsv(Mlim);
1248 Green[z].set_tol_compr(tol_compr);
1249
1250 Green[z].set_OxPhiFull(OxPhiCellKet[z]);
1251 Green[z].compute_one(j0, Hwork, OxPhiCellKet[z], Eg, TIME_DIR(spec), false); // COUNTERPROPAGATE=false
1252 Green[z].save(false); // IGNORE_TX=false
1253 }
1254}
1255
1256template<typename Hamiltonian>
1257void SpectralManager<Hamiltonian>::
1258compute_finiteCell (int Lcell, int x0, string wd, string label, double tmax, double dt, double wmin, double wmax, int wpoints, Q_RANGE QR, int qpoints, GREEN_INTEGRATION INT, size_t Mlim, double tol_DeltaS, double tol_compr)
1259{
1260 if (Green.size()==0)
1261 {
1262 resize_Green(wd,label,1,tmax,dt,wmin,wmax,wpoints,QR,qpoints,INT);
1263 }
1264
1265 for (int z=0; z<Nspec; ++z)
1266 {
1267 Green[z].set_Lcell(Lcell);
1268 }
1269
1270 // propagation
1271 #pragma omp parallel for
1272 for (int z=0; z<Nspec; ++z)
1273 {
1274 string spec = specs[z];
1275 Green[z].set_tol_DeltaS(tol_DeltaS);
1276 Green[z].set_lim_Nsv(Mlim);
1277 Green[z].set_tol_compr(tol_compr);
1278
1279 Green[z].set_OxPhiFull(OxPhiCellKet[z]);
1280 if (Oshift.size() != 0)
1281 {
1282 if (Oshift[z][0] != 0.)
1283 {
1284 Green[z].set_Oshift(Oshift[z]);
1285 }
1286 }
1287 Green[z].compute_cell(Hwork, OxPhiCellKet[z], Eg, TIME_DIR(spec), false, x0); // COUNTERPROPAGATE=false
1288 Green[z].save(false); // IGNORE_TX=false
1289 }
1290}
1291
1292template<typename Hamiltonian>
1293void SpectralManager<Hamiltonian>::
1294compute_thermal (string wd, string label, int dLphys, double tmax, double dt, double wmin, double wmax, int wpoints, Q_RANGE QR, int qpoints, GREEN_INTEGRATION INT, size_t Mlim, double tol_DeltaS, double tol_compr)
1295{
1296 if (Green.size()==0)
1297 {
1298 resize_Green(wd,label,1,tmax,dt,wmin,wmax,wpoints,QR,qpoints,INT);
1299 }
1300
1301 #pragma omp parallel for
1302 for (int z=0; z<Nspec; ++z)
1303 {
1304 string spec = specs[z];
1305 Green[z].set_tol_DeltaS(tol_DeltaS);
1306 Green[z].set_lim_Nsv(Mlim);
1307 Green[z].set_tol_compr(tol_compr);
1308
1309 Green[z].compute_thermal_cell(Hwork, Odag[z], PhiTt, OxPhiTt[z], TIME_DIR(spec));
1310// Green[z].FT_allSites();
1311 Green[z].save(false); // IGNORE_TX=false
1312 }
1313}
1314
1315template<typename Hamiltonian>
1316void SpectralManager<Hamiltonian>::
1317FTcell_xq()
1318{
1319 assert(Green.size() == Nspec);
1320 for (int z=0; z<Nspec; ++z)
1321 {
1322 Green[z].FTcell_xq();
1323 Green[z].save_GtqCell();
1324 }
1325}
1326
1327template<typename Hamiltonian>
1328void SpectralManager<Hamiltonian>::
1329make_A1P (GreenPropagator<Hamiltonian,Symmetry,Scalar,complex<double> > &Gfull, string wd, string label, int Ns, double tmax, double wmin, double wmax, int wpoints, Q_RANGE QR, int qpoints, GREEN_INTEGRATION INT, bool SAVE_N_MU)
1330{
1331 auto itPES = find(specs.begin(), specs.end(), "PES");
1332 auto itIPE = find(specs.begin(), specs.end(), "IPE");
1333
1334 if (itPES != specs.end() and itIPE != specs.end())
1335 {
1336 int iPES = distance(specs.begin(), itPES);
1337 int iIPE = distance(specs.begin(), itIPE);
1338
1339 // Add PES+IPE
1340 vector<vector<MatrixXcd>> GinA1P(L); for (int i=0; i<L; ++i) {GinA1P[i].resize(L);}
1341 for (int i=0; i<L; ++i)
1342 for (int j=0; j<L; ++j)
1343 {
1344 GinA1P[i][j].resize(Green[iPES].get_GtxCell()[0][0].rows(),
1345 Green[iPES].get_GtxCell()[0][0].cols());
1346 GinA1P[i][j].setZero();
1347 }
1348
1349 for (int i=0; i<L; ++i)
1350 for (int j=0; j<L; ++j)
1351 {
1352 GinA1P[i][j] += Green[iPES].get_GtxCell()[i][j] + Green[iIPE].get_GtxCell()[i][j];
1353 }
1354 Gfull = GreenPropagator<Hamiltonian,Symmetry,Scalar,complex<double> >(wd+"A1P"+"_"+label,L,Ncells,Ns,tmax,GinA1P,QR,qpoints,INT);
1355 Gfull.recalc_FTwCell(wmin,wmax,wpoints);
1356
1357 if (SAVE_N_MU)
1358 {
1359 IntervalIterator mu(wmin,wmax,501);
1360 for (mu=mu.begin(); mu!=mu.end(); ++mu)
1361 {
1362 mu << Gfull.integrate_Glocw_cell(*mu);
1363 mu.save(make_string(wd,"n(μ)_tmax=",tmax,"_",label,".dat"));
1364 }
1365 }
1366 }
1367}
1368
1369template<typename Hamiltonian>
1370void SpectralManager<Hamiltonian>::
1371make_A1P_finite (GreenPropagator<Hamiltonian,Symmetry,Scalar,complex<double> > &Gfull, string wd, string label, double tmax, double wmin, double wmax, int wpoints, GREEN_INTEGRATION INT)
1372{
1373 auto itPES = find(specs.begin(), specs.end(), "PES");
1374 auto itIPE = find(specs.begin(), specs.end(), "IPE");
1375
1376 if (itPES != specs.end() and itIPE != specs.end())
1377 {
1378 int iPES = distance(specs.begin(), itPES);
1379 int iIPE = distance(specs.begin(), itIPE);
1380
1381 // Add PES+IPE
1382 MatrixXcd GinA1P;
1383 GinA1P.resize(Green[iPES].get_Gtx().rows(), Green[iPES].get_Gtx().cols());
1384 GinA1P.setZero();
1385 GinA1P += Green[iPES].get_Gtx() + Green[iIPE].get_Gtx();
1386
1387 Gfull = GreenPropagator<Hamiltonian,Symmetry,Scalar,complex<double> >(wd+"A1P"+"_"+label,L,tmax,GinA1P,INT);
1388 Gfull.recalc_FTw(wmin,wmax,wpoints);
1389 }
1390}
1391
1392template<typename Hamiltonian>
1393void SpectralManager<Hamiltonian>::
1394reload (string wd, const vector<string> &specs_input, string label, int L_input, int Ncells_input, int Ns, double tmax, double wmin, double wmax, int wpoints, Q_RANGE QR, int qpoints, GREEN_INTEGRATION INT)
1395{
1396 L = L_input;
1397 Ncells = Ncells_input;
1398 Lhetero = L*Ncells;
1399 x0 = Lhetero/2;
1400 specs = specs_input;
1401 Nspec = specs.size();
1402
1403 for (const auto &spec:specs)
1404 {
1405 assert(CHECK_SPEC(spec) and "Wrong spectral abbreviation!");
1406 }
1407
1408 Green.resize(Nspec);
1409
1410 for (int z=0; z<Nspec; ++z)
1411 {
1412 string spec = specs[z];
1413 Green[z] = GreenPropagator<Hamiltonian,Symmetry,Scalar,complex<double>>
1414 (wd+spec+"_"+label,L,Ncells,Ns,tmax,{wd+spec+"_"+label+make_string("_L=",L,"x",Ncells,"_tmax=",tmax)},QR,qpoints,INT);
1415 Green[z].recalc_FTwCell(wmin,wmax,wpoints, TIME_DIR(spec));
1416 Green[z].save(true); // IGNORE_TX=true
1417 }
1418}
1419
1420template<typename Hamiltonian>
1421Mpo<typename Hamiltonian::Symmetry,typename Hamiltonian::Mpo::Scalar_> SpectralManager<Hamiltonian>::
1422get_Op (const Hamiltonian &H, size_t loc, std::string spec, double factor, size_t locy, int dLphys)
1423{
1424 Mpo<Symmetry,Scalar> Res;
1425 bool OPERATOR_IS_SET = false;
1426
1427 // spin structure factor
1428 if (spec == "SSF")
1429 {
1430 if constexpr (Symmetry::IS_SPIN_SU2())
1431 {
1432 Res = H.S(loc,locy,factor);
1433 OPERATOR_IS_SET = true;
1434 }
1435 else
1436 {
1437 Res = H.Scomp(SP,loc,locy);
1438 OPERATOR_IS_SET = true;
1439 }
1440 }
1441 else if (spec == "SDAGSF")
1442 {
1443 if constexpr (Symmetry::IS_SPIN_SU2())
1444 {
1445 Res = H.Sdag(loc,locy,factor);
1446 OPERATOR_IS_SET = true;
1447 }
1448 else
1449 {
1450 Res = H.Scomp(SM,loc,locy);
1451 OPERATOR_IS_SET = true;
1452 }
1453 }
1454 else if (spec == "SSZ")
1455 {
1456 if constexpr (!Symmetry::IS_SPIN_SU2())
1457 {
1458 Res = H.Scomp(SZ,loc,locy);
1459 OPERATOR_IS_SET = true;
1460 }
1461 else
1462 {
1463 throw;
1464 }
1465 }
1466 else if (spec == "SPM")
1467 {
1468 if constexpr (!Symmetry::IS_SPIN_SU2())
1469 {
1470 Res = H.Scomp(SM,loc,locy);
1471 OPERATOR_IS_SET = true;
1472 }
1473 else
1474 {
1475 throw;
1476 }
1477 }
1478 else if (spec == "SMP")
1479 {
1480 if constexpr (!Symmetry::IS_SPIN_SU2())
1481 {
1482 Res = H.Scomp(SP,loc,locy);
1483 OPERATOR_IS_SET = true;
1484 }
1485 else
1486 {
1487 throw;
1488 }
1489 }
1490 if constexpr (Hamiltonian::FAMILY == HEISENBERG)
1491 {
1492 if (spec == "QSF")
1493 {
1494 if constexpr (Symmetry::IS_SPIN_SU2())
1495 {
1496 Res = H.Q(loc,locy,factor);
1497 OPERATOR_IS_SET = true;
1498 }
1499 else
1500 {
1501 Res = H.Scomp(SP,loc,locy);
1502 OPERATOR_IS_SET = true;
1503 }
1504 }
1505 else if (spec == "QDAGSF")
1506 {
1507 if constexpr (Symmetry::IS_SPIN_SU2())
1508 {
1509 Res = H.Qdag(loc,locy,factor);
1510 OPERATOR_IS_SET = true;
1511 }
1512 else
1513 {
1514 Res = H.Scomp(SM,loc,locy);
1515 OPERATOR_IS_SET = true;
1516 }
1517 }
1518 }
1519 if constexpr (Hamiltonian::FAMILY == HUBBARD or Hamiltonian::FAMILY == KONDO)
1520 {
1521 // photemission
1522 if (spec == "PES" or spec == "PESUP")
1523 {
1524 if constexpr (Symmetry::IS_SPIN_SU2()) // or spinless
1525 {
1526 Res = H.c(loc,locy,factor);
1527 OPERATOR_IS_SET = true;
1528 }
1529 else
1530 {
1531 Res = H.template c<UP>(loc,locy);
1532 OPERATOR_IS_SET = true;
1533 }
1534 }
1535 else if (spec == "PESDN")
1536 {
1537 if constexpr (Symmetry::IS_SPIN_SU2()) // or spinless
1538 {
1539 Res = H.c(loc,locy,factor);
1540 OPERATOR_IS_SET = true;
1541 }
1542 else
1543 {
1544 Res = H.template c<DN>(loc,locy);
1545 OPERATOR_IS_SET = true;
1546 }
1547 }
1548 // inverse photoemission
1549 else if (spec == "IPE" or spec == "IPEUP")
1550 {
1551 if constexpr (Symmetry::IS_SPIN_SU2()) // or spinless
1552 {
1553 Res = H.cdag(loc,locy,factor);
1554 OPERATOR_IS_SET = true;
1555 }
1556 else
1557 {
1558 Res = H.template cdag<UP>(loc,locy,factor);
1559 OPERATOR_IS_SET = true;
1560 }
1561 }
1562 else if (spec == "IPEDN")
1563 {
1564 if constexpr (Symmetry::IS_SPIN_SU2()) // or spinless
1565 {
1566 Res = H.cdag(loc,locy,factor);
1567 OPERATOR_IS_SET = true;
1568 }
1569 else
1570 {
1571 Res = H.template cdag<DN>(loc,locy,factor);
1572 OPERATOR_IS_SET = true;
1573 }
1574 }
1575 // charge structure factor
1576 else if (spec == "CSF" or spec == "ICSF")
1577 {
1578 if constexpr (!Symmetry::IS_CHARGE_SU2())
1579 {
1580 Res = H.n(loc,locy);
1581 OPERATOR_IS_SET = true;
1582 }
1583 else
1584 {
1585 throw;
1586 }
1587 }
1588 // Auger electron spectroscopy
1589 else if (spec == "AES")
1590 {
1591 if constexpr (!Symmetry::IS_CHARGE_SU2())
1592 {
1593 Res = H.cc(loc,locy);
1594 OPERATOR_IS_SET = true;
1595 }
1596 else
1597 {
1598 throw;
1599 }
1600 }
1601 // Appearance potential spectroscopy
1602 else if (spec == "APS")
1603 {
1604 if constexpr (!Symmetry::IS_CHARGE_SU2())
1605 {
1606 Res = H.cdagcdag(loc,locy);
1607 OPERATOR_IS_SET = true;
1608 }
1609 else
1610 {
1611 throw;
1612 }
1613 }
1614 // pseudospin structure factor
1615 else if (spec == "PSF")
1616 {
1617 if constexpr (Symmetry::IS_CHARGE_SU2())
1618 {
1619 Res = H.T(loc,locy);
1620 OPERATOR_IS_SET = true;
1621 }
1622 else
1623 {
1624 Res = H.Tp(loc,locy);
1625 OPERATOR_IS_SET = true;
1626 }
1627 }
1628 // pseudospin structure factor
1629 else if (spec == "PDAGSF")
1630 {
1631 if constexpr (Symmetry::IS_CHARGE_SU2())
1632 {
1633 Res = H.Tdag(loc,locy);
1634 OPERATOR_IS_SET = true;
1635 }
1636 else
1637 {
1638 Res = H.Tm(loc,locy);
1639 OPERATOR_IS_SET = true;
1640 }
1641 }
1642 // pseudospin structure factor: z-component
1643 else if (spec == "PSZ" or spec == "IPSZ")
1644 {
1645 if constexpr (!Symmetry::IS_CHARGE_SU2())
1646 {
1647 Res = H.Tz(loc,locy);
1648 OPERATOR_IS_SET = true;
1649 }
1650 else
1651 {
1652 throw;
1653 }
1654 }
1655 // hybridization structure factor
1656 else if (spec == "HSF")
1657 {
1658 if constexpr (Symmetry::IS_SPIN_SU2())
1659 {
1660 if (loc<H.length()-dLphys)
1661 {
1662 Res = H.cdagc(loc,loc+dLphys,0,0);
1663 OPERATOR_IS_SET = true;
1664 }
1665 else
1666 {
1667 lout << termcolor::yellow << "HSF operator hit right edge! Returning zero." << termcolor::reset << endl;
1668 Res = Hamiltonian::Zero(H.qPhys);
1669 OPERATOR_IS_SET = true;
1670 }
1671 }
1672 else
1673 {
1674 if (loc<H.length()-dLphys)
1675 {
1676 Res = H.template cdagc<UP,UP>(loc,loc+dLphys,0,0);
1677 OPERATOR_IS_SET = true;
1678 }
1679 else
1680 {
1681 lout << termcolor::yellow << "HSF operator hit right edge! Returning zero." << termcolor::reset << endl;
1682 Res = Hamiltonian::Zero(H.qPhys);
1683 OPERATOR_IS_SET = true;
1684 }
1685 }
1686 }
1687 // inverse hybridization structure factor
1688 else if (spec == "IHSF")
1689 {
1690 if constexpr (Symmetry::IS_SPIN_SU2())
1691 {
1692 if (loc<H.length()-dLphys)
1693 {
1694 Res = H.cdagc(loc+dLphys,loc,0,0);
1695 OPERATOR_IS_SET = true;
1696 }
1697 else
1698 {
1699 lout << termcolor::yellow << "IHSF operator hit right edge! Returning zero." << termcolor::reset << endl;
1700 Res = Hamiltonian::Zero(H.qPhys);
1701 OPERATOR_IS_SET = true;
1702 }
1703 }
1704 else
1705 {
1706 if (loc<H.length()-dLphys)
1707 {
1708 Res = H.template cdagc<UP,UP>(loc+dLphys,loc,0,0);
1709 OPERATOR_IS_SET = true;
1710 }
1711 else
1712 {
1713 lout << termcolor::yellow << "IHSF operator hit right edge! Returning zero." << termcolor::reset << endl;
1714 Res = Hamiltonian::Zero(H.qPhys);
1715 OPERATOR_IS_SET = true;
1716 }
1717 }
1718 }
1719 // hybridization triplet structure factor
1720 else if (spec == "HTS")
1721 {
1722 if constexpr (Symmetry::IS_SPIN_SU2())
1723 {
1724 if (loc<H.length()-dLphys)
1725 {
1726 Res = H.cdagc3(loc,loc+dLphys,0,0);
1727 OPERATOR_IS_SET = true;
1728 }
1729 else
1730 {
1731 lout << termcolor::yellow << "HTS operator hit right edge! Returning zero." << termcolor::reset << endl;
1732 Res = Hamiltonian::Zero(H.qPhys);
1733 OPERATOR_IS_SET = true;
1734 }
1735 }
1736 else
1737 {
1738 if (loc<H.length()-dLphys)
1739 {
1740 Res = H.template cdagc<UP,DN>(loc,loc+dLphys,0,0);
1741 OPERATOR_IS_SET = true;
1742 }
1743 else
1744 {
1745 lout << termcolor::yellow << "HTS operator hit right edge! Returning zero." << termcolor::reset << endl;
1746 Res = Hamiltonian::Zero(H.qPhys);
1747 OPERATOR_IS_SET = true;
1748 }
1749 }
1750 }
1751 // inverse hybridization triplet structure factor
1752 else if (spec == "IHTS")
1753 {
1754 if constexpr (Symmetry::IS_SPIN_SU2())
1755 {
1756 if (loc<H.length()-dLphys)
1757 {
1758 Res = H.cdagc3(loc+dLphys,loc,0,0);
1759 OPERATOR_IS_SET = true;
1760 }
1761 else
1762 {
1763 lout << termcolor::yellow << "IHTS operator hit right edge! Returning zero." << termcolor::reset << endl;
1764 Res = Hamiltonian::Zero(H.qPhys);
1765 OPERATOR_IS_SET = true;
1766 }
1767 }
1768 else
1769 {
1770 if (loc<H.length()-dLphys)
1771 {
1772 Res = H.template cdagc<DN,UP>(loc+dLphys,loc,0,0);
1773 OPERATOR_IS_SET = true;
1774 }
1775 else
1776 {
1777 lout << termcolor::yellow << "IHTS operator hit right edge! Returning zero." << termcolor::reset << endl;
1778 Res = Hamiltonian::Zero(H.qPhys);
1779 OPERATOR_IS_SET = true;
1780 }
1781 }
1782 }
1783 }
1784
1785 Res.set_locality(loc);
1786 assert(OPERATOR_IS_SET);
1787 return Res;
1788}
1789
1790#endif
DmrgLinearAlgebra.h
External functions to manipulate Mps and Mpo objects.
dot
Scalar dot(const Mps< Symmetry, Scalar > &Vbra, const Mps< Symmetry, Scalar > &Vket)
Definition: DmrgLinearAlgebra.h:24
sum
Hamiltonian sum(const Hamiltonian &H1, const Hamiltonian &H2, DMRG::VERBOSITY::OPTION VERBOSITY=DMRG::VERBOSITY::SILENT)
Definition: DmrgLinearAlgebra.h:866
OxV
void OxV(const Mpo< Symmetry, MpoScalar > &H, const Mpo< Symmetry, MpoScalar > &Hdag, const Mps< Symmetry, Scalar > &Vin, Mps< Symmetry, Scalar > &Vout, bool VERBOSE=true, double tol_compr=1e-4, int Mincr=100, int Mlimit=10000, int max_halfsweeps=100, int min_halfsweeps=1, bool MEASURE_DISTANCE=true)
Definition: DmrgLinearAlgebra.h:593
prod
Hamiltonian prod(const Hamiltonian &H1, const Hamiltonian &H2, const qarray< Hamiltonian::Symmetry::Nq > &Qtot=Hamiltonian::Symmetry::qvacuum(), DMRG::VERBOSITY::OPTION VERBOSITY=DMRG::VERBOSITY::SILENT)
Definition: DmrgLinearAlgebra.h:880
avg_hetero
Scalar avg_hetero(const Mps< Symmetry, Scalar > &Vbra, const Mpo< Symmetry, MpoScalar > &O, const Mps< Symmetry, Scalar > &Vket, bool USE_BOUNDARY=false, size_t usePower=1ul, const qarray< Symmetry::Nq > &Qmid=Symmetry::qvacuum())
Definition: DmrgLinearAlgebra.h:287
avg
Scalar avg(const Mps< Symmetry, Scalar > &Vbra, const Mpo< Symmetry, MpoScalar > &O, const Mps< Symmetry, Scalar > &Vket, size_t power_of_O=1ul, DMRG::DIRECTION::OPTION DIR=DMRG::DIRECTION::RIGHT, DMRG::VERBOSITY::OPTION CHOSEN_VERBOSITY=DMRG::VERBOSITY::SILENT)
Definition: DmrgLinearAlgebra.h:149
OxV_exact
void OxV_exact(const Mpo< Symmetry, MpoScalar > &O, const Mps< Symmetry, Scalar > &Vin, Mps< Symmetry, Scalar > &Vout, double tol_compr=1e-7, DMRG::VERBOSITY::OPTION VERBOSITY=DMRG::VERBOSITY::HALFSWEEPWISE, int max_halfsweeps=200, int min_halfsweeps=1, int Minit=-1, DMRG::BROOM::OPTION BROOMOPTION=DMRG::BROOM::QR)
Definition: DmrgLinearAlgebra.h:722
norm
double norm(const PivotMatrix0< Symmetry, Scalar, MpoScalar > &H)
Definition: DmrgPivotMatrix0.h:114
KONDO
@ KONDO
Definition: DmrgTypedefs.h:96
HUBBARD
@ HUBBARD
Definition: DmrgTypedefs.h:96
HEISENBERG
@ HEISENBERG
Definition: DmrgTypedefs.h:96
SZ
@ SZ
Definition: DmrgTypedefs.h:60
SP
@ SP
Definition: DmrgTypedefs.h:60
SM
@ SM
Definition: DmrgTypedefs.h:60
isReal
double isReal(double x)
Definition: DmrgTypedefs.h:21
Q_RANGE
Q_RANGE
Definition: GreenPropagator.h:27
ZERO_2PI
@ ZERO_2PI
Definition: GreenPropagator.h:27
GREEN_INTEGRATION
GREEN_INTEGRATION
Definition: GreenPropagator.h:28
OOURA
@ OOURA
Definition: GreenPropagator.h:28
DmrgJanitor::eps_truncWeight
double eps_truncWeight
Definition: DmrgJanitor.h:97
GreenPropagator::recalc_FTw
void recalc_FTw(double wmin_new, double wmax_new, int Nw_new=501)
Definition: GreenPropagator.h:2729
GreenPropagator::set_verbosity
void set_verbosity(DMRG::VERBOSITY::OPTION VERBOSITY)
Definition: GreenPropagator.h:459
GreenPropagator::recalc_FTwCell
void recalc_FTwCell(double wmin_new, double wmax_new, int Nw_new=501, bool TRANSPOSE=false)
Definition: GreenPropagator.h:2759
MpoTerms::scale
void scale(const Scalar factor, const Scalar offset=0., const std::size_t power=0ul, const double tolerance=1.e-14)
Definition: MpoTerms.h:2857
Mpo
Definition: Mpo.h:40
Mpo::Identity
static Mpo< Symmetry, Scalar > Identity(const std::vector< std::vector< qType > > &qPhys, const qType &Q=Symmetry::qvacuum())
Mpo::set_locality
void set_locality(std::size_t LocalSite_input)
Definition: Mpo.h:125
Mps
Definition: Mps.h:39
SpectralManager::compute_finiteCell
void compute_finiteCell(int Lcell, int x0, string wd, string label, double tmax, double dt=0.1, double wmin=-10., double wmax=10., int wpoints=501, Q_RANGE QR=ZERO_2PI, int qpoints=501, GREEN_INTEGRATION INT=OOURA, size_t Mlim=500ul, double tol_DeltaS=1e-2, double tol_compr=1e-4)
Definition: SpectralManager.h:1258
SpectralManager::PhiT
Mps< Symmetry, Scalar > PhiT
Definition: SpectralManager.h:173
SpectralManager::Symmetry
Hamiltonian::Symmetry Symmetry
Definition: SpectralManager.h:15
SpectralManager::OxPhiCellKet
vector< vector< Mps< Symmetry, complex< double > > > > OxPhiCellKet
Definition: SpectralManager.h:166
SpectralManager::energy
const double & energy() const
Definition: SpectralManager.h:75
SpectralManager::DAG
static string DAG(std::string spec)
Definition: SpectralManager.h:100
SpectralManager::Oshift
vector< vector< Scalar > > Oshift
Definition: SpectralManager.h:179
SpectralManager::TIME_DIR
static bool TIME_DIR(std::string spec)
Definition: SpectralManager.h:93
SpectralManager::CHECK_SPEC
static bool CHECK_SPEC(string spec)
Definition: SpectralManager.h:130
SpectralManager::OxPhiCellBra
vector< vector< Mps< Symmetry, complex< double > > > > OxPhiCellBra
Definition: SpectralManager.h:165
SpectralManager::make_A1P_finite
void make_A1P_finite(GreenPropagator< Hamiltonian, Symmetry, Scalar, complex< double > > &Gfull, string wd, string label, double tmax, double wmin=-10., double wmax=10., int wpoints=501, GREEN_INTEGRATION INT=OOURA)
Definition: SpectralManager.h:1371
SpectralManager::Hwork
Hamiltonian Hwork
Definition: SpectralManager.h:162
SpectralManager::get_PhiT
const Mps< Symmetry, Scalar > & get_PhiT() const
Definition: SpectralManager.h:76
SpectralManager::continue_beta_propagation
void continue_beta_propagation(const Hamiltonian &Hprop, int Lcell, int dLphys, double s_betainit, double betainit, double betamax, double dbeta, double tol_compr_beta, size_t Mlim, qarray< Hamiltonian::Symmetry::Nq > Q, double betaswitch, string wd, string th_label, string LOAD_BETA, bool SAVE_BETA, DMRG::VERBOSITY::OPTION VERB, vector< double > stateSavePoints, vector< string > stateSaveLabels, bool CALC_C, bool CALC_CHI)
Definition: SpectralManager.h:887
SpectralManager::g
Eigenstate< Umps< Symmetry, Scalar > > g
Definition: SpectralManager.h:168
SpectralManager::PhiTt
Mps< Symmetry, complex< double > > PhiTt
Definition: SpectralManager.h:174
SpectralManager::resize_Green
void resize_Green(string wd, string label, int Ns, double tmax, double dt, double wmin, double wmax, int wpoints, Q_RANGE QR, int qpoints, GREEN_INTEGRATION INT)
Definition: SpectralManager.h:1193
SpectralManager::Green
vector< GreenPropagator< Hamiltonian, typename Hamiltonian::Symmetry, Scalar, complex< double > > > Green
Definition: SpectralManager.h:170
SpectralManager::compute_thermal
void compute_thermal(string wd, string label, int dLphys, double tmax, double dt=0.1, double wmin=-10., double wmax=10., int wpoints=501, Q_RANGE QR=ZERO_2PI, int qpoints=501, GREEN_INTEGRATION INT=OOURA, size_t Mlim=500ul, double tol_DeltaS=1e-2, double tol_compr=1e-4)
Definition: SpectralManager.h:1294
SpectralManager::get_Op
Mpo< Symmetry, Scalar > get_Op(const Hamiltonian &H, size_t loc, std::string spec, double factor=1., size_t locy=0, int dLphys=1)
Definition: SpectralManager.h:1422
SpectralManager::gfinite
Eigenstate< Mps< Symmetry, Scalar > > gfinite
Definition: SpectralManager.h:169
SpectralManager::get_GwqCell
vector< vector< MatrixXcd > > get_GwqCell(int z) const
Definition: SpectralManager.h:146
SpectralManager::compute
void compute(string wd, string label, int Ns, double tmax, double dt=0.2, double wmin=-10., double wmax=10., int wpoints=501, Q_RANGE QR=ZERO_2PI, int qpoints=501, GREEN_INTEGRATION INT=OOURA, size_t Mlim=500ul, double tol_DeltaS=1e-2, double tol_compr=1e-4)
Definition: SpectralManager.h:1211
SpectralManager::OxPhiTt
vector< vector< Mps< Symmetry, complex< double > > > > OxPhiTt
Definition: SpectralManager.h:175
SpectralManager::apply_operators_on_thermal_state
void apply_operators_on_thermal_state(int Lcell, int dLphys, bool CHECK=true)
Definition: SpectralManager.h:1090
SpectralManager::make_A1P
void make_A1P(GreenPropagator< Hamiltonian, Symmetry, Scalar, complex< double > > &Gfull, string wd, string label, int Ns, double tmax, double wmin=-10., double wmax=10., int wpoints=501, Q_RANGE QR=ZERO_2PI, int qpoints=501, GREEN_INTEGRATION INT=OOURA, bool SAVE_N_MU=true)
Definition: SpectralManager.h:1329
SpectralManager::set_Ncells
void set_Ncells(int Ncells_input)
Definition: SpectralManager.h:148
SpectralManager::Phi
Mps< Symmetry, Scalar > Phi
Definition: SpectralManager.h:163
SpectralManager::specs
vector< string > specs
Definition: SpectralManager.h:160
SpectralManager::Odag
vector< vector< Mpo< typename Hamiltonian::Symmetry, Scalar > > > Odag
Definition: SpectralManager.h:176
SpectralManager::SpectralManager
SpectralManager(const vector< string > &specs_input, const Hamiltonian &H, DMRG::VERBOSITY::OPTION VERB=DMRG::VERBOSITY::HALFSWEEPWISE)
Definition: SpectralManager.h:32
SpectralManager::set_measurement
void set_measurement(int iz, string spec, double factor, int dLphys, qarray< Symmetry::Nq > Q, int Lcell, int measure_interval_input=10, string measure_name_input="M", string measure_subfolder_input=".", bool TRANSFORM=false)
Definition: SpectralManager.h:1178
SpectralManager::beta_propagation
void beta_propagation(const Hamiltonian &Hprop, const HamiltonianThermal &Htherm, int Lcell, int dLphys, double betamax_input, double dbeta_input, double tol_compr_beta_input, size_t Mlim, qarray< Symmetry::Nq > Q, double s_betainit, double betaswitch, string wd, string th_label, bool LOAD_BETA=false, bool SAVE_BETA=true, DMRG::VERBOSITY::OPTION VERB=DMRG::VERBOSITY::HALFSWEEPWISE, vector< double > stateSavePoints={}, vector< string > stateSaveLabels={}, int Ntaylor=0, bool CALC_C=true, bool CALC_CHI=true, bool USE_PHIT=false)
SpectralManager::compute_finite
void compute_finite(size_t j0, string wd, string label, int Ns, double tmax, double dt=0.1, double wmin=-10., double wmax=10., int wpoints=501, GREEN_INTEGRATION INT=OOURA, size_t Mlim=500ul, double tol_DeltaS=1e-2, double tol_compr=1e-4)
Definition: SpectralManager.h:1234
SpectralManager::reload
void reload(string wd, const vector< string > &specs_input, string label, int L, int Ncells, int Ns, double tmax, double wmin=-10., double wmax=10., int wpoints=501, Q_RANGE QR=ZERO_2PI, int qpoints=501, GREEN_INTEGRATION INT=OOURA)
Definition: SpectralManager.h:1394
SpectralManager::CHOSEN_VERB
DMRG::VERBOSITY::OPTION CHOSEN_VERB
Definition: SpectralManager.h:178
SpectralManager::Scalar
Hamiltonian::Mpo::Scalar_ Scalar
Definition: SpectralManager.h:14
SpectralManager::ground
const Umps< Symmetry, Scalar > & ground() const
Definition: SpectralManager.h:74
SpectralManager::set_PhiT
void set_PhiT(const Mps< Symmetry, Scalar > &PhiT_input)
Definition: SpectralManager.h:150
Umps
Definition: Umps.h:42
DMRG::CONTROL::DYN::max_alpha_rsvd
function< double(size_t)> max_alpha_rsvd
Definition: DmrgTypedefs.h:397
DMRG::CONTROL::DYN::iteration
function< DMRG::ITERATION::OPTION(size_t)> iteration
Definition: DmrgTypedefs.h:407
DMRG::CONTROL::GLOB::INITDIR
DMRG::DIRECTION::OPTION INITDIR
Definition: DmrgTypedefs.h:390
qarray
Definition: qarray.h:26