TC2.h

Go to the documentation of this file.
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure
00002  * calculations.
00003  * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
00004  * 
00005  * This program is free software: you can redistribute it and/or modify
00006  * it under the terms of the GNU General Public License as published by
00007  * the Free Software Foundation, either version 3 of the License, or
00008  * (at your option) any later version.
00009  * 
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  * 
00015  * You should have received a copy of the GNU General Public License
00016  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
00017  * 
00018  * Primary academic reference:
00019  * Kohn−Sham Density Functional Theory Electronic Structure Calculations 
00020  * with Linearly Scaling Computational Time and Memory Usage,
00021  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
00022  * J. Chem. Theory Comput. 7, 340 (2011),
00023  * <http://dx.doi.org/10.1021/ct100611z>
00024  * 
00025  * For further information about Ergo, see <http://www.ergoscf.org>.
00026  */
00027 
00036 #ifndef MAT_TC2
00037 #define MAT_TC2
00038 #include <math.h>
00039 #include "bisection.h"
00040 namespace mat {
00061   template<typename Treal, typename Tmatrix>
00062     class TC2 {
00063   public:    
00064     TC2(Tmatrix& F, 
00065         Tmatrix& DM, 
00066         const int size, 
00067         const int noc, 
00068         const Treal trunc = 0,
00070         const int maxmm = 100 
00072         ); 
00076     Treal fermi_level(Treal tol = 1e-15 
00077                       ) const; 
00081     Treal homo(Treal tol = 1e-15 
00082                ) const;
00086     Treal lumo(Treal tol = 1e-15 
00087                ) const;
00092     inline int n_multiplies() const {
00093       return nmul;
00094     }
00096     void print_data(int const start, int const stop) const;
00097     virtual ~TC2() {
00098       delete[] idemerror;
00099       delete[] tracediff;
00100       delete[] polys;
00101     } 
00102   protected:
00103     Tmatrix& X; 
00107     Tmatrix& D; 
00108     const int n; 
00109     const int nocc;   
00110     const Treal frob_trunc; 
00111     const int maxmul; 
00112     Treal lmin; 
00113     Treal lmax; 
00114     int nmul; 
00115     int nmul_firstpart; 
00118     Treal* idemerror; 
00127     Treal* tracediff; 
00131     int* polys; 
00134     void purify(); 
00137   private:
00138     class Fun;
00139   };
00145   template<typename Treal, typename Tmatrix>
00146     class TC2<Treal, Tmatrix>::Fun {
00147   public:
00148     Fun(int const* const p, int const pl, Treal const s)
00149       :pol(p), pollength(pl), shift(s) {
00150       assert(shift <= 1 && shift >= 0);
00151       assert(pollength >= 0);
00152     }
00153       Treal eval(Treal const x) const {
00154         Treal y = x;
00155         for (int ind = 0; ind < pollength; ind++ )
00156           y = 2 * pol[ind] * y + (1 - 2 * pol[ind]) * y * y;
00157         /*
00158          * pol[ind] == 0 --> y = y * y
00159          * pol[ind] == 1 --> y = 2 * y - y * y
00160          */
00161         return y - shift;
00162       }
00163   protected:
00164   private:
00165       int const* const pol;
00166       int const pollength;
00167       Treal const shift;
00168   };
00169   
00170 
00171   template<typename Treal, typename Tmatrix>
00172     TC2<Treal, Tmatrix>::TC2(Tmatrix& F, Tmatrix& DM, const int size, 
00173                              const int noc, 
00174                              const Treal trunc, const int maxmm)
00175     :X(F), D(DM), n(size), nocc(noc), frob_trunc(trunc), maxmul(maxmm), 
00176     lmin(0), lmax(0), nmul(0), nmul_firstpart(0), 
00177     idemerror(0), tracediff(0), polys(0) {
00178     assert(frob_trunc >= 0);
00179     assert(nocc >= 0);
00180     assert(maxmul >= 0);
00181     X.gersgorin(lmin, lmax);    /* Find eigenvalue bounds              */
00182     X.add_identity(-lmax);      /* Scale to [0, 1] interval and negate */
00183     X *= ((Treal)1.0 / (lmin - lmax));
00184     D = X;
00185     idemerror = new Treal[maxmul];
00186     tracediff = new Treal[maxmul + 1];
00187     polys = new int[maxmul];
00188     tracediff[0] = X.trace() - nocc;
00189     purify(); 
00190   } 
00193   template<typename Treal, typename Tmatrix>
00194     void TC2<Treal, Tmatrix>::purify() {
00195     assert(nmul == 0);
00196     assert(nmul_firstpart == 0);
00197     Treal delta, beta, trD2;
00198     int ind;
00199     Treal const ONE = 1;
00200     Treal const TWO = 2;
00201     do {
00202       if (nmul >= maxmul) {
00203         print_data(0, nmul);
00204         throw AcceptableMaxIter("TC2<Treal, Tmatrix>::purify(): "
00205                                 "Purification reached maxmul"
00206                                 " without convergence", maxmul);
00207       }
00208       if (tracediff[nmul] > 0) {
00209         D = ONE * X * X;
00210         polys[nmul] = 0;
00211       }
00212       else {
00213         D = -ONE * X * X + TWO * D;
00214         polys[nmul] = 1;
00215       }
00216       D.frob_thresh(frob_trunc);
00217       idemerror[nmul] = Tmatrix::frob_diff(D, X);
00218       ++nmul;
00219       tracediff[nmul] = D.trace() - nocc;
00220       X = D;
00221       /* Setting up convergence criteria */
00222       beta = (3 - template_blas_sqrt(5)) / 2 - frob_trunc;
00223       if (idemerror[nmul - 1] < 1 / (Treal)4 && 
00224           (1 - template_blas_sqrt(1 - 4 * idemerror[nmul - 1])) / 2 < beta)
00225         beta = (1 + template_blas_sqrt(1 - 4 * idemerror[nmul - 1])) / 2;
00226       trD2  = (tracediff[nmul] + nocc - 
00227                2 * polys[nmul - 1] * (tracediff[nmul - 1] + nocc)) / 
00228         (1 - 2 * polys[nmul - 1]);
00229       delta = frob_trunc;
00230       ind = nmul - 1;
00231       while (ind > 0 && polys[ind] == polys[ind - 1]) {
00232         delta = delta + frob_trunc;
00233         ind--;
00234       }
00235       delta = delta < (template_blas_sqrt(1 + 4 * idemerror[nmul - 1]) - 1) / 2 ?
00236         delta : (template_blas_sqrt(1 + 4 * idemerror[nmul - 1]) - 1) / 2;
00237     } while((trD2 + beta * (1 + delta) * n - (1 + delta + beta) * 
00238              (tracediff[nmul - 1] + nocc)) / 
00239             ((1 + 2 * delta) * (delta + beta)) < n - nocc - 1 ||
00240             (trD2 - delta * (1 - beta) * n - (1 - delta - beta) * 
00241              (tracediff[nmul - 1] + nocc)) / 
00242             ((1 + 2 * delta) * (delta + beta)) < nocc - 1);
00243     
00244     /* Note that: */
00245     /* tracediff[i] - tracediff[i - 1] = trace(D[i]) - trace(D[i - 1]) */
00246     /* i.e. the change of the trace. */
00247     
00248     /* Take one step to make sure the eigenvalues stays in */
00249     /* { [ 0 , 2 * epsilon [ , ] 1 - 2 * epsilon , 1] }    */
00250     if (tracediff[nmul - 1] > 0) { 
00251       /* The same tracediff as in the last step is used since we want to */
00252       /* take a step with the other direction (with the other polynomial).*/
00253       D = -ONE * X * X + TWO * D; /* This is correct!! */
00254       polys[nmul] = 1;
00255     }
00256     else {
00257       D = ONE * X * X; /* This is correct!! */
00258       polys[nmul] = 0;
00259     }
00260     D.frob_thresh(frob_trunc);
00261     idemerror[nmul] = Tmatrix::frob_diff(D, X);
00262     ++nmul;
00263     tracediff[nmul] = D.trace() - nocc;
00264 
00265     nmul_firstpart = nmul; /* First part of purification finished. At this */
00266     /* point the eigenvalues are separated but have not yet converged.     */
00267     /* Use second order convergence polynomials to converge completely:    */
00268     do { 
00269       if (nmul + 1 >= maxmul) {
00270         print_data(0, nmul);
00271         throw AcceptableMaxIter("TC2<Treal, Tmatrix>::purify(): "
00272                                 "Purification reached maxmul"
00273                                 " without convergence", maxmul);
00274       }
00275       if (tracediff[nmul] > 0) {
00276         X = ONE * D * D;
00277         idemerror[nmul] = Tmatrix::frob_diff(D, X);
00278         D = X;
00279         polys[nmul] = 0;
00280         ++nmul;
00281         tracediff[nmul] = D.trace() - nocc;
00282         
00283         D = -ONE * X * X + TWO * D;
00284         idemerror[nmul] = Tmatrix::frob_diff(D, X);
00285         polys[nmul] = 1;
00286         ++nmul;
00287         tracediff[nmul] = D.trace() - nocc;
00288       }
00289       else {
00290         X = D;
00291         X = -ONE * D * D + TWO * X;
00292         idemerror[nmul] = Tmatrix::frob_diff(D, X);
00293         polys[nmul] = 1;
00294         ++nmul;
00295         tracediff[nmul] = X.trace() - nocc;
00296         
00297         D = ONE * X * X;
00298         idemerror[nmul] = Tmatrix::frob_diff(D, X);
00299         polys[nmul] = 0;
00300         ++nmul;
00301         tracediff[nmul] = D.trace() - nocc;
00302       }
00303       D.frob_thresh(frob_trunc);
00304 #if 0
00305     } while (idemerror[nmul - 1] > frob_trunc); /* FIXME Check conv. crit. */
00306 #else
00307   } while ((1 - template_blas_sqrt(1 - 4 * idemerror[nmul - 1])) / 2 > frob_trunc); 
00308 #endif
00309   X.clear();
00310 }
00311 
00312   template<typename Treal, typename Tmatrix>
00313     Treal TC2<Treal, Tmatrix>::fermi_level(Treal tol) const {
00314     Fun const fermifun(polys, nmul, 0.5);
00315     Treal chempot = bisection(fermifun, (Treal)0, (Treal)1, tol);
00316     return (lmin - lmax) * chempot + lmax;
00317   }
00318 
00319   template<typename Treal, typename Tmatrix>
00320     Treal TC2<Treal, Tmatrix>::homo(Treal tol) const {
00321     Treal homo = 0;
00322     Treal tmp;
00323     for (int mul = nmul_firstpart; mul < nmul; mul++) {
00324       if (idemerror[mul] < 1.0 / 4) {
00325         Fun const homofun(polys, mul, (1 + template_blas_sqrt(1 - 4 * idemerror[mul])) / 2);
00326         tmp = bisection(homofun, (Treal)0, (Treal)1, tol);
00327         /*
00328           std::cout <<  tmp << " ,   ";
00329           std::cout << (lmin - lmax) * tmp + lmax << std::endl;
00330         */
00331         homo = tmp > homo ? tmp : homo;
00332       }
00333     }
00334     return  (lmin - lmax) * homo + lmax;
00335   }
00336   
00337   template<typename Treal, typename Tmatrix>
00338     Treal  TC2<Treal, Tmatrix>::lumo(Treal tol) const {
00339     Treal lumo = 1;
00340     Treal tmp;
00341     for (int mul = nmul_firstpart; mul < nmul; mul++) {
00342       if (idemerror[mul] < 1.0 / 4) {
00343         Fun const lumofun(polys, mul, (1 - template_blas_sqrt(1 - 4 * idemerror[mul])) / 2);
00344         tmp = bisection(lumofun, (Treal)0, (Treal)1, tol);
00345         /*
00346           std::cout <<  tmp << " ,   ";
00347           std::cout << (lmin - lmax) * tmp + lmax << std::endl;
00348         */
00349         lumo = tmp < lumo ? tmp : lumo;
00350       }
00351     }
00352     return  (lmin - lmax) * lumo + lmax;
00353   }
00354    
00355 template<typename Treal, typename Tmatrix>
00356   void TC2<Treal, Tmatrix>::print_data(int const start, int const stop) const {
00357   for (int ind = start; ind < stop; ind ++) {
00358     std::cout << "Iteration: " << ind
00359               << "   Idempotency error: " << idemerror[ind]
00360               << "   Tracediff: " << tracediff[ind]
00361               << "   Poly: " << polys[ind] 
00362               << std::endl;
00363   } 
00364 }
00365 
00366 
00367 } /* end namespace mat */
00368 #endif

Generated on 21 Nov 2012 for ergo by  doxygen 1.6.1