-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathTheBlock.cpp
151 lines (141 loc) · 7.13 KB
/
TheBlock.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
#include "FreeFunctions.h"
using namespace Eigen;
TheBlock::TheBlock(int m, const MatrixX_t& hS, const std::vector<int>& qNumList,
const matPair newRhoBasisH2s, int l)
: m(m), hS(hS), qNumList(qNumList), off0RhoBasisH2(newRhoBasisH2s.first),
off1RhoBasisH2(newRhoBasisH2s.second), l(l) {};
TheBlock::TheBlock(const Hamiltonian& ham)
: m(d), hS(MatrixD_t::Zero()), qNumList(ham.oneSiteQNums), l(0)
{
off0RhoBasisH2.assign(ham.siteBasisH2.begin(),
ham.siteBasisH2.begin() + nIndepCouplingOperators);
};
TheBlock TheBlock::nextBlock(const stepData& data, rmMatrixX_t& psiGround,
double& cumulativeTruncationError)
{
std::vector<int> hSprimeQNumList;
MatrixX_t hSprime = createHprime(this, data.ham, hSprimeQNumList);
// expanded system block
int md = m * d;
if(data.exactDiag)
return TheBlock(md, hSprime, hSprimeQNumList,
createNewRhoBasisH2s(data.ham.siteBasisH2, true), l + 1);
// if near edge of system, no truncation necessary so skip DMRG algorithm
HamSolver hSuperSolver = createHSuperSolver(data, hSprime, hSprimeQNumList,
psiGround);
// calculate superblock ground state
int compm = data.compBlock -> m;
hSuperSolver.lowestEvec.resize(md, compm * d);
DMSolver rhoSolver(hSuperSolver, data.mMax);
// find density matrix eigenstates
cumulativeTruncationError += rhoSolver.truncationError;
primeToRhoBasis = rhoSolver.highestEvecs; // construct change-of-basis matrix
int nextBlockm = primeToRhoBasis.cols();
if(!data.infiniteStage) // modify psiGround to predict the next ground state
{
psiGround = hSuperSolver.lowestEvec;
for(int sPrimeIndex = 0; sPrimeIndex < md; sPrimeIndex++)
// transpose the environment block and the right-hand free site
{
rmMatrixX_t ePrime = psiGround.row(sPrimeIndex);
ePrime.resize(compm, d);
ePrime.transposeInPlace();
ePrime.resize(1, d * compm);
psiGround.row(sPrimeIndex) = ePrime;
};
psiGround = primeToRhoBasis.adjoint() * psiGround;
// change the expanded system block basis
psiGround.resize(nextBlockm * d, compm);
psiGround *= data.beforeCompBlock -> primeToRhoBasis.transpose();
// change the environment block basis
psiGround.resize(nextBlockm * d * data.beforeCompBlock -> m * d, 1);
};
return TheBlock(nextBlockm, changeBasis(hSprime), rhoSolver.highestEvecQNums,
createNewRhoBasisH2s(data.ham.siteBasisH2, false), l + 1);
// save expanded-block operators in new basis
};
MatrixX_t TheBlock::createHprime(const TheBlock* block, const Hamiltonian& ham,
std::vector<int>& hPrimeQNumList) const
{
MatrixX_t hPrime = kp(block -> hS, Id_d)
+ ham.blockAdjacentSiteJoin(1, block -> off0RhoBasisH2);
if(l != 0)
hPrime += ham.blockAdjacentSiteJoin(2, block -> off1RhoBasisH2);
// hSprime += (l == 0 ? -ham.blockAdjacentSiteJoin(1, off0RhoBasisH2) / 2
// : ham.blockAdjacentSiteJoin(2, off1RhoBasisH2));
hPrimeQNumList = vectorProductSum(block -> qNumList, ham.oneSiteQNums);
// add in quantum numbers of new site
return hPrime;
};
matPair TheBlock::createNewRhoBasisH2s(const vecMatD_t& siteBasisH2,
bool exactDiag) const
{
std::vector<MatrixX_t> newOff0RhoBasisH2,
newOff1RhoBasisH2;
newOff0RhoBasisH2.reserve(nIndepCouplingOperators);
newOff1RhoBasisH2.reserve(nIndepCouplingOperators);
for(int i = 0; i < nIndepCouplingOperators; i++)
{
newOff0RhoBasisH2.push_back(exactDiag ?
kp(Id(m), siteBasisH2[i]) :
changeBasis(kp(Id(m), siteBasisH2[i])));
newOff1RhoBasisH2.push_back(exactDiag ?
kp(off0RhoBasisH2[i], Id_d) :
changeBasis(kp(off0RhoBasisH2[i], Id_d)));
};
return std::make_pair(newOff0RhoBasisH2, newOff1RhoBasisH2);
};
HamSolver TheBlock::createHSuperSolver(const stepData data,
const MatrixX_t& hSprime,
const std::vector<int>& hSprimeQNumList,
rmMatrixX_t& psiGround) const
{
MatrixX_t hEprime; // expanded environment block
std::vector<int> hEprimeQNumList;
int scaledTargetQNum;
if(data.infiniteStage)
{
hEprime = hSprime;
hEprimeQNumList = hSprimeQNumList;
scaledTargetQNum = data.ham.targetQNum * (l + 2) / data.ham.lSys * 2;
// int automatically rounds down
// during iDMRG stage, targets correct quantum number
// per unit site by scaling to fit current system size
// - note: this will change if d != 2
}
else
{
hEprime = createHprime(data.compBlock, data.ham, hEprimeQNumList);
scaledTargetQNum = data.ham.targetQNum;
};
int compm = data.compBlock -> m;
MatrixX_t hSuper = kp(hSprime, Id(compm * d))
+ data.ham.lBlockrSiteJoin(off0RhoBasisH2, compm)
+ data.ham.siteSiteJoin(m, compm)
+ data.ham.lSiterBlockJoin(m, data.compBlock
-> off0RhoBasisH2)
+ kp(Id(m * d), hEprime); // superblock
return HamSolver(hSuper, hSprimeQNumList, hEprimeQNumList, scaledTargetQNum,
psiGround, data.lancTolerance);
};
FinalSuperblock TheBlock::createHSuperFinal(const stepData& data,
rmMatrixX_t& psiGround, int skips)
const
{
std::vector<int> hSprimeQNumList;
MatrixX_t hSprime = createHprime(this, data.ham, hSprimeQNumList);
// expanded system block
HamSolver hSuperSolver = createHSuperSolver(data, hSprime, hSprimeQNumList,
psiGround);
// calculate superblock ground state
return FinalSuperblock(hSuperSolver, data.ham.lSys, m, data.compBlock -> m,
skips);
};
obsMatrixX_t TheBlock::obsChangeBasis(const obsMatrixX_t& mat) const
{
return primeToRhoBasis.adjoint() * mat * primeToRhoBasis;
};
MatrixX_t TheBlock::changeBasis(const MatrixX_t& mat) const
{
return primeToRhoBasis.adjoint() * mat * primeToRhoBasis;
};