Feature: Symmetry for DFT+U

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\

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

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)

source/source_lcao/module_ri/module_exx_symmetry/irred_sec_bvk.cpp → 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;
     // }
 
-};
+};

source/source_lcao/module_ri/module_exx_symmetry/irreducible_sector.cpp → 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<TCdouble>(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<TC>& 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<TC>& 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);
     }
 }

source/source_lcao/module_ri/module_exx_symmetry/irreducible_sector.h → source/source_cell/module_symmetry/irreducible_sector.h

@@ -2,7 +2,7 @@
 #include <vector>
 #include <map>
 #include <set>
-#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<int, int>;
     using TC = std::array<int, 3>;
     using TapR = std::pair<Tap, TC>;
-    using TCdouble = Abfs::Vector3_Order<double>;
+    using TCdouble = ModuleBase::Vector3<double>;
 
     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
+        ///                    <output_dir>/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<TC>& Rs, const TC& period, const Lattice& lat);
+            const std::vector<TC>& Rs, const TC& period, const Lattice& lat, const std::string& output_dir = "");
         const std::map<Tap, std::set<TC>>& get_irreducible_sector()const { return this->irreducible_sector_; }
         // const std::map<int, std::set<std::pair<int, TC>>> 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;
     };
-}
+}