Merge branch 'develop' into ppcg-v2

Author: cheerly-pku

.github/workflows/build_test_cmake.yml

@@ -41,7 +41,7 @@ jobs:
             name: "Build without LCAO"
           - tag: gnu
             external_toolchain_args: ""
-            build_args: "-DUSE_ELPA=OFF "
+            build_args: "-DENABLE_ELPA=OFF "
             name: "Build without ELPA"
           - tag: gnu
             external_toolchain_args: ""

CMakeLists.txt

@@ -65,7 +65,7 @@ option(ENABLE_RAPIDJSON "Enable rapid-json usage" OFF)
 option(ENABLE_CNPY "Enable cnpy usage" OFF)
 
 # Options requiring MPI and LCAO
-cmake_dependent_option(ENABLE_ELPA "Enable ELPA for LCAO" ON "ENABLE_LCAO;ENABLE_MPI" OFF)
+cmake_dependent_option(ENABLE_ELPA "Enable ELPA eigensolver" ON "ENABLE_MPI" OFF)
 cmake_dependent_option(ENABLE_LIBRI "Enable LibRI for hybrid functional"
                        OFF "ENABLE_LCAO;ENABLE_MPI" OFF)
 cmake_dependent_option(ENABLE_PEXSI "Enable PEXSI for LCAO" OFF "ENABLE_LCAO;ENABLE_MPI" OFF)
@@ -304,17 +304,15 @@ if(CMAKE_CXX_COMPILER_ID MATCHES Intel)
   set(ENABLE_ABACUS_LIBM OFF)
 endif()
 
-if(ENABLE_LCAO)
-  if(ENABLE_ELPA)
-    find_package(ELPA REQUIRED)
-  endif()
-  if(ENABLE_PEXSI)
-    find_package(PEXSI REQUIRED CONFIG)
-    if(PEXSI_VERSION VERSION_LESS "2.0.0")
-      message(FATAL_ERROR "PEXSI >= 2.0.0 is required")
-    endif()
-    set(CMAKE_CXX_STANDARD 14)
+if(ENABLE_ELPA)
+  find_package(ELPA REQUIRED)
+endif()
+if(ENABLE_PEXSI)
+  find_package(PEXSI REQUIRED CONFIG)
+  if(PEXSI_VERSION VERSION_LESS "2.0.0")
+    message(FATAL_ERROR "PEXSI >= 2.0.0 is required")
   endif()
+  set(CMAKE_CXX_STANDARD 14)
 endif()
 
 if(ENABLE_MPI)

cmake/CollectBuildInfoVars.cmake

@@ -126,13 +126,6 @@ else()
   set(ABACUS_LCAO_ENABLED "False")
 endif()
 
-# Core Math Libraries
-if(ENABLE_LCAO AND ENABLE_ELPA)
-    set(ABACUS_ELPA_VERSION "yes (v${ELPA_VERSION})")
-else()
-    set(ABACUS_ELPA_VERSION "no")
-endif()
-
 # BLAS and FFTW libraries
 if(MKL_FOUND)
   set(ABACUS_BLAS_VENDOR "MKL")
@@ -179,6 +172,13 @@ else()
   endif()
 endif()
 
+# ELPA eigensolver
+if(ENABLE_ELPA)
+    set(ABACUS_ELPA_VERSION "yes (v${ELPA_VERSION})")
+else()
+    set(ABACUS_ELPA_VERSION "no")
+endif()
+
 if(ENABLE_LIBXC AND Libxc_VERSION)
     set(ABACUS_LIBXC_VERSION "yes (v${Libxc_VERSION})")
 elseif(ENABLE_LIBXC)

docs/advanced/input_files/input-main.md

@@ -673,13 +673,22 @@
 ### kpar
 
 - **Type**: Integer
-- **Description**: Divide all processors into kpar groups, and k points will be distributed among each group. The value taken should be less than or equal to the number of k points as well as the number of MPI processes.
+- **Description**: Controls k-point parallelism. The value must be positive and should not exceed either the number of k-points or the number of MPI processes.
+  - For PW calculations, divide all MPI processes into persistent k-point pools. Each pool stores and processes a subset of the k-points.
+  - For LCAO calculations with lapack, genelpa, elpa, or scalapack_gvx, divide the diagonalization work into temporary k-point pools. After diagonalization, the eigenvalues and distributed wavefunctions are restored for all k-points before occupations, density matrices, and output are evaluated.
+  - Multi-process LCAO cusolver uses its own active-GPU distribution and does not use this value to define its k-point layout. Other LCAO eigensolvers do not use the temporary k-point-pool implementation.
 - **Default**: 1
 
 ### bndpar
 
 - **Type**: Integer
-- **Description**: Divide all processors into bndpar groups for SDFT or the BPCG solver. bndpar must be positive, no greater than the number of MPI processes, and kpar * bndpar must divide the number of MPI processes exactly.
+- **Availability**: *([`basis_type`](#basis_type)==pw and [`esolver_type`](#esolver_type)==sdft) or ([`basis_type`](#basis_type)==pw and [`esolver_type`](#esolver_type)==ksdft and [`ks_solver`](#ks_solver)==bpcg)*
+- **Description**: Controls band-group parallelism for PW SDFT and PW KSDFT calculations using the BPCG eigensolver.
+  - Within each k-point pool, divide the MPI processes into bndpar band groups. Each group contains NPROC / (kpar * bndpar) processes when bndpar is greater than 1.
+  - With BPCG, distribute contiguous ranges of global Kohn-Sham bands among the band groups. nbands does not need to be divisible by bndpar, but bndpar cannot exceed a positive nbands. Groups with lower indices receive one additional band when necessary.
+  - In SDFT, distribute stochastic orbitals among the band groups. When the deterministic Kohn-Sham eigensolver is not BPCG, band group 0 calculates the deterministic orbitals and broadcasts them to the other groups.
+  - bndpar must be positive and no greater than the number of MPI processes. When bndpar is greater than 1, kpar * bndpar must divide the number of MPI processes exactly.
+  > Note: For PW calculations on GPU, if the input kpar * bndpar differs from the number of MPI processes, ABACUS automatically sets the effective kpar to NPROC / bndpar.
 - **Default**: 1
 
 ### latname
@@ -4516,9 +4525,9 @@
 ### ocp_set
 
 - **Type**: String
-- **Description**: Fixed occupation weights used when ocp is true. Values are assigned band by band for each k-point, following k-point order. In LCAO RT-TDDFT, the initial ground-state SCF uses its normally determined occupations, and this array is applied only during subsequent real-time propagation steps. The repetition syntax N*x expands to N copies of x.
+- **Description**: Fixed occupation weights used when ocp is true. Values are assigned in band order for each k-point, following k-point order. In LCAO RT-TDDFT, the initial ground-state SCF uses its normally determined occupations, and this array is applied only during subsequent real-time propagation steps. The repetition syntax N*x expands to N copies of x.
   - Example: 1 10*1 0 1 expands to 13 values, with the 12th value equal to 0 and all other values equal to 1.
-  - After expansion, the array length must equal nks * nbands.
+  - After expansion, provide one block of nbands values for each k-point. If nspin is 2, provide all k-point blocks for spin up followed by all k-point blocks for spin down; otherwise, provide one block per k-point.
   - The sum of all weights must equal nelec; otherwise the calculation terminates with an error.
 - **Default**: None
 

docs/parameters.yaml

@@ -100,18 +100,26 @@ parameters:
     category: System variables
     type: Integer
     description: |
-      Divide all processors into kpar groups, and k points will be distributed among each group. The value taken should be less than or equal to the number of k points as well as the number of MPI processes.
+      Controls k-point parallelism. The value must be positive and should not exceed either the number of k-points or the number of MPI processes.
+      * For PW calculations, divide all MPI processes into persistent k-point pools. Each pool stores and processes a subset of the k-points.
+      * For LCAO calculations with lapack, genelpa, elpa, or scalapack_gvx, divide the diagonalization work into temporary k-point pools. After diagonalization, the eigenvalues and distributed wavefunctions are restored for all k-points before occupations, density matrices, and output are evaluated.
+      * Multi-process LCAO cusolver uses its own active-GPU distribution and does not use this value to define its k-point layout. Other LCAO eigensolvers do not use the temporary k-point-pool implementation.
     default_value: "1"
     unit: ""
     availability: ""
   - name: bndpar
     category: System variables
     type: Integer
     description: |
-      Divide all processors into bndpar groups for SDFT or the BPCG solver. bndpar must be positive, no greater than the number of MPI processes, and kpar * bndpar must divide the number of MPI processes exactly.
+      Controls band-group parallelism for PW SDFT and PW KSDFT calculations using the BPCG eigensolver.
+      * Within each k-point pool, divide the MPI processes into bndpar band groups. Each group contains NPROC / (kpar * bndpar) processes when bndpar is greater than 1.
+      * With BPCG, distribute contiguous ranges of global Kohn-Sham bands among the band groups. nbands does not need to be divisible by bndpar, but bndpar cannot exceed a positive nbands. Groups with lower indices receive one additional band when necessary.
+      * In SDFT, distribute stochastic orbitals among the band groups. When the deterministic Kohn-Sham eigensolver is not BPCG, band group 0 calculates the deterministic orbitals and broadcasts them to the other groups.
+      * bndpar must be positive and no greater than the number of MPI processes. When bndpar is greater than 1, kpar * bndpar must divide the number of MPI processes exactly.
+      [NOTE] For PW calculations on GPU, if the input kpar * bndpar differs from the number of MPI processes, ABACUS automatically sets the effective kpar to NPROC / bndpar.
     default_value: "1"
     unit: ""
-    availability: ""
+    availability: (basis_type==pw and esolver_type==sdft) or (basis_type==pw and esolver_type==ksdft and ks_solver==bpcg)
   - name: latname
     category: System variables
     type: String
@@ -2745,9 +2753,9 @@ parameters:
     category: "RT-TDDFT: Real-Time Time-Dependent Density Functional Theory"
     type: String
     description: |
-      Fixed occupation weights used when ocp is true. Values are assigned band by band for each k-point, following k-point order. In LCAO RT-TDDFT, the initial ground-state SCF uses its normally determined occupations, and this array is applied only during subsequent real-time propagation steps. The repetition syntax N*x expands to N copies of x.
+      Fixed occupation weights used when ocp is true. Values are assigned in band order for each k-point, following k-point order. In LCAO RT-TDDFT, the initial ground-state SCF uses its normally determined occupations, and this array is applied only during subsequent real-time propagation steps. The repetition syntax N*x expands to N copies of x.
       * Example: 1 10*1 0 1 expands to 13 values, with the 12th value equal to 0 and all other values equal to 1.
-      * After expansion, the array length must equal nks * nbands.
+      * After expansion, provide one block of nbands values for each k-point. If nspin is 2, provide all k-point blocks for spin up followed by all k-point blocks for spin down; otherwise, provide one block per k-point.
       * The sum of all weights must equal nelec; otherwise the calculation terminates with an error.
     default_value: None
     unit: ""

source/CMakeLists.txt

@@ -687,6 +687,10 @@ target_link_libraries(
     device
     base)
 
+if(ENABLE_ELPA)
+  target_link_libraries(${ABACUS_BIN_NAME} PRIVATE genelpa)
+endif()
+
 if(ENABLE_LCAO)
   target_link_libraries(
     ${ABACUS_BIN_NAME}
@@ -702,9 +706,6 @@ if(ENABLE_LCAO)
   if(ENABLE_LIBRI)
     target_link_libraries(${ABACUS_BIN_NAME} PRIVATE bse)
   endif()
-  if(ENABLE_ELPA)
-    target_link_libraries(${ABACUS_BIN_NAME} PRIVATE genelpa)
-  endif()
   if(ENABLE_PEXSI)
     target_link_libraries(${ABACUS_BIN_NAME} PRIVATE pexsi)
   endif()

source/Makefile.Objects

@@ -552,6 +552,7 @@ OBJS_IO=module_parameter/input_conv.o\
     cal_test.o\
     module_dos/write_dos_pw.o\
     nscf_fermi_surf.o\
+    band_parallel_output.o\
     write_bands.o\
     module_dos/cal_dos.o\
     module_dos/cal_pdos_gamma.o\

source/source_base/parallel_grid.cpp

@@ -422,60 +422,81 @@ void Parallel_Grid::zpiece_to_stogroup(double* zpiece, const int& iz, double* rh
     // ofs_running << "\n iz = " << iz << " Done.";
     return;
 }
+
+// Taoni modified on 2026-08-21, fixed BPCG out_chg MPI_ERR_RANK
 void Parallel_Grid::reduce(double* rhotot, const double* const rhoin, const bool reduce_all_pool) const
 {
-    // ModuleBase::TITLE("Parallel_Grid","reduce");
-
-    // if not the first pool, wait here until processpr 0
-    // send the Barrier command.
+    // POOL_WORLD communicators are disjoint, so inactive pools may return
+    // without skipping a collective required by an active pool.
     if (!reduce_all_pool && GlobalV::MY_POOL != 0)
     {
         return;
     }
 
-    double* zpiece = new double[this->ncxy];
+    assert(rhoin != nullptr);
+    assert(this->nrxx == this->ncxy * this->nczp);
 
-    for (int iz = 0; iz < this->ncz; iz++)
-    {
-        const int znow = iz - this->startz[GlobalV::MY_POOL][GlobalV::RANK_IN_POOL];
-        const int proc = this->whichpro[GlobalV::MY_POOL][iz];
-        const int proc_loc = this->whichpro_loc[GlobalV::MY_POOL][iz]; // Obtain the local processor index in the pool
-        ModuleBase::GlobalFunc::ZEROS(zpiece, this->ncxy);
-        int tag = iz;
-        MPI_Status ierror;
-
-        // Local processor 0 collects data from all other processors in the pool
-        // proc = proc_loc if GlobalV::MY_POOL == 0
-        if (proc_loc == GlobalV::RANK_IN_POOL)
-        {
-            for (int ir = 0; ir < ncxy; ir++)
-            {
-                zpiece[ir] = rhoin[ir * this->nczp + znow];
-            }
-            // Send data to the root of the pool
-            if (GlobalV::RANK_IN_POOL != 0)
-            {
-                MPI_Send(zpiece, ncxy, MPI_DOUBLE, 0, tag, POOL_WORLD);
-            }
-        }
+    int pool_size = 0;
+    MPI_Comm_size(POOL_WORLD, &pool_size);
+
+    // Derive slab ownership from the communicator used below. Precomputed
+    // source ranks may belong to a larger BPCG layout and are not necessarily
+    // valid ranks in the current POOL_WORLD.
+    std::vector<int> local_z_counts(pool_size);
+    MPI_Allgather(&this->nczp, 1, MPI_INT, local_z_counts.data(), 1, MPI_INT, POOL_WORLD);
 
-        // The root of the pool receives data from other processors
-        if (GlobalV::RANK_IN_POOL == 0 && proc_loc != GlobalV::RANK_IN_POOL)
+    // MPI_Gatherv concatenates complete [xy][local_z] buffers by rank.
+    // These displacements therefore address rank blocks, not global z coordinates.
+    std::vector<int> receive_counts(pool_size);
+    std::vector<int> displacements(pool_size, 0);
+    int total_z = 0;
+    for (int rank = 0; rank < pool_size; ++rank)
+    {
+        receive_counts[rank] = local_z_counts[rank] * this->ncxy;
+        if (rank > 0)
         {
-            MPI_Recv(zpiece, ncxy, MPI_DOUBLE, proc_loc, tag, POOL_WORLD, &ierror);
+            displacements[rank] = displacements[rank - 1] + receive_counts[rank - 1];
         }
+        total_z += local_z_counts[rank];
+    }
+    assert(total_z == this->ncz);
 
-        if (GlobalV::RANK_IN_POOL == 0)
+    // Only the pool root needs storage for the gathered global grid.
+    // MPI ignores the receive buffer on all non-root ranks.
+    std::vector<double> gathered_data;
+    if (GlobalV::RANK_IN_POOL == 0)
+    {
+        assert(rhotot != nullptr);
+        gathered_data.resize(this->ncxyz);
+    }
+    MPI_Gatherv(rhoin,
+                this->nrxx,
+                MPI_DOUBLE,
+                gathered_data.data(),
+                receive_counts.data(),
+                displacements.data(),
+                MPI_DOUBLE,
+                0,
+                POOL_WORLD);
+
+    if (GlobalV::RANK_IN_POOL == 0)
+    {
+        // Rank blocks cannot be copied directly to rhotot: each block stores
+        // [xy][local_z], whereas Cube output expects [xy][global_z].
+        // The slab decomposition is contiguous and ordered by POOL_WORLD rank.
+        int global_z_start = 0;
+        for (int rank = 0; rank < pool_size; ++rank)
         {
+            const int local_z = local_z_counts[rank];
             for (int ixy = 0; ixy < this->ncxy; ++ixy)
             {
-                rhotot[ixy * ncz + iz] = zpiece[ixy];
+                for (int iz = 0; iz < local_z; ++iz)
+                {
+                    rhotot[ixy * this->ncz + global_z_start + iz] = gathered_data[displacements[rank] + ixy * local_z + iz];
+                }
             }
+            global_z_start += local_z;
         }
     }
-
-    delete[] zpiece;
-
-    return;
 }
 #endif

source/source_cell/cal_atoms_info.h

@@ -143,6 +143,11 @@ class CalAtomsInfo
         result.nbands = nbands;
         unitcell::cal_nbands(static_cast<int>(result.nelec), result.nlocal, nelec_spin, result.nbands,
                               esolver_type, lspinorb, nspin, basis_type, smearing_method);
+        // Taoni add check for bndpar > nbands on 2026-08-21
+        if (ks_solver == "bpcg" && result.nbands > 0 && bndpar > result.nbands)
+        {
+            ModuleBase::WARNING_QUIT("CalAtomsInfo", "bndpar cannot exceed nbands when ks_solver is bpcg");
+        }
 
         // calculate the number of nbands_local
         result.nbands_l = result.nbands;

source/source_cell/parallel_kpoints.cpp

@@ -113,16 +113,20 @@ void Parallel_Kpoints::set_startpro_pool()
     return;
 }
 
-
 // gather kpoints from all processor pools, only need to be called by the first processor of each pool.
 void Parallel_Kpoints::gatherkvec(const std::vector<ModuleBase::Vector3<double>>& vec_local,
                                   std::vector<ModuleBase::Vector3<double>>& vec_global) const
 {
+    int world_rank = 0;
+    MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
+    // Taoni fix bndpar on 2026-08-21
+    // rank_in_pool is repeated in every band group under BPCG. The first
+    // world rank of the k-point pool instead identifies one unique contributor.
+    const bool is_pool_root = (world_rank == this->startpro_pool[this->my_pool]);
     vec_global.resize(this->nkstot_np, ModuleBase::Vector3<double>(0.0, 0.0, 0.0));
     for (int i = 0; i < this->nks_np; ++i)
     {
-
-        if (this->rank_in_pool == 0)
+        if (is_pool_root)
         {
             vec_global[i + startk_pool[this->my_pool]] = vec_local[i];
         }
@@ -154,7 +158,6 @@ void Parallel_Kpoints::pool_collection(double& value, const double* wk, const in
             {
                 MPI_Status ierror;
                 MPI_Recv(&value, 1, MPI_DOUBLE, this->startpro_pool[pool], ik, MPI_COMM_WORLD, &ierror);
-
             }
         }
         else
@@ -215,7 +218,7 @@ void Parallel_Kpoints::pool_collection_aux(T* value, const V& w, const int& dim,
     T* p = &w.ptr[begin];
     // temprary restrict kpar=1 for NSPIN=2 case for generating_orbitals
     int pool = 0;
-    if (this->nspin != 2) 
+    if (this->nspin != 2)
     {
         pool = this->whichpool[ik];
     }