From e1f1dccd77f7a0315329f11696ead35f6e4abebd Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Fri, 14 Aug 2026 18:30:19 +0800 Subject: [PATCH 1/2] Feature: Enable correct Laplacian-dependent meta-GGA calculations --- source/source_hamilt/module_xc/libxc_pot.cpp | 33 ++++++++++--------- .../source_hamilt/module_xc/test/xc3_mock.h | 16 +++++++++ tests/01_PW/207_PW_SCANL/result.ref | 8 ++--- tests/01_PW/207_PW_SCANL/threshold | 10 +++--- 4 files changed, 43 insertions(+), 24 deletions(-) diff --git a/source/source_hamilt/module_xc/libxc_pot.cpp b/source/source_hamilt/module_xc/libxc_pot.cpp index 9252150b599..7e69a9a7804 100644 --- a/source/source_hamilt/module_xc/libxc_pot.cpp +++ b/source/source_hamilt/module_xc/libxc_pot.cpp @@ -462,32 +462,35 @@ std::tuple XC_Functional_Li voflapl(is,ir) += vlapl[ir*nspin+is] * sgn[ir*nspin+is]; } } -} + } // v_xc += nabla^2(vlapl) where vlapl = d(rho*eps_xc)/d(nabla^2 rho) if (need_laplacian) { const int ng = chr->rhopw->npw; const double tpiba2 = tpiba * tpiba; - std::vector> lapl_tmp(chr->rhopw->nmaxgr); - for(int is = 0; is < voflapl.nr; is++) + std::vector lapl_r(nrxx); + std::vector> lapl_g(ng); + for (int is = 0; is < voflapl.nr; is++) { - for(int ir = 0; ir < nrxx; ir++) - lapl_tmp[ir] = std::complex(voflapl(is, ir), 0.0); - for(int ig = ng; ig < chr->rhopw->nmaxgr; ig++) - lapl_tmp[ig] = std::complex(0.0, 0.0); - chr->rhopw->real2recip(lapl_tmp.data(), lapl_tmp.data()); - for(int ig = 0; ig < ng; ig++) + for (int ir = 0; ir < nrxx; ir++) + { + lapl_r[ir] = voflapl(is, ir); + } + chr->rhopw->real2recip(lapl_r.data(), lapl_g.data()); + for (int ig = 0; ig < ng; ig++) { double g2 = 0.0; - for(int i = 0; i < 3; i++) + for (int i = 0; i < 3; i++) + { g2 += chr->rhopw->gcar[ig][i] * chr->rhopw->gcar[ig][i]; - lapl_tmp[ig] *= -g2 * tpiba2; + } + lapl_g[ig] *= -g2 * tpiba2; } - chr->rhopw->recip2real(lapl_tmp.data(), lapl_tmp.data()); - for(int ir = 0; ir < nrxx; ir++) + chr->rhopw->recip2real(lapl_g.data(), lapl_r.data()); + for (int ir = 0; ir < nrxx; ir++) { - double vlapl_corr = ModuleBase::e2 * lapl_tmp[ir].real(); + double vlapl_corr = ModuleBase::e2 * lapl_r[ir]; v(is, ir) += vlapl_corr; vtxc += vlapl_corr * chr->rho[is][ir]; } @@ -511,4 +514,4 @@ std::tuple XC_Functional_Li return std::make_tuple( etxc, vtxc, std::move(v), std::move(vofk) ); } -#endif \ No newline at end of file +#endif diff --git a/source/source_hamilt/module_xc/test/xc3_mock.h b/source/source_hamilt/module_xc/test/xc3_mock.h index afa1017ee56..ab5dd9ea111 100644 --- a/source/source_hamilt/module_xc/test/xc3_mock.h +++ b/source/source_hamilt/module_xc/test/xc3_mock.h @@ -52,6 +52,22 @@ namespace ModulePW const bool add, const double factor) const; + template + void PW_Basis::recip2real(const std::complex* in, + FPTYPE* out, + const bool add, + const FPTYPE factor) const + { + for (int i = 0; i < nrxx; i++) + { + out[i] = (-ModuleBase::IMAG_UNIT * in[i]).real(); + } + } + template void PW_Basis::recip2real(const std::complex* in, + double* out, + const bool add, + const double factor) const; + template void PW_Basis_K::recip2real(const std::complex* in, std::complex* out, diff --git a/tests/01_PW/207_PW_SCANL/result.ref b/tests/01_PW/207_PW_SCANL/result.ref index 62907df9714..90b60e8a4e1 100644 --- a/tests/01_PW/207_PW_SCANL/result.ref +++ b/tests/01_PW/207_PW_SCANL/result.ref @@ -1,5 +1,5 @@ -etotref -205.28423347 -etotperatomref -102.64211673 -totalforceref 16.23426800 -totalstressref 1297.505181 +etotref -205.4194741913589439 +etotperatomref -102.7097370957 +totalforceref 16.622862 +totalstressref 1344.868111 totaltimeref 246.99 diff --git a/tests/01_PW/207_PW_SCANL/threshold b/tests/01_PW/207_PW_SCANL/threshold index 390d09ad3ca..42d33cecbb0 100644 --- a/tests/01_PW/207_PW_SCANL/threshold +++ b/tests/01_PW/207_PW_SCANL/threshold @@ -1,5 +1,5 @@ -# SCAN-L Laplacian is sensitive to the plane-wave cutoff, requiring -# relaxed thresholds compared to the default (1e-7/1e-4/1e-3) -threshold 0.0000003 -force_threshold 0.0005 -stress_threshold 0.02 \ No newline at end of file +# Libxc 5.1.7 and 7.1.x use different numerical implementations of SCAN-L. +# These thresholds cover both supported versions. +threshold 0.000002 +force_threshold 0.005 +stress_threshold 0.05 From d27363fa66bf0388c6fb5988da9426e589ae9fce Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Fri, 14 Aug 2026 19:01:59 +0800 Subject: [PATCH 2/2] Fix FFT in source/source_hamilt/module_xc/test/xc3_mock.h --- .../source_hamilt/module_xc/test/xc3_mock.h | 421 ++++++++++-------- 1 file changed, 225 insertions(+), 196 deletions(-) diff --git a/source/source_hamilt/module_xc/test/xc3_mock.h b/source/source_hamilt/module_xc/test/xc3_mock.h index ab5dd9ea111..83bc9ce08ef 100644 --- a/source/source_hamilt/module_xc/test/xc3_mock.h +++ b/source/source_hamilt/module_xc/test/xc3_mock.h @@ -4,250 +4,279 @@ namespace ModulePW { - PW_Basis::PW_Basis(){}; - PW_Basis::~PW_Basis(){}; - - template - void PW_Basis::real2recip(const FPTYPE* in, std::complex* out, const bool add, const FPTYPE factor) const - { - for (int i=0;i(const double* in, - std::complex* out, - bool add, - double factor) const; - - template - void PW_Basis::real2recip(const std::complex* in, - std::complex* out, - const bool add, - const FPTYPE factor) const +namespace +{ +// Preserve the shared prefix and use zero for mock outputs without a +// one-to-one counterpart between real- and reciprocal-space buffers. +template +void mock_real2recip(const InputType* in, std::complex* out, const int nrxx, const int nrecip, const bool add, const FPTYPE factor) +{ + for (int i = 0; i < nrecip; ++i) { - for (int i=0;i value = (i < nrxx) ? std::complex(in[i]) : std::complex(0.0, 0.0); + if (add) { - out[i] = in[i]; + out[i] += factor * value; } - } - template void PW_Basis::real2recip(const std::complex* in, - std::complex* out, - const bool add, - const double factor) const; - - template - void PW_Basis::recip2real(const std::complex* in, - std::complex* out, - const bool add, - const FPTYPE factor) const // in:(nz, ns) ; out(nplane,nx*ny) - { - for (int i=0;i* in, - std::complex* out, - const bool add, - const double factor) const; - - template - void PW_Basis::recip2real(const std::complex* in, - FPTYPE* out, - const bool add, - const FPTYPE factor) const +} + +template +void mock_recip2real(const std::complex* in, + std::complex* out, + const int nrecip, + const int nrxx, + const bool add, + const FPTYPE factor) +{ + for (int i = 0; i < nrxx; ++i) { - for (int i = 0; i < nrxx; i++) + const std::complex value = (i < nrecip) ? -ModuleBase::IMAG_UNIT * in[i] : std::complex(0.0, 0.0); + if (add) { - out[i] = (-ModuleBase::IMAG_UNIT * in[i]).real(); + out[i] += factor * value; } - } - template void PW_Basis::recip2real(const std::complex* in, - double* out, - const bool add, - const double factor) const; - - template - void PW_Basis_K::recip2real(const std::complex* in, - std::complex* out, - const int ik, - const bool add, - const FPTYPE factor) const // in:(nz, ns) ; out(nplane,nx*ny) - { - for (int i = 0; i < nrxx; i++) + else { - out[i] = -ModuleBase::IMAG_UNIT * in[i]; + out[i] = value; } } - template void PW_Basis_K::recip2real(const std::complex* in, - std::complex* out, - const int ik, - const bool add, - const double factor) const; - - ModuleBase::Vector3 PW_Basis_K::getgpluskcar(int, int) const - { - ModuleBase::Vector3 x = {1,2,3}; - return x; - } - +} - template - void PW_Basis_K::real_to_recip(const Device* ctx, - const std::complex* in, - std::complex* out, - const int ik, - const bool add, - const FPTYPE factor) const // in:(nplane,nx*ny) ; out(nz, ns) +template +void mock_recip2real(const std::complex* in, FPTYPE* out, const int nrecip, const int nrxx, const bool add, const FPTYPE factor) +{ + for (int i = 0; i < nrxx; ++i) { - for (int i=0;i - void PW_Basis_K::recip_to_real(const Device* ctx, - const std::complex* in, - std::complex* out, - const int ik, - const bool add, - const FPTYPE factor) const - { - for (int i = 0; i < nrxx; i++) + else { - out[i] = -ModuleBase::IMAG_UNIT * in[i]; + out[i] = value; } } +} +} // namespace - template void PW_Basis_K::real_to_recip(const base_device::DEVICE_CPU* ctx, - const std::complex* in, - std::complex* out, - const int ik, - const bool add, - const double factor) const; - template void PW_Basis_K::recip_to_real(const base_device::DEVICE_CPU* ctx, - const std::complex* in, - std::complex* out, - const int ik, - const bool add, - const double factor) const; -#if __CUDA || __ROCM - template void PW_Basis_K::real_to_recip(const base_device::DEVICE_GPU* ctx, - const std::complex* in, - std::complex* out, - const int ik, - const bool add, - const double factor) const; - - template void PW_Basis_K::recip_to_real(const base_device::DEVICE_GPU* ctx, - const std::complex* in, - std::complex* out, - const int ik, - const bool add, - const double factor) const; -#endif +PW_Basis::PW_Basis() {}; +PW_Basis::~PW_Basis() {}; +template +void PW_Basis::real2recip(const FPTYPE* in, std::complex* out, const bool add, const FPTYPE factor) const +{ + mock_real2recip(in, out, nrxx, npw, add, factor); +} +template void PW_Basis::real2recip(const double* in, std::complex* out, bool add, double factor) const; - void PW_Basis::initgrids(double, ModuleBase::Matrix3, double){}; - void PW_Basis::distribute_r(){}; - void PW_Basis::initgrids(double, ModuleBase::Matrix3, int, int, int){}; +template +void PW_Basis::real2recip(const std::complex* in, std::complex* out, const bool add, const FPTYPE factor) const +{ + mock_real2recip(in, out, nrxx, npw, add, factor); +} +template void PW_Basis::real2recip(const std::complex* in, + std::complex* out, + const bool add, + const double factor) const; + +template +void PW_Basis::recip2real(const std::complex* in, + std::complex* out, + const bool add, + const FPTYPE factor) const // in:(nz, ns) ; out(nplane,nx*ny) +{ + mock_recip2real(in, out, npw, nrxx, add, factor); +} +template void PW_Basis::recip2real(const std::complex* in, std::complex* out, const bool add, const double factor) const; - PW_Basis_K::PW_Basis_K(){}; - PW_Basis_K::~PW_Basis_K(){}; +template +void PW_Basis::recip2real(const std::complex* in, FPTYPE* out, const bool add, const FPTYPE factor) const +{ + mock_recip2real(in, out, npw, nrxx, add, factor); +} +template void PW_Basis::recip2real(const std::complex* in, double* out, const bool add, const double factor) const; + +template +void PW_Basis_K::recip2real(const std::complex* in, + std::complex* out, + const int ik, + const bool add, + const FPTYPE factor) const // in:(nz, ns) ; out(nplane,nx*ny) +{ + mock_recip2real(in, out, npwk[ik], nrxx, add, factor); } +template void PW_Basis_K::recip2real(const std::complex* in, + std::complex* out, + const int ik, + const bool add, + const double factor) const; -namespace ModuleBase +ModuleBase::Vector3 PW_Basis_K::getgpluskcar(int, int) const { - void WARNING_QUIT(const std::string &file,const std::string &description) - { - std::cout << " " << file <<" warning : "<< description< x = {1, 2, 3}; + return x; +} - void Matrix3::Identity(){}; +template +void PW_Basis_K::real_to_recip(const Device* ctx, + const std::complex* in, + std::complex* out, + const int ik, + const bool add, + const FPTYPE factor) const // in:(nplane,nx*ny) ; out(nz, ns) +{ + mock_real2recip(in, out, nrxx, npwk[ik], add, factor); +} +template +void PW_Basis_K::recip_to_real(const Device* ctx, + const std::complex* in, + std::complex* out, + const int ik, + const bool add, + const FPTYPE factor) const +{ + mock_recip2real(in, out, npwk[ik], nrxx, add, factor); +} + +template void PW_Basis_K::real_to_recip(const base_device::DEVICE_CPU* ctx, + const std::complex* in, + std::complex* out, + const int ik, + const bool add, + const double factor) const; +template void PW_Basis_K::recip_to_real(const base_device::DEVICE_CPU* ctx, + const std::complex* in, + std::complex* out, + const int ik, + const bool add, + const double factor) const; +#if __CUDA || __ROCM +template void PW_Basis_K::real_to_recip(const base_device::DEVICE_GPU* ctx, + const std::complex* in, + std::complex* out, + const int ik, + const bool add, + const double factor) const; + +template void PW_Basis_K::recip_to_real(const base_device::DEVICE_GPU* ctx, + const std::complex* in, + std::complex* out, + const int ik, + const bool add, + const double factor) const; +#endif - IntArray::IntArray(int,int){}; - IntArray::~IntArray(){}; +void PW_Basis::initgrids(double, ModuleBase::Matrix3, double) {}; +void PW_Basis::distribute_r() {}; +void PW_Basis::initgrids(double, ModuleBase::Matrix3, int, int, int) {}; - void TITLE(const std::string &class_function_name,bool disable){}; - void TITLE(const std::string &class_name,const std::string &function_name,bool disable){}; +PW_Basis_K::PW_Basis_K() {}; +PW_Basis_K::~PW_Basis_K() {}; +} // namespace ModulePW +namespace ModuleBase +{ +void WARNING_QUIT(const std::string& file, const std::string& description) +{ + std::cout << " " << file << " warning : " << description << std::endl; + exit(1); } +void WARNING(const std::string& file, const std::string& description) {}; + +void Matrix3::Identity() {}; + +IntArray::IntArray(int, int) {}; +IntArray::~IntArray() {}; + +void TITLE(const std::string& class_function_name, bool disable) {}; +void TITLE(const std::string& class_name, const std::string& function_name, bool disable) {}; + +} // namespace ModuleBase namespace GlobalV { - std::string BASIS_TYPE = ""; - bool CAL_STRESS = false; - int CAL_FORCE = 0; - int NSPIN; - int NPOL; - bool DOMAG; - bool DOMAG_Z; - std::ofstream ofs_device; - std::ofstream ofs_running; -} +std::string BASIS_TYPE = ""; +bool CAL_STRESS = false; +int CAL_FORCE = 0; +int NSPIN; +int NPOL; +bool DOMAG; +bool DOMAG_Z; +std::ofstream ofs_device; +std::ofstream ofs_running; +} // namespace GlobalV namespace GlobalC { - Exx_Info exx_info; +Exx_Info exx_info; } -UnitCell::UnitCell(){}; -UnitCell::~UnitCell(){}; +UnitCell::UnitCell() {}; +UnitCell::~UnitCell() {}; -Charge::Charge(){}; -Charge::~Charge(){}; +Charge::Charge() {}; +Charge::~Charge() {}; -Magnetism::Magnetism(){}; -Magnetism::~Magnetism(){}; +Magnetism::Magnetism() {}; +Magnetism::~Magnetism() {}; -SepPot::SepPot(){} -SepPot::~SepPot(){} -Sep_Cell::Sep_Cell() noexcept {} -Sep_Cell::~Sep_Cell() noexcept {} +SepPot::SepPot() +{ +} +SepPot::~SepPot() +{ +} +Sep_Cell::Sep_Cell() noexcept +{ +} +Sep_Cell::~Sep_Cell() noexcept +{ +} namespace unitcell { - void cal_ux(UnitCell& ucell, const int nspin) - { - ucell.magnet.lsign_ = false; - - ucell.magnet.ux_[0] = 0; - ucell.magnet.ux_[1] = 1; - ucell.magnet.ux_[2] = 2; - - ucell.magnet.lsign_ = true; - }; -} +void cal_ux(UnitCell& ucell, const int nspin) +{ + ucell.magnet.lsign_ = false; + ucell.magnet.ux_[0] = 0; + ucell.magnet.ux_[1] = 1; + ucell.magnet.ux_[2] = 2; + ucell.magnet.lsign_ = true; +}; +} // namespace unitcell namespace Parallel_Reduce { - /// reduce in all process - template - void reduce_all(T& object){}; - template - void reduce_all(T* object, const int n){}; - template - void reduce_pool(T& object){}; - template - void reduce_pool(T* object, const int n){}; - - template<> - void Parallel_Reduce::reduce_pool(double& object) - { - #ifdef __MPI - double swap = object; - MPI_Allreduce(&swap , &object , 1, MPI_DOUBLE , MPI_SUM , MPI_COMM_WORLD); - #endif - return; - } - template void reduce_all(double& object); - template void reduce_all(double* object, const int n); - template void reduce_pool(float& object); - template void reduce_pool(float* object, const int n); - template void reduce_pool(double* object, const int n); +/// reduce in all process +template +void reduce_all(T& object) {}; +template +void reduce_all(T* object, const int n) {}; +template +void reduce_pool(T& object) {}; +template +void reduce_pool(T* object, const int n) {}; + +template <> +void Parallel_Reduce::reduce_pool(double& object) +{ +#ifdef __MPI + double swap = object; + MPI_Allreduce(&swap, &object, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); +#endif + return; } +template void reduce_all(double& object); +template void reduce_all(double* object, const int n); +template void reduce_pool(float& object); +template void reduce_pool(float* object, const int n); +template void reduce_pool(double* object, const int n); +} // namespace Parallel_Reduce