Feature: add a parameter symmetry_autoclose to control the behavior when symmetry error occurs

docs/advanced/input_files/input-main.md

@@ -403,12 +403,19 @@ These variables are used to control general system parameters.
 
 - **Type**: Real
 - **Description**: The accuracy for symmetry judgment. Usually the default value is good enough, but if the lattice parameters or atom positions in STRU file is not accurate enough, this value should be enlarged. 
-
-  Note: if *[calculation](#calculation)=cell_relax*, this value can be dynamically enlarged corresponding to the accuracy loss of the lattice parameters and atom positions during the relaxation. There will be a warning message in that case.
-
+  > Note: if *[calculation](#calculation)==cell_relax*, this value can be dynamically changed corresponding to the variation of accuracy of the lattice parameters and atom positions during the relaxation. The new value will be printed in `OUT.${suffix}/running_cell-relax.log` in that case.
 - **Default**: 1.0e-5
 - **Unit**:  Bohr
 
+### symmetry_autoclose
+
+- **Type**: Boolean
+- **Availability**: *[symmetry](#symmetry)==1*
+- **Description**: Control how to deal with error in symmetry analysis due to inaccurate lattice parameters or atom positions in STRU file, especially useful when *[calculation](#calculation)==cell-relax*
+  - False: quit with an error message
+  - True: automatically set symmetry to 0 and continue running without symmetry analysis
+- **Default**: False
+
 ### kpar
 
 - **Type**: Integer

source/module_cell/klist.cpp

@@ -98,7 +98,16 @@ void K_Vectors::set(
         {
             std::cout<< "Optimized lattice type of reciprocal lattice cannot match the optimized real lattice. " <<std::endl;
             std::cout << "It is often because the inaccuracy of lattice parameters in STRU." << std::endl;
-            ModuleBase::WARNING_QUIT("K_Vectors::ibz_kpoint", "Refine the lattice parameters in STRU or use a different`symmetry_prec`. ");
+            if (ModuleSymmetry::Symmetry::symm_autoclose)
+            {
+                ModuleBase::WARNING("K_Vectors::ibz_kpoint", "Automatically set symmetry to 0 and continue ...");
+                std::cout << "Automatically set symmetry to 0 and continue ..." << std::endl;
+                ModuleSymmetry::Symmetry::symm_flag = 0;
+                match = true;
+                this->ibz_kpoint(symm, ModuleSymmetry::Symmetry::symm_flag, skpt1, GlobalC::ucell, match);
+            }
+            else
+                ModuleBase::WARNING_QUIT("K_Vectors::ibz_kpoint", "Refine the lattice parameters in STRU or use a different`symmetry_prec`. ");
         }
         if (ModuleSymmetry::Symmetry::symm_flag || is_mp)
         {

source/module_cell/module_symmetry/symmetry.cpp

@@ -20,7 +20,8 @@ Symmetry::~Symmetry()
 }
 
 
-int Symmetry::symm_flag=0;
+int Symmetry::symm_flag = 0;
+bool Symmetry::symm_autoclose = false;
 
 
 void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running)
@@ -199,16 +200,17 @@ void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running)
             if (valid_index > 0)//epsilon is set smaller
                 ofs_running << " Narrowing `symmetry_prec` from " << eps_current << " to " << epsilon << " ..." << std::endl;
         }
-        // final number of symmetry operations
-#ifdef __DEBUG
-        ofs_running << "symmetry_prec(epsilon) in current ion step: " << this->epsilon << std::endl;
-        ofs_running << "number of symmetry operations in current ion step: " << this->nrotk << std::endl;
-#endif
     }
     else
         lattice_to_group(this->nrot, this->nrotk, ofs_running);
 
-	this->pointgroup(this->nrot, this->pgnumber, this->pgname, this->gmatrix, ofs_running);
+    // final number of symmetry operations
+#ifdef __DEBUG
+    ofs_running << "symmetry_prec(epsilon) in current ion step: " << this->epsilon << std::endl;
+    ofs_running << "number of symmetry operations in current ion step: " << this->nrotk << std::endl;
+#endif
+
+    this->pointgroup(this->nrot, this->pgnumber, this->pgname, this->gmatrix, ofs_running);
 	ModuleBase::GlobalFunc::OUT(ofs_running,"POINT GROUP", this->pgname);
     this->pointgroup(this->nrotk, this->spgnumber, this->spgname, this->gmatrix, ofs_running);
 	ModuleBase::GlobalFunc::OUT(ofs_running,"POINT GROUP IN SPACE GROUP", this->spgname);

source/module_cell/module_symmetry/symmetry.h

@@ -17,7 +17,8 @@ class Symmetry : public Symmetry_Basic
 	//-1 : no symmetry at all, k points would be total nks in KPT
 	//0 : only basic time-reversal symmetry is considered, point k and -k would fold to k
 	//1 : point group symmetry is considered
-	static int symm_flag;
+    static int symm_flag;
+    static bool symm_autoclose;
 
 	void analy_sys(const UnitCell &ucell, std::ofstream &ofs_running);
 	bool available;

source/module_io/input.cpp

@@ -200,6 +200,7 @@ void Input::Default(void)
     init_vel = false;
     ref_cell_factor = 1.0;
     symmetry_prec = 1.0e-5; // LiuXh add 2021-08-12, accuracy for symmetry
+    symmetry_autoclose = false; // whether to close symmetry automatically when error occurs in symmetry analysis
     cal_force = 0;
     force_thr = 1.0e-3;
     force_thr_ev2 = 0;
@@ -872,6 +873,10 @@ bool Input::Read(const std::string &fn)
         {
             read_value(ifs, symmetry_prec);
         }
+        else if (strcmp("symmetry_autoclose", word) == 0)
+        {
+            read_value(ifs, symmetry_autoclose);
+        }
         else if (strcmp("cal_force", word) == 0)
         {
             read_bool(ifs, cal_force);
@@ -2840,6 +2845,7 @@ void Input::Bcast()
     Parallel_Common::bcast_bool(init_vel); // liuyu 2021-07-14
     Parallel_Common::bcast_double(ref_cell_factor);
     Parallel_Common::bcast_double(symmetry_prec); // LiuXh add 2021-08-12, accuracy for symmetry
+    Parallel_Common::bcast_bool(symmetry_autoclose);
     Parallel_Common::bcast_bool(cal_force);
     Parallel_Common::bcast_double(force_thr);
     Parallel_Common::bcast_double(force_thr_ev2);

source/module_io/input.h

@@ -58,6 +58,7 @@ class Input
       1, point group symmetry would be considered*/
     std::string symmetry; 
     double symmetry_prec; // LiuXh add 2021-08-12, accuracy for symmetry
+    bool symmetry_autoclose; // whether to close symmetry automatically when error occurs in symmetry analysis
     int kpar; // ecch pool is for one k point
 
     bool berry_phase; // berry phase calculation

source/module_io/input_conv.cpp

@@ -317,6 +317,7 @@ void Input_Conv::Convert(void)
     Ions_Move_CG::RELAX_CG_THR = INPUT.relax_cg_thr; // pengfei add 2013-09-09
 
     ModuleSymmetry::Symmetry::symm_flag = std::stoi(INPUT.symmetry);
+    ModuleSymmetry::Symmetry::symm_autoclose = INPUT.symmetry_autoclose;
     GlobalV::BASIS_TYPE = INPUT.basis_type;
     GlobalV::KS_SOLVER = INPUT.ks_solver;
     GlobalV::SEARCH_RADIUS = INPUT.search_radius;

source/module_io/test/for_testing_input_conv.h

@@ -114,6 +114,7 @@ int Ions_Move_Basic::out_stru = 0;
 double Ions_Move_CG::RELAX_CG_THR = -1.0;
 std::string Lattice_Change_Basic::fixed_axes = "None";
 int ModuleSymmetry::Symmetry::symm_flag = 0;
+bool ModuleSymmetry::Symmetry::symm_autoclose = false;
 
 Charge_Mixing::Charge_Mixing()
 {

source/module_io/test/input_test.cpp

@@ -82,7 +82,8 @@ TEST_F(InputTest, Default)
         EXPECT_EQ(INPUT.symmetry,"default");
         EXPECT_FALSE(INPUT.init_vel);
         EXPECT_DOUBLE_EQ(INPUT.ref_cell_factor,1.0);
-        EXPECT_DOUBLE_EQ(INPUT.symmetry_prec,1.0e-5);
+        EXPECT_DOUBLE_EQ(INPUT.symmetry_prec, 1.0e-5);
+        EXPECT_FALSE(INPUT.symmetry_autoclose);
         EXPECT_EQ(INPUT.cal_force, 0);
         EXPECT_DOUBLE_EQ(INPUT.force_thr,1.0e-3);
         EXPECT_DOUBLE_EQ(INPUT.force_thr_ev2,0);
@@ -417,7 +418,8 @@ TEST_F(InputTest, Read)
         EXPECT_TRUE(INPUT.search_pbc);
         EXPECT_EQ(INPUT.symmetry,"1");
         EXPECT_FALSE(INPUT.init_vel);
-        EXPECT_DOUBLE_EQ(INPUT.symmetry_prec,1.0e-5);
+        EXPECT_DOUBLE_EQ(INPUT.symmetry_prec, 1.0e-5);
+        EXPECT_FALSE(INPUT.symmetry_autoclose);
         EXPECT_EQ(INPUT.cal_force, 0);
         EXPECT_NEAR(INPUT.force_thr,1.0e-3,1.0e-7);
         EXPECT_DOUBLE_EQ(INPUT.force_thr_ev2,0);

source/module_io/test/input_test_para.cpp

@@ -88,7 +88,8 @@ TEST_F(InputParaTest,Bcast)
         EXPECT_TRUE(INPUT.search_pbc);
         EXPECT_EQ(INPUT.symmetry,"default");
         EXPECT_FALSE(INPUT.init_vel);
-        EXPECT_DOUBLE_EQ(INPUT.symmetry_prec,1.0e-5);
+        EXPECT_DOUBLE_EQ(INPUT.symmetry_prec, 1.0e-5);
+        EXPECT_FALSE(INPUT.symmetry_autoclose);
         EXPECT_EQ(INPUT.cal_force, 0);
         EXPECT_DOUBLE_EQ(INPUT.force_thr,1.0e-3);
         EXPECT_DOUBLE_EQ(INPUT.force_thr_ev2,0);

source/module_io/test/write_input_test.cpp

@@ -52,6 +52,7 @@ TEST_F(write_input,print)
         EXPECT_THAT(output,testing::HasSubstr("symmetry                       1 #the control of symmetry"));
         EXPECT_THAT(output,testing::HasSubstr("init_vel                       0 #read velocity from STRU or not"));
         EXPECT_THAT(output,testing::HasSubstr("symmetry_prec                  1e-05 #accuracy for symmetry"));
+        EXPECT_THAT(output, testing::HasSubstr("symmetry_autoclose             0 #whether to close symmetry automatically when error occurs in symmetry analysis"));
         EXPECT_THAT(output,testing::HasSubstr("nelec                          0 #input number of electrons"));
         EXPECT_THAT(output,testing::HasSubstr("out_mul                        0 # mulliken  charge or not"));
         EXPECT_THAT(output,testing::HasSubstr("noncolin                       0 #using non-collinear-spin"));

source/module_io/write_input.cpp

@@ -66,6 +66,7 @@ void Input::Print(const std::string &fn) const
                                  "symmetry_prec",
                                  symmetry_prec,
                                  "accuracy for symmetry"); // LiuXh add 2021-08-12, accuracy for symmetry
+    ModuleBase::GlobalFunc::OUTP(ofs, "symmetry_autoclose", symmetry_autoclose, "whether to close symmetry automatically when error occurs in symmetry analysis");
     ModuleBase::GlobalFunc::OUTP(ofs, "nelec", nelec, "input number of electrons");
     ModuleBase::GlobalFunc::OUTP(ofs, "out_mul", GlobalV::out_mul, " mulliken  charge or not"); // qifeng add 2019/9/10
     ModuleBase::GlobalFunc::OUTP(ofs, "noncolin", noncolin, "using non-collinear-spin");