diff --git a/CMakeLists.txt b/CMakeLists.txt index 2cc19db8f2..277f1924ec 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_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) @@ -281,6 +282,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_SW) + add_compile_definitions(__SW) + set(SW ON) + include_directories(${SW_MATH}/include) + include_directories(${SW_FFT}/include) + + 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) + target_link_libraries(${ABACUS_BIN_NAME} ${SW_MATH}/libswblas.a) +endif() find_package(Threads REQUIRED) target_link_libraries(${ABACUS_BIN_NAME} Threads::Threads) @@ -459,6 +472,9 @@ if(MKLROOT) if(CMAKE_CXX_COMPILER_ID MATCHES Intel) 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) find_package(Lapack 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; } 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_basis/module_pw/pw_gatherscatter.h b/source/source_basis/module_pw/pw_gatherscatter.h index 827182d997..e6b5998446 100644 --- a/source/source_basis/module_pw/pw_gatherscatter.h +++ b/source/source_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/source_cell/module_symmetry/symm_rho.cpp b/source/source_cell/module_symmetry/symm_rho.cpp index 05c4f322fe..754279ffbc 100644 --- a/source/source_cell/module_symmetry/symm_rho.cpp +++ b/source/source_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/source_cell/pseudo.cpp b/source/source_cell/pseudo.cpp index e3124bbf4d..233632069d 100644 --- a/source/source_cell/pseudo.cpp +++ b/source/source_cell/pseudo.cpp @@ -9,6 +9,34 @@ pseudo::~pseudo() { } +void pseudo::check_betar() +{ + bool min_flag = false; + for (int ib = 0; ib < nbeta; ib++) + { + for (int ir = 0; ir < mesh; ir++) + { + // 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)) + { + min_flag = true; + betar(ib, ir) = 0.0; + } + } + } + if (min_flag) + { + std::cout << "WARNING: some of potential function is set to zero cause of less than 1e-30.\n"; + } +} + void pseudo::print_pseudo(std::ofstream& ofs) { print_pseudo_vl(ofs); diff --git a/source/source_cell/pseudo.h b/source/source_cell/pseudo.h index 45f321a4eb..b654abcc8f 100644 --- a/source/source_cell/pseudo.h +++ b/source/source_cell/pseudo.h @@ -74,6 +74,12 @@ class pseudo // uspp ModuleBase::realArray qfuncl; // qfuncl(2*lmax+1,nbeta*(nbeta+1)/2,mesh) Q_{mu,nu}(|r|) function for |r|> 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..0e91dffccc 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_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 45351ca3cc..e0941664bf 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_betar(); // Read the hamiltonian terms D_ij if (ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "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/source_estate/module_charge/charge_mixing_residual.cpp b/source/source_estate/module_charge/charge_mixing_residual.cpp index 819dcbea92..9e9c5e0131 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,19 +200,19 @@ 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; } 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 @@ -298,14 +295,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 +318,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 +370,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,18 +407,19 @@ 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; - const int ig0 = this->rhopw->ig_gge0; if (ig0 > 0) { sum += fac2 diff --git a/source/source_estate/module_pot/H_Hartree_pw.cpp b/source/source_estate/module_pot/H_Hartree_pw.cpp index fab0230f2c..a9cecd0b64 100644 --- a/source/source_estate/module_pot/H_Hartree_pw.cpp +++ b/source/source_estate/module_pot/H_Hartree_pw.cpp @@ -44,18 +44,20 @@ 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) { - 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]; + 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/source_hamilt/module_ewald/H_Ewald_pw.cpp b/source/source_hamilt/module_ewald/H_Ewald_pw.cpp index 154fd0877d..20e6762cc0 100644 --- a/source/source_hamilt/module_ewald/H_Ewald_pw.cpp +++ b/source/source_hamilt/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/source_hamilt/module_surchem/cal_pseudo.cpp b/source/source_hamilt/module_surchem/cal_pseudo.cpp index c99324a5dc..1e0a64c5b5 100644 --- a/source/source_hamilt/module_surchem/cal_pseudo.cpp +++ b/source/source_hamilt/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/source_hamilt/module_surchem/cal_totn.cpp b/source/source_hamilt/module_surchem/cal_totn.cpp index b74e7e44d0..501046cb28 100644 --- a/source/source_hamilt/module_surchem/cal_totn.cpp +++ b/source/source_hamilt/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,17 +40,15 @@ void surchem::induced_charge(const UnitCell& cell, const ModulePW::PW_Basis* rho std::complex *induced_rhog = new std::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; } - 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_io/read_input.cpp b/source/source_io/read_input.cpp index 470f62903c..c59c31d8c6 100644 --- a/source/source_io/read_input.cpp +++ b/source/source_io/read_input.cpp @@ -196,6 +196,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 __SW ModuleBase::Global_File::make_dir_out(param.input.suffix, param.input.calculation, out_dir, @@ -203,13 +204,14 @@ void ReadInput::create_directory(const Parameter& param) 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 __SW if (system(ss.c_str())) { ModuleBase::WARNING_QUIT("ReadInput", "please set right files directory for reading in."); } - + #endif return; } 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/VNL_in_pw.cpp b/source/source_pw/module_pwdft/VNL_in_pw.cpp index 7aaf2dce45..e26e09fa3e 100644 --- a/source/source_pw/module_pwdft/VNL_in_pw.cpp +++ b/source/source_pw/module_pwdft/VNL_in_pw.cpp @@ -724,20 +724,18 @@ void pseudopot_cell_vnl::init_vnl(UnitCell& cell, const ModulePW::PW_Basis* rho_ double* jl = new double[kkbeta]; double* aux = new double[kkbeta]; - 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; + double vqint=0.0; ModuleBase::Integral::Simpson_Integral(kkbeta, aux, cell.atoms[it].ncpp.rab.data(), vqint); this->tab(it, ib, iq) = vqint * pref; } diff --git a/source/source_pw/module_pwdft/forces.cpp b/source/source_pw/module_pwdft/forces.cpp index 0fb9d63325..6aab09161f 100644 --- a/source/source_pw/module_pwdft/forces.cpp +++ b/source/source_pw/module_pwdft/forces.cpp @@ -472,14 +472,18 @@ 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 (ig== ig0) + { + continue; // skip G=0 + } aux[ig] *= ModuleBase::libm::exp(-1.0 * rho_basis->gg[ig] * ucell.tpiba2 / alpha / 4.0) - / (rho_basis->gg[ig] * ucell.tpiba2); + / (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 2e062bca8b..faa2de25b7 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,13 +159,15 @@ 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++) { + 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_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..b99ee7fb04 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,12 +68,11 @@ 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; } - + 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++) {