diff --git a/CMakeLists.txt b/CMakeLists.txt index fb0d9b38a2..671f2736d5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -49,6 +49,8 @@ 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) +cmake_dependent_option(ENABLE_LIBVDWXC "Enable using the libvdwxc package" OFF + "ENABLE_LIBXC" OFF) option(ENABLE_FLOAT_FFTW "Enable using single-precision FFTW library." OFF) option(ENABLE_LCAO "Enable LCAO algorithm" ON) @@ -285,6 +287,10 @@ if(ENABLE_DFTD4 OR ENABLE_PEXSI) enable_language(C) enable_language(Fortran) endif() +# The libvdwxc package configuration looks up MPI for C when libvdwxc was built with MPI. +if(ENABLE_LIBVDWXC) + enable_language(C) +endif() macro(set_if_higher VARIABLE VALUE) if(${VARIABLE} LESS ${VALUE}) @@ -643,6 +649,24 @@ if(ENABLE_LIBXC) endif() endif() +if(ENABLE_LIBVDWXC) + find_package(libvdwxc CONFIG REQUIRED) + # An MPI build calls vdwxc_init_mpi, which only an MPI-enabled libvdwxc provides. + if(ENABLE_MPI AND NOT VDWXC_ENABLE_MPI) + message( + FATAL_ERROR + "ENABLE_LIBVDWXC with ENABLE_MPI requires libvdwxc configured with -DVDWXC_ENABLE_MPI=ON" + ) + endif() + # vdwxc_init_mpi creates FFTW-MPI plans, which the FFTW interface of MKL or KML cannot execute. + if(ENABLE_MPI AND (MKL_FOUND OR USE_KML)) + message( + FATAL_ERROR + "ENABLE_LIBVDWXC with ENABLE_MPI requires ABACUS built with FFTW3, not with the FFTW interface of MKL or KML" + ) + endif() +endif() + if(DEFINED DeePMD_DIR) find_package(DeePMD REQUIRED) endif() diff --git a/cmake/CollectBuildInfoVars.cmake b/cmake/CollectBuildInfoVars.cmake index d74e8c0243..0894c06f90 100644 --- a/cmake/CollectBuildInfoVars.cmake +++ b/cmake/CollectBuildInfoVars.cmake @@ -192,6 +192,12 @@ else() set(ABACUS_DFTD4_VERSION "no") endif() +if(ENABLE_LIBVDWXC) + set(ABACUS_LIBVDWXC_VERSION "yes (v${libvdwxc_VERSION})") +else() + set(ABACUS_LIBVDWXC_VERSION "no") +endif() + # Accelerators if(USE_CUDA AND CUDAToolkit_VERSION) set(ABACUS_CUDA_VERSION "yes (v${CUDAToolkit_VERSION})") @@ -480,6 +486,7 @@ message(" NEP = ${ABACUS_NEP_VERSION}") message(" Features : Libxc = ${ABACUS_LIBXC_VERSION}") message(" LibRI = ${ABACUS_LIBRI_VERSION}") message(" DFT-D4 = ${ABACUS_DFTD4_VERSION}") +message(" libvdwxc = ${ABACUS_LIBVDWXC_VERSION}") message(" JSON = ${ABACUS_JSON_VERSION}") message(" PEXSI = ${ABACUS_PEXSI_VERSION}") message(" CNPY = ${ABACUS_CNPY_VERSION}") diff --git a/docs/CITATIONS.md b/docs/CITATIONS.md index 43176829a1..88a32d7fcc 100644 --- a/docs/CITATIONS.md +++ b/docs/CITATIONS.md @@ -55,3 +55,13 @@ 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 rVV10 nonlocal correlation (`vdw_method rvv10`) 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. + + R. 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. + + A. H. Larsen, M. Kuisma, J. Löfgren, Y. Pouillon, P. Erhart, and P. Hyldgaard. "libvdwxc: a library for exchange-correlation functionals in the vdW-DF family." Modelling and Simulation in Materials Science and Engineering 25, 065004 (2017). DOI: 10.1088/1361-651X/aa7320. + + For SCAN+rVV10 and r2SCAN+rVV10, also cite the parametrization: H. Peng, Z.-H. Yang, J. P. Perdew, and J. Sun. "Versatile van der Waals density functional based on a meta-generalized gradient approximation." Physical Review X 6, 041005 (2016). DOI: 10.1103/PhysRevX.6.041005; J. Ning, M. Kothakonda, J. W. Furness, A. D. Kaplan, S. Ehlert, J. G. Brandenburg, J. P. Perdew, and J. Sun. "Workhorse minimally empirical dispersion-corrected density functional with tests for weakly bound systems: r2SCAN+rVV10." Physical Review B 106, 075422 (2022). DOI: 10.1103/PhysRevB.106.075422. diff --git a/docs/advanced/install.md b/docs/advanced/install.md index 86d63cfd1e..ceb4fcf9ad 100644 --- a/docs/advanced/install.md +++ b/docs/advanced/install.md @@ -99,6 +99,25 @@ 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 support + +`vdw_method rvv10` evaluates the rVV10 nonlocal correlation with the external [libvdwxc](https://gitlab.com/libvdwxc/libvdwxc) library. Until a tagged libvdwxc release contains the rVV10 interface, build the branch of [libvdwxc merge request !32](https://gitlab.com/libvdwxc/libvdwxc/-/merge_requests/32) (tested at commit `da2d55f`) with libvdwxc's CMake build, which installs the package that ABACUS looks up, and against the FFTW that ABACUS uses. For an MPI build of ABACUS, configure libvdwxc with MPI, which needs FFTW with its MPI library: + +```bash +cmake -S libvdwxc -B build-libvdwxc -DBUILD_SHARED_LIBS=ON -DVDWXC_ENABLE_MPI=ON -DCMAKE_INSTALL_PREFIX=/path/to/libvdwxc +cmake --build build-libvdwxc && cmake --install build-libvdwxc +``` + +Then pass `-DENABLE_LIBXC=ON -DENABLE_LIBVDWXC=ON` to CMake for ABACUS and add `/path/to/libvdwxc` to `CMAKE_PREFIX_PATH`. ABACUS hands each k-point pool's density to libvdwxc's distributed interface, so the nonlocal convolution and its memory are shared by the processes of the pool. libvdwxc's MPI build plans its FFTs with FFTW-MPI, so an MPI build of ABACUS with libvdwxc must use FFTW3 rather than the FFTW interface of MKL or KML; CMake stops otherwise. The legacy Makefile does not build this feature. + +In the input file: + +```text +vdw_method rvv10 +dft_functional scan +vdw_rvv10_b 15.7 # b of SCAN+rVV10; see vdw_rvv10_b for other functionals +``` + ## 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/docs/advanced/scf/advanced.md b/docs/advanced/scf/advanced.md index a01280344b..afbd4dade7 100644 --- a/docs/advanced/scf/advanced.md +++ b/docs/advanced/scf/advanced.md @@ -96,6 +96,7 @@ To use VdW-correction, users need to supply value to the `vdw_method` keyword in - d2: DFT-D2 method - d3_0: DFT-D3(0) method - d3_bj: DFT-D3(BJ) method + - rvv10: rVV10 nonlocal correlation (requires libvdwxc). Unlike the DFT-D methods, it depends on the electron density: it is added self-consistently to the exchange-correlation energy and potential of `dft_functional`, and its energy is reported in the exchange-correlation term. Set `vdw_rvv10_b` for the functional; its description lists the published values. Furthermore, ABACUS also provides a [list of keywords](../input_files/input-main.md#vdw-correction) to control relevant parmeters used in calculating the VdW correction, such as the scale factor (s6) term. Recommended values of such parameters can be found on the [webpage](https://www.chemie.uni-bonn.de/pctc/mulliken-center/software/dft-d3). The default values of the parameters in ABACUS are set to be the recommended values for PBE. diff --git a/generate_build_info.sh b/generate_build_info.sh index f50f98ac72..afe6325ca5 100755 --- a/generate_build_info.sh +++ b/generate_build_info.sh @@ -279,6 +279,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@#no#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 d849396a7d..287206ef03 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -67,6 +67,7 @@ set(_abacus_feature_definitions $<$:TEST_EXX_RADIAL=1> $<$:__LIBXC> $<$:__DFTD4> + $<$:__LIBVDWXC> $<$:GATHER_INFO>) if(ABACUS_COMMIT_INFO_INCLUDE_DIR) @@ -407,6 +408,7 @@ foreach( ELPA::ELPA Libxc::xc dftd4::dftd4 + libvdwxc::vdwxc DeePMD::deepmd_c DeePMD::deepmd_cc NEP::nep diff --git a/source/source_base/parallel_device.cpp b/source/source_base/parallel_device.cpp index 8ba8c3ae8e..a66837d2b4 100644 --- a/source/source_base/parallel_device.cpp +++ b/source/source_base/parallel_device.cpp @@ -350,6 +350,16 @@ void gatherv_data(const std::complex* sendbuf, int sendcount, std::comple { MPI_Allgatherv(sendbuf, sendcount, MPI_COMPLEX, recvbuf, recvcounts, displs, MPI_COMPLEX, comm); } +void alltoallv_data(const double* sendbuf, + const int* sendcounts, + const int* sdispls, + double* recvbuf, + const int* recvcounts, + const int* rdispls, + const MPI_Comm& comm) +{ + MPI_Alltoallv(sendbuf, sendcounts, sdispls, MPI_DOUBLE, recvbuf, recvcounts, rdispls, MPI_DOUBLE, comm); +} #ifndef __CUDA_MPI template diff --git a/source/source_base/parallel_device.h b/source/source_base/parallel_device.h index 7826e4b465..f4d35e4e99 100644 --- a/source/source_base/parallel_device.h +++ b/source/source_base/parallel_device.h @@ -31,6 +31,13 @@ void gatherv_data(const double* sendbuf, int sendcount, double* recvbuf, const i void gatherv_data(const std::complex* sendbuf, int sendcount, std::complex* recvbuf, const int* recvcounts, const int* displs, MPI_Comm& comm); void gatherv_data(const float* sendbuf, int sendcount, float* recvbuf, const int* recvcounts, const int* displs, MPI_Comm& comm); void gatherv_data(const std::complex* sendbuf, int sendcount, std::complex* recvbuf, const int* recvcounts, const int* displs, MPI_Comm& comm); +void alltoallv_data(const double* sendbuf, + const int* sendcounts, + const int* sdispls, + double* recvbuf, + const int* recvcounts, + const int* rdispls, + const MPI_Comm& comm); #if defined(__NCCL_PARALLEL_DEVICE) void nccl_bcast_data(double* object, const int& n, MPI_Comm& comm, int root = 0); diff --git a/source/source_base/test_parallel/parallel_device_test.cpp b/source/source_base/test_parallel/parallel_device_test.cpp index ba4a4d90b8..118d73f183 100644 --- a/source/source_base/test_parallel/parallel_device_test.cpp +++ b/source/source_base/test_parallel/parallel_device_test.cpp @@ -131,5 +131,48 @@ TEST(ParallelDevice, CoversMpiTypeOverloads) exercise_mpi_wrappers>(domain); } +TEST(ParallelDevice, AlltoallvExchangesUnevenBlocks) +{ + const ModuleBase::CommunicationDomain domain = ModuleBase::world_comm_domain(); + MPI_Comm communicator = domain.communicator(); + const int rank = domain.rank(); + MPICommGroup world_group(communicator); + const int size = world_group.gsize; + // rank r sends (r + p) % 2 + 1 copies of 100 r + p to rank p + std::vector send_counts(size); + std::vector send_offsets(size); + std::vector recv_counts(size); + std::vector recv_offsets(size); + std::vector sent; + int received_total = 0; + for (int peer = 0; peer < size; ++peer) + { + const double value = 100.0 * rank + peer; + send_counts[peer] = (rank + peer) % 2 + 1; + send_offsets[peer] = static_cast(sent.size()); + sent.insert(sent.end(), send_counts[peer], value); + recv_counts[peer] = send_counts[peer]; + recv_offsets[peer] = received_total; + received_total += recv_counts[peer]; + } + std::vector received(received_total); + Parallel_Common::alltoallv_data(sent.data(), + send_counts.data(), + send_offsets.data(), + received.data(), + recv_counts.data(), + recv_offsets.data(), + communicator); + for (int peer = 0; peer < size; ++peer) + { + const double expected = 100.0 * peer + rank; + for (int index = 0; index < recv_counts[peer]; ++index) + { + const int position = recv_offsets[peer] + index; + EXPECT_EQ(received[position], expected); + } + } +} + } // namespace #endif diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp index 149f50dc8c..97badb8360 100644 --- a/source/source_esolver/esolver_fp.cpp +++ b/source/source_esolver/esolver_fp.cpp @@ -60,6 +60,9 @@ void ESolver_FP::before_all_runners(BaseCell& basecell, const Input_para& inp) xc_parameters.exx_erfc_alpha = inp.exx_erfc_alpha; xc_parameters.xc_exch_ext = inp.xc_exch_ext; xc_parameters.xc_corr_ext = inp.xc_corr_ext; + xc_parameters.rvv10 = inp.vdw_method == "rvv10"; + xc_parameters.rvv10_b = inp.vdw_rvv10_b; + xc_parameters.rvv10_c = inp.vdw_rvv10_c; XC_Functional::set_runtime_parameters(xc_parameters); ModuleBase::TITLE("ESolver_FP", "before_all_runners"); 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..4e582307ab --- /dev/null +++ b/source/source_hamilt/module_vdw/nonlocal_xc.cpp @@ -0,0 +1,427 @@ +#include "nonlocal_xc.h" + +#include "source_base/constants.h" +#include "source_base/matrix.h" +#include "source_base/matrix3.h" +#include "source_base/parallel_device.h" +#include "source_base/timer.h" +#include "source_base/tool_quit.h" +#include "source_base/tool_title.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_estate/module_charge/charge.h" +#include "source_hamilt/module_xc/xc_functional.h" + +#ifdef __MPI +// vdwxc_mpi.h includes inside an extern "C" block, so include it first here. +#include +#endif +// vdwxc.h includes inside an extern "C" block. Include the C++ header first so that +// libc++ does not parse its templates with C linkage, then drop the I macro. +#include +#include +#include +#ifdef __MPI +#include +#endif +#undef I + +#include +#include + +namespace vdw +{ +namespace +{ +/// Cartesian components of the stress tensor returned by vdwxc_stress, stored as [3 * a + b] +const int STRESS_COMPONENTS = 9; + +/** + * libvdwxc object for one FFT grid and pool, kept between calls and rebuilt when either changes. + * + * libvdwxc gets the lattice vectors in the cyclic order (a3, a1, a2), so that the first axis of its grid + * is the z axis that ABACUS distributes. It splits that axis into blocks of ceil(nz/P) planes, the FFTW-MPI + * default; vdwxc_init_mpi does not report the local slab, so plan_exchange computes it the same way. The + * slab of a process is stored as [z][x][y], the transpose of the ABACUS layout [x][y][z], and whole z + * planes are moved between the two distributions by all-to-all exchanges. + */ +class Rvv10Workspace +{ + public: + /// libvdwxc energy (Ha) and derivatives of the ABACUS-distributed fields rho and sigma; with + /// stress != nullptr also the explicit strain derivative (Ha, Cartesian, not divided by the volume). + double evaluate(const ModulePW::PW_Basis& pw, + const double b, + const double c, + const std::vector& rho, + const std::vector& sigma, + std::vector& dedn, + std::vector& dedsigma, + double* stress); + + private: + void prepare(const ModulePW::PW_Basis& pw, const double b, const double c); + void plan_exchange(const ModulePW::PW_Basis& pw); + void to_libvdwxc(const ModulePW::PW_Basis& pw, const std::vector& field, std::vector& slab); + void from_libvdwxc(const ModulePW::PW_Basis& pw, const std::vector& slab, std::vector& field); + + vdwxc_data data_ = nullptr; + int nx_ = 0; + int ny_ = 0; + int nz_ = 0; + int nproc_ = 0; +#ifdef __MPI + MPI_Comm comm_ = MPI_COMM_NULL; +#endif + // planes moved to (send) and from (recv) each process, in doubles + std::vector send_count_; + std::vector send_offset_; + std::vector recv_count_; + std::vector recv_offset_; + // the ABACUS planes of this process as [z][xy], and the libvdwxc slab + std::vector planes_; + std::vector rho_; + std::vector sigma_; + std::vector dedn_; + std::vector dedsigma_; +}; + +void Rvv10Workspace::prepare(const ModulePW::PW_Basis& pw, const double b, const double c) +{ + bool same_grid = data_ != nullptr && pw.nx == nx_ && pw.ny == ny_ && pw.nz == nz_ && pw.poolnproc == nproc_; +#ifdef __MPI + same_grid = same_grid && (pw.poolnproc == 1 || pw.pool_world == comm_); +#endif + if (!same_grid) + { + if (data_ != nullptr) + { + vdwxc_finalize(&data_); + } + data_ = vdwxc_new(VDWXC_RVV10); + nx_ = pw.nx; + ny_ = pw.ny; + nz_ = pw.nz; + nproc_ = pw.poolnproc; + } + // The cell is set on every call: a cell relaxation keeps the FFT grid but changes the cell. libvdwxc + // integrates with the signed cell volume, and inverting a left-handed cell leaves rVV10 unchanged. + ModuleBase::Matrix3 a = pw.latvec * pw.lat0; + if (a.Det() < 0.0) + { + a = a * -1.0; + } + vdwxc_set_unit_cell(data_, pw.nz, pw.nx, pw.ny, a.e31, a.e32, a.e33, a.e11, a.e12, a.e13, a.e21, a.e22, a.e23); + if (!same_grid) + { +#ifdef __MPI + comm_ = pw.pool_world; + if (pw.poolnproc > 1) + { + vdwxc_init_mpi(data_, comm_); + } + else + { + vdwxc_init_serial(data_); + } +#else + vdwxc_init_serial(data_); +#endif + } + // The libvdwxc slabs depend only on nz and the number of processes, but the ABACUS planes depend on + // the basis, so the exchange is planned on every call. + this->plan_exchange(pw); + if (vdwxc_set_rvv10_parameters(data_, b, c) != 0) + { + ModuleBase::WARNING_QUIT("vdw::add_rvv10", "libvdwxc rejected vdw_rvv10_b or vdw_rvv10_c"); + } +} + +void Rvv10Workspace::plan_exchange(const ModulePW::PW_Basis& pw) +{ + const int nproc = pw.poolnproc; + const int block = (pw.nz + nproc - 1) / nproc; + const int first_slab_plane = pw.poolrank * block; + const int slab_start = std::min(first_slab_plane, pw.nz); + const int planes_after_start = pw.nz - slab_start; + const int nslab = std::min(block, planes_after_start); + const int slab_end = slab_start + nslab; + const int plane_end = pw.startz_current + pw.nplane; + send_count_.assign(nproc, 0); + send_offset_.assign(nproc, 0); + recv_count_.assign(nproc, 0); + recv_offset_.assign(nproc, 0); + for (int ip = 0; nproc > 1 && ip < nproc; ++ip) + { + // ABACUS planes of this process that lie in the slab of ip + const int first_plane_of_ip = ip * block; + const int begin = std::min(first_plane_of_ip, pw.nz); + const int block_end = begin + block; + const int end = std::min(block_end, pw.nz); + int low = std::max(pw.startz_current, begin); + int high = std::min(plane_end, end); + if (high > low) + { + send_count_[ip] = (high - low) * pw.nxy; + send_offset_[ip] = (low - pw.startz_current) * pw.nxy; + } + // ABACUS planes of ip that lie in the slab of this process + const int plane_end_of_ip = pw.startz[ip] + pw.numz[ip]; + low = std::max(pw.startz[ip], slab_start); + high = std::min(plane_end_of_ip, slab_end); + if (high > low) + { + recv_count_[ip] = (high - low) * pw.nxy; + recv_offset_[ip] = (low - slab_start) * pw.nxy; + } + } + // libvdwxc expects valid pointers also for an empty slab + const int slab_points = nslab * pw.nxy; + const int slab_size = std::max(1, slab_points); + // the transposed ABACUS planes are only needed for the exchange between processes + const int planes_size = (nproc > 1) ? std::max(1, pw.nrxx) : 0; + planes_.resize(planes_size); + rho_.resize(slab_size); + sigma_.resize(slab_size); + dedn_.resize(slab_size); + dedsigma_.resize(slab_size); +} + +void Rvv10Workspace::to_libvdwxc(const ModulePW::PW_Basis& pw, + const std::vector& field, + std::vector& slab) +{ + std::vector& planes = (pw.poolnproc > 1) ? planes_ : slab; + for (int ixy = 0; ixy < pw.nxy; ++ixy) + { + for (int iz = 0; iz < pw.nplane; ++iz) + { + planes[iz * pw.nxy + ixy] = field[ixy * pw.nplane + iz]; + } + } +#ifdef __MPI + if (pw.poolnproc > 1) + { + Parallel_Common::alltoallv_data(planes_.data(), + send_count_.data(), + send_offset_.data(), + slab.data(), + recv_count_.data(), + recv_offset_.data(), + pw.pool_world); + } +#endif +} + +void Rvv10Workspace::from_libvdwxc(const ModulePW::PW_Basis& pw, + const std::vector& slab, + std::vector& field) +{ +#ifdef __MPI + if (pw.poolnproc > 1) + { + Parallel_Common::alltoallv_data(slab.data(), + recv_count_.data(), + recv_offset_.data(), + planes_.data(), + send_count_.data(), + send_offset_.data(), + pw.pool_world); + } +#endif + const std::vector& planes = (pw.poolnproc > 1) ? planes_ : slab; + for (int ixy = 0; ixy < pw.nxy; ++ixy) + { + for (int iz = 0; iz < pw.nplane; ++iz) + { + field[ixy * pw.nplane + iz] = planes[iz * pw.nxy + ixy]; + } + } +} + +double Rvv10Workspace::evaluate(const ModulePW::PW_Basis& pw, + const double b, + const double c, + const std::vector& rho, + const std::vector& sigma, + std::vector& dedn, + std::vector& dedsigma, + double* stress) +{ + this->prepare(pw, b, c); + this->to_libvdwxc(pw, rho, rho_); + this->to_libvdwxc(pw, sigma, sigma_); + double energy = 0.0; + if (stress == nullptr) + { + energy = vdwxc_calculate(data_, rho_.data(), sigma_.data(), dedn_.data(), dedsigma_.data()); + } + else + { + energy = vdwxc_stress(data_, stress, rho_.data(), sigma_.data(), dedn_.data(), dedsigma_.data()); + } + this->from_libvdwxc(pw, dedn_, dedn); + this->from_libvdwxc(pw, dedsigma_, dedsigma); +#ifdef __MPI + // the distributed backend returns the contribution of this process + if (pw.poolnproc > 1) + { + Parallel_Common::reduce_data(&energy, 1, pw.pool_world); + if (stress != nullptr) + { + Parallel_Common::reduce_data(stress, STRESS_COMPONENTS, pw.pool_world); + } + } +#endif + return energy; +} + +/// The workspace is never destroyed: a static object would be destroyed after MPI_Finalize and +/// fftw_cleanup_threads, which libvdwxc's FFT plans must not outlive. +Rvv10Workspace& workspace() +{ + static Rvv10Workspace* const instance = new Rvv10Workspace; + return *instance; +} + +/// Total density of chr (valence plus core), its plane-wave gradient and sigma = |grad n|^2. +void total_density(const Charge& chr, + std::vector& rho, + std::vector>& grad, + std::vector& sigma) +{ + const ModulePW::PW_Basis& pw = *chr.rhopw; + const int nrxx = chr.nrxx; + const double* core_begin = chr.rho_core; + const double* core_end = core_begin + nrxx; + // for nspin = 4, rho[0] is the total valence density + const int nspin0 = (chr.nspin == 2) ? 2 : 1; + rho.assign(core_begin, core_end); + for (int is = 0; is < nspin0; ++is) + { + for (int ir = 0; ir < nrxx; ++ir) + { + rho[ir] += chr.rho[is][ir]; + } + } + std::vector> rhog(pw.npw); + pw.real2recip(rho.data(), rhog.data()); + grad.resize(nrxx); + XC_Functional::grad_rho(rhog.data(), grad.data(), &pw, pw.tpiba); + sigma.resize(nrxx); + for (int ir = 0; ir < nrxx; ++ir) + { + sigma[ir] = grad[ir].norm2(); + } +} +} // namespace + +void add_rvv10(const XCFunctionalParameters& parameters, + const Charge& chr, + double& etxc, + double& vtxc, + ModuleBase::matrix& v) +{ + if (!parameters.rvv10) + { + return; + } + ModuleBase::TITLE("vdw", "add_rvv10"); + ModuleBase::timer::start("vdw", "add_rvv10"); + const ModulePW::PW_Basis& pw = *chr.rhopw; + const int nrxx = chr.nrxx; + std::vector rho; + std::vector> grad; + std::vector sigma; + total_density(chr, rho, grad, sigma); + std::vector dedn(nrxx); + std::vector dedsigma(nrxx); + const double energy + = workspace().evaluate(pw, parameters.rvv10_b, parameters.rvv10_c, rho, sigma, dedn, dedsigma, nullptr); + + // v = dE/dn - div(2 dE/dsigma grad n), converted from Hartree to Ry by e2 + for (int ir = 0; ir < nrxx; ++ir) + { + grad[ir] *= 2.0 * dedsigma[ir]; + } + std::vector divergence(nrxx); + XC_Functional::grad_dot(grad.data(), divergence.data(), &pw, pw.tpiba); + // for nspin = 4 the potential of the total density goes to v(0) only + const int nspin0 = (chr.nspin == 2) ? 2 : 1; + double vtxc_nl = 0.0; + for (int ir = 0; ir < nrxx; ++ir) + { + const double v_nl = ModuleBase::e2 * (dedn[ir] - divergence[ir]); + for (int is = 0; is < nspin0; ++is) + { + v(is, ir) += v_nl; + vtxc_nl += v_nl * chr.rho[is][ir]; + } + } +#ifdef __MPI + Parallel_Common::reduce_data(&vtxc_nl, 1, pw.pool_world); +#endif + etxc += ModuleBase::e2 * energy; + vtxc += vtxc_nl * pw.omega / pw.nxyz; + ModuleBase::timer::end("vdw", "add_rvv10"); +} + +void add_rvv10_stress(const XCFunctionalParameters& parameters, const Charge& chr, ModuleBase::matrix& sigma) +{ + if (!parameters.rvv10) + { + return; + } + ModuleBase::TITLE("vdw", "add_rvv10_stress"); + ModuleBase::timer::start("vdw", "add_rvv10_stress"); + const ModulePW::PW_Basis& pw = *chr.rhopw; + const int nrxx = chr.nrxx; + std::vector rho; + std::vector> grad; + std::vector density_sigma; + total_density(chr, rho, grad, density_sigma); + std::vector dedn(nrxx); + std::vector dedsigma(nrxx); + double explicit_stress[STRESS_COMPONENTS] = {0.0}; + const double energy = workspace().evaluate(pw, + parameters.rvv10_b, + parameters.rvv10_c, + rho, + density_sigma, + dedn, + dedsigma, + explicit_stress); + + // integral of dE/dsigma d_a n d_b n, in units of the volume element + double gradient_term[STRESS_COMPONENTS] = {0.0}; + for (int ir = 0; ir < nrxx; ++ir) + { + const double g[3] = {grad[ir].x, grad[ir].y, grad[ir].z}; + for (int a = 0; a < 3; ++a) + { + for (int b = 0; b < 3; ++b) + { + gradient_term[3 * a + b] += dedsigma[ir] * g[a] * g[b]; + } + } + } +#ifdef __MPI + Parallel_Common::reduce_data(gradient_term, STRESS_COMPONENTS, pw.pool_world); +#endif + + // The full strain derivative is T - delta_ab Int(n dE/dn + 2 sigma dE/dsigma) - 2 Int dE/dsigma d_a n d_b n, + // with T the explicit part from libvdwxc. The XC diagonal -(etxc - vtxc)/omega in PW (in LCAO, -etxc/omega + // and the Pulay term of the effective potential) already holds -e2 (E - Int(n dE/dn + 2 sigma dE/dsigma)) / omega, + // which leaves the terms below (ABACUS sign, Ry). + for (int a = 0; a < 3; ++a) + { + for (int b = 0; b < 3; ++b) + { + const double diagonal = (a == b) ? energy : 0.0; + sigma(a, b) += -ModuleBase::e2 / pw.omega * (explicit_stress[3 * a + b] - diagonal) + + 2.0 * ModuleBase::e2 * gradient_term[3 * a + b] / pw.nxyz; + } + } + ModuleBase::timer::end("vdw", "add_rvv10_stress"); +} + +} // 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..22eb4a56b3 --- /dev/null +++ b/source/source_hamilt/module_vdw/nonlocal_xc.h @@ -0,0 +1,31 @@ +#ifndef NONLOCAL_XC_H +#define NONLOCAL_XC_H + +class Charge; +struct XCFunctionalParameters; +namespace ModuleBase +{ +class matrix; +} + +namespace vdw +{ + +// rVV10 nonlocal correlation of vdw_method rvv10, evaluated by libvdwxc. Both functions do nothing +// unless parameters.rvv10 is set. They use the total density of chr (valence plus core) and are +// collective over the pool communicator of chr.rhopw. The libvdwxc object, its FFT plans and its +// buffers are kept between calls. + +/// Add the rVV10 energy, potential and vtxc (Ry) to the pool-reduced results of an XC evaluation. +void add_rvv10(const XCFunctionalParameters& parameters, + const Charge& chr, + double& etxc, + double& vtxc, + ModuleBase::matrix& v); + +/// Add the rVV10 stress (Ry/Bohr^3) that the XC diagonal and potential terms do not contain. +void add_rvv10_stress(const XCFunctionalParameters& parameters, const Charge& chr, ModuleBase::matrix& sigma); + +} // namespace vdw + +#endif diff --git a/source/source_hamilt/module_vdw/test/CMakeLists.txt b/source/source_hamilt/module_vdw/test/CMakeLists.txt index fddd34b681..2f056a6999 100644 --- a/source/source_hamilt/module_vdw/test/CMakeLists.txt +++ b/source/source_hamilt/module_vdw/test/CMakeLists.txt @@ -22,3 +22,35 @@ if(ENABLE_DFTD4 AND ENABLE_MPI) WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} ) endif() + +if(ENABLE_LIBVDWXC AND ENABLE_MPI) + addtest( + TARGET MODULE_HAMILT_nonlocal_xcTest + LIBS xc_ + planewave + parameter + MPI::MPI_CXX + Libxc::xc + psi + device + container + base + SOURCES test_nonlocal_xc.cpp) + set_tests_properties(MODULE_HAMILT_nonlocal_xcTest + PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=1") + add_test( + NAME MODULE_HAMILT_nonlocal_xcTest_2np + COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 2 ${MPIEXEC_PREFLAGS} + $ ${MPIEXEC_POSTFLAGS} + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + set_tests_properties(MODULE_HAMILT_nonlocal_xcTest_2np + PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=1" PROCESSORS 2) + # four processes give an uneven z distribution and an empty libvdwxc slab + add_test( + NAME MODULE_HAMILT_nonlocal_xcTest_4np + COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} + $ ${MPIEXEC_POSTFLAGS} + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + set_tests_properties(MODULE_HAMILT_nonlocal_xcTest_4np + PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=1" PROCESSORS 4) +endif() diff --git a/source/source_hamilt/module_vdw/test/test_nonlocal_xc.cpp b/source/source_hamilt/module_vdw/test/test_nonlocal_xc.cpp new file mode 100644 index 0000000000..59689d1693 --- /dev/null +++ b/source/source_hamilt/module_vdw/test/test_nonlocal_xc.cpp @@ -0,0 +1,633 @@ +#include "../nonlocal_xc.h" +#include "source_base/constants.h" +#include "source_base/matrix.h" +#include "source_base/parallel_comm.h" +#include "source_base/parallel_device.h" +#include "source_base/parallel_global.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_cell/unitcell.h" +#include "source_estate/module_charge/charge.h" +#include "source_hamilt/module_xc/libxc_abacus.h" +#include "source_hamilt/module_xc/xc_functional.h" + +#include +#include +#include +#include + +// The reference calls libvdwxc directly. Include its C header after the C++ headers. +#include +#include +#include +#undef I + +// The target links PW/XC objects only; the tests own the charge storage. +Charge::Charge() +{ +} +Charge::~Charge() +{ +} +UnitCell::UnitCell() +{ +} +UnitCell::~UnitCell() +{ +} +Magnetism::Magnetism() +{ +} +Magnetism::~Magnetism() +{ +} +SepPot::SepPot() +{ +} +SepPot::~SepPot() +{ +} +Sep_Cell::Sep_Cell() noexcept +{ +} +Sep_Cell::~Sep_Cell() noexcept +{ +} + +namespace +{ +int test_rank = 0; +int test_size = 1; +typedef std::tuple VxcResult; +typedef std::tuple MetaResult; +// rVV10 parameters of the rPW86-PBE pairing +const double RVV10_B = 6.3; +const double RVV10_C = 0.0093; +// Results that share the plane-wave gradient agree to rounding. +const double ROUNDING_TOLERANCE = 2.0e-12; +// Central differences on smooth densities: the truncation error is of order step^2. +const double FINITE_DIFFERENCE_STEP = 1.0e-5; +const double FINITE_DIFFERENCE_TOLERANCE = 2.0e-6; +const double STRAIN_STEP = 1.0e-4; +const double STRAIN_TOLERANCE = 1.0e-5; +const double BASIS_CUTOFF = 4.0; +// PW_Basis lowers a larger cutoff to the largest G sphere inside the FFT box. The test densities have +// only first-shell G components, well inside that sphere, so a small strain changes none of them. +const double LARGEST_SPHERE_CUTOFF = 1.0e6; + +typedef double (*GridFunction)(const ModulePW::PW_Basis&, const int, const int, const int); + +double valence_density(const ModulePW::PW_Basis& pw, const int ix, const int iy, const int iz) +{ + const double ax = ModuleBase::TWO_PI * ix / pw.nx; + const double ay = ModuleBase::TWO_PI * iy / pw.ny; + const double az = ModuleBase::TWO_PI * iz / pw.nz; + const double ayz = ay + az; + return 0.004 + 0.006 * (1.0 + 0.25 * std::sin(ax) + 0.15 * std::cos(ayz)); +} + +double core_density(const ModulePW::PW_Basis& pw, const int ix, const int iy, const int iz) +{ + const double az = ModuleBase::TWO_PI * iz / pw.nz; + return 0.002 * (1.0 + std::cos(az)); +} + +double kinetic_density(const ModulePW::PW_Basis& pw, const int ix, const int iy, const int iz) +{ + const double ay = ModuleBase::TWO_PI * iy / pw.ny; + return 0.003 * (1.0 + 0.2 * std::cos(ay)); +} + +double perturbation(const ModulePW::PW_Basis& pw, const int ix, const int iy, const int iz) +{ + const double ax = ModuleBase::TWO_PI * ix / pw.nx; + const double ay = ModuleBase::TWO_PI * iy / pw.ny; + const double axy = ax + ay; + const double ay2 = 2.0 * ay; + return std::sin(axy) + 0.3 * std::cos(ay2); +} + +/// Values of f on the z planes of this process, in the ABACUS layout [ixy * nplane + iz]. +std::vector local_grid(const ModulePW::PW_Basis& pw, GridFunction f) +{ + std::vector values(pw.nrxx); + for (int ix = 0; ix < pw.nx; ++ix) + { + for (int iy = 0; iy < pw.ny; ++iy) + { + for (int iz = 0; iz < pw.nplane; ++iz) + { + values[(ix * pw.ny + iy) * pw.nplane + iz] = f(pw, ix, iy, pw.startz_current + iz); + } + } + } + return values; +} + +/// Oblique lattice; z_sign = -1 gives its left-handed mirror image. +ModuleBase::Matrix3 oblique_lattice(const double z_sign) +{ + const double c3 = 2.0 * z_sign; + return ModuleBase::Matrix3(1.0, 0.0, 0.0, -0.5, std::sqrt(0.75), 0.0, 0.0, 0.0, c3); +} + +/// The lattice after the strain eps_ab = h: every lattice vector gets a_i[a] += h * a_i[b]. +ModuleBase::Matrix3 strained_lattice(const ModuleBase::Matrix3& lattice, const int a, const int b, const double h) +{ + double m[3][3] = {{lattice.e11, lattice.e12, lattice.e13}, + {lattice.e21, lattice.e22, lattice.e23}, + {lattice.e31, lattice.e32, lattice.e33}}; + for (int i = 0; i < 3; ++i) + { + m[i][a] += h * m[i][b]; + } + return ModuleBase::Matrix3(m[0][0], m[0][1], m[0][2], m[1][0], m[1][1], m[1][2], m[2][0], m[2][1], m[2][2]); +} + +/// Basis on the nx x ny x nz FFT grid of the given cell. +void setup_grid(ModulePW::PW_Basis& basis, + const ModuleBase::Matrix3& lattice, + const double cutoff, + const int nx, + const int ny, + const int nz) +{ + basis.initgrids(10.0, lattice, nx, ny, nz); + basis.initparameters(false, cutoff, 1, true); + basis.setuptransform(); + basis.collect_local_pw(); +} + +/// 5x7x9 grid: on two processes the z planes split 5 + 4, on four processes 3 + 2 + 2 + 2 against +/// libvdwxc's 3 + 3 + 3 + 0, which exercises the redistribution and an empty slab. Odd nx and ny leave +/// no Nyquist plane, so libvdwxc's result does not depend on the order in which it gets the axes. +void setup(ModulePW::PW_Basis& basis, const ModuleBase::Matrix3& lattice, const double cutoff) +{ + setup_grid(basis, lattice, cutoff, 5, 7, 9); +} + +/// XC runtime parameters with the rVV10 switch set as given and the parameters of the rPW86-PBE pairing. +XCFunctionalParameters rvv10_parameters(const bool rvv10) +{ + XCFunctionalParameters parameters = XC_Functional::get_runtime_parameters(); + parameters.rvv10 = rvv10; + parameters.rvv10_b = RVV10_B; + parameters.rvv10_c = RVV10_C; + return parameters; +} + +/// Sets the XC runtime parameters for the lifetime of the object, as ESolver_FP does. +class RuntimeParameters +{ + public: + explicit RuntimeParameters(const bool rvv10) + { + saved_ = XC_Functional::get_runtime_parameters(); + const XCFunctionalParameters parameters = rvv10_parameters(rvv10); + XC_Functional::set_runtime_parameters(parameters); + } + ~RuntimeParameters() + { + XC_Functional::set_runtime_parameters(saved_); + } + + private: + XCFunctionalParameters saved_; +}; + +/// A Charge on pw that owns its density components, core density and the G-space buffers of gradcorr. +class ChargeOnGrid +{ + public: + ChargeOnGrid(ModulePW::PW_Basis& pw, + const std::vector>& channels, + const std::vector& core) + { + channels_ = channels; + core_ = core; + const std::size_t ncomponents = channels.size(); + const std::vector> reciprocal(pw.npw); + rhog_.assign(ncomponents, reciprocal); + core_g_ = reciprocal; + for (std::size_t is = 0; is < ncomponents; ++is) + { + rows_.push_back(channels_[is].data()); + rhog_rows_.push_back(rhog_[is].data()); + } + chr.rhopw = &pw; + chr.nrxx = pw.nrxx; + chr.nspin = static_cast(ncomponents); + chr.rho = rows_.data(); + chr.rhog = rhog_rows_.data(); + chr.rho_core = core_.data(); + chr.rhog_core = core_g_.data(); + } + ChargeOnGrid(const ChargeOnGrid&) = delete; + ChargeOnGrid& operator=(const ChargeOnGrid&) = delete; + + Charge chr; + + private: + std::vector> channels_; + std::vector rows_; + std::vector core_; + std::vector>> rhog_; + std::vector*> rhog_rows_; + std::vector> core_g_; +}; + +/// rVV10 contribution for the given valence channels and core density on pw. +VxcResult rvv10_result(ModulePW::PW_Basis& pw, + const std::vector>& channels, + const std::vector& core) +{ + const ChargeOnGrid charge(pw, channels, core); + const XCFunctionalParameters parameters = rvv10_parameters(true); + double etxc = 0.0; + double vtxc = 0.0; + ModuleBase::matrix v(charge.chr.nspin, pw.nrxx); + vdw::add_rvv10(parameters, charge.chr, etxc, vtxc, v); + return VxcResult(etxc, vtxc, v); +} + +/// rVV10 stress beyond the XC diagonal term for the given density on pw. +ModuleBase::matrix rvv10_stress(ModulePW::PW_Basis& pw, + const std::vector>& channels, + const std::vector& core) +{ + const ChargeOnGrid charge(pw, channels, core); + const XCFunctionalParameters parameters = rvv10_parameters(true); + ModuleBase::matrix sigma(3, 3); + vdw::add_rvv10_stress(parameters, charge.chr, sigma); + return sigma; +} + +void expect_same_result(const VxcResult& actual, const VxcResult& expected, const double tolerance) +{ + EXPECT_NEAR(std::get<0>(actual), std::get<0>(expected), tolerance); + EXPECT_NEAR(std::get<1>(actual), std::get<1>(expected), tolerance); + const ModuleBase::matrix& v = std::get<2>(actual); + const ModuleBase::matrix& w = std::get<2>(expected); + ASSERT_EQ(v.nr, w.nr); + ASSERT_EQ(v.nc, w.nc); + for (int is = 0; is < v.nr; ++is) + { + for (int ir = 0; ir < v.nc; ++ir) + { + EXPECT_NEAR(v(is, ir), w(is, ir), tolerance); + } + } +} + +/// A functional evaluated with rVV10 equals the evaluation without it plus add_rvv10. +void expect_rvv10_added(const VxcResult& total, VxcResult semilocal, const Charge& chr) +{ + const XCFunctionalParameters parameters = rvv10_parameters(true); + double& etxc = std::get<0>(semilocal); + double& vtxc = std::get<1>(semilocal); + ModuleBase::matrix& v = std::get<2>(semilocal); + vdw::add_rvv10(parameters, chr, etxc, vtxc, v); + expect_same_result(total, semilocal, ROUNDING_TOLERANCE); +} + +/// XC_Functional::v_xc for the given functional and density with and without rVV10. +void expect_vxc_adds_rvv10(const char* functional, + const Charge& chr, + const UnitCell& ucell, + const bool domag, + const int gga_grad) +{ + XC_Functional::ScopedType scope(functional); + const int nrxx = chr.nrxx; + const int nspin = chr.nspin; + const bool domag_z = false; + const double hybrid_alpha = 0.0; + const double hse_omega = 0.0; + VxcResult total; + VxcResult semilocal; + for (int with_rvv10 = 0; with_rvv10 < 2; ++with_rvv10) + { + const bool rvv10 = with_rvv10 == 1; + RuntimeParameters parameters(rvv10); + VxcResult& result = rvv10 ? total : semilocal; + result = XC_Functional::v_xc(nrxx, &chr, &ucell, nspin, domag, domag_z, gga_grad, hybrid_alpha, hse_omega); + } + expect_rvv10_added(total, semilocal, chr); +} + +class NonlocalXCTest : public testing::Test +{ + protected: + /// grid on an oblique cell, distributed over the processes of the test + ModulePW::PW_Basis pw{"cpu", "double"}; + /// the same grid, complete on every process + ModulePW::PW_Basis serial{"cpu", "double"}; + + void SetUp() override + { + // The oblique cell catches row/column mistakes in the cell handed to libvdwxc. + const ModuleBase::Matrix3 lattice = oblique_lattice(1.0); + pw.initmpi(test_size, test_rank, MPI_COMM_WORLD); + serial.initmpi(1, 0, MPI_COMM_SELF); + setup(pw, lattice, BASIS_CUTOFF); + setup(serial, lattice, BASIS_CUTOFF); + } +}; + +TEST_F(NonlocalXCTest, MatchesDirectLibvdwxc) +{ + // unpolarized and spin-polarized layouts against libvdwxc called directly in the ABACUS axis order + const std::vector valence = local_grid(serial, valence_density); + const std::vector core = local_grid(serial, core_density); + const VxcResult result = rvv10_result(serial, {valence}, core); + + // reference: libvdwxc on the total density in the ABACUS axis order, Hartree units + std::vector total(serial.nrxx); + for (int ir = 0; ir < serial.nrxx; ++ir) + { + total[ir] = valence[ir] + core[ir]; + } + std::vector> rhog(serial.npw); + serial.real2recip(total.data(), rhog.data()); + std::vector> grad(serial.nrxx); + XC_Functional::grad_rho(rhog.data(), grad.data(), &serial, serial.tpiba); + std::vector sigma(serial.nrxx); + std::vector dedn(serial.nrxx); + std::vector dedsigma(serial.nrxx); + for (int ir = 0; ir < serial.nrxx; ++ir) + { + sigma[ir] = grad[ir].norm2(); + } + vdwxc_data direct = vdwxc_new(VDWXC_RVV10); + ASSERT_NE(direct, nullptr); + ASSERT_EQ(vdwxc_set_rvv10_parameters(direct, RVV10_B, RVV10_C), 0); + const ModuleBase::Matrix3 a = serial.latvec * serial.lat0; + const int nx = serial.nx; + const int ny = serial.ny; + const int nz = serial.nz; + vdwxc_set_unit_cell(direct, nx, ny, nz, a.e11, a.e12, a.e13, a.e21, a.e22, a.e23, a.e31, a.e32, a.e33); + vdwxc_init_serial(direct); + const double energy_ha = vdwxc_calculate(direct, total.data(), sigma.data(), dedn.data(), dedsigma.data()); + vdwxc_finalize(&direct); + for (int ir = 0; ir < serial.nrxx; ++ir) + { + grad[ir] *= 2.0 * dedsigma[ir]; + } + std::vector divergence(serial.nrxx); + XC_Functional::grad_dot(grad.data(), divergence.data(), &serial, serial.tpiba); + + ModuleBase::matrix v(1, serial.nrxx); + double vtxc = 0.0; + for (int ir = 0; ir < serial.nrxx; ++ir) + { + const double potential = ModuleBase::e2 * (dedn[ir] - divergence[ir]); + v(0, ir) = potential; + vtxc += potential * valence[ir]; + } + const double etxc = ModuleBase::e2 * energy_ha; + vtxc *= serial.omega / serial.nxyz; + const VxcResult expected(etxc, vtxc, v); + expect_same_result(result, expected, ROUNDING_TOLERANCE); + + // spin channels share the total density: same energy and vtxc, and each channel gets the potential + std::vector up(serial.nrxx); + std::vector down(serial.nrxx); + for (int ir = 0; ir < serial.nrxx; ++ir) + { + up[ir] = 0.7 * valence[ir]; + down[ir] = valence[ir] - up[ir]; + } + const VxcResult polarized = rvv10_result(serial, {up, down}, core); + EXPECT_NEAR(std::get<0>(polarized), etxc, ROUNDING_TOLERANCE); + EXPECT_NEAR(std::get<1>(polarized), vtxc, ROUNDING_TOLERANCE); + for (int ir = 0; ir < serial.nrxx; ++ir) + { + EXPECT_NEAR(std::get<2>(polarized)(0, ir), v(0, ir), ROUNDING_TOLERANCE); + EXPECT_NEAR(std::get<2>(polarized)(1, ir), v(0, ir), ROUNDING_TOLERANCE); + } +} + +TEST_F(NonlocalXCTest, DistributedGridMatchesCompleteGrid) +{ + const std::vector valence = local_grid(pw, valence_density); + const std::vector core = local_grid(pw, core_density); + const std::vector valence_complete = local_grid(serial, valence_density); + const std::vector core_complete = local_grid(serial, core_density); + const VxcResult distributed = rvv10_result(pw, {valence}, core); + const VxcResult complete = rvv10_result(serial, {valence_complete}, core_complete); + EXPECT_NEAR(std::get<0>(distributed), std::get<0>(complete), ROUNDING_TOLERANCE); + EXPECT_NEAR(std::get<1>(distributed), std::get<1>(complete), ROUNDING_TOLERANCE); + for (int ixy = 0; ixy < pw.nxy; ++ixy) + { + for (int iz = 0; iz < pw.nplane; ++iz) + { + const int local = ixy * pw.nplane + iz; + const int global = ixy * pw.nz + pw.startz_current + iz; + EXPECT_NEAR(std::get<2>(distributed)(0, local), std::get<2>(complete)(0, global), ROUNDING_TOLERANCE); + } + } + const ModuleBase::matrix stress = rvv10_stress(pw, {valence}, core); + const ModuleBase::matrix stress_complete = rvv10_stress(serial, {valence_complete}, core_complete); + EXPECT_GT(std::abs(stress_complete(0, 0)), 1.0e-10); + for (int a = 0; a < 3; ++a) + { + for (int b = 0; b < 3; ++b) + { + EXPECT_NEAR(stress(a, b), stress_complete(a, b), ROUNDING_TOLERANCE); + } + } +} + +TEST_F(NonlocalXCTest, PotentialIsTheDerivativeOfTheEnergy) +{ + const std::vector valence = local_grid(pw, valence_density); + const std::vector core = local_grid(pw, core_density); + const std::vector direction = local_grid(pw, perturbation); + std::vector plus = valence; + std::vector minus = valence; + for (int ir = 0; ir < pw.nrxx; ++ir) + { + plus[ir] += FINITE_DIFFERENCE_STEP * direction[ir]; + minus[ir] -= FINITE_DIFFERENCE_STEP * direction[ir]; + } + const VxcResult reference = rvv10_result(pw, {valence}, core); + const VxcResult result_plus = rvv10_result(pw, {plus}, core); + const VxcResult result_minus = rvv10_result(pw, {minus}, core); + const double energy_difference = std::get<0>(result_plus) - std::get<0>(result_minus); + const double finite_difference = energy_difference / (2.0 * FINITE_DIFFERENCE_STEP); + + double contraction = 0.0; + for (int ir = 0; ir < pw.nrxx; ++ir) + { + contraction += std::get<2>(reference)(0, ir) * direction[ir]; + } + Parallel_Common::reduce_data(&contraction, 1, pw.pool_world); + contraction *= pw.omega / pw.nxyz; + EXPECT_NEAR(finite_difference, contraction, FINITE_DIFFERENCE_TOLERANCE); +} + +TEST_F(NonlocalXCTest, LeftHandedCellGivesTheMirrorImageResult) +{ + // A mirror image has the same result; libvdwxc's signed cell volume would flip the sign of the energy. + ModulePW::PW_Basis mirror{"cpu", "double"}; + mirror.initmpi(1, 0, MPI_COMM_SELF); + const ModuleBase::Matrix3 mirrored = oblique_lattice(-1.0); + setup(mirror, mirrored, BASIS_CUTOFF); + ASSERT_LT(mirror.latvec.Det(), 0.0); + const std::vector valence = local_grid(serial, valence_density); + const std::vector core = local_grid(serial, core_density); + const VxcResult right_handed = rvv10_result(serial, {valence}, core); + const VxcResult left_handed = rvv10_result(mirror, {valence}, core); + expect_same_result(left_handed, right_handed, ROUNDING_TOLERANCE); +} + +TEST_F(NonlocalXCTest, WorkspaceFollowsGridCellAndProcessChanges) +{ + // The libvdwxc object is reused between calls; switching the process layout, the cell, the grid and + // back must reproduce the earlier results. + ModulePW::PW_Basis orthogonal{"cpu", "double"}; + orthogonal.initmpi(1, 0, MPI_COMM_SELF); + const ModuleBase::Matrix3 rectangle(std::sqrt(0.75), 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 2.0); + setup(orthogonal, rectangle, BASIS_CUTOFF); + ModulePW::PW_Basis finer{"cpu", "double"}; + finer.initmpi(1, 0, MPI_COMM_SELF); + const ModuleBase::Matrix3 lattice = oblique_lattice(1.0); + setup_grid(finer, lattice, BASIS_CUTOFF, 5, 7, 11); + const std::vector valence_finer = local_grid(finer, valence_density); + const std::vector core_finer = local_grid(finer, core_density); + const std::vector valence = local_grid(serial, valence_density); + const std::vector core = local_grid(serial, core_density); + const std::vector valence_local = local_grid(pw, valence_density); + const std::vector core_local = local_grid(pw, core_density); + const VxcResult first = rvv10_result(serial, {valence}, core); + const VxcResult repeated = rvv10_result(serial, {valence}, core); + const VxcResult distributed = rvv10_result(pw, {valence_local}, core_local); + const VxcResult other_cell = rvv10_result(orthogonal, {valence}, core); + const VxcResult other_grid = rvv10_result(finer, {valence_finer}, core_finer); + const VxcResult again = rvv10_result(serial, {valence}, core); + const VxcResult other_grid_again = rvv10_result(finer, {valence_finer}, core_finer); + expect_same_result(repeated, first, 0.0); + expect_same_result(again, first, ROUNDING_TOLERANCE); + expect_same_result(other_grid_again, other_grid, ROUNDING_TOLERANCE); + EXPECT_NEAR(std::get<0>(distributed), std::get<0>(first), ROUNDING_TOLERANCE); + EXPECT_GT(std::abs(std::get<0>(other_cell) - std::get<0>(first)), 1.0e-8); +} + +TEST_F(NonlocalXCTest, StressIsTheStrainDerivativeOfTheEnergy) +{ + // A strain keeps the grid and every G component of the density, so the stress is the exact derivative + // of the discrete energy: sigma_ab = -(1/omega) dE/d(eps_ab) at a fixed electron number. + ModulePW::PW_Basis box{"cpu", "double"}; + box.initmpi(1, 0, MPI_COMM_SELF); + const ModuleBase::Matrix3 lattice = oblique_lattice(1.0); + setup(box, lattice, LARGEST_SPHERE_CUTOFF); + const std::vector density = local_grid(box, valence_density); + const std::vector no_core(box.nrxx, 0.0); + const VxcResult reference = rvv10_result(box, {density}, no_core); + const ModuleBase::matrix extra = rvv10_stress(box, {density}, no_core); + // the XC diagonal term that ABACUS adds from etxc and vtxc + const double diagonal = -(std::get<0>(reference) - std::get<1>(reference)) / box.omega; + + const int components[4][2] = {{0, 0}, {2, 2}, {0, 1}, {1, 2}}; + for (const int* ab: components) + { + const int a = ab[0]; + const int b = ab[1]; + double energy[2] = {0.0, 0.0}; + for (int side = 0; side < 2; ++side) + { + const double h = (side == 0) ? STRAIN_STEP : -STRAIN_STEP; + ModulePW::PW_Basis strained{"cpu", "double"}; + strained.initmpi(1, 0, MPI_COMM_SELF); + const ModuleBase::Matrix3 strained_cell = strained_lattice(lattice, a, b, h); + setup(strained, strained_cell, LARGEST_SPHERE_CUTOFF); + // keep the electron number: the density scales with the inverse volume ratio + const double volume_ratio = strained.omega / box.omega; + std::vector scaled = density; + for (double& value: scaled) + { + value /= volume_ratio; + } + const VxcResult strained_result = rvv10_result(strained, {scaled}, no_core); + energy[side] = std::get<0>(strained_result); + } + const double finite_difference = -(energy[0] - energy[1]) / (2.0 * STRAIN_STEP * box.omega); + const double analytic = extra(a, b) + ((a == b) ? diagonal : 0.0); + const double tolerance = STRAIN_TOLERANCE * std::abs(analytic) + 1.0e-12; + SCOPED_TRACE(testing::Message() << "component " << a << b); + EXPECT_NEAR(analytic, finite_difference, tolerance); + } +} + +TEST_F(NonlocalXCTest, VxcAddsRvv10ToEveryFunctionalPath) +{ + // v_xc and v_xc_meta with rVV10 equal the semilocal result plus add_rvv10, for native PBE (final + // return), a libxc GGA (v_xc_libxc), noncollinear PBE (NCGGA return) and SCAN (v_xc_meta). + const std::vector valence = local_grid(pw, valence_density); + const std::vector core = local_grid(pw, core_density); + const std::vector zeros(pw.nrxx); + std::vector> collinear(2, zeros); + std::vector> noncollinear(4, zeros); + for (int ir = 0; ir < pw.nrxx; ++ir) + { + collinear[0][ir] = 0.7 * valence[ir]; + collinear[1][ir] = valence[ir] - collinear[0][ir]; + noncollinear[0][ir] = valence[ir]; + noncollinear[1][ir] = 0.1 * valence[ir]; + noncollinear[2][ir] = -0.05 * valence[ir]; + noncollinear[3][ir] = 0.2 * valence[ir]; + } + UnitCell ucell; + ucell.omega = pw.omega; + ucell.tpiba = pw.tpiba; + ucell.magnet.lsign_ = false; + const ChargeOnGrid polarized(pw, collinear, core); + const ChargeOnGrid magnetized(pw, noncollinear, core); + expect_vxc_adds_rvv10("PBE", polarized.chr, ucell, false, 0); + expect_vxc_adds_rvv10("GGA_X_RPW86+GGA_C_PBE", polarized.chr, ucell, false, 0); + expect_vxc_adds_rvv10("PBE", magnetized.chr, ucell, true, 2); + + XC_Functional::ScopedType scan("SCAN"); + const std::vector> unpolarized(1, valence); + std::vector kinetic = local_grid(pw, kinetic_density); + std::vector kinetic_rows(1, kinetic.data()); + ChargeOnGrid charge(pw, unpolarized, core); + charge.chr.kin_r = kinetic_rows.data(); + const std::vector func_id = XC_Functional::get_func_id(); + const int nspin = 1; + const double hybrid_alpha = 0.0; + const double hse_omega = 0.0; + VxcResult results[2]; + for (int with_rvv10 = 0; with_rvv10 < 2; ++with_rvv10) + { + const bool rvv10 = with_rvv10 == 1; + RuntimeParameters parameters(rvv10); + const MetaResult meta = XC_Functional_Libxc::v_xc_meta(func_id, + pw.nrxx, + pw.omega, + pw.tpiba, + &charge.chr, + nspin, + hybrid_alpha, + hse_omega); + results[with_rvv10] = VxcResult(std::get<0>(meta), std::get<1>(meta), std::get<2>(meta)); + } + expect_rvv10_added(results[1], results[0], charge.chr); +} + +} // namespace + +int main(int argc, char** argv) +{ + int threads = 1; + Parallel_Global::read_pal_param(argc, argv, test_size, threads, test_rank); + POOL_WORLD = MPI_COMM_WORLD; + KP_WORLD = MPI_COMM_NULL; + INT_BGROUP = MPI_COMM_NULL; + BP_WORLD = MPI_COMM_NULL; + GRID_WORLD = MPI_COMM_NULL; + DIAG_WORLD = MPI_COMM_NULL; + testing::InitGoogleTest(&argc, argv); + const int result = RUN_ALL_TESTS(); + Parallel_Global::finalize_mpi(); + return result; +} diff --git a/source/source_hamilt/module_vdw/test/vdw_test.cpp b/source/source_hamilt/module_vdw/test/vdw_test.cpp index 71e3682ce6..5e0105e198 100644 --- a/source/source_hamilt/module_vdw/test/vdw_test.cpp +++ b/source/source_hamilt/module_vdw/test/vdw_test.cpp @@ -264,6 +264,18 @@ TEST_F(vdwd2Test, WrongVdwType) std::string output = testing::internal::GetCapturedStdout(); } +TEST_F(vdwd2Test, Rvv10HasNoPairwiseTerm) +{ + input.vdw_method = "rvv10"; +#ifdef __LIBVDWXC + EXPECT_TRUE(vdw::make_vdw(ucell, input) == nullptr); +#else + testing::internal::CaptureStdout(); + EXPECT_EXIT(vdw::make_vdw(ucell, input);, ::testing::ExitedWithCode(1), ""); + std::string output = testing::internal::GetCapturedStdout(); +#endif +} + TEST_F(vdwd2Test, D2ReadFile) { input.vdw_C6_file = "c6.txt"; diff --git a/source/source_hamilt/module_vdw/vdw.cpp b/source/source_hamilt/module_vdw/vdw.cpp index 2bcad9de68..ac8f6441a4 100644 --- a/source/source_hamilt/module_vdw/vdw.cpp +++ b/source/source_hamilt/module_vdw/vdw.cpp @@ -94,6 +94,21 @@ std::unique_ptr make_vdw(const UnitCell& ucell, const Input_para& input, st "Rebuild ABACUS with -DENABLE_DFTD4=ON."); #endif } + else if (input.vdw_method == "rvv10") + { +#ifndef __LIBVDWXC + ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::make_vdw", + "rVV10 support was not enabled at build time. " + "Rebuild ABACUS with -DENABLE_LIBXC=ON -DENABLE_LIBVDWXC=ON."); +#endif + // rVV10 is evaluated inside XC_Functional, so there is no pairwise term here + if (plog != nullptr) + { + *plog << "\nrVV10 nonlocal correlation (libvdwxc): b=" << input.vdw_rvv10_b << " C=" << input.vdw_rvv10_c + << "\nincluded in the exchange-correlation energy, potential, forces and stress" << std::endl; + } + return nullptr; + } else if (input.vdw_method != "none") { ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::make_vdw", diff --git a/source/source_hamilt/module_xc/CMakeLists.txt b/source/source_hamilt/module_xc/CMakeLists.txt index 495c5ae199..8496738f86 100644 --- a/source/source_hamilt/module_xc/CMakeLists.txt +++ b/source/source_hamilt/module_xc/CMakeLists.txt @@ -31,6 +31,12 @@ if(ENABLE_COVERAGE) add_coverage(xc_) endif() +if(ENABLE_LIBVDWXC) + # The XC evaluators (v_xc, v_xc_libxc, v_xc_meta) and Stress_Func::stress_gga add rVV10 through this adapter. + target_sources(xc_ PRIVATE ../module_vdw/nonlocal_xc.cpp) + target_link_libraries(xc_ PUBLIC libvdwxc::vdwxc) +endif() + if(BUILD_TESTING) if(ENABLE_MPI) if(ENABLE_LIBXC) diff --git a/source/source_hamilt/module_xc/libxc_pot.cpp b/source/source_hamilt/module_xc/libxc_pot.cpp index eb4275cf5e..6ef84c6099 100644 --- a/source/source_hamilt/module_xc/libxc_pot.cpp +++ b/source/source_hamilt/module_xc/libxc_pot.cpp @@ -9,6 +9,9 @@ #ifdef __EXX #include "source_hamilt/module_xc/exx_info.h" #endif +#ifdef __LIBVDWXC +#include "source_hamilt/module_vdw/nonlocal_xc.h" +#endif #include @@ -449,6 +452,10 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( etxc *= omega / chr->rhopw->nxyz; vtxc *= omega / chr->rhopw->nxyz; +#ifdef __LIBVDWXC + const XCFunctionalParameters& xc_parameters = XC_Functional::get_runtime_parameters(); + vdw::add_rvv10(xc_parameters, *chr, etxc, vtxc, v); +#endif XC_Functional_Libxc::finish_func(funcs); @@ -768,6 +775,10 @@ std::tuple XC_Functional_Li etxc *= omega / chr->rhopw->nxyz; vtxc *= omega / chr->rhopw->nxyz; +#ifdef __LIBVDWXC + const XCFunctionalParameters& xc_parameters = XC_Functional::get_runtime_parameters(); + vdw::add_rvv10(xc_parameters, *chr, etxc, vtxc, v); +#endif XC_Functional_Libxc::finish_func(funcs); diff --git a/source/source_hamilt/module_xc/test/CMakeLists.txt b/source/source_hamilt/module_xc/test/CMakeLists.txt index 82e8188923..1b639bac29 100644 --- a/source/source_hamilt/module_xc/test/CMakeLists.txt +++ b/source/source_hamilt/module_xc/test/CMakeLists.txt @@ -164,3 +164,7 @@ foreach(xc_gradient_test MODULE_HAMILT_XCTest_GRADCORR MODULE_HAMILT_XCTest_LAPL ../xc_functional_ncgga_sf.cpp ../xc_ncgga_radial.cpp ../libxc_pot.cpp ../libxc_tools.cpp ../../../source_base/timer.cpp) endforeach() + +# The targets in this directory compile xc_pot.cpp or libxc_pot.cpp against a mocked PW basis, without the +# libvdwxc adapter. +abacus_disable_feature_definitions(__LIBVDWXC) diff --git a/source/source_hamilt/module_xc/xc_functional.cpp b/source/source_hamilt/module_xc/xc_functional.cpp index afb4aff039..79a7ecfdb6 100644 --- a/source/source_hamilt/module_xc/xc_functional.cpp +++ b/source/source_hamilt/module_xc/xc_functional.cpp @@ -59,6 +59,9 @@ void XC_Functional::ScopedType::restore() noexcept XC_Functional::runtime_parameters.exx_erfc_alpha.swap(runtime_parameters_.exx_erfc_alpha); XC_Functional::runtime_parameters.xc_exch_ext.swap(runtime_parameters_.xc_exch_ext); XC_Functional::runtime_parameters.xc_corr_ext.swap(runtime_parameters_.xc_corr_ext); + XC_Functional::runtime_parameters.rvv10 = runtime_parameters_.rvv10; + XC_Functional::runtime_parameters.rvv10_b = runtime_parameters_.rvv10_b; + XC_Functional::runtime_parameters.rvv10_c = runtime_parameters_.rvv10_c; XC_Functional::scaling_factor_xc.swap(scaling_factor_xc_); } diff --git a/source/source_hamilt/module_xc/xc_functional.h b/source/source_hamilt/module_xc/xc_functional.h index a8b57f0046..cce26f3c80 100644 --- a/source/source_hamilt/module_xc/xc_functional.h +++ b/source/source_hamilt/module_xc/xc_functional.h @@ -39,6 +39,10 @@ struct XCFunctionalParameters std::vector exx_erfc_alpha = {"default"}; std::vector xc_exch_ext; std::vector xc_corr_ext; + /// rVV10 nonlocal correlation (vdw_method rvv10) added to any functional, and its parameters b and C + bool rvv10 = false; + double rvv10_b = 0.0; + double rvv10_c = 0.0; }; class XC_Functional diff --git a/source/source_hamilt/module_xc/xc_pot.cpp b/source/source_hamilt/module_xc/xc_pot.cpp index 278de82754..21475732da 100644 --- a/source/source_hamilt/module_xc/xc_pot.cpp +++ b/source/source_hamilt/module_xc/xc_pot.cpp @@ -16,6 +16,9 @@ #include "source_hamilt/module_xc/exx_info.h" #endif #endif +#ifdef __LIBVDWXC +#include "source_hamilt/module_vdw/nonlocal_xc.h" +#endif // [etxc, vtxc, v] = XC_Functional::v_xc(...) @@ -56,7 +59,15 @@ std::tuple XC_Functional::v_xc( // Evaluate the regularized mapped-density local-collinear graph. if (nspin == 4 && (domag || domag_z) && gga_grad == 2) { - return ModuleXC::NCGGA_SF_Builtin::v_xc_ncgga_sf_builtin(nrxx, ucell->omega, ucell->tpiba, chr); + std::tuple etxc_vtxc_v + = ModuleXC::NCGGA_SF_Builtin::v_xc_ncgga_sf_builtin(nrxx, ucell->omega, ucell->tpiba, chr); +#ifdef __LIBVDWXC + double& etxc = std::get<0>(etxc_vtxc_v); + double& vtxc = std::get<1>(etxc_vtxc_v); + ModuleBase::matrix& v = std::get<2>(etxc_vtxc_v); + vdw::add_rvv10(runtime_parameters, *chr, etxc, vtxc, v); +#endif + return etxc_vtxc_v; } ModuleBase::timer::start("XC_Functional", "v_xc"); @@ -205,6 +216,9 @@ std::tuple XC_Functional::v_xc( #endif etxc *= ucell->omega / chr->rhopw->nxyz; vtxc *= ucell->omega / chr->rhopw->nxyz; +#ifdef __LIBVDWXC + vdw::add_rvv10(runtime_parameters, *chr, etxc, vtxc, v); +#endif ModuleBase::timer::end("XC_Functional", "v_xc"); return std::make_tuple(etxc, vtxc, std::move(v)); diff --git a/source/source_io/build_info.h.in b/source/source_io/build_info.h.in index 543ef3e1b4..a7eb5e8995 100644 --- a/source/source_io/build_info.h.in +++ b/source/source_io/build_info.h.in @@ -52,6 +52,7 @@ // --- Other Libraries --- #define ABACUS_DFTD4_VERSION "@ABACUS_DFTD4_VERSION@" +#define ABACUS_LIBVDWXC_VERSION "@ABACUS_LIBVDWXC_VERSION@" #define ABACUS_JSON_VERSION "@ABACUS_JSON_VERSION@" #define ABACUS_PEXSI_VERSION "@ABACUS_PEXSI_VERSION@" #define ABACUS_CNPY_VERSION "@ABACUS_CNPY_VERSION@" diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 1c624aa5a0..155415eb03 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -625,6 +625,8 @@ struct Input_para std::string vdw_cn_thr_unit = "Bohr"; ///< unit of cn_thr, Bohr or Angstrom std::string vdw_d4_xc = "default"; ///< functional name passed to DFT-D4 std::string vdw_d4_model = "d4"; ///< DFT-D4 dispersion model (d4 or d4s) + double vdw_rvv10_b = 6.3; ///< rVV10 parameter b + double vdw_rvv10_c = 0.0093; ///< rVV10 parameter C ModuleBase::Vector3 vdw_cutoff_period = {3, 3, 3}; ///< periods of periodic structure // ============== #Parameters (15.exx) ==================== diff --git a/source/source_io/module_parameter/read_inp_model.cpp b/source/source_io/module_parameter/read_inp_model.cpp index 961002946c..6d68178184 100644 --- a/source/source_io/module_parameter/read_inp_model.cpp +++ b/source/source_io/module_parameter/read_inp_model.cpp @@ -263,8 +263,12 @@ void ReadInput::item_model() // vdW Correction { Input_Item item("vdw_method"); - item.domain.set_choices({"none", "d2", "d3_0", "d3_bj", "d4"}); - item.annotation = "the method of calculating vdw (none ; d2 ; d3_0 ; d3_bj ; d4)"; + item.domain.set_choices({"none", "d2", "d3_0", "d3_bj", "d4", "rvv10"}); + // The rVV10 response kernel needed by DFPT and LR-TDDFT forces is not implemented. + item.domain.add_requirement( + "vdw_method==rvv10", + "esolver_type!=dfpt and (cal_force==false or (esolver_type!=lr and esolver_type!=ks-lr))"); + item.annotation = "the method of calculating vdw (none ; d2 ; d3_0 ; d3_bj ; d4 ; rvv10)"; item.category = "vdW correction"; item.type = "String"; item.description = R"(Specifies the method used for Van der Waals (VdW) correction. Available options are: @@ -272,9 +276,12 @@ void ReadInput::item_model() * d3_0: Grimme's DFT-D3(0) dispersion correction method (zero-damping) * d3_bj: Grimme's DFTD3(BJ) dispersion correction method (BJ-damping) * d4: Grimme's DFT-D4 dispersion correction method using the external DFT-D4 library +* rvv10: rVV10 nonlocal correlation evaluated self-consistently with the external libvdwxc library and added to the exchange-correlation energy, potential, forces and stress of dft_functional; see vdw_rvv10_b and vdw_rvv10_c * none: no vdW correction -[NOTE] ABACUS automatically loads DFT-D3 parameters for supported functionals according to dft_functional setting. Individual user values overwrite the corresponding tabulated values. Setting all four of vdw_s6, vdw_s8, vdw_a1 and vdw_a2 defines a fully custom set and bypasses functional lookup.)"; +[NOTE] ABACUS automatically loads DFT-D3 parameters for supported functionals according to dft_functional setting. Individual user values overwrite the corresponding tabulated values. Setting all four of vdw_s6, vdw_s8, vdw_a1 and vdw_a2 defines a fully custom set and bypasses functional lookup. + +[NOTE] vdw_method rvv10 does not support DFPT or LR-TDDFT forces. LR-TDDFT spectra use the kernel selected by xc_kernel, without the rVV10 response.)"; item.default_value = "none"; item.unit = ""; read_sync_string(input.vdw_method); @@ -309,6 +316,34 @@ Available options are: read_sync_string(input.vdw_d4_model); this->add_item(item); } + { + Input_Item item("vdw_rvv10_b"); + item.domain.set_rule("value>0"); + item.annotation = "rVV10 parameter b"; + item.category = "vdW correction"; + item.type = "Real"; + item.description + = R"(Parameter b of rVV10, which sets the short-range damping and depends on dft_functional. Published values are 6.3 for rPW86 exchange with PBE correlation (dft_functional GGA_X_RPW86+GGA_C_PBE, the original rVV10), 15.7 for SCAN and 11.95 for r2SCAN.)"; + item.default_value = "6.3"; + item.unit = ""; + item.set_availability("vdw_method==rvv10"); + read_sync_double(input.vdw_rvv10_b); + this->add_item(item); + } + { + Input_Item item("vdw_rvv10_c"); + item.domain.set_rule("value>=0"); + item.annotation = "rVV10 parameter C"; + item.category = "vdW correction"; + item.type = "Real"; + item.description = "Parameter C of rVV10, which sets the gradient dependence of the local polarizability; " + "0.0093 in the original rVV10 and in the SCAN and r2SCAN parametrizations."; + item.default_value = "0.0093"; + item.unit = ""; + item.set_availability("vdw_method==rvv10"); + read_sync_double(input.vdw_rvv10_c); + this->add_item(item); + } { Input_Item item("vdw_s6"); item.annotation = "scale parameter of d2/d3_0/d3_bj"; @@ -547,7 +582,8 @@ Namely, each line contains the element name and the corresponding parameter.)"; } }; item.check_value = [](const Input_Item& item, const Parameter& para) { - if (std::stod(para.input.vdw_cutoff_radius) <= 0 && para.input.vdw_method != "none") + if (std::stod(para.input.vdw_cutoff_radius) <= 0 && para.input.vdw_method != "none" + && para.input.vdw_method != "rvv10") { ModuleBase::WARNING_QUIT("ReadInput", "vdw_cutoff_radius <= 0 is not allowd"); } diff --git a/source/source_io/parse_args.cpp b/source/source_io/parse_args.cpp index 046bc52702..2115f20d21 100644 --- a/source/source_io/parse_args.cpp +++ b/source/source_io/parse_args.cpp @@ -165,6 +165,7 @@ void print_build_info() // --- 11. Other Libraries --- print_section("Other Libraries"); print_info("DFT-D4 Support", ABACUS_DFTD4_VERSION); + print_info("libvdwxc Support", ABACUS_LIBVDWXC_VERSION); print_info("JSON Support", ABACUS_JSON_VERSION); print_info("PEXSI Support", ABACUS_PEXSI_VERSION); print_info("cnpy Support", ABACUS_CNPY_VERSION); diff --git a/source/source_io/test_serial/read_input_test.cpp b/source/source_io/test_serial/read_input_test.cpp index 059e52c462..badce05377 100644 --- a/source/source_io/test_serial/read_input_test.cpp +++ b/source/source_io/test_serial/read_input_test.cpp @@ -212,6 +212,26 @@ TEST_F(InputTest, RejectAutoDevice) EXPECT_TRUE(std::remove("./auto_device_INPUT") == 0); } +TEST_F(InputTest, VdwMethodRvv10) +{ + const std::string input_file = testing::TempDir() + "rvv10_INPUT"; + Parameter param; + EXPECT_NO_THROW(read_parameters(input_file, + "vdw_method rvv10\ndft_functional scan\nvdw_rvv10_b 15.7\ncalculation cell-relax\n", + param)); + EXPECT_EQ(param.inp.vdw_method, "rvv10"); + EXPECT_DOUBLE_EQ(param.inp.vdw_rvv10_b, 15.7); + EXPECT_DOUBLE_EQ(param.inp.vdw_rvv10_c, 0.0093); + const std::string reason = "Requires esolver_type!=dfpt"; + expect_invalid_input(input_file, "vdw_method rvv10\nesolver_type dfpt\n", reason); + expect_invalid_input(input_file, "vdw_method rvv10\nesolver_type lr\ncalculation nscf\ncal_force 1\n", reason); + // relax resets cal_force to true, so LR-TDDFT relaxation is rejected without an explicit cal_force + expect_invalid_input(input_file, "vdw_method rvv10\nesolver_type ks-lr\ncalculation relax\ngamma_only 1\n", reason); + expect_invalid_input(input_file, "vdw_method rvv10\nvdw_rvv10_b 0\n", "vdw_rvv10_b"); + expect_invalid_input(input_file, "vdw_method rvv10\nvdw_rvv10_c -0.1\n", "vdw_rvv10_c"); + expect_invalid_input(input_file, "vdw_method RVV10\n", "vdw_method"); +} + TEST_F(InputTest, ValidateRelaxMethodScalar) { for (const std::string method : {"cg", "cg_trad", "bfgs", "bfgs_trad", "sd", "cg_bfgs", "lbfgs"}) diff --git a/source/source_pw/module_pwdft/stress_gga.cpp b/source/source_pw/module_pwdft/stress_gga.cpp index 1716bd093d..7a8a78a450 100644 --- a/source/source_pw/module_pwdft/stress_gga.cpp +++ b/source/source_pw/module_pwdft/stress_gga.cpp @@ -1,6 +1,9 @@ #include "stress_func.h" #include "source_base/parallel_reduce.h" #include "source_hamilt/module_xc/xc_functional.h" +#ifdef __LIBVDWXC +#include "source_hamilt/module_vdw/nonlocal_xc.h" +#endif //calculate the GGA stress correction in PW and LCAO template @@ -15,6 +18,11 @@ void Stress_Func::stress_gga(const UnitCell& ucell, { ModuleBase::TITLE("Stress","stress_gga"); ModuleBase::timer::start("Stress","stress_gga"); +#ifdef __LIBVDWXC + // rVV10 is added before the LDA return below; stress_gga is the XC stress hook of PW, LCAO, SDFT and OFDFT + const XCFunctionalParameters& xc_parameters = XC_Functional::get_runtime_parameters(); + vdw::add_rvv10_stress(xc_parameters, *chr, sigma); +#endif int func_type = XC_Functional::get_func_type(); if (func_type == 0 || func_type == 1)