Feature: Space group symmetry for EXX based on LibRI v0.2.1.0

CMakeLists.txt

@@ -535,7 +535,7 @@ if(ENABLE_LIBRI)
   else()
     message(FATAL_ERROR "Must provide LIBRI_DIR for RI related features.")
   endif()
-  target_link_libraries(${ABACUS_BIN_NAME} ri)
+  target_link_libraries(${ABACUS_BIN_NAME} ri module_exx_symmetry)
   add_compile_definitions(__EXX EXX_DM=3 EXX_H_COMM=2 TEST_EXX_LCAO=0
                           TEST_EXX_RADIAL=1)
 endif()

docs/advanced/input_files/input-main.md

@@ -243,6 +243,7 @@
     - [exx\_opt\_orb\_ecut](#exx_opt_orb_ecut)
     - [exx\_opt\_orb\_tolerence](#exx_opt_orb_tolerence)
     - [exx\_real\_number](#exx_real_number)
+    - [exx\_symmetry\_realspace](#exx_symmetry_realspace)
     - [rpa\_ccp\_rmesh\_times](#rpa_ccp_rmesh_times)
   - [Molecular dynamics](#molecular-dynamics)
     - [md\_type](#md_type)
@@ -2425,6 +2426,15 @@ These variables are relevant when using hybrid functionals.
 - **Description**: How many times larger the radial mesh required is to that of atomic orbitals in the postprocess calculation of the **bare** Coulomb matrix for RPA, GW, etc.
 - **Default**: 10
 
+### exx_symmetry_realspace
+
+- **Type**: Boolean
+- **Availability**: *[symmetry](#symmetry)==1* and exx calculation  (*[dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0/opt_orb* or *[rpa](#rpa)==True*)
+- **Description**: 
+  - False: only rotate k-space density matrix D(k) from irreducible k-points to accelerate diagonalization
+  - True: rotate both D(k) and Hexx(R) to accelerate both diagonalization and EXX calculation
+- **Default**: True
+
 [back to top](#full-list-of-input-keywords)
 
 ## Molecular dynamics

examples/hse/lcao_Si2/INPUT

@@ -10,7 +10,7 @@ ecutwfc                 100
 scf_thr                 1e-8
 scf_nmax                100
 gamma_only              0
-symmetry                -1
+symmetry                1
 smearing_method         gauss
 mixing_type             broyden
 mixing_beta             0.4

source/Makefile.Objects

@@ -602,6 +602,8 @@ OBJS_MODULE_RI=conv_coulomb_pot_k.o\
       RI_2D_Comm.o\
       Mix_DMk_2D.o\
       Mix_Matrix.o\
+      symmetry_rotation.o\
+      symmetry_irreducible_sector.o\
 
 OBJS_PARALLEL=parallel_common.o\
     parallel_global.o\

source/module_base/scalapack_connector.h

@@ -208,7 +208,7 @@ class ScalapackConnector
 			B, &IB, &JB, DESCB, &beta, C, &IC, &JC, DESCC);
 	}
 
-	static inline
+    static inline
 	void getrf(
 		const int M, const int N, 
 		std::complex<double> *A, const int IA, const int JA, const int *DESCA,

source/module_cell/klist.cpp

@@ -876,9 +876,9 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry& symm,
         // if really there is no equivalent k point in the list, then add it.
         if (!already_exist)
         {
-            // if it's a new ibz kpoint.
-            // nkstot_ibz indicate the index of ibz kpoint.
-            kvec_d_ibz[nkstot_ibz] = kvec_rot;
+			//if it's a new ibz kpoint.
+			//nkstot_ibz indicate the index of ibz kpoint.
+            kvec_d_ibz[nkstot_ibz] = kvec_d[i];
             // output in kpoints file
             ibz_index[i] = nkstot_ibz;
 
@@ -922,6 +922,38 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry& symm,
     }
 
     delete[] kkmatrix;
+
+#ifdef __EXX
+    // setup kstars according to the final (max-norm) kvec_d_ibz
+    this->kstars.resize(nkstot_ibz);
+    if (ModuleSymmetry::Symmetry::symm_flag == 1)
+    {
+        for (int i = 0; i < nkstot; ++i)
+        {
+            int exist_number = -1;
+            int isym = 0;
+            for (int j = 0; j < nrotkm; ++j)
+            {
+                kvec_rot = kvec_d[i] * kgmatrix[j];
+                restrict_kpt(kvec_rot);
+                for (int k = 0; k < nkstot_ibz; ++k)
+                {
+                    if (symm.equal(kvec_rot.x, kvec_d_ibz[k].x) &&
+                        symm.equal(kvec_rot.y, kvec_d_ibz[k].y) &&
+                        symm.equal(kvec_rot.z, kvec_d_ibz[k].z))
+                    {
+                        isym = j;
+                        exist_number = k;
+                        break;
+                    }
+                }
+                if (exist_number != -1) break;
+            }
+            this->kstars[exist_number].insert(std::make_pair(isym, kvec_d[i]));
+        }
+    }
+#endif
+
     // output in kpoints file
     std::stringstream ss;
     ss << " " << std::setw(40) << "nkstot"
@@ -1194,6 +1226,40 @@ void K_Vectors::mpi_k()
         wk[i] = wk_aux[k_index];
         isk[i] = isk_aux[k_index];
     }
+
+#ifdef __EXX
+    if (ModuleSymmetry::Symmetry::symm_flag == 1)
+    {//bcast kstars
+        this->kstars.resize(nkstot);
+        for (int ikibz = 0;ikibz < nkstot;++ikibz)
+        {
+            int starsize = this->kstars[ikibz].size();
+            Parallel_Common::bcast_int(starsize);
+            GlobalV::ofs_running << "starsize: " << starsize << std::endl;
+            auto ks = this->kstars[ikibz].begin();
+            for (int ik = 0;ik < starsize;++ik)
+            {
+                int isym = 0;
+                ModuleBase::Vector3<double> ks_vec(0, 0, 0);
+                if (GlobalV::MY_RANK == 0)
+                {
+                    isym = ks->first;
+                    ks_vec = ks->second;
+                    ++ks;
+                }
+                Parallel_Common::bcast_int(isym);
+                Parallel_Common::bcast_double(ks_vec.x);
+                Parallel_Common::bcast_double(ks_vec.y);
+                Parallel_Common::bcast_double(ks_vec.z);
+                GlobalV::ofs_running << "isym: " << isym << " ks_vec: " << ks_vec.x << " " << ks_vec.y << " " << ks_vec.z << std::endl;
+                if (GlobalV::MY_RANK != 0)
+                {
+                    kstars[ikibz].insert(std::make_pair(isym, ks_vec));
+                }
+            }
+        }
+    }
+#endif
 } // END SUBROUTINE
 #endif
 

source/module_cell/klist.h

@@ -10,7 +10,7 @@
 
 class K_Vectors
 {
-  public:
+public:
     std::vector<ModuleBase::Vector3<double>> kvec_c; /// Cartesian coordinates of k points
     std::vector<ModuleBase::Vector3<double>> kvec_d; /// Direct coordinates of k points
 
@@ -22,6 +22,10 @@ class K_Vectors
     int nmp[3];                 /// Number of Monhorst-Pack
     std::vector<int> kl_segids; /// index of kline segment
 
+    /// @brief equal k points to each ibz-kpont, corresponding to a certain symmetry operations. 
+    /// dim: [iks_ibz][(isym, kvec_d)]
+    std::vector<std::map<int, ModuleBase::Vector3<double>>> kstars;
+
     K_Vectors();
     ~K_Vectors();
     K_Vectors& operator=(const K_Vectors&) = default;
@@ -47,11 +51,11 @@ class K_Vectors
      * @note Only available for nspin = 1 or 2 or 4.
      */
     void set(const ModuleSymmetry::Symmetry& symm,
-             const std::string& k_file_name,
-             const int& nspin,
-             const ModuleBase::Matrix3& reciprocal_vec,
-             const ModuleBase::Matrix3& latvec,
-             std::ofstream& ofs);
+        const std::string& k_file_name,
+        const int& nspin,
+        const ModuleBase::Matrix3& reciprocal_vec,
+        const ModuleBase::Matrix3& latvec,
+        std::ofstream& ofs);
 
     /**
      * @brief Generates irreducible k-points in the Brillouin zone considering symmetry operations.
@@ -67,10 +71,10 @@ class K_Vectors
      * @param match A boolean flag that indicates if the results matches the real condition.
      */
     void ibz_kpoint(const ModuleSymmetry::Symmetry& symm,
-                    bool use_symm,
-                    std::string& skpt,
-                    const UnitCell& ucell,
-                    bool& match);
+        bool use_symm,
+        std::string& skpt,
+        const UnitCell& ucell,
+        bool& match);
     // LiuXh add 20180515
 
     /**
@@ -149,7 +153,7 @@ class K_Vectors
         this->nkstot_full = value;
     }
 
-  private:
+private:
     int nks;         // number of symmetry-reduced k points in this pool(processor, up+dw)
     int nkstot;      /// number of symmetry-reduced k points in full k mesh
     int nkstot_full; /// number of k points before symmetry reduction in full k mesh
@@ -278,8 +282,8 @@ class K_Vectors
      * be recalculated.
      */
     void update_use_ibz(const int& nkstot_ibz,
-                        const std::vector<ModuleBase::Vector3<double>>& kvec_d_ibz,
-                        const std::vector<double>& wk_ibz);
+        const std::vector<ModuleBase::Vector3<double>>& kvec_d_ibz,
+        const std::vector<double>& wk_ibz);
 
     /**
      * @brief Sets both the direct and Cartesian k-vectors.
@@ -400,4 +404,4 @@ inline int K_Vectors::getik_global(const int& ik) const
     }
 }
 
-#endif // KVECT_H
+#endif // KVECT_H

source/module_cell/module_symmetry/symmetry.cpp

@@ -155,7 +155,8 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms,
                     this->itmin_start = istart[it];
                 }
             }
-            this->getgroup(nrot_out, nrotk_out, ofs_running, pos_spinup.data());
+            this->getgroup(nrot_out, nrotk_out, ofs_running, this->nop, this->symop, this->gmatrix, this->gtrans,
+                pos_spinup.data(), this->rotpos, this->index, this->itmin_type, this->itmin_start, this->istart, this->na);
             // recover na and istart
             for (int it = 0;it < ntype;++it)
             {
@@ -166,11 +167,14 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms,
             }
             // For AFM analysis End
             //------------------------------------------------------------
-        } else {
+        }
+        else
+        {
             // get the real symmetry operations according to the input structure
             // nrot_out: the number of pure point group rotations
             // nrotk_out: the number of all space group operations
-            this->getgroup(nrot_out, nrotk_out, ofs_running, this->newpos);
+            this->getgroup(nrot_out, nrotk_out, ofs_running, this->nop, this->symop, this->gmatrix, this->gtrans,
+                this->newpos, this->rotpos, this->index, this->itmin_type, this->itmin_start, this->istart, this->na);
         }
         };
 
@@ -289,16 +293,16 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms,
     // 3. output to running.log
     //----------------------------------
     // output the point group
-    this->pointgroup(this->nrot, this->pgnumber, this->pgname, this->gmatrix, ofs_running);
+    bool valid_group = this->pointgroup(this->nrot, this->pgnumber, this->pgname, this->gmatrix, ofs_running);
 	ModuleBase::GlobalFunc::OUT(ofs_running,"POINT GROUP", this->pgname);
     // output the space group
-    this->pointgroup(this->nrotk, this->spgnumber, this->spgname, this->gmatrix, ofs_running);
+    valid_group = this->pointgroup(this->nrotk, this->spgnumber, this->spgname, this->gmatrix, ofs_running);
     ModuleBase::GlobalFunc::OUT(ofs_running, "POINT GROUP IN SPACE GROUP", this->spgname);
 
     //-----------------------------
     // 4. For the case where point group is not complete due to symmetry_prec
     //-----------------------------
-    if (!this->valid_group)
+    if (!valid_group)
     {   // select the operations that have the inverse
         std::vector<int>invmap(this->nrotk, -1);
         this->gmatrix_invmap(this->gmatrix, this->nrotk, invmap.data());
@@ -883,7 +887,9 @@ void Symmetry::lattice_type(
 }
 
 
-void Symmetry::getgroup(int& nrot, int& nrotk, std::ofstream& ofs_running, double* pos)
+void Symmetry::getgroup(int& nrot, int& nrotk, std::ofstream& ofs_running, const int& nop,
+    const ModuleBase::Matrix3* symop, ModuleBase::Matrix3* gmatrix, ModuleBase::Vector3<double>* gtrans,
+    double* pos, double* rotpos, int* index, const int itmin_type, const int itmin_start, int* istart, int* na)const
 {
     ModuleBase::TITLE("Symmetry", "getgroup");
 
@@ -907,9 +913,7 @@ void Symmetry::getgroup(int& nrot, int& nrotk, std::ofstream& ofs_running, doubl
     //std::cout << "nop = " <<nop <<std::endl;
     for (int i = 0; i < nop; ++i)
     {
-    //    std::cout << "symop = " << symop[i].e11 <<" "<< symop[i].e12 <<" "<< symop[i].e13 <<" "<< symop[i].e21 <<" "<< symop[i].e22 <<" "<< symop[i].e23 <<" "<< symop[i].e31 <<" "<< symop[i].e32 <<" "<< symop[i].e33 << std::endl;
-        this->checksym(this->symop[i], this->gtrans[i], pos);
-      //  std::cout << "s_flag =" <<s_flag<<std::endl;
+        bool s_flag = this->checksym(symop[i], gtrans[i], pos, rotpos, index, itmin_type, itmin_start, istart, na);
         if (s_flag == 1)
         {
 			//------------------------------
@@ -987,7 +991,8 @@ void Symmetry::getgroup(int& nrot, int& nrotk, std::ofstream& ofs_running, doubl
     return;
 }
 
-void Symmetry::checksym(ModuleBase::Matrix3 &s, ModuleBase::Vector3<double> &gtrans, double* pos)
+bool Symmetry::checksym(const ModuleBase::Matrix3& s, ModuleBase::Vector3<double>& gtrans,
+    double* pos, double* rotpos, int* index, const int itmin_type, const int itmin_start, int* istart, int* na)const
 {
 	//----------------------------------------------
     // checks whether a point group symmetry element 
@@ -996,7 +1001,7 @@ void Symmetry::checksym(ModuleBase::Matrix3 &s, ModuleBase::Vector3<double> &gtr
     // the start atom index.
     bool no_diff = false;
     ModuleBase::Vector3<double> trans(2.0, 2.0, 2.0);
-    s_flag = 0;
+    bool s_flag = 0;
 
     for (int it = 0; it < ntype; it++)
     {
@@ -1051,9 +1056,9 @@ void Symmetry::checksym(ModuleBase::Matrix3 &s, ModuleBase::Vector3<double> &gtr
 	//---------------------------------------------------------
     // itmin_start = the start atom positions of species itmin
 	//---------------------------------------------------------
-    sptmin.x = rotpos[itmin_start*3];
-    sptmin.y = rotpos[itmin_start*3+1];
-    sptmin.z = rotpos[itmin_start*3+2];
+    // (s)tart (p)osition of atom (t)ype which has (min)inal number.
+    ModuleBase::Vector3<double> sptmin(rotpos[itmin_start * 3], rotpos[itmin_start * 3 + 1], rotpos[itmin_start * 3 + 2]);
+
     for (int i = itmin_start; i < itmin_start + na[itmin_type]; ++i)
     {
         //set up the current test std::vector "gtrans"
@@ -1143,13 +1148,12 @@ void Symmetry::checksym(ModuleBase::Matrix3 &s, ModuleBase::Vector3<double> &gtr
         gtrans.y = trans.y;
         gtrans.z = trans.z;
     }
-    return;
+    return s_flag;
 }
 
 void Symmetry::pricell(double* pos, const Atom* atoms)
 {
     bool no_diff = false;
-    s_flag = 0;
     ptrans.clear();
 
     for (int it = 0; it < ntype; it++)
@@ -1184,10 +1188,10 @@ void Symmetry::pricell(double* pos, const Atom* atoms)
 
 	//---------------------------------------------------------
     // itmin_start = the start atom positions of species itmin
-	//---------------------------------------------------------
-    sptmin.x = pos[itmin_start*3];
-    sptmin.y = pos[itmin_start*3+1];
-    sptmin.z = pos[itmin_start*3+2];
+    //---------------------------------------------------------
+    // (s)tart (p)osition of atom (t)ype which has (min)inal number.
+    ModuleBase::Vector3<double> sptmin(pos[itmin_start * 3], pos[itmin_start * 3 + 1], pos[itmin_start * 3 + 2]);
+
     for (int i = itmin_start; i < itmin_start + na[itmin_type]; ++i)
     {
         //set up the current test std::vector "gtrans"
@@ -1906,7 +1910,7 @@ void Symmetry::gtrans_convert(const ModuleBase::Vector3<double>* va, ModuleBase:
           vb[i]=va[i]*a*bi;
     }
 }
-void Symmetry::gmatrix_invmap(const ModuleBase::Matrix3* s, const int n, int* invmap)
+void Symmetry::gmatrix_invmap(const ModuleBase::Matrix3* s, const int n, int* invmap) const
 {
     ModuleBase::Matrix3 eig(1, 0, 0, 0, 1, 0, 0, 0, 1);
     ModuleBase::Matrix3 tmp;