-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathSpinDensityWave.cpp
266 lines (245 loc) · 11 KB
/
SpinDensityWave.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
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
#include "TheBlock.h"
#define j1 couplingConstants[0]
#define j2 couplingConstants[1]
#define jprime couplingConstants[2]
#define hIn couplingConstantsIn[3]
#define splus siteBasisH2[0]
#define sz siteBasisH2[1]
#define sminus siteBasisH2[2]
#define rhoBasissplus rhoBasisH2[0]
#define rhoBasissz rhoBasisH2[1]
#define off0RhoBasissplus off0RhoBasisH2[0]
#define off0RhoBasissz off0RhoBasisH2[1]
#define off1RhoBasissplus off1RhoBasisH2[0]
#define off1RhoBasissz off1RhoBasisH2[1]
#define lS blockParts.l
#define mS blockParts.m
#define lE compBlockParts.l
#define mE compBlockParts.m
using namespace Eigen;
Hamiltonian::Hamiltonian()
{
siteBasisH2.resize(nCouplingOperators);
splus << 0., 1.,
0., 0.; // define one-site operators
sminus << 0., 0.,
1., 0.;
sz << .5, 0.,
0., -.5;
};
void Hamiltonian::setParams(const std::vector<double>& couplingConstantsIn,
int lSysIn)
{
couplingConstants.assign(couplingConstantsIn.begin(),
couplingConstantsIn.end() - 1);
lSys = lSysIn;
h = hIn;
};
void Hamiltonian::calcEffectiveH(const std::vector<double>& intSpins)
{
hFromIntSpins.clear();
hFromIntSpins.reserve(lSys);
hFromIntSpins.push_back(-jprime * intSpins[0]);
// interstitial spins have hard-wall boundary conditions
for(int i = 1; i < lSys; i++)
hFromIntSpins.push_back(-jprime * (intSpins[i - 1] + intSpins[i]));
};
rmMatrixX_t
Hamiltonian::EDBASCoupling(int j, const std::vector<MatrixX_t>& rhoBasisH2)
const
{
MatrixX_t plusMinus = kp(rhoBasissplus, sminus);
return couplingConstants[j - 1] * ( kp(rhoBasissz, sz)
+ (plusMinus + plusMinus.adjoint()) / 2);
};
MatrixD_t Hamiltonian::h1(int sitesFromWestEnd) const
{
return -(hFromIntSpins[sitesFromWestEnd] + h) * sz;
};
Matrix<hamScalar, Dynamic, 1>
Hamiltonian::act(const effectiveHams& blockParts,
const effectiveHams& compBlockParts, rmMatrixX_t x,
bool sweepingEast) const
{
int mSd = mS * d,
mEd = mE * d,
totalDimension = mSd * mEd,
// dimension of state on which Hamiltonian acts
lFreeSiteDistFromWestEnd,
rFreeSiteDistFromWestEnd;
if(sweepingEast)
{
lFreeSiteDistFromWestEnd = lS + 1;
rFreeSiteDistFromWestEnd = lSys - 2 - lE;
}
else
{
lFreeSiteDistFromWestEnd = lSys - 2 - lS;
rFreeSiteDistFromWestEnd = lE + 1;
};
// act with system-block-right-hand-free-site bond:
x.resize(mS, d * mEd);
rmMatrixX_t y1234
= j2 * ( sysBlockRSiteCoupling(blockParts.off0RhoBasissz, sz, x, mS, mE)
+ ( sysBlockRSiteCoupling(blockParts.off0RhoBasissplus, sminus,
x, mS, mE)
+ sysBlockRSiteCoupling(blockParts.off0RhoBasissplus
.adjoint(), splus, x, mS, mE)) / 2);
y1234.resize(totalDimension, 1);
// act with left-hand side of system:
rmMatrixX_t y2341 = (blockParts.hS * x).transpose(), // act with system block
cycledx = x.transpose();
y2341.resize(d, mEd * mS);
y2341.noalias()
+= j1 * ( sysBlockLSiteCoupling(blockParts.off0RhoBasissz, sz,
cycledx, mS, mE)
+ ( sysBlockLSiteCoupling(blockParts.off0RhoBasissplus,
sminus, cycledx, mS, mE)
+ sysBlockLSiteCoupling(blockParts.off0RhoBasissplus
.adjoint(), splus, cycledx, mS, mE))
/ 2)
+ j2 * ( sysBlockLSiteCoupling(blockParts.off1RhoBasissz, sz,
cycledx, mS, mE)
+ ( sysBlockLSiteCoupling(blockParts.off1RhoBasissplus,
sminus, cycledx, mS, mE)
+ sysBlockLSiteCoupling(blockParts.off1RhoBasissplus
.adjoint(),
splus, cycledx, mS, mE)) / 2);
// act with system-block-left-hand-free-site bonds
cycledx.resize(d, mEd * mS);
y2341.noalias() += h1(lFreeSiteDistFromWestEnd) * cycledx;
// act with left-hand free site
y2341.resize(d * mEd, mS);
y2341.transposeInPlace();
y2341.resize(totalDimension, 1);
// act with free sites' bond:
rmMatrixX_t yFSBond(mS, d * mEd);
for(int i = 0; i < mS; i++)
{
rmMatrixX_t RHS = x.row(i);
RHS.resize(d, mEd);
yFSBond.row(i)
= j1 * ( freeSiteCouplingTerm(sz, sz, RHS, mE)
+ ( freeSiteCouplingTerm(splus, sminus, RHS, mE)
+ freeSiteCouplingTerm(sminus, splus, RHS, mE)) / 2);
};
yFSBond.resize(totalDimension, 1);
// act with environment-block-left-hand-free-site bond:
x.resize(mSd, mEd);
cycledx = x.transpose();
cycledx.resize(mE, d * mSd);
rmMatrixX_t y3412
= j2 * ( envBlockLSiteCoupling(compBlockParts.off0RhoBasissz, sz,
cycledx, mS, mE)
+ ( envBlockLSiteCoupling(compBlockParts.off0RhoBasissplus,
sminus, cycledx, mS, mE)
+ envBlockLSiteCoupling(compBlockParts.off0RhoBasissplus
.adjoint(),
splus, cycledx, mS, mE)) / 2);
y3412.resize(mEd, mSd);
y3412.transposeInPlace();
y3412.resize(totalDimension, 1);
// act with right-hand half of system:
x.resize(mSd * mE, d);
cycledx = x.transpose();
rmMatrixX_t y4123 = h1(rFreeSiteDistFromWestEnd) * cycledx;
// act with right-hand free site
y4123.resize(d * mSd, mE);
y4123.noalias()
+= j1 * ( envBlockRSiteCoupling(compBlockParts.off0RhoBasissz, sz,
cycledx, mS, mE)
+ ( envBlockRSiteCoupling(compBlockParts.off0RhoBasissplus,
sminus, cycledx, mS, mE)
+ envBlockRSiteCoupling(compBlockParts.off0RhoBasissplus
.adjoint(),
splus, cycledx, mS, mE)) / 2)
+ j2 * ( envBlockRSiteCoupling(compBlockParts.off1RhoBasissz, sz,
cycledx, mS, mE)
+ ( envBlockRSiteCoupling(compBlockParts.off1RhoBasissplus,
sminus, cycledx, mS, mE)
+ envBlockRSiteCoupling(compBlockParts.off1RhoBasissplus
.adjoint(),
splus, cycledx, mS, mE)) / 2);
// act with environment-block-right-hand-free-site bonds
cycledx.resize(d * mSd, mE);
y4123.noalias() += cycledx * compBlockParts.hS.transpose();
// act with environment block
y4123.resize(d, mSd * mE);
y4123.transposeInPlace();
y4123.resize(totalDimension, 1);
return y1234 + y2341 + yFSBond + y3412 + y4123;
};
rmMatrixX_t Hamiltonian::sysBlockRSiteCoupling(const MatrixX_t& hSBond,
const MatrixD_t& hSiteBond,
const rmMatrixX_t& x, int mSIn,
int mEIn) const
{
rmMatrixX_t hBondx = hSBond * x;
hBondx.resize(mSIn * d * mEIn, d);
return hBondx * hSiteBond.transpose();
};
rmMatrixX_t Hamiltonian::sysBlockLSiteCoupling(const MatrixX_t& hSBond,
const MatrixD_t& hSiteBond,
const rmMatrixX_t& cycledx,
int mSIn, int mEIn) const
{
rmMatrixX_t hBondx = cycledx * hSBond.transpose();
// act with system-block half of bond term
hBondx.resize(d, mEIn * d * mSIn);
return hSiteBond * hBondx; // act with left-hand-free-site half of bond term
};
Matrix<hamScalar, 1, Dynamic>
Hamiltonian::freeSiteCouplingTerm(const MatrixD_t& lFreeSiteBond,
const MatrixD_t& rFreeSiteBond,
const rmMatrixX_t& RHS, int mEIn) const
{
rmMatrixX_t hSiteBondRHS = lFreeSiteBond * RHS;
// act with left-hand-free-site half of bond
hSiteBondRHS.resize(mEIn * d, d);
rmMatrixX_t hSiteBondRHShSiteBond
= hSiteBondRHS * rFreeSiteBond.transpose();
// act with right-hand-free-site half of bond
hSiteBondRHShSiteBond.resize(1, d * mEIn * d);
return hSiteBondRHShSiteBond;
};
rmMatrixX_t Hamiltonian::envBlockLSiteCoupling(const MatrixX_t& hEBond,
const MatrixD_t& hSiteBond,
const rmMatrixX_t& cycledx,
int mSIn, int mEIn) const
{
rmMatrixX_t hBondx = hEBond * cycledx;
hBondx.resize(mEIn * d * mSIn, d);
return hBondx * hSiteBond.transpose();
};
rmMatrixX_t Hamiltonian::envBlockRSiteCoupling(const MatrixX_t& hEBond,
const MatrixD_t& hSiteBond,
const rmMatrixX_t& cycledx,
int mSIn, int mEIn) const
{
rmMatrixX_t hBondx = hSiteBond * cycledx;
// act with right-hand free site half of bond term
hBondx.resize(d * mSIn * d, mEIn);
return hBondx * hEBond.transpose();
// act with environment-block half of bond term
};
rmMatrixX_t Hamiltonian::projectedBASCouplings(TheBlock * block) const
{
return j1 * ( block
-> projectBASCoupling(block -> blockParts.off0RhoBasissz, sz)
+ ( block
-> projectBASCoupling(block -> blockParts.off0RhoBasissplus,
sminus)
+ block
-> projectBASCoupling(block -> blockParts
.off0RhoBasissplus.adjoint(),
splus)) / 2)
+ j2 * ( block
-> projectBASCoupling(block -> blockParts.off1RhoBasissz, sz)
+ ( block
-> projectBASCoupling(block -> blockParts.off1RhoBasissplus,
sminus)
+ block
-> projectBASCoupling(block -> blockParts.off1RhoBasissplus
.adjoint(),
splus)) / 2);
};