Change the file names of output charge densities and kinetic energy densities (#6254)

* add notes for the modified output file names, specify the names in the 3.10-LTS version

* update the CHG file names in input-main.md and change the output file names of charge density

* add explanations for taus1.cube, etc.

* change the file names for 203_PW_OK example

* update read in charge density codes and tests

* fix example 050_PW_CHG_mismatch

* update the file names of charge densities

* update CHG files in support directories

* fix file names in rho_io_test.cpp

* update 203_PW_OK

* update result.ref in 23_NO_KP_OK

Author: mohanchen

docs/advanced/input_files/input-main.md

@@ -608,7 +608,7 @@ These variables are used to control general system parameters.
 - **Description**: This variable is used for both plane wave set and localized orbitals set. It indicates the type of starting density.
 
   - atomic: the density is starting from the summation of the atomic density of single atoms.
-  - file: the density will be read in from a binary file `charge-density.dat` first. If it does not exist, the charge density will be read in from cube files. Besides, when you do `nspin=1` calculation, you only need the density file SPIN1_CHG.cube. However, if you do `nspin=2` calculation, you also need the density file SPIN2_CHG.cube. The density file should be output with these names if you set out_chg = 1 in INPUT file.
+  - file: the density will be read in from a binary file `charge-density.dat` first. If it does not exist, the charge density will be read in from cube files. Besides, when you do `nspin=1` calculation, you only need the density file chgs1.cube. However, if you do `nspin=2` calculation, you also need the density file chgs2.cube. The density file should be output with these names if you set out_chg = 1 in INPUT file.
   - wfc: the density will be calculated by wavefunctions and occupations. Wavefunctions are read in from binary files `wf*.dat` (see [out_wfc_pw](#out_wfc_pw)) while occupations are read in from file `istate.info`.
   - auto: Abacus first attempts to read the density from a file; if not found, it defaults to using atomic density.
 - **Default**: atomic
@@ -742,7 +742,7 @@ These variables are used to control parameters related to input files.
 ### read_file_dir
 
 - **Type**: String
-- **Description**: Indicates the location of files, such as electron density (`SPIN1_CHG.cube`), required as a starting point.
+- **Description**: Location of files, such as the electron density (`chgs1.cube`), required as a starting point.
   - Example: './' implies the files to be read are located in the working directory.
 - **Default**: OUT.$suffix
 
@@ -1602,24 +1602,26 @@ These variables are used to control the output of properties.
 - **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:
-    - 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.
+    - nspin = 1: `chgs1.cube`;
+    - nspin = 2: `chgs1.cube`, and `chgs2.cube`;
+    - nspin = 4: `chgs1.cube`, `chgs2.cube`, `chgs3.cube`, and `chgs4.cube`;
+  Note that by using the Meta-GGA functional, additional files containing the kinetic energy density will be output with the following names:
+    - nspin = 1: `taus1.cube`;
+    - nspin = 2: `taus1.cube`, and `taus2.cube`;
+    - nspin = 4: `taus1.cube`, `taus2.cube`, `taus3.cube`, and `taus4.cube`;
+  - 2: On top of 1, also output the initial charge density files with a suffix name as '_ini', such as `taus1_ini.cube`, etc.
   - -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`.
 
   ---
   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 EXX(exact exchange) calculations, (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).
+  In molecular dynamics simulations, the output frequency is controlled by [out_interval](#out_interval).
 - **Default**: 0 3
+- **NOTE**: In the 3.10-LTS version, the file names are SPIN1_CHG.cube and SPIN1_CHG_INI.cube, etc. 
 
 ### out_pot
 
@@ -1638,6 +1640,7 @@ These variables are used to control the output of properties.
 
   In molecular dynamics calculations, the output frequency is controlled by [out_interval](#out_interval).
 - **Default**: 0
+- **NOTE**: In the 3.10-LTS version, the file names are SPIN1_POT.cube and SPIN1_POT_INI.cube, etc. 
 
 ### out_dm
 
@@ -1651,6 +1654,7 @@ These variables are used to control the output of properties.
     - nspin = 1: `dms1k1_nao.csr`, `dms1k2_nao.csr`, ...;
     - nspin = 2: `dms1k1_nao.csr`... and `dms2k1_nao.csr`... for the two spin channels. 
 - **Default**: False
+- **NOTE**: In the 3.10-LTS version, the file names are SPIN1_DM and SPIN2_DM, etc.
 
 ### out_dm1
 
@@ -1660,6 +1664,7 @@ These variables are used to control the output of properties.
   - nspin = 1: `dmrs1_nao.csr`;
   - nspin = 2: `dmrs1_nao.csr` and `dmrs2_nao.csr` for the two spin channels.
 - **Default**: False
+- **NOTE**: In the 3.10-LTS version, the file names are data-DMR-sparse_SPIN0.csr and data-DMR-sparse_SPIN1.csr, etc.
 
 ### out_wfc_pw
 
@@ -1672,6 +1677,7 @@ These variables are used to control the output of properties.
   - 2: (binary format)
     - non-gamma-only: `wfs1k1_pw.dat` or `wfs1k2_pw.dat`, ....
 - **Default**: 0
+- **NOTE**: In the 3.10-LTS version, the file names are OUT.${suffix}/WAVEFUNC${K}.dat, etc. 
 
 ### out_wfc_lcao
 
@@ -1689,7 +1695,8 @@ These variables are used to control the output of properties.
   The corresponding sequence of the orbitals can be seen in [Basis Set](../pp_orb.md#basis-set).
 
   Also controled by [out_interval](#out_interval) and [out_app_flag](#out_app_flag).
-- **Default**: Flase
+- **Default**: False
+- **NOTE**: In the 3.10-LTS version, the file names are WFC_NAO_GAMMA1_ION1.txt and WFC_NAO_K1_ION1.txt, etc.
 
 ### out_dos
 

source/module_elecstate/module_charge/charge_init.cpp

@@ -28,8 +28,12 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
 {
     ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "init_chg", PARAM.inp.init_chg);
 
+    const int nspin = PARAM.inp.nspin;
+    assert(nspin>0);
+
     std::cout << " START CHARGE      : " << PARAM.inp.init_chg << std::endl;
-    //here we need to set the omega for the charge density
+
+    // we need to set the omega for the charge density
     set_omega(&ucell.omega);
     this->pgrid = &pgrid;
 
@@ -46,17 +50,17 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
         if (ModuleIO::read_rhog(binary.str(), rhopw, rhog))
         {
             GlobalV::ofs_running << " Read electron density from file: " << binary.str() << std::endl;
-            for (int is = 0; is < PARAM.inp.nspin; ++is)
+            for (int is = 0; is < nspin; ++is)
             {
                 rhopw->recip2real(rhog[is], rho[is]);
             }
         }
         else
         {
-            for (int is = 0; is < PARAM.inp.nspin; ++is)
+            for (int is = 0; is < nspin; ++is)
             {
                 std::stringstream ssc;
-                ssc << PARAM.globalv.global_readin_dir << "SPIN" << is + 1 << "_CHG.cube";
+                ssc << PARAM.globalv.global_readin_dir << "chgs" << is + 1 << ".cube";
                 if (ModuleIO::read_vdata_palgrid(pgrid,
                     (PARAM.inp.esolver_type == "sdft" ? GlobalV::RANK_IN_BPGROUP : GlobalV::MY_RANK),
                     GlobalV::ofs_running,
@@ -103,7 +107,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
         {
             const std::string warn_msg
                 = " WARNING: \"init_chg\" is enabled but ABACUS failed to read\n charge density from file.\n"
-                  " Please check if there is SPINX_CHG.cube (X=1,...) or\n {suffix}-CHARGE-DENSITY.restart in the "
+                  " Please check if there is chgsx.cube (x=1,2,etc.) or\n {suffix}-CHARGE-DENSITY.restart in the "
                   "directory.\n";
             std::cout << warn_msg;
             if (PARAM.inp.init_chg == "file")
@@ -121,9 +125,9 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
             {
                 GlobalV::ofs_running << " try to read kinetic energy density from file" << std::endl;
                 // try to read charge from binary file first, which is the same as QE
-                std::vector<std::complex<double>> kin_g_space(PARAM.inp.nspin * this->ngmc, {0.0, 0.0});
+                std::vector<std::complex<double>> kin_g_space(nspin * this->ngmc, {0.0, 0.0});
                 std::vector<std::complex<double>*> kin_g;
-                for (int is = 0; is < PARAM.inp.nspin; is++)
+                for (int is = 0; is < nspin; is++)
                 {
                     kin_g.push_back(kin_g_space.data() + is * this->ngmc);
                 }
@@ -133,14 +137,14 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
                 if (ModuleIO::read_rhog(binary.str(), rhopw, kin_g.data()))
                 {
                     GlobalV::ofs_running << " Read in the kinetic energy density: " << binary.str() << std::endl;
-                    for (int is = 0; is < PARAM.inp.nspin; ++is)
+                    for (int is = 0; is < nspin; ++is)
                     {
                         rhopw->recip2real(kin_g[is], this->kin_r[is]);
                     }
                 }
                 else
                 {
-                    for (int is = 0; is < PARAM.inp.nspin; is++)
+                    for (int is = 0; is < nspin; is++)
                     {
                         std::stringstream ssc;
                         ssc << PARAM.globalv.global_readin_dir << "SPIN" << is + 1 << "_TAU.cube";
@@ -175,16 +179,16 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
         }
     }
 
-    if (PARAM.inp.init_chg == "atomic" || read_error) // mohan add 2007-10-17
+    if (PARAM.inp.init_chg == "atomic" || read_error)
     {
         if (read_error)
         {
             std::cout << " Charge::init_rho: use atomic initialization instead." << std::endl;
         }
-        this->atomic_rho(PARAM.inp.nspin, ucell.omega, rho, strucFac, ucell);
+        this->atomic_rho(nspin, ucell.omega, rho, strucFac, ucell);
     }
 
-    // wenfei 2021-7-29 : initial tau = 3/5 rho^2/3, Thomas-Fermi
+    // initial tau = 3/5 rho^2/3, Thomas-Fermi
     if (XC_Functional::get_ked_flag())
     {
         if (PARAM.inp.init_chg == "atomic" || read_kin_error)
@@ -194,11 +198,11 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
                 std::cout << " Charge::init_rho: init kinetic energy density from rho." << std::endl;
             }
             const double fact = (3.0 / 5.0) * pow(3.0 * ModuleBase::PI * ModuleBase::PI, 2.0 / 3.0);
-            for (int is = 0; is < PARAM.inp.nspin; ++is)
+            for (int is = 0; is < nspin; ++is)
             {
                 for (int ir = 0; ir < this->rhopw->nrxx; ++ir)
                 {
-                    kin_r[is][ir] = fact * pow(std::abs(rho[is][ir]) * PARAM.inp.nspin, 5.0 / 3.0) / PARAM.inp.nspin;
+                    kin_r[is][ir] = fact * pow(std::abs(rho[is][ir]) * nspin, 5.0 / 3.0) / nspin;
                 }
             }
         }
@@ -207,7 +211,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
     // Peize Lin add 2020.04.04
     if (GlobalC::restart.info_load.load_charge && !GlobalC::restart.info_load.load_charge_finish)
     {
-        for (int is = 0; is < PARAM.inp.nspin; ++is)
+        for (int is = 0; is < nspin; ++is)
         {
             try
             {
@@ -217,20 +221,21 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
             {
                 // try to load from the output of `out_chg` 
                 std::stringstream ssc;
-                ssc << PARAM.globalv.global_readin_dir << "SPIN" << is + 1 << "_CHG.cube";
+                ssc << PARAM.globalv.global_readin_dir << "chgs" << is + 1 << ".cube";
                 if (ModuleIO::read_vdata_palgrid(pgrid,
                     (PARAM.inp.esolver_type == "sdft" ? GlobalV::RANK_IN_BPGROUP : GlobalV::MY_RANK),
                     GlobalV::ofs_running,
                     ssc.str(),
                     this->rho[is],
                     ucell.nat))
                 {
-                    GlobalV::ofs_running << " Read in the electron density: " << ssc.str() << std::endl;
+                    GlobalV::ofs_running << " Read in electron density: " << ssc.str() << std::endl;
                 }
             }
         }
         GlobalC::restart.info_load.load_charge_finish = true;
     }
+
 #ifdef __MPI
     this->init_chgmpi();
 #endif
@@ -248,7 +253,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
                 PARAM.globalv.global_readin_dir,
 				GlobalV::KPAR, GlobalV::MY_POOL, GlobalV::MY_RANK, 
                 GlobalV::NPROC_IN_POOL, GlobalV::RANK_IN_POOL,
-				PARAM.inp.nbands, PARAM.inp.nspin, PARAM.globalv.npol,
+				PARAM.inp.nbands, nspin, PARAM.globalv.npol,
 				kv->get_nkstot(),kv->ik2iktot,kv->isk,GlobalV::ofs_running);
     }
 }

source/module_elecstate/test/support/SPIN1_CHG.cube renamed to source/module_elecstate/test/support/chgs1.cube

@@ -165,7 +165,7 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
             for (int is = 0; is < PARAM.inp.nspin; is++)
             {
                 this->pw_rhod->real2recip(this->chr.rho_save[is], this->chr.rhog_save[is]);
-                std::string fn =PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_CHG.cube";
+                std::string fn =PARAM.globalv.global_out_dir + "/chgs" + std::to_string(is + 1) + ".cube";
                 ModuleIO::write_vdata_palgrid(Pgrid,
                                               this->chr.rho_save[is],
                                               is,
@@ -176,9 +176,10 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
                                               &(ucell),
                                               PARAM.inp.out_chg[1],
                                               1);
+
                 if (XC_Functional::get_ked_flag())
                 {
-                    fn =PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_TAU.cube";
+                    fn =PARAM.globalv.global_out_dir + "/taus" + std::to_string(is + 1) + ".cube";
                     ModuleIO::write_vdata_palgrid(Pgrid,
                                                   this->chr.kin_r_save[is],
                                                   is,
@@ -339,7 +340,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
         for (int is = 0; is < PARAM.inp.nspin; is++)
         {
             std::stringstream ss;
-            ss << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_CHG_INI.cube";
+            ss << PARAM.globalv.global_out_dir << "/chgs" << is + 1 << "_ini.cube";
             ModuleIO::write_vdata_palgrid(this->Pgrid,
                                           this->chr.rho[is],
                                           is,

source/module_esolver/esolver_fp.cpp

@@ -91,7 +91,7 @@ class RhoIOTest : public ::testing::Test
 TEST_F(RhoIOTest, Read)
 {
     int is = 0;
-    std::string fn = "./support/SPIN1_CHG.cube";
+    std::string fn = "./support/chgs1.cube";
     int nx = 36;
     int ny = 36;
     int nz = 36;
@@ -117,9 +117,9 @@ TEST_F(RhoIOTest, Write)
     ucell->atoms[0].tau[1] = ModuleBase::Vector3<double>(-0.75, 0.75, 0.75);
     ucell->atoms[0].ncpp.zv = 4;
     Parallel_Grid pgrid(nx, ny, nz, nz, nrxx, nz, 1);
-    ModuleIO::read_vdata_palgrid(pgrid, my_rank, ofs_running, "support/SPIN1_CHG.cube", rho[0], ucell->nat);
+    ModuleIO::read_vdata_palgrid(pgrid, my_rank, ofs_running, "support/chgs1.cube", rho[0], ucell->nat);
     ModuleIO::write_vdata_palgrid(pgrid, rho[0], 0, nspin, 0, "test_write_vdata_palgrid.cube", 0.461002, ucell, 11, 1);
-    EXPECT_EQ(system("diff -q test_write_vdata_palgrid.cube support/SPIN1_CHG.cube"), 0);
+    EXPECT_EQ(system("diff -q test_write_vdata_palgrid.cube support/chgs1.cube"), 0);
 }
 
 TEST_F(RhoIOTest, TrilinearInterpolate)
@@ -130,7 +130,7 @@ TEST_F(RhoIOTest, TrilinearInterpolate)
     int nx_read = 36;
     int ny_read = 36;
     int nz_read = 36;
-    std::ifstream ifs("./support/SPIN1_CHG.cube");
+    std::ifstream ifs("./support/chgs1.cube");
     for (int i = 0; i < 8; ++i)
     {
         ifs.ignore(300, '\n');
@@ -190,7 +190,7 @@ struct CubeIOTest : public ::testing::Test
     std::vector<double> atom_charge;
     std::vector<std::vector<double>> atom_pos;
     std::vector<double> data_read;
-    const std::string fn = "./support/SPIN1_CHG.cube";
+    const std::string fn = "./support/chgs1.cube";
 };
 
 
@@ -236,5 +236,5 @@ TEST_F(CubeIOTest, WriteCube)
 			dx, dy, dz, atom_type, 
 			atom_charge, atom_pos, data_read, 11);
 
-	EXPECT_EQ(system("diff -q test_write.cube ./support/SPIN1_CHG.cube"), 0);
+	EXPECT_EQ(system("diff -q test_write.cube ./support/chgs1.cube"), 0);
 }