VMPS++
Loading...
Searching...
No Matches
NuclearManager.h
Go to the documentation of this file.
1#ifndef NUCLEAR_MANAGER
2#define NUCLEAR_MANAGER
3
4#include <boost/algorithm/string.hpp>
5
6#include <Eigen/SparseCore>
7#include <boost/rational.hpp>
8
9#include "solvers/DmrgSolver.h"
10#include "HDF5Interface.h"
11#include "ParamHandler.h"
12#include "DmrgLinearAlgebra.h"
13
14template<typename VectorType>
15struct jEigenstate
16{
17 Eigenstate<VectorType> eigenstate;
18 string label;
19 int index;
20};
21
22struct orbital
23{
24 orbital (const int &r)
25 {
26 l = r/2;
27 s = (r%2==0)? UP:DN;
28 }
29
30 bool operator == (const orbital &b) const
31 {
32 return (l==b.l and s==b.s)? true:false;
33 }
34
35 int l;
36 SPIN_INDEX s;
37};
38
39struct orbinfo
40{
41 bool operator == (const orbinfo &b) const
42 {
43 return (lev==b.lev and jx2==b.jx2 and label==b.label)? true:false;
44 }
45
46 std::array<int,4> lev;
47 std::array<int,4> jx2;
48 std::array<string,4> label;
49};
50
51double cgc (int j1x2, int j2x2, int Jx2, int m1x2, int m2x2, int Mx2)
52{
53 double Wigner3j = gsl_sf_coupling_3j(j1x2,j2x2,Jx2,m1x2,m2x2,-Mx2);
54// return pow(-1,-j1x2/2+j2x2/2-Mx2/2) * sqrt(Jx2+1.) * Wigner3j;
55 return pow(-1,+j1x2/2-j2x2/2+Mx2/2) * sqrt(Jx2+1.) * Wigner3j;
56}
57
58struct NuclearInfo
59{
60 string info() const
61 {
62 stringstream ss;
63 ss << "el=" << el << ", Zsingle=" << Zsingle << ", Nsingle=" << Nsingle << ", Nclosed=" << Nclosed << ", Nmax=" << Nmax << endl;
64 return ss.str();
65 }
66
67 string el;
68 int Nclosed;
69 int Nmax;
70 int Nsingle;
71 int Zsingle;
72 Nuclear::PARTICLE P;
73 string PARAM;
74 vector<string> levels;
75 string levelstring;
76 int Nlev;
77};
78
79void load_nuclearData (string el, string set, string rootfolder, NuclearInfo &info, double &V0, MatrixXd &Vij, vector<MatrixXd> &Vijnocgc, VectorXi &deg, VectorXd &eps0, bool &REF)
80{
81 REF = false;
82
83 ifstream DataLoader(rootfolder+"data_nuclear/sets.sav");
84 string line;
85 while (getline(DataLoader, line))
86 {
87 if (line[0] == '#') continue;
88
89 vector<string> split;
90 boost::split(split, line, boost::is_any_of("\t"));
91 if (split[0] == el and split[1] == set)
92 {
93 //info.Z = Nuclear::Ztable(split[0]);
94 info.el = split[0];
95 info.Nclosed = stoi(split[2]);
96 info.Nmax = stoi(split[3]);
97 info.Nsingle = stoi(split[4]);
98 info.Zsingle = stoi(split[5]);
99 info.P = static_cast<Nuclear::PARTICLE>(stoi(split[9]));
100 info.PARAM = split[10];
101 lout << info.info() << endl;
102
103 info.levelstring = split[6];
104 boost::split(info.levels, info.levelstring, boost::is_any_of(","));
105 info.Nlev = info.levels.size();
106 Nuclear::replace_halves(info.levelstring);
107 break;
108 }
109 }
110
111 if (el == "Sn" and set == "benchmark")
112 {
113 // benchmark Sn values
114 V0 = 1.;
115 info.Nclosed = 50;
116 info.Zsingle = 50;
117
118 info.Nlev = 5;
119 deg.resize(info.Nlev); deg << 4, 3, 2, 1, 6;
120
121 eps0.resize(info.Nlev); eps0 << -6.121, -5.508, -3.749, -3.891, -3.778;
122
123 Vij.resize(info.Nlev,info.Nlev);
124 Vij.setZero();
125 Vij(0,0) = 0.9850;
126 Vij(0,1) = 0.5711;
127 Vij(0,2) = 0.5184;
128 Vij(0,3) = 0.2920;
129 Vij(0,4) = 1.1454;
130 Vij(1,1) = 0.7063;
131 Vij(1,2) = 0.9056;
132 Vij(1,3) = 0.3456;
133 Vij(1,4) = 0.9546;
134 Vij(2,2) = 0.4063;
135 Vij(2,3) = 0.3515;
136 Vij(2,4) = 0.6102;
137 Vij(3,3) = 0.7244;
138 Vij(3,4) = 0.4265;
139 Vij(4,4) = 1.0599;
140
141 for (int i=0; i<info.Nlev; ++i)
142 for (int j=0; j<i; ++j)
143 {
144 Vij(i,j) = Vij(j,i);
145 }
146
147 int Jmax = 12;
148 Vijnocgc.resize(Jmax/2+1);
149 for (int J=0; J<=Jmax; J+=2)
150 {
151 Vijnocgc[J/2] = Vij;
152 }
153
154 for (int i=0; i<info.Nlev; ++i)
155 for (int j=0; j<info.Nlev; ++j)
156 {
157 Vij(i,j) /= sqrt(deg(i)*deg(j));
158 }
159
160 REF = true;
161 }
162 else if (el == "h11shell" and set == "benchmark")
163 {
164 V0 = 1.;
165 info.Nclosed = 0;
166 info.Zsingle = 0;
167
168 info.Nlev = 1;
169 deg.resize(info.Nlev); deg << 6;
170
171 eps0.resize(info.Nlev); eps0 << 0.;
172
173 Vij.resize(info.Nlev,info.Nlev);
174 Vij.setConstant(0.25);
175
176 int Jmax = 10;
177 Vijnocgc.resize(Jmax/2+1);
178 for (int J=0; J<=Jmax; J+=2)
179 {
180 Vijnocgc[J/2] = Vij*6.;
181 }
182 }
183 else if (el == "Richardson50" and set == "benchmark")
184 {
185 V0 = 3.;
186 info.Nclosed = 0;
187 info.Zsingle = 0;
188
189 info.Nlev = 50;
190 deg.resize(info.Nlev); deg.setConstant(1);
191
192 eps0.resize(info.Nlev);
193 for (int j=0; j<info.Nlev; ++j)
194 {
195 eps0[j] = j;
196 }
197
198 Vij.resize(info.Nlev,info.Nlev);
199 Vij.setConstant(1.);
200
201 int Jmax = 0;
202 Vijnocgc.resize(Jmax/2+1);
203 for (int J=0; J<=Jmax; J+=2)
204 {
205 Vijnocgc[J/2] = Vij;
206 }
207 }
208 // From: Zelevinsky, Volya: Physics of Atomic Nuclei (2017), pp. 218-220
209 else if (el == "f7shell" and set == "benchmark")
210 {
211 V0 = 2.535;
212 std::array<double,4> V_J = {1., 0.981/V0, -0.140/V0, -0.664/V0};
213 info.Nclosed = 0;
214 info.Zsingle = 0;
215
216 info.Nlev = 1;
217 deg.resize(info.Nlev); deg << 4;
218
219 eps0.resize(info.Nlev); eps0 << -9.626;
220
221 Vij.resize(info.Nlev,info.Nlev);
222 Vij.setConstant(1.);
223
224 int Jmax = 6;
225 Vijnocgc.resize(Jmax/2+1);
226 for (int J=0; J<=Jmax; J+=2)
227 {
228 Vijnocgc[J/2] = Vij * V_J[J/2];
229 }
230 }
231 else if (el == "testshell" and set == "benchmark")
232 {
233 V0 = 1.;
234 std::array<double,3> V_J = {1., 0.5, 0.25};
235 info.Nclosed = 0;
236 info.Zsingle = 0;
237
238 info.Nlev = 3;
239 deg.resize(info.Nlev); deg << 1, 2, 3;
240
241 eps0.resize(info.Nlev); eps0 << 0., 0., 0.;
242
243 double lower_bound = 0.;
244 double upper_bound = 1.;
245 std::uniform_real_distribution<double> unif(lower_bound,upper_bound);
246 std::default_random_engine re;
247
248 Vij.resize(info.Nlev,info.Nlev);
249 Vij.setRandom();
250 for (int i=0; i<info.Nlev; ++i)
251 for (int j=0; j<=i; ++j)
252 {
253 Vij(i,j) = unif(re);
254 Vij(j,i) = Vij(i,j);
255 }
256
257 int Jmax = 4;
258 Vijnocgc.resize(Jmax/2+1);
259 for (int J=0; J<=Jmax; J+=2)
260 {
261 Vijnocgc[J/2] = Vij * V_J[J/2];
262 }
263 }
264
265 if (set != "benchmark")
266 {
267 HDF5Interface target(make_string(rootfolder+"NuclearLevels",info.PARAM,"/hdf5/NuclearData_Z=",info.Zsingle,"_N=",info.Nsingle,"_P=",info.P,"_j=",info.levelstring,".h5"), READ);
268 target.load_matrix(Vij,"Vij","");
269 target.load_vector(eps0,"eps0","");
270 target.load_vector(deg,"deg","");
271
272 int Jmax;
273 target.load_scalar(Jmax,"Jmax","");
274 Vijnocgc.resize(Jmax/2+1);
275 for (int J=0; J<=Jmax; J+=2)
276 {
277 target.load_matrix(Vijnocgc[J/2],"Vijnocgc",make_string("J=",J));
278 }
279
280// if (RANDOM)
281// {
282// static thread_local mt19937 generatorUniformReal(random_device{}());
283// uniform_real_distribution<double> distribution(0.0015, 0.015);
284//
285// for (int i=0; i<Vij.rows(); ++i)
286// for (int j=0; j<=i; ++j)
287// {
288// Vij(i,j) = distribution(generatorUniformReal);
289// Vij(j,i) = Vij(i,j);
290// }
291// lout << "Vij random=" << endl << Vij << endl;
292// }
293 }
294}
295
296template<typename MODEL>
297class NuclearManager
298{
299public:
300
301 NuclearManager(){};
302
303 NuclearManager (int Nclosed_input, int Nsingle_input, int Z_input, const ArrayXi &deg_input, const vector<string> &labels_input,
304 const ArrayXd &eps0_input, const ArrayXXd &Vij_input, double V0_input=1.,
305 bool REF_input=false, string PARAM_input="Seminole")
306 :Nclosed(Nclosed_input), Nsingle(Nsingle_input), Z(Z_input), deg(deg_input), labels(labels_input), eps0(eps0_input), Vij(Vij_input), V0(V0_input), REF(REF_input), PARAM(PARAM_input)
307 {
308 Nlev = deg.rows();
309 L = deg.sum();
310 lout << "Nlev=" << Nlev << ", L=" << L << endl;
311 Vij.triangularView<Lower>() = Vij.transpose();
312 levelstring = str(labels_input);
313
314 construct();
315 };
316
317 NuclearManager (int Nclosed_input, int Nsingle_input, int Z_input,const string &filename, double V0_input, bool REF_input=false)
318 :Nclosed(Nclosed_input), Nsingle(Nsingle_input), Z(Z_input), V0(V0_input), REF(REF_input)
319 {
320 HDF5Interface source(filename+".h5", READ);
321
322 source.load_vector<int>(deg, "deg", "");
323 source.load_vector<double>(eps0, "eps0", "");
324 source.load_matrix<double>(Vij, "Vij", "");
325
326 Nlev = deg.rows();
327 L = deg.sum();
328 lout << "Nlev=" << Nlev << ", L=" << L << endl;
329 Vij.triangularView<Lower>() = Vij.transpose();
330
331 labels.resize(Nlev);
332 for (int i=0; i<Nlev; ++i)
333 {
334 source.load_char(labels[i], make_string("lev",i));
335 }
336 levelstring = str(labels);
337
338 construct();
339 };
340
341 void construct();
342
343 void make_Hamiltonian (bool LOAD=false, bool SAVE=false, string wd="./", bool ODD=true);
344 void make_fullHamiltonian (bool LOAD=false, bool SAVE=false, string wd="./");
345
346 tuple<Eigenstate<typename MODEL::StateXd>,string,int>
347 calc_gs (int Nshell, LANCZOS::EDGE::OPTION EDGE=LANCZOS::EDGE::GROUND, bool CALC_VAR=true, DMRG::VERBOSITY::OPTION VERB = DMRG::VERBOSITY::ON_EXIT) const;
348
349 Eigenstate<typename MODEL::StateXd>
350 calc_gs_full (int Nshell, LANCZOS::EDGE::OPTION EDGE=LANCZOS::EDGE::GROUND, bool CALC_VAR=true, DMRG::VERBOSITY::OPTION VERB = DMRG::VERBOSITY::ON_EXIT) const;
351
352 void compute (bool LOAD_MPO=false, bool SAVE_MPO=false, string wd="./", int minNshell=0, int maxNshell=-1, bool SAVE=true);
353 void compute_parallel (bool LOADv=false, bool SAVE_MPO=false, string wd="./", int minNshell=0, int maxNshell=-1, bool SAVE=true);
354
355 ArrayXd get_occ() const;
356 inline ArrayXd get_eps0() const {return eps0;}
357 inline ArrayXi get_deg() const {return deg;}
358 inline ArrayXi get_offset() const {return offset;}
359 inline size_t length() const {return L;}
360 inline ArrayXd get_sign() {return sign;}
361 inline ArrayXd get_levelsign (int level) {return sign_per_level[level];}
362 inline VectorXd get_onsite_free() const {return onsite_free;}
363 inline MODEL get_H() const {return H;}
364 inline MODEL get_Hfull() const {return Hfull;}
365
366 double Sn_Eref (int N) const;
367 ArrayXd Sn_nref (int Nshell) const;
368 ArrayXd Sn_Nref (int Nshell) const;
369 ArrayXd Sn_Pref (int Nshell) const;
370 ArrayXd Sn_Sref14() const;
371
372 const Eigenstate<typename MODEL::StateXd> &get_g (int Nshell) const {return g[Nshell];}
373 MatrixXd get_Ghop() const {return Ghop;}
374 VectorXd get_onsite() const {return onsite;}
375 VectorXd get_Gloc() const {return Gloc;}
376 VectorXi get_mfactor() const {return mfactor;}
377 std::array<MatrixXd,2> get_FlipChart() const {return FlipChart;}
378 ArrayXXd get_sign2d() const {return sign2d;}
379
380 vector<tuple<orbital,orbital,orbital,orbital,double,orbinfo>> generate_couplingList (int J) const;
381
382private:
383
384 void save (string wd, int Nfrst=0, int Nlast=-1, int dNshell=1) const;
385
386 bool REF = false;
387 string PARAM = "Seminole";
388
389 int Nclosed, Nsingle, Z;
390 int Nlev;
391 size_t L;
392
393 string levelstring;
394 vector<string> labels;
395 VectorXi deg, offset, mfactor;
396 VectorXd eps0;
397 ArrayXd sign;
398 std::array<MatrixXd,2> FlipChart;
399 vector<boost::rational<int>> mlist, jlist;
400 vector<int> orblist;
401 ArrayXXd sign2d;
402
403 vector<ArrayXd> sign_per_level;
404
405 MatrixXd Vij, Ghop;
406 double V0;
407 VectorXd onsite, Gloc, onsite_free;
408 VectorXd energies, energies_free;
409
410 string outfile;
411
412 vector<Eigenstate<typename MODEL::StateXd>> g;
413 vector<double> var;
414 vector<int> Mmax;
415 vector<VectorXd> n, avgN;
416 vector<string> Jlabel;
417 vector<int> Jindex;
418
419 MODEL H;
420 vector<MODEL> Hodd;
421 MODEL Hfull;
422
423 vector<MODEL> Hperturb;
424 vector<double> perturb_sweeps;
425 vector<Mpo<typename MODEL::Symmetry>> make_cdagj (const MODEL &H_input) const;
426 Mpo<typename MODEL::Symmetry> make_cdagj (int j, const MODEL &H_input) const;
427 vector<Mpo<typename MODEL::Symmetry>> cdagj;
428};
429
430template<typename MODEL>
431void NuclearManager<MODEL>::
432construct()
433{
434 // offsets
435 offset.resize(Nlev);
436 offset(0) = 0;
437 for (int i=1; i<Nlev; ++i)
438 {
439 offset(i) = deg.head(i).sum();
440 }
441 lout << "offset=" << offset.transpose() << endl;
442
443 // G hopping
444 Ghop.resize(L,L); Ghop.setZero();
445 for (int i=0; i<Nlev; ++i)
446 for (int j=0; j<Nlev; ++j)
447 {
448 Ghop.block(offset(i),offset(j),deg(i),deg(j)).setConstant(-V0*Vij(i,j));
449 }
450
451 if constexpr (MODEL::FAMILY == HUBBARD)
452 {
453 for (int i=0; i<L; ++i)
454 for (int j=0; j<L; ++j)
455 {
456 Ghop(i,j) *= pow(-1,i+j); // cancels sign, so that one can use Vxy
457 }
458
459 // Convention: start with highest projection:
460 // ...5/2,3/2,1/2
461 // Careful when changing!
462 sign2d.resize(L,L);
463 sign2d.setZero();
464 for (int j1=0; j1<Nlev; ++j1)
465 for (int j2=0; j2<Nlev; ++j2)
466 {
467 boost::rational<int> J1 = deg(j1)-boost::rational<int>(1,2);
468 boost::rational<int> J2 = deg(j2)-boost::rational<int>(1,2);
469
470 ArrayXXd signblock(deg(j1),deg(j2));
471
472 for (int m1=0; m1<deg(j1); ++m1)
473 for (int m2=0; m2<deg(j2); ++m2)
474 {
475 boost::rational<int> M1 = J1-m1;
476 boost::rational<int> M2 = J2-m2;
477 assert((J1+J2+M1+M2).denominator() == 1);
478 signblock(m1,m2) = pow(-1,(J1+J2-M1-M2).numerator());
479 }
480
481 sign2d.block(offset(j1),offset(j2),deg(j1),deg(j2)) = signblock;
482 }
483
484 sign.resize(L);
485 sign_per_level.resize(Nlev);
486 sign.setZero();
487 for (int j=0; j<Nlev; ++j)
488 {
489 boost::rational<int> J = deg(j)-boost::rational<int>(1,2);
490
491 ArrayXd signblock(deg(j));
492
493 for (int m=0; m<deg(j); ++m)
494 {
495 boost::rational<int> M = J-m;
496 signblock(m) = pow(-1,(J-M).numerator());
497 }
498
499 sign_per_level[j] = signblock;
500 lout << signblock.transpose() << endl;
501
502 sign.segment(offset(j),deg(j)) = signblock;
503 }
504// lout << "sign2d=" << endl << sign2d << endl;
505
506 Ghop = (Ghop.array()*sign2d).matrix();
507 // lout << "Ghop=" << endl << Ghop << endl;
508 }
509
510 // G local
511 Gloc = Ghop.diagonal();
512 Ghop.diagonal().setZero();
513
514 // onsite
515 onsite.resize(L); onsite.setZero();
516 for (int i=0; i<Nlev; ++i)
517 {
518 onsite.segment(offset(i),deg(i)).setConstant(eps0(i));
519 }
520
521 onsite_free.resize(2*L);
522 for (int i=0; i<L; ++i)
523 {
524 onsite_free(2*i) = onsite(i);
525 onsite_free(2*i+1) = onsite(i);
526 }
527
528 outfile = make_string("PairingResult_Z=",Z,"_Nclosed=",Nclosed,"_Nsingle=",Nsingle,"_V0=",V0);
529 outfile += make_string("_j=",levelstring);
530 if (PARAM!="Seminole") outfile += make_string("_PARAM=",PARAM);
531 outfile += ".h5";
532 lout << "outfile=" << outfile << endl;
533
534 /*PermutationMatrix<Dynamic,Dynamic> P(L);
535 P.setIdentity();
536 std::random_device rd;
537 std::mt19937 g(rd());
538 std::shuffle(P.indices().data(), P.indices().data()+P.indices().size(), g);
539 Ghop = P.inverse()*Ghop*P;
540 onsite = P.inverse()*onsite;
541 Gloc = P.inverse()*Gloc;*/
542
543 mfactor.resize(L); mfactor.setZero();
544 int l = 0;
545 for (int j=0; j<Nlev; ++j)
546 {
547 boost::rational<int> jval = deg(j)-boost::rational<int>(1,2);
548
549 // Convention: start with highest projection:
550 // ...5/2,3/2,1/2
551 // Careful when changing!
552 for (boost::rational<int> mval=jval; mval>0; --mval)
553 {
554 mfactor(l) = (boost::rational<int>(2,1)*mval).numerator();
555 mlist.push_back(+mval);
556 mlist.push_back(-mval);
557 jlist.push_back(jval);
558 jlist.push_back(jval);
559 orblist.push_back(j);
560 orblist.push_back(j);
561 ++l;
562 }
563 }
564// lout << "mfactor=" << mfactor.transpose() << endl;
565// for (int i=0; i<mlist.size(); ++i)
566// {
567// lout << "m=" << mlist[i] << ", j=" << jlist[i] << ", orb=" << orblist[i] << endl;
568// }
569
570 FlipChart[UP].resize(2*L,2*L); FlipChart[UP].setZero();
571 FlipChart[DN].resize(2*L,2*L); FlipChart[DN].setZero();
572 for (int i=0; i<mlist.size(); ++i)
573 for (int j=0; j<mlist.size(); ++j)
574 {
575 double jval = boost::rational_cast<double>(jlist[i]);
576 double m1 = boost::rational_cast<double>(mlist[i]);
577 double m2 = boost::rational_cast<double>(mlist[j]);
578
579 if (mlist[j] == mlist[i]+boost::rational<int>(1,1) and jlist[i]==jlist[j])
580 {
581 FlipChart[DN](i,j) = sqrt(jval*(jval+1.)-m1*m2);
582 }
583 if (mlist[j] == mlist[i]-boost::rational<int>(1,1) and jlist[i]==jlist[j])
584 {
585 FlipChart[UP](i,j) = sqrt(jval*(jval+1.)-m1*m2);
586 }
587 }
588// cout << "UP=" << endl << FlipChart[UP] << endl;
589// cout << "DN=" << endl << FlipChart[DN] << endl;
590}
591
592template<typename MODEL>
593ArrayXd NuclearManager<MODEL>::
594get_occ() const
595{
596 ArrayXd res(Nlev);
597 for (int i=0; i<deg.rows(); ++i)
598 {
599 res(i) = 2.*deg(i);
600 }
601 return res;
602}
603
604template<typename MODEL>
605double NuclearManager<MODEL>::
606Sn_Eref (int Nshell) const
607{
608 double res = std::nan("0");
609 /*if (Nshell == 12)
610 {
611 res = -75.831;
612 }
613 else if (Nshell == 14)
614 {
615 res = -86.308;
616 }
617 else if (Nshell == 16)
618 {
619 res = -95.942;
620 }*/
621 vector<double> EDenergies =
622 {
623 -0, //0
624 -6.121, //1
625 -14.0158183065649, //2
626 -19.81090128197716, //3
627 -27.49864249287803, //4
628 -32.94931684903962, //5
629 -40.44176691839276, //6
630 -45.5231949893755, //7
631 -52.83210266861026, //8
632 -57.52194523311658, //9
633 -64.64401299541129, //10
634 -69.01055649707473, //11
635 -75.83057632660048, //12
636 -79.75643390536675, //13
637 -86.30836166475844, //14
638 -89.69722801960962, //15
639 -95.94172568796624, //16
640 -99.06316035101028, //17 -99.0631603509355 // 10 digits
641 -104.9348539406987, //18
642 -107.8672420394976, //19
643 -113.3725141279236, //20
644 -116.108976116819, //21
645 -121.3002558685274, //22
646 -123.8366779467813, //23
647 -128.7462271575389, //24
648 -131.0794704088455, //25
649 -135.7291949708218, //26
650 -137.8562329243708, //27
651 -142.2621883681673, //28
652 -144.1794965107825, //29
653 -148.3543512502746, //30
654 -150.0597499999999, //31
655 -154.0119 //32
656 };
657 return EDenergies[Nshell];
658}
659
660template<typename MODEL>
661ArrayXd NuclearManager<MODEL>::
662Sn_nref (int Nshell) const
663{
664 ArrayXd res = Sn_Nref(Nshell);
665 for (int i=0; i<5; ++i) res(i) /= deg(i)*2;
666 /*res(0) = 0.Z70;
667 res(1) = 0.744;
668 res(2) = 0.157;
669 res(3) = 0.178;
670 res(4) = 0.133;*/
671 return res;
672}
673
674template<typename MODEL>
675ArrayXd NuclearManager<MODEL>::
676Sn_Nref (int Nshell) const
677{
678 ArrayXd res(5);
679 if (Nshell == 12)
680 {
681 res(0) = 6.4551;
682 res(1) = 3.6458;
683 res(2) = 0.4757;
684 res(3) = 0.2358;
685 res(4) = 1.1877;
686 }
687 else if (Nshell == 14)
688 {
689 res(0) = 6.9562;
690 res(1) = 4.4630;
691 res(2) = 0.6265;
692 res(3) = 0.3558;
693 res(4) = 1.5985;
694 }
695 else if (Nshell == 16)
696 {
697 res(0) = 7.1413;
698 res(1) = 4.7697;
699 res(2) = 0.9280;
700 res(3) = 0.6463;
701 res(4) = 2.5147;
702 }
703 else if (Nshell == 18)
704 {
705 res(0) = 7.2727;
706 res(1) = 4.9743;
707 res(2) = 1.2526;
708 res(3) = 0.9171;
709 res(4) = 3.5833;
710 }
711 return res;
712}
713
714template<typename MODEL>
715ArrayXd NuclearManager<MODEL>::
716Sn_Pref (int Nshell) const
717{
718 ArrayXd res(5);
719 if (Nshell == 12)
720 {
721 res(0) = 1.8870;
722 res(1) = 1.7040;
723 res(2) = 0.6569;
724 res(3) = 0.3287;
725 res(4) = 1.8089;
726 }
727 else if (Nshell == 14)
728 {
729 res(0) = 1.6197;
730 res(1) = 1.6107;
731 res(2) = 0.7411;
732 res(3) = 0.3958;
733 res(4) = 2.0687;
734 }
735 else if (Nshell == 16)
736 {
737 res(0) = 1.3592;
738 res(1) = 1.3494;
739 res(2) = 0.8722;
740 res(3) = 0.5138;
741 res(4) = 2.5256;
742 }
743 else if (Nshell == 18)
744 {
745 res(0) = 1.2466;
746 res(1) = 1.2336;
747 res(2) = 0.9720;
748 res(3) = 0.5563;
749 res(4) = 2.8951;
750 }
751 for (int i=0; i<5; ++i) res(i) /= sqrt(deg(i));
752 /*res(0) = 0.81;
753 res(1) = 0.93;
754 res(2) = 0.524;
755 res(3) = 0.396;
756 res(4) = 0.845;*/
757 return res;
758}
759
760template<typename MODEL>
761ArrayXd NuclearManager<MODEL>::
762Sn_Sref14() const
763{
764 ArrayXd res(5);
765 res(0) = 6.86;
766 res(1) = 6.55;
767 res(2) = 7.25;
768 res(3) = 6.98;
769 res(4) = 7.12;
770 return res;
771}
772
773template<typename MODEL>
774void NuclearManager<MODEL>::
775make_Hamiltonian (bool LOAD, bool SAVE, string wd, bool ODD)
776{
777 string filename = make_string(wd,"H_Z=",Z,"_Nclosed=",Nclosed,"_Nsingle=",Nsingle,"_V0=",V0);
778 filename += make_string("_j=",levelstring);
779 if (PARAM!="Seminole") filename += make_string("_PARAM=",PARAM);
780
781 if (LOAD)
782 {
783 MODEL Hres(L,{{"maxPower",1ul}});
784 Hres.load(filename);
785 H = Hres;
786 lout << H.info() << endl;
787 }
788 else
789 {
790 ArrayXXd Ghopx2 = 2.*Ghop; // 2 compensates 0.5 in our definition
791
792 vector<Param> params_loc;
793 params_loc.push_back({"t",0.});
794 params_loc.push_back({"J",0.});
795 params_loc.push_back({"maxPower",1ul});
796 if constexpr (MODEL::FAMILY == HUBBARD) params_loc.push_back({"REMOVE_SINGLE",true});
797
798 for (size_t l=0; l<L; ++l)
799 {
800 if constexpr (MODEL::FAMILY == HUBBARD)
801 {
802 params_loc.push_back({"t0",onsite(l),l});
803 params_loc.push_back({"U",Gloc(l),l});
804 params_loc.push_back({"mfactor",mfactor(l),l});
805 }
806 else if constexpr (MODEL::FAMILY == HEISENBERG)
807 {
808 params_loc.push_back({"nu",2.*onsite(l)+Gloc(l),l});
809 }
810 }
811
812 MODEL Hloc(L,params_loc,BC::OPEN,DMRG::VERBOSITY::SILENT);
813 MpoTerms<typename MODEL::Symmetry,double> Terms = Hloc;
814
815 Stopwatch<> CompressionTimer;
816 for (int i=0; i<Ghopx2.rows(); ++i)
817 {
818// cout << "compressing: site " << i+1 << "/" << Ghopx2.rows() << endl;
819 vector<Param> params_tmp;
820 params_tmp.push_back({"t",0.});
821 params_tmp.push_back({"J",0.});
822 params_tmp.push_back({"maxPower",1ul});
823 if constexpr (MODEL::FAMILY == HUBBARD) params_tmp.push_back({"REMOVE_SINGLE",true});
824
825 ArrayXXd Ghopx2row(Ghop.rows(),Ghop.cols());
826 Ghopx2row.setZero();
827 Ghopx2row.matrix().row(i) = Ghopx2.matrix().row(i);
828 if constexpr (MODEL::FAMILY == HUBBARD)
829 {
830 params_tmp.push_back({"Vxyfull",Ghopx2row});
831 for (size_t l=0; l<L; ++l) params_tmp.push_back({"mfactor",mfactor(l),l});
832 }
833 else if constexpr (MODEL::FAMILY == HEISENBERG)
834 {
835 params_tmp.push_back({"Jxyfull",Ghopx2row});
836 }
837
838 MODEL Hterm = MODEL(L, params_tmp, BC::OPEN, DMRG::VERBOSITY::SILENT);
839 Terms.set_verbosity(DMRG::VERBOSITY::SILENT);
840 Hterm.set_verbosity(DMRG::VERBOSITY::SILENT);
841 Terms = MODEL::sum(Terms,Hterm);
842 }
843
844 vector<Param> params_full = params_loc;
845 if constexpr (MODEL::FAMILY == HUBBARD)
846 {
847 params_full.push_back({"Vxyfull",Ghopx2});
848 }
849 else if constexpr (MODEL::FAMILY == HEISENBERG)
850 {
851 params_full.push_back({"Jxyfull",Ghopx2});
852 }
853
854// MODEL Hfull(L,params_full);
855// lout << "Hfull:" << Hfull.info() << endl;
856
857 Mpo<typename MODEL::Symmetry,double> Hmpo(Terms);
858 H = MODEL(Hmpo,params_full);
859// H.calc(2ul); // takes too long
860 lout << H.info() << endl;
861 lout << CompressionTimer.info("MPO compression even") << endl;
862
863 if (SAVE) H.save(filename);
864
865 // ODD
866 if constexpr (MODEL::FAMILY == HUBBARD)
867 {
868 if (ODD)
869 {
870 Hodd.resize(Nlev);
871
872 Stopwatch<> CompressionTimerOdd;
873 #pragma omp parallel for
874 for (int j=0; j<Nlev; ++j)
875 {
876 vector<Param> params_loc_odd;
877 params_loc_odd.push_back({"t",0.});
878 params_loc_odd.push_back({"J",0.});
879 params_loc_odd.push_back({"maxPower",1ul});
880 for (size_t l=0; l<L; ++l)
881 {
882 if (l<offset(j) or l>=offset(j)+deg(j)) params_loc_odd.push_back({"REMOVE_SINGLE",true,l});
883 else params_loc_odd.push_back({"REMOVE_SINGLE",false,l});
884
885 params_loc_odd.push_back({"t0",onsite(l),l});
886 params_loc_odd.push_back({"U",Gloc(l),l});
887 params_loc_odd.push_back({"mfactor",mfactor(l),l});
888 }
889
890 MODEL Hloc_odd(L,params_loc_odd,BC::OPEN,DMRG::VERBOSITY::SILENT);
891 MpoTerms<typename MODEL::Symmetry,double> Terms_odd = Hloc_odd;
892
893 for (int i=0; i<Ghopx2.rows(); ++i)
894 {
895 vector<Param> params_tmp_odd;
896 params_tmp_odd.push_back({"t",0.});
897 params_tmp_odd.push_back({"J",0.});
898 params_tmp_odd.push_back({"maxPower",1ul});
899 for (size_t l=0; l<L; ++l)
900 {
901 if (l<offset(j) or l>=offset(j)+deg(j)) params_tmp_odd.push_back({"REMOVE_SINGLE",true,l});
902 else params_tmp_odd.push_back({"REMOVE_SINGLE",false,l});
903
904 params_tmp_odd.push_back({"mfactor",mfactor(l),l});
905 }
906
907 ArrayXXd Ghopx2row(Ghop.rows(),Ghop.cols());
908 Ghopx2row.setZero();
909 Ghopx2row.matrix().row(i) = Ghopx2.row(i);
910 params_tmp_odd.push_back({"Vxyfull",Ghopx2row});
911
912 MODEL Hterm = MODEL(L, params_tmp_odd, BC::OPEN, DMRG::VERBOSITY::SILENT);
913 Terms_odd.set_verbosity(DMRG::VERBOSITY::SILENT);
914 Hterm.set_verbosity(DMRG::VERBOSITY::SILENT);
915 Terms_odd = MODEL::sum(Terms_odd,Hterm);
916 }
917
918 vector<Param> params_full_odd = params_loc_odd;
919 params_loc_odd.push_back({"Vxyfull",Ghopx2});
920
921 Mpo<typename MODEL::Symmetry,double> Hmpo(Terms_odd);
922 Hodd[j] = MODEL(Hmpo,params_full_odd);
923 }
924 lout << CompressionTimerOdd.info("MPO compression odd") << endl;
925 }
926 }
927 }
928
929 /*if constexpr (MODEL::FAMILY == HUBBARD)
930 {
931 // good:
932 vector<double> tvals = {1., 0.5, 0.25, 0.125, 0.05};
933 perturb_sweeps = {4, 4, 4, 4, 4};
934 // bad:
935// vector<double> tvals = {4., 3., 2., 1., 0.5, 0.25, 0.125, 0.05, 0.01};
936// perturb_sweeps = {2, 2, 2, 2, 2, 2, 2, 2, 2};
937 Hperturb.resize(tvals.size());
938 for (int i=0; i<tvals.size(); ++i)
939 {
940 MODEL Hperturb(L,{{"t",tvals[i]},{"Bx",tvals[i]},{"maxPower",1ul}}, BC::OPEN, DMRG::VERBOSITY::SILENT);
941 Hperturb = sum<MODEL,double>(H, Hperturb, DMRG::VERBOSITY::SILENT);
942 Hperturb[i] = Hperturb;
943 }
944 }*/
945
946// cdagj = make_cdagj(H);
947}
948
949template<typename MODEL>
950void NuclearManager<MODEL>::
951make_fullHamiltonian (bool LOAD, bool SAVE, string wd)
952{
953 string filename = make_string(wd,"H_Z=",Z,"_Nclosed=",Nclosed,"_Nsingle=",Nsingle,"_V0=",V0);
954 filename += make_string("_j=",levelstring);
955 if (PARAM!="Seminole") filename += make_string("_PARAM=",PARAM);
956
957 if (LOAD)
958 {
959 MODEL Hres(L,{{"maxPower",1ul}});
960 Hres.load(filename);
961 H = Hres;
962 lout << H.info() << endl;
963 }
964 else
965 {
966 ArrayXXd Ghopx2 = 2.*Ghop; // 2 compensates 0.5 in our definition
967
968 vector<Param> params_loc;
969 params_loc.push_back({"t",0.});
970 params_loc.push_back({"J",0.});
971 params_loc.push_back({"maxPower",1ul});
972
973 for (size_t l=0; l<L; ++l)
974 {
975 if constexpr (MODEL::FAMILY == HUBBARD)
976 {
977 params_loc.push_back({"t0",onsite(l),l});
978 params_loc.push_back({"U",Gloc(l),l});
979 params_loc.push_back({"mfactor",mfactor(l),l});
980 }
981 else if constexpr (MODEL::FAMILY == HEISENBERG)
982 {
983 params_loc.push_back({"nu",2.*onsite(l)+Gloc(l),l});
984 }
985 }
986
987 MODEL Hloc(L,params_loc,BC::OPEN,DMRG::VERBOSITY::SILENT);
988 MpoTerms<typename MODEL::Symmetry,double> Terms = Hloc;
989
990 Stopwatch<> CompressionTimer;
991 for (int i=0; i<Ghopx2.rows(); ++i)
992 {
993// cout << "compressing: site " << i+1 << "/" << Ghopx2.rows() << endl;
994 vector<Param> params_tmp;
995 params_tmp.push_back({"t",0.});
996 params_tmp.push_back({"J",0.});
997 params_tmp.push_back({"maxPower",1ul});
998
999 ArrayXXd Ghopx2row(Ghop.rows(),Ghop.cols());
1000 Ghopx2row.setZero();
1001 Ghopx2row.matrix().row(i) = Ghopx2.matrix().row(i);
1002 if constexpr (MODEL::FAMILY == HUBBARD)
1003 {
1004 params_tmp.push_back({"Vxyfull",Ghopx2row});
1005 for (size_t l=0; l<L; ++l) params_tmp.push_back({"mfactor",mfactor(l),l});
1006 }
1007 else if constexpr (MODEL::FAMILY == HEISENBERG)
1008 {
1009 params_tmp.push_back({"Jxyfull",Ghopx2row});
1010 }
1011
1012 MODEL Hterm = MODEL(L, params_tmp, BC::OPEN, DMRG::VERBOSITY::SILENT);
1013 Terms.set_verbosity(DMRG::VERBOSITY::SILENT);
1014 Hterm.set_verbosity(DMRG::VERBOSITY::SILENT);
1015 Terms = MODEL::sum(Terms,Hterm);
1016 }
1017
1018 vector<Param> params_full = params_loc;
1019 if constexpr (MODEL::FAMILY == HUBBARD)
1020 {
1021 params_full.push_back({"Vxyfull",Ghopx2});
1022 }
1023 else if constexpr (MODEL::FAMILY == HEISENBERG)
1024 {
1025 params_full.push_back({"Jxyfull",Ghopx2});
1026 }
1027
1028 Mpo<typename MODEL::Symmetry,double> Hmpo(Terms);
1029 Hfull = MODEL(Hmpo,params_full);
1030 lout << CompressionTimer.info("MPO compression full") << endl;
1031
1032 if (SAVE) H.save(filename);
1033 }
1034}
1035
1036template<>
1037vector<Mpo<typename VMPS::HubbardU1::Symmetry>> NuclearManager<VMPS::HubbardU1>::
1038make_cdagj (const VMPS::HubbardU1 &H_input) const
1039{
1040 vector<Mpo<typename VMPS::HubbardU1::Symmetry>> res;
1041 res.resize(Nlev);
1042 for (int j=0; j<Nlev; ++j)
1043 {
1044 for (int i=0; i<deg[j]; ++i)
1045 {
1046 auto cdagj1 = H_input.template cdag<UP>(offset(j)+i);
1047 auto cdagj2 = H_input.template cdag<DN>(offset(j)+i);
1048 if (i==0)
1049 {
1050 res[j] = sum(cdagj1,cdagj2);
1051 }
1052 else
1053 {
1054 res[j] = sum(res[j],cdagj1);
1055 res[j] = sum(res[j],cdagj2);
1056 }
1057 }
1058 }
1059 return res;
1060}
1061
1062template<>
1063Mpo<typename VMPS::HubbardU1::Symmetry> NuclearManager<VMPS::HubbardU1>::
1064make_cdagj (int j, const VMPS::HubbardU1 &H_input) const
1065{
1066 Mpo<typename VMPS::HubbardU1::Symmetry> res;
1067 for (int i=0; i<deg[j]; ++i)
1068 {
1069 auto cdagj1 = H_input.template cdag<UP>(offset(j)+i);
1070 auto cdagj2 = H_input.template cdag<DN>(offset(j)+i);
1071 if (i==0)
1072 {
1073 res = sum(cdagj1,cdagj2);
1074 }
1075 else
1076 {
1077 res = sum(res,cdagj1);
1078 res = sum(res,cdagj2);
1079 }
1080 }
1081 return res;
1082}
1083
1084template<>
1085vector<Mpo<typename VMPS::HubbardU1xU1::Symmetry>> NuclearManager<VMPS::HubbardU1xU1>::
1086make_cdagj (const VMPS::HubbardU1xU1 &H_input) const
1087{
1088 vector<Mpo<typename VMPS::HubbardU1xU1::Symmetry>> res;
1089 res.resize(Nlev);
1090 for (int j=0; j<Nlev; ++j)
1091 {
1092 for (int i=0; i<deg[j]; ++i)
1093 {
1094 auto cdagj1 = H_input.template cdag<UP>(offset(j)+i);
1095// auto cdagj2 = H_input.template cdag<DN>(offset(j)+i);
1096 if (i==0)
1097 {
1098 res[j] = cdagj1; //sum(cdagj1,cdagj2);
1099 }
1100 else
1101 {
1102 res[j] = sum(res[j],cdagj1);
1103// res[j] = sum(res[j],cdagj2);
1104 }
1105 }
1106 }
1107 return res;
1108}
1109
1110template<>
1111Mpo<typename VMPS::HubbardU1xU1::Symmetry> NuclearManager<VMPS::HubbardU1xU1>::
1112make_cdagj (int j, const VMPS::HubbardU1xU1 &H_input) const
1113{
1114 Mpo<typename VMPS::HubbardU1xU1::Symmetry> res;
1115 for (int i=0; i<deg[j]; ++i)
1116 {
1117 auto cdagj1 = H_input.template cdag<UP>(offset(j)+i);
1118// auto cdagj2 = H_input.template cdag<DN>(offset(j)+i);
1119 if (i==0)
1120 {
1121 res = cdagj1; //sum(cdagj1,cdagj2);
1122 }
1123 else
1124 {
1125 res = sum(res,cdagj1);
1126// res = sum(res,cdagj2);
1127 }
1128 }
1129 return res;
1130}
1131
1132template<>
1133vector<Mpo<typename VMPS::HeisenbergU1::Symmetry>> NuclearManager<VMPS::HeisenbergU1>::
1134make_cdagj (const VMPS::HeisenbergU1 &H_input) const
1135{
1136 vector<Mpo<typename VMPS::HeisenbergU1::Symmetry>> res;
1137 return res;
1138}
1139
1140/*template<>
1141Mpo<typename VMPS::HeisenbergU1::Symmetry> NuclearManager<VMPS::HeisenbergU1>::
1142make_cdagj (int j, const VMPS::HeisenbergU1 &H_input) const
1143{
1144 Mpo<typename VMPS::HeisenbergU1::Symmetry> res;
1145 return res;
1146}*/
1147
1148template<typename MODEL>
1149vector<tuple<orbital,orbital,orbital,orbital,double,orbinfo>> NuclearManager<MODEL>::
1150generate_couplingList (int J) const
1151{
1152 vector<tuple<orbital,orbital,orbital,orbital,double,orbinfo>> res;
1153
1154 for (int r1=0; r1<2*L; ++r1)
1155 for (int r2=0; r2<2*L; ++r2)
1156 for (int r3=0; r3<2*L; ++r3)
1157 for (int r4=0; r4<2*L; ++r4)
1158 {
1159 orbital o1(r1);
1160 orbital o2(r2);
1161 orbital o3(r3);
1162 orbital o4(r4);
1163
1164 int j1x2 = (2*jlist[r1]).numerator();
1165 int j2x2 = (2*jlist[r2]).numerator();
1166 int j3x2 = (2*jlist[r3]).numerator();
1167 int j4x2 = (2*jlist[r4]).numerator();
1168
1169 int m1x2 = (2*mlist[r1]).numerator();
1170 int m2x2 = (2*mlist[r2]).numerator();
1171 int m3x2 = (2*mlist[r3]).numerator();
1172 int m4x2 = (2*mlist[r4]).numerator();
1173
1174 int lev1 = orblist[r1];
1175 int lev2 = orblist[r2];
1176 int lev3 = orblist[r3];
1177 int lev4 = orblist[r4];
1178
1179 orbinfo oinfo;
1180 oinfo.jx2[0] = j1x2;
1181 oinfo.jx2[1] = j2x2;
1182 oinfo.jx2[2] = j3x2;
1183 oinfo.jx2[3] = j4x2;
1184 oinfo.lev[0] = lev1;
1185 oinfo.lev[1] = lev2;
1186 oinfo.lev[2] = lev3;
1187 oinfo.lev[3] = lev4;
1188 oinfo.label[0] = labels[lev1];
1189 oinfo.label[1] = labels[lev2];
1190 oinfo.label[2] = labels[lev3];
1191 oinfo.label[3] = labels[lev4];
1192
1193 double factor = 0.;
1194 if (lev1==lev2 and lev3==lev4)
1195 //if (lev1==lev2 and lev3==lev4 and m1x2==-m2x2 and m3x2==-m4x2)
1196 {
1197 for (int M=-J; M<=J; ++M)
1198 {
1199 factor += cgc(j1x2,j2x2,2*J, m1x2,m2x2,2*M) * cgc(j3x2,j4x2,2*J, m3x2,m4x2,2*M);
1200 }
1201 //factor += cgc(j1x2,j2x2,2*J, m1x2,m2x2,2*0) * cgc(j3x2,j4x2,2*J, m3x2,m4x2,2*0);
1202 }
1203 if (factor != 0.)
1204 {
1205 bool ADDED = false;
1206 for (int i=0; i<res.size(); ++i)
1207 {
1208 if (
1209 (get<0>(res[i])==o2 and get<1>(res[i])==o1 and // exchanged -> minus
1210 get<2>(res[i])==o4 and get<3>(res[i])==o3 and
1211 get<5>(res[i])==oinfo)
1212 or
1213 (get<0>(res[i])==o1 and get<1>(res[i])==o2 and
1214 get<2>(res[i])==o3 and get<3>(res[i])==o4 and // exchanged -> minus
1215 get<5>(res[i])==oinfo)
1216 )
1217 {
1218 get<4>(res[i]) -= factor;
1219 ADDED = true;
1220 }
1221 else if (get<0>(res[i])==o2 and get<1>(res[i])==o1 and // exchanged -> minus
1222 get<2>(res[i])==o3 and get<3>(res[i])==o4 and // doubly exchanged -> plus
1223 get<5>(res[i])==oinfo)
1224 {
1225 get<4>(res[i]) += factor;
1226 ADDED = true;
1227 }
1228 if (ADDED) break;
1229 }
1230 if (!ADDED and abs(factor)>1e-8)
1231 {
1232 // permute o3 & o4 -> minus
1233 // a†(jm)a†(jm)a(j'n')a(j'm') → -a†(jm)a†(jn)a(j'm')a(j'n')
1234 // in order to have ordering a†(↑)a†(↓) * a(↓)a(↑) = L^+ * L^-
1235 res.push_back(make_tuple(o1,o2,o3,o4, -factor, oinfo));
1236 }
1237 }
1238 }
1239 return res;
1240}
1241
1242template<>
1243tuple<Eigenstate<typename VMPS::HubbardU1::StateXd>,string,int> NuclearManager<VMPS::HubbardU1>::
1244calc_gs (int Nshell, LANCZOS::EDGE::OPTION EDGE, bool CALC_VAR, DMRG::VERBOSITY::OPTION VERB) const
1245{
1246 qarray<VMPS::HubbardU1::Symmetry::Nq> Q = VMPS::HubbardU1::singlet(Nshell);
1247
1248 DMRG::CONTROL::GLOB GlobParam;
1249 GlobParam.min_halfsweeps = 12ul;
1250 GlobParam.max_halfsweeps = 24ul;
1251 GlobParam.Minit = 100ul;
1252 GlobParam.Qinit = 100ul;
1253 GlobParam.CONVTEST = DMRG::CONVTEST::VAR_2SITE; // DMRG::CONVTEST::VAR_HSQ;
1254 GlobParam.CALC_S_ON_EXIT = false;
1255 GlobParam.Mlimit = 500ul;
1256// GlobParam.CONVTEST = DMRG::CONVTEST::VAR_HSQ;
1257 GlobParam.tol_eigval = 1e-10;
1258 GlobParam.tol_state = 1e-8;
1259
1260 DMRG::CONTROL::DYN DynParam;
1261 size_t lim2site = 24ul;
1262 DynParam.iteration = [lim2site] (size_t i) {return (i<lim2site and i%2==0)? DMRG::ITERATION::TWO_SITE : DMRG::ITERATION::ONE_SITE;};
1263 DynParam.max_alpha_rsvd = [lim2site] (size_t i) {return (i<=lim2site+4ul)? 1e4:0.;};
1264 size_t Mincr_per = 2ul;
1265 DynParam.Mincr_per = [Mincr_per] (size_t i) {return Mincr_per;};
1266// size_t min_Nsv_val = (Nshell<=4 or Nshell>=min(2*static_cast<int>(L)-4,0))? 1ul:0ul;
1267// DynParam.min_Nsv = [min_Nsv_val] (size_t i) {return min_Nsv_val;};
1268
1269 Stopwatch<> Timer;
1270
1271 string Jlabel_res = "0";
1272 int Jindex_res = -1;
1273 Eigenstate<typename VMPS::HubbardU1::StateXd> gres;
1274 size_t Nguess = 4;
1275 vector<jEigenstate<typename VMPS::HubbardU1::StateXd>> godd;
1276 vector<jEigenstate<typename VMPS::HubbardU1::StateXd>> ginit;
1277
1278 typename VMPS::HubbardU1::Solver DMRGs;
1279 DMRGs.userSetGlobParam();
1280 DMRGs.userSetDynParam();
1281 DMRGs.GlobParam = GlobParam;
1282 DMRGs.DynParam = DynParam;
1283
1284 if (Nshell%2 == 0)
1285 {
1286 DMRGs.edgeState(H, gres, Q, EDGE);
1287 }
1288 else
1289 {
1290 //#pragma omp parallel for
1291 for (int j=0; j<Nlev; ++j)
1292 {
1293 double Ediff = 1.;
1294 int Niter = 0;
1295 Eigenstate<typename VMPS::HubbardU1::StateXd> gprev;
1296 while (Ediff > 1e-6)
1297 {
1298 auto DMRGt = DMRGs;
1299 DMRGt.GlobParam.min_halfsweeps = 1ul;
1300 DMRGt.set_verbosity(DMRG::VERBOSITY::SILENT);
1301 DMRGt.GlobParam.INITDIR = (Niter%2==1)? DMRG::DIRECTION::LEFT:DMRG::DIRECTION::RIGHT;
1302 DMRGt.edgeState(Hodd[j], gprev, {Nshell-1}, EDGE);
1303 Ediff = abs(gprev.energy-g[Nshell-1].energy);
1304 if (VERB > DMRG::VERBOSITY::SILENT) lout << "Ediff=" << Ediff << endl;
1305 ++Niter;
1306 if (Niter == 4) break;
1307 }
1308 jEigenstate<typename VMPS::HubbardU1::StateXd> gtmp;
1309 bool INCLUDE = true;
1310 /*if (abs(gprev.energy)>1e-7 and Nshell>L)
1311 {
1312 INCLUDE = false;
1313 }*/
1314 if (2.*deg(j)-avgN[Nshell-1](j) > 0.01 and INCLUDE)
1315 {
1316 double compr_tol = (Nshell<=10 or Nshell>=max(2*static_cast<int>(L)-10,0))? 2.:1e-8;
1317 OxV_exact(make_cdagj(j,Hodd[j]), gprev.state, gtmp.eigenstate.state, compr_tol, DMRG::VERBOSITY::SILENT);
1318 gtmp.eigenstate.state /= sqrt(dot(gtmp.eigenstate.state,gtmp.eigenstate.state));
1319 gtmp.eigenstate.energy = avg(gtmp.eigenstate.state, Hodd[j], gtmp.eigenstate.state);
1320 gtmp.label = labels[j];
1321 gtmp.index = j;
1322 //#pragma omp critical
1323 {
1324 ginit.push_back(gtmp);
1325 }
1326 }
1327 if (VERB > DMRG::VERBOSITY::SILENT)
1328 {
1329 #pragma omp critical
1330 {
1331 lout << "guess: j=" << j
1332 << ", label=" << labels[j]
1333 << ", energy=" << gtmp.eigenstate.energy
1334 << ", Niter=" << Niter << endl;
1335 }
1336 }
1337 }
1338
1339 assert(ginit.size()>0);
1340
1341 sort(ginit.begin(), ginit.end(), [] (const auto& lhs, const auto& rhs) {return lhs.eigenstate.energy < rhs.eigenstate.energy;});
1342
1343 // Calculate the ground state for the lowest guesses
1344 //#pragma omp parallel for
1345 for (int i=0; i<min(Nguess,ginit.size()); ++i)
1346 {
1347 int j = ginit[i].index;
1348 godd.push_back(ginit[i]);
1349
1350 double var = 1.;
1351 int Niter = 0;
1352 while (var > 1e-5)
1353 {
1354 if (Niter>0 and VERB > DMRG::VERBOSITY::SILENT)
1355 {
1356 #pragma omp critical
1357 {
1358 lout << termcolor::yellow << "Nshell=" << Nshell
1359 << ": must restart j=" << j << " with var=" << var << termcolor::reset << endl;
1360 }
1361 }
1362
1363 auto DMRGr = DMRGs;
1364 DMRGr.GlobParam.min_halfsweeps = 2ul;
1365 DMRGr.set_verbosity(DMRG::VERBOSITY::SILENT);
1366 DMRGr.edgeState(Hodd[j], godd[i].eigenstate, Q, EDGE, (Niter==0)?true:false);
1367 var = abs(avg(godd[i].eigenstate.state,Hodd[j],Hodd[j],godd[i].eigenstate.state)-pow(godd[i].eigenstate.energy,2))/L;
1368 if (Niter>0 and var<1e-5)
1369 {
1370 #pragma omp critical
1371 {
1372 lout << "Nshell=" << Nshell << ", j=" << j << ", new try brought: var=" << var << endl;
1373 }
1374 }
1375 ++Niter;
1376 if (Niter == 21) break;
1377 }
1378 }
1379
1380 sort(godd.begin(), godd.end(), [] (const auto& lhs, const auto& rhs) {return lhs.eigenstate.energy < rhs.eigenstate.energy;});
1381
1382 gres = godd[0].eigenstate;
1383 Jlabel_res = godd[0].label;
1384 Jindex_res = godd[0].index;
1385
1386 if (VERB > DMRG::VERBOSITY::SILENT)
1387 {
1388 for (int i=0; i<min(Nguess,ginit.size()); ++i)
1389 {
1390 lout << "j=" << godd[i].index << ", label=" << godd[i].label << ", energy=" << godd[i].eigenstate.energy << endl;
1391// godd[i].eigenstate.state.graph(make_string("j=",godd[i].index));
1392 }
1393 }
1394 }
1395
1396 if (VERB > DMRG::VERBOSITY::SILENT)
1397 {
1398 lout << termcolor::blue << setprecision(16) << gres.energy << setprecision(6) << termcolor::reset << endl;
1399 }
1400
1401 if (CALC_VAR)
1402 {
1403 double var;
1404 if (Nshell%2 == 0)
1405 {
1406 var = abs(avg(gres.state,H,H,gres.state)-pow(gres.energy,2))/L;
1407 }
1408 else
1409 {
1410 var = abs(avg(gres.state,Hodd[Jindex_res],Hodd[Jindex_res],gres.state)-pow(gres.energy,2))/L;
1411 }
1412 if (VERB > DMRG::VERBOSITY::SILENT) lout << "var=" << var << endl;
1413 }
1414
1415 if (VERB > DMRG::VERBOSITY::SILENT) lout << Timer.info("edge state") << endl;
1416
1417 return {gres, Jlabel_res, Jindex_res};
1418}
1419
1420template<>
1421tuple<Eigenstate<typename VMPS::HubbardU1xU1::StateXd>,string,int> NuclearManager<VMPS::HubbardU1xU1>::
1422calc_gs (int Nshell, LANCZOS::EDGE::OPTION EDGE, bool CALC_VAR, DMRG::VERBOSITY::OPTION VERB) const
1423{
1424 DMRG::CONTROL::GLOB GlobParam;
1425 GlobParam.min_halfsweeps = 12ul;
1426 GlobParam.max_halfsweeps = 24ul;
1427 GlobParam.Minit = 100ul;
1428 GlobParam.Qinit = 100ul;
1429 GlobParam.CONVTEST = DMRG::CONVTEST::VAR_2SITE; // DMRG::CONVTEST::VAR_HSQ;
1430 GlobParam.CALC_S_ON_EXIT = false;
1431 GlobParam.Mlimit = 500ul;
1432// GlobParam.CONVTEST = DMRG::CONVTEST::VAR_HSQ;
1433 GlobParam.tol_eigval = 1e-10;
1434 GlobParam.tol_state = 1e-8;
1435
1436 DMRG::CONTROL::DYN DynParam;
1437 size_t lim2site = 24ul;
1438 DynParam.iteration = [lim2site] (size_t i) {return (i<lim2site and i%2==0)? DMRG::ITERATION::TWO_SITE : DMRG::ITERATION::ONE_SITE;};
1439 DynParam.max_alpha_rsvd = [lim2site] (size_t i) {return (i<=lim2site+4ul)? 1e4:0.;};
1440 size_t Mincr_per = 2ul;
1441 DynParam.Mincr_per = [Mincr_per] (size_t i) {return Mincr_per;};
1442// size_t min_Nsv_val = (Nshell<=4 or Nshell>=min(2*static_cast<int>(L)-4,0))? 1ul:0ul;
1443// DynParam.min_Nsv = [min_Nsv_val] (size_t i) {return min_Nsv_val;};
1444
1445 Stopwatch<> Timer;
1446
1447 string Jlabel_res = "0";
1448 int Jindex_res = -1;
1449 Eigenstate<typename VMPS::HubbardU1xU1::StateXd> gres;
1450 vector<jEigenstate<typename VMPS::HubbardU1xU1::StateXd>> godd;
1451
1452 typename VMPS::HubbardU1xU1::Solver DMRGs;
1453 DMRGs.userSetGlobParam();
1454 DMRGs.userSetDynParam();
1455 DMRGs.GlobParam = GlobParam;
1456 DMRGs.DynParam = DynParam;
1457
1458 if (Nshell%2 == 0)
1459 {
1460 qarray<VMPS::HubbardU1xU1::Symmetry::Nq> Q = VMPS::HubbardU1xU1::singlet(Nshell);
1461 DMRGs.edgeState(H, gres, Q, EDGE);
1462 }
1463 else
1464 {
1465 //#pragma omp parallel for
1466 for (int j=0; j<Nlev; ++j)
1467 {
1468 qarray<VMPS::HubbardU1xU1::Symmetry::Nq> Q = {2*deg(j)-1,Nshell};
1469
1470 typename VMPS::HubbardU1xU1::Solver DMRGr;
1471 DMRGr.userSetGlobParam();
1472 DMRGr.userSetDynParam();
1473 DMRGr.GlobParam = GlobParam;
1474 DMRGr.DynParam = DynParam;
1475 jEigenstate<typename VMPS::HubbardU1xU1::StateXd> gtry;
1476 gtry.label = labels[j];
1477 gtry.index = j;
1478 DMRGr.edgeState(Hfull, gtry.eigenstate, Q, EDGE);
1479 godd.push_back(gtry);
1480 cout << "j=" << j << ", label=" << labels[j] << ", energy=" << setprecision(16) << gtry.eigenstate.energy << setprecision(6) << ", Q=" << Q << endl;
1481 }
1482 gres = godd[0].eigenstate;
1483 Jlabel_res = godd[0].label;
1484 Jindex_res = godd[0].index;
1485 for (int j=1; j<Nlev; ++j)
1486 {
1487 if (godd[j].eigenstate.energy < gres.energy)
1488 {
1489 gres = godd[j].eigenstate;
1490 Jlabel_res = godd[j].label;
1491 Jindex_res = godd[j].index;
1492 }
1493 }
1494 }
1495
1496 if (VERB > DMRG::VERBOSITY::SILENT)
1497 {
1498 lout << termcolor::blue << setprecision(16) << gres.energy << setprecision(6) << termcolor::reset << endl;
1499 }
1500
1501 if (CALC_VAR)
1502 {
1503 double var;
1504 if (Nshell%2 == 0)
1505 {
1506 var = abs(avg(gres.state,H,H,gres.state)-pow(gres.energy,2))/L;
1507 }
1508 else
1509 {
1510// var = abs(avg(gres.state,Hodd[Jindex_res],Hodd[Jindex_res],gres.state)-pow(gres.energy,2))/L;
1511 var = abs(avg(gres.state,Hfull,Hfull,gres.state)-pow(gres.energy,2))/L;
1512 }
1513 if (VERB > DMRG::VERBOSITY::SILENT) lout << "var=" << var << endl;
1514 }
1515
1516 if (VERB > DMRG::VERBOSITY::SILENT) lout << Timer.info("edge state") << endl;
1517
1518 return {gres, Jlabel_res, Jindex_res};
1519}
1520
1521template<typename MODEL>
1522Eigenstate<typename MODEL::StateXd> NuclearManager<MODEL>::
1523calc_gs_full (int Nshell, LANCZOS::EDGE::OPTION EDGE, bool CALC_VAR, DMRG::VERBOSITY::OPTION VERB) const
1524{
1525 assert(Nshell%2 == 0);
1526 assert(Hfull.length() != 0 and "Calculate full Hamiltonian first!");
1527
1528 DMRG::CONTROL::GLOB GlobParam;
1529 GlobParam.min_halfsweeps = 12ul;
1530 GlobParam.max_halfsweeps = 24ul;
1531 GlobParam.Minit = 100ul;
1532 GlobParam.Qinit = 100ul;
1533 GlobParam.CONVTEST = DMRG::CONVTEST::VAR_2SITE; // DMRG::CONVTEST::VAR_HSQ;
1534 GlobParam.CALC_S_ON_EXIT = false;
1535 GlobParam.Mlimit = 500ul;
1536// GlobParam.CONVTEST = DMRG::CONVTEST::VAR_HSQ;
1537 GlobParam.tol_eigval = 1e-10;
1538 GlobParam.tol_state = 1e-8;
1539
1540 DMRG::CONTROL::DYN DynParam;
1541 size_t lim2site = 24ul;
1542 DynParam.iteration = [lim2site] (size_t i) {return (i<lim2site and i%2==0)? DMRG::ITERATION::TWO_SITE : DMRG::ITERATION::ONE_SITE;};
1543 DynParam.max_alpha_rsvd = [lim2site] (size_t i) {return (i<=lim2site+4ul)? 1e4:0.;};
1544 size_t Mincr_per = 2ul;
1545 DynParam.Mincr_per = [Mincr_per] (size_t i) {return Mincr_per;};
1546// size_t min_Nsv_val = (Nshell<=4 or Nshell>=min(2*static_cast<int>(L)-4,0))? 1ul:0ul;
1547// DynParam.min_Nsv = [min_Nsv_val] (size_t i) {return min_Nsv_val;};
1548
1549 Eigenstate<typename MODEL::StateXd> gres;
1550
1551 typename MODEL::Solver DMRGs;
1552 DMRGs.userSetGlobParam();
1553 DMRGs.userSetDynParam();
1554 DMRGs.GlobParam = GlobParam;
1555 DMRGs.DynParam = DynParam;
1556
1557 if (Nshell%2 == 0)
1558 {
1559 qarray<MODEL::Symmetry::Nq> Q = MODEL::singlet(Nshell);
1560 DMRGs.edgeState(Hfull, gres, Q, EDGE);
1561 }
1562
1563 if (VERB > DMRG::VERBOSITY::SILENT)
1564 {
1565 lout << termcolor::blue << setprecision(16) << gres.energy << setprecision(6) << termcolor::reset << endl;
1566 }
1567
1568 return gres;
1569}
1570
1571/*template<>
1572Eigenstate<typename VMPS::HubbardU1::StateXd> NuclearManager<VMPS::HubbardU1>::
1573calc_gs_full (int Nshell, LANCZOS::EDGE::OPTION EDGE, int Nruns, bool CALC_VAR, DMRG::VERBOSITY::OPTION VERB) const
1574{
1575 qarray<VMPS::HubbardU1::Symmetry::Nq> Q = VMPS::HubbardU1::singlet(Nshell);
1576
1577 DMRG::CONTROL::GLOB GlobParam;
1578 GlobParam.min_halfsweeps = 12ul;
1579 GlobParam.max_halfsweeps = 36ul;
1580 GlobParam.Minit = 100ul;
1581 GlobParam.Qinit = 100ul;
1582 GlobParam.CONVTEST = DMRG::CONVTEST::VAR_2SITE; // DMRG::CONVTEST::VAR_HSQ;
1583 GlobParam.CALC_S_ON_EXIT = false;
1584 GlobParam.Mlimit = 500ul;
1585// GlobParam.CONVTEST = DMRG::CONVTEST::VAR_HSQ;
1586 GlobParam.tol_eigval = 1e-9;
1587 GlobParam.tol_state = 1e-7;
1588
1589 DMRG::CONTROL::DYN DynParam;
1590 size_t lim2site = 24ul;
1591 DynParam.iteration = [lim2site] (size_t i) {return (i<lim2site and i%2==0)? DMRG::ITERATION::TWO_SITE : DMRG::ITERATION::ONE_SITE;};
1592 DynParam.max_alpha_rsvd = [lim2site] (size_t i) {return (i<=lim2site+4ul)? 1e4:0.;};
1593 size_t Mincr_per = 2ul;
1594 DynParam.Mincr_per = [Mincr_per] (size_t i) {return Mincr_per;};
1595
1596 Stopwatch<> Timer;
1597
1598 int Nruns_adjusted = (Nshell%2==0)? Nruns:1;
1599 if (Nshell==0 or Nshell==2*L) Nruns_adjusted = 1;
1600// if (Nshell==42 and Nlev==13) Nruns_adjusted = 10;
1601
1602 vector<Eigenstate<typename VMPS::HubbardU1::StateXd>> gL(Nruns_adjusted);
1603 vector<Eigenstate<typename VMPS::HubbardU1::StateXd>> gR(Nruns_adjusted);
1604 ArrayXd energiesL(Nruns_adjusted);
1605 ArrayXd energiesR(Nruns_adjusted);
1606 string JgsL, JgsR;
1607
1608 vector<jEigenstate<typename VMPS::HubbardU1::StateXd>> ginit;
1609 if (Nshell%2==1)
1610 {
1611 jEigenstate<typename VMPS::HubbardU1::StateXd> gtmp;
1612 for (int j=0; j<Nlev; ++j)
1613 {
1614// #pragma omp critical
1615// {
1616// lout << "j=" << j << ", label=" << labels[j] << ", free space=" << 2.*deg(j)-avgN[Nshell-1](j) << endl;
1617// }
1618 if (2.*deg(j)-avgN[Nshell-1](j) > 0.01)
1619 {
1620 double compr_tol = (Nshell<=10 or Nshell>=max(2*static_cast<int>(L)-10,0))? 2.:1e-8;
1621// double compr_tol = 2.;
1622 OxV_exact(cdagj[j], g[Nshell-1].state, gtmp.eigenstate.state, compr_tol, DMRG::VERBOSITY::SILENT);
1623 gtmp.eigenstate.state /= sqrt(dot(gtmp.eigenstate.state,gtmp.eigenstate.state));
1624 gtmp.eigenstate.energy = avg(gtmp.eigenstate.state, H, gtmp.eigenstate.state);
1625 ginit.push_back(gtmp);
1626 ginit[ginit.size()-1].label = labels[j];
1627 }
1628 }
1629
1630 sort(ginit.begin(), ginit.end(), [] (const auto& lhs, const auto& rhs) {return lhs.eigenstate.energy < rhs.eigenstate.energy;});
1631
1632 if (VERB > DMRG::VERBOSITY::SILENT)
1633 {
1634 for (int j=0; j<ginit.size(); ++j)
1635 {
1636 lout << "E(" << ginit[j].label << ")=" << ginit[j].eigenstate.energy << endl;
1637 }
1638 }
1639 }
1640
1641 for (int r=0; r<Nruns_adjusted; ++r)
1642 {
1643 if (VERB > DMRG::VERBOSITY::SILENT) lout << "------r=" << r << "------" << endl;
1644
1645 std::array<typename VMPS::HubbardU1::Solver,2> DMRG_LR;
1646 for (int i=0; i<2; ++i)
1647 {
1648 DMRG_LR[i] = typename VMPS::HubbardU1::Solver(DMRG::VERBOSITY::SILENT);
1649 DMRG_LR[i].userSetGlobParam();
1650 DMRG_LR[i].userSetDynParam();
1651 DMRG_LR[i].GlobParam = GlobParam;
1652 DMRG_LR[i].DynParam = DynParam;
1653 DMRG_LR[i].Epenalty = (EDGE==LANCZOS::EDGE::GROUND)? 1e4:-1e4;
1654 }
1655
1656 #pragma omp parallel sections
1657 {
1658 #pragma omp section
1659 {
1660 for (int s=0; s<r; ++s)
1661 {
1662 DMRG_LR[0].push_back(gL[s].state);
1663 DMRG_LR[0].push_back(gR[s].state);
1664 }
1665
1666 if (Nshell%2 == 0)
1667 {
1668 for (int i=0; i<Hperturb.size(); ++i)
1669 {
1670 DMRG_LR[0].GlobParam.min_halfsweeps = perturb_sweeps[i];
1671 DMRG_LR[0].GlobParam.max_halfsweeps = perturb_sweeps[i];
1672 DMRG_LR[0].edgeState(Hperturb[i], gL[r], Q, EDGE, (i==0)?false:true);
1673 }
1674 DMRG_LR[0].GlobParam.min_halfsweeps = 10ul;
1675 DMRG_LR[0].GlobParam.max_halfsweeps = 32ul;
1676 DMRG_LR[0].set_verbosity(VERB);
1677 DMRG_LR[0].edgeState(H, gL[r], Q, EDGE, true);
1678 DMRG_LR[0].edgeState(H, gL[r], Q, EDGE);
1679 }
1680 else
1681 {
1682 DMRG_LR[0].GlobParam.min_halfsweeps = 12ul;
1683 DMRG_LR[0].GlobParam.max_halfsweeps = 36ul;
1684 vector<jEigenstate<typename VMPS::HubbardU1::StateXd>> gfinalL = ginit;
1685 for (int j=0; j<min(gfinalL.size(),4ul); ++j)
1686 {
1687 DMRG_LR[0].edgeState(H, gfinalL[j].eigenstate, Q, EDGE, true);
1688 }
1689 sort(gfinalL.begin(), gfinalL.end(), [] (const auto& lhs, const auto& rhs) {return lhs.eigenstate.energy < rhs.eigenstate.energy;});
1690 gL[r] = gfinalL[0].eigenstate;
1691 JgsL = gfinalL[0].label;
1692 }
1693
1694 energiesL® = gL[r].energy;
1695 }
1696 #pragma omp section
1697 {
1698 DMRG_LR[1].GlobParam.INITDIR = DMRG::DIRECTION::LEFT;
1699 for (int s=0; s<r; ++s)
1700 {
1701 DMRG_LR[1].push_back(gL[s].state);
1702 DMRG_LR[1].push_back(gR[s].state);
1703 }
1704
1705 if (Nshell%2 == 0)
1706 {
1707 for (int i=0; i<Hperturb.size(); ++i)
1708 {
1709 DMRG_LR[1].GlobParam.min_halfsweeps = perturb_sweeps[i];
1710 DMRG_LR[1].GlobParam.max_halfsweeps = perturb_sweeps[i];
1711 DMRG_LR[1].edgeState(Hperturb[i], gR[r], Q, EDGE, (i==0)?false:true);
1712 }
1713 DMRG_LR[1].GlobParam.min_halfsweeps = 10ul;
1714 DMRG_LR[1].GlobParam.max_halfsweeps = 32ul;
1715 DMRG_LR[1].set_verbosity(VERB);
1716 DMRG_LR[1].edgeState(H, gR[r], Q, EDGE, true);
1717 }
1718 else
1719 {
1720 DMRG_LR[1].GlobParam.min_halfsweeps = 12ul;
1721 DMRG_LR[1].GlobParam.max_halfsweeps = 36ul;
1722 vector<jEigenstate<typename VMPS::HubbardU1::StateXd>> gfinalR = ginit;
1723 for (int j=0; j<min(gfinalR.size(),4ul); ++j)
1724 {
1725 DMRG_LR[1].edgeState(H, gfinalR[j].eigenstate, Q, EDGE, true);
1726 }
1727 sort(gfinalR.begin(), gfinalR.end(), [] (const auto& lhs, const auto& rhs) {return lhs.eigenstate.energy < rhs.eigenstate.energy;});
1728 gR[r] = gfinalR[0].eigenstate;
1729 JgsR = gfinalR[0].label;
1730 }
1731
1732 energiesR® = gR[r].energy;
1733 }
1734 }
1735 }
1736
1737 if (VERB > DMRG::VERBOSITY::SILENT)
1738 {
1739 for (int r=0; r<Nruns_adjusted; ++r)
1740 {
1741 lout << termcolor::blue << setprecision(9) << gL[r].energy << "\t" << gR[r].energy << setprecision(6) << termcolor::reset << endl;
1742 }
1743 }
1744
1745 double valL, valR;
1746 int indexL, indexR;
1747 string edgeStateLabel;
1748 if (EDGE == LANCZOS::EDGE::GROUND)
1749 {
1750 valL = energiesL.minCoeff(&indexL);
1751 valR = energiesR.minCoeff(&indexR);
1752 edgeStateLabel = "ground";
1753 }
1754 else
1755 {
1756 valL = energiesL.maxCoeff(&indexL);
1757 valR = energiesR.maxCoeff(&indexR);
1758 edgeStateLabel = "roof";
1759 }
1760 if (VERB > DMRG::VERBOSITY::SILENT)
1761 {
1762 if (indexL == 1 and valL<valR)
1763 {
1764 lout << termcolor::yellow << "Found " << edgeStateLabel << " state at r=" << indexL << " (L→R)" << termcolor::reset << endl;
1765 }
1766 else if (indexL > 1 and valL<valR)
1767 {
1768 lout << termcolor::red << "Found " << edgeStateLabel << " state at r=" << indexL << " (L→R)" << termcolor::reset << endl;
1769 }
1770 if (indexR == 1 and valL>valR)
1771 {
1772 lout << termcolor::yellow << "Found " << edgeStateLabel << " state at r=" << indexR << " (R→L)" << termcolor::reset << endl;
1773 }
1774 else if (indexR > 1 and valL>valR)
1775 {
1776 lout << termcolor::red << "Found " << edgeStateLabel << " state at r=" << indexR << " (R→L)" << termcolor::reset << endl;
1777 }
1778 }
1779
1780 Eigenstate<typename VMPS::HubbardU1::StateXd> res;
1781 string Jres;
1782
1783 if (EDGE == LANCZOS::EDGE::GROUND)
1784 {
1785 res = (gL[0].energy<=gR[0].energy)? gL[0]:gR[0];
1786 if (Nshell%2 == 1)
1787 {
1788 Jres = (gL[0].energy<=gR[0].energy)? JgsL:JgsR;
1789 }
1790 else
1791 {
1792 Jres = "0";
1793 }
1794 for (int r=1; r<Nruns_adjusted; ++r)
1795 {
1796 if (gL[r].energy < res.energy) res = gL[r];
1797 if (gR[r].energy < res.energy) res = gR[r];
1798 }
1799 }
1800 else
1801 {
1802 res = (gL[0].energy>=gR[0].energy)? gL[0]:gR[0];
1803 if (Nshell%2 == 1)
1804 {
1805 Jres = (gL[0].energy>=gR[0].energy)? JgsL:JgsR;
1806 }
1807 else
1808 {
1809 Jres = "0";
1810 }
1811 for (int r=1; r<Nruns_adjusted; ++r)
1812 {
1813 if (gL[r].energy > res.energy) res = gL[r];
1814 if (gR[r].energy > res.energy) res = gR[r];
1815 }
1816 }
1817
1818 if (CALC_VAR)
1819 {
1820 double var = abs(avg(res.state,H,H,res.state)-pow(res.energy,2))/L;
1821 if (VERB > DMRG::VERBOSITY::SILENT) lout << "var=" << var << endl;
1822 }
1823
1824 if (VERB > DMRG::VERBOSITY::SILENT) lout << Timer.info("edge state") << endl;
1825
1826 return res;
1827}*/
1828
1829template<>
1830tuple<Eigenstate<typename VMPS::HeisenbergU1::StateXd>,string,int> NuclearManager<VMPS::HeisenbergU1>::
1831calc_gs (int Nshell, LANCZOS::EDGE::OPTION EDGE, bool CALC_VAR, DMRG::VERBOSITY::OPTION VERB) const
1832{
1833 qarray<VMPS::HeisenbergU1::Symmetry::Nq> Q = {Nshell-static_cast<int>(L)};
1834
1835 DMRG::CONTROL::GLOB GlobParam;
1836 GlobParam.min_halfsweeps = 16ul;
1837 GlobParam.max_halfsweeps = 36ul;
1838 GlobParam.Minit = 100ul;
1839 GlobParam.Qinit = 100ul;
1840 GlobParam.CONVTEST = DMRG::CONVTEST::VAR_2SITE; // DMRG::CONVTEST::VAR_HSQ;
1841 GlobParam.CALC_S_ON_EXIT = false;
1842 GlobParam.Mlimit = 500ul;
1843// GlobParam.CONVTEST = DMRG::CONVTEST::VAR_HSQ;
1844 GlobParam.tol_eigval = 1e-9;
1845 GlobParam.tol_state = 1e-7;
1846
1847 DMRG::CONTROL::DYN DynParam;
1848 size_t lim2site = 24ul;
1849 DynParam.iteration = [lim2site] (size_t i) {return (i<lim2site and i%2==0)? DMRG::ITERATION::TWO_SITE : DMRG::ITERATION::ONE_SITE;};
1850 DynParam.max_alpha_rsvd = [lim2site] (size_t i) {return (i<=lim2site+4ul)? 1e4:0.;};
1851 size_t Mincr_per = 2ul;
1852 DynParam.Mincr_per = [Mincr_per] (size_t i) {return Mincr_per;};
1853// size_t min_Nsv_val = (Nshell<=4 or Nshell>=min(2*static_cast<int>(L)-4,0))? 1ul:0ul;
1854// DynParam.min_Nsv = [min_Nsv_val] (size_t i) {return min_Nsv_val;};
1855
1856 Stopwatch<> Timer;
1857
1858 Eigenstate<typename VMPS::HeisenbergU1::StateXd> gres;
1859
1860 typename VMPS::HeisenbergU1::Solver DMRGs;
1861 DMRGs.userSetGlobParam();
1862 DMRGs.userSetDynParam();
1863 DMRGs.GlobParam = GlobParam;
1864 DMRGs.DynParam = DynParam;
1865
1866 DMRGs.edgeState(H, gres, Q, EDGE);
1867
1868 if (VERB > DMRG::VERBOSITY::SILENT)
1869 {
1870 lout << termcolor::blue << setprecision(16) << gres.energy << setprecision(6) << termcolor::reset << endl;
1871 }
1872
1873 if (CALC_VAR)
1874 {
1875 double var = abs(avg(gres.state,H,H,gres.state)-pow(gres.energy,2))/L;
1876 if (VERB > DMRG::VERBOSITY::SILENT) lout << "var=" << var << endl;
1877 }
1878
1879 if (VERB > DMRG::VERBOSITY::SILENT) lout << Timer.info("edge state") << endl;
1880
1881 return {gres,"0",-1};
1882}
1883
1884template<>
1885void NuclearManager<VMPS::HubbardU1xU1>::
1886compute (bool LOAD_MPO, bool SAVE_MPO, string wd, int minNshell, int maxNshell, bool SAVE)
1887{
1888 make_Hamiltonian(LOAD_MPO,SAVE_MPO,wd,false);
1889 make_fullHamiltonian(LOAD_MPO,SAVE_MPO,wd);
1890
1891 g.resize(2*L+1);
1892 var.resize(2*L+1);
1893 Mmax.resize(2*L+1);
1894 avgN.resize(2*L+1);
1895 n.resize(2*L+1);
1896 Jlabel.resize(2*L+1);
1897 Jindex.resize(2*L+1);
1898 energies.resize(2*L+1);
1899 energies_free.resize(2*L+1);
1900
1901 int Nshelllast = (maxNshell==-1)? 2*L : maxNshell;
1902
1903 for (int Nshell=minNshell; Nshell<=Nshelllast; Nshell+=1)
1904 {
1905 int A = Z+Nclosed+Nshell;
1906 lout << termcolor::bold << "A=" << A << ", Z=" << Z << ", N=" << Nclosed+Nshell << ", Nshell=" << Nshell << "/" << 2*L << ", progress=" << round(Nshell*1./(2*L)*100,1) << "%" << termcolor::reset << endl;
1907
1908 auto [gres, Jlabel_res, Jindex_res] = calc_gs(Nshell, LANCZOS::EDGE::GROUND, false);
1909 g[Nshell] = gres;
1910 Jlabel[Nshell] = Jlabel_res;
1911 Jindex[Nshell] = Jindex_res;
1912
1913 lout << g[Nshell].state.info() << endl;
1914
1915 lout << "J of ground state: " << Jlabel[Nshell] << ", index=" << Jindex[Nshell] << endl;
1916
1917 if (Nshell != 0 and Nshell != 2*L)
1918 {
1919 if (Nshell%2 == 0)
1920 {
1921 var[Nshell] = abs(avg(g[Nshell].state,H,H,g[Nshell].state)-pow(g[Nshell].energy,2))/L;
1922 }
1923 else
1924 {
1925// var[Nshell] = abs(avg(g[Nshell].state,Hodd[Jindex[Nshell]],Hodd[Jindex[Nshell]],g[Nshell].state)-pow(g[Nshell].energy,2))/L;
1926 var[Nshell] = abs(avg(g[Nshell].state,Hfull,Hfull,g[Nshell].state)-pow(g[Nshell].energy,2))/L;
1927 }
1928 }
1929 else
1930 {
1931 var[Nshell] = 0.;
1932 }
1933 lout << termcolor::blue << "var=" << var[Nshell] << termcolor::reset << endl;
1934
1935 Mmax[Nshell] = g[Nshell].state.calc_Mmax();
1936
1937 energies(Nshell) = g[Nshell].energy;
1938 energies_free(Nshell) = onsite_free.head(Nshell).sum();
1939 lout << "noninteracting energy=" << energies_free(Nshell) << ", interacting energy=" << energies(Nshell) << endl;
1940 if (Nshell > 1 and Nshell<2*L-1) assert(g[Nshell].energy < energies_free(Nshell));
1941
1942 if (Z==50 and Nclosed==50 and REF)
1943 {
1944 double diff = abs(g[Nshell].energy-Sn_Eref(Nshell));
1945 if (diff<1e-2) lout << termcolor::green;
1946 else lout << termcolor::red;
1947 lout << "ref=" << Sn_Eref(Nshell) << ", diff=" << diff << endl;
1948 lout << termcolor::reset;
1949 }
1950 lout << "E_B/A=" << abs(g[Nshell].energy)/Nshell << " MeV" << endl;
1951
1952 n[Nshell].resize(Nlev);
1953 avgN[Nshell].resize(Nlev);
1954
1955 for (int j=0; j<Nlev; ++j)
1956 {
1957 avgN[Nshell](j) = 0.;
1958 for (int i=0; i<deg(j); ++i)
1959 {
1960 if (Nshell%2==0)
1961 {
1962 avgN[Nshell](j) += avg(g[Nshell].state, H.n(offset(j)+i), g[Nshell].state);
1963 }
1964 else
1965 {
1966// avgN[Nshell](j) += avg(g[Nshell].state, Hodd[Jindex[Nshell]].n(offset(j)+i), g[Nshell].state);
1967 avgN[Nshell](j) += avg(g[Nshell].state, Hfull.n(offset(j)+i), g[Nshell].state);
1968 }
1969 }
1970 n[Nshell](j) = avgN[Nshell](j)/(2.*deg(j));
1971
1972 double Nplot = (avgN[Nshell](j)<1e-14)? 0 : avgN[Nshell](j);
1973 double nplot = (n[Nshell](j) <1e-14)? 0 : n[Nshell](j);
1974 lout << "N(" << labels[j] << ")=" << Nplot << "/" << 2*deg(j) << ", " << "n=" << nplot;
1975 //<< ", " << "P=" << resP;
1976
1977 if (Z==50 and Nclosed==50 and REF and (Nshell==12 or Nshell==14 or Nshell==16 or Nshell==18))
1978 {
1979 auto Sn_nref_ = Sn_nref(Nshell);
1980 auto Sn_Nref_ = Sn_Nref(Nshell);
1981 auto Sn_Pref_ = Sn_Pref(Nshell);
1982
1983 lout << ", ref(N)=";
1984 if (abs(Sn_Nref_(j)-avgN[Nshell](j))<1e-3)
1985 {
1986 lout << termcolor::green;
1987 }
1988 else
1989 {
1990 lout << termcolor::red;
1991 }
1992 lout << Sn_Nref_(j) << termcolor::reset;
1993
1994 lout << ", ref(n)=";
1995 if (abs(Sn_nref_(j)-n[Nshell](j))<1e-3)
1996 {
1997 lout << termcolor::green;
1998 }
1999 else
2000 {
2001 lout << termcolor::red;
2002 }
2003 lout << Sn_nref_(j) << termcolor::reset;
2004 }
2005 lout << endl;
2006 }
2007 lout << endl;
2008 }
2009
2010 if (SAVE) save(wd, minNshell, Nshelllast);
2011}
2012
2013template<>
2014void NuclearManager<VMPS::HubbardU1>::
2015compute (bool LOAD_MPO, bool SAVE_MPO, string wd, int minNshell, int maxNshell, bool SAVE)
2016{
2017 make_Hamiltonian(LOAD_MPO,SAVE_MPO,wd);
2018
2019 g.resize(2*L+1);
2020 var.resize(2*L+1);
2021 Mmax.resize(2*L+1);
2022 avgN.resize(2*L+1);
2023 n.resize(2*L+1);
2024 Jlabel.resize(2*L+1);
2025 Jindex.resize(2*L+1);
2026 energies.resize(2*L+1);
2027 energies_free.resize(2*L+1);
2028
2029 int Nshelllast = (maxNshell==-1)? 2*L : maxNshell;
2030
2031 for (int Nshell=minNshell; Nshell<=Nshelllast; ++Nshell)
2032 {
2033 int A = Z+Nclosed+Nshell;
2034 lout << termcolor::bold << "A=" << A << ", Z=" << Z << ", N=" << Nclosed+Nshell << ", Nshell=" << Nshell << "/" << 2*L << ", progress=" << round(Nshell*1./(2*L)*100,1) << "%" << termcolor::reset << endl;
2035
2036 auto [gres, Jlabel_res, Jindex_res] = calc_gs(Nshell, LANCZOS::EDGE::GROUND, false);
2037 g[Nshell] = gres;
2038 Jlabel[Nshell] = Jlabel_res;
2039 Jindex[Nshell] = Jindex_res;
2040
2041 lout << g[Nshell].state.info() << endl;
2042
2043 lout << "J of ground state: " << Jlabel[Nshell] << ", index=" << Jindex[Nshell] << endl;
2044
2045 if (Nshell != 0 and Nshell != 2*L)
2046 {
2047 if (Nshell%2 == 0)
2048 {
2049 var[Nshell] = abs(avg(g[Nshell].state,H,H,g[Nshell].state)-pow(g[Nshell].energy,2))/L;
2050 }
2051 else
2052 {
2053 var[Nshell] = abs(avg(g[Nshell].state,Hodd[Jindex[Nshell]],Hodd[Jindex[Nshell]],g[Nshell].state)-pow(g[Nshell].energy,2))/L;
2054 }
2055 }
2056 else
2057 {
2058 var[Nshell] = 0.;
2059 }
2060 lout << termcolor::blue << "var=" << var[Nshell] << termcolor::reset << endl;
2061
2062 Mmax[Nshell] = g[Nshell].state.calc_Mmax();
2063
2064 energies(Nshell) = g[Nshell].energy;
2065 energies_free(Nshell) = onsite_free.head(Nshell).sum();
2066 lout << "noninteracting energy=" << energies_free(Nshell) << ", interacting energy=" << energies(Nshell) << endl;
2067 if (Nshell > 1 and Nshell<2*L-1) assert(g[Nshell].energy < energies_free(Nshell));
2068
2069 if (Z==50 and Nclosed==50 and REF)
2070 {
2071 double diff = abs(g[Nshell].energy-Sn_Eref(Nshell));
2072 if (diff<1e-2) lout << termcolor::green;
2073 else lout << termcolor::red;
2074 lout << "ref=" << Sn_Eref(Nshell) << ", diff=" << diff << endl;
2075 lout << termcolor::reset;
2076 }
2077 lout << "E_B/A=" << abs(g[Nshell].energy)/Nshell << " MeV" << endl;
2078
2079 n[Nshell].resize(Nlev);
2080 avgN[Nshell].resize(Nlev);
2081
2082 for (int j=0; j<Nlev; ++j)
2083 {
2084 avgN[Nshell](j) = 0.;
2085 for (int i=0; i<deg(j); ++i)
2086 {
2087 if (Nshell%2==0)
2088 {
2089 avgN[Nshell](j) += avg(g[Nshell].state, H.n(offset(j)+i), g[Nshell].state);
2090 }
2091 else
2092 {
2093 avgN[Nshell](j) += avg(g[Nshell].state, Hodd[Jindex[Nshell]].n(offset(j)+i), g[Nshell].state);
2094 }
2095 }
2096 n[Nshell](j) = avgN[Nshell](j)/(2.*deg(j));
2097
2098 double Nplot = (avgN[Nshell](j)<1e-14)? 0 : avgN[Nshell](j);
2099 double nplot = (n[Nshell](j) <1e-14)? 0 : n[Nshell](j);
2100 lout << "N(" << labels[j] << ")=" << Nplot << "/" << 2*deg(j) << ", " << "n=" << nplot;
2101 //<< ", " << "P=" << resP;
2102
2103 if (Z==50 and Nclosed==50 and REF and (Nshell==12 or Nshell==14 or Nshell==16 or Nshell==18))
2104 {
2105 auto Sn_nref_ = Sn_nref(Nshell);
2106 auto Sn_Nref_ = Sn_Nref(Nshell);
2107 auto Sn_Pref_ = Sn_Pref(Nshell);
2108
2109 lout << ", ref(N)=";
2110 if (abs(Sn_Nref_(j)-avgN[Nshell](j))<1e-3)
2111 {
2112 lout << termcolor::green;
2113 }
2114 else
2115 {
2116 lout << termcolor::red;
2117 }
2118 lout << Sn_Nref_(j) << termcolor::reset;
2119
2120 lout << ", ref(n)=";
2121 if (abs(Sn_nref_(j)-n[Nshell](j))<1e-3)
2122 {
2123 lout << termcolor::green;
2124 }
2125 else
2126 {
2127 lout << termcolor::red;
2128 }
2129 lout << Sn_nref_(j) << termcolor::reset;
2130 }
2131 lout << endl;
2132 }
2133 lout << endl;
2134 }
2135
2136 // calc P
2137// double resP = 0.;
2138// if (Nshell>=minNshell+2)
2139// {
2140// for (int i=0; i<deg(j); ++i)
2141// {
2142// double Pcontrib = sign(offset(j)+i) * avg(g[Nshell-2].state, H.cc(offset(j)+i), g[Nshell].state);
2144// resP += Pcontrib;
2145// }
2146// }
2147// resP /= sqrt(deg(j));
2148// resP = abs(resP);
2149
2150 // Pref
2151// lout << ", refâ„—=";
2152// if (abs(Sn_Pref_(j)-resP)<1e-3)
2153// {
2154// lout << termcolor::green;
2155// }
2156// else
2157// {
2158// lout << termcolor::red;
2159// }
2160// lout << Sn_Pref_(j) << termcolor::reset;
2161
2162 if (SAVE) save(wd, minNshell, Nshelllast);
2163}
2164
2165template<>
2166void NuclearManager<VMPS::HeisenbergU1>::
2167compute (bool LOAD_MPO, bool SAVE_MPO, string wd, int minNshell, int maxNshell, bool SAVE)
2168{
2169 make_Hamiltonian(LOAD_MPO,SAVE_MPO,wd);
2170
2171 g.resize(2*L+1);
2172 var.resize(2*L+1);
2173 Mmax.resize(2*L+1);
2174 avgN.resize(2*L+1);
2175 n.resize(2*L+1);
2176 Jlabel.resize(2*L+1);
2177 Jindex.resize(2*L+1);
2178 energies.resize(2*L+1);
2179 energies_free.resize(2*L+1);
2180
2181 int Nshelllast = (maxNshell==-1)? 2*L : maxNshell;
2182
2183 for (int Nshell=minNshell; Nshell<=Nshelllast; Nshell+=2)
2184 {
2185 int A = Z+Nclosed+Nshell;
2186 lout << termcolor::bold << "A=" << A << ", Z=" << Z << ", N=" << Nclosed+Nshell << ", Nshell=" << Nshell << "/" << 2*L << ", progress=" << round(Nshell*1./(2*L)*100,1) << "%" << termcolor::reset << endl;
2187
2188 auto [gres, Jlabel_res, Jindex_res] = calc_gs(Nshell, LANCZOS::EDGE::GROUND, false);
2189 g[Nshell] = gres;
2190 Jlabel[Nshell] = Jlabel_res;
2191 Jindex[Nshell] = Jindex_res;
2192
2193 lout << g[Nshell].state.info() << endl;
2194
2195 if (Nshell != 0 and Nshell != 2*L)
2196 {
2197 var[Nshell] = abs(avg(g[Nshell].state,H,H,g[Nshell].state)-pow(g[Nshell].energy,2))/L;
2198 }
2199 else
2200 {
2201 var[Nshell] = 0.;
2202 }
2203 lout << termcolor::blue << "var=" << var[Nshell] << termcolor::reset << endl;
2204
2205 Mmax[Nshell] = g[Nshell].state.calc_Mmax();
2206
2207 energies(Nshell) = g[Nshell].energy;
2208 energies_free(Nshell) = onsite_free.head(Nshell).sum();
2209 lout << "noninteracting energy=" << energies_free(Nshell) << ", interacting energy=" << energies(Nshell) << endl;
2210 if (Nshell > 1 and Nshell<2*L-1) assert(g[Nshell].energy < energies_free(Nshell));
2211
2212 if (Z==50 and Nclosed==50 and REF and (Nshell==12 or Nshell==14 or Nshell==16))
2213 {
2214 double diff = abs(g[Nshell].energy-Sn_Eref(Nshell));
2215 if (diff<1e-2)
2216 {
2217 lout << termcolor::green;
2218 }
2219 else
2220 {
2221 lout << termcolor::red;
2222 }
2223 lout << "ref=" << Sn_Eref(Nshell) << endl;
2224 lout << termcolor::reset;
2225 }
2226 lout << "E_B/A=" << abs(g[Nshell].energy)/Nshell << " MeV" << endl;
2227
2228 n[Nshell].resize(Nlev);
2229 avgN[Nshell].resize(Nlev);
2230
2231 for (int j=0; j<Nlev; ++j)
2232 {
2233 avgN[Nshell](j) = 0.;
2234 for (int i=0; i<deg(j); ++i)
2235 {
2236 avgN[Nshell](j) += 2.*avg(g[Nshell].state, H.Sz(offset(j)+i), g[Nshell].state)+1.;
2237 }
2238 n[Nshell](j) = avgN[Nshell](j)/(2.*deg(j));
2239
2240 double Nplot = (avgN[Nshell](j)<1e-14)? 0 : avgN[Nshell](j);
2241 double nplot = (n[Nshell](j)<1e-14)? 0 : n[Nshell](j);
2242 lout << "N(" << labels[j] << ")=" << Nplot << "/" << 2*deg(j) << ", " << "n=" << nplot;
2243 //<< ", " << "P=" << resP;
2244
2245 if (Z==50 and Nclosed==50 and REF and (Nshell==12 or Nshell==14 or Nshell==16 or Nshell==18))
2246 {
2247 auto Sn_nref_ = Sn_nref(Nshell);
2248 auto Sn_Nref_ = Sn_Nref(Nshell);
2249 auto Sn_Pref_ = Sn_Pref(Nshell);
2250
2251 lout << ", ref(N)=";
2252 if (abs(Sn_Nref_(j)-avgN[Nshell](j))<1e-3)
2253 {
2254 lout << termcolor::green;
2255 }
2256 else
2257 {
2258 lout << termcolor::red;
2259 }
2260 lout << Sn_Nref_(j) << termcolor::reset;
2261
2262 lout << ", ref(n)=";
2263 if (abs(Sn_nref_(j)-n[Nshell](j))<1e-3)
2264 {
2265 lout << termcolor::green;
2266 }
2267 else
2268 {
2269 lout << termcolor::red;
2270 }
2271 lout << Sn_nref_(j) << termcolor::reset;
2272 }
2273 lout << endl;
2274 }
2275 lout << endl;
2276 }
2277
2278 if (SAVE) save(wd, minNshell, Nshelllast, 2);
2279}
2280
2281template<>
2282void NuclearManager<VMPS::HubbardU1>::
2283compute_parallel (bool LOAD_MPO, bool SAVE_MPO, string wd, int minNshell, int maxNshell, bool SAVE)
2284{
2285 make_Hamiltonian(LOAD_MPO,SAVE_MPO,wd);
2286
2287 g.resize(2*L+1);
2288 var.resize(2*L+1);
2289 Mmax.resize(2*L+1);
2290 avgN.resize(2*L+1);
2291 n.resize(2*L+1);
2292 Jlabel.resize(2*L+1);
2293 Jindex.resize(2*L+1);
2294 energies.resize(2*L+1);
2295 energies_free.resize(2*L+1);
2296
2297 int Nshelllast = (maxNshell==-1)? 2*L : maxNshell;
2298
2299 #pragma omp parallel for schedule(dynamic)
2300 for (int iNshell=minNshell; iNshell<=Nshelllast; iNshell+=2)
2301 {
2302 int jmaxNshell = (iNshell==Nshelllast)? 0:1;
2303 for (int jNshell=0; jNshell<=jmaxNshell; ++jNshell)
2304 {
2305 int Nshell = iNshell+jNshell;
2306 int A = Z+Nclosed+Nshell;
2307
2308 #pragma omp critical
2309 {
2310 lout << termcolor::bold << "A=" << A << ", Z=" << Z << ", N=" << Nclosed+Nshell << ", Nshell=" << Nshell << "/" << 2*L << termcolor::reset << endl;
2311 }
2312
2313 auto [gres, Jlabel_res, Jindex_res] = calc_gs(Nshell, LANCZOS::EDGE::GROUND, false, DMRG::VERBOSITY::SILENT);
2314 g[Nshell] = gres;
2315 Jlabel[Nshell] = Jlabel_res;
2316 Jindex[Nshell] = Jindex_res;
2317
2318 #pragma omp critical
2319 {
2320 lout << "Nshell=" << Nshell << " done!, E0=" << g[Nshell].energy << endl;
2321 }
2322
2323 if (Nshell != 0 and Nshell != 2*L)
2324 {
2325 if (Nshell%2 == 0)
2326 {
2327 var[Nshell] = abs(avg(g[Nshell].state,H,H,g[Nshell].state)-pow(g[Nshell].energy,2))/L;
2328 }
2329 else
2330 {
2331 var[Nshell] = abs(avg(g[Nshell].state,Hodd[Jindex[Nshell]],Hodd[Jindex[Nshell]],g[Nshell].state)-pow(g[Nshell].energy,2))/L;
2332 }
2333 }
2334 else
2335 {
2336 var[Nshell] = 0.;
2337 }
2338
2339 Mmax[Nshell] = g[Nshell].state.calc_Mmax();
2340
2341 energies(Nshell) = g[Nshell].energy;
2342 energies_free(Nshell) = onsite_free.head(Nshell).sum();
2343 if (Nshell > 1 and Nshell<2*L-1) assert(g[Nshell].energy < energies_free(Nshell));
2344
2345 n[Nshell].resize(Nlev);
2346 avgN[Nshell].resize(Nlev);
2347
2348 for (int j=0; j<Nlev; ++j)
2349 {
2350 avgN[Nshell](j) = 0.;
2351 for (int i=0; i<deg(j); ++i)
2352 {
2353 if (Nshell%2==0)
2354 {
2355 avgN[Nshell](j) += avg(g[Nshell].state, H.n(offset(j)+i), g[Nshell].state);
2356 }
2357 else
2358 {
2359 avgN[Nshell](j) += avg(g[Nshell].state, Hodd[Jindex[Nshell]].n(offset(j)+i), g[Nshell].state);
2360 }
2361 }
2362 n[Nshell](j) = avgN[Nshell](j)/(2.*deg(j));
2363 }
2364 }
2365 }
2366
2367 if (SAVE) save(wd, minNshell, Nshelllast);
2368}
2369
2370template<>
2371void NuclearManager<VMPS::HubbardU1xU1>::
2372compute_parallel (bool LOAD_MPO, bool SAVE_MPO, string wd, int minNshell, int maxNshell, bool SAVE)
2373{
2374 make_Hamiltonian(LOAD_MPO,SAVE_MPO,wd);
2375
2376 g.resize(2*L+1);
2377 var.resize(2*L+1);
2378 Mmax.resize(2*L+1);
2379 avgN.resize(2*L+1);
2380 n.resize(2*L+1);
2381 Jlabel.resize(2*L+1);
2382 Jindex.resize(2*L+1);
2383 energies.resize(2*L+1);
2384 energies_free.resize(2*L+1);
2385
2386 int Nshelllast = (maxNshell==-1)? 2*L : maxNshell;
2387
2388 #pragma omp parallel for schedule(dynamic)
2389 for (int iNshell=minNshell; iNshell<=Nshelllast; iNshell+=2)
2390 {
2391 int jmaxNshell = (iNshell==Nshelllast)? 0:1;
2392 for (int jNshell=0; jNshell<=jmaxNshell; ++jNshell)
2393 {
2394 int Nshell = iNshell+jNshell;
2395 int A = Z+Nclosed+Nshell;
2396
2397 #pragma omp critical
2398 {
2399 lout << termcolor::bold << "A=" << A << ", Z=" << Z << ", N=" << Nclosed+Nshell << ", Nshell=" << Nshell << "/" << 2*L << termcolor::reset << endl;
2400 }
2401
2402 auto [gres, Jlabel_res, Jindex_res] = calc_gs(Nshell, LANCZOS::EDGE::GROUND, false, DMRG::VERBOSITY::SILENT);
2403 g[Nshell] = gres;
2404 Jlabel[Nshell] = Jlabel_res;
2405 Jindex[Nshell] = Jindex_res;
2406
2407 #pragma omp critical
2408 {
2409 lout << "Nshell=" << Nshell << " done!, E0=" << g[Nshell].energy << endl;
2410 }
2411
2412 if (Nshell != 0 and Nshell != 2*L)
2413 {
2414 if (Nshell%2 == 0)
2415 {
2416 var[Nshell] = abs(avg(g[Nshell].state,H,H,g[Nshell].state)-pow(g[Nshell].energy,2))/L;
2417 }
2418 else
2419 {
2420 var[Nshell] = abs(avg(g[Nshell].state,Hodd[Jindex[Nshell]],Hodd[Jindex[Nshell]],g[Nshell].state)-pow(g[Nshell].energy,2))/L;
2421 }
2422 }
2423 else
2424 {
2425 var[Nshell] = 0.;
2426 }
2427
2428 Mmax[Nshell] = g[Nshell].state.calc_Mmax();
2429
2430 energies(Nshell) = g[Nshell].energy;
2431 energies_free(Nshell) = onsite_free.head(Nshell).sum();
2432 if (Nshell > 1 and Nshell<2*L-1) assert(g[Nshell].energy < energies_free(Nshell));
2433
2434 n[Nshell].resize(Nlev);
2435 avgN[Nshell].resize(Nlev);
2436
2437 for (int j=0; j<Nlev; ++j)
2438 {
2439 avgN[Nshell](j) = 0.;
2440 for (int i=0; i<deg(j); ++i)
2441 {
2442 if (Nshell%2==0)
2443 {
2444 avgN[Nshell](j) += avg(g[Nshell].state, H.n(offset(j)+i), g[Nshell].state);
2445 }
2446 else
2447 {
2448 avgN[Nshell](j) += avg(g[Nshell].state, Hodd[Jindex[Nshell]].n(offset(j)+i), g[Nshell].state);
2449 }
2450 }
2451 n[Nshell](j) = avgN[Nshell](j)/(2.*deg(j));
2452 }
2453 }
2454 }
2455
2456 if (SAVE) save(wd, minNshell, Nshelllast);
2457}
2458
2459template<>
2460void NuclearManager<VMPS::HeisenbergU1>::
2461compute_parallel (bool LOAD_MPO, bool SAVE_MPO, string wd, int minNshell, int maxNshell, bool SAVE)
2462{
2463 make_Hamiltonian(LOAD_MPO,SAVE_MPO,wd);
2464
2465 g.resize(2*L+1);
2466 var.resize(2*L+1);
2467 Mmax.resize(2*L+1);
2468 avgN.resize(2*L+1);
2469 n.resize(2*L+1);
2470 Jlabel.resize(2*L+1);
2471 Jindex.resize(2*L+1);
2472 energies.resize(2*L+1);
2473 energies_free.resize(2*L+1);
2474
2475 int Nshelllast = (maxNshell==-1)? 2*L : maxNshell;
2476
2477 #pragma omp parallel for schedule(dynamic)
2478 for (int Nshell=minNshell; Nshell<=Nshelllast; Nshell+=2)
2479 {
2480 int A = Z+Nclosed+Nshell;
2481
2482 #pragma omp critical
2483 {
2484 lout << termcolor::bold << "A=" << A << ", Z=" << Z << ", N=" << Nclosed+Nshell << ", Nshell=" << Nshell << "/" << 2*L << termcolor::reset << endl;
2485 }
2486
2487 auto [gres, Jlabel_res, Jindex_res] = calc_gs(Nshell, LANCZOS::EDGE::GROUND, false, DMRG::VERBOSITY::SILENT);
2488 g[Nshell] = gres;
2489 Jlabel[Nshell] = Jlabel_res;
2490 Jindex[Nshell] = Jindex_res;
2491
2492 if (Nshell != 0 and Nshell != 2*L)
2493 {
2494 var[Nshell] = abs(avg(g[Nshell].state,H,H,g[Nshell].state)-pow(g[Nshell].energy,2))/L;
2495 }
2496 else
2497 {
2498 var[Nshell] = 0.;
2499 }
2500
2501 Mmax[Nshell] = g[Nshell].state.calc_Mmax();
2502
2503 energies(Nshell) = g[Nshell].energy;
2504 energies_free(Nshell) = onsite_free.head(Nshell).sum();
2505 if (Nshell > 1 and Nshell<2*L-1) assert(g[Nshell].energy < energies_free(Nshell));
2506
2507 n[Nshell].resize(Nlev);
2508 avgN[Nshell].resize(Nlev);
2509
2510 for (int j=0; j<Nlev; ++j)
2511 {
2512 avgN[Nshell](j) = 0.;
2513 for (int i=0; i<deg(j); ++i)
2514 {
2515 avgN[Nshell](j) += 2.*avg(g[Nshell].state, H.Sz(offset(j)+i), g[Nshell].state)+1.;
2516 }
2517 n[Nshell](j) = avgN[Nshell](j)/(2.*deg(j));
2518 }
2519 }
2520
2521 if (SAVE) save(wd, minNshell, Nshelllast, 2);
2522}
2523
2524template<typename MODEL>
2525void NuclearManager<MODEL>::
2526save (string wd, int Nfrst, int Nlast, int dNshell) const
2527{
2528 HDF5Interface target(wd+outfile, WRITE);
2529 target.save_vector(eps0,"eps0","");
2530 target.save_vector(deg,"deg","");
2531
2532 for (int Nshell=Nfrst; Nshell<=Nlast; Nshell+=dNshell)
2533 {
2534 target.create_group(make_string(Nshell));
2535 target.save_scalar(g[Nshell].energy,"E0",make_string(Nshell));
2536 target.save_scalar(energies_free(Nshell),"E0free",make_string(Nshell));
2537 target.save_vector(avgN[Nshell],"avgN",make_string(Nshell));
2538 target.save_vector(n[Nshell],"n",make_string(Nshell));
2539 target.save_scalar(var[Nshell],"var",make_string(Nshell));
2540 target.save_scalar(Mmax[Nshell],"Mmax",make_string(Nshell));
2541 target.save_scalar(Jlabel[Nshell],"Jlabel",make_string(Nshell));
2542 target.save_scalar(Jindex[Nshell],"Jindex",make_string(Nshell));
2543 }
2544
2545 MatrixXd Delta3(0,2);
2546 for (int Nshell=1, i=0; Nshell<=2*L-1; ++Nshell, ++i)
2547 {
2548 double DeltaVal = 0.5*(2.*energies(Nshell)-energies(Nshell-1)-energies(Nshell+1));
2549
2550// lout << "Nn=" << Nshell+Nclosed << ", Delta3=" << DeltaVal << endl;
2551 Delta3.conservativeResize(Delta3.rows()+1,Delta3.cols());
2552 Delta3(i,0) = Nshell;
2553 Delta3(i,1) = DeltaVal;
2554 }
2555 target.save_matrix(Delta3,"Delta3");
2556
2557 MatrixXd Delta3free(0,2);
2558 for (int Nshell=1, i=0; Nshell<=2*L-1; ++Nshell, ++i)
2559 {
2560 double DeltaVal = 0.5*(2.*energies_free(Nshell)-energies_free(Nshell-1)-energies_free(Nshell+1));
2561
2562 Delta3free.conservativeResize(Delta3free.rows()+1,Delta3free.cols());
2563 Delta3free(i,0) = Nshell;
2564 Delta3free(i,1) = DeltaVal;
2565 }
2566 target.save_matrix(Delta3free,"Delta3free");
2567
2568 MatrixXd Delta3b(0,2);
2569 for (int Nshell=2, i=0; Nshell<=2*L; ++Nshell, ++i)
2570 {
2571 double DeltaVal = 0.5*(energies(Nshell)-2.*energies(Nshell-1)+energies(Nshell-2));
2572
2573// lout << "Nn=" << Nshell+Nclosed << ", Delta3b=" << DeltaVal << endl;
2574 Delta3b.conservativeResize(Delta3b.rows()+1,Delta3b.cols());
2575 Delta3b(i,0) = Nshell;
2576 Delta3b(i,1) = DeltaVal;
2577 }
2578 target.save_matrix(Delta3b,"Delta3b");
2579
2580 MatrixXd Delta4(0,2);
2581 for (int Nshell=2, i=0; Nshell<=2*L-1; ++Nshell, ++i)
2582 {
2583 double DeltaVal = 0.25*(energies(Nshell-2)-3.*energies(Nshell-1)+3.*energies(Nshell)-energies(Nshell+1));
2584
2585// lout << "Nn=" << Nshell+Nclosed << ", Delta4=" << DeltaVal << endl;
2586 Delta4.conservativeResize(Delta4.rows()+1,Delta4.cols());
2587 Delta4(i,0) = Nshell;
2588 Delta4(i,1) = DeltaVal;
2589 }
2590 target.save_matrix(Delta4,"Delta4");
2591
2592 MatrixXd Delta5(0,2);
2593 for (int Nshell=2, i=0; Nshell<=2*L-2; ++Nshell, ++i)
2594 {
2595 double DeltaVal = -0.125*(6.*energies(Nshell-2)-4.*energies(Nshell+1)-4.*energies(Nshell-1)+energies(Nshell+2)+energies(Nshell-2));
2596
2597// lout << "Nn=" << Nshell+Nclosed << ", Delta5=" << DeltaVal << endl;
2598 Delta5.conservativeResize(Delta5.rows()+1,Delta5.cols());
2599 Delta5(i,0) = Nshell;
2600 Delta5(i,1) = DeltaVal;
2601 }
2602 target.save_matrix(Delta5,"Delta5");
2603
2604 lout << termcolor::green << "saved: " << wd+outfile << termcolor::reset << endl;
2605
2606 target.close();
2607}
2608
2609#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
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
diff
Mpo< Symmetry, Scalar > diff(const Mpo< Symmetry, Scalar > &H1, const Mpo< Symmetry, Scalar > &H2, DMRG::VERBOSITY::OPTION VERBOSITY=DMRG::VERBOSITY::SILENT)
Definition: DmrgLinearAlgebra.h:914
SPIN_INDEX
SPIN_INDEX
Definition: DmrgTypedefs.h:36
DN
@ DN
Definition: DmrgTypedefs.h:38
UP
@ UP
Definition: DmrgTypedefs.h:37
HUBBARD
@ HUBBARD
Definition: DmrgTypedefs.h:96
HEISENBERG
@ HEISENBERG
Definition: DmrgTypedefs.h:96
ODD
@ ODD
Definition: DmrgTypedefs.h:128
OPEN
@ OPEN
Definition: DmrgTypedefs.h:163
A
@ A
Definition: DmrgTypedefs.h:130
load_nuclearData
void load_nuclearData(string el, string set, string rootfolder, NuclearInfo &info, double &V0, MatrixXd &Vij, vector< MatrixXd > &Vijnocgc, VectorXi &deg, VectorXd &eps0, bool &REF)
Definition: NuclearManager.h:79
cgc
double cgc(int j1x2, int j2x2, int Jx2, int m1x2, int m2x2, int Mx2)
Definition: NuclearManager.h:51
MpoTerms::set_verbosity
void set_verbosity(const DMRG::VERBOSITY::OPTION VERB_in)
Definition: MpoTerms.h:881
Mpo
Definition: Mpo.h:40
NuclearManager::get_Gloc
VectorXd get_Gloc() const
Definition: NuclearManager.h:375
NuclearManager::make_cdagj
Mpo< typename MODEL::Symmetry > make_cdagj(int j, const MODEL &H_input) const
NuclearManager::g
vector< Eigenstate< typename MODEL::StateXd > > g
Definition: NuclearManager.h:412
NuclearManager::FlipChart
std::array< MatrixXd, 2 > FlipChart
Definition: NuclearManager.h:398
NuclearManager::make_fullHamiltonian
void make_fullHamiltonian(bool LOAD=false, bool SAVE=false, string wd="./")
Definition: NuclearManager.h:951
NuclearManager::NuclearManager
NuclearManager(int Nclosed_input, int Nsingle_input, int Z_input, const string &filename, double V0_input, bool REF_input=false)
Definition: NuclearManager.h:317
NuclearManager::get_occ
ArrayXd get_occ() const
Definition: NuclearManager.h:594
NuclearManager::get_levelsign
ArrayXd get_levelsign(int level)
Definition: NuclearManager.h:361
NuclearManager::sign_per_level
vector< ArrayXd > sign_per_level
Definition: NuclearManager.h:403
NuclearManager::get_offset
ArrayXi get_offset() const
Definition: NuclearManager.h:358
NuclearManager::length
size_t length() const
Definition: NuclearManager.h:359
NuclearManager::get_sign2d
ArrayXXd get_sign2d() const
Definition: NuclearManager.h:378
NuclearManager::get_onsite_free
VectorXd get_onsite_free() const
Definition: NuclearManager.h:362
NuclearManager::get_FlipChart
std::array< MatrixXd, 2 > get_FlipChart() const
Definition: NuclearManager.h:377
NuclearManager::get_eps0
ArrayXd get_eps0() const
Definition: NuclearManager.h:356
NuclearManager::generate_couplingList
vector< tuple< orbital, orbital, orbital, orbital, double, orbinfo > > generate_couplingList(int J) const
Definition: NuclearManager.h:1150
NuclearManager::get_Ghop
MatrixXd get_Ghop() const
Definition: NuclearManager.h:373
NuclearManager::jlist
vector< boost::rational< int > > jlist
Definition: NuclearManager.h:399
NuclearManager::avgN
vector< VectorXd > avgN
Definition: NuclearManager.h:415
NuclearManager::get_sign
ArrayXd get_sign()
Definition: NuclearManager.h:360
NuclearManager::make_Hamiltonian
void make_Hamiltonian(bool LOAD=false, bool SAVE=false, string wd="./", bool ODD=true)
Definition: NuclearManager.h:775
NuclearManager::mlist
vector< boost::rational< int > > mlist
Definition: NuclearManager.h:399
NuclearManager::cdagj
vector< Mpo< typename MODEL::Symmetry > > cdagj
Definition: NuclearManager.h:427
NuclearManager::Sn_nref
ArrayXd Sn_nref(int Nshell) const
Definition: NuclearManager.h:662
NuclearManager::NuclearManager
NuclearManager(int Nclosed_input, int Nsingle_input, int Z_input, const ArrayXi &deg_input, const vector< string > &labels_input, const ArrayXd &eps0_input, const ArrayXXd &Vij_input, double V0_input=1., bool REF_input=false, string PARAM_input="Seminole")
Definition: NuclearManager.h:303
NuclearManager::Sn_Nref
ArrayXd Sn_Nref(int Nshell) const
Definition: NuclearManager.h:676
NuclearManager::save
void save(string wd, int Nfrst=0, int Nlast=-1, int dNshell=1) const
Definition: NuclearManager.h:2526
NuclearManager::Sn_Sref14
ArrayXd Sn_Sref14() const
Definition: NuclearManager.h:762
NuclearManager::get_onsite
VectorXd get_onsite() const
Definition: NuclearManager.h:374
NuclearManager::onsite_free
VectorXd onsite_free
Definition: NuclearManager.h:407
NuclearManager::n
vector< VectorXd > n
Definition: NuclearManager.h:415
NuclearManager::get_g
const Eigenstate< typename MODEL::StateXd > & get_g(int Nshell) const
Definition: NuclearManager.h:372
NuclearManager::energies_free
VectorXd energies_free
Definition: NuclearManager.h:408
NuclearManager::compute_parallel
void compute_parallel(bool LOADv=false, bool SAVE_MPO=false, string wd="./", int minNshell=0, int maxNshell=-1, bool SAVE=true)
NuclearManager::Sn_Pref
ArrayXd Sn_Pref(int Nshell) const
Definition: NuclearManager.h:716
NuclearManager::Mmax
vector< int > Mmax
Definition: NuclearManager.h:414
NuclearManager::get_mfactor
VectorXi get_mfactor() const
Definition: NuclearManager.h:376
NuclearManager::orblist
vector< int > orblist
Definition: NuclearManager.h:400
NuclearManager::Jindex
vector< int > Jindex
Definition: NuclearManager.h:417
NuclearManager::var
vector< double > var
Definition: NuclearManager.h:413
NuclearManager::Sn_Eref
double Sn_Eref(int N) const
Definition: NuclearManager.h:606
NuclearManager::calc_gs_full
Eigenstate< typename MODEL::StateXd > calc_gs_full(int Nshell, LANCZOS::EDGE::OPTION EDGE=LANCZOS::EDGE::GROUND, bool CALC_VAR=true, DMRG::VERBOSITY::OPTION VERB=DMRG::VERBOSITY::ON_EXIT) const
Definition: NuclearManager.h:1523
NuclearManager::get_Hfull
MODEL get_Hfull() const
Definition: NuclearManager.h:364
NuclearManager::get_H
MODEL get_H() const
Definition: NuclearManager.h:363
NuclearManager::Hodd
vector< MODEL > Hodd
Definition: NuclearManager.h:420
NuclearManager::Hperturb
vector< MODEL > Hperturb
Definition: NuclearManager.h:423
NuclearManager::calc_gs
tuple< Eigenstate< typename MODEL::StateXd >, string, int > calc_gs(int Nshell, LANCZOS::EDGE::OPTION EDGE=LANCZOS::EDGE::GROUND, bool CALC_VAR=true, DMRG::VERBOSITY::OPTION VERB=DMRG::VERBOSITY::ON_EXIT) const
NuclearManager::labels
vector< string > labels
Definition: NuclearManager.h:394
NuclearManager::compute
void compute(bool LOAD_MPO=false, bool SAVE_MPO=false, string wd="./", int minNshell=0, int maxNshell=-1, bool SAVE=true)
NuclearManager::Jlabel
vector< string > Jlabel
Definition: NuclearManager.h:416
NuclearManager::perturb_sweeps
vector< double > perturb_sweeps
Definition: NuclearManager.h:424
NuclearManager::get_deg
ArrayXi get_deg() const
Definition: NuclearManager.h:357
NuclearManager::make_cdagj
vector< Mpo< typename MODEL::Symmetry > > make_cdagj(const MODEL &H_input) const
VMPS::HeisenbergU1
Heisenberg Model.
Definition: HeisenbergU1.h:41
VMPS::HubbardU1::singlet
static qarray< 1 > singlet(int N)
Definition: HubbardU1.h:41
VMPS::HubbardU1xU1
Hubbard model with U(1) symmetries. MPO representation of the Hubbard model.
Definition: HubbardU1xU1.h:52
VMPS::HubbardU1xU1::singlet
static qarray< 2 > singlet(int N=0, int L=0)
Definition: HubbardU1xU1.h:78
DMRG::CONTROL::DYN::Mincr_per
function< size_t(size_t)> Mincr_per
Definition: DmrgTypedefs.h:403
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::CONVTEST
DMRG::CONVTEST::OPTION CONVTEST
Definition: DmrgTypedefs.h:387
NuclearInfo::info
string info() const
Definition: NuclearManager.h:60
NuclearInfo::P
Nuclear::PARTICLE P
Definition: NuclearManager.h:72
NuclearInfo::levels
vector< string > levels
Definition: NuclearManager.h:74
NuclearInfo::levelstring
string levelstring
Definition: NuclearManager.h:75
jEigenstate::eigenstate
Eigenstate< VectorType > eigenstate
Definition: NuclearManager.h:17
orbinfo::operator==
bool operator==(const orbinfo &b) const
Definition: NuclearManager.h:41
orbinfo::label
std::array< string, 4 > label
Definition: NuclearManager.h:48
orbinfo::jx2
std::array< int, 4 > jx2
Definition: NuclearManager.h:47
orbinfo::lev
std::array< int, 4 > lev
Definition: NuclearManager.h:46
orbital::orbital
orbital(const int &r)
Definition: NuclearManager.h:24
orbital::s
SPIN_INDEX s
Definition: NuclearManager.h:36
orbital::operator==
bool operator==(const orbital &b) const
Definition: NuclearManager.h:30
qarray
Definition: qarray.h:26