functions.hpp

Go to the documentation of this file.

#ifndef XPED_FUNCTIONS_H_
#define XPED_FUNCTIONS_H_
 
#include "Xped/Symmetry/SU2Wrappers.hpp"
#include "Xped/Symmetry/qarray.hpp"
#include "Xped/Util/Constfct.hpp"
 
#include <boost/rational.hpp>
 
namespace Xped {
 
namespace util {
 
inline std::string print_frac_nice(boost::rational<int> r)
{
    std::stringstream ss;
    if(r.denominator() == 1) {
        ss << r.numerator();
    } else {
        ss << r;
    }
    return ss.str();
}
 
} // namespace util
 
namespace Sym {
#ifndef KIND_ENUM
#    define KIND_ENUM
enum KIND
{
    S,
    Salt,
    T,
    FT,
    N,
    FN,
    M,
    Nup,
    Ndn,
    Z2
};
#endif
 
template <typename Symmetry>
std::string format(qarray<Symmetry::Nq> qnum)
{
    std::stringstream ss;
    for(int q = 0; q < Symmetry::Nq; ++q) {
        if(Symmetry::kind()[q] == KIND::S or Symmetry::kind()[q] == KIND::Salt or Symmetry::kind()[q] == KIND::T) {
            ss << util::print_frac_nice(boost::rational<int>(qnum[q] - 1, 2));
        } else if(Symmetry::kind()[q] == KIND::M) {
            if(Symmetry::IS_MODULAR[q]) {
                int q_nice = qnum[q];
                if(std::abs(qnum[q]) > std::abs(qnum[q] - Symmetry::MOD_N[q])) { q_nice = q_nice - Symmetry::MOD_N[q]; }
                ss << util::print_frac_nice(boost::rational<int>(q_nice, 2));
            } else {
                ss << util::print_frac_nice(boost::rational<int>(qnum[q], 2));
            }
        } else if(Symmetry::kind()[q] == KIND::Z2) {
            std::string parity = util::constFct::posmod<2>(qnum[q] == 0) ? "evn" : "odd";
            ss << parity;
        } else {
            ss << qnum[q];
        }
        if(q != Symmetry::Nq - 1) { ss << ","; }
    }
    if(Symmetry::Nq == 0) { ss << "()"; }
    return ss.str();
}
 
template <typename Scalar>
Scalar phase(int q)
{
    if(q % 2) { return Scalar(-1.); }
    return Scalar(1.);
}
 
template <typename Symmetry>
std::vector<std::pair<typename Symmetry::qType, typename Symmetry::qType>>
split(const typename Symmetry::qType Q, const std::vector<typename Symmetry::qType>& ql, const std::vector<typename Symmetry::qType> qr)
{
    std::vector<std::pair<typename Symmetry::qType, typename Symmetry::qType>> vout;
    for(std::size_t q1 = 0; q1 < ql.size(); q1++)
        for(std::size_t q2 = 0; q2 < qr.size(); q2++) {
            auto Qs = Symmetry::reduceSilent(ql[q1], qr[q2]);
            if(auto it = std::find(Qs.begin(), Qs.end(), Q) != Qs.end()) { vout.push_back({ql[q1], qr[q2]}); }
        }
    return vout;
}
 
template <typename Symmetry>
std::vector<std::pair<std::size_t, std::size_t>>
split(const typename Symmetry::qType Q, const std::vector<typename Symmetry::qType>& ql, const std::vector<typename Symmetry::qType> qr)
{
    std::vector<std::pair<std::size_t, std::size_t>> vout;
    for(std::size_t q1 = 0; q1 < ql.size(); q1++)
        for(std::size_t q2 = 0; q2 < qr.size(); q2++) {
            auto Qs = Symmetry::reduceSilent(ql[q1], qr[q2]);
            if(auto it = std::find(Qs.begin(), Qs.end(), Q) != Qs.end()) { vout.push_back({q1, q2}); }
        }
    return vout;
}
 
inline void initialize(int maxJ = 1, std::string f_3j = "", std::string f_6j = "", std::string f_9j = "")
{
    std::ignore = maxJ;
    std::ignore = f_3j;
    std::ignore = f_6j;
    std::ignore = f_9j;
#ifdef USE_WIG_SU2_COEFFS
    wig_table_init(2 * maxJ, 9);
    wig_temp_init(2 * maxJ);
#endif
 
#ifdef USE_FAST_WIG_SU2_COEFFS
    fastwigxj_load(f_3j.c_str(), 3, NULL);
    fastwigxj_load(f_6j.c_str(), 6, NULL);
    fastwigxj_load(f_9j.c_str(), 9, NULL);
 
    wig_table_init(2 * maxJ, 9);
    wig_temp_init(2 * maxJ);
#endif
}
 
inline void finalize(bool PRINT_STATS = false)
{
    std::ignore = PRINT_STATS;
#ifdef USE_WIG_SU2_COEFFS
    wig_temp_free();
    wig_table_free();
#endif
 
#ifdef USE_FAST_WIG_SU2_COEFFS
    if(PRINT_STATS) {
        std::cout << std::endl;
        fastwigxj_print_stats();
    }
 
    fastwigxj_unload(3);
    fastwigxj_unload(6);
    fastwigxj_unload(9);
 
    wig_temp_free();
    wig_table_free();
#endif
}
 
} // namespace Sym
 
template <typename Symmetry>
qarray<Symmetry::Nq> adjustQN(const qarray<Symmetry::Nq>& qin, const size_t number_cells, bool BACK = false)
{
    qarray<Symmetry::Nq> out;
    for(size_t q = 0; q < Symmetry::Nq; ++q) {
        if(Symmetry::kind()[q] != Sym::KIND::S and Symmetry::kind()[q] != Sym::KIND::Salt and
           Symmetry::kind()[q] != Sym::KIND::T) // Do not transform the base for non-Abelian symmetries
        {
            if(BACK) {
                out[q] = qin[q] / number_cells;
            } else {
                out[q] = qin[q] * number_cells;
            }
        } else {
            out[q] = qin[q];
        }
    }
    return out;
};
 
} // namespace Xped
#endif