Feature: output symmetry information for LibRPA (#6435)

* output symmetry information for RPA and GW

* .dat -> .txt

Author: maki49

docs/advanced/input_files/input-main.md

@@ -1978,6 +1978,7 @@ These variables are used to control the output of properties.
 
 - **Type**: Boolean
 - **Description**: Generate output files used in rpa calculations.
+- **Note**: If [`symmetry`](#symmetry) is set to 1, additional files containing the necessary information for exploiting symmetry in the subsequent rpa calculation will be output:  `irreducible_sector.txt`, `symrot_k.txt` and `symrot_R.txt`
 - **Default**: False
 
 ### out_pchg

source/Makefile.Objects

@@ -677,6 +677,7 @@ OBJS_MODULE_RI=conv_coulomb_pot_k.o\
       Mix_DMk_2D.o\
       Mix_Matrix.o\
       symmetry_rotation.o\
+      symmetry_rotation_output.o\
       symmetry_irreducible_sector.o\
 
 OBJS_PARALLEL=parallel_common.o\

source/source_lcao/module_ri/Exx_LRI_interface.hpp

@@ -139,7 +139,7 @@ void Exx_LRI_Interface<T, Tdata>::exx_before_all_runners(
         this->symrot_.find_irreducible_sector(
             ucell.symm, ucell.atoms, ucell.st,
             RI_Util::get_Born_von_Karmen_cells(period), period, ucell.lat);
-        // this->symrot_.set_Cs_rotation(this->exx_ptr->get_abfs_nchis());
+        this->symrot_.set_abfs_Lmax(GlobalC::exx_info.info_ri.abfs_Lmax);
         this->symrot_.cal_Ms(kv, ucell, pv);
     }
 }

source/source_lcao/module_ri/RPA_LRI.hpp

@@ -90,9 +90,15 @@ void RPA_LRI<T, Tdata>::cal_postSCF_exx(const elecstate::DensityMatrix<T, Tdata>
     if (exx_spacegroup_symmetry)
     {
         const std::array<Tcell, Ndim> 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,
-            RI_Util::get_Born_von_Karmen_cells(period), period, ucell.lat);
+            Rs, period, ucell.lat);
+        // set Lmax of the rotation matrices to max(l_ao, l_abf), to support rotation under ABF
+        symrot.set_abfs_Lmax(GlobalC::exx_info.info_ri.abfs_Lmax);
         symrot.cal_Ms(kv, ucell, *dm.get_paraV_pointer());
+        // 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);
         mix_DMk_2D.mix(symrot.restore_dm(kv, dm.get_DMK_vector(), *dm.get_paraV_pointer()), true);
     }
     else { mix_DMk_2D.mix(dm.get_DMK_vector(), true); }

source/source_lcao/module_ri/module_exx_symmetry/CMakeLists.txt

@@ -4,6 +4,7 @@ if (ENABLE_LIBRI)
         irreducible_sector.cpp
         irreducible_sector_bvk.cpp
         symmetry_rotation.cpp
+        symmetry_rotation_output.cpp
     )
     add_library(
         module_exx_symmetry

source/source_lcao/module_ri/module_exx_symmetry/irreducible_sector.cpp

@@ -148,7 +148,7 @@ namespace ModuleSymmetry
         if(GlobalV::MY_RANK == 0)
         {
             std::ofstream ofs;
-            ofs.open(PARAM.globalv.global_out_dir + "irreducible_sector.dat");
+            ofs.open(PARAM.globalv.global_out_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";}

source/source_lcao/module_ri/module_exx_symmetry/irreducible_sector.h

@@ -68,6 +68,17 @@ namespace ModuleSymmetry
             const ModuleBase::Matrix3& gmatd, const TCdouble gtransd,
             const TCdouble& posd_a1, const TCdouble& posd_a2)const;
 
+        // Getting calculated return lattice
+        TCdouble get_return_lattice(const int iat, const int isym) const
+        {
+            if (iat < 0 || iat >= static_cast<int>(this->return_lattice_.size())) {
+                throw std::out_of_range("Invalid atom index in get_return_lattice");
+            }
+            if (isym < 0 || isym >= static_cast<int>(this->return_lattice_[iat].size())) {
+                throw std::out_of_range("Invalid symmetry index in get_return_lattice");
+            }
+            return this->return_lattice_[iat][isym];
+        }
     protected:
         void cal_return_lattice_all(const Symmetry& symm, const Atom* atoms, const Statistics& st);
 

source/source_lcao/module_ri/module_exx_symmetry/symmetry_rotation.cpp

@@ -16,9 +16,7 @@ namespace ModuleSymmetry
     {
         this->reduce_Cs_ = true;
         this->abfs_l_nchi_ = abfs_l_nchi;
-        this->abfs_Lmax_ = 0;
-        for (auto& abfs_T : abfs_l_nchi) { this->abfs_Lmax_ = std::max(this->abfs_Lmax_, static_cast<int>(abfs_T.size()) - 1);
-}
+        for (auto& abfs_T : abfs_l_nchi) { this->abfs_Lmax_ = std::max(this->abfs_Lmax_, static_cast<int>(abfs_T.size()) - 1); }
     }
     void Symmetry_rotation::cal_Ms(const K_Vectors& kv,
         //const std::vector<std::map<int, TCdouble>>& kstars,
@@ -36,9 +34,8 @@ namespace ModuleSymmetry
         }
         // 1. calculate the rotation matrix in real spherical harmonics representation for each symmetry operation: [T_l (isym)]_mm'
         std::vector<ModuleBase::Matrix3> gmatc(nsym_);
-        for (int i = 0;i < nsym_;++i) { gmatc[i] = this->irs_.direct_to_cartesian(ucell.symm.gmatrix[i], ucell.latvec);
-}
-        this->cal_rotmat_Slm(gmatc.data(), reduce_Cs_ ? std::max(this->abfs_Lmax_, ucell.lmax) : ucell.lmax);
+        for (int i = 0;i < nsym_;++i) { gmatc[i] = this->irs_.direct_to_cartesian(ucell.symm.gmatrix[i], ucell.latvec); }
+        this->cal_rotmat_Slm(gmatc.data(), std::max(this->abfs_Lmax_, ucell.lmax));
 
         // 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
@@ -47,12 +44,9 @@ namespace ModuleSymmetry
                 kvec_res.x = fmod(kvec.x + 100.5 - 0.5 * symm_prec, 1) - 0.5 + 0.5 * symm_prec;
                 kvec_res.y = fmod(kvec.y + 100.5 - 0.5 * symm_prec, 1) - 0.5 + 0.5 * symm_prec;
                 kvec_res.z = fmod(kvec.z + 100.5 - 0.5 * symm_prec, 1) - 0.5 + 0.5 * symm_prec;
-                if (std::abs(kvec_res.x) < symm_prec) { kvec_res.x = 0.0;
-}
-                if (std::abs(kvec_res.y) < symm_prec) { kvec_res.y = 0.0;
-}
-                if (std::abs(kvec_res.z) < symm_prec) { kvec_res.z = 0.0;
-}
+                if (std::abs(kvec_res.x) < symm_prec) { kvec_res.x = 0.0; }
+                if (std::abs(kvec_res.y) < symm_prec) { kvec_res.y = 0.0; }
+                if (std::abs(kvec_res.z) < symm_prec) { kvec_res.z = 0.0; }
                 return kvec_res;
             };
         int nks_ibz = kv.kstars.size(); // kv.nks = 2 * kv.nks_ibz when nspin=2

source/source_lcao/module_ri/module_exx_symmetry/symmetry_rotation.h

@@ -27,13 +27,21 @@ namespace ModuleSymmetry
         {
             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<std::vector<RI::Tensor<std::complex<double>>>>& rotmat_Slm = this->rotmat_Slm_;
+        const int& abfs_Lmax = this->abfs_Lmax_;
         //--------------------------------------------------------------------------------
         // setters
         void find_irreducible_sector(const Symmetry& symm, const Atom* atoms, const Statistics& st,
             const std::vector<TC>& Rs, const TC& period, const Lattice& lat)
         {
             this->irs_.find_irreducible_sector(symm, atoms, st, Rs, period, lat);
         }
+        void set_abfs_Lmax(const int l) { this->abfs_Lmax_ = l; }
         void set_Cs_rotation(const std::vector<std::vector<int>>& abfs_l_nchi);
         //--------------------------------------------------------------------------------
         /// functions  to contruct rotation matrix in AO-representation
@@ -167,6 +175,20 @@ namespace ModuleSymmetry
         Irreducible_Sector irs_;
 
     };
+
+    template<typename T>  std::string vec3_fmt(const T& x, const T& y, const T& z)
+    {
+        return  "(" + std::to_string(x) + " " + std::to_string(y) + " " + std::to_string(z) + ")";
+    }
+    template<typename T>  std::string vec3_fmt(const ModuleBase::Vector3<T>& v)
+    {
+        return vec3_fmt(v.x, v.y, v.z);
+    }
+    // output k stars and the rotation matrices of Bloch orbitals
+    void print_symrot_info_k(const ModuleSymmetry::Symmetry_rotation& symrot,
+        const K_Vectors& kv, const UnitCell& ucell);
+    void print_symrot_info_R(const Symmetry_rotation& symrot, const Symmetry& symm,
+        const int lmax_ao, const std::vector<TC>& Rs);
 }
 
 #include "symmetry_rotation_R.hpp"

source/source_lcao/module_ri/module_exx_symmetry/symmetry_rotation_output.cpp

@@ -0,0 +1,100 @@
+#include "./symmetry_rotation.h"
+namespace ModuleSymmetry
+{
+    std::string mat3_fmt(const ModuleBase::Matrix3& m)
+    {
+        auto s = [](auto x) { return std::to_string(x); };
+        return s(m.e11) + " " + s(m.e12) + " " + s(m.e13) + "\n" +
+               s(m.e21) + " " + s(m.e22) + " " + s(m.e23) + "\n" +
+               s(m.e31) + " " + s(m.e32) + " " + s(m.e33);
+    }
+
+    // needs to calculate Ts from l=0 to l=max(l_ao,l_abf) before
+
+    void print_symrot_info_R(const Symmetry_rotation& symrot, const Symmetry& symm,
+        const int lmax_ao, const std::vector<TC>& Rs)
+    {
+        ModuleBase::TITLE("ModuleSymmetry", "print_symrot_info_R");
+        std::ofstream ofs(PARAM.globalv.global_out_dir + "symrot_R.txt");
+        // Print the irreducible sector (to be optimized)
+        ofs << "Number of irreducible sector: " << symrot.get_irreducible_sector().size() << std::endl;
+        ofs << "Lmax of AOs: " << lmax_ao << "\n";
+        ofs << "Lmax of ABFs: " << symrot.abfs_Lmax << "\n";
+        // print AO rotation matrix T
+        ofs << "Format:\n"
+            << "The index of the symmetry operation\n"
+            << "The rotation matrix of this symmetry operation (3*3)\n"
+            << "(The translation vector of this symmetry operation)\n"
+            << "Orbital rotation matrix (T) of each angular momentum with size ((2l + 1) * (2l + 1)) \n\n";
+        const int lmax = std::max(lmax_ao, symrot.abfs_Lmax);
+        for (int isym = 0;isym < symm.nrotk;++isym)
+        {
+            ofs << isym << "\n" << mat3_fmt(symm.gmatrix[isym]) << "\n"
+                << vec3_fmt(symm.gtrans[isym]) << "\n";
+            for (int l=0;l <= lmax;++l)
+            {
+                const int nm = 2 * l + 1;
+                // ofs << "l = " << l << ", nm = " << nm << "\n";
+                const auto& T_block = symrot.rotmat_Slm[isym][l];
+                for (int m1 = 0;m1 < nm;++m1)
+                {
+                    for (int m2 = 0;m2 < nm;++m2)
+                    {
+                        //note: the order of m in orbitals may be different from increasing
+                        //note: is Ts row- or col-major ?
+                        ofs << T_block(m1, m2);
+                    }
+                    ofs << "\n";
+                }
+            }
+            }
+        ofs.close();
+    }
+
+    void print_symrot_info_k(const Symmetry_rotation& symrot, const K_Vectors& kv, const UnitCell& ucell)
+    {
+        ModuleBase::TITLE("Symmetry_rotation", "print_symrot_info_k");
+        std::ofstream ofs(PARAM.globalv.global_out_dir + "symrot_k.txt");
+        ofs << "Number of IBZ k-points (k stars): " << kv.kstars.size() << std::endl;
+        ofs << "Format:\n" << "The symmetry operation index to the irreducible k-point. For the irreducible k-points, isym=0.\n\n"
+            << "(The direct coordinate of the original k-point)\n"
+            << "For each atom: \n"
+            << "- Original index->transformed index, type and the Lmax\n"
+            << "- Bloch orbital rotation matrix (M) of the given operation and atom, for each angular momentum\n\n";
+        for (int istar = 0;istar < kv.kstars.size();++istar)
+        {
+            ofs << "Star " << istar + 1 << " of IBZ k-point " << vec3_fmt(kv.kstars[istar].at(0)) << ":\n";
+            for (const auto& isym_kvd : kv.kstars[istar])
+            {
+                const int& isym = isym_kvd.first;
+                ofs << isym << "\n" << vec3_fmt(isym_kvd.second) << "\n";
+                for (int iat1 =0;iat1 < ucell.nat;++iat1)
+                {
+                    const int it = ucell.iat2it[iat1];  // it1=it2
+                    const int lmax = ucell.atoms[it].nwl;
+                    const int iat2 = ucell.symm.get_rotated_atom(isym, iat1);
+                    const double arg = 2 * ModuleBase::PI * isym_kvd.second * symrot.get_return_lattice(iat1,isym);
+                    std::complex<double>phase_factor = std::complex<double>(std::cos(arg), std::sin(arg));
+                    ofs << "atom " << iat1 + 1 << " -> " << iat2 + 1 << " of type " << it + 1 << " with Lmax= " << lmax << "\n";
+                    for (int l = 0;l < lmax + 1;++l)
+                    {
+                        const int nm = 2 * l + 1;
+                        const auto& m_block = symrot.rotmat_Slm[isym][l];
+                        for (int m1 = 0;m1 < nm;++m1)
+                        {
+                            // const int m1_start = m2 * nm;
+                            for (int m2 = 0;m2 < nm;++m2)
+                            {
+                                ofs << phase_factor * m_block(m1, m2);    // row-major
+                            }
+                            ofs << "\n";
+                        }
+                    }// end l
+                }   // end iat
+            }   // end (k, op)
+            ofs << "\n";
+        }   // end star
+        ofs.close();
+        ModuleBase::timer::tick("Symmetry_rotation", "print_symrot_info_k");
+    }
+}