From 8c713022b44356cc167da72faf59bbeb1ea51443 Mon Sep 17 00:00:00 2001 From: ubuntu <3158793232@qq.com> Date: Mon, 7 Jul 2025 16:56:34 +0800 Subject: [PATCH 01/12] remove float error in sunway --- CMakeLists.txt | 15 ++++++++ .../module_basis/module_pw/pw_gatherscatter.h | 2 +- .../module_cell/module_symmetry/symm_rho.cpp | 2 +- .../module_elecstate/module_charge/charge.cpp | 7 ++-- .../hamilt_pwdft/VNL_in_pw.cpp | 35 +++++++++++++++---- .../module_hamilt_pw/hamilt_pwdft/forces.cpp | 7 ++-- source/module_io/read_input.cpp | 6 ++-- 7 files changed, 59 insertions(+), 15 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e5afbe4793..3b25f6c3e8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -17,6 +17,7 @@ option(USE_CUDA_MPI "Enable CUDA-aware MPI" OFF) option(USE_CUDA_ON_DCU "Enable CUDA on DCU" OFF) option(USE_ROCM "Enable ROCm" OFF) option(USE_DSP "Enable DSP" OFF) +option(USE_SUNWAY "Enable Sunway Architecture" OFF) option(USE_ABACUS_LIBM "Build libmath from source to speed up" OFF) option(ENABLE_LIBXC "Enable using the LibXC package" OFF) @@ -278,6 +279,18 @@ if (USE_DSP) target_link_libraries(${ABACUS_BIN_NAME} ${MT_HOST_DIR}/hthreads/lib/libhthread_device.a) target_link_libraries(${ABACUS_BIN_NAME} ${MT_HOST_DIR}/hthreads/lib/libhthread_host.a) endif() +if (USE_SUNWAY) + add_compile_definitions(__SUNWAY) + set(SUNWAY ON) + include_directories(${SUNWAY_MATH}/include) + include_directories(${SUNWAY_FFT}/include) + include_directories(${SUNWAY_FFT}/lib/libfftw3.a) + + target_link_libraries(${ABACUS_BIN_NAME} ${SUNWAY_MATH}/libswfft.a) + target_link_libraries(${ABACUS_BIN_NAME} ${SUNWAY_MATH}/libswscalapack.a) + target_link_libraries(${ABACUS_BIN_NAME} ${SUNWAY_MATH}/libswlapack.a) + target_link_libraries(${ABACUS_BIN_NAME} ${SUNWAY_MATH}/libswblas.a) +endif() find_package(Threads REQUIRED) target_link_libraries(${ABACUS_BIN_NAME} Threads::Threads) @@ -452,6 +465,8 @@ if(MKLROOT) if(CMAKE_CXX_COMPILER_ID MATCHES Intel) list(APPEND math_libs ifcore) endif() +elseif(USE_SUNWAY) + # Sunway architecture can only use its own math library else() find_package(FFTW3 REQUIRED) find_package(Lapack REQUIRED) diff --git a/source/module_basis/module_pw/pw_gatherscatter.h b/source/module_basis/module_pw/pw_gatherscatter.h index 97be6e5c23..3488ec2ddb 100644 --- a/source/module_basis/module_pw/pw_gatherscatter.h +++ b/source/module_basis/module_pw/pw_gatherscatter.h @@ -123,7 +123,7 @@ void PW_Basis::gathers_scatterp(std::complex* in, std::complex* out) const outp[iz] = inp[iz]; } } - ModuleBase::timer::tick(this->classname, "gathers_scatterp"); + // ModuleBase::timer::tick(this->classname, "gathers_scatterp"); return; } #ifdef __MPI diff --git a/source/module_cell/module_symmetry/symm_rho.cpp b/source/module_cell/module_symmetry/symm_rho.cpp index 858054df23..3b369e9592 100644 --- a/source/module_cell/module_symmetry/symm_rho.cpp +++ b/source/module_cell/module_symmetry/symm_rho.cpp @@ -285,7 +285,7 @@ void Symmetry::rhog_symmetry(std::complex *rhogtot, //assert(rot_count<=nrotk); }//end if section }//end c_index loop - sum /= rot_count; + if (rot_count!=0) sum/= rot_count; for (int isym = 0; isym < rot_count; ++isym) { rhogtot[ipw_record[isym]] = sum/gphase_record[isym]; diff --git a/source/module_elecstate/module_charge/charge.cpp b/source/module_elecstate/module_charge/charge.cpp index 252b336af6..30e020e420 100644 --- a/source/module_elecstate/module_charge/charge.cpp +++ b/source/module_elecstate/module_charge/charge.cpp @@ -279,8 +279,11 @@ void Charge::atomic_rho(const int spin_number_need, for (int ir = 0; ir < mesh; ++ir) { double r2 = atom->ncpp.r[ir] * atom->ncpp.r[ir]; - rhoatm[ir] = atom->ncpp.rho_at[ir] / ModuleBase::FOUR_PI / r2; - } + if (r2!=0) + { + rhoatm[ir] = atom->ncpp.rho_at[ir] / ModuleBase::FOUR_PI / r2; + } + } rhoatm[0] = pow((rhoatm[2] / rhoatm[1]), atom->ncpp.r[1] / (atom->ncpp.r[2] - atom->ncpp.r[1])); // zws add, sunliang updated 2024-03-04 if (rhoatm[0] < 1e-12) diff --git a/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp index ee9b2bb659..751a5c25e8 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp @@ -724,21 +724,42 @@ void pseudopot_cell_vnl::init_vnl(UnitCell& cell, const ModulePW::PW_Basis* rho_ double* jl = new double[kkbeta]; double* aux = new double[kkbeta]; - + #ifdef __SUNWAY + for (int ib = 0; ib < nbeta; ib++) + { + int flag = 0; + for (int ir = 0; ir < kkbeta; ir++) + { + if (std::abs(cell.atoms[it].ncpp.betar(ib, ir)) < 1e-280) + { + std::cout << "ir is " << ir << " " << std::abs(cell.atoms[it].ncpp.betar(ib, ir)) << std::endl; + flag = ir; + break; + } + } + if (flag!=0) + { + for (int ir = flag; ir < kkbeta; ir++) + { + cell.atoms[it].ncpp.betar(ib, ir) = 0; + } + } + } + #endif for (int ib = 0; ib < nbeta; ib++) { const int l = cell.atoms[it].ncpp.lll[ib]; for (int iq = 0; iq < PARAM.globalv.nqx; iq++) { const double q = iq * PARAM.globalv.dq; - ModuleBase::Sphbes::Spherical_Bessel(kkbeta, cell.atoms[it].ncpp.r.data(), q, l, jl); - + ModuleBase::Sphbes::Spherical_Bessel(kkbeta, cell.atoms[it].ncpp.r.data(), q, l, jl); for (int ir = 0; ir < kkbeta; ir++) - { - aux[ir] = cell.atoms[it].ncpp.betar(ib, ir) * jl[ir] * cell.atoms[it].ncpp.r[ir]; + { + aux[ir] = cell.atoms[it].ncpp.betar(ib, ir) * jl[ir] * cell.atoms[it].ncpp.r[ir]; } - double vqint; - ModuleBase::Integral::Simpson_Integral(kkbeta, aux, cell.atoms[it].ncpp.rab.data(), vqint); + + double vqint; + ModuleBase::Integral::Simpson_Integral(kkbeta, aux, cell.atoms[it].ncpp.rab.data(), vqint); this->tab(it, ib, iq) = vqint * pref; } } diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces.cpp index a722bc5449..de9343ff4d 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/forces.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/forces.cpp @@ -478,8 +478,11 @@ void Forces::cal_force_ew(const UnitCell& ucell, #endif for (int ig = 0; ig < rho_basis->npw; ig++) { - aux[ig] *= ModuleBase::libm::exp(-1.0 * rho_basis->gg[ig] * ucell.tpiba2 / alpha / 4.0) - / (rho_basis->gg[ig] * ucell.tpiba2); + if (rho_basis->gg[ig] !=0) + { + aux[ig] *= ModuleBase::libm::exp(-1.0 * rho_basis->gg[ig] * ucell.tpiba2 / alpha / 4.0) + / (rho_basis->gg[ig] * ucell.tpiba2); + } } // set pos rho_basis->ig_gge0 to zero diff --git a/source/module_io/read_input.cpp b/source/module_io/read_input.cpp index 14d6438010..09111824ff 100644 --- a/source/module_io/read_input.cpp +++ b/source/module_io/read_input.cpp @@ -188,19 +188,21 @@ void ReadInput::create_directory(const Parameter& param) } // NOTE: "make_dir_out" must be called by all processes!!! // Maybe it is not good, because only rank 0 can create the directory. + #ifndef __SUNWAY ModuleBase::Global_File::make_dir_out(param.input.suffix, param.input.calculation, out_dir, this->rank, param.input.mdp.md_restart, param.input.out_alllog); // xiaohui add 2013-09-01 - + #endif const std::string ss = "test -d " + PARAM.inp.read_file_dir; + #ifndef __SUNWAY if (system(ss.c_str())) { ModuleBase::WARNING_QUIT("ReadInput", "please set right files directory for reading in."); } - + #endif return; } From 03c1172159bda22340507ac3c2a5c7f15adb048b Mon Sep 17 00:00:00 2001 From: ubuntu <3158793232@qq.com> Date: Mon, 7 Jul 2025 17:03:27 +0800 Subject: [PATCH 02/12] fix ig=0 --- source/module_elecstate/module_pot/H_Hartree_pw.cpp | 6 +++--- source/module_hamilt_general/module_ewald/H_Ewald_pw.cpp | 4 ++-- .../module_hamilt_general/module_surchem/cal_pseudo.cpp | 1 + source/module_hamilt_general/module_surchem/cal_totn.cpp | 6 ++++-- source/module_hamilt_pw/hamilt_pwdft/forces.cpp | 8 ++++++-- source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp | 3 +-- 6 files changed, 17 insertions(+), 11 deletions(-) diff --git a/source/module_elecstate/module_pot/H_Hartree_pw.cpp b/source/module_elecstate/module_pot/H_Hartree_pw.cpp index 898817f736..f3c4ffe70e 100644 --- a/source/module_elecstate/module_pot/H_Hartree_pw.cpp +++ b/source/module_elecstate/module_pot/H_Hartree_pw.cpp @@ -44,18 +44,18 @@ ModuleBase::matrix H_Hartree_pw::v_hartree(const UnitCell &cell, double ehart = 0.0; std::vector> vh_g(rho_basis->npw); + const int ig0 = rho_basis->ig_gge0; #ifdef _OPENMP #pragma omp parallel for reduction(+:ehart) #endif for (int ig = 0; ig < rho_basis->npw; ig++) { - if (rho_basis->gg[ig] >= 1.0e-8) // LiuXh 20180410 - { + if (ig == ig0) continue; // skip G=0 const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI / (cell.tpiba2 * rho_basis->gg[ig]); ehart += (conj(Porter[ig]) * Porter[ig]).real() * fac; vh_g[ig] = fac * Porter[ig]; - } + } Parallel_Reduce::reduce_pool(ehart); diff --git a/source/module_hamilt_general/module_ewald/H_Ewald_pw.cpp b/source/module_hamilt_general/module_ewald/H_Ewald_pw.cpp index 2c8cc60f9d..fea3368a06 100644 --- a/source/module_hamilt_general/module_ewald/H_Ewald_pw.cpp +++ b/source/module_hamilt_general/module_ewald/H_Ewald_pw.cpp @@ -124,10 +124,10 @@ double H_Ewald_pw::compute_ewald(const UnitCell& cell, fact = 1.0; //GlobalV::ofs_running << "\n pwb.gstart = " << pwb.gstart << std::endl; - + const int ig0 = rho_basis->ig_gge0; for (int ig = 0; ig < rho_basis->npw; ig++) { - if(ig == rho_basis->ig_gge0) + if(ig == ig0) { continue; } diff --git a/source/module_hamilt_general/module_surchem/cal_pseudo.cpp b/source/module_hamilt_general/module_surchem/cal_pseudo.cpp index 7ff5bb507c..81ef7eb6fc 100644 --- a/source/module_hamilt_general/module_surchem/cal_pseudo.cpp +++ b/source/module_hamilt_general/module_surchem/cal_pseudo.cpp @@ -10,6 +10,7 @@ void surchem::gauss_charge(const UnitCell& cell, Structure_Factor* sf) { sf->setup_structure_factor(&cell, pgrid, rho_basis); // call strucFac(ntype,ngmc) + const int ig0 = rho_basis->ig_gge0; // G=0 index for (int it = 0; it < cell.ntype; it++) { double RCS = GetAtom.atom_RCS[cell.atoms[it].ncpp.psd]; diff --git a/source/module_hamilt_general/module_surchem/cal_totn.cpp b/source/module_hamilt_general/module_surchem/cal_totn.cpp index 3650a67dfa..34b4a4da39 100644 --- a/source/module_hamilt_general/module_surchem/cal_totn.cpp +++ b/source/module_hamilt_general/module_surchem/cal_totn.cpp @@ -12,8 +12,9 @@ void surchem::cal_totn(const UnitCell& cell, ModuleBase::GlobalFunc::ZEROS(vloc_g, rho_basis->npw); rho_basis->real2recip(vlocal, vloc_g); // now n is vloc in Recispace + const int ig0 = rho_basis->ig_gge0; // G=0 index for (int ig = 0; ig < rho_basis->npw; ig++) { - if(ig==rho_basis->ig_gge0) + if(ig==ig0) { N[ig] = Porter_g[ig]; continue; @@ -39,9 +40,10 @@ void surchem::induced_charge(const UnitCell& cell, const ModulePW::PW_Basis* rho std::complex *induced_rhog = new complex[rho_basis->npw]; ModuleBase::GlobalFunc::ZEROS(induced_rhog, rho_basis->npw); rho_basis->real2recip(delta_phi, delta_phig); + const int ig0 = rho_basis->ig_gge0; // G=0 index for (int ig = 0; ig < rho_basis->npw; ig++) { - if(rho_basis->ig_gge0 == ig) + if(ig==ig0) { continue; } diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces.cpp index de9343ff4d..89d06af65b 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/forces.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/forces.cpp @@ -472,13 +472,17 @@ void Forces::cal_force_ew(const UnitCell& ucell, upperbound = 2.0 * charge * charge * sqrt(2.0 * alpha / ModuleBase::TWO_PI) * erfc(sqrt(ucell.tpiba2 * rho_basis->ggecut / 4.0 / alpha)); } while (upperbound > 1.0e-6); - + const int ig0 = rho_basis->ig_gge0; #ifdef _OPENMP #pragma omp parallel for #endif for (int ig = 0; ig < rho_basis->npw; ig++) { - if (rho_basis->gg[ig] !=0) + if (ig== ig0) + { + continue; // skip G=0 + } + else { aux[ig] *= ModuleBase::libm::exp(-1.0 * rho_basis->gg[ig] * ucell.tpiba2 / alpha / 4.0) / (rho_basis->gg[ig] * ucell.tpiba2); diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp index 32d5c2339d..c87231149c 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp @@ -159,8 +159,7 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, for(int ig = 0;ig< rho_basis->npw;ig++) { const FPTYPE norm_g = sqrt(rho_basis->gg[ig]); - if(norm_g < 1e-4) { continue; -} + if(norm_g < 1e-4) { continue;} for (int l = 0; l < 3; l++) { for (int m = 0;m< 3;m++) From cf96e9398433bcb6192f0d580ce0ae081ab1e7ef Mon Sep 17 00:00:00 2001 From: ubuntu <3158793232@qq.com> Date: Mon, 7 Jul 2025 17:30:12 +0800 Subject: [PATCH 03/12] add the sw --- CMakeLists.txt | 26 ++++++++++----------- source/source_io/read_input.cpp | 4 ++-- source/source_pw/module_pwdft/VNL_in_pw.cpp | 2 +- 3 files changed, 16 insertions(+), 16 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 4505306534..97f5d06203 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -17,7 +17,7 @@ option(USE_CUDA_MPI "Enable CUDA-aware MPI" OFF) option(USE_CUDA_ON_DCU "Enable CUDA on DCU" OFF) option(USE_ROCM "Enable ROCm" OFF) option(USE_DSP "Enable DSP" OFF) -option(USE_SUNWAY "Enable Sunway Architecture" OFF) +option(USE_SW "Enable SW Architecture" OFF) option(USE_ABACUS_LIBM "Build libmath from source to speed up" OFF) option(ENABLE_LIBXC "Enable using the LibXC package" OFF) @@ -282,17 +282,17 @@ if (USE_DSP) target_link_libraries(${ABACUS_BIN_NAME} ${MT_HOST_DIR}/hthreads/lib/libhthread_device.a) target_link_libraries(${ABACUS_BIN_NAME} ${MT_HOST_DIR}/hthreads/lib/libhthread_host.a) endif() -if (USE_SUNWAY) - add_compile_definitions(__SUNWAY) - set(SUNWAY ON) - include_directories(${SUNWAY_MATH}/include) - include_directories(${SUNWAY_FFT}/include) - include_directories(${SUNWAY_FFT}/lib/libfftw3.a) +if (USE_SW) + add_compile_definitions(__SW) + set(SW ON) + include_directories(${SW_MATH}/include) + include_directories(${SW_FFT}/include) + include_directories(${SW_FFT}/lib/libfftw3.a) - target_link_libraries(${ABACUS_BIN_NAME} ${SUNWAY_MATH}/libswfft.a) - target_link_libraries(${ABACUS_BIN_NAME} ${SUNWAY_MATH}/libswscalapack.a) - target_link_libraries(${ABACUS_BIN_NAME} ${SUNWAY_MATH}/libswlapack.a) - target_link_libraries(${ABACUS_BIN_NAME} ${SUNWAY_MATH}/libswblas.a) + target_link_libraries(${ABACUS_BIN_NAME} ${SW_MATH}/libswfft.a) + target_link_libraries(${ABACUS_BIN_NAME} ${SW_MATH}/libswscalapack.a) + target_link_libraries(${ABACUS_BIN_NAME} ${SW_MATH}/libswlapack.a) + target_link_libraries(${ABACUS_BIN_NAME} ${SW_MATH}/libswblas.a) endif() find_package(Threads REQUIRED) @@ -472,8 +472,8 @@ if(MKLROOT) if(CMAKE_CXX_COMPILER_ID MATCHES Intel) list(APPEND math_libs ifcore) endif() -elseif(USE_SUNWAY) - # Sunway architecture can only use its own math library +elseif(USE_SW) + # SW architecture can only use its own math library else() find_package(FFTW3 REQUIRED) find_package(Lapack REQUIRED) diff --git a/source/source_io/read_input.cpp b/source/source_io/read_input.cpp index 17502e3f86..d08e1ac4de 100644 --- a/source/source_io/read_input.cpp +++ b/source/source_io/read_input.cpp @@ -191,7 +191,7 @@ void ReadInput::create_directory(const Parameter& param) } // NOTE: "make_dir_out" must be called by all processes!!! // Maybe it is not good, because only rank 0 can create the directory. - #ifndef __SUNWAY + #ifndef __SW ModuleBase::Global_File::make_dir_out(param.input.suffix, param.input.calculation, out_dir, @@ -200,7 +200,7 @@ void ReadInput::create_directory(const Parameter& param) param.input.out_alllog); // xiaohui add 2013-09-01 #endif const std::string ss = "test -d " + PARAM.inp.read_file_dir; - #ifndef __SUNWAY + #ifndef __SW if (system(ss.c_str())) { ModuleBase::WARNING_QUIT("ReadInput", "please set right files directory for reading in."); diff --git a/source/source_pw/module_pwdft/VNL_in_pw.cpp b/source/source_pw/module_pwdft/VNL_in_pw.cpp index e874634e3b..02911c5918 100644 --- a/source/source_pw/module_pwdft/VNL_in_pw.cpp +++ b/source/source_pw/module_pwdft/VNL_in_pw.cpp @@ -724,7 +724,7 @@ void pseudopot_cell_vnl::init_vnl(UnitCell& cell, const ModulePW::PW_Basis* rho_ double* jl = new double[kkbeta]; double* aux = new double[kkbeta]; - #ifdef __SUNWAY + #ifdef __SW for (int ib = 0; ib < nbeta; ib++) { int flag = 0; From c65408b4a78cb7367101af347c3fc644930d9922 Mon Sep 17 00:00:00 2001 From: liiutao <3158793232@qq.com> Date: Thu, 10 Jul 2025 11:45:45 +0800 Subject: [PATCH 04/12] change the make_dir --- CMakeLists.txt | 5 +++-- source/source_base/global_function.cpp | 2 ++ 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 97f5d06203..277f1924ec 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -287,8 +287,8 @@ if (USE_SW) set(SW ON) include_directories(${SW_MATH}/include) include_directories(${SW_FFT}/include) - include_directories(${SW_FFT}/lib/libfftw3.a) - + + target_link_libraries(${ABACUS_BIN_NAME} ${SW_FFT}/lib/libfftw3.a) target_link_libraries(${ABACUS_BIN_NAME} ${SW_MATH}/libswfft.a) target_link_libraries(${ABACUS_BIN_NAME} ${SW_MATH}/libswscalapack.a) target_link_libraries(${ABACUS_BIN_NAME} ${SW_MATH}/libswlapack.a) @@ -473,6 +473,7 @@ if(MKLROOT) list(APPEND math_libs ifcore) endif() elseif(USE_SW) + list(APPEND math_libs gfortran) # SW architecture can only use its own math library else() find_package(FFTW3 REQUIRED) diff --git a/source/source_base/global_function.cpp b/source/source_base/global_function.cpp index b51315e158..6124714bb1 100644 --- a/source/source_base/global_function.cpp +++ b/source/source_base/global_function.cpp @@ -61,6 +61,7 @@ void OUT(std::ofstream &ofs, const std::string &name) void MAKE_DIR(const std::string &fn) { // ModuleBase::TITLE("global_function","MAKE_DIR"); + #ifndef __SW if (GlobalV::MY_RANK == 0) { std::stringstream ss; @@ -73,6 +74,7 @@ void MAKE_DIR(const std::string &fn) ModuleBase::WARNING_QUIT("MAKE_DIR", fn); } } + #endif return; } From 89d8344ec2ab9a3650f9fe5528459089019ef577 Mon Sep 17 00:00:00 2001 From: liiutao <3158793232@qq.com> Date: Thu, 10 Jul 2025 17:44:42 +0800 Subject: [PATCH 05/12] unify the gg use --- source/source_basis/module_pw/pw_basis.cpp | 7 ++++ .../module_charge/charge_mixing_residual.cpp | 38 ++++++++++--------- .../source_estate/module_pot/H_Hartree_pw.cpp | 12 +++--- .../source_hamilt/module_surchem/cal_totn.cpp | 7 +--- .../module_surchem/sol_force.cpp | 15 +++----- source/source_pw/module_ofdft/kedf_wt.cpp | 19 +++++----- source/source_pw/module_pwdft/forces.cpp | 7 +--- .../source_pw/module_pwdft/stress_func_cc.cpp | 7 +++- .../module_pwdft/stress_func_ewa.cpp | 1 - .../module_pwdft/stress_func_har.cpp | 7 ++-- 10 files changed, 61 insertions(+), 59 deletions(-) diff --git a/source/source_basis/module_pw/pw_basis.cpp b/source/source_basis/module_pw/pw_basis.cpp index 6a27cda397..a7f2b16ff7 100644 --- a/source/source_basis/module_pw/pw_basis.cpp +++ b/source/source_basis/module_pw/pw_basis.cpp @@ -144,6 +144,7 @@ void PW_Basis::collect_local_pw() delete[] this->gcar; this->gcar = new ModuleBase::Vector3[this->npw]; ModuleBase::Vector3 f; + int gamma_num = 0; for(int ig = 0 ; ig < this-> npw ; ++ig) { int isz = this->ig2isz[ig]; @@ -173,6 +174,12 @@ void PW_Basis::collect_local_pw() if(this->gg[ig] < 1e-8) { this->ig_gge0 = ig; + ++gamma_num; + if (gamma_num > 1) + { + ModuleBase::WARNING_QUIT("PW_Basis::collect_local_pw", + "More than one gamma point found in the plane wave basis set.\n"); + } } } return; diff --git a/source/source_estate/module_charge/charge_mixing_residual.cpp b/source/source_estate/module_charge/charge_mixing_residual.cpp index 819dcbea92..4cc38b2a88 100644 --- a/source/source_estate/module_charge/charge_mixing_residual.cpp +++ b/source/source_estate/module_charge/charge_mixing_residual.cpp @@ -129,15 +129,13 @@ double Charge_Mixing::inner_product_recip_rho(std::complex* rho1, std::c auto part_of_noncolin = [&]() { double sum = 0.0; + const int ig0 = this->rhopw->ig_gge0; #ifdef _OPENMP #pragma omp parallel for reduction(+ : sum) #endif for (int ig = 0; ig < this->rhopw->npw; ++ig) { - if (this->rhopw->gg[ig] < 1e-8) - { - continue; - } + if (ig == ig0) {continue;} sum += (conj(rhog1[0][ig]) * rhog2[0][ig]).real() / this->rhopw->gg[ig]; } sum *= fac; @@ -152,14 +150,13 @@ double Charge_Mixing::inner_product_recip_rho(std::complex* rho1, std::c case 2: { // (1) First part of density error. + const int ig0 = this->rhopw->ig_gge0; #ifdef _OPENMP #pragma omp parallel for reduction(+ : sum) #endif for (int ig = 0; ig < this->rhopw->npw; ++ig) { - if (this->rhopw->gg[ig] < 1e-8) { - continue; -} + if (ig == ig0) {continue;} sum += (conj(rhog1[0][ig] + rhog1[1][ig]) * (rhog2[0][ig] + rhog2[1][ig])).real() / this->rhopw->gg[ig]; } sum *= fac; @@ -203,12 +200,13 @@ double Charge_Mixing::inner_product_recip_rho(std::complex* rho1, std::c } else { // another part with magnetization + const int ig0 = this->rhopw->ig_gge0; #ifdef _OPENMP #pragma omp parallel for reduction(+ : sum) #endif for (int ig = 0; ig < this->rhopw->npw; ig++) { - if (ig == this->rhopw->ig_gge0) + if (ig == ig0) { continue; } @@ -298,14 +296,16 @@ double Charge_Mixing::inner_product_recip_hartree(std::complex* rhog1, s auto part_of_rho = [&]() { double sum = 0.0; + const int ig0 = this->rhopw->ig_gge0; #ifdef _OPENMP #pragma omp parallel for reduction(+ : sum) #endif for (int ig = 0; ig < this->rhopw->npw; ++ig) { - if (this->rhopw->gg[ig] < 1e-8) { + if (ig == ig0) + { continue; -} + } sum += (conj(rhog1[ig]) * rhog2[ig]).real() / this->rhopw->gg[ig]; } sum *= fac; @@ -319,14 +319,16 @@ double Charge_Mixing::inner_product_recip_hartree(std::complex* rhog1, s else if (PARAM.inp.nspin==2) { // charge density part + const int ig0 = this->rhopw->ig_gge0; #ifdef _OPENMP #pragma omp parallel for reduction(+ : sum) #endif for (int ig = 0; ig < this->rhopw->npw; ++ig) { - if (this->rhopw->gg[ig] < 1e-8) { + if (ig == ig0) + { continue; -} + } sum += (conj(rhog1[ig]) * (rhog2[ig])).real() / this->rhopw->gg[ig]; } sum *= fac; @@ -369,18 +371,16 @@ double Charge_Mixing::inner_product_recip_hartree(std::complex* rhog1, s else if (this->mixing_angle <= 0) { // sum for tradtional mixing + const int ig0 = this->rhopw->ig_gge0; #ifdef _OPENMP #pragma omp parallel for reduction(+ : sum) #endif for (int ig = 0; ig < this->rhopw->npw; ig++) { - if (ig == this->rhopw->ig_gge0) { - continue; -} + if (ig == ig0) {continue;} sum += (conj(rhog1[ig]) * rhog2[ig]).real() / this->rhopw->gg[ig]; } sum *= fac; - const int ig0 = this->rhopw->ig_gge0; if (ig0 > 0) { sum += fac2 @@ -408,14 +408,16 @@ double Charge_Mixing::inner_product_recip_hartree(std::complex* rhog1, s else if (this->mixing_angle > 0) { // sum for angle mixing + const int ig0 = this->rhopw->ig_gge0; #ifdef _OPENMP #pragma omp parallel for reduction(+ : sum) #endif for (int ig = 0; ig < this->rhopw->npw; ig++) { - if (ig == this->rhopw->ig_gge0) { + if (ig == ig0) + { continue; -} + } sum += (conj(rhog1[ig]) * rhog2[ig]).real() / this->rhopw->gg[ig]; } sum *= fac; diff --git a/source/source_estate/module_pot/H_Hartree_pw.cpp b/source/source_estate/module_pot/H_Hartree_pw.cpp index 30f8d80484..a9cecd0b64 100644 --- a/source/source_estate/module_pot/H_Hartree_pw.cpp +++ b/source/source_estate/module_pot/H_Hartree_pw.cpp @@ -50,11 +50,13 @@ ModuleBase::matrix H_Hartree_pw::v_hartree(const UnitCell &cell, #endif for (int ig = 0; ig < rho_basis->npw; ig++) { - if (ig == ig0) continue; // skip G=0 - const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI / (cell.tpiba2 * rho_basis->gg[ig]); - - ehart += (conj(Porter[ig]) * Porter[ig]).real() * fac; - vh_g[ig] = fac * Porter[ig]; + if (ig == ig0) + { + continue; // skip G=0 + } + const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI / (cell.tpiba2 * rho_basis->gg[ig]); + ehart += (conj(Porter[ig]) * Porter[ig]).real() * fac; + vh_g[ig] = fac * Porter[ig]; } diff --git a/source/source_hamilt/module_surchem/cal_totn.cpp b/source/source_hamilt/module_surchem/cal_totn.cpp index 1e014c33f4..501046cb28 100644 --- a/source/source_hamilt/module_surchem/cal_totn.cpp +++ b/source/source_hamilt/module_surchem/cal_totn.cpp @@ -47,11 +47,8 @@ void surchem::induced_charge(const UnitCell& cell, const ModulePW::PW_Basis* rho { continue; } - else - { - const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI /(cell.tpiba2 * rho_basis->gg[ig]); - induced_rhog[ig] = -delta_phig[ig] / fac; - } + const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI /(cell.tpiba2 * rho_basis->gg[ig]); + induced_rhog[ig] = -delta_phig[ig] / fac; } rho_basis->recip2real(induced_rhog, induced_rho); diff --git a/source/source_hamilt/module_surchem/sol_force.cpp b/source/source_hamilt/module_surchem/sol_force.cpp index 89db985468..8ce4220b9e 100644 --- a/source/source_hamilt/module_surchem/sol_force.cpp +++ b/source/source_hamilt/module_surchem/sol_force.cpp @@ -21,7 +21,7 @@ void surchem::force_cor_one(const UnitCell& cell, // double Ael1 = 0; // ModuleBase::GlobalFunc::ZEROS(vg, ngmc); int iat = 0; - + const int ig0 = rho_basis->ig_gge0; for (int it = 0;it < cell.ntype;it++) { for (int ia = 0;ia < cell.atoms[it].na ; ia++) @@ -31,17 +31,14 @@ void surchem::force_cor_one(const UnitCell& cell, std::complex phase = exp( ModuleBase::NEG_IMAG_UNIT *ModuleBase::TWO_PI * ( rho_basis->gcar[ig] * cell.atoms[it].tau[ia])); //vloc for each atom vloc_at[ig] = vloc(it, rho_basis->ig2igg[ig]) * phase; - if(rho_basis->ig_gge0 == ig) + if(ig==ig0) { N[ig] = cell.atoms[it].ncpp.zv / cell.omega; + continue; // skip G=0 } - else - { - const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI / - (cell.tpiba2 * rho_basis->gg[ig]); - - N[ig] = -vloc_at[ig] / fac; - } + const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI / + (cell.tpiba2 * rho_basis->gg[ig]); + N[ig] = -vloc_at[ig] / fac; //force for each atom forcesol(iat, 0) += rho_basis->gcar[ig][0] * imag(conj(delta_phi_g[ig]) * N[ig]); diff --git a/source/source_pw/module_ofdft/kedf_wt.cpp b/source/source_pw/module_ofdft/kedf_wt.cpp index 460c515d3d..437a199593 100644 --- a/source/source_pw/module_ofdft/kedf_wt.cpp +++ b/source/source_pw/module_ofdft/kedf_wt.cpp @@ -292,6 +292,7 @@ void KEDF_WT::get_stress(const double* const* prho, ModulePW::PW_Basis* pw_rho, double eta = 0.; double diff = 0.; this->stress.zero_out(); + const int ig0 = pw_rho->ig_gge0; for (int is = 0; is < PARAM.inp.nspin; ++is) { for (int ip = 0; ip < pw_rho->npw; ++ip) @@ -299,22 +300,20 @@ void KEDF_WT::get_stress(const double* const* prho, ModulePW::PW_Basis* pw_rho, eta = sqrt(pw_rho->gg[ip]) * pw_rho->tpiba / this->tkf_; diff = this->diff_linhard(eta, vw_weight); diff *= eta * (recipRhoAlpha[is][ip] * std::conj(recipRhoBeta[is][ip])).real(); - // cout << "diff " << diff << endl; + if (ip == ig0) + { + continue; + } for (int a = 0; a < 3; ++a) { for (int b = a; b < 3; ++b) { - if (pw_rho->gg[ip] == 0.) - { - continue; - } - else + this->stress(a, b) += -diff * pw_rho->gcar[ip][a] * pw_rho->gcar[ip][b] / pw_rho->gg[ip]; + if (a == b) { - this->stress(a, b) += -diff * pw_rho->gcar[ip][a] * pw_rho->gcar[ip][b] / pw_rho->gg[ip]; - if (a == b) { - this->stress(a, b) += diff * coef; -} + this->stress(a, b) += diff * coef; } + } } } diff --git a/source/source_pw/module_pwdft/forces.cpp b/source/source_pw/module_pwdft/forces.cpp index 155ffa0484..6aab09161f 100644 --- a/source/source_pw/module_pwdft/forces.cpp +++ b/source/source_pw/module_pwdft/forces.cpp @@ -482,11 +482,8 @@ void Forces::cal_force_ew(const UnitCell& ucell, { continue; // skip G=0 } - else - { - aux[ig] *= ModuleBase::libm::exp(-1.0 * rho_basis->gg[ig] * ucell.tpiba2 / alpha / 4.0) - / (rho_basis->gg[ig] * ucell.tpiba2); - } + aux[ig] *= ModuleBase::libm::exp(-1.0 * rho_basis->gg[ig] * ucell.tpiba2 / alpha / 4.0) + / (rho_basis->gg[ig] * ucell.tpiba2); } // set pos rho_basis->ig_gge0 to zero diff --git a/source/source_pw/module_pwdft/stress_func_cc.cpp b/source/source_pw/module_pwdft/stress_func_cc.cpp index b58eb41295..500be82e2f 100644 --- a/source/source_pw/module_pwdft/stress_func_cc.cpp +++ b/source/source_pw/module_pwdft/stress_func_cc.cpp @@ -148,6 +148,7 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, rho_basis, 0); // non diagonal term (g=0 contribution missing) + const int ig0 = rho_basis->ig_gge0; #ifdef _OPENMP #pragma omp parallel { @@ -158,8 +159,10 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, #endif for(int ig = 0;ig< rho_basis->npw;ig++) { - const FPTYPE norm_g = sqrt(rho_basis->gg[ig]); - if(norm_g < 1e-4) { continue;} + if (ig == ig0) + { + continue; // skip G=0 + } for (int l = 0; l < 3; l++) { for (int m = 0;m< 3;m++) diff --git a/source/source_pw/module_pwdft/stress_func_ewa.cpp b/source/source_pw/module_pwdft/stress_func_ewa.cpp index c7c8eb1a3e..c3297cfd72 100644 --- a/source/source_pw/module_pwdft/stress_func_ewa.cpp +++ b/source/source_pw/module_pwdft/stress_func_ewa.cpp @@ -95,7 +95,6 @@ void Stress_Func::stress_ewa(const UnitCell& ucell, { continue; } - g2 = rho_basis->gg[ig]* ucell.tpiba2; g2a = g2 /4.0/alpha; rhostar=std::complex(0.0,0.0); diff --git a/source/source_pw/module_pwdft/stress_func_har.cpp b/source/source_pw/module_pwdft/stress_func_har.cpp index d4902dcb36..9478ee6575 100644 --- a/source/source_pw/module_pwdft/stress_func_har.cpp +++ b/source/source_pw/module_pwdft/stress_func_har.cpp @@ -57,7 +57,7 @@ void Stress_Func::stress_har(const UnitCell& ucell, //============================= rho_basis->real2recip(aux, aux); - + const int ig0 = rho_basis->ig_gge0; #ifndef _OPENMP ModuleBase::matrix& local_sigma = sigma; #else @@ -68,9 +68,8 @@ void Stress_Func::stress_har(const UnitCell& ucell, #endif for (int ig = 0 ; ig < rho_basis->npw ; ++ig) { - const FPTYPE g2 = rho_basis->gg[ig]; - if(g2 < 1e-8) - { + if (ig == ig0) + { continue; } From 2eab680d0b42e32cafa9b29a2a4df4ebdd8c08f6 Mon Sep 17 00:00:00 2001 From: liiutao <3158793232@qq.com> Date: Thu, 10 Jul 2025 18:01:10 +0800 Subject: [PATCH 06/12] fix compile bug --- source/source_estate/module_charge/charge_mixing_residual.cpp | 2 -- source/source_pw/module_pwdft/stress_func_cc.cpp | 1 + source/source_pw/module_pwdft/stress_func_har.cpp | 2 +- 3 files changed, 2 insertions(+), 3 deletions(-) diff --git a/source/source_estate/module_charge/charge_mixing_residual.cpp b/source/source_estate/module_charge/charge_mixing_residual.cpp index 4cc38b2a88..9e9c5e0131 100644 --- a/source/source_estate/module_charge/charge_mixing_residual.cpp +++ b/source/source_estate/module_charge/charge_mixing_residual.cpp @@ -213,7 +213,6 @@ double Charge_Mixing::inner_product_recip_rho(std::complex* rho1, std::c sum += (conj(rhog1[0][ig]) * rhog2[0][ig]).real() / this->rhopw->gg[ig]; } sum *= fac; - const int ig0 = this->rhopw->ig_gge0; if (ig0 > 0) { sum += fac2 @@ -421,7 +420,6 @@ double Charge_Mixing::inner_product_recip_hartree(std::complex* rhog1, s sum += (conj(rhog1[ig]) * rhog2[ig]).real() / this->rhopw->gg[ig]; } sum *= fac; - const int ig0 = this->rhopw->ig_gge0; if (ig0 > 0) { sum += fac2 diff --git a/source/source_pw/module_pwdft/stress_func_cc.cpp b/source/source_pw/module_pwdft/stress_func_cc.cpp index 500be82e2f..faa2de25b7 100644 --- a/source/source_pw/module_pwdft/stress_func_cc.cpp +++ b/source/source_pw/module_pwdft/stress_func_cc.cpp @@ -167,6 +167,7 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, { for (int m = 0;m< 3;m++) { + const FPTYPE norm_g = sqrt(rho_basis->gg[ig]); const std::complex t = conj(psic[ig]) * p_sf->strucFac(nt, ig) * rhocg[rho_basis->ig2igg[ig]] * ucell.tpiba * rho_basis->gcar[ig][l] * rho_basis->gcar[ig][m] / norm_g * fact; diff --git a/source/source_pw/module_pwdft/stress_func_har.cpp b/source/source_pw/module_pwdft/stress_func_har.cpp index 9478ee6575..b99ee7fb04 100644 --- a/source/source_pw/module_pwdft/stress_func_har.cpp +++ b/source/source_pw/module_pwdft/stress_func_har.cpp @@ -72,7 +72,7 @@ void Stress_Func::stress_har(const UnitCell& ucell, { continue; } - + const FPTYPE g2 = rho_basis->gg[ig]; FPTYPE shart= ( conj( aux[ig] ) * aux[ig] ).real()/(ucell.tpiba2 * g2); for(int l=0;l<3;l++) { From 4c4a87a112caf52c2fe01b43846cad3cf2746c7b Mon Sep 17 00:00:00 2001 From: liiutao <3158793232@qq.com> Date: Thu, 10 Jul 2025 22:13:20 +0800 Subject: [PATCH 07/12] add init --- source/source_pw/module_pwdft/VNL_in_pw.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/source/source_pw/module_pwdft/VNL_in_pw.cpp b/source/source_pw/module_pwdft/VNL_in_pw.cpp index 02911c5918..dc3752d0c3 100644 --- a/source/source_pw/module_pwdft/VNL_in_pw.cpp +++ b/source/source_pw/module_pwdft/VNL_in_pw.cpp @@ -757,9 +757,8 @@ void pseudopot_cell_vnl::init_vnl(UnitCell& cell, const ModulePW::PW_Basis* rho_ { aux[ir] = cell.atoms[it].ncpp.betar(ib, ir) * jl[ir] * cell.atoms[it].ncpp.r[ir]; } - - double vqint; - ModuleBase::Integral::Simpson_Integral(kkbeta, aux, cell.atoms[it].ncpp.rab.data(), vqint); + double vqint=0.0; + ModuleBase::Integral::Simpson_Integral(kkbeta, aux, cell.atoms[it].ncpp.rab.data(), vqint); this->tab(it, ib, iq) = vqint * pref; } } From 0ca0462fa4a414b080795dee22ea9f5e662b690e Mon Sep 17 00:00:00 2001 From: liiutao <3158793232@qq.com> Date: Fri, 11 Jul 2025 09:42:35 +0800 Subject: [PATCH 08/12] temporarily remove the sunway define --- source/source_pw/module_pwdft/VNL_in_pw.cpp | 22 --------------------- 1 file changed, 22 deletions(-) diff --git a/source/source_pw/module_pwdft/VNL_in_pw.cpp b/source/source_pw/module_pwdft/VNL_in_pw.cpp index dc3752d0c3..e26e09fa3e 100644 --- a/source/source_pw/module_pwdft/VNL_in_pw.cpp +++ b/source/source_pw/module_pwdft/VNL_in_pw.cpp @@ -724,28 +724,6 @@ void pseudopot_cell_vnl::init_vnl(UnitCell& cell, const ModulePW::PW_Basis* rho_ double* jl = new double[kkbeta]; double* aux = new double[kkbeta]; - #ifdef __SW - for (int ib = 0; ib < nbeta; ib++) - { - int flag = 0; - for (int ir = 0; ir < kkbeta; ir++) - { - if (std::abs(cell.atoms[it].ncpp.betar(ib, ir)) < 1e-280) - { - std::cout << "ir is " << ir << " " << std::abs(cell.atoms[it].ncpp.betar(ib, ir)) << std::endl; - flag = ir; - break; - } - } - if (flag!=0) - { - for (int ir = flag; ir < kkbeta; ir++) - { - cell.atoms[it].ncpp.betar(ib, ir) = 0; - } - } - } - #endif for (int ib = 0; ib < nbeta; ib++) { const int l = cell.atoms[it].ncpp.lll[ib]; From e245b26f5483363362b8a5b11a0371ca84f24cb0 Mon Sep 17 00:00:00 2001 From: liiutao <3158793232@qq.com> Date: Fri, 11 Jul 2025 12:49:15 +0800 Subject: [PATCH 09/12] add the pesduo --- source/source_cell/pseudo.cpp | 25 +++++++++++++++++++++++++ source/source_cell/pseudo.h | 6 ++++++ source/source_cell/read_pp_upf100.cpp | 3 ++- source/source_cell/read_pp_upf201.cpp | 3 ++- 4 files changed, 35 insertions(+), 2 deletions(-) diff --git a/source/source_cell/pseudo.cpp b/source/source_cell/pseudo.cpp index e3124bbf4d..fa165968c0 100644 --- a/source/source_cell/pseudo.cpp +++ b/source/source_cell/pseudo.cpp @@ -9,6 +9,31 @@ pseudo::~pseudo() { } +void pseudo::check_pseudo() +{ + for (int ib = 0; ib < nbeta; ib++) + { + int flag = mesh+1; + for (int ir = 0; ir < mesh; ir++) + { + if (std::abs(betar(ib, ir)) < 1.0e-30) + { + flag = ir; + break; + } + } + if (flag r_L ModuleBase::matrix qqq; // qqq(nbeta,nbeta) q_{mu,nu} + /** + * @brief Check the input data for non-normal numbers in the betar. + * Subsequent values following non-normal numbers will be reset to zero + * to prevent potential computational issues arising from invalid data. + */ + void check_betar(); void print_pseudo_h(std::ofstream& ofs); void print_pseudo_atom(std::ofstream& ofs); diff --git a/source/source_cell/read_pp_upf100.cpp b/source/source_cell/read_pp_upf100.cpp index bbac3f08aa..bc093b329e 100644 --- a/source/source_cell/read_pp_upf100.cpp +++ b/source/source_cell/read_pp_upf100.cpp @@ -319,7 +319,8 @@ void Pseudopot_upf::read_pseudo_nl(std::ifstream& ifs, Atom_pseudo& pp) ModuleBase::GlobalFunc::SCAN_END(ifs, ""); pp.kkbeta = (this->kbeta[i] > pp.kkbeta) ? this->kbeta[i] : pp.kkbeta; } - + //check the betar for non-normal number + pp.check_pseudo(); // DIJ ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "", false); ModuleBase::GlobalFunc::READ_VALUE(ifs, this->nd); // nl_4 diff --git a/source/source_cell/read_pp_upf201.cpp b/source/source_cell/read_pp_upf201.cpp index 45351ca3cc..68cc712f5d 100644 --- a/source/source_cell/read_pp_upf201.cpp +++ b/source/source_cell/read_pp_upf201.cpp @@ -530,7 +530,8 @@ void Pseudopot_upf::read_pseudo_upf201_nonlocal(std::ifstream& ifs, Atom_pseudo& word = ""; ModuleBase::GlobalFunc::SCAN_END(ifs, word); } - + //check the betar for non-normal number + pp.check_pseudo(); // Read the hamiltonian terms D_ij if (ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, " Date: Fri, 11 Jul 2025 15:17:04 +0800 Subject: [PATCH 10/12] fix compile bug --- source/source_cell/pseudo.cpp | 6 +++--- source/source_cell/read_pp_upf100.cpp | 2 +- source/source_cell/read_pp_upf201.cpp | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/source/source_cell/pseudo.cpp b/source/source_cell/pseudo.cpp index fa165968c0..39393cae8c 100644 --- a/source/source_cell/pseudo.cpp +++ b/source/source_cell/pseudo.cpp @@ -9,7 +9,7 @@ pseudo::~pseudo() { } -void pseudo::check_pseudo() +void pseudo::check_betar() { for (int ib = 0; ib < nbeta; ib++) { @@ -24,8 +24,8 @@ void pseudo::check_pseudo() } if (flagkbeta[i] > pp.kkbeta) ? this->kbeta[i] : pp.kkbeta; } //check the betar for non-normal number - pp.check_pseudo(); + pp.check_betar(); // DIJ ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "", false); ModuleBase::GlobalFunc::READ_VALUE(ifs, this->nd); // nl_4 diff --git a/source/source_cell/read_pp_upf201.cpp b/source/source_cell/read_pp_upf201.cpp index 68cc712f5d..e0941664bf 100644 --- a/source/source_cell/read_pp_upf201.cpp +++ b/source/source_cell/read_pp_upf201.cpp @@ -531,7 +531,7 @@ void Pseudopot_upf::read_pseudo_upf201_nonlocal(std::ifstream& ifs, Atom_pseudo& ModuleBase::GlobalFunc::SCAN_END(ifs, word); } //check the betar for non-normal number - pp.check_pseudo(); + pp.check_betar(); // Read the hamiltonian terms D_ij if (ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, " Date: Sat, 12 Jul 2025 11:56:19 +0800 Subject: [PATCH 11/12] fix bug in the betar --- source/source_cell/pseudo.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/source/source_cell/pseudo.cpp b/source/source_cell/pseudo.cpp index 39393cae8c..fa267b4e9b 100644 --- a/source/source_cell/pseudo.cpp +++ b/source/source_cell/pseudo.cpp @@ -16,7 +16,8 @@ void pseudo::check_betar() int flag = mesh+1; for (int ir = 0; ir < mesh; ir++) { - if (std::abs(betar(ib, ir)) < 1.0e-30) + // to check is non-normal number and it shoule not be zero + if ((std::abs(betar(ib, ir)) < 1.0e-30) && (std::abs(betar(ib, ir))!=0)) { flag = ir; break; From 97e45a4551c50c11a649029f18416f7ffc701375 Mon Sep 17 00:00:00 2001 From: ubuntu <3158793232@qq.com> Date: Sat, 12 Jul 2025 17:08:42 +0800 Subject: [PATCH 12/12] modify the test --- source/source_cell/pseudo.cpp | 30 ++++++++++++++++-------------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/source/source_cell/pseudo.cpp b/source/source_cell/pseudo.cpp index fa267b4e9b..233632069d 100644 --- a/source/source_cell/pseudo.cpp +++ b/source/source_cell/pseudo.cpp @@ -11,27 +11,29 @@ pseudo::~pseudo() void pseudo::check_betar() { + bool min_flag = false; for (int ib = 0; ib < nbeta; ib++) { - int flag = mesh+1; for (int ir = 0; ir < mesh; ir++) { - // to check is non-normal number and it shoule not be zero - if ((std::abs(betar(ib, ir)) < 1.0e-30) && (std::abs(betar(ib, ir))!=0)) + // Get the bit representation of the double + uint64_t bits = *(uint64_t*)&betar(ib, ir); + // Extract exponent field (bits 52-62) + uint64_t exponent = (bits >> 52) & 0x7FF; + // Define exponent threshold for 1e-30 + // Calculated as: bias + floor(log2(1e-30)) + // Where bias = 1023 and log2(1e-30) ≈ -99.657 + // Thus threshold is approximately 923 + if ((exponent <= 923)) { - flag = ir; - break; + min_flag = true; + betar(ib, ir) = 0.0; } } - if (flag