diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 0f7c22bfe22..428f38fd01d 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -3419,6 +3419,7 @@ - **Availability**: *[`symmetry`](#symmetry)==1 and ([`dft_functional`](#dft_functional) in [hse, hf, pbe0, scan0] or ([`basis_type`](#basis_type)==lcao and [`rpa`](#rpa)==true))* - **Description**: - False: only rotate k-space density matrix D(k) from irreducible k-points to accelerate diagonalization - True: rotate both D(k) and Hexx(R) to accelerate both diagonalization and EXX calculation + For multi-k calculations, D(k) is averaged over the unitary little group of each irreducible k point before star expansion, for either setting. - **Default**: True ### out_ri_cv diff --git a/docs/parameters.yaml b/docs/parameters.yaml index 1c4693afbd7..a0f76236e86 100644 --- a/docs/parameters.yaml +++ b/docs/parameters.yaml @@ -4681,6 +4681,7 @@ parameters: description: | * False: only rotate k-space density matrix D(k) from irreducible k-points to accelerate diagonalization * True: rotate both D(k) and Hexx(R) to accelerate both diagonalization and EXX calculation + For multi-k calculations, D(k) is averaged over the unitary little group of each irreducible k point before star expansion, for either setting. default_value: "True" unit: "" availability: "symmetry==1 and (dft_functional in [hse, hf, pbe0, scan0] or (basis_type==lcao and rpa==true))" diff --git a/source/source_hsolver/diago_elpa_native.cpp b/source/source_hsolver/diago_elpa_native.cpp index a07f8303202..9d30afa590b 100644 --- a/source/source_hsolver/diago_elpa_native.cpp +++ b/source/source_hsolver/diago_elpa_native.cpp @@ -3,6 +3,7 @@ #include "source_base/global_function.h" #include "source_base/module_external/blas_connector.h" #include "source_base/module_external/blacs_connector.h" +#include "source_base/module_external/scalapack_connector.h" #include "source_base/timer.h" #include "source_base/tool_quit.h" #include "source_hsolver/module_genelpa/elpa_new.h" @@ -71,6 +72,79 @@ void DiagoElpaNative::diag_pool(ModuleBase::MatrixBlock& h_mat, std::vector eigen(this->nlocal, 0.0); std::vector eigenvectors(narows * nacols); + // The complex LCAO matrices follow the LAPACK UPLO='U' convention: only + // their upper triangles are guaranteed to contain the Hermitian matrix. + // ELPA's native generalized solver consumes both triangles, so complete + // private Hermitian copies before the solve. + std::vector h_work; + std::vector s_work; + T* h_elpa = h_mat.p; + T* s_elpa = s_mat.p; + int decomposed_state = this->DecomposedState; + if (!std::is_same::value) + { + h_work.resize(narows * nacols); + s_work.resize(narows * nacols); + const int one = 1; + ScalapackConnector::tranc(nFull, + nFull, + T(1.0), + h_mat.p, + one, + one, + h_mat.desc, + T(0.0), + h_work.data(), + one, + one, + h_mat.desc); + ScalapackConnector::tranc(nFull, + nFull, + T(1.0), + s_mat.p, + one, + one, + s_mat.desc, + T(0.0), + s_work.data(), + one, + one, + s_mat.desc); + const auto local_to_global = [](const int local_index, + const int block_size, + const int process_coordinate, + const int source_coordinate, + const int process_count) { + if (source_coordinate < 0) + { + return local_index; + } + const int process_offset + = (process_coordinate - source_coordinate + process_count) % process_count; + return ((local_index / block_size) * process_count + process_offset) * block_size + + local_index % block_size; + }; + for (int local_col = 0; local_col < nacols; ++local_col) + { + const int global_col + = local_to_global(local_col, h_mat.desc[5], mypcol, h_mat.desc[7], npcols); + for (int local_row = 0; local_row < narows; ++local_row) + { + const int global_row + = local_to_global(local_row, h_mat.desc[4], myprow, h_mat.desc[6], nprows); + if (global_row <= global_col) + { + const int local_index = local_row + local_col * narows; + h_work[local_index] = h_mat.p[local_index]; + s_work[local_index] = s_mat.p[local_index]; + } + } + } + h_elpa = h_work.data(); + s_elpa = s_work.data(); + decomposed_state = 0; + } + if (elpa_init(20210430) != ELPA_OK) { fprintf(stderr, "Error: ELPA API version not supported"); @@ -114,11 +188,11 @@ void DiagoElpaNative::diag_pool(ModuleBase::MatrixBlock& h_mat, #endif elpa_generalized_eigenvectors(handle, - h_mat.p, - s_mat.p, + h_elpa, + s_elpa, eigen.data(), eigenvectors.data(), - this->DecomposedState, + decomposed_state, &success); elpa_deallocate(handle, &success); elpa_uninit(&success); diff --git a/source/source_hsolver/kernels/cuda/diag_cusolver.cu b/source/source_hsolver/kernels/cuda/diag_cusolver.cu index 642198cfd92..52e206483a2 100644 --- a/source/source_hsolver/kernels/cuda/diag_cusolver.cu +++ b/source/source_hsolver/kernels/cuda/diag_cusolver.cu @@ -9,7 +9,8 @@ Diag_Cusolver_gvd::Diag_Cusolver_gvd(){ itype = CUSOLVER_EIG_TYPE_1; // A*x = (lambda)*B*x jobz = CUSOLVER_EIG_MODE_VECTOR; // compute eigenvalues and eigenvectors. - uplo = CUBLAS_FILL_MODE_LOWER; + // LCAO supplies the authoritative upper triangle of the Hermitian H/S matrices. + uplo = CUBLAS_FILL_MODE_UPPER; d_A = NULL; d_B = NULL; diff --git a/source/source_hsolver/test/CMakeLists.txt b/source/source_hsolver/test/CMakeLists.txt index 96895a504ed..043c0596c1f 100644 --- a/source/source_hsolver/test/CMakeLists.txt +++ b/source/source_hsolver/test/CMakeLists.txt @@ -93,7 +93,7 @@ if (ENABLE_MPI) AddTest( TARGET MODULE_HSOLVER_LCAO LIBS parameter ELPA::ELPA base genelpa psi device - SOURCES diago_lcao_test.cpp ../diago_elpa.cpp ../diago_scalapack.cpp ../diago_lapack.cpp + SOURCES diago_lcao_test.cpp ../diago_elpa.cpp ../diago_elpa_native.cpp ../diago_scalapack.cpp ../diago_lapack.cpp ) else() AddTest( diff --git a/source/source_hsolver/test/diago_elpa_utils.h b/source/source_hsolver/test/diago_elpa_utils.h index 22e986935d9..843c53b6a51 100644 --- a/source/source_hsolver/test/diago_elpa_utils.h +++ b/source/source_hsolver/test/diago_elpa_utils.h @@ -167,7 +167,7 @@ void lapack_diago(double *hmatrix, double *smatrix, double *e, int &nFull) const char jobz = 'V'; // 'N':only calc eigenvalue, 'V': eigenvalues and eigenvectors const char uplo = 'U'; // Upper triangles int lwork = (nFull + 2) * nFull, info = 0; - double *ev = new double[nFull * nFull]; + double *ev = new double[lwork]; double *a = new double[nFull * nFull]; double *b = new double[nFull * nFull]; @@ -196,7 +196,7 @@ void lapack_diago(std::complex *hmatrix, std::complex *smatrix, const char uplo = 'U'; // Upper triangles int lwork = (nFull + 1) * nFull, info = 0; double *rwork = new double[3 * nFull - 2]; - std::complex *ev = new std::complex[nFull * nFull]; + std::complex *ev = new std::complex[lwork]; std::complex *a = new std::complex[nFull * nFull]; std::complex *b = new std::complex[nFull * nFull]; diff --git a/source/source_hsolver/test/diago_lcao_cusolver_test.cpp b/source/source_hsolver/test/diago_lcao_cusolver_test.cpp index 88170826de7..bab49060153 100644 --- a/source/source_hsolver/test/diago_lcao_cusolver_test.cpp +++ b/source/source_hsolver/test/diago_lcao_cusolver_test.cpp @@ -134,6 +134,21 @@ class DiagoPrepare return ok; } + void poison_lower_triangle() + { + // The distributed solver buffers are column-major. Keep the original + // row-major fixtures intact for the independent LAPACK reference. + for (int col = 0; col < nlocal; ++col) + { + for (int row = col + 1; row < nlocal; ++row) + { + const int index = row + col * nlocal; + this->h_local[index] = T(123.0 + row + col); + this->s_local[index] = T(0.0); + } + } + } + void print_hs() { if (!PRINT_HS) @@ -203,6 +218,10 @@ class DiagoPrepare { this->pb2d(); this->distribute_data(); + if (ks_solver == "cusolver") + { + this->poison_lower_triangle(); + } this->print_hs(); this->set_env(); diff --git a/source/source_hsolver/test/diago_lcao_test.cpp b/source/source_hsolver/test/diago_lcao_test.cpp index ecf91b5afaa..18b4a30fb0c 100644 --- a/source/source_hsolver/test/diago_lcao_test.cpp +++ b/source/source_hsolver/test/diago_lcao_test.cpp @@ -9,6 +9,7 @@ #include #ifdef __ELPA #include "source_hsolver/diago_elpa.h" +#include "source_hsolver/diago_elpa_native.h" #include "source_hsolver/module_genelpa/elpa_solver.h" #endif #include "source_base/module_external/scalapack_connector.h" @@ -136,6 +137,47 @@ class DiagoPrepare return ok; } +#ifdef __ELPA + void diago_native_with_poisoned_lower() + { + this->pb2d(); + this->distribute_data(); + int nprows; + int npcols; + int myprow; + int mypcol; + Cblacs_gridinfo(icontxt, &nprows, &npcols, &myprow, &mypcol); + for (int col = 0; col < hmtest.ncol; ++col) + { + const int global_col = (col / nb2d * npcols + mypcol) * nb2d + col % nb2d; + for (int row = 0; row < hmtest.nrow; ++row) + { + const int global_row = (row / nb2d * nprows + myprow) * nb2d + row % nb2d; + if (global_row > global_col) + { + const int index = row + col * hmtest.nrow; + h_local[index] = T(123.0 + global_row + global_col); + s_local[index] = T(0.0); + } + } + } + hmtest.h_local = h_local; + hmtest.s_local = s_local; + ModuleBase::MatrixBlock h_mat; + ModuleBase::MatrixBlock s_mat; + hmtest.matrix(h_mat, s_mat); + hsolver::DiagoElpaNative solver(nlocal, nbands, false); + solver.diag(h_mat, s_mat, psi, e_solver.data()); + EXPECT_EQ(hmtest.h_local, h_local); + EXPECT_EQ(hmtest.s_local, s_local); + // A second solve must not reuse an in-place decomposition of the + // private overlap copy or alter the caller's matrix storage. + solver.diag(h_mat, s_mat, psi, e_solver.data()); + EXPECT_EQ(hmtest.h_local, h_local); + EXPECT_EQ(hmtest.s_local, s_local); + } +#endif + void print_hs() { if (!PRINT_HS) @@ -371,6 +413,26 @@ INSTANTIATE_TEST_SUITE_P( DiagoPrepare>(0, 0, 32, 0, "scalapack_gvx", "H-KPoints-Si64.dat", "S-KPoints-Si64.dat"))); #ifdef __ELPA +class DiagoElpaNativeUpperTest : public ::testing::TestWithParam +{ +}; + +TEST_P(DiagoElpaNativeUpperTest, PreservesInputsAndIgnoresLowerTriangle) +{ + DiagoPrepare> dp(0, 0, GetParam(), 0, "genelpa", + "H-KPoints-Si2.dat", "S-KPoints-Si2.dat"); + ASSERT_TRUE(dp.produce_HS()); + dp.diago_native_with_poisoned_lower(); + if (dp.myrank == 0) + { + dp.diago_lapack(); + std::stringstream out_info; + EXPECT_TRUE(dp.compare_eigen(out_info)) << out_info.str(); + } +} + +INSTANTIATE_TEST_SUITE_P(BlockSizes, DiagoElpaNativeUpperTest, ::testing::Values(1, 2, 3)); + TEST(DiagoElpaComplexTest, UsesAuthoritativeUpperTriangle) { std::stringstream out_info; diff --git a/source/source_io/module_parameter/read_inp_exx_dftu.cpp b/source/source_io/module_parameter/read_inp_exx_dftu.cpp index a0f6e728604..7eff18f20a8 100644 --- a/source/source_io/module_parameter/read_inp_exx_dftu.cpp +++ b/source/source_io/module_parameter/read_inp_exx_dftu.cpp @@ -525,7 +525,8 @@ void ReadInput::item_exx() item.category = "Exact Exchange (LCAO)"; item.type = "Boolean"; item.description = R"(* False: only rotate k-space density matrix D(k) from irreducible k-points to accelerate diagonalization -* True: rotate both D(k) and Hexx(R) to accelerate both diagonalization and EXX calculation)"; +* True: rotate both D(k) and Hexx(R) to accelerate both diagonalization and EXX calculation +For multi-k calculations, D(k) is averaged over the unitary little group of each irreducible k point before star expansion, for either setting.)"; item.default_value = "True"; item.unit = ""; item.set_availability("symmetry==1 and (dft_functional in [hse, hf, pbe0, scan0] or (basis_type==lcao and rpa==true))"); diff --git a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.cpp b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.cpp index 18eb826888c..a1942918f70 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.cpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.cpp @@ -52,28 +52,37 @@ namespace ModuleSymmetry } this->spin_U_ = spin_U; // keep for restore_HR_nspin4 (real-space EXX H(R) spin mixing) - // 2. calculate the rotation matrix in AO-representation for each ibz_kpoint and symmetry operation: M(k, isym) - auto restrict_kpt = [](const TCdouble& kvec, const double& symm_prec) -> TCdouble - {// in (-0.5, 0.5] - TCdouble kvec_res; - kvec_res.x = fmod(kvec.x + 100.5 - 0.5 * symm_prec, 1) - 0.5 + 0.5 * symm_prec; - kvec_res.y = fmod(kvec.y + 100.5 - 0.5 * symm_prec, 1) - 0.5 + 0.5 * symm_prec; - kvec_res.z = fmod(kvec.z + 100.5 - 0.5 * symm_prec, 1) - 0.5 + 0.5 * symm_prec; - if (std::abs(kvec_res.x) < symm_prec) { kvec_res.x = 0.0; } - if (std::abs(kvec_res.y) < symm_prec) { kvec_res.y = 0.0; } - if (std::abs(kvec_res.z) < symm_prec) { kvec_res.z = 0.0; } - return kvec_res; - }; - int nks_ibz = kv.kstars.size(); // kv.nks = 2 * kv.nks_ibz when nspin=2 - this->Ms_.resize(nks_ibz); - for (int ik_ibz = 0;ik_ibz < nks_ibz;++ik_ibz) + // A k-star contains only one operation per distinct k point. The other + // operations fixing k (modulo a reciprocal lattice vector) must still + // be averaged: a finite-grid SCF density need not respect this little group. + const int nks_ibz = kv.kstars.size(); + this->Ms_.assign(nks_ibz, {}); + this->little_groups_.assign(nks_ibz, {}); + for (int ik_ibz = 0; ik_ibz < nks_ibz; ++ik_ibz) { - // const TCdouble& kvec_d_ibz = restrict_kpt((*kstars[ik_ibz].begin()).second * ucell.symm.kgmatrix[(*kstars[ik_ibz].begin()).first], ucell.symm.epsilon); - for (auto& isym_kvd : kv.kstars[ik_ibz]) { - if (isym_kvd.first < nop_tot) { - this->Ms_[ik_ibz][isym_kvd.first] = this->contruct_2d_rot_mat_ao(ucell.symm, ucell.atoms, ucell.st, kv.kvec_d[ik_ibz], isym_kvd.first, pv, spin_U[isym_kvd.first]); -} -} + std::set needed; + for (const auto& member : kv.kstars[ik_ibz]) + { + const int op = (!this->magnetic_nspin4_ && member.first >= nsym_) + ? member.first - nsym_ : member.first; + needed.insert(op); + } + for (int op = 0; op < nsym_; ++op) + { + const auto delta = kv.kvec_d[ik_ibz] * ucell.symm.kgmatrix[op] - kv.kvec_d[ik_ibz]; + if (std::abs(delta.x - std::round(delta.x)) < this->eps_ + && std::abs(delta.y - std::round(delta.y)) < this->eps_ + && std::abs(delta.z - std::round(delta.z)) < this->eps_) + { + this->little_groups_[ik_ibz].push_back(op); + needed.insert(op); + } + } + for (const int op : needed) + { + this->Ms_[ik_ibz][op] = this->contruct_2d_rot_mat_ao( + ucell.symm, ucell.atoms, ucell.st, kv.kvec_d[ik_ibz], op, pv, spin_U[op]); + } } // output Ms of isym=1 // std::ofstream ofs("Ms_kibz7_sym7.dat"); @@ -109,18 +118,37 @@ namespace ModuleSymmetry { for (int ik_ibz = 0;ik_ibz < nk;++ik_ibz) { + // P_k D = |G_k|^{-1} sum_g M_g^T D M_g^*. This preserves + // Hermiticity and makes restoration independent of the chosen + // star representative; rotating just one arbitrary D does not. + const auto& little_group = this->little_groups_.at(ik_ibz); + assert(!little_group.empty()); + std::vector> projected = dm_k_ibz[ik_ibz + is * nk]; + if (little_group.size() > 1) + { + std::fill(projected.begin(), projected.end(), 0.0); + for (const int op : little_group) + { + const auto rotated = this->rot_matrix_ao( + dm_k_ibz[ik_ibz + is * nk], ik_ibz, little_group.size(), op, pv); + for (size_t i = 0; i < projected.size(); ++i) + { + projected[i] += rotated[i]; + } + } + } for (auto& isym_kvd : kv.kstars[ik_ibz]) { if (isym_kvd.first == 0) { double factor = 1.0 / static_cast(kv.kstars[ik_ibz].size()); std::vector> dm_scaled(pv.get_local_size()); - for (int i = 0;i < pv.get_local_size();++i) { dm_scaled[i] = factor * dm_k_ibz[ik_ibz + is * nk][i]; } + for (int i = 0;i < pv.get_local_size();++i) { dm_scaled[i] = factor * projected[i]; } dm_k_full.push_back(dm_scaled); } else if (isym_kvd.first < nsym_) { //space group operations - dm_k_full.push_back(this->rot_matrix_ao(dm_k_ibz[ik_ibz + is * nk], ik_ibz, kv.kstars[ik_ibz].size(), isym_kvd.first, pv)); + dm_k_full.push_back(this->rot_matrix_ao(projected, ik_ibz, kv.kstars[ik_ibz].size(), isym_kvd.first, pv)); } else { // antiunitary elements: Theta * (spatial operation) @@ -140,12 +168,12 @@ namespace ModuleSymmetry // m=0: gray group: the space-group part of anti-unitary elements are the same of the unitary elements, isym_M < nsym_ // m!=0: Shubnikov group: using different space-group part of anti-unitary elements stored in gmatrix_anti with isym_M >= nsym_ dm_k_full.push_back(this->trs_spin_rotate( - this->rot_matrix_ao(dm_k_ibz[ik_ibz + is * nk], ik_ibz, kv.kstars[ik_ibz].size(), isym_M, pv, false), + this->rot_matrix_ao(projected, ik_ibz, kv.kstars[ik_ibz].size(), isym_M, pv, false), sigma_y, pv, 1.0)); } else { - dm_k_full.push_back(this->rot_matrix_ao(dm_k_ibz[ik_ibz + is * nk], ik_ibz, kv.kstars[ik_ibz].size(), isym_M, pv, true)); + dm_k_full.push_back(this->rot_matrix_ao(projected, ik_ibz, kv.kstars[ik_ibz].size(), isym_M, pv, true)); } } } @@ -492,7 +520,7 @@ namespace ModuleSymmetry const char notrans = 'N'; std::complex alpha(1.0, 0.0); const std::complex beta(0.0, 0.0); - const int nbasis = PARAM.globalv.nlocal; + const int nbasis = pv.get_global_row_size(); const int i1 = 1; if (TRS_conj) { diff --git a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h index 3ff91c7523e..1d27edfef2f 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h +++ b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h @@ -74,6 +74,17 @@ namespace ModuleSymmetry std::vector> trs_spin_rotate(const std::vector>& X, const std::vector>& sigma_y, const Parallel_2D& pv, const double scale) const; + /// Inject synthetic AO rotations for density-restoration regression tests. + void set_density_rotations_for_testing( + const std::vector>>>& rotations, + const std::vector>& little_groups, + const int nrot) + { + this->Ms_ = rotations; + this->little_groups_ = little_groups; + this->nsym_ = nrot; + } + /// calculate Wigner D matrix double wigner_d(const double beta, const int l, const int m1, const int m2) const; std::complex wigner_D(const TCdouble& euler_angle, const int l, const int m1, const int m2, const bool inv) const; @@ -209,6 +220,10 @@ namespace ModuleSymmetry /// size: [nks_ibz][nsym][nbasis*nbasis], only need to calculate once. std::vector>>> Ms_; + /// Unitary operations fixing each IBZ k point modulo reciprocal lattice vectors. + /// Geometry data built with Ms_ in cal_Ms, not an SCF workflow switch. + std::vector> little_groups_; + /// (nspin=4) the SU(2) spin-1/2 rotation U(isym) for each symmetry operation, size [nsym]. /// The spinor AO rotation is T(isym) (x) U(isym); restore_HR_nspin4 uses it to mix the 4 spin /// channels of the real-space EXX H(R). Filled in cal_Ms (identity for nspin<4). diff --git a/source/source_lcao/module_ri/module_exx_symmetry/test/CMakeLists.txt b/source/source_lcao/module_ri/module_exx_symmetry/test/CMakeLists.txt index cf9ad10d1af..ba45e581124 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/test/CMakeLists.txt +++ b/source/source_lcao/module_ri/module_exx_symmetry/test/CMakeLists.txt @@ -4,7 +4,7 @@ abacus_disable_feature_definitions(__ROCM) AddTest( TARGET MODULE_RI_EXX_SYMMETRY_rotation LIBS base device symmetry neighbor parameter - SOURCES symm_rotation_test.cpp ../symm_rotation.cpp ../symm_rot_out.cpp ../irreducible_sector.cpp ../irred_sec_bvk.cpp + SOURCES symm_rotation_test.cpp test_symm_rotation.cpp ../symm_rotation.cpp ../symm_rot_out.cpp ../irreducible_sector.cpp ../irred_sec_bvk.cpp ../../../../source_basis/module_ao/parallel_orbitals.cpp ) \ No newline at end of file diff --git a/source/source_lcao/module_ri/module_exx_symmetry/test/test_symm_rotation.cpp b/source/source_lcao/module_ri/module_exx_symmetry/test/test_symm_rotation.cpp new file mode 100644 index 00000000000..01ec4938bfb --- /dev/null +++ b/source/source_lcao/module_ri/module_exx_symmetry/test/test_symm_rotation.cpp @@ -0,0 +1,161 @@ +#include "../symm_rotation.h" +#include "source_io/module_parameter/parameter.h" +#include "gtest/gtest.h" + +class TestParameters +{ + public: + TestParameters(Parameter& parameters, const int nspin) + : parameters_(parameters), original_nspin_(parameters.inp.nspin) + { + parameters_.input.nspin = nspin; + } + ~TestParameters() { parameters_.input.nspin = original_nspin_; } + + private: + Parameter& parameters_; + const int original_nspin_; +}; + +// K-point generation is outside this test: use explicit stars, but provide +// the virtual symbols needed by the existing lightweight rotation test target. +void ModuleCell::ReciprocalGrid::renew(const int&) +{ + ADD_FAILURE() << "Unexpected k-point generation"; +} +void K_Vectors::renew(const int&) +{ + ADD_FAILURE() << "Unexpected k-point generation"; +} +void K_Vectors::reduce_by_symmetry(const UnitCell&, const ModuleSymmetry::Symmetry&, + bool, std::string&, bool&, const int, std::ofstream&) +{ + ADD_FAILURE() << "Unexpected k-point reduction"; +} + +namespace +{ +using Complex = std::complex; + +// Independent dense product in the stored (transposed density) convention. +std::vector rotate_reference(const std::vector& density, + const std::vector& rotation, + const int n) +{ + std::vector result(n * n, 0.0); + for (int i = 0; i < n; ++i) + { + for (int j = 0; j < n; ++j) + { + for (int a = 0; a < n; ++a) + { + for (int b = 0; b < n; ++b) + { + result[i + j * n] += rotation[a + i * n] * density[a + b * n] + * std::conj(rotation[b + j * n]); + } + } + } + } + return result; +} + +void check_little_group_restoration(const int nspin) +{ + const TestParameters parameters(PARAM, nspin); + const int channels = nspin == 2 ? 2 : 1; + const int n = 4; + Parallel_2D pv; + pv.init(n, n, 1, MPI_COMM_WORLD); + ModuleSymmetry::Symmetry_rotation rotation; + std::vector identity(n * n, 0.0); + std::vector little(n * n, 0.0); + std::vector representative(n * n, 0.0); + std::vector alternate(n * n, 0.0); + const int sign[n] = {1, 1, -1, -1}; + for (int i = 0; i < n; ++i) + { + identity[i + i * n] = 1.0; + little[i + i * n] = sign[i]; + const int row = (i + 1) % n; + representative[row + i * n] = std::polar(1.0, 0.3 * i); + alternate[row + i * n] = double(sign[row]) * representative[row + i * n]; + } + auto local = [&pv, n](const std::vector& dense) { + std::vector result(pv.get_local_size()); + for (int i = 0; i < n; ++i) + { + for (int j = 0; j < n; ++j) + { + if (pv.in_this_processor(i, j)) + { + result[pv.global2local_row(i) + pv.global2local_col(j) * pv.get_row_size()] + = dense[i + j * n]; + } + } + } + return result; + }; + rotation.set_density_rotations_for_testing( + {{{0, local(identity)}, {1, local(little)}, {2, local(representative)}, {3, local(alternate)}}}, + {{0, 1}}, 4); + K_Vectors kv; + kv.set_nkstot(channels); + kv.set_nkstot_nospin(2); + kv.kstars = {{{0, {0.25, 0.0, 0.0}}, {2, {0.0, 0.25, 0.0}}}}; + std::vector> inputs; + std::vector> expected; + for (int spin = 0; spin < channels; ++spin) + { + std::vector density(n * n); + std::vector projected(n * n); + for (int i = 0; i < n; ++i) + { + for (int j = 0; j < n; ++j) + { + const Complex value((spin + 1) * (2.0 + i + j), 0.2 * (i - j)); + density[i + j * n] = value; + // For this C2 little group, averaging removes exactly the odd blocks. + projected[i + j * n] = sign[i] == sign[j] ? 0.5 * value : Complex(0.0); + } + } + inputs.push_back(local(density)); + expected.push_back(local(projected)); + expected.push_back(local(rotate_reference(projected, representative, n))); + } + const auto restored = rotation.restore_dm(kv, inputs, pv); + ASSERT_EQ(restored.size(), expected.size()); + for (size_t k = 0; k < expected.size(); ++k) + { + for (size_t i = 0; i < expected[k].size(); ++i) + { + EXPECT_NEAR(std::abs(restored[k][i] - expected[k][i]), 0.0, 1e-12); + } + } + // A different representative of the same star must give the same density. + kv.kstars = {{{0, {0.25, 0.0, 0.0}}, {3, {0.0, 0.25, 0.0}}}}; + const auto changed_representative = rotation.restore_dm(kv, inputs, pv); + for (size_t k = 0; k < expected.size(); ++k) + { + for (size_t i = 0; i < expected[k].size(); ++i) + { + EXPECT_NEAR(std::abs(changed_representative[k][i] - restored[k][i]), 0.0, 1e-12); + } + } +} +} // namespace + +TEST(SymmetryDensityRestoration, LittleGroupAndStarWeight) +{ + check_little_group_restoration(1); +} + +TEST(SymmetryDensityRestoration, IndependentSpinChannels) +{ + check_little_group_restoration(2); +} + +TEST(SymmetryDensityRestoration, SpinorDensity) +{ + check_little_group_restoration(4); +} diff --git a/tests/08_EXX/08_KP_HSE_symm/README b/tests/08_EXX/08_KP_HSE_symm/README index 88ac4f7c94d..4dc65661511 100644 --- a/tests/08_EXX/08_KP_HSE_symm/README +++ b/tests/08_EXX/08_KP_HSE_symm/README @@ -1 +1,36 @@ -HSE calculation on Si, multiple k-points, cal force and stress, exx real number, symmetry=1 \ No newline at end of file +HSE calculation on Si, multiple k-points, force and stress, real EXX, symmetry=1. + +The density at each irreducible k point must be averaged over its unitary +little group before expanding its star. A star stores one representative per +unique k point; it does not contain all operations fixing that point. On this +coarse finite grid, omitting that average leaves the EXX density outside the +space-group-invariant subspace, so rotating irreducible H(R) blocks can produce +a non-Hermitian Hamiltonian. + +In ABACUS's transposed density storage, the average is + P_k(D) = sum_{g in G_k} M_g^T D M_g^* / |G_k|, +where g k = k modulo a reciprocal lattice vector. P_k is a group projector: +it preserves Hermiticity and removes dependence on the chosen unitary star +representative. The existing star-size factor is applied only after averaging. +The unit tests in module_exx_symmetry/test/test_symm_rotation.cpp exercise this +contract with complex matrices and separate spin channels. + +The reference was regenerated after merging upstream develop 54aaffefb, +including its ELPA fallback correction, on 2026-09-14. It uses the unchanged +case with the default genelpa solver on one MPI rank. Its total energy is +-189.4406826032011395 eV. ScaLAPACK runs with 1 and 4 MPI ranks and a full EXX +contraction run (exx_symmetry_realspace=0, with the same projected density) +agree within 1.4e-9 eV. All four report totalstressref=2143.792554 and +totalforceref=0.000000. INPUT, thresholds and the historical timing reference +are unchanged. This is a fixed-input regression reference, not a claim of +cutoff or outer-loop convergence. + +Verification used GCC, OpenMPI 5.0.10, ELPA 2026.02.001, LibXC 7.0.0 and the +repository-pinned LibRI on Sai DSPRHBM, with OMP_NUM_THREADS=1. The 12 EXX +symmetry tests and, on one rank, both complex ELPA upper-triangle tests and +all three native ELPA block-size tests passed. The complete LCAO solver test +still crashes inside the installed two-stage ELPA real routine; the four-rank +complex ELPA tests and HSE case also fail in that library after disabling +XPMEM to avoid a UCX communication error. These failures remain unresolved, +and four-rank ELPA is not included in the agreement claim. No solver-stage +or production INPUT default change is part of this merge. diff --git a/tests/08_EXX/08_KP_HSE_symm/result.ref b/tests/08_EXX/08_KP_HSE_symm/result.ref index eb0739164ed..c8757c65e8d 100644 --- a/tests/08_EXX/08_KP_HSE_symm/result.ref +++ b/tests/08_EXX/08_KP_HSE_symm/result.ref @@ -1,7 +1,7 @@ -etotref -189.41351672 -etotperatomref -94.70675836 +etotref -189.4406826032011395 +etotperatomref -94.7203413016 totalforceref 0.000000 -totalstressref 2153.17347300 +totalstressref 2143.792554 pointgroupref D_3d spacegroupref O_h nksibzref 3