Feature: support control on precison of charge cube output (#4958)

* Feature: support control on precison of charge cube output * [pre-commit.ci lite] apply automatic fixes * change to avoid some by-pass allocation caused segfault cases * update doc * correct the doc --------- Co-authored-by: pre-commit-ci-lite[bot] <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com>
deepmodeling · Aug 16, 2024 · 3d44f1f · 3d44f1f
1 parent 3751dc4
commit 3d44f1f
Show file tree

Hide file tree

Showing 10 changed files with 101 additions and 55 deletions.
diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
@@ -422,8 +422,9 @@
     - [lr\_nstates](#lr_nstates)
     - [abs\_wavelen\_range](#abs_wavelen_range)
     - [out\_wfc\_lr](#out_wfc_lr)
-[back to top](#full-list-of-input-keywords)
+    - [abs\_broadening](#abs_broadening)
 
+[back to top](#full-list-of-input-keywords)
 ## System variables
 
 These variables are used to control general system parameters.
@@ -1511,8 +1512,10 @@ These variables are used to control the output of properties.
 
 ### out_chg
 
-- **Type**: Integer
+- **Type**: Integer \[Integer\](optional)
 - **Description**: 
+
+  The first integer controls whether to output the charge density on real space grids:
   - 1. Output the charge density (in Bohr^-3) on real space grids into the density files in the folder `OUT.${suffix}`. The files are named as:
     - nspin = 1: SPIN1_CHG.cube;
     - nspin = 2: SPIN1_CHG.cube, and SPIN2_CHG.cube;
@@ -1521,13 +1524,16 @@ These variables are used to control the output of properties.
     - nspin = 1: SPIN1_CHG_INI.cube
     - nspin = 2: SPIN1_CHG_INI.cube, and SPIN2_CHG_INI.cube;
     - nspin = 4: SPIN1_CHG_INI.cube, SPIN2_CHG_INI.cube, SPIN3_CHG_INI.cube, and SPIN4_CHG_INI.cube.
+
+  The second integer controls the precision of the charge density output, if not given, will use `3` as default. For purpose restarting from this file and other high-precision involved calculation, recommend to use `10`.
 
+  ---
   The circle order of the charge density on real space grids is: x is the outer loop, then y and finally z (z is moving fastest).
 
   If EXX(exact exchange) is calculated, (i.e. *[dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0/opt_orb* or *[rpa](#rpa)==True*), the Hexx(R) files will be output in the folder `OUT.${suffix}` too, which can be read in NSCF calculation.
 
   In molecular dynamics calculations, the output frequency is controlled by [out_interval](#out_interval).
-- **Default**: 0
+- **Default**: 0 3
 
 ### out_pot
 
@@ -1612,7 +1618,7 @@ These variables are used to control the output of properties.
 
 ### out_band
 
-- **Type**: Boolean Integer(optional)
+- **Type**: Boolean \[Integer\](optional)
 - **Description**: Whether to output the band structure (in eV), optionally output precision can be set by a second parameter, default is 8. For more information, refer to the [band.md](../elec_properties/band.md)
 - **Default**: False
 
@@ -1656,7 +1662,7 @@ These variables are used to control the output of properties.
 
 ### out_mat_hs
 
-- **Type**: Boolean Integer(optional)
+- **Type**: Boolean \[Integer\](optional)
 - **Availability**: Numerical atomic orbital basis
 - **Description**: Whether to print the upper triangular part of the Hamiltonian matrices (in Ry) and overlap matrices for each k point into files in the directory `OUT.${suffix}`. The second number controls precision. For more information, please refer to [hs_matrix.md](../elec_properties/hs_matrix.md#out_mat_hs). Also controled by [out_interval](#out_interval) and [out_app_flag](#out_app_flag).
 - **Default**: False 8

diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp
@@ -122,7 +122,7 @@ void ESolver_FP::after_scf(const int istep)
     if (istep % PARAM.inp.out_interval == 0)
     {
         // 3) write charge density
-        if (PARAM.inp.out_chg)
+        if (PARAM.inp.out_chg[0])
         {
             for (int is = 0; is < GlobalV::NSPIN; is++)
             {
@@ -153,7 +153,7 @@ void ESolver_FP::after_scf(const int istep)
                     this->pw_rhod->nz,
                     this->pelec->eferm.get_efval(is),
                     &(GlobalC::ucell),
-                    3,
+                    PARAM.inp.out_chg[1],
                     1);
                 if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
                 {

diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp
@@ -230,7 +230,7 @@ void ESolver_KS_PW<T, Device>::before_scf(const int istep)
     this->pelec->init_scf(istep, this->sf.strucFac);
 
     //! output the initial charge density
-    if (PARAM.inp.out_chg == 2)
+    if (PARAM.inp.out_chg[0] == 2)
     {
         for (int is = 0; is < GlobalV::NSPIN; is++)
         {
@@ -472,7 +472,7 @@ void ESolver_KS_PW<T, Device>::iter_finish(int& iter)
 
     if (this->out_freq_elec && iter % this->out_freq_elec == 0)
     {
-        if (PARAM.inp.out_chg > 0)
+        if (PARAM.inp.out_chg[0] > 0)
         {
             for (int is = 0; is < GlobalV::NSPIN; is++)
             {

diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp
@@ -220,7 +220,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(const int istep)
     this->pelec->init_scf(istep, this->sf.strucFac);
 
     //! output the initial charge density
-    if (PARAM.inp.out_chg == 2)
+    if (PARAM.inp.out_chg[0] == 2)
     {
         for (int is = 0; is < GlobalV::NSPIN; is++)
         {

diff --git a/source/module_hamilt_lcao/module_dftu/dftu_io.cpp b/source/module_hamilt_lcao/module_dftu/dftu_io.cpp
@@ -52,7 +52,7 @@ void DFTU::output()
 
     //Write onsite.dm
     std::ofstream ofdftu;
-    if(PARAM.inp.out_chg){
+    if(PARAM.inp.out_chg[0]){
       if(GlobalV::MY_RANK == 0){
         ofdftu.open(GlobalV::global_out_dir + "onsite.dm");
       }

diff --git a/source/module_io/read_input_item_output.cpp b/source/module_io/read_input_item_output.cpp
@@ -37,14 +37,25 @@ void ReadInput::item_output()
     }
     {
         Input_Item item("out_chg");
-        item.annotation = ">0 output charge density for selected electron steps";
-        item.reset_value = [](const Input_Item& item, Parameter& para) {
-            if (para.input.calculation == "get_wf" || para.input.calculation == "get_pchg")
+        item.annotation = ">0 output charge density for selected electron steps"
+                          ", second parameter controls the precision, default is 3.";
+        item.read_value = [](const Input_Item& item, Parameter& para) {
+            size_t count = item.get_size();
+            std::vector<int> out_chg(count, -1);
+            std::transform(item.str_values.begin(), item.str_values.end(), out_chg.begin(), [](std::string s) { return std::stoi(s); });
+            // assign non-negative values to para.input.out_chg
+            for (size_t i = 0; i < count; ++i)
             {
-                para.input.out_chg = 1;
+                if (out_chg[i] >= 0)
+                {
+                    para.input.out_chg[i] = out_chg[i];
+                }
             }
         };
-        read_sync_int(input.out_chg);
+        item.reset_value = [](const Input_Item& item, Parameter& para) {
+            para.input.out_chg[0] = (para.input.calculation == "get_wf" || para.input.calculation == "get_pchg") ? 1 : para.input.out_chg[0];
+        };
+        sync_intvec(input.out_chg, 2, 0);
         this->add_item(item);
     }
     {

diff --git a/source/module_io/test/read_input_ptest.cpp b/source/module_io/test/read_input_ptest.cpp
@@ -179,7 +179,8 @@ TEST_F(InputParaTest, ParaRead)
     EXPECT_EQ(param.inp.chg_extrap, "atomic");
     EXPECT_EQ(param.inp.out_freq_elec, 0);
     EXPECT_EQ(param.inp.out_freq_ion, 0);
-    EXPECT_EQ(param.inp.out_chg, 0);
+    EXPECT_EQ(param.inp.out_chg[0], 0);
+    EXPECT_EQ(param.inp.out_chg[1], 3);
     EXPECT_EQ(param.inp.out_dm, 0);
     EXPECT_EQ(param.inp.out_dm1, 0);
     EXPECT_EQ(param.inp.deepks_out_labels, 0);

diff --git a/source/module_io/test_serial/read_input_item_test.cpp b/source/module_io/test_serial/read_input_item_test.cpp
@@ -347,9 +347,9 @@ TEST_F(InputTest, Item_test)
     { // out_chg
         auto it = find_label("out_chg", readinput.input_lists);
         param.input.calculation = "get_wf";
-        param.input.out_chg = 0;
+        param.input.out_chg = {0};
         it->second.reset_value(it->second, param);
-        EXPECT_EQ(param.input.out_chg, 1);
+        EXPECT_EQ(param.input.out_chg[0], 1);
     }
     { // out_pot
         auto it = find_label("out_pot", readinput.input_lists);

diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h
@@ -310,7 +310,7 @@ struct Input_para
                                           ///< 0: output only when converged
     int out_freq_ion = 0;                 ///< the frequency ( >= 0 ) of ionic step to output charge density;
                                           ///< 0: output only when ion steps are finished
-    int out_chg = 0;                      ///< output charge density. 0: no; 1: yes
+    std::vector<int> out_chg = {0, 3};    ///< output charge density. 0: no; 1: yes
     int out_pot = 0;                      ///< yes or no
     int out_wfc_pw = 0;                   ///< 0: no; 1: txt; 2: dat
     bool out_wfc_r = false;               ///< 0: no; 1: yes