Merge branch 'develop' into elpa_gpu

CMakeLists.txt

@@ -492,7 +492,10 @@ if (ENABLE_CNPY)
     include_directories(${cnpy_INCLUDE_DIR})
   endif()
   include_directories(${cnpy_SOURCE_DIR})
-  target_link_libraries(${ABACUS_BIN_NAME} cnpy)
+
+  # find ZLIB and link
+  find_package(ZLIB REQUIRED)
+  target_link_libraries(${ABACUS_BIN_NAME} cnpy ZLIB::ZLIB)
   add_compile_definitions(__USECNPY)
 endif()
 

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

@@ -477,7 +477,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_extra.cpp

@@ -20,11 +20,7 @@ Charge_Extra::~Charge_Extra()
     }
 }
 
-void Charge_Extra::Init_CE(const int& nspin,
-                           const int& natom,
-                           const double& volume,
-                           const int& nrxx,
-                           const std::string chg_extrap)
+void Charge_Extra::Init_CE(const int& nspin, const int& natom, const int& nrxx, const std::string chg_extrap)
 {
     if (chg_extrap == "none")
     {
@@ -47,13 +43,13 @@ void Charge_Extra::Init_CE(const int& nspin,
         ModuleBase::WARNING_QUIT("Charge_Extra","charge extrapolation method is not available !");
     }
 
-    this->omega_old = volume;
     this->nspin = nspin;
 
-    if (pot_order > 1)
+    if (pot_order > 0)
     {
         delta_rho1.resize(this->nspin, std::vector<double>(nrxx, 0.0));
         delta_rho2.resize(this->nspin, std::vector<double>(nrxx, 0.0));
+        delta_rho3.resize(this->nspin, std::vector<double>(nrxx, 0.0));
     }
 
     if(pot_order == 3)
@@ -110,39 +106,18 @@ void Charge_Extra::extrapolate_charge(
 
     // if(lsda || noncolin) rho2zeta();
 
-    double** rho_atom = new double*[this->nspin];
-    for (int is = 0; is < this->nspin; is++)
-    {
-        rho_atom[is] = new double[chr->rhopw->nrxx];
-    }
-    chr->atomic_rho(this->nspin, omega_old, rho_atom, sf->strucFac, ucell);
-#ifdef _OPENMP
-#pragma omp parallel for collapse(2) schedule(static, 512)
-#endif
-    for (int is = 0; is < this->nspin; is++)
-    {
-        for (int ir = 0; ir < chr->rhopw->nrxx; ir++)
-        {
-            chr->rho[is][ir] -= rho_atom[is][ir];
-            chr->rho[is][ir] *= omega_old;
-        }
-    }
-
     if(rho_extr == 1)
     {
         ofs_running << " NEW-OLD atomic charge density approx. for the potential !" << std::endl;
 
-        if (pot_order > 1)
-        {
 #ifdef _OPENMP
-#pragma omp parallel for collapse(2) schedule(static, 512)
+#pragma omp parallel for collapse(2) schedule(static, 128)
 #endif
-            for (int is = 0; is < this->nspin; is++)
+        for (int is = 0; is < this->nspin; is++)
+        {
+            for (int ir = 0; ir < chr->rhopw->nrxx; ir++)
             {
-                for (int ir = 0; ir < chr->rhopw->nrxx; ir++)
-                {
-                    delta_rho1[is][ir] = chr->rho[is][ir];
-                }
+                chr->rho[is][ir] = delta_rho1[is][ir];
             }
         }
     }
@@ -158,8 +133,6 @@ void Charge_Extra::extrapolate_charge(
         {
             for (int ir = 0; ir < chr->rhopw->nrxx; ir++)
             {
-                delta_rho2[is][ir] = delta_rho1[is][ir];
-                delta_rho1[is][ir] = chr->rho[is][ir];
                 chr->rho[is][ir] = 2 * delta_rho1[is][ir] - delta_rho2[is][ir];
             }
         }
@@ -174,24 +147,25 @@ void Charge_Extra::extrapolate_charge(
         const double one_add_alpha = 1 + alpha;
         const double beta_alpha = beta - alpha;
 
-        std::vector<std::vector<double>> delta_rho3(this->nspin, std::vector<double>(chr->rhopw->nrxx, 0.0));
 #ifdef _OPENMP
 #pragma omp parallel for collapse(2) schedule(static, 64)
 #endif
         for (int is = 0; is < this->nspin; is++)
         {
             for (int ir = 0; ir < chr->rhopw->nrxx; ir++)
             {
-                delta_rho3[is][ir] = delta_rho2[is][ir];
-                delta_rho2[is][ir] = delta_rho1[is][ir];
-                delta_rho1[is][ir] = chr->rho[is][ir];
                 chr->rho[is][ir]
                     = one_add_alpha * delta_rho1[is][ir] + beta_alpha * delta_rho2[is][ir] - beta * delta_rho3[is][ir];
             }
         }
     }
 
     sf->setup_structure_factor(&ucell, chr->rhopw);
+    double** rho_atom = new double*[this->nspin];
+    for (int is = 0; is < this->nspin; is++)
+    {
+        rho_atom[is] = new double[chr->rhopw->nrxx];
+    }
     chr->atomic_rho(this->nspin, ucell.omega, rho_atom, sf->strucFac, ucell);
 #ifdef _OPENMP
 #pragma omp parallel for collapse(2) schedule(static, 512)
@@ -205,8 +179,6 @@ void Charge_Extra::extrapolate_charge(
         }
     }
 
-    omega_old = ucell.omega;
-
     for (int is = 0; is < this->nspin; is++)
     {
         delete[] rho_atom[is];
@@ -296,4 +268,41 @@ void Charge_Extra::update_all_dis(const UnitCell& ucell)
         assert(iat == ucell.nat);
     }
     return;
+}
+
+void Charge_Extra::update_delta_rho(const UnitCell& ucell, const Charge* chr, const Structure_Factor* sf)
+{
+    if (pot_order == 0)
+    {
+        return;
+    }
+
+    // obtain the difference between chr->rho and atomic_rho
+    double** rho_atom = new double*[this->nspin];
+    for (int is = 0; is < this->nspin; is++)
+    {
+        rho_atom[is] = new double[chr->rhopw->nrxx];
+    }
+    chr->atomic_rho(this->nspin, ucell.omega, rho_atom, sf->strucFac, ucell);
+
+#ifdef _OPENMP
+#pragma omp parallel for collapse(2) schedule(static, 512)
+#endif
+    for (int is = 0; is < this->nspin; is++)
+    {
+        for (int ir = 0; ir < chr->rhopw->nrxx; ir++)
+        {
+            delta_rho3[is][ir] = delta_rho2[is][ir];
+            delta_rho2[is][ir] = delta_rho1[is][ir];
+            delta_rho1[is][ir] = chr->rho[is][ir] - rho_atom[is][ir];
+            delta_rho1[is][ir] *= ucell.omega;
+        }
+    }
+
+    for (int is = 0; is < this->nspin; is++)
+    {
+        delete[] rho_atom[is];
+    }
+    delete[] rho_atom;
+    return;
 }

source/module_elecstate/module_charge/charge_extra.h

@@ -47,15 +47,10 @@ class Charge_Extra
      *
      * @param nspin the number of spins
      * @param natom the number of atoms
-     * @param volume the volume of the cell
      * @param nrxx the number of grids
      * @param chg_extrap the charge extrapolation method
      */
-    void Init_CE(const int& nspin,
-                 const int& natom,
-                 const double& volume,
-                 const int& nrxx,
-                 const std::string chg_extrap);
+    void Init_CE(const int& nspin, const int& natom, const int& nrxx, const std::string chg_extrap);
 
     /**
      * @brief charge extrapolation method
@@ -87,11 +82,19 @@ class Charge_Extra
      */
     void update_all_dis(const UnitCell& ucell);
 
+    /**
+     * @brief update the difference of charge density
+     *
+     * @param ucell the cell information
+     * @param chr the charge density
+     * @param sf the structure factor
+     */
+    void update_delta_rho(const UnitCell& ucell, const Charge* chr, const Structure_Factor* sf);
+
   private:
     int istep = 0; ///< the current step
     int pot_order; ///< the specified charge extrapolation method
     int rho_extr;  ///< the actually used method
-    double omega_old; ///< the old volume of the last step
     int nspin;        ///< the number of spins
 
     ModuleBase::Vector3<double>* dis_old1 = nullptr; ///< dis_old2 = pos_old1 - pos_old2
@@ -100,6 +103,7 @@ class Charge_Extra
 
     std::vector<std::vector<double>> delta_rho1; ///< the last step difference of rho and atomic_rho
     std::vector<std::vector<double>> delta_rho2; ///< the second last step difference of rho and atomic_rho
+    std::vector<std::vector<double>> delta_rho3; ///< the third last step difference of rho and atomic_rho
 
     double alpha; ///< parameter used in the second order extrapolation
     double beta;  ///< parameter used in the second order extrapolation

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_elecstate/test/charge_extra_test.cpp

@@ -136,7 +136,7 @@ TEST_F(ChargeExtraTest, InitCEWarningQuit)
 {
     GlobalV::chg_extrap = "wwww";
     testing::internal::CaptureStdout();
-    EXPECT_EXIT(CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap),
+    EXPECT_EXIT(CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap),
                 ::testing::ExitedWithCode(0),
                 "");
     std::string output = testing::internal::GetCapturedStdout();
@@ -146,21 +146,21 @@ TEST_F(ChargeExtraTest, InitCEWarningQuit)
 TEST_F(ChargeExtraTest, InitCECase1)
 {
     GlobalV::chg_extrap = "none";
-    CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap);
+    CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap);
     EXPECT_EQ(CE.pot_order, 0);
 }
 
 TEST_F(ChargeExtraTest, InitCECase2)
 {
     GlobalV::chg_extrap = "atomic";
-    CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap);
+    CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap);
     EXPECT_EQ(CE.pot_order, 1);
 }
 
 TEST_F(ChargeExtraTest, InitCECase3)
 {
     GlobalV::chg_extrap = "first-order";
-    CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap);
+    CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap);
     EXPECT_EQ(CE.pot_order, 2);
     EXPECT_NE(CE.delta_rho1.size(), 0);
     EXPECT_NE(CE.delta_rho2.size(), 0);
@@ -169,7 +169,7 @@ TEST_F(ChargeExtraTest, InitCECase3)
 TEST_F(ChargeExtraTest, InitCECase4)
 {
     GlobalV::chg_extrap = "second-order";
-    CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap);
+    CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap);
     EXPECT_EQ(CE.pot_order, 3);
     EXPECT_DOUBLE_EQ(CE.alpha, 1.0);
     EXPECT_DOUBLE_EQ(CE.beta, 0.0);
@@ -183,7 +183,7 @@ TEST_F(ChargeExtraTest, InitCECase4)
 TEST_F(ChargeExtraTest, ExtrapolateChargeCase1)
 {
     GlobalV::chg_extrap = "second-order";
-    CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap);
+    CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap);
     CE.istep = 0;
     CE.pot_order = 3;
 
@@ -205,7 +205,7 @@ TEST_F(ChargeExtraTest, ExtrapolateChargeCase1)
 TEST_F(ChargeExtraTest, ExtrapolateChargeCase2)
 {
     GlobalV::chg_extrap = "second-order";
-    CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap);
+    CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap);
     CE.istep = 1;
     CE.pot_order = 3;
 
@@ -227,7 +227,7 @@ TEST_F(ChargeExtraTest, ExtrapolateChargeCase2)
 TEST_F(ChargeExtraTest, ExtrapolateChargeCase3)
 {
     GlobalV::chg_extrap = "second-order";
-    CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap);
+    CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap);
     CE.istep = 2;
     CE.pot_order = 3;
 
@@ -249,7 +249,7 @@ TEST_F(ChargeExtraTest, ExtrapolateChargeCase3)
 TEST_F(ChargeExtraTest, ExtrapolateChargeCase4)
 {
     GlobalV::chg_extrap = "second-order";
-    CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap);
+    CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap);
     CE.istep = 3;
 
     GlobalV::ofs_running.open("log");
@@ -271,7 +271,7 @@ TEST_F(ChargeExtraTest, ExtrapolateChargeCase4)
 TEST_F(ChargeExtraTest, UpdateAllDis)
 {
     GlobalV::chg_extrap = "second-order";
-    CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap);
+    CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap);
     CE.istep = 3;
     for (int i = 0; i < ucell->nat; ++i)
     {
@@ -293,7 +293,7 @@ TEST_F(ChargeExtraTest, UpdateAllDis)
 TEST_F(ChargeExtraTest, FindAlphaAndBeta)
 {
     GlobalV::chg_extrap = "second-order";
-    CE.Init_CE(GlobalV::NSPIN, ucell->nat, ucell->omega, charge.rhopw->nrxx, GlobalV::chg_extrap);
+    CE.Init_CE(GlobalV::NSPIN, ucell->nat, charge.rhopw->nrxx, GlobalV::chg_extrap);
     CE.istep = 3;
     for (int i = 0; i < ucell->nat; ++i)
     {