Kernel_Elec_Utils.h

Go to the documentation of this file.

 
#ifndef Kernel_Elec_Utils_H
#define Kernel_Elec_Utils_H
 
#include "kids/Kernel.h"
#include "kids/Kernel_Elec.h"
#include "kids/Kernel_Random.h"
#include "kids/linalg.h"
#include "kids/vars_list.h"
 
namespace PROJECT_NS {
 
DEFINE_POLICY(SQCPolicy,
              SQR,  // square window
              TRI,  // triangle window
              SPX,  // simplex window
              BIG   // simplex window
);
 
class elec_utils {
   public:
    static inline double gamma_wigner(int fdim) { return (sqrt(fdim + 1.0e0) - 1) / fdim; }
 
    static inline double gamma_opt(int fdim) {
        double sum = 0.0e0;
        for (int i = 0; i < fdim; ++i) sum += 1.0e0 / (i + 1);
        return ((fdim - sum) / (sum - 1.0e0)) / fdim;
    }
 
    static int max_choose(kids_complex* rho) {
        int       imax = 0;
        kids_real vmax = 0.0f;
        for (int i = 0, ii = 0; i < Dimension::F; ++i, ii += Dimension::Fadd1) {
            if (std::real(rho[ii]) > vmax) {
                vmax = std::real(rho[ii]);
                imax = i;
            }
        }
        return imax;
    }
 
    static int pop_choose(kids_complex* rho) {
        kids_real rand_tmp;
        kids_real sum = 0.0f;
        Kernel_Random::rand_uniform(&rand_tmp);
        for (int i = 0, ii = 0; i < Dimension::F; ++i, ii += Dimension::Fadd1) {
            sum += std::min({std::max({std::real(rho[ii]), 0.0e0}), 1.0e0});
            if (rand_tmp < sum) return i;
        }
        return 0;
    }
 
    static int pop_neg_choose(kids_complex* rho) {
        kids_real rand_tmp;
        kids_real sum = 0.0f;
        for (int i = 0, ii = 0; i < Dimension::F; ++i, ii += Dimension::Fadd1) sum += std::abs(rho[ii]);
        Kernel_Random::rand_uniform(&rand_tmp);
        rand_tmp *= sum;
        sum = 0.0e0;
        for (int i = 0, ii = 0; i < Dimension::F; ++i, ii += Dimension::Fadd1) {
            sum += std::abs(rho[ii]);
            if (rand_tmp < sum) return i;
        }
        return 0;
    }
 
    static int hopping_choose(kids_complex* rho, kids_complex* H, int from, kids_real dt) {
        int       to = from;
        kids_real rand_tmp, sumprob = 0.0f;
        kids_real rhoii = std::real(rho[from * Dimension::Fadd1]);
        Kernel_Random::rand_uniform(&rand_tmp);
 
        for (int n = 0; n < Dimension::F; ++n) {
            kids_real prob =
                (n == from) ? 0.0f
                            : -2.0f * std::imag(rho[n * Dimension::F + from] * H[from * Dimension::F + n]) / rhoii * dt;
            prob = (prob > 1.0f) ? 1.0f : ((prob < 0.0f) ? 0.0f : prob);  // hopping cut-off
            sumprob += prob;
            if (rand_tmp < sumprob) {
                to = n;
                break;
            }
        }
        return to;
    }
 
    static void hopping_direction(kids_real* direction, kids_real* dE, int from, int to) {
        if (to == from) return;
        for (int i = 0; i < Dimension::N; ++i) { direction[i] = dE[i * Dimension::FF + from * Dimension::F + to]; }
    }
 
    static void hopping_direction(kids_real* direction, kids_real* E, kids_real* dE, kids_complex* rho, int from,
                                  int to) {
        if (to == from) return;
 
        for (int i = 0; i < Dimension::N; ++i) {
            direction[i] = 0.0f;
            for (int k = 0; k < Dimension::F; ++k)
                if (k != from)
                    direction[i] +=
                        std::real(rho[from * Dimension::F + k] * dE[i * Dimension::FF + k * Dimension::F + from]) /
                        (E[from] - E[k]);
            for (int k = 0; k < Dimension::F; ++k)
                if (k != to)
                    direction[i] -=
                        std::real(rho[to * Dimension::F + k] * dE[i * Dimension::FF + k * Dimension::F + to]) /
                        (E[to] - E[k]);
        }
    }
 
    static int hopping_impulse(kids_real* direction, kids_real* np, kids_real* nm, kids_real* E,  //
                               int from, int to, bool reflect) {
        if (to == from) return from;
 
        // solve x: Ef + P**2 / (2*M) = Et + (P + direction*x)**2 / (2*M)
        kids_real coeffa = 0.0f, coeffb = 0.0f, coeffc = E[to] - E[from];
        for (int i = 0; i < Dimension::N; ++i) {
            coeffa += 0.5f * direction[i] * direction[i] / nm[i];
            coeffb += np[i] / nm[i] * direction[i];
        }
        coeffb /= coeffa, coeffc /= coeffa;  // normalization for safety
 
        kids_real coeffd = coeffb * coeffb - 4 * coeffc;
        if (coeffd > 0) {
            kids_real x1 = 0.5f * (-coeffb + sqrt(coeffd)), x2 = 0.5f * (-coeffb - sqrt(coeffd));
            kids_real xx = (std::abs(x1) < std::abs(x2)) ? x1 : x2;
            for (int i = 0; i < Dimension::N; ++i) np[i] += xx * direction[i];
            return to;
        } else if (reflect) {  // 2008Algorithm
            kids_real xx = -coeffb;
            for (int i = 0; i < Dimension::N; ++i) np[i] += xx * direction[i];
            return from;
        } else {  // 1990Algorithm, do nothing
            return from;
        }
        return from;
    }
 
    static int c_sphere(kids_complex* c, int fdim) {
        Kernel_Random::rand_gaussian(reinterpret_cast<double*>(c), 2 * fdim);
        double l = std::sqrt(std::real(ARRAY_INNER_TRANS1(c, c, fdim)));
        for (int i = 0; i < fdim; ++i) c[i] /= l;
        return 0;
    }
 
    static int c_focus(kids_complex* c, double xi, double gamma, int occ, int fdim) {
        for (int i = 0; i < fdim; ++i) c[i] = ((i == occ ? 1.0e0 : 0.0e0) + gamma) / xi;
        double randu;
        for (int i = 0; i < fdim; ++i) {
            Kernel_Random::rand_uniform(&randu);
            randu *= phys::math::twopi;
            c[i] = std::sqrt(std::abs(c[i])) * (cos(randu) + phys::math::im * sin(randu));
        }
        return 0;
    };
 
    static int rho_focus(kids_complex* rho, int iocc, double gamma_ou, double gamma_uu, int fdim, bool rand_act,
                         bool pure_phase, bool cont_phase) {
        double randu = 1.0e0;
        if (pure_phase) {
            for (int j = 0; j < fdim; ++j) {
                if (j == iocc) {
                    rho[iocc * fdim + j] = 1.0e0;
                    continue;
                }
                Kernel_Random::rand_uniform(&randu);
                randu = (cont_phase) ? phys::math::twopi * randu : phys::math::halfpi * (int(randu / 0.25f) + 1);
                rho[iocc * fdim + j] = cos(randu) + phys::math::im * sin(randu);
                rho[j * fdim + iocc] = std::conj(rho[iocc * fdim + j]);
            }
            for (int i = 0, ij = 0; i < fdim; ++i) {
                for (int j = 0; j < fdim; ++j, ++ij) {
                    if (i == iocc || j == iocc) continue;
                    rho[ij] = rho[iocc * fdim + j] / rho[iocc * fdim + i];
                }
            }
        } else {
            for (int i = 0; i < fdim; ++i) {
                for (int j = i + 1; j < fdim; ++j) {
                    Kernel_Random::rand_uniform(&randu);
                    randu = (cont_phase) ? phys::math::twopi * randu : phys::math::halfpi * (int(randu / 0.25f) + 1);
                    rho[i * fdim + j] = cos(randu) + phys::math::im * sin(randu);
                    rho[j * fdim + i] = std::conj(rho[i * fdim + j]);
                }
            }
        }
        Kernel_Random::rand_uniform(&randu);
        double occrand = (rand_act) ? int(randu * fdim) : iocc;
        for (int i = 0, ij = 0; i < fdim; ++i) {
            for (int j = 0; j < fdim; ++j, ++ij) {
                if (i == j) {
                    rho[ij] = (i == occrand) ? phys::math::iu : phys::math::iz;
                } else if (i == occrand || j == occrand) {
                    rho[ij] *= gamma_ou;
                } else {
                    rho[ij] *= gamma_uu;
                }
            }
        }
        return 0;
    }
 
    static int c_window(kids_complex* c, int iocc, int type, int fdim) {
        switch (type) {
            case SQCPolicy::TRI: {
                kids_real tmp2[2];
                Kernel_Random::rand_uniform(tmp2, 2);
                while (tmp2[0] + tmp2[1] > 1.0f) Kernel_Random::rand_uniform(tmp2, 2);
                c[iocc] = tmp2[0];
                tmp2[1] = 1.0f - std::real(c[iocc]);
                for (int i = 0; i < fdim; ++i) {
                    if (i != iocc) {
                        Kernel_Random::rand_uniform(tmp2, 1);
                        c[i] = tmp2[0] * tmp2[1];
                    }
                }
                c[iocc] += 1.0e0;
                break;
            }
            case SQCPolicy::SPX: {
                c_sphere(c, Dimension::F);
                for (int i = 0; i < Dimension::F; ++i) c[i] = std::abs(c[i] * c[i]);
                c[iocc] += 1.0e0;
                break;
            }
            case SQCPolicy::BIG: {
                kids_complex* cadd1 = new kids_complex[Dimension::Fadd1];
                c_sphere(cadd1, Dimension::Fadd1);
                for (int i = 0; i < Dimension::F; ++i) c[i] = std::abs(cadd1[i] * cadd1[i]);
                c[iocc] += 1.0e0;
                delete[] cadd1;
                break;
            }
            case SQCPolicy::SQR: {
                const kids_real gm0 = gamma_wigner(2.0f);
                for (int i = 0; i < fdim; ++i) {
                    kids_real randu;
                    Kernel_Random::rand_uniform(&randu);
                    c[i] = 2.0 * gm0 * randu;
                }
                c[iocc] += 1.0e0;
                break;
            }
        }
        for (int i = 0; i < fdim; ++i) {
            kids_real randu;
            Kernel_Random::rand_uniform(&randu);
            randu *= phys::math::twopi;
            c[i] = sqrt(c[i]);
            c[i] *= (cos(randu) + phys::math::im * sin(randu));
        }
        return 0;
    };
 
    static int ker_binning(kids_complex* ker, kids_complex* rho, int sqc_type) {
        const kids_real gm0 = gamma_wigner(2.0f), gm1 = 1 + gm0, gmh = 0.5f + gm0;
 
        // set all elements to 1
        for (int i = 0; i < Dimension::FF; ++i) ker[i] = rho[i] / std::abs(rho[i]);
 
        // for ker[ij], loop k-index to find if ker[ij] should be set to 0
        for (int i = 0, ij = 0; i < Dimension::F; ++i) {
            for (int j = 0; j < Dimension::F; ++j, ++ij) {
                for (int k = 0, kk = 0; k < Dimension::F; ++k, kk += Dimension::Fadd1) {
                    double vk      = std::abs(rho[kk]);
                    bool   Outlier = false;
                    switch (sqc_type) {
                        case SQCPolicy::TRI:
                        case SQCPolicy::SPX:
                        case SQCPolicy::BIG:
                            Outlier = (i == j) ? ((k != i && vk > 1) || (k == i && vk < 1))
                                               : ((k != i && k != j && vk > 1) || ((k == i || k == j) && vk < 0.5f));
 
                            if (i != j) Outlier = false;
                            break;
                        case SQCPolicy::SQR:  // @bug?
                            Outlier = (i == j) ? ((k != i && std::abs(vk - gm0) < gm0) ||
                                                  (k == i && std::abs(vk - gm1) < gm0))
                                               : ((k != i && std::abs(vk - gm0) > gm0) ||  //
                                                  (k == i && std::abs(vk - gmh) > gm0) ||  //
                                                  (k == j && std::abs(vk - gmh) > gm0));
                            break;
                    }
                    if (Outlier) {
                        ker[ij] = phys::math::iz;
                        break;
                    }
                }
            }
        }
        return 0;
    };
};
 
};  // namespace PROJECT_NS
 
#endif  // Kernel_Elec_Utils_H