feat(exx-symmetry): SSG spin-flip coset support for H(R) rotation (nspin=2)

Author: maki49

source/source_io/module_parameter/read_inp_sys.cpp

@@ -278,7 +278,7 @@ void ReadInput::item_system()
                           "* 0: disabled (default)\n"
                           "* 1: enabled";
         item.default_value = "0";
-        item.availability = "symmetry==1";
+        item.set_availability("symmetry==1");
         read_sync_int(input.symmetry_ssg);
         this->add_item(item);
     }

source/source_lcao/module_ri/exx_lri.hpp

@@ -803,14 +803,20 @@ void Exx_LRI<Tdata>::cal_exx_elec(const std::vector<std::map<TA, std::map<TAC, R
 		return;
 	}
 
+	// (nspin=2 SSG) the up-sublattice-to-down-sublattice spin-flip coset makes the down channel's
+	// H(R) a single flip-op rotation of the up channel's, so we compute cal_Hs for is=0 only and
+	// derive is=1 from the stored up-channel full H(R); see restore_HR_flip_nspin2.
+	const bool ssg_flip = (p_symrot && PARAM.inp.nspin == 2 && ucell.symm.spin_flip_nspin2);
+	std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> Hs_full_up;
+
 	this->Hexxs.resize(PARAM.inp.nspin);
 	this->Eexx = 0;
 	for(int is=0; is<PARAM.inp.nspin; ++is)
 	{
 		const std::string suffix = ((PARAM.inp.cal_force || PARAM.inp.cal_stress) ? std::to_string(is) : "");
 
 		this->exx_lri.set_Ds(Ds[is], this->info.dm_threshold, suffix);
-		this->exx_lri.cal_Hs({ "","",suffix });
+		if (!(ssg_flip && is == 1)) { this->exx_lri.cal_Hs({ "","",suffix }); }
 
 		if (!p_symrot)
 		{
@@ -819,10 +825,20 @@ void Exx_LRI<Tdata>::cal_exx_elec(const std::vector<std::map<TA, std::map<TAC, R
 		}
 		else
 		{
-			// reduce but not repeat
-			auto Hs_a2D = this->exx_lri.post_2D.set_tensors_map2(this->exx_lri.Hs);
-			// rotate locally without repeat
-			Hs_a2D = p_symrot->restore_HR(ucell.symm, ucell.atoms, ucell.st, 'H', Hs_a2D);
+			std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> Hs_a2D;
+			if (ssg_flip && is == 1)
+			{
+				// down channel: derive the full H(R) from the stored up channel by one spin-flip op
+				Hs_a2D = p_symrot->restore_HR_flip_nspin2(ucell.symm, ucell.atoms, ucell.st, 'H', Hs_full_up);
+			}
+			else
+			{
+				// reduce but not repeat
+				Hs_a2D = this->exx_lri.post_2D.set_tensors_map2(this->exx_lri.Hs);
+				// rotate locally without repeat
+				Hs_a2D = p_symrot->restore_HR(ucell.symm, ucell.atoms, ucell.st, 'H', Hs_a2D);
+				if (ssg_flip && is == 0) { Hs_full_up = Hs_a2D; }   // keep up channel for the flip derivation
+			}
 			// cal energy using full Hs without repeat
 			this->exx_lri.energy = this->exx_lri.post_2D.cal_energy(
 				this->exx_lri.post_2D.saves["Ds_" + suffix],

source/source_lcao/module_ri/module_exx_symmetry/irreducible_sector.cpp

@@ -9,12 +9,17 @@ namespace ModuleSymmetry
     // 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];
+        if (isym < symm.nrotk) { return symm.gmatrix[isym]; }
+        // (nspin=2 SSG) the spin-flip coset reuses the [nrotk, nrotk+nrotk_flip) index range
+        // (nrotk_anti==0 in that case); its spatial part is gmatrix_flip[isym-nrotk].
+        return symm.spin_flip_nspin2 ? symm.gmatrix_flip[isym - symm.nrotk]
+                                     : 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);
+        if (isym < symm.nrotk) { return symm.get_rotated_atom(isym, iat); }
+        return symm.spin_flip_nspin2 ? symm.get_rotated_atom_flip(isym - symm.nrotk, iat)
+                                     : symm.get_rotated_atom_anti(isym - symm.nrotk, iat);
     }
 
     TC Irreducible_Sector::rotate_R(const Symmetry& symm,

source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h

@@ -127,6 +127,16 @@ namespace ModuleSymmetry
         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;
+        /// (nspin=2 SSG) derive the FULL H(R) of the OPPOSITE spin channel from a full H(R), using a
+        /// single spin-flip coset operation f=[C2_perp||g]. With SOC off the spin channels are
+        /// independent real matrices and f's spatial/orbital rotation is identical to a unitary op,
+        /// so H_down(f.apR)=rotate_f(H_up(apR)); f is a bijection on the full atom-pair set, so one
+        /// forward pass over HR_full_this yields the whole opposite channel. This lets cal_exx_elec
+        /// skip the second spin's cal_Hs entirely. Uses flip op j=0 (raw isym = nrotk).
+        template<typename Tdata>    // RI::Tensor type
+        std::map<int, std::map<std::pair<int, TC>, RI::Tensor<Tdata>>> restore_HR_flip_nspin2(
+            const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode,
+            const std::map<int, std::map<std::pair<int, TC>, RI::Tensor<Tdata>>>& HR_full_this)const;
 
         //--------------------------------------------------------------------------------
         /// test functions

source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r.hpp

@@ -248,6 +248,48 @@ namespace ModuleSymmetry
         return HR_full;
     }
 
+    // (nspin=2 SSG) derive the opposite spin channel's FULL H(R) from a full H(R) by one spin-flip
+    // coset op f=[C2_perp||g] (raw isym = nrotk, i.e. gmatrix_flip[0]). With SOC off the two spin
+    // channels are independent real matrices, and f's spatial/orbital rotation is identical to a
+    // unitary op (only the target spin channel differs), so the covariance is exactly the one used
+    // by restore_HR:  H_other[Z] = rotate_f( H_this[ f.Z ] ),  with f.Z = rotate_apR_by_formula(f, Z).
+    // f is a bijection on the full atom-pair set, so a single forward pass fills the whole channel.
+    template<typename Tdata>
+    std::map<int, std::map<std::pair<int, TC>, RI::Tensor<Tdata>>> Symmetry_rotation::restore_HR_flip_nspin2(
+        const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode,
+        const std::map<int, std::map<std::pair<int, TC>, RI::Tensor<Tdata>>>& HR_full_this) const
+    {
+        ModuleBase::TITLE("Symmetry_rotation", "restore_HR_flip_nspin2");
+        ModuleBase::timer::start("Symmetry_rotation", "restore_HR_flip_nspin2");
+        assert(symm.spin_flip_nspin2);
+        assert(symm.nrotk_flip > 0);
+        const int isym_flip = symm.nrotk;   // flip op j=0, raw sector index nrotk (nrotk_anti==0 here)
+        std::map<int, std::map<std::pair<int, TC>, RI::Tensor<Tdata>>> HR_other;
+        for (auto& tmp1 : HR_full_this)
+        {
+            const int& z1 = tmp1.first;
+            for (auto& tmp2 : tmp1.second)
+            {
+                const int& z2 = tmp2.first.first;
+                const TC& Rz = tmp2.first.second;
+                const TapR& src = this->irs_.rotate_apR_by_formula(symm, isym_flip, { { z1, z2 }, Rz });
+                const int& s1 = src.first.first;
+                const int& s2 = src.first.second;
+                const TC& Rs = src.second;
+                // f is a bijection on the full apR set, so H_this[src] should exist; a missing entry
+                // can only be a below-threshold drop, treated as zero (skip).
+                auto it1 = HR_full_this.find(s1);
+                if (it1 == HR_full_this.end()) { continue; }
+                auto it2 = it1->second.find({ s2, Rs });
+                if (it2 == it1->second.end()) { continue; }
+                HR_other[z1][{z2, Rz}] = rotate_atompair_serial(it2->second, isym_flip,
+                    atoms[st.iat2it[s1]], atoms[st.iat2it[s2]], mode);
+            }
+        }
+        ModuleBase::timer::end("Symmetry_rotation", "restore_HR_flip_nspin2");
+        return HR_other;
+    }
+
     template<typename Tdata>
     inline void set_block(const int starti, const int startj, const RI::Tensor<std::complex<double>>& block,
         RI::Tensor<Tdata>& obj_tensor)