Merge branch 'develop' of https://github.com/deepmodeling/abacus-develop into develop

Author: lyb9812

docs/advanced/elec_properties/Mulliken.md

@@ -1,7 +1,7 @@
 # Mulliken Charge Analysis
 
 From version 2.1.0, ABACUS has the function of Mulliken population analysis. The example can be found in [examples/mulliken](https://github.com/deepmodeling/abacus-develop/tree/develop/examples/mulliken). \
-To use this function, set [out_mul](./input-main.md#out_mul) to `1` in the INPUT file. After calculation, there will be an output file named `mulliken.txt` in the output directory. In MD calculations, the output interval is controlled by the keyword [out_interval](./input-main.md#out_interval). In the file, there are contents like (`nspin 1`):
+To use this function, set [out_mul](./input-main.md#out_mul) to `1` in the INPUT file. After calculation, there will be an output file named `mulliken.txt` in the output directory. In MD calculations, the output interval is controlled by the keyword [out_freq_ion](./input-main.md#out_freq_ion). In the file, there are contents like (`nspin 1`):
 
 ```
 STEP: 0

docs/advanced/elec_properties/hs_matrix.md

@@ -45,7 +45,7 @@ The CSR format stores a sparse m × n matrix M in row form using three (one-dime
   - The arrays V and COL_INDEX are of length NNZ, and contain the non-zero values and the column indices of those values respectively.
   - The array ROW_INDEX is of length m + 1 and encodes the index in V and COL_INDEX where the given row starts. This is equivalent to ROW_INDEX[j] encoding the total number of nonzeros above row j. The last element is NNZ , i.e., the fictitious index in V immediately after the last valid index NNZ - 1.
 
-For calculations involving ionic movements, the output frequency of the matrix is controlled by [out_interval](../input_files/input-main.md#out_interval) and [out_app_flag](../input_files/input-main.md#out_app_flag). 
+For calculations involving ionic movements, the output frequency of the matrix is controlled by [out_freq_ion](../input_files/input-main.md#out_freq_ion) and [out_app_flag](../input_files/input-main.md#out_app_flag). 
 
 ## get_s
 We also offer the option of only calculating the overlap matrix without running SCF. For that purpose, in `INPUT` file we need to set the value keyword [calculation](../input_files/input-main.md#calculation) to be `get_s`.

docs/advanced/elec_properties/position_matrix.md

@@ -18,7 +18,7 @@ Each file or each section of the appended file starts with "STEP: " followed by
 
 Each block here contains the matrix for the corresponding cell. There are three columns in each block, giving the matrix elements in x, y, z directions, respectively. There are altogether nbasis * nbasis lines in each block, which emulates the matrix elements.
 
-In molecular dynamics (MD) calculations, if [out_app_flag](../input_files/input-main.md#out_app_flag) is set to true, then `data-rR-tr` is written in an append manner. Otherwise, output files will be put in a separate directory, `matrix`, and named as `$x`_data-rR-tr, where `$x` is the number of MD step. In addition, the output frequency is controlled by [out_interval](../input_files/input-main.md#out_interval). For example, if we are running a 10-step MD with out_interval = 3, then `$x` will be 0, 3, 6, and 9.
+In molecular dynamics (MD) calculations, if [out_app_flag](../input_files/input-main.md#out_app_flag) is set to true, then `data-rR-tr` is written in an append manner. Otherwise, output files will be put in a separate directory, `matrix`, and named as `$x`_data-rR-tr, where `$x` is the number of MD step. In addition, the output frequency is controlled by [out_freq_ion](../input_files/input-main.md#out_freq_ion). For example, if we are running a 10-step MD with out_freq_ion = 3, then `$x` will be 0, 3, 6, and 9.
 
 ## get_s
 We also offer the option of only calculating the position matrix without running SCF. For that purpose, in `INPUT` file we need to set the keyword [calculation](../input_files/input-main.md#calculation) to `get_s`, and [out_mat_r](../input_files/input-main.md#out_mat_r) to `true`.

docs/advanced/input_files/input-main.md

@@ -8,7 +8,7 @@ ecutwfc			10
 scf_nmax		20
 
 basis_type		pw
-relax_nmax      5
+relax_nmax      4
 
 cal_stress      1
 stress_thr      1e-6
@@ -25,18 +25,23 @@ relax_method    lbfgs
 pseudo_dir      ../../../tests/PP_ORB	
 orbital_dir	    ../../../tests/PP_ORB
 
-nspin           2
+nspin           4
+
+out_freq_ion    1
+#out_freq_elec   0 
+
+out_chg         1  # chg.txt g
+out_pot         1  # pot.txt g
+out_wfc_pw      1  # wfkx_pw.txt g
+out_dos         1  # dos.txt g
+out_elf         1  # elf.txt
+out_band        1  # eig.txt
+out_stru        1  # g
+out_bandgap     1
 
-out_chg         1
-out_pot         1
-out_wfc_pw      1
-out_dos         1
-out_stru        1
 out_app_flag    0
 
-out_interval    1
-kpar            2
+kpar            1
 symmetry        -1
 
-#out_freq_elec   2
-#out_band        1
+

examples/relax/pw_output/INPUT

@@ -2,3 +2,4 @@ K_POINTS
 0
 Gamma
 2 2 2 0 0 0
+

examples/relax/pw_output/KPT

@@ -202,6 +202,8 @@ OBJS_CELL=atom_pseudo.o\
     bcast_cell.o\
     read_stru.o\
     read_atom_species.o\
+    sep.o\
+    sep_cell.o\
 
 OBJS_DEEPKS=LCAO_deepks.o\
         deepks_basic.o\
@@ -220,7 +222,7 @@ OBJS_DEEPKS=LCAO_deepks.o\
         deepks_phialpha.o\
         LCAO_deepks_io.o\
         LCAO_deepks_interface.o\
-        
+
 
 OBJS_ELECSTAT=elecstate.o\
     elecstate_energy_terms.o\
@@ -246,6 +248,7 @@ OBJS_ELECSTAT=elecstate.o\
     cal_nelec_nband.o\
     read_pseudo.o\
     cal_wfc.o\
+    pot_sep.o\
 
 OBJS_ELECSTAT_LCAO=elecstate_lcao.o\
       elecstate_lcao_cal_tau.o\
@@ -398,7 +401,7 @@ OBJS_HSOLVER=diago_cg.o\
     diag_const_nums.o\
     diag_hs_para.o\
     diago_pxxxgvx.o\
-    
+
 OBJS_HSOLVER_LCAO=hsolver_lcao.o\
       diago_scalapack.o\
       diago_lapack.o\
@@ -416,7 +419,7 @@ OBJS_HSOLVER_PEXSI=diago_pexsi.o\
       dist_bcd_matrix.o\
       dist_ccs_matrix.o\
       dist_matrix_transformer.o\
-      
+
 OBJS_MD=fire.o\
     langevin.o\
     md_base.o\
@@ -491,7 +494,7 @@ OBJS_RELAXATION=bfgs_basic.o\
     lbfgs.o\
     matrix_methods.o\
     line_search.o\
-    
+
 
 OBJS_SURCHEM=surchem.o\
     H_correction_pw.o\
@@ -691,7 +694,7 @@ OBJS_PARALLEL=parallel_common.o\
     parallel_kpoints.o\
     parallel_reduce.o\
     parallel_device.o
-      
+
 OBJS_SRCPW=H_Ewald_pw.o\
     dnrm2.o\
     VL_in_pw.o\
@@ -755,7 +758,8 @@ OBJS_SRCPW=H_Ewald_pw.o\
     sto_elecond.o\
     sto_dos.o\
     onsite_projector.o\
-    onsite_proj_tools.o
+    onsite_proj_tools.o\
+    VSep_in_pw.o
 
 OBJS_VDW=vdw.o\
     vdwd2_parameters.o\

source/Makefile.Objects

@@ -26,6 +26,8 @@ add_library(
     print_cell.cpp
     read_atom_species.cpp
     k_vector_utils.cpp
+    sep.cpp
+    sep_cell.cpp
 )
 
 if(ENABLE_COVERAGE)

source/source_cell/CMakeLists.txt

@@ -6,15 +6,15 @@
  * 1-1. test if the bravis lattice analysis is right;
  * 1-2. check if matrix-type and vector3-type;
  *     input and optimized lattice vectors are right;
- * 1-3. double-check for  if `gtrans_convert` 
+ * 1-3. double-check for  if `gtrans_convert`
  *     gives the same results as `veccon`;
  * 1-4. check `invmap` function gives the right result;
  * 1-5 test if `gmatrix_convert` and `gmatrix_convert_int`
  *     gives the right result;
  * 2. function: `atom_ordering_new:
  * test the new atom-sort algorithm gives the right result;
  *3. function: `pricell`:
- * test if the number of primitive cells are right, 
+ * test if the number of primitive cells are right,
  * using cases whose space group
  * is different from its point group.
  ***********************************************/
@@ -34,6 +34,10 @@ UnitCell::UnitCell(){}
 UnitCell::~UnitCell(){}
 Magnetism::Magnetism(){}
 Magnetism::~Magnetism() {}
+SepPot::SepPot(){}
+SepPot::~SepPot(){}
+Sep_Cell::Sep_Cell() noexcept {}
+Sep_Cell::~Sep_Cell() noexcept {}
 
 TEST_F(SymmetryTest, AnalySys)
 {
@@ -50,8 +54,8 @@ TEST_F(SymmetryTest, AnalySys)
         int cal_ibrav = symm.real_brav;
         EXPECT_EQ(cal_ibrav, ref_ibrav);
         EXPECT_EQ(cal_point_group, ref_point_group) << "ibrav=" << stru_lib[stru].ibrav;
-        
-        //2. input and optimized lattice, gtrans_convert and veccon 
+
+        //2. input and optimized lattice, gtrans_convert and veccon
         //input lattice
         EXPECT_EQ(symm.s1, ucell.a1);
         EXPECT_EQ(symm.s2, ucell.a2);
@@ -73,7 +77,7 @@ TEST_F(SymmetryTest, AnalySys)
         symm.gtrans_convert(symm.gtrans, gtrans_optconf.data(), symm.nrotk, ucell.latvec, symm.optlat);
         symm.veccon(gtrans_veccon, gtrans_optconf_veccon, symm.nrotk, symm.s1, symm.s2, symm.s3, symm.a1, symm.a2, symm.a3);
         for(int i=0;i<symm.nrotk;++i)
-            EXPECT_EQ(gtrans_optconf[i], ModuleBase::Vector3<double>(gtrans_optconf_veccon[i*3], 
+            EXPECT_EQ(gtrans_optconf[i], ModuleBase::Vector3<double>(gtrans_optconf_veccon[i*3],
             gtrans_optconf_veccon[i*3+1], gtrans_optconf_veccon[i*3+2]));
         delete[] gtrans_veccon;
         delete[] gtrans_optconf_veccon;
@@ -108,7 +112,7 @@ TEST_F(SymmetryTest, AnalySys)
         symm.gmatrix_convert_int(symm.gmatrix, gmatrix_opt, symm.nrotk, ucell.latvec, symm.optlat); //1->2
         symm.gmatrix_convert_int(gmatrix_opt, gmatrix_input_back, symm.nrotk, symm.optlat, ucell.latvec);   //2->3
         symm.gmatrix_convert_int(gmatrix_input_back, gmatrix_opt_back, symm.nrotk, ucell.latvec, symm.optlat); //3->4
-        
+
         symm.gmatrix_convert(symm.gmatrix, kgmatrix_nonint, symm.nrotk, ucell.latvec, ucell.G);
         for (int i=0;i<symm.nrotk;++i)
         {
@@ -130,7 +134,7 @@ TEST_F(SymmetryTest, AnalySys)
             EXPECT_NEAR(gmatrix_opt[i].e31, gmatrix_opt_back[i].e31, DOUBLETHRESHOLD);
             EXPECT_NEAR(gmatrix_opt[i].e23, gmatrix_opt_back[i].e23, DOUBLETHRESHOLD);
             EXPECT_NEAR(gmatrix_opt[i].e32, gmatrix_opt_back[i].e32, DOUBLETHRESHOLD);
-            
+
             ModuleBase::Matrix3 tmpA=symm.optlat.Inverse()*gmatrix_opt[i]*symm.optlat; //A^-1*SA*A
             ModuleBase::Matrix3 tmpB=ucell.latvec.Inverse()*symm.gmatrix[i]*ucell.latvec;//B^-1*SB*B
             ModuleBase::Matrix3 tmpG_int=ucell.G.Inverse()*symm.kgmatrix[i]*ucell.G;//G^-1*SG*G
@@ -168,7 +172,7 @@ TEST_F(SymmetryTest, AnalySys)
             EXPECT_NEAR(tmpA.e23, tmpG.e23, DOUBLETHRESHOLD);
             EXPECT_NEAR(tmpA.e32, tmpG.e32, DOUBLETHRESHOLD);
         }
-            
+
         delete[] gmatrix_input_back;
         delete[] gmatrix_opt;
         delete[] gmatrix_opt_back;
@@ -180,7 +184,7 @@ TEST_F(SymmetryTest, AnalySys)
 
 TEST_F(SymmetryTest, AtomOrderingNew)
 {
-    // the old function `atom_ordering` has bugs 
+    // the old function `atom_ordering` has bugs
     // so here I do not compare with its results
     ModuleSymmetry::Symmetry symm;
     symm.epsilon=1e-5;
@@ -199,7 +203,7 @@ TEST_F(SymmetryTest, AtomOrderingNew)
     }
     //ordering
     symm.test_atom_ordering(new_pos, nat, subindex);
-    //check 
+    //check
     for (int i=0;i<nat-1;++i)
     {
         //x[i]<=x[i+1]
@@ -250,7 +254,7 @@ TEST_F(SymmetryTest, SG_Pricell)
         std::string cal_point_group = symm.pgname;
         std::string ref_space_group = supercell_lib[stru].space_group;
         std::string cal_space_group = symm.spgname;
-        
+
         int ref_ncells = supercell_lib[stru].ibrav;
         EXPECT_EQ(symm.ncell, ref_ncells);
         EXPECT_EQ(cal_point_group, ref_point_group);

source/source_cell/module_symmetry/test/symmetry_test_analysis.cpp

@@ -19,6 +19,10 @@ UnitCell::UnitCell() {}
 UnitCell::~UnitCell() {}
 Magnetism::Magnetism() {}
 Magnetism::~Magnetism() {}
+SepPot::SepPot(){}
+SepPot::~SepPot(){}
+Sep_Cell::Sep_Cell() noexcept {}
+Sep_Cell::~Sep_Cell() noexcept {}
 
 inline std::vector<double> allocate_pos(ModuleSymmetry::Symmetry& symm, UnitCell& ucell)
 {
@@ -46,7 +50,7 @@ TEST_F(SymmetryTest, ForceSymmetry)
 {
     auto check_force = [](stru_& conf, ModuleBase::matrix& force)
     {
-        // 1. check zeros  
+        // 1. check zeros
         for (auto iat : conf.force_zero_iat)
             for (int j = 0; j < 3; ++j)
                 EXPECT_NEAR(force(iat, j), 0.0, DOUBLETHRESHOLD);