symmetry_rotation.h

Go to the documentation of this file.

#pragma once
#include "irreducible_sector.h"
#include "source_basis/module_ao/parallel_orbitals.h"
#include <RI/global/Tensor.h>
#include "source_lcao/module_hcontainer/hcontainer.h"
#include "source_cell/module_neighbor/sltk_grid_driver.h"
 
namespace ModuleSymmetry
{
    using Tap = std::pair<int, int>;
    using TC = std::array<int, 3>;
    using TapR = std::pair<Tap, TC>;
    using TCdouble = Abfs::Vector3_Order<double>;
 
    class Symmetry_rotation
    {
    public:
        Symmetry_rotation() {};
        ~Symmetry_rotation() {};
 
        //--------------------------------------------------------------------------------
        // getters
        const std::map<Tap, std::set<TC>>& get_irreducible_sector()const { return this->irs_.get_irreducible_sector(); }
        TCdouble get_return_lattice(const Symmetry& symm,
            const ModuleBase::Matrix3& gmatd, const TCdouble gtransd,
            const TCdouble& posd_a1, const TCdouble& posd_a2)const
        {
            return this->irs_.get_return_lattice(symm, gmatd, gtransd, posd_a1, posd_a2);
        }
        TCdouble get_return_lattice(const int iat, const int isym) const
        {
            return this->irs_.get_return_lattice(iat, isym);
        }
        const std::vector<std::vector<RI::Tensor<std::complex<double>>>>& rotmat_Slm = this->rotmat_Slm_;
        const int& abfs_Lmax = this->abfs_Lmax_;
        //--------------------------------------------------------------------------------
        // setters
        void find_irreducible_sector(const Symmetry& symm, const Atom* atoms, const Statistics& st,
            const std::vector<TC>& Rs, const TC& period, const Lattice& lat)
        {
            this->irs_.find_irreducible_sector(symm, atoms, st, Rs, period, lat);
        }
        void set_abfs_Lmax(const int l) { this->abfs_Lmax_ = l; }
        void set_Cs_rotation(const std::vector<std::vector<int>>& abfs_l_nchi);
        //--------------------------------------------------------------------------------
 
        void cal_Ms(const K_Vectors& kv,
            //const std::vector<std::map<int, TCdouble>>& kstars,
            const UnitCell& ucell, const Parallel_2D& pv);
 
        std::vector<std::vector<std::complex<double>>>restore_dm(const K_Vectors& kv,
            const std::vector<std::vector<std::complex<double>>>& dm_k_ibz,
            const Parallel_2D& pv)const;
        std::vector<std::vector<double>>restore_dm(const K_Vectors& kv,
            const std::vector<std::vector<double>>& dm_k_ibz,
            const Parallel_2D& pv)const;
        std::vector<std::complex<double>> rot_matrix_ao(const std::vector<std::complex<double>>& DMkibz,
            const int ik_ibz, const int kstar_size, const int isym, const Parallel_2D& pv, const bool TRS_conj = false) const;
 
        double wigner_d(const double beta, const int l, const int m1, const int m2) const;
        std::complex<double> wigner_D(const TCdouble& euler_angle, const int l, const int m1, const int m2, const bool inv) const;
 
        std::complex<double> ovlp_Ylm_Slm(const int l, const int m1, const int m2) const;
 
        TCdouble get_euler_angle(const ModuleBase::Matrix3& gmatc) const;
 
        void cal_rotmat_Slm(const ModuleBase::Matrix3* gmatc, const int lmax);
 
        void set_block_to_mat2d(const int starti, const int startj, const RI::Tensor<std::complex<double>>& block,
            std::vector<std::complex<double>>& obj_mat, const Parallel_2D& pv, const bool trans = false) const;
        void set_block_to_mat2d(const int starti, const int startj, const RI::Tensor<std::complex<double>>& block,
            std::vector<double>& obj_mat, const Parallel_2D& pv, const bool trans = false) const;
 
        std::vector<std::complex<double>> contruct_2d_rot_mat_ao(const Symmetry& symm, const Atom* atoms, const Statistics& cell_st,
            const TCdouble& kvec_d_ibz, int isym, const Parallel_2D& pv) const;
 
        std::vector<std::vector<RI::Tensor<std::complex<double>>>>& get_rotmat_Slm() { return this->rotmat_Slm_; }
 
        //--------------------------------------------------------------------------------
        template<typename Tdata>    // RI::Tensor type
        std::map<int, std::map<std::pair<int, TC>, RI::Tensor<Tdata>>> restore_HR(
            const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode,
            const std::map<int, std::map<std::pair<int, TC>, RI::Tensor<Tdata>>>& HR_irreduceble)const;
        template<typename TR>   // HContainer type
        void restore_HR(
            const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode,
            const hamilt::HContainer<TR>& HR_irreduceble, hamilt::HContainer<TR>& HR_rotated)const;
 
        //--------------------------------------------------------------------------------
        template<typename Tdata>    // RI::Tensor type, using col-major implementation
        void test_HR_rotation(const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode,
            const std::map<int, std::map<std::pair<int, TC>, RI::Tensor<Tdata>>>& HR_full);
        template<typename Tdata>    // test the rotation of RI coefficients 
        void test_Cs_rotation(const Symmetry& symm, const Atom* atoms, const Statistics& st,
            const std::map<int, std::map<std::pair<int, TC>, RI::Tensor<Tdata>>>& Cs_full)const;
        template<typename TR>   // HContainer type, using row-major implementation
        void test_HR_rotation(const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode,
            const hamilt::HContainer<TR>& HR_full);
        template<typename Tdata>    // HContainer type
        void print_HR(const std::map<int, std::map<std::pair<int, TC>, RI::Tensor<Tdata>>>& HR, const std::string name, const double& threshold = 0.0);
        //--------------------------------------------------------------------------------
 
    private:
        //--------------------------------------------------------------------------------
        std::vector<TC> get_Rs_from_BvK(const K_Vectors& kv)const;
        std::vector<TC> get_Rs_from_adjacent_list(const UnitCell& ucell,
                                                  const Grid_Driver& gd,
                                                  const Parallel_Orbitals& pv) const;
        //--------------------------------------------------------------------------------
 
        template<typename Tdata>    // RI::Tensor type, blas
        RI::Tensor<Tdata> rotate_atompair_serial(const RI::Tensor<Tdata>& A, const int isym,
            const Atom& a1, const Atom& a2, const char mode, bool output = false)const;
        template<typename Tdata>    // pointer type, blas
        void rotate_atompair_serial(Tdata* TAT, const Tdata* A, const int& nw1, const int& nw2, const int isym,
            const Atom& a1, const Atom& a2, const char mode)const;
        template<typename TR>    // HContainer type, pblas
        void rotate_atompair_parallel(const TR* Alocal_in, const int isym, const Atom* atoms, const Statistics& st,
            const Tap& ap_in, const Tap& ap_out, const char mode, const Parallel_Orbitals& pv, TR* Alocal_out, const bool output = false)const;
 
        template<typename Tdata>
        RI::Tensor<Tdata> rotate_singleC_serial(const RI::Tensor<Tdata>& C, const int isym,
            const Atom& a1, const Atom& a2, const int& type1, bool output = false)const;
 
        template<typename Tdata>
        RI::Tensor<Tdata> set_rotation_matrix(const Atom& a, const int& isym)const;
        template<typename Tdata>
        RI::Tensor<Tdata> set_rotation_matrix_abf(const int& type, const int& isym)const;
        //--------------------------------------------------------------------------------
 
        int nsym_ = 1;
 
        double eps_ = 1e-6;
 
        bool TRS_first_ = true; //if R(k)=-k, firstly use TRS to restore D(k) from D(R(k)), i.e conjugate D(R(k)).
 
        bool reduce_Cs_ = false;
        int abfs_Lmax_ = 0;
 
        std::vector<std::vector<int>> abfs_l_nchi_;
 
        std::vector<std::vector<RI::Tensor<std::complex<double>>>> rotmat_Slm_;
        // [natom][nsym], phase factor corresponding to a certain kvec_d_ibz
        // std::vector<std::vector<std::complex<double>>> phase_factor_;
 
        std::vector<std::map<int, std::vector<std::complex<double>>>> Ms_;
 
        Irreducible_Sector irs_;
 
    };
 
    template<typename T>  std::string vec3_fmt(const T& x, const T& y, const T& z)
    {
        return  "(" + std::to_string(x) + " " + std::to_string(y) + " " + std::to_string(z) + ")";
    }
    template<typename T>  std::string vec3_fmt(const ModuleBase::Vector3<T>& v)
    {
        return vec3_fmt(v.x, v.y, v.z);
    }
    // output k stars and the rotation matrices of Bloch orbitals
    void print_symrot_info_k(const ModuleSymmetry::Symmetry_rotation& symrot,
        const K_Vectors& kv, const UnitCell& ucell);
    void print_symrot_info_R(const Symmetry_rotation& symrot, const Symmetry& symm,
        const int lmax_ao, const std::vector<TC>& Rs);
}
 
#include "symmetry_rotation_R.hpp"
#include "symmetry_rotation_R_hcontainer.hpp"