diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 87bf922499..5af84d0711 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -115,6 +115,11 @@ jobs: cmake --build build -j8 cmake --install build + - name: rVV10 host contract + run: | + source toolchain/install/setup + ctest --test-dir build --output-on-failure --no-tests=error -R RVV10_LIBVDWXC_CONTRACT + - name: Report compiler cache if: always() run: | diff --git a/CMakeLists.txt b/CMakeLists.txt index 62843068f6..dd8b64261d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -46,6 +46,7 @@ option(USE_SW "Enable SW Architecture" OFF) option(ENABLE_ABACUS_LIBM "Build libmath from source to speed up" OFF) option(ENABLE_LIBXC "Enable using the LibXC package" OFF) +option(ENABLE_LIBVDWXC "Enable libvdwxc-backed nonlocal XC" OFF) option(ENABLE_FLOAT_FFTW "Enable using single-precision FFTW library." OFF) option(ENABLE_LCAO "Enable LCAO algorithm" ON) @@ -618,6 +619,17 @@ if(ENABLE_LIBXC) endif() endif() +if(ENABLE_LIBVDWXC) + if(NOT ENABLE_LIBXC) + message(FATAL_ERROR + "ENABLE_LIBVDWXC currently requires ENABLE_LIBXC=ON because the first " + "rVV10 interface uses Libxc for the semilocal RPW86+PBE part.") + endif() + # The config target supplies headers and transitive FFT dependencies. + # rVV10 currently requires the development revision documented in install.md. + find_package(libvdwxc CONFIG REQUIRED) +endif() + if(DEFINED DeePMD_DIR) find_package(DeePMD REQUIRED) endif() diff --git a/cmake/CollectBuildInfoVars.cmake b/cmake/CollectBuildInfoVars.cmake index 740c50ca4b..15ec9e43e8 100644 --- a/cmake/CollectBuildInfoVars.cmake +++ b/cmake/CollectBuildInfoVars.cmake @@ -181,6 +181,14 @@ else() set(ABACUS_DFTD4_VERSION "no") endif() +if(ENABLE_LIBVDWXC AND libvdwxc_VERSION) + set(ABACUS_LIBVDWXC_VERSION "yes (v${libvdwxc_VERSION})") +elseif(ENABLE_LIBVDWXC) + set(ABACUS_LIBVDWXC_VERSION "yes (version unknown)") +else() + set(ABACUS_LIBVDWXC_VERSION "no") +endif() + # Accelerators if(USE_CUDA AND CUDAToolkit_VERSION) set(ABACUS_CUDA_VERSION "yes (v${CUDAToolkit_VERSION})") @@ -382,6 +390,7 @@ list(APPEND CMAKE_OPTIONS_LIST " USE_ROCM=${USE_ROCM}") list(APPEND CMAKE_OPTIONS_LIST " ENABLE_LCAO=${ENABLE_LCAO}") list(APPEND CMAKE_OPTIONS_LIST " ENABLE_ELPA=${ENABLE_ELPA}") list(APPEND CMAKE_OPTIONS_LIST " ENABLE_LIBXC=${ENABLE_LIBXC}") +list(APPEND CMAKE_OPTIONS_LIST " ENABLE_LIBVDWXC=${ENABLE_LIBVDWXC}") list(APPEND CMAKE_OPTIONS_LIST " ENABLE_MLALGO=${ENABLE_MLALGO}") list(APPEND CMAKE_OPTIONS_LIST " ENABLE_ASAN=${ENABLE_ASAN}") list(APPEND CMAKE_OPTIONS_LIST " BUILD_TESTING=${BUILD_TESTING}") @@ -396,6 +405,7 @@ list(APPEND CMAKE_FIND_PACKAGES_LIST " CUDA Found=${CUDA_FOUND}") list(APPEND CMAKE_FIND_PACKAGES_LIST " ROCm Found=${ROCM_FOUND}") list(APPEND CMAKE_FIND_PACKAGES_LIST " ELPA Found=${ELPA_FOUND}") list(APPEND CMAKE_FIND_PACKAGES_LIST " Libxc Found=${Libxc_FOUND}") +list(APPEND CMAKE_FIND_PACKAGES_LIST " libvdwxc Found=${libvdwxc_FOUND}") list(APPEND CMAKE_FIND_PACKAGES_LIST " PEXSI Found=${PEXSI_FOUND}") list(APPEND CMAKE_FIND_PACKAGES_LIST " LibRI Found=${LIBRI_FOUND}") foreach(package_line ${CMAKE_FIND_PACKAGES_LIST}) @@ -439,6 +449,7 @@ message(" Torch = ${ABACUS_LIBTORCH_VERSION}") message(" DeePMD = ${ABACUS_DEEPMD_VERSION}") message(" NEP = ${ABACUS_NEP_VERSION}") message(" Features : Libxc = ${ABACUS_LIBXC_VERSION}") +message(" libvdwxc = ${ABACUS_LIBVDWXC_VERSION}") message(" LibRI = ${ABACUS_LIBRI_VERSION}") message(" DFT-D4 = ${ABACUS_DFTD4_VERSION}") message(" JSON = ${ABACUS_JSON_VERSION}") diff --git a/docs/CITATIONS.md b/docs/CITATIONS.md index 43176829a1..da0e68dadc 100644 --- a/docs/CITATIONS.md +++ b/docs/CITATIONS.md @@ -55,3 +55,11 @@ The following references are required to be cited when using ABACUS. Specificall Eike Caldeweyher, Jan-Michael Mewes, Sebastian Ehlert, and Stefan Grimme. "Extension and evaluation of the D4 London-dispersion model for periodic systems." Physical Chemistry Chemical Physics 22, 8499-8512 (2020). DOI: 10.1039/D0CP00502A. Nikolay V. Tkachenko, Linus B. Dittmer, Rebecca Tomann, and Martin Head-Gordon. "Smooth D4S model." The Journal of Physical Chemistry Letters 15, 10629-10637 (2024). DOI: 10.1021/acs.jpclett.4c02653. + +- **If the libvdwxc-backed rVV10 nonlocal correction is used:** + + O. A. Vydrov and T. Van Voorhis. "Nonlocal van der Waals density functional: The simpler the better." The Journal of Chemical Physics 133, 244103 (2010). DOI: 10.1063/1.3521275. + + P. Sabatini, T. Gorni, and S. de Gironcoli. "Nonlocal van der Waals density functional made simple and efficient." Physical Review B 87, 041108(R) (2013). DOI: 10.1103/PhysRevB.87.041108. + + Cite the exact [libvdwxc](https://gitlab.com/libvdwxc/libvdwxc) revision used for the calculation, because the rVV10 CMake interface is version-dependent during development. diff --git a/docs/advanced/install.md b/docs/advanced/install.md index ab1939181b..6c2606ab9a 100644 --- a/docs/advanced/install.md +++ b/docs/advanced/install.md @@ -99,6 +99,44 @@ vdw_d4_model d4 # or d4s for the smooth D4S model If `vdw_d4_xc` is set to `default`, ABACUS will infer the functional name from `dft_functional` or pseudopotential metadata and pass it to the DFT-D4 library. The `vdw_d4_model` keyword selects the dispersion model inside the DFT-D4 library; the default is `d4`, while `d4s` enables the smooth D4S model. +## Build with libvdwxc-backed rVV10 support + +The first rVV10 integration uses the external [libvdwxc](https://gitlab.com/libvdwxc/libvdwxc) library for the nonlocal kernel. The current rVV10 API is under review in Draft MR !32, so this option is intended for development and reproducible integration testing rather than a stable release build. + +Build and install the pinned libvdwxc development branch first (MR !32, +`rvv10-cmake`, commit +`94410ebdf3dc9615a40f40c3a5b59dbe20a433d2`), then pass its CMake prefix to +ABACUS together with Libxc. This revision fixes Hermitian symmetry on stored +FFT boundary planes in oblique cells. With GCC 11.2.1, its eight library tests +pass using either FFTW 3.3.11 or the MKL 2022.2.1 compatibility backend. +The ABACUS adapter is tested with standalone FFTW for both programs; the +library's MKL tests alone do not validate an MKL-based ABACUS build. Do not +mix the two FFT providers in one process. To use the validated serial +configuration, leave `MKLROOT`/`MKL_ROOT` unset and select OpenBLAS/LAPACK plus +FFTW through ABACUS's existing dependency options: + +```bash +cmake -B build-libvdwxc -DBUILD_TESTING=ON \ + -DENABLE_MPI=OFF -DENABLE_LCAO=OFF \ + -DENABLE_LIBXC=ON -DENABLE_LIBVDWXC=ON \ + -DCMAKE_PREFIX_PATH="/path/to/libvdwxc-prefix;/path/to/libxc-prefix" \ + -DFFTW3_DIR=/path/to/fftw-prefix \ + -DBLA_VENDOR=OpenBLAS \ + -DBLAS_DIR=/path/to/openblas-prefix \ + -DLAPACK_DIR=/path/to/openblas-prefix +``` + +Enable the correction with `xc_nonlocal rvv10` and explicitly select +`dft_functional GGA_X_RPW86+GGA_C_PBE`; this first parameterization uses +`rvv10_b 6.3` and `rvv10_c 0.0093`. + +The first adapter supports CPU double-precision, full-complex PW, one process, +`nspin=1`, a single density grid, a right-handed cell and norm-conserving +pseudopotentials for SCF energy/potential evaluation. Other spin/layout modes, +forces, stress, relaxation and pairwise D2/D3/D4 combinations are rejected. +The legacy Makefile does not build this backend. See `tests/rvv10/README.md` +for the numerical tests and host/backend accounting. + ## Build with JSON support ABACUS uses [nlohmann-json](https://github.com/nlohmann/json) for optional JSON output. To build ABACUS with JSON support, pass `-DENABLE_JSON=ON` to CMake. diff --git a/generate_build_info.sh b/generate_build_info.sh index f50f98ac72..48f46416f6 100755 --- a/generate_build_info.sh +++ b/generate_build_info.sh @@ -165,6 +165,12 @@ if [ "${USE_CUDA}" == "ON" ]; then fi fi +# 6b. libvdwxc +# +# The legacy Makefile does not compile the adapter, even if an environment +# variable requests it. Report the compiled feature, not a discovered library. +LIBVDWXC_VERSION="no" + # 7. ROCm/HIP ROCM_VERSION="no" if [ "${USE_ROCM}" == "ON" ]; then @@ -279,6 +285,7 @@ sed \ -e "s#@ABACUS_OPENMP_VERSION@#$OPENMP_VERSION#g" \ -e "s#@ABACUS_MKL_SUPPORT@#$MKL_SUPPORT#g" \ -e "s#@ABACUS_LIBXC_VERSION@#$LIBXC_VERSION#g" \ + -e "s#@ABACUS_LIBVDWXC_VERSION@#$LIBVDWXC_VERSION#g" \ -e "s#@ABACUS_FFTW_VERSION@#$FFTW3_VERSION#g" \ -e "s#@ABACUS_CUDA_VERSION@#$CUDA_VERSION#g" \ -e "s#@ABACUS_ROCM_VERSION@#$ROCM_VERSION#g" \ diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index 2322dda497..588f39c317 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -66,6 +66,7 @@ set(_abacus_feature_definitions $<$:TEST_EXX_LCAO=0> $<$:TEST_EXX_RADIAL=1> $<$:__LIBXC> + $<$:__LIBVDWXC> $<$:__DFTD4> $<$:GATHER_INFO>) @@ -270,7 +271,8 @@ set(_abacus_feature_libs $<$:ZLIB::ZLIB> $<$:cereal::cereal> $<$:Libxc::xc> - $<$:PLUMED::plumed>) + $<$:PLUMED::plumed> + $<$:libvdwxc::vdwxc>) set(_abacus_feature_include_dirs) set(_abacus_feature_compile_options) @@ -392,6 +394,7 @@ foreach(_abacus_feature_target IN ITEMS cereal::cereal ELPA::ELPA Libxc::xc + libvdwxc::vdwxc dftd4::dftd4 DeePMD::deepmd_c DeePMD::deepmd_cc @@ -714,6 +717,34 @@ endif() # base uses symbols from device (memory_op and math_ylm_op). target_link_libraries(base PUBLIC device) +if(BUILD_TESTING) + # Register after all object-library dependencies exist; the optional adapter + # test links the existing PW, INPUT and device targets. + add_subdirectory(source_hamilt/module_vdw/test) + if(ENABLE_LIBVDWXC AND NOT ENABLE_MPI AND NOT USE_CUDA AND NOT USE_ROCM) + # Test the additive host with real Charge and UnitCell lifecycles after + # its dependencies exist, without linking the full estate/cell libraries. + AddTest( + TARGET MODULE_ESTATE_pot_nonlocal_xc_LIBVDWXC + LIBS xc_ planewave base parameter device nonlocal_xc + SOURCES + source_estate/test/test_pot_nonlocal_xc.cpp + source_estate/module_pot/pot_nonlocal_xc.cpp + source_estate/module_charge/charge.cpp + source_estate/module_charge/chg_tools.cpp + source_cell/unitcell.cpp + source_cell/atom_spec.cpp + source_cell/atom_pseudo.cpp + source_cell/pseudo.cpp + source_cell/magnetism.cpp + source_cell/sep.cpp + source_cell/sep_cell.cpp + ) + target_link_libraries(MODULE_ESTATE_pot_nonlocal_xc_LIBVDWXC PRIVATE Libxc::xc libvdwxc::vdwxc) + set_tests_properties(MODULE_ESTATE_pot_nonlocal_xc_LIBVDWXC PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=1") + endif() +endif() + if(ENABLE_COVERAGE) add_coverage(driver) endif() @@ -738,6 +769,7 @@ target_link_libraries( module_stodft module_dfpt xc_ + nonlocal_xc vdw relax io_advanced diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 34b39f96da..12c2a4bd0b 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -289,6 +289,7 @@ OBJS_ELECSTAT=elecstate.o\ h_hartree_pw.o\ h_tddft_pw.o\ pot_xc.o\ + pot_nonlocal_xc.o\ pot_ml_exx.o\ pot_ml_exx_label.o\ pot_xc_fdm.o\ @@ -1014,6 +1015,7 @@ OBJS_SRCPW=h_ewald_pw.o\ vsep_pw.o OBJS_VDW=vdw.o\ + nonlocal_xc.o\ vdwd2_parameters.o\ vdwd3_parameters.o\ vdwd2.o\ diff --git a/source/source_estate/CMakeLists.txt b/source/source_estate/CMakeLists.txt index 15c09ac3ce..dc2a7a9977 100644 --- a/source/source_estate/CMakeLists.txt +++ b/source/source_estate/CMakeLists.txt @@ -19,6 +19,7 @@ list(APPEND objects module_pot/h_hartree_pw.cpp module_pot/h_tddft_pw.cpp module_pot/pot_xc.cpp + module_pot/pot_nonlocal_xc.cpp module_pot/pot_local.cpp module_pot/pot_pcc.cpp module_pot/potential_new.cpp diff --git a/source/source_estate/module_pot/pot_nonlocal_xc.cpp b/source/source_estate/module_pot/pot_nonlocal_xc.cpp new file mode 100644 index 0000000000..c41f033a01 --- /dev/null +++ b/source/source_estate/module_pot/pot_nonlocal_xc.cpp @@ -0,0 +1,59 @@ +#include "pot_nonlocal_xc.h" + +#include "source_hamilt/module_vdw/nonlocal_xc.h" +#include "source_base/timer.h" +#include "source_base/tool_quit.h" + +#include + +namespace elecstate +{ +PotNonlocalXC::PotNonlocalXC(const ModulePW::PW_Basis* density_basis, + double* etxc, + double* vtxc, + const std::string& method, + double b, + double c) + : etxc_(etxc), vtxc_(vtxc) +{ + backend_ = vdw::make_nonlocal_xc(density_basis, method, b, c); + rho_basis_ = density_basis; + dynamic_mode = true; + fixed_mode = false; +} + +PotNonlocalXC::~PotNonlocalXC() = default; + +void PotNonlocalXC::cal_v_eff(const Charge* chg, const UnitCell* cell, ModuleBase::matrix& v_eff) +{ + ModuleBase::TITLE("PotNonlocalXC", "cal_v_eff"); + ModuleBase::timer::start("PotNonlocalXC", "cal_v_eff"); + if (!chg || !cell || !rho_basis_ || chg->rhopw != rho_basis_ || chg->nspin != 1 + || chg->nrxx != rho_basis_->nrxx || !chg->rho || !chg->rho_core || v_eff.nr != chg->nspin + || v_eff.nc != chg->nrxx) + { + ModuleBase::WARNING_QUIT("PotNonlocalXC", "the first rVV10 interface requires a consistent nspin=1 density basis"); + } + std::vector valence(chg->nrxx, 0.0); + for (int is = 0; is < chg->nspin; ++is) + for (int r = 0; r < chg->nrxx; ++r) + valence[r] += chg->rho[is][r]; + std::vector total = valence; + for (int r = 0; r < chg->nrxx; ++r) + total[r] += chg->rho_core[r]; + const vdw::NonlocalXCResult nl = backend_->evaluate(*rho_basis_, total, valence); + if (nl.potential.size() != static_cast(chg->nrxx)) + { + ModuleBase::WARNING_QUIT("PotNonlocalXC", "nonlocal XC backend returned an invalid potential size"); + } + + // PotXC is registered before this component and owns the semilocal + // contribution. This component only adds the selected nonlocal correction. + *etxc_ += nl.energy; + *vtxc_ += nl.vtxc; + for (int is = 0; is < chg->nspin; ++is) + for (int r = 0; r < chg->nrxx; ++r) + v_eff(is, r) += nl.potential[r]; + ModuleBase::timer::end("PotNonlocalXC", "cal_v_eff"); +} +} // namespace elecstate diff --git a/source/source_estate/module_pot/pot_nonlocal_xc.h b/source/source_estate/module_pot/pot_nonlocal_xc.h new file mode 100644 index 0000000000..feafc824bd --- /dev/null +++ b/source/source_estate/module_pot/pot_nonlocal_xc.h @@ -0,0 +1,34 @@ +#ifndef ABACUS_POT_NONLOCAL_XC_H +#define ABACUS_POT_NONLOCAL_XC_H +#include "pot_base.h" + +#include +#include + +namespace vdw +{ +class NonlocalXC; +} + +namespace elecstate +{ +// Independent rVV10 nonlocal correction for the CPU/PW SCF path. +class PotNonlocalXC final : public PotBase +{ + public: + PotNonlocalXC(const ModulePW::PW_Basis* density_basis, + double* etxc, + double* vtxc, + const std::string& method, + double b, + double c); + ~PotNonlocalXC() override; + void cal_v_eff(const Charge* chg, const UnitCell* cell, ModuleBase::matrix& v_eff) override; + + private: + double* etxc_; + double* vtxc_; + std::unique_ptr backend_; +}; +} // namespace elecstate +#endif // ABACUS_POT_NONLOCAL_XC_H diff --git a/source/source_estate/module_pot/potential_new.cpp b/source/source_estate/module_pot/potential_new.cpp index b23b8710bb..73cb7f242c 100644 --- a/source/source_estate/module_pot/potential_new.cpp +++ b/source/source_estate/module_pot/potential_new.cpp @@ -24,8 +24,36 @@ Potential::Potential(const ModulePW::PW_Basis* rho_basis_in, double* etxc_in, double* vtxc_in, VSep* vsep_cell_in) + : Potential(rho_basis_in, + rho_basis_smooth_in, + ucell_in, + vloc_in, + structure_factors_in, + solvent_in, + etxc_in, + vtxc_in, + vsep_cell_in, + 6.3, + 0.0093, + "none") +{ +} + +Potential::Potential(const ModulePW::PW_Basis* rho_basis_in, + const ModulePW::PW_Basis* rho_basis_smooth_in, + const UnitCell* ucell_in, + const ModuleBase::matrix* vloc_in, + Structure_Factor* structure_factors_in, + surchem* solvent_in, + double* etxc_in, + double* vtxc_in, + VSep* vsep_cell_in, + double nonlocal_xc_b_in, + double nonlocal_xc_c_in, + const std::string& nonlocal_xc_in) : ucell_(ucell_in), vloc_(vloc_in), structure_factors_(structure_factors_in), solvent_(solvent_in), vsep_cell(vsep_cell_in), - etxc_(etxc_in), vtxc_(vtxc_in) + etxc_(etxc_in), vtxc_(vtxc_in), nonlocal_xc_b_(nonlocal_xc_b_in), nonlocal_xc_c_(nonlocal_xc_c_in), + nonlocal_xc_(nonlocal_xc_in) { this->rho_basis_ = rho_basis_in; this->rho_basis_smooth_ = rho_basis_smooth_in; diff --git a/source/source_estate/module_pot/potential_new.h b/source/source_estate/module_pot/potential_new.h index e1af826981..177365c8e1 100644 --- a/source/source_estate/module_pot/potential_new.h +++ b/source/source_estate/module_pot/potential_new.h @@ -8,6 +8,7 @@ #include "source_pw/module_pwdft/vsep_pw.h" #include +#include #include namespace elecstate @@ -67,6 +68,21 @@ class Potential : public PotBase double* etxc_in, double* vtxc_in, VSep* vsep_cell_in = nullptr); + // Constructor used when a density-based nonlocal XC correction is + // selected. Keeping the method in the immutable Potential configuration + // avoids a mutable workflow switch after construction. + Potential(const ModulePW::PW_Basis* rho_basis_in, + const ModulePW::PW_Basis* rho_basis_smooth_in, + const UnitCell* ucell_in, + const ModuleBase::matrix* vloc_in, + Structure_Factor* structure_factors_in, + surchem* solvent_in, + double* etxc_in, + double* vtxc_in, + VSep* vsep_cell_in, + double nonlocal_xc_b_in, + double nonlocal_xc_c_in, + const std::string& nonlocal_xc_in); ~Potential(); // initialize potential when SCF begin @@ -189,6 +205,12 @@ class Potential : public PotBase { return this->rho_basis_; } + + bool has_nonlocal_xc() const + { + return !nonlocal_xc_.empty() && nonlocal_xc_ != "none"; + } + const UnitCell* get_ucell() const { return this->ucell_; @@ -259,6 +281,9 @@ class Potential : public PotBase Structure_Factor* structure_factors_ = nullptr; surchem* solvent_ = nullptr; VSep* vsep_cell = nullptr; + double nonlocal_xc_b_ = 6.3; + double nonlocal_xc_c_ = 0.0093; + std::string nonlocal_xc_ = "none"; bool use_gpu_ = false; std::shared_ptr td_field_manager_; }; diff --git a/source/source_estate/module_pot/potential_types.cpp b/source/source_estate/module_pot/potential_types.cpp index 4280bb8210..f154767ce5 100644 --- a/source/source_estate/module_pot/potential_types.cpp +++ b/source/source_estate/module_pot/potential_types.cpp @@ -3,6 +3,7 @@ #include "gatefield.h" #include "pot_local.h" #include "pot_pcc.h" +#include "pot_nonlocal_xc.h" #include "pot_sep.h" #include "pot_surchem.h" #include "pot_xc.h" @@ -14,6 +15,7 @@ #include "source_base/tool_title.h" #include "source_io/module_parameter/parameter.h" #include "h_tddft_pw.h" + #ifdef __MLALGO #include "pot_ml_exx.h" #endif @@ -47,6 +49,28 @@ PotBase* Potential::get_pot_type(const std::string& pot_type) const int open_axis = input.pcc_2d_axis; return new PotPcc(this->rho_basis_, PotPcc::Dimension::slab, open_axis, electron_count); } + else if (pot_type == "xc_nonlocal") + { + // USPP can use a single density grid, so pointer equality does not + // determine the pseudopotential type. + for (int it = 0; it < this->ucell_->ntype; ++it) + { + if (this->ucell_->atoms[it].ncpp.tvanp) + { + ModuleBase::WARNING_QUIT("Potential", "nonlocal XC currently supports norm-conserving pseudopotentials only"); + } + } + if (this->rho_basis_ != this->rho_basis_smooth_) + { + ModuleBase::WARNING_QUIT("Potential", "nonlocal XC currently requires the same smooth and dense density basis"); + } + return new PotNonlocalXC(this->rho_basis_, + this->etxc_, + this->vtxc_, + this->nonlocal_xc_, + this->nonlocal_xc_b_, + this->nonlocal_xc_c_); + } else if (pot_type == "surchem") { return new PotSurChem(this->rho_basis_, this->structure_factors_, this->v_eff_fixed.data(), this->solvent_); diff --git a/source/source_estate/setup_estate_pw.cpp b/source/source_estate/setup_estate_pw.cpp index 9e045377cf..a134715f88 100644 --- a/source/source_estate/setup_estate_pw.cpp +++ b/source/source_estate/setup_estate_pw.cpp @@ -103,7 +103,8 @@ void setup_estate_pw_impl( { pelec->pot = new elecstate::Potential(pw_rhod, pw_rho, &ucell, &locpp.vloc, &sf, - &solvent, &(pelec->f_en.etxc), &(pelec->f_en.vtxc), vsep_cell); + &solvent, &(pelec->f_en.etxc), &(pelec->f_en.vtxc), vsep_cell, + inp.rvv10_b, inp.rvv10_c, inp.xc_nonlocal); } locpp.init_vloc(ucell, pw_rhod); diff --git a/source/source_estate/test/test_pot_nonlocal_xc.cpp b/source/source_estate/test/test_pot_nonlocal_xc.cpp new file mode 100644 index 0000000000..2df6f7ba14 --- /dev/null +++ b/source/source_estate/test/test_pot_nonlocal_xc.cpp @@ -0,0 +1,70 @@ +#include "source_estate/module_pot/pot_nonlocal_xc.h" +#include "source_hamilt/module_vdw/nonlocal_xc_libvdwxc.h" + +#include +#include +#include + +TEST(PotNonlocalXCTest, AddsCorrectionFromValenceAndCoreWithoutReplacingExistingXc) +{ + ModulePW::PW_Basis pw("cpu", "double"); + const ModuleBase::Matrix3 lattice(1.0, 0.0, 0.0, -0.5, std::sqrt(0.75), 0.0, 0.0, 0.0, 2.0); + pw.initgrids(10.0, lattice, 6, 8, 12); + pw.initparameters(false, 4.0, 1, true); + pw.setuptransform(); + pw.collect_local_pw(); + + UnitCell cell; + cell.lat0 = pw.lat0; + cell.latvec = lattice; + cell.omega = pw.omega; + Charge charge; + charge.set_rhopw(&pw); + charge.allocate(1, false, false, 0); + + const double pi = std::acos(-1.0); + std::vector valence(pw.nrxx); + std::vector core(pw.nrxx); + std::vector total(pw.nrxx); + ModuleBase::matrix potential(1, pw.nrxx); + std::vector initial_potential(pw.nrxx); + for (int x = 0; x < pw.nx; ++x) + for (int y = 0; y < pw.ny; ++y) + for (int z = 0; z < pw.nz; ++z) + { + const int r = (x * pw.ny + y) * pw.nz + z; + const double ax = 2.0 * pi * x / pw.nx; + const double ay = 2.0 * pi * y / pw.ny; + const double az = 2.0 * pi * z / pw.nz; + valence[r] = 0.007 + 0.001 * std::cos(ax) + 0.0005 * std::sin(ay + az); + core[r] = 0.003 + 0.0008 * std::cos(ay); + total[r] = valence[r] + core[r]; + charge.rho[0][r] = valence[r]; + charge.rho_core[r] = core[r]; + initial_potential[r] = 0.25 + 0.0001 * r; + potential(0, r) = initial_potential[r]; + } + + // The separately tested adapter supplies the correction; this test checks + // the production host's density construction and additive XC bookkeeping. + vdw::LibvdwxcNonlocalXC reference(pw, 6.3, 0.0093); + const vdw::NonlocalXCResult correction = reference.evaluate(pw, total, valence); + const double initial_energy = -1.25; + const double initial_vtxc = 0.375; + double energy = initial_energy; + double vtxc = initial_vtxc; + elecstate::PotNonlocalXC host(&pw, &energy, &vtxc, "rvv10", 6.3, 0.0093); + host.cal_v_eff(&charge, &cell, potential); + + const double expected_energy = initial_energy + correction.energy; + const double expected_vtxc = initial_vtxc + correction.vtxc; + EXPECT_NEAR(energy, expected_energy, 2.0e-12); + EXPECT_NEAR(vtxc, expected_vtxc, 2.0e-12); + for (int r = 0; r < pw.nrxx; ++r) + { + const double expected_potential = initial_potential[r] + correction.potential[r]; + EXPECT_NEAR(potential(0, r), expected_potential, 2.0e-12); + EXPECT_DOUBLE_EQ(charge.rho[0][r], valence[r]); + EXPECT_DOUBLE_EQ(charge.rho_core[r], core[r]); + } +} diff --git a/source/source_hamilt/module_vdw/CMakeLists.txt b/source/source_hamilt/module_vdw/CMakeLists.txt index e51299da04..dc5f128683 100644 --- a/source/source_hamilt/module_vdw/CMakeLists.txt +++ b/source/source_hamilt/module_vdw/CMakeLists.txt @@ -15,6 +15,17 @@ endif() add_library(vdw OBJECT ${vdw_sources}) +# Density-grid nonlocal XC has a different dependency and data model from the +# atom-pair D2/D3/D4 correction. Keep it in a separate object target so the +# existing pairwise-vdW unit tests do not acquire PW/XC link dependencies. +set(nonlocal_xc_sources nonlocal_xc.cpp) + +if(ENABLE_LIBVDWXC) + list(APPEND nonlocal_xc_sources nonlocal_xc_libvdwxc.cpp) +endif() + +add_library(nonlocal_xc OBJECT ${nonlocal_xc_sources}) + if(ENABLE_DFTD4) target_link_libraries(vdw PUBLIC dftd4::dftd4) endif() @@ -22,9 +33,3 @@ endif() if(ENABLE_COVERAGE) add_coverage(vdw) endif() - -if(BUILD_TESTING) - if(ENABLE_MPI) - add_subdirectory(test) - endif() -endif() diff --git a/source/source_hamilt/module_vdw/nonlocal_xc.cpp b/source/source_hamilt/module_vdw/nonlocal_xc.cpp new file mode 100644 index 0000000000..8ae9f54447 --- /dev/null +++ b/source/source_hamilt/module_vdw/nonlocal_xc.cpp @@ -0,0 +1,33 @@ +#include "nonlocal_xc.h" +#ifdef __LIBVDWXC +#include "nonlocal_xc_libvdwxc.h" +#endif + +#include "source_base/tool_quit.h" + +namespace vdw +{ + +std::unique_ptr make_nonlocal_xc(const ModulePW::PW_Basis* pw, + const std::string& method, + double b, + double c) +{ + if (method == "rvv10") + { +#ifdef __LIBVDWXC + if (pw == nullptr) + { + ModuleBase::WARNING_QUIT("NonlocalXC", "libvdwxc rVV10 requires a density PW basis"); + } + return std::unique_ptr(new LibvdwxcNonlocalXC(*pw, b, c)); +#else + ModuleBase::WARNING_QUIT("NonlocalXC", "xc_nonlocal=rvv10 requires a CMake build with ENABLE_LIBVDWXC=ON"); +#endif + } + + ModuleBase::WARNING_QUIT("NonlocalXC", "Unsupported nonlocal XC method: " + method); + return nullptr; +} + +} // namespace vdw diff --git a/source/source_hamilt/module_vdw/nonlocal_xc.h b/source/source_hamilt/module_vdw/nonlocal_xc.h new file mode 100644 index 0000000000..4941edfbc1 --- /dev/null +++ b/source/source_hamilt/module_vdw/nonlocal_xc.h @@ -0,0 +1,46 @@ +#ifndef ABACUS_NONLOCAL_XC_H +#define ABACUS_NONLOCAL_XC_H + +#include +#include +#include + +namespace ModulePW +{ +class PW_Basis; +} + +namespace vdw +{ + +struct NonlocalXCResult +{ + double energy; + double vtxc; + std::vector potential; +}; + +/** + * @brief Density-grid interface for a nonlocal exchange-correlation backend. + * + * The host potential code only knows this interface. Concrete methods such as + * rVV10, and future libvdwxc-backed methods, remain in the vdW module. + */ +class NonlocalXC +{ + public: + virtual ~NonlocalXC() = default; + + virtual NonlocalXCResult evaluate(const ModulePW::PW_Basis& pw, + const std::vector& total_density, + const std::vector& valence_density) const = 0; +}; + +std::unique_ptr make_nonlocal_xc(const ModulePW::PW_Basis* pw, + const std::string& method, + double b, + double c); + +} // namespace vdw + +#endif // ABACUS_NONLOCAL_XC_H diff --git a/source/source_hamilt/module_vdw/nonlocal_xc_libvdwxc.cpp b/source/source_hamilt/module_vdw/nonlocal_xc_libvdwxc.cpp new file mode 100644 index 0000000000..f7529c0733 --- /dev/null +++ b/source/source_hamilt/module_vdw/nonlocal_xc_libvdwxc.cpp @@ -0,0 +1,220 @@ +#include "nonlocal_xc_libvdwxc.h" + +#ifdef __LIBVDWXC + +#include "source_basis/module_pw/pw_basis.h" +#include "source_base/tool_quit.h" +#include "source_hamilt/module_xc/xc_functional.h" + +#include +#include +#include + +namespace vdw +{ +namespace +{ +void check_serial_basis(const ModulePW::PW_Basis& pw) +{ + const double determinant = pw.latvec.Det(); + if (pw.gamma_only || pw.get_device() != "cpu" || pw.get_precision() != "double" || pw.poolnproc != 1 + || pw.nrxx != pw.nxyz || pw.nx <= 0 || pw.ny <= 0 || pw.nz <= 0 || pw.nxyz <= 0 || pw.npw <= 0 + || !std::isfinite(pw.omega) || pw.omega <= 0.0 || !std::isfinite(determinant) || determinant <= 0.0 + || !std::isfinite(pw.tpiba) || pw.tpiba <= 0.0) + { + ModuleBase::WARNING_QUIT("LibvdwxcNonlocalXC", + "libvdwxc rVV10 requires one-process CPU/double full-complex PW density grid and a right-handed cell"); + } +} + +std::array cell_vectors(const ModulePW::PW_Basis& pw) +{ + const auto& a = pw.latvec; + // ABACUS stores a1, a2, a3 as rows. The C API groups C00,C10,C20 as + // one vector, so pass the ABACUS rows consecutively in those groups. + return {a.e11 * pw.lat0, + a.e12 * pw.lat0, + a.e13 * pw.lat0, + a.e21 * pw.lat0, + a.e22 * pw.lat0, + a.e23 * pw.lat0, + a.e31 * pw.lat0, + a.e32 * pw.lat0, + a.e33 * pw.lat0}; +} + +bool same_cell(const std::array& lhs, const std::array& rhs) +{ + for (std::size_t i = 0; i < lhs.size(); ++i) + { + const double lhs_abs = std::abs(lhs[i]); + const double rhs_abs = std::abs(rhs[i]); + double scale = 1.0; + if (lhs_abs > scale) + { + scale = lhs_abs; + } + if (rhs_abs > scale) + { + scale = rhs_abs; + } + if (std::abs(lhs[i] - rhs[i]) > 1.0e-12 * scale) + { + return false; + } + } + return true; +} +} // namespace + +LibvdwxcNonlocalXC::LibvdwxcNonlocalXC(const ModulePW::PW_Basis& pw, double b, double c) +{ + check_serial_basis(pw); + data_ = vdwxc_new(VDWXC_RVV10); + if (data_ == nullptr) + { + ModuleBase::WARNING_QUIT("LibvdwxcNonlocalXC", "libvdwxc failed to create an rVV10 evaluator"); + } + if (vdwxc_set_rvv10_parameters(data_, b, c) != 0) + { + vdwxc_finalize(&data_); + ModuleBase::WARNING_QUIT("LibvdwxcNonlocalXC", "libvdwxc rejected the rVV10 b/C parameters"); + } + + cell_ = cell_vectors(pw); + // libvdwxc names the arguments by matrix columns (C00,C10,C20,...). + vdwxc_set_unit_cell(data_, + pw.nx, + pw.ny, + pw.nz, + cell_[0], + cell_[1], + cell_[2], + cell_[3], + cell_[4], + cell_[5], + cell_[6], + cell_[7], + cell_[8]); + vdwxc_init_serial(data_); + nx_ = pw.nx; + ny_ = pw.ny; + nz_ = pw.nz; + omega_ = pw.omega; +} + +LibvdwxcNonlocalXC::~LibvdwxcNonlocalXC() +{ + if (data_ != nullptr) + { + vdwxc_finalize(&data_); + } +} + +NonlocalXCResult LibvdwxcNonlocalXC::evaluate(const ModulePW::PW_Basis& pw, + const std::vector& total_density, + const std::vector& valence_density) const +{ + check_serial_basis(pw); + const std::array current_cell = cell_vectors(pw); + const double omega_abs = std::abs(pw.omega); + const double initialized_omega_abs = std::abs(omega_); + double volume_scale = 1.0; + if (omega_abs > volume_scale) + { + volume_scale = omega_abs; + } + if (initialized_omega_abs > volume_scale) + { + volume_scale = initialized_omega_abs; + } + const bool same_grid = pw.nx == nx_ && pw.ny == ny_ && pw.nz == nz_; + const bool same_volume = std::abs(pw.omega - omega_) <= 1.0e-12 * volume_scale; + if (!same_grid || !same_cell(cell_, current_cell) || !same_volume + || total_density.size() != static_cast(pw.nrxx) + || valence_density.size() != total_density.size()) + { + ModuleBase::WARNING_QUIT("LibvdwxcNonlocalXC", + "libvdwxc rVV10 received a density grid different from its initialized cell"); + } + + using Complex = std::complex; + std::vector density_g(pw.npw); + std::vector> gradient(pw.nrxx); + std::vector sigma(pw.nrxx, 0.0); + std::vector dedn(pw.nrxx, 0.0); + std::vector dedsigma(pw.nrxx, 0.0); + std::vector> field(pw.nrxx); + std::vector divergence(pw.nrxx, 0.0); + + for (double value: total_density) + { + if (!std::isfinite(value)) + { + ModuleBase::WARNING_QUIT("LibvdwxcNonlocalXC", "libvdwxc rVV10 received a nonfinite total density"); + } + } + for (double value: valence_density) + { + if (!std::isfinite(value)) + { + ModuleBase::WARNING_QUIT("LibvdwxcNonlocalXC", "libvdwxc rVV10 received a nonfinite valence density"); + } + } + + pw.real2recip(total_density.data(), density_g.data()); + XC_Functional::grad_rho(density_g.data(), gradient.data(), &pw, pw.tpiba); + for (int r = 0; r < pw.nrxx; ++r) + { + sigma[r] = gradient[r].x * gradient[r].x + gradient[r].y * gradient[r].y + gradient[r].z * gradient[r].z; + field[r] = ModuleBase::Vector3(0.0, 0.0, 0.0); + } + + std::vector density = total_density; + const double energy_ha = vdwxc_calculate(data_, + density.data(), + sigma.data(), + dedn.data(), + dedsigma.data()); + if (!std::isfinite(energy_ha)) + { + ModuleBase::WARNING_QUIT("LibvdwxcNonlocalXC", "libvdwxc returned a nonfinite rVV10 energy"); + } + + for (int r = 0; r < pw.nrxx; ++r) + { + if (!std::isfinite(dedn[r]) || !std::isfinite(dedsigma[r])) + { + ModuleBase::WARNING_QUIT("LibvdwxcNonlocalXC", "libvdwxc returned a nonfinite rVV10 derivative"); + } + const double coefficient = 2.0 * dedsigma[r]; + field[r].x = coefficient * gradient[r].x; + field[r].y = coefficient * gradient[r].y; + field[r].z = coefficient * gradient[r].z; + } + XC_Functional::grad_dot(field.data(), divergence.data(), &pw, pw.tpiba); + + const double energy_ry = 2.0 * energy_ha; + std::vector potential(pw.nrxx, 0.0); + NonlocalXCResult result{energy_ry, 0.0, potential}; + const double dv = pw.omega / pw.nxyz; + for (int r = 0; r < pw.nrxx; ++r) + { + // libvdwxc returns Hartree derivatives; ABACUS stores XC quantities in Ry. + result.potential[r] = 2.0 * (dedn[r] - divergence[r]); + if (!std::isfinite(divergence[r]) || !std::isfinite(result.potential[r])) + { + ModuleBase::WARNING_QUIT("LibvdwxcNonlocalXC", "libvdwxc produced a nonfinite rVV10 potential"); + } + result.vtxc += dv * valence_density[r] * result.potential[r]; + } + if (!std::isfinite(result.vtxc)) + { + ModuleBase::WARNING_QUIT("LibvdwxcNonlocalXC", "libvdwxc produced a nonfinite rVV10 vtxc"); + } + return result; +} + +} // namespace vdw + +#endif // __LIBVDWXC diff --git a/source/source_hamilt/module_vdw/nonlocal_xc_libvdwxc.h b/source/source_hamilt/module_vdw/nonlocal_xc_libvdwxc.h new file mode 100644 index 0000000000..4b460c4765 --- /dev/null +++ b/source/source_hamilt/module_vdw/nonlocal_xc_libvdwxc.h @@ -0,0 +1,57 @@ +#ifndef ABACUS_NONLOCAL_XC_LIBVDWXC_H +#define ABACUS_NONLOCAL_XC_LIBVDWXC_H + +#include "nonlocal_xc.h" + +#include + +#ifdef __LIBVDWXC +// libvdwxc's public C header includes inside an extern "C" +// block. Pre-include the standard header here so libc++ does not parse its +// C++ templates with C linkage. +#include +#include +// complex.h defines I as a macro; it collides with template identifiers in +// GoogleTest and other C++ headers, but is not needed after vdwxc.h parsing. +#ifdef I +#undef I +#endif +#endif + +namespace vdw +{ + +#ifdef __LIBVDWXC +/** + * @brief Thin serial/PW adapter for libvdwxc's rVV10 evaluator. + * + * libvdwxc owns the nonlocal kernel, interpolation and convolution. ABACUS + * owns density construction, FFT gradients/divergence, units and XC + * bookkeeping. The first adapter intentionally rejects distributed layouts; + * the libvdwxc MPI grid/communicator contract must be validated separately. + */ +class LibvdwxcNonlocalXC final : public NonlocalXC +{ + public: + LibvdwxcNonlocalXC(const ModulePW::PW_Basis& pw, double b, double c); + ~LibvdwxcNonlocalXC() override; + LibvdwxcNonlocalXC(const LibvdwxcNonlocalXC&) = delete; + LibvdwxcNonlocalXC& operator=(const LibvdwxcNonlocalXC&) = delete; + + NonlocalXCResult evaluate(const ModulePW::PW_Basis& pw, + const std::vector& total_density, + const std::vector& valence_density) const override; + + private: + vdwxc_data data_ = nullptr; + int nx_ = 0; + int ny_ = 0; + int nz_ = 0; + double omega_ = 0.0; + std::array cell_{}; +}; +#endif + +} // namespace vdw + +#endif // ABACUS_NONLOCAL_XC_LIBVDWXC_H diff --git a/source/source_hamilt/module_vdw/test/CMakeLists.txt b/source/source_hamilt/module_vdw/test/CMakeLists.txt index 195a5312ce..c1d2d56517 100644 --- a/source/source_hamilt/module_vdw/test/CMakeLists.txt +++ b/source/source_hamilt/module_vdw/test/CMakeLists.txt @@ -6,11 +6,13 @@ install(FILES c6.txt DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) install(FILES r0.txt DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) install(FILES r0_zero.txt DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) -AddTest( - TARGET MODULE_HAMILT_vdwTest - LIBS parameter base device vdw - SOURCES vdw_test.cpp vdwd3_evaluator_test.cpp -) +if(ENABLE_MPI) + AddTest( + TARGET MODULE_HAMILT_vdwTest + LIBS parameter base device vdw + SOURCES vdw_test.cpp vdwd3_evaluator_test.cpp + ) +endif() if(ENABLE_DFTD4 AND ENABLE_MPI) add_test( @@ -21,3 +23,19 @@ if(ENABLE_DFTD4 AND ENABLE_MPI) WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} ) endif() +if(ENABLE_LIBVDWXC AND NOT ENABLE_MPI AND NOT USE_CUDA AND NOT USE_ROCM) + AddTest( + TARGET MODULE_HAMILT_vdwTest_LIBVDWXC + LIBS xc_ planewave base parameter device nonlocal_xc + SOURCES test_nonlocal_xc_libvdwxc.cpp + ) + if(ENABLE_LIBXC) + target_link_libraries(MODULE_HAMILT_vdwTest_LIBVDWXC PRIVATE Libxc::xc) + endif() + target_link_libraries(MODULE_HAMILT_vdwTest_LIBVDWXC PRIVATE libvdwxc::vdwxc) + set_tests_properties(MODULE_HAMILT_vdwTest_LIBVDWXC PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=1") +endif() + +if(ENABLE_DFTD4 AND ENABLE_MPI) + target_link_libraries(MODULE_HAMILT_vdwTest PRIVATE dftd4::dftd4) +endif() diff --git a/source/source_hamilt/module_vdw/test/test_nonlocal_xc_libvdwxc.cpp b/source/source_hamilt/module_vdw/test/test_nonlocal_xc_libvdwxc.cpp new file mode 100644 index 0000000000..ea27dbc132 --- /dev/null +++ b/source/source_hamilt/module_vdw/test/test_nonlocal_xc_libvdwxc.cpp @@ -0,0 +1,240 @@ +#include "../nonlocal_xc_libvdwxc.h" + +#include "source_basis/module_pw/pw_basis.h" +#include "source_hamilt/module_xc/xc_functional.h" + +#include +#include +#include +#include +#include +#include + +namespace +{ +const double pi = std::acos(-1.0); + +class LibvdwxcAdapterTest : public testing::Test +{ + protected: + ModulePW::PW_Basis pw{"cpu", "double"}; + + void SetUp() override + { + // The oblique cell catches row/column-major cell mapping mistakes. + const ModuleBase::Matrix3 lattice(1.0, 0.0, 0.0, -0.5, std::sqrt(0.75), 0.0, 0.0, 0.0, 2.0); + pw.initgrids(10.0, lattice, 6, 8, 12); + pw.initparameters(false, 4.0, 1, true); + pw.setuptransform(); + pw.collect_local_pw(); + } + + std::vector density(double offset, double amplitude) const + { + std::vector values(pw.nrxx); + for (int x = 0; x < pw.nx; ++x) + for (int y = 0; y < pw.ny; ++y) + for (int z = 0; z < pw.nz; ++z) + { + const int index = (x * pw.ny + y) * pw.nz + z; + const double ax = 2.0 * pi * x / pw.nx; + const double ay = 2.0 * pi * y / pw.ny; + const double az = 2.0 * pi * z / pw.nz; + values[index] = offset + amplitude * (1.0 + 0.25 * std::sin(ax) + 0.15 * std::cos(ay + az)); + } + return values; + } + + std::vector perturbation() const + { + std::vector values(pw.nrxx, 0.0); + for (int x = 0; x < pw.nx; ++x) + for (int y = 0; y < pw.ny; ++y) + for (int z = 0; z < pw.nz; ++z) + { + const int index = (x * pw.ny + y) * pw.nz + z; + const double ax = 2.0 * pi * x / pw.nx; + const double ay = 2.0 * pi * y / pw.ny; + values[index] = std::sin(ax + ay) + 0.3 * std::cos(2.0 * ay); + } + return values; + } +}; + +#ifdef __LIBVDWXC +TEST(LibvdwxcOwnershipTest, BackendHandleCannotBeCopied) +{ + EXPECT_FALSE(std::is_copy_constructible::value); + EXPECT_FALSE(std::is_copy_assignable::value); +} + +TEST_F(LibvdwxcAdapterTest, MatchesDirectSerialApiForNonuniformDensity) +{ + const double b = 6.3; + const double c = 0.0093; + const std::vector total = density(0.004, 0.006); + const std::vector valence = density(0.002, 0.003); + vdw::LibvdwxcNonlocalXC adapter(pw, b, c); + const vdw::NonlocalXCResult adapted = adapter.evaluate(pw, total, valence); + + vdwxc_data direct = vdwxc_new(VDWXC_RVV10); + ASSERT_NE(direct, nullptr); + ASSERT_EQ(vdwxc_set_rvv10_parameters(direct, b, c), 0); + const auto& a = pw.latvec; + vdwxc_set_unit_cell(direct, pw.nx, pw.ny, pw.nz, + // libvdwxc stores each C00,C10,C20 group as one + // lattice vector; ABACUS stores lattice vectors in rows. + a.e11 * pw.lat0, a.e12 * pw.lat0, a.e13 * pw.lat0, + a.e21 * pw.lat0, a.e22 * pw.lat0, a.e23 * pw.lat0, + a.e31 * pw.lat0, a.e32 * pw.lat0, a.e33 * pw.lat0); + vdwxc_init_serial(direct); + + std::vector> total_g(pw.npw); + std::vector> gradient(pw.nrxx); + std::vector sigma(pw.nrxx, 0.0); + std::vector dedn(pw.nrxx, 0.0); + std::vector dedsigma(pw.nrxx, 0.0); + std::vector> field(pw.nrxx); + std::vector divergence(pw.nrxx, 0.0); + pw.real2recip(total.data(), total_g.data()); + XC_Functional::grad_rho(total_g.data(), gradient.data(), &pw, pw.tpiba); + for (int r = 0; r < pw.nrxx; ++r) + { + sigma[r] = gradient[r].x * gradient[r].x + gradient[r].y * gradient[r].y + gradient[r].z * gradient[r].z; + } + std::vector direct_density = total; + const double direct_energy = vdwxc_calculate(direct, direct_density.data(), sigma.data(), dedn.data(), dedsigma.data()); + for (int r = 0; r < pw.nrxx; ++r) + { + const double coefficient = 2.0 * dedsigma[r]; + field[r].x = coefficient * gradient[r].x; + field[r].y = coefficient * gradient[r].y; + field[r].z = coefficient * gradient[r].z; + } + XC_Functional::grad_dot(field.data(), divergence.data(), &pw, pw.tpiba); + + const double dv = pw.omega / pw.nxyz; + double expected_vtxc = 0.0; + EXPECT_NEAR(adapted.energy, 2.0 * direct_energy, 2.0e-12); + for (int r = 0; r < pw.nrxx; ++r) + { + const double expected_potential = 2.0 * (dedn[r] - divergence[r]); + expected_vtxc += dv * valence[r] * expected_potential; + EXPECT_NEAR(adapted.potential[r], expected_potential, 2.0e-12); + } + EXPECT_NEAR(adapted.vtxc, expected_vtxc, 2.0e-12); + vdwxc_finalize(&direct); +} + +TEST_F(LibvdwxcAdapterTest, RejectsLeftHandedCell) +{ + ModulePW::PW_Basis reflected("cpu", "double"); + const ModuleBase::Matrix3 lattice(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, -2.0); + reflected.initgrids(10.0, lattice, 6, 8, 12); + reflected.initparameters(false, 4.0, 1, true); + reflected.setuptransform(); + reflected.collect_local_pw(); + EXPECT_EXIT(vdw::LibvdwxcNonlocalXC(reflected, 6.3, 0.0093), testing::ExitedWithCode(1), ""); +} + +TEST_F(LibvdwxcAdapterTest, RejectsInvalidParametersWithStandardDiagnostic) +{ + EXPECT_EXIT(vdw::LibvdwxcNonlocalXC(pw, 0.0, 0.0093), testing::ExitedWithCode(1), ""); +} + +TEST_F(LibvdwxcAdapterTest, RejectsInvalidDensityWithStandardDiagnostic) +{ + vdw::LibvdwxcNonlocalXC adapter(pw, 6.3, 0.0093); + const std::vector total = density(0.004, 0.006); + std::vector invalid = total; + invalid.pop_back(); + EXPECT_EXIT(adapter.evaluate(pw, invalid, total), testing::ExitedWithCode(1), ""); + invalid = total; + invalid[0] = std::numeric_limits::quiet_NaN(); + EXPECT_EXIT(adapter.evaluate(pw, invalid, total), testing::ExitedWithCode(1), ""); +} + +TEST_F(LibvdwxcAdapterTest, CellMappingPreservesIndependentReciprocalGeometry) +{ + // Both cells have volume sqrt(3)*lat0^3 and |G100|=4*pi/(sqrt(3)*lat0). + // A density depending only on fractional x must therefore give identical + // isotropic nonlocal energy; transposing the oblique cell breaks this equality. + ModulePW::PW_Basis orthogonal("cpu", "double"); + const ModuleBase::Matrix3 lattice(std::sqrt(0.75), 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 2.0); + orthogonal.initgrids(10.0, lattice, pw.nx, pw.ny, pw.nz); + orthogonal.initparameters(false, 4.0, 1, true); + orthogonal.setuptransform(); + orthogonal.collect_local_pw(); + + const double amplitude = 0.003; + std::vector total(pw.nrxx); + for (int x = 0; x < pw.nx; ++x) + for (int y = 0; y < pw.ny; ++y) + for (int z = 0; z < pw.nz; ++z) + { + const int index = (x * pw.ny + y) * pw.nz + z; + const double phase = 2.0 * pi * x / pw.nx; + total[index] = 0.01 + amplitude * std::cos(phase); + } + + // The hand-derived gradient of cos(2*pi*s1) uses b1=(1,1/sqrt(3),0). + std::vector> density_g(pw.npw); + std::vector> gradient(pw.nrxx); + pw.real2recip(total.data(), density_g.data()); + XC_Functional::grad_rho(density_g.data(), gradient.data(), &pw, pw.tpiba); + for (int x = 0; x < pw.nx; ++x) + for (int y = 0; y < pw.ny; ++y) + for (int z = 0; z < pw.nz; ++z) + { + const int index = (x * pw.ny + y) * pw.nz + z; + const double phase = 2.0 * pi * x / pw.nx; + const double derivative = -amplitude * (2.0 * pi / pw.lat0) * std::sin(phase); + EXPECT_NEAR(gradient[index].x, derivative, 2.0e-12); + EXPECT_NEAR(gradient[index].y, derivative / std::sqrt(3.0), 2.0e-12); + EXPECT_NEAR(gradient[index].z, 0.0, 2.0e-12); + } + + vdw::LibvdwxcNonlocalXC oblique_adapter(pw, 6.3, 0.0093); + vdw::LibvdwxcNonlocalXC orthogonal_adapter(orthogonal, 6.3, 0.0093); + const vdw::NonlocalXCResult oblique_result = oblique_adapter.evaluate(pw, total, total); + const vdw::NonlocalXCResult orthogonal_result = orthogonal_adapter.evaluate(orthogonal, total, total); + EXPECT_NEAR(pw.omega, orthogonal.omega, 2.0e-12); + ASSERT_NEAR(oblique_result.energy, orthogonal_result.energy, 2.0e-12); + EXPECT_NEAR(oblique_result.vtxc, orthogonal_result.vtxc, 2.0e-12); + for (int r = 0; r < pw.nrxx; ++r) + { + EXPECT_NEAR(oblique_result.potential[r], orthogonal_result.potential[r], 2.0e-12); + } +} + +TEST_F(LibvdwxcAdapterTest, PotentialMatchesFiniteDifferenceOfEnergy) +{ + const double b = 6.3; + const double c = 0.0093; + const double step = 1.0e-5; + const std::vector total = density(0.004, 0.006); + const std::vector valence = total; + const std::vector direction = perturbation(); + vdw::LibvdwxcNonlocalXC adapter(pw, b, c); + const vdw::NonlocalXCResult reference = adapter.evaluate(pw, total, valence); + std::vector plus = total; + std::vector minus = total; + for (int r = 0; r < pw.nrxx; ++r) + { + plus[r] += step * direction[r]; + minus[r] -= step * direction[r]; + } + const vdw::NonlocalXCResult forward = adapter.evaluate(pw, plus, valence); + const vdw::NonlocalXCResult backward = adapter.evaluate(pw, minus, valence); + const double finite_difference = (forward.energy - backward.energy) / (2.0 * step); + const double dv = pw.omega / pw.nxyz; + double potential_contraction = 0.0; + for (int r = 0; r < pw.nrxx; ++r) + { + potential_contraction += dv * reference.potential[r] * direction[r]; + } + EXPECT_NEAR(finite_difference, potential_contraction, 2.0e-6); +} +#endif + +} // namespace diff --git a/source/source_io/build_info.h.in b/source/source_io/build_info.h.in index 543ef3e1b4..aa9ba4c1f7 100644 --- a/source/source_io/build_info.h.in +++ b/source/source_io/build_info.h.in @@ -31,6 +31,7 @@ #define ABACUS_ELPA_VERSION "@ABACUS_ELPA_VERSION@" #define ABACUS_MKL_SUPPORT "@ABACUS_MKL_SUPPORT@" #define ABACUS_LIBXC_VERSION "@ABACUS_LIBXC_VERSION@" +#define ABACUS_LIBVDWXC_VERSION "@ABACUS_LIBVDWXC_VERSION@" #define ABACUS_FFTW_VERSION "@ABACUS_FFTW_VERSION@" // --- Accelerators & Specific Hardware --- diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index d86993b63d..07cd7bcbd1 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -88,6 +88,9 @@ struct Input_para double nupdown = 0.0; std::string dft_functional = "default"; ///< input DFT functional. int gga_grad = 0; ///< Noncollinear GGA: 0 original, 1 local axis, 2 regularized mapped-density LCA. + std::string xc_nonlocal = "none"; ///< nonlocal XC correction (none; rvv10) + double rvv10_b = 6.3; ///< rVV10 damping parameter b + double rvv10_c = 0.0093; ///< rVV10 parameter C double xc_temperature = 0.0; ///< only relevant if finite temperature functional is used double pseudo_rcut = 15.0; ///< cut-off radius for calculating msh bool pseudo_mesh = false; ///< 0: use msh to normalize radial wave functions; 1: diff --git a/source/source_io/module_parameter/read_inp_estruc.cpp b/source/source_io/module_parameter/read_inp_estruc.cpp index b5841a81e3..2762b3e02a 100644 --- a/source/source_io/module_parameter/read_inp_estruc.cpp +++ b/source/source_io/module_parameter/read_inp_estruc.cpp @@ -3,6 +3,10 @@ #include "read_input.h" #include "read_input_tool.h" +#include +#include +#include + namespace ModuleIO { void ReadInput::item_elec_stru() @@ -333,6 +337,107 @@ The other way is only available when compiling with LIBXC, and it allows for sup item.default_value = "Used the same as DFT functional as specified in the pseudopotential files."; item.unit = ""; read_sync_string(input.dft_functional); + item.check_value = [](const Input_Item&, const Parameter& para) { + std::string name = para.input.dft_functional; + std::transform(name.begin(), name.end(), name.begin(), [](unsigned char c) { + return static_cast(std::toupper(c)); + }); + if (name == "RVV10") + ModuleBase::WARNING_QUIT( + "ReadInput", "Use xc_nonlocal=rvv10 with an explicit semilocal dft_functional"); + }; + this->add_item(item); + } + { + Input_Item item("xc_nonlocal"); + item.annotation = "none; rvv10"; + item.category = "Electronic structure"; + item.type = "String"; + item.description = R"(Select an optional libvdwxc-backed nonlocal exchange-correlation correction. `none` disables the correction. The current rVV10 interface requires a CMake build with `ENABLE_LIBXC=ON` and `ENABLE_LIBVDWXC=ON`, `dft_functional=GGA_X_RPW86+GGA_C_PBE`, fixed b=6.3 and C=0.0093, and a single-process, CPU double-precision, plane-wave, single-density-grid, right-handed-cell, norm-conserving, spin-unpolarized SCF calculation. Forces, stress, geometry optimization, molecular dynamics, LCAO, GPU, reduced-gamma layouts, ultrasoft/PAW pseudopotentials, spin-polarized calculations and pairwise vdW corrections are not part of this interface.)"; + item.default_value = "none"; + item.unit = ""; + read_sync_string(input.xc_nonlocal); + item.reset_value = [](const Input_Item&, Parameter& para) { + std::transform(para.input.xc_nonlocal.begin(), para.input.xc_nonlocal.end(), para.input.xc_nonlocal.begin(), + [](unsigned char c) { return static_cast(std::tolower(c)); }); + }; + item.check_value = [](const Input_Item&, const Parameter& para) { + std::string nonlocal = para.input.xc_nonlocal; + std::transform(nonlocal.begin(), nonlocal.end(), nonlocal.begin(), [](unsigned char c) { + return static_cast(std::tolower(c)); + }); + if (nonlocal != "none" && nonlocal != "rvv10") + ModuleBase::WARNING_QUIT("ReadInput", "xc_nonlocal must be either none or rvv10"); + if (nonlocal != "rvv10") + return; +#if !defined(__LIBXC) || !defined(__LIBVDWXC) + ModuleBase::WARNING_QUIT("ReadInput", "xc_nonlocal=rvv10 requires Libxc and libvdwxc-enabled CMake builds"); +#endif + if (!std::isfinite(para.input.rvv10_b) || para.input.rvv10_b <= 0.0 || !std::isfinite(para.input.rvv10_c) + || para.input.rvv10_c < 0.0) + { + ModuleBase::WARNING_QUIT( + "ReadInput", "rvv10_b must be finite and positive, and rvv10_c must be finite and non-negative"); + } + std::string functional = para.input.dft_functional; + std::transform(functional.begin(), functional.end(), functional.begin(), [](unsigned char c) { + return static_cast(std::toupper(c)); + }); + functional.erase(std::remove_if(functional.begin(), functional.end(), [](unsigned char c) { + return std::isspace(c); + }), functional.end()); + if (functional.empty() || functional == "DEFAULT") + ModuleBase::WARNING_QUIT( + "ReadInput", "xc_nonlocal=rvv10 requires an explicit semilocal dft_functional"); + if (functional != "GGA_X_RPW86+GGA_C_PBE") + ModuleBase::WARNING_QUIT( + "ReadInput", + "xc_nonlocal=rvv10 requires dft_functional=GGA_X_RPW86+GGA_C_PBE for the libvdwxc backend"); + if (std::abs(para.input.rvv10_b - 6.3) > 1.0e-12 || std::abs(para.input.rvv10_c - 0.0093) > 1.0e-12) + ModuleBase::WARNING_QUIT( + "ReadInput", "xc_nonlocal=rvv10 uses fixed parameters rvv10_b=6.3 and rvv10_c=0.0093"); + const auto& p = para.input; + if (para.globalv.nproc != 1 || p.calculation != "scf" || p.esolver_type != "ksdft" || p.basis_type != "pw" || p.device != "cpu" + || p.precision != "double" || p.nspin != 1 || p.cal_force || p.cal_stress || para.sys.double_grid + || p.gamma_only || p.vdw_method != "none" || !p.vl_in_h) + ModuleBase::WARNING_QUIT("ReadInput", + "xc_nonlocal=rvv10 supports only one MPI rank with CPU double PW nspin=1 " + "serial SCF with a single density grid, without forces, stress, gamma_only or other vdW corrections"); + }; + this->add_item(item); + } + { + Input_Item item("rvv10_b"); + item.annotation = "positive real"; + item.category = "Electronic structure"; + item.type = "Real"; + item.description = "Fixed damping parameter b=6.3 for the current libvdwxc-backed `GGA_X_RPW86+GGA_C_PBE` rVV10 interface."; + item.default_value = "6.3"; + item.unit = "dimensionless"; + item.set_availability("xc_nonlocal==rvv10"); + read_sync_double(input.rvv10_b); + item.check_value = [](const Input_Item&, const Parameter& para) { + if (para.input.xc_nonlocal == "rvv10" + && (!std::isfinite(para.input.rvv10_b) || para.input.rvv10_b <= 0.0)) + ModuleBase::WARNING_QUIT("ReadInput", "rvv10_b must be finite and positive"); + }; + this->add_item(item); + } + { + Input_Item item("rvv10_c"); + item.annotation = "non-negative real"; + item.category = "Electronic structure"; + item.type = "Real"; + item.description = "Fixed parameter C=0.0093 for the current libvdwxc-backed `GGA_X_RPW86+GGA_C_PBE` rVV10 interface."; + item.default_value = "0.0093"; + item.unit = "dimensionless"; + item.set_availability("xc_nonlocal==rvv10"); + read_sync_double(input.rvv10_c); + item.check_value = [](const Input_Item&, const Parameter& para) { + if (para.input.xc_nonlocal == "rvv10" + && (!std::isfinite(para.input.rvv10_c) || para.input.rvv10_c < 0.0)) + ModuleBase::WARNING_QUIT("ReadInput", "rvv10_c must be finite and non-negative"); + }; this->add_item(item); } { @@ -805,6 +910,18 @@ Note: If gamma_only is set to 1, the KPT file will be overwritten. So make sure item.reset_value = [](const Input_Item& item, Parameter& para) { if (para.input.basis_type == "pw" && para.input.gamma_only) { + // Resets precede validation. Reject before the PW fallback + // overwrites KPT; xc_nonlocal may not yet be lowercased. + std::string nonlocal = para.input.xc_nonlocal; + std::transform(nonlocal.begin(), nonlocal.end(), nonlocal.begin(), [](unsigned char c) { + return static_cast(std::tolower(c)); + }); + if (nonlocal == "rvv10") + { + ModuleBase::WARNING_QUIT("ReadInput", + "xc_nonlocal=rvv10 supports only one MPI rank with CPU double PW nspin=1 " + "serial SCF without forces, stress, gamma_only or other vdW corrections"); + } para.input.gamma_only = false; GlobalV::ofs_warning << " WARNING : gamma_only has not been implemented for pw yet" << std::endl; GlobalV::ofs_warning << "gamma_only is not supported in the pw model" << std::endl; diff --git a/source/source_io/parse_args.cpp b/source/source_io/parse_args.cpp index 046bc52702..41cf5054c6 100644 --- a/source/source_io/parse_args.cpp +++ b/source/source_io/parse_args.cpp @@ -140,6 +140,7 @@ void print_build_info() print_info("ELPA Support", ABACUS_ELPA_VERSION); print_info("MKL Support", ABACUS_MKL_SUPPORT); print_info("LibXC Support", ABACUS_LIBXC_VERSION); + print_info("libvdwxc Support", ABACUS_LIBVDWXC_VERSION); print_info("FFTW Support", ABACUS_FFTW_VERSION); // --- 8. Accelerators & Specific Hardware --- diff --git a/source/source_pw/module_pwdft/Makefile b/source/source_pw/module_pwdft/Makefile index ad634ae49e..8091a70670 100644 --- a/source/source_pw/module_pwdft/Makefile +++ b/source/source_pw/module_pwdft/Makefile @@ -139,4 +139,4 @@ ${OBJ_DIR}/%.o:%.cu .PHONY:clean clean: @ if [ -d $(OBJ_DIR) ]; then rm -rf $(OBJ_DIR); fi - @ if [ -e ${VERSION}-PW.x ]; then rm -f ${VERSION}-PW.x; fi \ No newline at end of file + @ if [ -e ${VERSION}-PW.x ]; then rm -f ${VERSION}-PW.x; fi diff --git a/source/source_pw/module_pwdft/hamilt_pw.cpp b/source/source_pw/module_pwdft/hamilt_pw.cpp index 13977cd9f6..fdc189c464 100644 --- a/source/source_pw/module_pwdft/hamilt_pw.cpp +++ b/source/source_pw/module_pwdft/hamilt_pw.cpp @@ -83,6 +83,11 @@ HamiltPW::HamiltPW(elecstate::Potential* pot_in, } // no variable can choose xc, maybe it is necessary pot_register_in.push_back("xc"); + // Keep the selected nonlocal XC backend after the canonical semilocal XC. + if (pot_in->has_nonlocal_xc()) + { + pot_register_in.push_back("xc_nonlocal"); + } if (input.imp_sol) { pot_register_in.push_back("surchem"); diff --git a/source/source_pw/module_pwdft/hamilt_pw.h b/source/source_pw/module_pwdft/hamilt_pw.h index 72158eadaf..c4117810b3 100644 --- a/source/source_pw/module_pwdft/hamilt_pw.h +++ b/source/source_pw/module_pwdft/hamilt_pw.h @@ -27,6 +27,7 @@ class HamiltPW : public Hamilt using syncmem_complex_h2d_op = base_device::memory::synchronize_memory_op; public: + // Explicit component selection keeps nonlocal XC out of global state. HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index bcbc620a1a..b0976b7eea 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1,4 +1,5 @@ add_subdirectory(integrate) +add_subdirectory(rvv10) add_subdirectory(01_PW) add_subdirectory(02_NAO_Gamma) add_subdirectory(03_NAO_multik) diff --git a/tests/rvv10/CMakeLists.txt b/tests/rvv10/CMakeLists.txt new file mode 100644 index 0000000000..a9ee1775d3 --- /dev/null +++ b/tests/rvv10/CMakeLists.txt @@ -0,0 +1,10 @@ +find_package(Python3 QUIET COMPONENTS Interpreter) + +if(Python3_Interpreter_FOUND) + add_test(NAME RVV10_LIBVDWXC_CONTRACT + COMMAND ${Python3_EXECUTABLE} -m unittest discover + -s ${CMAKE_CURRENT_SOURCE_DIR} -p test_libvdwxc_contract.py) + + set_tests_properties(RVV10_LIBVDWXC_CONTRACT PROPERTIES + ENVIRONMENT "OMP_NUM_THREADS=1") +endif() diff --git a/tests/rvv10/README.md b/tests/rvv10/README.md new file mode 100644 index 0000000000..826dff73c6 --- /dev/null +++ b/tests/rvv10/README.md @@ -0,0 +1,54 @@ +# libvdwxc-backed rVV10 checks + +The numerical kernel is provided by libvdwxc MR !32, revision +`94410ebdf3dc9615a40f40c3a5b59dbe20a433d2`. This revision fixes FFT boundary +symmetry in oblique cells. Its eight library tests pass with both standalone +FFTW 3.3.11 and the MKL 2022.2.1 compatibility backend (GCC 11.2.1). +The ABACUS adapter and host tests are verified with standalone FFTW, not an +MKL-based ABACUS build. Do not mix both FFT providers in one process. +Installation and supported inputs are documented in +`docs/advanced/install.md`. + +The host adapter is intentionally limited to the current API contract: + +* CPU, double precision, full-complex plane-wave density grid; +* one process and `nspin=1` SCF; +* Libxc semilocal `GGA_X_RPW86+GGA_C_PBE` plus `xc_nonlocal=rvv10`; +* norm-conserving pseudopotentials, energy, potential and `vtxc` only. + +MPI, spin-polarized calculations, GPU/LCAO, gamma-only layouts, ultrasoft/PAW +pseudopotentials, forces and stress are rejected by the input contract and are +follow-up work. The initial integration remains Draft pending independent QE +reference data and maintainer agreement on this staged scope. + +The dependency-free test executes the Makefile build-info generator and +checks that an installed library cannot be mistaken for a compiled feature: + +```bash +cmake -B build -DBUILD_TESTING=ON +ctest --test-dir build -R RVV10_LIBVDWXC_CONTRACT --output-on-failure +``` + +With the pinned libvdwxc installation available, enable the real adapter test: + +```bash +# Configure build-libvdwxc with the standalone-FFTW command in install.md. +cmake --build build-libvdwxc --target MODULE_HAMILT_vdwTest_LIBVDWXC \ + MODULE_ESTATE_pot_nonlocal_xc_LIBVDWXC +ctest --test-dir build-libvdwxc -R LIBVDWXC --output-on-failure +``` + +The C++ tests check energy and potential against the serial API, the potential +by a finite difference of energy, and geometry using equivalent oblique and +orthogonal cells with an analytically matched reciprocal vector. Negative +tests exercise the standard ABACUS error diagnostics. +The adapter cannot be copied because it owns its library handle. A separate +host test uses real `Charge` and `UnitCell` objects to check valence/core density +construction, additive energy/potential bookkeeping and unchanged input data. + +The host reuses `PW_Basis`, `grad_rho` and `grad_dot` to supply +`n = n_valence + n_core` and `sigma = |grad(n)|^2` to libvdwxc. The library +includes the beta compensation. The returned derivatives form +`v = 2*(dedn - div(2*dedsigma*grad(n)))` in Ry; `vtxc` contracts this potential +with the valence density. The generic host component adds these quantities +after `PotXC`, while the numerical backend stays in `module_vdw`. diff --git a/tests/rvv10/test_libvdwxc_contract.py b/tests/rvv10/test_libvdwxc_contract.py new file mode 100644 index 0000000000..77612a50c4 --- /dev/null +++ b/tests/rvv10/test_libvdwxc_contract.py @@ -0,0 +1,27 @@ +"""Exercise the optional-backend build report without external dependencies.""" + +import os +from pathlib import Path +import subprocess +import tempfile +import unittest + + +ROOT = Path(__file__).resolve().parents[2] + + +class LibvdwxcContract(unittest.TestCase): + def test_makefile_reports_only_compiled_libvdwxc_support(self): + environment = os.environ.copy() + environment["ENABLE_LIBVDWXC"] = "ON" + with tempfile.TemporaryDirectory() as directory: + header = Path(directory) / "build_info.h" + subprocess.run( + ["bash", str(ROOT / "generate_build_info.sh"), str(header)], + env=environment, check=True, capture_output=True, text=True, + ) + self.assertIn('#define ABACUS_LIBVDWXC_VERSION "no"', header.read_text()) + + +if __name__ == "__main__": + unittest.main()