From 0d0f29bf2d9e26174b44c3ad17679d8aa85ccdd4 Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Thu, 15 Oct 2026 23:56:03 +0800 Subject: [PATCH 01/24] experiment(lcao): connect device selection to LibRI DDLA GPU contractions --- CMakeLists.txt | 9 ++++++ cmake/modules/FindLibRI.cmake | 6 +++- source/CMakeLists.txt | 8 +++++ .../source_lcao/module_ri/exx_lri_interface.h | 10 +++++++ source/source_lcao/setup_exx.cpp | 30 +++++++++++++++++++ 5 files changed, 62 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 62843068f6d..92d1c6cbf80 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -606,6 +606,15 @@ endfunction() if(ENABLE_LIBRI) set_if_higher(CMAKE_CXX_STANDARD 14) + option(ENABLE_LIBRI_GPU "Enable experimental LibRI CUDA contractions through DDLA" OFF) + if(ENABLE_LIBRI_GPU) + if(NOT USE_CUDA) + message(FATAL_ERROR "ENABLE_LIBRI_GPU requires USE_CUDA") + endif() + set_if_higher(CMAKE_CXX_STANDARD 17) + find_path(DDLA_INCLUDE_DIR ddla/ddla_connector.h REQUIRED) + find_library(DDLA_LIBRARY NAMES ddla REQUIRED) + endif() find_package(LibRI REQUIRED) find_package(LibComm REQUIRED) find_package(cereal REQUIRED CONFIG) diff --git a/cmake/modules/FindLibRI.cmake b/cmake/modules/FindLibRI.cmake index fa7e8833f41..49f053f7c55 100644 --- a/cmake/modules/FindLibRI.cmake +++ b/cmake/modules/FindLibRI.cmake @@ -14,9 +14,13 @@ find_path(LIBRI_DIR if(NOT LIBRI_DIR) include(FetchContent) + set(_libri_revision b0eff7a) + if(ENABLE_LIBRI_GPU) + set(_libri_revision 6ce3c8480ae38a085746dc81cc4dbb59b8074430) + endif() FetchContent_Declare( LibRI - URL https://codeload.github.com/abacusmodeling/LibRI/tar.gz/b0eff7a + URL https://codeload.github.com/abacusmodeling/LibRI/tar.gz/${_libri_revision} ) FetchContent_Populate(LibRI) set(LIBRI_DIR ${libri_SOURCE_DIR}) diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index 9f7ab0196cc..c6d92c9fd22 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -60,6 +60,9 @@ set(_abacus_feature_definitions $<$:__FLOAT_FFTW> $<$:__MLALGO> $<$:__CNPY> + $<$:__GPU_RI> + $<$:__DDLA_RI> + $<$:DDLA_USE_CUDA> $<$:__EXX> $<$:EXX_DM=3> $<$:EXX_H_COMM=2> @@ -307,6 +310,11 @@ if(ENABLE_CNPY) endif() endif() +if(ENABLE_LIBRI_GPU) + list(APPEND _abacus_feature_include_dirs ${DDLA_INCLUDE_DIR}) + list(APPEND _abacus_feature_libs ${DDLA_LIBRARY}) +endif() + if(ENABLE_LIBRI) list(APPEND _abacus_feature_include_dirs ${LIBRI_DIR}/include diff --git a/source/source_lcao/module_ri/exx_lri_interface.h b/source/source_lcao/module_ri/exx_lri_interface.h index 8b5bc61cf64..4476d03e9e5 100644 --- a/source/source_lcao/module_ri/exx_lri_interface.h +++ b/source/source_lcao/module_ri/exx_lri_interface.h @@ -53,6 +53,16 @@ class Exx_LRI_Interface } Exx_LRI_Interface() = delete; + void set_gpu(const bool use_gpu) + { +#ifdef __GPU_RI + auto& mode = this->exx_ptr->exx_lri.lri.cal_mode; + mode = use_gpu ? RI::LRI_Cal_Mode::GPU : RI::LRI_Cal_Mode::CPU; +#else + (void)use_gpu; +#endif + } + ///// read and write Hexxs using cereal //void write_Hexxs_cereal(const std::string& file_name) const; //void read_Hexxs_cereal(const std::string& file_name); diff --git a/source/source_lcao/setup_exx.cpp b/source/source_lcao/setup_exx.cpp index 0c17b3d2d9b..fe4df875fb3 100644 --- a/source/source_lcao/setup_exx.cpp +++ b/source/source_lcao/setup_exx.cpp @@ -50,6 +50,21 @@ void Exx_NAO::before_runner( Exx_Info& exx_info) { #ifdef __EXX + const bool use_gpu = inp.device == "gpu"; + if (this->exd) + { + this->exd->set_gpu(use_gpu); + } + if (this->exc) + { + this->exc->set_gpu(use_gpu); + } +#ifdef __GPU_RI + if (exx_info.info_global.cal_exx) + { + GlobalV::ofs_running << "LibRI contractions: " << inp.device << std::endl; + } +#endif if (inp.calculation == "scf" || inp.calculation == "relax" || inp.calculation == "cell-relax" || inp.calculation == "md") { @@ -104,6 +119,21 @@ void Exx_NAO::before_scf( Exx_Info& exx_info) { #ifdef __EXX + const bool use_gpu = inp.device == "gpu"; + if (this->exd) + { + this->exd->set_gpu(use_gpu); + } + if (this->exc) + { + this->exc->set_gpu(use_gpu); + } +#ifdef __GPU_RI + if (exx_info.info_global.cal_exx) + { + GlobalV::ofs_running << "LibRI contractions: " << inp.device << std::endl; + } +#endif if (inp.calculation != "nscf") { const bool update_symmetry = inp.calculation == "cell-relax" From b408a31374fa6134618b60d075cb50d8cf1b1c0c Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Fri, 16 Oct 2026 23:17:01 +0800 Subject: [PATCH 02/24] refactor(lcao): select RI backend once before runner initialization --- source/source_lcao/setup_exx.cpp | 21 --------------------- 1 file changed, 21 deletions(-) diff --git a/source/source_lcao/setup_exx.cpp b/source/source_lcao/setup_exx.cpp index fe4df875fb3..6599470964a 100644 --- a/source/source_lcao/setup_exx.cpp +++ b/source/source_lcao/setup_exx.cpp @@ -59,12 +59,6 @@ void Exx_NAO::before_runner( { this->exc->set_gpu(use_gpu); } -#ifdef __GPU_RI - if (exx_info.info_global.cal_exx) - { - GlobalV::ofs_running << "LibRI contractions: " << inp.device << std::endl; - } -#endif if (inp.calculation == "scf" || inp.calculation == "relax" || inp.calculation == "cell-relax" || inp.calculation == "md") { @@ -119,21 +113,6 @@ void Exx_NAO::before_scf( Exx_Info& exx_info) { #ifdef __EXX - const bool use_gpu = inp.device == "gpu"; - if (this->exd) - { - this->exd->set_gpu(use_gpu); - } - if (this->exc) - { - this->exc->set_gpu(use_gpu); - } -#ifdef __GPU_RI - if (exx_info.info_global.cal_exx) - { - GlobalV::ofs_running << "LibRI contractions: " << inp.device << std::endl; - } -#endif if (inp.calculation != "nscf") { const bool update_symmetry = inp.calculation == "cell-relax" From d20dde74ba2bf3333ada64770f812514646aaf99 Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Sat, 17 Oct 2026 00:52:41 +0800 Subject: [PATCH 03/24] test(lcao): record HSE GPU backend parity performance and memory --- cmake/modules/FindLibRI.cmake | 7 +- tests/lcao_hse_gpu/README.md | 72 +++++++++++++ tests/lcao_hse_gpu/RESULTS.md | 140 ++++++++++++++++++++++++++ tests/lcao_hse_gpu/benchmark.py | 126 +++++++++++++++++++++++ tests/lcao_hse_gpu/build.sbatch.in | 88 ++++++++++++++++ tests/lcao_hse_gpu/config.ini | 50 +++++++++ tests/lcao_hse_gpu/libri-compat.patch | 82 +++++++++++++++ 7 files changed, 561 insertions(+), 4 deletions(-) create mode 100644 tests/lcao_hse_gpu/README.md create mode 100644 tests/lcao_hse_gpu/RESULTS.md create mode 100644 tests/lcao_hse_gpu/benchmark.py create mode 100644 tests/lcao_hse_gpu/build.sbatch.in create mode 100644 tests/lcao_hse_gpu/config.ini create mode 100644 tests/lcao_hse_gpu/libri-compat.patch diff --git a/cmake/modules/FindLibRI.cmake b/cmake/modules/FindLibRI.cmake index 49f053f7c55..fe41ddaf760 100644 --- a/cmake/modules/FindLibRI.cmake +++ b/cmake/modules/FindLibRI.cmake @@ -13,14 +13,13 @@ find_path(LIBRI_DIR ) if(NOT LIBRI_DIR) - include(FetchContent) - set(_libri_revision b0eff7a) if(ENABLE_LIBRI_GPU) - set(_libri_revision 6ce3c8480ae38a085746dc81cc4dbb59b8074430) + message(FATAL_ERROR "ENABLE_LIBRI_GPU requires an explicitly supplied, patched LIBRI_DIR; see tests/lcao_hse_gpu/README.md") endif() + include(FetchContent) FetchContent_Declare( LibRI - URL https://codeload.github.com/abacusmodeling/LibRI/tar.gz/${_libri_revision} + URL https://codeload.github.com/abacusmodeling/LibRI/tar.gz/b0eff7a ) FetchContent_Populate(LibRI) set(LIBRI_DIR ${libri_SOURCE_DIR}) diff --git a/tests/lcao_hse_gpu/README.md b/tests/lcao_hse_gpu/README.md new file mode 100644 index 00000000000..859756d3a76 --- /dev/null +++ b/tests/lcao_hse_gpu/README.md @@ -0,0 +1,72 @@ +# Experimental LCAO HSE GPU integration + +Fork-only experiment; not production-ready. ABACUS base: 8ae37fcd6. +Runtime selection uses the existing `device cpu/gpu`. Build with +`ENABLE_LIBRI=ON`, `USE_CUDA=ON`, `ENABLE_LIBRI_GPU=ON`, and explicit +`LIBRI_DIR`, `DDLA_INCLUDE_DIR`, `DDLA_LIBRARY`. + +Dependencies used in the SAI experiment: + +- LibRI 6ce3c8480ae38a085746dc81cc4dbb59b8074430; archive SHA256 + 09c8f766687e5c2236bfe580f372ea5ecef38e935a56c6a24a786d0430ff8fc4. +- LibDDLA 9c404ac65148ec1b48977b846f38ba6ef9e134f7; archive SHA256 + bc349e64c0a308aa9cd0f766e9a766e942bbc974f557522604bd4af9e1ce9baa. +- LibComm and cereal use the existing ABACUS toolchain pins. + +Apply `libri-compat.patch` to LibRI after normalizing its two affected files +to LF, then use `git apply --unidiff-zero /path/to/libri-compat.patch` +from that dependency root. The patch omits context to preserve upstream tab +indentation without introducing whitespace diagnostics in this repository. It repairs the DDLA diagnostic helper name, the split cleanup API, +and GPU contraction accumulation that otherwise ignores `fac_add_Ds`. +The original GPU force regression fails while energy passes; patched force +and stress must be verified against the unchanged references. +The build template records equivalent edits used in the actual experiment. + +Use the existing `.ci/slurm/runner.py` transport, build, worker and artifact +machinery. Experiment controls and full raw results are retained locally at +`/mnt/g/work/agent/abacus/lcao-hse-results/attempt*`. +The copied control runner supports `autotest_gpu` and `hse_benchmark`, as in +the earlier fork experiment; these extensions are not in the base runner. +Set OMP_NUM_THREADS=1 after loading cluster modules. + +Limitations: LibRI's 5 GiB temporary limit is not a total GPU memory budget. +Multi-node GPU mapping and all GPU contraction labels are not validated. +Do not infer production MD readiness from SCF smoke tests. + +## Benchmark interpretation + +The initial Si8/Si32 benchmark compares CPU ScaLAPACK with GPU cuSOLVER; +wall time therefore includes a solver change. Use the additional GPU-device +build with `ENABLE_LIBRI_GPU=OFF` to isolate RI offload. +`exx_hybrid_step=50` is an outer-loop cap, not a convergence guarantee. +Inspect the printed outer deltaE before calling a result converged. +GPU process memory is sampled at 100 ms; it is not a strict allocator peak +or a proof of an upper bound for other systems. + +## Running a backend comparison + +On a GPU allocation with the cluster modules loaded, set `OMP_NUM_THREADS=1` +and make DDLA's `lib` directory available in `LD_LIBRARY_PATH`. Build two +executables with the same patched LibRI and identical ABACUS options, changing +only `ENABLE_LIBRI_GPU`. Both runs use `device gpu` and `ks_solver cusolver`. + +```bash +python3 tests/lcao_hse_gpu/benchmark.py \ + --abacus /path/to/gpu-ri/abacus --baseline /path/to/cpu-ri/abacus \ + --source /path/to/abacus-source --output /path/to/new-results \ + --orbital Si_dojo_6au.orb --repeats 2 +``` + +The source must contain the test PP_ORB assets. `ABACUS_CPU_RI` may supply the +baseline path when invoked by the local `.ci/slurm` experiment controls. +The benchmark alternates backend order, checks outer energy convergence and +iteration count, and compares total energy, every force component and every +stress component. Thresholds are 1e-6 eV/atom, 1e-4 eV/Angstrom and 1e-3 kbar. +These are backend parity tests, not basis/cutoff convergence studies; the +small SZ system can break symmetry under prolonged outer iteration. + +GPU mode deliberately requires an externally supplied, patched LibRI tree. +The automatic download still uses ABACUS's original CPU dependency. Downloading +an unpatched newer revision is insufficient to reproduce the working build. +The build option defaults to OFF; the HSE approximation is unchanged. +See RESULTS.md for measured outcomes and remaining validation limits. diff --git a/tests/lcao_hse_gpu/RESULTS.md b/tests/lcao_hse_gpu/RESULTS.md new file mode 100644 index 00000000000..775b75ec11c --- /dev/null +++ b/tests/lcao_hse_gpu/RESULTS.md @@ -0,0 +1,140 @@ +# SAI LCAO HSE GPU experiment — 2026-10-17 + +Experimental integration only; no upstream PR. Computing resources were +provided by the Open Source Supercomputing Center of SAI (https://www.open-sai.com/). + +ABACUS base: `8ae37fcd6cb9043a4a0c7ebde2c062d97f393bdb`. +Runtime source: `b408a31374fa6134618b60d075cb50d8cf1b1c0c`. +Subsequent changes are test/documentation and a configure-time guard against +fetching an unpatched GPU dependency; runtime C++ is unchanged. +LibRI/DDLA revisions and patches are in README.md and libri-compat.patch. +Hardware: one V100 16 GB, one MPI rank, OMP_NUM_THREADS=1, CUDA 12.9.86, +GNU 13.3.0, OpenMPI. Both comparison executables use GPU cuSOLVER. +The CPU RI baseline uses the same newer LibRI, not the original b0eff7a. + +## Build and regression evidence + +| Job | Result | +| --- | --- | +| 1818882 | Patched GPU RI build passed | +| 1818916_0 / 1818916_2 | Existing HSE Gamma / spin-polarized Gamma energy, force, stress references passed | +| 1819337_0 | Existing complex multi-k HSE references passed with GPU cuSOLVER | +| 1819559 | CPU RI control build passed, CUDA and same newer LibRI retained | +| 1819665 | SZ Si8/Si32 matched-backend comparison passed | +| 1819761 | DZP Si8/Si32, two alternating-order repeats per backend: all 8 calculations and 4 paired comparisons passed | + +Failures retained in raw artifacts: unpatched GPU force accumulation failed +reference tests despite passing energy; incompatible DDLA helper and LibRI +cleanup APIs required fixes. The installed ELPA/GenELPA paths were unsuitable +for these GPU comparisons (including an illegal-instruction exit); cuSOLVER +is used consistently in the final matched comparison. Earlier failed jobs +are not counted as all-pass suites. + +With default `scf_ene_thr=-1`, the Si32 CPU RI and GPU RI trajectories had +matching early EXX energies but different outer stopping times. GPU RI reached +the 50-step cap and the final states differed. Explicit `scf_ene_thr=1e-8` eV +produced equal outer counts and matching results. This supports convergence +sensitivity as the cause of that discrepancy, but does not establish default +parameter robustness or production MD readiness. + +## SZ matched comparison + +| System | CPU RI wall s | GPU RI wall s | CPU/GPU sampled MiB | Outer steps | +| --- | ---: | ---: | ---: | ---: | +| Si8 SZ | 8.91 | 9.19 | 422 / 498 | 4 / 4 | +| Si32 SZ | 24.31 | 24.13 | 438 / 530 | 5 / 5 | + +Single paired measurements: no demonstrated speedup. Total energy differences +were below 4e-11 eV, forces identical at printed precision, maximum stress +component difference 9.98e-7 kbar. Si32 electronic EXX time increased from +0.50 to 0.80 s, force EXX from 0.80 to 1.17 s, stress EXX from 0.96 to 1.57 s. +Small contractions can be slower after offload. + +## DZP matched comparison + +Si_dojo_6au.orb is a 2s2p1d basis. Two repeats per backend, alternating +CPU/GPU order; means below include SCF, force and stress calculation. + +| System | CPU RI wall s | GPU RI wall s | Wall speedup | CPU/GPU sampled MiB | +| --- | ---: | ---: | ---: | ---: | +| Si8 DZP | 30.99 | 29.64 | 1.05x | 442 / 558–564 | +| Si32 DZP | 250.79 | 73.81 | 3.40x | 450 / 2756 | + +Si32 mean electronic EXX contraction time: 50.20 -> 5.36 s (9.36x). +EXX force: 64.95 -> 6.89 s; EXX stress: 85.09 -> 9.68 s. Both backends +used 7 outer iterations for both sizes. Maximum paired energy difference: +1.03e-12 eV/atom; force difference zero at printed precision; stress difference +1.11e-7 kbar. Symmetric silicon forces are zero, so nonzero-force coverage +comes from the existing HSE regression references, not these benchmark cells. + +These speedups are relative to single-thread CPU RI; a tuned multi-core CPU +comparison has not been performed. They are not a universal HSE speedup. +All eight memory logs contain a single GPU process PID. Si32 GPU RI added +2306 MiB (2.25 GiB) over the GPU-solver/CPU-RI baseline: a substantial increase, +not a guarantee of small memory overhead. If the requirement is less than +2x *total* GPU memory relative to that baseline, this implementation fails it. + +## Scope and limitations + +Memory is sampled every 100 ms and includes the GPU eigensolver. It is not +a strict allocation peak. LibRI's hardcoded 5 GiB temporary batching limit +excludes other simultaneous allocations, and is not a total memory bound. +No guarantee that large-system memory stays below twice a baseline is made. +Multi-GPU/multi-node, long MD, sustained allocation behavior, SOC and larger +production cells remain unvalidated. Do not enable this path by default. + +Local raw results: `/mnt/g/work/agent/abacus/lcao-hse-results/`. +The corresponding remote Slurm archives are under +`~/abacus_hse_gpu_experiment/archives/lcao-hse/`. +The experiment reused `.ci/slurm` transport/worker/artifact handling with local +controls; benchmark and dependency fixes are preserved here without importing +unrelated runner changes. + +## Follow-up priorities suggested by the evidence + +1. Keep the small SZ contraction workload on CPU unless profiling shows a + crossover. GPU dispatch/packing/transfer overhead can dominate these blocks. +2. Profile and optimize `refresh_coulomb`/initialization separately. Offloading + contraction GEMMs cannot accelerate this measured CPU-dominated phase. +3. LibRI currently constructs a GPU context per `cal_loop3` call and creates, + packs and uploads input/temporary/output objects per contraction label. + Reuse a single explicitly capped allocation arena and queue where lifetimes + allow it; do not keep two full workspaces or unbounded per-shape caches. +4. Investigate batching by shape and selective reuse of invariant tensors, + measuring transfer/packing/kernel time separately. Persistent tensors must + share a total memory budget with temporary and output buffers. + +These are proposed next experiments, not demonstrated speedups. The current +change connects an existing double-precision GPU contraction implementation; +it does not introduce mixed precision, new screening thresholds, or a different +HSE approximation. Floating-point accumulation order can differ (DDLA uses +atomic accumulation in relevant batched kernels), so bitwise identity is not +expected. Stronger screening or mixed precision would need separate accuracy +validation and are outside this experiment. + +## CP2K comparison + +The relevant transferable ideas are block batching and explicit memory-aware +contraction scheduling, not replacing an ABACUS matrix call with DBCSR directly. +CP2K's [RI HFX implementation](https://github.com/cp2k/cp2k/blob/0adc3f08c83fcdafa7d545917ddb4e901470c070/src/hfx_ri.F) +uses DBCSR operations and explicitly batches occupied-MO/AO indices to save +memory. DBCSR's [libsmm_acc](https://github.com/cp2k/dbcsr/blob/8f63c8f49ce807b5127cf8023b4684f579d19225/src/acc/libsmm_acc/README.md) +executes stacks of small matrix products and tunes kernels for matrix dimensions +and GPU models. LibRI already has variable-size batched contractions, so merely +adding batching is not a new optimization; block-shape distribution, packing, +transfer and kernel efficiency should determine the next change. Different +RI formulations/data structures preclude assuming CP2K's speedup transfers +unchanged to ABACUS. + +## Local checks + +- Repository `.ci/slurm/test_runner.py`: 45 tests passed. Running that suite + from the copied experiment controls initially failed because it assumes the + original 117-case configuration and repository layout; the repository suite + was rerun in its intended location without modifying expectations. +- `agent_governance_check.py --staged`: no findings. +- Python benchmark syntax/CLI and shell template syntax checks passed. +- LibRI patch apply check passed after LF normalization (`--unidiff-zero`). +- Missing GPU LibRI dependency produces the expected CMake diagnostic; this + guard was checked remotely with CMake 3.31.6. +- Final staged whitespace check passed. diff --git a/tests/lcao_hse_gpu/benchmark.py b/tests/lcao_hse_gpu/benchmark.py new file mode 100644 index 00000000000..4527e527406 --- /dev/null +++ b/tests/lcao_hse_gpu/benchmark.py @@ -0,0 +1,126 @@ +#!/usr/bin/env python3 +"""Compare two RI backends using identical GPU LCAO HSE inputs.""" +import argparse +import json +import os +import re +import subprocess +import time +from pathlib import Path + + +def derivatives(log, size): + """Read the final total force (eV/Angstrom) and stress (kbar).""" + force_block = log.rsplit('#TOTAL-FORCE (eV/Angstrom)#', 1)[1] + forces = re.findall(r'^\s*Si\d+\s+([-+\d.eE]+)\s+([-+\d.eE]+)\s+([-+\d.eE]+)', force_block, re.M) + stress_block = log.rsplit('#TOTAL-STRESS (kbar)#', 1)[1] + stress = re.findall(r'^\s*([-+\d.eE]+)\s+([-+\d.eE]+)\s+([-+\d.eE]+)\s*$', stress_block, re.M)[:3] + if len(forces) != size or len(stress) != 3: + raise RuntimeError('Missing force/stress output') + return {key: [float(x) for row in values for x in row] + for key, values in [('force', forces), ('stress', stress)]} + + +def write_inputs(work, source, size, orbital): + points = [(0, 0, 0), (0, .5, .5), (.5, 0, .5), (.5, .5, 0), + (.25, .25, .25), (.25, .75, .75), (.75, .25, .75), (.75, .75, .25)] + nx = 2 if size == 32 else 1 + positions = [((x+i)/nx, (y+j)/nx, z) for i in range(nx) + for j in range(nx) for x, y, z in points] + work.mkdir() + (work / 'STRU').write_text( + 'ATOMIC_SPECIES\nSi 28.0855 Si_dojo_nsoc.upf\nNUMERICAL_ORBITAL\n' + f'{orbital}\nLATTICE_CONSTANT\n10.26\nLATTICE_VECTORS\n' + f'{nx} 0 0\n0 {nx} 0\n0 0 1\nATOMIC_POSITIONS\nDirect\nSi\n0\n{size}\n' + + ''.join('%g %g %g 0 0 0\n' % p for p in positions)) + (work / 'KPT').write_text('K_POINTS\n0\nGamma\n1 1 1 0 0 0\n') + values = dict( + suffix='bench', calculation='scf', basis_type='lcao', device='gpu', + ks_solver='cusolver', gamma_only=1, nspin=1, ntype=1, nbands=int(size*2.5), + pseudo_dir=str(source/'tests/PP_ORB'), orbital_dir=str(source/'tests/PP_ORB'), + ecutwfc=30, scf_thr='1e-9', scf_ene_thr='1e-8', scf_nmax=200, symmetry=0, + smearing_method='gaussian', smearing_sigma=.002, mixing_beta=.1, + dft_functional='hse', exx_real_number=1, exx_pca_threshold='1e-3', + exx_ccp_rmesh_times=1, exx_separate_loop=1, exx_hybrid_step=50, + cal_force=1, cal_stress=1) + (work/'INPUT').write_text('INPUT_PARAMETERS\n' + ''.join(f'{k} {v}\n' for k, v in values.items())) + + +def run_case(binary, work, size): + with (work/'gpu-memory.csv').open('w') as memory: + monitor = subprocess.Popen( + ['nvidia-smi', '--query-compute-apps=pid,used_gpu_memory', + '--format=csv,noheader,nounits', '-lms', '100'], + stdout=memory, stderr=subprocess.STDOUT) + start = time.monotonic() + try: + with (work/'stdout.log').open('w') as stdout: + result = subprocess.run( + ['timeout', '--kill-after=30s', '1200s', 'mpirun', '-np', '1', str(binary)], + cwd=work, stdout=stdout, stderr=subprocess.STDOUT, check=False) + finally: + monitor.terminate() + monitor.wait() + seconds = time.monotonic() - start + if result.returncode: + raise RuntimeError(f'{work}: exit code {result.returncode}') + log = (work/'OUT.bench/running_scf.log').read_text() + stdout = (work/'stdout.log').read_text() + energy = re.findall(r'!FINAL_ETOT_IS\s+([-+\d.eE]+)', log+'\n'+stdout) + deltas = [float(x) for x in re.findall(r'deltaE \(eV\) from outer loop:\s*([-+\d.eE]+)', stdout)] + # An inner SCF marker alone does not establish outer EXX convergence. + if not energy or '#SCF IS CONVERGED#' not in log or not deltas or deltas[-1] >= 1e-8 or len(deltas) >= 50: + raise RuntimeError(f'{work}: outer energy convergence not established') + timing = {cls+'.'+fn: float(value) for cls, fn, value in re.findall( + r'^\s*(Exx_LRI|RI_2D_Comm|HSolverLCAO)\s+(\w+)\s+([\d.]+)', stdout, re.M)} + memory = [int(x) for x in re.findall(r'^\s*\d+,\s*(\d+)\s*$', (work/'gpu-memory.csv').read_text(), re.M)] + return dict(seconds=seconds, energy_ev=float(energy[-1]), outer_steps=len(deltas), + outer_delta_ev=deltas[-1], timing=timing, + sampled_max_process_mib=max(memory) if memory else None, + **derivatives(log, size)) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument('--abacus', type=Path, required=True, help='GPU RI executable') + parser.add_argument('--baseline', type=Path, default=os.environ.get('ABACUS_CPU_RI'), + help='CPU RI executable, built with the same LibRI and GPU solver') + parser.add_argument('--source', type=Path, required=True) + parser.add_argument('--output', type=Path, required=True) + parser.add_argument('--orbital', default='Si_dojo_6au_sz.orb') + parser.add_argument('--repeats', type=int, default=2) + args = parser.parse_args() + if args.baseline is None or args.repeats < 1: + parser.error('A baseline executable and positive repeat count are required') + args.output.mkdir(parents=True, exist_ok=True) + rows = [] + comparisons = [] + for size in (8, 32): + for repeat in range(args.repeats): + pair = {} + order = ('cpu_ri', 'gpu_ri') if repeat % 2 == 0 else ('gpu_ri', 'cpu_ri') + for backend in order: + work = args.output/f'si{size}_{backend}_{repeat}' + write_inputs(work, args.source, size, args.orbital) + binary = args.baseline if backend == 'cpu_ri' else args.abacus + row = dict(size=size, repeat=repeat, backend=backend, + orbital=args.orbital, **run_case(binary, work, size)) + rows.append(row) + pair[backend] = row + (args.output/'measurements.json').write_text(json.dumps(rows, indent=2)+'\n') + print(json.dumps({k: v for k, v in row.items() if k not in ('force', 'stress')}), flush=True) + cpu, gpu = pair['cpu_ri'], pair['gpu_ri'] + comparison = dict(size=size, repeat=repeat, + energy_ev_per_atom=abs(cpu['energy_ev']-gpu['energy_ev'])/size, + force_ev_angstrom=max(abs(x-y) for x, y in zip(cpu['force'], gpu['force'])), + stress_kbar=max(abs(x-y) for x, y in zip(cpu['stress'], gpu['stress'])), + same_outer_steps=cpu['outer_steps'] == gpu['outer_steps']) + comparisons.append(comparison) + (args.output/'comparison.json').write_text(json.dumps(comparisons, indent=2)+'\n') + if (comparison['energy_ev_per_atom'] > 1e-6 or comparison['force_ev_angstrom'] > 1e-4 + or comparison['stress_kbar'] > 1e-3 or not comparison['same_outer_steps']): + raise RuntimeError(f'Backend comparison failed: {comparison}') + + +if __name__ == '__main__': + main() diff --git a/tests/lcao_hse_gpu/build.sbatch.in b/tests/lcao_hse_gpu/build.sbatch.in new file mode 100644 index 00000000000..a8b90abc51b --- /dev/null +++ b/tests/lcao_hse_gpu/build.sbatch.in @@ -0,0 +1,88 @@ +#!/usr/bin/env bash +#SBATCH --job-name=@JOB_NAME@ +#SBATCH --partition=@PARTITION@ +#SBATCH --qos=@QOS@ +#SBATCH --account=abacus-group +#SBATCH --nodes=@NODES@ +#SBATCH --ntasks=@TASKS@ +#SBATCH --ntasks-per-node=@TASKS_PER_NODE@ +#SBATCH --gpus-per-node=@GPUS_PER_NODE@ +#SBATCH --time=@TIME@ +#SBATCH --output=@OUTPUT@ +#SBATCH --export=NIL + +set -euo pipefail + +export HOME=@HOME@ +export USER=${SLURM_JOB_USER:?} +export LOGNAME=$USER +export PATH=/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin + +source @CONTROL@/modules.sh +module -t list 2>&1 | tee @RESULTS@/modules.txt +command -v cmake make mpicxx nvcc | tee @RESULTS@/tools.txt +nvidia-smi dmon -s pucvmte -o T > @RESULTS@/dmon-${SLURM_JOB_ID}.log & + +# Header-only dependencies pinned by toolchain/scripts/package_versions.sh. +# Reuse the cereal installation settings from stage4/install_cereal.sh. +for package in libri libcomm cereal; do + mkdir -p @BUILD@/hse-deps/$package + tar -xzf @CONTROL@/$package.tar.gz --strip-components=1 -C @BUILD@/hse-deps/$package +done +cmake -S @BUILD@/hse-deps/cereal -B @BUILD@/hse-deps/cereal-build \ + -DCMAKE_INSTALL_PREFIX=@BUILD@/hse-deps/cereal-install \ + -DCMAKE_INSTALL_LIBDIR=lib -DBUILD_DOC=OFF -DBUILD_SANDBOX=OFF \ + -DBUILD_TESTS=OFF -DTHREAD_SAFE=ON -DSKIP_PORTABILITY_TEST=ON \ + -DSKIP_PERFORMANCE_COMPARISON=ON -DJUST_INSTALL_CEREAL=ON +cmake --install @BUILD@/hse-deps/cereal-build + +mkdir -p @BUILD@/ddla-source + tar -xzf @CONTROL@/ddla.tar.gz --strip-components=1 -C @BUILD@/ddla-source +cmake -S @BUILD@/ddla-source -B @BUILD@/ddla-build \ + -DCMAKE_BUILD_TYPE=Release -DCMAKE_CUDA_ARCHITECTURES=70 \ + -DCUBLAS_LIBRARY=/opt/devtools/nvidia/cuda-12.9.1/lib64/libcublas.so \ + -DCURAND_LIBRARY=/opt/devtools/nvidia/cuda-12.9.1/lib64/libcurand.so \ + -DCUSOLVER_LIBRARY=/opt/devtools/nvidia/cuda-12.9.1/lib64/libcusolver.so \ + -DDDLA_USE_CUDA=ON -DDDLA_USE_CCL=OFF -DBUILD_TESTS=OFF \ + -DCMAKE_INSTALL_PREFIX=@INSTALL@/ddla +cmake --build @BUILD@/ddla-build --parallel 16 +cmake --install @BUILD@/ddla-build +# LibRI adapter uses an obsolete DDLA diagnostic helper name. +sed -i 's/DEVICE_CHECK_RI/DEVICE_CHECK/g' @BUILD@/hse-deps/libri/include/RI/global/gpu/GPU_Backend.h +sed -i 's/LRI_Cal_Aux::destroy_lock_result(lock_Ds_result_add_map, Ds_result);/LRI_Cal_Aux::destroy_lock(lock_Ds_result_add_map); LRI_Cal_Aux::destroy_Ds_result(Ds_result);/' @BUILD@/hse-deps/libri/include/RI/ri/LRI-cal_loop3_GPU.hpp +# Accumulate into the existing local map, then apply fac_add_Ds once at the final merge. +sed -i 's/rDs_output.download(Ds_result, queue);/rDs_output.download(Ds_result_thread, queue);/g' @BUILD@/hse-deps/libri/include/RI/ri/LRI-cal_loop3_GPU.hpp +cmake -S @SOURCE@ -B @BUILD@ -G "Unix Makefiles" \ + -DCMAKE_BUILD_TYPE=Release \ + -DCMAKE_INSTALL_PREFIX=@INSTALL@ \ + -DCMAKE_CUDA_ARCHITECTURES=70 \ + -DENABLE_LIBRI_GPU=ON \ + -DDDLA_INCLUDE_DIR=@INSTALL@/ddla/include \ + -DDDLA_LIBRARY=@INSTALL@/ddla/lib/libddla.so \ + -DENABLE_LIBXC=ON -DENABLE_LIBRI=ON \ + -DLIBRI_DIR=@BUILD@/hse-deps/libri \ + -DLIBCOMM_DIR=@BUILD@/hse-deps/libcomm \ + -Dcereal_DIR=@BUILD@/hse-deps/cereal-install/lib/cmake/cereal \ + -DUSE_CUDA=ON \ + -DUSE_CUDA_MPI=ON \ + -DENABLE_CUSOLVERMP=ON \ + -DENABLE_CUBLASMP=ON \ + -DENABLE_NCCL_PARALLEL_DEVICE=ON \ + -DBUILD_TESTING=OFF \ + -DGIT_SUBMODULE=OFF 2>&1 | tee @RESULTS@/configure.log +cmake --build @BUILD@ --parallel 32 2>&1 | tee @RESULTS@/build.log +cmake --install @BUILD@ 2>&1 | tee @RESULTS@/install.log + +@INSTALL@/bin/abacus --info | tee @RESULTS@/abacus-info.txt +ldd @INSTALL@/bin/abacus | tee @RESULTS@/ldd.txt +if grep -q 'not found' @RESULTS@/ldd.txt; then + echo "ERROR: unresolved dynamic dependencies" >&2 + sed -n 's/.*not found.*/&/p' @RESULTS@/ldd.txt >&2 + exit 1 +fi +cp @BUILD@/CMakeCache.txt @RESULTS@/CMakeCache.txt + +mkdir -p "$HOME/abacus_hse_gpu_experiment/reuse/lcao" +cp @INSTALL@/bin/abacus "$HOME/abacus_hse_gpu_experiment/reuse/lcao/abacus" +cp -a @INSTALL@/ddla "$HOME/abacus_hse_gpu_experiment/reuse/lcao/" +exit 0 diff --git a/tests/lcao_hse_gpu/config.ini b/tests/lcao_hse_gpu/config.ini new file mode 100644 index 00000000000..a47520fee7d --- /dev/null +++ b/tests/lcao_hse_gpu/config.ini @@ -0,0 +1,50 @@ +[site] +name = Open Source Supercomputing Center of SAI +url = https://www.open-sai.com/ +acknowledgement = Computing resources were provided by + +[remote] +host = c0.sai.ai-4s.com +port = 12022 +user = abacususer01 +project_root = ~/abacus_hse_gpu_experiment +allowed_project_roots = /home, /org + +[cluster] +partition = 16V100 +mapping_root = /opt/sai_config/mps_mapping.d +disable_nccl_ib = false +poll_seconds = 10 + +[build] +qos = huge-gpu +nodes = 1 +tasks_per_node = 4 +gpus_per_node = 4 +time_seconds = 3600 + +[resource.gpu1] +qos = flood-1o2gpu +nodes = 1 +tasks_per_node = 1 +gpus_per_node = 1 +time_seconds = 1800 +parallelism = 1 + +[case.001] +suite = 08_RI +name = scf_hse_gamma +resource = gpu1 +runner = autotest_gpu + +[case.002] +suite = 08_RI +name = scf_hse_complex_multik +resource = gpu1 +runner = autotest_gpu + +[case.003] +suite = 08_RI +name = scf_hse_spin2_gamma +resource = gpu1 +runner = autotest_gpu diff --git a/tests/lcao_hse_gpu/libri-compat.patch b/tests/lcao_hse_gpu/libri-compat.patch new file mode 100644 index 00000000000..3e17c0c3a93 --- /dev/null +++ b/tests/lcao_hse_gpu/libri-compat.patch @@ -0,0 +1,82 @@ +--- a/include/RI/global/gpu/GPU_Backend.h ++++ b/include/RI/global/gpu/GPU_Backend.h +@@ -110 +110 @@ +- ddla::DEVICE_CHECK_RI(ddla::deviceGetDeviceCount(&device_count)); ++ ddla::DEVICE_CHECK(ddla::deviceGetDeviceCount(&device_count)); +@@ -177 +177 @@ +- ddla::DEVICE_CHECK_RI(ddla::deviceMalloc(pointer, count * sizeof(T))); ++ ddla::DEVICE_CHECK(ddla::deviceMalloc(pointer, count * sizeof(T))); +@@ -189 +189 @@ +- ddla::DEVICE_CHECK_RI(ddla::deviceFree(pointer)); ++ ddla::DEVICE_CHECK(ddla::deviceFree(pointer)); +@@ -202 +202 @@ +- ddla::DEVICE_CHECK_RI(deviceMemcpyAsync( ++ ddla::DEVICE_CHECK(deviceMemcpyAsync( +@@ -218 +218 @@ +- ddla::DEVICE_CHECK_RI(deviceMemcpyAsync( ++ ddla::DEVICE_CHECK(deviceMemcpyAsync( +@@ -232 +232 @@ +- ddla::DEVICE_CHECK_RI(deviceMemsetAsync( ++ ddla::DEVICE_CHECK(deviceMemsetAsync( +@@ -242 +242 @@ +- ddla::DEVICE_CHECK_RI( ++ ddla::DEVICE_CHECK( +--- a/include/RI/ri/LRI-cal_loop3_GPU.hpp ++++ b/include/RI/ri/LRI-cal_loop3_GPU.hpp +@@ -205 +205 @@ +- rDs_output.download(Ds_result, queue); ++ rDs_output.download(Ds_result_thread, queue); +@@ -341 +341 @@ +- rDs_output.download(Ds_result, queue); ++ rDs_output.download(Ds_result_thread, queue); +@@ -483 +483 @@ +- rDs_output.download(Ds_result, queue); ++ rDs_output.download(Ds_result_thread, queue); +@@ -623 +623 @@ +- rDs_output.download(Ds_result, queue); ++ rDs_output.download(Ds_result_thread, queue); +@@ -763 +763 @@ +- rDs_output.download(Ds_result, queue); ++ rDs_output.download(Ds_result_thread, queue); +@@ -903 +903 @@ +- rDs_output.download(Ds_result, queue); ++ rDs_output.download(Ds_result_thread, queue); +@@ -1045 +1045 @@ +- rDs_output.download(Ds_result, queue); ++ rDs_output.download(Ds_result_thread, queue); +@@ -1185 +1185 @@ +- rDs_output.download(Ds_result, queue); ++ rDs_output.download(Ds_result_thread, queue); +@@ -1325 +1325 @@ +- rDs_output.download(Ds_result, queue); ++ rDs_output.download(Ds_result_thread, queue); +@@ -1465 +1465 @@ +- rDs_output.download(Ds_result, queue); ++ rDs_output.download(Ds_result_thread, queue); +@@ -1606 +1606 @@ +- rDs_output.download(Ds_result, queue); ++ rDs_output.download(Ds_result_thread, queue); +@@ -1745 +1745 @@ +- rDs_output.download(Ds_result, queue); ++ rDs_output.download(Ds_result_thread, queue); +@@ -1884 +1884 @@ +- rDs_output.download(Ds_result, queue); ++ rDs_output.download(Ds_result_thread, queue); +@@ -2023 +2023 @@ +- rDs_output.download(Ds_result, queue); ++ rDs_output.download(Ds_result_thread, queue); +@@ -2189 +2189 @@ +- rDs_output.download(Ds_result, queue); ++ rDs_output.download(Ds_result_thread, queue); +@@ -2353 +2353 @@ +- rDs_output.download(Ds_result, queue); ++ rDs_output.download(Ds_result_thread, queue); +@@ -2517 +2517 @@ +- rDs_output.download(Ds_result, queue); ++ rDs_output.download(Ds_result_thread, queue); +@@ -2681 +2681 @@ +- rDs_output.download(Ds_result, queue); ++ rDs_output.download(Ds_result_thread, queue); +@@ -2693 +2693 @@ +- LRI_Cal_Aux::destroy_lock_result(lock_Ds_result_add_map, Ds_result); ++ LRI_Cal_Aux::destroy_lock(lock_Ds_result_add_map); LRI_Cal_Aux::destroy_Ds_result(Ds_result); From 784610dfadc952bfbefd98c298e250814c0e9684 Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Sat, 17 Oct 2026 22:31:40 +0800 Subject: [PATCH 04/24] docs(lcao): record successful multi-GPU RI memory scaling --- tests/lcao_hse_gpu/RESULTS.md | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/tests/lcao_hse_gpu/RESULTS.md b/tests/lcao_hse_gpu/RESULTS.md index 775b75ec11c..9caf680bd02 100644 --- a/tests/lcao_hse_gpu/RESULTS.md +++ b/tests/lcao_hse_gpu/RESULTS.md @@ -138,3 +138,37 @@ unchanged to ABACUS. - Missing GPU LibRI dependency produces the expected CMake diagnostic; this guard was checked remotely with CMake 3.31.6. - Final staged whitespace check passed. + +## Single-node multi-GPU memory scaling + +Follow-up Slurm job 1832072 passed on node 16v100n48. The same saved GPU RI +binary and Si32 DZP input ran sequentially with 1, 2 and 4 MPI ranks, one rank +per GPU, OMP_NUM_THREADS=1, GPU cuSOLVER. All runs used 7 outer iterations. +One measurement per configuration, on the same four-V100 allocation: + +| Active GPUs / ranks | Per-process sampled peak MiB | Wall seconds | +| --- | --- | ---: | +| 1 | 2756 | 63.78 | +| 2 | 1682, 1682 | 57.23 | +| 4 | 1158, 1106, 1106, 1158 | 52.71 | + +A live nvidia-smi UUID/PID snapshot confirmed the four listed processes were +on four distinct GPUs; a separate live query confirmed distinct GPUs for the +two-rank run. Peaks are sampled at 100 ms, not allocator-level bounds. +Maximum total-energy difference from one rank: 5.10e-11 eV. Forces agree at +printed precision; maximum stress-component difference: 1.30e-7 kbar. + +Thus distributing RI work *does* lower per-card memory in this tested case: +39% reduction with two cards and 58% with four. Total replicated memory grows, +and the slowest card does not follow ideal 1/N scaling. Performance improved +only 1.21x from one to four cards on this small system. These measurements +support using more GPUs to relieve capacity pressure, but do not establish +large-system or multi-node scaling or guarantee arbitrary capacity growth. +The CPU-RI baseline's small GPU footprint is not an appropriate standalone +rejection criterion for GPU offload; the relevant capacity metric is maximum +per-card memory as system size and GPU count increase. + +Controls, full input/output, per-PID memory traces and UUID mapping snapshot: +`/mnt/g/work/agent/abacus/lcao-hse-results/multigpu/`. +Remote archive: `archives/lcao-hse/lcao-multigpu-20261017-1.tar.gz` beneath the +experiment project root. No runtime code was changed for this test. From fa9c8531b08d13853157cf595d426f2ce072480c Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Sat, 17 Oct 2026 22:47:06 +0800 Subject: [PATCH 05/24] docs(lcao): verify size-dependent GPU RI memory and strong scaling --- tests/lcao_hse_gpu/RESULTS.md | 5 ++ tests/lcao_hse_gpu/SCALING.md | 130 ++++++++++++++++++++++++++++++++++ 2 files changed, 135 insertions(+) create mode 100644 tests/lcao_hse_gpu/SCALING.md diff --git a/tests/lcao_hse_gpu/RESULTS.md b/tests/lcao_hse_gpu/RESULTS.md index 9caf680bd02..667793bd1d2 100644 --- a/tests/lcao_hse_gpu/RESULTS.md +++ b/tests/lcao_hse_gpu/RESULTS.md @@ -172,3 +172,8 @@ Controls, full input/output, per-PID memory traces and UUID mapping snapshot: `/mnt/g/work/agent/abacus/lcao-hse-results/multigpu/`. Remote archive: `archives/lcao-hse/lcao-multigpu-20261017-1.tar.gz` beneath the experiment project root. No runtime code was changed for this test. + +Further size-dependent memory and multi-GPU checks (Si16/Si64, job 1832165) +are documented in [SCALING.md](SCALING.md). They reveal a memory plateau +between one and two cards for Si64; do not extrapolate the Si32 reduction +as a universal 1/P law. diff --git a/tests/lcao_hse_gpu/SCALING.md b/tests/lcao_hse_gpu/SCALING.md new file mode 100644 index 00000000000..9bed36b355b --- /dev/null +++ b/tests/lcao_hse_gpu/SCALING.md @@ -0,0 +1,130 @@ +# LCAO HSE GPU memory and strong-scaling verification + +## Measurements + +Slurm job 1832165 passed on SAI node 16v100n41 (four V100 16 GB GPUs). +Si16 and Si64 each ran with 1, 2 and 4 MPI ranks, one rank per GPU and +OMP_NUM_THREADS=1. Saved GPU RI binary and dependency patches are unchanged. +All configurations use GPU cuSOLVER, DZP Si_dojo_6au.orb, Gamma sampling, +SCF energy threshold 1e-8 eV and density threshold 1e-9. Forces/stress included. +Si32 is the earlier same-node comparison, job 1832072 on 16v100n48. +One run per point: timings are observations, not statistically established +speedups. Compare timings within each row; nodes differ between Si32 and the +new measurements. CPU resources also grow with MPI/GPU count. + +| Atoms | Max per-card sampled MiB, 1 / 2 / 4 cards | Wall seconds, 1 / 2 / 4 cards | Four-card speedup | Four-card parallel efficiency | +| --- | --- | --- | ---: | ---: | +| 16 | 860 / 708 / 616 | 40.37 / 39.16 / 38.02 | 1.06x | 26.5% | +| 32 | 2756 / 1682 / 1158 | 63.78 / 57.23 / 52.71 | 1.21x | 30.3% | +| 64 | 6568 / 6262 / 3980 | 178.30 / 137.05 / 101.07 | 1.76x | 44.1% | + +Memory is nvidia-smi per-process usage sampled every 100 ms. Mapping snapshots +confirm the Si64 two-rank processes occupy two distinct GPUs, and four-rank +processes four distinct GPUs. Si64 per-process peaks on four cards were +3980, 3426, 3426 and 3962 MiB. These individual peaks need not occur at the same +time; their sum is not a measured simultaneous aggregate peak. This includes +non-RI allocations; an allocator trace separating RI components was not taken. + +Si16 used 8 outer steps in all configurations, with 74 / 74 / 73 total inner +iterations. Si64 used 7 outer steps, with 63 / 62 / 62 inner iterations. +Maximum energy difference within a size: 2.05e-10 eV total; force differences +zero at printed precision; stress-component differences <= 1.88e-7 kbar. +These symmetric cells have zero total forces; existing regression cases provide +nonzero-force validation. Inner counts differ slightly, so these wall times +are time-to-solution scaling, not perfectly fixed-work kernel benchmarks. + +| Si64 EXX timer | 1 card s | 2 cards s | 4 cards s | 1-to-4 speedup | +| --- | ---: | ---: | ---: | ---: | +| cal_exx_elec | 25.42 | 22.65 | 12.30 | 2.07x | +| cal_exx_force | 33.62 | 24.95 | 17.46 | 1.93x | +| cal_exx_stress | 43.92 | 30.74 | 21.72 | 2.02x | +| init | 11.92 | 12.73 | 13.71 | 0.87x | + +Timers include their enclosing routine's work, not just CUDA kernel duration. +They do not isolate MPI communication, packing or synchronization. Si16's +cal_exx_elec instead rises from 1.74 to 2.53 s with four cards. Small workloads +are unsuitable for assuming efficient strong scaling. + +## What the data establish + +- Multi-GPU distribution relieves per-card capacity pressure, but not as 1/P + universally. Si64 improves only 4.7% with two cards and 39.4% with four. +- A plateau is compatible with the implementation's per-rank temporary cap: + min(T/P, 5 GiB) stays near 5 GiB while each rank still exceeds the cap. This + is a code-supported explanation, not an allocation-level attribution of + the Si64 peak. Other allocations and load imbalance also matter. +- At fixed card count, memory grows substantially with size: 860 -> 2756 -> + 6568 MiB for one card and 616 -> 1158 -> 3980 MiB for four cards. These finite + points do not establish a universal power law. Supercell shapes change + (2x1x1, 2x2x1, 2x2x2), affecting periodic neighborhoods and screening. +- Multi-GPU speed improves for the larger workload, but four-card efficiency + remains about 44% end-to-end and 52% for electronic EXX. Communication scaling + or a single bottleneck cannot be inferred without profiling. +- Increasing cards proportionally with atoms does not guarantee fixed per-card + memory. Unified memory-budgeted batching remains useful even with MPI. +- No multi-node or weak-scaling benchmark has been completed. No host-RAM peak + has been measured. The results do not establish arbitrary-size scalability. + +## Source-level allocation model + +For one rank, one contraction label, the allocator-level components are: + +M_RI = s (I + U + O + min(T, 5 GiB / s)) + B + C, + +where s is the scalar size (8 bytes real double, 16 bytes complex double), +I is the total element count of unique uploaded inputs, U of unique accumulated +intermediates, O of unique outputs, T of unsegmented temporary products, B of +per-task pointer/dimension arrays, and C of backend context/workspace allocations. +Take the maximum over simultaneously live objects/phases; add other ABACUS GPU +allocations to obtain whole-process memory. This is an accounting model, not a +prediction from atom count alone. Certain labels can have different live objects. + +Input, Mul and Output deduplicate by atom-pair key before allocating full local +data. Tmp allocates min(totalSize, limit/sizeof(T)) scalar elements, but its +pointer array is still allocated for the whole task list. In the examined +three-index contraction, a tmp entry and dimension/pointer entries are added +per surviving triple; two-index inputs and accumulated products are deduplicated. +The hardcoded budget is per temporary data buffer, not for total allocations. + +For fixed species/basis and fixed periodic-image range, the dense unscreened +pair storage can grow quadratically with atom count; the triple task list and +unsegmented temp volume can grow cubically. This is an upper-envelope counting +argument for these contraction stages, not a claim that screened silicon follows +N^3, nor a universal complexity of ABACUS HSE. Fixed-density finite-range +screening can make the surviving neighborhoods bounded and substantially lower +the growth. Auxiliary basis size, image count, cutoff and screening matter. + +The default Parallel_LRI_Equally partitions two index ranges with MPI. The GPU +loop builds its batches from rank-local lists and requests corresponding tensor +maps. Tasks can divide approximately as 1/P when balanced. Shared inputs/outputs +can be replicated across one process-grid direction; whole-process memory need +not follow 1/P. In an ideal balanced dense two-dimensional grid, a matrix slice +shared along a row/column can follow N^2/sqrt(P) rather than N^2/P. This is an +illustration of replication, not an empirically fitted scaling law here. + +CPU-side original maps, packed host input copies, output downloads and task +metadata coexist. Host memory cannot be inferred from nvidia-smi measurements; +no host-memory asymptotic scaling has been measured in this experiment. + +## Source locations and reproducibility + +The analysis refers to LibRI revision +`6ce3c8480ae38a085746dc81cc4dbb59b8074430`, under `include/RI/`: + +- `ri/gpu/GPU_Data_Input.h`: deduplicated input packing and upload. +- `ri/gpu/GPU_Data_Tmp.h`: capped temporary data, full-batch pointer allocation. +- `ri/gpu/GPU_Data_Mul.h`, `GPU_Data_Output.h`: full local accumulated/output data. +- `ri/LRI-cal_loop3_GPU.hpp`: rank-local loop lists, per-label lifetime, uploads, + explicit synchronization and downloads. +- `parallel/Parallel_LRI_Equally.hpp`: two-index task split and tensor requests. +- `distribute/Distribute_Equally.hpp`, `Split_Processes.hpp`: index partitioning. + +Actual tested controls, generated inputs, outputs, memory traces and summarized +comparisons are retained in `/mnt/g/work/agent/abacus/lcao-hse-results/scaling/`. +`control/hse_benchmark.py` is the tested matrix driver (Si16 and Si64, 1/2/4 +ranks). Input generation was checked for cell dimensions, atom counts and unique +positions. No production code changes were required for this experiment. +The runner shipped source HEAD `784610dfadc952bfbefd98c298e250814c0e9684` and +reused the previously built GPU RI executable (runtime C++ unchanged). +Remote archive beneath the experiment project root: +`archives/lcao-hse/lcao-scaling-20261017-1.tar.gz`. From 33d400fd6012262c781697716f1a0e19950843d0 Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Sun, 18 Oct 2026 00:28:03 +0800 Subject: [PATCH 06/24] build(lcao): encapsulate LibRI GPU dependencies in CMake targets --- CMakeLists.txt | 12 ++-- cmake/modules/FindDDLA.cmake | 29 +++++++++ cmake/modules/FindLibRI.cmake | 60 ++++++++++++++++++- source/CMakeLists.txt | 18 +----- tests/lcao_hse_gpu/README.md | 32 +++++++++- tests/lcao_hse_gpu/cmake/CMakeLists.txt | 23 +++++++ .../cmake/test_libri_consumer.cpp | 27 +++++++++ 7 files changed, 176 insertions(+), 25 deletions(-) create mode 100644 cmake/modules/FindDDLA.cmake create mode 100644 tests/lcao_hse_gpu/cmake/CMakeLists.txt create mode 100644 tests/lcao_hse_gpu/cmake/test_libri_consumer.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 92d1c6cbf80..f1168f8bc98 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -72,6 +72,10 @@ option(ENABLE_CNPY "Enable cnpy usage" OFF) cmake_dependent_option(ENABLE_ELPA "Enable ELPA eigensolver" ON "ENABLE_MPI" OFF) cmake_dependent_option(ENABLE_LIBRI "Enable LibRI for hybrid functional" OFF "ENABLE_LCAO;ENABLE_MPI" OFF) +option(ENABLE_LIBRI_GPU "Enable experimental LibRI CUDA contractions through DDLA" OFF) +if(ENABLE_LIBRI_GPU AND (NOT ENABLE_LIBRI OR NOT USE_CUDA)) + message(FATAL_ERROR "ENABLE_LIBRI_GPU requires ENABLE_LIBRI and USE_CUDA") +endif() cmake_dependent_option(ENABLE_PEXSI "Enable PEXSI for LCAO" OFF "ENABLE_LCAO;ENABLE_MPI" OFF) cmake_dependent_option(ENABLE_MLALGO "Enable the machine learning algorithms" OFF "ENABLE_LCAO;ENABLE_MPI" OFF) @@ -606,18 +610,10 @@ endfunction() if(ENABLE_LIBRI) set_if_higher(CMAKE_CXX_STANDARD 14) - option(ENABLE_LIBRI_GPU "Enable experimental LibRI CUDA contractions through DDLA" OFF) if(ENABLE_LIBRI_GPU) - if(NOT USE_CUDA) - message(FATAL_ERROR "ENABLE_LIBRI_GPU requires USE_CUDA") - endif() set_if_higher(CMAKE_CXX_STANDARD 17) - find_path(DDLA_INCLUDE_DIR ddla/ddla_connector.h REQUIRED) - find_library(DDLA_LIBRARY NAMES ddla REQUIRED) endif() find_package(LibRI REQUIRED) - find_package(LibComm REQUIRED) - find_package(cereal REQUIRED CONFIG) endif() if(ENABLE_LIBXC) diff --git a/cmake/modules/FindDDLA.cmake b/cmake/modules/FindDDLA.cmake new file mode 100644 index 00000000000..79024169a4c --- /dev/null +++ b/cmake/modules/FindDDLA.cmake @@ -0,0 +1,29 @@ +# CUDA adapter for the DDLA interface currently consumed by LibRI. +# DDLA_ROOT or CMAKE_PREFIX_PATH may name an installation prefix; the explicit +# DDLA_INCLUDE_DIR and DDLA_LIBRARY overrides remain supported. + +find_path(DDLA_INCLUDE_DIR + NAMES ddla/ddla_connector.h + HINTS ${DDLA_ROOT} + PATH_SUFFIXES include) +find_library(DDLA_LIBRARY + NAMES ddla + HINTS ${DDLA_ROOT} + PATH_SUFFIXES lib lib64) + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(DDLA + REQUIRED_VARS DDLA_INCLUDE_DIR DDLA_LIBRARY) + +if(DDLA_FOUND AND NOT TARGET DDLA::DDLA) + add_library(DDLA::DDLA UNKNOWN IMPORTED) + set_target_properties(DDLA::DDLA PROPERTIES + IMPORTED_LOCATION "${DDLA_LIBRARY}" + INTERFACE_INCLUDE_DIRECTORIES "${DDLA_INCLUDE_DIR}" + INTERFACE_COMPILE_DEFINITIONS DDLA_USE_CUDA + INTERFACE_COMPILE_FEATURES cxx_std_17 + INTERFACE_LINK_LIBRARIES + "MPI::MPI_CXX;CUDA::cudart;CUDA::cublas;CUDA::curand;CUDA::cusolver") +endif() + +mark_as_advanced(DDLA_INCLUDE_DIR DDLA_LIBRARY) diff --git a/cmake/modules/FindLibRI.cmake b/cmake/modules/FindLibRI.cmake index fe41ddaf760..c0793fd6740 100644 --- a/cmake/modules/FindLibRI.cmake +++ b/cmake/modules/FindLibRI.cmake @@ -4,6 +4,7 @@ # # LIBRI_FOUND - True if LibRI is found. # LIBRI_DIR - Where to find LibRI files. +# LibRI::LibRI - Header-only target, including optional DDLA requirements. find_path(LIBRI_DIR include/RI/version.h @@ -29,5 +30,62 @@ endif() include(FindPackageHandleStandardArgs) find_package_handle_standard_args(LibRI DEFAULT_MSG LIBRI_DIR) -# Copy the results to the output variables and target. +if(LibRI_FOUND AND NOT TARGET LibRI::LibRI) + find_package(LibComm REQUIRED) + find_package(cereal REQUIRED CONFIG) + + add_library(LibRI::LibRI INTERFACE IMPORTED) + set_target_properties(LibRI::LibRI PROPERTIES + INTERFACE_INCLUDE_DIRECTORIES "${LIBRI_DIR}/include;${LIBCOMM_DIR}/include" + INTERFACE_COMPILE_FEATURES cxx_std_14 + INTERFACE_LINK_LIBRARIES "MPI::MPI_CXX;cereal::cereal") + if(ENABLE_OPENMP) + target_link_libraries(LibRI::LibRI INTERFACE OpenMP::OpenMP_CXX) + endif() + # Reuse ABACUS's selected math backend. Special platforms without these + # targets retain their existing numerical libraries in the final link. + foreach(_libri_math_target IN ITEMS LAPACK::LAPACK BLAS::BLAS) + if(TARGET ${_libri_math_target}) + target_link_libraries(LibRI::LibRI INTERFACE ${_libri_math_target}) + endif() + endforeach() + + if(ENABLE_LIBRI_GPU) + find_package(DDLA REQUIRED) + target_compile_definitions(LibRI::LibRI INTERFACE __GPU_RI __DDLA_RI) + target_compile_features(LibRI::LibRI INTERFACE cxx_std_17) + target_link_libraries(LibRI::LibRI INTERFACE DDLA::DDLA) + + include(CMakePushCheckState) + include(CheckCXXSourceCompiles) + cmake_push_check_state(RESET) + set(CMAKE_REQUIRED_LIBRARIES LibRI::LibRI) + # Recheck when users change dependency paths or apply compatibility fixes. + unset(LIBRI_GPU_API_WORKS CACHE) + check_cxx_source_compiles([=[ + #include + #include + template void check_gpu() + { + RI::LRI lri; + lri.cal_mode = RI::LRI_Cal_Mode::GPU; + std::map>, RI::Tensor>> result; + lri.cal_loop3({}, result, 1.0); + } + int main() + { + check_gpu(); + check_gpu>(); + } + ]=] LIBRI_GPU_API_WORKS) + cmake_pop_check_state() + if(NOT LIBRI_GPU_API_WORKS) + message(FATAL_ERROR + "LibRI/DDLA GPU compile/link check failed. Supply compatible, patched " + "LibRI and CUDA DDLA installations; see tests/lcao_hse_gpu/README.md " + "and the CMake configure log.") + endif() + endif() +endif() + mark_as_advanced(LIBRI_DIR) diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index c6d92c9fd22..be0a36eef61 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -60,9 +60,6 @@ set(_abacus_feature_definitions $<$:__FLOAT_FFTW> $<$:__MLALGO> $<$:__CNPY> - $<$:__GPU_RI> - $<$:__DDLA_RI> - $<$:DDLA_USE_CUDA> $<$:__EXX> $<$:EXX_DM=3> $<$:EXX_H_COMM=2> @@ -277,7 +274,7 @@ set(_abacus_feature_libs $<$:PEXSI::PEXSI> $<$:cnpy> $<$:ZLIB::ZLIB> - $<$:cereal::cereal> + $<$:LibRI::LibRI> $<$:Libxc::xc> $<$:PLUMED::plumed>) set(_abacus_feature_include_dirs) @@ -310,17 +307,6 @@ if(ENABLE_CNPY) endif() endif() -if(ENABLE_LIBRI_GPU) - list(APPEND _abacus_feature_include_dirs ${DDLA_INCLUDE_DIR}) - list(APPEND _abacus_feature_libs ${DDLA_LIBRARY}) -endif() - -if(ENABLE_LIBRI) - list(APPEND _abacus_feature_include_dirs - ${LIBRI_DIR}/include - ${LIBCOMM_DIR}/include) -endif() - if(DEFINED DeePMD_DIR) if(DeePMDC_FOUND) list(APPEND _abacus_feature_libs DeePMD::deepmd_c) @@ -403,6 +389,8 @@ target_compile_options(abacus_compile_requirements INTERFACE foreach(_abacus_feature_target IN ITEMS nlohmann_json::nlohmann_json + LibRI::LibRI + DDLA::DDLA cereal::cereal ELPA::ELPA Libxc::xc diff --git a/tests/lcao_hse_gpu/README.md b/tests/lcao_hse_gpu/README.md index 859756d3a76..7000dea9755 100644 --- a/tests/lcao_hse_gpu/README.md +++ b/tests/lcao_hse_gpu/README.md @@ -3,7 +3,21 @@ Fork-only experiment; not production-ready. ABACUS base: 8ae37fcd6. Runtime selection uses the existing `device cpu/gpu`. Build with `ENABLE_LIBRI=ON`, `USE_CUDA=ON`, `ENABLE_LIBRI_GPU=ON`, and explicit -`LIBRI_DIR`, `DDLA_INCLUDE_DIR`, `DDLA_LIBRARY`. +`LIBRI_DIR`, and a CUDA DDLA installation selected with `DDLA_ROOT` or +`CMAKE_PREFIX_PATH`. Explicit `DDLA_INCLUDE_DIR` and `DDLA_LIBRARY` overrides +remain supported. + +ABACUS's Find modules provide `LibRI::LibRI` and `DDLA::DDLA`; no CMake changes +in either dependency are needed. The LibRI target carries its headers, C++ +requirements, LibComm/cereal and available MPI/OpenMP/math targets. GPU builds +add DDLA and the backend macros through those targets. ABACUS's existing +compile-requirements collector keeps these header-only definitions consistent +across its translation units, with final linking in `source/CMakeLists.txt`. +CPU builds do not discover or link DDLA. Enabling GPU RI without LibRI/CUDA +is an error. A configure-time compile/link probe instantiates real and complex +GPU contractions to reject incompatible dependency APIs; it does not execute +GPU code or verify numerical correctness. The dependency patches below are +still required. Dependencies used in the SAI experiment: @@ -70,3 +84,19 @@ The automatic download still uses ABACUS's original CPU dependency. Downloading an unpatched newer revision is insufficient to reproduce the working build. The build option defaults to OFF; the HSE approximation is unchanged. See RESULTS.md for measured outcomes and remaining validation limits. + +## Checking the CMake dependency interface + +The small consumer below links only `LibRI::LibRI`, without ABACUS's common +include paths or link flags. It checks headers, backend macros, the CPU default +in GPU-enabled LibRI, and transitive BLAS linkage. Run with GPU support OFF +and ON; the ON variant requires CUDA and the patched dependencies above. + +```bash +cmake -S tests/lcao_hse_gpu/cmake -B /path/to/consumer-build \ + -DLIBRI_DIR=/path/to/libri -DLIBCOMM_DIR=/path/to/libcomm \ + -Dcereal_DIR=/path/to/cereal/lib/cmake/cereal \ + -DENABLE_LIBRI_GPU=ON -DDDLA_ROOT=/path/to/ddla +cmake --build /path/to/consumer-build +ctest --test-dir /path/to/consumer-build --output-on-failure +``` diff --git a/tests/lcao_hse_gpu/cmake/CMakeLists.txt b/tests/lcao_hse_gpu/cmake/CMakeLists.txt new file mode 100644 index 00000000000..32ce121adcb --- /dev/null +++ b/tests/lcao_hse_gpu/cmake/CMakeLists.txt @@ -0,0 +1,23 @@ +cmake_minimum_required(VERSION 3.18) +project(libri_consumer_check LANGUAGES CXX) + +# Exercise the dependency target without ABACUS's global includes or link flags. +list(PREPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_LIST_DIR}/../../../cmake/modules") +option(ENABLE_OPENMP "Enable OpenMP" ON) +option(ENABLE_LIBRI_GPU "Check the LibRI CUDA dependency interface" OFF) +find_package(MPI REQUIRED COMPONENTS CXX) +if(ENABLE_OPENMP) + find_package(OpenMP REQUIRED COMPONENTS CXX) +endif() +find_package(Lapack REQUIRED) +if(ENABLE_LIBRI_GPU) + find_package(CUDAToolkit REQUIRED) +endif() +find_package(LibRI REQUIRED) + +add_executable(libri_consumer test_libri_consumer.cpp) +target_link_libraries(libri_consumer PRIVATE LibRI::LibRI) +target_compile_definitions(libri_consumer PRIVATE + $<$:EXPECT_LIBRI_GPU>) +enable_testing() +add_test(NAME libri_consumer COMMAND libri_consumer) diff --git a/tests/lcao_hse_gpu/cmake/test_libri_consumer.cpp b/tests/lcao_hse_gpu/cmake/test_libri_consumer.cpp new file mode 100644 index 00000000000..6ab8e5a5b09 --- /dev/null +++ b/tests/lcao_hse_gpu/cmake/test_libri_consumer.cpp @@ -0,0 +1,27 @@ +#include + +#if defined(EXPECT_LIBRI_GPU) +#if !defined(__GPU_RI) || !defined(__DDLA_RI) || !defined(DDLA_USE_CUDA) +#error "LibRI GPU usage requirements did not reach the consumer" +#endif +#else +#if defined(__GPU_RI) || defined(__DDLA_RI) || defined(DDLA_USE_CUDA) +#error "CPU LibRI consumer unexpectedly requires a GPU backend" +#endif +#endif + +int main() +{ + RI::LRI lri; +#ifdef EXPECT_LIBRI_GPU + if (lri.cal_mode != RI::LRI_Cal_Mode::CPU) + { + return 1; + } + lri.cal_mode = RI::LRI_Cal_Mode::GPU; +#endif + // The consumer must receive the BLAS dependency through LibRI::LibRI. + const double data[] = {3.0, 4.0}; + const double norm = RI::Blas_Interface::nrm2(2, data); + return norm == 5.0 ? 0 : 1; +} From 8dd69254063302fb70b08b01b08cc284a2d9e37e Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Sun, 18 Oct 2026 00:41:15 +0800 Subject: [PATCH 07/24] build: make MKL adapters usable in dependency probes --- cmake/modules/FindMKL.cmake | 17 ++++++++--------- tests/lcao_hse_gpu/README.md | 4 +++- tests/lcao_hse_gpu/cmake/CMakeLists.txt | 6 +++++- 3 files changed, 16 insertions(+), 11 deletions(-) diff --git a/cmake/modules/FindMKL.cmake b/cmake/modules/FindMKL.cmake index 8d324666bc1..3296bc475a2 100644 --- a/cmake/modules/FindMKL.cmake +++ b/cmake/modules/FindMKL.cmake @@ -241,30 +241,29 @@ if(MKL_FOUND) set(_mkl_base_archives ${MKL_INTERFACE_LIB} ${MKL_THREAD} ${MKL_CORE}) _mkl_link_group(_mkl_base ${_mkl_base_archives}) - add_library(abacus_mkl INTERFACE) - add_library(abacus::mkl ALIAS abacus_mkl) + # Imported interfaces also propagate into dependency compile/link probes. + add_library(abacus::mkl INTERFACE IMPORTED GLOBAL) target_include_directories( - abacus_mkl + abacus::mkl INTERFACE "${MKL_INCLUDE}" "${MKL_FFTW_INCLUDE}") - target_link_libraries(abacus_mkl INTERFACE ${_mkl_base} ${_mkl_runtime}) + target_link_libraries(abacus::mkl INTERFACE ${_mkl_base} ${_mkl_runtime}) if(ENABLE_MPI) - add_library(abacus_mkl_scalapack INTERFACE) - add_library(abacus::mkl_scalapack ALIAS abacus_mkl_scalapack) + add_library(abacus::mkl_scalapack INTERFACE IMPORTED GLOBAL) target_include_directories( - abacus_mkl_scalapack + abacus::mkl_scalapack INTERFACE "${MKL_INCLUDE}" "${MKL_FFTW_INCLUDE}") if(_mkl_all_static AND CMAKE_SYSTEM_NAME STREQUAL "Linux") _mkl_link_group(_mkl_cluster ${MKL_SCALAPACK} ${MKL_BLACS} ${_mkl_base_archives}) target_link_libraries( - abacus_mkl_scalapack INTERFACE ${_mkl_cluster} MPI::MPI_CXX ${_mkl_runtime}) + abacus::mkl_scalapack INTERFACE ${_mkl_cluster} MPI::MPI_CXX ${_mkl_runtime}) else() target_link_libraries( - abacus_mkl_scalapack INTERFACE + abacus::mkl_scalapack INTERFACE ${MKL_SCALAPACK} ${MKL_BLACS} abacus::mkl MPI::MPI_CXX) endif() set(MKL_LIBRARIES abacus::mkl_scalapack) diff --git a/tests/lcao_hse_gpu/README.md b/tests/lcao_hse_gpu/README.md index 7000dea9755..fd0cac4c839 100644 --- a/tests/lcao_hse_gpu/README.md +++ b/tests/lcao_hse_gpu/README.md @@ -2,7 +2,7 @@ Fork-only experiment; not production-ready. ABACUS base: 8ae37fcd6. Runtime selection uses the existing `device cpu/gpu`. Build with -`ENABLE_LIBRI=ON`, `USE_CUDA=ON`, `ENABLE_LIBRI_GPU=ON`, and explicit +`ENABLE_LIBRI=ON`, `USE_CUDA=ON`, `ENABLE_LIBRI_GPU=ON`, an explicit `LIBRI_DIR`, and a CUDA DDLA installation selected with `DDLA_ROOT` or `CMAKE_PREFIX_PATH`. Explicit `DDLA_INCLUDE_DIR` and `DDLA_LIBRARY` overrides remain supported. @@ -91,6 +91,8 @@ The small consumer below links only `LibRI::LibRI`, without ABACUS's common include paths or link flags. It checks headers, backend macros, the CPU default in GPU-enabled LibRI, and transitive BLAS linkage. Run with GPU support OFF and ON; the ON variant requires CUDA and the patched dependencies above. +Pass `MKL_ROOT` and optionally `MKL_LINK=static` or `dynamic` to exercise +ABACUS's MKL adapter instead of the default BLAS/LAPACK discovery. ```bash cmake -S tests/lcao_hse_gpu/cmake -B /path/to/consumer-build \ diff --git a/tests/lcao_hse_gpu/cmake/CMakeLists.txt b/tests/lcao_hse_gpu/cmake/CMakeLists.txt index 32ce121adcb..500d5fde3f0 100644 --- a/tests/lcao_hse_gpu/cmake/CMakeLists.txt +++ b/tests/lcao_hse_gpu/cmake/CMakeLists.txt @@ -9,7 +9,11 @@ find_package(MPI REQUIRED COMPONENTS CXX) if(ENABLE_OPENMP) find_package(OpenMP REQUIRED COMPONENTS CXX) endif() -find_package(Lapack REQUIRED) +if(MKLROOT OR MKL_ROOT OR DEFINED ENV{MKLROOT}) + find_package(MKL REQUIRED) +else() + find_package(Lapack REQUIRED) +endif() if(ENABLE_LIBRI_GPU) find_package(CUDAToolkit REQUIRED) endif() From 9d404068145e5bd506a814843a641bde2dca4c3c Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Sun, 18 Oct 2026 00:52:07 +0800 Subject: [PATCH 08/24] test(lcao): record CMake dependency and HSE regressions --- tests/lcao_hse_gpu/RESULTS.md | 61 +++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/tests/lcao_hse_gpu/RESULTS.md b/tests/lcao_hse_gpu/RESULTS.md index 667793bd1d2..fbc126743bd 100644 --- a/tests/lcao_hse_gpu/RESULTS.md +++ b/tests/lcao_hse_gpu/RESULTS.md @@ -177,3 +177,64 @@ Further size-dependent memory and multi-GPU checks (Si16/Si64, job 1832165) are documented in [SCALING.md](SCALING.md). They reveal a memory plateau between one and two cards for Si64; do not extrapolate the Si32 reduction as a universal 1/P law. + +## ABACUS-side CMake integration — 2026-10-18 + +The dependency adapter now lives in ABACUS's Find modules. `LibRI::LibRI` +and `DDLA::DDLA` carry headers, language requirements, backend definitions +and link dependencies. These are ABACUS-provided targets, not upstream +LibRI package exports. `ENABLE_LIBRI_GPU` remains OFF by default; runtime +selection, numerical kernels, memory allocation and dependency patches are +unchanged by this CMake work. + +Full CPU/GPU builds used source `33d400fd6012262c781697716f1a0e19950843d0`, +GNU 13.3.0, CMake 3.31.6 and CUDA 12.9.86 on SAI. The MKL follow-up used +`8dd69254063302fb70b08b01b08cc284a2d9e37e` and oneMKL 2024.2. It changes +the MKL adapters to named imported interfaces so CMake can export their +transitive requirements into a compile/link probe. Existing target names +and numerical link ordering are retained. No runtime C++ changed. + +| Verification | Result | +| --- | --- | +| Job 1833459: standalone CPU/GPU consumers linking only `LibRI::LibRI` | Both built; both CTest runs passed | +| Job 1833459: real/complex GPU contraction compile/link probe | Passed with patched dependencies; rejected unpatched LibRI as expected | +| Job 1833459: full CPU build, CUDA OFF, DDLA discovery disabled | Built; three HSE cases passed all 10 existing reference checks | +| Job 1833459: full GPU build, GPU RI ON | Built; all dynamic dependencies resolved | +| Jobs 1833493_0/1/2: GPU HSE Gamma, complex multi-k and spin-polarized Gamma | All three passed, with 10 existing reference checks in total | +| Job 1833601: MKL static/dynamic, CPU/GPU standalone consumers | All four configure/build/CTest combinations passed | +| Job 1833601: original CPU LibRI `b0eff7a`, DDLA discovery disabled | Configure, build and CTest passed | +| Job 1833601: previous MKL adapter with the new GPU probe | Reproduced the missing-target error; corrected adapter passed | + +HSE cases use the repository's unchanged references. Gamma and spin-polarized +Gamma check energy, force and stress; the complex multi-k reference checks +energy. Copied CPU inputs use ScaLAPACK where GenELPA was explicitly selected; +copied GPU inputs use cuSOLVER. All runtime tests use one MPI rank and +`OMP_NUM_THREADS=1`. This is build/regression validation, not a timing +comparison or a new test of multi-GPU scaling. + +The CPU executable's `ldd` has no DDLA/CUDA libraries. The RI and base object +compile flags contain the three GPU RI definitions only in the GPU build. +MKL consumer link lines select the requested static/shared MKL libraries; +CPU consumers do not link DDLA. The MKL checks validate dependency interfaces, +not a full MKL ABACUS build. The configure-time probe cannot detect numerical +errors such as the force-accumulation issue; the existing LibRI patches and +numerical regressions remain necessary. + +Three local invalid configurations (GPU RI with LibRI OFF, CUDA OFF or MPI +OFF) fail with the intended diagnostic. The repository command +`python3 -m unittest discover -s .ci/slurm -p 'test_runner.py' -q` passed all +45 tests. Governance and whitespace checks passed. + +Two experiment-control issues are retained in the raw evidence: job 1833432 +stopped before compilation because its DDLA reuse path was wrong; the corrected +run passed. Result collection then rejected two resource-directory symlinks +introduced by the local CPU test control. Only those symlink entries were +omitted from a copy of the archive, which passed the original safe extractor; +calculation logs and numerical outputs were retained. No reference or runner +safety check was relaxed. + +Controls and logs: `/mnt/g/work/agent/abacus/lcao-hse-results/cmake-integration2/` +and `/mnt/g/work/agent/abacus/lcao-hse-results/cmake-mkl/`. +Full-build remote archive: +`archives/lcao-hse/cmake-targets-20261018b-1.tar.gz` under the experiment root. +MKL controls and results are under `checks/cmake-mkl-20261018/` there. From d8d32046f74fa713081f0012302ddc03ded08b83 Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Sun, 18 Oct 2026 07:54:11 +0800 Subject: [PATCH 09/24] build(libri): minimize GPU integration and follow CMake conventions --- CMakeLists.txt | 6 +- cmake/CollectBuildInfoVars.cmake | 2 + cmake/modules/FindLibRI.cmake | 4 +- cmake/modules/FindMKL.cmake | 9 +- docs/advanced/install.md | 25 +++ tests/lcao_hse_gpu/README.md | 104 ------------- tests/lcao_hse_gpu/RESULTS.md | 240 ----------------------------- tests/lcao_hse_gpu/SCALING.md | 130 ---------------- tests/lcao_hse_gpu/benchmark.py | 126 --------------- tests/lcao_hse_gpu/build.sbatch.in | 88 ----------- tests/lcao_hse_gpu/config.ini | 50 ------ 11 files changed, 36 insertions(+), 748 deletions(-) delete mode 100644 tests/lcao_hse_gpu/README.md delete mode 100644 tests/lcao_hse_gpu/RESULTS.md delete mode 100644 tests/lcao_hse_gpu/SCALING.md delete mode 100644 tests/lcao_hse_gpu/benchmark.py delete mode 100644 tests/lcao_hse_gpu/build.sbatch.in delete mode 100644 tests/lcao_hse_gpu/config.ini diff --git a/CMakeLists.txt b/CMakeLists.txt index f1168f8bc98..0646ba0ff15 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -72,10 +72,8 @@ option(ENABLE_CNPY "Enable cnpy usage" OFF) cmake_dependent_option(ENABLE_ELPA "Enable ELPA eigensolver" ON "ENABLE_MPI" OFF) cmake_dependent_option(ENABLE_LIBRI "Enable LibRI for hybrid functional" OFF "ENABLE_LCAO;ENABLE_MPI" OFF) -option(ENABLE_LIBRI_GPU "Enable experimental LibRI CUDA contractions through DDLA" OFF) -if(ENABLE_LIBRI_GPU AND (NOT ENABLE_LIBRI OR NOT USE_CUDA)) - message(FATAL_ERROR "ENABLE_LIBRI_GPU requires ENABLE_LIBRI and USE_CUDA") -endif() +cmake_dependent_option(ENABLE_LIBRI_GPU "Enable experimental LibRI CUDA contractions through DDLA" + OFF "ENABLE_LIBRI;USE_CUDA" OFF) cmake_dependent_option(ENABLE_PEXSI "Enable PEXSI for LCAO" OFF "ENABLE_LCAO;ENABLE_MPI" OFF) cmake_dependent_option(ENABLE_MLALGO "Enable the machine learning algorithms" OFF "ENABLE_LCAO;ENABLE_MPI" OFF) diff --git a/cmake/CollectBuildInfoVars.cmake b/cmake/CollectBuildInfoVars.cmake index 740c50ca4b0..e72b2cebbc7 100644 --- a/cmake/CollectBuildInfoVars.cmake +++ b/cmake/CollectBuildInfoVars.cmake @@ -382,6 +382,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_LIBRI_GPU=${ENABLE_LIBRI_GPU}") 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}") @@ -440,6 +441,7 @@ message(" DeePMD = ${ABACUS_DEEPMD_VERSION}") message(" NEP = ${ABACUS_NEP_VERSION}") message(" Features : Libxc = ${ABACUS_LIBXC_VERSION}") message(" LibRI = ${ABACUS_LIBRI_VERSION}") +message(" LibRI GPU = ${ENABLE_LIBRI_GPU}") message(" DFT-D4 = ${ABACUS_DFTD4_VERSION}") message(" JSON = ${ABACUS_JSON_VERSION}") message(" PEXSI = ${ABACUS_PEXSI_VERSION}") diff --git a/cmake/modules/FindLibRI.cmake b/cmake/modules/FindLibRI.cmake index c0793fd6740..0de205d9d31 100644 --- a/cmake/modules/FindLibRI.cmake +++ b/cmake/modules/FindLibRI.cmake @@ -15,7 +15,7 @@ find_path(LIBRI_DIR if(NOT LIBRI_DIR) if(ENABLE_LIBRI_GPU) - message(FATAL_ERROR "ENABLE_LIBRI_GPU requires an explicitly supplied, patched LIBRI_DIR; see tests/lcao_hse_gpu/README.md") + message(FATAL_ERROR "ENABLE_LIBRI_GPU requires an explicitly supplied, patched LIBRI_DIR; see docs/advanced/install.md") endif() include(FetchContent) FetchContent_Declare( @@ -82,7 +82,7 @@ if(LibRI_FOUND AND NOT TARGET LibRI::LibRI) if(NOT LIBRI_GPU_API_WORKS) message(FATAL_ERROR "LibRI/DDLA GPU compile/link check failed. Supply compatible, patched " - "LibRI and CUDA DDLA installations; see tests/lcao_hse_gpu/README.md " + "LibRI and CUDA DDLA installations; see docs/advanced/install.md " "and the CMake configure log.") endif() endif() diff --git a/cmake/modules/FindMKL.cmake b/cmake/modules/FindMKL.cmake index 3296bc475a2..d5e4486a7b0 100644 --- a/cmake/modules/FindMKL.cmake +++ b/cmake/modules/FindMKL.cmake @@ -251,19 +251,20 @@ if(MKL_FOUND) target_link_libraries(abacus::mkl INTERFACE ${_mkl_base} ${_mkl_runtime}) if(ENABLE_MPI) - add_library(abacus::mkl_scalapack INTERFACE IMPORTED GLOBAL) + add_library(abacus_mkl_scalapack INTERFACE) + add_library(abacus::mkl_scalapack ALIAS abacus_mkl_scalapack) target_include_directories( - abacus::mkl_scalapack + abacus_mkl_scalapack INTERFACE "${MKL_INCLUDE}" "${MKL_FFTW_INCLUDE}") if(_mkl_all_static AND CMAKE_SYSTEM_NAME STREQUAL "Linux") _mkl_link_group(_mkl_cluster ${MKL_SCALAPACK} ${MKL_BLACS} ${_mkl_base_archives}) target_link_libraries( - abacus::mkl_scalapack INTERFACE ${_mkl_cluster} MPI::MPI_CXX ${_mkl_runtime}) + abacus_mkl_scalapack INTERFACE ${_mkl_cluster} MPI::MPI_CXX ${_mkl_runtime}) else() target_link_libraries( - abacus::mkl_scalapack INTERFACE + abacus_mkl_scalapack INTERFACE ${MKL_SCALAPACK} ${MKL_BLACS} abacus::mkl MPI::MPI_CXX) endif() set(MKL_LIBRARIES abacus::mkl_scalapack) diff --git a/docs/advanced/install.md b/docs/advanced/install.md index ab1939181be..726b15bc5e5 100644 --- a/docs/advanced/install.md +++ b/docs/advanced/install.md @@ -81,6 +81,31 @@ The new EXX implementation directly depends on [LibRI](https://github.com/abacus `ENABLE_LIBCOMM` is deprecated because LibComm is not a standalone ABACUS feature. CMake locates it automatically as a dependency of LibRI. If you prefer using manually downloaded libraries, enable LibRI and provide their locations via `-DLIBRI_DIR=/path/to/LibRI` and `-DLIBCOMM_DIR=/path/to/LibComm`. ``` +Experimental LCAO GPU contractions additionally require `ENABLE_LIBRI_GPU=ON` +(default OFF), `USE_CUDA=ON`, and a CUDA build of DDLA. Like other dependent +options, GPU RI is disabled when its LibRI/CUDA prerequisites are disabled. +With GPU RI compiled in, INPUT `device gpu` selects GPU contractions and +`device cpu` retains CPU contractions. Otherwise, RI contractions remain on CPU. + +The tested combination is LibRI `6ce3c8480ae38a085746dc81cc4dbb59b8074430` +and DDLA `9c404ac65148ec1b48977b846f38ba6ef9e134f7`. LibRI requires the +compatibility and force-accumulation fixes in `tests/lcao_hse_gpu/libri-compat.patch`: +normalize `include/RI/global/gpu/GPU_Backend.h` and +`include/RI/ri/LRI-cal_loop3_GPU.hpp` to LF, then run +`git apply --unidiff-zero /path/to/abacus/tests/lcao_hse_gpu/libri-compat.patch` +from the LibRI source directory. GPU mode requires this explicitly supplied +LibRI tree; the default CPU dependency download is unchanged. + +```bash +cmake -B build -DENABLE_LIBRI=ON -DUSE_CUDA=ON -DENABLE_LIBRI_GPU=ON \ + -DLIBRI_DIR=/path/to/patched/LibRI -DDDLA_ROOT=/path/to/DDLA +``` + +`CMAKE_PREFIX_PATH` can also locate DDLA. A configure-time compile/link check +rejects incompatible GPU APIs, but does not validate numerical correctness. +This experimental backend has additional GPU memory requirements; its temporary +buffer limit is not a bound on total GPU memory. + ## Build with DFT-D4 support ABACUS can use the external [DFT-D4](https://github.com/dftd4/dftd4) library for Grimme's DFT-D4 dispersion correction. DFT-D4 support is optional and disabled by default. diff --git a/tests/lcao_hse_gpu/README.md b/tests/lcao_hse_gpu/README.md deleted file mode 100644 index fd0cac4c839..00000000000 --- a/tests/lcao_hse_gpu/README.md +++ /dev/null @@ -1,104 +0,0 @@ -# Experimental LCAO HSE GPU integration - -Fork-only experiment; not production-ready. ABACUS base: 8ae37fcd6. -Runtime selection uses the existing `device cpu/gpu`. Build with -`ENABLE_LIBRI=ON`, `USE_CUDA=ON`, `ENABLE_LIBRI_GPU=ON`, an explicit -`LIBRI_DIR`, and a CUDA DDLA installation selected with `DDLA_ROOT` or -`CMAKE_PREFIX_PATH`. Explicit `DDLA_INCLUDE_DIR` and `DDLA_LIBRARY` overrides -remain supported. - -ABACUS's Find modules provide `LibRI::LibRI` and `DDLA::DDLA`; no CMake changes -in either dependency are needed. The LibRI target carries its headers, C++ -requirements, LibComm/cereal and available MPI/OpenMP/math targets. GPU builds -add DDLA and the backend macros through those targets. ABACUS's existing -compile-requirements collector keeps these header-only definitions consistent -across its translation units, with final linking in `source/CMakeLists.txt`. -CPU builds do not discover or link DDLA. Enabling GPU RI without LibRI/CUDA -is an error. A configure-time compile/link probe instantiates real and complex -GPU contractions to reject incompatible dependency APIs; it does not execute -GPU code or verify numerical correctness. The dependency patches below are -still required. - -Dependencies used in the SAI experiment: - -- LibRI 6ce3c8480ae38a085746dc81cc4dbb59b8074430; archive SHA256 - 09c8f766687e5c2236bfe580f372ea5ecef38e935a56c6a24a786d0430ff8fc4. -- LibDDLA 9c404ac65148ec1b48977b846f38ba6ef9e134f7; archive SHA256 - bc349e64c0a308aa9cd0f766e9a766e942bbc974f557522604bd4af9e1ce9baa. -- LibComm and cereal use the existing ABACUS toolchain pins. - -Apply `libri-compat.patch` to LibRI after normalizing its two affected files -to LF, then use `git apply --unidiff-zero /path/to/libri-compat.patch` -from that dependency root. The patch omits context to preserve upstream tab -indentation without introducing whitespace diagnostics in this repository. It repairs the DDLA diagnostic helper name, the split cleanup API, -and GPU contraction accumulation that otherwise ignores `fac_add_Ds`. -The original GPU force regression fails while energy passes; patched force -and stress must be verified against the unchanged references. -The build template records equivalent edits used in the actual experiment. - -Use the existing `.ci/slurm/runner.py` transport, build, worker and artifact -machinery. Experiment controls and full raw results are retained locally at -`/mnt/g/work/agent/abacus/lcao-hse-results/attempt*`. -The copied control runner supports `autotest_gpu` and `hse_benchmark`, as in -the earlier fork experiment; these extensions are not in the base runner. -Set OMP_NUM_THREADS=1 after loading cluster modules. - -Limitations: LibRI's 5 GiB temporary limit is not a total GPU memory budget. -Multi-node GPU mapping and all GPU contraction labels are not validated. -Do not infer production MD readiness from SCF smoke tests. - -## Benchmark interpretation - -The initial Si8/Si32 benchmark compares CPU ScaLAPACK with GPU cuSOLVER; -wall time therefore includes a solver change. Use the additional GPU-device -build with `ENABLE_LIBRI_GPU=OFF` to isolate RI offload. -`exx_hybrid_step=50` is an outer-loop cap, not a convergence guarantee. -Inspect the printed outer deltaE before calling a result converged. -GPU process memory is sampled at 100 ms; it is not a strict allocator peak -or a proof of an upper bound for other systems. - -## Running a backend comparison - -On a GPU allocation with the cluster modules loaded, set `OMP_NUM_THREADS=1` -and make DDLA's `lib` directory available in `LD_LIBRARY_PATH`. Build two -executables with the same patched LibRI and identical ABACUS options, changing -only `ENABLE_LIBRI_GPU`. Both runs use `device gpu` and `ks_solver cusolver`. - -```bash -python3 tests/lcao_hse_gpu/benchmark.py \ - --abacus /path/to/gpu-ri/abacus --baseline /path/to/cpu-ri/abacus \ - --source /path/to/abacus-source --output /path/to/new-results \ - --orbital Si_dojo_6au.orb --repeats 2 -``` - -The source must contain the test PP_ORB assets. `ABACUS_CPU_RI` may supply the -baseline path when invoked by the local `.ci/slurm` experiment controls. -The benchmark alternates backend order, checks outer energy convergence and -iteration count, and compares total energy, every force component and every -stress component. Thresholds are 1e-6 eV/atom, 1e-4 eV/Angstrom and 1e-3 kbar. -These are backend parity tests, not basis/cutoff convergence studies; the -small SZ system can break symmetry under prolonged outer iteration. - -GPU mode deliberately requires an externally supplied, patched LibRI tree. -The automatic download still uses ABACUS's original CPU dependency. Downloading -an unpatched newer revision is insufficient to reproduce the working build. -The build option defaults to OFF; the HSE approximation is unchanged. -See RESULTS.md for measured outcomes and remaining validation limits. - -## Checking the CMake dependency interface - -The small consumer below links only `LibRI::LibRI`, without ABACUS's common -include paths or link flags. It checks headers, backend macros, the CPU default -in GPU-enabled LibRI, and transitive BLAS linkage. Run with GPU support OFF -and ON; the ON variant requires CUDA and the patched dependencies above. -Pass `MKL_ROOT` and optionally `MKL_LINK=static` or `dynamic` to exercise -ABACUS's MKL adapter instead of the default BLAS/LAPACK discovery. - -```bash -cmake -S tests/lcao_hse_gpu/cmake -B /path/to/consumer-build \ - -DLIBRI_DIR=/path/to/libri -DLIBCOMM_DIR=/path/to/libcomm \ - -Dcereal_DIR=/path/to/cereal/lib/cmake/cereal \ - -DENABLE_LIBRI_GPU=ON -DDDLA_ROOT=/path/to/ddla -cmake --build /path/to/consumer-build -ctest --test-dir /path/to/consumer-build --output-on-failure -``` diff --git a/tests/lcao_hse_gpu/RESULTS.md b/tests/lcao_hse_gpu/RESULTS.md deleted file mode 100644 index fbc126743bd..00000000000 --- a/tests/lcao_hse_gpu/RESULTS.md +++ /dev/null @@ -1,240 +0,0 @@ -# SAI LCAO HSE GPU experiment — 2026-10-17 - -Experimental integration only; no upstream PR. Computing resources were -provided by the Open Source Supercomputing Center of SAI (https://www.open-sai.com/). - -ABACUS base: `8ae37fcd6cb9043a4a0c7ebde2c062d97f393bdb`. -Runtime source: `b408a31374fa6134618b60d075cb50d8cf1b1c0c`. -Subsequent changes are test/documentation and a configure-time guard against -fetching an unpatched GPU dependency; runtime C++ is unchanged. -LibRI/DDLA revisions and patches are in README.md and libri-compat.patch. -Hardware: one V100 16 GB, one MPI rank, OMP_NUM_THREADS=1, CUDA 12.9.86, -GNU 13.3.0, OpenMPI. Both comparison executables use GPU cuSOLVER. -The CPU RI baseline uses the same newer LibRI, not the original b0eff7a. - -## Build and regression evidence - -| Job | Result | -| --- | --- | -| 1818882 | Patched GPU RI build passed | -| 1818916_0 / 1818916_2 | Existing HSE Gamma / spin-polarized Gamma energy, force, stress references passed | -| 1819337_0 | Existing complex multi-k HSE references passed with GPU cuSOLVER | -| 1819559 | CPU RI control build passed, CUDA and same newer LibRI retained | -| 1819665 | SZ Si8/Si32 matched-backend comparison passed | -| 1819761 | DZP Si8/Si32, two alternating-order repeats per backend: all 8 calculations and 4 paired comparisons passed | - -Failures retained in raw artifacts: unpatched GPU force accumulation failed -reference tests despite passing energy; incompatible DDLA helper and LibRI -cleanup APIs required fixes. The installed ELPA/GenELPA paths were unsuitable -for these GPU comparisons (including an illegal-instruction exit); cuSOLVER -is used consistently in the final matched comparison. Earlier failed jobs -are not counted as all-pass suites. - -With default `scf_ene_thr=-1`, the Si32 CPU RI and GPU RI trajectories had -matching early EXX energies but different outer stopping times. GPU RI reached -the 50-step cap and the final states differed. Explicit `scf_ene_thr=1e-8` eV -produced equal outer counts and matching results. This supports convergence -sensitivity as the cause of that discrepancy, but does not establish default -parameter robustness or production MD readiness. - -## SZ matched comparison - -| System | CPU RI wall s | GPU RI wall s | CPU/GPU sampled MiB | Outer steps | -| --- | ---: | ---: | ---: | ---: | -| Si8 SZ | 8.91 | 9.19 | 422 / 498 | 4 / 4 | -| Si32 SZ | 24.31 | 24.13 | 438 / 530 | 5 / 5 | - -Single paired measurements: no demonstrated speedup. Total energy differences -were below 4e-11 eV, forces identical at printed precision, maximum stress -component difference 9.98e-7 kbar. Si32 electronic EXX time increased from -0.50 to 0.80 s, force EXX from 0.80 to 1.17 s, stress EXX from 0.96 to 1.57 s. -Small contractions can be slower after offload. - -## DZP matched comparison - -Si_dojo_6au.orb is a 2s2p1d basis. Two repeats per backend, alternating -CPU/GPU order; means below include SCF, force and stress calculation. - -| System | CPU RI wall s | GPU RI wall s | Wall speedup | CPU/GPU sampled MiB | -| --- | ---: | ---: | ---: | ---: | -| Si8 DZP | 30.99 | 29.64 | 1.05x | 442 / 558–564 | -| Si32 DZP | 250.79 | 73.81 | 3.40x | 450 / 2756 | - -Si32 mean electronic EXX contraction time: 50.20 -> 5.36 s (9.36x). -EXX force: 64.95 -> 6.89 s; EXX stress: 85.09 -> 9.68 s. Both backends -used 7 outer iterations for both sizes. Maximum paired energy difference: -1.03e-12 eV/atom; force difference zero at printed precision; stress difference -1.11e-7 kbar. Symmetric silicon forces are zero, so nonzero-force coverage -comes from the existing HSE regression references, not these benchmark cells. - -These speedups are relative to single-thread CPU RI; a tuned multi-core CPU -comparison has not been performed. They are not a universal HSE speedup. -All eight memory logs contain a single GPU process PID. Si32 GPU RI added -2306 MiB (2.25 GiB) over the GPU-solver/CPU-RI baseline: a substantial increase, -not a guarantee of small memory overhead. If the requirement is less than -2x *total* GPU memory relative to that baseline, this implementation fails it. - -## Scope and limitations - -Memory is sampled every 100 ms and includes the GPU eigensolver. It is not -a strict allocation peak. LibRI's hardcoded 5 GiB temporary batching limit -excludes other simultaneous allocations, and is not a total memory bound. -No guarantee that large-system memory stays below twice a baseline is made. -Multi-GPU/multi-node, long MD, sustained allocation behavior, SOC and larger -production cells remain unvalidated. Do not enable this path by default. - -Local raw results: `/mnt/g/work/agent/abacus/lcao-hse-results/`. -The corresponding remote Slurm archives are under -`~/abacus_hse_gpu_experiment/archives/lcao-hse/`. -The experiment reused `.ci/slurm` transport/worker/artifact handling with local -controls; benchmark and dependency fixes are preserved here without importing -unrelated runner changes. - -## Follow-up priorities suggested by the evidence - -1. Keep the small SZ contraction workload on CPU unless profiling shows a - crossover. GPU dispatch/packing/transfer overhead can dominate these blocks. -2. Profile and optimize `refresh_coulomb`/initialization separately. Offloading - contraction GEMMs cannot accelerate this measured CPU-dominated phase. -3. LibRI currently constructs a GPU context per `cal_loop3` call and creates, - packs and uploads input/temporary/output objects per contraction label. - Reuse a single explicitly capped allocation arena and queue where lifetimes - allow it; do not keep two full workspaces or unbounded per-shape caches. -4. Investigate batching by shape and selective reuse of invariant tensors, - measuring transfer/packing/kernel time separately. Persistent tensors must - share a total memory budget with temporary and output buffers. - -These are proposed next experiments, not demonstrated speedups. The current -change connects an existing double-precision GPU contraction implementation; -it does not introduce mixed precision, new screening thresholds, or a different -HSE approximation. Floating-point accumulation order can differ (DDLA uses -atomic accumulation in relevant batched kernels), so bitwise identity is not -expected. Stronger screening or mixed precision would need separate accuracy -validation and are outside this experiment. - -## CP2K comparison - -The relevant transferable ideas are block batching and explicit memory-aware -contraction scheduling, not replacing an ABACUS matrix call with DBCSR directly. -CP2K's [RI HFX implementation](https://github.com/cp2k/cp2k/blob/0adc3f08c83fcdafa7d545917ddb4e901470c070/src/hfx_ri.F) -uses DBCSR operations and explicitly batches occupied-MO/AO indices to save -memory. DBCSR's [libsmm_acc](https://github.com/cp2k/dbcsr/blob/8f63c8f49ce807b5127cf8023b4684f579d19225/src/acc/libsmm_acc/README.md) -executes stacks of small matrix products and tunes kernels for matrix dimensions -and GPU models. LibRI already has variable-size batched contractions, so merely -adding batching is not a new optimization; block-shape distribution, packing, -transfer and kernel efficiency should determine the next change. Different -RI formulations/data structures preclude assuming CP2K's speedup transfers -unchanged to ABACUS. - -## Local checks - -- Repository `.ci/slurm/test_runner.py`: 45 tests passed. Running that suite - from the copied experiment controls initially failed because it assumes the - original 117-case configuration and repository layout; the repository suite - was rerun in its intended location without modifying expectations. -- `agent_governance_check.py --staged`: no findings. -- Python benchmark syntax/CLI and shell template syntax checks passed. -- LibRI patch apply check passed after LF normalization (`--unidiff-zero`). -- Missing GPU LibRI dependency produces the expected CMake diagnostic; this - guard was checked remotely with CMake 3.31.6. -- Final staged whitespace check passed. - -## Single-node multi-GPU memory scaling - -Follow-up Slurm job 1832072 passed on node 16v100n48. The same saved GPU RI -binary and Si32 DZP input ran sequentially with 1, 2 and 4 MPI ranks, one rank -per GPU, OMP_NUM_THREADS=1, GPU cuSOLVER. All runs used 7 outer iterations. -One measurement per configuration, on the same four-V100 allocation: - -| Active GPUs / ranks | Per-process sampled peak MiB | Wall seconds | -| --- | --- | ---: | -| 1 | 2756 | 63.78 | -| 2 | 1682, 1682 | 57.23 | -| 4 | 1158, 1106, 1106, 1158 | 52.71 | - -A live nvidia-smi UUID/PID snapshot confirmed the four listed processes were -on four distinct GPUs; a separate live query confirmed distinct GPUs for the -two-rank run. Peaks are sampled at 100 ms, not allocator-level bounds. -Maximum total-energy difference from one rank: 5.10e-11 eV. Forces agree at -printed precision; maximum stress-component difference: 1.30e-7 kbar. - -Thus distributing RI work *does* lower per-card memory in this tested case: -39% reduction with two cards and 58% with four. Total replicated memory grows, -and the slowest card does not follow ideal 1/N scaling. Performance improved -only 1.21x from one to four cards on this small system. These measurements -support using more GPUs to relieve capacity pressure, but do not establish -large-system or multi-node scaling or guarantee arbitrary capacity growth. -The CPU-RI baseline's small GPU footprint is not an appropriate standalone -rejection criterion for GPU offload; the relevant capacity metric is maximum -per-card memory as system size and GPU count increase. - -Controls, full input/output, per-PID memory traces and UUID mapping snapshot: -`/mnt/g/work/agent/abacus/lcao-hse-results/multigpu/`. -Remote archive: `archives/lcao-hse/lcao-multigpu-20261017-1.tar.gz` beneath the -experiment project root. No runtime code was changed for this test. - -Further size-dependent memory and multi-GPU checks (Si16/Si64, job 1832165) -are documented in [SCALING.md](SCALING.md). They reveal a memory plateau -between one and two cards for Si64; do not extrapolate the Si32 reduction -as a universal 1/P law. - -## ABACUS-side CMake integration — 2026-10-18 - -The dependency adapter now lives in ABACUS's Find modules. `LibRI::LibRI` -and `DDLA::DDLA` carry headers, language requirements, backend definitions -and link dependencies. These are ABACUS-provided targets, not upstream -LibRI package exports. `ENABLE_LIBRI_GPU` remains OFF by default; runtime -selection, numerical kernels, memory allocation and dependency patches are -unchanged by this CMake work. - -Full CPU/GPU builds used source `33d400fd6012262c781697716f1a0e19950843d0`, -GNU 13.3.0, CMake 3.31.6 and CUDA 12.9.86 on SAI. The MKL follow-up used -`8dd69254063302fb70b08b01b08cc284a2d9e37e` and oneMKL 2024.2. It changes -the MKL adapters to named imported interfaces so CMake can export their -transitive requirements into a compile/link probe. Existing target names -and numerical link ordering are retained. No runtime C++ changed. - -| Verification | Result | -| --- | --- | -| Job 1833459: standalone CPU/GPU consumers linking only `LibRI::LibRI` | Both built; both CTest runs passed | -| Job 1833459: real/complex GPU contraction compile/link probe | Passed with patched dependencies; rejected unpatched LibRI as expected | -| Job 1833459: full CPU build, CUDA OFF, DDLA discovery disabled | Built; three HSE cases passed all 10 existing reference checks | -| Job 1833459: full GPU build, GPU RI ON | Built; all dynamic dependencies resolved | -| Jobs 1833493_0/1/2: GPU HSE Gamma, complex multi-k and spin-polarized Gamma | All three passed, with 10 existing reference checks in total | -| Job 1833601: MKL static/dynamic, CPU/GPU standalone consumers | All four configure/build/CTest combinations passed | -| Job 1833601: original CPU LibRI `b0eff7a`, DDLA discovery disabled | Configure, build and CTest passed | -| Job 1833601: previous MKL adapter with the new GPU probe | Reproduced the missing-target error; corrected adapter passed | - -HSE cases use the repository's unchanged references. Gamma and spin-polarized -Gamma check energy, force and stress; the complex multi-k reference checks -energy. Copied CPU inputs use ScaLAPACK where GenELPA was explicitly selected; -copied GPU inputs use cuSOLVER. All runtime tests use one MPI rank and -`OMP_NUM_THREADS=1`. This is build/regression validation, not a timing -comparison or a new test of multi-GPU scaling. - -The CPU executable's `ldd` has no DDLA/CUDA libraries. The RI and base object -compile flags contain the three GPU RI definitions only in the GPU build. -MKL consumer link lines select the requested static/shared MKL libraries; -CPU consumers do not link DDLA. The MKL checks validate dependency interfaces, -not a full MKL ABACUS build. The configure-time probe cannot detect numerical -errors such as the force-accumulation issue; the existing LibRI patches and -numerical regressions remain necessary. - -Three local invalid configurations (GPU RI with LibRI OFF, CUDA OFF or MPI -OFF) fail with the intended diagnostic. The repository command -`python3 -m unittest discover -s .ci/slurm -p 'test_runner.py' -q` passed all -45 tests. Governance and whitespace checks passed. - -Two experiment-control issues are retained in the raw evidence: job 1833432 -stopped before compilation because its DDLA reuse path was wrong; the corrected -run passed. Result collection then rejected two resource-directory symlinks -introduced by the local CPU test control. Only those symlink entries were -omitted from a copy of the archive, which passed the original safe extractor; -calculation logs and numerical outputs were retained. No reference or runner -safety check was relaxed. - -Controls and logs: `/mnt/g/work/agent/abacus/lcao-hse-results/cmake-integration2/` -and `/mnt/g/work/agent/abacus/lcao-hse-results/cmake-mkl/`. -Full-build remote archive: -`archives/lcao-hse/cmake-targets-20261018b-1.tar.gz` under the experiment root. -MKL controls and results are under `checks/cmake-mkl-20261018/` there. diff --git a/tests/lcao_hse_gpu/SCALING.md b/tests/lcao_hse_gpu/SCALING.md deleted file mode 100644 index 9bed36b355b..00000000000 --- a/tests/lcao_hse_gpu/SCALING.md +++ /dev/null @@ -1,130 +0,0 @@ -# LCAO HSE GPU memory and strong-scaling verification - -## Measurements - -Slurm job 1832165 passed on SAI node 16v100n41 (four V100 16 GB GPUs). -Si16 and Si64 each ran with 1, 2 and 4 MPI ranks, one rank per GPU and -OMP_NUM_THREADS=1. Saved GPU RI binary and dependency patches are unchanged. -All configurations use GPU cuSOLVER, DZP Si_dojo_6au.orb, Gamma sampling, -SCF energy threshold 1e-8 eV and density threshold 1e-9. Forces/stress included. -Si32 is the earlier same-node comparison, job 1832072 on 16v100n48. -One run per point: timings are observations, not statistically established -speedups. Compare timings within each row; nodes differ between Si32 and the -new measurements. CPU resources also grow with MPI/GPU count. - -| Atoms | Max per-card sampled MiB, 1 / 2 / 4 cards | Wall seconds, 1 / 2 / 4 cards | Four-card speedup | Four-card parallel efficiency | -| --- | --- | --- | ---: | ---: | -| 16 | 860 / 708 / 616 | 40.37 / 39.16 / 38.02 | 1.06x | 26.5% | -| 32 | 2756 / 1682 / 1158 | 63.78 / 57.23 / 52.71 | 1.21x | 30.3% | -| 64 | 6568 / 6262 / 3980 | 178.30 / 137.05 / 101.07 | 1.76x | 44.1% | - -Memory is nvidia-smi per-process usage sampled every 100 ms. Mapping snapshots -confirm the Si64 two-rank processes occupy two distinct GPUs, and four-rank -processes four distinct GPUs. Si64 per-process peaks on four cards were -3980, 3426, 3426 and 3962 MiB. These individual peaks need not occur at the same -time; their sum is not a measured simultaneous aggregate peak. This includes -non-RI allocations; an allocator trace separating RI components was not taken. - -Si16 used 8 outer steps in all configurations, with 74 / 74 / 73 total inner -iterations. Si64 used 7 outer steps, with 63 / 62 / 62 inner iterations. -Maximum energy difference within a size: 2.05e-10 eV total; force differences -zero at printed precision; stress-component differences <= 1.88e-7 kbar. -These symmetric cells have zero total forces; existing regression cases provide -nonzero-force validation. Inner counts differ slightly, so these wall times -are time-to-solution scaling, not perfectly fixed-work kernel benchmarks. - -| Si64 EXX timer | 1 card s | 2 cards s | 4 cards s | 1-to-4 speedup | -| --- | ---: | ---: | ---: | ---: | -| cal_exx_elec | 25.42 | 22.65 | 12.30 | 2.07x | -| cal_exx_force | 33.62 | 24.95 | 17.46 | 1.93x | -| cal_exx_stress | 43.92 | 30.74 | 21.72 | 2.02x | -| init | 11.92 | 12.73 | 13.71 | 0.87x | - -Timers include their enclosing routine's work, not just CUDA kernel duration. -They do not isolate MPI communication, packing or synchronization. Si16's -cal_exx_elec instead rises from 1.74 to 2.53 s with four cards. Small workloads -are unsuitable for assuming efficient strong scaling. - -## What the data establish - -- Multi-GPU distribution relieves per-card capacity pressure, but not as 1/P - universally. Si64 improves only 4.7% with two cards and 39.4% with four. -- A plateau is compatible with the implementation's per-rank temporary cap: - min(T/P, 5 GiB) stays near 5 GiB while each rank still exceeds the cap. This - is a code-supported explanation, not an allocation-level attribution of - the Si64 peak. Other allocations and load imbalance also matter. -- At fixed card count, memory grows substantially with size: 860 -> 2756 -> - 6568 MiB for one card and 616 -> 1158 -> 3980 MiB for four cards. These finite - points do not establish a universal power law. Supercell shapes change - (2x1x1, 2x2x1, 2x2x2), affecting periodic neighborhoods and screening. -- Multi-GPU speed improves for the larger workload, but four-card efficiency - remains about 44% end-to-end and 52% for electronic EXX. Communication scaling - or a single bottleneck cannot be inferred without profiling. -- Increasing cards proportionally with atoms does not guarantee fixed per-card - memory. Unified memory-budgeted batching remains useful even with MPI. -- No multi-node or weak-scaling benchmark has been completed. No host-RAM peak - has been measured. The results do not establish arbitrary-size scalability. - -## Source-level allocation model - -For one rank, one contraction label, the allocator-level components are: - -M_RI = s (I + U + O + min(T, 5 GiB / s)) + B + C, - -where s is the scalar size (8 bytes real double, 16 bytes complex double), -I is the total element count of unique uploaded inputs, U of unique accumulated -intermediates, O of unique outputs, T of unsegmented temporary products, B of -per-task pointer/dimension arrays, and C of backend context/workspace allocations. -Take the maximum over simultaneously live objects/phases; add other ABACUS GPU -allocations to obtain whole-process memory. This is an accounting model, not a -prediction from atom count alone. Certain labels can have different live objects. - -Input, Mul and Output deduplicate by atom-pair key before allocating full local -data. Tmp allocates min(totalSize, limit/sizeof(T)) scalar elements, but its -pointer array is still allocated for the whole task list. In the examined -three-index contraction, a tmp entry and dimension/pointer entries are added -per surviving triple; two-index inputs and accumulated products are deduplicated. -The hardcoded budget is per temporary data buffer, not for total allocations. - -For fixed species/basis and fixed periodic-image range, the dense unscreened -pair storage can grow quadratically with atom count; the triple task list and -unsegmented temp volume can grow cubically. This is an upper-envelope counting -argument for these contraction stages, not a claim that screened silicon follows -N^3, nor a universal complexity of ABACUS HSE. Fixed-density finite-range -screening can make the surviving neighborhoods bounded and substantially lower -the growth. Auxiliary basis size, image count, cutoff and screening matter. - -The default Parallel_LRI_Equally partitions two index ranges with MPI. The GPU -loop builds its batches from rank-local lists and requests corresponding tensor -maps. Tasks can divide approximately as 1/P when balanced. Shared inputs/outputs -can be replicated across one process-grid direction; whole-process memory need -not follow 1/P. In an ideal balanced dense two-dimensional grid, a matrix slice -shared along a row/column can follow N^2/sqrt(P) rather than N^2/P. This is an -illustration of replication, not an empirically fitted scaling law here. - -CPU-side original maps, packed host input copies, output downloads and task -metadata coexist. Host memory cannot be inferred from nvidia-smi measurements; -no host-memory asymptotic scaling has been measured in this experiment. - -## Source locations and reproducibility - -The analysis refers to LibRI revision -`6ce3c8480ae38a085746dc81cc4dbb59b8074430`, under `include/RI/`: - -- `ri/gpu/GPU_Data_Input.h`: deduplicated input packing and upload. -- `ri/gpu/GPU_Data_Tmp.h`: capped temporary data, full-batch pointer allocation. -- `ri/gpu/GPU_Data_Mul.h`, `GPU_Data_Output.h`: full local accumulated/output data. -- `ri/LRI-cal_loop3_GPU.hpp`: rank-local loop lists, per-label lifetime, uploads, - explicit synchronization and downloads. -- `parallel/Parallel_LRI_Equally.hpp`: two-index task split and tensor requests. -- `distribute/Distribute_Equally.hpp`, `Split_Processes.hpp`: index partitioning. - -Actual tested controls, generated inputs, outputs, memory traces and summarized -comparisons are retained in `/mnt/g/work/agent/abacus/lcao-hse-results/scaling/`. -`control/hse_benchmark.py` is the tested matrix driver (Si16 and Si64, 1/2/4 -ranks). Input generation was checked for cell dimensions, atom counts and unique -positions. No production code changes were required for this experiment. -The runner shipped source HEAD `784610dfadc952bfbefd98c298e250814c0e9684` and -reused the previously built GPU RI executable (runtime C++ unchanged). -Remote archive beneath the experiment project root: -`archives/lcao-hse/lcao-scaling-20261017-1.tar.gz`. diff --git a/tests/lcao_hse_gpu/benchmark.py b/tests/lcao_hse_gpu/benchmark.py deleted file mode 100644 index 4527e527406..00000000000 --- a/tests/lcao_hse_gpu/benchmark.py +++ /dev/null @@ -1,126 +0,0 @@ -#!/usr/bin/env python3 -"""Compare two RI backends using identical GPU LCAO HSE inputs.""" -import argparse -import json -import os -import re -import subprocess -import time -from pathlib import Path - - -def derivatives(log, size): - """Read the final total force (eV/Angstrom) and stress (kbar).""" - force_block = log.rsplit('#TOTAL-FORCE (eV/Angstrom)#', 1)[1] - forces = re.findall(r'^\s*Si\d+\s+([-+\d.eE]+)\s+([-+\d.eE]+)\s+([-+\d.eE]+)', force_block, re.M) - stress_block = log.rsplit('#TOTAL-STRESS (kbar)#', 1)[1] - stress = re.findall(r'^\s*([-+\d.eE]+)\s+([-+\d.eE]+)\s+([-+\d.eE]+)\s*$', stress_block, re.M)[:3] - if len(forces) != size or len(stress) != 3: - raise RuntimeError('Missing force/stress output') - return {key: [float(x) for row in values for x in row] - for key, values in [('force', forces), ('stress', stress)]} - - -def write_inputs(work, source, size, orbital): - points = [(0, 0, 0), (0, .5, .5), (.5, 0, .5), (.5, .5, 0), - (.25, .25, .25), (.25, .75, .75), (.75, .25, .75), (.75, .75, .25)] - nx = 2 if size == 32 else 1 - positions = [((x+i)/nx, (y+j)/nx, z) for i in range(nx) - for j in range(nx) for x, y, z in points] - work.mkdir() - (work / 'STRU').write_text( - 'ATOMIC_SPECIES\nSi 28.0855 Si_dojo_nsoc.upf\nNUMERICAL_ORBITAL\n' - f'{orbital}\nLATTICE_CONSTANT\n10.26\nLATTICE_VECTORS\n' - f'{nx} 0 0\n0 {nx} 0\n0 0 1\nATOMIC_POSITIONS\nDirect\nSi\n0\n{size}\n' - + ''.join('%g %g %g 0 0 0\n' % p for p in positions)) - (work / 'KPT').write_text('K_POINTS\n0\nGamma\n1 1 1 0 0 0\n') - values = dict( - suffix='bench', calculation='scf', basis_type='lcao', device='gpu', - ks_solver='cusolver', gamma_only=1, nspin=1, ntype=1, nbands=int(size*2.5), - pseudo_dir=str(source/'tests/PP_ORB'), orbital_dir=str(source/'tests/PP_ORB'), - ecutwfc=30, scf_thr='1e-9', scf_ene_thr='1e-8', scf_nmax=200, symmetry=0, - smearing_method='gaussian', smearing_sigma=.002, mixing_beta=.1, - dft_functional='hse', exx_real_number=1, exx_pca_threshold='1e-3', - exx_ccp_rmesh_times=1, exx_separate_loop=1, exx_hybrid_step=50, - cal_force=1, cal_stress=1) - (work/'INPUT').write_text('INPUT_PARAMETERS\n' + ''.join(f'{k} {v}\n' for k, v in values.items())) - - -def run_case(binary, work, size): - with (work/'gpu-memory.csv').open('w') as memory: - monitor = subprocess.Popen( - ['nvidia-smi', '--query-compute-apps=pid,used_gpu_memory', - '--format=csv,noheader,nounits', '-lms', '100'], - stdout=memory, stderr=subprocess.STDOUT) - start = time.monotonic() - try: - with (work/'stdout.log').open('w') as stdout: - result = subprocess.run( - ['timeout', '--kill-after=30s', '1200s', 'mpirun', '-np', '1', str(binary)], - cwd=work, stdout=stdout, stderr=subprocess.STDOUT, check=False) - finally: - monitor.terminate() - monitor.wait() - seconds = time.monotonic() - start - if result.returncode: - raise RuntimeError(f'{work}: exit code {result.returncode}') - log = (work/'OUT.bench/running_scf.log').read_text() - stdout = (work/'stdout.log').read_text() - energy = re.findall(r'!FINAL_ETOT_IS\s+([-+\d.eE]+)', log+'\n'+stdout) - deltas = [float(x) for x in re.findall(r'deltaE \(eV\) from outer loop:\s*([-+\d.eE]+)', stdout)] - # An inner SCF marker alone does not establish outer EXX convergence. - if not energy or '#SCF IS CONVERGED#' not in log or not deltas or deltas[-1] >= 1e-8 or len(deltas) >= 50: - raise RuntimeError(f'{work}: outer energy convergence not established') - timing = {cls+'.'+fn: float(value) for cls, fn, value in re.findall( - r'^\s*(Exx_LRI|RI_2D_Comm|HSolverLCAO)\s+(\w+)\s+([\d.]+)', stdout, re.M)} - memory = [int(x) for x in re.findall(r'^\s*\d+,\s*(\d+)\s*$', (work/'gpu-memory.csv').read_text(), re.M)] - return dict(seconds=seconds, energy_ev=float(energy[-1]), outer_steps=len(deltas), - outer_delta_ev=deltas[-1], timing=timing, - sampled_max_process_mib=max(memory) if memory else None, - **derivatives(log, size)) - - -def main(): - parser = argparse.ArgumentParser(description=__doc__) - parser.add_argument('--abacus', type=Path, required=True, help='GPU RI executable') - parser.add_argument('--baseline', type=Path, default=os.environ.get('ABACUS_CPU_RI'), - help='CPU RI executable, built with the same LibRI and GPU solver') - parser.add_argument('--source', type=Path, required=True) - parser.add_argument('--output', type=Path, required=True) - parser.add_argument('--orbital', default='Si_dojo_6au_sz.orb') - parser.add_argument('--repeats', type=int, default=2) - args = parser.parse_args() - if args.baseline is None or args.repeats < 1: - parser.error('A baseline executable and positive repeat count are required') - args.output.mkdir(parents=True, exist_ok=True) - rows = [] - comparisons = [] - for size in (8, 32): - for repeat in range(args.repeats): - pair = {} - order = ('cpu_ri', 'gpu_ri') if repeat % 2 == 0 else ('gpu_ri', 'cpu_ri') - for backend in order: - work = args.output/f'si{size}_{backend}_{repeat}' - write_inputs(work, args.source, size, args.orbital) - binary = args.baseline if backend == 'cpu_ri' else args.abacus - row = dict(size=size, repeat=repeat, backend=backend, - orbital=args.orbital, **run_case(binary, work, size)) - rows.append(row) - pair[backend] = row - (args.output/'measurements.json').write_text(json.dumps(rows, indent=2)+'\n') - print(json.dumps({k: v for k, v in row.items() if k not in ('force', 'stress')}), flush=True) - cpu, gpu = pair['cpu_ri'], pair['gpu_ri'] - comparison = dict(size=size, repeat=repeat, - energy_ev_per_atom=abs(cpu['energy_ev']-gpu['energy_ev'])/size, - force_ev_angstrom=max(abs(x-y) for x, y in zip(cpu['force'], gpu['force'])), - stress_kbar=max(abs(x-y) for x, y in zip(cpu['stress'], gpu['stress'])), - same_outer_steps=cpu['outer_steps'] == gpu['outer_steps']) - comparisons.append(comparison) - (args.output/'comparison.json').write_text(json.dumps(comparisons, indent=2)+'\n') - if (comparison['energy_ev_per_atom'] > 1e-6 or comparison['force_ev_angstrom'] > 1e-4 - or comparison['stress_kbar'] > 1e-3 or not comparison['same_outer_steps']): - raise RuntimeError(f'Backend comparison failed: {comparison}') - - -if __name__ == '__main__': - main() diff --git a/tests/lcao_hse_gpu/build.sbatch.in b/tests/lcao_hse_gpu/build.sbatch.in deleted file mode 100644 index a8b90abc51b..00000000000 --- a/tests/lcao_hse_gpu/build.sbatch.in +++ /dev/null @@ -1,88 +0,0 @@ -#!/usr/bin/env bash -#SBATCH --job-name=@JOB_NAME@ -#SBATCH --partition=@PARTITION@ -#SBATCH --qos=@QOS@ -#SBATCH --account=abacus-group -#SBATCH --nodes=@NODES@ -#SBATCH --ntasks=@TASKS@ -#SBATCH --ntasks-per-node=@TASKS_PER_NODE@ -#SBATCH --gpus-per-node=@GPUS_PER_NODE@ -#SBATCH --time=@TIME@ -#SBATCH --output=@OUTPUT@ -#SBATCH --export=NIL - -set -euo pipefail - -export HOME=@HOME@ -export USER=${SLURM_JOB_USER:?} -export LOGNAME=$USER -export PATH=/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin - -source @CONTROL@/modules.sh -module -t list 2>&1 | tee @RESULTS@/modules.txt -command -v cmake make mpicxx nvcc | tee @RESULTS@/tools.txt -nvidia-smi dmon -s pucvmte -o T > @RESULTS@/dmon-${SLURM_JOB_ID}.log & - -# Header-only dependencies pinned by toolchain/scripts/package_versions.sh. -# Reuse the cereal installation settings from stage4/install_cereal.sh. -for package in libri libcomm cereal; do - mkdir -p @BUILD@/hse-deps/$package - tar -xzf @CONTROL@/$package.tar.gz --strip-components=1 -C @BUILD@/hse-deps/$package -done -cmake -S @BUILD@/hse-deps/cereal -B @BUILD@/hse-deps/cereal-build \ - -DCMAKE_INSTALL_PREFIX=@BUILD@/hse-deps/cereal-install \ - -DCMAKE_INSTALL_LIBDIR=lib -DBUILD_DOC=OFF -DBUILD_SANDBOX=OFF \ - -DBUILD_TESTS=OFF -DTHREAD_SAFE=ON -DSKIP_PORTABILITY_TEST=ON \ - -DSKIP_PERFORMANCE_COMPARISON=ON -DJUST_INSTALL_CEREAL=ON -cmake --install @BUILD@/hse-deps/cereal-build - -mkdir -p @BUILD@/ddla-source - tar -xzf @CONTROL@/ddla.tar.gz --strip-components=1 -C @BUILD@/ddla-source -cmake -S @BUILD@/ddla-source -B @BUILD@/ddla-build \ - -DCMAKE_BUILD_TYPE=Release -DCMAKE_CUDA_ARCHITECTURES=70 \ - -DCUBLAS_LIBRARY=/opt/devtools/nvidia/cuda-12.9.1/lib64/libcublas.so \ - -DCURAND_LIBRARY=/opt/devtools/nvidia/cuda-12.9.1/lib64/libcurand.so \ - -DCUSOLVER_LIBRARY=/opt/devtools/nvidia/cuda-12.9.1/lib64/libcusolver.so \ - -DDDLA_USE_CUDA=ON -DDDLA_USE_CCL=OFF -DBUILD_TESTS=OFF \ - -DCMAKE_INSTALL_PREFIX=@INSTALL@/ddla -cmake --build @BUILD@/ddla-build --parallel 16 -cmake --install @BUILD@/ddla-build -# LibRI adapter uses an obsolete DDLA diagnostic helper name. -sed -i 's/DEVICE_CHECK_RI/DEVICE_CHECK/g' @BUILD@/hse-deps/libri/include/RI/global/gpu/GPU_Backend.h -sed -i 's/LRI_Cal_Aux::destroy_lock_result(lock_Ds_result_add_map, Ds_result);/LRI_Cal_Aux::destroy_lock(lock_Ds_result_add_map); LRI_Cal_Aux::destroy_Ds_result(Ds_result);/' @BUILD@/hse-deps/libri/include/RI/ri/LRI-cal_loop3_GPU.hpp -# Accumulate into the existing local map, then apply fac_add_Ds once at the final merge. -sed -i 's/rDs_output.download(Ds_result, queue);/rDs_output.download(Ds_result_thread, queue);/g' @BUILD@/hse-deps/libri/include/RI/ri/LRI-cal_loop3_GPU.hpp -cmake -S @SOURCE@ -B @BUILD@ -G "Unix Makefiles" \ - -DCMAKE_BUILD_TYPE=Release \ - -DCMAKE_INSTALL_PREFIX=@INSTALL@ \ - -DCMAKE_CUDA_ARCHITECTURES=70 \ - -DENABLE_LIBRI_GPU=ON \ - -DDDLA_INCLUDE_DIR=@INSTALL@/ddla/include \ - -DDDLA_LIBRARY=@INSTALL@/ddla/lib/libddla.so \ - -DENABLE_LIBXC=ON -DENABLE_LIBRI=ON \ - -DLIBRI_DIR=@BUILD@/hse-deps/libri \ - -DLIBCOMM_DIR=@BUILD@/hse-deps/libcomm \ - -Dcereal_DIR=@BUILD@/hse-deps/cereal-install/lib/cmake/cereal \ - -DUSE_CUDA=ON \ - -DUSE_CUDA_MPI=ON \ - -DENABLE_CUSOLVERMP=ON \ - -DENABLE_CUBLASMP=ON \ - -DENABLE_NCCL_PARALLEL_DEVICE=ON \ - -DBUILD_TESTING=OFF \ - -DGIT_SUBMODULE=OFF 2>&1 | tee @RESULTS@/configure.log -cmake --build @BUILD@ --parallel 32 2>&1 | tee @RESULTS@/build.log -cmake --install @BUILD@ 2>&1 | tee @RESULTS@/install.log - -@INSTALL@/bin/abacus --info | tee @RESULTS@/abacus-info.txt -ldd @INSTALL@/bin/abacus | tee @RESULTS@/ldd.txt -if grep -q 'not found' @RESULTS@/ldd.txt; then - echo "ERROR: unresolved dynamic dependencies" >&2 - sed -n 's/.*not found.*/&/p' @RESULTS@/ldd.txt >&2 - exit 1 -fi -cp @BUILD@/CMakeCache.txt @RESULTS@/CMakeCache.txt - -mkdir -p "$HOME/abacus_hse_gpu_experiment/reuse/lcao" -cp @INSTALL@/bin/abacus "$HOME/abacus_hse_gpu_experiment/reuse/lcao/abacus" -cp -a @INSTALL@/ddla "$HOME/abacus_hse_gpu_experiment/reuse/lcao/" -exit 0 diff --git a/tests/lcao_hse_gpu/config.ini b/tests/lcao_hse_gpu/config.ini deleted file mode 100644 index a47520fee7d..00000000000 --- a/tests/lcao_hse_gpu/config.ini +++ /dev/null @@ -1,50 +0,0 @@ -[site] -name = Open Source Supercomputing Center of SAI -url = https://www.open-sai.com/ -acknowledgement = Computing resources were provided by - -[remote] -host = c0.sai.ai-4s.com -port = 12022 -user = abacususer01 -project_root = ~/abacus_hse_gpu_experiment -allowed_project_roots = /home, /org - -[cluster] -partition = 16V100 -mapping_root = /opt/sai_config/mps_mapping.d -disable_nccl_ib = false -poll_seconds = 10 - -[build] -qos = huge-gpu -nodes = 1 -tasks_per_node = 4 -gpus_per_node = 4 -time_seconds = 3600 - -[resource.gpu1] -qos = flood-1o2gpu -nodes = 1 -tasks_per_node = 1 -gpus_per_node = 1 -time_seconds = 1800 -parallelism = 1 - -[case.001] -suite = 08_RI -name = scf_hse_gamma -resource = gpu1 -runner = autotest_gpu - -[case.002] -suite = 08_RI -name = scf_hse_complex_multik -resource = gpu1 -runner = autotest_gpu - -[case.003] -suite = 08_RI -name = scf_hse_spin2_gamma -resource = gpu1 -runner = autotest_gpu From 2ca0421cde1d7ca9a83ba3efb46a36074300dd6d Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Sun, 18 Oct 2026 09:33:37 +0800 Subject: [PATCH 10/24] style(libri): align GPU dependency adapter formatting --- cmake/modules/FindLibRI.cmake | 4 +++- tests/lcao_hse_gpu/libri-compat.patch | 5 +++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/cmake/modules/FindLibRI.cmake b/cmake/modules/FindLibRI.cmake index 0de205d9d31..5417d387544 100644 --- a/cmake/modules/FindLibRI.cmake +++ b/cmake/modules/FindLibRI.cmake @@ -15,7 +15,9 @@ find_path(LIBRI_DIR if(NOT LIBRI_DIR) if(ENABLE_LIBRI_GPU) - message(FATAL_ERROR "ENABLE_LIBRI_GPU requires an explicitly supplied, patched LIBRI_DIR; see docs/advanced/install.md") + message(FATAL_ERROR + "ENABLE_LIBRI_GPU requires an explicitly supplied, patched LIBRI_DIR; " + "see docs/advanced/install.md") endif() include(FetchContent) FetchContent_Declare( diff --git a/tests/lcao_hse_gpu/libri-compat.patch b/tests/lcao_hse_gpu/libri-compat.patch index 3e17c0c3a93..5ad3bcd365c 100644 --- a/tests/lcao_hse_gpu/libri-compat.patch +++ b/tests/lcao_hse_gpu/libri-compat.patch @@ -77,6 +77,7 @@ @@ -2681 +2681 @@ - rDs_output.download(Ds_result, queue); + rDs_output.download(Ds_result_thread, queue); -@@ -2693 +2693 @@ +@@ -2693 +2693,2 @@ - LRI_Cal_Aux::destroy_lock_result(lock_Ds_result_add_map, Ds_result); -+ LRI_Cal_Aux::destroy_lock(lock_Ds_result_add_map); LRI_Cal_Aux::destroy_Ds_result(Ds_result); ++ LRI_Cal_Aux::destroy_lock(lock_Ds_result_add_map); ++ LRI_Cal_Aux::destroy_Ds_result(Ds_result); From f07d8d14f39d26782490cc4bc15d0a58c5700e62 Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Sun, 18 Oct 2026 11:04:09 +0800 Subject: [PATCH 11/24] feat(lcao): select LibRI contractions with exx_device --- docs/advanced/install.md | 9 ++- .../module_parameter/input_parameter.h | 1 + .../module_parameter/read_inp_exx_dftu.cpp | 37 ++++++++++++ .../test_serial/read_input_item_test.cpp | 58 +++++++++++++++++++ source/source_lcao/setup_exx.cpp | 2 +- 5 files changed, 103 insertions(+), 4 deletions(-) diff --git a/docs/advanced/install.md b/docs/advanced/install.md index 726b15bc5e5..755cbaed969 100644 --- a/docs/advanced/install.md +++ b/docs/advanced/install.md @@ -84,11 +84,14 @@ The new EXX implementation directly depends on [LibRI](https://github.com/abacus Experimental LCAO GPU contractions additionally require `ENABLE_LIBRI_GPU=ON` (default OFF), `USE_CUDA=ON`, and a CUDA build of DDLA. Like other dependent options, GPU RI is disabled when its LibRI/CUDA prerequisites are disabled. -With GPU RI compiled in, INPUT `device gpu` selects GPU contractions and -`device cpu` retains CPU contractions. Otherwise, RI contractions remain on CPU. +INPUT `exx_device` defaults to `cpu`, retaining CPU contractions even with +`device gpu`. To use GPU contractions, set both `device gpu` and `exx_device gpu`. +The latter requires GPU RI to be compiled in. Real-space EXX symmetry reduction +is unsupported by this backend: with `symmetry 1`, use `exx_device cpu` or +`exx_symmetry_realspace 0`. The tested combination is LibRI `6ce3c8480ae38a085746dc81cc4dbb59b8074430` -and DDLA `9c404ac65148ec1b48977b846f38ba6ef9e134f7`. LibRI requires the +and DDLA `v0.0.4` (`9c404ac65148ec1b48977b846f38ba6ef9e134f7`). LibRI requires the compatibility and force-accumulation fixes in `tests/lcao_hse_gpu/libri-compat.patch`: normalize `include/RI/global/gpu/GPU_Backend.h` and `include/RI/ri/LRI-cal_loop3_GPU.hpp` to LF, then run diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index d86993b63d9..d393d60291c 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -628,6 +628,7 @@ struct Input_para ///< the evaluation of Fock exchange double exx_mixing_beta = 1.0; ///< mixing_beta for outer-loop when exx_separate_loop=1 std::string exx_real_number = "default"; ///< exx calculated in real or complex + std::string exx_device = "cpu"; ///< computing device for LCAO EXX tensor contractions double exx_pca_threshold = 0.0001; ///< threshold to screen on-site ABFs in exx double exx_c_threshold = 0.0001; ///< threshold to screen C matrix in exx double exx_v_threshold = 0.1; ///< threshold to screen C matrix in exx diff --git a/source/source_io/module_parameter/read_inp_exx_dftu.cpp b/source/source_io/module_parameter/read_inp_exx_dftu.cpp index bbecf5cace4..cfc72707342 100644 --- a/source/source_io/module_parameter/read_inp_exx_dftu.cpp +++ b/source/source_io/module_parameter/read_inp_exx_dftu.cpp @@ -237,6 +237,43 @@ void ReadInput::item_exx() sync_stringvec(input.exx_fock_lambda, para.input.exx_fock_lambda.size(), ""); this->add_item(item); } + { + Input_Item item("exx_device"); + item.annotation = "computing device for LCAO EXX tensor contractions"; + item.category = "Exact Exchange (LCAO)"; + item.type = "String"; + item.description = R"(Selects the LibRI tensor-contraction device for LCAO exact exchange. +* cpu: retain CPU contractions, including when device=gpu. +* gpu: use the experimental GPU backend; requires ENABLE_LIBRI_GPU=ON and device=gpu. +GPU contractions do not support real-space EXX symmetry reduction. With symmetry=1, +use exx_device=cpu or disable exx_symmetry_realspace.)"; + item.default_value = "cpu"; + item.set_availability("basis_type==lcao"); + read_sync_string(input.exx_device); + item.check_value = [](const Input_Item& item, const Parameter& para) { + if (para.input.exx_device != "cpu" && para.input.exx_device != "gpu") + { + ModuleBase::WARNING_QUIT("ReadInput", "exx_device must be cpu or gpu"); + } + if (para.input.exx_device == "gpu") + { + if (para.input.basis_type != "lcao" || para.input.device != "gpu") + { + ModuleBase::WARNING_QUIT("ReadInput", "exx_device=gpu requires basis_type=lcao and device=gpu"); + } +#ifndef __GPU_RI + ModuleBase::WARNING_QUIT("ReadInput", "exx_device=gpu requires a build with ENABLE_LIBRI_GPU=ON"); +#endif + if (para.input.symmetry == "1" && para.input.exx_symmetry_realspace) + { + ModuleBase::WARNING_QUIT("ReadInput", + "LibRI GPU contractions do not support real-space EXX symmetry; " + "set exx_device=cpu or exx_symmetry_realspace=0"); + } + } + }; + this->add_item(item); + } { Input_Item item("exx_pca_threshold"); item.annotation = "threshold to screen on-site ABFs in exx"; diff --git a/source/source_io/test_serial/read_input_item_test.cpp b/source/source_io/test_serial/read_input_item_test.cpp index 6107585e733..1e3c694caf8 100644 --- a/source/source_io/test_serial/read_input_item_test.cpp +++ b/source/source_io/test_serial/read_input_item_test.cpp @@ -145,6 +145,64 @@ TEST_F(InputTest, PccSelectionAndUnsupportedCombinations) EXPECT_EXIT(it->second.check_value(it->second, param), testing::ExitedWithCode(1), ""); } +TEST_F(InputTest, ExxDevice) +{ + ModuleIO::ReadInput readinput(0); + Parameter param; + auto it = find_label("exx_device", readinput.input_lists); + ASSERT_NE(it, readinput.input_lists.end()); + EXPECT_EQ(param.inp.exx_device, "cpu"); + EXPECT_EQ(it->second.default_value, "cpu"); + + auto& input = TestParameters::input(param); + input.basis_type = "lcao"; + input.device = "gpu"; + input.symmetry = "1"; + input.exx_symmetry_realspace = true; + it->second.check_value(it->second, param); + it->second.get_final_value(it->second, param); + EXPECT_EQ(it->second.final_value.str(), "cpu"); + + const auto expect_rejected = [&](const std::string& message) { + testing::internal::CaptureStdout(); + EXPECT_EXIT(it->second.check_value(it->second, param), ::testing::ExitedWithCode(1), ""); + const std::string output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, testing::HasSubstr(message)); + }; + + it->second.str_values = {"ddla"}; + it->second.read_value(it->second, param); + expect_rejected("exx_device must be cpu or gpu"); + + it->second.str_values = {"gpu"}; + it->second.read_value(it->second, param); + EXPECT_EQ(param.inp.exx_device, "gpu"); + input.basis_type = "pw"; + expect_rejected("exx_device=gpu requires basis_type=lcao and device=gpu"); + input.basis_type = "lcao"; + input.device = "cpu"; + expect_rejected("exx_device=gpu requires basis_type=lcao and device=gpu"); + input.device = "gpu"; + +#ifdef __GPU_RI + expect_rejected("LibRI GPU contractions do not support real-space EXX symmetry"); + input.exx_symmetry_realspace = false; + it->second.check_value(it->second, param); + input.exx_symmetry_realspace = true; + input.symmetry = "0"; + it->second.check_value(it->second, param); +#else + expect_rejected("exx_device=gpu requires a build with ENABLE_LIBRI_GPU=ON"); +#endif + + it->second.str_values = {"cpu"}; + it->second.read_value(it->second, param); + EXPECT_EQ(param.inp.exx_device, "cpu"); + input.symmetry = "1"; + input.exx_symmetry_realspace = true; + it->second.check_value(it->second, param); +} + TEST_F(InputTest, Item_test) { ModuleIO::ReadInput readinput(0); diff --git a/source/source_lcao/setup_exx.cpp b/source/source_lcao/setup_exx.cpp index 6599470964a..219a22356f0 100644 --- a/source/source_lcao/setup_exx.cpp +++ b/source/source_lcao/setup_exx.cpp @@ -50,7 +50,7 @@ void Exx_NAO::before_runner( Exx_Info& exx_info) { #ifdef __EXX - const bool use_gpu = inp.device == "gpu"; + const bool use_gpu = inp.exx_device == "gpu"; if (this->exd) { this->exd->set_gpu(use_gpu); From 63ca5ef073427565e2835e0e7583cb6f88ad766d Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Sun, 18 Oct 2026 22:31:10 +0800 Subject: [PATCH 12/24] test(libri): cover GPU accumulation and clarify runtime limits --- docs/advanced/install.md | 16 +++ .../module_parameter/read_inp_exx_dftu.cpp | 9 +- tests/lcao_hse_gpu/cmake/CMakeLists.txt | 16 +++ .../cmake/test_libri_gpu_accumulation.cpp | 122 ++++++++++++++++++ 4 files changed, 161 insertions(+), 2 deletions(-) create mode 100644 tests/lcao_hse_gpu/cmake/test_libri_gpu_accumulation.cpp diff --git a/docs/advanced/install.md b/docs/advanced/install.md index 755cbaed969..a4aafa470c6 100644 --- a/docs/advanced/install.md +++ b/docs/advanced/install.md @@ -90,6 +90,19 @@ The latter requires GPU RI to be compiled in. Real-space EXX symmetry reduction is unsupported by this backend: with `symmetry 1`, use `exx_device cpu` or `exx_symmetry_realspace 0`. +The selector applies to the LCAO Kohn-Sham exchange engine. Independent LR and +RDMFT exchange engines retain CPU contractions; LR inherits the selected backend +only when it reuses the Kohn-Sham exchange object. + +With the LibRI version below, all MPI ranks on a node must see the same allocated +GPU list, with at least as many visible GPUs as local ranks. Giving each rank its +own single-GPU `CUDA_VISIBLE_DEVICES` list fails LibRI's GPU-count check when +multiple ranks share a node. For example, on one node allocated GPUs 0 and 1: + +```bash +CUDA_VISIBLE_DEVICES=0,1 mpirun -np 2 abacus +``` + The tested combination is LibRI `6ce3c8480ae38a085746dc81cc4dbb59b8074430` and DDLA `v0.0.4` (`9c404ac65148ec1b48977b846f38ba6ef9e134f7`). LibRI requires the compatibility and force-accumulation fixes in `tests/lcao_hse_gpu/libri-compat.patch`: @@ -106,6 +119,9 @@ cmake -B build -DENABLE_LIBRI=ON -DUSE_CUDA=ON -DENABLE_LIBRI_GPU=ON \ `CMAKE_PREFIX_PATH` can also locate DDLA. A configure-time compile/link check rejects incompatible GPU APIs, but does not validate numerical correctness. +The standalone tests in `tests/lcao_hse_gpu/cmake` additionally accept +`ENABLE_LIBRI_GPU_TESTS=ON` with `ENABLE_LIBRI_GPU=ON` to check real and complex +GPU contractions with signed accumulation. This optional CTest requires one GPU. This experimental backend has additional GPU memory requirements; its temporary buffer limit is not a bound on total GPU memory. diff --git a/source/source_io/module_parameter/read_inp_exx_dftu.cpp b/source/source_io/module_parameter/read_inp_exx_dftu.cpp index cfc72707342..159787b1135 100644 --- a/source/source_io/module_parameter/read_inp_exx_dftu.cpp +++ b/source/source_io/module_parameter/read_inp_exx_dftu.cpp @@ -242,11 +242,16 @@ void ReadInput::item_exx() item.annotation = "computing device for LCAO EXX tensor contractions"; item.category = "Exact Exchange (LCAO)"; item.type = "String"; - item.description = R"(Selects the LibRI tensor-contraction device for LCAO exact exchange. + item.description = R"(Selects the LibRI tensor-contraction device for LCAO Kohn-Sham exact exchange. * cpu: retain CPU contractions, including when device=gpu. * gpu: use the experimental GPU backend; requires ENABLE_LIBRI_GPU=ON and device=gpu. GPU contractions do not support real-space EXX symmetry reduction. With symmetry=1, -use exx_device=cpu or disable exx_symmetry_realspace.)"; +use exx_device=cpu or disable exx_symmetry_realspace. +Independent LR and RDMFT exchange engines retain CPU contractions. LR inherits +this selection only when it reuses the Kohn-Sham exchange object. +For GPU contractions, all MPI ranks on a node must see the same allocated GPU +list, with at least as many visible GPUs as local ranks; per-rank single-GPU +masking is unsupported.)"; item.default_value = "cpu"; item.set_availability("basis_type==lcao"); read_sync_string(input.exx_device); diff --git a/tests/lcao_hse_gpu/cmake/CMakeLists.txt b/tests/lcao_hse_gpu/cmake/CMakeLists.txt index 500d5fde3f0..3bf7726bbf5 100644 --- a/tests/lcao_hse_gpu/cmake/CMakeLists.txt +++ b/tests/lcao_hse_gpu/cmake/CMakeLists.txt @@ -5,6 +5,10 @@ project(libri_consumer_check LANGUAGES CXX) list(PREPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_LIST_DIR}/../../../cmake/modules") option(ENABLE_OPENMP "Enable OpenMP" ON) option(ENABLE_LIBRI_GPU "Check the LibRI CUDA dependency interface" OFF) +option(ENABLE_LIBRI_GPU_TESTS "Run numerical LibRI tests requiring a GPU" OFF) +if(ENABLE_LIBRI_GPU_TESTS AND NOT ENABLE_LIBRI_GPU) + message(FATAL_ERROR "ENABLE_LIBRI_GPU_TESTS requires ENABLE_LIBRI_GPU=ON") +endif() find_package(MPI REQUIRED COMPONENTS CXX) if(ENABLE_OPENMP) find_package(OpenMP REQUIRED COMPONENTS CXX) @@ -25,3 +29,15 @@ target_compile_definitions(libri_consumer PRIVATE $<$:EXPECT_LIBRI_GPU>) enable_testing() add_test(NAME libri_consumer COMMAND libri_consumer) + +if(ENABLE_LIBRI_GPU_TESTS) + add_executable(libri_gpu_accumulation test_libri_gpu_accumulation.cpp) + target_link_libraries(libri_gpu_accumulation PRIVATE LibRI::LibRI) + add_test(NAME libri_gpu_accumulation + COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 1 ${MPIEXEC_PREFLAGS} + $ ${MPIEXEC_POSTFLAGS}) + set_tests_properties(libri_gpu_accumulation PROPERTIES + LABELS gpu + TIMEOUT 120 + ENVIRONMENT "OMP_NUM_THREADS=1;OPENBLAS_NUM_THREADS=1;MKL_NUM_THREADS=1") +endif() diff --git a/tests/lcao_hse_gpu/cmake/test_libri_gpu_accumulation.cpp b/tests/lcao_hse_gpu/cmake/test_libri_gpu_accumulation.cpp new file mode 100644 index 00000000000..137247c260f --- /dev/null +++ b/tests/lcao_hse_gpu/cmake/test_libri_gpu_accumulation.cpp @@ -0,0 +1,122 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace +{ +template +RI::Tensor make_tensor(const RI::Shape_Vector& shape, const T imaginary) +{ + RI::Tensor tensor(shape); + for (std::size_t i = 0; i < tensor.data->size(); ++i) + { + (*tensor.data)[i] = T(0.5 + 0.07 * i) + imaginary * T(0.1 + 0.03 * i); + } + return tensor; +} + +template +void check_accumulation(const char* scalar_name, const T imaginary) +{ + using LRI = RI::LRI; + using Tensors = std::map>>; + const typename LRI::TAC key = {0, {0, 0, 0}}; + const std::map> atoms = {{0, {0.0, 0.0, 0.0}}}; + const std::array, 3> lattice = {{{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}}; + const std::vector labels(RI::Label::array_ab_ab.begin(), RI::Label::array_ab_ab.end()); + LRI lri; + lri.set_parallel(MPI_COMM_WORLD, atoms, lattice, {1, 1, 1}, labels); + for (const auto label: RI::Label::array_ab) + { + const bool coefficient = label == RI::Label::ab::a || label == RI::Label::ab::b; + const RI::Shape_Vector shape = coefficient ? RI::Shape_Vector{2, 2, 2} : RI::Shape_Vector{2, 2}; + Tensors tensors; + tensors[0][key] = make_tensor(shape, imaginary); + lri.set_tensors_map2(tensors, {label}); + } + + const RI::Tensor seed = make_tensor({2, 2}, imaginary); + for (const auto label: labels) + { + Tensors reference; + lri.cal_mode = RI::LRI_Cal_Mode::CPU; + lri.cal_loop3({label}, reference, 1.0); + const RI::Tensor contribution = reference.at(0).at(key).copy(); + const double norm = contribution.norm(2); + if (!std::isfinite(norm) || norm < 1e-8) + { + throw std::runtime_error("Empty or invalid CPU reference contraction"); + } + + for (const auto mode: {RI::LRI_Cal_Mode::CPU, RI::LRI_Cal_Mode::GPU}) + { + lri.cal_mode = mode; + Tensors result; + // Tensor copies share storage unless copy() is used explicitly. + result[0][key] = seed.copy(); + double accumulated_factor = 0.0; + for (const double factor: {-1.0, 1.0, 1.0, -1.0}) + { + lri.cal_loop3({label}, result, factor); + accumulated_factor += factor; + const RI::Tensor& actual = result.at(0).at(key); + if (result.size() != 1 || result.at(0).size() != 1 || actual.shape.size() != 2 || actual.shape[0] != 2 + || actual.shape[1] != 2) + { + throw std::runtime_error("Unexpected contraction output shape"); + } + for (std::size_t i = 0; i < seed.data->size(); ++i) + { + const T expected = (*seed.data)[i] + accumulated_factor * (*contribution.data)[i]; + const double error = std::abs((*actual.data)[i] - expected); + const double tolerance = 1e-11 * (1.0 + std::abs(expected)); + if (!std::isfinite(error) || error > tolerance) + { + std::cerr << scalar_name << " label=" << static_cast(label) + << " mode=" << static_cast(mode) << " factor=" << factor << " element=" << i + << " error=" << error << '\n'; + throw std::runtime_error("Signed accumulation differs from the CPU reference"); + } + } + } + } + } + std::cout << scalar_name << ": all 18 contractions preserve signed accumulation\n"; +} +} // namespace + +int main(int argc, char** argv) +{ + // This standalone dependency test has no ABACUS runtime to initialize MPI. + int provided = 0; + if (MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &provided) != MPI_SUCCESS) + { + return 1; + } + int status = 0; + try + { + const int ranks = RI::MPI_Wrapper::mpi_get_size(MPI_COMM_WORLD); + if (ranks != 1 || provided < MPI_THREAD_MULTIPLE) + { + throw std::runtime_error("This test requires one MPI rank and MPI_THREAD_MULTIPLE"); + } + check_accumulation("real", 0.0); + const std::complex imaginary(0.0, 1.0); + check_accumulation>("complex", imaginary); + } + catch (const std::exception& error) + { + std::cerr << error.what() << '\n'; + status = 1; + } + MPI_Finalize(); + return status; +} From 7acdd18ec3cf3c05953554d774732a923ba27e89 Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Mon, 19 Oct 2026 00:21:51 +0800 Subject: [PATCH 13/24] perf(xc): distribute Libxc batches across OpenMP threads --- source/source_hamilt/module_xc/libxc_pot.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/source/source_hamilt/module_xc/libxc_pot.cpp b/source/source_hamilt/module_xc/libxc_pot.cpp index eb4275cf5e3..361cdca5914 100644 --- a/source/source_hamilt/module_xc/libxc_pot.cpp +++ b/source/source_hamilt/module_xc/libxc_pot.cpp @@ -65,7 +65,7 @@ void XC_Functional_Libxc::gradcorr_ncgga_sf_libxc(const std::vector& func_i std::vector vsigma(3 * nrxx); constexpr int nr_batch_size = 1024; #ifdef _OPENMP -#pragma omp parallel for schedule(static, nr_batch_size) +#pragma omp parallel for schedule(static) #endif for (int ir_start = 0; ir_start < static_cast(nrxx); ir_start += nr_batch_size) { @@ -319,7 +319,7 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( case XC_FAMILY_LDA: { constexpr int nr_batch_size = 1024; #ifdef _OPENMP -#pragma omp parallel for schedule(static, nr_batch_size) +#pragma omp parallel for schedule(static) #endif for (int ir_start = 0; ir_start < nrxx; ir_start += nr_batch_size) { @@ -337,7 +337,7 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( case XC_FAMILY_HYB_GGA: { constexpr int nr_batch_size = 1024; #ifdef _OPENMP -#pragma omp parallel for schedule(static, nr_batch_size) +#pragma omp parallel for schedule(static) #endif for (int ir_start = 0; ir_start < nrxx; ir_start += nr_batch_size) { @@ -574,7 +574,7 @@ std::tuple XC_Functional_Li ModuleBase::timer::start("Libxc","xc_mgga_exc_vxc"); constexpr int nr_batch_size = 1024; #ifdef _OPENMP - #pragma omp parallel for schedule(static, nr_batch_size) + #pragma omp parallel for schedule(static) #endif for( int ir_start = 0; ir_start < nrxx; ir_start += nr_batch_size ) { From 4cbebe880cd6bf3d9c29f2cf676dae34f4a1cea4 Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Mon, 19 Oct 2026 20:38:01 +0800 Subject: [PATCH 14/24] perf(ri): batch GPU radial table quadrature with bounded workspace --- source/source_lcao/center2orb_orb11.cpp | 244 ++++++++++++++++++ source/source_lcao/center2orb_orb11.h | 3 + source/source_lcao/center2orb_orb21.cpp | 55 ++-- source/source_lcao/center2orb_orb21.h | 4 + source/source_lcao/module_ri/ewald_vq.hpp | 2 +- source/source_lcao/module_ri/exx_lri.hpp | 9 +- source/source_lcao/module_ri/lri_cv.h | 3 +- source/source_lcao/module_ri/lri_cv.hpp | 7 +- .../source_lcao/module_ri/matrix_orbs11.cpp | 19 +- source/source_lcao/module_ri/matrix_orbs11.h | 1 + .../source_lcao/module_ri/matrix_orbs21.cpp | 19 +- source/source_lcao/module_ri/matrix_orbs21.h | 1 + .../source_lcao/module_ri/rpa_lri_overlap.cpp | 2 +- source/source_lcao/test/CMakeLists.txt | 18 ++ .../test/test_center2orb_orb11.cpp | 123 +++++++++ 15 files changed, 468 insertions(+), 42 deletions(-) create mode 100644 source/source_lcao/test/test_center2orb_orb11.cpp diff --git a/source/source_lcao/center2orb_orb11.cpp b/source/source_lcao/center2orb_orb11.cpp index 7241aaba141..7568d8dae3e 100644 --- a/source/source_lcao/center2orb_orb11.cpp +++ b/source/source_lcao/center2orb_orb11.cpp @@ -7,11 +7,26 @@ #include "center2orb.h" #include "source_base/constants.h" +#include "source_base/math_integral.h" +#include "source_base/mathzone_add1.h" +#include "source_base/module_external/blas_connector.h" +#include "source_base/timer.h" + +#ifdef __CUDA +#include "source_base/module_device/memory_op.h" + +#include +#include +#endif #include "source_base/math_polyint.h" #include "source_base/sph_bessel_recursive.h" #include "source_base/ylm.h" +#include #include +#include +#include +#include Center2_Orb::Orb11::Orb11(const Numerical_Orbital_Lm& nA_in, const Numerical_Orbital_Lm& nB_in, @@ -86,6 +101,235 @@ void Center2_Orb::Orb11::init_radial_table(const std::set& radials) } } +namespace +{ +double radial_product(const Numerical_Orbital_Lm& left, const Numerical_Orbital_Lm& right, const int ik) +{ + if (!left.get_psif().empty() && !right.get_psi_k2().empty()) + { + return left.getPsif(ik) * right.getPsi_k2(ik); + } + if (!left.get_psi_k().empty() && !right.get_psi_k().empty()) + { + return left.getPsi_k(ik) * right.getPsi_k(ik); + } + if (!left.get_psi_k2().empty() && !right.get_psif().empty()) + { + return left.getPsi_k2(ik) * right.getPsif(ik); + } + return 0.0; +} +} // namespace + +void Center2_Orb::Orb11::init_radial_table(const std::vector& orbitals, const bool use_gpu) +{ + ModuleBase::timer::start("Center2_Orb", "init_radial_table_batch"); +#ifndef __CUDA + if (use_gpu) + { + throw std::invalid_argument("GPU radial tables require a CUDA build"); + } +#else + if (use_gpu) + { + container::kernels::createGpuBlasHandle(); + } +#endif + // Group by quadrature grid and angular momentum. All orbitals in a batch + // belong to one Matrix_Orbs object and share its immutable Bessel table. + using Grid = std::tuple; + std::map> groups; + for (Orb11* orbital: orbitals) + { + if (orbital->psb_ != orbitals.front()->psb_) + { + throw std::invalid_argument("A radial batch must share its Bessel table"); + } + const int rmesh = Center2_Orb::get_rmesh(orbital->nA.getRcut(), orbital->nB.getRcut(), orbital->dr_); + const int kmesh = orbital->nA.getNk(); + const double dk = orbital->nA.getDk(); + const int first_l = std::abs(orbital->nA.getL() - orbital->nB.getL()); + const int last_l = orbital->nA.getL() + orbital->nB.getL(); + for (int l = first_l; l <= last_l; l += 2) + { + orbital->Table_r[l].resize(rmesh); + orbital->Table_dr[l].resize(rmesh); + const Grid grid(l, kmesh, rmesh, dk); + groups[grid].push_back(orbital); + } + } + for (const auto& group: groups) + { + const int l = std::get<0>(group.first); + const int kmesh = std::get<1>(group.first); + const int rmesh = std::get<2>(group.first); + const double dk = std::get<3>(group.first); + const auto& tasks = group.second; + const auto& bessel = tasks.front()->psb_->get_jlx(); + const double* kpoint = tasks.front()->nA.getKpoint(); + // Bound all three device buffers together: at most 160 MiB, released + // before returning to the SCF/LibRI contraction stages. + constexpr std::size_t input_limit = 16 * 1024 * 1024 / sizeof(double); + constexpr std::size_t table_limit = 128 * 1024 * 1024 / sizeof(double); + constexpr std::size_t output_limit = 16 * 1024 * 1024 / sizeof(double); + if (kmesh < 3 || kmesh % 2 == 0 || static_cast(kmesh) > input_limit) + { + throw std::invalid_argument("Unsupported radial quadrature grid for batched tables"); + } + const int radial_block = std::min(rmesh, table_limit / (2 * kmesh)); + const std::size_t batch_limit = std::min(input_limit / kmesh, output_limit / (2 * radial_block)); + const int batch_size = std::min(tasks.size(), std::min(128, batch_limit)); + const std::size_t input_size = static_cast(batch_size) * kmesh; + const std::size_t table_size = static_cast(2 * radial_block) * kmesh; + const std::size_t output_size = static_cast(batch_size) * 2 * radial_block; + std::vector input(input_size); + std::vector table(table_size); + std::vector output(output_size); + double* input_ptr = input.data(); + double* table_ptr = table.data(); + double* output_ptr = output.data(); +#ifdef __CUDA + std::unique_ptr workspace; + if (use_gpu) + { + const int64_t workspace_size = input_size + table_size + output_size; + const container::TensorShape shape({workspace_size}); + workspace.reset( + new container::Tensor(container::DataType::DT_DOUBLE, container::DeviceType::GpuDevice, shape)); + input_ptr = workspace->data(); + table_ptr = input_ptr + input_size; + output_ptr = table_ptr + table_size; + } +#endif + for (int radial_start = 0; radial_start < rmesh; radial_start += radial_block) + { + const int rows = std::min(radial_block, rmesh - radial_start); + const int columns = 2 * rows; + const double derivative_factor = -(l + 1.0) / (2 * l + 1.0); + const double lower_factor = l / (l + 1.0); + for (int ir = 0; ir < rows; ++ir) + { + const int radial = radial_start + ir; + for (int ik = 0; ik < kmesh; ++ik) + { + const double lower = l == 0 ? 0.0 : bessel[l - 1][radial][ik]; + table[(2 * ir) * kmesh + ik] = bessel[l][radial][ik]; + table[(2 * ir + 1) * kmesh + ik] + = derivative_factor * kpoint[ik] * (bessel[l + 1][radial][ik] - lower_factor * lower); + } + } +#ifdef __CUDA + if (use_gpu) + { + const std::size_t count = static_cast(columns) * kmesh; + base_device::memory::synchronize_memory_op()( + table_ptr, + table.data(), + count); + } +#endif + for (std::size_t start = 0; start < tasks.size(); start += batch_size) + { + const int count = std::min(batch_size, tasks.size() - start); + for (int it = 0; it < count; ++it) + { + const Orb11& orbital = *tasks[start + it]; + for (int ik = 0; ik < kmesh; ++ik) + { + const int weight = (ik == 0 || ik == kmesh - 1) ? 1 : (ik % 2 == 1 ? 4 : 2); + const double product = radial_product(orbital.nA, orbital.nB, ik); + input[it * kmesh + ik] = product * weight * dk * ModuleBase::FOUR_PI / 3.0; + } + } +#ifdef __CUDA + if (use_gpu) + { + const std::size_t size = static_cast(count) * kmesh; + base_device::memory:: + synchronize_memory_op()(input_ptr, + input.data(), + size); + } +#endif +#ifdef __CUDA + if (use_gpu) + { + const double one = 1.0; + const double zero = 0.0; + container::kernels::blas_gemm()('T', + 'N', + columns, + count, + kmesh, + &one, + table_ptr, + kmesh, + input_ptr, + kmesh, + &zero, + output_ptr, + columns); + } + else +#endif + { + BlasConnector::gemm_cm('T', + 'N', + columns, + count, + kmesh, + 1.0, + table_ptr, + kmesh, + input_ptr, + kmesh, + 0.0, + output_ptr, + columns, + base_device::AbacusDevice_t::CpuDevice); + } +#ifdef __CUDA + if (use_gpu) + { + const std::size_t size = static_cast(count) * columns; + base_device::memory:: + synchronize_memory_op()(output.data(), + output_ptr, + size); + } +#endif + for (int it = 0; it < count; ++it) + { + auto& values = tasks[start + it]->Table_r.at(l); + auto& derivatives = tasks[start + it]->Table_dr.at(l); + for (int ir = 0; ir < rows; ++ir) + { + values[radial_start + ir] = output[it * columns + 2 * ir]; + derivatives[radial_start + ir] = output[it * columns + 2 * ir + 1]; + } + } + } + } + // Preserve the analytic R=0 convention used by cal_overlap. + if (l > 0) + { + std::vector integrand(kmesh); + for (Orb11* orbital: tasks) + { + for (int ik = 0; ik < kmesh; ++ik) + { + integrand[ik] = radial_product(orbital->nA, orbital->nB, ik) * std::pow(kpoint[ik], l); + } + double integral = 0.0; + ModuleBase::Integral::Simpson_Integral(kmesh, integrand.data(), dk, integral); + orbital->Table_r.at(l)[0] + = ModuleBase::FOUR_PI / ModuleBase::Mathzone_Add1::dualfac(2 * l + 1) * integral; + } + } + } + ModuleBase::timer::end("Center2_Orb", "init_radial_table_batch"); +} + double Center2_Orb::Orb11::cal_overlap(const ModuleBase::Vector3& RA, const ModuleBase::Vector3& RB, const int& mA, diff --git a/source/source_lcao/center2orb_orb11.h b/source/source_lcao/center2orb_orb11.h index a7871eb1728..2b718eee529 100644 --- a/source/source_lcao/center2orb_orb11.h +++ b/source/source_lcao/center2orb_orb11.h @@ -27,6 +27,9 @@ class Center2_Orb::Orb11 void init_radial_table(); void init_radial_table(const std::set& radials); // unit: Bohr/MOT.dr + /// Build full tables in bounded batches, using the same Simpson quadrature. + static void init_radial_table(const std::vector& orbitals, const bool use_gpu); + double cal_overlap(const ModuleBase::Vector3& RA, const ModuleBase::Vector3& RB, // unit: Bohr const int& mA, diff --git a/source/source_lcao/center2orb_orb21.cpp b/source/source_lcao/center2orb_orb21.cpp index 261c7402526..3acc44fdf91 100644 --- a/source/source_lcao/center2orb_orb21.cpp +++ b/source/source_lcao/center2orb_orb21.cpp @@ -20,7 +20,7 @@ Center2_Orb::Orb21::Orb21(const Numerical_Orbital_Lm& nA1_in, { } -void Center2_Orb::Orb21::init_radial_table() +void Center2_Orb::Orb21::prepare_radial_table() { const Numerical_Orbital_Lm& nA_short = (this->nA1.getNr() <= this->nA2.getNr()) ? this->nA1 : this->nA2; @@ -57,49 +57,40 @@ void Center2_Orb::Orb21::init_radial_table() this->orb11s.insert(std::make_pair(LA, Center2_Orb::Orb11(this->nA[LA], nB, psb_, this->MGT))); - this->orb11s.at(LA).init_radial_table(); } } -void Center2_Orb::Orb21::init_radial_table(const std::set& radials) +void Center2_Orb::Orb21::init_radial_table() { - const Numerical_Orbital_Lm& nA_short = (this->nA1.getNr() <= this->nA2.getNr()) ? this->nA1 : this->nA2; + this->prepare_radial_table(); + for (auto& orbital: this->orb11s) + { + orbital.second.init_radial_table(); + } +} - std::vector nA_tmp(nA_short.getNr()); - for (size_t ir = 0; ir != nA_tmp.size(); ++ir) +void Center2_Orb::Orb21::init_radial_table(const std::set& radials) +{ + this->prepare_radial_table(); + for (auto& orbital: this->orb11s) { - nA_tmp[ir] = this->nA1.getPsi(ir) * this->nA2.getPsi(ir); + orbital.second.init_radial_table(radials); } +} - const int LA1 = this->nA1.getL(); - const int LA2 = this->nA2.getL(); - for (int LA = std::abs(LA1 - LA2); LA <= LA1 + LA2; ++LA) +void Center2_Orb::Orb21::init_radial_table(const std::vector& orbitals, const bool use_gpu) +{ + std::vector pairs; + for (Orb21* orbital: orbitals) { - if ((LA - std::abs(LA1 - LA2)) % 2 == 1) // if LA+LB-LAB == odd, then Gaunt_Coefficients = 0 + // Orbital transforms grow the shared Bessel cache, so prepare them serially. + orbital->prepare_radial_table(); + for (auto& pair: orbital->orb11s) { - continue; + pairs.push_back(&pair.second); } - - this->nA[LA].set_orbital_info(nA_short.getLabel(), - nA_short.getType(), - LA, - 1, // N? - nA_short.getNr(), - nA_short.getRab(), - nA_short.getRadial(), - Numerical_Orbital_Lm::Psi_Type::Psi, - nA_tmp.data(), - nA_short.getNk(), - nA_short.getDk(), - nA_short.getDruniform(), - false, - true, - PARAM.inp.cal_force); // mohan add 2021-05-07 - - this->orb11s.insert(std::make_pair(LA, Center2_Orb::Orb11(this->nA[LA], nB, psb_, this->MGT))); - - this->orb11s.at(LA).init_radial_table(radials); } + Orb11::init_radial_table(pairs, use_gpu); } double Center2_Orb::Orb21::cal_overlap(const ModuleBase::Vector3& RA, diff --git a/source/source_lcao/center2orb_orb21.h b/source/source_lcao/center2orb_orb21.h index bc7cb5c7fa1..330aebb4197 100644 --- a/source/source_lcao/center2orb_orb21.h +++ b/source/source_lcao/center2orb_orb21.h @@ -28,6 +28,8 @@ class Center2_Orb::Orb21 void init_radial_table(); void init_radial_table(const std::set& radials); // unit: Bohr/MOT.dr + static void init_radial_table(const std::vector& orbitals, const bool use_gpu); + double cal_overlap(const ModuleBase::Vector3& RA, const ModuleBase::Vector3& RB, // unit: Bohr const int& mA1, @@ -41,6 +43,8 @@ class Center2_Orb::Orb21 const int& mB) const; private: + void prepare_radial_table(); + const Numerical_Orbital_Lm& nA1; const Numerical_Orbital_Lm& nA2; const Numerical_Orbital_Lm& nB; diff --git a/source/source_lcao/module_ri/ewald_vq.hpp b/source/source_lcao/module_ri/ewald_vq.hpp index d368b5565c7..d707d8a3152 100644 --- a/source/source_lcao/module_ri/ewald_vq.hpp +++ b/source/source_lcao/module_ri/ewald_vq.hpp @@ -68,7 +68,7 @@ void Ewald_Vq::init(const UnitCell& ucell, this->index_abfs = ModuleBase::Element_Basis_Index::construct_index(range_abfs); this->cv - .set_orbitals(ucell, orb, this->g_lcaos, this->g_abfs, this->g_abfs_ccp, kmesh_times_in, MGT_in, false); + .set_orbitals(ucell, orb, this->g_lcaos, this->g_abfs, this->g_abfs_ccp, kmesh_times_in, MGT_in, false, false); this->gaunt.create(MGT_in->Gaunt_Coefficients.getBound1(), MGT_in->Gaunt_Coefficients.getBound2(), MGT_in->Gaunt_Coefficients.getBound3()); diff --git a/source/source_lcao/module_ri/exx_lri.hpp b/source/source_lcao/module_ri/exx_lri.hpp index 2ba07c692f0..6950a1a7162 100644 --- a/source/source_lcao/module_ri/exx_lri.hpp +++ b/source/source_lcao/module_ri/exx_lri.hpp @@ -96,6 +96,11 @@ void Exx_LRI::refresh_coulomb(const UnitCell& ucell, ++object; } } +#ifdef __GPU_RI + const bool use_gpu = this->exx_lri.lri.cal_mode == RI::LRI_Cal_Mode::GPU; +#else + const bool use_gpu = false; +#endif for (const auto& settings_list : settings) { const auto previous = this->coulomb_settings.find(settings_list.first); @@ -113,7 +118,7 @@ void Exx_LRI::refresh_coulomb(const UnitCell& ucell, object.abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp( this->abfs, settings_list.second.second, this->info.ccp_rmesh_times); object.cv.set_orbitals(ucell, orb, this->lcaos, this->abfs, object.abfs_ccp, - this->info.kmesh_times, this->MGT, settings_list.second.first); + this->info.kmesh_times, this->MGT, settings_list.second.first, use_gpu); object.cv.set_info_ri(&this->info); if (settings_list.first == Conv_Coulomb_Pot_K::Coulomb_Method::Ewald) { @@ -190,7 +195,7 @@ void Exx_LRI::init_spencer( this->exx_objs[Conv_Coulomb_Pot_K::Coulomb_Method::Center2].abfs_ccp, this->info.kmesh_times, this->MGT, - center2_settings->second.first); + center2_settings->second.first, false); this->exx_objs[Conv_Coulomb_Pot_K::Coulomb_Method::Center2].cv.set_info_ri(&this->info); ModuleBase::timer::end("Exx_LRI", "init_spencer"); diff --git a/source/source_lcao/module_ri/lri_cv.h b/source/source_lcao/module_ri/lri_cv.h index 12359d3b5e3..0df1662a274 100644 --- a/source/source_lcao/module_ri/lri_cv.h +++ b/source/source_lcao/module_ri/lri_cv.h @@ -44,7 +44,8 @@ class LRI_CV const std::vector>> &abfs_ccp_in, const double &kmesh_times, std::shared_ptr MGT, - const bool& init_C); + const bool& init_C, + const bool use_gpu); inline std::map>> cal_Vs( const UnitCell &ucell, diff --git a/source/source_lcao/module_ri/lri_cv.hpp b/source/source_lcao/module_ri/lri_cv.hpp index 49f197816c0..8fc0221f68f 100644 --- a/source/source_lcao/module_ri/lri_cv.hpp +++ b/source/source_lcao/module_ri/lri_cv.hpp @@ -45,7 +45,8 @@ void LRI_CV::set_orbitals( const std::vector>> &abfs_ccp_in, const double &kmesh_times, std::shared_ptr MGT, - const bool& init_C) + const bool& init_C, + const bool use_gpu) { ModuleBase::TITLE("LRI_CV", "set_orbitals"); ModuleBase::timer::start("LRI_CV", "set_orbitals"); @@ -74,9 +75,9 @@ void LRI_CV::set_orbitals( this->abfs_ccp, this->lcaos, this->lcaos, ucell, orb, kmesh_times); - this->m_abfs_abfs.init_radial_table(); + this->m_abfs_abfs.init_radial_table(use_gpu); if (init_C) { - this->m_abfslcaos_lcaos.init_radial_table(); + this->m_abfslcaos_lcaos.init_radial_table(use_gpu); } ModuleBase::timer::end("LRI_CV", "set_orbitals"); diff --git a/source/source_lcao/module_ri/matrix_orbs11.cpp b/source/source_lcao/module_ri/matrix_orbs11.cpp index 2673535cddb..a7b90200ab2 100644 --- a/source/source_lcao/module_ri/matrix_orbs11.cpp +++ b/source/source_lcao/module_ri/matrix_orbs11.cpp @@ -91,22 +91,39 @@ void Matrix_Orbs11::init_radial(const LCAO_Orbitals& orb_A, const LCAO_Orbitals& */ void Matrix_Orbs11::init_radial_table() +{ + this->init_radial_table(false); +} + +void Matrix_Orbs11::init_radial_table(const bool use_gpu) { ModuleBase::TITLE("Matrix_Orbs11", "init_radial_table"); ModuleBase::timer::start("Matrix_Orbs11", "init_radial_table"); + std::vector orbitals; for (auto& coA: center2_orb11_s) { for (auto& coB: coA.second) { for (auto& coC: coB.second) { for (auto& coD: coC.second) { for (auto& coE: coD.second) { for (auto& coF: coE.second) { - coF.second.init_radial_table(); + if (use_gpu) + { + orbitals.push_back(&coF.second); + } + else + { + coF.second.init_radial_table(); + } } } } } } } + if (use_gpu) + { + Center2_Orb::Orb11::init_radial_table(orbitals, true); + } ModuleBase::timer::end("Matrix_Orbs11", "init_radial_table"); } diff --git a/source/source_lcao/module_ri/matrix_orbs11.h b/source/source_lcao/module_ri/matrix_orbs11.h index 0fbb6ca1cc7..87dd7208191 100644 --- a/source/source_lcao/module_ri/matrix_orbs11.h +++ b/source/source_lcao/module_ri/matrix_orbs11.h @@ -29,6 +29,7 @@ class Matrix_Orbs11 const double kmesh_times); // extend Kcut, keep dK void init_radial_table(); + void init_radial_table(const bool use_gpu); void init_radial_table(const std::map>>& Rs); // unit: ucell.lat0 enum class Matrix_Order diff --git a/source/source_lcao/module_ri/matrix_orbs21.cpp b/source/source_lcao/module_ri/matrix_orbs21.cpp index 963d7786b42..31491a04d0c 100644 --- a/source/source_lcao/module_ri/matrix_orbs21.cpp +++ b/source/source_lcao/module_ri/matrix_orbs21.cpp @@ -117,9 +117,15 @@ void Matrix_Orbs21::init_radial(const std::vectorinit_radial_table(false); +} + +void Matrix_Orbs21::init_radial_table(const bool use_gpu) { ModuleBase::TITLE("Matrix_Orbs21", "init_radial_table"); ModuleBase::timer::start("Matrix_Orbs21", "init_radial_table"); + std::vector orbitals; for (auto& coA: center2_orb21_s) { for (auto& coB: coA.second) @@ -136,7 +142,14 @@ void Matrix_Orbs21::init_radial_table() { for (auto& coH: coG.second) { - coH.second.init_radial_table(); + if (use_gpu) + { + orbitals.push_back(&coH.second); + } + else + { + coH.second.init_radial_table(); + } } } } @@ -145,6 +158,10 @@ void Matrix_Orbs21::init_radial_table() } } } + if (use_gpu) + { + Center2_Orb::Orb21::init_radial_table(orbitals, true); + } ModuleBase::timer::end("Matrix_Orbs21", "init_radial_table"); } diff --git a/source/source_lcao/module_ri/matrix_orbs21.h b/source/source_lcao/module_ri/matrix_orbs21.h index 10991a4cade..707eca20c3d 100644 --- a/source/source_lcao/module_ri/matrix_orbs21.h +++ b/source/source_lcao/module_ri/matrix_orbs21.h @@ -29,6 +29,7 @@ class Matrix_Orbs21 const double kmesh_times); // extend Kcut, keep dK void init_radial_table(); + void init_radial_table(const bool use_gpu); void init_radial_table(const std::map>>& Rs); // unit: ucell.lat0 enum class Matrix_Order diff --git a/source/source_lcao/module_ri/rpa_lri_overlap.cpp b/source/source_lcao/module_ri/rpa_lri_overlap.cpp index 1bb5e27ad6f..1409bacc966 100644 --- a/source/source_lcao/module_ri/rpa_lri_overlap.cpp +++ b/source/source_lcao/module_ri/rpa_lri_overlap.cpp @@ -64,7 +64,7 @@ void RPA_LRI::cal_large_Cs(const UnitCell& ucell, const LCAO_Orbitals& center2_obj_it->second.abfs_ccp, this->info.kmesh_times, this->MGT, // get MGT from exx_cut_coulomb and used in `cal_abfs_overlap` - true); + true, false); const std::array period_Vs = LRI_CV_Tools::cal_latvec_range(1 + this->ccp_rmesh_times_cut, ucell, orb_cutoff_); diff --git a/source/source_lcao/test/CMakeLists.txt b/source/source_lcao/test/CMakeLists.txt index 4c993f3023b..1885e7c907b 100644 --- a/source/source_lcao/test/CMakeLists.txt +++ b/source/source_lcao/test/CMakeLists.txt @@ -3,6 +3,24 @@ abacus_disable_feature_definitions(__CUDA) abacus_disable_feature_definitions(__ROCM) if(ENABLE_LCAO) +AddTest( + TARGET MODULE_LCAO_center2orb_table + LIBS parameter base device container + SOURCES test_center2orb_orb11.cpp ../center2orb.cpp ../center2orb_orb11.cpp + ${ABACUS_SOURCE_DIR}/source_basis/module_ao/orb_atomic_lm.cpp + ${ABACUS_SOURCE_DIR}/source_basis/module_ao/orb_gaunt_table.cpp +) +if(USE_CUDA AND ENABLE_LIBRI_GPU) + AddTest( + TARGET MODULE_LCAO_center2orb_table_gpu + KEEP_FEATURE_DEFINITIONS __CUDA + LIBS parameter base device container + SOURCES test_center2orb_orb11.cpp ../center2orb.cpp ../center2orb_orb11.cpp + ${ABACUS_SOURCE_DIR}/source_basis/module_ao/orb_atomic_lm.cpp + ${ABACUS_SOURCE_DIR}/source_basis/module_ao/orb_gaunt_table.cpp + ) + set_tests_properties(MODULE_LCAO_center2orb_table_gpu PROPERTIES LABELS gpu) +endif() AddTest( TARGET MODULE_LCAO_init_dm_from_file_test LIBS parameter base device symmetry diff --git a/source/source_lcao/test/test_center2orb_orb11.cpp b/source/source_lcao/test/test_center2orb_orb11.cpp new file mode 100644 index 00000000000..617445731dd --- /dev/null +++ b/source/source_lcao/test/test_center2orb_orb11.cpp @@ -0,0 +1,123 @@ +#include "source_base/module_device/device.h" +#include "source_lcao/center2orb_orb11.h" + +#include "gtest/gtest.h" +#include +#include +#include + +namespace +{ +Numerical_Orbital_Lm make_orbital(const int l, const int index, const int nr) +{ + const double dr = 0.02; + const double cutoff = (nr - 1) * dr; + std::vector radial(nr); + std::vector rab(nr, dr); + std::vector values(nr); + for (int ir = 0; ir < nr; ++ir) + { + const double r = ir * dr; + radial[ir] = r; + values[ir] = std::pow(r, l) * std::pow(cutoff - r, 2) * std::exp(-(0.7 + 0.03 * index) * r); + } + Numerical_Orbital_Lm orbital; + const double* rab_ptr = rab.data(); + const double* radial_ptr = radial.data(); + const double* values_ptr = values.data(); + orbital.set_orbital_info("X", + 0, + l, + index, + nr, + rab_ptr, + radial_ptr, + Numerical_Orbital_Lm::Psi_Type::Psi, + values_ptr, + 65, + 0.1, + 0.01, + false, + true, + true); + return orbital; +} + +void compare_pairs(const bool use_gpu) +{ + // More than 128 same-grid pairs exercise full and partial batches. + std::vector basis; + for (int l = 0; l <= 3; ++l) + { + basis.push_back(make_orbital(l, l, 101)); + } + basis.push_back(make_orbital(2, 4, 81)); + ORB_gaunt_table gaunt; + gaunt.init_Gaunt_CH(3); + gaunt.init_Gaunt(3); + ModuleBase::Sph_Bessel_Recursive::D2 bessel; + bessel.set_dx(0.001); + bessel.cal_jlx(7, 405, 65); + std::vector> reference; + std::vector> candidate; + std::vector tasks; + for (int i = 0; i < 140; ++i) + { + const int left = i < 130 ? 1 : i % 5; + const int right = i < 130 ? 2 : (i + 2) % 5; + reference.emplace_back(new Center2_Orb::Orb11(basis[left], basis[right], &bessel, gaunt)); + candidate.emplace_back(new Center2_Orb::Orb11(basis[left], basis[right], &bessel, gaunt)); + reference.back()->init_radial_table(); + tasks.push_back(candidate.back().get()); + } + Center2_Orb::Orb11::init_radial_table(tasks, use_gpu); + const ModuleBase::Vector3 origin(0.0, 0.0, 0.0); + for (int i = 0; i < 140; ++i) + { + SCOPED_TRACE(i); + const int left = i < 130 ? 1 : i % 5; + const int right = i < 130 ? 2 : (i + 2) % 5; + const int last_left = 2 * basis[left].getL(); + const int last_right = 2 * basis[right].getL(); + for (double distance: {0.0, 0.037, 0.73, 1.91, 3.87, 4.2}) + { + const double x = 0.8 * distance; + const double y = 0.36 * distance; + const double z = 0.48 * distance; + const ModuleBase::Vector3 position(x, y, z); + for (int ma: {0, last_left}) + { + for (int mb: {0, last_right}) + { + const double expected = reference[i]->cal_overlap(origin, position, ma, mb); + const double actual = candidate[i]->cal_overlap(origin, position, ma, mb); + const double tolerance = 2e-11 * (1.0 + std::abs(expected)); + EXPECT_NEAR(actual, expected, tolerance); + const auto grad_expected = reference[i]->cal_grad_overlap(origin, position, ma, mb); + const auto grad_actual = candidate[i]->cal_grad_overlap(origin, position, ma, mb); + const double error = (grad_actual - grad_expected).norm(); + const double grad_tolerance = 2e-10 * (1.0 + grad_expected.norm()); + EXPECT_LT(error, grad_tolerance); + } + } + } + } +} +} // namespace + +TEST(Center2Orb11, BatchedCpuMatchesIndividualTables) +{ + compare_pairs(false); +} + +#ifdef __CUDA +TEST(Center2Orb11, BatchedGpuMatchesIndividualTables) +{ + const bool available = base_device::information::probe_gpu_availability(); + if (!available) + { + GTEST_SKIP() << "A CUDA device is required"; + } + compare_pairs(true); +} +#endif From d11afe60bac432cfd38a9f92fc16d746b1f72dd3 Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Mon, 19 Oct 2026 20:46:15 +0800 Subject: [PATCH 15/24] perf(ri): skip unused product-orbital origin derivatives --- source/source_lcao/center2orb_orb21.cpp | 3 +- source/source_lcao/test/CMakeLists.txt | 4 +- .../test/test_center2orb_orb11.cpp | 38 +++++++++++++++++++ 3 files changed, 41 insertions(+), 4 deletions(-) diff --git a/source/source_lcao/center2orb_orb21.cpp b/source/source_lcao/center2orb_orb21.cpp index 3acc44fdf91..8454e6c86cb 100644 --- a/source/source_lcao/center2orb_orb21.cpp +++ b/source/source_lcao/center2orb_orb21.cpp @@ -1,6 +1,5 @@ //========================================================= // AUTHOR : Peize Lin -#include "source_io/module_parameter/parameter.h" // DATE : 2016-01-24 //========================================================= @@ -53,7 +52,7 @@ void Center2_Orb::Orb21::prepare_radial_table() nA_short.getDruniform(), false, true, - PARAM.inp.cal_force); // mohan add 2021-05-07 + false); // Product orbitals use reciprocal-space gradients, not zty. this->orb11s.insert(std::make_pair(LA, Center2_Orb::Orb11(this->nA[LA], nB, psb_, this->MGT))); diff --git a/source/source_lcao/test/CMakeLists.txt b/source/source_lcao/test/CMakeLists.txt index 1885e7c907b..f7e97b7e18e 100644 --- a/source/source_lcao/test/CMakeLists.txt +++ b/source/source_lcao/test/CMakeLists.txt @@ -6,7 +6,7 @@ if(ENABLE_LCAO) AddTest( TARGET MODULE_LCAO_center2orb_table LIBS parameter base device container - SOURCES test_center2orb_orb11.cpp ../center2orb.cpp ../center2orb_orb11.cpp + SOURCES test_center2orb_orb11.cpp ../center2orb.cpp ../center2orb_orb11.cpp ../center2orb_orb21.cpp ${ABACUS_SOURCE_DIR}/source_basis/module_ao/orb_atomic_lm.cpp ${ABACUS_SOURCE_DIR}/source_basis/module_ao/orb_gaunt_table.cpp ) @@ -15,7 +15,7 @@ if(USE_CUDA AND ENABLE_LIBRI_GPU) TARGET MODULE_LCAO_center2orb_table_gpu KEEP_FEATURE_DEFINITIONS __CUDA LIBS parameter base device container - SOURCES test_center2orb_orb11.cpp ../center2orb.cpp ../center2orb_orb11.cpp + SOURCES test_center2orb_orb11.cpp ../center2orb.cpp ../center2orb_orb11.cpp ../center2orb_orb21.cpp ${ABACUS_SOURCE_DIR}/source_basis/module_ao/orb_atomic_lm.cpp ${ABACUS_SOURCE_DIR}/source_basis/module_ao/orb_gaunt_table.cpp ) diff --git a/source/source_lcao/test/test_center2orb_orb11.cpp b/source/source_lcao/test/test_center2orb_orb11.cpp index 617445731dd..560e6de565f 100644 --- a/source/source_lcao/test/test_center2orb_orb11.cpp +++ b/source/source_lcao/test/test_center2orb_orb11.cpp @@ -1,5 +1,6 @@ #include "source_base/module_device/device.h" #include "source_lcao/center2orb_orb11.h" +#include "source_lcao/center2orb_orb21.h" #include "gtest/gtest.h" #include @@ -110,6 +111,43 @@ TEST(Center2Orb11, BatchedCpuMatchesIndividualTables) compare_pairs(false); } +TEST(Center2Orb21, ProductGradientMatchesFiniteDifference) +{ + const auto left = make_orbital(2, 0, 101); + const auto middle = make_orbital(1, 1, 101); + const auto right = make_orbital(3, 2, 101); + ORB_gaunt_table gaunt; + gaunt.init_Gaunt_CH(3); + gaunt.init_Gaunt(3); + ModuleBase::Sph_Bessel_Recursive::D2 bessel; + bessel.set_dx(0.001); + bessel.cal_jlx(7, 405, 65); + Center2_Orb::Orb21 orbital(left, middle, right, &bessel, gaunt); + orbital.init_radial_table(); + const ModuleBase::Vector3 origin(0.0, 0.0, 0.0); + const ModuleBase::Vector3 position(0.813, 0.321, 0.473); + const double step = 1e-5; + for (int ma: {0, 2, 4}) + { + for (int mb: {0, 1, 2}) + { + for (int mc: {0, 3, 6}) + { + const auto gradient = orbital.cal_grad_overlap(origin, position, ma, mb, mc); + auto plus = position; + auto minus = position; + plus.x += step; + minus.x -= step; + const double upper = orbital.cal_overlap(origin, plus, ma, mb, mc); + const double lower = orbital.cal_overlap(origin, minus, ma, mb, mc); + const double finite_difference = (upper - lower) / (2 * step); + const double tolerance = 2e-7 * (1.0 + std::abs(finite_difference)); + EXPECT_NEAR(gradient.x, finite_difference, tolerance); + } + } + } +} + #ifdef __CUDA TEST(Center2Orb11, BatchedGpuMatchesIndividualTables) { From 522118a63ce35a99dbc309192d06a61705d598e2 Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Mon, 19 Oct 2026 20:51:22 +0800 Subject: [PATCH 16/24] perf(ri): reuse product orbitals across batched three-center tables --- source/source_lcao/center2orb_orb21.cpp | 59 +++++++++++++------ source/source_lcao/center2orb_orb21.h | 5 +- .../test/test_center2orb_orb11.cpp | 48 +++++++++++++++ 3 files changed, 93 insertions(+), 19 deletions(-) diff --git a/source/source_lcao/center2orb_orb21.cpp b/source/source_lcao/center2orb_orb21.cpp index 8454e6c86cb..75f2177458b 100644 --- a/source/source_lcao/center2orb_orb21.cpp +++ b/source/source_lcao/center2orb_orb21.cpp @@ -21,6 +21,7 @@ Center2_Orb::Orb21::Orb21(const Numerical_Orbital_Lm& nA1_in, void Center2_Orb::Orb21::prepare_radial_table() { + auto products = std::make_shared(); const Numerical_Orbital_Lm& nA_short = (this->nA1.getNr() <= this->nA2.getNr()) ? this->nA1 : this->nA2; std::vector nA_tmp(nA_short.getNr()); @@ -38,24 +39,32 @@ void Center2_Orb::Orb21::prepare_radial_table() continue; } - this->nA[LA].set_orbital_info(nA_short.getLabel(), - nA_short.getType(), - LA, - 1, // N? - nA_short.getNr(), - nA_short.getRab(), - nA_short.getRadial(), - Numerical_Orbital_Lm::Psi_Type::Psi, - nA_tmp.data(), - nA_short.getNk(), - nA_short.getDk(), - nA_short.getDruniform(), - false, - true, - false); // Product orbitals use reciprocal-space gradients, not zty. - - this->orb11s.insert(std::make_pair(LA, Center2_Orb::Orb11(this->nA[LA], nB, psb_, this->MGT))); + (*products)[LA].set_orbital_info(nA_short.getLabel(), + nA_short.getType(), + LA, + 1, // N? + nA_short.getNr(), + nA_short.getRab(), + nA_short.getRadial(), + Numerical_Orbital_Lm::Psi_Type::Psi, + nA_tmp.data(), + nA_short.getNk(), + nA_short.getDk(), + nA_short.getDruniform(), + false, + true, + false); // Product orbitals use reciprocal-space gradients, not zty. + } + this->nA = std::move(products); + this->bind_product_orbitals(); +} +void Center2_Orb::Orb21::bind_product_orbitals() +{ + this->orb11s.clear(); + for (const auto& product: *this->nA) + { + this->orb11s.emplace(product.first, Orb11(product.second, this->nB, this->psb_, this->MGT)); } } @@ -80,10 +89,24 @@ void Center2_Orb::Orb21::init_radial_table(const std::set& radials) void Center2_Orb::Orb21::init_radial_table(const std::vector& orbitals, const bool use_gpu) { std::vector pairs; + // A product depends on its two orbitals on center A, not the orbital on + // center B. Share it within this build and retain ownership with each user. + std::map>> + products; for (Orb21* orbital: orbitals) { // Orbital transforms grow the shared Bessel cache, so prepare them serially. - orbital->prepare_radial_table(); + auto& product = products[&orbital->nA1][&orbital->nA2]; + if (!product) + { + orbital->prepare_radial_table(); + product = orbital->nA; + } + else + { + orbital->nA = product; + orbital->bind_product_orbitals(); + } for (auto& pair: orbital->orb11s) { pairs.push_back(&pair.second); diff --git a/source/source_lcao/center2orb_orb21.h b/source/source_lcao/center2orb_orb21.h index 330aebb4197..d2c6fd5eccb 100644 --- a/source/source_lcao/center2orb_orb21.h +++ b/source/source_lcao/center2orb_orb21.h @@ -13,6 +13,7 @@ #include "source_basis/module_ao/orb_gaunt_table.h" #include +#include #include #include @@ -44,6 +45,7 @@ class Center2_Orb::Orb21 private: void prepare_radial_table(); + void bind_product_orbitals(); const Numerical_Orbital_Lm& nA1; const Numerical_Orbital_Lm& nA2; @@ -52,7 +54,8 @@ class Center2_Orb::Orb21 const ModuleBase::Sph_Bessel_Recursive::D2* psb_ = nullptr; const ORB_gaunt_table& MGT; - std::map nA; + using ProductOrbitals = std::map; + std::shared_ptr nA; std::map orb11s; }; diff --git a/source/source_lcao/test/test_center2orb_orb11.cpp b/source/source_lcao/test/test_center2orb_orb11.cpp index 560e6de565f..7d509c9b48e 100644 --- a/source/source_lcao/test/test_center2orb_orb11.cpp +++ b/source/source_lcao/test/test_center2orb_orb11.cpp @@ -104,11 +104,58 @@ void compare_pairs(const bool use_gpu) } } } + +void compare_products(const bool use_gpu) +{ + std::vector basis; + for (int l = 0; l <= 3; ++l) + { + basis.push_back(make_orbital(l, l, 101)); + } + ORB_gaunt_table gaunt; + gaunt.init_Gaunt_CH(3); + gaunt.init_Gaunt(3); + ModuleBase::Sph_Bessel_Recursive::D2 bessel; + bessel.set_dx(0.001); + bessel.cal_jlx(7, 405, 65); + std::vector> reference; + std::vector> candidate; + std::vector tasks; + for (int i = 0; i < 4; ++i) + { + reference.emplace_back(new Center2_Orb::Orb21(basis[1], basis[2], basis[i], &bessel, gaunt)); + candidate.emplace_back(new Center2_Orb::Orb21(basis[1], basis[2], basis[i], &bessel, gaunt)); + reference.back()->init_radial_table(); + tasks.push_back(candidate.back().get()); + } + Center2_Orb::Orb21::init_radial_table(tasks, use_gpu); + const ModuleBase::Vector3 origin(0.0, 0.0, 0.0); + for (int i = 0; i < 4; ++i) + { + for (double distance: {0.0, 0.037, 0.73, 1.91, 3.87}) + { + const double x = 0.8 * distance; + const double y = 0.36 * distance; + const double z = 0.48 * distance; + const ModuleBase::Vector3 position(x, y, z); + const double expected = reference[i]->cal_overlap(origin, position, 1, 2, i); + const double actual = candidate[i]->cal_overlap(origin, position, 1, 2, i); + const double tolerance = 2e-11 * (1.0 + std::abs(expected)); + EXPECT_NEAR(actual, expected, tolerance); + const auto grad_expected = reference[i]->cal_grad_overlap(origin, position, 1, 2, i); + const auto grad_actual = candidate[i]->cal_grad_overlap(origin, position, 1, 2, i); + const double error = (grad_actual - grad_expected).norm(); + const double grad_tolerance = 2e-10 * (1.0 + grad_expected.norm()); + EXPECT_LT(error, grad_tolerance); + } + } +} } // namespace TEST(Center2Orb11, BatchedCpuMatchesIndividualTables) { compare_pairs(false); + compare_products(false); } TEST(Center2Orb21, ProductGradientMatchesFiniteDifference) @@ -157,5 +204,6 @@ TEST(Center2Orb11, BatchedGpuMatchesIndividualTables) GTEST_SKIP() << "A CUDA device is required"; } compare_pairs(true); + compare_products(true); } #endif From f399b63c861fdcd660cc341414110e8d2abe8570 Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Mon, 19 Oct 2026 21:00:49 +0800 Subject: [PATCH 17/24] test(ri): cover high angular momenta and bounded radial tiles --- .../test/test_center2orb_orb11.cpp | 98 +++++++++++++++++-- 1 file changed, 90 insertions(+), 8 deletions(-) diff --git a/source/source_lcao/test/test_center2orb_orb11.cpp b/source/source_lcao/test/test_center2orb_orb11.cpp index 7d509c9b48e..2b6583c3747 100644 --- a/source/source_lcao/test/test_center2orb_orb11.cpp +++ b/source/source_lcao/test/test_center2orb_orb11.cpp @@ -9,7 +9,7 @@ namespace { -Numerical_Orbital_Lm make_orbital(const int l, const int index, const int nr) +Numerical_Orbital_Lm make_orbital(const int l, const int index, const int nr, const int nk) { const double dr = 0.02; const double cutoff = (nr - 1) * dr; @@ -35,7 +35,7 @@ Numerical_Orbital_Lm make_orbital(const int l, const int index, const int nr) radial_ptr, Numerical_Orbital_Lm::Psi_Type::Psi, values_ptr, - 65, + nk, 0.1, 0.01, false, @@ -50,9 +50,9 @@ void compare_pairs(const bool use_gpu) std::vector basis; for (int l = 0; l <= 3; ++l) { - basis.push_back(make_orbital(l, l, 101)); + basis.push_back(make_orbital(l, l, 101, 65)); } - basis.push_back(make_orbital(2, 4, 81)); + basis.push_back(make_orbital(2, 4, 81, 65)); ORB_gaunt_table gaunt; gaunt.init_Gaunt_CH(3); gaunt.init_Gaunt(3); @@ -110,7 +110,7 @@ void compare_products(const bool use_gpu) std::vector basis; for (int l = 0; l <= 3; ++l) { - basis.push_back(make_orbital(l, l, 101)); + basis.push_back(make_orbital(l, l, 101, 65)); } ORB_gaunt_table gaunt; gaunt.init_Gaunt_CH(3); @@ -150,19 +150,99 @@ void compare_products(const bool use_gpu) } } } + +void compare_high_l_products(const bool use_gpu) +{ + // Heavy-element auxiliary bases include l=6; multiplying by an f orbital + // requires product angular momenta up to 9 and Bessel functions up to 13. + const auto left = make_orbital(6, 0, 101, 65); + const auto middle = make_orbital(3, 1, 101, 65); + const auto right = make_orbital(3, 2, 81, 65); + ORB_gaunt_table gaunt; + gaunt.init_Gaunt_CH(9); + gaunt.init_Gaunt(9); + ModuleBase::Sph_Bessel_Recursive::D2 bessel; + bessel.set_dx(0.001); + bessel.cal_jlx(13, 405, 65); + Center2_Orb::Orb21 reference(left, middle, right, &bessel, gaunt); + Center2_Orb::Orb21 candidate(left, middle, right, &bessel, gaunt); + reference.init_radial_table(); + const std::vector tasks{&candidate}; + Center2_Orb::Orb21::init_radial_table(tasks, use_gpu); + const ModuleBase::Vector3 origin(0.0, 0.0, 0.0); + for (double distance: {0.0, 0.037, 0.73, 1.91, 3.57, 3.8}) + { + const double x = 0.8 * distance; + const double y = 0.36 * distance; + const double z = 0.48 * distance; + const ModuleBase::Vector3 position(x, y, z); + for (int ma: {0, 6, 12}) + { + for (int mb: {0, 3, 6}) + { + for (int mc: {0, 3, 6}) + { + const double expected = reference.cal_overlap(origin, position, ma, mb, mc); + const double actual = candidate.cal_overlap(origin, position, ma, mb, mc); + const double tolerance = 2e-11 * (1.0 + std::abs(expected)); + EXPECT_NEAR(actual, expected, tolerance); + const auto grad_expected = reference.cal_grad_overlap(origin, position, ma, mb, mc); + const auto grad_actual = candidate.cal_grad_overlap(origin, position, ma, mb, mc); + const double error = (grad_actual - grad_expected).norm(); + const double grad_tolerance = 2e-10 * (1.0 + grad_expected.norm()); + EXPECT_LT(error, grad_tolerance); + } + } + } + } +} + +void compare_radial_tiles(const bool use_gpu) +{ + // The paired value/derivative Bessel matrix exceeds the 128 MiB workspace + // limit. Check interpolation across its radial tile boundary and tail. + const auto orbital = make_orbital(0, 0, 527, 4097); + ORB_gaunt_table gaunt; + gaunt.init_Gaunt_CH(0); + gaunt.init_Gaunt(0); + ModuleBase::Sph_Bessel_Recursive::D2 bessel; + bessel.set_dx(0.001); + bessel.cal_jlx(1, 2109, 4097); + Center2_Orb::Orb11 reference(orbital, orbital, &bessel, gaunt); + Center2_Orb::Orb11 candidate(orbital, orbital, &bessel, gaunt); + reference.init_radial_table(); + const std::vector tasks{&candidate}; + Center2_Orb::Orb11::init_radial_table(tasks, use_gpu); + const ModuleBase::Vector3 origin(0.0, 0.0, 0.0); + for (double distance: {0.0, 0.037, 1.91, 20.465, 20.475, 21.035, 21.1}) + { + const ModuleBase::Vector3 position(distance, 0.0, 0.0); + const double expected = reference.cal_overlap(origin, position, 0, 0); + const double actual = candidate.cal_overlap(origin, position, 0, 0); + const double tolerance = 2e-11 * (1.0 + std::abs(expected)); + EXPECT_NEAR(actual, expected, tolerance); + const auto grad_expected = reference.cal_grad_overlap(origin, position, 0, 0); + const auto grad_actual = candidate.cal_grad_overlap(origin, position, 0, 0); + const double error = (grad_actual - grad_expected).norm(); + const double grad_tolerance = 2e-10 * (1.0 + grad_expected.norm()); + EXPECT_LT(error, grad_tolerance); + } +} } // namespace TEST(Center2Orb11, BatchedCpuMatchesIndividualTables) { compare_pairs(false); compare_products(false); + compare_high_l_products(false); + compare_radial_tiles(false); } TEST(Center2Orb21, ProductGradientMatchesFiniteDifference) { - const auto left = make_orbital(2, 0, 101); - const auto middle = make_orbital(1, 1, 101); - const auto right = make_orbital(3, 2, 101); + const auto left = make_orbital(2, 0, 101, 65); + const auto middle = make_orbital(1, 1, 101, 65); + const auto right = make_orbital(3, 2, 101, 65); ORB_gaunt_table gaunt; gaunt.init_Gaunt_CH(3); gaunt.init_Gaunt(3); @@ -205,5 +285,7 @@ TEST(Center2Orb11, BatchedGpuMatchesIndividualTables) } compare_pairs(true); compare_products(true); + compare_high_l_products(true); + compare_radial_tiles(true); } #endif From 80b2059e2b20477b588cca56350bc78080e1efc1 Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Mon, 19 Oct 2026 21:06:39 +0800 Subject: [PATCH 18/24] perf(ri): reuse product orbitals in selective radial tables --- source/source_lcao/center2orb_orb21.cpp | 17 ++++++++++++++++- source/source_lcao/center2orb_orb21.h | 2 ++ source/source_lcao/module_ri/matrix_orbs21.cpp | 4 +++- .../source_lcao/test/test_center2orb_orb11.cpp | 18 ++++++++++++++++++ 4 files changed, 39 insertions(+), 2 deletions(-) diff --git a/source/source_lcao/center2orb_orb21.cpp b/source/source_lcao/center2orb_orb21.cpp index 75f2177458b..f980ca7ed99 100644 --- a/source/source_lcao/center2orb_orb21.cpp +++ b/source/source_lcao/center2orb_orb21.cpp @@ -87,6 +87,21 @@ void Center2_Orb::Orb21::init_radial_table(const std::set& radials) } void Center2_Orb::Orb21::init_radial_table(const std::vector& orbitals, const bool use_gpu) +{ + const auto pairs = prepare_radial_tables(orbitals); + Orb11::init_radial_table(pairs, use_gpu); +} + +void Center2_Orb::Orb21::init_radial_table(const std::vector& orbitals, const std::set& radials) +{ + const auto pairs = prepare_radial_tables(orbitals); + for (Orb11* pair: pairs) + { + pair->init_radial_table(radials); + } +} + +std::vector Center2_Orb::Orb21::prepare_radial_tables(const std::vector& orbitals) { std::vector pairs; // A product depends on its two orbitals on center A, not the orbital on @@ -112,7 +127,7 @@ void Center2_Orb::Orb21::init_radial_table(const std::vector& orbitals, pairs.push_back(&pair.second); } } - Orb11::init_radial_table(pairs, use_gpu); + return pairs; } double Center2_Orb::Orb21::cal_overlap(const ModuleBase::Vector3& RA, diff --git a/source/source_lcao/center2orb_orb21.h b/source/source_lcao/center2orb_orb21.h index d2c6fd5eccb..dab39312ee8 100644 --- a/source/source_lcao/center2orb_orb21.h +++ b/source/source_lcao/center2orb_orb21.h @@ -30,6 +30,7 @@ class Center2_Orb::Orb21 void init_radial_table(const std::set& radials); // unit: Bohr/MOT.dr static void init_radial_table(const std::vector& orbitals, const bool use_gpu); + static void init_radial_table(const std::vector& orbitals, const std::set& radials); double cal_overlap(const ModuleBase::Vector3& RA, const ModuleBase::Vector3& RB, // unit: Bohr @@ -46,6 +47,7 @@ class Center2_Orb::Orb21 private: void prepare_radial_table(); void bind_product_orbitals(); + static std::vector prepare_radial_tables(const std::vector& orbitals); const Numerical_Orbital_Lm& nA1; const Numerical_Orbital_Lm& nA2; diff --git a/source/source_lcao/module_ri/matrix_orbs21.cpp b/source/source_lcao/module_ri/matrix_orbs21.cpp index 31491a04d0c..0a90a64409c 100644 --- a/source/source_lcao/module_ri/matrix_orbs21.cpp +++ b/source/source_lcao/module_ri/matrix_orbs21.cpp @@ -190,6 +190,7 @@ void Matrix_Orbs21::init_radial_table(const std::map orbitals; for (auto& coC: *center2_orb21_sAB) { for (auto& coD: coC.second) @@ -202,13 +203,14 @@ void Matrix_Orbs21::init_radial_table(const std::map> reference; std::vector> candidate; + std::vector> selective; std::vector tasks; + std::vector selective_tasks; for (int i = 0; i < 4; ++i) { reference.emplace_back(new Center2_Orb::Orb21(basis[1], basis[2], basis[i], &bessel, gaunt)); candidate.emplace_back(new Center2_Orb::Orb21(basis[1], basis[2], basis[i], &bessel, gaunt)); + selective.emplace_back(new Center2_Orb::Orb21(basis[1], basis[2], basis[i], &bessel, gaunt)); reference.back()->init_radial_table(); tasks.push_back(candidate.back().get()); + selective_tasks.push_back(selective.back().get()); } Center2_Orb::Orb21::init_radial_table(tasks, use_gpu); + std::set radials{10000}; + for (size_t first: {0, 3, 72, 190, 386}) + { + for (size_t offset = 0; offset < 6; ++offset) + { + radials.insert(first + offset); + } + } + Center2_Orb::Orb21::init_radial_table(selective_tasks, radials); const ModuleBase::Vector3 origin(0.0, 0.0, 0.0); for (int i = 0; i < 4; ++i) { @@ -142,11 +155,16 @@ void compare_products(const bool use_gpu) const double actual = candidate[i]->cal_overlap(origin, position, 1, 2, i); const double tolerance = 2e-11 * (1.0 + std::abs(expected)); EXPECT_NEAR(actual, expected, tolerance); + const double selective_actual = selective[i]->cal_overlap(origin, position, 1, 2, i); + EXPECT_NEAR(selective_actual, expected, tolerance); const auto grad_expected = reference[i]->cal_grad_overlap(origin, position, 1, 2, i); const auto grad_actual = candidate[i]->cal_grad_overlap(origin, position, 1, 2, i); const double error = (grad_actual - grad_expected).norm(); const double grad_tolerance = 2e-10 * (1.0 + grad_expected.norm()); EXPECT_LT(error, grad_tolerance); + const auto selective_grad = selective[i]->cal_grad_overlap(origin, position, 1, 2, i); + const double selective_error = (selective_grad - grad_expected).norm(); + EXPECT_LT(selective_error, grad_tolerance); } } } From 1ce1fdc410f03f5c6324234dc2d51243956556c3 Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Mon, 19 Oct 2026 21:16:59 +0800 Subject: [PATCH 19/24] perf(ri): reuse radial interpolation across magnetic components --- source/source_lcao/center2orb_orb11.cpp | 57 +++++++++++-------------- 1 file changed, 24 insertions(+), 33 deletions(-) diff --git a/source/source_lcao/center2orb_orb11.cpp b/source/source_lcao/center2orb_orb11.cpp index 7568d8dae3e..5ea9f8e8592 100644 --- a/source/source_lcao/center2orb_orb11.cpp +++ b/source/source_lcao/center2orb_orb11.cpp @@ -363,9 +363,17 @@ double Center2_Orb::Orb11::cal_overlap(const ModuleBase::Vector3& RA, const int idx2 = this->MGT.get_lm_index(LB, mB); const double* Gaunt_Coefficients_ptr = &(this->MGT.Gaunt_Coefficients(idx1, idx2, 0)); + const int rmesh = Center2_Orb::get_rmesh(RcutA, RcutB, dr_); for (const auto& tb_r: this->Table_r) { const int LAB = tb_r.first; + const double distance_power = std::pow(distance, LAB); + const int phase = (LA - LB - LAB) / 2; + const double i_exp = std::pow(-1.0, phase); + const double Interp_Tlm + = (distance > tiny2) + ? ModuleBase::PolyInt::Polynomial_Interpolation(tb_r.second.data(), rmesh, dr_, distance) + : tb_r.second.at(0); for (int mAB = 0; mAB != 2 * LAB + 1; ++mAB) // const int mAB = mA + mB; @@ -382,17 +390,7 @@ double Center2_Orb::Orb11::cal_overlap(const ModuleBase::Vector3& RA, { continue; } - const double ylm_real = (distance > tiny2) ? ylm_solid / pow(distance, LAB) : ylm_solid; - - const double i_exp = std::pow(-1.0, (LA - LB - LAB) / 2); - - const double Interp_Tlm - = (distance > tiny2) - ? ModuleBase::PolyInt::Polynomial_Interpolation(tb_r.second.data(), - Center2_Orb::get_rmesh(RcutA, RcutB, dr_), - dr_, - distance) - : tb_r.second.at(0); + const double ylm_real = (distance > tiny2) ? ylm_solid / distance_power : ylm_solid; overlap += // pow(2*PI,1.5) @@ -441,9 +439,22 @@ ModuleBase::Vector3 Center2_Orb::Orb11::cal_grad_overlap( // caoyu add 2 ModuleBase::Vector3 grad_overlap(0.0, 0.0, 0.0); + const int rmesh = Center2_Orb::get_rmesh(RcutA, RcutB, dr_); for (const auto& tb_r: this->Table_r) { const int LAB = tb_r.first; + const double distance_power = std::pow(distance, LAB); + const int phase = (LA - LB - LAB) / 2; + const double i_exp = std::pow(-1.0, phase); + const double Interp_Tlm + = (distance > tiny2) + ? ModuleBase::PolyInt::Polynomial_Interpolation(tb_r.second.data(), rmesh, dr_, distance) + : tb_r.second.at(0); + const double grad_Interp_Tlm + = (distance > tiny2) + ? ModuleBase::PolyInt::Polynomial_Interpolation(this->Table_dr.at(LAB).data(), rmesh, dr_, distance) + - Interp_Tlm * LAB / distance + : 0.0; for (int mAB = 0; mAB != 2 * LAB + 1; ++mAB) // const int mAB = mA + mB; { @@ -455,30 +466,10 @@ ModuleBase::Vector3 Center2_Orb::Orb11::cal_grad_overlap( // caoyu add 2 } const double ylm_solid = rly[idx3]; - const double ylm_real = (distance > tiny2) ? ylm_solid / pow(distance, LAB) : ylm_solid; + const double ylm_real = (distance > tiny2) ? ylm_solid / distance_power : ylm_solid; const ModuleBase::Vector3 gylm_solid = grly[idx3]; - const ModuleBase::Vector3 gylm_real - = (distance > tiny2) ? gylm_solid / pow(distance, LAB) : gylm_solid; - - const double i_exp = std::pow(-1.0, (LA - LB - LAB) / 2); - - const double Interp_Tlm - = (distance > tiny2) - ? ModuleBase::PolyInt::Polynomial_Interpolation(tb_r.second.data(), - Center2_Orb::get_rmesh(RcutA, RcutB, dr_), - dr_, - distance) - : tb_r.second.at(0); - - const double grad_Interp_Tlm - = (distance > tiny2) - ? ModuleBase::PolyInt::Polynomial_Interpolation(this->Table_dr.at(LAB).data(), - Center2_Orb::get_rmesh(RcutA, RcutB, dr_), - dr_, - distance) // Interp(Table_dr) - - Interp_Tlm * LAB / distance - : 0.0; + const ModuleBase::Vector3 gylm_real = (distance > tiny2) ? gylm_solid / distance_power : gylm_solid; grad_overlap += i_exp // pow(2*PI,1.5) * Gaunt_real_A_B_AB From 6283aea20810c90b70056bb62ccc8948679e7070 Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Mon, 19 Oct 2026 21:26:54 +0800 Subject: [PATCH 20/24] perf(ri): share overlaps across three-center magnetic components --- source/source_lcao/center2orb_orb21.cpp | 68 +++++++++++++++++++ source/source_lcao/center2orb_orb21.h | 13 ++++ .../source_lcao/module_ri/matrix_orbs21.hpp | 53 +++++++++------ .../test/test_center2orb_orb11.cpp | 25 +++++++ 4 files changed, 137 insertions(+), 22 deletions(-) diff --git a/source/source_lcao/center2orb_orb21.cpp b/source/source_lcao/center2orb_orb21.cpp index f980ca7ed99..421acb79d41 100644 --- a/source/source_lcao/center2orb_orb21.cpp +++ b/source/source_lcao/center2orb_orb21.cpp @@ -197,3 +197,71 @@ ModuleBase::Vector3 Center2_Orb::Orb21::cal_grad_overlap(const ModuleBas return grad_overlap; } + +template +std::vector Center2_Orb::Orb21::cal_magnetic_components( + const ModuleBase::Vector3& RA, + const ModuleBase::Vector3& RB, + T (Orb11::*evaluate)(const ModuleBase::Vector3&, const ModuleBase::Vector3&, const int&, const int&) + const) const +{ + const int LA1 = this->nA1.getL(); + const int LA2 = this->nA2.getL(); + const int count_A1 = 2 * LA1 + 1; + const int count_A2 = 2 * LA2 + 1; + const int count_B = 2 * this->nB.getL() + 1; + const int count = count_A1 * count_A2 * count_B; + std::vector result(count); + for (const auto& orbital: this->orb11s) + { + const int LA = orbital.first; + const int count_A = 2 * LA + 1; + const int pair_count = count_A * count_B; + std::vector overlaps(pair_count); + // These two-center overlaps are shared by all ways to couple the + // magnetic components of the two original orbitals on center A. + for (int ma = 0; ma < count_A; ++ma) + { + for (int mb = 0; mb < count_B; ++mb) + { + overlaps[ma * count_B + mb] = (orbital.second.*evaluate)(RA, RB, ma, mb); + } + } + for (int ma1 = 0; ma1 < count_A1; ++ma1) + { + const int idx1 = this->MGT.get_lm_index(LA1, ma1); + for (int ma2 = 0; ma2 < count_A2; ++ma2) + { + const int idx2 = this->MGT.get_lm_index(LA2, ma2); + const int offset = (ma1 * count_A2 + ma2) * count_B; + for (int ma = 0; ma < count_A; ++ma) + { + const int idx = this->MGT.get_lm_index(LA, ma); + const double coefficient = this->MGT.Gaunt_Coefficients(idx1, idx2, idx); + if (coefficient == 0.0) + { + continue; + } + for (int mb = 0; mb < count_B; ++mb) + { + result[offset + mb] += coefficient * overlaps[ma * count_B + mb]; + } + } + } + } + } + return result; +} + +std::vector Center2_Orb::Orb21::cal_overlap(const ModuleBase::Vector3& RA, + const ModuleBase::Vector3& RB) const +{ + return this->cal_magnetic_components(RA, RB, &Orb11::cal_overlap); +} + +std::vector> Center2_Orb::Orb21::cal_grad_overlap( + const ModuleBase::Vector3& RA, + const ModuleBase::Vector3& RB) const +{ + return this->cal_magnetic_components(RA, RB, &Orb11::cal_grad_overlap); +} diff --git a/source/source_lcao/center2orb_orb21.h b/source/source_lcao/center2orb_orb21.h index dab39312ee8..f5cb77d1658 100644 --- a/source/source_lcao/center2orb_orb21.h +++ b/source/source_lcao/center2orb_orb21.h @@ -44,11 +44,24 @@ class Center2_Orb::Orb21 const int& mA2, const int& mB) const; + /// All magnetic components, ordered as [mA1][mA2][mB]. + std::vector cal_overlap(const ModuleBase::Vector3& RA, const ModuleBase::Vector3& RB) const; + std::vector> cal_grad_overlap(const ModuleBase::Vector3& RA, + const ModuleBase::Vector3& RB) const; + private: void prepare_radial_table(); void bind_product_orbitals(); static std::vector prepare_radial_tables(const std::vector& orbitals); + template + std::vector cal_magnetic_components(const ModuleBase::Vector3& RA, + const ModuleBase::Vector3& RB, + T (Orb11::*evaluate)(const ModuleBase::Vector3&, + const ModuleBase::Vector3&, + const int&, + const int&) const) const; + const Numerical_Orbital_Lm& nA1; const Numerical_Orbital_Lm& nA2; const Numerical_Orbital_Lm& nB; diff --git a/source/source_lcao/module_ri/matrix_orbs21.hpp b/source/source_lcao/module_ri/matrix_orbs21.hpp index b6f464d7440..47d3bfa0f51 100644 --- a/source/source_lcao/module_ri/matrix_orbs21.hpp +++ b/source/source_lcao/module_ri/matrix_orbs21.hpp @@ -36,31 +36,35 @@ RI::Tensor Matrix_Orbs21::cal_overlap_matrix( default: throw std::invalid_argument(std::string(__FILE__)+" line "+std::to_string(__LINE__)); } + const ModuleBase::Vector3 position_A = tauA * lat0; + const ModuleBase::Vector3 position_B = tauB * lat0; for( const auto &co3 : center2_orb21_s.at(TA).at(TB) ) { const int LA1 = co3.first; for( const auto &co4 : co3.second ) { const size_t NA1 = co4.first; - for( size_t MA1=0; MA1!=2*LA1+1; ++MA1 ) + for (const auto& co5: co4.second) { - for( const auto &co5 : co4.second ) + const int LA2 = co5.first; + for (const auto& co6: co5.second) { - const int LA2 = co5.first; - for( const auto &co6 : co5.second ) + const size_t NA2 = co6.first; + for (const auto& co7: co6.second) { - const size_t NA2 = co6.first; - for( size_t MA2=0; MA2!=2*LA2+1; ++MA2 ) + const int LB = co7.first; + for (const auto& co8: co7.second) { - for( const auto &co7 : co6.second ) + const size_t NB = co8.first; + const auto overlaps = co8.second.cal_overlap(position_A, position_B); + for (size_t MA1 = 0; MA1 != 2 * LA1 + 1; ++MA1) { - const int LB = co7.first; - for( const auto &co8 : co7.second ) + for (size_t MA2 = 0; MA2 != 2 * LA2 + 1; ++MA2) { - const size_t NB = co8.first; for( size_t MB=0; MB!=2*LB+1; ++MB ) { - const Tdata overlap = co8.second.cal_overlap( tauA*lat0, tauB*lat0, MA1, MA2, MB ); + const size_t offset = (MA1 * (2 * LA2 + 1) + MA2) * (2 * LB + 1) + MB; + const Tdata overlap = overlaps[offset]; const size_t iA1 = index_A1[TA][LA1][NA1][MA1]; const size_t iA2 = index_A2[TA][LA2][NA2][MA2]; const size_t iB = index_B[TB][LB][NB][MB]; @@ -116,31 +120,36 @@ std::array,3> Matrix_Orbs21::cal_grad_overlap_matrix( } } + const ModuleBase::Vector3 position_A = tauA * lat0; + const ModuleBase::Vector3 position_B = tauB * lat0; for( const auto &co3 : center2_orb21_s.at(TA).at(TB) ) { const int LA1 = co3.first; for( const auto &co4 : co3.second ) { const size_t NA1 = co4.first; - for( size_t MA1=0; MA1!=2*LA1+1; ++MA1 ) + for (const auto& co5: co4.second) { - for( const auto &co5 : co4.second ) + const int LA2 = co5.first; + for (const auto& co6: co5.second) { - const int LA2 = co5.first; - for( const auto &co6 : co5.second ) + const size_t NA2 = co6.first; + for (const auto& co7: co6.second) { - const size_t NA2 = co6.first; - for( size_t MA2=0; MA2!=2*LA2+1; ++MA2 ) + const int LB = co7.first; + for (const auto& co8: co7.second) { - for( const auto &co7 : co6.second ) + const size_t NB = co8.first; + const auto overlaps = co8.second.cal_grad_overlap(position_A, position_B); + for (size_t MA1 = 0; MA1 != 2 * LA1 + 1; ++MA1) { - const int LB = co7.first; - for( const auto &co8 : co7.second ) + for (size_t MA2 = 0; MA2 != 2 * LA2 + 1; ++MA2) { - const size_t NB = co8.first; for( size_t MB=0; MB!=2*LB+1; ++MB ) { - const std::array grad_overlap = RI_Util::Vector3_to_array3(co8.second.cal_grad_overlap( tauA*lat0, tauB*lat0, MA1, MA2, MB )); + const size_t offset = (MA1 * (2 * LA2 + 1) + MA2) * (2 * LB + 1) + MB; + const std::array grad_overlap + = RI_Util::Vector3_to_array3(overlaps[offset]); const size_t iA1 = index_A1[TA][LA1][NA1][MA1]; const size_t iA2 = index_A2[TA][LA2][NA2][MA2]; const size_t iB = index_B[TB][LB][NB][MB]; diff --git a/source/source_lcao/test/test_center2orb_orb11.cpp b/source/source_lcao/test/test_center2orb_orb11.cpp index a3f2116851d..42fd6c8b835 100644 --- a/source/source_lcao/test/test_center2orb_orb11.cpp +++ b/source/source_lcao/test/test_center2orb_orb11.cpp @@ -165,6 +165,25 @@ void compare_products(const bool use_gpu) const auto selective_grad = selective[i]->cal_grad_overlap(origin, position, 1, 2, i); const double selective_error = (selective_grad - grad_expected).norm(); EXPECT_LT(selective_error, grad_tolerance); + const auto overlaps = candidate[i]->cal_overlap(origin, position); + const auto gradients = candidate[i]->cal_grad_overlap(origin, position); + for (int ma = 0; ma < 3; ++ma) + { + for (int mb = 0; mb < 5; ++mb) + { + for (int mc = 0; mc < 2 * i + 1; ++mc) + { + const int offset = (ma * 5 + mb) * (2 * i + 1) + mc; + const double value = reference[i]->cal_overlap(origin, position, ma, mb, mc); + const double value_tolerance = 2e-11 * (1.0 + std::abs(value)); + EXPECT_NEAR(overlaps[offset], value, value_tolerance); + const auto gradient = reference[i]->cal_grad_overlap(origin, position, ma, mb, mc); + const double component_error = (gradients[offset] - gradient).norm(); + const double component_tolerance = 2e-10 * (1.0 + gradient.norm()); + EXPECT_LT(component_error, component_tolerance); + } + } + } } } } @@ -194,6 +213,8 @@ void compare_high_l_products(const bool use_gpu) const double y = 0.36 * distance; const double z = 0.48 * distance; const ModuleBase::Vector3 position(x, y, z); + const auto overlaps = candidate.cal_overlap(origin, position); + const auto gradients = candidate.cal_grad_overlap(origin, position); for (int ma: {0, 6, 12}) { for (int mb: {0, 3, 6}) @@ -204,11 +225,15 @@ void compare_high_l_products(const bool use_gpu) const double actual = candidate.cal_overlap(origin, position, ma, mb, mc); const double tolerance = 2e-11 * (1.0 + std::abs(expected)); EXPECT_NEAR(actual, expected, tolerance); + const int offset = (ma * 7 + mb) * 7 + mc; + EXPECT_NEAR(overlaps[offset], expected, tolerance); const auto grad_expected = reference.cal_grad_overlap(origin, position, ma, mb, mc); const auto grad_actual = candidate.cal_grad_overlap(origin, position, ma, mb, mc); const double error = (grad_actual - grad_expected).norm(); const double grad_tolerance = 2e-10 * (1.0 + grad_expected.norm()); EXPECT_LT(error, grad_tolerance); + const double component_error = (gradients[offset] - grad_expected).norm(); + EXPECT_LT(component_error, grad_tolerance); } } } From 26d01a67dffb858af5db20ef329b55f3e657b7da Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Mon, 19 Oct 2026 21:44:01 +0800 Subject: [PATCH 21/24] fix(ri): keep radial tiles within the output workspace limit --- source/source_lcao/center2orb_orb11.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/source/source_lcao/center2orb_orb11.cpp b/source/source_lcao/center2orb_orb11.cpp index 5ea9f8e8592..bbb139c0224 100644 --- a/source/source_lcao/center2orb_orb11.cpp +++ b/source/source_lcao/center2orb_orb11.cpp @@ -176,7 +176,8 @@ void Center2_Orb::Orb11::init_radial_table(const std::vector& orbitals, { throw std::invalid_argument("Unsupported radial quadrature grid for batched tables"); } - const int radial_block = std::min(rmesh, table_limit / (2 * kmesh)); + const std::size_t radial_limit = std::min(table_limit / (2 * kmesh), output_limit / 2); + const int radial_block = std::min(rmesh, radial_limit); const std::size_t batch_limit = std::min(input_limit / kmesh, output_limit / (2 * radial_block)); const int batch_size = std::min(tasks.size(), std::min(128, batch_limit)); const std::size_t input_size = static_cast(batch_size) * kmesh; From 5156a5d4ebf3190108d0d87f5a817c62ee9dfe96 Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Mon, 19 Oct 2026 22:33:26 +0800 Subject: [PATCH 22/24] style(ri): use the function name for radial table timing --- source/source_lcao/center2orb_orb11.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/source_lcao/center2orb_orb11.cpp b/source/source_lcao/center2orb_orb11.cpp index bbb139c0224..d53af9a7b45 100644 --- a/source/source_lcao/center2orb_orb11.cpp +++ b/source/source_lcao/center2orb_orb11.cpp @@ -123,7 +123,7 @@ double radial_product(const Numerical_Orbital_Lm& left, const Numerical_Orbital_ void Center2_Orb::Orb11::init_radial_table(const std::vector& orbitals, const bool use_gpu) { - ModuleBase::timer::start("Center2_Orb", "init_radial_table_batch"); + ModuleBase::timer::start("Center2_Orb", "init_radial_table"); #ifndef __CUDA if (use_gpu) { @@ -328,7 +328,7 @@ void Center2_Orb::Orb11::init_radial_table(const std::vector& orbitals, } } } - ModuleBase::timer::end("Center2_Orb", "init_radial_table_batch"); + ModuleBase::timer::end("Center2_Orb", "init_radial_table"); } double Center2_Orb::Orb11::cal_overlap(const ModuleBase::Vector3& RA, From 3fa1f972b60b95767c10f2a14b30fe237f96dedb Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Tue, 20 Oct 2026 10:33:14 +0800 Subject: [PATCH 23/24] refactor(input): rename exx_device to device_exx --- docs/advanced/install.md | 6 +++--- .../module_parameter/input_parameter.h | 2 +- .../module_parameter/read_inp_exx_dftu.cpp | 18 +++++++++--------- .../test_serial/read_input_item_test.cpp | 18 +++++++++--------- source/source_lcao/setup_exx.cpp | 2 +- 5 files changed, 23 insertions(+), 23 deletions(-) diff --git a/docs/advanced/install.md b/docs/advanced/install.md index a4aafa470c6..9c3959932f3 100644 --- a/docs/advanced/install.md +++ b/docs/advanced/install.md @@ -84,10 +84,10 @@ The new EXX implementation directly depends on [LibRI](https://github.com/abacus Experimental LCAO GPU contractions additionally require `ENABLE_LIBRI_GPU=ON` (default OFF), `USE_CUDA=ON`, and a CUDA build of DDLA. Like other dependent options, GPU RI is disabled when its LibRI/CUDA prerequisites are disabled. -INPUT `exx_device` defaults to `cpu`, retaining CPU contractions even with -`device gpu`. To use GPU contractions, set both `device gpu` and `exx_device gpu`. +INPUT `device_exx` defaults to `cpu`, retaining CPU contractions even with +`device gpu`. To use GPU contractions, set both `device gpu` and `device_exx gpu`. The latter requires GPU RI to be compiled in. Real-space EXX symmetry reduction -is unsupported by this backend: with `symmetry 1`, use `exx_device cpu` or +is unsupported by this backend: with `symmetry 1`, use `device_exx cpu` or `exx_symmetry_realspace 0`. The selector applies to the LCAO Kohn-Sham exchange engine. Independent LR and diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 0e45d797dfd..d2f4d013958 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -634,7 +634,7 @@ struct Input_para ///< the evaluation of Fock exchange double exx_mixing_beta = 1.0; ///< mixing_beta for outer-loop when exx_separate_loop=1 std::string exx_real_number = "default"; ///< exx calculated in real or complex - std::string exx_device = "cpu"; ///< computing device for LCAO EXX tensor contractions + std::string device_exx = "cpu"; ///< computing device for LCAO EXX tensor contractions double exx_pca_threshold = 0.0001; ///< threshold to screen on-site ABFs in exx double exx_c_threshold = 0.0001; ///< threshold to screen C matrix in exx double exx_v_threshold = 0.1; ///< threshold to screen C matrix in exx diff --git a/source/source_io/module_parameter/read_inp_exx_dftu.cpp b/source/source_io/module_parameter/read_inp_exx_dftu.cpp index c3ce4d961d4..f6f546dc440 100644 --- a/source/source_io/module_parameter/read_inp_exx_dftu.cpp +++ b/source/source_io/module_parameter/read_inp_exx_dftu.cpp @@ -238,7 +238,7 @@ void ReadInput::item_exx() this->add_item(item); } { - Input_Item item("exx_device"); + Input_Item item("device_exx"); item.annotation = "computing device for LCAO EXX tensor contractions"; item.category = "Exact Exchange (LCAO)"; item.type = "String"; @@ -246,7 +246,7 @@ void ReadInput::item_exx() * cpu: retain CPU contractions, including when device=gpu. * gpu: use the experimental GPU backend; requires ENABLE_LIBRI_GPU=ON and device=gpu. GPU contractions do not support real-space EXX symmetry reduction. With symmetry=1, -use exx_device=cpu or disable exx_symmetry_realspace. +use device_exx=cpu or disable exx_symmetry_realspace. Independent LR and RDMFT exchange engines retain CPU contractions. LR inherits this selection only when it reuses the Kohn-Sham exchange object. For GPU contractions, all MPI ranks on a node must see the same allocated GPU @@ -254,26 +254,26 @@ list, with at least as many visible GPUs as local ranks; per-rank single-GPU masking is unsupported.)"; item.default_value = "cpu"; item.set_availability("basis_type==lcao"); - read_sync_string(input.exx_device); + read_sync_string(input.device_exx); item.check_value = [](const Input_Item& item, const Parameter& para) { - if (para.input.exx_device != "cpu" && para.input.exx_device != "gpu") + if (para.input.device_exx != "cpu" && para.input.device_exx != "gpu") { - ModuleBase::WARNING_QUIT("ReadInput", "exx_device must be cpu or gpu"); + ModuleBase::WARNING_QUIT("ReadInput", "device_exx must be cpu or gpu"); } - if (para.input.exx_device == "gpu") + if (para.input.device_exx == "gpu") { if (para.input.basis_type != "lcao" || para.input.device != "gpu") { - ModuleBase::WARNING_QUIT("ReadInput", "exx_device=gpu requires basis_type=lcao and device=gpu"); + ModuleBase::WARNING_QUIT("ReadInput", "device_exx=gpu requires basis_type=lcao and device=gpu"); } #ifndef __GPU_RI - ModuleBase::WARNING_QUIT("ReadInput", "exx_device=gpu requires a build with ENABLE_LIBRI_GPU=ON"); + ModuleBase::WARNING_QUIT("ReadInput", "device_exx=gpu requires a build with ENABLE_LIBRI_GPU=ON"); #endif if (para.input.symmetry == "1" && para.input.exx_symmetry_realspace) { ModuleBase::WARNING_QUIT("ReadInput", "LibRI GPU contractions do not support real-space EXX symmetry; " - "set exx_device=cpu or exx_symmetry_realspace=0"); + "set device_exx=cpu or exx_symmetry_realspace=0"); } } }; diff --git a/source/source_io/test_serial/read_input_item_test.cpp b/source/source_io/test_serial/read_input_item_test.cpp index 9d555aaf0bd..f91202f9bfe 100644 --- a/source/source_io/test_serial/read_input_item_test.cpp +++ b/source/source_io/test_serial/read_input_item_test.cpp @@ -134,13 +134,13 @@ TEST_F(InputTest, PccSelectionAndUnsupportedCombinations) EXPECT_EXIT(it->second.check_value(it->second, param), testing::ExitedWithCode(1), ""); } -TEST_F(InputTest, ExxDevice) +TEST_F(InputTest, DeviceExx) { ModuleIO::ReadInput readinput(0); Parameter param; - auto it = find_label("exx_device", readinput.input_lists); + auto it = find_label("device_exx", readinput.input_lists); ASSERT_NE(it, readinput.input_lists.end()); - EXPECT_EQ(param.inp.exx_device, "cpu"); + EXPECT_EQ(param.inp.device_exx, "cpu"); EXPECT_EQ(it->second.default_value, "cpu"); auto& input = TestParameters::input(param); @@ -161,16 +161,16 @@ TEST_F(InputTest, ExxDevice) it->second.str_values = {"ddla"}; it->second.read_value(it->second, param); - expect_rejected("exx_device must be cpu or gpu"); + expect_rejected("device_exx must be cpu or gpu"); it->second.str_values = {"gpu"}; it->second.read_value(it->second, param); - EXPECT_EQ(param.inp.exx_device, "gpu"); + EXPECT_EQ(param.inp.device_exx, "gpu"); input.basis_type = "pw"; - expect_rejected("exx_device=gpu requires basis_type=lcao and device=gpu"); + expect_rejected("device_exx=gpu requires basis_type=lcao and device=gpu"); input.basis_type = "lcao"; input.device = "cpu"; - expect_rejected("exx_device=gpu requires basis_type=lcao and device=gpu"); + expect_rejected("device_exx=gpu requires basis_type=lcao and device=gpu"); input.device = "gpu"; #ifdef __GPU_RI @@ -181,12 +181,12 @@ TEST_F(InputTest, ExxDevice) input.symmetry = "0"; it->second.check_value(it->second, param); #else - expect_rejected("exx_device=gpu requires a build with ENABLE_LIBRI_GPU=ON"); + expect_rejected("device_exx=gpu requires a build with ENABLE_LIBRI_GPU=ON"); #endif it->second.str_values = {"cpu"}; it->second.read_value(it->second, param); - EXPECT_EQ(param.inp.exx_device, "cpu"); + EXPECT_EQ(param.inp.device_exx, "cpu"); input.symmetry = "1"; input.exx_symmetry_realspace = true; it->second.check_value(it->second, param); diff --git a/source/source_lcao/setup_exx.cpp b/source/source_lcao/setup_exx.cpp index 219a22356f0..8be35a6149b 100644 --- a/source/source_lcao/setup_exx.cpp +++ b/source/source_lcao/setup_exx.cpp @@ -50,7 +50,7 @@ void Exx_NAO::before_runner( Exx_Info& exx_info) { #ifdef __EXX - const bool use_gpu = inp.exx_device == "gpu"; + const bool use_gpu = inp.device_exx == "gpu"; if (this->exd) { this->exd->set_gpu(use_gpu); From 810d4909c7b5392fe6830a0a663044381bd832f1 Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Tue, 20 Oct 2026 11:01:44 +0800 Subject: [PATCH 24/24] perf(ri): reuse angular terms across magnetic components Evaluate displacement-dependent harmonics and radial interpolation once for each Orb11 magnetic block, and reuse the vector path in Orb21 assembly. Extend scalar-reference and finite-difference coverage without changing numerical thresholds. --- source/source_lcao/center2orb_orb11.cpp | 138 ++++++++++++++++++ source/source_lcao/center2orb_orb11.h | 9 ++ source/source_lcao/center2orb_orb21.cpp | 18 +-- source/source_lcao/center2orb_orb21.h | 11 +- .../test/test_center2orb_orb11.cpp | 44 +++--- 5 files changed, 183 insertions(+), 37 deletions(-) diff --git a/source/source_lcao/center2orb_orb11.cpp b/source/source_lcao/center2orb_orb11.cpp index d53af9a7b45..58bcf2f1e06 100644 --- a/source/source_lcao/center2orb_orb11.cpp +++ b/source/source_lcao/center2orb_orb11.cpp @@ -480,3 +480,141 @@ ModuleBase::Vector3 Center2_Orb::Orb11::cal_grad_overlap( // caoyu add 2 return grad_overlap; } + +template +std::vector Center2_Orb::Orb11::contract_angular(const std::vector& components) const +{ + const int LA = this->nA.getL(); + const int LB = this->nB.getL(); + const int count_A = 2 * LA + 1; + const int count_B = 2 * LB + 1; + const int count = count_A * count_B; + std::vector result(count); + for (int ma = 0; ma < count_A; ++ma) + { + const int idx_A = this->MGT.get_lm_index(LA, ma); + for (int mb = 0; mb < count_B; ++mb) + { + const int idx_B = this->MGT.get_lm_index(LB, mb); + const double* coefficients = &this->MGT.Gaunt_Coefficients(idx_A, idx_B, 0); + T& value = result[ma * count_B + mb]; + // Keep the scalar path's angular-momentum accumulation order. + for (const std::pair>& table: this->Table_r) + { + const int l = table.first; + for (int m = 0; m <= 2 * l; ++m) + { + const int idx = this->MGT.get_lm_index(l, m); + if (coefficients[idx] != 0.0) + { + value += coefficients[idx] * components[idx]; + } + } + } + } + } + return result; +} + +std::vector Center2_Orb::Orb11::cal_overlap(const ModuleBase::Vector3& RA, + const ModuleBase::Vector3& RB) const +{ + const double tiny1 = 1e-12; + const double tiny2 = 1e-10; + const ModuleBase::Vector3 delta_R = RB - RA; + const double distance_true = delta_R.norm(); + const double distance = distance_true >= tiny1 ? distance_true : distance_true + tiny1; + const int LA = this->nA.getL(); + const int LB = this->nB.getL(); + const double RcutA = this->nA.getRcut(); + const double RcutB = this->nB.getRcut(); + if (distance > RcutA + RcutB) + { + const int count = (2 * LA + 1) * (2 * LB + 1); + return std::vector(count); + } + + const int lmax = LA + LB; + std::vector rly; + ModuleBase::Ylm::rl_sph_harm(lmax, delta_R.x, delta_R.y, delta_R.z, rly); + const int angular_count = (lmax + 1) * (lmax + 1); + std::vector components(angular_count); + const int rmesh = Center2_Orb::get_rmesh(RcutA, RcutB, this->dr_); + for (const std::pair>& table: this->Table_r) + { + const int l = table.first; + const double distance_power = std::pow(distance, l); + const int phase = (LA - LB - l) / 2; + const double i_exp = std::pow(-1.0, phase); + const double radial + = distance > tiny2 + ? ModuleBase::PolyInt::Polynomial_Interpolation(table.second.data(), rmesh, dr_, distance) + : table.second.at(0); + for (int m = 0; m <= 2 * l; ++m) + { + const int idx = this->MGT.get_lm_index(l, m); + const double angular = distance > tiny2 ? rly[idx] / distance_power : rly[idx]; + components[idx] = i_exp * radial * angular; + } + } + return this->contract_angular(components); +} + +std::vector> Center2_Orb::Orb11::cal_grad_overlap( + const ModuleBase::Vector3& RA, + const ModuleBase::Vector3& RB) const +{ + const double tiny1 = 1e-12; + const double tiny2 = 1e-10; + const ModuleBase::Vector3 delta_R = RB - RA; + const double distance_true = delta_R.norm(); + const double distance = distance_true >= tiny1 ? distance_true : distance_true + tiny1; + const int LA = this->nA.getL(); + const int LB = this->nB.getL(); + const double RcutA = this->nA.getRcut(); + const double RcutB = this->nB.getRcut(); + if (distance > RcutA + RcutB) + { + const int count = (2 * LA + 1) * (2 * LB + 1); + return std::vector>(count); + } + + const int lmax = LA + LB; + const int angular_count = (lmax + 1) * (lmax + 1); + const int gradient_count = 3 * angular_count; + std::vector rly(angular_count); + std::vector grly(gradient_count); + ModuleBase::Ylm::grad_rl_sph_harm(lmax, delta_R.x, delta_R.y, delta_R.z, rly.data(), grly.data()); + std::vector> components(angular_count); + const int rmesh = Center2_Orb::get_rmesh(RcutA, RcutB, this->dr_); + for (const std::pair>& table: this->Table_r) + { + const int l = table.first; + const double distance_power = std::pow(distance, l); + const int phase = (LA - LB - l) / 2; + const double i_exp = std::pow(-1.0, phase); + const double radial + = distance > tiny2 + ? ModuleBase::PolyInt::Polynomial_Interpolation(table.second.data(), rmesh, dr_, distance) + : table.second.at(0); + const double derivative + = distance > tiny2 + ? ModuleBase::PolyInt::Polynomial_Interpolation(this->Table_dr.at(l).data(), rmesh, dr_, distance) + - radial * l / distance + : 0.0; + for (int m = 0; m <= 2 * l; ++m) + { + const int idx = this->MGT.get_lm_index(l, m); + const double angular = distance > tiny2 ? rly[idx] / distance_power : rly[idx]; + const int offset = 3 * idx; + const double gx = grly[offset]; + const double gy = grly[offset + 1]; + const double gz = grly[offset + 2]; + const ModuleBase::Vector3 solid_gradient(gx, gy, gz); + const ModuleBase::Vector3 angular_gradient + = distance > tiny2 ? solid_gradient / distance_power : solid_gradient; + components[idx] = i_exp * (radial * angular_gradient + derivative * angular * delta_R / distance); + } + } + return this->contract_angular(components); +} diff --git a/source/source_lcao/center2orb_orb11.h b/source/source_lcao/center2orb_orb11.h index 2b718eee529..6cb2af3828c 100644 --- a/source/source_lcao/center2orb_orb11.h +++ b/source/source_lcao/center2orb_orb11.h @@ -41,7 +41,16 @@ class Center2_Orb::Orb11 const int& mA, const int& mB) const; + /// All magnetic components, ordered as [mA][mB]. + std::vector cal_overlap(const ModuleBase::Vector3& RA, + const ModuleBase::Vector3& RB) const; + std::vector> cal_grad_overlap(const ModuleBase::Vector3& RA, + const ModuleBase::Vector3& RB) const; + private: + template + std::vector contract_angular(const std::vector& components) const; + const Numerical_Orbital_Lm& nA; const Numerical_Orbital_Lm& nB; diff --git a/source/source_lcao/center2orb_orb21.cpp b/source/source_lcao/center2orb_orb21.cpp index 421acb79d41..3f1242f13e6 100644 --- a/source/source_lcao/center2orb_orb21.cpp +++ b/source/source_lcao/center2orb_orb21.cpp @@ -202,8 +202,8 @@ template std::vector Center2_Orb::Orb21::cal_magnetic_components( const ModuleBase::Vector3& RA, const ModuleBase::Vector3& RB, - T (Orb11::*evaluate)(const ModuleBase::Vector3&, const ModuleBase::Vector3&, const int&, const int&) - const) const + std::vector (Orb11::*evaluate)(const ModuleBase::Vector3&, + const ModuleBase::Vector3&) const) const { const int LA1 = this->nA1.getL(); const int LA2 = this->nA2.getL(); @@ -216,17 +216,9 @@ std::vector Center2_Orb::Orb21::cal_magnetic_components( { const int LA = orbital.first; const int count_A = 2 * LA + 1; - const int pair_count = count_A * count_B; - std::vector overlaps(pair_count); // These two-center overlaps are shared by all ways to couple the // magnetic components of the two original orbitals on center A. - for (int ma = 0; ma < count_A; ++ma) - { - for (int mb = 0; mb < count_B; ++mb) - { - overlaps[ma * count_B + mb] = (orbital.second.*evaluate)(RA, RB, ma, mb); - } - } + const std::vector overlaps = (orbital.second.*evaluate)(RA, RB); for (int ma1 = 0; ma1 < count_A1; ++ma1) { const int idx1 = this->MGT.get_lm_index(LA1, ma1); @@ -256,12 +248,12 @@ std::vector Center2_Orb::Orb21::cal_magnetic_components( std::vector Center2_Orb::Orb21::cal_overlap(const ModuleBase::Vector3& RA, const ModuleBase::Vector3& RB) const { - return this->cal_magnetic_components(RA, RB, &Orb11::cal_overlap); + return this->cal_magnetic_components(RA, RB, &Orb11::cal_overlap); } std::vector> Center2_Orb::Orb21::cal_grad_overlap( const ModuleBase::Vector3& RA, const ModuleBase::Vector3& RB) const { - return this->cal_magnetic_components(RA, RB, &Orb11::cal_grad_overlap); + return this->cal_magnetic_components>(RA, RB, &Orb11::cal_grad_overlap); } diff --git a/source/source_lcao/center2orb_orb21.h b/source/source_lcao/center2orb_orb21.h index f5cb77d1658..8ca689e4e07 100644 --- a/source/source_lcao/center2orb_orb21.h +++ b/source/source_lcao/center2orb_orb21.h @@ -55,12 +55,11 @@ class Center2_Orb::Orb21 static std::vector prepare_radial_tables(const std::vector& orbitals); template - std::vector cal_magnetic_components(const ModuleBase::Vector3& RA, - const ModuleBase::Vector3& RB, - T (Orb11::*evaluate)(const ModuleBase::Vector3&, - const ModuleBase::Vector3&, - const int&, - const int&) const) const; + std::vector cal_magnetic_components( + const ModuleBase::Vector3& RA, + const ModuleBase::Vector3& RB, + std::vector (Orb11::*evaluate)(const ModuleBase::Vector3&, + const ModuleBase::Vector3&) const) const; const Numerical_Orbital_Lm& nA1; const Numerical_Orbital_Lm& nA2; diff --git a/source/source_lcao/test/test_center2orb_orb11.cpp b/source/source_lcao/test/test_center2orb_orb11.cpp index 42fd6c8b835..3c475badd7e 100644 --- a/source/source_lcao/test/test_center2orb_orb11.cpp +++ b/source/source_lcao/test/test_center2orb_orb11.cpp @@ -80,25 +80,34 @@ void compare_pairs(const bool use_gpu) const int right = i < 130 ? 2 : (i + 2) % 5; const int last_left = 2 * basis[left].getL(); const int last_right = 2 * basis[right].getL(); - for (double distance: {0.0, 0.037, 0.73, 1.91, 3.87, 4.2}) + for (double distance: {0.0, 1e-12, 5e-11, 2e-10, 0.037, 0.73, 1.91, 3.87, 4.0, 4.2}) { const double x = 0.8 * distance; const double y = 0.36 * distance; const double z = 0.48 * distance; const ModuleBase::Vector3 position(x, y, z); - for (int ma: {0, last_left}) + const auto all_values = candidate[i]->cal_overlap(origin, position); + const auto all_gradients = candidate[i]->cal_grad_overlap(origin, position); + const int count = (last_left + 1) * (last_right + 1); + ASSERT_EQ(all_values.size(), count); + ASSERT_EQ(all_gradients.size(), count); + for (int ma = 0; ma <= last_left; ++ma) { - for (int mb: {0, last_right}) + for (int mb = 0; mb <= last_right; ++mb) { const double expected = reference[i]->cal_overlap(origin, position, ma, mb); const double actual = candidate[i]->cal_overlap(origin, position, ma, mb); const double tolerance = 2e-11 * (1.0 + std::abs(expected)); EXPECT_NEAR(actual, expected, tolerance); + const int offset = ma * (last_right + 1) + mb; + EXPECT_NEAR(all_values[offset], expected, tolerance); const auto grad_expected = reference[i]->cal_grad_overlap(origin, position, ma, mb); const auto grad_actual = candidate[i]->cal_grad_overlap(origin, position, ma, mb); const double error = (grad_actual - grad_expected).norm(); const double grad_tolerance = 2e-10 * (1.0 + grad_expected.norm()); EXPECT_LT(error, grad_tolerance); + const double component_error = (all_gradients[offset] - grad_expected).norm(); + EXPECT_LT(component_error, grad_tolerance); } } } @@ -297,23 +306,22 @@ TEST(Center2Orb21, ProductGradientMatchesFiniteDifference) const ModuleBase::Vector3 origin(0.0, 0.0, 0.0); const ModuleBase::Vector3 position(0.813, 0.321, 0.473); const double step = 1e-5; - for (int ma: {0, 2, 4}) + const auto gradients = orbital.cal_grad_overlap(origin, position); + for (int axis = 0; axis < 3; ++axis) { - for (int mb: {0, 1, 2}) + auto plus = position; + auto minus = position; + plus[axis] += step; + minus[axis] -= step; + const auto upper = orbital.cal_overlap(origin, plus); + const auto lower = orbital.cal_overlap(origin, minus); + ASSERT_EQ(upper.size(), gradients.size()); + ASSERT_EQ(lower.size(), gradients.size()); + for (std::size_t i = 0; i < gradients.size(); ++i) { - for (int mc: {0, 3, 6}) - { - const auto gradient = orbital.cal_grad_overlap(origin, position, ma, mb, mc); - auto plus = position; - auto minus = position; - plus.x += step; - minus.x -= step; - const double upper = orbital.cal_overlap(origin, plus, ma, mb, mc); - const double lower = orbital.cal_overlap(origin, minus, ma, mb, mc); - const double finite_difference = (upper - lower) / (2 * step); - const double tolerance = 2e-7 * (1.0 + std::abs(finite_difference)); - EXPECT_NEAR(gradient.x, finite_difference, tolerance); - } + const double finite_difference = (upper[i] - lower[i]) / (2 * step); + const double tolerance = 2e-7 * (1.0 + std::abs(finite_difference)); + EXPECT_NEAR(gradients[i][axis], finite_difference, tolerance); } } }