diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 03cfdc35e96..51735fe0d61 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -231,6 +231,10 @@ OBJS_CELL=atom_pseudo.o\ mdcell.o\ cif_io.o\ ucell_io.o\ + read_cube.o\ + write_cube.o\ + write_pao.o\ + output_log.o\ OBJS_DEEPKS=lcao_deepks.o\ deepks_basic.o\ @@ -626,7 +630,6 @@ OBJS_IO=module_parameter/input_conv.o\ module_bessel/numerical_basis_output.o\ output.o\ module_output/print_info.o\ - module_output/read_cube.o\ module_wf/read_wfc_pw.o\ module_wf/read_wf2rho_pw.o\ module_restart/restart.o\ @@ -645,9 +648,7 @@ OBJS_IO=module_parameter/input_conv.o\ module_wannier/to_w90_pw_setup.o\ module_wannier/fr_overlap.o\ module_unk/unk_overlap_pw.o\ - module_output/write_pao.o\ module_wf/write_wfc_pw.o\ - module_output/write_cube.o\ module_elf/write_elf.o\ module_dipole/write_dipole.o\ module_current/td_current_io.o\ @@ -655,7 +656,6 @@ OBJS_IO=module_parameter/input_conv.o\ td_efield_io.o\ td_vector_pot_io.o\ module_chgpot/write_libxc_r.o\ - module_output/output_log.o\ module_hs/output_mat_sparse.o\ module_ctrl/ctrl_scf_lcao.o\ module_ctrl/ctrl_runner_lcao.o\ @@ -692,7 +692,7 @@ OBJS_IO=module_parameter/input_conv.o\ module_hs/cal_plpr.o\ OBJS_IO_LCAO=module_hs/cal_r_overlap_r.o\ - module_output/write_orb_info.o\ + write_orb_info.o\ module_dos/write_dos_lcao.o\ module_energy/write_proj_band_lcao.o\ module_energy/write_eig_occ.o\ @@ -797,7 +797,7 @@ OBJS_PARALLEL=parallel_common.o\ para_pw_world.o\ para_diag_world.o\ para_rgrid_world.o\ - para_bgroup_world.o\ + para_bdiff_ksame_world.o\ para_matrix_world.o\ para_mpi_func.o\ para_setup.o\ diff --git a/source/source_base/CMakeLists.txt b/source/source_base/CMakeLists.txt index e2ccaaf1732..79bf625ae42 100644 --- a/source/source_base/CMakeLists.txt +++ b/source/source_base/CMakeLists.txt @@ -79,7 +79,7 @@ add_library( module_parallel/para_pw_world.cpp module_parallel/para_diag_world.cpp module_parallel/para_rgrid_world.cpp - module_parallel/para_bgroup_world.cpp + module_parallel/para_bdiff_ksame_world.cpp module_parallel/para_matrix_world.cpp module_parallel/para_mpi_func.cpp module_parallel/para_setup.cpp diff --git a/source/source_base/module_parallel/para_bdiff_ksame_world.cpp b/source/source_base/module_parallel/para_bdiff_ksame_world.cpp new file mode 100644 index 00000000000..147ae3f347c --- /dev/null +++ b/source/source_base/module_parallel/para_bdiff_ksame_world.cpp @@ -0,0 +1,33 @@ +#include "para_bdiff_ksame_world.h" + +namespace Parallel +{ + +ParaBdiffKsameWorld::ParaBdiffKsameWorld() + : ParaWorld("bdiff_ksame"), my_bndgroup_(0), nbndgroup_(1) +{ +} + +#ifdef __MPI +ParaBdiffKsameWorld::ParaBdiffKsameWorld(const MPI_Comm& intra_comm, const MPI_Comm& inter_comm, int nbndgroup) + : ParaWorld("bdiff_ksame", intra_comm), inter_comm_(inter_comm), nbndgroup_(nbndgroup) +{ + if (inter_comm != MPI_COMM_NULL) + { + MPI_Comm_rank(inter_comm, &my_bndgroup_); + } +} +#endif + +void ParaBdiffKsameWorld::reduce_across_bdiff_ksame(double& value) const +{ +#ifdef __MPI + if (inter_comm_ == MPI_COMM_NULL || nbndgroup_ <= 1) + { + return; + } + MPI_Allreduce(MPI_IN_PLACE, &value, 1, MPI_DOUBLE, MPI_SUM, inter_comm_); +#endif +} + +} // namespace Parallel diff --git a/source/source_base/module_parallel/para_bdiff_ksame_world.h b/source/source_base/module_parallel/para_bdiff_ksame_world.h new file mode 100644 index 00000000000..2584a491534 --- /dev/null +++ b/source/source_base/module_parallel/para_bdiff_ksame_world.h @@ -0,0 +1,87 @@ +#ifndef PARA_BDIFF_KSAME_WORLD_H +#define PARA_BDIFF_KSAME_WORLD_H + +#include "para_world.h" + +namespace Parallel +{ + +/** + * @brief bdiff_ksame parallel domain: band-group communication topology + * inside one k-pool. + * + * Self-contained replacement for INT_BGROUP + BP_WORLD + + * GlobalV::MY_BNDGROUP/NPROC_IN_BNDGROUP/RANK_IN_BPGROUP. + * + * The domain has two communicators: + * - intra: INT_BGROUP (bsame_kdiff; same band group, different k/pw) + * - inter: BP_WORLD (bdiff_ksame; different band groups, same k) + * + * Tests only need this header. + */ +class ParaBdiffKsameWorld : public ParaWorld +{ +public: + /** + * @brief Construct a serial domain (single band group). + */ + ParaBdiffKsameWorld(); + +#ifdef __MPI + /** + * @brief Construct a domain from intra and inter communicators. + * + * @param[in] intra_comm intra-group communicator (e.g. INT_BGROUP) + * @param[in] inter_comm inter-group communicator (e.g. BP_WORLD) + * @param[in] nbndgroup number of band groups + */ + ParaBdiffKsameWorld(const MPI_Comm& intra_comm, const MPI_Comm& inter_comm, int nbndgroup); +#endif + + /// Band group index of this process. + int my_bndgroup() const { return my_bndgroup_; } + + /// Number of band groups. + int nbndgroup() const { return nbndgroup_; } + + /// Rank within the band group (alias for rank()). + int rank_in_bpgroup() const { return rank(); } + + /// Number of processes in the band group (alias for size()). + int nproc_in_bndgroup() const { return size(); } + +#ifdef __MPI + /// Inter-group communicator (BP_WORLD / bdiff_ksame equivalent). + MPI_Comm inter_comm() const { return inter_comm_; } +#endif + + /** + * @brief Sum a scalar across the band groups of this k-pool. + * + * Band-parallel eigensolvers (bpcg) shard the band range across the + * BNDPAR band groups of a k-pool: every process only accumulates the + * partial sum over its own band window. This reduction combines those + * partial sums on the bdiff_ksame (BP_WORLD) communicator, which links + * the same rank position of every band group inside one k-pool, so each + * band window contributes exactly once. + * + * It must run BEFORE ParaKmeshWorld::reduce_across_pools so that the + * k-pool reduction receives one complete per-k-pool partial sum. + * No-op when there is only a single band group. + * + * @param[in,out] value local partial sum, overwritten with the + * k-pool-wide total + */ + void reduce_across_bdiff_ksame(double& value) const; + +private: + int my_bndgroup_ = 0; + int nbndgroup_ = 1; +#ifdef __MPI + MPI_Comm inter_comm_ = MPI_COMM_NULL; +#endif +}; + +} // namespace Parallel + +#endif // PARA_BDIFF_KSAME_WORLD_H diff --git a/source/source_base/module_parallel/para_bgroup_world.cpp b/source/source_base/module_parallel/para_bgroup_world.cpp deleted file mode 100644 index 83493253757..00000000000 --- a/source/source_base/module_parallel/para_bgroup_world.cpp +++ /dev/null @@ -1,22 +0,0 @@ -#include "para_bgroup_world.h" - -namespace Parallel -{ - -ParaBgroupWorld::ParaBgroupWorld() - : ParaWorld("bdiff_ksame"), my_bndgroup_(0), nbndgroup_(1) -{ -} - -#ifdef __MPI -ParaBgroupWorld::ParaBgroupWorld(const MPI_Comm& intra_comm, const MPI_Comm& inter_comm, int nbndgroup) - : ParaWorld("bdiff_ksame", intra_comm), inter_comm_(inter_comm), nbndgroup_(nbndgroup) -{ - if (inter_comm != MPI_COMM_NULL) - { - MPI_Comm_rank(inter_comm, &my_bndgroup_); - } -} -#endif - -} // namespace Parallel diff --git a/source/source_base/module_parallel/para_bgroup_world.h b/source/source_base/module_parallel/para_bgroup_world.h deleted file mode 100644 index e2d82cc99b9..00000000000 --- a/source/source_base/module_parallel/para_bgroup_world.h +++ /dev/null @@ -1,67 +0,0 @@ -#ifndef PARA_BGROUP_WORLD_H -#define PARA_BGROUP_WORLD_H - -#include "para_world.h" - -namespace Parallel -{ - -/** - * @brief bgroup parallel domain: band group communication topology. - * - * Self-contained replacement for INT_BGROUP + BP_WORLD + - * GlobalV::MY_BNDGROUP/NPROC_IN_BNDGROUP/RANK_IN_BPGROUP. - * - * The band group domain has two communicators: - * - intra: INT_BGROUP (same band group, different k/pw) - * - inter: BP_WORLD (different band groups, same k) - * - * Tests only need this header. - */ -class ParaBgroupWorld : public ParaWorld -{ -public: - /** - * @brief Construct a serial bgroup domain (single band group). - */ - ParaBgroupWorld(); - -#ifdef __MPI - /** - * @brief Construct a bgroup domain from intra and inter communicators. - * - * @param[in] intra_comm intra-group communicator (e.g. INT_BGROUP) - * @param[in] inter_comm inter-group communicator (e.g. BP_WORLD) - * @param[in] nbndgroup number of band groups - */ - ParaBgroupWorld(const MPI_Comm& intra_comm, const MPI_Comm& inter_comm, int nbndgroup); -#endif - - /// Band group index of this process. - int my_bndgroup() const { return my_bndgroup_; } - - /// Number of band groups. - int nbndgroup() const { return nbndgroup_; } - - /// Rank within the band group (alias for rank()). - int rank_in_bpgroup() const { return rank(); } - - /// Number of processes in the band group (alias for size()). - int nproc_in_bndgroup() const { return size(); } - -#ifdef __MPI - /// Inter-group communicator (BP_WORLD equivalent). - MPI_Comm inter_comm() const { return inter_comm_; } -#endif - -private: - int my_bndgroup_ = 0; - int nbndgroup_ = 1; -#ifdef __MPI - MPI_Comm inter_comm_ = MPI_COMM_NULL; -#endif -}; - -} // namespace Parallel - -#endif // PARA_BGROUP_WORLD_H diff --git a/source/source_base/module_parallel/para_bridge.cpp b/source/source_base/module_parallel/para_bridge.cpp index c2c009b1580..37369202361 100644 --- a/source/source_base/module_parallel/para_bridge.cpp +++ b/source/source_base/module_parallel/para_bridge.cpp @@ -2,6 +2,7 @@ #include "para_tag.h" #ifdef __MPI +#include "source_base/global_variable.h" #include "source_base/parallel_comm.h" #endif @@ -19,4 +20,63 @@ ParaWorld make_pw_world() #endif } +// Reduce-only overload: no k-point distribution data. +ParaKmeshWorld make_kmesh_world() +{ +#ifdef __MPI + int mpi_initialized = 0; + MPI_Initialized(&mpi_initialized); + // Any distributed layout (k pools or band groups) may need the + // world-wide max/min reductions, so build the MPI domain whenever + // more than one process is running. The sum reduction no-ops for + // kpar <= 1 on its own. + if (mpi_initialized && GlobalV::NPROC > 1) + { + // Build from globals but skip distribute_kpoints (nkstot=0). + return ParaKmeshWorld(MPI_COMM_WORLD, GlobalV::KPAR, GlobalV::MY_POOL, 0, 1); + } +#endif + return ParaKmeshWorld(); +} + +// Temporary bridge: construct a kmesh-domain ParaKmeshWorld from the old +// globals. Delete this file once ParaCollection is wired into driver init. +ParaKmeshWorld make_kmesh_world(int nkstot, int nspin) +{ +#ifdef __MPI + // Fall back to a serial single-pool domain when MPI is not initialized + // (e.g. unit tests linked against the MPI-compiled base library) or when + // there is only one k-pool, so that no MPI call is made on an unset + // communicator. + int mpi_initialized = 0; + MPI_Initialized(&mpi_initialized); + if (mpi_initialized && GlobalV::KPAR > 1) + { + return ParaKmeshWorld(MPI_COMM_WORLD, GlobalV::KPAR, GlobalV::MY_POOL, + nkstot, nspin); + } +#endif + return ParaKmeshWorld(nkstot, nspin); +} + +// Temporary bridge: construct a bdiff_ksame-domain ParaBdiffKsameWorld from +// the old globals. Delete this file once ParaCollection is wired into driver init. +ParaBdiffKsameWorld make_bdiff_ksame_world() +{ +#ifdef __MPI + int mpi_initialized = 0; + MPI_Initialized(&mpi_initialized); + // NPROC_IN_BNDGROUP stays 0 until divide_pools has run, which also + // guards unit tests that link the MPI base library without a layout. + if (mpi_initialized && INT_BGROUP != MPI_COMM_NULL && BP_WORLD != MPI_COMM_NULL + && GlobalV::NPROC_IN_BNDGROUP > 1) + { + int nbndgroup = 1; + MPI_Comm_size(BP_WORLD, &nbndgroup); + return ParaBdiffKsameWorld(INT_BGROUP, BP_WORLD, nbndgroup); + } +#endif + return ParaBdiffKsameWorld(); +} + } // namespace Parallel diff --git a/source/source_base/module_parallel/para_bridge.h b/source/source_base/module_parallel/para_bridge.h index c0df2a61946..73319de8f54 100644 --- a/source/source_base/module_parallel/para_bridge.h +++ b/source/source_base/module_parallel/para_bridge.h @@ -1,6 +1,8 @@ #ifndef PARA_BRIDGE_H #define PARA_BRIDGE_H +#include "para_bdiff_ksame_world.h" +#include "para_kmesh_world.h" #include "para_world.h" namespace Parallel @@ -16,6 +18,43 @@ namespace Parallel */ ParaWorld make_pw_world(); +/** + * @brief Temporary bridge: construct a kmesh-domain ParaKmeshWorld from + * the old globals KP_WORLD / GlobalV::KPAR (MPI) or as a serial domain + * (non-MPI). + * + * Falls back to a serial single-pool domain when MPI is not initialized + * (e.g. unit tests linked against the MPI-compiled base library) or when + * there is only one k-point pool, so that no MPI call is made on an + * unset communicator. + * + * @param[in] nkstot total number of k-points (without spin) + * @param[in] nspin number of spin components + */ +ParaKmeshWorld make_kmesh_world(int nkstot, int nspin); + +/** + * @brief Reduce-only overload: construct a kmesh domain for call sites + * that only need cross-pool reduction (reduce_across_pools etc.) and + * have no k-point information to pass. + * + * The k-point distribution data (nks_pool_, whichpool_, ...) is left + * empty; calling pool_collection / gather_kvec on the returned object + * is invalid. Use the (nkstot, nspin) overload when those are needed. + */ +ParaKmeshWorld make_kmesh_world(); + +/** + * @brief Temporary bridge: construct a bdiff_ksame-domain + * ParaBdiffKsameWorld from the old globals INT_BGROUP / BP_WORLD (MPI) or + * as a serial domain. + * + * Falls back to a serial single-band-group domain when MPI is not + * initialized or the pool layout has not been set up yet (e.g. unit + * tests), so that no MPI call is made on an unset communicator. + */ +ParaBdiffKsameWorld make_bdiff_ksame_world(); + } // namespace Parallel #endif // PARA_BRIDGE_H diff --git a/source/source_base/module_parallel/para_kmesh_world.cpp b/source/source_base/module_parallel/para_kmesh_world.cpp index 319654df6a3..8e58dc4f34f 100644 --- a/source/source_base/module_parallel/para_kmesh_world.cpp +++ b/source/source_base/module_parallel/para_kmesh_world.cpp @@ -7,28 +7,40 @@ namespace Parallel { ParaKmeshWorld::ParaKmeshWorld(int nkstot, int nspin) - : ParaWorld("kmesh"), kpar_(1), my_pool_(0), rank_in_pool_(0), - nproc_(1), nspin_(nspin), nkstot_(nkstot) + : ParaWorld("kmesh"), nspin_(nspin), nkstot_(nkstot) { distribute_kpoints(); nks_local_ = nkstot_; startk_global_ = 0; + nproc_ = size(); +} + +ParaKmeshWorld::ParaKmeshWorld() + : ParaWorld("kmesh"), nspin_(1) +{ + // Intentionally empty: no k-point distribution data. + // Only kpar_ / comm() are valid for reduce_across_pools. + nproc_ = size(); } #ifdef __MPI -ParaKmeshWorld::ParaKmeshWorld(const MPI_Comm& comm, int kpar, int my_pool, int nproc, int nkstot, int nspin) +ParaKmeshWorld::ParaKmeshWorld(const MPI_Comm& comm, int kpar, int my_pool, int nkstot, int nspin) : ParaWorld("kmesh", comm), kpar_(kpar), my_pool_(my_pool), - rank_in_pool_(rank()), nproc_(nproc), nspin_(nspin), nkstot_(nkstot) + rank_in_pool_(rank()), nspin_(nspin), nkstot_(nkstot) { + // nproc_ must be known before distribute_kpoints(), which derives the + // first rank of every k-pool from it. + nproc_ = size(); distribute_kpoints(); nks_local_ = nks_pool_[my_pool_]; startk_global_ = startk_pool_[my_pool_]; + kpool_root_ = (rank_in_pool_ == startpro_pool_[my_pool_]); } #endif void ParaKmeshWorld::distribute_kpoints() { - // k-points per pool (evenly divided, remainder to front) + // k-points per k-pool (evenly divided, remainder to front) nks_pool_.resize(kpar_, 0); const int nks_ave = nkstot_ / kpar_; const int nks_rem = nkstot_ % kpar_; @@ -37,14 +49,14 @@ void ParaKmeshWorld::distribute_kpoints() nks_pool_[i] = nks_ave + (i < nks_rem ? 1 : 0); } - // global start index per pool + // global start index per k-pool startk_pool_.resize(kpar_, 0); for (int i = 1; i < kpar_; ++i) { startk_pool_[i] = startk_pool_[i - 1] + nks_pool_[i - 1]; } - // pool index per k-point + // k-pool index per k-point whichpool_.resize(nkstot_, 0); for (int p = 0; p < kpar_; ++p) { @@ -54,7 +66,9 @@ void ParaKmeshWorld::distribute_kpoints() } } - // first world rank per pool + // first communicator rank per k-pool (processes are split into + // consecutive rank blocks, remainder to the front k-pools; this + // mirrors Parallel_Global::divide_mpi_groups) startpro_pool_.resize(kpar_, 0); const int nproc_ave = nproc_ / kpar_; const int nproc_rem = nproc_ % kpar_; @@ -90,9 +104,51 @@ int ParaKmeshWorld::startpro_pool(int pool) const int ParaKmeshWorld::max_nks_pool() const { + // Reduce-only domains carry no distribution arrays (see the default + // constructor); querying them here would dereference an empty vector. + assert(!nks_pool_.empty()); return *std::max_element(nks_pool_.begin(), nks_pool_.end()); } +void ParaKmeshWorld::reduce_across_pools(double& value) const +{ + if (kpar_ <= 1) + { + return; + } +#ifdef __MPI + // Exactly one contribution per k-pool: the first process of each + // k-pool injects the partial sum, all other processes inject zero. + // A single world-wide Allreduce therefore returns the sum of the + // per-k-pool partial sums, with no normalization division and with + // uneven k-pool sizes handled naturally. + const double local = kpool_root_ ? value : 0.0; + MPI_Allreduce(&local, &value, 1, MPI_DOUBLE, MPI_SUM, comm()); +#endif +} + +void ParaKmeshWorld::reduce_max_across_pools(double& value) const +{ + if (nproc_ <= 1) + { + return; + } +#ifdef __MPI + MPI_Allreduce(MPI_IN_PLACE, &value, 1, MPI_DOUBLE, MPI_MAX, comm()); +#endif +} + +void ParaKmeshWorld::reduce_min_across_pools(double& value) const +{ + if (nproc_ <= 1) + { + return; + } +#ifdef __MPI + MPI_Allreduce(MPI_IN_PLACE, &value, 1, MPI_DOUBLE, MPI_MIN, comm()); +#endif +} + void ParaKmeshWorld::pool_collection(double& value, const double* wk, int ik) const { #ifdef __MPI diff --git a/source/source_base/module_parallel/para_kmesh_world.h b/source/source_base/module_parallel/para_kmesh_world.h index f78607112d7..e65f7d0eba3 100644 --- a/source/source_base/module_parallel/para_kmesh_world.h +++ b/source/source_base/module_parallel/para_kmesh_world.h @@ -10,11 +10,17 @@ namespace Parallel { /** - * @brief k-mesh parallel domain: k-point distribution across pools. + * @brief k-mesh parallel domain: pure k-point pool (k-pool) topology. * * Self-contained replacement for Parallel_Kpoints + KP_WORLD + - * GlobalV::KPAR / MY_POOL / RANK_IN_POOL. Owns all k-point pool - * topology data and provides query / collection operations. + * GlobalV::KPAR / MY_POOL. Owns the k-point pool layout and provides + * query / collection operations. + * + * The k-pool split (kpar) is independent of bndpar: the domain knows + * nothing about band groups. Reductions that must span both dimensions + * (e.g. the total electron count under BPCG) therefore combine + * ParaBdiffKsameWorld::reduce_across_bdiff_ksame (band dimension, run FIRST) + * with reduce_across_pools (k dimension, run SECOND). * * In serial builds all operations degenerate to single-pool behavior. * Tests only need this header; no GlobalV, no parallel_comm.h. @@ -30,27 +36,38 @@ class ParaKmeshWorld : public ParaWorld */ ParaKmeshWorld(int nkstot, int nspin); + /** + * @brief Construct a reduce-only k-mesh domain with no k-point + * distribution data. + * + * kpar_/comm are set from the bridge globals so that + * reduce_across_pools / reduce_max/min_across_pools work correctly. + * The distribution arrays (nks_pool_, whichpool_, ...) are left + * empty; calling pool_collection / gather_kvec is invalid. + */ + ParaKmeshWorld(); + #ifdef __MPI /** * @brief Construct a k-mesh domain on an existing communicator. * - * @param[in] comm k-point pool communicator (e.g. KP_WORLD) - * @param[in] kpar number of pools - * @param[in] my_pool pool index of this process - * @param[in] nproc total number of processes (MPI_COMM_WORLD size) + * @param[in] comm communicator spanning every process of every + * k-pool (MPI_COMM_WORLD in the current bridge) + * @param[in] kpar number of k-pools + * @param[in] my_pool k-pool index of this process * @param[in] nkstot total number of k-points (without spin) * @param[in] nspin number of spin components */ - ParaKmeshWorld(const MPI_Comm& comm, int kpar, int my_pool, int nproc, int nkstot, int nspin); + ParaKmeshWorld(const MPI_Comm& comm, int kpar, int my_pool, int nkstot, int nspin); #endif - /// Number of pools. + /// Number of k-pools. int kpar() const { return kpar_; } - /// Pool index of this process. + /// k-pool index of this process. int my_pool() const { return my_pool_; } - /// Rank within the pool. + /// Rank within the communicator (world rank in the bridge layout). int rank_in_pool() const { return rank_in_pool_; } /// Total number of processes. @@ -62,27 +79,76 @@ class ParaKmeshWorld : public ParaWorld /// Total number of k-points (without spin). int nkstot() const { return nkstot_; } - /// Number of k-points in this pool. + /// Number of k-points in this k-pool. int nks_local() const { return nks_local_; } - /// Global start index of this pool's k-points. + /// Global start index of this k-pool's k-points. int startk_global() const { return startk_global_; } - /// Number of k-points in the given pool. + /// Number of k-points in the given k-pool. int nks_pool(int pool) const; - /// Global start index of the given pool's k-points. + /// Global start index of the given k-pool's k-points. int startk_pool(int pool) const; - /// Which pool owns the given global k-point index. + /// Which k-pool owns the given global k-point index. int which_pool(int ik_global) const; - /// First MPI_COMM_WORLD rank of the given pool. + /// First communicator rank of the given k-pool. int startpro_pool(int pool) const; - /// Maximum number of k-points across all pools. + /// Maximum number of k-points across all k-pools. int max_nks_pool() const; + /// Whether this process is the first process of its k-pool. + bool kpool_root() const { return kpool_root_; } + + // ===== Cross-pool reductions ===== + + /** + * @brief Sum a scalar across the k-pools: one contribution per pool. + * + * Replaces Parallel_Reduce::reduce_double_allpool. The first process + * of each k-pool injects the partial sum, all other processes inject + * zero, so a single world-wide MPI_Allreduce yields the sum of the + * per-k-pool partial sums. Correct for uneven k-pool sizes and free + * of the legacy normalization division (which divided by an average + * pool size and double-counted pools of uneven layouts). + * + * Precondition: with band parallelism (bndpar > 1) the caller must + * first combine the band-group partial sums (e.g. + * ParaBdiffKsameWorld::reduce_across_bdiff_ksame) so that every process of a + * k-pool holds one complete per-pool partial sum. + * + * No-op when kpar() <= 1. + * + * @param[in,out] value local partial sum, overwritten with global total + */ + void reduce_across_pools(double& value) const; + + /** + * @brief Global max across all k-pools and band groups. + * + * Max/min are idempotent, so a plain world-wide Allreduce is correct + * for every pool layout and covers both the k-pool and the band-group + * dimension (band-parallel shards see different eigenvalue windows, + * so the Fermi-level bounds must be extremized across both). + * Replaces Parallel_Reduce::reduce_max (all of MPI_COMM_WORLD). + * No-op when this domain spans a single process. + * + * @param[in,out] value local value, overwritten with global max + */ + void reduce_max_across_pools(double& value) const; + + /** + * @brief Global min across all k-pools and band groups. + * + * @param[in,out] value local value, overwritten with global min + */ + void reduce_min_across_pools(double& value) const; + + // ===== Cross-domain operations ===== + /** * @brief Collect a scalar value from the pool that owns k-point ik. * @@ -128,11 +194,12 @@ class ParaKmeshWorld : public ParaWorld int nkstot_ = 0; int nks_local_ = 0; int startk_global_ = 0; + bool kpool_root_ = true; ///< first process of my k-pool (reduction contributor) - std::vector nks_pool_; ///< k-points per pool - std::vector startk_pool_; ///< global start index per pool - std::vector whichpool_; ///< pool index per k-point - std::vector startpro_pool_; ///< first world rank per pool + std::vector nks_pool_; ///< k-points per k-pool + std::vector startk_pool_; ///< global start index per k-pool + std::vector whichpool_; ///< k-pool index per k-point + std::vector startpro_pool_; ///< first communicator rank per k-pool }; } // namespace Parallel diff --git a/source/source_base/module_parallel/para_tag.h b/source/source_base/module_parallel/para_tag.h index 7d2ff9543ce..2573cf2fd2e 100644 --- a/source/source_base/module_parallel/para_tag.h +++ b/source/source_base/module_parallel/para_tag.h @@ -9,14 +9,19 @@ namespace Parallel /** * @brief Domain tag constants for the parallel communication domains. * - * These tags replace raw string literals to avoid typo-induced runtime - * failures. They map to the legacy global communicators as follows: + * Pool terminology (see parallel_comm.cpp): + * - k-pool: one of the KPAR groups of processes that share one subset of + * k-points. This split happens first and is independent of bndpar. + * - band-pool: one of the BNDPAR sub-groups of a k-pool, holding one + * band window ("band group"). + * + * The tags map to the legacy global communicators as follows: * - esolver -> one esolver instance (intra-image communicator) * - images -> cross-image communicator (same rank_in_esolver) - * - pw -> POOL_WORLD - * - kmesh -> KP_WORLD - * - bsame_kdiff -> INT_BGROUP - * - bdiff_ksame -> BP_WORLD + * - pw -> POOL_WORLD (one band-pool) + * - kmesh -> KP_WORLD (links k-pools; only valid when the k-pool split is even) + * - bsame_kdiff -> INT_BGROUP (same band group across k-pools) + * - bdiff_ksame -> BP_WORLD (different band groups inside one k-pool) * - rgrid -> GRID_WORLD * - diag -> DIAG_WORLD * - matrix -> matrix domain diff --git a/source/source_base/module_parallel/para_world.cpp b/source/source_base/module_parallel/para_world.cpp index 0a4ca51748e..63138d0e322 100644 --- a/source/source_base/module_parallel/para_world.cpp +++ b/source/source_base/module_parallel/para_world.cpp @@ -3,22 +3,23 @@ namespace Parallel { -ParaWorld::ParaWorld(const std::string& tag) : tag_(tag), rank_(0), size_(1) +ParaWorld::ParaWorld(const std::string& tag) : tag_(tag), rank_(0), size_(1), comm_(nullptr) { #ifdef __MPI if (!tag.empty()) { - comm_ = MPI_COMM_SELF; + comm_ = handle_from_comm(MPI_COMM_SELF); } else { - comm_ = MPI_COMM_NULL; + comm_ = handle_from_comm(MPI_COMM_NULL); } #endif } #ifdef __MPI -ParaWorld::ParaWorld(const std::string& tag, const MPI_Comm& comm) : tag_(tag), comm_(comm) +ParaWorld::ParaWorld(const std::string& tag, const MPI_Comm& comm) + : tag_(tag), comm_(handle_from_comm(comm)) { if (comm == MPI_COMM_NULL) { @@ -34,7 +35,7 @@ ParaWorld::ParaWorld(const std::string& tag, const MPI_Comm& comm) : tag_(tag), bool ParaWorld::valid() const { #ifdef __MPI - return comm_ != MPI_COMM_NULL; + return comm() != MPI_COMM_NULL; #else return !tag_.empty(); #endif diff --git a/source/source_base/module_parallel/para_world.h b/source/source_base/module_parallel/para_world.h index a8291fad8bd..975af2a7f78 100644 --- a/source/source_base/module_parallel/para_world.h +++ b/source/source_base/module_parallel/para_world.h @@ -1,6 +1,7 @@ #ifndef PARA_WORLD_H #define PARA_WORLD_H +#include #include #include @@ -20,9 +21,13 @@ namespace Parallel * GlobalV::RANK_IN_POOL / POOL_WORLD by an object that functions * receive explicitly. * - * In serial builds (no __MPI) the communicator member does not - * exist; rank() always returns 0 and size() always returns 1, so - * call sites compile unchanged in both serial and MPI builds. + * The communicator is stored as an opaque handle so that the class + * layout is identical in serial and MPI builds. Binaries that mix + * translation units compiled with different __MPI settings (e.g. unit + * tests linked against the MPI-compiled base library) would otherwise + * be an ODR violation with undefined behavior. In serial builds rank() + * always returns 0 and size() always returns 1, so call sites compile + * unchanged in both serial and MPI builds. */ class ParaWorld { @@ -59,7 +64,11 @@ class ParaWorld /// Underlying MPI communicator (MPI builds only). MPI_Comm comm() const { - return comm_; + MPI_Comm comm = MPI_COMM_NULL; + static_assert(sizeof(MPI_Comm) <= sizeof(comm_), + "MPI_Comm does not fit into the opaque handle"); + std::memcpy(&comm, &comm_, sizeof(MPI_Comm)); + return comm; } #endif @@ -127,12 +136,23 @@ class ParaWorld #endif private: +#ifdef __MPI + /// Wrap an MPI communicator into the opaque handle storage. + static void* handle_from_comm(const MPI_Comm& comm) + { + void* handle = nullptr; + std::memcpy(&handle, &comm, sizeof(MPI_Comm)); + return handle; + } +#endif + std::string tag_; ///< domain tag int rank_; ///< rank inside domain int size_; ///< number of processes in domain -#ifdef __MPI - MPI_Comm comm_; ///< wrapped communicator (never owned/freed here) -#endif + // Opaque communicator handle, present in both serial and MPI builds + // so that the class layout never depends on the __MPI macro (see the + // class comment). Never owned/freed here. + void* comm_; }; } // namespace Parallel diff --git a/source/source_base/module_parallel/test/CMakeLists.txt b/source/source_base/module_parallel/test/CMakeLists.txt index 6f3a2169a7c..f0dc5174190 100644 --- a/source/source_base/module_parallel/test/CMakeLists.txt +++ b/source/source_base/module_parallel/test/CMakeLists.txt @@ -30,8 +30,8 @@ AddTest( ) AddTest( - TARGET MODULE_BASE_para_bgroup_world - SOURCES para_bgroup_world_test.cpp ../para_bgroup_world.cpp ../para_world.cpp + TARGET MODULE_BASE_para_bdiff_ksame_world + SOURCES para_bdiff_ksame_world_test.cpp ../para_bdiff_ksame_world.cpp ../para_world.cpp ) AddTest( @@ -77,10 +77,18 @@ AddTest( ) target_compile_definitions(MODULE_BASE_para_setup_mpi PRIVATE __MPI) +# Built with add_executable (not AddTest) so that no direct-run CTest entry is +# created; the binary is only exercised through mpirun by the .sh test below, +# matching the multi-process requirement of these cases. +add_executable(MODULE_BASE_para_kmesh_world_mpi test_para_kmesh_world_mpi.cpp ../para_bdiff_ksame_world.cpp ../para_kmesh_world.cpp ../para_world.cpp) +target_link_libraries(MODULE_BASE_para_kmesh_world_mpi PRIVATE MPI::MPI_CXX GTest::gtest GTest::gtest_main abacus::linalg_libs) +target_compile_definitions(MODULE_BASE_para_kmesh_world_mpi PRIVATE __MPI) + file(COPY para_world_mpi_test.sh DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) file(COPY para_collection_mpi_test.sh DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) file(COPY para_mpi_func_mpi_test.sh DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) file(COPY para_setup_mpi_test.sh DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) +file(COPY test_para_kmesh_world_mpi.sh DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) find_program(BASH bash) add_test(NAME MODULE_BASE_para_world_mpi_test COMMAND ${BASH} para_world_mpi_test.sh @@ -98,3 +106,7 @@ add_test(NAME MODULE_BASE_para_setup_mpi_test COMMAND ${BASH} para_setup_mpi_test.sh WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} ) +add_test(NAME MODULE_BASE_para_kmesh_world_mpi_test + COMMAND ${BASH} test_para_kmesh_world_mpi.sh + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} +) diff --git a/source/source_base/module_parallel/test/para_bgroup_world_test.cpp b/source/source_base/module_parallel/test/para_bdiff_ksame_world_test.cpp similarity index 63% rename from source/source_base/module_parallel/test/para_bgroup_world_test.cpp rename to source/source_base/module_parallel/test/para_bdiff_ksame_world_test.cpp index 8dceccf4ad9..4fde3d5c6c8 100644 --- a/source/source_base/module_parallel/test/para_bgroup_world_test.cpp +++ b/source/source_base/module_parallel/test/para_bdiff_ksame_world_test.cpp @@ -1,10 +1,10 @@ #include "gtest/gtest.h" -#include "../para_bgroup_world.h" +#include "../para_bdiff_ksame_world.h" -TEST(ParaBgroupWorldTest, SerialMode) +TEST(ParaBdiffKsameWorldTest, SerialMode) { - const Parallel::ParaBgroupWorld world; + const Parallel::ParaBdiffKsameWorld world; EXPECT_EQ(world.tag(), "bdiff_ksame"); EXPECT_EQ(world.my_bndgroup(), 0); EXPECT_EQ(world.nbndgroup(), 1); @@ -13,9 +13,9 @@ TEST(ParaBgroupWorldTest, SerialMode) EXPECT_TRUE(world.valid()); } -TEST(ParaBgroupWorldTest, AliasesMatchBase) +TEST(ParaBdiffKsameWorldTest, AliasesMatchBase) { - const Parallel::ParaBgroupWorld world; + const Parallel::ParaBdiffKsameWorld world; EXPECT_EQ(world.rank_in_bpgroup(), world.rank()); EXPECT_EQ(world.nproc_in_bndgroup(), world.size()); } diff --git a/source/source_base/module_parallel/test/para_collection_mpi_test.cpp b/source/source_base/module_parallel/test/para_collection_mpi_test.cpp index a166ff05322..6c06c4c4dd6 100644 --- a/source/source_base/module_parallel/test/para_collection_mpi_test.cpp +++ b/source/source_base/module_parallel/test/para_collection_mpi_test.cpp @@ -8,7 +8,7 @@ TEST(ParaCollectionMpiTest, AssembleAndFind) { Parallel::ParaCollection coll; coll.add(std::unique_ptr( - new Parallel::ParaKmeshWorld(MPI_COMM_WORLD, 1, 0, 1, 4, 1))); + new Parallel::ParaKmeshWorld(MPI_COMM_WORLD, 1, 0, 4, 1))); coll.add(Parallel::ParaWorld::make_serial(Parallel::ParaTag::pw)); EXPECT_EQ(coll.size(), 2u); @@ -25,7 +25,7 @@ TEST(ParaCollectionMpiTest, FindMissingReturnsInvalid) { Parallel::ParaCollection coll; coll.add(std::unique_ptr( - new Parallel::ParaKmeshWorld(MPI_COMM_WORLD, 1, 0, 1, 4, 1))); + new Parallel::ParaKmeshWorld(MPI_COMM_WORLD, 1, 0, 4, 1))); const Parallel::ParaWorld& missing = coll.find("nonexistent"); EXPECT_FALSE(missing.valid()); @@ -35,7 +35,7 @@ TEST(ParaCollectionMpiTest, FindAsSubclass) { Parallel::ParaCollection coll; coll.add(std::unique_ptr( - new Parallel::ParaKmeshWorld(MPI_COMM_WORLD, 1, 0, 1, 8, 1))); + new Parallel::ParaKmeshWorld(MPI_COMM_WORLD, 1, 0, 8, 1))); const Parallel::ParaKmeshWorld* kmesh = coll.find_as(Parallel::ParaTag::kmesh); ASSERT_NE(kmesh, nullptr); diff --git a/source/source_base/module_parallel/test/test_para_kmesh_world_mpi.cpp b/source/source_base/module_parallel/test/test_para_kmesh_world_mpi.cpp new file mode 100644 index 00000000000..20258140e03 --- /dev/null +++ b/source/source_base/module_parallel/test/test_para_kmesh_world_mpi.cpp @@ -0,0 +1,131 @@ +#include "gtest/gtest.h" + +#include "../para_bdiff_ksame_world.h" +#include "../para_kmesh_world.h" + +// Run with: mpirun -np 4 ./MODULE_BASE_para_kmesh_world_mpi +// +// The sum reduction protocol has two layers: +// 1. ParaBdiffKsameWorld::reduce_across_bdiff_ksame (band dimension, BPCG shards) +// 2. ParaKmeshWorld::reduce_across_pools (k dimension, one +// contribution per k-pool: the first rank of each k-pool injects the +// partial sum, everyone else injects zero) +// The band layer must run first so that the k layer receives one complete +// per-k-pool partial sum. + +TEST(ParaKmeshWorldMpiTest, ReduceAcrossBandGroupsBndpar2) +{ + int nprocs = 0; + int myrank = 0; + MPI_Comm_size(MPI_COMM_WORLD, &nprocs); + MPI_Comm_rank(MPI_COMM_WORLD, &myrank); + ASSERT_EQ(nprocs, 4); + + // kpar=1, bndpar=2, 4 ranks: band group = myrank/2, rank position + // inside the band group = myrank%2. Reproduce the BP_WORLD layout, + // which links the same rank position of every band group. + MPI_Comm bp_world = MPI_COMM_NULL; + MPI_Comm_split(MPI_COMM_WORLD, myrank % 2, myrank / 2, &bp_world); + Parallel::ParaBdiffKsameWorld bdiff(MPI_COMM_WORLD, bp_world, 2); + + // Each band group holds a partial occupation sum of 14 (28 electrons + // split into two band windows). + double sumk = 14.0; + bdiff.reduce_across_bdiff_ksame(sumk); + EXPECT_DOUBLE_EQ(sumk, 28.0); + + // max/min stay world-wide (idempotent) and must span the band groups + // as well: the two shards see different eigenvalue windows. + Parallel::ParaKmeshWorld kmesh(MPI_COMM_WORLD, 1, 0, 0, 1); + double eup = (myrank < 2) ? 40.0 : 45.0; + kmesh.reduce_max_across_pools(eup); + EXPECT_DOUBLE_EQ(eup, 45.0); + + double elw = (myrank < 2) ? -1.0 : -5.0; + kmesh.reduce_min_across_pools(elw); + EXPECT_DOUBLE_EQ(elw, -5.0); + + MPI_Comm_free(&bp_world); +} + +TEST(ParaKmeshWorldMpiTest, ReduceAcrossKpoolsKpar2) +{ + int nprocs = 0; + int myrank = 0; + MPI_Comm_size(MPI_COMM_WORLD, &nprocs); + MPI_Comm_rank(MPI_COMM_WORLD, &myrank); + ASSERT_EQ(nprocs, 4); + + // kpar=2, bndpar=1: k-pool 0 = ranks {0,1}, k-pool 1 = ranks {2,3} + // (consecutive rank blocks, divide_mpi_groups layout). + const int my_pool = myrank / 2; + Parallel::ParaKmeshWorld kmesh(MPI_COMM_WORLD, 2, my_pool, 4, 1); + EXPECT_EQ(kmesh.startpro_pool(0), 0); + EXPECT_EQ(kmesh.startpro_pool(1), 2); + EXPECT_EQ(kmesh.kpool_root(), (myrank % 2 == 0)); + + // Every process holds its k-pool's partial sum; the reduction must + // count each pool exactly once. + double sumk = 3.5; + kmesh.reduce_across_pools(sumk); + EXPECT_DOUBLE_EQ(sumk, 7.0); +} + +TEST(ParaKmeshWorldMpiTest, UnevenKpoolsKpar3) +{ + int nprocs = 0; + int myrank = 0; + MPI_Comm_size(MPI_COMM_WORLD, &nprocs); + MPI_Comm_rank(MPI_COMM_WORLD, &myrank); + ASSERT_EQ(nprocs, 4); + + // nproc=4, kpar=3 (the 007_PW_UPF201_USPP_Fe layout): k-pool sizes + // are [2,1,1]. divide_mpi_groups puts ranks {0,1} in pool 0, rank 2 + // in pool 1 and rank 3 in pool 2. + int my_pool = 0; + if (myrank >= 3) + { + my_pool = 2; + } + else if (myrank >= 2) + { + my_pool = 1; + } + Parallel::ParaKmeshWorld kmesh(MPI_COMM_WORLD, 3, my_pool, 0, 1); + EXPECT_EQ(kmesh.startpro_pool(0), 0); + EXPECT_EQ(kmesh.startpro_pool(1), 2); + EXPECT_EQ(kmesh.startpro_pool(2), 3); + EXPECT_EQ(kmesh.kpool_root(), (myrank != 1)); + + // Every process holds its k-pool's partial sum. The reduction must + // count each pool exactly once even though the pools are uneven: + // the legacy average-pool-size division (4/3 = 1) double-counted + // pool 0 here and corrupted the electron count / Fermi level. + const double pool_sum = (my_pool == 0) ? 10.0 : ((my_pool == 1) ? 20.0 : 30.0); + double sumk = pool_sum; + kmesh.reduce_across_pools(sumk); + EXPECT_DOUBLE_EQ(sumk, 60.0); +} + +TEST(ParaKmeshWorldMpiTest, SinglePoolIsNoOp) +{ + int nprocs = 0; + MPI_Comm_size(MPI_COMM_WORLD, &nprocs); + ASSERT_EQ(nprocs, 4); + + // kpar == 1: the sum reduction must be a no-op regardless of the + // world size (the band dimension is handled by ParaBdiffKsameWorld). + Parallel::ParaKmeshWorld kmesh(MPI_COMM_WORLD, 1, 0, 0, 1); + double sumk = 42.0; + kmesh.reduce_across_pools(sumk); + EXPECT_DOUBLE_EQ(sumk, 42.0); +} + +int main(int argc, char** argv) +{ + MPI_Init(&argc, &argv); + testing::InitGoogleTest(&argc, argv); + const int result = RUN_ALL_TESTS(); + MPI_Finalize(); + return result; +} diff --git a/source/source_base/module_parallel/test/test_para_kmesh_world_mpi.sh b/source/source_base/module_parallel/test/test_para_kmesh_world_mpi.sh new file mode 100644 index 00000000000..dd934e2ad4a --- /dev/null +++ b/source/source_base/module_parallel/test/test_para_kmesh_world_mpi.sh @@ -0,0 +1,18 @@ +#!/bin/bash -e + +np=`cat /proc/cpuinfo | grep "cpu cores" | uniq| awk '{print $NF}'` +echo "nprocs in this machine is $np" + +for i in 4;do + if [[ $i -gt $np ]];then + continue + fi + echo "TEST in parallel, nprocs=$i" + mpirun -np $i ./MODULE_BASE_para_kmesh_world_mpi + if [[ $? -ne 0 ]]; then + echo -e "\e[1;33m [ FAILED ] \e[0m"\ + "execute UT with $i cores error." + exit 1 + fi + break +done diff --git a/source/source_base/parallel_comm.cpp b/source/source_base/parallel_comm.cpp index 5d03447b5aa..27eac402b1a 100644 --- a/source/source_base/parallel_comm.cpp +++ b/source/source_base/parallel_comm.cpp @@ -3,10 +3,18 @@ #include "mpi.h" #include "parallel_global.h" -MPI_Comm POOL_WORLD; //groups for different plane waves. In this group, only plane waves are different. K-points and bands are the same. -MPI_Comm KP_WORLD; // groups for differnt k. In this group, only k-points are different. Bands and plane waves are the same. -MPI_Comm BP_WORLD; // groups for differnt bands. In this group, only bands are different. K-points and plane waves are the same. -MPI_Comm INT_BGROUP; // internal comm groups for same bands. In this group, only bands are the same. K-points and plane waves are different. +// Two-level pool terminology used across the parallel layer: +// - k-pool: a group of processes that share one subset of k-points. The +// processes are split into KPAR k-pools first (divide_pools); this split +// is independent of bndpar. MY_POOL / KP_WORLD refer to this level. +// - band-pool: a sub-group of one k-pool, created afterwards by dividing +// the k-pool into BNDPAR band groups. NPROC_IN_POOL / RANK_IN_POOL / +// POOL_WORLD refer to this level, i.e. the term "pool" in those globals +// means the (k-pool, band-group) cell, NOT the k-pool itself. +MPI_Comm POOL_WORLD; // one band-pool (k-pool x band-group cell): plane waves are distributed, k-points and the band window are shared. +MPI_Comm KP_WORLD; // links k-pools: only k-points differ; same rank_in_pool position in every k-pool. Valid ONLY when k-pools are equally sized (NPROC % KPAR == 0), otherwise MPI_COMM_NULL. +MPI_Comm BP_WORLD; // links band groups inside one k-pool: only the band window differs; k-points and plane-wave slab are the same. One communicator per rank position. +MPI_Comm INT_BGROUP; // same band-group index across all k-pools (plus the plane-wave ranks of that band group): k-points differ, the band window is the same. Always valid, also for uneven k-pools. MPI_Comm GRID_WORLD; // mohan add 2012-01-13 MPI_Comm DIAG_WORLD; // mohan add 2012-01-13 diff --git a/source/source_base/parallel_comm.h b/source/source_base/parallel_comm.h index 2243aea729f..14222728630 100644 --- a/source/source_base/parallel_comm.h +++ b/source/source_base/parallel_comm.h @@ -3,10 +3,10 @@ #ifdef __MPI #include "mpi.h" -extern MPI_Comm POOL_WORLD; -extern MPI_Comm KP_WORLD; // communicator among different pools -extern MPI_Comm INT_BGROUP; -extern MPI_Comm BP_WORLD; +extern MPI_Comm POOL_WORLD; // one band-pool (k-pool x band-group cell): only plane waves are distributed +extern MPI_Comm KP_WORLD; // links k-pools at the same rank_in_pool position; MPI_COMM_NULL when k-pools are uneven +extern MPI_Comm INT_BGROUP; // same band-group index across all k-pools +extern MPI_Comm BP_WORLD; // links band groups inside one k-pool (same k, different band windows) extern MPI_Comm GRID_WORLD; // mohan add 2012-01-13 extern MPI_Comm DIAG_WORLD; // mohan add 2012-01-13 diff --git a/source/source_base/parallel_global.cpp b/source/source_base/parallel_global.cpp index 697b7f1f702..fa7a5aa40ae 100644 --- a/source/source_base/parallel_global.cpp +++ b/source/source_base/parallel_global.cpp @@ -230,9 +230,14 @@ void Parallel_Global::divide_pools(const int& NPROC, int& RANK_IN_POOL, int& MY_POOL) { - // note: the order of k-point parallelization and band parallelization is important - // The order will not change the behavior of KP_WORLD or BP_WORLD, and MY_POOL - // and MY_BNDGROUP will be the same as well. + // Two-level split, order matters: + // 1. k-point parallelization: NPROC processes are divided into KPAR + // k-pools FIRST, independent of BNDPAR. MY_POOL is the k-pool index. + // Uneven k-pool sizes (NPROC % KPAR != 0) are allowed here; in that + // case KP_WORLD is MPI_COMM_NULL (see MPICommGroup::divide_group_comm). + // 2. band parallelization: each k-pool is divided into BNDPAR band + // groups ("band-pools"). NPROC_IN_POOL / RANK_IN_POOL / POOL_WORLD + // belong to this (k-pool x band-group) cell, NOT to the k-pool. if(BNDPAR > 1 && NPROC %(BNDPAR * KPAR) != 0) { std::cout << "Error: When BNDPAR = " << BNDPAR << " > 1, number of processes (" << NPROC diff --git a/source/source_base/test_parallel/parallel_device_test.cpp b/source/source_base/test_parallel/parallel_device_test.cpp index 0c6746267eb..ba4a4d90b86 100644 --- a/source/source_base/test_parallel/parallel_device_test.cpp +++ b/source/source_base/test_parallel/parallel_device_test.cpp @@ -124,7 +124,7 @@ TEST(ParallelDevice, CoversGpuStagingWithoutAccelerator) TEST(ParallelDevice, CoversMpiTypeOverloads) { - const ModuleBase::CommunicationDomain domain = ModuleBase::world_communication_domain(); + const ModuleBase::CommunicationDomain domain = ModuleBase::world_comm_domain(); exercise_mpi_wrappers(domain); exercise_mpi_wrappers(domain); exercise_mpi_wrappers>(domain); diff --git a/source/source_base/test_parallel/parallel_domain_grid_test.cpp b/source/source_base/test_parallel/parallel_domain_grid_test.cpp index f71581e74e2..5a066c6e0ab 100644 --- a/source/source_base/test_parallel/parallel_domain_grid_test.cpp +++ b/source/source_base/test_parallel/parallel_domain_grid_test.cpp @@ -14,7 +14,7 @@ TEST(CommunicationDomainTest, ReportsDefaultAndWorldDomains) EXPECT_EQ(local_domain.rank(), 0); EXPECT_EQ(local_domain.communicator(), MPI_COMM_NULL); - const ModuleBase::CommunicationDomain world_domain = ModuleBase::world_communication_domain(); + const ModuleBase::CommunicationDomain world_domain = ModuleBase::world_comm_domain(); MPICommGroup world_group(world_domain.communicator()); EXPECT_EQ(world_domain.communicator(), MPI_COMM_WORLD); EXPECT_GE(world_domain.rank(), 0); @@ -28,7 +28,7 @@ TEST(CommunicationDomainTest, ReportsDefaultAndWorldDomains) TEST(MPICommGroupTest, DividesWorldIntoEvenGroups) { - const ModuleBase::CommunicationDomain world_domain = ModuleBase::world_communication_domain(); + const ModuleBase::CommunicationDomain world_domain = ModuleBase::world_comm_domain(); MPICommGroup group(MPI_COMM_WORLD); EXPECT_EQ(group.grank, world_domain.rank()); @@ -46,7 +46,7 @@ TEST(MPICommGroupTest, DividesWorldIntoEvenGroups) TEST(ParallelGridTest, BroadcastsAndReducesDistributedGrid) { - const ModuleBase::CommunicationDomain world_domain = ModuleBase::world_communication_domain(); + const ModuleBase::CommunicationDomain world_domain = ModuleBase::world_comm_domain(); MPICommGroup world_group(world_domain.communicator()); const int nx = 2; const int ny = 1; diff --git a/source/source_base/test_parallel/test_para_gemm.cpp b/source/source_base/test_parallel/test_para_gemm.cpp index f116507e7f6..0e62823c3be 100644 --- a/source/source_base/test_parallel/test_para_gemm.cpp +++ b/source/source_base/test_parallel/test_para_gemm.cpp @@ -74,7 +74,7 @@ void expect_near_value(const T& actual, const T& expected) template void test_additional_type_paths() { - const ModuleBase::CommunicationDomain domain = ModuleBase::world_communication_domain(); + const ModuleBase::CommunicationDomain domain = ModuleBase::world_comm_domain(); MPI_Comm world = domain.communicator(); const int rank = domain.rank(); MPICommGroup world_group(world); diff --git a/source/source_cell/CMakeLists.txt b/source/source_cell/CMakeLists.txt index 46eacaec8f0..06f70560853 100644 --- a/source/source_cell/CMakeLists.txt +++ b/source/source_cell/CMakeLists.txt @@ -44,8 +44,16 @@ add_library( cal_ux.cpp cif_io.cpp ucell_io.cpp + read_cube.cpp + write_cube.cpp + write_pao.cpp + output_log.cpp ) +if(ENABLE_LCAO) + target_sources(cell PRIVATE write_orb_info.cpp) +endif() + if(ENABLE_COVERAGE) add_coverage(cell) endif() diff --git a/source/source_io/module_output/cube_io.h b/source/source_cell/cube_io.h similarity index 100% rename from source/source_io/module_output/cube_io.h rename to source/source_cell/cube_io.h diff --git a/source/source_io/module_output/output_log.cpp b/source/source_cell/output_log.cpp similarity index 100% rename from source/source_io/module_output/output_log.cpp rename to source/source_cell/output_log.cpp diff --git a/source/source_io/module_output/output_log.h b/source/source_cell/output_log.h similarity index 100% rename from source/source_io/module_output/output_log.h rename to source/source_cell/output_log.h diff --git a/source/source_io/module_output/read_cube.cpp b/source/source_cell/read_cube.cpp similarity index 99% rename from source/source_io/module_output/read_cube.cpp rename to source/source_cell/read_cube.cpp index a4155f78f89..cbae974ed04 100644 --- a/source/source_io/module_output/read_cube.cpp +++ b/source/source_cell/read_cube.cpp @@ -1,4 +1,4 @@ -#include "source_io/module_output/cube_io.h" +#include "source_cell/cube_io.h" #include #include "source_base/parallel_grid.h" #include "source_io/module_parameter/parameter.h" diff --git a/source/source_io/module_output/write_cube.cpp b/source/source_cell/write_cube.cpp similarity index 99% rename from source/source_io/module_output/write_cube.cpp rename to source/source_cell/write_cube.cpp index 1771a28026b..e328338f8f0 100644 --- a/source/source_io/module_output/write_cube.cpp +++ b/source/source_cell/write_cube.cpp @@ -1,7 +1,7 @@ #include "source_base/element_name.h" #include "source_base/parallel_comm.h" #include "source_base/parallel_grid.h" -#include "source_io/module_output/cube_io.h" +#include "source_cell/cube_io.h" #include diff --git a/source/source_io/module_output/write_orb_info.cpp b/source/source_cell/write_orb_info.cpp similarity index 100% rename from source/source_io/module_output/write_orb_info.cpp rename to source/source_cell/write_orb_info.cpp diff --git a/source/source_io/module_output/write_orb_info.h b/source/source_cell/write_orb_info.h similarity index 100% rename from source/source_io/module_output/write_orb_info.h rename to source/source_cell/write_orb_info.h diff --git a/source/source_io/module_output/write_pao.cpp b/source/source_cell/write_pao.cpp similarity index 100% rename from source/source_io/module_output/write_pao.cpp rename to source/source_cell/write_pao.cpp diff --git a/source/source_io/module_output/write_pao.h b/source/source_cell/write_pao.h similarity index 100% rename from source/source_io/module_output/write_pao.h rename to source/source_cell/write_pao.h diff --git a/source/source_esolver/esolver_dm2rho.cpp b/source/source_esolver/esolver_dm2rho.cpp index 4f1c042deb7..8a84a37b81e 100644 --- a/source/source_esolver/esolver_dm2rho.cpp +++ b/source/source_esolver/esolver_dm2rho.cpp @@ -5,7 +5,7 @@ #include "source_cell/read_pp_ucell.h" #include "source_estate/elecstate_lcao.h" #include "source_io/module_ml/io_npz.h" -#include "source_io/module_output/cube_io.h" +#include "source_cell/cube_io.h" #include "source_lcao/lcao_domain.h" #include "source_lcao/hamilt_lcao.h" #include "source_lcao/module_operator_lcao/operator_lcao.h" diff --git a/source/source_esolver/esolver_dp.cpp b/source/source_esolver/esolver_dp.cpp index 0a641ef462e..127d8f7cbee 100644 --- a/source/source_esolver/esolver_dp.cpp +++ b/source/source_esolver/esolver_dp.cpp @@ -23,7 +23,7 @@ #include "source_cell/mdcell.h" #include "source_cell/module_neighlist/neighbor_search.h" #include "source_cell/cif_io.h" -#include "source_io/module_output/output_log.h" +#include "source_cell/output_log.h" #include "source_io/module_parameter/parameter.h" #include diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp index 0d824a84244..05eb897a8dc 100644 --- a/source/source_esolver/esolver_fp.cpp +++ b/source/source_esolver/esolver_fp.cpp @@ -7,7 +7,7 @@ #include "source_estate/param_update.h" #include "source_hamilt/module_ewald/h_ewald_pw.h" #include "source_hamilt/module_vdw/vdw.h" -#include "source_io/module_output/output_log.h" +#include "source_cell/output_log.h" #include "source_io/module_output/print_info.h" #include "source_estate/rhog_io.h" #include "source_io/module_parameter/parameter.h" @@ -153,8 +153,8 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso CE.update_delta_rho(ucell, &(this->chr), &(this->sf)); //! print out charge density, potential, elf, etc. - ModuleIO::ctrl_output_fp(ucell, this->pelec, this->pw_big, this->pw_rhod, - this->chr, this->solvent, this->Pgrid, istep); + ModuleIO::ctrl_output_fp(ucell, this->pelec, this->pw_big, this->pw_rhod, + this->chr, this->solvent, this->Pgrid, istep, PARAM.inp); } diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp index dec3fdd5770..d0264e70e5a 100644 --- a/source/source_esolver/esolver_ks.cpp +++ b/source/source_esolver/esolver_ks.cpp @@ -11,7 +11,7 @@ #include "source_io/module_energy/write_eig_occ.h" #include "source_io/module_energy/write_bands.h" #include "source_hamilt/module_xc/xc_functional.h" -#include "source_io/module_output/output_log.h" // use write_head +#include "source_cell/output_log.h" // use write_head #include "source_estate/elecstate_print.h" // print_etot #include "source_lcao/module_dftu/dftu_nao.h" // mohan add 2025-11-07 #include "source_hamilt/module_xc/general_exx_info.h" // for init_general_exx_info diff --git a/source/source_esolver/esolver_ks_lcao_tddft.cpp b/source/source_esolver/esolver_ks_lcao_tddft.cpp index 0b3c63ce16c..67a9f12aa9b 100644 --- a/source/source_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/source_esolver/esolver_ks_lcao_tddft.cpp @@ -8,7 +8,7 @@ #include "source_io/module_ctrl/ctrl_output_td.h" #include "source_io/module_efield/td_efield_io.h" #include "source_io/module_efield/td_vector_pot_io.h" -#include "source_io/module_output/output_log.h" +#include "source_cell/output_log.h" #include "source_io/module_parameter/parameter.h" #include "source_io/module_wf/read_wfc_nao.h" //------LCAO HSolver ElecState------- diff --git a/source/source_esolver/esolver_lj.cpp b/source/source_esolver/esolver_lj.cpp index 4daee953a5f..5b085d30302 100644 --- a/source/source_esolver/esolver_lj.cpp +++ b/source/source_esolver/esolver_lj.cpp @@ -6,7 +6,7 @@ #include "source_cell/module_neighlist/neighbor_types.h" #include "source_io/module_parameter/parameter.h" #include "source_cell/cif_io.h" -#include "source_io/module_output/output_log.h" +#include "source_cell/output_log.h" #ifdef __MPI #include #endif diff --git a/source/source_esolver/esolver_lr_lcao_tddft.cpp b/source/source_esolver/esolver_lr_lcao_tddft.cpp index d53902d3131..688eeac89de 100644 --- a/source/source_esolver/esolver_lr_lcao_tddft.cpp +++ b/source/source_esolver/esolver_lr_lcao_tddft.cpp @@ -11,7 +11,7 @@ #include #include "source_lcao/hamilt_lcao.h" #include "source_io/module_wf/read_wfc_nao.h" -#include "source_io/module_output/cube_io.h" +#include "source_cell/cube_io.h" #include "source_io/module_output/print_info.h" #include "source_cell/module_neighbor/sltk_atom_arrange.h" #include "source_lcao/module_lr/utils/lr_util_print.h" diff --git a/source/source_esolver/esolver_nep.cpp b/source/source_esolver/esolver_nep.cpp index e9728020035..6622f3c592d 100644 --- a/source/source_esolver/esolver_nep.cpp +++ b/source/source_esolver/esolver_nep.cpp @@ -21,7 +21,7 @@ #include "source_cell/mdcell.h" #include "source_cell/module_neighlist/neighbor_search.h" #include "source_cell/cif_io.h" -#include "source_io/module_output/output_log.h" +#include "source_cell/output_log.h" #include "source_io/module_parameter/parameter.h" #include diff --git a/source/source_esolver/test/CMakeLists.txt b/source/source_esolver/test/CMakeLists.txt index 80ca8f061b1..38c823ea67c 100644 --- a/source/source_esolver/test/CMakeLists.txt +++ b/source/source_esolver/test/CMakeLists.txt @@ -24,5 +24,5 @@ AddTest( ../esolver_dp.cpp ../../source_cell/basecell.cpp ../../source_cell/cif_io.cpp - ../../source_io/module_output/output_log.cpp + ../../source_cell/output_log.cpp ) diff --git a/source/source_estate/elecstate_tools.cpp b/source/source_estate/elecstate_tools.cpp index e4e6c2930a5..3f9147d29e2 100644 --- a/source/source_estate/elecstate_tools.cpp +++ b/source/source_estate/elecstate_tools.cpp @@ -1,8 +1,8 @@ #include "elecstate_tools.h" #include "occupy.h" +#include "source_base/module_parallel/para_bridge.h" #include "source_base/parallel_comm.h" -#include "source_base/parallel_reduce.h" #include #include @@ -74,11 +74,16 @@ void calEBand(const ModuleBase::matrix& ekb, const ModuleBase::matrix& wg, fener } f_en.eband = eband; -#ifdef __MPI - // Combine contributions distributed by both KPAR and BNDPAR. - const int npool = GlobalV::KPAR * PARAM.inp.bndpar; - Parallel_Reduce::reduce_double_allpool(npool, GlobalV::NPROC_IN_POOL, f_en.eband); -#endif + // Two-step reduction, order matters: + // 1. Band dimension: BPCG shards the band range across the band + // groups of this k-pool, so combine the per-window partial sums + // first (no-op with a single band group). + // 2. k dimension: exactly one contribution per k-pool. + // Reduce-only domains: no k-point distribution data needed. + Parallel::ParaKmeshWorld kmesh = Parallel::make_kmesh_world(); + Parallel::ParaBdiffKsameWorld bdiff = Parallel::make_bdiff_ksame_world(); + bdiff.reduce_across_bdiff_ksame(f_en.eband); + kmesh.reduce_across_pools(f_en.eband); return; } @@ -98,27 +103,41 @@ void calculate_weights(const ModuleBase::matrix& ekb, const int nbands = ekb.nc; const int nks = ekb.nr; + const int nspin = PARAM.inp.nspin; if (!(Occupy::use_gaussian_broadening || Occupy::fixed_occupations)) { // Taoni fix smearing_method=fixed for BPCG on 2026-08-21 // Integer occupations use global band indices even when ekb is a local // contiguous BPCG shard. const int band_offset = get_band_offset(nbands, global_nbands); + // The kmesh domain is only built in the branches that dereference + // klist, so that callers passing no k-list (fixed occupations) are + // not affected. Only the reduction is needed here, so use the + // reduce-only overload and skip building the distribution arrays + // on every SCF iteration. + Parallel::ParaKmeshWorld kmesh = Parallel::make_kmesh_world(); if (PARAM.globalv.two_fermi) { - Occupy::iweights(nks, klist->wk, nbands, band_offset, nelec_spin[0], ekb, eferm.ef_up, wg, 0, klist->isk); - Occupy::iweights(nks, klist->wk, nbands, band_offset, nelec_spin[1], ekb, eferm.ef_dw, wg, 1, klist->isk); + Occupy::iweights(nks, klist->wk, nbands, band_offset, nelec_spin[0], ekb, eferm.ef_up, wg, + nspin, 0, klist->isk, kmesh); + Occupy::iweights(nks, klist->wk, nbands, band_offset, nelec_spin[1], ekb, eferm.ef_dw, wg, + nspin, 1, klist->isk, kmesh); // ef = ( ef_up + ef_dw ) / 2.0_dp need??? mohan add 2012-04-16 // Keep independent Fermi levels for the two spin channels. } else { // A spin selector of -1 requests the combined-spin occupation path. - Occupy::iweights(nks, klist->wk, nbands, band_offset, PARAM.inp.nelec, ekb, eferm.ef, wg, -1, klist->isk); + Occupy::iweights(nks, klist->wk, nbands, band_offset, PARAM.inp.nelec, ekb, eferm.ef, wg, + nspin, -1, klist->isk, kmesh); } } else if (Occupy::use_gaussian_broadening) { + // The kmesh domain is only built in the branches that dereference + // klist, so that callers passing no k-list (fixed occupations) are + // not affected. Reduce-only overload: see the iweights branch above. + Parallel::ParaKmeshWorld kmesh = Parallel::make_kmesh_world(); if (PARAM.globalv.two_fermi) { double demet_up = 0.0; @@ -134,7 +153,8 @@ void calculate_weights(const ModuleBase::matrix& ekb, demet_up, wg, 0, - klist->isk); + klist->isk, + kmesh); Occupy::gweights(nks, klist->wk, nbands, @@ -146,7 +166,8 @@ void calculate_weights(const ModuleBase::matrix& ekb, demet_dw, wg, 1, - klist->isk); + klist->isk, + kmesh); f_en.demet = demet_up + demet_dw; } else @@ -163,13 +184,15 @@ void calculate_weights(const ModuleBase::matrix& ekb, f_en.demet, wg, -1, - klist->isk); + klist->isk, + kmesh); } -#ifdef __MPI // demet is accumulated independently on every k-point and band partition. - const int npool = GlobalV::KPAR * PARAM.inp.bndpar; - Parallel_Reduce::reduce_double_allpool(npool, GlobalV::NPROC_IN_POOL, f_en.demet); -#endif + // Band dimension first (BPCG band shards), then one contribution + // per k-pool; see calEBand for the ordering rationale. + Parallel::ParaBdiffKsameWorld bdiff = Parallel::make_bdiff_ksame_world(); + bdiff.reduce_across_bdiff_ksame(f_en.demet); + kmesh.reduce_across_pools(f_en.demet); } else if (Occupy::fixed_occupations) { diff --git a/source/source_estate/module_charge/charge_extra.cpp b/source/source_estate/module_charge/charge_extra.cpp index 7513469a31b..1b64ee5c80d 100644 --- a/source/source_estate/module_charge/charge_extra.cpp +++ b/source/source_estate/module_charge/charge_extra.cpp @@ -4,7 +4,7 @@ #include "source_base/global_variable.h" #include "source_base/timer.h" #include "source_base/tool_threading.h" -#include "source_io/module_output/cube_io.h" +#include "source_cell/cube_io.h" Charge_Extra::Charge_Extra() { diff --git a/source/source_estate/module_charge/charge_init.cpp b/source/source_estate/module_charge/charge_init.cpp index 672d800f718..d3ce9534e51 100644 --- a/source/source_estate/module_charge/charge_init.cpp +++ b/source/source_estate/module_charge/charge_init.cpp @@ -13,7 +13,7 @@ #include "source_base/tool_threading.h" #include "source_cell/magnetism.h" #include "source_base/parallel_grid.h" -#include "source_io/module_output/cube_io.h" +#include "source_cell/cube_io.h" #include "source_estate/rhog_io.h" #include "source_io/module_wf/read_wf2rho_pw.h" #include "source_io/module_restart/restart.h" diff --git a/source/source_estate/occupy.cpp b/source/source_estate/occupy.cpp index 0a492735c6d..9699d407a70 100644 --- a/source/source_estate/occupy.cpp +++ b/source/source_estate/occupy.cpp @@ -2,8 +2,8 @@ #include "source_base/constants.h" #include "source_base/mymath.h" -#include "source_base/parallel_reduce.h" -#include "source_io/module_parameter/parameter.h" +#include "source_base/module_parallel/para_bridge.h" +#include "source_base/module_parallel/para_kmesh_world.h" Occupy::Occupy() { @@ -126,8 +126,10 @@ void Occupy::decision(const std::string& name, const std::string& smearing_metho * @param ekb the array save the band energy. * @param ef output: the highest occupied Kohn-Sham level. * @param wg output: weight for each k, each band. + * @param nspin number of spin components: 1 (spin-degenerate), 2 (collinear) or 4 (non-collinear). * @param is the spin index now. * @param isk distinguish k point belong to which spin. + * @param kmesh k-point parallel domain for cross-pool reduction. */ void Occupy::iweights( const int nks, @@ -138,17 +140,22 @@ void Occupy::iweights( const ModuleBase::matrix& ekb, double& ef, ModuleBase::matrix& wg, + const int nspin, const int& is, //<- is should be -1, 0, or 1. -1 means set all spins, and 0 means spin up, 1 means spin down. - const std::vector& isk) + const std::vector& isk, + const Parallel::ParaKmeshWorld& kmesh) { - assert(is < 2); + assert(nspin == 1 || nspin == 2 || nspin == 4); + assert(is >= -1 && is < 2); double degspin = 2.0; - if (PARAM.inp.nspin == 4) { + if (nspin == 4) + { degspin = 1.0; -} - if (is != -1) { + } + if (is != -1) + { degspin = 1.0; -} + } double ib_mind = nelec / degspin; int ib_min = std::ceil(ib_mind); @@ -163,7 +170,7 @@ void Occupy::iweights( for (int ik = 0; ik < nks; ++ik) { // when NSPIN=2, only calculate spin up or spin down with TWO_FERMI mode(nupdown != 0) - if (PARAM.inp.nspin == 2 && isk[ik] != is && is != -1) + if (nspin == 2 && isk[ik] != is && is != -1) { continue; } @@ -182,9 +189,7 @@ void Occupy::iweights( } } } - #ifdef __MPI - Parallel_Reduce::reduce_max(ef); - #endif + kmesh.reduce_max_across_pools(ef); return; } @@ -203,6 +208,7 @@ void Occupy::iweights( * @param wg output: weight of each band at each k point * @param is spin * @param isk array to point out each k belong to which spin + * @param kmesh k-point parallel domain for cross-pool reduction. */ void Occupy::gweights(const int nks, const std::vector& wk, @@ -215,24 +221,26 @@ void Occupy::gweights(const int nks, double& demet, ModuleBase::matrix& wg, const int& is, - const std::vector& isk) + const std::vector& isk, + const Parallel::ParaKmeshWorld& kmesh) { // ModuleBase::TITLE("Occupy","gweights"); //=============================== // Calculate the Fermi energy ef //=============================== // call efermig - Occupy::efermig(ekb, nband, nks, nelec, wk, smearing_sigma, ngauss, ef, is, isk); + Occupy::efermig(ekb, nband, nks, nelec, wk, smearing_sigma, ngauss, ef, is, isk, kmesh); demet = 0.0; for (int ik = 0; ik < nks; ik++) { // mohan add 2011-04-03 - if (is != -1 && is != isk[ik]) { + if (is != -1 && is != isk[ik]) + { continue; -} + } - for (int ib = 0; ib < PARAM.globalv.nbands_l; ib++) + for (int ib = 0; ib < nband; ib++) { //================================ // Calculate the gaussian weights @@ -266,6 +274,7 @@ void Occupy::gweights(const int nks, * @param ef output: fermi level * @param is spin * @param isk array to point out each k belong to which spin + * @param kmesh k-point parallel domain for cross-pool reduction. */ void Occupy::efermig(const ModuleBase::matrix& ekb, const int nband, @@ -276,7 +285,8 @@ void Occupy::efermig(const ModuleBase::matrix& ekb, const int ngauss, double& ef, const int& is, - const std::vector& isk) + const std::vector& isk, + const Parallel::ParaKmeshWorld& kmesh) { // ModuleBase::TITLE("Occupy","efermig"); //================================================================== @@ -309,10 +319,8 @@ void Occupy::efermig(const ModuleBase::matrix& ekb, eup += 2 * smearing_sigma; elw -= 2 * smearing_sigma; // find min and max across pools - #ifdef __MPI - Parallel_Reduce::reduce_max(eup); - Parallel_Reduce::reduce_min(elw); - #endif + kmesh.reduce_max_across_pools(eup); + kmesh.reduce_min_across_pools(elw); //================= // Bisection method //================= @@ -320,8 +328,8 @@ void Occupy::efermig(const ModuleBase::matrix& ekb, int changetime = 0; while (true) { - const double sumkup = Occupy::sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, eup, is, isk); - const double sumklw = Occupy::sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, elw, is, isk); + const double sumkup = Occupy::sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, eup, is, isk, kmesh); + const double sumklw = Occupy::sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, elw, is, isk, kmesh); if (changetime > 1000) { @@ -360,7 +368,7 @@ void Occupy::efermig(const ModuleBase::matrix& ekb, // change ef value //====================== ef = (eup + elw) / 2.0; - const double sumkmid = sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, ef, is, isk); + const double sumkmid = sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, ef, is, isk, kmesh); if (std::abs(sumkmid - nelec) < eps) { @@ -390,6 +398,7 @@ void Occupy::efermig(const ModuleBase::matrix& ekb, * @param e a givern energy * @param is spin * @param isk array to point out each k belong to which spin + * @param kmesh k-point parallel domain for cross-pool reduction. * @return (double) the number of states */ double Occupy::sumkg(const ModuleBase::matrix& ekb, @@ -400,15 +409,17 @@ double Occupy::sumkg(const ModuleBase::matrix& ekb, const int ngauss, const double& e, const int& is, - const std::vector& isk) + const std::vector& isk, + const Parallel::ParaKmeshWorld& kmesh) { // ModuleBase::TITLE("Occupy","sumkg"); double sum2 = 0.0; for (int ik = 0; ik < nks; ik++) { - if (is != -1 && is != isk[ik]) { + if (is != -1 && is != isk[ik]) + { continue; -} + } double sum1 = 0.0; for (int ib = 0; ib < nband; ib++) @@ -421,14 +432,14 @@ double Occupy::sumkg(const ModuleBase::matrix& ekb, sum2 += wk[ik] * sum1; } - // GlobalV::ofs_running << "\n sum2 before reduce = " << sum2 << std::endl; - -#ifdef __MPI - const int npool = GlobalV::KPAR * PARAM.inp.bndpar; - Parallel_Reduce::reduce_double_allpool(npool, GlobalV::NPROC_IN_POOL, sum2); -#endif - - // GlobalV::ofs_running << "\n sum2 after reduce = " << sum2 << std::endl; + // Two-step reduction, order matters: + // 1. Band dimension: BPCG shards the band range across the band + // groups of this k-pool, so combine the per-window partial sums + // first (no-op with a single band group). + Parallel::ParaBdiffKsameWorld bdiff = Parallel::make_bdiff_ksame_world(); + bdiff.reduce_across_bdiff_ksame(sum2); + // 2. k dimension: exactly one contribution per k-pool. + kmesh.reduce_across_pools(sum2); return sum2; } @@ -487,7 +498,7 @@ double Occupy::wgauss(const double& x, const int n) //==================== wga = 0.5 * (1 - erf(-x)); // wga = gauss_freq(x * ModuleBase::SQRT2); - // std::cout<<"\n x="<& isk); + const std::vector& isk, + const Parallel::ParaKmeshWorld& kmesh); static void gweights(const int nks, const std::vector& wk, @@ -57,43 +64,43 @@ class Occupy double& demet, ModuleBase::matrix& wg, const int& is, - const std::vector& isk); + const std::vector& isk, + const Parallel::ParaKmeshWorld& kmesh); - static void tweights(const int nks,const int nspin,const int nband, - const double &nelec,const int ntetra, - const ModuleBase::matrix &tetra,double **ekb,double &ef,ModuleBase::matrix &wg); + static void tweights(const int nks, const int nspin, const int nband, + const double& nelec, const int ntetra, + const ModuleBase::matrix& tetra, double** ekb, double& ef, ModuleBase::matrix& wg); static double wsweight(const ModuleBase::Vector3 &r, ModuleBase::Vector3 *rws,const int nrws); -private: - static void efermig(const ModuleBase::matrix& ekb, - const int nbnd, - const int nks, - const double& nelec, - const std::vector& wk, - const double& smearing_sigma, - const int ngauss, - double& ef, - const int& is, - const std::vector& isk); - - static double sumkg(const ModuleBase::matrix& ekb, - const int nband, - const int nks, - const std::vector& wk, - const double& smearing_sigma, - const int ngauss, - const double& e, - const int& is, - const std::vector& isk); - - static double wgauss(const double& x, const int n); + static void efermig(const ModuleBase::matrix& ekb, + const int nbnd, + const int nks, + const double& nelec, + const std::vector& wk, + const double& smearing_sigma, + const int ngauss, + double& ef, + const int& is, + const std::vector& isk, + const Parallel::ParaKmeshWorld& kmesh); + + static double sumkg(const ModuleBase::matrix& ekb, + const int nband, + const int nks, + const std::vector& wk, + const double& smearing_sigma, + const int ngauss, + const double& e, + const int& is, + const std::vector& isk, + const Parallel::ParaKmeshWorld& kmesh); + + static double wgauss(const double& x, const int n); + + static double w1gauss(const double& x, const int n); - static double w1gauss(const double& x, const int n); - - //============================ - // Needed in tweights - //============================ +private: static void efermit(double** ekb, const int nband, const int nks, diff --git a/source/source_estate/test/CMakeLists.txt b/source/source_estate/test/CMakeLists.txt index 13be13db5cc..53ec125444d 100644 --- a/source/source_estate/test/CMakeLists.txt +++ b/source/source_estate/test/CMakeLists.txt @@ -28,8 +28,8 @@ endif() AddTest( TARGET MODULE_ESTATE_elecstate_occupy - LIBS parameter base device - SOURCES elecstate_occupy_test.cpp ../occupy.cpp + LIBS base device + SOURCES test_occupy.cpp ../occupy.cpp ) AddTest( @@ -111,7 +111,7 @@ AddTest( AddTest( TARGET MODULE_ESTATE_charge_extra LIBS parameter base device cell_info - SOURCES charge_extra_test.cpp ../module_charge/charge_extra.cpp ../../source_io/module_output/read_cube.cpp ../../source_io/module_output/write_cube.cpp + SOURCES charge_extra_test.cpp ../module_charge/charge_extra.cpp ../../source_cell/read_cube.cpp ../../source_cell/write_cube.cpp ../../source_base/module_fft/fft_bundle.cpp ../../source_base/module_fft/fft_cpu.cpp ) diff --git a/source/source_estate/test/elecstate_occupy_test.cpp b/source/source_estate/test/test_occupy.cpp similarity index 88% rename from source/source_estate/test/elecstate_occupy_test.cpp rename to source/source_estate/test/test_occupy.cpp index 8f034ee93b6..ee14467858b 100644 --- a/source/source_estate/test/elecstate_occupy_test.cpp +++ b/source/source_estate/test/test_occupy.cpp @@ -1,9 +1,8 @@ #include #include "gmock/gmock.h" #include "gtest/gtest.h" -#define private public -#include "source_io/module_parameter/parameter.h" -#undef private +#include "source_base/module_parallel/para_kmesh_world.h" +#include "source_estate/occupy.h" /*************************************************************** * unit test of class Occupy @@ -14,9 +13,6 @@ * - Occupy::Occupy() * - Occupy::decision() */ -#define private public -#include "source_estate/occupy.h" -#undef private class OccupyTest : public ::testing::Test { protected: @@ -182,21 +178,20 @@ TEST_F(OccupyTest, DecisionArbitrary) TEST_F(OccupyTest, IweightsNOSPIN) { - PARAM.input.nspin = 1; double ef = 0.0; ModuleBase::matrix wg(1, 1); std::vector wk(1, 2.0); ModuleBase::matrix ekb(1, 1); std::vector isk(1); ekb(0, 0) = 0.1; - occupy.iweights(1, wk, 1, 0, 2.0, ekb, ef, wg, 0, isk); + Parallel::ParaKmeshWorld kmesh(1, 1); + occupy.iweights(1, wk, 1, 0, 2.0, ekb, ef, wg, 1, 0, isk, kmesh); EXPECT_DOUBLE_EQ(wg(0, 0), 2.0); EXPECT_DOUBLE_EQ(ef, 0.1); } TEST_F(OccupyTest, IweightsSPIN) { - PARAM.input.nspin = 2; double ef_up = 0.0; double ef_dw = 0.0; ModuleBase::matrix wg(2, 1); @@ -207,8 +202,9 @@ TEST_F(OccupyTest, IweightsSPIN) isk[1] = 1; ekb(0, 0) = 0.1; ekb(1, 0) = 0.2; - occupy.iweights(2, wk, 1, 0, 1.0, ekb, ef_up, wg, 0, isk); - occupy.iweights(2, wk, 1, 0, 1.0, ekb, ef_dw, wg, 1, isk); + Parallel::ParaKmeshWorld kmesh(2, 2); + occupy.iweights(2, wk, 1, 0, 1.0, ekb, ef_up, wg, 2, 0, isk, kmesh); + occupy.iweights(2, wk, 1, 0, 1.0, ekb, ef_dw, wg, 2, 1, isk, kmesh); EXPECT_DOUBLE_EQ(wg(0, 0), 1.0); EXPECT_DOUBLE_EQ(wg(1, 0), 1.0); EXPECT_DOUBLE_EQ(ef_up, 0.1); @@ -217,7 +213,6 @@ TEST_F(OccupyTest, IweightsSPIN) TEST_F(OccupyTest, IweightsWarning) { - PARAM.input.nspin = 1; double ef = 0.0; ModuleBase::matrix wg(1, 1); std::vector wk(1, 2.0); @@ -225,8 +220,9 @@ TEST_F(OccupyTest, IweightsWarning) std::vector isk(1); ekb(0, 0) = 0.1; + Parallel::ParaKmeshWorld kmesh(1, 1); testing::internal::CaptureStdout(); - EXPECT_EXIT(occupy.iweights(1, wk, 1, 0, 1.0, ekb, ef, wg, -1, isk);, ::testing::ExitedWithCode(1), ""); + EXPECT_EXIT(occupy.iweights(1, wk, 1, 0, 1.0, ekb, ef, wg, 1, -1, isk, kmesh);, ::testing::ExitedWithCode(1), ""); output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("It is not a semiconductor or insulator. Please do not set 'smearing_method=fixed', and try other options.")); } @@ -260,7 +256,8 @@ TEST_F(OccupyTest, Sumkg) double e = 0.0; int is = 0; std::vector isk = {0, 0}; - EXPECT_DOUBLE_EQ(occupy.sumkg(ekb, 1, 1, wk, smearing_sigma, ngauss, e, is, isk), 1.0); + Parallel::ParaKmeshWorld kmesh(1, 1); + EXPECT_DOUBLE_EQ(occupy.sumkg(ekb, 1, 1, wk, smearing_sigma, ngauss, e, is, isk, kmesh), 1.0); } TEST_F(OccupyTest, Efermig) @@ -274,7 +271,8 @@ TEST_F(OccupyTest, Efermig) int is = 0; std::vector isk = {0, 0}; double ef = 0.0; - occupy.efermig(ekb, 1, 1, 1.0, wk, smearing_sigma, ngauss, ef, is, isk); + Parallel::ParaKmeshWorld kmesh(1, 1); + occupy.efermig(ekb, 1, 1, 1.0, wk, smearing_sigma, ngauss, ef, is, isk, kmesh); EXPECT_NEAR(ef, -0.5, 1e-13); } @@ -290,10 +288,12 @@ TEST_F(OccupyTest, Gweights) std::vector isk = {0, 0}; double ef = 0.0; ModuleBase::matrix wg(1, 1); - wg(0, 0) = 1.0; double demet = 0.0; - occupy.gweights(1, wk, 1, 1.0, smearing_sigma, ngauss, ekb, ef, demet, wg, is, isk); - EXPECT_NEAR(ef, -0.5, 1e-13); - EXPECT_NEAR(demet, 0.0, 1e-13); - EXPECT_NEAR(wg(0, 0), 1.0, 1e-13); + // Half-filled single band: the Fermi energy stays at the band energy, the + // occupation is 1/2 and demet equals sigma * w1gauss(0, 0). + Parallel::ParaKmeshWorld kmesh(1, 1); + occupy.gweights(1, wk, 1, 0.5, smearing_sigma, ngauss, ekb, ef, demet, wg, is, isk, kmesh); + EXPECT_NEAR(ef, -1.0, 1e-13); + EXPECT_NEAR(wg(0, 0), 0.5, 1e-13); + EXPECT_NEAR(demet, smearing_sigma * (-0.28209479177387814), 1e-13); } diff --git a/source/source_estate/write_elecstat_pot.cpp b/source/source_estate/write_elecstat_pot.cpp index 15a5363ed62..2ef5a2b9f69 100644 --- a/source/source_estate/write_elecstat_pot.cpp +++ b/source/source_estate/write_elecstat_pot.cpp @@ -1,12 +1,13 @@ #include "source_base/element_name.h" #include "source_base/timer.h" -#include "source_io/module_parameter/parameter.h" #include "source_estate/module_pot/h_hartree_pw.h" #include "source_estate/module_pot/efield.h" -#include "source_io/module_output/cube_io.h" -#include "source_io/module_output/output_log.h" +#include "source_cell/cube_io.h" +#include "source_cell/output_log.h" #include "write_elecstat_pot.h" +#include + namespace ModuleIO { @@ -22,17 +23,19 @@ void write_elecstat_pot( const UnitCell* ucell, const double* v_eff, const surchem& solvent, - const int precision) + const int precision, + const int nspin, + const bool efield_flag, + const bool dip_cor_flag, + const bool imp_sol, + const bool two_fermi) { ModuleBase::TITLE("ModuleIO", "write_elecstat_pot"); ModuleBase::timer::start("ModuleIO", "write_elecstat_pot"); - std::vector v_elecstat(rho_basis->nrxx, 0.0); + assert(nspin == 1 || nspin == 2 || nspin == 4); - const int nspin = PARAM.inp.nspin; - const int efield = PARAM.inp.efield_flag; - const int dip_corr = PARAM.inp.dip_cor_flag; - const bool imp_sol = PARAM.inp.imp_sol; + std::vector v_elecstat(rho_basis->nrxx, 0.0); //========================================== // Hartree potential @@ -44,7 +47,7 @@ void write_elecstat_pot( //! Dipole correction //========================================== ModuleBase::matrix v_efield; - if (efield>0 && dip_corr>0) + if (efield_flag && dip_cor_flag) { v_efield.create(nspin, rho_basis->nrxx); v_efield = elecstate::Efield::add_efield(*ucell, @@ -62,11 +65,11 @@ void write_elecstat_pot( // the spin index is 0 v_elecstat[ir] = vh(0, ir) + v_eff[ir]; - if (efield>0 && dip_corr>0) + if (efield_flag && dip_cor_flag) { v_elecstat[ir] += v_efield(0, ir); } - if(imp_sol == true) + if(imp_sol) { v_elecstat[ir] += solvent.delta_phi[ir]; } @@ -103,7 +106,7 @@ void write_elecstat_pot( ucell, precision, out_fermi, - PARAM.globalv.two_fermi, + two_fermi, false); ModuleBase::timer::end("ModuleIO", "write_elecstat_pot"); diff --git a/source/source_estate/write_elecstat_pot.h b/source/source_estate/write_elecstat_pot.h index bee575b95c3..967a6200e8b 100644 --- a/source/source_estate/write_elecstat_pot.h +++ b/source/source_estate/write_elecstat_pot.h @@ -20,7 +20,12 @@ namespace ModuleIO /// @param ucell_ /// @param v_eff_fixed /// @param solvent: for solvation model -/// #param precision: output precision +/// @param precision: output precision +/// @param nspin: number of spin channels (1, 2, or 4) +/// @param efield_flag: whether electric field is applied +/// @param dip_cor_flag: whether dipole correction is applied +/// @param imp_sol: whether implicit solvation model is used +/// @param two_fermi: whether two Fermi levels are used void write_elecstat_pot( #ifdef __MPI const int& bz, @@ -33,7 +38,12 @@ void write_elecstat_pot( const UnitCell* ucell_, const double* v_eff_fixed, const surchem& solvent, - const int precision); + const int precision, + const int nspin, + const bool efield_flag, + const bool dip_cor_flag, + const bool imp_sol, + const bool two_fermi); } // namespace ModuleIO diff --git a/source/source_estate/write_init.cpp b/source/source_estate/write_init.cpp index 58cdb534287..ccd8e125a98 100644 --- a/source/source_estate/write_init.cpp +++ b/source/source_estate/write_init.cpp @@ -15,7 +15,7 @@ // ===================================================================== #include "source_estate/write_init.h" -#include "source_io/module_output/cube_io.h" +#include "source_cell/cube_io.h" #include "source_base/tool_quit.h" #include diff --git a/source/source_io/CMakeLists.txt b/source/source_io/CMakeLists.txt index 85ffd82d545..33a757c5836 100644 --- a/source/source_io/CMakeLists.txt +++ b/source/source_io/CMakeLists.txt @@ -22,13 +22,10 @@ list(APPEND objects module_bessel/numerical_basis_jyjy.cpp module_bessel/numerical_descriptor.cpp module_output/print_info.cpp - module_output/read_cube.cpp module_wf/read_wfc_pw.cpp module_wf/read_wf2rho_pw.cpp module_restart/restart.cpp module_wf/write_wfc_pw.cpp - module_output/write_pao.cpp - module_output/write_cube.cpp module_elf/write_elf.cpp module_dipole/write_dipole.cpp module_ml/write_mlkedf_desc.cpp @@ -37,7 +34,6 @@ list(APPEND objects module_efield/td_efield_io.cpp module_efield/td_vector_pot_io.cpp module_chgpot/write_libxc_r.cpp - module_output/output_log.cpp module_json/para_json.cpp parse_args.cpp input_help.cpp @@ -64,7 +60,6 @@ if(ENABLE_LCAO) module_dos/write_dos_lcao.cpp module_dos/cal_pdos_gamma.cpp module_dos/cal_pdos_multik.cpp - module_output/write_orb_info.cpp module_energy/write_proj_band_lcao.cpp module_chgpot/get_pchg_lcao.cpp module_wf/get_wf_lcao.cpp diff --git a/source/source_io/module_chgpot/get_pchg_lcao.cpp b/source/source_io/module_chgpot/get_pchg_lcao.cpp index ae11f7b6e74..50e7d09d5a1 100644 --- a/source/source_io/module_chgpot/get_pchg_lcao.cpp +++ b/source/source_io/module_chgpot/get_pchg_lcao.cpp @@ -3,7 +3,7 @@ #include "source_estate/module_charge/symm_rho.h" #include "source_estate/module_dm/cal_dm_psi.h" #include "source_hamilt/module_gint/gint_interface.h" -#include "source_io/module_output/cube_io.h" +#include "source_cell/cube_io.h" #include #include diff --git a/source/source_io/module_chgpot/get_pchg_pw.h b/source/source_io/module_chgpot/get_pchg_pw.h index d7ad0b7cc80..56bd208a22b 100644 --- a/source/source_io/module_chgpot/get_pchg_pw.h +++ b/source/source_io/module_chgpot/get_pchg_pw.h @@ -5,7 +5,7 @@ #include "source_base/parallel_comm.h" #include "source_estate/module_charge/symm_rho.h" #include "source_io/module_output/band_parallel_output.h" -#include "source_io/module_output/cube_io.h" +#include "source_cell/cube_io.h" namespace ModuleIO { diff --git a/source/source_io/module_chgpot/write_libxc_r.cpp b/source/source_io/module_chgpot/write_libxc_r.cpp index d13e7563850..9f49474b031 100644 --- a/source/source_io/module_chgpot/write_libxc_r.cpp +++ b/source/source_io/module_chgpot/write_libxc_r.cpp @@ -12,7 +12,7 @@ #include "source_estate/module_charge/charge.h" #include "source_basis/module_pw/pw_basis_big.h" #include "source_basis/module_pw/pw_basis.h" -#include "source_io/module_output/cube_io.h" +#include "source_cell/cube_io.h" #include "source_base/global_variable.h" #include "source_io/module_parameter/parameter.h" #include "source_base/timer.h" diff --git a/source/source_io/module_ctrl/ctrl_output_fp.cpp b/source/source_io/module_ctrl/ctrl_output_fp.cpp index 318c571e3c4..faaaab9eba4 100644 --- a/source/source_io/module_ctrl/ctrl_output_fp.cpp +++ b/source/source_io/module_ctrl/ctrl_output_fp.cpp @@ -1,10 +1,11 @@ #include "ctrl_output_fp.h" // use ctrl_output_fp() -#include "../module_output/cube_io.h" // use write_vdata_palgrid +#include "source_cell/cube_io.h" // use write_vdata_palgrid #include "../module_dipole/dipole_io.h" // use write_dipole #include "source_estate/module_charge/symm_rho.h" // use Symmetry_rho #include "source_hamilt/module_xc/xc_functional.h" // use XC_Functional #include "source_estate/write_elecstat_pot.h" // use write_elecstat_pot #include "source_io/module_elf/write_elf.h" +#include "source_io/module_parameter/input_parameter.h" #ifdef __LIBXC #include "source_io/module_chgpot/write_libxc_r.h" @@ -20,35 +21,36 @@ void ctrl_output_fp(UnitCell& ucell, Charge& chr, surchem& solvent, Parallel_Grid& para_grid, - const int istep) + const int istep, + const Input_para& inp) { ModuleBase::TITLE("ModuleIO", "ctrl_output_fp"); ModuleBase::timer::start("ModuleIO", "ctrl_output_fp"); - const bool out_app_flag = PARAM.inp.out_app_flag; + const bool out_app_flag = inp.out_app_flag; const bool gamma_only = PARAM.globalv.gamma_only_local; - const int nspin = PARAM.inp.nspin; + const int nspin = inp.nspin; const std::string global_out_dir = PARAM.globalv.global_out_dir; // print out the 'g' index when istep_in != -1 int istep_in = -1; - if (PARAM.inp.esolver_type != "tddft" && PARAM.inp.out_freq_ion > 0) // default value of out_freq_ion is 0 + if (inp.esolver_type != "tddft" && inp.out_freq_ion > 0) // default value of out_freq_ion is 0 { - if (istep % PARAM.inp.out_freq_ion == 0) + if (istep % inp.out_freq_ion == 0) { istep_in = istep; } } - else if (PARAM.inp.esolver_type == "tddft" && PARAM.inp.out_freq_td > 0) // default value of out_freq_td is 0 + else if (inp.esolver_type == "tddft" && inp.out_freq_td > 0) // default value of out_freq_td is 0 { - if (istep % PARAM.inp.out_freq_td == 0) + if (istep % inp.out_freq_td == 0) { istep_in = istep; } } std::string geom_block; - bool should_output = (PARAM.inp.out_freq_ion == 0); + bool should_output = (inp.out_freq_ion == 0); if (istep_in >= 0) { geom_block = "g" + std::to_string(istep + 1); @@ -56,7 +58,7 @@ void ctrl_output_fp(UnitCell& ucell, } // 4) write charge density - if (PARAM.inp.out_chg[0] > 0 && should_output) + if (inp.out_chg[0] > 0 && should_output) { for (int is = 0; is < nspin; ++is) { @@ -82,7 +84,7 @@ void ctrl_output_fp(UnitCell& ucell, fn, pelec->eferm.get_efval(is), &(ucell), - PARAM.inp.out_chg[1], + inp.out_chg[1], 1, PARAM.globalv.two_fermi, false); @@ -110,7 +112,7 @@ void ctrl_output_fp(UnitCell& ucell, } // 5) write potential - if ((PARAM.inp.out_pot[0] == 1 || PARAM.inp.out_pot[0] == 3) && should_output) + if ((inp.out_pot[0] == 1 || inp.out_pot[0] == 3) && should_output) { for (int is = 0; is < nspin; is++) { @@ -136,13 +138,13 @@ void ctrl_output_fp(UnitCell& ucell, fn, 0.0, // efermi &(ucell), - PARAM.inp.out_pot[1], // precision + inp.out_pot[1], // precision 0, // out_fermi PARAM.globalv.two_fermi, false); } } - else if (PARAM.inp.out_pot[0] == 2 && should_output) + else if (inp.out_pot[0] == 2 && should_output) { std::string fn = PARAM.globalv.global_out_dir + "potes"; fn += geom_block + ".cube"; @@ -159,11 +161,16 @@ void ctrl_output_fp(UnitCell& ucell, &(ucell), pelec->pot->get_fixed_v(), solvent, - PARAM.inp.out_pot[1]); + inp.out_pot[1], + nspin, + inp.efield_flag, + inp.dip_cor_flag, + inp.imp_sol, + PARAM.globalv.two_fermi); } // 6) write ELF - if (PARAM.inp.out_elf[0] > 0 && should_output) + if (inp.out_elf[0] > 0 && should_output) { chr.cal_elf = true; Symmetry_rho srho; @@ -181,16 +188,16 @@ void ctrl_output_fp(UnitCell& ucell, pw_rhod, para_grid, &(ucell), - PARAM.inp.out_elf[1], + inp.out_elf[1], geom_block, PARAM.globalv.two_fermi); } #ifdef __LIBXC // 7) write xc(r) - if (PARAM.inp.out_xc_r[0] >= 0 && should_output) + if (inp.out_xc_r[0] >= 0 && should_output) { - ModuleIO::write_libxc_r(PARAM.inp.out_xc_r[0], + ModuleIO::write_libxc_r(inp.out_xc_r[0], XC_Functional::get_func_id(), pw_rhod->nrxx, // number of real-space grid ucell.omega, // volume of cell @@ -202,7 +209,7 @@ void ctrl_output_fp(UnitCell& ucell, #endif // 8) write dipole moment - if (PARAM.inp.out_dipole == 1 && should_output) + if (inp.out_dipole == 1 && should_output) { for (int is = 0; is < nspin; ++is) { diff --git a/source/source_io/module_ctrl/ctrl_output_fp.h b/source/source_io/module_ctrl/ctrl_output_fp.h index ec9f4e20e78..a1203cd107d 100644 --- a/source/source_io/module_ctrl/ctrl_output_fp.h +++ b/source/source_io/module_ctrl/ctrl_output_fp.h @@ -3,6 +3,8 @@ #include "source_estate/elecstate_lcao.h" +struct Input_para; + namespace ModuleIO { @@ -13,7 +15,8 @@ void ctrl_output_fp(UnitCell& ucell, Charge& chr, surchem& solvent, Parallel_Grid& para_grid, - const int istep); + const int istep, + const Input_para& inp); } #endif diff --git a/source/source_io/module_dos/cal_ldos.cpp b/source/source_io/module_dos/cal_ldos.cpp index b6ae29fba84..195ff51c97d 100644 --- a/source/source_io/module_dos/cal_ldos.cpp +++ b/source/source_io/module_dos/cal_ldos.cpp @@ -1,7 +1,7 @@ #include "cal_ldos.h" #include "cal_dos.h" -#include "../module_output/cube_io.h" +#include "source_cell/cube_io.h" #include "source_estate/module_dm/cal_dm_psi.h" #include "source_hamilt/module_gint/gint_interface.h" #include "source_base/module_device/memory_op.h" diff --git a/source/source_io/module_dos/cal_pdos_gamma.cpp b/source/source_io/module_dos/cal_pdos_gamma.cpp index 50acad4f265..ad7ee387ce0 100644 --- a/source/source_io/module_dos/cal_pdos_gamma.cpp +++ b/source/source_io/module_dos/cal_pdos_gamma.cpp @@ -6,7 +6,7 @@ #include "source_base/global_function.h" #include "source_base/global_variable.h" #include "source_lcao/hamilt_lcao.h" -#include "source_io/module_output/write_orb_info.h" +#include "source_cell/write_orb_info.h" void ModuleIO::cal_pdos( diff --git a/source/source_io/module_dos/cal_pdos_multik.cpp b/source/source_io/module_dos/cal_pdos_multik.cpp index b7766b22f4a..87bcc0aac22 100644 --- a/source/source_io/module_dos/cal_pdos_multik.cpp +++ b/source/source_io/module_dos/cal_pdos_multik.cpp @@ -3,7 +3,7 @@ #include "source_base/parallel_reduce.h" #include "source_base/module_external/blas_connector.h" #include "source_base/module_external/scalapack_connector.h" -#include "source_io/module_output/write_orb_info.h" +#include "source_cell/write_orb_info.h" #include "source_base/global_function.h" #include "source_base/global_variable.h" #include "source_lcao/hamilt_lcao.h" diff --git a/source/source_io/module_elf/write_elf.cpp b/source/source_io/module_elf/write_elf.cpp index ea6e9b88050..8c9eba0126f 100644 --- a/source/source_io/module_elf/write_elf.cpp +++ b/source/source_io/module_elf/write_elf.cpp @@ -1,5 +1,5 @@ #include "write_elf.h" -#include "source_io/module_output/cube_io.h" +#include "source_cell/cube_io.h" #ifdef _OPENMP #include #endif diff --git a/source/source_io/module_energy/write_proj_band_lcao.cpp b/source/source_io/module_energy/write_proj_band_lcao.cpp index ecf83e6a4d0..80c4e38dbaf 100644 --- a/source/source_io/module_energy/write_proj_band_lcao.cpp +++ b/source/source_io/module_energy/write_proj_band_lcao.cpp @@ -6,7 +6,7 @@ #include "source_base/module_external/scalapack_connector.h" #include "source_base/timer.h" #include "source_cell/module_neighbor/sltk_atom_arrange.h" -#include "source_io/module_output/write_orb_info.h" +#include "source_cell/write_orb_info.h" #include "source_lcao/hamilt_lcao.h" template<> diff --git a/source/source_io/module_wf/get_wf_lcao.cpp b/source/source_io/module_wf/get_wf_lcao.cpp index ea1189acaee..54167861985 100644 --- a/source/source_io/module_wf/get_wf_lcao.cpp +++ b/source/source_io/module_wf/get_wf_lcao.cpp @@ -2,7 +2,7 @@ #include "source_hamilt/module_gint/gint_env_gamma.h" #include "source_hamilt/module_gint/gint_env_k.h" -#include "source_io/module_output/cube_io.h" +#include "source_cell/cube_io.h" #include #include diff --git a/source/source_io/test/CMakeLists.txt b/source/source_io/test/CMakeLists.txt index 309aba37986..eca660d7e82 100644 --- a/source/source_io/test/CMakeLists.txt +++ b/source/source_io/test/CMakeLists.txt @@ -113,7 +113,7 @@ add_test(NAME MODULE_IO_write_wfc_nao_para AddTest( TARGET MODULE_IO_write_orb_info LIBS parameter base device cell_info - SOURCES write_orb_info_test.cpp ../module_output/write_orb_info.cpp + SOURCES write_orb_info_test.cpp ../../source_cell/write_orb_info.cpp ) AddTest( @@ -137,7 +137,7 @@ AddTest( AddTest( TARGET MODULE_IO_output_log_test LIBS parameter base device - SOURCES ../module_output/output_log.cpp outputlog_test.cpp ../../source_basis/module_pw/test/test_tool.cpp + SOURCES ../../source_cell/output_log.cpp outputlog_test.cpp ../../source_basis/module_pw/test/test_tool.cpp ) if(ENABLE_LCAO) diff --git a/source/source_io/test/outputlog_test.cpp b/source/source_io/test/outputlog_test.cpp index b84a9ef4a0e..ef91ae4a482 100644 --- a/source/source_io/test/outputlog_test.cpp +++ b/source/source_io/test/outputlog_test.cpp @@ -10,7 +10,7 @@ #include "source_base/constants.h" #include "source_base/global_variable.h" -#include "source_io/module_output/output_log.h" +#include "source_cell/output_log.h" #ifdef __MPI #include "source_basis/module_pw/test/test_tool.h" diff --git a/source/source_io/test/write_orb_info_test.cpp b/source/source_io/test/write_orb_info_test.cpp index 1829cf69fc3..5509c40c37e 100644 --- a/source/source_io/test/write_orb_info_test.cpp +++ b/source/source_io/test/write_orb_info_test.cpp @@ -3,7 +3,7 @@ #define private public #include "source_io/module_parameter/parameter.h" #undef private -#include "source_io/module_output/write_orb_info.h" +#include "source_cell/write_orb_info.h" #include "source_cell/unitcell.h" #include "prepare_unitcell.h" #include "source_cell/read_pp_ucell.h" diff --git a/source/source_io/test_serial/CMakeLists.txt b/source/source_io/test_serial/CMakeLists.txt index dd87df580c5..43a339320a4 100644 --- a/source/source_io/test_serial/CMakeLists.txt +++ b/source/source_io/test_serial/CMakeLists.txt @@ -51,7 +51,7 @@ AddTest( AddTest( TARGET MODULE_IO_rho_io LIBS parameter base device cell_info - SOURCES rho_io_test.cpp ../module_output/read_cube.cpp ../module_output/write_cube.cpp + SOURCES rho_io_test.cpp ../../source_cell/read_cube.cpp ../../source_cell/write_cube.cpp ) AddTest( diff --git a/source/source_io/test_serial/rho_io_test.cpp b/source/source_io/test_serial/rho_io_test.cpp index 5939eecaef1..f3b600b7f8d 100644 --- a/source/source_io/test_serial/rho_io_test.cpp +++ b/source/source_io/test_serial/rho_io_test.cpp @@ -1,9 +1,9 @@ -#include "source_io/module_output/cube_io.h" +#include "source_cell/cube_io.h" #include "gmock/gmock.h" #include "gtest/gtest.h" #include "source_base/global_variable.h" -#include "source_io/module_output/cube_io.h" +#include "source_cell/cube_io.h" #include "prepare_unitcell.h" #include "source_base/parallel_grid.h" diff --git a/source/source_lcao/force_stress_lcao.cpp b/source/source_lcao/force_stress_lcao.cpp index 895e03c1378..c5ce4f2d5db 100644 --- a/source/source_lcao/force_stress_lcao.cpp +++ b/source/source_lcao/force_stress_lcao.cpp @@ -3,7 +3,7 @@ #include "source_base/parallel_reduce.h" #include "source_lcao/module_dftu/dftu_nao.h" //Quxin add for DFT+U on 20201029 #include "source_lcao/module_dftu/dftu_nao_fs_k.h" -#include "source_io/module_output/output_log.h" +#include "source_cell/output_log.h" #include "source_io/module_parameter/parameter.h" // new #include "source_base/timer.h" diff --git a/source/source_lcao/module_lr/potentials/xc_kernel.cpp b/source/source_lcao/module_lr/potentials/xc_kernel.cpp index de7146db670..aa02bfca7e5 100644 --- a/source/source_lcao/module_lr/potentials/xc_kernel.cpp +++ b/source/source_lcao/module_lr/potentials/xc_kernel.cpp @@ -6,7 +6,7 @@ #include "source_lcao/module_lr/utils/lr_util_xc.hpp" #include #include -#include "source_io/module_output/cube_io.h" +#include "source_cell/cube_io.h" #ifdef __LIBXC #include #include "source_hamilt/module_xc/libxc_abacus.h" diff --git a/source/source_lcao/module_lr/utils/exciton_plotter.h b/source/source_lcao/module_lr/utils/exciton_plotter.h index 18ad7521233..53113519d4b 100644 --- a/source/source_lcao/module_lr/utils/exciton_plotter.h +++ b/source/source_lcao/module_lr/utils/exciton_plotter.h @@ -5,7 +5,7 @@ #include "source_cell/atom_spec.h" #include "source_cell/klist.h" #include "source_estate/module_dm/density_matrix.h" -#include "source_io/module_output/cube_io.h" +#include "source_cell/cube_io.h" #include "source_hamilt/module_gint/gint_interface.h" #include "source_lcao/module_lr/dm_trans/dm_trans.h" #include "source_lcao/module_lr/utils/lr_util.h" diff --git a/source/source_md/md_func.cpp b/source/source_md/md_func.cpp index 31d8111e24c..7e2b35178a2 100644 --- a/source/source_md/md_func.cpp +++ b/source/source_md/md_func.cpp @@ -2,7 +2,7 @@ #include "source_base/global_variable.h" #include "source_base/timer.h" -#include "source_io/module_output/output_log.h" +#include "source_cell/output_log.h" #include "source_io/module_parameter/parameter.h" #include diff --git a/source/source_md/test/CMakeLists.txt b/source/source_md/test/CMakeLists.txt index 6d2821ff91d..84e05122ff0 100644 --- a/source/source_md/test/CMakeLists.txt +++ b/source/source_md/test/CMakeLists.txt @@ -53,7 +53,7 @@ list(APPEND depend_files ../../source_cell/module_neighlist/domain_decomposition.cpp ../../source_cell/mdcell.cpp ../../source_base/output.cpp - ../../source_io/module_output/output_log.cpp + ../../source_cell/output_log.cpp ../../source_io/module_output/print_info.cpp ../../source_cell/cif_io.cpp ../../source_esolver/esolver_lj.cpp diff --git a/source/source_psi/psi_init_atomic.cpp b/source/source_psi/psi_init_atomic.cpp index 39db75dbbaf..04592d79a5c 100644 --- a/source/source_psi/psi_init_atomic.cpp +++ b/source/source_psi/psi_init_atomic.cpp @@ -10,7 +10,7 @@ #include "source_base/tool_quit.h" #include "source_base/timer.h" #include "source_base/global_variable.h" -#include "source_io/module_output/write_pao.h" +#include "source_cell/write_pao.h" template void psi_init_atomic::allocate_ps_table() diff --git a/source/source_psi/test/CMakeLists.txt b/source/source_psi/test/CMakeLists.txt index 79438154d99..69af9a2687f 100644 --- a/source/source_psi/test/CMakeLists.txt +++ b/source/source_psi/test/CMakeLists.txt @@ -16,7 +16,7 @@ AddTest( ../../source_cell/atom_spec.cpp ../../source_cell/test/support/mock_unitcell.cpp - ../../source_io/module_output/write_pao.cpp + ../../source_cell/write_pao.cpp ../../source_io/module_wf/read_wfc_pw.cpp ) endif() diff --git a/source/source_pw/module_ofdft/of_stress_pw.cpp b/source/source_pw/module_ofdft/of_stress_pw.cpp index 6b988baf786..4dade34559c 100644 --- a/source/source_pw/module_ofdft/of_stress_pw.cpp +++ b/source/source_pw/module_ofdft/of_stress_pw.cpp @@ -3,7 +3,7 @@ #include "source_base/timer.h" #include "source_base/tool_quit.h" #include "source_hamilt/module_vdw/vdw.h" -#include "source_io/module_output/output_log.h" +#include "source_cell/output_log.h" // Since the kinetic stress of OFDFT is calculated by kinetic functionals in esolver_of.cpp, here we regard it as an // input variable. diff --git a/source/source_pw/module_pwdft/force_pw.cpp b/source/source_pw/module_pwdft/force_pw.cpp index f02afc017b6..a87b30614ac 100644 --- a/source/source_pw/module_pwdft/force_pw.cpp +++ b/source/source_pw/module_pwdft/force_pw.cpp @@ -4,7 +4,7 @@ #include "source_base/parallel_reduce.h" #include "source_pw/module_pwdft/kernels/force_op.h" #include "source_io/module_parameter/parameter.h" -#include "source_io/module_output/output_log.h" +#include "source_cell/output_log.h" // new #include "source_base/complexmatrix.h" #include "source_base/libm/libm.h" diff --git a/source/source_pw/module_pwdft/stress_pw.cpp b/source/source_pw/module_pwdft/stress_pw.cpp index bcf6aff8db2..b0b6890ef36 100644 --- a/source/source_pw/module_pwdft/stress_pw.cpp +++ b/source/source_pw/module_pwdft/stress_pw.cpp @@ -4,7 +4,7 @@ #include "source_base/tool_quit.h" #include "source_base/global_variable.h" // use GlobalC #include "source_hamilt/module_vdw/vdw.h" -#include "source_io/module_output/output_log.h" +#include "source_cell/output_log.h" #include "source_hamilt/module_xc/xc_functional.h" #include "source_hamilt/module_xc/general_exx_info.h" // for General_Exx_Info type diff --git a/source/source_pw/module_stodft/sto_forces.cpp b/source/source_pw/module_stodft/sto_forces.cpp index e092b8f9327..2e0d6ef2864 100644 --- a/source/source_pw/module_stodft/sto_forces.cpp +++ b/source/source_pw/module_stodft/sto_forces.cpp @@ -5,7 +5,7 @@ #include "source_estate/elecstate.h" #include "source_estate/module_pot/efield.h" #include "source_estate/module_pot/gatefield.h" -#include "source_io/module_output/output_log.h" +#include "source_cell/output_log.h" #include "source_io/module_parameter/parameter.h" #include "source_pw/module_pwdft/fs_nonlocal_tools.h" diff --git a/source/source_pw/module_stodft/sto_stress_pw.cpp b/source/source_pw/module_stodft/sto_stress_pw.cpp index b85ab58b503..887629e2e63 100644 --- a/source/source_pw/module_stodft/sto_stress_pw.cpp +++ b/source/source_pw/module_stodft/sto_stress_pw.cpp @@ -5,7 +5,7 @@ #include "source_pw/module_pwdft/fs_kin_tools.h" #include "source_pw/module_pwdft/fs_nonlocal_tools.h" #include "source_pw/module_pwdft/stru_fac.h" -#include "source_io/module_output/output_log.h" +#include "source_cell/output_log.h" #include "source_io/module_parameter/parameter.h" template diff --git a/source/source_relax/relax_driver.cpp b/source/source_relax/relax_driver.cpp index 1d2fd132eed..8fbfc6002a2 100644 --- a/source/source_relax/relax_driver.cpp +++ b/source/source_relax/relax_driver.cpp @@ -4,7 +4,7 @@ #include "source_base/version.h" #include "source_cell/cif_io.h" #include "source_io/module_json/output_info.h" -#include "source_io/module_output/output_log.h" +#include "source_cell/output_log.h" #include "source_io/module_output/print_info.h" #include "source_base/module_out/read_exit_file.h" #include "source_io/module_parameter/parameter.h"