VMPS++
Loading...
Searching...
No Matches
FermionBase.h
Go to the documentation of this file.
1#ifndef FERMIONBASE_H_
2#define FERMIONBASE_H_
3
4#include "DmrgTypedefs.h"
5
6#include "tensors/SiteOperatorQ.h"
7#include "sites/FermionSite.h"
8#include <boost/dynamic_bitset.hpp>
9#include <Eigen/Core>
10#include <unsupported/Eigen/MatrixFunctions>
12
13#include "symmetry/kind_dummies.h"
14#include "DmrgTypedefs.h" // for SPIN_INDEX, SPINOP_LABEL
15#include "tensors/SiteOperator.h"
16#include "DmrgExternal.h" // for posmod
17
18//Note: Don't put a name in this documentation with \class .. because doxygen gets confused with template symbols
25template<typename Symmetry_>
26class FermionBase : public FermionSite<Symmetry_>
27{
28 typedef Eigen::Index Index;
29 typedef double Scalar;
30
31public:
32
33 typedef Symmetry_ Symmetry;
34 typedef SiteOperatorQ<Symmetry,Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> > OperatorType;
35 typedef SiteOperatorQ<Symmetry,Eigen::Matrix<complex<Scalar>,Eigen::Dynamic,Eigen::Dynamic> > ComplexOperatorType;
36 typedef typename Symmetry::qType qType;
37
38 FermionBase(){};
39
44 FermionBase (std::size_t L_input, bool REMOVE_DOUBLE=false, bool REMVOVE_EMPTY=false, bool REMOVE_UP=false, bool REMOVE_DN=false, int mfactor=1, int k_input=0);
45
47 inline Index dim() const {return static_cast<Index>(N_states);}
48
50 inline std::size_t orbitals() const {return N_orbitals;}
51
53
65 template<class Dummy = Symmetry>
66 typename std::enable_if<Dummy::IS_SPIN_SU2() and !Dummy::IS_CHARGE_SU2(),OperatorType>::type c (size_t orbital=0) const;
67
68 template<class Dummy = Symmetry>
69 typename std::enable_if<!Dummy::IS_SPIN_SU2() and !Dummy::IS_CHARGE_SU2(),OperatorType>::type c (SPIN_INDEX sigma, size_t orbital=0) const;
70
71 template<class Dummy = Symmetry>
72 typename std::enable_if<Dummy::IS_CHARGE_SU2() and !Dummy::IS_SPIN_SU2(),OperatorType>::type c (SPIN_INDEX sigma, SUB_LATTICE G, size_t orbital=0) const;
73
74 template<class Dummy = Symmetry>
75 typename std::enable_if<Dummy::IS_CHARGE_SU2() and Dummy::IS_SPIN_SU2(),OperatorType>::type c (SUB_LATTICE G, size_t orbital=0) const;
76
90 template<class Dummy = Symmetry>
91 typename std::enable_if<Dummy::IS_SPIN_SU2() and !Dummy::IS_CHARGE_SU2(),OperatorType>::type cdag (size_t orbital=0) const;
92
93 template<class Dummy = Symmetry>
94 typename std::enable_if<!Dummy::IS_SPIN_SU2() and !Dummy::IS_CHARGE_SU2(),OperatorType>::type cdag (SPIN_INDEX sigma, size_t orbital=0) const;
95
96 template<class Dummy = Symmetry>
97 typename std::enable_if<Dummy::IS_CHARGE_SU2() and !Dummy::IS_SPIN_SU2(),OperatorType>::type cdag (SPIN_INDEX sigma, SUB_LATTICE G, size_t orbital=0) const;
98
99 template<class Dummy = Symmetry>
100 typename std::enable_if<Dummy::IS_CHARGE_SU2() and Dummy::IS_SPIN_SU2(),OperatorType>::type cdag (SUB_LATTICE G, size_t orbital=0) const;
101
107 OperatorType sign (std::size_t orb1=0, std::size_t orb2=0) const;
108
113 OperatorType sign_local (std::size_t orbital=0) const;
114
119 template<class Dummy = Symmetry>
120 typename std::enable_if<true,OperatorType>::type n (std::size_t orbital=0) const;
121
122 template<class Dummy = Symmetry>
123 typename std::enable_if<!Dummy::IS_SPIN_SU2(),OperatorType>::type n (SPIN_INDEX sigma, size_t orbital=0) const;
124
129 template<class Dummy = Symmetry>
130 typename std::enable_if<!Dummy::IS_CHARGE_SU2(),OperatorType>::type d (std::size_t orbital=0) const;
131
136 OperatorType ns (std::size_t orbital=0) const;
137
142 OperatorType nh (std::size_t orbital=0) const;
144
146
150 template<class Dummy = Symmetry>
151 typename std::enable_if<Dummy::IS_SPIN_SU2(),OperatorType>::type S (size_t orbital=0) const;
152
157 template<class Dummy = Symmetry>
158 typename std::enable_if<Dummy::IS_SPIN_SU2(),OperatorType>::type Sdag (size_t orbital=0) const;
159
160 template<class Dummy = Symmetry>
161 typename std::enable_if<!Dummy::IS_SPIN_SU2(),OperatorType>::type Sz (size_t orbital=0) const;
162
163 template<class Dummy = Symmetry>
164 typename std::enable_if<!Dummy::IS_SPIN_SU2(),ComplexOperatorType>::type exp_ipiSz (size_t orbital=0) const;
165
166 template<class Dummy = Symmetry>
167 typename std::enable_if<!Dummy::IS_SPIN_SU2(),OperatorType>::type Sp (size_t orbital=0) const;
168
169 template<class Dummy = Symmetry>
170 typename std::enable_if<!Dummy::IS_SPIN_SU2(),OperatorType>::type Sm (size_t orbital=0) const;
171
172 template<class Dummy = Symmetry>
173 typename std::enable_if<Dummy::NO_SPIN_SYM(),OperatorType>::type Sx (size_t orbital=0) const;
174
175 template<class Dummy = Symmetry>
176 typename std::enable_if<Dummy::NO_SPIN_SYM(),OperatorType>::type iSy (size_t orbital=0) const;
177
178 template<class Dummy = Symmetry>
179 typename std::enable_if<!Dummy::IS_SPIN_SU2(),OperatorType>::type Scomp (SPINOP_LABEL Sa, int orbital) const
180 {
181 assert(Sa != SY);
182 OperatorType out;
183 if constexpr (Dummy::NO_SPIN_SYM())
184 {
185 if (Sa==SX) {out = Sx(orbital);}
186 else if (Sa==iSY) {out = iSy(orbital);}
187 else if (Sa==SZ) {out = Sz(orbital);}
188 else if (Sa==SP) {out = Sp(orbital);}
189 else if (Sa==SM) {out = Sm(orbital);}
190 }
191 else
192 {
193 if (Sa==SZ) {out = Sz(orbital);}
194 else if (Sa==SP) {out = Sp(orbital);}
195 else if (Sa==SM) {out = Sm(orbital);}
196 }
197 return out;
198 };
200
201 template<class Dummy = Symmetry>
202 typename std::enable_if<!Dummy::IS_SPIN_SU2(),SiteOperatorQ<Symmetry,Eigen::MatrixXcd> >::type Rcomp (SPINOP_LABEL Sa, int orbital) const;
203
205
209 template<class Dummy = Symmetry>
210 typename std::enable_if<Dummy::IS_CHARGE_SU2(),OperatorType>::type T (size_t orbital=0) const;
211
216 template<class Dummy = Symmetry>
217 typename std::enable_if<Dummy::IS_CHARGE_SU2(),OperatorType>::type Tdag (size_t orbital=0) const;
218
223 template<class Dummy = Symmetry>
224 typename std::enable_if<!Dummy::IS_CHARGE_SU2(),OperatorType>::type Tz (size_t orbital=0) const;
225
226 template<class Dummy = Symmetry>
227 typename std::enable_if<!Dummy::IS_CHARGE_SU2(),OperatorType>::type tz (size_t orbital=0) const;
228
233 template<class Dummy = Symmetry>
234 typename std::enable_if<Dummy::NO_CHARGE_SYM(),OperatorType>::type Tx (size_t orbital=0, SUB_LATTICE G=A) const;
235
240 template<class Dummy = Symmetry>
241 typename std::enable_if<Dummy::NO_CHARGE_SYM(),OperatorType>::type iTy (size_t orbital=0, SUB_LATTICE G=A) const;
242
243 template<class Dummy = Symmetry>
244 typename std::enable_if<!Dummy::IS_CHARGE_SU2(),OperatorType>::type Tp (size_t orbital=0, SUB_LATTICE G=A) const
245 {
246 double factor = static_cast<double>(static_cast<int>(G));
247 return factor*cc(orbital);
248 };
249
250 template<class Dummy = Symmetry>
251 typename std::enable_if<!Dummy::IS_CHARGE_SU2(),OperatorType>::type Tm (size_t orbital=0, SUB_LATTICE G=A) const
252 {
253 double factor = static_cast<double>(static_cast<int>(G));
254 return factor*cdagcdag(orbital);
255 };
256
258
262 template<class Dummy = Symmetry>
263 typename std::enable_if<!Dummy::IS_CHARGE_SU2(),OperatorType>::type cc (size_t orbital=0) const;
264
269 template<class Dummy = Symmetry>
270 typename std::enable_if<!Dummy::IS_CHARGE_SU2(),OperatorType>::type cdagcdag (size_t orbital=0) const;
272
274 OperatorType Id (std::size_t orbital=0) const;
275
277 ArrayXd ZeroField() const { return ArrayXd::Zero(N_orbitals); }
278
280 ArrayXXd ZeroHopping() const { return ArrayXXd::Zero(N_orbitals,N_orbitals); }
281
293 template<typename Scalar_>
294 SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> > HubbardHamiltonian (const Array<Scalar_,Dynamic,1> &Uph,
295 const Array<Scalar_,Dynamic,Dynamic> &t,
296 const Array<Scalar_,Dynamic,Dynamic> &V,
297 const Array<Scalar_,Dynamic,Dynamic> &J) const;
298
299 template<typename Scalar_, typename Dummy = Symmetry>
300 typename std::enable_if<Dummy::ABELIAN,SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> > >::type
301 HubbardHamiltonian (const Array<Scalar_,Dynamic,1> &U,
302 const Array<Scalar_,Dynamic,1> &Uph,
303 const Array<Scalar_,Dynamic,1> &Eorb,
304 const Array<Scalar_,Dynamic,1> &Bz,
305 const Array<Scalar_,Dynamic,Dynamic> &t,
306 const Array<Scalar_,Dynamic,Dynamic> &V,
307 const Array<Scalar_,Dynamic,Dynamic> &Vz,
308 const Array<Scalar_,Dynamic,Dynamic> &Vxy,
309 const Array<Scalar_,Dynamic,Dynamic> &Jz,
310 const Array<Scalar_,Dynamic,Dynamic> &Jxy,
311 const Array<Scalar_,Dynamic,Dynamic> &C) const;
312
313 template<typename Scalar_, typename Dummy = Symmetry>
314 typename std::enable_if<Dummy::IS_SPIN_SU2() and !Dummy::IS_CHARGE_SU2(),SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> > >::type
315 HubbardHamiltonian (const Array<Scalar_,Dynamic,1> &U,
316 const Array<Scalar_,Dynamic,1> &Uph,
317 const Array<Scalar_,Dynamic,1> &Eorb,
318 const Array<Scalar_,Dynamic,Dynamic> &t,
319 const Array<Scalar_,Dynamic,Dynamic> &V,
320 const Array<Scalar_,Dynamic,Dynamic> &Vz,
321 const Array<Scalar_,Dynamic,Dynamic> &Vxy,
322 const Array<Scalar_,Dynamic,Dynamic> &J) const;
323
324 template<typename Scalar_, typename Dummy = Symmetry>
325 typename std::enable_if<Dummy::NO_SPIN_SYM() and Dummy::IS_CHARGE_SU2(),SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> > >::type
326 HubbardHamiltonian (const Array<Scalar_,Dynamic,1> &Uph,
327 const Array<Scalar_,Dynamic,Dynamic> &t,
328 const Array<Scalar_,Dynamic,Dynamic> &V,
329 const Array<Scalar_,Dynamic,Dynamic> &Jz,
330 const Array<Scalar_,Dynamic,Dynamic> &Jxy,
331 const Array<Scalar_,Dynamic,1> &Bz,
332 const Array<Scalar_,Dynamic,1> &Bx) const;
333
334 template<typename Scalar_, typename Dummy = Symmetry>
335 typename std::enable_if<Dummy::IS_SPIN_U1() and Dummy::IS_CHARGE_SU2(),SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> > >::type
336 HubbardHamiltonian (const Array<Scalar_,Dynamic,1> &Uph,
337 const Array<Scalar_,Dynamic,Dynamic> &t,
338 const Array<Scalar_,Dynamic,Dynamic> &V,
339 const Array<Scalar_,Dynamic,Dynamic> &Jz,
340 const Array<Scalar_,Dynamic,Dynamic> &Jxy,
341 const Array<Scalar_,Dynamic,1> &Bz) const;
342
343 template<typename Scalar_, typename Dummy = Symmetry>
344 typename std::enable_if<Dummy::NO_SPIN_SYM(),SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> >>::type coupling_Bx (const Array<double,Dynamic,1> &Bx) const;
345
346 template<typename Dummy = Symmetry>
347 typename std::enable_if<Dummy::NO_SPIN_SYM(),SiteOperatorQ<Symmetry_,Eigen::Matrix<complex<double>,-1,-1> >>::type coupling_By (const Array<double,Dynamic,1> &By) const;
348
349 template<typename Scalar_, typename Dummy = Symmetry>
350 typename std::enable_if<Dummy::IS_TRIVIAL,SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> >>::type coupling_singleFermion (const Array<double,Dynamic,1> &Fp) const;
351
352 template<typename Scalar_, typename Dummy = Symmetry>
353 typename std::enable_if<Dummy::NO_SPIN_SYM(),SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> >>::type coupling_XYZspin (const Array<double,Dynamic,Dynamic> &Jx,
354 const Array<double,Dynamic,Dynamic> &Jy,
355 const Array<double,Dynamic,Dynamic> &Jz) const;
356
358 Qbasis<Symmetry> get_basis() const { return TensorBasis; }
359
360private:
361
362 OperatorType make_operator (const OperatorType &Op_1s, size_t orbital=0, bool FERMIONIC = false, string label="") const;
363 std::size_t N_orbitals;
364 std::size_t N_states;
365
366 Qbasis<Symmetry> TensorBasis; //Final basis for N_orbital sites
367
368 //operators defined on zero orbitals
369 OperatorType Id_vac, Zero_vac;
370};
371
372template <typename Symmetry_>
373FermionBase<Symmetry_>::
374FermionBase (std::size_t L_input, bool REMOVE_DOUBLE, bool REMVOVE_EMPTY, bool REMOVE_UP, bool REMOVE_DN, int mfactor, int k_input)
375:FermionSite<Symmetry>(REMOVE_DOUBLE, REMVOVE_EMPTY, REMOVE_UP, REMOVE_DN, mfactor, k_input), N_orbitals(L_input)
376{
377 //create basis for zero orbitals
378 typename Symmetry::qType Q=Symmetry::qvacuum();
379 Eigen::Index inner_dim = 1;
380 Qbasis<Symmetry_> vacuum;
381 vacuum.push_back(Q, inner_dim);
382
383 // create operators for zero orbitals
384 Zero_vac = OperatorType(Symmetry::qvacuum(), vacuum);
385 Zero_vac.setZero();
386 Id_vac = OperatorType(Symmetry::qvacuum(), vacuum);
387 Id_vac.setIdentity();
388
389 // create basis for N_orbitals fermionic sites
390 if (N_orbitals == 1) {TensorBasis = this->basis_1s();}
391 else if (N_orbitals == 0) {TensorBasis = vacuum;}
392 else
393 {
394 TensorBasis = this->basis_1s().combine(this->basis_1s());
395 for (std::size_t o=2; o<N_orbitals; o++)
396 {
397 TensorBasis = TensorBasis.combine(this->basis_1s());
398 }
399 }
400
401 N_states = TensorBasis.size();
402}
403
404template<typename Symmetry_>
405SiteOperatorQ<Symmetry_, Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > FermionBase<Symmetry_>::
406make_operator (const OperatorType &Op_1s, size_t orbital, bool FERMIONIC, string label) const
407{
408 OperatorType out;
409 if (N_orbitals == 1) {out = Op_1s; out.label() = label; return out;}
410 else if (N_orbitals == 0) {return Zero_vac;}
411 else
412 {
413 OperatorType stringOp;
414 if (FERMIONIC) {stringOp = this->F_1s();}
415 else {stringOp = this->Id_1s();}
416 bool TOGGLE=false;
417 if (orbital == 0) {out = OperatorType::outerprod(Op_1s,this->Id_1s(),Op_1s.Q()); TOGGLE=true;}
418 else
419 {
420 if (orbital == 1) {out = OperatorType::outerprod(stringOp,Op_1s,Op_1s.Q()); TOGGLE=true;}
421 else {out = OperatorType::outerprod(stringOp,stringOp,Symmetry_::qvacuum());}
422 }
423 for(std::size_t o=2; o<N_orbitals; o++)
424 {
425 if (orbital == o) {out = OperatorType::outerprod(out,Op_1s,Op_1s.Q()); TOGGLE=true; }
426 else if (TOGGLE==false) {out = OperatorType::outerprod(out,stringOp,Symmetry_::qvacuum());}
427 else if (TOGGLE==true) {out = OperatorType::outerprod(out,this->Id_1s(),Op_1s.Q());}
428 }
429 out.label() = label;
430 return out;
431 }
432}
433
434template <typename Symmetry_>
435template <typename Dummy>
436typename std::enable_if<Dummy::IS_SPIN_SU2() and !Dummy::IS_CHARGE_SU2(),SiteOperatorQ<Symmetry_, Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
437c (std::size_t orbital) const
438{
439 return make_operator(this->c_1s(),orbital,PROP::FERMIONIC, "c");
440}
441
442template <typename Symmetry_>
443template <typename Dummy>
444typename std::enable_if<Dummy::IS_SPIN_SU2() and !Dummy::IS_CHARGE_SU2(),SiteOperatorQ<Symmetry_, Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
445cdag (std::size_t orbital) const
446{
447 //return c(orbital).adjoint();
448 return make_operator(this->cdag_1s(),orbital,PROP::FERMIONIC, "c†");
449}
450
451template <typename Symmetry_>
452template <typename Dummy>
453typename std::enable_if<!Dummy::IS_SPIN_SU2() and !Dummy::IS_CHARGE_SU2(),SiteOperatorQ<Symmetry_, Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
454c (SPIN_INDEX sigma, std::size_t orbital) const
455{
456 stringstream ss;
457 ss << "c" << sigma;
458 return make_operator(this->c_1s(sigma),orbital,PROP::FERMIONIC, ss.str());
459}
460
461template <typename Symmetry_>
462template <typename Dummy>
463typename std::enable_if<!Dummy::IS_SPIN_SU2() and !Dummy::IS_CHARGE_SU2(),SiteOperatorQ<Symmetry_, Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
464cdag (SPIN_INDEX sigma, std::size_t orbital) const
465{
466// return c(sigma,orbital).adjoint();
467 stringstream ss;
468 ss << "c†" << sigma;
469 return make_operator(this->cdag_1s(sigma),orbital,PROP::FERMIONIC, ss.str());
470}
471
472template <typename Symmetry_>
473template <typename Dummy>
474typename std::enable_if<Dummy::IS_CHARGE_SU2() and !Dummy::IS_SPIN_SU2(),SiteOperatorQ<Symmetry_, Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
475c (SPIN_INDEX sigma, SUB_LATTICE G, std::size_t orbital) const
476{
477 stringstream ss;
478 ss << "c" << sigma << G;
479 return make_operator(this->c_1s(sigma,G),orbital,PROP::FERMIONIC, ss.str());
480}
481
482template <typename Symmetry_>
483template <typename Dummy>
484typename std::enable_if<Dummy::IS_CHARGE_SU2() and !Dummy::IS_SPIN_SU2(),SiteOperatorQ<Symmetry_, Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
485cdag (SPIN_INDEX sigma, SUB_LATTICE G, std::size_t orbital) const
486{
487 //return c(sigma,G,orbital).adjoint();
488 stringstream ss;
489 ss << "c†" << sigma << G;
490 return make_operator(this->cdag_1s(sigma,G),orbital,PROP::FERMIONIC, ss.str());
491}
492
493template <typename Symmetry_>
494template <typename Dummy>
495typename std::enable_if<Dummy::IS_CHARGE_SU2() and Dummy::IS_SPIN_SU2(),SiteOperatorQ<Symmetry_, Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
496c (SUB_LATTICE G, std::size_t orbital) const
497{
498 stringstream ss;
499 ss << "c" << G;
500 return make_operator(this->c_1s(G),orbital,PROP::FERMIONIC,ss.str());
501}
502
503template <typename Symmetry_>
504template <typename Dummy>
505typename std::enable_if<Dummy::IS_CHARGE_SU2() and Dummy::IS_SPIN_SU2(),SiteOperatorQ<Symmetry_, Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
506cdag (SUB_LATTICE G, std::size_t orbital) const
507{
508// return c(G,orbital).adjoint();
509 stringstream ss;
510 ss << "c†" << G;
511 return make_operator(this->cdag_1s(G),orbital,PROP::FERMIONIC,ss.str());
512}
513
514template <typename Symmetry_>
515SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > FermionBase<Symmetry_>::
516sign_local (std::size_t orbital) const
517{
518 return make_operator(this->F_1s(), orbital, PROP::FERMIONIC, "F");
519}
520
521template <typename Symmetry_>
522SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > FermionBase<Symmetry_>::
523sign (std::size_t orb1, std::size_t orb2) const
524{
525 OperatorType Oout;
526 if (N_orbitals == 1) {Oout = this->F_1s(); Oout.label()="sign"; return Oout;}
527 else if (N_orbitals == 0) {return Zero_vac;}
528 else
529 {
530 Oout = Id();
531 for (int i=orb1; i<N_orbitals; ++i)
532 {
533 Oout = Oout * (nh(i) - ns(i));
534 }
535 for (int i=0; i<orb2; ++i)
536 {
537 Oout = Oout * (nh(i) - ns(i));
538 }
539 Oout.label() = "sign";
540 return Oout;
541 }
542}
543
544template <typename Symmetry_>
545template <typename Dummy>
546typename std::enable_if<true, SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
547n (std::size_t orbital) const
548{
549 return make_operator(this->n_1s(), orbital, PROP::NON_FERMIONIC,"n");
550}
551
552template <typename Symmetry_>
553template <typename Dummy>
554typename std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
555n (SPIN_INDEX sigma, std::size_t orbital) const
556{
557 stringstream ss;
558 ss << "n" << sigma;
559 return make_operator(this->n_1s(sigma), orbital, PROP::NON_FERMIONIC,ss.str());
560}
561
562template <typename Symmetry_>
563SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > FermionBase<Symmetry_>::
564ns (std::size_t orbital) const
565{
566 return make_operator(this->ns_1s(), orbital, PROP::NON_FERMIONIC, "ns");
567}
568
569template <typename Symmetry_>
570SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > FermionBase<Symmetry_>::
571nh (std::size_t orbital) const
572{
573 return make_operator(this->nh_1s(), orbital, PROP::NON_FERMIONIC, "nh");
574}
575
576template <typename Symmetry_>
577template <typename Dummy>
578typename std::enable_if<!Dummy::IS_CHARGE_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
579d (std::size_t orbital) const
580{
581 return make_operator(this->d_1s(), orbital, PROP::NON_FERMIONIC, "d");
582}
583
584template<typename Symmetry_>
585template <typename Dummy>
586typename std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<complex<double>,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
587Rcomp (SPINOP_LABEL Sa, int orbital) const
588{
589 assert(orbital<N_orbitals);
590 SiteOperatorQ<Symmetry_,Eigen::Matrix<complex<double>,Eigen::Dynamic,Eigen::Dynamic> > Oout;
591 SiteOperatorQ<Symmetry_,Eigen::Matrix<complex<double>,Eigen::Dynamic,Eigen::Dynamic> > Scomp_cmplx = Scomp(Sa,orbital).template cast<complex<double> >();
592 if (Sa==iSY)
593 {
594 Oout = 2.*M_PI*Scomp_cmplx;
595 }
596 else
597 {
598 Oout = 2.*1.i*M_PI*Scomp_cmplx;
599 }
600 Oout.data() = Oout.data().exp(1.);
601
602 cout << "Rcomp=" << Oout << endl << endl;
603 // cout << "Re=" << Mtmp.exp().real() << endl << endl;
604 // cout << "Im=" << Mtmp.exp().imag() << endl << endl;
605 // cout << "Op=" << Op << endl << endl;
606
607 return Oout; //SiteOperator<Symmetry,complex<double>>(Op,getQ(Sa));
608}
609
610template <typename Symmetry_>
611template <typename Dummy>
612typename std::enable_if<Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
613S (std::size_t orbital) const
614{
615 return make_operator(this->S_1s(), orbital, PROP::NON_FERMIONIC,"S");
616}
617
618template <typename Symmetry_>
619template <typename Dummy>
620typename std::enable_if<Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
621Sdag (std::size_t orbital) const
622{
623 return S(orbital).adjoint();
624}
625
626template <typename Symmetry_>
627template <typename Dummy>
628typename std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
629Sz (std::size_t orbital) const
630{
631 return make_operator(this->Sz_1s(), orbital, PROP::NON_FERMIONIC,"Sz");
632}
633
634template <typename Symmetry_>
635template <typename Dummy>
636typename std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<complex<double>,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
637exp_ipiSz (std::size_t orbital) const
638{
639 //return make_operator(this->exp_ipiSz(), orbital, PROP::NON_FERMIONIC,"exp(i*Ï€*Sz)");
640
641 ComplexOperatorType out;
642 if (N_orbitals == 1)
643 {
644 out = this->exp_ipiSz_1s(); out.label() = "exp(i*Ï€*Sz)";
645 return out;
646 }
647 else
648 {
649 throw; // not implemented
650 }
651}
652
653template <typename Symmetry_>
654template <typename Dummy>
655typename std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
656Sp (std::size_t orbital) const
657{
658 return make_operator(this->Sp_1s(), orbital, PROP::NON_FERMIONIC,"Sp");
659}
660
661template <typename Symmetry_>
662template <typename Dummy>
663typename std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
664Sm (std::size_t orbital) const
665{
666 return make_operator(this->Sm_1s(), orbital, PROP::NON_FERMIONIC,"Sm");
667}
668
669template <typename Symmetry_>
670template <typename Dummy>
671typename std::enable_if<Dummy::NO_SPIN_SYM(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
672Sx (std::size_t orbital) const
673{
674 OperatorType out = 0.5 * (Sp(orbital) + Sm(orbital));
675 out.label() = "Sx";
676 return out;
677}
678
679template <typename Symmetry_>
680template <typename Dummy>
681typename std::enable_if<Dummy::NO_SPIN_SYM(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
682iSy (std::size_t orbital) const
683{
684 OperatorType out = 0.5 * (Sp(orbital) - Sm(orbital));
685 out.label() = "iSy";
686 return out;
687}
688
689template <typename Symmetry_>
690template <typename Dummy>
691typename std::enable_if<Dummy::IS_CHARGE_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
692T (std::size_t orbital) const
693{
694 return make_operator(this->T_1s(), orbital, PROP::NON_FERMIONIC,"T");
695}
696
697template <typename Symmetry_>
698template <typename Dummy>
699typename std::enable_if<Dummy::IS_CHARGE_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
700Tdag (std::size_t orbital) const
701{
702 return T(orbital).adjoint();
703}
704
705template <typename Symmetry_>
706template <typename Dummy>
707typename std::enable_if<!Dummy::IS_CHARGE_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
708cc (std::size_t orbital) const
709{
710 return make_operator(this->cc_1s(), orbital, PROP::NON_FERMIONIC, "cc");
711}
712
713template <typename Symmetry_>
714template <typename Dummy>
715typename std::enable_if<!Dummy::IS_CHARGE_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
716cdagcdag (std::size_t orbital) const
717{
718// return cc(orbital).adjoint();
719 return make_operator(this->cdagcdag_1s(), orbital, PROP::NON_FERMIONIC, "c†c†");
720}
721
722//SiteOperatorQ<Sym::S1xS2<Sym::SU2<Sym::SpinSU2>,Sym::U1<Sym::ChargeU1> >,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > FermionBase<Sym::S1xS2<Sym::SU2<Sym::SpinSU2>,Sym::U1<Sym::ChargeU1> > >::
723//Tp (std::size_t orbital) const
724//{
725// return Eta(orbital);
726//}
727
728//SiteOperatorQ<Sym::S1xS2<Sym::SU2<Sym::SpinSU2>,Sym::U1<Sym::ChargeU1> >,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > FermionBase<Sym::S1xS2<Sym::SU2<Sym::SpinSU2>,Sym::U1<Sym::ChargeU1> > >::
729//Tm (std::size_t orbital) const
730//{
731// return Eta(orbital).adjoint();
732//}
733
734template <typename Symmetry_>
735template <typename Dummy>
736typename std::enable_if<!Dummy::IS_CHARGE_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
737Tz (std::size_t orbital) const
738{
739 OperatorType out = 0.5*(n(orbital)-Id());
740 out.label() = "Tz=1/2*(n-Id)";
741 return out;
742}
743
744template <typename Symmetry_>
745template <typename Dummy>
746typename std::enable_if<!Dummy::IS_CHARGE_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
747tz (std::size_t orbital) const
748{
749 OperatorType out = n(orbital)-0.5*Id();
750 out.label() = "tz=n-0.5*Id";
751 return out;
752}
753
754template <typename Symmetry_>
755template <typename Dummy>
756typename std::enable_if<Dummy::NO_CHARGE_SYM(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
757Tx (std::size_t orbital, SUB_LATTICE G) const
758{
759 OperatorType out = 0.5 * (Tp(orbital,G) + Tm(orbital, G));
760 out.label() = "Tx";
761 return out;
762}
763
764template <typename Symmetry_>
765template <typename Dummy>
766typename std::enable_if<Dummy::NO_CHARGE_SYM(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type FermionBase<Symmetry_>::
767iTy (std::size_t orbital, SUB_LATTICE G) const
768{
769 OperatorType out = 0.5 * (Tp(orbital,G) - Tm(orbital,G));
770 out.label() = "iTy";
771 return out;
772}
773
774template <typename Symmetry_>
775SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > FermionBase<Symmetry_>::
776Id (std::size_t orbital) const
777{
778 return make_operator(this->Id_1s(), orbital, PROP::NON_FERMIONIC,"Id");
779}
780
781template<typename Symmetry_>
782template<typename Scalar_>
783SiteOperatorQ<Symmetry_, Eigen::Matrix<Scalar_,Eigen::Dynamic,Eigen::Dynamic> > FermionBase<Symmetry_>::
784HubbardHamiltonian (const Array<Scalar_,Dynamic,1> &Uph, const Array<Scalar_,Dynamic,Dynamic> &t,
785 const Array<Scalar_,Dynamic,Dynamic> &V, const Array<Scalar_,Dynamic,Dynamic> &J) const
786{
787 //lout << "HubbardHamiltonian (const Array<Scalar_,Dynamic,1> &Uph, const Array<Scalar_,Dynamic,Dynamic> &t, const Array<Scalar_,Dynamic,Dynamic> &V, const Array<Scalar_,Dynamic,Dynamic> &J)" << endl;
788 //lout << "J=" << J << endl;
789 SiteOperatorQ<Symmetry_, Eigen::Matrix<Scalar_,Eigen::Dynamic,Eigen::Dynamic> > Oout(Symmetry::qvacuum(),TensorBasis);
790
791 for (int i=0; i<N_orbitals; ++i)
792 for (int j=0; j<N_orbitals; ++j)
793 {
794 // Attention: Should be changed to a parameter given to the function:
795 auto G1 = static_cast<SUB_LATTICE>(static_cast<int>(pow(-1,i)));
796 auto G2 = static_cast<SUB_LATTICE>(static_cast<int>(pow(-1,j)));
797 if (t(i,j) != 0.)
798 {
799 if constexpr (Symmetry_::IS_SPIN_SU2() and Symmetry::IS_CHARGE_SU2())
800 {
801 Oout += -t(i,j) * std::sqrt(2.)*std::sqrt(2.) * OperatorType::prod(cdag(G1,i),c(G2,j),Symmetry::qvacuum()).template cast<Scalar_>();
802 }
803 else if constexpr (Symmetry_::IS_SPIN_SU2() and !Symmetry_::IS_CHARGE_SU2())
804 {
805 Oout += -t(i,j)*std::sqrt(2.) * (OperatorType::prod(cdag(i),c(j),Symmetry::qvacuum())).template cast<Scalar_>();
806 Oout += -t(j,i)*std::sqrt(2.) * (OperatorType::prod(c(i),cdag(j),Symmetry::qvacuum())).template cast<Scalar_>();
807 }
808 else if constexpr (!Symmetry_::IS_SPIN_SU2() and Symmetry_::IS_CHARGE_SU2())
809 {
810 Oout += -t(i,j) * std::sqrt(2.) * (OperatorType::prod(cdag(UP,G1,i),c(UP,G2,j),Symmetry::qvacuum()) + OperatorType::prod(cdag(DN,G1,i),c(DN,G2,j),Symmetry::qvacuum())).template cast<Scalar_>();
811 }
812 else if constexpr (!Symmetry_::IS_SPIN_SU2() and !Symmetry_::IS_CHARGE_SU2())
813 {
814 Oout += -t(i,j)*(cdag(UP,i) * c(UP,j) + cdag(DN,i) * c(DN,j)).template cast<Scalar_>();
815 Oout += -t(j,i)*(cdag(UP,j) * c(UP,i) + cdag(DN,j) * c(DN,i)).template cast<Scalar_>();
816 }
817 else
818 {
819 // static_assert(false, "You use a symmetry combination for which there is no implementation of the hopping part in FermionBase::HubbardHamiltonian()");
820 }
821 }
822 if (V(i,j) != 0.)
823 {
824 if constexpr (Symmetry::IS_CHARGE_SU2())
825 {
826 Oout += V(i,j)*std::sqrt(3.) * (OperatorType::prod(Tdag(i),T(j),Symmetry::qvacuum())).template cast<Scalar_>();
827 }
828 else
829 {
830 Oout += V(i,j) * (OperatorType::prod(Tz(i),Tz(j),Symmetry::qvacuum())).template cast<Scalar_>();
831 Oout += 0.5*V(i,j) * (OperatorType::prod(Tp(i,G1),Tm(j,G2),Symmetry::qvacuum()) +
832 OperatorType::prod(Tm(i,G1),Tp(j,G2),Symmetry::qvacuum())).template cast<Scalar_>();
833 }
834 }
835 if (J(i,j) != 0.)
836 {
837 if constexpr (Symmetry::IS_SPIN_SU2())
838 {
839 lout << "Setting SdagS" << endl;
840 Oout += J(i,j)*std::sqrt(3.) * (OperatorType::prod(Sdag(i),S(j),Symmetry::qvacuum())).template cast<Scalar_>();
841 }
842 else
843 {
844 Oout += J(i,j) * (OperatorType::prod(Sz(i),Sz(j),Symmetry::qvacuum())).template cast<Scalar_>();
845 Oout += 0.5*J(i,j) * (OperatorType::prod(Sp(i),Sm(j),Symmetry::qvacuum()) +
846 OperatorType::prod(Sm(i),Sp(j),Symmetry::qvacuum())).template cast<Scalar_>();
847 }
848 }
849 }
850
851 for (int i=0; i<N_orbitals; ++i)
852 {
853 if (Uph(i) != 0. and Uph(i) != std::numeric_limits<double>::infinity())
854 {
855 Oout += 0.5*Uph(i) * nh(i).template cast<Scalar_>();
856 }
857 }
858 Oout.label() = "Hloc";
859 return Oout;
860}
861
862template<typename Symmetry_>
863template<typename Scalar_, typename Dummy>
864typename std::enable_if<Dummy::ABELIAN,SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> > >::type FermionBase<Symmetry_>::
865HubbardHamiltonian (const Array<Scalar_,Dynamic,1> &U,
866 const Array<Scalar_,Dynamic,1> &Uph,
867 const Array<Scalar_,Dynamic,1> &Eorb,
868 const Array<Scalar_,Dynamic,1> &Bz,
869 const Array<Scalar_,Dynamic,Dynamic> &t,
870 const Array<Scalar_,Dynamic,Dynamic> &V,
871 const Array<Scalar_,Dynamic,Dynamic> &Vz,
872 const Array<Scalar_,Dynamic,Dynamic> &Vxy,
873 const Array<Scalar_,Dynamic,Dynamic> &Jz,
874 const Array<Scalar_,Dynamic,Dynamic> &Jxy,
875 const Array<Scalar_,Dynamic,Dynamic> &C) const
876{
877 //lout << "HubbardHamiltonian (const Array<Scalar_,Dynamic,1> &U, const Array<Scalar_,Dynamic,1> &Uph, const Array<Scalar_,Dynamic,1> &Eorb, const Array<Scalar_,Dynamic,1> &Bz, const Array<Scalar_,Dynamic,Dynamic> &t, const Array<Scalar_,Dynamic,Dynamic> &V, const Array<Scalar_,Dynamic,Dynamic> &Vz, const Array<Scalar_,Dynamic,Dynamic> &Vxy, const Array<Scalar_,Dynamic,Dynamic> &Jz, const Array<Scalar_,Dynamic,Dynamic> &Jxy, const Array<Scalar_,Dynamic,Dynamic> &C)" << endl;
878 auto Oout = HubbardHamiltonian<Scalar_>(Uph, t, ZeroHopping(), ZeroHopping());
879
880 for (int i=0; i<N_orbitals; ++i)
881 for (int j=0; j<N_orbitals; ++j)
882 {
883 auto G1 = static_cast<SUB_LATTICE>(static_cast<int>(pow(-1,i)));
884 auto G2 = static_cast<SUB_LATTICE>(static_cast<int>(pow(-1,j)));
885
886 if (V(i,j) != 0.)
887 {
888 Oout += V(i,j) * (n(i)*n(j)).template cast<Scalar_>();
889 }
890 if (Vz(i,j) != 0.)
891 {
892 Oout += Vz(i,j) * (Tz(i)*Tz(j)).template cast<Scalar_>();
893 }
894 if (Vxy(i,j) != 0.)
895 {
896 Oout += 0.5*Vxy(i,j) * (OperatorType::prod(Tp(i,G1),Tm(j,G2),Symmetry::qvacuum()) +
897 OperatorType::prod(Tm(i,G1),Tp(j,G2),Symmetry::qvacuum())).template cast<Scalar_>();
898 }
899 if (Jz(i,j) != 0.)
900 {
901 Oout += Jz(i,j) * (Sz(i)*Sz(j)).template cast<Scalar_>();
902 }
903 if (Jxy(i,j) != 0.)
904 {
905 Oout += 0.5*Jxy(i,j) * (Sp(i)*Sm(j) + Sm(i)*Sp(j)).template cast<Scalar_>();
906 }
907 }
908
909 for (int i=0; i<N_orbitals; ++i)
910 {
911 if (U(i) != 0. and U(i) != numeric_limits<double>::infinity())
912 {
913 Oout += U(i) * d(i).template cast<Scalar_>();
914 }
915 if (Eorb(i) != 0.)
916 {
917 Oout += Eorb(i) * n(i).template cast<Scalar_>();
918 }
919 if (Bz(i) != 0.)
920 {
921 Oout -= Bz(i) * Sz(i).template cast<Scalar_>();
922 }
923 if (C(i) != 0.)
924 {
925 // convention: cUP*cDN + cdagDN*cdagUP
926 // convention for cc, cdagcdag is opposite, therefore needs one commutation
927// Oout -= C(i) * cc(i).template cast<Scalar_>();
928// Oout -= C(i) * cdagcdag(i).template cast<Scalar_>();
929 Oout += C(i) * (c(UP,i)*c(DN,i)).template cast<Scalar_>();
930 Oout += C(i) * (cdag(DN,i)*cdag(UP,i)).template cast<Scalar_>();
931 }
932 }
933 Oout.label() = "Hloc";
934 return Oout;
935}
936
937template<typename Symmetry_>
938template<typename Scalar_, typename Dummy>
939typename std::enable_if<Dummy::IS_SPIN_SU2() and !Dummy::IS_CHARGE_SU2(),SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> > >::type FermionBase<Symmetry_>::
940HubbardHamiltonian (const Array<Scalar_,Dynamic,1> &U,
941 const Array<Scalar_,Dynamic,1> &Uph,
942 const Array<Scalar_,Dynamic,1> &Eorb,
943 const Array<Scalar_,Dynamic,Dynamic> &t,
944 const Array<Scalar_,Dynamic,Dynamic> &V,
945 const Array<Scalar_,Dynamic,Dynamic> &Vz,
946 const Array<Scalar_,Dynamic,Dynamic> &Vxy,
947 const Array<Scalar_,Dynamic,Dynamic> &J) const
948{
949 //lout << "HubbardHamiltonian (const Array<Scalar_,Dynamic,1> &U, const Array<Scalar_,Dynamic,1> &Uph, const Array<Scalar_,Dynamic,1> &Eorb, const Array<Scalar_,Dynamic,Dynamic> &t, const Array<Scalar_,Dynamic,Dynamic> &V, const Array<Scalar_,Dynamic,Dynamic> &Vz, const Array<Scalar_,Dynamic,Dynamic> &Vxy, const Array<Scalar_,Dynamic,Dynamic> &J)" << endl;
950 //lout << "J=" << J << endl;
951 auto Oout = HubbardHamiltonian<Scalar_>(Uph, t, ZeroHopping(), J);
952
953 for (int i=0; i<N_orbitals; ++i)
954 for (int j=0; j<N_orbitals; ++j)
955 {
956 auto G1 = static_cast<SUB_LATTICE>(static_cast<int>(pow(-1,i)));
957 auto G2 = static_cast<SUB_LATTICE>(static_cast<int>(pow(-1,j)));
958 if (V(i,j) != 0.)
959 {
960 Oout += V(i,j) * (n(i)*n(j)).template cast<Scalar_>();
961 }
962 if (Vz(i,j) != 0.)
963 {
964 Oout += Vz(i,j) * (Tz(i)*Tz(j)).template cast<Scalar_>();
965 }
966 if (Vxy(i,j) != 0.)
967 {
968 Oout += 0.5*Vxy(i,j) * (OperatorType::prod(Tp(i,G1),Tm(j,G2),Symmetry::qvacuum()) +
969 OperatorType::prod(Tm(i,G1),Tp(j,G2),Symmetry::qvacuum())).template cast<Scalar_>();
970 }
971 }
972
973 for (int i=0; i<N_orbitals; ++i)
974 {
975 if (U(i) != 0. and U(i) != numeric_limits<double>::infinity())
976 {
977 Oout += U(i) * d(i).template cast<Scalar_>();
978 }
979 if (Eorb(i) != 0.)
980 {
981 Oout += Eorb(i) * n(i).template cast<Scalar_>();
982 }
983 }
984 Oout.label() = "Hloc";
985 return Oout;
986}
987
988template<typename Symmetry_>
989template<typename Scalar_, typename Dummy>
990typename std::enable_if<Dummy::NO_SPIN_SYM() and Dummy::IS_CHARGE_SU2(),SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> > >::type FermionBase<Symmetry_>::
991HubbardHamiltonian (const Array<Scalar_,Dynamic,1> &Uph,
992 const Array<Scalar_,Dynamic,Dynamic> &t,
993 const Array<Scalar_,Dynamic,Dynamic> &V,
994 const Array<Scalar_,Dynamic,Dynamic> &Jz,
995 const Array<Scalar_,Dynamic,Dynamic> &Jxy,
996 const Array<Scalar_,Dynamic,1> &Bz,
997 const Array<Scalar_,Dynamic,1> &Bx) const
998{
999 //lout << "HubbardHamiltonian (const Array<Scalar_,Dynamic,1> &Uph, const Array<Scalar_,Dynamic,Dynamic> &t, const Array<Scalar_,Dynamic,Dynamic> &V, const Array<Scalar_,Dynamic,Dynamic> &Jz, const Array<Scalar_,Dynamic,Dynamic> &Jxy, const Array<Scalar_,Dynamic,1> &Bz, const Array<Scalar_,Dynamic,1> &Bx)" << endl;
1000 auto Oout = HubbardHamiltonian<Scalar_>(Uph, t, V, ZeroHopping());
1001
1002 for (int i=0; i<N_orbitals; ++i)
1003 for (int j=0; j<N_orbitals; ++j)
1004 {
1005 if (Jz(i,j) != 0.)
1006 {
1007 Oout += Jz(i,j) * (Sz(i)*Sz(j)).template cast<Scalar_>();
1008 }
1009 if (Jxy(i,j) != 0.)
1010 {
1011 Oout += 0.5*Jxy(i,j) * (OperatorType::prod(Sp(i),Sm(j),Symmetry::qvacuum()) +
1012 OperatorType::prod(Sm(i),Sp(j),Symmetry::qvacuum())).template cast<Scalar_>();
1013 }
1014 }
1015
1016 for (int i=0; i<N_orbitals; ++i)
1017 {
1018 if (Bz(i) != 0.)
1019 {
1020 Oout += -1. * Bz(i) * Sz(i).template cast<Scalar_>();
1021 }
1022 if (Bx(i) != 0.)
1023 {
1024 Oout += -1. * Bx(i) * Sx(i).template cast<Scalar_>();
1025 }
1026 }
1027 Oout.label() = "Hloc";
1028 return Oout;
1029}
1030
1031template<typename Symmetry_>
1032template<typename Scalar_, typename Dummy>
1033typename std::enable_if<Dummy::IS_SPIN_U1() and Dummy::IS_CHARGE_SU2(),SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> > >::type FermionBase<Symmetry_>::
1034HubbardHamiltonian (const Array<Scalar_,Dynamic,1> &Uph,
1035 const Array<Scalar_,Dynamic,Dynamic> &t,
1036 const Array<Scalar_,Dynamic,Dynamic> &V,
1037 const Array<Scalar_,Dynamic,Dynamic> &Jz,
1038 const Array<Scalar_,Dynamic,Dynamic> &Jxy,
1039 const Array<Scalar_,Dynamic,1> &Bz) const
1040{
1041 //lout << "HubbardHamiltonian (const Array<Scalar_,Dynamic,1> &Uph, const Array<Scalar_,Dynamic,Dynamic> &t, const Array<Scalar_,Dynamic,Dynamic> &V, const Array<Scalar_,Dynamic,Dynamic> &Jz, const Array<Scalar_,Dynamic,Dynamic> &Jxy, const Array<Scalar_,Dynamic,1> &Bz)" << endl;
1042 auto Oout = HubbardHamiltonian<Scalar_>(Uph, t, V, ZeroHopping());
1043
1044 for (int i=0; i<N_orbitals; ++i)
1045 for (int j=0; j<N_orbitals; ++j)
1046 {
1047 if (Jz(i,j) != 0.)
1048 {
1049 Oout += Jz(i,j) * (Sz(i)*Sz(j)).template cast<Scalar_>();
1050 }
1051 if (Jxy(i,j) != 0.)
1052 {
1053 Oout += 0.5*Jxy(i,j) * (OperatorType::prod(Sp(i),Sm(j),Symmetry::qvacuum()) +
1054 OperatorType::prod(Sm(i),Sp(j),Symmetry::qvacuum())).template cast<Scalar_>();
1055 }
1056 }
1057
1058 for (int i=0; i<N_orbitals; ++i)
1059 {
1060 if (Bz(i) != 0.)
1061 {
1062 Oout += -1. * Bz(i) * Sz(i).template cast<Scalar_>();
1063 }
1064 }
1065 Oout.label() = "Hloc";
1066 return Oout;
1067}
1068
1069template<typename Symmetry_>
1070template<typename Scalar_, typename Dummy>
1071typename std::enable_if<Dummy::NO_SPIN_SYM(),SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> > >::type FermionBase<Symmetry_>::
1072coupling_Bx (const Array<double,Dynamic,1> &Bx) const
1073{
1074 SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> > Mout(Symmetry::qvacuum(), TensorBasis);
1075 for (int i=0; i<N_orbitals; ++i)
1076 {
1077 if (Bx(i) != 0.)
1078 {
1079 Mout -= Bx(i) * Sx(i).template cast<Scalar_>();
1080 }
1081 }
1082 return Mout;
1083}
1084
1085template<typename Symmetry_>
1086template<typename Dummy>
1087typename std::enable_if<Dummy::NO_SPIN_SYM(),SiteOperatorQ<Symmetry_,Eigen::Matrix<complex<double>,-1,-1> >>::type FermionBase<Symmetry_>::
1088coupling_By (const Array<double,Dynamic,1> &By) const
1089{
1090 SiteOperatorQ<Symmetry_,Eigen::Matrix<complex<double>,-1,-1> > Mout(Symmetry::qvacuum(), TensorBasis);
1091 for (int i=0; i<N_orbitals; ++i)
1092 {
1093 if (By(i) != 0.)
1094 {
1095 Mout -= -1i*By(i) * iSy(i).template cast<complex<double> >();
1096 }
1097 }
1098 return Mout;
1099}
1100
1101template<typename Symmetry_>
1102template<typename Scalar_, typename Dummy>
1103typename std::enable_if<Dummy::IS_TRIVIAL,SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> > >::type FermionBase<Symmetry_>::
1104coupling_singleFermion (const Array<double,Dynamic,1> &Fp) const
1105{
1106 SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> > Mout(Symmetry::qvacuum(), TensorBasis);
1107 for (int i=0; i<N_orbitals; ++i)
1108 {
1109 if (Fp(i) != 0.)
1110 {
1111 Mout += Fp(i) * (c(UP,i) + cdag(UP,i) + c(DN,i) + cdag(DN,i)).template cast<Scalar_>();
1112 }
1113 }
1114 return Mout;
1115}
1116
1117template<typename Symmetry_>
1118template<typename Scalar_, typename Dummy>
1119typename std::enable_if<Dummy::NO_SPIN_SYM(),SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> > >::type FermionBase<Symmetry_>::
1120coupling_XYZspin (const Array<double,Dynamic,Dynamic> &Jx, const Array<double,Dynamic,Dynamic> &Jy, const Array<double,Dynamic,Dynamic> &Jz) const
1121{
1122 SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> > Mout(Symmetry::qvacuum(), TensorBasis);
1123 for (int i=0; i<N_orbitals; ++i)
1124 for (int j=0; j<N_orbitals; ++j)
1125 {
1126 if (Jx(i,j) != 0.)
1127 {
1128 Mout += Jx(i,j) * OperatorType::prod(Sx(i), Sx(j), Symmetry::qvacuum()).template cast<Scalar_>();
1129 }
1130 if (Jy(i,j) != 0.)
1131 {
1132 Mout += -Jy(i,j) * OperatorType::prod(iSy(i), iSy(j), Symmetry::qvacuum()).template cast<Scalar_>();
1133 }
1134 if (Jz(i,j) != 0.)
1135 {
1136 Mout += Jz(i,j) * OperatorType::prod(Sz(i), Sz(j), Symmetry::qvacuum()).template cast<Scalar_>();
1137 }
1138 }
1139 return Mout;
1140}
1141
1142#endif
SPIN_INDEX
SPIN_INDEX
Definition: DmrgTypedefs.h:36
DN
@ DN
Definition: DmrgTypedefs.h:38
UP
@ UP
Definition: DmrgTypedefs.h:37
SPINOP_LABEL
SPINOP_LABEL
Definition: DmrgTypedefs.h:60
iSY
@ iSY
Definition: DmrgTypedefs.h:60
SZ
@ SZ
Definition: DmrgTypedefs.h:60
SP
@ SP
Definition: DmrgTypedefs.h:60
SX
@ SX
Definition: DmrgTypedefs.h:60
SM
@ SM
Definition: DmrgTypedefs.h:60
SY
@ SY
Definition: DmrgTypedefs.h:60
SUB_LATTICE
SUB_LATTICE
Definition: DmrgTypedefs.h:130
A
@ A
Definition: DmrgTypedefs.h:130
FermionBase::sign_local
OperatorType sign_local(std::size_t orbital=0) const
Definition: FermionBase.h:516
FermionBase::Sx
std::enable_if< Dummy::NO_SPIN_SYM(), OperatorType >::type Sx(size_t orbital=0) const
FermionBase::qType
Symmetry::qType qType
Definition: FermionBase.h:36
FermionBase::Tp
std::enable_if<!Dummy::IS_CHARGE_SU2(), OperatorType >::type Tp(size_t orbital=0, SUB_LATTICE G=A) const
Definition: FermionBase.h:244
FermionBase::Index
Eigen::Index Index
Definition: FermionBase.h:28
FermionBase::N_states
std::size_t N_states
Definition: FermionBase.h:364
FermionBase::cdag
std::enable_if< Dummy::IS_SPIN_SU2() and!Dummy::IS_CHARGE_SU2(), OperatorType >::type cdag(size_t orbital=0) const
FermionBase::make_operator
OperatorType make_operator(const OperatorType &Op_1s, size_t orbital=0, bool FERMIONIC=false, string label="") const
Definition: FermionBase.h:406
FermionBase::coupling_Bx
std::enable_if< Dummy::NO_SPIN_SYM(), SiteOperatorQ< Symmetry_, Eigen::Matrix< Scalar_,-1,-1 > > >::type coupling_Bx(const Array< double, Dynamic, 1 > &Bx) const
Definition: FermionBase.h:1072
FermionBase::n
std::enable_if< true, OperatorType >::type n(std::size_t orbital=0) const
Definition: FermionBase.h:547
FermionBase::OperatorType
SiteOperatorQ< Symmetry, Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > > OperatorType
Definition: FermionBase.h:34
FermionBase::cdagcdag
std::enable_if<!Dummy::IS_CHARGE_SU2(), OperatorType >::type cdagcdag(size_t orbital=0) const
FermionBase::iTy
std::enable_if< Dummy::NO_CHARGE_SYM(), OperatorType >::type iTy(size_t orbital=0, SUB_LATTICE G=A) const
FermionBase::coupling_XYZspin
std::enable_if< Dummy::NO_SPIN_SYM(), SiteOperatorQ< Symmetry_, Eigen::Matrix< Scalar_,-1,-1 > > >::type coupling_XYZspin(const Array< double, Dynamic, Dynamic > &Jx, const Array< double, Dynamic, Dynamic > &Jy, const Array< double, Dynamic, Dynamic > &Jz) const
Definition: FermionBase.h:1120
FermionBase::Id
OperatorType Id(std::size_t orbital=0) const
Definition: FermionBase.h:776
FermionBase::dim
Index dim() const
Definition: FermionBase.h:47
FermionBase::S
std::enable_if< Dummy::IS_SPIN_SU2(), OperatorType >::type S(size_t orbital=0) const
FermionBase::orbitals
std::size_t orbitals() const
Definition: FermionBase.h:50
FermionBase::nh
OperatorType nh(std::size_t orbital=0) const
Definition: FermionBase.h:571
FermionBase::Sp
std::enable_if<!Dummy::IS_SPIN_SU2(), OperatorType >::type Sp(size_t orbital=0) const
FermionBase::TensorBasis
Qbasis< Symmetry > TensorBasis
Definition: FermionBase.h:366
FermionBase::Tdag
std::enable_if< Dummy::IS_CHARGE_SU2(), OperatorType >::type Tdag(size_t orbital=0) const
FermionBase::ComplexOperatorType
SiteOperatorQ< Symmetry, Eigen::Matrix< complex< Scalar >, Eigen::Dynamic, Eigen::Dynamic > > ComplexOperatorType
Definition: FermionBase.h:35
FermionBase::Rcomp
std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ< Symmetry, Eigen::MatrixXcd > >::type Rcomp(SPINOP_LABEL Sa, int orbital) const
Definition: FermionBase.h:587
FermionBase::ZeroField
ArrayXd ZeroField() const
Definition: FermionBase.h:277
FermionBase::coupling_singleFermion
std::enable_if< Dummy::IS_TRIVIAL, SiteOperatorQ< Symmetry_, Eigen::Matrix< Scalar_,-1,-1 > > >::type coupling_singleFermion(const Array< double, Dynamic, 1 > &Fp) const
Definition: FermionBase.h:1104
FermionBase::Tm
std::enable_if<!Dummy::IS_CHARGE_SU2(), OperatorType >::type Tm(size_t orbital=0, SUB_LATTICE G=A) const
Definition: FermionBase.h:251
FermionBase::ns
OperatorType ns(std::size_t orbital=0) const
Definition: FermionBase.h:564
FermionBase::Sz
std::enable_if<!Dummy::IS_SPIN_SU2(), OperatorType >::type Sz(size_t orbital=0) const
FermionBase::cc
std::enable_if<!Dummy::IS_CHARGE_SU2(), OperatorType >::type cc(size_t orbital=0) const
FermionBase::Zero_vac
OperatorType Zero_vac
Definition: FermionBase.h:369
FermionBase::Tz
std::enable_if<!Dummy::IS_CHARGE_SU2(), OperatorType >::type Tz(size_t orbital=0) const
FermionBase::Id_vac
OperatorType Id_vac
Definition: FermionBase.h:369
FermionBase::d
std::enable_if<!Dummy::IS_CHARGE_SU2(), OperatorType >::type d(std::size_t orbital=0) const
Definition: FermionBase.h:579
FermionBase::sign
OperatorType sign(std::size_t orb1=0, std::size_t orb2=0) const
Definition: FermionBase.h:523
FermionBase::Symmetry
Symmetry_ Symmetry
Definition: FermionBase.h:33
FermionBase::exp_ipiSz
std::enable_if<!Dummy::IS_SPIN_SU2(), ComplexOperatorType >::type exp_ipiSz(size_t orbital=0) const
FermionBase::iSy
std::enable_if< Dummy::NO_SPIN_SYM(), OperatorType >::type iSy(size_t orbital=0) const
FermionBase::HubbardHamiltonian
SiteOperatorQ< Symmetry_, Eigen::Matrix< Scalar_,-1,-1 > > HubbardHamiltonian(const Array< Scalar_, Dynamic, 1 > &Uph, const Array< Scalar_, Dynamic, Dynamic > &t, const Array< Scalar_, Dynamic, Dynamic > &V, const Array< Scalar_, Dynamic, Dynamic > &J) const
FermionBase::ZeroHopping
ArrayXXd ZeroHopping() const
Definition: FermionBase.h:280
FermionBase::Scomp
std::enable_if<!Dummy::IS_SPIN_SU2(), OperatorType >::type Scomp(SPINOP_LABEL Sa, int orbital) const
Definition: FermionBase.h:179
FermionBase::Scalar
double Scalar
Definition: FermionBase.h:29
FermionBase::c
std::enable_if< Dummy::IS_SPIN_SU2() and!Dummy::IS_CHARGE_SU2(), OperatorType >::type c(size_t orbital=0) const
FermionBase::T
std::enable_if< Dummy::IS_CHARGE_SU2(), OperatorType >::type T(size_t orbital=0) const
FermionBase::N_orbitals
std::size_t N_orbitals
Definition: FermionBase.h:363
FermionBase::Tx
std::enable_if< Dummy::NO_CHARGE_SYM(), OperatorType >::type Tx(size_t orbital=0, SUB_LATTICE G=A) const
FermionBase::Sdag
std::enable_if< Dummy::IS_SPIN_SU2(), OperatorType >::type Sdag(size_t orbital=0) const
FermionBase::Sm
std::enable_if<!Dummy::IS_SPIN_SU2(), OperatorType >::type Sm(size_t orbital=0) const
FermionBase::get_basis
Qbasis< Symmetry > get_basis() const
Definition: FermionBase.h:358
FermionBase::tz
std::enable_if<!Dummy::IS_CHARGE_SU2(), OperatorType >::type tz(size_t orbital=0) const
FermionBase::coupling_By
std::enable_if< Dummy::NO_SPIN_SYM(), SiteOperatorQ< Symmetry_, Eigen::Matrix< complex< double >,-1,-1 > > >::type coupling_By(const Array< double, Dynamic, 1 > &By) const
Definition: FermionBase.h:1088
FermionSite::basis_1s
Qbasis< Symmetry > basis_1s() const
Definition: FermionSite.h:48
Qbasis
Definition: Qbasis.h:39
Qbasis::push_back
void push_back(const std::tuple< qType, Eigen::Index, std::vector< std::string > > &state)
Definition: Qbasis.h:214
SiteOperatorQ::label
std::string & label()
Definition: SiteOperatorQ.h:120
SiteOperatorQ::setIdentity
void setIdentity()
Definition: SiteOperatorQ.h:403
kind_dummies.h
PROP::FERMIONIC
const bool FERMIONIC
Definition: DmrgTypedefs.h:496
PROP::NON_FERMIONIC
const bool NON_FERMIONIC
Definition: DmrgTypedefs.h:497
std
STL namespace.
Biped::exp
Biped< Symmetry, MatrixType_ > exp(const expScalar x) const
Definition: Biped.h:1069