ergo
|
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