Skip to content
Snippets Groups Projects
Commit 936bca0c authored by Alvarez, Gonzalo's avatar Alvarez, Gonzalo
Browse files

HubbardOneBand uses HubbardHelper

parent 3a510515
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/Engine/ModelSelector.h
+ 1
− 0
View file @ 936bca0c
......@@ -90,6 +90,7 @@ public:
~ModelSelector()
{
delete modelPtr_;
modelPtr_ = 0;
}
const ModelBaseType& operator()() const { return *modelPtr_; }
......
src/Models/HubbardOneOrbital/BasisOneSpin.h
+ 0
− 12
View file @ 936bca0c
......@@ -183,17 +183,6 @@ namespace LanczosPlusPlus {
return sum;
}
// SizeType getNbyKet(SizeType ket) const
// {
// SizeType sum = 0;
// WordType ketCopy = ket;
// while(ketCopy) {
// if (ketCopy & 1) sum++;
// ketCopy <<= 1;
// }
// return sum;
// }
void doCombinatorial()
{
/* look-up table for binomial coefficients */
......@@ -217,7 +206,6 @@ namespace LanczosPlusPlus {
bitmask_[i] = bitmask_[i-1]<<1;
}
SizeType size_;
SizeType npart_;
PsimagLite::Vector<WordType>::Type data_;
......
src/Models/HubbardOneOrbital/HubbardOneOrbital.h
+ 14
− 242
View file @ 936bca0c
......@@ -7,20 +7,16 @@
#include "BasisHubbardLanczos.h"
#include "BitManip.h"
#include "TypeToString.h"
#include "SparseRow.h"
#include "ParametersModelHubbard.h"
#include "ProgramGlobals.h"
#include "../../Engine/ModelBase.h"
#include "HubbardHelper.h"
namespace LanczosPlusPlus {
template<typename ComplexOrRealType,typename GeometryType,typename InputType>
class HubbardOneOrbital : public ModelBase<ComplexOrRealType,GeometryType,InputType> {
typedef typename PsimagLite::Real<ComplexOrRealType>::Type RealType;
typedef PsimagLite::Matrix<ComplexOrRealType> MatrixType;
typedef ModelBase<ComplexOrRealType,GeometryType,InputType> BaseType;
enum {TERM_HOPPING=0,TERM_NINJ=1,TERM_SUPER=2};
enum {SO_HOPPING_TERM = 1};
......@@ -29,6 +25,9 @@ class HubbardOneOrbital : public ModelBase<ComplexOrRealType,GeometryType,InputT
public:
typedef typename PsimagLite::Real<ComplexOrRealType>::Type RealType;
typedef PsimagLite::Matrix<ComplexOrRealType> MatrixType;
typedef ModelBase<ComplexOrRealType,GeometryType,InputType> BaseType;
typedef ParametersModelHubbard<RealType,InputType> ParametersModelType;
typedef BasisHubbardLanczos<GeometryType> BasisType;
typedef typename BasisType::PairIntType PairIntType;
......@@ -38,20 +37,21 @@ public:
typedef typename BaseType::SparseMatrixType SparseMatrixType;
typedef typename BaseType::VectorType VectorType;
typedef typename BasisType::LabeledOperatorType LabeledOperatorType;
typedef PsimagLite::SparseRow<SparseMatrixType> SparseRowType;
typedef HubbardHelper<BaseType, BasisType, ParametersModelType> HubbardHelperType;
static int const FERMION_SIGN = BasisType::FERMION_SIGN;
HubbardOneOrbital(SizeType nup,
SizeType ndown,
const ParametersModelType& mp,
InputType& io,
const GeometryType& geometry)
: mp_(mp),
: mp_(io),
geometry_(geometry),
basis_(geometry,nup,ndown),
hoppings_(geometry_.numberOfSites(),geometry_.numberOfSites()),
hasJcoupling_(false),
hasCoulombCoupling_(false)
hasCoulombCoupling_(false),
helper_(geometry, mp_, hoppings_, hasJcoupling_, hasCoulombCoupling_)
{
if (mp_.model == "HubbardOneBandExtended" ||
mp_.model == "SuperHubbardExtended") hasCoulombCoupling_ = true;
......@@ -111,61 +111,19 @@ public:
void setupHamiltonian(SparseMatrixType& matrix,
const BasisBaseType& basis) const
{
SizeType hilbert=basis.size();
typename PsimagLite::Vector<RealType>::Type diag(hilbert);
calcDiagonalElements(diag,basis);
SizeType nsite = geometry_.numberOfSites();
matrix.resize(hilbert,hilbert);
// Calculate off-diagonal elements AND store matrix
SizeType nCounter=0;
for (SizeType ispace=0;ispace<hilbert;ispace++) {
SparseRowType sparseRow;
matrix.setRow(ispace,nCounter);
WordType ket1 = basis(ispace,SPIN_UP);
WordType ket2 = basis(ispace,SPIN_DOWN);
// Save diagonal
sparseRow.add(ispace,diag[ispace]);
for (SizeType i=0;i<nsite;i++) {
setHoppingTerm(sparseRow,ket1,ket2,i,basis);
setJTermOffDiagonal(sparseRow,ket1,ket2,i,basis);
}
nCounter += sparseRow.finalize(matrix);
}
matrix.setRow(hilbert,nCounter);
helper_.setupHamiltonian(matrix, basis);
}
void matrixVectorProduct(VectorType &x,VectorType const &y) const
void matrixVectorProduct(VectorType& x,const VectorType& y) const
{
matrixVectorProduct(x,y,basis_);
matrixVectorProduct(x, y, basis_);
}
void matrixVectorProduct(VectorType &x,
VectorType const &y,
const BasisBaseType& basis) const
{
SizeType hilbert=basis.size();
typename PsimagLite::Vector<RealType>::Type diag(hilbert);
calcDiagonalElements(diag,basis);
for (SizeType ispace=0;ispace<hilbert;ispace++)
x[ispace] += diag[ispace]*y[ispace];
diag.clear();
SizeType nsite = geometry_.numberOfSites();
// Calculate off-diagonal elements AND store matrix
for (SizeType ispace=0;ispace<hilbert;ispace++) {
SparseRowType sparseRow;
WordType ket1 = basis(ispace,SPIN_UP);
WordType ket2 = basis(ispace,SPIN_DOWN);
for (SizeType i=0;i<nsite;i++) {
setHoppingTerm(sparseRow,ket1,ket2,i,basis);
setJTermOffDiagonal(sparseRow,ket1,ket2,i,basis);
}
x[ispace] += sparseRow.finalize(y);
}
helper_.matrixVectorProduct(x, y, basis);
}
bool hasNewParts(std::pair<SizeType,SizeType>& newParts,
......@@ -336,193 +294,6 @@ private:
return true;
}
void calcDiagonalElements(typename PsimagLite::Vector<RealType>::Type& diag,
const BasisBaseType& basis) const
{
const RealType zeroPointFive = 0.5;
SizeType hilbert=basis.size();
SizeType nsite = geometry_.numberOfSites();
SizeType orb = 0;
// Calculate diagonal elements
for (SizeType ispace=0;ispace<hilbert;ispace++) {
WordType ket1 = basis(ispace,SPIN_UP);
WordType ket2 = basis(ispace,SPIN_DOWN);
ComplexOrRealType s=0;
for (SizeType i=0;i<nsite;i++) {
// Hubbard term U0
s += mp_.hubbardU[i] *
basis.isThereAnElectronAt(ket1,ket2,i,SPIN_UP,orb) *
basis.isThereAnElectronAt(ket1,ket2,i,SPIN_DOWN,orb);
// SzSz
for (SizeType j=0;j<nsite;j++) {
ComplexOrRealType value = jCoupling(i,j);
if (PsimagLite::real(value) == 0 && PsimagLite::imag(value) == 0) continue;
s += value*zeroPointFive* // double counting i,j
szTerm(ket1,ket2,i,basis)*
szTerm(ket1,ket2,j,basis);
}
// Coulomb
RealType ne = (basis.getN(ket1,ket2,i,SPIN_UP,orb) +
basis.getN(ket1,ket2,i,SPIN_DOWN,orb));
for (SizeType j=0;j<nsite;j++) {
ComplexOrRealType value = coulombCoupling(i,j);
if (PsimagLite::real(value) == 0 && PsimagLite::imag(value) == 0) continue;
RealType tmp2 = basis.getN(ket1,ket2,j,SPIN_UP,orb) +
basis.getN(ket1,ket2,j,SPIN_DOWN,orb);
s += value * ne * tmp2;
}
// Potential term
RealType tmp = mp_.potentialV[i];
if (mp_.potentialT.size()>0)
tmp += mp_.potentialT[i]*mp_.timeFactor;
if (tmp!=0) s += tmp * ne;
}
assert(fabs(PsimagLite::imag(s))<1e-12);
diag[ispace] = PsimagLite::real(s);
}
}
void setHoppingTerm(SparseRowType& sparseRow,
const WordType& ket1,
const WordType& ket2,
SizeType i,
const BasisBaseType& basis) const
{
WordType s1i=(ket1 & BasisType::bitmask(i));
if (s1i>0) s1i=1;
WordType s2i=(ket2 & BasisType::bitmask(i));
if (s2i>0) s2i=1;
SizeType nsite = geometry_.numberOfSites();
SizeType orb = 0;
// Hopping term
for (SizeType j=0;j<nsite;j++) {
if (j<i) continue;
ComplexOrRealType h = hoppings_(i,j);
if (PsimagLite::real(h) == 0 && PsimagLite::imag(h) == 0) continue;
WordType s1j= (ket1 & BasisType::bitmask(j));
if (s1j>0) s1j=1;
WordType s2j= (ket2 & BasisType::bitmask(j));
if (s2j>0) s2j=1;
if (s1i+s1j==1) {
WordType bra1= ket1 ^(BasisType::bitmask(i)|BasisType::bitmask(j));
SizeType temp = basis.perfectIndex(bra1,ket2);
RealType extraSign = (s1i==1) ? FERMION_SIGN : 1;
RealType tmp2 = basis.doSign(ket1,ket2,i,orb,j,orb,SPIN_UP);
ComplexOrRealType cTemp = h*extraSign*tmp2;
if (s1i == 0) cTemp = PsimagLite::conj(cTemp);
assert(temp<basis.size());
sparseRow.add(temp,cTemp);
}
if (s2i+s2j==1) {
WordType bra2= ket2 ^(BasisType::bitmask(i)|BasisType::bitmask(j));
SizeType temp = basis.perfectIndex(ket1,bra2);
RealType extraSign = (s2i==1) ? FERMION_SIGN : 1;
RealType tmp2 = basis.doSign(ket1,ket2,i,orb,j,orb,SPIN_DOWN);
ComplexOrRealType cTemp = h*extraSign*tmp2;
if (s2j == 0) cTemp = PsimagLite::conj(cTemp);
assert(temp<basis.size());
sparseRow.add(temp,cTemp);
}
}
}
void setJTermOffDiagonal(SparseRowType& sparseRow,
const WordType& ket1,
const WordType& ket2,
SizeType i,
const BasisBaseType& basis) const
{
const RealType zeroPointFive = 0.5;
for (SizeType j=0;j<geometry_.numberOfSites();j++) {
ComplexOrRealType value = jCoupling(i,j)*zeroPointFive;
if (PsimagLite::real(value) == 0 && PsimagLite::imag(value) == 0) continue;
value *= 0.5; // double counting i,j
assert(i!=j);
RealType sign = jTermSign(ket1,ket2,i,j,basis);
setSplusSminus(sparseRow,ket1,ket2,i,j,value*sign,basis);
setSplusSminus(sparseRow,ket1,ket2,j,i,value*sign,basis);
}
}
void setSplusSminus(SparseRowType &sparseRow,
const WordType& ket1,
const WordType& ket2,
SizeType i,
SizeType j,
ComplexOrRealType value,
const BasisBaseType &basis) const
{
if (splusSminusNonZero(ket1,ket2,i,j,basis)==0) return;
assert(i!=j);
WordType bra1 = ket1 ^ (BasisType::bitmask(i)|BasisType::bitmask(j));
WordType bra2 = ket2 ^ (BasisType::bitmask(i)|BasisType::bitmask(j));
SizeType temp = basis.perfectIndex(bra1,bra2);
sparseRow.add(temp,value);
}
SizeType splusSminusNonZero(const WordType& ket1,
const WordType& ket2,
SizeType i,
SizeType j,
const BasisBaseType &basis) const
{
SizeType orb = 0;
if (basis.isThereAnElectronAt(ket1,ket2,j,SPIN_UP,orb)==0) return 0;
if (basis.isThereAnElectronAt(ket1,ket2,i,SPIN_UP,orb)==1) return 0;
if (basis.isThereAnElectronAt(ket1,ket2,i,SPIN_DOWN,orb)==0) return 0;
if (basis.isThereAnElectronAt(ket1,ket2,j,SPIN_DOWN,orb)==1) return 0;
return 1;
}
int jTermSign(const WordType& ket1,
const WordType& ket2,
SizeType i,
SizeType j,
const BasisBaseType& basis) const
{
if (i>j) return jTermSign(ket1,ket2,j,i,basis);
SizeType orb = 0;
int x = basis.doSign(ket1,ket2,i,orb,j,orb,SPIN_UP);
x *= basis.doSign(ket1,ket2,i,orb,j,orb,SPIN_DOWN);
return x;
}
RealType szTerm(const WordType& ket1,
const WordType& ket2,
SizeType i,
const BasisBaseType& basis) const
{
SizeType orb = 0;
RealType sz = basis.isThereAnElectronAt(ket1,ket2,i,SPIN_UP,orb);
sz -= basis.isThereAnElectronAt(ket1,ket2,i,SPIN_DOWN,orb);
return 0.5*sz;
}
ComplexOrRealType jCoupling(SizeType i,SizeType j) const
{
if (!hasJcoupling_) return 0;
return geometry_(i,0,j,0,TERM_SUPER);
}
ComplexOrRealType coulombCoupling(SizeType i,SizeType j) const
{
if (!hasCoulombCoupling_) return 0;
return geometry_(i,0,j,0,TERM_NINJ);
}
const ParametersModelType mp_;
const GeometryType& geometry_;
BasisType basis_;
......@@ -530,6 +301,7 @@ private:
bool hasJcoupling_;
bool hasCoulombCoupling_;
mutable typename PsimagLite::Vector<BasisType*>::Type garbage_;
HubbardHelperType helper_;
}; // class HubbardOneOrbital
} // namespace LanczosPlusPlus
#endif
......
src/Models/HubbardOneOrbital/ParametersModelHubbard.h
+ 4
− 0
View file @ 936bca0c
......@@ -99,6 +99,10 @@ struct ParametersModelHubbard {
} catch (std::exception& e) {}
}
ParametersModelHubbard(const ParametersModelHubbard&) = delete;
ParametersModelHubbard& operator=(const ParametersModelHubbard &) = delete;
PsimagLite::String model;
// Do not include here connection parameters
// those are handled by the Geometry
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment