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
00028
00029
00038 #ifndef MAT_TC2
00039 #define MAT_TC2
00040 #include <math.h>
00041 #include "bisection.h"
00042 namespace mat {
00063 template<typename Treal, typename Tmatrix>
00064 class TC2 {
00065 public:
00066 TC2(Tmatrix& F,
00067 Tmatrix& DM,
00068 const int size,
00069 const int noc,
00070 const Treal trunc = 0,
00072 const int maxmm = 100
00074 );
00078 Treal fermi_level(Treal tol = 1e-15
00079 ) const;
00083 Treal homo(Treal tol = 1e-15
00084 ) const;
00088 Treal lumo(Treal tol = 1e-15
00089 ) const;
00094 inline int n_multiplies() const {
00095 return nmul;
00096 }
00098 void print_data(int const start, int const stop) const;
00099 virtual ~TC2() {
00100 delete[] idemerror;
00101 delete[] tracediff;
00102 delete[] polys;
00103 }
00104 protected:
00105 Tmatrix& X;
00109 Tmatrix& D;
00110 const int n;
00111 const int nocc;
00112 const Treal frob_trunc;
00113 const int maxmul;
00114 Treal lmin;
00115 Treal lmax;
00116 int nmul;
00117 int nmul_firstpart;
00120 Treal* idemerror;
00129 Treal* tracediff;
00133 int* polys;
00136 void purify();
00139 private:
00140 class Fun;
00141 };
00147 template<typename Treal, typename Tmatrix>
00148 class TC2<Treal, Tmatrix>::Fun {
00149 public:
00150 Fun(int const* const p, int const pl, Treal const s)
00151 :pol(p), pollength(pl), shift(s) {
00152 assert(shift <= 1 && shift >= 0);
00153 assert(pollength >= 0);
00154 }
00155 Treal eval(Treal const x) const {
00156 Treal y = x;
00157 for (int ind = 0; ind < pollength; ind++ )
00158 y = 2 * pol[ind] * y + (1 - 2 * pol[ind]) * y * y;
00159
00160
00161
00162
00163 return y - shift;
00164 }
00165 protected:
00166 private:
00167 int const* const pol;
00168 int const pollength;
00169 Treal const shift;
00170 };
00171
00172
00173 template<typename Treal, typename Tmatrix>
00174 TC2<Treal, Tmatrix>::TC2(Tmatrix& F, Tmatrix& DM, const int size,
00175 const int noc,
00176 const Treal trunc, const int maxmm)
00177 :X(F), D(DM), n(size), nocc(noc), frob_trunc(trunc), maxmul(maxmm),
00178 lmin(0), lmax(0), nmul(0), nmul_firstpart(0),
00179 idemerror(0), tracediff(0), polys(0) {
00180 assert(frob_trunc >= 0);
00181 assert(nocc >= 0);
00182 assert(maxmul >= 0);
00183 X.gershgorin(lmin, lmax);
00184 X.add_identity(-lmax);
00185 X *= ((Treal)1.0 / (lmin - lmax));
00186 D = X;
00187 idemerror = new Treal[maxmul];
00188 tracediff = new Treal[maxmul + 1];
00189 polys = new int[maxmul];
00190 tracediff[0] = X.trace() - nocc;
00191 purify();
00192 }
00195 template<typename Treal, typename Tmatrix>
00196 void TC2<Treal, Tmatrix>::purify() {
00197 assert(nmul == 0);
00198 assert(nmul_firstpart == 0);
00199 Treal delta, beta, trD2;
00200 int ind;
00201 Treal const ONE = 1;
00202 Treal const TWO = 2;
00203 do {
00204 if (nmul >= maxmul) {
00205 print_data(0, nmul);
00206 throw AcceptableMaxIter("TC2<Treal, Tmatrix>::purify(): "
00207 "Purification reached maxmul"
00208 " without convergence", maxmul);
00209 }
00210 if (tracediff[nmul] > 0) {
00211 D = ONE * X * X;
00212 polys[nmul] = 0;
00213 }
00214 else {
00215 D = -ONE * X * X + TWO * D;
00216 polys[nmul] = 1;
00217 }
00218 D.frob_thresh(frob_trunc);
00219 idemerror[nmul] = Tmatrix::frob_diff(D, X);
00220 ++nmul;
00221 tracediff[nmul] = D.trace() - nocc;
00222 X = D;
00223
00224 beta = (3 - template_blas_sqrt(5)) / 2 - frob_trunc;
00225 if (idemerror[nmul - 1] < 1 / (Treal)4 &&
00226 (1 - template_blas_sqrt(1 - 4 * idemerror[nmul - 1])) / 2 < beta)
00227 beta = (1 + template_blas_sqrt(1 - 4 * idemerror[nmul - 1])) / 2;
00228 trD2 = (tracediff[nmul] + nocc -
00229 2 * polys[nmul - 1] * (tracediff[nmul - 1] + nocc)) /
00230 (1 - 2 * polys[nmul - 1]);
00231 delta = frob_trunc;
00232 ind = nmul - 1;
00233 while (ind > 0 && polys[ind] == polys[ind - 1]) {
00234 delta = delta + frob_trunc;
00235 ind--;
00236 }
00237 delta = delta < (template_blas_sqrt(1 + 4 * idemerror[nmul - 1]) - 1) / 2 ?
00238 delta : (template_blas_sqrt(1 + 4 * idemerror[nmul - 1]) - 1) / 2;
00239 } while((trD2 + beta * (1 + delta) * n - (1 + delta + beta) *
00240 (tracediff[nmul - 1] + nocc)) /
00241 ((1 + 2 * delta) * (delta + beta)) < n - nocc - 1 ||
00242 (trD2 - delta * (1 - beta) * n - (1 - delta - beta) *
00243 (tracediff[nmul - 1] + nocc)) /
00244 ((1 + 2 * delta) * (delta + beta)) < nocc - 1);
00245
00246
00247
00248
00249
00250
00251
00252 if (tracediff[nmul - 1] > 0) {
00253
00254
00255 D = -ONE * X * X + TWO * D;
00256 polys[nmul] = 1;
00257 }
00258 else {
00259 D = ONE * X * X;
00260 polys[nmul] = 0;
00261 }
00262 D.frob_thresh(frob_trunc);
00263 idemerror[nmul] = Tmatrix::frob_diff(D, X);
00264 ++nmul;
00265 tracediff[nmul] = D.trace() - nocc;
00266
00267 nmul_firstpart = nmul;
00268
00269
00270 do {
00271 if (nmul + 1 >= maxmul) {
00272 print_data(0, nmul);
00273 throw AcceptableMaxIter("TC2<Treal, Tmatrix>::purify(): "
00274 "Purification reached maxmul"
00275 " without convergence", maxmul);
00276 }
00277 if (tracediff[nmul] > 0) {
00278 X = ONE * D * D;
00279 idemerror[nmul] = Tmatrix::frob_diff(D, X);
00280 D = X;
00281 polys[nmul] = 0;
00282 ++nmul;
00283 tracediff[nmul] = D.trace() - nocc;
00284
00285 D = -ONE * X * X + TWO * D;
00286 idemerror[nmul] = Tmatrix::frob_diff(D, X);
00287 polys[nmul] = 1;
00288 ++nmul;
00289 tracediff[nmul] = D.trace() - nocc;
00290 }
00291 else {
00292 X = D;
00293 X = -ONE * D * D + TWO * X;
00294 idemerror[nmul] = Tmatrix::frob_diff(D, X);
00295 polys[nmul] = 1;
00296 ++nmul;
00297 tracediff[nmul] = X.trace() - nocc;
00298
00299 D = ONE * X * X;
00300 idemerror[nmul] = Tmatrix::frob_diff(D, X);
00301 polys[nmul] = 0;
00302 ++nmul;
00303 tracediff[nmul] = D.trace() - nocc;
00304 }
00305 D.frob_thresh(frob_trunc);
00306 #if 0
00307 } while (idemerror[nmul - 1] > frob_trunc);
00308 #else
00309 } while ((1 - template_blas_sqrt(1 - 4 * idemerror[nmul - 1])) / 2 > frob_trunc);
00310 #endif
00311 X.clear();
00312 }
00313
00314 template<typename Treal, typename Tmatrix>
00315 Treal TC2<Treal, Tmatrix>::fermi_level(Treal tol) const {
00316 Fun const fermifun(polys, nmul, 0.5);
00317 Treal chempot = bisection(fermifun, (Treal)0, (Treal)1, tol);
00318 return (lmin - lmax) * chempot + lmax;
00319 }
00320
00321 template<typename Treal, typename Tmatrix>
00322 Treal TC2<Treal, Tmatrix>::homo(Treal tol) const {
00323 Treal homo = 0;
00324 Treal tmp;
00325 for (int mul = nmul_firstpart; mul < nmul; mul++) {
00326 if (idemerror[mul] < 1.0 / 4) {
00327 Fun const homofun(polys, mul, (1 + template_blas_sqrt(1 - 4 * idemerror[mul])) / 2);
00328 tmp = bisection(homofun, (Treal)0, (Treal)1, tol);
00329
00330
00331
00332
00333 homo = tmp > homo ? tmp : homo;
00334 }
00335 }
00336 return (lmin - lmax) * homo + lmax;
00337 }
00338
00339 template<typename Treal, typename Tmatrix>
00340 Treal TC2<Treal, Tmatrix>::lumo(Treal tol) const {
00341 Treal lumo = 1;
00342 Treal tmp;
00343 for (int mul = nmul_firstpart; mul < nmul; mul++) {
00344 if (idemerror[mul] < 1.0 / 4) {
00345 Fun const lumofun(polys, mul, (1 - template_blas_sqrt(1 - 4 * idemerror[mul])) / 2);
00346 tmp = bisection(lumofun, (Treal)0, (Treal)1, tol);
00347
00348
00349
00350
00351 lumo = tmp < lumo ? tmp : lumo;
00352 }
00353 }
00354 return (lmin - lmax) * lumo + lmax;
00355 }
00356
00357 template<typename Treal, typename Tmatrix>
00358 void TC2<Treal, Tmatrix>::print_data(int const start, int const stop) const {
00359 for (int ind = start; ind < stop; ind ++) {
00360 std::cout << "Iteration: " << ind
00361 << " Idempotency error: " << idemerror[ind]
00362 << " Tracediff: " << tracediff[ind]
00363 << " Poly: " << polys[ind]
00364 << std::endl;
00365 }
00366 }
00367
00368
00369 }
00370 #endif