Merge branch 'develop' into new

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_hamilt_lcao/hamilt_lcaodft/LCAO_gen_fixedH.cpp

@@ -169,7 +169,7 @@ void LCAO_gen_fixedH::build_ST_new(const char& dtype, const bool& calc_deri, con
 									// according to global2local_row and global2local_col
 									// the last paramete: 1 for Sloc, 2 for Hloc
 									// and 3 for Hloc_fixed.
-									this->LM->set_HSgamma(iw1_all, iw2_all, olm[0], dtype, HSloc);
+                                    this->LM->set_HSgamma(iw1_all, iw2_all, olm[0], HSloc);
 								}
 								else // k point algorithm
 								{
@@ -651,7 +651,7 @@ void LCAO_gen_fixedH::build_Nonlocal_mu_new(double* NLloc, const bool &calc_deri
 											// mohan add 2010-12-20
 											if( nlm_tmp!=0.0 )
 											{
-												this->LM->set_HSgamma(iw1_all,iw2_all,nlm_tmp,'N', NLloc);//N stands for nonlocal.
+                                                this->LM->set_HSgamma(iw1_all, iw2_all, nlm_tmp, NLloc);//N stands for nonlocal.
 											}
 										}
 										else

source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.cpp

@@ -141,8 +141,7 @@ void LCAO_Matrix::allocate_HS_R(const int &nnR)
 void LCAO_Matrix::set_HSgamma(
     const int &iw1_all, // index i for atomic orbital (row)
     const int &iw2_all, // index j for atomic orbital (column)
-    const double &v, // value for matrix element (i,j) 
-    const char &dtype, // type of the matrix
+    const double& v, // value for matrix element (i,j) 
     double* HSloc) //input pointer for store the matrix
 {
     // use iw1_all and iw2_all to set Hloc

source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h

@@ -190,7 +190,7 @@ class LCAO_Matrix
     double* DHloc_fixed_33;
 
 
-    void set_HSgamma(const int &iw1_all, const int &iw2_all, const double &v, const char &dtype, double* HSloc);
+    void set_HSgamma(const int& iw1_all, const int& iw2_all, const double& v, double* HSloc);
     void set_HSk(const int &iw1_all, const int &iw2_all, const std::complex<double> &v, const char &dtype, const int spin = 0);
 
     void set_force (const int& iw1_all, const int& iw2_all, const double& vx, const double& vy, 

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp

@@ -59,13 +59,13 @@ void Veff<OperatorLCAO<double>>::contributeHk(int ik)
 
     if(XC_Functional::get_func_type()==3 || XC_Functional::get_func_type()==5)
     {
-        Gint_inout inout(vr_eff1, vofk_eff1, this->LM, Gint_Tools::job_type::vlocal_meta);
-        this->GG->cal_vlocal(&inout, new_e_iteration);
+        Gint_inout inout(vr_eff1, vofk_eff1, Gint_Tools::job_type::vlocal_meta);
+        this->GG->cal_vlocal(&inout, this->LM, new_e_iteration);
     }
     else
     {
-        Gint_inout inout(vr_eff1, this->LM, Gint_Tools::job_type::vlocal);
-        this->GG->cal_vlocal(&inout, new_e_iteration);
+        Gint_inout inout(vr_eff1, Gint_Tools::job_type::vlocal);
+        this->GG->cal_vlocal(&inout, this->LM, new_e_iteration);
     }
 
     this->new_e_iteration = false;

source/module_hamilt_lcao/module_gint/gint_gamma.h

@@ -33,7 +33,7 @@ class Gint_Gamma : public Gint
     // there is an additional step in calculating vlocal for gamma point
     // namely the redistribution of Hamiltonian from grid to 2D block format
     // hence we have an additional layer outside the unified interface
-	void cal_vlocal(Gint_inout *inout, const bool new_e_iteration);
+    void cal_vlocal(Gint_inout* inout, LCAO_Matrix* lm, const bool new_e_iteration);
 
     //------------------------------------------------------
     // in gint_gamma_env.cpp 
@@ -45,12 +45,18 @@ class Gint_Gamma : public Gint
 
     double***  DM;   //pointer to LOC.DM
 
-    //------------------------------------------------------
-    // in gint_gamma_vl.cpp 
-    //------------------------------------------------------
-    // method for redistributing the Hamiltonian
-    // from grid to 2D format
-    void vl_grid_to_2D(const int lgd, LCAO_Matrix& lm, const bool new_e_iteration);
+    ///------------------------------------------------------
+    /// in gint_gamma_vl.cpp 
+    ///------------------------------------------------------
+    /// method for redistributing the Hamiltonian
+    /// from grid to 2D format
+    /// pass a setter function to customize row/col major and outputs
+    void vl_grid_to_2D(const double* vl_grid,
+        const Parallel_2D& p2d,
+        const int loc_grid_dim,
+        const bool new_e_iteration,
+        double* vl_2d,
+        std::function<void(const int&, const int&, const double&, double*)> setfunc);
 
     ///===============================
     /// Use MPI_Alltoallv to convert a grid distributed matrix

source/module_hamilt_lcao/module_gint/gint_gamma_vl.cpp

@@ -26,7 +26,7 @@ extern "C"
     void Cblacs_pcoord(int icontxt, int pnum, int *prow, int *pcol);
 }
 
-void Gint_Gamma::cal_vlocal(Gint_inout *inout, const bool new_e_iteration)
+void Gint_Gamma::cal_vlocal(Gint_inout* inout, LCAO_Matrix* lm, bool new_e_iteration)
 {
 	const int max_size = this->gridt->max_atom;
 	const int lgd = this->gridt->lgd;
@@ -40,10 +40,10 @@ void Gint_Gamma::cal_vlocal(Gint_inout *inout, const bool new_e_iteration)
         }
 
         this->cal_gint(inout);
+        this->vl_grid_to_2D(pvpR_grid, *lm->ParaV, lgd, new_e_iteration, lm->Hloc.data(),
+            std::bind(&LCAO_Matrix::set_HSgamma, lm, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4));
 
-		this->vl_grid_to_2D(lgd,inout->lm[0], new_e_iteration);
-
-		if (max_size >0 && lgd > 0) delete[] pvpR_grid;
+        if (max_size > 0 && lgd > 0) delete[] pvpR_grid;
 	}
 }
 
@@ -192,7 +192,8 @@ inline int setBufferParameter(
 }
 #endif
 
-void Gint_Gamma::vl_grid_to_2D(const int lgd_now, LCAO_Matrix &lm, const bool new_e_iteration)
+void Gint_Gamma::vl_grid_to_2D(const double* vl_grid, const Parallel_2D& p2d, const int loc_grid_dim, const bool new_e_iteration, double* vl_2d,
+    std::function<void(const int&, const int&, const double&, double*)> setfunc)
 {
     ModuleBase::timer::tick("Gint_Gamma","distri_vl");
     // setup send buffer and receive buffer size
@@ -201,7 +202,7 @@ void Gint_Gamma::vl_grid_to_2D(const int lgd_now, LCAO_Matrix &lm, const bool ne
     {
         ModuleBase::timer::tick("Gint_Gamma","distri_vl_index");
         #ifdef __MPI
-        setBufferParameter(*this->gridt, lm.ParaV->comm_2D, lm.ParaV->blacs_ctxt, lm.ParaV->nb,
+        setBufferParameter(*this->gridt, p2d.comm_2D, p2d.blacs_ctxt, p2d.nb,
                            this->sender_index_size, this->sender_local_index,
                            this->sender_size_process, this->sender_displacement_process,
                            this->sender_size, this->sender_buffer,
@@ -226,11 +227,11 @@ void Gint_Gamma::vl_grid_to_2D(const int lgd_now, LCAO_Matrix &lm, const bool ne
         const int icol=this->sender_local_index[i+1];
         if(irow<=icol)
 		{
-            this->sender_buffer[i/2]=pvpR_grid[irow*lgd_now+icol];
+            this->sender_buffer[i / 2] = vl_grid[irow * loc_grid_dim + icol];
 		}
         else
 		{
-            this->sender_buffer[i/2]=pvpR_grid[icol*lgd_now+irow];
+            this->sender_buffer[i / 2] = vl_grid[icol * loc_grid_dim + irow];
 		}
     }
 
@@ -242,7 +243,7 @@ void Gint_Gamma::vl_grid_to_2D(const int lgd_now, LCAO_Matrix &lm, const bool ne
     #ifdef __MPI
     MPI_Alltoallv(this->sender_buffer, this->sender_size_process, this->sender_displacement_process, MPI_DOUBLE,
                   this->receiver_buffer, this->receiver_size_process,
-					this->receiver_displacement_process, MPI_DOUBLE, lm.ParaV->comm_2D);
+        this->receiver_displacement_process, MPI_DOUBLE, p2d.comm_2D);
     #endif
 
 #ifdef __DEBUG
@@ -252,9 +253,9 @@ void Gint_Gamma::vl_grid_to_2D(const int lgd_now, LCAO_Matrix &lm, const bool ne
     // put local data to H matrix
     for(int i=0; i<this->receiver_index_size; i+=2)
     {
-        const int g_row=this->receiver_global_index[i];
-        const int g_col=this->receiver_global_index[i+1];
-        lm.set_HSgamma(g_row,g_col,this->receiver_buffer[i/2],'L', lm.Hloc.data());
+        const int g_row = this->receiver_global_index[i];
+        const int g_col = this->receiver_global_index[i + 1];
+        setfunc(g_row, g_col, this->receiver_buffer[i / 2], vl_2d);
     }
 
     ModuleBase::timer::tick("Gint_Gamma","distri_vl_value");

source/module_hamilt_lcao/module_gint/gint_tools.h

@@ -6,7 +6,6 @@
 #include "grid_technique.h"
 #include <cstdlib>
 #include "module_elecstate/module_charge/charge.h"
-#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h"
 
 namespace Gint_Tools
 {
@@ -28,7 +27,6 @@ class Gint_inout
         bool isforce;
         bool isstress;
         int ispin;
-		LCAO_Matrix *lm;
 
     //output
         double** rho;
@@ -129,19 +127,17 @@ class Gint_inout
         }
 
 	// vlocal, gamma point
-		Gint_inout(const double* vl_in, LCAO_Matrix *lm_in, Gint_Tools::job_type job_in)
+        Gint_inout(const double* vl_in, Gint_Tools::job_type job_in)
         {
             vl = vl_in;
-            lm = lm_in;
             job = job_in;
         }
 
 	// mGGA vlocal, gamma point
-		Gint_inout(const double* vl_in, const double* vofk_in, LCAO_Matrix *lm_in, Gint_Tools::job_type job_in)
+        Gint_inout(const double* vl_in, const double* vofk_in, Gint_Tools::job_type job_in)
         {
             vl = vl_in;
-			vofk = vofk_in;
-            lm = lm_in;
+            vofk = vofk_in;
             job = job_in;
         }
 };

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");

source/module_ri/RI_2D_Comm.hpp

@@ -122,7 +122,7 @@ void RI_2D_Comm::add_Hexx(
 						if(GlobalV::GAMMA_ONLY_LOCAL)
 							lm.set_HSgamma(iwt0, iwt1,
 								RI::Global_Func::convert<double>(H(iw0_b, iw1_b)) * RI::Global_Func::convert<double>(frac),
-								'L', lm.Hloc.data());
+                                lm.Hloc.data());
 						else
 							lm.set_HSk(iwt0, iwt1,
 								RI::Global_Func::convert<std::complex<double>>(H(iw0_b, iw1_b)) * frac,