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