Merge branch 'develop' into fft15

Author: A-006

docs/advanced/input_files/input-main.md

@@ -599,7 +599,7 @@ These variables are used to control general system parameters.
   - nao: from numerical atomic orbitals. If they are not enough, other wave functions are initialized with random numbers.
   - nao+random: add small random numbers on numerical atomic orbitals
 
-  > Only the `file` option is useful for the lcao basis set, which is mostly used when [calculation](#calculation) is set to `set_wf` and `get_pchg`. See more details in [out_wfc_lcao](#out_wfc_lcao).
+  > Only the `file` option is useful for the lcao basis set, which is mostly used when [calculation](#calculation) is set to `get_wf` and `get_pchg`. See more details in [out_wfc_lcao](#out_wfc_lcao).
 - **Default**: atomic
 
 ### init_chg
@@ -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
 
@@ -1921,7 +1928,7 @@ The band (KS orbital) energy for each (k-point, spin, band) will be printed in t
 ### if_separate_k
 
 - **Type**: Boolean
-- **Availability**: Only for LCAO, used only when `calculation = get_pchg` and `gamma_only` is turned off.
+- **Availability**: For both PW and LCAO. When `basis_type = pw`, used if `out_pchg` is set. When `basis_type = lcao`, used only when `calculation = get_pchg` and `gamma_only` is turned off.
 - **Description**: Specifies whether to write the partial charge densities for all k-points to individual files or merge them. **Warning**: Enabling symmetry may produce incorrect results due to incorrect k-point weights. Therefore, when calculating partial charge densities, it is strongly recommended to set `symmetry = -1`.
 - **Default**: false
 

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/module_charge/symmetry_rho.cpp

@@ -97,7 +97,7 @@ void Symmetry_rho::psymm(double* rho_part,
         rhotot.resize(rho_basis->nxyz);
         ModuleBase::GlobalFunc::ZEROS(rhotot.data(), rho_basis->nxyz);
     }
-    Pgrid.reduce(rhotot.data(), rho_part);
+    Pgrid.reduce(rhotot.data(), rho_part, false);
 
     // (2)
     if (GlobalV::MY_RANK == 0)

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

@@ -3,6 +3,10 @@
 
 #include "esolver.h"
 
+#ifndef __MPI
+#include <chrono>
+#endif
+
 //! plane wave basis
 #include "module_basis/module_pw/pw_basis.h"
 
@@ -98,6 +102,12 @@ class ESolver_FP: public ESolver
 
     int pw_rho_flag  = false; ///< flag for pw_rho, 0: not initialized, 1: initialized
 
+    //! the start time of scf iteration
+    #ifdef __MPI
+        double iter_time;
+    #else
+        std::chrono::system_clock::time_point iter_time;
+    #endif
 };
 } // namespace ModuleESolver
 

source/module_esolver/esolver_fp.h

@@ -1,12 +1,6 @@
 #ifndef ESOLVER_KS_H
 #define ESOLVER_KS_H
 
-#ifdef __MPI
-#include <mpi.h>
-#else
-#include <chrono>
-#endif
-
 #include <cstring>
 //#include <fstream>
 
@@ -79,13 +73,6 @@ class ESolver_KS : public ESolver_FP
     //! Electronic wavefunctions
     psi::Psi<T>* psi = nullptr;
 
-    //! the start time of scf iteration
-#ifdef __MPI
-    double iter_time;
-#else
-    std::chrono::system_clock::time_point iter_time;
-#endif
-
     std::string basisname;      //! esolver_ks_lcao.cpp
     double esolver_KS_ne = 0.0; //! number of electrons
     double diag_ethr;           //! the threshold for diagonalization