Merge branch 'develop' into new-new-setdiag

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()
 

Dockerfile.cuda

@@ -1,9 +1,9 @@
 FROM nvidia/cuda:12.2.0-devel-ubuntu22.04
 
 RUN apt update && apt install -y --no-install-recommends \
-    libopenblas-openmp-dev liblapack-dev libscalapack-mpi-dev libelpa-dev libfftw3-dev libcereal-dev \
+    libopenblas-openmp-dev liblapack-dev libscalapack-mpi-dev libfftw3-dev libcereal-dev \
     libxc-dev libgtest-dev libgmock-dev libbenchmark-dev python3-numpy \
-    bc cmake git g++ make bc time sudo unzip vim wget
+    bc cmake git g++ make bc time sudo unzip vim wget libopenmpi-dev gfortran libtool-bin
 
 ENV GIT_SSL_NO_VERIFY=true TERM=xterm-256color \
     OMPI_ALLOW_RUN_AS_ROOT=1 OMPI_ALLOW_RUN_AS_ROOT_CONFIRM=1 \
@@ -13,6 +13,17 @@ RUN git clone https://github.com/llohse/libnpy.git && \
     cp libnpy/include/npy.hpp /usr/local/include && \
     rm -r libnpy
 
+RUN cd /tmp && \
+    ELPA_VER=2024.05.001 && \
+    wget -q https://elpa.mpcdf.mpg.de/software/tarball-archive/Releases/$ELPA_VER/elpa-$ELPA_VER.tar.gz && \
+    tar xzf elpa-$ELPA_VER.tar.gz  && rm elpa-$ELPA_VER.tar.gz && \
+    cd elpa-$ELPA_VER && \
+    ./configure CXX=mpic++ CFLAGS="-O3 -march=native" FCFLAGS="-O3" LDFLAGS="-L/usr/local/cuda/lib64 -lstdc++" NVCCFLAGS="-arch sm_75 -arch sm_80" --enable-openmp --enable-nvidia-gpu --with-NVIDIA-GPU-compute-capability="sm_70" --with-cuda-path=/usr/local/cuda/ && \
+    make -j`nproc` && \
+    make PREFIX=/usr/local install && \
+    ln -s /usr/local/include/elpa_openmp-$ELPA_VER/elpa /usr/local/include/ && \
+    cd /tmp && rm -rf elpa-$ELPA_VER
+
 ADD https://api.github.com/repos/deepmodeling/abacus-develop/git/refs/heads/develop /dev/null
 
 RUN git clone https://github.com/deepmodeling/abacus-develop.git --depth 1 && \

docs/advanced/input_files/input-main.md

@@ -143,6 +143,7 @@
     - [out\_level](#out_level)
     - [out\_alllog](#out_alllog)
     - [out\_mat\_hs](#out_mat_hs)
+    - [out\_mat\_tk](#out_mat_tk)
     - [out\_mat\_r](#out_mat_r)
     - [out\_mat\_hs2](#out_mat_hs2)
     - [out\_mat\_t](#out_mat_t)
@@ -1514,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`.
 
   ---
@@ -1667,6 +1668,13 @@ These variables are used to control the output of properties.
 - **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
 
+### out_mat_tk
+
+- **Type**: Boolean \[Integer\](optional)
+- **Availability**: Numerical atomic orbital basis
+- **Description**: Whether to print the upper triangular part of the kinetic matrices (in Ry) for each k point into `OUT.${suffix}/data-i-T`, where i is the index of k points (see `OUT.${suffix}/kpoints`). One may optionally provide a second parameter to specify the precision. 
+- **Default**: False \[8\]
+
 ### out_mat_r
 
 - **Type**: Boolean

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_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")
         {