ergo
TC2.h
Go to the documentation of this file.
1/* Ergo, version 3.8, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2019 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
38#ifndef MAT_TC2
39#define MAT_TC2
40#include <math.h>
41#include "bisection.h"
42namespace mat {
63 template<typename Treal, typename Tmatrix>
64 class TC2 {
65 public:
66 TC2(Tmatrix& F,
67 Tmatrix& DM,
68 const int size,
69 const int noc,
70 const Treal trunc = 0,
72 const int maxmm = 100
74 );
78 Treal fermi_level(Treal tol = 1e-15
79 ) const;
83 Treal homo(Treal tol = 1e-15
84 ) const;
88 Treal lumo(Treal tol = 1e-15
89 ) const;
94 inline int n_multiplies() const {
95 return nmul;
96 }
98 void print_data(int const start, int const stop) const;
99 virtual ~TC2() {
100 delete[] idemerror;
101 delete[] tracediff;
102 delete[] polys;
103 }
104 protected:
105 Tmatrix& X;
109 Tmatrix& D;
110 const int n;
111 const int nocc;
112 const Treal frob_trunc;
113 const int maxmul;
114 Treal lmin;
115 Treal lmax;
116 int nmul;
120 Treal* idemerror;
129 Treal* tracediff;
133 int* polys;
136 void purify();
139 private:
140 class Fun;
141 };
147 template<typename Treal, typename Tmatrix>
148 class TC2<Treal, Tmatrix>::Fun {
149 public:
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);
154 }
155 Treal eval(Treal const x) const {
156 Treal y = x;
157 for (int ind = 0; ind < pollength; ind++ )
158 y = 2 * pol[ind] * y + (1 - 2 * pol[ind]) * y * y;
159 /*
160 * pol[ind] == 0 --> y = y * y
161 * pol[ind] == 1 --> y = 2 * y - y * y
162 */
163 return y - shift;
164 }
165 protected:
166 private:
167 int const* const pol;
168 int const pollength;
169 Treal const shift;
170 };
171
172
173 template<typename Treal, typename Tmatrix>
174 TC2<Treal, Tmatrix>::TC2(Tmatrix& F, Tmatrix& DM, const int size,
175 const int noc,
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) {
180 assert(frob_trunc >= 0);
181 assert(nocc >= 0);
182 assert(maxmul >= 0);
183 X.gershgorin(lmin, lmax); /* Find eigenvalue bounds */
184 X.add_identity(-lmax); /* Scale to [0, 1] interval and negate */
185 X *= ((Treal)1.0 / (lmin - lmax));
186 D = X;
187 idemerror = new Treal[maxmul];
188 tracediff = new Treal[maxmul + 1];
189 polys = new int[maxmul];
190 tracediff[0] = X.trace() - nocc;
191 purify();
192 }
195 template<typename Treal, typename Tmatrix>
197 assert(nmul == 0);
198 assert(nmul_firstpart == 0);
199 Treal delta, beta, trD2;
200 int ind;
201 Treal const ONE = 1;
202 Treal const TWO = 2;
203 do {
204 if (nmul >= maxmul) {
205 print_data(0, nmul);
206 throw AcceptableMaxIter("TC2<Treal, Tmatrix>::purify(): "
207 "Purification reached maxmul"
208 " without convergence", maxmul);
209 }
210 if (tracediff[nmul] > 0) {
211 D = ONE * X * X;
212 polys[nmul] = 0;
213 }
214 else {
215 D = -ONE * X * X + TWO * D;
216 polys[nmul] = 1;
217 }
218 D.frob_thresh(frob_trunc);
219 idemerror[nmul] = Tmatrix::frob_diff(D, X);
220 ++nmul;
221 tracediff[nmul] = D.trace() - nocc;
222 X = D;
223 /* Setting up convergence criteria */
224 beta = (3 - template_blas_sqrt(5)) / 2 - frob_trunc;
225 if (idemerror[nmul - 1] < 1 / (Treal)4 &&
226 (1 - template_blas_sqrt(1 - 4 * idemerror[nmul - 1])) / 2 < beta)
227 beta = (1 + template_blas_sqrt(1 - 4 * idemerror[nmul - 1])) / 2;
228 trD2 = (tracediff[nmul] + nocc -
229 2 * polys[nmul - 1] * (tracediff[nmul - 1] + nocc)) /
230 (1 - 2 * polys[nmul - 1]);
231 delta = frob_trunc;
232 ind = nmul - 1;
233 while (ind > 0 && polys[ind] == polys[ind - 1]) {
234 delta = delta + frob_trunc;
235 ind--;
236 }
237 delta = delta < (template_blas_sqrt(1 + 4 * idemerror[nmul - 1]) - 1) / 2 ?
238 delta : (template_blas_sqrt(1 + 4 * idemerror[nmul - 1]) - 1) / 2;
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);
245
246 /* Note that: */
247 /* tracediff[i] - tracediff[i - 1] = trace(D[i]) - trace(D[i - 1]) */
248 /* i.e. the change of the trace. */
249
250 /* Take one step to make sure the eigenvalues stays in */
251 /* { [ 0 , 2 * epsilon [ , ] 1 - 2 * epsilon , 1] } */
252 if (tracediff[nmul - 1] > 0) {
253 /* The same tracediff as in the last step is used since we want to */
254 /* take a step with the other direction (with the other polynomial).*/
255 D = -ONE * X * X + TWO * D; /* This is correct!! */
256 polys[nmul] = 1;
257 }
258 else {
259 D = ONE * X * X; /* This is correct!! */
260 polys[nmul] = 0;
261 }
262 D.frob_thresh(frob_trunc);
263 idemerror[nmul] = Tmatrix::frob_diff(D, X);
264 ++nmul;
265 tracediff[nmul] = D.trace() - nocc;
266
267 nmul_firstpart = nmul; /* First part of purification finished. At this */
268 /* point the eigenvalues are separated but have not yet converged. */
269 /* Use second order convergence polynomials to converge completely: */
270 do {
271 if (nmul + 1 >= maxmul) {
272 print_data(0, nmul);
273 throw AcceptableMaxIter("TC2<Treal, Tmatrix>::purify(): "
274 "Purification reached maxmul"
275 " without convergence", maxmul);
276 }
277 if (tracediff[nmul] > 0) {
278 X = ONE * D * D;
279 idemerror[nmul] = Tmatrix::frob_diff(D, X);
280 D = X;
281 polys[nmul] = 0;
282 ++nmul;
283 tracediff[nmul] = D.trace() - nocc;
284
285 D = -ONE * X * X + TWO * D;
286 idemerror[nmul] = Tmatrix::frob_diff(D, X);
287 polys[nmul] = 1;
288 ++nmul;
289 tracediff[nmul] = D.trace() - nocc;
290 }
291 else {
292 X = D;
293 X = -ONE * D * D + TWO * X;
294 idemerror[nmul] = Tmatrix::frob_diff(D, X);
295 polys[nmul] = 1;
296 ++nmul;
297 tracediff[nmul] = X.trace() - nocc;
298
299 D = ONE * X * X;
300 idemerror[nmul] = Tmatrix::frob_diff(D, X);
301 polys[nmul] = 0;
302 ++nmul;
303 tracediff[nmul] = D.trace() - nocc;
304 }
305 D.frob_thresh(frob_trunc);
306#if 0
307 } while (idemerror[nmul - 1] > frob_trunc); /* FIXME Check conv. crit. */
308#else
309 } while ((1 - template_blas_sqrt(1 - 4 * idemerror[nmul - 1])) / 2 > frob_trunc);
310#endif
311 X.clear();
312}
313
314 template<typename Treal, typename Tmatrix>
315 Treal TC2<Treal, Tmatrix>::fermi_level(Treal tol) const {
316 Fun const fermifun(polys, nmul, 0.5);
317 Treal chempot = bisection(fermifun, (Treal)0, (Treal)1, tol);
318 return (lmin - lmax) * chempot + lmax;
319 }
320
321 template<typename Treal, typename Tmatrix>
322 Treal TC2<Treal, Tmatrix>::homo(Treal tol) const {
323 Treal homo = 0;
324 Treal tmp;
325 for (int mul = nmul_firstpart; mul < nmul; mul++) {
326 if (idemerror[mul] < 1.0 / 4) {
327 Fun const homofun(polys, mul, (1 + template_blas_sqrt(1 - 4 * idemerror[mul])) / 2);
328 tmp = bisection(homofun, (Treal)0, (Treal)1, tol);
329 /*
330 std::cout << tmp << " , ";
331 std::cout << (lmin - lmax) * tmp + lmax << std::endl;
332 */
333 homo = tmp > homo ? tmp : homo;
334 }
335 }
336 return (lmin - lmax) * homo + lmax;
337 }
338
339 template<typename Treal, typename Tmatrix>
340 Treal TC2<Treal, Tmatrix>::lumo(Treal tol) const {
341 Treal lumo = 1;
342 Treal tmp;
343 for (int mul = nmul_firstpart; mul < nmul; mul++) {
344 if (idemerror[mul] < 1.0 / 4) {
345 Fun const lumofun(polys, mul, (1 - template_blas_sqrt(1 - 4 * idemerror[mul])) / 2);
346 tmp = bisection(lumofun, (Treal)0, (Treal)1, tol);
347 /*
348 std::cout << tmp << " , ";
349 std::cout << (lmin - lmax) * tmp + lmax << std::endl;
350 */
351 lumo = tmp < lumo ? tmp : lumo;
352 }
353 }
354 return (lmin - lmax) * lumo + lmax;
355 }
356
357template<typename Treal, typename Tmatrix>
358 void TC2<Treal, Tmatrix>::print_data(int const start, int const stop) const {
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]
364 << std::endl;
365 }
366}
367
368
369} /* end namespace mat */
370#endif
Bisection method.
Definition: Failure.h:81
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
Treal template_blas_sqrt(Treal x)