Feature: support automatic binary output of charge density

docs/advanced/input_files/input-main.md

@@ -1515,17 +1515,17 @@ These variables are used to control the output of properties.
 
 - **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;
     - nspin = 4: SPIN1_CHG.cube, SPIN2_CHG.cube, SPIN3_CHG.cube, and SPIN4_CHG.cube.
-  - 2. On top of 1, also output the initial charge density. The files are named as:
+  - 2: On top of 1, also output the initial charge density. The files are named as:
     - 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.
-
+  - -1: disable the charge density auto-back-up file `{suffix}-CHARGE-DENSITY.restart`, useful for large systems.
+
   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`.
 
   ---

source/Makefile.Objects

@@ -476,7 +476,7 @@ OBJS_IO=input_conv.o\
     print_info.o\
     read_cube.o\
     read_rho.o\
-    read_rhog.o\
+    rhog_io.o\
     read_exit_file.o\
     read_wfc_pw.o\
     restart.o\

source/module_basis/module_pw/pw_distributeg.cpp

@@ -199,4 +199,4 @@ void PW_Basis::get_ig2isz_is2fftixy(
 #endif
     return;
 }
-}
+} // namespace ModulePW

source/module_elecstate/module_charge/charge_init.cpp

@@ -32,7 +32,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout, const ModuleBase::ComplexMa
         // try to read charge from binary file first, which is the same as QE
         // liuyu 2023-12-05
         std::stringstream binary;
-        binary << GlobalV::global_readin_dir << "charge-density.dat";
+        binary << GlobalV::global_readin_dir << PARAM.inp.suffix + "-CHARGE-DENSITY.restart";
         if (ModuleIO::read_rhog(binary.str(), rhopw, rhog))
         {
             GlobalV::ofs_running << " Read in the charge density: " << binary.str() << std::endl;
@@ -150,13 +150,12 @@ void Charge::init_rho(elecstate::efermi& eferm_iout, const ModuleBase::ComplexMa
     }
     if (read_error)
     {
-        std::cout << " WARNING: Failed to read charge density from file. Possible reasons: " << std::endl;
-        std::cout << " - not found: The default directory of SPIN1_CHG.cube is OUT.suffix, \n"
-            "or you must set read_file_dir to a specific directory. " << std::endl;
-        std::cout << " - parameter mismatch: check the previous warning." << std::endl;
+        const std::string warn_msg = " WARNING: \"init_chg\" is enabled but ABACUS failed to read charge density from file.\n"
+                                     " Please check if there is SPINX_CHG.cube (X=1,...) or {suffix}-CHARGE-DENSITY.restart in the directory.\n";
+        std::cout << std::endl << warn_msg;
         if (GlobalV::init_chg == "auto")
         {
-            std::cout << " Use atomic initialization instead." << std::endl;
+            std::cout << " Charge::init_rho: use atomic initialization instead." << std::endl << std::endl;
         }
         else if (GlobalV::init_chg == "file")
         {

source/module_esolver/esolver_fp.cpp

@@ -6,6 +6,7 @@
 #include "module_io/output_log.h"
 #include "module_io/write_elecstat_pot.h"
 #include "module_parameter/parameter.h"
+#include "module_io/rhog_io.h"
 namespace ModuleESolver
 {
 
@@ -122,18 +123,20 @@ void ESolver_FP::after_scf(const int istep)
     if (istep % PARAM.inp.out_interval == 0)
     {
         // 3) write charge density
-        if (PARAM.inp.out_chg[0])
+        if (PARAM.inp.out_chg[0] > 0)
         {
             for (int is = 0; is < GlobalV::NSPIN; is++)
             {
                 double* data = nullptr;
                 if (PARAM.inp.dm_to_rho)
                 {
                     data = this->pelec->charge->rho[is];
+                    this->pw_rhod->real2recip(this->pelec->charge->rho[is], this->pelec->charge->rhog[is]);
                 }
                 else
                 {
                     data = this->pelec->charge->rho_save[is];
+                    this->pw_rhod->real2recip(this->pelec->charge->rho_save[is], this->pelec->charge->rhog_save[is]);
                 }
                 std::string fn = GlobalV::global_out_dir + "/SPIN" + std::to_string(is + 1) + "_CHG.cube";
                 ModuleIO::write_cube(
@@ -178,6 +181,23 @@ void ESolver_FP::after_scf(const int istep)
                 }
             }
         }
+        if (PARAM.inp.out_chg[0] != -1)
+        {
+            std::complex<double>** rhog_tot = (PARAM.inp.dm_to_rho)? this->pelec->charge->rhog : this->pelec->charge->rhog_save;
+            double** rhor_tot = (PARAM.inp.dm_to_rho)? this->pelec->charge->rho : this->pelec->charge->rho_save;
+            for (int is = 0; is < GlobalV::NSPIN; is++)
+            {
+                this->pw_rhod->real2recip(rhor_tot[is], rhog_tot[is]);
+            }
+            ModuleIO::write_rhog(GlobalV::global_out_dir + PARAM.inp.suffix + "-CHARGE-DENSITY.restart",
+                                 GlobalV::GAMMA_ONLY_PW || GlobalV::GAMMA_ONLY_LOCAL,
+                                 this->pw_rhod,
+                                 GlobalV::NSPIN,
+                                 GlobalC::ucell.GT,
+                                 rhog_tot,
+                                 GlobalV::RANK_IN_POOL,
+                                 GlobalV::NPROC_IN_POOL);
+        }
 
         // 4) write potential
         if (PARAM.inp.out_pot == 1 || PARAM.inp.out_pot == 3)

source/module_io/CMakeLists.txt

@@ -13,7 +13,7 @@ list(APPEND objects
     print_info.cpp
     read_cube.cpp
     read_rho.cpp
-    read_rhog.cpp
+    rhog_io.cpp
     read_exit_file.cpp
     read_wfc_pw.cpp
     restart.cpp