00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
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
00159
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);
00182 X.add_identity(-lmax);
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
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
00245
00246
00247
00248
00249
00250 if (tracediff[nmul - 1] > 0) {
00251
00252
00253 D = -ONE * X * X + TWO * D;
00254 polys[nmul] = 1;
00255 }
00256 else {
00257 D = ONE * X * X;
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;
00266
00267
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);
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
00329
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
00347
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 }
00368 #endif