63 template<
typename Treal,
typename Tmatrix>
70 const Treal trunc = 0,
83 Treal
homo(Treal tol = 1e-15
88 Treal
lumo(Treal tol = 1e-15
98 void print_data(
int const start,
int const stop)
const;
147 template<
typename Treal,
typename Tmatrix>
150 Fun(
int const*
const p,
int const pl, Treal
const s)
151 :pol(p), pollength(pl), shift(s) {
152 assert(shift <= 1 && shift >= 0);
153 assert(pollength >= 0);
155 Treal
eval(Treal
const x)
const {
157 for (
int ind = 0; ind < pollength; ind++ )
158 y = 2 * pol[ind] * y + (1 - 2 * pol[ind]) * y * y;
173 template<
typename Treal,
typename Tmatrix>
176 const Treal trunc,
const int maxmm)
177 :X(F), D(DM), n(size), nocc(noc), frob_trunc(trunc), maxmul(maxmm),
178 lmin(0), lmax(0), nmul(0), nmul_firstpart(0),
179 idemerror(0), tracediff(0), polys(0) {
184 X.add_identity(-
lmax);
195 template<
typename Treal,
typename Tmatrix>
198 assert(nmul_firstpart == 0);
199 Treal delta, beta, trD2;
204 if (nmul >= maxmul) {
207 "Purification reached maxmul" 208 " without convergence", maxmul);
210 if (tracediff[nmul] > 0) {
215 D = -ONE * X * X + TWO * D;
218 D.frob_thresh(frob_trunc);
219 idemerror[nmul] = Tmatrix::frob_diff(D, X);
221 tracediff[nmul] = D.trace() - nocc;
225 if (idemerror[nmul - 1] < 1 / (Treal)4 &&
228 trD2 = (tracediff[nmul] + nocc -
229 2 * polys[nmul - 1] * (tracediff[nmul - 1] + nocc)) /
230 (1 - 2 * polys[nmul - 1]);
233 while (ind > 0 && polys[ind] == polys[ind - 1]) {
234 delta = delta + frob_trunc;
239 }
while((trD2 + beta * (1 + delta) * n - (1 + delta + beta) *
240 (tracediff[nmul - 1] + nocc)) /
241 ((1 + 2 * delta) * (delta + beta)) < n - nocc - 1 ||
242 (trD2 - delta * (1 - beta) * n - (1 - delta - beta) *
243 (tracediff[nmul - 1] + nocc)) /
244 ((1 + 2 * delta) * (delta + beta)) < nocc - 1);
252 if (tracediff[nmul - 1] > 0) {
255 D = -ONE * X * X + TWO * D;
262 D.frob_thresh(frob_trunc);
263 idemerror[nmul] = Tmatrix::frob_diff(D, X);
265 tracediff[nmul] = D.trace() - nocc;
267 nmul_firstpart = nmul;
271 if (nmul + 1 >= maxmul) {
274 "Purification reached maxmul" 275 " without convergence", maxmul);
277 if (tracediff[nmul] > 0) {
279 idemerror[nmul] = Tmatrix::frob_diff(D, X);
283 tracediff[nmul] = D.trace() - nocc;
285 D = -ONE * X * X + TWO * D;
286 idemerror[nmul] = Tmatrix::frob_diff(D, X);
289 tracediff[nmul] = D.trace() - nocc;
293 X = -ONE * D * D + TWO * X;
294 idemerror[nmul] = Tmatrix::frob_diff(D, X);
297 tracediff[nmul] = X.trace() - nocc;
300 idemerror[nmul] = Tmatrix::frob_diff(D, X);
303 tracediff[nmul] = D.trace() - nocc;
305 D.frob_thresh(frob_trunc);
307 }
while (idemerror[nmul - 1] > frob_trunc);
314 template<
typename Treal,
typename Tmatrix>
316 Fun const fermifun(polys, nmul, 0.5);
317 Treal chempot =
bisection(fermifun, (Treal)0, (Treal)1, tol);
318 return (lmin - lmax) * chempot + lmax;
321 template<
typename Treal,
typename Tmatrix>
325 for (
int mul = nmul_firstpart; mul < nmul; mul++) {
326 if (idemerror[mul] < 1.0 / 4) {
328 tmp =
bisection(homofun, (Treal)0, (Treal)1, tol);
333 homo = tmp > homo ? tmp : homo;
336 return (lmin - lmax) * homo + lmax;
339 template<
typename Treal,
typename Tmatrix>
343 for (
int mul = nmul_firstpart; mul < nmul; mul++) {
344 if (idemerror[mul] < 1.0 / 4) {
346 tmp =
bisection(lumofun, (Treal)0, (Treal)1, tol);
351 lumo = tmp < lumo ? tmp : lumo;
354 return (lmin - lmax) * lumo + lmax;
357 template<
typename Treal,
typename Tmatrix>
359 for (
int ind = start; ind < stop; ind ++) {
360 std::cout <<
"Iteration: " << ind
361 <<
" Idempotency error: " << idemerror[ind]
362 <<
" Tracediff: " << tracediff[ind]
363 <<
" Poly: " << polys[ind]
Treal const shift
Definition: TC2.h:169
const Treal frob_trunc
Threshold for the truncation.
Definition: TC2.h:112
Treal bisection(Tfun const &fun, Treal min, Treal max, Treal const tol)
Bisection algorithm for root finding.
Definition: bisection.h:70
Treal eval(Treal const x) const
Definition: TC2.h:155
Treal fermi_level(Treal tol=1e-15) const
Returns the Fermi level.
Definition: TC2.h:315
Tmatrix & D
Density matrix after purification.
Definition: TC2.h:109
int nmul
Number of used matrix multiplications.
Definition: TC2.h:116
Treal * idemerror
Upper bound of euclidean norm ||D-D^2||_2 before each step.
Definition: TC2.h:120
Treal * tracediff
The difference between the trace of the matrix and the number of occupied orbitals before each step...
Definition: TC2.h:129
Treal lmax
Upper bound for eigenvalue spectrum.
Definition: TC2.h:115
Definition: allocate.cc:39
Treal lmin
Lower bound for eigenvalue spectrum.
Definition: TC2.h:114
Fun(int const *const p, int const pl, Treal const s)
Definition: TC2.h:150
Treal lumo(Treal tol=1e-15) const
Returns lower bound of the LUMO eigenvalue.
Definition: TC2.h:340
int n_multiplies() const
Returns the number of used matrix matrix multiplications.
Definition: TC2.h:94
void print_data(int const start, int const stop) const
Definition: TC2.h:358
const int maxmul
Number of tolerated matrix multiplications.
Definition: TC2.h:113
int nmul_firstpart
Number of used matrix multiplications in the first part of the purification.
Definition: TC2.h:117
const int nocc
Number of occupied orbitals.
Definition: TC2.h:111
int const pollength
Definition: TC2.h:168
TC2(Tmatrix &F, Tmatrix &DM, const int size, const int noc, const Treal trunc=0, const int maxmm=100)
Constructor Initializes everything.
Definition: TC2.h:174
int const *const pol
Definition: TC2.h:167
const int n
System size.
Definition: TC2.h:110
void purify()
Runs purification.
Definition: TC2.h:196
int * polys
Choices of polynomials 0 for x^2 and 1 for 2x-x^2 Length: nmul.
Definition: TC2.h:133
Treal homo(Treal tol=1e-15) const
Returns upper bound of the HOMO eigenvalue.
Definition: TC2.h:322
Help class for bisection root finding calls.
Definition: TC2.h:148
virtual ~TC2()
Destructor.
Definition: TC2.h:99
Treal template_blas_sqrt(Treal x)
Tmatrix & X
Fock / Kohn-Sham matrix at initialization.
Definition: TC2.h:105
Trace correcting purification.
Definition: TC2.h:64