ParitySector.hpp

Go to the documentation of this file.

#pragma once
 
#include "Headers/mytypes.hpp"
#include "HilbertSpace.hpp"
#include "Eigen/SparseCore"
 
template<class ManyBodySpace, typename RealType = Real_t>
class ParitySector : public SubSpace<ManyBodySpace, RealType> {
    private:
        int m_parity;
        int m_NumParityEigens = 0;
 
        __host__ __device__ int numParityPairs() {
            assert((this->dim() - m_NumParityEigens) % 2 == 0);
            return (this->dim() - m_NumParityEigens) / 2;
        }
 
    public:
        using Base = SubSpace<ManyBodySpace, RealType>;
        template<typename std::enable_if_t< std::is_convertible_v<
                     ManyBodySpace, ManyBodySpaceBase<ManyBodySpace>> >* = nullptr>
        __host__ __device__ ParitySector(int parity, ManyBodySpace const& mbHSpace)
            : Base(mbHSpace), m_parity(parity) {
            debug_constructor_printf(1);
            assert(parity == +1 || parity == -1);
            if(this->totalSpace().sysSize() == 0) { return; }
 
            Eigen::ArrayX<bool> calculated = Eigen::ArrayX<bool>::Constant(mbHSpace.dim(), false);
            for(int state = 0; state != mbHSpace.dim(); ++state) {
                if(calculated[state]) continue;
                calculated[state] = true;
 
                int reversed = mbHSpace.reverse(state);
                if(reversed == state)
                    m_NumParityEigens += 1;
                else
                    calculated[reversed] = true;
            }
            int dim = (mbHSpace.dim() - m_NumParityEigens) / 2;
            if(m_parity == +1) dim += m_NumParityEigens;
            debug_print("mbHSpace.dim() = " << mbHSpace.dim() << ", dim = " << dim
                                            << ", m_NumParityEigens = " << m_NumParityEigens);
 
            this->basis().resize(mbHSpace.dim(), dim);
            this->basis().reserve(Eigen::VectorXi::Constant(dim, 2));
            dim        = 0;
            calculated = Eigen::ArrayX<bool>::Constant(mbHSpace.dim(), false);
            Eigen::ArrayXi parityEigen(m_NumParityEigens);
            m_NumParityEigens = 0;
            for(int state = 0; state != mbHSpace.dim(); ++state) {
                if(calculated[state]) continue;
                calculated[state] = true;
 
                int reversed = mbHSpace.reverse(state);
                if(reversed == state) { parityEigen(m_NumParityEigens++) = state; }
                else {
                    this->basis().insert(state, dim)    = 1.0 / sqrt(2.0);
                    this->basis().insert(reversed, dim) = 1.0 / sqrt(2.0);
                    calculated[reversed]                = true;
                    dim += 1;
                }
            }
            if(m_parity == 1) {
                for(int j = 0; j != parityEigen.size(); ++j)
                    this->basis().insert(parityEigen[j], dim++) = 1;
            }
 
            debug_print("mbHSpace.dim() = " << mbHSpace.dim() << ", dim = " << dim
                                            << ", this->dim() = " << this->dim());
            this->basis().makeCompressed();
            assert(dim == this->dim());
        };
        template<typename... Args>
        __host__ __device__ ParitySector(int parity, int systemSize, Args... args)
            : ParitySector(parity, ManyBodySpace(systemSize, args...)) {
            debug_constructor_printf((Default));
        }
 
        __host__ __device__ int parity() const { return m_parity; };
 
        __host__ __device__ int rep(int j, int trans = 0) const {
            assert(j < this->dim());
            int const period
                = this->basis().outerIndexPtr()[j + 1] - this->basis().outerIndexPtr()[j];
            int const innerId = this->basis().outerIndexPtr()[j] + (trans % period);
            assert(innerId < this->basis().nonZeros());
            return this->basis().innerIndexPtr()[innerId];
        };
 
        __host__ __device__ int period(int j) const {
            assert(j < this->dim() && "j < this->dim()");
            return this->basis().outerIndexPtr()[j + 1] - this->basis().outerIndexPtr()[j];
        };
};