diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 27158e05069..444640004c0 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -588,6 +588,9 @@ OBJS_SYMMETRY=symm_other.o\ symm_rho.o\ little_group.o\ symmetry.o\ + irreducible_sector.o\ + irred_sec_bvk.o\ + symm_rotation_k.o\ OBJS_XC=xc_functional.o\ xc_functional_op.o\ diff --git a/source/source_cell/klist.cpp b/source/source_cell/klist.cpp index fd8080076ad..c0be2ebc946 100644 --- a/source/source_cell/klist.cpp +++ b/source/source_cell/klist.cpp @@ -421,8 +421,9 @@ void K_Vectors::update_use_ibz(const int& nkstot_ibz, ModuleBase::GlobalFunc::OUT(ofs_running, "nkstot now", nkstot); - // qianrui fix a bug 2021-7-13: size for the spin_mult=2 doubling in set_kup_and_kdw() - this->kvec_d.resize(this->nkstot * this->spin_mult); + // qianrui fix a bug 2021-7-13: shrink kvec_d to the (now smaller) ibz count; + // the spin_mult=2 doubling for set_kup_and_kdw() is reserved there, on demand. + this->kvec_d.resize(this->nkstot); for (int i = 0; i < this->nkstot; ++i) { @@ -445,6 +446,18 @@ void K_Vectors::set_kup_and_kdw(std::ofstream& ofs_running) { ModuleBase::TITLE("K_Vectors", "setup_kup_and_kdw"); + // grow the containers expand_spin_kpoints() is about to append the down-spin + // copy into (indices [nkstot, 2*nkstot)). Resize only these; NOT via renew(), + // which would also resize kvec_c_full -- that one must keep holding the + // original, un-doubled, un-symmetry-reduced full-BZ mesh for later consumers + // (e.g. Ewald_Vq) regardless of what nkstot has become by this point. + const int nkstot_spin = this->nkstot * this->spin_mult; + this->kvec_c.resize(nkstot_spin); + this->kvec_d.resize(nkstot_spin); + this->wk.resize(nkstot_spin); + this->ngk.resize(nkstot_spin); + this->isk.resize(nkstot_spin); + KListIO::expand_spin_kpoints(this->spin_mult, this->kvec_c, this->kvec_d, @@ -523,8 +536,9 @@ void K_Vectors::reduce_by_symmetry(const UnitCell& ucell, ibz2bz); const int nkstot_ibz = kvec_d_ibz.size(); -#ifdef __EXX // setup kstars according to the final (max-norm) kvec_d_ibz + // (used by both EXX and DFT+U's crystal-symmetry density-matrix restoration; + // no LibRI dependency, so this must not be gated behind __EXX) if (ModuleSymmetry::Symmetry::symm_flag == 1) { KListIO::build_kstars(this->kvec_d, @@ -535,7 +549,6 @@ void K_Vectors::reduce_by_symmetry(const UnitCell& ucell, [&symm](double a, double b) { return symm.equal(a, b); }, this->kstars); } -#endif // output in kpoints file skpt = KListIO::ibz_kpt_table(this->nkstot, this->kvec_d, this->ibz_index, kvec_d_ibz); @@ -665,12 +678,12 @@ void K_Vectors::mpi_k(std::ofstream& ofs_running, const int my_rank, const int m this->kvec_d, this->kvec_c_full); -#ifdef __EXX // bcast kstars (rank 0 holds the filled maps; other ranks rebuild them) + // (no LibRI dependency; needed by DFT+U's symmetry restoration too, so this + // must not be gated behind __EXX) if (ModuleSymmetry::Symmetry::symm_flag == 1) { KListIO::bcast_kstars(this->kstars, this->nkstot, my_rank); } -#endif } // END SUBROUTINE mpi_k #endif diff --git a/source/source_cell/module_symmetry/CMakeLists.txt b/source/source_cell/module_symmetry/CMakeLists.txt index 7cdd95d129b..5b2a03cd967 100644 --- a/source/source_cell/module_symmetry/CMakeLists.txt +++ b/source/source_cell/module_symmetry/CMakeLists.txt @@ -17,6 +17,9 @@ add_library( symmetry.cpp symm_rot_spin.cpp little_group.cpp + irreducible_sector.cpp + irred_sec_bvk.cpp + symm_rotation_k.cpp ) if(ENABLE_COVERAGE) diff --git a/source/source_lcao/module_ri/module_exx_symmetry/irred_sec_bvk.cpp b/source/source_cell/module_symmetry/irred_sec_bvk.cpp similarity index 99% rename from source/source_lcao/module_ri/module_exx_symmetry/irred_sec_bvk.cpp rename to source/source_cell/module_symmetry/irred_sec_bvk.cpp index b97f95a990b..5d506cab434 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/irred_sec_bvk.cpp +++ b/source/source_cell/module_symmetry/irred_sec_bvk.cpp @@ -75,7 +75,7 @@ namespace ModuleSymmetry return; } - // extern lattice to minimal BvK lattice, and set direct coordinates in min BvK lattice + // extern lattice to minimal BvK lattice, and set direct coordinates in min BvK lattice int bvk_gcd = gcd(bvk_period[0], gcd(bvk_period[1], bvk_period[2])); const TC bvk_min_period = TC({ bvk_period[0] / bvk_gcd, bvk_period[1] / bvk_gcd, bvk_period[2] / bvk_gcd }); const int bvk_nat = st.nat * bvk_min_period[0] * bvk_min_period[1] * bvk_min_period[2]; @@ -199,4 +199,4 @@ namespace ModuleSymmetry // return in_plain; // } -}; \ No newline at end of file +}; diff --git a/source/source_lcao/module_ri/module_exx_symmetry/irreducible_sector.cpp b/source/source_cell/module_symmetry/irreducible_sector.cpp similarity index 97% rename from source/source_lcao/module_ri/module_exx_symmetry/irreducible_sector.cpp rename to source/source_cell/module_symmetry/irreducible_sector.cpp index 26157460c3c..c6521c7ab52 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/irreducible_sector.cpp +++ b/source/source_cell/module_symmetry/irreducible_sector.cpp @@ -1,5 +1,4 @@ -#include "source_lcao/module_ri/module_exx_symmetry/irreducible_sector.h" -#include "source_io/module_parameter/parameter.h" +#include "source_cell/module_symmetry/irreducible_sector.h" namespace ModuleSymmetry { // Raw-index dispatch shared by the real-space sector helpers, matching the convention used @@ -62,7 +61,7 @@ namespace ModuleSymmetry } } - // Perfoming {R|t} to atom position r in the R=0 lattice, we get Rr+t, which may get out of R=0 lattice, + // Perfoming {R|t} to atom position r in the R=0 lattice, we get Rr+t, which may get out of R=0 lattice, // whose image in R=0 lattice is r'=Rr+t-O. This function is to get O for each atom and each symmetry operation. // the range of direct position is [-0.5, 0.5). TCdouble Irreducible_Sector::get_return_lattice(const Symmetry& symm, @@ -109,7 +108,7 @@ namespace ModuleSymmetry ModuleBase::TITLE("Symmetry_rotation", "cal_return_lattice_all"); // Columns [0, nrotk) are the unitary operations; columns [nrotk, nrotk+nrotk_anti) are the // spatial parts of the antiunitary elements Theta*g of the Shubnikov group (nspin=4 magnetic), - // so that Symmetry_rotation can address both with one raw index. + // so that Symmetry_rotation can address both with one raw index. this->return_lattice_.resize(st.nat, std::vector(symm.nrotk + symm.nrotk_anti)); for (int iat1 = 0;iat1 < st.nat;++iat1) { @@ -170,12 +169,12 @@ namespace ModuleSymmetry std::cout << std::endl; } } - void Irreducible_Sector::write_irreducible_sector() + void Irreducible_Sector::write_irreducible_sector(const std::string& output_dir) { - if(GlobalV::MY_RANK == 0) + if(GlobalV::MY_RANK == 0 && !output_dir.empty()) { std::ofstream ofs; - ofs.open(PARAM.globalv.global_out_dir + "irreducible_sector.txt"); + ofs.open(output_dir + "irreducible_sector.txt"); for (auto& irap_irR : this->irreducible_sector_) { for (auto& irR : irap_irR.second){ofs << "atompair (" << irap_irR.first.first << ", " << irap_irR.first.second << "), R = (" << irR[0] << ", " << irR[1] << ", " << irR[2] << ") \n";} @@ -184,7 +183,7 @@ namespace ModuleSymmetry } } - void Irreducible_Sector::find_irreducible_sector(const Symmetry& symm, const Atom* atoms, const Statistics& st, const std::vector& Rs, const TC& period, const Lattice& lat) + void Irreducible_Sector::find_irreducible_sector(const Symmetry& symm, const Atom* atoms, const Statistics& st, const std::vector& Rs, const TC& period, const Lattice& lat, const std::string& output_dir) { this->full_map_to_irreducible_sector_.clear(); this->irreducible_sector_.clear(); @@ -277,6 +276,6 @@ namespace ModuleSymmetry assert(total_apR_in_star == this->full_map_to_irreducible_sector_.size()); // this->output_full_map_to_irreducible_sector(st.nat); // this->output_sector_star(); - this->write_irreducible_sector(); + this->write_irreducible_sector(output_dir); } } diff --git a/source/source_lcao/module_ri/module_exx_symmetry/irreducible_sector.h b/source/source_cell/module_symmetry/irreducible_sector.h similarity index 89% rename from source/source_lcao/module_ri/module_exx_symmetry/irreducible_sector.h rename to source/source_cell/module_symmetry/irreducible_sector.h index 7bfdd1d3eae..2b091a7e42c 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/irreducible_sector.h +++ b/source/source_cell/module_symmetry/irreducible_sector.h @@ -2,7 +2,7 @@ #include #include #include -#include "source_lcao/module_ri/abfs_vector3_order.h" +#include "source_base/vector3.h" #include "source_base/matrix3.h" #include "source_cell/module_symmetry/symmetry.h" #include "source_cell/klist.h" @@ -12,7 +12,7 @@ namespace ModuleSymmetry using Tap = std::pair; using TC = std::array; using TapR = std::pair; - using TCdouble = Abfs::Vector3_Order; + using TCdouble = ModuleBase::Vector3; class Irreducible_Sector { @@ -54,13 +54,18 @@ namespace ModuleSymmetry //-------------------------------------------------------------------------------- /// The main function to find irreducible sector: {abR} + /// @param output_dir if non-empty, write_irreducible_sector() dumps a debug listing to + /// /irreducible_sector.txt (callers pass the run's global + /// output directory setting); kept as an explicit argument rather than + /// reading the global config directly so this LibRI-free class has no + /// module_parameter link dependency. void find_irreducible_sector(const Symmetry& symm, const Atom* atoms, const Statistics& st, - const std::vector& Rs, const TC& period, const Lattice& lat); + const std::vector& Rs, const TC& period, const Lattice& lat, const std::string& output_dir = ""); const std::map>& get_irreducible_sector()const { return this->irreducible_sector_; } // const std::map>> convirt_irreducible_sector() {}; //-------------------------------------------------------------------------------- - /// Perfoming {R|t} to atom position r in the R=0 lattice, we get Rr+t, which may get out of R=0 lattice, + /// Perfoming {R|t} to atom position r in the R=0 lattice, we get Rr+t, which may get out of R=0 lattice, /// whose image in R=0 lattice is r'=Rr+t-O. This function is to get O for each atom and each symmetry operation. /// the range of direct position is [-0.5, 0.5). TCdouble get_return_lattice(const Symmetry& symm, @@ -104,7 +109,7 @@ namespace ModuleSymmetry void output_full_map_to_irreducible_sector(const int nat); void output_sector_star(); - void write_irreducible_sector(); + void write_irreducible_sector(const std::string& output_dir); //-------------------------------------------------------------------------------- /// The sub functions judge special symmetry @@ -140,5 +145,6 @@ namespace ModuleSymmetry int bvk_nsym_; friend class Symmetry_rotation; + friend class Symmetry_rotation_k; }; -} \ No newline at end of file +} diff --git a/source/source_cell/module_symmetry/symm_rotation_k.cpp b/source/source_cell/module_symmetry/symm_rotation_k.cpp new file mode 100644 index 00000000000..3de1d1b47b2 --- /dev/null +++ b/source/source_cell/module_symmetry/symm_rotation_k.cpp @@ -0,0 +1,602 @@ +#include "symm_rotation_k.h" +#include "source_base/constants.h" +#include +#include "source_base/parallel_reduce.h" +#include "source_base/parallel_global.h" +#include "source_base/module_external/scalapack_connector.h" +#include "source_base/module_external/blas_connector.h" +#include "source_base/tool_title.h" +#include "source_base/timer.h" + +namespace ModuleSymmetry +{ + std::vector Symmetry_rotation_k::get_bvk_cells(const TC& period) + { + std::vector cells; + cells.reserve(static_cast(period[0]) * period[1] * period[2]); + for (int ix = 0; ix < period[0]; ++ix) { + for (int iy = 0; iy < period[1]; ++iy) { + for (int iz = 0; iz < period[2]; ++iz) { + cells.push_back({ix, iy, iz}); + } } } + return cells; + } + + void Symmetry_rotation_k::cal_Ms(const K_Vectors& kv, + const UnitCell& ucell, const Parallel_2D& pv, const int nspin) + { + ModuleBase::TITLE("Symmetry_rotation_k", "cal_Ms"); + ModuleBase::timer::start("Symmetry_rotation_k", "cal_Ms"); + + this->nspin_ = nspin; + this->nsym_ = ucell.symm.nrotk; + this->nanti_ = ucell.symm.nrotk_anti; + this->magnetic_nspin4_ = ucell.symm.magnetic_nspin4; + this->eps_ = ucell.symm.epsilon; + if (this->irs_.invmap_.empty()) + { + this->irs_.invmap_.resize(ucell.symm.nrotk); + ucell.symm.gmatrix_invmap(ucell.symm.gmatrix, ucell.symm.nrotk, this->irs_.invmap_.data()); + } + // 1. calculate the rotation matrix in real spherical harmonics representation for each symmetry operation: [T_l (isym)]_mm' + const int nop_tot = this->nsym_ + this->nanti_; + std::vector gmatc(nop_tot); + for (int i = 0;i < nsym_;++i) { gmatc[i] = this->irs_.direct_to_cartesian(ucell.symm.gmatrix[i], ucell.latvec); } + for (int j = 0;j < this->nanti_;++j) + { gmatc[nsym_ + j] = this->irs_.direct_to_cartesian(ucell.symm.gmatrix_anti[j], ucell.latvec); } + this->cal_rotmat_Slm(gmatc.data(), std::max(this->abfs_Lmax_, ucell.lmax), nop_tot); + + // 1.5 (nspin=4) the SU(2) spin-1/2 rotation U(isym) for each symmetry operation. The AO + // rotation matrix M becomes the spinor operator T(isym) (x) U(isym) so that the same + // gemm D(k)=M^dagger D(k_ibz) M rotates both the orbital and the spin part at once. + // For an antiunitary element Theta*g only the spatial part g enters M here; the Theta + // (sigma_y (.)^* sigma_y) is applied afterwards in restore_dm. + std::vector spin_U(nop_tot, SpinRotation::Su2{ 1.0, 0.0, 0.0, 1.0 }); + if (this->nspin_ == 4) + { + for (int i = 0;i < nop_tot;++i) { spin_U[i] = SpinRotation::so3_to_su2(gmatc[i]); } + } + this->spin_U_ = spin_U; // keep for restore_HR_nspin4 (real-space EXX H(R) spin mixing) + + // 2. calculate the rotation matrix in AO-representation for each ibz_kpoint and symmetry operation: M(k, isym) + int nks_ibz = kv.kstars.size(); // kv.nks = 2 * kv.nks_ibz when nspin=2 + this->Ms_.assign(nks_ibz, {}); + this->little_groups_.assign(nks_ibz, {}); + + // (k-point pools, KPAR>1) kv.kvec_d only holds the k-points owned by this pool, so + // kv.kvec_d[ik_ibz] is only valid for ik_ibz < kv.para_k.nks_np and is otherwise either + // out of range or (after a caller's spin-doubling resize) a meaningless zero placeholder. + // Gather the (small, size nks_ibz) global ibz-representative k-vector list once so every + // pool builds the correct rotation matrix for every ibz k, not just the ones it owns. +#ifdef __MPI + // inlined equivalent of Parallel_Kpoints::gatherkvec (avoided as a direct call so this + // class does not pull in a link dependency on parallel_kpoints.cpp for every target that + // links the "symmetry" library): every rank in the owning pool holds the same local + // k-vectors, so only the pool root contributes to the MPI_Allreduce, matching gatherkvec. + int world_rank = 0; + MPI_Comm_rank(MPI_COMM_WORLD, &world_rank); + const bool is_pool_root = (world_rank == kv.para_k.get_startpro_pool(kv.para_k.my_pool)); + std::vector> kvec_d_ibz_global(nks_ibz, ModuleBase::Vector3(0.0, 0.0, 0.0)); + for (int i = 0; i < kv.para_k.nks_np; ++i) + { + if (is_pool_root) { kvec_d_ibz_global[i + kv.para_k.startk_pool[kv.para_k.my_pool]] = kv.kvec_d[i]; } + } + MPI_Allreduce(MPI_IN_PLACE, kvec_d_ibz_global.data(), 3 * nks_ibz, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); +#else + const std::vector>& kvec_d_ibz_global = kv.kvec_d; +#endif + + // A k-star contains only one operation per distinct k point. The other + // operations fixing k (modulo a reciprocal lattice vector) must still + // be averaged: a finite-grid SCF density need not respect this little group. + for (int ik_ibz = 0; ik_ibz < nks_ibz; ++ik_ibz) + { + std::set needed; + for (const std::pair>& member : kv.kstars[ik_ibz]) + { + const int op = (!this->magnetic_nspin4_ && member.first >= nsym_) + ? member.first - nsym_ : member.first; + needed.insert(op); + } + for (int op = 0; op < nsym_; ++op) + { + const ModuleBase::Vector3 delta = kvec_d_ibz_global[ik_ibz] * ucell.symm.kgmatrix[op] - kvec_d_ibz_global[ik_ibz]; + if (std::abs(delta.x - std::round(delta.x)) < this->eps_ + && std::abs(delta.y - std::round(delta.y)) < this->eps_ + && std::abs(delta.z - std::round(delta.z)) < this->eps_) + { + this->little_groups_[ik_ibz].push_back(op); + needed.insert(op); + } + } + for (const int op : needed) + { + this->Ms_[ik_ibz][op] = this->contruct_2d_rot_mat_ao( + ucell.symm, ucell.atoms, ucell.st, kvec_d_ibz_global[ik_ibz], op, pv, spin_U[op]); + } + } + + ModuleBase::timer::end("Symmetry_rotation_k", "cal_Ms"); + } + + std::vector>> Symmetry_rotation_k::restore_dm(const K_Vectors& kv, + const std::vector>>& dm_k_ibz, const Parallel_2D& pv)const + { + ModuleBase::TITLE("Symmetry_rotation_k", "restore_dm"); + ModuleBase::timer::start("Symmetry_rotation_k", "restore_dm"); + std::vector>> dm_k_full; + int nspin0 = this->nspin_ == 2 ? 2 : 1; + // (k-point pools, KPAR>1) dm_k_ibz (elecstate::DensityMatrix::_DMK) only ever holds + // the irreducible k-points owned by THIS pool (_nk = kv.get_nks()/nspin, see + // setup_dm.cpp), never the global set -- so nk here must be the local count, and + // kv.kstars (which is global, identical on every pool) must be indexed via the + // local-to-global map kv.ik2iktot, not via the local loop variable directly. + // This is safe: D(k) -> D(R) (or, for DFT+U, the occupation matrix built from it) + // is a linear sum over k, so each pool returning only the stars of its own local + // irreducible k-points, to be combined by the caller's existing cross-pool + // reduction (e.g. compute_occ_from_dmr's Parallel_Reduce::reduce_all), gives the + // exact same total as if every pool held the full global k-set -- no pool needs + // (or has to pay for gathering) the complete global D(k) at any point. + int nk = kv.get_nks() / nspin0; + const int nks_ibz_global = static_cast(this->little_groups_.size()); + + // (nspin=4) Sigma_y = I (x) sigma_y for the time-reversal spin flip; k-independent, build once. + std::vector> sigma_y; + if (this->nspin_ == 4) { sigma_y = this->set_sigma_y_2d(pv); } + + for (int is = 0;is < nspin0;++is) + { + for (int ik_local = 0;ik_local < nk;++ik_local) + { + const int ik_ibz = kv.ik2iktot[ik_local + is * nk] % nks_ibz_global; + // P_k D = |G_k|^{-1} sum_g M_g^T D M_g^*. This preserves + // Hermiticity and makes restoration independent of the chosen + // star representative; rotating just one arbitrary D does not. + const std::vector& little_group = this->little_groups_.at(ik_ibz); + assert(!little_group.empty()); + std::vector> projected = dm_k_ibz[ik_local + is * nk]; + if (little_group.size() > 1) + { + std::fill(projected.begin(), projected.end(), 0.0); + for (const int op : little_group) + { + const std::vector> rotated = this->rot_matrix_ao( + dm_k_ibz[ik_local + is * nk], ik_ibz, little_group.size(), op, pv); + for (size_t i = 0; i < projected.size(); ++i) + { + projected[i] += rotated[i]; + } + } + } + for (const std::pair>& isym_kvd : kv.kstars[ik_ibz]) + { + if (isym_kvd.first == 0) + { + double factor = 1.0 / static_cast(kv.kstars[ik_ibz].size()); + std::vector> dm_scaled(pv.get_local_size()); + for (int i = 0;i < pv.get_local_size();++i) { dm_scaled[i] = factor * projected[i]; } + dm_k_full.push_back(dm_scaled); + } + else if (isym_kvd.first < nsym_) + { //space group operations + dm_k_full.push_back(this->rot_matrix_ao(projected, ik_ibz, kv.kstars[ik_ibz].size(), isym_kvd.first, pv)); + } + else + { // antiunitary elements: Theta * (spatial operation) + // D(Theta*g k_ibz) = sigma_y [D(g k_ibz)]^* sigma_y with D(g k_ibz) = M^dagger D M. + // For nspin=4, first do the (non-conjugated) spatial rotation, then the spin flip; + // for nspin<4 (Theta=K) the original TRS_conj path already gives the conjugate. + // + // Which spatial operation the index denotes depends on the regime, matching + // how the k-reduction filled kgmatrix[] (see K_Vectors::reduce_by_symmetry): + // - nspin=4 magnetic (Shubnikov): index j+nsym_ is the antiunitary element + // Theta*gmatrix_anti[j]; its Ms is stored under the RAW key j+nsym_. + // - otherwise (grey group / nspin<4): index i+nsym_ is Theta*gmatrix[i], + // i.e. the unitary operation i, whose Ms is stored under key i. + const int isym_M = this->magnetic_nspin4_ ? isym_kvd.first : (isym_kvd.first - nsym_); + if (this->nspin_ == 4) + { + // m=0: gray group: the space-group part of anti-unitary elements are the same of the unitary elements, isym_M < nsym_ + // m!=0: Shubnikov group: using different space-group part of anti-unitary elements stored in gmatrix_anti with isym_M >= nsym_ + dm_k_full.push_back(this->trs_spin_rotate( + this->rot_matrix_ao(projected, ik_ibz, kv.kstars[ik_ibz].size(), isym_M, pv, false), + sigma_y, pv, 1.0)); + } + else + { + dm_k_full.push_back(this->rot_matrix_ao(projected, ik_ibz, kv.kstars[ik_ibz].size(), isym_M, pv, true)); + } + } + } + } + } + ModuleBase::timer::end("Symmetry_rotation_k", "restore_dm"); + return dm_k_full; + } + std::vector> Symmetry_rotation_k::restore_dm(const K_Vectors& kv, + const std::vector>& dm_k_ibz, const Parallel_2D& pv)const + { + return dm_k_ibz;// do nothing for gamma_only + } + + // calculate Wigner D matrix + double Symmetry_rotation_k::wigner_d(const double beta, const int l, const int m1, const int m2) const + { + auto factorial = [](int n) -> int { + int result = 1; + for (int i = 1;i <= n;++i) { result *= i; +} + return result; + }; + double result = 0.0; + for (int i = std::max(0, m2 - m1);i <= std::min(l - m1, l + m2);++i) { + result += std::pow(-1, i) * std::sqrt(factorial(l + m1) * factorial(l - m1) * factorial(l + m2) * factorial(l - m2)) + * std::pow(std::cos(beta / 2), 2 * l + m2 - m1 - 2 * i) * std::pow(-std::sin(beta / 2), m1 - m2 + 2 * i) + / (factorial(i) * factorial(l - m1 - i) * factorial(l + m2 - i) * factorial(i - m2 + m1)); +} + return result; + } + + std::complex Symmetry_rotation_k::wigner_D(const TCdouble& euler_angle, const int l, const int m1, const int m2, const bool inv) const + { + std::complex prefac(inv ? std::pow(-1, l) : 1, 0); + return std::exp(-ModuleBase::IMAG_UNIT * static_cast(m1) * euler_angle.x) + * std::exp(-ModuleBase::IMAG_UNIT * static_cast(m2) * euler_angle.z) + * wigner_d(euler_angle.y, l, m1, m2) * prefac; + } + + // c^l_{m1, m2}= + std::complex Symmetry_rotation_k::ovlp_Ylm_Slm(const int l, const int m1, const int m2) const + { + if (m1 == m2) + { + if (m1 == 0) { return 1.0; +} + if (m1 > 0) { return 1 / std::sqrt(2); +} + if (m1 < 0) { return std::pow(-1, m1) * ModuleBase::IMAG_UNIT / std::sqrt(2); +} + } + else if (m1 == -m2) + { + if (m1 > 0) { return -ModuleBase::IMAG_UNIT / std::sqrt(2); +} + if (m1 < 0) { return std::pow(-1, m1) / std::sqrt(2); +} + } + return 0.0; + } + + // reference: https://github.com/minyez/abf_trans/blob/f9e68e68069a94610d89e077bfe6e8ffac0b097d/src/rotate.cpp#L118 + // because the atom position here is row vector, the original gmatrix(eular angle) is transposed. + // gmatc: the rotation matrix under the basis of cartesian coordinates + // gmatc should be a rotation matrix, i.e. det(gmatc)=1 + TCdouble Symmetry_rotation_k::get_euler_angle(const ModuleBase::Matrix3& gmatc) const + { + double threshold = this->eps_; + double alpha = 0.0, beta = 0.0, gamma = 0.0; + if (std::fabs(gmatc.e32) > threshold || std::fabs(gmatc.e31) > threshold) // sin(beta) is not zero + { + // use the 2-angle elements to get alpha and gamma + alpha = std::atan2(gmatc.e32, gmatc.e31); + if (alpha < 0) { alpha += 2 * ModuleBase::PI; +} + gamma = std::atan2(gmatc.e23, -gmatc.e13); + if (gamma < 0) { gamma += 2 * ModuleBase::PI; +} + // use the larger one of 2-angle elements to calculate beta + if (std::fabs(gmatc.e32) > std::fabs(gmatc.e31)) { + beta = std::atan2(gmatc.e32 / std::sin(alpha), gmatc.e33); + } else { + beta = std::atan2(gmatc.e31 / std::cos(alpha), gmatc.e33); +} + } + else + {//sin(beta)=0, beta = 0 or pi, only (alpha+gamma) or (alpha-gamma) is important. now assign this to alpha. + alpha = std::atan2(gmatc.e12, gmatc.e11); + if (alpha < 0) { alpha += 2 * ModuleBase::PI; +} + // if beta=0, gmatc.e11=cos(alpha+gamma), gmatc.e21=sin(alpha+gamma) + // if beta=pi, gmatc.e11=cos(pi+alpha-gamma), gmatc.e21=sin(pi+alpha-gamma) + if (gmatc.e33 > 0) + { + beta = 0; + gamma = 0; //alpha+gamma=alpha => gamma=0 + } + else + { + beta = ModuleBase::PI; + gamma = ModuleBase::PI;// pi+alpha-gamma=alpha => gamma=pi + } + } + return TCdouble(alpha, beta, gamma); + } + + // in: the real value of m in range {-l, -l+1, ..., 0, ..., l-1, l} + // out: the index of the orbital in a fixed {n， l}, i.e. the index in array [0, 1, -1, 2, -2, ...] + inline int m2im_k(int m) + { + return (m > 0 ? 2 * m - 1 : -2 * m); + } + + /// T_mm' = [c^\dagger D c]_mm' + void Symmetry_rotation_k::cal_rotmat_Slm(const ModuleBase::Matrix3* gmatc, const int lmax, const int nop) + { + ++this->rotmat_Slm_version_; + const int nop_tot = (nop < 0) ? this->nsym_ : nop; + this->rotmat_Slm_.resize(nop_tot); + // c matrix is independent on isym + std::vector c_mm(lmax + 1); + for (int l = 0;l <= lmax;++l) { + c_mm[l].create(2 * l + 1, 2 * l + 1); +} + for (int l = 0;l <= lmax;++l) { + for (int m1 = -l;m1 <= l;++m1) { + for (int m2 = -l;m2 <= l;++m2) { + c_mm[l](m2im_k(m1), m2im_k(m2)) = ovlp_Ylm_Slm(l, m1, m2); +} +} +} + + for (int isym = 0;isym < nop_tot;++isym) + { + // if R is a reflection operation, calculate D^l(R)=(-1)^l*D^l(IR), so the euler angle of (IR) is needed. + TCdouble euler_angle = get_euler_angle(gmatc[isym].Det() > 0 ? + gmatc[isym] : gmatc[isym] * ModuleBase::Matrix3(-1, 0, 0, 0, -1, 0, 0, 0, -1)); + + this->rotmat_Slm_[isym].resize(lmax + 1); + for (int l = 0;l <= lmax;++l) + {// wigner D matrix + ModuleBase::ComplexMatrix D_mm(2 * l + 1, 2 * l + 1); + for (int m1 = -l;m1 <= l;++m1) { + for (int m2 = -l;m2 <= l;++m2) { + D_mm(m2im_k(m1), m2im_k(m2)) = wigner_D(euler_angle, l, m1, m2, (gmatc[isym].Det() < 0)); +} +} + this->rotmat_Slm_[isym][l] = transpose(c_mm[l], true) * D_mm * c_mm[l]; + } + } + } + + void Symmetry_rotation_k::set_block_to_mat2d(const int starti, const int startj, const ModuleBase::ComplexMatrix& block, + std::vector>& obj_mat, const Parallel_2D& pv, const bool trans) const + { // caution: ComplaxMatrix is row-major(col-continuous), but obj_mat is col-major(row-continuous) + for (int j = 0;j < block.nr;++j) {//outside dimension + for (int i = 0;i < block.nc;++i) { //inside dimension + if (pv.in_this_processor(starti + i, startj + j)) + { + int index = pv.global2local_col(startj + j) * pv.get_row_size() + pv.global2local_row(starti + i); + obj_mat[index] = trans ? block(i, j) : block(j, i); + } +} +} + } + + void Symmetry_rotation_k::set_block_to_mat2d(const int starti, const int startj, const ModuleBase::ComplexMatrix& block, + std::vector& obj_mat, const Parallel_2D& pv, const bool trans) const + { // caution: ComplaxMatrix is row-major(col-continuous), but obj_mat is col-major(row-continuous) + for (int j = 0;j < block.nr;++j) {//outside dimension + for (int i = 0;i < block.nc;++i) { //inside dimension + if (pv.in_this_processor(starti + i, startj + j)) + { + int index = pv.global2local_col(startj + j) * pv.get_row_size() + pv.global2local_row(starti + i); + obj_mat[index] = trans ? block(i, j).real() : block(j, i).real(); + } +} +} + } + + // 2d-block parallized rotation matrix in AO-representation, denoted as M. + // finally we will use D(k)=M(R, k)^\dagger*D(Rk)*M(R, k) to D(k) from D(Rk) in cal_Ms. + std::vector> Symmetry_rotation_k::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 SpinRotation::Su2& spin_U) const + { + const bool soc = (this->nspin_ == 4); + const int npol = soc ? 2 : 1; // spinor: global AO index is spin-fast interleaved, I = npol*iw_orb + s + std::vector> M_isym(pv.get_local_size(), 0.0); + // isym >= symm.nrotk addresses the antiunitary coset (spatial part gmatrix_anti[isym-nrotk]), + // whose atom map lives in a separate table. + const int nrotk_u = symm.nrotk; + auto rotated_atom = [&symm, nrotk_u](const int is, const int iat) -> int + { + return (is < nrotk_u) ? symm.get_rotated_atom(is, iat) + : symm.get_rotated_atom_anti(is - nrotk_u, iat); + }; + for (int iat1 = 0;iat1 < cell_st.nat;++iat1) + { + int it = cell_st.iat2it[iat1]; // it1=it2 + int ia1 = cell_st.iat2ia[iat1]; + int iat2 = rotated_atom(isym, iat1); //iat2=rot(iat1) + int ia2 = cell_st.iat2ia[iat2]; + // cal phase factor from return lattice: exp(-ik_ibz*O) + double arg = -2 * ModuleBase::PI * kvec_d_ibz * this->irs_.return_lattice_[iat1][isym]; + std::complexphase_factor = std::complex(std::cos(arg), std::sin(arg)); + int iw1start = atoms[it].stapos_wf + ia1 * atoms[it].nw; + int iw2start = atoms[it].stapos_wf + ia2 * atoms[it].nw; + int iw = 0; + while (iw < atoms[it].nw) + { + int l = atoms[it].iw2l[iw]; + int nm = 2 * l + 1; + //caution: the order of m in orbitals may be different from increasing + if (!soc) + { + set_block_to_mat2d(iw2start + iw, iw1start + iw, + phase_factor * this->rotmat_Slm_[isym][l], M_isym, pv, true); + } + else + { + // M = T(isym) (x) U(isym): scatter phase * T_l(m,m') * U(a,b) to the interleaved + // spinor positions (row = rotated atom/spin, col = original atom/spin). For nspin=4 + // stapos_wf already carries the npol factor, so the per-atom offset is ia*nw*npol + // and the within-atom spinor index is (iw_orb)*npol + spin (spin is the fast index). + const int base2 = atoms[it].stapos_wf + ia2 * atoms[it].nw * npol; + const int base1 = atoms[it].stapos_wf + ia1 * atoms[it].nw * npol; + const ModuleBase::ComplexMatrix& Tl = this->rotmat_Slm_[isym][l]; + for (int m = 0;m < nm;++m) + { + for (int mp = 0;mp < nm;++mp) + { + const std::complex t = phase_factor * Tl(m, mp); + for (int a = 0;a < npol;++a) + { + for (int b = 0;b < npol;++b) + { + const int gi = base2 + (iw + m) * npol + a; + const int gj = base1 + (iw + mp) * npol + b; + if (pv.in_this_processor(gi, gj)) + { + const int index = pv.global2local_col(gj) * pv.get_row_size() + pv.global2local_row(gi); + // M(isym) = T_l (x) U is the spinor rep, with U = so3_to_su2 placed as-is: + // M[(m,a),(m',b)] = phase * T_l(m,m') * U_{ab}, U_{ab} = spin_U[a*npol + b]. + // Both T_l (rotmat_Slm) and U are ANTI-homomorphisms here (row-vector / R^T convention: + // rotmat_Slm(g)=R_orb(g)^{-1}, so3_to_su2 likewise), so this M is a consistent rep + // and rot_matrix_ao's stored-DM rotation M^T D M^* is exact for ALL ops. + M_isym[index] = t * spin_U[a * npol + b]; + } + } + } + } + } + } + iw += nm; + } + } + return M_isym; + } + + // D(k) = M^T(R, k) D(k_ibz) M^*(R, k), if D(k) is col-maj + // D^T(k) = M^\dagger(R, k) D^T(k_ibz) M(R, k), if D(k) is row-maj + // Ds from RI_2D_Comm are row-maj + // the link ik_ibz-isym-ik can be found in kstars. + std::vector> Symmetry_rotation_k::rot_matrix_ao(const std::vector>& DMkibz, + const int ik_ibz, const int kstar_size, const int isym, const Parallel_2D& pv, const bool TRS_conj) const + { + std::vector> DMk(pv.nloc, 0.0); + std::vector> DMkibz_M(pv.nloc, 0.0); // intermediate result + const char dagger = 'C'; + const char transpose = 'T'; + const char notrans = 'N'; + std::complex alpha(1.0, 0.0); + const std::complex beta(0.0, 0.0); + const int nbasis = pv.get_global_row_size(); + const int i1 = 1; + if (TRS_conj) + { + // D^T* = M^T [M^T (D^T)^T]^\dagger +#ifdef __MPI + ScalapackConnector::gemm(transpose, transpose, nbasis, nbasis, nbasis, + alpha, this->Ms_[ik_ibz].at(isym).data(), i1, i1, pv.desc, DMkibz.data(), i1, i1, pv.desc, + beta, DMkibz_M.data(), i1, i1, pv.desc); +#else + // without MPI, pv holds the whole (non-block-cyclic) dense matrix locally, + // so the 2D-block-cyclic pdgemm/pzgemm degenerates to a plain col-major gemm. + BlasConnector::gemm_cm(transpose, transpose, nbasis, nbasis, nbasis, + alpha, this->Ms_[ik_ibz].at(isym).data(), nbasis, DMkibz.data(), nbasis, + beta, DMkibz_M.data(), nbasis); +#endif + alpha.real(1.0 / static_cast(kstar_size)); +#ifdef __MPI + ScalapackConnector::gemm(transpose, dagger, nbasis, nbasis, nbasis, + alpha, this->Ms_[ik_ibz].at(isym).data(), i1, i1, pv.desc, DMkibz_M.data(), i1, i1, pv.desc, + beta, DMk.data(), i1, i1, pv.desc); +#else + BlasConnector::gemm_cm(transpose, dagger, nbasis, nbasis, nbasis, + alpha, this->Ms_[ik_ibz].at(isym).data(), nbasis, DMkibz_M.data(), nbasis, + beta, DMk.data(), nbasis); +#endif + } + else + { + // Physical DM rotation D(k) = M^dagger D(k_ibz) M, with M = T (x) U is the anti-homomorphism rep in row-major convention. + // ABACUS stores the DM transposed (S = D^T), for which this becomes S(gk) = M^T S(k_ibz) M^* = (conj M)^dagger S (conj M) + // For nspin<4 the orbital-only M is real, so Mc = M and this is bit-identical to the old M^dagger D M. + const std::vector>& Mref = this->Ms_[ik_ibz].at(isym); + std::vector> Mc(Mref.size()); + for (size_t i = 0; i < Mref.size(); ++i) { Mc[i] = std::conj(Mref[i]); } +#ifdef __MPI + ScalapackConnector::gemm(dagger, notrans, nbasis, nbasis, nbasis, + alpha, Mc.data(), i1, i1, pv.desc, DMkibz.data(), i1, i1, pv.desc, + beta, DMkibz_M.data(), i1, i1, pv.desc); +#else + BlasConnector::gemm_cm(dagger, notrans, nbasis, nbasis, nbasis, + alpha, Mc.data(), nbasis, DMkibz.data(), nbasis, + beta, DMkibz_M.data(), nbasis); +#endif + alpha.real(1.0 / static_cast(kstar_size)); +#ifdef __MPI + ScalapackConnector::gemm(notrans, notrans, nbasis, nbasis, nbasis, + alpha, DMkibz_M.data(), i1, i1, pv.desc, Mc.data(), i1, i1, pv.desc, + beta, DMk.data(), i1, i1, pv.desc); +#else + BlasConnector::gemm_cm(notrans, notrans, nbasis, nbasis, nbasis, + alpha, DMkibz_M.data(), nbasis, Mc.data(), nbasis, + beta, DMk.data(), nbasis); +#endif + } + return DMk; + } + + std::vector> Symmetry_rotation_k::set_sigma_y_2d(const Parallel_2D& pv) const + { + std::vector> sigma_y(pv.get_local_size(), 0.0); + const int nlocal = pv.get_global_row_size(); // = 2*nao for nspin=4 + // sigma_y = [[0, -i], [i, 0]] on the interleaved spin index (I = 2*iorb + spin) + const std::complex sy[2][2] = { {std::complex(0.0, 0.0), std::complex(0.0, -1.0)}, + {std::complex(0.0, 1.0), std::complex(0.0, 0.0)} }; + for (int iorb = 0; 2 * iorb < nlocal; ++iorb) + { + for (int a = 0; a < 2; ++a) + { + const int b = 1 - a; // only the off-diagonal spin entries are non-zero + const int gi = 2 * iorb + a; + const int gj = 2 * iorb + b; + if (pv.in_this_processor(gi, gj)) + { + const int index = pv.global2local_col(gj) * pv.get_row_size() + pv.global2local_row(gi); + sigma_y[index] = sy[a][b]; + } + } + } + return sigma_y; + } + + std::vector> Symmetry_rotation_k::trs_spin_rotate(const std::vector>& X, + const std::vector>& sigma_y, const Parallel_2D& pv, const double scale) const + { + // stored (transposed 2d-block) form of D_new = sigma_y * conj(D) * sigma_y is + // Sigma_y * conj(X) * Sigma_y (Sigma_y^T = -Sigma_y, the two minus signs cancel). + const char notrans = 'N'; + const int nbasis = pv.get_global_row_size(); + const int i1 = 1; + const std::complex one(1.0, 0.0); + const std::complex beta(0.0, 0.0); + std::vector> Xc(X.size()); + for (size_t i = 0; i < X.size(); ++i) { Xc[i] = std::conj(X[i]); } + std::vector> tmp(pv.get_local_size(), 0.0); + std::vector> out(pv.get_local_size(), 0.0); + // tmp = Sigma_y * conj(X) +#ifdef __MPI + ScalapackConnector::gemm(notrans, notrans, nbasis, nbasis, nbasis, + one, sigma_y.data(), i1, i1, pv.desc, Xc.data(), i1, i1, pv.desc, + beta, tmp.data(), i1, i1, pv.desc); +#else + // without MPI, pv holds the whole (non-block-cyclic) dense matrix locally, + // so the 2D-block-cyclic pzgemm degenerates to a plain col-major gemm. + BlasConnector::gemm_cm(notrans, notrans, nbasis, nbasis, nbasis, + one, sigma_y.data(), nbasis, Xc.data(), nbasis, + beta, tmp.data(), nbasis); +#endif + // out = scale * tmp * Sigma_y +#ifdef __MPI + ScalapackConnector::gemm(notrans, notrans, nbasis, nbasis, nbasis, + std::complex(scale, 0.0), tmp.data(), i1, i1, pv.desc, sigma_y.data(), i1, i1, pv.desc, + beta, out.data(), i1, i1, pv.desc); +#else + BlasConnector::gemm_cm(notrans, notrans, nbasis, nbasis, nbasis, + std::complex(scale, 0.0), tmp.data(), nbasis, sigma_y.data(), nbasis, + beta, out.data(), nbasis); +#endif + return out; + } +} diff --git a/source/source_cell/module_symmetry/symm_rotation_k.h b/source/source_cell/module_symmetry/symm_rotation_k.h new file mode 100644 index 00000000000..48ee839f4ff --- /dev/null +++ b/source/source_cell/module_symmetry/symm_rotation_k.h @@ -0,0 +1,177 @@ +#ifndef SYMM_ROTATION_K_H +#define SYMM_ROTATION_K_H +#include "irreducible_sector.h" +#include "source_base/parallel_2d.h" +#include "source_base/complexmatrix.h" +#include "source_cell/module_symmetry/symm_rot_spin.h" + +namespace ModuleSymmetry +{ + /// @brief k-space AO-representation symmetry restoration: reconstructs D(k) at every + /// k-star member from D(k_ibz), for crystal-symmetry-reduced BZ sampling. + /// + /// This is the LibRI-independent subset of what used to be a single + /// source_lcao/module_ri/module_exx_symmetry/Symmetry_rotation class: everything needed + /// to go from an irreducible-k-point density matrix to the full-BZ one, so callers that + /// only need that (e.g. DFT+U's occupation-matrix/DMR restoration) do not have to depend + /// on module_ri or LibRI. EXX/RPA's own real-space H(R)/RI-coefficient restoration (which + /// does need RI::Tensor) is built on top of this class in + /// source_lcao/module_ri/module_exx_symmetry/symm_rotation.h (ModuleSymmetry::Symmetry_rotation, + /// which inherits from this one). + class Symmetry_rotation_k + { + public: + Symmetry_rotation_k() {}; + virtual ~Symmetry_rotation_k() {}; + + //-------------------------------------------------------------------------------- + // getters + const std::map>& 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); + } + /// the rotation matrix under the basis of S_l^m. size: [nsym][lmax][nm*nm] + const std::vector>& 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& Rs, const TC& period, const Lattice& lat, const std::string& output_dir = "") + { + this->irs_.find_irreducible_sector(symm, atoms, st, Rs, period, lat, output_dir); + } + void set_abfs_Lmax(const int l) { this->abfs_Lmax_ = l; } + //-------------------------------------------------------------------------------- + /// functions to contruct rotation matrix in AO-representation + + /// The top-level calculation interface of this class. calculate the rotation matrix in AO representation: M + /// only need once call in each ion step (decided by the configuration) + /// @param kstars equal k points to each ibz-kpont, corresponding to a certain symmetry operations. + /// @param nspin stored as a member so restore_dm()/contruct_2d_rot_mat_ao() do not each + /// need to read the global nspin config setting (keeps this LibRI-free class + /// free of a module_parameter link dependency; every existing caller already + /// has nspin in scope). + void cal_Ms(const K_Vectors& kv, + const UnitCell& ucell, const Parallel_2D& pv, const int nspin); + + /// Use calculated M matrix to recover D(k) from D(k_ibz): D(k) = M(R, k)^\dagger D(k_ibz) M(R, k) + /// the link "ik_ibz-isym-ik" can be found in kstars: k_bz = gmat[isym](k) + std::vector>>restore_dm(const K_Vectors& kv, + const std::vector>>& dm_k_ibz, + const Parallel_2D& pv)const; + std::vector>restore_dm(const K_Vectors& kv, + const std::vector>& dm_k_ibz, + const Parallel_2D& pv)const; + std::vector> rot_matrix_ao(const std::vector>& DMkibz, + const int ik_ibz, const int kstar_size, const int isym, const Parallel_2D& pv, const bool TRS_conj = false) const; + + /// (nspin=4) build the 2*nao spin operator Sigma_y = I_nao (x) sigma_y in 2d-block layout. + std::vector> set_sigma_y_2d(const Parallel_2D& pv) const; + + /// (nspin=4) time-reversal on the spin density matrix: D(k) = sigma_y D^*(-k) sigma_y, + /// realized distribution-safely as scale * Sigma_y * conj(X) * Sigma_y (X is the already + /// space-group-rotated D(-k) stored in the transposed 2d-block convention). + std::vector> trs_spin_rotate(const std::vector>& X, + const std::vector>& sigma_y, const Parallel_2D& pv, const double scale) const; + + /// calculate Wigner D matrix + double wigner_d(const double beta, const int l, const int m1, const int m2) const; + std::complex wigner_D(const TCdouble& euler_angle, const int l, const int m1, const int m2, const bool inv) const; + + /// c^l_{m1, m2}= + std::complex ovlp_Ylm_Slm(const int l, const int m1, const int m2) const; + + /// calculate euler angle from rotation matrix + TCdouble get_euler_angle(const ModuleBase::Matrix3& gmatc) const; + + /// T_mm' = [c^\dagger D c]_mm', the rotation matrix in the representation of real sphere harmonics + /// @param nop number of operations in gmatc; <0 means nsym_ (the unitary ones only). + /// Pass nsym_+nanti_ to also build the antiunitary operations' T_l. + void cal_rotmat_Slm(const ModuleBase::Matrix3* gmatc, const int lmax, const int nop); + + /// set a block matrix onto a 2d-parallelized matrix(col-maj), at the position (starti, startj) + /// if trans=true, the block matrix is transposed before setting + void set_block_to_mat2d(const int starti, const int startj, const ModuleBase::ComplexMatrix& block, + std::vector>& obj_mat, const Parallel_2D& pv, const bool trans = false) const; + void set_block_to_mat2d(const int starti, const int startj, const ModuleBase::ComplexMatrix& block, + std::vector& obj_mat, const Parallel_2D& pv, const bool trans = false) const; + + /// 2d-block parallized rotation matrix in AO-representation, denoted as M. + /// finally we will use D(k)=M(R, k)^\dagger*D(Rk)*M(R, k) to recover D(k) from D(Rk). + std::vector> 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 SpinRotation::Su2& spin_U /*= SpinRotation::Su2{ 1.0, 0.0, 0.0, 1.0 }*/) const; + + std::vector>& get_rotmat_Slm() { return this->rotmat_Slm_; } + + /// test-only: inject Ms_/little_groups_/nsym_ directly, bypassing cal_Ms(), so restore_dm() + /// can be unit-tested against synthetic k-stars without a real UnitCell/K_Vectors setup. + void set_density_rotations_for_testing(const std::vector>>>& Ms, + const std::vector>& little_groups, const int nsym, const int nspin) + { + this->Ms_ = Ms; + this->little_groups_ = little_groups; + this->nsym_ = nsym; + this->nspin_ = nspin; + } + + //-------------------------------------------------------------------------------- + /// list all cells in a Born-von-Karman supercell of the given period (no LibRI dependency, + /// unlike RI_Util::get_Born_von_Karmen_cells which this mirrors for 3D periods). + static std::vector get_bvk_cells(const TC& period); + + protected: + /// set by cal_Ms() (or set_density_rotations_for_testing()); avoids reading the global + /// nspin config setting in restore_dm()/contruct_2d_rot_mat_ao(), which would otherwise + /// pull a module_parameter link dependency into every target that links this LibRI-free class. + int nspin_ = 1; + + int nsym_ = 1; + /// (nspin=4, magnetic) number of ANTIUNITARY elements Theta*g of the Shubnikov group. + /// Their orbital rotations / return lattices / Ms are appended after the nsym_ unitary + /// ones, so the raw index isym in [nsym_, nsym_+nanti_) addresses gmatrix_anti[isym-nsym_]. + int nanti_ = 0; + /// (nspin=4) true when the configuration carries a non-zero local moment. Then pure time + /// reversal is NOT a symmetry (it reverses m) and the k-star must be restored with the + /// Shubnikov elements Theta*gmatrix_anti[] instead of the generic -k shortcut. + bool magnetic_nspin4_ = false; + + double eps_ = 1e-6; + + int abfs_Lmax_ = 0; + + /// the rotation matrix under the basis of S_l^m. size: [nsym][lmax][nm*nm] + std::vector> rotmat_Slm_; + + /// bumped every time cal_rotmat_Slm() (re)fills rotmat_Slm_, so derived classes caching a + /// converted copy of rotmat_Slm_ (e.g. EXX's RI::Tensor mirror) can detect staleness + /// without recomparing the whole matrix. + int rotmat_Slm_version_ = 0; + + /// The unitary matrix associate D(Rk) with D(k) for each ibz-kpoint Rk and each symmetry operation. + /// size: [nks_ibz][nsym][nbasis*nbasis], only need to calculate once. + std::vector>>> Ms_; + + /// The little group of each ibz-kpoint: the subset of unitary space-group operations that + /// fix kvec_d_ibz modulo a reciprocal lattice vector. D(k_ibz) is averaged over this group + /// before star-expansion, since a finite-grid SCF density need not exactly respect it. + /// size: [nks_ibz][<=nsym_], always non-empty (identity is always a member). + std::vector> little_groups_; + + /// (nspin=4) the SU(2) spin-1/2 rotation U(isym) for each symmetry operation, size [nsym]. + /// The spinor AO rotation is T(isym) (x) U(isym); restore_HR_nspin4 (EXX) uses it to mix + /// the 4 spin channels of the real-space H(R). Filled in cal_Ms (identity for nspin<4). + std::vector spin_U_; + + /// irreducible sector + Irreducible_Sector irs_; + }; +} +#endif // SYMM_ROTATION_K_H diff --git a/source/source_cell/reciprocal_grid.cpp b/source/source_cell/reciprocal_grid.cpp index 43015c4f265..10f5e1a32cc 100644 --- a/source/source_cell/reciprocal_grid.cpp +++ b/source/source_cell/reciprocal_grid.cpp @@ -71,8 +71,10 @@ void ReciprocalGrid::Monkhorst_Pack(const int* nmp_in, const double* koffset_in, const int mpnz = nmp_in[2]; this->nkstot = mpnx * mpny * mpnz; - // only can renew after nkstot is estimated. - this->renew(nkstot * spin_factor()); + // only can renew after nkstot is estimated. The spin_factor() doubling (for + // K_Vectors) is applied later, on demand, by set_kup_and_kdw() itself, so it + // covers the post-symmetry-reduction k-count rather than this pre-reduction one. + this->renew(nkstot); for (int x = 1; x <= mpnx; x++) { diff --git a/source/source_io/module_parameter/input_conv.cpp b/source/source_io/module_parameter/input_conv.cpp index adeac3353c1..94ddd807d68 100644 --- a/source/source_io/module_parameter/input_conv.cpp +++ b/source/source_io/module_parameter/input_conv.cpp @@ -256,6 +256,14 @@ void Input_Conv::Convert() ModuleBase::WARNING_QUIT("Input_Conv", "EXX stress is not supported for basis_type = lcao_in_pw"); } + if (cal_exx && inp.basis_type == "lcao" && inp.kpar != 1) + { + // module_ri's real-space D(R)/Ds construction (RI_2D_Comm::split_m2D_ktoR_k) only + // sums each pool's local k-points, with no reduce_pool/reduce_all afterwards, so + // E_exx/H_exx(R)/forces would silently be built from a single pool's k-points. + ModuleBase::WARNING_QUIT("Input_Conv", "EXX (basis_type = lcao) does not support k-point parallelism (kpar > 1)"); + } + //---------------------------------------------------------- // reset symmetry flag to avoid error //---------------------------------------------------------- diff --git a/source/source_lcao/hamilt_lcao_factory.cpp b/source/source_lcao/hamilt_lcao_factory.cpp index 1cbc40c96f8..4cdbf9f966d 100644 --- a/source/source_lcao/hamilt_lcao_factory.cpp +++ b/source/source_lcao/hamilt_lcao_factory.cpp @@ -58,6 +58,7 @@ void add_dftu_op(Operator*& ops, two_center_bundle.overlap_orb_onsite.get(), orb.cutoffs(), p_dftu, inp.nspin, inp.onsite_radius, DM_in); + static_cast>*>(plus_u)->set_kv(kv); } else if (inp.dft_plus_u == 2) { diff --git a/source/source_lcao/module_dftu/dftu_nao_occ.cpp b/source/source_lcao/module_dftu/dftu_nao_occ.cpp index f0702f67c70..4f7511b507c 100644 --- a/source/source_lcao/module_dftu/dftu_nao_occ.cpp +++ b/source/source_lcao/module_dftu/dftu_nao_occ.cpp @@ -6,11 +6,67 @@ #include "source_base/parallel_reduce.h" #include "source_estate/occ_matrix.h" #include "source_lcao/hamilt_lcao.h" +#include "source_cell/module_symmetry/symmetry.h" +#include "source_cell/module_symmetry/symm_rotation_k.h" // cal_occ_mat_k / cal_occ_mat_gamma take Plus_U_Base& dftu directly and read all // occupation-matrix state (occ/save arrays, lookup table, nspin/npol, and the // occmat_ready flag) from dftu.occmat() and the Plus_U_Base accessors. +namespace +{ +// (symmetry) lazily-built rotation machinery, shared across SCF iterations +// of one run: process-lifetime static since cal_occ_mat_k has no natural +// per-ion-step owning object to hang this off (unlike the dft_plus_u=1 +// operator path, which owns its own copy). +ModuleSymmetry::Symmetry_rotation_k dftu_occ_symrot; +bool dftu_occ_symrot_built = false; + +/// @brief accumulate one k-star member's rotated S*DM product into occmat, +/// redistributing the ibz k-point's full weight (already baked into +/// srho_ibz) across all kstar_size members via Symmetry_rotation's +/// built-in 1/kstar_size scaling (see restore_dm/rot_matrix_ao). +void accumulate_occ_over_kstar(OccupationMatrix& occmat, + const UnitCell& ucell, + const Parallel_Orbitals& pv, + const K_Vectors& kv, + const std::vector>& srho_ibz, + const int ik_ibz, + const int spin, + const int nspin, + const std::vector& l_channel) +{ + const int nsym = ucell.symm.nrotk; + const size_t kstar_size = kv.kstars[ik_ibz].size(); + std::vector> sigma_y; + for (const std::pair>& isym_kvd : kv.kstars[ik_ibz]) + { + const int isym = isym_kvd.first; + std::vector> srho_rot; + if (isym < nsym) + { // unitary space-group operation (isym==0 is the identity) + srho_rot = dftu_occ_symrot.rot_matrix_ao(srho_ibz, ik_ibz, kstar_size, isym, pv); + } + else + { // antiunitary element: TRS * (spatial operation), see restore_dm + const int isym_M = ucell.symm.magnetic_nspin4 ? isym : (isym - nsym); + if (nspin == 4) + { + if (sigma_y.empty()) { sigma_y = dftu_occ_symrot.set_sigma_y_2d(pv); } + srho_rot = dftu_occ_symrot.trs_spin_rotate( + dftu_occ_symrot.rot_matrix_ao(srho_ibz, ik_ibz, kstar_size, isym_M, pv, false), + sigma_y, pv, 1.0); + } + else + { + srho_rot = dftu_occ_symrot.rot_matrix_ao(srho_ibz, ik_ibz, kstar_size, isym_M, pv, true); + } + } + DFTU_LCAO::accumulate_occ_k_for_ik(occmat, ucell, pv, srho_rot.data(), spin, l_channel); + } +} +} // namespace + void DFTU_LCAO::cal_occ_mat_k(const Parallel_Orbitals* pv, const UnitCell& ucell, @@ -42,6 +98,22 @@ void DFTU_LCAO::cal_occ_mat_k(const Parallel_Orbitals* pv, std::vector> srho(pv->nloc); + // (symmetry) when crystal symmetry reduces the k-mesh, each ik below is only + // the irreducible representative; build the AO rotation machinery once so + // its k-star can be correctly re-expanded (see accumulate_occ_over_kstar). + // Symmetry is analyzed once at the beginning and preserved by symmetrization. + // Accordingly, symrot_, dftu_occ_symrot, and the cached Ms_ remain valid and + // do not need to be rebuilt each ionic step. + const bool dftu_spacegroup_symmetry = (ModuleSymmetry::Symmetry::symm_flag == 1) && !kv.kstars.empty(); + if (dftu_spacegroup_symmetry && !dftu_occ_symrot_built) + { + const std::array period{ kv.nmp[0], kv.nmp[1], kv.nmp[2] }; + dftu_occ_symrot.find_irreducible_sector(ucell.symm, ucell.atoms, ucell.st, + ModuleSymmetry::Symmetry_rotation_k::get_bvk_cells(period), period, ucell.lat); + dftu_occ_symrot.cal_Ms(kv, ucell, *pv, nspin); + dftu_occ_symrot_built = true; + } + for (int ik = 0; ik < kv.get_nks(); ik++) { // srho(mu,nu) = \sum_{iw} S(mu,iw)*dm_k(iw,nu) @@ -82,7 +154,20 @@ void DFTU_LCAO::cal_occ_mat_k(const Parallel_Orbitals* pv, const int spin = kv.isk[ik]; // Walk (it, ia, l, n=0) and accumulate each qualifying channel - accumulate_occ_k_for_ik(dftu.occmat(), ucell, *pv, srho.data(), spin, l_channel); + if (dftu_spacegroup_symmetry) + { + // kv.kstars/Ms_ are sized per spin and indexed by GLOBAL ibz position + // (kv.kstars.size() == nks_ibz); ik is local to this k-point pool, so + // map it to the global k index first (kv.ik2iktot), then wrap into the + // per-spin ibz range (mirrors RI_2D_Comm::split_m2D_ktoR_k's + // "ik % ik_list.size()", but on the global index rather than the local one). + const int ik_ibz = kv.ik2iktot[ik] % static_cast(kv.kstars.size()); + accumulate_occ_over_kstar(dftu.occmat(), ucell, *pv, kv, srho, ik_ibz, spin, nspin, l_channel); + } + else + { + accumulate_occ_k_for_ik(dftu.occmat(), ucell, *pv, srho.data(), spin, l_channel); + } } // ik // MPI Allreduce + symmetrize per (iat, l, n=0) channel across all ranks diff --git a/source/source_lcao/module_dftu/dftu_nao_op.cpp b/source/source_lcao/module_dftu/dftu_nao_op.cpp index 282aad6aa2d..dd518b44c70 100644 --- a/source/source_lcao/module_dftu/dftu_nao_op.cpp +++ b/source/source_lcao/module_dftu/dftu_nao_op.cpp @@ -9,6 +9,10 @@ #include "source_lcao/module_operator_lcao/operator_lcao.h" #include "source_pw/module_pwdft/dftu_base.h" #include "source_base/parallel_reduce.h" +#include "source_cell/klist.h" +#include "source_cell/module_symmetry/symmetry.h" + +#include #include "dftu_nao_adj.h" #include "dftu_nao_fs_r.h" @@ -124,6 +128,46 @@ void hamilt::DFTU_onsite>::contributeHR() const Parallel_Orbitals* pv = this->hR->get_atom_pair(0).get_paraV(); // nlm_tot is precomputed in the constructor (structure snapshot) + // (symmetry) when crystal symmetry reduces the k-mesh, dm_->get_DMR_pointer() + // was Fourier-transformed from the irreducible k-points only and is not + // actually symmetric; reconstruct the full-BZ DMR once here (reused by every + // atom below) via the same D(k) restoration EXX already uses for its own + // real-space density matrix (ModuleSymmetry::Symmetry_rotation::restore_dm). + std::unique_ptr> dmr_sym; + if (!this->dftu->is_occmat_ready() && this->kv_ != nullptr && ModuleSymmetry::Symmetry::symm_flag == 1 + && !this->kv_->kstars.empty()) + { + if (!this->symrot_built_) + { + const std::array period{ this->kv_->nmp[0], this->kv_->nmp[1], this->kv_->nmp[2] }; + // for return_lattice to calculate Ms + this->symrot_.find_irreducible_sector(this->ucell->symm, this->ucell->atoms, this->ucell->st, + ModuleSymmetry::Symmetry_rotation_k::get_bvk_cells(period), period, this->ucell->lat); + this->symrot_.cal_Ms(*this->kv_, *this->ucell, *pv, this->nspin); + this->symrot_built_ = true; + } + const int nspin0 = (this->nspin == 2) ? 2 : 1; + // (k-point pools, KPAR>1) restore_dm() now returns only the stars of THIS pool's own + // local irreducible k-points (see its definition for why that's enough); kvec_d_full + // must be built the same way -- one entry per star member of each local ibz-k, + // enumerated in the same order restore_dm uses for its spin-0 block (kv.ik2iktot maps + // the spin-0 and spin-1 blocks to the same sequence of global ibz indices, so a single + // list built from the spin-0 mapping is valid for the whole nspin0-block DensityMatrix). + const int nk_local = this->kv_->get_nks() / nspin0; + const int nks_ibz_global = static_cast(this->kv_->kstars.size()); + std::vector> kvec_d_full; + for (int ik_local = 0; ik_local < nk_local; ++ik_local) + { + const int ik_ibz = this->kv_->ik2iktot[ik_local] % nks_ibz_global; + for (const std::pair>& isym_kvd : this->kv_->kstars[ik_ibz]) { kvec_d_full.push_back(isym_kvd.second); } + } + const std::vector> dmk_full = this->symrot_.restore_dm(*this->kv_, this->dm_->get_DMK_vector(), *pv); + dmr_sym.reset(new elecstate::DensityMatrix(pv, nspin0, kvec_d_full, static_cast(kvec_d_full.size()))); + dmr_sym->init_DMR(*this->dm_->get_DMR_pointer(1)); + dmr_sym->get_DMK_vector() = dmk_full; + dmr_sym->cal_DMR(); + } + // loop over all Hubbard-projector center atoms (iat0) int atom_index = 0; for (int iat0 = 0; iat0 < this->ucell->nat; iat0++) @@ -148,6 +192,7 @@ void hamilt::DFTU_onsite>::contributeHR() // DMR is guaranteed ready here: otherwise the early exit above // would have returned. DMR index is 1-based, hence +1. const hamilt::HContainer* dmr = this->dm_->get_DMR_pointer(this->current_spin + 1); + if (dmr_sym) { dmr = dmr_sym->get_DMR_pointer(this->current_spin + 1); } DFTU_LCAO::compute_occ_from_dmr(*this->ucell, *this->dftu, iat0, diff --git a/source/source_lcao/module_dftu/dftu_nao_op.h b/source/source_lcao/module_dftu/dftu_nao_op.h index 51d919d3166..8b3a6c28d90 100644 --- a/source/source_lcao/module_dftu/dftu_nao_op.h +++ b/source/source_lcao/module_dftu/dftu_nao_op.h @@ -5,12 +5,14 @@ #include "source_lcao/module_operator_lcao/operator_lcao.h" #include "source_hamilt/module_hcontainer/hcontainer.h" #include "source_lcao/module_dftu/dftu_nao_adj.h" +#include "source_cell/module_symmetry/symm_rotation_k.h" // ModuleSymmetry::Symmetry_rotation_k (value member); LibRI-free #include class Plus_U_Base; class TwoCenterIntegrator; class UnitCell; +class K_Vectors; namespace elecstate { @@ -62,6 +64,12 @@ class DFTU_onsite> : public OperatorLCAO */ void contributeHR() override; + /// @brief provide the K_Vectors needed to restore the full-BZ density matrix + /// under crystal symmetry (kv.kstars); set once after construction from + /// HamiltLCAO. When left null, contributeHR() falls back to the + /// unsymmetrized DMR (previous behavior). + void set_kv(const K_Vectors* kv_in) { this->kv_ = kv_in; } + private: const UnitCell* ucell = nullptr; @@ -83,6 +91,17 @@ class DFTU_onsite> : public OperatorLCAO /// @brief cached overlap values; structure snapshot computed /// once in the constructor, reused across SCF iterations of one ionic step DFTU_LCAO::NlmTot nlm_tot; + + /// @brief K_Vectors of the current run; only used (when non-null) to restore + /// the full-BZ density matrix under crystal symmetry before it is read into + /// the DFT+U occupation matrix. Set via set_kv() after construction. + const K_Vectors* kv_ = nullptr; + + /// @brief LibRI-independent k-space rotation machinery (shared base of EXX's own + /// Symmetry_rotation) to reconstruct D(k) at every k-star member from D(k_ibz); + /// built lazily (once) the first time it is needed (rebuilt every ionic step). + ModuleSymmetry::Symmetry_rotation_k symrot_; + bool symrot_built_ = false; }; } // namespace hamilt diff --git a/source/source_lcao/module_dftu/test/CMakeLists.txt b/source/source_lcao/module_dftu/test/CMakeLists.txt index 8ca8651135f..ba9d848ae74 100644 --- a/source/source_lcao/module_dftu/test/CMakeLists.txt +++ b/source/source_lcao/module_dftu/test/CMakeLists.txt @@ -15,7 +15,7 @@ AddTest( if(ENABLE_LCAO AND ENABLE_MPI) AddTest( TARGET dftu_lcao_test - LIBS parameter psi base device container + LIBS parameter psi base device container symmetry SOURCES dftu_lcao_test.cpp ../dftu_nao_op.cpp ../dftu_nao_adj.cpp ../dftu_nao_pots.cpp ../dftu_nao_fs_r.cpp ../dftu_nao_for_r.cpp ../dftu_nao_str_r.cpp ../../../source_estate/module_dm/density_matrix.cpp ../../../source_estate/module_dm/density_matrix_io.cpp diff --git a/source/source_lcao/module_rdmft/rdmft.cpp b/source/source_lcao/module_rdmft/rdmft.cpp index 5a27d9d86f8..094e71dbe1b 100644 --- a/source/source_lcao/module_rdmft/rdmft.cpp +++ b/source/source_lcao/module_rdmft/rdmft.cpp @@ -158,8 +158,8 @@ void RDMFT::init(Parallel_Orbitals& ParaV_in, { const std::array& period = RI_Util::get_Born_vonKarmen_period(*kv); this->symrot_exx.find_irreducible_sector(ucell->symm, ucell->atoms, ucell->st, - RI_Util::get_Born_von_Karmen_cells(period), period, ucell->lat); - this->symrot_exx.cal_Ms(*kv, *ucell, *ParaV); + RI_Util::get_Born_von_Karmen_cells(period), period, ucell->lat, PARAM.globalv.global_out_dir); + this->symrot_exx.cal_Ms(*kv, *ucell, *ParaV, nspin); } if (this->exx_info_->info_ri.real_number) diff --git a/source/source_lcao/module_ri/exx_lri_interface.hpp b/source/source_lcao/module_ri/exx_lri_interface.hpp index f00011ab833..dce0a78d463 100644 --- a/source/source_lcao/module_ri/exx_lri_interface.hpp +++ b/source/source_lcao/module_ri/exx_lri_interface.hpp @@ -121,9 +121,9 @@ void Exx_LRI_Interface::exx_before_all_runners( const std::array& period = RI_Util::get_Born_vonKarmen_period(kv); this->symrot_.find_irreducible_sector( ucell.symm, ucell.atoms, ucell.st, - RI_Util::get_Born_von_Karmen_cells(period), period, ucell.lat); + RI_Util::get_Born_von_Karmen_cells(period), period, ucell.lat, PARAM.globalv.global_out_dir); this->symrot_.set_abfs_Lmax(Exx_Abfs::Construct_Orbs::get_Lmax(this->exx_ptr->abfs)); - this->symrot_.cal_Ms(kv, ucell, pv); + this->symrot_.cal_Ms(kv, ucell, pv, PARAM.inp.nspin); } } diff --git a/source/source_lcao/module_ri/module_exx_symmetry/CMakeLists.txt b/source/source_lcao/module_ri/module_exx_symmetry/CMakeLists.txt index 2af8a5bf7fa..e5b6e32ce24 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/CMakeLists.txt +++ b/source/source_lcao/module_ri/module_exx_symmetry/CMakeLists.txt @@ -1,8 +1,6 @@ if (ENABLE_LIBRI) list(APPEND objects - irreducible_sector.cpp - irred_sec_bvk.cpp symm_rotation.cpp symm_rot_out.cpp ) diff --git a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.cpp b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.cpp index a1942918f70..b8261ad6fe1 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.cpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.cpp @@ -7,6 +7,7 @@ #include "source_base/tool_title.h" #include "source_base/timer.h" #include "source_base/mathzone.h" +#include "source_lcao/module_ri/ri_util.h" namespace ModuleSymmetry { @@ -16,590 +17,6 @@ namespace ModuleSymmetry this->abfs_l_nchi_ = abfs_l_nchi; for (auto& abfs_T : abfs_l_nchi) { this->abfs_Lmax_ = std::max(this->abfs_Lmax_, static_cast(abfs_T.size()) - 1); } } - void Symmetry_rotation::cal_Ms(const K_Vectors& kv, - //const std::vector>& kstars, - const UnitCell& ucell, const Parallel_2D& pv) - { - ModuleBase::TITLE("Symmetry_rotation", "cal_Ms"); - ModuleBase::timer::start("Symmetry_rotation", "cal_Ms"); - - this->nsym_ = ucell.symm.nrotk; - this->nanti_ = ucell.symm.nrotk_anti; - this->magnetic_nspin4_ = ucell.symm.magnetic_nspin4; - this->eps_ = ucell.symm.epsilon; - if (this->irs_.invmap_.empty()) - { - this->irs_.invmap_.resize(ucell.symm.nrotk); - ucell.symm.gmatrix_invmap(ucell.symm.gmatrix, ucell.symm.nrotk, this->irs_.invmap_.data()); - } - // 1. calculate the rotation matrix in real spherical harmonics representation for each symmetry operation: [T_l (isym)]_mm' - const int nop_tot = this->nsym_ + this->nanti_; - std::vector gmatc(nop_tot); - for (int i = 0;i < nsym_;++i) { gmatc[i] = this->irs_.direct_to_cartesian(ucell.symm.gmatrix[i], ucell.latvec); } - for (int j = 0;j < this->nanti_;++j) - { gmatc[nsym_ + j] = this->irs_.direct_to_cartesian(ucell.symm.gmatrix_anti[j], ucell.latvec); } - this->cal_rotmat_Slm(gmatc.data(), std::max(this->abfs_Lmax_, ucell.lmax), nop_tot); - - // 1.5 (nspin=4) the SU(2) spin-1/2 rotation U(isym) for each symmetry operation. The AO - // rotation matrix M becomes the spinor operator T(isym) (x) U(isym) so that the same - // gemm D(k)=M^dagger D(k_ibz) M rotates both the orbital and the spin part at once. - // For an antiunitary element Theta*g only the spatial part g enters M here; the Theta - // (sigma_y (.)^* sigma_y) is applied afterwards in restore_dm. - std::vector spin_U(nop_tot, SpinRotation::Su2{ 1.0, 0.0, 0.0, 1.0 }); - if (PARAM.inp.nspin == 4) - { - for (int i = 0;i < nop_tot;++i) { spin_U[i] = SpinRotation::so3_to_su2(gmatc[i]); } - } - this->spin_U_ = spin_U; // keep for restore_HR_nspin4 (real-space EXX H(R) spin mixing) - - // A k-star contains only one operation per distinct k point. The other - // operations fixing k (modulo a reciprocal lattice vector) must still - // be averaged: a finite-grid SCF density need not respect this little group. - const int nks_ibz = kv.kstars.size(); - this->Ms_.assign(nks_ibz, {}); - this->little_groups_.assign(nks_ibz, {}); - for (int ik_ibz = 0; ik_ibz < nks_ibz; ++ik_ibz) - { - std::set needed; - for (const auto& member : kv.kstars[ik_ibz]) - { - const int op = (!this->magnetic_nspin4_ && member.first >= nsym_) - ? member.first - nsym_ : member.first; - needed.insert(op); - } - for (int op = 0; op < nsym_; ++op) - { - const auto delta = kv.kvec_d[ik_ibz] * ucell.symm.kgmatrix[op] - kv.kvec_d[ik_ibz]; - if (std::abs(delta.x - std::round(delta.x)) < this->eps_ - && std::abs(delta.y - std::round(delta.y)) < this->eps_ - && std::abs(delta.z - std::round(delta.z)) < this->eps_) - { - this->little_groups_[ik_ibz].push_back(op); - needed.insert(op); - } - } - for (const int op : needed) - { - this->Ms_[ik_ibz][op] = this->contruct_2d_rot_mat_ao( - ucell.symm, ucell.atoms, ucell.st, kv.kvec_d[ik_ibz], op, pv, spin_U[op]); - } - } - // output Ms of isym=1 - // std::ofstream ofs("Ms_kibz7_sym7.dat"); - // for (int i = 0;i < pv.get_row_size();++i) - // { - // for (int j = 0;j < pv.get_col_size();++j) - // { - // ofs << std::setprecision(10) << this->Ms_[7][7][j * pv.get_col_size() + i] << " "; - // } - // ofs << std::endl; - // } - // ofs << std::endl; - // ofs.close(); - - ModuleBase::timer::end("Symmetry_rotation", "cal_Ms"); - } - - std::vector>> Symmetry_rotation::restore_dm(const K_Vectors& kv, - const std::vector>>& dm_k_ibz, const Parallel_2D& pv)const - { - ModuleBase::TITLE("Symmetry_rotation", "restore_dm"); - ModuleBase::timer::start("Symmetry_rotation", "restore_dm"); - std::vector>> dm_k_full; - int nspin0 = PARAM.inp.nspin == 2 ? 2 : 1; - dm_k_full.reserve(kv.get_nkstot_nospin() * nspin0); //nkstot_nospin didn't doubled by spin - int nk = kv.get_nkstot() / nspin0; - - // (nspin=4) Sigma_y = I (x) sigma_y for the time-reversal spin flip; k-independent, build once. - std::vector> sigma_y; - if (PARAM.inp.nspin == 4) { sigma_y = this->set_sigma_y_2d(pv); } - - for (int is = 0;is < nspin0;++is) - { - for (int ik_ibz = 0;ik_ibz < nk;++ik_ibz) - { - // P_k D = |G_k|^{-1} sum_g M_g^T D M_g^*. This preserves - // Hermiticity and makes restoration independent of the chosen - // star representative; rotating just one arbitrary D does not. - const auto& little_group = this->little_groups_.at(ik_ibz); - assert(!little_group.empty()); - std::vector> projected = dm_k_ibz[ik_ibz + is * nk]; - if (little_group.size() > 1) - { - std::fill(projected.begin(), projected.end(), 0.0); - for (const int op : little_group) - { - const auto rotated = this->rot_matrix_ao( - dm_k_ibz[ik_ibz + is * nk], ik_ibz, little_group.size(), op, pv); - for (size_t i = 0; i < projected.size(); ++i) - { - projected[i] += rotated[i]; - } - } - } - for (auto& isym_kvd : kv.kstars[ik_ibz]) - { - if (isym_kvd.first == 0) - { - double factor = 1.0 / static_cast(kv.kstars[ik_ibz].size()); - std::vector> dm_scaled(pv.get_local_size()); - for (int i = 0;i < pv.get_local_size();++i) { dm_scaled[i] = factor * projected[i]; } - dm_k_full.push_back(dm_scaled); - } - else if (isym_kvd.first < nsym_) - { //space group operations - dm_k_full.push_back(this->rot_matrix_ao(projected, ik_ibz, kv.kstars[ik_ibz].size(), isym_kvd.first, pv)); - } - else - { // antiunitary elements: Theta * (spatial operation) - // D(Theta*g k_ibz) = sigma_y [D(g k_ibz)]^* sigma_y with D(g k_ibz) = M^dagger D M. - // For nspin=4, first do the (non-conjugated) spatial rotation, then the spin flip; - // for nspin<4 (Theta=K) the original TRS_conj path already gives the conjugate. - // - // Which spatial operation the index denotes depends on the regime, matching - // how the k-reduction filled kgmatrix[] (see K_Vectors::reduce_by_symmetry): - // - nspin=4 magnetic (Shubnikov): index j+nsym_ is the antiunitary element - // Theta*gmatrix_anti[j]; its Ms is stored under the RAW key j+nsym_. - // - otherwise (grey group / nspin<4): index i+nsym_ is Theta*gmatrix[i], - // i.e. the unitary operation i, whose Ms is stored under key i. - const int isym_M = this->magnetic_nspin4_ ? isym_kvd.first : (isym_kvd.first - nsym_); - if (PARAM.inp.nspin == 4) - { - // m=0: gray group: the space-group part of anti-unitary elements are the same of the unitary elements, isym_M < nsym_ - // m!=0: Shubnikov group: using different space-group part of anti-unitary elements stored in gmatrix_anti with isym_M >= nsym_ - dm_k_full.push_back(this->trs_spin_rotate( - this->rot_matrix_ao(projected, ik_ibz, kv.kstars[ik_ibz].size(), isym_M, pv, false), - sigma_y, pv, 1.0)); - } - else - { - dm_k_full.push_back(this->rot_matrix_ao(projected, ik_ibz, kv.kstars[ik_ibz].size(), isym_M, pv, true)); - } - } - } - } - } - // test for output -/* - std::ofstream ofs("DM.dat"); - int ik = 0; - for (int ikibz = 0;ikibz < kv.get_nkstot() / nspin0;++ikibz) - for (auto& isym_kvd : kv.kstars[ikibz]) - { - ofs << "isym=" << isym_kvd.first << std::endl; - ofs << " k = " << isym_kvd.second.x << " " << isym_kvd.second.y << " " << isym_kvd.second.z << std::endl; - ofs << "DM(k):" << std::endl; - for (int i = 0;i < pv.get_row_size();++i) - { - for (int j = 0;j < pv.get_col_size();++j) - { - ofs << dm_k_full[ik][j * pv.get_row_size() + i] << " "; - } - ofs << std::endl; - } - ++ik; - ofs << std::endl; - } - ofs.close(); -*/ - ModuleBase::timer::end("Symmetry_rotation", "restore_dm"); - return dm_k_full; - } - std::vector> Symmetry_rotation::restore_dm(const K_Vectors& kv, - const std::vector>& dm_k_ibz, const Parallel_2D& pv)const - { - return dm_k_ibz;// do nothing for gamma_only - } - - // calculate Wigner D matrix - double Symmetry_rotation::wigner_d(const double beta, const int l, const int m1, const int m2) const - { - auto factorial = [](int n) -> int { - int result = 1; - for (int i = 1;i <= n;++i) { result *= i; -} - return result; - }; - double result = 0.0; - for (int i = std::max(0, m2 - m1);i <= std::min(l - m1, l + m2);++i) { - result += std::pow(-1, i) * std::sqrt(factorial(l + m1) * factorial(l - m1) * factorial(l + m2) * factorial(l - m2)) - * std::pow(std::cos(beta / 2), 2 * l + m2 - m1 - 2 * i) * std::pow(-std::sin(beta / 2), m1 - m2 + 2 * i) - / (factorial(i) * factorial(l - m1 - i) * factorial(l + m2 - i) * factorial(i - m2 + m1)); -} - return result; - } - - std::complex Symmetry_rotation::wigner_D(const TCdouble& euler_angle, const int l, const int m1, const int m2, const bool inv) const - { - std::complex prefac(inv ? std::pow(-1, l) : 1, 0); - return std::exp(-ModuleBase::IMAG_UNIT * static_cast(m1) * euler_angle.x) - * std::exp(-ModuleBase::IMAG_UNIT * static_cast(m2) * euler_angle.z) - * wigner_d(euler_angle.y, l, m1, m2) * prefac; - } - - // c^l_{m1, m2}= - std::complex Symmetry_rotation::ovlp_Ylm_Slm(const int l, const int m1, const int m2) const - { - if (m1 == m2) - { - if (m1 == 0) { return 1.0; -} - if (m1 > 0) { return 1 / std::sqrt(2); -} - if (m1 < 0) { return std::pow(-1, m1) * ModuleBase::IMAG_UNIT / std::sqrt(2); -} - } - else if (m1 == -m2) - { - if (m1 > 0) { return -ModuleBase::IMAG_UNIT / std::sqrt(2); -} - if (m1 < 0) { return std::pow(-1, m1) / std::sqrt(2); -} - } - return 0.0; - } - - // reference: https://github.com/minyez/abf_trans/blob/f9e68e68069a94610d89e077bfe6e8ffac0b097d/src/rotate.cpp#L118 - // because the atom position here is row vector, the original gmatrix(eular angle) is transposed. - // gmatc: the rotation matrix under the basis of cartesian coordinates - // gmatc should be a rotation matrix, i.e. det(gmatc)=1 - TCdouble Symmetry_rotation::get_euler_angle(const ModuleBase::Matrix3& gmatc) const - { - double threshold = this->eps_; - double alpha = 0.0, beta = 0.0, gamma = 0.0; - if (std::fabs(gmatc.e32) > threshold || std::fabs(gmatc.e31) > threshold) // sin(beta) is not zero - { - // use the 2-angle elements to get alpha and gamma - alpha = std::atan2(gmatc.e32, gmatc.e31); - if (alpha < 0) { alpha += 2 * ModuleBase::PI; -} - gamma = std::atan2(gmatc.e23, -gmatc.e13); - if (gamma < 0) { gamma += 2 * ModuleBase::PI; -} - // use the larger one of 2-angle elements to calculate beta - if (std::fabs(gmatc.e32) > std::fabs(gmatc.e31)) { - beta = std::atan2(gmatc.e32 / std::sin(alpha), gmatc.e33); - } else { - beta = std::atan2(gmatc.e31 / std::cos(alpha), gmatc.e33); -} - } - else - {//sin(beta)=0, beta = 0 or pi, only (alpha+gamma) or (alpha-gamma) is important. now assign this to alpha. - alpha = std::atan2(gmatc.e12, gmatc.e11); - if (alpha < 0) { alpha += 2 * ModuleBase::PI; -} - // if beta=0, gmatc.e11=cos(alpha+gamma), gmatc.e21=sin(alpha+gamma) - // if beta=pi, gmatc.e11=cos(pi+alpha-gamma), gmatc.e21=sin(pi+alpha-gamma) - if (gmatc.e33 > 0) - { - beta = 0; - gamma = 0; //alpha+gamma=alpha => gamma=0 - } - else - { - beta = ModuleBase::PI; - gamma = ModuleBase::PI;// pi+alpha-gamma=alpha => gamma=pi - } - } - return TCdouble(alpha, beta, gamma); - } - - // in: the real value of m in range {-l, -l+1, ..., 0, ..., l-1, l} - // out: the index of the orbital in a fixed {n， l}, i.e. the index in array [0, 1, -1, 2, -2, ...] - inline int m2im(int m) - { - return (m > 0 ? 2 * m - 1 : -2 * m); - } - - /// T_mm' = [c^\dagger D c]_mm' - void Symmetry_rotation::cal_rotmat_Slm(const ModuleBase::Matrix3* gmatc, const int lmax, const int nop) - { - const int nop_tot = (nop < 0) ? this->nsym_ : nop; - auto set_integer = [](RI::Tensor>& mat) -> void - { - double zero_thres = 1e-10; - for (int i = 0;i < mat.shape[0];++i) { - for (int j = 0;j < mat.shape[1];++j) - { - if (std::abs(mat(i, j).real() - std::round(mat(i, j).real())) < zero_thres) { mat(i, j).real(std::round(mat(i, j).real())); -} - if (std::abs(mat(i, j).imag() - std::round(mat(i, j).imag())) < zero_thres) { mat(i, j).imag(std::round(mat(i, j).imag())); -} - } -} - }; - this->rotmat_Slm_.resize(nop_tot); - // c matrix is independent on isym - std::vector>> c_mm(lmax + 1); - for (int l = 0;l <= lmax;++l) { - c_mm[l] = RI::Tensor>({ size_t(2 * l + 1), size_t(2 * l + 1) }); -} - for (int l = 0;l <= lmax;++l) { - for (int m1 = -l;m1 <= l;++m1) { - for (int m2 = -l;m2 <= l;++m2) { - c_mm[l](m2im(m1), m2im(m2)) = ovlp_Ylm_Slm(l, m1, m2); -} -} -} - - for (int isym = 0;isym < nop_tot;++isym) - { - // if R is a reflection operation, calculate D^l(R)=(-1)^l*D^l(IR), so the euler angle of (IR) is needed. - TCdouble euler_angle = get_euler_angle(gmatc[isym].Det() > 0 ? - gmatc[isym] : gmatc[isym] * ModuleBase::Matrix3(-1, 0, 0, 0, -1, 0, 0, 0, -1)); - - this->rotmat_Slm_[isym].resize(lmax + 1); - for (int l = 0;l <= lmax;++l) - {// wigner D matrix - RI::Tensor> D_mm({ size_t(2 * l + 1), size_t(2 * l + 1) }); - for (int m1 = -l;m1 <= l;++m1) { - for (int m2 = -l;m2 <= l;++m2) { - D_mm(m2im(m1), m2im(m2)) = wigner_D(euler_angle, l, m1, m2, (gmatc[isym].Det() < 0)); -} -} - this->rotmat_Slm_[isym][l] = c_mm[l].dagger() * D_mm * c_mm[l]; - // set_integer(this->rotmat_Slm_[isym][l]); - } - } - /* - std::vector euler_angles_test(nsym_); - for (int isym = 0;isym < nsym_;++isym) euler_angles_test[isym] = - get_euler_angle(gmatc[isym].Det() > 0 ? gmatc[isym] : gmatc[isym] * ModuleBase::Matrix3(-1, 0, 0, 0, -1, 0, 0, 0, -1)); - - auto test_Tmm = [&]()-> void - { - std::ofstream ofs("Tlm.dat"); - for (int isym = 0;isym < nsym_;++isym) - { - ofs << "isym=" << isym << std::endl; - ofs << "gmatrix_cart=" << std::endl; - ofs << gmatc[isym].e11 << " " << gmatc[isym].e12 << " " << gmatc[isym].e13 << std::endl; - ofs << gmatc[isym].e21 << " " << gmatc[isym].e22 << " " << gmatc[isym].e23 << std::endl; - ofs << gmatc[isym].e31 << " " << gmatc[isym].e32 << " " << gmatc[isym].e33 << std::endl; - ofs << "gmatrix_direct=" << std::endl; - ofs << ucell.symm.gmatrix[isym].e11 << " " << ucell.symm.gmatrix[isym].e12 << " " << ucell.symm.gmatrix[isym].e13 << std::endl; - ofs << ucell.symm.gmatrix[isym].e21 << " " << ucell.symm.gmatrix[isym].e22 << " " << ucell.symm.gmatrix[isym].e23 << std::endl; - ofs << ucell.symm.gmatrix[isym].e31 << " " << ucell.symm.gmatrix[isym].e32 << " " << ucell.symm.gmatrix[isym].e33 << std::endl; - ofs << "kgmatrix_direct=" << std::endl; - ofs << ucell.symm.kgmatrix[isym].e11 << " " << ucell.symm.kgmatrix[isym].e12 << " " << ucell.symm.kgmatrix[isym].e13 << std::endl; - ofs << ucell.symm.kgmatrix[isym].e21 << " " << ucell.symm.kgmatrix[isym].e22 << " " << ucell.symm.kgmatrix[isym].e23 << std::endl; - ofs << ucell.symm.kgmatrix[isym].e31 << " " << ucell.symm.kgmatrix[isym].e32 << " " << ucell.symm.kgmatrix[isym].e33 << std::endl; - ofs << "euler_angle/pi: " << euler_angles_test[isym].x / ModuleBase::PI << " " - << euler_angles_test[isym].y / ModuleBase::PI << " " << euler_angles_test[isym].z / ModuleBase::PI << std::endl; - for (int l = 0;l <= lmax;++l) - for (int i = 0;i < 2 * l + 1;++i) - { - for (int j = 0;j < 2 * l + 1;++j) ofs << this->rotmat_Slm_[isym][l](i, j) << " "; - ofs << std::endl; - } - } - ofs.close(); - }; - test_Tmm(); - */ - } - - void Symmetry_rotation::set_block_to_mat2d(const int starti, const int startj, const RI::Tensor>& block, - std::vector>& obj_mat, const Parallel_2D& pv, const bool trans) const - { // caution: ComplaxMatrix is row-major(col-continuous), but obj_mat is col-major(row-continuous) - for (int j = 0;j < block.shape[0];++j) {//outside dimension - for (int i = 0;i < block.shape[1];++i) { //inside dimension - if (pv.in_this_processor(starti + i, startj + j)) - { - int index = pv.global2local_col(startj + j) * pv.get_row_size() + pv.global2local_row(starti + i); - obj_mat[index] = trans ? block(i, j) : block(j, i); - } -} -} - } - - void Symmetry_rotation::set_block_to_mat2d(const int starti, const int startj, const RI::Tensor>& block, - std::vector& obj_mat, const Parallel_2D& pv, const bool trans) const - { // caution: ComplaxMatrix is row-major(col-continuous), but obj_mat is col-major(row-continuous) - for (int j = 0;j < block.shape[0];++j) {//outside dimension - for (int i = 0;i < block.shape[1];++i) { //inside dimension - if (pv.in_this_processor(starti + i, startj + j)) - { - int index = pv.global2local_col(startj + j) * pv.get_row_size() + pv.global2local_row(starti + i); - obj_mat[index] = trans ? block(i, j).real() : block(j, i).real(); - } -} -} - } - - // 2d-block parallized rotation matrix in AO-representation, denoted as M. - // finally we will use D(k)=M(R, k)^\dagger*D(Rk)*M(R, k) to D(k) from D(Rk) in cal_Ms. - std::vector> Symmetry_rotation::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 SpinRotation::Su2& spin_U) const - { - const bool soc = (PARAM.inp.nspin == 4); - const int npol = soc ? 2 : 1; // spinor: global AO index is spin-fast interleaved, I = npol*iw_orb + s - std::vector> M_isym(pv.get_local_size(), 0.0); - // isym >= symm.nrotk addresses the antiunitary coset (spatial part gmatrix_anti[isym-nrotk]), - // whose atom map lives in a separate table. - const int nrotk_u = symm.nrotk; - auto rotated_atom = [&symm, nrotk_u](const int is, const int iat) -> int - { - return (is < nrotk_u) ? symm.get_rotated_atom(is, iat) - : symm.get_rotated_atom_anti(is - nrotk_u, iat); - }; - for (int iat1 = 0;iat1 < cell_st.nat;++iat1) - { - int it = cell_st.iat2it[iat1]; // it1=it2 - int ia1 = cell_st.iat2ia[iat1]; - int iat2 = rotated_atom(isym, iat1); //iat2=rot(iat1) - int ia2 = cell_st.iat2ia[iat2]; - // cal phase factor from return lattice: exp(-ik_ibz*O) - double arg = -2 * ModuleBase::PI * kvec_d_ibz * this->irs_.return_lattice_[iat1][isym]; - std::complexphase_factor = std::complex(std::cos(arg), std::sin(arg)); - int iw1start = atoms[it].stapos_wf + ia1 * atoms[it].nw; - int iw2start = atoms[it].stapos_wf + ia2 * atoms[it].nw; - int iw = 0; - while (iw < atoms[it].nw) - { - int l = atoms[it].iw2l[iw]; - int nm = 2 * l + 1; - //caution: the order of m in orbitals may be different from increasing - if (!soc) - { - set_block_to_mat2d(iw2start + iw, iw1start + iw, - phase_factor * this->rotmat_Slm_[isym][l], M_isym, pv, true); - } - else - { - // M = T(isym) (x) U(isym): scatter phase * T_l(m,m') * U(a,b) to the interleaved - // spinor positions (row = rotated atom/spin, col = original atom/spin). For nspin=4 - // stapos_wf already carries the npol factor, so the per-atom offset is ia*nw*npol - // and the within-atom spinor index is (iw_orb)*npol + spin (spin is the fast index). - const int base2 = atoms[it].stapos_wf + ia2 * atoms[it].nw * npol; - const int base1 = atoms[it].stapos_wf + ia1 * atoms[it].nw * npol; - const RI::Tensor>& Tl = this->rotmat_Slm_[isym][l]; - for (int m = 0;m < nm;++m) - { - for (int mp = 0;mp < nm;++mp) - { - const std::complex t = phase_factor * Tl(m, mp); - for (int a = 0;a < npol;++a) - { - for (int b = 0;b < npol;++b) - { - const int gi = base2 + (iw + m) * npol + a; - const int gj = base1 + (iw + mp) * npol + b; - if (pv.in_this_processor(gi, gj)) - { - const int index = pv.global2local_col(gj) * pv.get_row_size() + pv.global2local_row(gi); - // M(isym) = T_l (x) U is the spinor rep, with U = so3_to_su2 placed as-is: - // M[(m,a),(m',b)] = phase * T_l(m,m') * U_{ab}, U_{ab} = spin_U[a*npol + b]. - // Both T_l (rotmat_Slm) and U are ANTI-homomorphisms here (row-vector / R^T convention: - // rotmat_Slm(g)=R_orb(g)^{-1}, so3_to_su2 likewise), so this M is a consistent rep - // and rot_matrix_ao's stored-DM rotation M^T D M^* is exact for ALL ops. - M_isym[index] = t * spin_U[a * npol + b]; - } - } - } - } - } - } - iw += nm; - } - } - return M_isym; - } - - // D(k) = M^T(R, k) D(k_ibz) M^*(R, k), if D(k) is col-maj - // D^T(k) = M^\dagger(R, k) D^T(k_ibz) M(R, k), if D(k) is row-maj - // Ds from RI_2D_Comm are row-maj - // the link ik_ibz-isym-ik can be found in kstars. - std::vector> Symmetry_rotation::rot_matrix_ao(const std::vector>& DMkibz, - const int ik_ibz, const int kstar_size, const int isym, const Parallel_2D& pv, const bool TRS_conj) const - { - std::vector> DMk(pv.nloc, 0.0); - std::vector> DMkibz_M(pv.nloc, 0.0); // intermediate result - const char dagger = 'C'; - const char transpose = 'T'; - const char notrans = 'N'; - std::complex alpha(1.0, 0.0); - const std::complex beta(0.0, 0.0); - const int nbasis = pv.get_global_row_size(); - const int i1 = 1; - if (TRS_conj) - { - // D^T* = M^T [M^T (D^T)^T]^\dagger - ScalapackConnector::gemm(transpose, transpose, nbasis, nbasis, nbasis, - alpha, this->Ms_[ik_ibz].at(isym).data(), i1, i1, pv.desc, DMkibz.data(), i1, i1, pv.desc, - beta, DMkibz_M.data(), i1, i1, pv.desc); - alpha.real(1.0 / static_cast(kstar_size)); - ScalapackConnector::gemm(transpose, dagger, nbasis, nbasis, nbasis, - alpha, this->Ms_[ik_ibz].at(isym).data(), i1, i1, pv.desc, DMkibz_M.data(), i1, i1, pv.desc, - beta, DMk.data(), i1, i1, pv.desc); - } - else - { - // Physical DM rotation D(k) = M^dagger D(k_ibz) M, with M = T (x) U is the anti-homomorphism rep in row-major convention. - // ABACUS stores the DM transposed (S = D^T), for which this becomes S(gk) = M^T S(k_ibz) M^* = (conj M)^dagger S (conj M) - // For nspin<4 the orbital-only M is real, so Mc = M and this is bit-identical to the old M^dagger D M. - const std::vector>& Mref = this->Ms_[ik_ibz].at(isym); - std::vector> Mc(Mref.size()); - for (size_t i = 0; i < Mref.size(); ++i) { Mc[i] = std::conj(Mref[i]); } - ScalapackConnector::gemm(dagger, notrans, nbasis, nbasis, nbasis, - alpha, Mc.data(), i1, i1, pv.desc, DMkibz.data(), i1, i1, pv.desc, - beta, DMkibz_M.data(), i1, i1, pv.desc); - alpha.real(1.0 / static_cast(kstar_size)); - ScalapackConnector::gemm(notrans, notrans, nbasis, nbasis, nbasis, - alpha, DMkibz_M.data(), i1, i1, pv.desc, Mc.data(), i1, i1, pv.desc, - beta, DMk.data(), i1, i1, pv.desc); - } - return DMk; - } - - std::vector> Symmetry_rotation::set_sigma_y_2d(const Parallel_2D& pv) const - { - std::vector> sigma_y(pv.get_local_size(), 0.0); - const int nlocal = pv.get_global_row_size(); // = 2*nao for nspin=4 - // sigma_y = [[0, -i], [i, 0]] on the interleaved spin index (I = 2*iorb + spin) - const std::complex sy[2][2] = { {std::complex(0.0, 0.0), std::complex(0.0, -1.0)}, - {std::complex(0.0, 1.0), std::complex(0.0, 0.0)} }; - for (int iorb = 0; 2 * iorb < nlocal; ++iorb) - { - for (int a = 0; a < 2; ++a) - { - const int b = 1 - a; // only the off-diagonal spin entries are non-zero - const int gi = 2 * iorb + a; - const int gj = 2 * iorb + b; - if (pv.in_this_processor(gi, gj)) - { - const int index = pv.global2local_col(gj) * pv.get_row_size() + pv.global2local_row(gi); - sigma_y[index] = sy[a][b]; - } - } - } - return sigma_y; - } - - std::vector> Symmetry_rotation::trs_spin_rotate(const std::vector>& X, - const std::vector>& sigma_y, const Parallel_2D& pv, const double scale) const - { - // stored (transposed 2d-block) form of D_new = sigma_y * conj(D) * sigma_y is - // Sigma_y * conj(X) * Sigma_y (Sigma_y^T = -Sigma_y, the two minus signs cancel). - const char notrans = 'N'; - const int nbasis = pv.get_global_row_size(); - const int i1 = 1; - const std::complex one(1.0, 0.0); - const std::complex beta(0.0, 0.0); - std::vector> Xc(X.size()); - for (size_t i = 0; i < X.size(); ++i) { Xc[i] = std::conj(X[i]); } - std::vector> tmp(pv.get_local_size(), 0.0); - std::vector> out(pv.get_local_size(), 0.0); - // tmp = Sigma_y * conj(X) - ScalapackConnector::gemm(notrans, notrans, nbasis, nbasis, nbasis, - one, sigma_y.data(), i1, i1, pv.desc, Xc.data(), i1, i1, pv.desc, - beta, tmp.data(), i1, i1, pv.desc); - // out = scale * tmp * Sigma_y - ScalapackConnector::gemm(notrans, notrans, nbasis, nbasis, nbasis, - std::complex(scale, 0.0), tmp.data(), i1, i1, pv.desc, sigma_y.data(), i1, i1, pv.desc, - beta, out.data(), i1, i1, pv.desc); - return out; - } std::vector Symmetry_rotation::get_Rs_from_adjacent_list(const UnitCell& ucell, const Grid_Driver& gd, diff --git a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h index 1d27edfef2f..8c830c15974 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h +++ b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h @@ -1,119 +1,28 @@ #pragma once -#include "irreducible_sector.h" +#include "source_cell/module_symmetry/symm_rotation_k.h" #include "source_basis/module_ao/parallel_orbitals.h" #include #include "source_hamilt/module_hcontainer/hcontainer.h" #include "source_cell/module_neighbor/sltk_grid_driver.h" -#include "source_cell/module_symmetry/symm_rot_spin.h" namespace ModuleSymmetry { - using Tap = std::pair; - using TC = std::array; - using TapR = std::pair; - using TCdouble = Abfs::Vector3_Order; - - class Symmetry_rotation + /// Real-space (RI::Tensor / HContainer) H(R) and RI-coefficient symmetry restoration for + /// EXX/RPA, built on top of the LibRI-independent k-space restoration in + /// ModuleSymmetry::Symmetry_rotation_k (source_cell/module_symmetry/symm_rotation_k.h), + /// which provides cal_Ms/restore_dm/rot_matrix_ao and the shared rotation-matrix machinery + /// (rotmat_Slm_, irs_, Ms_, spin_U_, ...). Only the parts that genuinely need LibRI (RI::Tensor + /// atom-pair maps, HContainer real-space rotation) live here. + class Symmetry_rotation : public Symmetry_rotation_k { public: Symmetry_rotation() {}; ~Symmetry_rotation() {}; - //-------------------------------------------------------------------------------- - // getters - const std::map>& 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); - } - /// the rotation matrix under the basis of S_l^m. size: [nsym][lmax][nm*nm] - const std::vector>>>& 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& 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>& abfs_l_nchi); //-------------------------------------------------------------------------------- - /// functions to contruct rotation matrix in AO-representation - - /// The top-level calculation interface of this class. calculate the rotation matrix in AO representation: M - /// only need once call in each ion step (decided by the configuration) - /// @param kstars equal k points to each ibz-kpont, corresponding to a certain symmetry operations. - void cal_Ms(const K_Vectors& kv, - //const std::vector>& kstars, - const UnitCell& ucell, const Parallel_2D& pv); - - /// Use calculated M matrix to recover D(k) from D(k_ibz): D(k) = M(R, k)^\dagger D(k_ibz) M(R, k) - /// the link "ik_ibz-isym-ik" can be found in kstars: k_bz = gmat[isym](k) - std::vector>>restore_dm(const K_Vectors& kv, - const std::vector>>& dm_k_ibz, - const Parallel_2D& pv)const; - std::vector>restore_dm(const K_Vectors& kv, - const std::vector>& dm_k_ibz, - const Parallel_2D& pv)const; - std::vector> rot_matrix_ao(const std::vector>& DMkibz, - const int ik_ibz, const int kstar_size, const int isym, const Parallel_2D& pv, const bool TRS_conj = false) const; - - /// (nspin=4) build the 2*nao spin operator Sigma_y = I_nao (x) sigma_y in 2d-block layout. - std::vector> set_sigma_y_2d(const Parallel_2D& pv) const; - - /// (nspin=4) time-reversal on the spin density matrix: D(k) = sigma_y D^*(-k) sigma_y, - /// realized distribution-safely as scale * Sigma_y * conj(X) * Sigma_y (X is the already - /// space-group-rotated D(-k) stored in the transposed 2d-block convention). - std::vector> trs_spin_rotate(const std::vector>& X, - const std::vector>& sigma_y, const Parallel_2D& pv, const double scale) const; - - /// Inject synthetic AO rotations for density-restoration regression tests. - void set_density_rotations_for_testing( - const std::vector>>>& rotations, - const std::vector>& little_groups, - const int nrot) - { - this->Ms_ = rotations; - this->little_groups_ = little_groups; - this->nsym_ = nrot; - } - - /// calculate Wigner D matrix - double wigner_d(const double beta, const int l, const int m1, const int m2) const; - std::complex wigner_D(const TCdouble& euler_angle, const int l, const int m1, const int m2, const bool inv) const; - - /// c^l_{m1, m2}= - std::complex ovlp_Ylm_Slm(const int l, const int m1, const int m2) const; - - /// calculate euler angle from rotation matrix - TCdouble get_euler_angle(const ModuleBase::Matrix3& gmatc) const; - - /// T_mm' = [c^\dagger D c]_mm', the rotation matrix in the representation of real sphere harmonics - /// @param nop number of operations in gmatc; <0 means nsym_ (the unitary ones only). - /// Pass nsym_+nanti_ to also build the antiunitary operations' T_l. - void cal_rotmat_Slm(const ModuleBase::Matrix3* gmatc, const int lmax, const int nop); - - /// set a block matrix onto a 2d-parallelized matrix(col-maj), at the position (starti, startj) - /// if trans=true, the block matrix is transposed before setting - void set_block_to_mat2d(const int starti, const int startj, const RI::Tensor>& block, - std::vector>& 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>& block, - std::vector& obj_mat, const Parallel_2D& pv, const bool trans = false) const; - - /// 2d-block parallized rotation matrix in AO-representation, denoted as M. - /// finally we will use D(k)=M(R, k)^\dagger*D(Rk)*M(R, k) to recover D(k) from D(Rk). - std::vector> 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 SpinRotation::Su2& spin_U /*= SpinRotation::Su2{ 1.0, 0.0, 0.0, 1.0 }*/) const; - - std::vector>>>& get_rotmat_Slm() { return this->rotmat_Slm_; } //-------------------------------------------------------------------------------- /// The main functions to rotate matrices @@ -145,7 +54,7 @@ namespace ModuleSymmetry template // 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, RI::Tensor>>& HR_full); - template // test the rotation of RI coefficients + template // test the rotation of RI coefficients void test_Cs_rotation(const Symmetry& symm, const Atom* atoms, const Statistics& st, const std::map, RI::Tensor>>& Cs_full)const; template // HContainer type, using row-major implementation @@ -185,53 +94,19 @@ namespace ModuleSymmetry RI::Tensor set_rotation_matrix(const Atom& a, const int& isym)const; template RI::Tensor set_rotation_matrix_abf(const int& type, const int& isym)const; - //-------------------------------------------------------------------------------- - int nsym_ = 1; - /// (nspin=4, magnetic) number of ANTIUNITARY elements Theta*g of the Shubnikov group. - /// Their orbital rotations / return lattices / Ms are appended after the nsym_ unitary - /// ones, so the raw index isym in [nsym_, nsym_+nanti_) addresses gmatrix_anti[isym-nsym_]. - int nanti_ = 0; - /// (nspin=4) true when the configuration carries a non-zero local moment. Then pure time - /// reversal is NOT a symmetry (it reverses m) and the k-star must be restored with the - /// Shubnikov elements Theta*gmatrix_anti[] instead of the generic -k shortcut. - bool magnetic_nspin4_ = false; - - double eps_ = 1e-6; - - // (removed, not needed) TRS_first_: - // it used to short-circuit any star member equal to -k to pure time reversal, - // which silently pre-empted the genuine space-group operation that produced it. - // The operation is now decided by the index alone: isym=nsym_ antiunitary. - // A -k member reached through the TRS doubling lands on the antiunitary branch with M=I, - // which reduces exactly to the direct conjugation. + /// RI::Tensor mirror of rotmat_Slm_ (which is stored as ModuleBase::ComplexMatrix, shared + /// with the LibRI-free k-space code), rebuilt lazily and cached across the many + /// set_rotation_matrix/set_rotation_matrix_abf calls within one ion step (one per atom + /// pair/cell), instead of reconverting the same small block every time. + const RI::Tensor>& get_rotmat_Slm_tensor(const int isym, const int l)const; + mutable std::vector>>> rotmat_Slm_tensor_; + mutable int rotmat_Slm_tensor_version_ = -1; + //-------------------------------------------------------------------------------- bool reduce_Cs_ = false; - int abfs_Lmax_ = 0; std::vector> abfs_l_nchi_;///< number of abfs for each angular momentum - - /// the rotation matrix under the basis of S_l^m. size: [nsym][lmax][nm*nm] - std::vector>>> rotmat_Slm_; - // [natom][nsym], phase factor corresponding to a certain kvec_d_ibz - // std::vector>> phase_factor_; - - /// The unitary matrix associate D(Rk) with D(k) for each ibz-kpoint Rk and each symmetry operation. - /// size: [nks_ibz][nsym][nbasis*nbasis], only need to calculate once. - std::vector>>> Ms_; - - /// Unitary operations fixing each IBZ k point modulo reciprocal lattice vectors. - /// Geometry data built with Ms_ in cal_Ms, not an SCF workflow switch. - std::vector> little_groups_; - - /// (nspin=4) the SU(2) spin-1/2 rotation U(isym) for each symmetry operation, size [nsym]. - /// The spinor AO rotation is T(isym) (x) U(isym); restore_HR_nspin4 uses it to mix the 4 spin - /// channels of the real-space EXX H(R). Filled in cal_Ms (identity for nspin<4). - std::vector spin_U_; - - /// irreducible sector - Irreducible_Sector irs_; - }; template std::string vec3_fmt(const T& x, const T& y, const T& z) @@ -250,4 +125,4 @@ namespace ModuleSymmetry } #include "symm_rotation_r.hpp" -#include "symm_rotation_r_hcontainer.hpp" \ No newline at end of file +#include "symm_rotation_r_hcontainer.hpp" diff --git a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r.hpp b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r.hpp index 6c7cd64f8a3..381aaf5f435 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r.hpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r.hpp @@ -6,6 +6,38 @@ #include namespace ModuleSymmetry { + /// rotmat_Slm_ (inherited from Symmetry_rotation_k, shared with the LibRI-free k-space + /// restoration) is stored as ModuleBase::ComplexMatrix; convert to RI::Tensor here, at the + /// boundary of the real-space (RI::Tensor atom-pair) rotation code that needs that type. + inline RI::Tensor> rotmat_Slm_to_tensor(const ModuleBase::ComplexMatrix& block) + { + RI::Tensor> t({ static_cast(block.nr), static_cast(block.nc) }); + for (int i = 0; i < block.nr; ++i) { + for (int j = 0; j < block.nc; ++j) { + t(i, j) = block(i, j); + } } + return t; + } + + /// Cached (isym, l) lookup into rotmat_Slm_tensor_: rebuilt only when rotmat_Slm_version_ + /// (bumped by cal_rotmat_Slm, i.e. once per ion step) changes, instead of reconverting the + /// same block on every set_rotation_matrix/set_rotation_matrix_abf call (once per atom pair). + inline const RI::Tensor>& Symmetry_rotation::get_rotmat_Slm_tensor(const int isym, const int l)const + { + if (this->rotmat_Slm_tensor_version_ != this->rotmat_Slm_version_) + { + this->rotmat_Slm_tensor_.resize(this->rotmat_Slm_.size()); + for (size_t is = 0; is < this->rotmat_Slm_.size(); ++is) + { + this->rotmat_Slm_tensor_[is].resize(this->rotmat_Slm_[is].size()); + for (size_t il = 0; il < this->rotmat_Slm_[is].size(); ++il) { + this->rotmat_Slm_tensor_[is][il] = rotmat_Slm_to_tensor(this->rotmat_Slm_[is][il]); + } } + this->rotmat_Slm_tensor_version_ = this->rotmat_Slm_version_; + } + return this->rotmat_Slm_tensor_[isym][l]; + } + /// Elementwise complex conjugation used by the time-reversal branch of restore_HR_nspin4. /// Overloaded (not specialized) so a real Tdata compiles to the identity. inline float conj_elem(const float v) { return v; } @@ -268,7 +300,7 @@ namespace ModuleSymmetry { int l = a.iw2l[iw]; int nm = 2 * l + 1; - set_block(iw, iw, this->rotmat_Slm_[isym][l], T); + set_block(iw, iw, this->get_rotmat_Slm_tensor(isym, l), T); iw += nm; } return T; @@ -328,7 +360,7 @@ namespace ModuleSymmetry int nm = 2 * L + 1; for (int N = 0;N < this->abfs_l_nchi_[type][L];++N) { - set_block(iw, iw, this->rotmat_Slm_[isym][L], T); + set_block(iw, iw, this->get_rotmat_Slm_tensor(isym, L), T); iw += nm; // std::cout << "L=" << L << ", N=" << N << ", iw=" << iw << "\n"; } diff --git a/source/source_lcao/module_ri/module_exx_symmetry/test/CMakeLists.txt b/source/source_lcao/module_ri/module_exx_symmetry/test/CMakeLists.txt index ba45e581124..46fe581b419 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/test/CMakeLists.txt +++ b/source/source_lcao/module_ri/module_exx_symmetry/test/CMakeLists.txt @@ -4,7 +4,7 @@ abacus_disable_feature_definitions(__ROCM) AddTest( TARGET MODULE_RI_EXX_SYMMETRY_rotation LIBS base device symmetry neighbor parameter - SOURCES symm_rotation_test.cpp test_symm_rotation.cpp ../symm_rotation.cpp ../symm_rot_out.cpp ../irreducible_sector.cpp ../irred_sec_bvk.cpp + SOURCES symm_rotation_test.cpp test_symm_rotation.cpp ../symm_rotation.cpp ../symm_rot_out.cpp ../../../../source_basis/module_ao/parallel_orbitals.cpp ) \ No newline at end of file diff --git a/source/source_lcao/module_ri/module_exx_symmetry/test/symm_rotation_test.cpp b/source/source_lcao/module_ri/module_exx_symmetry/test/symm_rotation_test.cpp index e3dc3ef3db5..483df5d894a 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/test/symm_rotation_test.cpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/test/symm_rotation_test.cpp @@ -116,7 +116,7 @@ TEST_F(SymmetryRotationTest, OvlpYS) TEST_F(SymmetryRotationTest, RotMat) { symrot.cal_rotmat_Slm(&C41, 1, -1); - RI::Tensor>& rotmat = symrot.get_rotmat_Slm()[0][1]; + ModuleBase::ComplexMatrix& rotmat = symrot.get_rotmat_Slm()[0][1]; int l = 1; for (int m1 = -l;m1 <= l;++m1) for (int m2 = -l;m2 <= l;++m2) @@ -145,7 +145,7 @@ TEST_F(SymmetryRotationTest, SetBlockToMat2d) for (int j = 0;j < pv.get_col_size();++j) for (int i = 0;i < pv.get_row_size();++i) obj_mat[j * pv.get_row_size() + i] = std::complex(static_cast(pv.local2global_row(i)), static_cast(pv.local2global_col(j))); - RI::Tensor> block({ 2, 2 }); + ModuleBase::ComplexMatrix block(2, 2); block(0, 0) = 0; block(0, 1) = -1; block(1, 0) = -2; block(1, 1) = -3; symrot.set_block_to_mat2d(2, 3, block, obj_mat, pv); for (int i = 2;i < 4;++i) diff --git a/source/source_lcao/module_ri/module_exx_symmetry/test/test_symm_rotation.cpp b/source/source_lcao/module_ri/module_exx_symmetry/test/test_symm_rotation.cpp index 01ec4938bfb..64b8ad020ca 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/test/test_symm_rotation.cpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/test/test_symm_rotation.cpp @@ -1,22 +1,6 @@ #include "../symm_rotation.h" -#include "source_io/module_parameter/parameter.h" #include "gtest/gtest.h" -class TestParameters -{ - public: - TestParameters(Parameter& parameters, const int nspin) - : parameters_(parameters), original_nspin_(parameters.inp.nspin) - { - parameters_.input.nspin = nspin; - } - ~TestParameters() { parameters_.input.nspin = original_nspin_; } - - private: - Parameter& parameters_; - const int original_nspin_; -}; - // K-point generation is outside this test: use explicit stars, but provide // the virtual symbols needed by the existing lightweight rotation test target. void ModuleCell::ReciprocalGrid::renew(const int&) @@ -62,7 +46,6 @@ std::vector rotate_reference(const std::vector& density, void check_little_group_restoration(const int nspin) { - const TestParameters parameters(PARAM, nspin); const int channels = nspin == 2 ? 2 : 1; const int n = 4; Parallel_2D pv; @@ -98,10 +81,15 @@ void check_little_group_restoration(const int nspin) }; rotation.set_density_rotations_for_testing( {{{0, local(identity)}, {1, local(little)}, {2, local(representative)}, {3, local(alternate)}}}, - {{0, 1}}, 4); + {{0, 1}}, 4, nspin); K_Vectors kv; kv.set_nkstot(channels); kv.set_nkstot_nospin(2); + // Single-pool (KPAR=1) scenario: kv.get_nks() (local) equals kv.get_nkstot() (global), + // and ik2iktot is the identity map. With only one global ibz-k here, any value mod + // kv.kstars.size()==1 is 0, so the exact ik2iktot values don't matter, only its size. + kv.set_nks(channels); + kv.ik2iktot.assign(channels, 0); kv.kstars = {{{0, {0.25, 0.0, 0.0}}, {2, {0.0, 0.25, 0.0}}}}; std::vector> inputs; std::vector> expected; diff --git a/source/source_lcao/module_ri/rpa_lri.hpp b/source/source_lcao/module_ri/rpa_lri.hpp index b67d8163061..7a776992541 100644 --- a/source/source_lcao/module_ri/rpa_lri.hpp +++ b/source/source_lcao/module_ri/rpa_lri.hpp @@ -160,13 +160,13 @@ void RPA_LRI::cal_postSCF_exx(const elecstate::DensityMatrix { const std::array period = RI_Util::get_Born_vonKarmen_period(kv); const auto& Rs = RI_Util::get_Born_von_Karmen_cells(period); - symrot.find_irreducible_sector(ucell.symm, ucell.atoms, ucell.st, Rs, period, ucell.lat); + symrot.find_irreducible_sector(ucell.symm, ucell.atoms, ucell.st, Rs, period, ucell.lat, PARAM.globalv.global_out_dir); // set Lmax of the rotation matrices to max(l_ao, l_abf), to support rotation under ABF // NOTE: Using Exx_Abfs::Construct_Orbs::get_Lmax() to compute Lmax from the actual ABFs // instead of relying on exx_cut_coulomb->abfs_Lmax() (not yet initialized) or // this->info.abfs_Lmax (defaults to 0). This ensures correct Lmax for symmetry rotation. symrot.set_abfs_Lmax(Exx_Abfs::Construct_Orbs::get_Lmax(abfs_for_lmax)); - symrot.cal_Ms(kv, ucell, *dm.get_paraV_pointer()); + symrot.cal_Ms(kv, ucell, *dm.get_paraV_pointer(), PARAM.inp.nspin); // output Ts (symrot_R.txt) and Ms (symrot_k.txt) ModuleSymmetry::print_symrot_info_R(symrot, ucell.symm, ucell.lmax, Rs); ModuleSymmetry::print_symrot_info_k(symrot, kv, ucell); diff --git a/tests/17_DS_DFTU/03_LCAO_DFTU_S2_Z/INPUT b/tests/17_DS_DFTU/03_LCAO_DFTU_S2_Z/INPUT index 1eb50a84479..edd5eda7ef8 100644 --- a/tests/17_DS_DFTU/03_LCAO_DFTU_S2_Z/INPUT +++ b/tests/17_DS_DFTU/03_LCAO_DFTU_S2_Z/INPUT @@ -15,7 +15,8 @@ smearing_sigma 0.01 mixing_type broyden mixing_beta 0.4 ks_solver genelpa -symmetry 0 +symmetry 1 +kpar 4 # DFT+U parameters dft_plus_u 1 diff --git a/tests/17_DS_DFTU/03_LCAO_DFTU_S2_Z/result.ref b/tests/17_DS_DFTU/03_LCAO_DFTU_S2_Z/result.ref index bb1d7aa5ce8..2c12a092ecd 100644 --- a/tests/17_DS_DFTU/03_LCAO_DFTU_S2_Z/result.ref +++ b/tests/17_DS_DFTU/03_LCAO_DFTU_S2_Z/result.ref @@ -1,3 +1,6 @@ -etotref -6772.0999515218118177 -etotperatomref -3386.0499757609 +etotref -6771.6902262248995612 +etotperatomref -3385.8451131124 +pointgroupref C_3v +spacegroupref C_3v +nksibzref 2 totaltimeref 6.11 diff --git a/tests/17_DS_DFTU/04_LCAO_DFTU_S4_XY/INPUT b/tests/17_DS_DFTU/04_LCAO_DFTU_S4_XY/INPUT index a0de7786ca2..d1a025a67b7 100644 --- a/tests/17_DS_DFTU/04_LCAO_DFTU_S4_XY/INPUT +++ b/tests/17_DS_DFTU/04_LCAO_DFTU_S4_XY/INPUT @@ -16,7 +16,8 @@ smearing_sigma 0.01 mixing_type broyden mixing_beta 0.4 ks_solver genelpa -symmetry 0 +symmetry 1 +kpar 4 # DFT+U parameters dft_plus_u 1 diff --git a/tests/17_DS_DFTU/04_LCAO_DFTU_S4_XY/result.ref b/tests/17_DS_DFTU/04_LCAO_DFTU_S4_XY/result.ref index 1a258220cfb..3f7ce113b38 100644 --- a/tests/17_DS_DFTU/04_LCAO_DFTU_S4_XY/result.ref +++ b/tests/17_DS_DFTU/04_LCAO_DFTU_S4_XY/result.ref @@ -1,3 +1,7 @@ -etotref -6772.1004497577005168 -etotperatomref -3386.0500835053 +etotref -6771.5947932947956360 +etotperatomref -3385.7973966474 +pointgroupref C_3v +spacegroupref D_3d +nksibzref 4 +magpointgroupref C_2 totaltimeref 10.33