Symmetry: debug and new features (#1904)

* add a new atom_ordering algorithm in module symmetry

* fix wrong transpose

* add mat_rotate to ibrav=4,12,13 according to euler angles

* fix a gtrans=0 bug

* add space group analysis

* change euler angles according to their defination

* delete the duplicated call of point-group functions

* rotate the matrix back in k-space

* revert the transpose in force_symmetry

* remove unused gmatrix_fft

* revert to transpose with a special transpose for 60-hexagonal in mat_rotate

* fix a typo

* revert the double-call with more detailed output

* reverte back-rotation

* output bravais and point group analysis for input configuration

* simplify some variables

* add gmat_con and gtrans_con to bring optimized-symm into input config

* fix unset-newpos bug

* fix inverted-sign bug of gtrans

* implement rhog_symm and invmap

* add some unittests for trans-funcs in module_symmetry

* **implement rhog_symmetry both parallel and serial**

* fix some bugs

* modify the parallel of rhog_symmetry

* add round option in gmatrix_convert

* divide the phase factor back to fix space-group bug

* modify UT; divide gmatrix_convet_int from gmatrix_convert

* checksym: add gtrans firstly and then check_boundary

* add an imcomplete pricell

* complete pricell

* fix pricell bugs

* modify algorithm in pricell but still have bug

* fix a gtrans-missing bug in checksym

* modify pcell-lattice-vector-searching algorithm

* fix a bug: change 0-init into -1-init

* use pricell onto rhog_symmetry

* add get_shortest_latvec

* extract get_optlat as a subfunction from lattice_type

* add plat_type

* call the get_shortest_latvec and plat_type

* modify plat-algorithm and deal with collineation

* pricell: exclude collinear and coplane when searching b1b2b3

* pcell in rhosym: fix a nan-bug resulting from zero mean-gphase

* fix a compile error

* pricell: fix the condition searching jplane and iplane

* pricel in rhosym: add the missing ncell-division back

* comment out some couts

* pricell: fix collinear-condition

* rhosym: modify rotation

* add a space to output for convenient grep

* setgroup: change r2yp into r2zp in ibrav=12, fitting to mat_rotate

* setgroup: modify cel_const and use mat_rotate in ibrav=11

* improve cel_const for ibrav=12, 13

* standard_lat: strengthen conditions to ibrav-setting and remove mat_rotate

* remove redundant output and fix zero-celconst for triclinc

* convert ptrans to input configuration

* use a new rotate_recip without scaling gmatrix by the number of FFT-grids

* retain s1, s2, s3 as input lat-vec

* delete a exclamation mark

* restrict kvec_d into [0, 1)

* rhogsym: deal with error in gphase

* fix unittest

* fix a type error

* remove cel_const dependence from setgroup

* ibzkpt: use kgmatrix and row-vector-rotate (vec*mat)

* move the output out of lattice_type function

* symmetry: change some function into const and remove some useless code

* fix the dependence of the test

* fix fco-wrong-kgmatrix by changing the condition to modify optlat

* ibzkpt: add group-analysis in reciprocal and k-lattice, check and restrict kpoints

* pricell: add hermite_normal_form and remove plat_type

* fix the include in symmetry_rhog

* modify result.ref of some scf cases

* fix result.ref

* fix ibz_kpoint in unit tests

* set OMP_NUM_THREADS=1 in test workflow

* set scf_nmax=8 in case 121_PW_KPAR

* fix a MPI_Send bug sending to other pools (local to local, rather than global)

* rewrite  and fix KPAR-bug

* do rhog_symmetry respectively in each pool

* Revert "set OMP_NUM_THREADS=1 in test workflow"

This reverts commit c4719b6.

* Revert "set scf_nmax=8 in case 121_PW_KPAR"

This reverts commit a1305d2.

* fix Makefile

* clear ptrans befor every ion step

* modify result.ref

* add unittest to atom_ordering_new

* use EXPECT_LE

* delete some useless comments

* add introduction to unit test

Author: maki49

source/Makefile.Objects

@@ -479,6 +479,7 @@ OBJS_SRCPW=H_Ewald_pw.o\
     stress_pw.o\
     of_stress_pw.o\
     symmetry_rho.o\
+    symmetry_rhog.o\
     wavefunc.o\
     wf_atomic.o\
     wf_igk.o\

source/module_cell/klist.cpp

@@ -89,7 +89,7 @@ void K_Vectors::set(
     //if symm_flag is not set, only time-reversal symmetry would be considered.
     if(!berryphase::berry_phase_flag && ModuleSymmetry::Symmetry::symm_flag != -1)
     {
-        this->ibz_kpoint(symm, ModuleSymmetry::Symmetry::symm_flag,skpt1);
+        this->ibz_kpoint(symm, ModuleSymmetry::Symmetry::symm_flag, skpt1, GlobalC::ucell);
         if(ModuleSymmetry::Symmetry::symm_flag || is_mp)
         {
             this->update_use_ibz();
@@ -544,10 +544,21 @@ void K_Vectors::update_use_ibz( void )
     return;
 }
 
-void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,std::string& skpt)
+void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,std::string& skpt, const UnitCell &ucell)
 {
     if (GlobalV::MY_RANK!=0) return;
     ModuleBase::TITLE("K_Vectors", "ibz_kpoint");
+
+    // k-lattice: "pricell" of reciprocal space
+    // CAUTION: should fit into all k-input method, not only MP  !!!
+    ModuleBase::Vector3<double> gb1(ucell.G.e11, ucell.G.e12, ucell.G.e13);
+    ModuleBase::Vector3<double> gb2(ucell.G.e21, ucell.G.e22, ucell.G.e23);
+    ModuleBase::Vector3<double> gb3(ucell.G.e31, ucell.G.e32, ucell.G.e33);
+    ModuleBase::Vector3<double> gk1(gb1.x / nmp[0], gb1.y / nmp[0], gb1.z / nmp[0]);
+    ModuleBase::Vector3<double> gk2(gb2.x / nmp[1], gb2.y / nmp[1], gb2.z / nmp[1]);
+    ModuleBase::Vector3<double> gk3(gb3.x / nmp[2], gb3.y / nmp[2], gb3.z / nmp[2]);
+    ModuleBase::Matrix3 gk(gk1.x, gk1.y, gk1.z, gk2.x, gk2.y, gk2.z, gk3.x, gk3.y, gk3.z);
+
     //===============================================
     // search in all space group operations
     // if the operations does not already included
@@ -561,21 +572,71 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s
     int nrotkm;
     if(use_symm)
     {
+        // bravais type of reciprocal lattice and k-lattice
+        double b_const[6];
+        double b0_const[6];
+        double bk_const[6];
+        double bk0_const[6];
+        int bbrav=15;
+        int bkbrav=15;
+        std::string bbrav_name;
+        std::string bkbrav_name;
+        ModuleBase::Vector3<double> gb01=gb1, gb02=gb2, gb03=gb3, gk01=gk1, gk02=gk2, gk03=gk3;
+        symm.lattice_type(gb1, gb2, gb3, gb01, gb02, gb03, b_const, b0_const, bbrav, bbrav_name, ucell, false, nullptr);
+        GlobalV::ofs_running<<"(for reciprocal lattice: )"<<std::endl;
+        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"BRAVAIS TYPE", bbrav);
+        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"BRAVAIS LATTICE NAME", bbrav_name);
+        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"ibrav", bbrav);
+
+        symm.lattice_type(gk1, gk2, gk3, gk01, gk02, gk03, bk_const, bk0_const, bkbrav, bkbrav_name, ucell, false, nullptr);
+        GlobalV::ofs_running<<"(for k-lattice: )"<<std::endl;
+        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"BRAVAIS TYPE", bkbrav);
+        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"BRAVAIS LATTICE NAME", bkbrav_name);
+        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"ibrav", bkbrav);
+
+        // point-group analysis of reciprocal lattice
+        ModuleBase::Matrix3 bsymop[48];
+        int bnop=0;
+        symm.setgroup(bsymop, bnop, bbrav);
+        ModuleBase::Matrix3 b_optlat(gb1.x, gb1.y, gb1.z, gb2.x, gb2.y, gb2.z, gb3.x, gb3.y, gb3.z);
+        //symm.gmatrix_convert_int(bsymop, bsymop, bnop, b_optlat, ucell.G);
+        symm.gmatrix_convert(bsymop, bsymop, bnop, b_optlat, ucell.G);
+        //check if all the kgmatrix are in bsymop
+        auto matequal = [&symm] (ModuleBase::Matrix3 a, ModuleBase::Matrix3 b)
+        {
+            return (symm.equal(a.e11, b.e11) && symm.equal(a.e12, b.e12) && symm.equal(a.e13, b.e13) &&
+            symm.equal(a.e21, b.e21) && symm.equal(a.e22, b.e22) && symm.equal(a.e23, b.e23) &&
+            symm.equal(a.e31, b.e31) && symm.equal(a.e23, b.e23) && symm.equal(a.e33, b.e33));
+        };
+        for(int i=0;i<symm.nrotk;++i)
+        {
+            bool match = false;
+            for(int j=0;j<bnop;++j) 
+            {
+                if (matequal(symm.kgmatrix[i], bsymop[j])) {match=true; break;}
+            }
+            if(!match) 
+            {
+                std::cout<<"match failed:"<<std::endl;
+                ModuleBase::WARNING_QUIT("K_Vectors:ibz_kpoint", 
+                "symmetry operation of reciprocal lattice is wrong! ");
+            }
+        }
         nrotkm = symm.nrotk;// change if inv not included
         for (int i = 0; i < nrotkm; ++i)
         {
-            if (symm.gmatrix[i] == inv)
+            if (symm.kgmatrix[i] == inv)
             {
                 include_inv = true;
             }
-            kgmatrix[i] = symm.gmatrix[i];
+            kgmatrix[i] = symm.kgmatrix[i];
         }
 
         if (!include_inv)
         {
             for (int i = 0; i<symm.nrotk; ++i)
             {
-                kgmatrix[i + symm.nrotk] = inv * symm.gmatrix[i];
+                kgmatrix[i + symm.nrotk] = inv * symm.kgmatrix[i];
             }
             nrotkm = 2 * symm.nrotk;
         }
@@ -591,6 +652,13 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s
         return;
     }
 
+    // convert kgmatrix to k-lattice
+    ModuleBase::Matrix3* kkmatrix = new ModuleBase::Matrix3 [nrotkm];
+    symm.gmatrix_convert(kgmatrix.data(), kkmatrix, nrotkm, ucell.G, gk);
+    // direct coordinates of k-points in k-lattice
+    std::vector<ModuleBase::Vector3<double>> kvec_d_k(nkstot);
+    for (int i=0;i<nkstot;++i) kvec_d_k[i]=kvec_d[i]*ucell.G*gk.Inverse();
+
     // use operation : kgmatrix to find
     // the new set kvec_d : ir_kpt
     this->nkstot_ibz = 0;
@@ -604,18 +672,35 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s
     const double weight = 1.0 / static_cast<double>(nkstot);
 
     ModuleBase::Vector3<double> kvec_rot;
+    ModuleBase::Vector3<double> kvec_rot_k;
 
 
 //	for(int i=0; i<nrotkm; i++)
 //	{
 //		out.printM3("rot matrix",kgmatrix[i]);
 //	}
-
+    auto restrict_kpt = [&symm](ModuleBase::Vector3<double> &kvec){
+        // in (-0.5, 0.5]
+        kvec.x = fmod(kvec.x + 100.5-0.5*symm.epsilon, 1)-0.5+0.5*symm.epsilon;
+        kvec.y = fmod(kvec.y + 100.5-0.5*symm.epsilon, 1)-0.5+0.5*symm.epsilon;
+        kvec.z = fmod(kvec.z + 100.5-0.5*symm.epsilon, 1)-0.5+0.5*symm.epsilon;
+        // in [0, 1)
+        // kvec.x = fmod(kvec.x + 100 + symm.epsilon, 1) - symm.epsilon;
+        // kvec.y = fmod(kvec.y + 100 + symm.epsilon, 1) - symm.epsilon;
+        // kvec.z = fmod(kvec.z + 100 + symm.epsilon, 1) - symm.epsilon;
+        if(abs(kvec.x) < symm.epsilon) kvec.x = 0.0;
+        if(abs(kvec.y) < symm.epsilon) kvec.y = 0.0;
+        if(abs(kvec.z) < symm.epsilon) kvec.z = 0.0;
+        return;
+    };
     // for output in kpoints file
     int ibz_index[nkstot];
 	// search in all k-poins.
     for (int i = 0; i < nkstot; ++i)
     {
+        //restrict to [0, 1)
+        restrict_kpt(kvec_d[i]);
+
 		//std::cout << "\n kpoint = " << i << std::endl;
 		//std::cout << "\n kvec_d = " << kvec_d[i].x << " " << kvec_d[i].y << " " << kvec_d[i].z;
         bool already_exist = false;
@@ -631,13 +716,20 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s
 				// mohan modify 2010-01-30.
 				// mohan modify again 2010-01-31
 				// fix the bug like kvec_d * G; is wrong
-				//kvec_rot = kvec_d[i] * kgmatrix[j]; //wrong for total energy, but correct for nonlocal force.
-				kvec_rot = kgmatrix[j] * kvec_d[i]; //correct for total energy, but wrong for nonlocal force.
+				kvec_rot = kvec_d[i] * kgmatrix[j]; //wrong for total energy, but correct for nonlocal force.
+				//kvec_rot = kgmatrix[j] * kvec_d[i]; //correct for total energy, but wrong for nonlocal force.
+                restrict_kpt(kvec_rot);
 
-                kvec_rot.x = fmod(kvec_rot.x + 100, 1);
-                kvec_rot.y = fmod(kvec_rot.y + 100, 1);
-                kvec_rot.z = fmod(kvec_rot.z + 100, 1);
+                kvec_rot_k = kvec_d_k[i] * kkmatrix[j];   //k-lattice rotation
+                kvec_rot_k = kvec_rot_k * gk * ucell.G.Inverse(); //convert to recip lattice
+                restrict_kpt(kvec_rot_k);
 
+                assert(symm.equal(kvec_rot.x, kvec_rot_k.x));
+                assert(symm.equal(kvec_rot.y, kvec_rot_k.y));
+                assert(symm.equal(kvec_rot.z, kvec_rot_k.z));
+                // std::cout << "\n kvec_rot (in recip) = " << kvec_rot.x << " " << kvec_rot.y << " " << kvec_rot.z;
+                // std::cout << "\n kvec_rot(k to recip)= " << kvec_rot_k.x << " " << kvec_rot_k.y << " " << kvec_rot_k.z;
+                kvec_rot_k = kvec_rot_k * ucell.G * gk.Inverse(); //convert back to k-latice
 
 //				std::cout << "\n kvec_rot = " << kvec_rot.x << " " << kvec_rot.y << " " << kvec_rot.z;
 
@@ -708,6 +800,7 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s
 //		BLOCK_HERE("check k point");
     }
 
+    delete[] kkmatrix;
     // output in kpoints file
     std::stringstream ss;
     ss << " " << std::setw(40) <<"nkstot" << " = " << nkstot

source/module_cell/klist.h

@@ -38,7 +38,7 @@ class K_Vectors
         const ModuleBase::Matrix3 &reciprocal_vec,
         const ModuleBase::Matrix3 &latvec);
 
-    void ibz_kpoint( const ModuleSymmetry::Symmetry &symm, bool use_symm,std::string& skpt);
+    void ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,std::string& skpt, const UnitCell &ucell);
     //LiuXh add 20180515
     void set_after_vc(
             const ModuleSymmetry::Symmetry &symm,

source/module_cell/module_symmetry/symm_other.cpp

@@ -88,7 +88,7 @@ void Symm_Other::print1(const int &ibrav, const double *cel_const, std::ofstream
     	ModuleBase::GlobalFunc::OUT(ofs_running,"LATTICE CONSTANT A",cel_const[0]);
 		ModuleBase::GlobalFunc::OUT(ofs_running,"B/A RATIO",cel_const[1]);
 		ModuleBase::GlobalFunc::OUT(ofs_running,"C/A RATIO",cel_const[2]);
-		ModuleBase::GlobalFunc::OUT(ofs_running,"COS(BETA)",cel_const[4]);	
+		ModuleBase::GlobalFunc::OUT(ofs_running,"COS(BETA)",cel_const[4]);
 	}
 	else if(ibrav==14)
 	{