Fix: symmetrize charge density in reciprocal space in any case

source/module_cell/module_symmetry/symmetry.cpp

@@ -22,7 +22,7 @@ Symmetry::~Symmetry()
 
 int Symmetry::symm_flag = 0;
 bool Symmetry::symm_autoclose = false;
-
+bool Symmetry::pricell_loop = true;
 
 void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running)
 {
@@ -244,6 +244,8 @@ void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running)
 
     this->set_atom_map(ucell);
 
+    if (GlobalV::NSPIN > 1) pricell_loop = this->magmom_same_check(ucell);
+
 	delete[] newpos;
     delete[] na;
     delete[] rotpos;
@@ -1428,8 +1430,9 @@ void Symmetry::rhog_symmetry(std::complex<double> *rhogtot,
 	// allocate flag for each FFT grid.
     int* symflag = new int[fftnx*fftny*fftnz];// which group the grid belongs to
     int(*isymflag)[48] = new int[fftnx*fftny*fftnz][48];//which rotration operation the grid corresponds to
-    int(*table_xyz)[48] = new int[fftnx*fftny*fftnz][48];// group information
-    for (int i=0; i<fftnx*fftny*fftnz; i++)
+    int(*table_xyz)[48] = new int[fftnx * fftny * fftnz][48];// group information
+    int* count_xyz = new int[fftnx * fftny * fftnz];// how many symmetry operations has been covered
+    for (int i = 0; i < fftnx * fftny * fftnz; i++)
     {
         symflag[i] = -1;
     }
@@ -1529,8 +1532,9 @@ ModuleBase::timer::tick("Symmetry","group fft grids");
                         symflag[ixyz] = group_index;
                         isymflag[group_index][rot_count] = invmap[isym];
                         table_xyz[group_index][rot_count] = ixyz;
-                         ++rot_count;
-                        assert(rot_count<=nrotk);
+                        ++rot_count;
+                        assert(rot_count <= nrotk);
+                        count_xyz[group_index] = rot_count;
                     }
                 group_index++;
                 }
@@ -1560,10 +1564,10 @@ for (int g_index = 0; g_index < group_index; g_index++)
     std::complex<double> sum(0, 0);
     int rot_count=0;
 
-    for (int c_index = 0; c_index < nrotk; ++c_index)
+    for (int c_index = 0; c_index < count_xyz[g_index]; ++c_index)
     {
-                int ixyz0=table_xyz[g_index][c_index];
-                int ipw0=ixyz2ipw[ixyz0];
+        int ixyz0 = table_xyz[g_index][c_index];
+        int ipw0 = ixyz2ipw[ixyz0];
                 if (symflag[ixyz0] == g_index)
                 {
                     // note : do not use PBC after rotation. 
@@ -1582,7 +1586,7 @@ for (int g_index = 0; g_index < group_index; g_index++)
                     double arg_gtrans = tmp_gdirect_double * gtrans[isymflag[g_index][c_index]];
                     std::complex<double> phase_gtrans (ModuleBase::libm::cos(arg_gtrans), ModuleBase::libm::sin(arg_gtrans));
                     // for each pricell in supercell:
-                    for (int ipt = 0;ipt < this->ncell;++ipt)
+                    for (int ipt = 0;ipt < ((ModuleSymmetry::Symmetry::pricell_loop) ? this->ncell : 1);++ipt)
                     {
                         double arg = tmp_gdirect_double * ptrans[ipt];
                         double tmp_cos = 0.0, tmp_sin = 0.0;
@@ -2083,4 +2087,26 @@ void Symmetry::hermite_normal_form(const ModuleBase::Matrix3 &s3, ModuleBase::Ma
 #endif
     return;
 }
+
+bool Symmetry::magmom_same_check(const UnitCell& ucell) const
+{
+    ModuleBase::TITLE("Symmetry", "magmom_same_check");
+    bool pricell_loop = true;
+    for (int it = 0;it < ntype;++it)
+    {
+        if (pricell_loop)
+            for (int ia = 1;ia < ucell.atoms[it].na;++ia)
+            {
+                if (!equal(ucell.atoms[it].m_loc_[ia].x, ucell.atoms[it].m_loc_[0].x) ||
+                    !equal(ucell.atoms[it].m_loc_[ia].y, ucell.atoms[it].m_loc_[0].y) ||
+                    !equal(ucell.atoms[it].m_loc_[ia].z, ucell.atoms[it].m_loc_[0].z))
+                {
+                    pricell_loop = false;
+                    break;
+                }
+            }
+    }
+    return pricell_loop;
 }
+
+}

source/module_cell/module_symmetry/symmetry.h

@@ -19,6 +19,7 @@ class Symmetry : public Symmetry_Basic
 	//1 : point group symmetry is considered
     static int symm_flag;
     static bool symm_autoclose;
+    static bool pricell_loop;   ///< whether to loop primitive cell in rhog_symmetry
 
 	void analy_sys(const UnitCell &ucell, std::ofstream &ofs_running);
 	bool available;
@@ -140,7 +141,10 @@ class Symmetry : public Symmetry_Basic
 	void get_optlat(ModuleBase::Vector3<double> &v1, ModuleBase::Vector3<double> &v2, 
 			ModuleBase::Vector3<double> &v3, ModuleBase::Vector3<double> &w1, 
 			ModuleBase::Vector3<double> &w2, ModuleBase::Vector3<double> &w3, 
-			int& real_brav, double* cel_const, double* tmp_const)const;
+        int& real_brav, double* cel_const, double* tmp_const)const;
+
+    /// Loop the magmom of each atoms in its type when NSPIN>1. If not all the same, primitive cells should not be looped in rhog_symmetry.
+    bool magmom_same_check(const UnitCell& ucell)const;
 };
 }
 

source/module_elecstate/module_charge/symmetry_rho.cpp

@@ -21,12 +21,12 @@ void Symmetry_rho::begin(const int& spin_now,
 
 	if(ModuleSymmetry::Symmetry::symm_flag != 1) return;
 	// both parallel and serial
-	if(symm.nrot==symm.nrotk) //pure point-group, do rho_symm in real space
-	{
-		psymm(CHR.rho[spin_now], rho_basis, Pgrid, symm);
-		if(XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) psymm(CHR.kin_r[spin_now], rho_basis,Pgrid,symm);
-	}
-	else	//space group, do rho_symm in reciprocal space
+    // if(symm.nrot==symm.nrotk) //pure point-group, do rho_symm in real space
+    // {
+    // 	psymm(CHR.rho[spin_now], rho_basis, Pgrid, symm);
+    // 	if(XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) psymm(CHR.kin_r[spin_now], rho_basis,Pgrid,symm);
+    // }
+    // else	//space group, do rho_symm in reciprocal space
 	{
 		rho_basis->real2recip(CHR.rho[spin_now], CHR.rhog[spin_now]);
 		psymmg(CHR.rhog[spin_now], rho_basis, Pgrid, symm);	//need to modify