Perf: antiunitary ops are also used to reduce Hexx(R)

Author: maki49

source/source_io/module_parameter/input_conv.cpp

@@ -509,7 +509,7 @@ void Input_Conv::Convert()
         GlobalC::exx_info.info_opt_abfs.tolerence = PARAM.inp.exx_opt_orb_tolerence;
 
         // Space-group symmetry is supported for LCAO EXX (nspin=1,2 via restore_dm/restore_HR;
-        // nspin=4/SOC via restore_dm + restore_HR_soc), so symmetry=1 is honored here.
+        // nspin=4/SOC via restore_dm + restore_HR_nspin4), so symmetry=1 is honored here.
 
         GlobalC::exx_info.sync_from_global();
     }

source/source_lcao/module_ri/Exx_LRI.h

@@ -93,7 +93,7 @@ class Exx_LRI
 		const ModuleSymmetry::Symmetry_rotation* p_symrot = nullptr);
 	// (nspin=4) real-space symmetry EXX: the spinor H(R) rotation couples the 4 spin channels via
 	// the SU(2) part U(isym), so the 4 channels must be rotated together (not one-per-outer-loop).
-	// Gathers the irreducible Hs of all 4 channels, calls Symmetry_rotation::restore_HR_soc, then
+	// Gathers the irreducible Hs of all 4 channels, calls Symmetry_rotation::restore_HR_nspin4, then
 	// finishes energy/gather per channel. Called from cal_exx_elec when p_symrot && nspin==4.
 	void cal_exx_elec_soc(
 		const std::vector<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>>& Ds,

source/source_lcao/module_ri/Exx_LRI.hpp

@@ -859,7 +859,7 @@ void Exx_LRI<Tdata>::cal_exx_elec_soc(
 
 	// pass 2: spinor-coupled rotation of the 4 channels from the irreducible sector to the full BZ
 	std::array<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>, 4> Hs_full =
-		p_symrot->restore_HR_soc(ucell.symm, ucell.atoms, ucell.st, 'H', Hs_irr);
+		p_symrot->restore_HR_nspin4(ucell.symm, ucell.atoms, ucell.st, 'H', Hs_irr);
 
 	// pass 3: per-channel energy (full Hs, no repeat), then gather the repeated full Hs for abacus
 	for (int is = 0; is < 4; ++is)

source/source_lcao/module_ri/module_exx_symmetry/irreducible_sector.cpp

@@ -2,20 +2,38 @@
 #include "source_io/module_parameter/parameter.h"
 namespace ModuleSymmetry
 {
+    // Raw-index dispatch shared by the real-space sector helpers, matching the convention used
+    // everywhere else (symmetry_rotation.h): isym < nrotk -> unitary gmatrix[isym];
+    // isym >= nrotk -> spatial part of the antiunitary element Theta*gmatrix_anti[isym-nrotk].
+    // Only the SPATIAL part is needed here: Theta acts on H(R) as sigma_y (.)^* sigma_y and
+    // leaves R and the atom pair untouched, so the sector bookkeeping is identical for both kinds.
+    static inline const ModuleBase::Matrix3& sector_gmatrix(const Symmetry& symm, const int isym)
+    {
+        return (isym < symm.nrotk) ? symm.gmatrix[isym] : symm.gmatrix_anti[isym - symm.nrotk];
+    }
+    static inline int sector_rotated_atom(const Symmetry& symm, const int isym, const int iat)
+    {
+        return (isym < symm.nrotk) ? symm.get_rotated_atom(isym, iat)
+                                   : symm.get_rotated_atom_anti(isym - symm.nrotk, iat);
+    }
+
     TC Irreducible_Sector::rotate_R(const Symmetry& symm,
         const int isym, const int iat1, const int iat2, const TC& R, const char gauge) const
     {
         auto round2int = [symm](const double x) -> int { return x > 0 ? static_cast<int>(x + symm.epsilon) : static_cast<int>(x - symm.epsilon); };
         const TCdouble R_double(static_cast<double>(R[0]), static_cast<double>(R[1]), static_cast<double>(R[2]));
+        // return_lattice_ is already sized nrotk+nrotk_anti and indexed by the same raw isym.
+        const ModuleBase::Matrix3& gmat = sector_gmatrix(symm, isym);
         const TCdouble Rrot_double = (gauge == 'L')
-            ? R_double * symm.gmatrix[isym] + this->return_lattice_[iat1][isym] - this->return_lattice_[iat2][isym]
-            : R_double * symm.gmatrix[isym] + this->return_lattice_[iat2][isym] - this->return_lattice_[iat1][isym];
+            ? R_double * gmat + this->return_lattice_[iat1][isym] - this->return_lattice_[iat2][isym]
+            : R_double * gmat + this->return_lattice_[iat2][isym] - this->return_lattice_[iat1][isym];
         return { round2int(Rrot_double.x), round2int(Rrot_double.y), round2int(Rrot_double.z) };
     }
     TapR Irreducible_Sector::rotate_apR_by_formula(const Symmetry& symm,
         const int isym, const TapR& apR, const char gauge) const
     {
-        const Tap& aprot = { symm.get_rotated_atom(isym, apR.first.first), symm.get_rotated_atom(isym, apR.first.second) };
+        const Tap& aprot = { sector_rotated_atom(symm, isym, apR.first.first),
+                             sector_rotated_atom(symm, isym, apR.first.second) };
         return { aprot, this->rotate_R(symm, isym, apR.first.first, apR.first.second, apR.second, gauge) };
     }
 
@@ -183,11 +201,21 @@ namespace ModuleSymmetry
                 for (auto& R : Rs)
                     apR_all[{iat1, iat2}].insert(R);
 
-        // get invmap
+        // get invmap over the operation set actually used by the sector search.
+        // For nspin=4 magnetic that is the full Shubnikov group H (union) A, laid out as
+        // [gmatrix[0..nrotk) | gmatrix_anti[0..nrotk_anti)]. gmatrix_invmap needs no change:
+        // it searches the whole array for s[i]*s[j]==I, and A is closed under inversion
+        // (if g in A had g^-1 in H then g = (g^-1)^-1 would be in H, contradicting H n A = {}),
+        // so the concatenated array is exactly the parent group and every inverse is found,
+        // with the inverse of a coset member landing inside the coset block.
         if (this->invmap_.empty())
         {
-            this->invmap_.resize(symm.nrotk);
-            symm.gmatrix_invmap(symm.gmatrix, symm.nrotk, invmap_.data());
+            const int nop = symm.nrotk + symm.nrotk_anti;
+            std::vector<ModuleBase::Matrix3> gmat_all(nop);
+            for (int i = 0; i < symm.nrotk; ++i) { gmat_all[i] = symm.gmatrix[i]; }
+            for (int j = 0; j < symm.nrotk_anti; ++j) { gmat_all[symm.nrotk + j] = symm.gmatrix_anti[j]; }
+            this->invmap_.resize(nop);
+            symm.gmatrix_invmap(gmat_all.data(), nop, invmap_.data());
         }
 
         // get symmetry of BvK supercell
@@ -213,6 +241,13 @@ namespace ModuleSymmetry
                 //     if (!in_2d_plain[isym]) continue;
                 // }
                 const int& isym = this->isymbvk_to_isym_[isymbvk];
+                // A BvK operation with no counterpart in the unit-cell operation set is marked -1.
+                // For a magnetic nspin=4 system the unit-cell set is the Shubnikov group H (union) A,
+                // which is generally a PROPER subset of the crystallographic group (operations that
+                // merely tilt the moment belong to neither), so an unmatched BvK operation is a
+                // normal outcome, not an error: it is simply not a symmetry of the magnetic system.
+                // Skipping it only costs reduction, never correctness.
+                if (isym < 0) { continue; }
                 const TapR& apRrot = this->rotate_apR_by_formula(symm, this->invmap_[isym], irapR);
                 const Tap& aprot = apRrot.first;
                 const TC& Rrot = apRrot.second;

source/source_lcao/module_ri/module_exx_symmetry/irreducible_sector_bvk.cpp

@@ -25,6 +25,26 @@ namespace ModuleSymmetry
                     break;
                 }
             }
+            // (nspin=4 magnetic) second pass over the antiunitary coset: the spatial part of
+            // Theta*g is a genuine crystallographic operation of the BvK supercell too
+            // (time-reversal factor is only needed when the data is rotated.)
+            // Recorded after unitary ops so downstream (rotate_R, restore_HR_nspin4) can tell the two apart by isym vs. nrotk.
+            if (isymbvk2isym[isymbvk] < 0)
+            {
+                for (int j = 0;j < symm.nrotk_anti;++j)
+                {
+                    if (matequal(bvkgmat[isymbvk], symm.gmatrix_anti[j]))
+                    {
+                        isymbvk2isym[isymbvk] = symm.nrotk + j;
+                        break;
+                    }
+                }
+            }
+            // Unmatched stays -1. That is legitimate for nspin=4 magnetic: the unit-cell set is the
+            // Shubnikov group H (union) A, generally a PROPER subset of the crystallographic group
+            // (operations that merely tilt the moment belong to neither), so a BvK operation may
+            // have no counterpart. The consumer in find_irreducible_sector skips negative entries;
+            // it must never use one as an index.
         }
         return isymbvk2isym;
     }
@@ -42,9 +62,12 @@ namespace ModuleSymmetry
             -> ModuleBase::Matrix3 {return ModuleBase::Matrix3(a1.x, a1.y, a1.z, a2.x, a2.y, a2.z, a3.x, a3.y, a3.z);};
         auto set_bvk_same_as_ucell = [&symm, this]()->void
             {
-                this->bvk_nsym_ = symm.nrotk;
-                this->isymbvk_to_isym_.resize(symm.nrotk);
-                for (int isym = 0;isym < symm.nrotk;++isym) { this->isymbvk_to_isym_[isym] = isym; }
+                // include the antiunitary coset (nspin=4 magnetic); nrotk_anti is 0 otherwise,
+                // so this is unchanged for every other case.
+                const int nop = symm.nrotk + symm.nrotk_anti;
+                this->bvk_nsym_ = nop;
+                this->isymbvk_to_isym_.resize(nop);
+                for (int isym = 0;isym < nop;++isym) { this->isymbvk_to_isym_[isym] = isym; }
             };
         if (bvk_period[0] == bvk_period[1] && bvk_period[0] == bvk_period[2])
         {   //the BvK supercell has the same symmetry as the original cell
@@ -140,8 +163,10 @@ namespace ModuleSymmetry
         bvk_gmatrix.resize(bvk_nsg);
         bvk_gtrans.resize(bvk_nsg);
         this->bvk_nsym_ = bvk_nsg;
-        // bvk suppercell cannot have higher symmetry than the original cell
-        if (this->bvk_nsym_ > symm.nrotk)
+        // bvk suppercell cannot have higher symmetry than the original cell.
+        // The comparison is against the FULL operation set the sector search may use, i.e. the
+        // Shubnikov group H (union) A for nspin=4 magnetic (nrotk_anti is 0 in every other case).
+        if (this->bvk_nsym_ > symm.nrotk + symm.nrotk_anti)
         {
             std::cout << "reset bvk symmetry to the same as the original cell" << std::endl;
             set_bvk_same_as_ucell();

source/source_lcao/module_ri/module_exx_symmetry/symmetry_rotation.cpp

@@ -50,7 +50,7 @@ namespace ModuleSymmetry
         {
             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_soc (real-space EXX H(R) spin mixing)
+        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)
         auto restrict_kpt = [](const TCdouble& kvec, const double& symm_prec) -> TCdouble

source/source_lcao/module_ri/module_exx_symmetry/symmetry_rotation.h

@@ -119,9 +119,12 @@ namespace ModuleSymmetry
         /// orbital rotation T1^dagger(.)T2 (mode 'H') / T1^T(.)T2^* (mode 'D') applied to every
         /// channel, the SU(2) spin part U(isym) mixes them:  H'^{ab}=sum_{cd} conj(U_{ca}) U_{db} [T1^dagger H^{cd} T2].
         /// The 4 channels are ordered is=a*2+b (a=row spin, b=col spin), matching RI_2D_Comm::split_is_block.
-        /// Real-space (atom-pair) reduction uses space-group operations only, so no time-reversal branch is needed.
+        /// (nspin=4 magnetic) The atom-pair reduction may also use the ANTIUNITARY elements of the
+        /// Shubnikov group, flagged by isym >= nsym_. In real space time reversal acts as
+        /// H(R) -> sigma_y H^*(R) sigma_y (R and the orbital indices untouched), which becomes a
+        /// remap of the 4 channels applied after the SU(2) mixing; see symmetry_rotation_R.hpp.
         template<typename Tdata>    // RI::Tensor type
-        std::array<std::map<int, std::map<std::pair<int, TC>, RI::Tensor<Tdata>>>, 4> restore_HR_soc(
+        std::array<std::map<int, std::map<std::pair<int, TC>, RI::Tensor<Tdata>>>, 4> restore_HR_nspin4(
             const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode,
             const std::array<std::map<int, std::map<std::pair<int, TC>, RI::Tensor<Tdata>>>, 4>& HR_irreducible_soc)const;
 
@@ -207,7 +210,7 @@ namespace ModuleSymmetry
         std::vector<std::map<int, std::vector<std::complex<double>>>> Ms_;
 
         /// (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_soc uses it to mix the 4 spin
+        /// 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<SpinRotation::Su2> spin_U_;
 

source/source_lcao/module_ri/module_exx_symmetry/symmetry_rotation_R.hpp

@@ -6,6 +6,13 @@
 #include <RI/physics/symmetry/Symmetry_Rotation.h>
 namespace ModuleSymmetry
 {
+    /// 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; }
+    inline double conj_elem(const double v) { return v; }
+    inline std::complex<float> conj_elem(const std::complex<float>& v) { return std::conj(v); }
+    inline std::complex<double> conj_elem(const std::complex<double>& v) { return std::conj(v); }
+
     template<typename Tdata>
     inline void print_tensor(const RI::Tensor<Tdata>& t, const std::string& name, const double& threshold = 0.0)
     {
@@ -119,14 +126,24 @@ namespace ModuleSymmetry
     // which factorizes into the per-channel orbital rotation G^{cd}=T1^dagger H^{cd} T2 followed by the SU(2) spin mixing
     //     H'^{ab} = sum_{cd} conj(U_{ca}) U_{db} G^{cd}    (mode 'H')
     //     H'^{ab} = sum_{cd}      U_{ca} conj(U_{db}) G^{cd}    (mode 'D')
-    // Only space-group operations enter the atom-pair reduction, so no time-reversal branch is needed.
+    //
+    // (nspin=4 magnetic) The atom-pair reduction may also use the ANTIUNITARY elements Theta*g of
+    // the Shubnikov group, flagged by isym >= nsym_. In real space time reversal acts as
+    //     H(R) -> sigma_y H^*(R) sigma_y
+    // with R and the orbital indices untouched (in a real AO basis Theta = -i sigma_y K, and
+    // H(R) = sum_k H(k) e^{-ikR} turns H(k) -> sigma_y H^*(-k) sigma_y into exactly this).
+    // So the antiunitary case is the unitary result Y followed by one channel remap:
+    //     H'^{00} =  conj(Y^{11}),  H'^{01} = -conj(Y^{10})
+    //     H'^{10} = -conj(Y^{01}),  H'^{11} =  conj(Y^{00})
+    // The map is an involution (sigma_y^* = -sigma_y, sigma_y^2 = I), so no extra sign is needed
+    // when it is applied in either direction.
     template<typename Tdata>
-    std::array<std::map<int, std::map<std::pair<int, TC>, RI::Tensor<Tdata>>>, 4> Symmetry_rotation::restore_HR_soc(
+    std::array<std::map<int, std::map<std::pair<int, TC>, RI::Tensor<Tdata>>>, 4> Symmetry_rotation::restore_HR_nspin4(
         const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode,
         const std::array<std::map<int, std::map<std::pair<int, TC>, RI::Tensor<Tdata>>>, 4>& HR_irreducible_soc)const
     {
-        ModuleBase::TITLE("Symmetry_rotation", "restore_HR_soc");
-        ModuleBase::timer::start("Symmetry_rotation", "restore_HR_soc");
+        ModuleBase::TITLE("Symmetry_rotation", "restore_HR_nspin4");
+        ModuleBase::timer::start("Symmetry_rotation", "restore_HR_nspin4");
         assert(mode == 'H' || mode == 'D');
         std::array<std::map<int, std::map<std::pair<int, TC>, RI::Tensor<Tdata>>>, 4> HR_full;
 
@@ -181,6 +198,7 @@ namespace ModuleSymmetry
                     for (int is = 0;is < 4;++is) { G[is] = this->rotate_atompair_serial(Hir[is], isym, a1, a2, mode); }
                     // step 2: SU(2) spin mixing of the 4 rotated channels into the output channels
                     const SpinRotation::Su2& U = this->spin_U_[isym];
+                    std::array<RI::Tensor<Tdata>, 4> Hout_ch;
                     for (int a = 0;a < 2;++a) {
                         for (int b = 0;b < 2;++b)
                         {
@@ -194,13 +212,39 @@ namespace ModuleSymmetry
                                     Hout += RI::Global_Func::convert<Tdata>(coeff) * G[c * 2 + d];
                                 }
                             }
-                            HR_full[a * 2 + b][ap1][{ap2, R}] = Hout;
+                            Hout_ch[a * 2 + b] = Hout;
+                        }
+                    }
+                    // step 3 (antiunitary elements of the Shubnikov group): apply time reversal
+                    // sigma_y (.)^* sigma_y, i.e. the channel remap documented above.
+                    if (isym >= this->nsym_)
+                    {
+                        // NOTE: antiunitary elements only ever exist for nspin=4 (nrotk_anti is 0 otherwise), where Tdata is complex. 
+                        // conj_elem() is the identity for a
+                        // real Tdata, so this branch must not be reached with one -- it would
+                        // silently degrade into a bare channel swap.
+                        static const int src[4] = { 3, 2, 1, 0 };   // 00<-11, 01<-10, 10<-01, 11<-00
+                        static const bool neg[4] = { false, true, true, false };
+                        std::array<RI::Tensor<Tdata>, 4> Y = Hout_ch;
+                        for (int is = 0;is < 4;++is)
+                        {
+                            const RI::Tensor<Tdata>& s = Y[src[is]];
+                            RI::Tensor<Tdata> t({ static_cast<size_t>(a1.nw), static_cast<size_t>(a2.nw) });
+                            for (size_t i = 0;i < t.shape[0];++i) {
+                                for (size_t j = 0;j < t.shape[1];++j)
+                                {
+                                    const Tdata v = ModuleSymmetry::conj_elem(s(i, j));
+                                    t(i, j) = neg[is] ? -v : v;
+                                }
+                            }
+                            Hout_ch[is] = t;
                         }
                     }
+                    for (int is = 0;is < 4;++is) { HR_full[is][ap1][{ap2, R}] = Hout_ch[is]; }
                 }
             }
         }
-        ModuleBase::timer::end("Symmetry_rotation", "restore_HR_soc");
+        ModuleBase::timer::end("Symmetry_rotation", "restore_HR_nspin4");
         return HR_full;
     }