ergo
TC2.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 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"
42 namespace 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 
357 template<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
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
Definition: Failure.h:81
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
Bisection method.
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