VMPS++
Loading...
Searching...
No Matches
SpinBase.h
Go to the documentation of this file.
1#ifndef SPINBASE_H_
2#define SPINBASE_H_
3
4#include "DmrgTypedefs.h"
5
6#include "tensors/SiteOperatorQ.h"
7#include "sites/SpinSite.h"
8#include <Eigen/Core>
9#include <unsupported/Eigen/MatrixFunctions>
11
12#include "symmetry/kind_dummies.h"
13#include "DmrgTypedefs.h" // for SPIN_INDEX, SPINOP_LABEL
14#include "tensors/SiteOperator.h"
15#include "DmrgExternal.h" // for posmod
16
17//Note: Don't put a name in this documentation with \class .. because doxygen gets confused with template symbols
24template<typename Symmetry_, size_t order=0>
25class SpinBase : public SpinSite<Symmetry_,order>
26{
27 typedef Eigen::Index Index;
28 typedef double Scalar;
29
30public:
31
32 typedef Symmetry_ Symmetry;
33 typedef SiteOperatorQ<Symmetry,Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> > OperatorType;
34 typedef SiteOperatorQ<Symmetry,Eigen::Matrix<complex<Scalar>,Eigen::Dynamic,Eigen::Dynamic> > ComplexOperatorType;
35 typedef typename Symmetry::qType qType;
36
37 SpinBase(){};
38
43 SpinBase (std::size_t L_input, std::size_t D_input, int mfactor=1);
44
46 inline Index dim() const {return static_cast<Index>(N_states);}
47
49 inline std::size_t orbitals() const {return N_orbitals;}
50
52 inline size_t get_D() const {return this->D;}
53
59 OperatorType sign (std::size_t orb1=0, std::size_t orb2=0) const;
60
65 template<class Dummy = Symmetry>
66 typename std::enable_if<!Dummy::IS_SPIN_SU2(),OperatorType>::type n (std::size_t orbital=0) const;
67
69
73 template<class Dummy = Symmetry>
74 typename std::enable_if<Dummy::IS_SPIN_SU2(),OperatorType>::type S (size_t orbital=0) const;
75
80 template<class Dummy = Symmetry>
81 typename std::enable_if<Dummy::IS_SPIN_SU2(),OperatorType>::type Sdag (size_t orbital=0) const;
82
83 template<class Dummy = Symmetry>
84 typename std::enable_if<Dummy::IS_SPIN_SU2(),OperatorType>::type Q (size_t orbital=0) const;
85
86 template<class Dummy = Symmetry>
87 typename std::enable_if<Dummy::IS_SPIN_SU2(),OperatorType>::type Qdag (size_t orbital=0) const;
88
89 template<class Dummy = Symmetry>
90 typename std::enable_if<!Dummy::IS_SPIN_SU2(),OperatorType>::type Qz (size_t orbital=0) const;
91
92 template<class Dummy = Symmetry>
93 typename std::enable_if<!Dummy::IS_SPIN_SU2(),OperatorType>::type Qp (size_t orbital=0) const;
94
95 template<class Dummy = Symmetry>
96 typename std::enable_if<!Dummy::IS_SPIN_SU2(),OperatorType>::type Qm (size_t orbital=0) const;
97
98 template<class Dummy = Symmetry>
99 typename std::enable_if<!Dummy::IS_SPIN_SU2(),OperatorType>::type Qpz (size_t orbital=0) const;
100
101 template<class Dummy = Symmetry>
102 typename std::enable_if<!Dummy::IS_SPIN_SU2(),OperatorType>::type Qmz (size_t orbital=0) const;
103
104 template<class Dummy = Symmetry>
105 typename std::enable_if<!Dummy::IS_SPIN_SU2(),OperatorType>::type Sz (size_t orbital=0) const;
106
107 template<class Dummy = Symmetry>
108 typename std::enable_if<!Dummy::IS_SPIN_SU2(),OperatorType>::type Sp (size_t orbital=0) const;
109
110 template<class Dummy = Symmetry>
111 typename std::enable_if<!Dummy::IS_SPIN_SU2(),OperatorType>::type Sm (size_t orbital=0) const;
112
113 template<class Dummy = Symmetry>
114 typename std::enable_if<Dummy::NO_SPIN_SYM(),OperatorType>::type Sx (size_t orbital=0) const;
115
116 template<class Dummy = Symmetry>
117 typename std::enable_if<Dummy::NO_SPIN_SYM(),OperatorType>::type iSy (size_t orbital=0) const;
118
119 template<class Dummy = Symmetry>
120 typename std::enable_if<!Dummy::IS_SPIN_SU2(),OperatorType>::type Scomp (SPINOP_LABEL Sa, int orbital=0) const
121 {
122 assert(Sa != SY and Sa != QP and Sa != QM and Sa != QPZ and Sa != QMZ);
123 OperatorType out = Zero();
124 if constexpr (Dummy::NO_SPIN_SYM())
125 {
126 if (Sa==SX) { out = Sx(orbital); }
127 else if (Sa==iSY) { out = iSy(orbital); }
128 else if (Sa==SZ) { out = Sz(orbital); }
129 else if (Sa==SP) { out = Sp(orbital); }
130 else if (Sa==SM) { out = Sm(orbital); }
131 }
132 else
133 {
134 assert(Sa != SX and Sa != iSY);
135 if (Sa==SZ) { out = Sz(orbital); }
136 else if (Sa==SP) { out = Sp(orbital); }
137 else if (Sa==SM) { out = Sm(orbital); }
138 }
139 return out;
140 };
141
142 template<class Dummy = Symmetry>
143 typename std::enable_if<!Dummy::IS_SPIN_SU2(),OperatorType>::type Qcomp (SPINOP_LABEL Sa, int orbital=0) const
144 {
145 assert(Sa != SX and Sa != SY and Sa != iSY and Sa != SZ and Sa != SP and Sa != SM);
146 OperatorType out;
147 if (Sa==QZ) {out = Qz(orbital);}
148 else if (Sa==QP) {out = Qp(orbital);}
149 else if (Sa==QM) {out = Qm(orbital);}
150 else if (Sa==QPZ) {out = Qpz(orbital);}
151 else if (Sa==QMZ) {out = Qmz(orbital);}
152 return out;
153 };
155
156 template<class Dummy = Symmetry>
157 typename std::enable_if<!Dummy::IS_SPIN_SU2(),SiteOperatorQ<Symmetry,Eigen::MatrixXcd> >::type Rcomp (SPINOP_LABEL Sa, int orbital=0) const;
158
159 template<class Dummy = Symmetry>
160 typename std::enable_if<Dummy::IS_SPIN_SU2(),OperatorType>::type bead (STRING STR, size_t orbital=0) const;
161
162 template<class Dummy = Symmetry>
163 typename std::enable_if<!Dummy::IS_SPIN_SU2(),OperatorType>::type bead (STRING STR, size_t orbital=0) const;
164
165 template<class Dummy = Symmetry>
166 typename std::enable_if<!Dummy::IS_SPIN_SU2(),ComplexOperatorType>::type exp_ipiSz (size_t orbital=0) const;
167
169 OperatorType Id (std::size_t orbital=0) const;
170
172 OperatorType Zero (std::size_t orbital=0) const;
173
175 ArrayXd ZeroField() const { return ArrayXd::Zero(N_orbitals); }
176
178 ArrayXXd ZeroHopping() const { return ArrayXXd::Zero(N_orbitals,N_orbitals); }
179
191 template<typename Scalar_>
192 SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> > HeisenbergHamiltonian (const Array<Scalar_,Dynamic,Dynamic> &J, const ArrayXd &Offset) const;
193
194 template<typename Dummy = Symmetry>
195 typename std::enable_if<!Dummy::IS_SPIN_SU2(), OperatorType>::type HeisenbergHamiltonian (const ArrayXXd &Jxy, const ArrayXXd &Jz,
196 const ArrayXd &Bz, const ArrayXd &mu, const ArrayXd &nu,
197 const ArrayXd &Kz) const;
198
199 template<typename Dummy = Symmetry>
200 typename std::enable_if<!Dummy::IS_SPIN_SU2(), OperatorType>::type HeisenbergHamiltonian (const ArrayXXcd &Jxy, const ArrayXXcd &Jz) const;
201
213 template<typename Dummy = Symmetry>
214 typename std::enable_if<Dummy::NO_SPIN_SYM(), OperatorType>::type HeisenbergHamiltonian (const ArrayXXd &Jxy, const ArrayXXd &Jz,
215 const ArrayXd &Bz, const ArrayXd &Bx, const ArrayXd &mu, const ArrayXd &nu,
216 const ArrayXd &Kz, const ArrayXd &Kx,
217 const ArrayXXd &Dy) const;
225 template<typename Dummy = Symmetry>
226 typename std::enable_if<Dummy::NO_SPIN_SYM(),SiteOperatorQ<Symmetry_,Eigen::Matrix<complex<double>,-1,-1> > >::type HeisenbergHamiltonian (const std::array<ArrayXXd,3> &J,
227 const std::array<ArrayXd,3> &B,
228 const std::array<ArrayXd,3> &K,
229 const std::array<ArrayXXd,3> &D) const;
230 template<typename Scalar_, typename Dummy = Symmetry>
231 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;
232
233 template<typename Dummy = Symmetry>
234 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;
235
237 Qbasis<Symmetry> get_basis() const { return TensorBasis; }
238
239private:
240 OperatorType make_operator(const OperatorType &Op_1s, size_t orbital=0, string label="") const;
241 std::size_t N_orbitals;
242 std::size_t N_states;
243
244 Qbasis<Symmetry> TensorBasis; //Final basis for N_orbital sites
245};
246
247template <typename Symmetry_, size_t order>
248SpinBase<Symmetry_,order>::
249SpinBase (std::size_t L_input, std::size_t D_input, int mfactor)
250:SpinSite<Symmetry,order>(D_input,mfactor), N_orbitals(L_input)
251{
252 //create basis for spin 0
253 typename Symmetry::qType Q=Symmetry::qvacuum();
254 Eigen::Index inner_dim = 1;
255 Qbasis<Symmetry_> vacuum;
256 vacuum.push_back(Q, inner_dim);
257
258 // // create operators for zero orbitals
259 // Zero_vac = OperatorType(Symmetry::qvacuum(), vacuum);
260 // Zero_vac.setZero();
261 // Id_vac = OperatorType(Symmetry::qvacuum(), vacuum);
262 // Id_vac.setIdentity();
263
264 // create basis for N_orbitals fermionic sites
265 if (N_orbitals == 1) {TensorBasis = this->basis_1s();}
266 else if (N_orbitals == 0) {TensorBasis = vacuum;}
267 else
268 {
269 TensorBasis = this->basis_1s().combine(this->basis_1s());
270 for(std::size_t o=2; o<N_orbitals; o++)
271 {
272 TensorBasis = TensorBasis.combine(this->basis_1s());
273 }
274 }
275
276 N_states = TensorBasis.size();
277}
278
279template<typename Symmetry_, size_t order>
280SiteOperatorQ<Symmetry_, Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > SpinBase<Symmetry_,order>::
281make_operator (const OperatorType &Op_1s, size_t orbital, string label) const
282{
283 OperatorType out;
284 if (N_orbitals == 1) {out = Op_1s; out.label() = label; return out;}
285 else
286 {
287 OperatorType stringOp = this->Id_1s();;
288 bool TOGGLE=false;
289 if (orbital == 0) {out = OperatorType::outerprod(Op_1s,this->Id_1s(),Op_1s.Q()); TOGGLE=true;}
290 else
291 {
292 if (orbital == 1) {out = OperatorType::outerprod(stringOp,Op_1s,Op_1s.Q()); TOGGLE=true;}
293 else {out = OperatorType::outerprod(stringOp,stringOp,Symmetry_::qvacuum());}
294 }
295 for(std::size_t o=2; o<N_orbitals; o++)
296 {
297 if (orbital == o) {out = OperatorType::outerprod(out,Op_1s,Op_1s.Q()); TOGGLE=true; }
298 else if (TOGGLE==false) {out = OperatorType::outerprod(out,stringOp,Symmetry_::qvacuum());}
299 else if (TOGGLE==true) {out = OperatorType::outerprod(out,this->Id_1s(),Op_1s.Q());}
300 }
301 out.label() = label;
302 return out;
303 }
304}
305
306template<typename Symmetry_, size_t order>
307SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > SpinBase<Symmetry_,order>::
308sign (std::size_t orb1, std::size_t orb2) const
309{
310 OperatorType Oout;
311 if (N_orbitals == 1) {Oout = this->F_1s(); Oout.label()="sign"; return Oout;}
312 else
313 {
314 Oout = Id();
315 for (int i=orb1; i<N_orbitals; ++i)
316 {
317 Oout = Oout * (2.*Sz(i));
318 }
319 for (int i=0; i<orb2; ++i)
320 {
321 Oout = Oout * (2.*Sz(i));
322 }
323 Oout.label() = "sign";
324 return Oout;
325 }
326}
327
328template<typename Symmetry_, size_t order>
329template <typename Dummy>
330typename std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type SpinBase<Symmetry_,order>::
331n (std::size_t orbital) const
332{
333 OperatorType Oout = 0.5*Id() - Sz(orbital);
334 Oout.label() = "n";
335 return Oout;
336}
337
338/*template<typename Symmetry_, size_t order>*/
339/*template <typename Dummy>*/
340/*typename std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type SpinBase<Symmetry_,order>::*/
341/*bead (STRING STR, std::size_t orbital) const*/
342/*{*/
343/* if (STR == STRINGZ)*/
344/* {*/
345/* if (this->D%2==0)*/
346/* {*/
347/* lout << termcolor::red << "Warning: exp(iπSz) is not correct for half-integer S!" << termcolor::reset << endl;*/
348/* }*/
349/* return make_operator(this->exp_i_pi_Sz(), orbital, "exp(iπSz)");*/
350/* }*/
351/* else if (STR == STRINGX)*/
352/* {*/
353/* if (this->D%2==0)*/
354/* {*/
355/* lout << termcolor::red << "Warning: exp(iπSx) is not correct for half-integer S!" << termcolor::reset << endl;*/
356/* }*/
357/* return make_operator(this->exp_i_pi_Sx(), orbital, "exp(iπSx)");*/
358/* }*/
359/* else*/
360/* {*/
361/* if (this->D%2==0)*/
362/* {*/
363/* lout << termcolor::red << "Warning: exp(iπSy) is not correct for half-integer S!" << termcolor::reset << endl;*/
364/* }*/
365/* return make_operator(this->exp_i_pi_Sy(), orbital, "exp(iπSy)");*/
366/* }*/
367/*}*/
368
369/*template<typename Symmetry_, size_t order>*/
370/*template <typename Dummy>*/
371/*typename std::enable_if<Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type SpinBase<Symmetry_,order>::*/
372/*bead (STRING STR, std::size_t orbital) const*/
373/*{*/
374/* lout << termcolor::red << "Warning: returning Id instead of exp(iπSa) with SU(2) symmetry!" << termcolor::reset << endl;*/
375/* return make_operator(this->Id_1s(), orbital,"Id");*/
376/*}*/
377
378template<typename Symmetry_, size_t order>
379template <typename Dummy>
380typename std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<complex<double>,Eigen::Dynamic,Eigen::Dynamic> > >::type SpinBase<Symmetry_,order>::
381exp_ipiSz (std::size_t orbital) const
382{
383 //return make_operator(this->exp_ipiSz(), orbital, PROP::NON_FERMIONIC,"exp(i*Ï€*Sz)");
384 ComplexOperatorType out;
385 if (N_orbitals == 1)
386 {
387 out = this->exp_i_pi_Sz(); out.label() = "exp(i*Ï€*Sz)";
388 return out;
389 }
390 else
391 {
392 throw; // not implemented
393 }
394}
395
396template<typename Symmetry_, size_t order>
397template <typename Dummy>
398typename std::enable_if<Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type SpinBase<Symmetry_,order>::
399S (std::size_t orbital) const
400{
401 return make_operator(this->S_1s(), orbital,"S");
402}
403
404template<typename Symmetry_, size_t order>
405template <typename Dummy>
406typename std::enable_if<Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type SpinBase<Symmetry_,order>::
407Sdag (std::size_t orbital) const
408{
409 return S(orbital).adjoint();
410}
411
412template<typename Symmetry_, size_t order>
413template <typename Dummy>
414typename std::enable_if<Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type SpinBase<Symmetry_,order>::
415Q (std::size_t orbital) const
416{
417 return make_operator(this->Q_1s(), orbital,"Q");
418}
419
420template<typename Symmetry_, size_t order>
421template <typename Dummy>
422typename std::enable_if<Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type SpinBase<Symmetry_,order>::
423Qdag (std::size_t orbital) const
424{
425 return Q(orbital).adjoint();
426}
427
428template<typename Symmetry_, size_t order>
429template <typename Dummy>
430typename std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type SpinBase<Symmetry_,order>::
431Qz (std::size_t orbital) const
432{
433 return make_operator(this->Qz_1s(), orbital,"Qz");
434}
435
436template<typename Symmetry_, size_t order>
437template <typename Dummy>
438typename std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type SpinBase<Symmetry_,order>::
439Qp (std::size_t orbital) const
440{
441 return make_operator(this->Qp_1s(), orbital,"Qp");
442}
443
444template<typename Symmetry_, size_t order>
445template <typename Dummy>
446typename std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type SpinBase<Symmetry_,order>::
447Qm (std::size_t orbital) const
448{
449 return make_operator(this->Qm_1s(), orbital,"Qm");
450}
451
452template<typename Symmetry_, size_t order>
453template <typename Dummy>
454typename std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type SpinBase<Symmetry_,order>::
455Qpz (std::size_t orbital) const
456{
457 return make_operator(this->Qpz_1s(), orbital,"Qpz");
458}
459
460template<typename Symmetry_, size_t order>
461template <typename Dummy>
462typename std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type SpinBase<Symmetry_,order>::
463Qmz (std::size_t orbital) const
464{
465 return make_operator(this->Qmz_1s(), orbital,"Qmz");
466}
467
468template<typename Symmetry_, size_t order>
469template <typename Dummy>
470typename std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type SpinBase<Symmetry_,order>::
471Sz (std::size_t orbital) const
472{
473 return make_operator(this->Sz_1s(), orbital,"Sz");
474}
475
476template<typename Symmetry_, size_t order>
477template <typename Dummy>
478typename std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type SpinBase<Symmetry_,order>::
479Sp (std::size_t orbital) const
480{
481 return make_operator(this->Sp_1s(), orbital,"Sp");
482}
483
484template<typename Symmetry_, size_t order>
485template <typename Dummy>
486typename std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type SpinBase<Symmetry_,order>::
487Sm (std::size_t orbital) const
488{
489 return make_operator(this->Sm_1s(), orbital,"Sm");
490}
491
492template<typename Symmetry_, size_t order>
493template <typename Dummy>
494typename std::enable_if<Dummy::NO_SPIN_SYM(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type SpinBase<Symmetry_,order>::
495Sx (std::size_t orbital) const
496{
497 OperatorType out = 0.5 * (Sp(orbital) + Sm(orbital));
498 out.label() = "Sx";
499 return out;
500}
501
502template<typename Symmetry_, size_t order>
503template <typename Dummy>
504typename std::enable_if<Dummy::NO_SPIN_SYM(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > >::type SpinBase<Symmetry_,order>::
505iSy (std::size_t orbital) const
506{
507 OperatorType out = 0.5 * (Sp(orbital) - Sm(orbital));
508 out.label() = "iSy";
509 return out;
510}
511
512template<typename Symmetry_, size_t order>
513SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > SpinBase<Symmetry_,order>::
514Id (std::size_t orbital) const
515{
516 return make_operator(this->Id_1s(), orbital,"Id");
517}
518
519template<typename Symmetry_, size_t order>
520SiteOperatorQ<Symmetry_,Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > SpinBase<Symmetry_,order>::
521Zero (std::size_t orbital) const
522{
523 return make_operator(this->Zero_1s(), orbital, "Zero");
524}
525
526template<typename Symmetry_, size_t order>
527template<typename Scalar_>
528SiteOperatorQ<Symmetry_, Eigen::Matrix<Scalar_,Eigen::Dynamic,Eigen::Dynamic> > SpinBase<Symmetry_,order>::
529HeisenbergHamiltonian (const Array<Scalar_,Dynamic,Dynamic> &J, const ArrayXd &Offset) const
530{
531 SiteOperatorQ<Symmetry_, Eigen::Matrix<Scalar_,Eigen::Dynamic,Eigen::Dynamic> > Oout(Symmetry::qvacuum(),TensorBasis);
532
533 for (int i=0; i<N_orbitals; ++i)
534 for (int j=0; j<N_orbitals; ++j)
535 {
536 if (J(i,j) != 0.)
537 {
538 if constexpr (Symmetry::IS_SPIN_SU2())
539 {
540 Oout += J(i,j)*std::sqrt(3.) * (OperatorType::prod(Sdag(i),S(j),Symmetry::qvacuum())).template cast<Scalar_>();
541 }
542 else
543 {
544 Oout += J(i,j) * (OperatorType::prod(Sz(i),Sz(j),Symmetry::qvacuum())).template cast<Scalar_>();
545 Oout += 0.5*J(i,j) * (OperatorType::prod(Sp(i),Sm(j),Symmetry::qvacuum()) + OperatorType::prod(Sm(i),Sp(j),Symmetry::qvacuum())).template cast<Scalar_>();
546 }
547 }
548 }
549
550 for (int i=0; i<N_orbitals; ++i)
551 {
552 if (Offset(i) != 0.) {Oout += Offset(i) * Id(i);}
553 }
554
555 Oout.label() = "Hloc";
556 return Oout;
557}
558
559template<typename Symmetry_, size_t order>
560template<typename Dummy>
561typename std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,-1,-1> > >::type SpinBase<Symmetry_,order>::
562HeisenbergHamiltonian (const ArrayXXd &Jxy, const ArrayXXd &Jz, const ArrayXd &Bz, const ArrayXd &mu, const ArrayXd &nu, const ArrayXd &Kz) const
563{
564 assert(Bz.rows() == N_orbitals and Kz.rows() == N_orbitals);
565
566 OperatorType Mout(Symmetry::qvacuum(),TensorBasis);
567
568 for (int i=0; i<N_orbitals; ++i)
569 for (int j=0; j<N_orbitals; ++j)
570 {
571 if (Jxy(i,j) != 0.)
572 {
573 Mout += 0.5*Jxy(i,j) * (Scomp(SP,i) * Scomp(SM,j) + Scomp(SM,i) * Scomp(SP,j));
574 }
575 if (Jz(i,j) != 0.)
576 {
577 Mout += Jz(i,j) * Scomp(SZ,i)*Scomp(SZ,j);
578 }
579 }
580
581 for (int i=0; i<N_orbitals; ++i)
582 {
583 if (Bz(i) != 0.) {Mout -= Bz(i) * Scomp(SZ,i);}
584 }
585 for (int i=0; i<N_orbitals; ++i)
586 {
587 if (mu(i) != 0.) {Mout += mu(i) * (Scomp(SZ,i)-0.5*Id());} // for Kitaev chain: -mu*n = -mu*(1/2-Sz) = mu*(Sz-1/2)
588 }
589 for (int i=0; i<N_orbitals; ++i)
590 {
591 if (nu(i) != 0.) {Mout += nu(i) * (Scomp(SZ,i)+0.5*Id());}
592 }
593 for (int i=0; i<N_orbitals; ++i)
594 {
595 if (Kz(i)!=0.) {Mout += Kz(i) * Scomp(SZ,i) * Scomp(SZ,i);}
596 }
597 Mout.label() = "Hloc";
598 return Mout;
599}
600
601template<typename Symmetry_, size_t order>
602template<typename Dummy>
603typename std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,-1,-1> > >::type SpinBase<Symmetry_,order>::
604HeisenbergHamiltonian (const ArrayXXcd &Jxy, const ArrayXXcd &Jz) const
605{
606 OperatorType Mout(Symmetry::qvacuum(),TensorBasis);
607
608 for (int i=0; i<N_orbitals; ++i)
609 for (int j=0; j<N_orbitals; ++j)
610 {
611 if (abs(Jxy(i,j)) != 0.)
612 {
613 Mout += 0.5 * (Jxy(i,j) * Scomp(SP,i) * Scomp(SM,j) + conj(Jxy(i,j)) * Scomp(SM,i) * Scomp(SP,j));
614 }
615 if (abs(Jz(i,j)) != 0.)
616 {
617 Mout += Jz(i,j) * Scomp(SZ,i) * Scomp(SZ,j);
618 }
619 }
620 Mout.label() = "Hloc";
621 return Mout;
622}
623
624template<typename Symmetry_, size_t order>
625template<typename Dummy>
626typename std::enable_if<Dummy::NO_SPIN_SYM(), SiteOperatorQ<Symmetry_,Eigen::Matrix<double,-1,-1> > >::type SpinBase<Symmetry_,order>::
627HeisenbergHamiltonian (const ArrayXXd &Jxy, const ArrayXXd &Jz,
628 const ArrayXd &Bz, const ArrayXd &Bx,
629 const ArrayXd &mu, const ArrayXd &nu,
630 const ArrayXd &Kz, const ArrayXd &Kx,
631 const ArrayXXd &Dy) const
632{
633 assert(Bz.rows() == N_orbitals and Bx.rows() == N_orbitals and Kz.rows() == N_orbitals and Kx.rows() == N_orbitals);
634
635 OperatorType Mout = HeisenbergHamiltonian(Jxy, Jz, Bz, mu, nu, Kz);
636
637 for (int i=0; i<N_orbitals; ++i)
638 for (int j=0; j<i; ++j)
639 {
640 if (Dy(i,j) != 0.) {Mout += Dy(i,j) * (Scomp(SX,i)*Scomp(SZ,j) - Scomp(SZ,i)*Scomp(SX,j));}
641 }
642 for (int i=0; i<N_orbitals; ++i)
643 {
644 if (Bx(i) != 0.) {Mout -= Bx(i) * Scomp(SX,i);}
645 }
646 for (int i=0; i<N_orbitals; ++i)
647 {
648 if (Kx(i)!=0.) {Mout += Kx(i) * Scomp(SX,i) * Scomp(SX,i);}
649 }
650
651 return Mout;
652}
653
654template<typename Symmetry_, size_t order>
655template<typename Scalar_, typename Dummy>
656typename std::enable_if<Dummy::NO_SPIN_SYM(),SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> > >::type SpinBase<Symmetry_,order>::
657coupling_Bx (const Array<double,Dynamic,1> &Bx) const
658{
659 SiteOperatorQ<Symmetry_,Eigen::Matrix<Scalar_,-1,-1> > Mout(Symmetry::qvacuum(), TensorBasis);
660 for (int i=0; i<N_orbitals; ++i)
661 {
662 if (Bx(i) != 0.)
663 {
664 Mout -= Bx(i) * Sx(i).template cast<Scalar_>();
665 }
666 }
667 return Mout;
668}
669
670template<typename Symmetry_, size_t order>
671template<typename Dummy>
672typename std::enable_if<Dummy::NO_SPIN_SYM(),SiteOperatorQ<Symmetry_,Eigen::Matrix<complex<double>,-1,-1> >>::type SpinBase<Symmetry_,order>::
673coupling_By (const Array<double,Dynamic,1> &By) const
674{
675 SiteOperatorQ<Symmetry_,Eigen::Matrix<complex<double>,-1,-1> > Mout(Symmetry::qvacuum(), TensorBasis);
676 for (int i=0; i<N_orbitals; ++i)
677 {
678 if (By(i) != 0.)
679 {
680 Mout -= -1i*By(i) * iSy(i).template cast<complex<double> >();
681 }
682 }
683 return Mout;
684}
685
686template<typename Symmetry_, size_t order>
687template<typename Dummy>
688typename std::enable_if<Dummy::NO_SPIN_SYM(),SiteOperatorQ<Symmetry_,Eigen::Matrix<complex<double>,-1,-1> > >::type SpinBase<Symmetry_,order>::
689HeisenbergHamiltonian (const std::array<ArrayXXd,3> &J, const std::array<ArrayXd,3> &B, const std::array<ArrayXd,3> &K, const std::array<ArrayXXd,3> &D) const
690{
691 SiteOperatorQ<Symmetry_,Eigen::Matrix<complex<double>,-1,-1> > Mout(Symmetry::qvacuum(), TensorBasis);
692 for (int i=0; i<N_orbitals; ++i)
693 for (int j=0; j<N_orbitals; ++j)
694 {
695 //J
696 if (J[0](i,j) != 0.) {Mout += J[0](i,j) * (Sx(i) * Sx(j)).template cast<complex<double>>();}
697 if (J[1](i,j) != 0.) {Mout += -J[1](i,j) * (iSy(i) * iSy(j)).template cast<complex<double>>();}
698 if (J[2](i,j) != 0.) {Mout += J[2](i,j) * (Sz(i) * Sz(j)).template cast<complex<double>>();}
699
700 //D
701 if (D[0](i,j) != 0.) {Mout += D[0](i,j) * (-1.i) * (iSy(i) * Sz(j) - Sz(i) * iSy(j)).template cast<complex<double> >();}
702 if (D[1](i,j) != 0.) {Mout += D[1](i,j) * (Sx(i) * Sz(j) - Sz(i) * Sx(j)).template cast<complex<double> >();}
703 if (D[2](i,j) != 0.) {Mout += D[2](i,j) * (-1.i) * (Sx(i) * iSy(j) - iSy(i) * Sx(j)).template cast<complex<double> >();}
704 }
705
706 for (int i=0; i<N_orbitals; ++i)
707 {
708 // B
709 if (B[0](i) != 0.) {Mout -= B[0](i) * Sx(i).template cast<complex<double> >();}
710 if (B[1](i) != 0.) {Mout -= B[1](i) * (-1.i) * iSy(i).template cast<complex<double> >();}
711 if (B[2](i) != 0.) {Mout -= B[2](i) * Sz(i).template cast<complex<double> >();}
712
713 // K
714 if (K[0](i) != 0.) {Mout += K[0](i) * ( Sx(i) * Sx(i)).template cast<complex<double> >();}
715 if (K[1](i) != 0.) {Mout -= K[1](i) * (iSy(i) * iSy(i)).template cast<complex<double> >();}
716 if (K[2](i) != 0.) {Mout += K[2](i) * ( Sz(i) * Sz(i)).template cast<complex<double> >();}
717 }
718 Mout.label() = "Hloc";
719 return Mout;
720}
721#endif
SPINOP_LABEL
SPINOP_LABEL
Definition: DmrgTypedefs.h:60
QZ
@ QZ
Definition: DmrgTypedefs.h:60
iSY
@ iSY
Definition: DmrgTypedefs.h:60
QM
@ QM
Definition: DmrgTypedefs.h:60
QP
@ QP
Definition: DmrgTypedefs.h:60
SZ
@ SZ
Definition: DmrgTypedefs.h:60
SP
@ SP
Definition: DmrgTypedefs.h:60
SX
@ SX
Definition: DmrgTypedefs.h:60
QPZ
@ QPZ
Definition: DmrgTypedefs.h:60
SM
@ SM
Definition: DmrgTypedefs.h:60
SY
@ SY
Definition: DmrgTypedefs.h:60
QMZ
@ QMZ
Definition: DmrgTypedefs.h:60
B
@ B
Definition: DmrgTypedefs.h:130
STRING
STRING
Definition: DmrgTypedefs.h:78
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
SpinBase::Sdag
std::enable_if< Dummy::IS_SPIN_SU2(), OperatorType >::type Sdag(size_t orbital=0) const
SpinBase::OperatorType
SiteOperatorQ< Symmetry, Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > > OperatorType
Definition: SpinBase.h:33
SpinBase::Qpz
std::enable_if<!Dummy::IS_SPIN_SU2(), OperatorType >::type Qpz(size_t orbital=0) const
SpinBase::Sm
std::enable_if<!Dummy::IS_SPIN_SU2(), OperatorType >::type Sm(size_t orbital=0) const
SpinBase::HeisenbergHamiltonian
SiteOperatorQ< Symmetry_, Eigen::Matrix< Scalar_,-1,-1 > > HeisenbergHamiltonian(const Array< Scalar_, Dynamic, Dynamic > &J, const ArrayXd &Offset) const
SpinBase::Sp
std::enable_if<!Dummy::IS_SPIN_SU2(), OperatorType >::type Sp(size_t orbital=0) const
SpinBase::SpinBase
SpinBase()
Definition: SpinBase.h:37
SpinBase::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: SpinBase.h:657
SpinBase::orbitals
std::size_t orbitals() const
Definition: SpinBase.h:49
SpinBase::N_states
std::size_t N_states
Definition: SpinBase.h:242
SpinBase::Qp
std::enable_if<!Dummy::IS_SPIN_SU2(), OperatorType >::type Qp(size_t orbital=0) const
SpinBase::Qm
std::enable_if<!Dummy::IS_SPIN_SU2(), OperatorType >::type Qm(size_t orbital=0) const
SpinBase::bead
std::enable_if< Dummy::IS_SPIN_SU2(), OperatorType >::type bead(STRING STR, size_t orbital=0) const
SpinBase::ComplexOperatorType
SiteOperatorQ< Symmetry, Eigen::Matrix< complex< Scalar >, Eigen::Dynamic, Eigen::Dynamic > > ComplexOperatorType
Definition: SpinBase.h:34
SpinBase::Index
Eigen::Index Index
Definition: SpinBase.h:27
SpinBase::Zero
OperatorType Zero(std::size_t orbital=0) const
Definition: SpinBase.h:521
SpinBase::get_basis
Qbasis< Symmetry > get_basis() const
Definition: SpinBase.h:237
SpinBase::ZeroHopping
ArrayXXd ZeroHopping() const
Definition: SpinBase.h:178
SpinBase::get_D
size_t get_D() const
Definition: SpinBase.h:52
SpinBase::bead
std::enable_if<!Dummy::IS_SPIN_SU2(), OperatorType >::type bead(STRING STR, size_t orbital=0) const
SpinBase::Qdag
std::enable_if< Dummy::IS_SPIN_SU2(), OperatorType >::type Qdag(size_t orbital=0) const
SpinBase::Q
std::enable_if< Dummy::IS_SPIN_SU2(), OperatorType >::type Q(size_t orbital=0) const
SpinBase::Id
OperatorType Id(std::size_t orbital=0) const
Definition: SpinBase.h:514
SpinBase::Qcomp
std::enable_if<!Dummy::IS_SPIN_SU2(), OperatorType >::type Qcomp(SPINOP_LABEL Sa, int orbital=0) const
Definition: SpinBase.h:143
SpinBase::sign
OperatorType sign(std::size_t orb1=0, std::size_t orb2=0) const
Definition: SpinBase.h:308
SpinBase::TensorBasis
Qbasis< Symmetry > TensorBasis
Definition: SpinBase.h:244
SpinBase::dim
Index dim() const
Definition: SpinBase.h:46
SpinBase::Symmetry
Symmetry_ Symmetry
Definition: SpinBase.h:32
SpinBase::Sx
std::enable_if< Dummy::NO_SPIN_SYM(), OperatorType >::type Sx(size_t orbital=0) const
SpinBase::S
std::enable_if< Dummy::IS_SPIN_SU2(), OperatorType >::type S(size_t orbital=0) const
SpinBase::exp_ipiSz
std::enable_if<!Dummy::IS_SPIN_SU2(), ComplexOperatorType >::type exp_ipiSz(size_t orbital=0) const
SpinBase::n
std::enable_if<!Dummy::IS_SPIN_SU2(), OperatorType >::type n(std::size_t orbital=0) const
Definition: SpinBase.h:331
SpinBase::iSy
std::enable_if< Dummy::NO_SPIN_SYM(), OperatorType >::type iSy(size_t orbital=0) const
SpinBase::N_orbitals
std::size_t N_orbitals
Definition: SpinBase.h:241
SpinBase::Qmz
std::enable_if<!Dummy::IS_SPIN_SU2(), OperatorType >::type Qmz(size_t orbital=0) const
SpinBase::ZeroField
ArrayXd ZeroField() const
Definition: SpinBase.h:175
SpinBase::Scomp
std::enable_if<!Dummy::IS_SPIN_SU2(), OperatorType >::type Scomp(SPINOP_LABEL Sa, int orbital=0) const
Definition: SpinBase.h:120
SpinBase::Qz
std::enable_if<!Dummy::IS_SPIN_SU2(), OperatorType >::type Qz(size_t orbital=0) const
SpinBase::make_operator
OperatorType make_operator(const OperatorType &Op_1s, size_t orbital=0, string label="") const
Definition: SpinBase.h:281
SpinBase::qType
Symmetry::qType qType
Definition: SpinBase.h:35
SpinBase::Scalar
double Scalar
Definition: SpinBase.h:28
SpinBase::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: SpinBase.h:673
SpinBase::Rcomp
std::enable_if<!Dummy::IS_SPIN_SU2(), SiteOperatorQ< Symmetry, Eigen::MatrixXcd > >::type Rcomp(SPINOP_LABEL Sa, int orbital=0) const
SpinBase::Sz
std::enable_if<!Dummy::IS_SPIN_SU2(), OperatorType >::type Sz(size_t orbital=0) const
SpinSite::D
std::size_t D
Definition: SpinSite.h:54
SpinSite::basis_1s
Qbasis< Symmetry > basis_1s() const
Definition: SpinSite.h:50
SpinSite::mfactor
int mfactor
Definition: SpinSite.h:55
kind_dummies.h