61 template<
typename Treal,
typename Tmatrix>
68 const Treal trunc = 0,
81 Treal
homo(Treal tol = 1e-15
86 Treal
lumo(Treal tol = 1e-15
96 void print_data(
int const start,
int const stop)
const;
145 template<
typename Treal,
typename Tmatrix>
148 Fun(
int const*
const p,
int const pl, Treal
const s)
149 :pol(p), pollength(pl), shift(s) {
150 assert(shift <= 1 && shift >= 0);
151 assert(pollength >= 0);
153 Treal
eval(Treal
const x)
const {
155 for (
int ind = 0; ind < pollength; ind++ )
156 y = 2 * pol[ind] * y + (1 - 2 * pol[ind]) * y * y;
171 template<
typename Treal,
typename Tmatrix>
174 const Treal trunc,
const int maxmm)
175 :X(F), D(DM), n(size), nocc(noc), frob_trunc(trunc), maxmul(maxmm),
176 lmin(0), lmax(0), nmul(0), nmul_firstpart(0),
177 idemerror(0), tracediff(0), polys(0) {
182 X.add_identity(-
lmax);
193 template<
typename Treal,
typename Tmatrix>
196 assert(nmul_firstpart == 0);
197 Treal delta, beta, trD2;
202 if (nmul >= maxmul) {
205 "Purification reached maxmul"
206 " without convergence", maxmul);
208 if (tracediff[nmul] > 0) {
213 D = -ONE * X * X + TWO * D;
216 D.frob_thresh(frob_trunc);
217 idemerror[nmul] = Tmatrix::frob_diff(D, X);
219 tracediff[nmul] = D.trace() - nocc;
223 if (idemerror[nmul - 1] < 1 / (Treal)4 &&
226 trD2 = (tracediff[nmul] + nocc -
227 2 * polys[nmul - 1] * (tracediff[nmul - 1] + nocc)) /
228 (1 - 2 * polys[nmul - 1]);
231 while (ind > 0 && polys[ind] == polys[ind - 1]) {
232 delta = delta + frob_trunc;
237 }
while((trD2 + beta * (1 + delta) * n - (1 + delta + beta) *
238 (tracediff[nmul - 1] + nocc)) /
239 ((1 + 2 * delta) * (delta + beta)) < n - nocc - 1 ||
240 (trD2 - delta * (1 - beta) * n - (1 - delta - beta) *
241 (tracediff[nmul - 1] + nocc)) /
242 ((1 + 2 * delta) * (delta + beta)) < nocc - 1);
250 if (tracediff[nmul - 1] > 0) {
253 D = -ONE * X * X + TWO * D;
260 D.frob_thresh(frob_trunc);
261 idemerror[nmul] = Tmatrix::frob_diff(D, X);
263 tracediff[nmul] = D.trace() - nocc;
265 nmul_firstpart = nmul;
269 if (nmul + 1 >= maxmul) {
272 "Purification reached maxmul"
273 " without convergence", maxmul);
275 if (tracediff[nmul] > 0) {
277 idemerror[nmul] = Tmatrix::frob_diff(D, X);
281 tracediff[nmul] = D.trace() - nocc;
283 D = -ONE * X * X + TWO * D;
284 idemerror[nmul] = Tmatrix::frob_diff(D, X);
287 tracediff[nmul] = D.trace() - nocc;
291 X = -ONE * D * D + TWO * X;
292 idemerror[nmul] = Tmatrix::frob_diff(D, X);
295 tracediff[nmul] = X.trace() - nocc;
298 idemerror[nmul] = Tmatrix::frob_diff(D, X);
301 tracediff[nmul] = D.trace() - nocc;
303 D.frob_thresh(frob_trunc);
305 }
while (idemerror[nmul - 1] > frob_trunc);
312 template<
typename Treal,
typename Tmatrix>
314 Fun const fermifun(polys, nmul, 0.5);
315 Treal chempot =
bisection(fermifun, (Treal)0, (Treal)1, tol);
316 return (lmin - lmax) * chempot + lmax;
319 template<
typename Treal,
typename Tmatrix>
323 for (
int mul = nmul_firstpart; mul < nmul; mul++) {
324 if (idemerror[mul] < 1.0 / 4) {
326 tmp =
bisection(homofun, (Treal)0, (Treal)1, tol);
331 homo = tmp > homo ? tmp : homo;
334 return (lmin - lmax) * homo + lmax;
337 template<
typename Treal,
typename Tmatrix>
341 for (
int mul = nmul_firstpart; mul < nmul; mul++) {
342 if (idemerror[mul] < 1.0 / 4) {
344 tmp =
bisection(lumofun, (Treal)0, (Treal)1, tol);
349 lumo = tmp < lumo ? tmp : lumo;
352 return (lmin - lmax) * lumo + lmax;
355 template<
typename Treal,
typename Tmatrix>
357 for (
int ind = start; ind < stop; ind ++) {
358 std::cout <<
"Iteration: " << ind
359 <<
" Idempotency error: " << idemerror[ind]
360 <<
" Tracediff: " << tracediff[ind]
361 <<
" Poly: " << polys[ind]
Treal const shift
Definition: TC2.h:167
const Treal frob_trunc
Threshold for the truncation.
Definition: TC2.h:110
Treal bisection(Tfun const &fun, Treal min, Treal max, Treal const tol)
Bisection algorithm for root finding.
Definition: bisection.h:68
Tmatrix & D
Density matrix after purification.
Definition: TC2.h:107
int nmul
Number of used matrix multiplications.
Definition: TC2.h:114
Treal * idemerror
Upper bound of euclidean norm ||D-D^2||_2 before each step.
Definition: TC2.h:118
Treal * tracediff
The difference between the trace of the matrix and the number of occupied orbitals before each step...
Definition: TC2.h:127
Treal homo(Treal tol=1e-15) const
Returns upper bound of the HOMO eigenvalue.
Definition: TC2.h:320
Treal lmax
Upper bound for eigenvalue spectrum.
Definition: TC2.h:113
Treal lumo(Treal tol=1e-15) const
Returns lower bound of the LUMO eigenvalue.
Definition: TC2.h:338
Treal lmin
Lower bound for eigenvalue spectrum.
Definition: TC2.h:112
int n_multiplies() const
Returns the number of used matrix matrix multiplications.
Definition: TC2.h:92
Fun(int const *const p, int const pl, Treal const s)
Definition: TC2.h:148
Treal fermi_level(Treal tol=1e-15) const
Returns the Fermi level.
Definition: TC2.h:313
const int maxmul
Number of tolerated matrix multiplications.
Definition: TC2.h:111
Treal eval(Treal const x) const
Definition: TC2.h:153
int nmul_firstpart
Number of used matrix multiplications in the first part of the purification.
Definition: TC2.h:115
const int nocc
Number of occupied orbitals.
Definition: TC2.h:109
int const pollength
Definition: TC2.h:166
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:172
int const *const pol
Definition: TC2.h:165
const int n
System size.
Definition: TC2.h:108
void purify()
Runs purification.
Definition: TC2.h:194
int * polys
Choices of polynomials 0 for x^2 and 1 for 2x-x^2 Length: nmul.
Definition: TC2.h:131
void print_data(int const start, int const stop) const
Definition: TC2.h:356
Help class for bisection root finding calls.
Definition: TC2.h:146
virtual ~TC2()
Destructor.
Definition: TC2.h:97
Treal template_blas_sqrt(Treal x)
Tmatrix & X
Fock / Kohn-Sham matrix at initialization.
Definition: TC2.h:103
Trace correcting purification.
Definition: TC2.h:62