63 template<
typename Treal,
typename Tmatrix>
147 template<
typename Treal,
typename Tmatrix>
151 :pol(p), pollength(
pl), shift(s) {
153 assert(pollength >= 0);
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>
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);
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]) {
240 (tracediff[nmul - 1] + nocc)) /
241 ((1 + 2 *
delta) * (
delta + beta)) < n - nocc - 1 ||
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;
357template<
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]
Help class for bisection root finding calls.
Definition TC2.h:148
Treal const shift
Definition TC2.h:169
int const *const pol
Definition TC2.h:167
Fun(int const *const p, int const pl, Treal const s)
Definition TC2.h:150
Treal eval(Treal const x) const
Definition TC2.h:155
int const pollength
Definition TC2.h:168
Trace correcting purification.
Definition TC2.h:64
Treal lmax
Upper bound for eigenvalue spectrum.
Definition TC2.h:115
const Treal frob_trunc
Threshold for the truncation.
Definition TC2.h:112
const int nocc
Number of occupied orbitals.
Definition TC2.h:111
Tmatrix & X
Fock / Kohn-Sham matrix at initialization.
Definition TC2.h:105
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
const int maxmul
Number of tolerated matrix multiplications.
Definition TC2.h:113
int n_multiplies() const
Returns the number of used matrix matrix multiplications.
Definition TC2.h:94
int nmul_firstpart
Number of used matrix multiplications in the first part of the purification.
Definition TC2.h:117
virtual ~TC2()
Destructor.
Definition TC2.h:99
Treal * idemerror
Upper bound of euclidean norm ||D-D^2||_2 before each step.
Definition TC2.h:120
void print_data(int const start, int const stop) const
Definition TC2.h:358
Treal fermi_level(Treal tol=1e-15) const
Returns the Fermi level.
Definition TC2.h:315
Treal * tracediff
The difference between the trace of the matrix and the number of occupied orbitals before each step.
Definition TC2.h:129
Treal homo(Treal tol=1e-15) const
Returns upper bound of the HOMO eigenvalue.
Definition TC2.h:322
const int n
System size.
Definition TC2.h:110
int * polys
Choices of polynomials 0 for x^2 and 1 for 2x-x^2 Length: nmul.
Definition TC2.h:133
Tmatrix & D
Density matrix after purification.
Definition TC2.h:109
int nmul
Number of used matrix multiplications.
Definition TC2.h:116
Treal lumo(Treal tol=1e-15) const
Returns lower bound of the LUMO eigenvalue.
Definition TC2.h:340
void purify()
Runs purification.
Definition TC2.h:196
Treal lmin
Lower bound for eigenvalue spectrum.
Definition TC2.h:114
Definition allocate.cc:39
Treal bisection(Tfun const &fun, Treal min, Treal max, Treal const tol)
Bisection algorithm for root finding.
Definition bisection.h:70
static Treal getMachineEpsilon()
Definition matInclude.h:147
Treal template_blas_sqrt(Treal x)