From de7340f8bef98bc0457d1ab5616ab9d9da2f91e9 Mon Sep 17 00:00:00 2001 From: Marcello Puligheddu <33704456+marcelloPuligheddu@users.noreply.github.com> Date: Thu, 2 Jul 2026 13:30:11 +0100 Subject: [PATCH] PR for adding libGint to CP2K (#5446) --- CMakeLists.txt | 18 + cmake/modules/FindLibGint.cmake | 39 + docs/technologies/libraries.md | 9 + src/CMakeLists.txt | 3 + src/common/t_c_g0.F | 6 +- src/cp2k_info.F | 3 + src/hfx_energy_potential.F | 611 +++++++++------ src/hfx_pair_list_methods.F | 131 +++- src/hfx_types.F | 7 + src/input_constants.F | 5 + src/input_cp2k_hfx.F | 21 +- src/libgint_wrapper.F | 696 ++++++++++++++++++ src/offload/offload_library.c | 2 +- tests/QS/regtest-libgint/H1.inp | 67 ++ tests/QS/regtest-libgint/H2O-HFX-6.inp | 89 +++ tests/QS/regtest-libgint/TEST_FILES.toml | 4 + tests/QS/regtest-libgint/basis.dat | 10 + tests/QS/regtest-libgint/tidy | 18 + tests/TEST_DIRS | 1 + tools/docker/Dockerfile.test_cuda_A100 | 1 + tools/docker/Dockerfile.test_cuda_P100 | 1 + tools/docker/Dockerfile.test_cuda_V100 | 1 + .../Dockerfile.test_performance_cuda_A100 | 1 + .../Dockerfile.test_performance_cuda_P100 | 1 + .../Dockerfile.test_performance_cuda_V100 | 1 + tools/docker/generate_dockerfiles.py | 1 + tools/toolchain/install_cp2k_toolchain.sh | 15 +- .../scripts/stage4/install_libGint.sh | 77 ++ .../scripts/stage4/install_stage4.sh | 1 + 29 files changed, 1608 insertions(+), 232 deletions(-) create mode 100644 cmake/modules/FindLibGint.cmake create mode 100644 src/libgint_wrapper.F create mode 100644 tests/QS/regtest-libgint/H1.inp create mode 100644 tests/QS/regtest-libgint/H2O-HFX-6.inp create mode 100644 tests/QS/regtest-libgint/TEST_FILES.toml create mode 100644 tests/QS/regtest-libgint/basis.dat create mode 100644 tests/QS/regtest-libgint/tidy create mode 100755 tools/toolchain/scripts/stage4/install_libGint.sh diff --git a/CMakeLists.txt b/CMakeLists.txt index 0c1cf77296..0c6acfdee0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -239,6 +239,10 @@ cmake_dependent_option( "Enable CRAY power management framework with gpu support" "(NOT CP2K_USE_ACCEL MATCHES \"NONE\") AND (CP2K_USE_CRAY_PM_ENERGY)" OFF) +cmake_dependent_option( + CP2K_USE_LIBGINT "Enable LibGint support" ${CP2K_USE_EVERYTHING} + "CP2K_USE_ACCEL MATCHES \"CUDA\"" OFF) + set(CP2K_BLAS_VENDOR "auto" CACHE STRING "blas vendor/generic libraries") @@ -842,6 +846,10 @@ if(CP2K_USE_LIBINT2) include_directories(BEFORE ${_libint2_include_dirs} "${_libint2_mod_dir}") endif() +if(CP2K_USE_LIBGINT) + find_package(LibGint REQUIRED) +endif() + if(CP2K_USE_LIBFCI) find_package(libfci REQUIRED) endif() @@ -1214,6 +1222,12 @@ if(CP2K_USE_LIBINT2) " - Libraries: ${CP2K_LIBINT2_LINK_LIBRARIES}\n") endif() +if(CP2K_USE_LIBGINT) + message(" - LibGint\n" + " - include directories: ${CP2K_LIBGINT_INCLUDE_DIRS}\n" + " - libraries: ${CP2K_LIBGINT_LINK_LIBRARIES}\n\n") +endif() + if(CP2K_USE_LIBFCI) message(" - LibFCI\n" " - Include directories: ${CP2K_LIBFCI_INCLUDE_DIRS}\n" @@ -1334,6 +1348,10 @@ if(NOT CP2K_USE_LIBINT2) message(" - Libint2") endif() +if(NOT CP2K_USE_LIBGINT) + message(" - LibGint") +endif() + if(NOT CP2K_USE_LIBFCI) message(" - LibFCI") endif() diff --git a/cmake/modules/FindLibGint.cmake b/cmake/modules/FindLibGint.cmake new file mode 100644 index 0000000000..9170990afa --- /dev/null +++ b/cmake/modules/FindLibGint.cmake @@ -0,0 +1,39 @@ +#!-------------------------------------------------------------------------------------------------! +#! CP2K: A general program to perform molecular dynamics simulations ! +#! Copyright 2000-2026 CP2K developers group ! +#! ! +#! SPDX-License-Identifier: GPL-2.0-or-later ! +#!-------------------------------------------------------------------------------------------------! + +# Copyright (c) 2025- ETH Zurich +# +# author: Marcello Puligheddu, Sergey Chulkov + +include(FindPackageHandleStandardArgs) +include(cp2k_utils) + +cp2k_find_libraries(LIBGINT "cp2kGint") +cp2k_include_dirs(LIBGINT "libgint.mod") + +if(CP2K_LIBGINT_INCLUDE_DIRS) + find_package_handle_standard_args( + LibGint DEFAULT_MSG CP2K_LIBGINT_LINK_LIBRARIES CP2K_LIBGINT_INCLUDE_DIRS) +else() + find_package_handle_standard_args(LibGint DEFAULT_MSG + CP2K_LIBGINT_LINK_LIBRARIES) +endif() + +if(NOT TARGET cp2k::LibGint::libGint) + add_library(cp2k::LibGint::libGint INTERFACE IMPORTED) + set_target_properties( + cp2k::LibGint::libGint PROPERTIES INTERFACE_LINK_LIBRARIES + "${CP2K_LIBGINT_LINK_LIBRARIES}") + if(CP2K_LIBGINT_INCLUDE_DIRS) + set_target_properties( + cp2k::LibGint::libGint PROPERTIES INTERFACE_INCLUDE_DIRECTORIES + "${CP2K_LIBGINT_INCLUDE_DIRS}") + endif() +endif() + +mark_as_advanced(CP2K_LIBGINT_ROOT CP2K_LIBGINT_INCLUDE_DIRS + CP2K_LIBGINT_LINK_LIBRARIES CP2K_LIBGINT_LIBRARIES) diff --git a/docs/technologies/libraries.md b/docs/technologies/libraries.md index 878aa47edc..d32949b522 100644 --- a/docs/technologies/libraries.md +++ b/docs/technologies/libraries.md @@ -346,3 +346,12 @@ MiMiC - Multiscale simulation framework - For more information about the framework and supported programs see - Pass `-DCP2K_USE_MIMIC=ON` to CMake + +## libGint + +libGint - A library for the calculation of the Hartree Fock exchange on GPUs + +- Compared to a regular XF calculation, the changes needed in the input file are : +- &FORCE_EVAL &DFT &XC &HF HFX_LIBRARY libGint +- &FORCE_EVAL &DFT &XC &HF &MEMORY MAX_MEMORY X +- pass -DCP2K_USE_LIBGINT=ON\` to CMake. diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 27f110f3b7..e0e028ba77 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -336,6 +336,7 @@ list( libgrpp_integrals.F libint_2c_3c.F libint_wrapper.F + libgint_wrapper.F library_tests.F linesearch.F localization_tb.F @@ -1789,6 +1790,7 @@ target_link_libraries( $<$:CP2K::MIMIC::mcl> $<$:Libint2::int2> $<$:libfci::fci> + $<$:cp2k::LibGint::libGint> $<$:cp2k::cosma> $<$:${CP2K_DLAF_FORTRAN_TARGET}> $<$:cp2k::greenx> @@ -1861,6 +1863,7 @@ target_compile_definitions( $<$:__GREENX> $<$,$>:__FFTW3> $<$:__LIBINT> + $<$:__LIBGINT> $<$:__LIBFCI> $<$:__LIBTORCH> $<$:__COSMA> diff --git a/src/common/t_c_g0.F b/src/common/t_c_g0.F index ab7b77d6ec..db471aed8b 100644 --- a/src/common/t_c_g0.F +++ b/src/common/t_c_g0.F @@ -53,6 +53,7 @@ !> Nov 2008, 2009 Joost VandeVondele and Manuel Guidon !> May 2019 A. Bussy: Added a get_maxl_init function to get current status of nderiv_init and !> moved the file to common (made it accessible from aobasis, same place as gamma.F). +!> Oct 2025 M. Puligheddu: Added public qualifier to C0 to simplify reuse ! ************************************************************************************************** MODULE t_c_g0 USE kinds, ONLY: dp @@ -60,10 +61,13 @@ MODULE t_c_g0 #include "../base/base_uses.f90" IMPLICIT NONE + + REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE, SAVE, PUBLIC :: C0 + PRIVATE PUBLIC :: t_c_g0_n, init, free_C0, get_lmax_init - REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE, SAVE :: C0 + INTEGER, PARAMETER :: degree = 13 REAL(KIND=dp), PARAMETER :: target_error = 0.100000E-08 INTEGER, PARAMETER :: nderiv_max = 21 diff --git a/src/cp2k_info.F b/src/cp2k_info.F index d23f2c8960..afd451a01d 100644 --- a/src/cp2k_info.F +++ b/src/cp2k_info.F @@ -335,6 +335,9 @@ CONTAINS flags = TRIM(flags)//" openpmd" #endif +#if defined(__LIBGINT) + flags = TRIM(flags)//" libGint" +#endif END FUNCTION cp2k_flags ! ************************************************************************************************** diff --git a/src/hfx_energy_potential.F b/src/hfx_energy_potential.F index 200f378b6a..679a821530 100644 --- a/src/hfx_energy_potential.F +++ b/src/hfx_energy_potential.F @@ -9,6 +9,7 @@ !> \brief Routines to calculate HFX energy and potential !> \par History !> 11.2006 created [Manuel Guidon] +!> 10.2025 added libGint interface [M Puligheddu] !> \author Manuel Guidon ! ************************************************************************************************** MODULE hfx_energy_potential @@ -66,9 +67,15 @@ MODULE hfx_energy_potential hfx_pgf_list, hfx_pgf_product_list, hfx_potential_type, hfx_reset_memory_usage_counter, & hfx_screen_coeff_type, hfx_screening_type, hfx_task_list_type, hfx_type, init_t_c_g0_lmax, & log_zero, pair_list_type, pair_set_list_type + USE input_constants, ONLY: do_potential_mix_cl_trunc, & do_potential_truncated, & - hfx_do_eval_energy + do_potential_TShPSC, & + do_potential_short, & + do_potential_coulomb, & + hfx_do_eval_energy, & + hfx_library_is_libint, hfx_library_is_libGint, hfx_library_is_both + USE input_section_types, ONLY: section_vals_type USE kinds, ONLY: default_string_length, & dp, & @@ -90,6 +97,12 @@ MODULE hfx_energy_potential qs_ks_env_type USE t_c_g0, ONLY: init USE util, ONLY: sort + use, intrinsic :: iso_c_binding, only: c_ptr + + use cp2k_info, only: cp2k_flags + use offload_api, only: offload_get_chosen_device, offload_set_chosen_device, offload_activate_chosen_device + use libGint_wrapper, only: cp_libGint_init, libGint_update_env, libGint_set_density, libGint_coulomb4, & + libGint_update_fock_matrix, libGint_get_fock_matrix !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads @@ -154,7 +167,7 @@ CONTAINS INTEGER :: act_atomic_block_offset, act_set_offset, atomic_offset_ac, atomic_offset_ad, & atomic_offset_bc, atomic_offset_bd, bin, bits_max_val, buffer_left, buffer_size, & buffer_start, cache_size, current_counter, handle, handle_bin, handle_dist_ks, & - handle_getP, handle_load, handle_main, i, i_list_ij, i_list_kl, i_set_list_ij, & + handle_getP, handle_load, handle_main, i, inla, nla, i_list_ij, i_list_kl, i_set_list_ij, & i_set_list_ij_start, i_set_list_ij_stop, i_set_list_kl, i_set_list_kl_start, & i_set_list_kl_stop, i_thread, iatom, iatom_block, iatom_end, iatom_start, ikind, img, & iset, iw, j, jatom, jatom_block, jatom_end, jatom_start, jkind, jset, katom, katom_block, & @@ -177,7 +190,7 @@ CONTAINS tmp_i8(8) INTEGER(int_8), ALLOCATABLE, DIMENSION(:) :: tmp_task_list_cost INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, last_sgf_global, nimages, & - tmp_index + tmp_index, first_set_of_atom INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, lc_max, lc_min, ld_max, & ld_min, npgfa, npgfb, npgfc, npgfd, nsgfa, nsgfb, nsgfc, nsgfd, shm_block_offset INTEGER, DIMENSION(:, :), POINTER :: first_sgfb, nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d, & @@ -202,11 +215,13 @@ CONTAINS REAL(dp), DIMENSION(:), POINTER :: p_work REAL(dp), DIMENSION(:, :), POINTER :: full_density_alpha, full_density_beta, full_ks_alpha, & full_ks_beta, max_contraction, ptr_p_1, ptr_p_2, ptr_p_3, ptr_p_4, shm_pmax_atom, & - shm_pmax_block, sphi_b, zeta, zetb, zetc, zetd + shm_pmax_block, sphi_b, zeta, zetb, zetc, zetd, & + full_ks_alpha_from_gpu, full_ks_beta_from_gpu REAL(dp), DIMENSION(:, :, :), POINTER :: sphi_a_ext_set, sphi_b_ext_set, & sphi_c_ext_set, sphi_d_ext_set REAL(dp), DIMENSION(:, :, :, :), POINTER :: sphi_a_ext, sphi_b_ext, sphi_c_ext, & sphi_d_ext + REAL(dp), DIMENSION(:, :, :), POINTER :: gcc REAL(KIND=dp) :: coeffs_kind_max0 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(cell_type), POINTER :: cell @@ -247,6 +262,13 @@ CONTAINS TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(qs_ks_env_type), POINTER :: ks_env + logical :: screened + real(dp) :: max_abs_delta_KS + integer :: iatom_set, jatom_set, katom_set, latom_set, shared_dev + logical :: use_libint, use_libgint, use_only_alpha_spin + logical :: cp2k_was_compiled_with_cuda_offload, cp2k_was_compiled_with_libint + logical :: cp2k_was_compiled_with_libGint + NULLIFY (dft_control, matrix_ks_aux_fit_hfx) CALL timeset(routineN, handle) @@ -344,6 +366,8 @@ CONTAINS !$OMP nso,& !$OMP nco,& !$OMP full_ks_alpha,& +!$OMP full_ks_alpha_from_gpu,& +!$OMP full_ks_beta_from_gpu,& !$OMP full_ks_beta,& !$OMP n_threads,& !$OMP full_density_alpha,& @@ -379,7 +403,11 @@ CONTAINS !$OMP shm_pmax_block,& !$OMP shm_atomic_pair_list,& !$OMP shm_mem_compression_counter,& -!$OMP do_print_load_balance_info,& +!$OMP shared_dev, & +!$OMP do_print_load_balance_info,use_libgint,use_libint, & +!$OMP cp2k_was_compiled_with_cuda_offload, & +!$OMP cp2k_was_compiled_with_libint, & +!$OMP cp2k_was_compiled_with_libGint, & !$OMP my_nspins) & !$OMP PRIVATE(ln_10,i_thread,actual_x_data,do_periodic,screening_parameter,potential_parameter,& !$OMP general_parameter,load_balance_parameter,memory_parameter,cache_size,bits_max_val,& @@ -410,7 +438,9 @@ CONTAINS !$OMP tmp_screen_pgf1,tmp_screen_pgf2,cartesian_estimate,bintime_stop,iw,memsize_after,storage_counter_integrals,& !$OMP stor_count_int_disk,stor_count_max_val,ene_x_aa,ene_x_bb,mb_size_p,mb_size_f,mb_size_buffers,afac,ene_x_aa_diag,& !$OMP ene_x_bb_diag,act_atomic_block_offset,act_set_offset,j,handle_dist_ks,tmp_i8,tmp_i4,dft_control,& -!$OMP etmp,nkimages,img,bin,eps_scaling_str,eps_schwarz_min_str) +!$OMP etmp,nkimages,img,bin,eps_scaling_str,eps_schwarz_min_str, & +!$OMP i,inla,nla,gcc,first_set_of_atom,max_abs_delta_KS, screened, & +!$OMP iatom_set,jatom_set,katom_set,latom_set,use_only_alpha_spin) ln_10 = LOG(10.0_dp) i_thread = 0 @@ -463,6 +493,20 @@ CONTAINS neris_tmp = 0_int_8 max_val_memory = 1_int_8 + ! Set the integral engine + if (general_parameter%hfx_library == hfx_library_is_libint) then + use_libint = .true. + use_libGint = .false. + end if + if (general_parameter%hfx_library == hfx_library_is_libGint) then + use_libint = .false. + use_libGint = .true. + end if + if (general_parameter%hfx_library == hfx_library_is_both) then + use_libint = .true. + use_libGint = .true. + end if + max_am = basis_info%max_am CALL get_qs_env(qs_env=qs_env, & @@ -811,14 +855,14 @@ CONTAINS END IF END IF +!$OMP END MASTER +!$OMP BARRIER + !! Initialize schwarz screening matrices for near field estimates and boxing screening matrices !! for far field estimates. The update is only performed if the geomtry of the system changed. !! If the system is periodic, then the corresponding routines are called and some variables !! are initialized -!$OMP END MASTER -!$OMP BARRIER - IF (.NOT. shm_master_x_data%screen_funct_is_initialized) THEN CALL calc_pair_dist_radii(qs_env, basis_parameter, & shm_master_x_data%pair_dist_radii_pgf, max_set, max_pgf, eps_schwarz, & @@ -880,7 +924,30 @@ CONTAINS !$OMP END MASTER !$OMP BARRIER - !! Start caluclating integrals of the form (ab|cd) or (ij|kl) + use_only_alpha_spin = treat_lsd_in_core .or. .not. my_nspins == 2 + + ! Set up the gpu enviroment for libGint + if (use_libGint) then + + call offload_activate_chosen_device() + + ! At least once, and before any other call from all OMP threads + call cp_libGint_init(actual_x_data) + + ! At least at geo change + call libGint_update_env(fac, memory_parameter, do_periodic, cell, actual_x_data, & + nneighbors, max_pgf, natom, kind_of, particle_set, basis_parameter) + + ! Every cycle, comunicate the updated density to libGint + if (use_only_alpha_spin) then + call libGint_set_density(full_density_alpha) + else + call libGint_set_density(full_density_alpha, full_density_beta) + end if + end if + !!! + + !! Start calculating integrals of the form (ab|cd) or (ij|kl) !! In order to do so, there is a main four-loop structure that takes into account the two symmetries !! !! (ab|cd) = (ba|cd) = (ab|dc) = (ba|dc) @@ -918,6 +985,11 @@ CONTAINS !! lstart only the first time the loop is executed. All subsequent loops have to start with one or !! iatom and katom respectively. Therefore, we use flags like first_j_loop etc. + !! LibGint follows the same patter, distributing integrals over mpi and openmp with a shared-over-mpi + !! density and fock matrix. The main differences are that the fock matrix is not in a valid state + !! until after the libGint_get_fock_matrix has returned -in approximately a thousand lines-, + !! and the integrals are not saved in memory, but recomputed each cycle + do_dynamic_load_balancing = .TRUE. IF (n_threads == 1 .OR. do_disk_storage) do_dynamic_load_balancing = .FALSE. @@ -1088,7 +1160,9 @@ CONTAINS bins_left = .TRUE. do_it = .TRUE. bin = 0 + ! Main cycle starts here DO WHILE (bins_left) + ! Check if this thread can run this bin IF (.NOT. do_dynamic_load_balancing) THEN bin = bin + 1 IF (bin > my_bin_size) THEN @@ -1124,11 +1198,12 @@ CONTAINS END IF !$OMP END CRITICAL(hfxenergy_critical) END IF - IF (.NOT. do_it) CYCLE + !$OMP MASTER CALL timeset(routineN//"_bin", handle_bin) !$OMP END MASTER + bintime_start = m_walltime() my_istart = distribution_energy%istart my_current_counter = 0 @@ -1280,6 +1355,8 @@ CONTAINS sphi_d_u2 = UBOUND(sphi_d_ext, 2) sphi_d_u3 = UBOUND(sphi_d_ext, 3) + ! Note: shm_atomic_block_offset is already symmetric + ! see get_atomic_block_maps in hfx_communication.F atomic_offset_bd = shm_atomic_block_offset(jatom, latom) atomic_offset_bc = shm_atomic_block_offset(jatom, katom) atomic_offset_ad = shm_atomic_block_offset(iatom, latom) @@ -1355,7 +1432,6 @@ CONTAINS END IF !! At this stage, check for memory used in compression - IF (my_geo_change) THEN IF (.NOT. memory_parameter%do_all_on_the_fly) THEN ! ** We know the maximum amount of integrals that we can store per MPI-process @@ -1456,189 +1532,77 @@ CONTAINS IF (max_val2 < log10_eps_schwarz) CYCLE pmax_entry = EXP(log10_pmax*ln_10) - !! store current number of integrals, update total number and number of integrals in buffer - current_counter = nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset) - IF (buffer_overflow) THEN - neris_onthefly = neris_onthefly + current_counter - END IF - - !! Get integrals from buffer and update Kohn-Sham matrix - IF (.NOT. buffer_overflow .AND. .NOT. my_geo_change) THEN - nints = current_counter - IF (.NOT. use_disk_storage) THEN - CALL hfx_get_single_cache_element( & - estimate_to_store_int, 6, & - maxval_cache, maxval_container, memory_parameter%actual_memory_usage, & - use_disk_storage) - ELSE - CALL hfx_get_single_cache_element( & - estimate_to_store_int, 6, & - maxval_cache_disk, maxval_container_disk, memory_parameter%actual_memory_usage_disk, & - use_disk_storage) - END IF - spherical_estimate = SET_EXPONENT(1.0_dp, estimate_to_store_int + 1) - IF (spherical_estimate*pmax_entry < eps_schwarz) CYCLE - nbits = EXPONENT(ANINT(spherical_estimate*pmax_entry/eps_storage)) + 1 - buffer_left = nints - buffer_start = 1 - IF (.NOT. use_disk_storage) THEN - neris_incore = neris_incore + INT(nints, int_8) - ELSE - neris_disk = neris_disk + INT(nints, int_8) - END IF - DO WHILE (buffer_left > 0) - buffer_size = MIN(buffer_left, cache_size) - IF (.NOT. use_disk_storage) THEN - CALL hfx_get_mult_cache_elements(primitive_integrals(buffer_start), & - buffer_size, nbits, & - integral_caches(nbits), & - integral_containers(nbits), & - eps_storage, pmax_entry, & - memory_parameter%actual_memory_usage, & - use_disk_storage) - ELSE - CALL hfx_get_mult_cache_elements(primitive_integrals(buffer_start), & - buffer_size, nbits, & - integral_caches_disk(nbits), & - integral_containers_disk(nbits), & - eps_storage, pmax_entry, & - memory_parameter%actual_memory_usage_disk, & - use_disk_storage) - END IF - buffer_left = buffer_left - buffer_size - buffer_start = buffer_start + buffer_size - END DO - END IF - !! Calculate integrals if we run out of buffer or the geometry did change - IF (my_geo_change .OR. buffer_overflow) THEN - max_contraction_val = max_contraction(iset, iatom)* & - max_contraction(jset, jatom)* & - max_contraction(kset, katom)* & - max_contraction(lset, latom)*pmax_entry - tmp_R_1 => radii_pgf(:, :, jset, iset, jkind, ikind) - tmp_R_2 => radii_pgf(:, :, lset, kset, lkind, kkind) + ! libGint can only run in direct mode at the moment + ! so there is no need to check the buffers + if (use_libGint) then tmp_screen_pgf1 => screen_coeffs_pgf(:, :, jset, iset, jkind, ikind) tmp_screen_pgf2 => screen_coeffs_pgf(:, :, lset, kset, lkind, kkind) + ! TODO a more modern approach in which the call is more like + ! call coulom_4_gpu( HFX_data, iset, jset, kset, lset ) + call libGint_coulomb4( & + iatom, jatom, katom, latom, iset, jset, kset, lset, & + ra, rb, rc, rd, & + npgfa(iset), npgfb(jset), npgfc(kset), npgfd(lset), & + potential_parameter, & + screen_coeffs_set(jset, iset, jkind, ikind)%x, & + screen_coeffs_set(lset, kset, lkind, kkind)%x, & + log10_pmax, log10_eps_schwarz, & + tmp_screen_pgf1, tmp_screen_pgf2, & + actual_x_data%neighbor_cells, cell, do_periodic, screened) - CALL coulomb4(private_lib, ra, rb, rc, rd, npgfa(iset), npgfb(jset), npgfc(kset), npgfd(lset), & - la_min(iset), la_max(iset), lb_min(jset), lb_max(jset), & - lc_min(kset), lc_max(kset), ld_min(lset), ld_max(lset), & - nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), & - sphi_a_u1, sphi_a_u2, sphi_a_u3, & - sphi_b_u1, sphi_b_u2, sphi_b_u3, & - sphi_c_u1, sphi_c_u2, sphi_c_u3, & - sphi_d_u1, sphi_d_u2, sphi_d_u3, & - zeta(1:npgfa(iset), iset), zetb(1:npgfb(jset), jset), & - zetc(1:npgfc(kset), kset), zetd(1:npgfd(lset), lset), & - primitive_integrals, & - potential_parameter, & - actual_x_data%neighbor_cells, screen_coeffs_set(jset, iset, jkind, ikind)%x, & - screen_coeffs_set(lset, kset, lkind, kkind)%x, eps_schwarz, & - max_contraction_val, cartesian_estimate, cell, neris_tmp, & - log10_pmax, log10_eps_schwarz, & - tmp_R_1, tmp_R_2, tmp_screen_pgf1, tmp_screen_pgf2, & - pgf_list_ij, pgf_list_kl, pgf_product_list, & - nsgfl_a(:, iset), nsgfl_b(:, jset), & - nsgfl_c(:, kset), nsgfl_d(:, lset), & - sphi_a_ext_set, & - sphi_b_ext_set, & - sphi_c_ext_set, & - sphi_d_ext_set, & - ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp, & - nimages, do_periodic, p_work) + if (.not. screened) then - nints = nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset) - neris_total = neris_total + nints - nprim_ints = nprim_ints + neris_tmp -! IF (cartesian_estimate == 0.0_dp) cartesian_estimate = TINY(cartesian_estimate) -! estimate_to_store_int = EXPONENT(cartesian_estimate) -! estimate_to_store_int = MAX(estimate_to_store_int, -15_int_8) -! cartesian_estimate = SET_EXPONENT(1.0_dp, estimate_to_store_int+1) -! IF (.NOT. buffer_overflow .AND. my_geo_change) THEN -! IF (cartesian_estimate < eps_schwarz) THEN -! IF (.NOT. use_disk_storage) THEN -! CALL hfx_add_single_cache_element( & -! estimate_to_store_int, 6, & -! maxval_cache, maxval_container, memory_parameter%actual_memory_usage, & -! use_disk_storage, max_val_memory) -! ELSE -! CALL hfx_add_single_cache_element( & -! estimate_to_store_int, 6, & -! maxval_cache_disk, maxval_container_disk, memory_parameter%actual_memory_usage_disk, & -! use_disk_storage) -! END IF -! END IF -! END IF -! -! IF (cartesian_estimate < eps_schwarz) CYCLE + ! TODO a more modern approach in which the call is more like + ! call update_fock_matrix_gpu( HFX_data, iset, jset, kset, lset ) + call libGint_update_fock_matrix( & + symm_fac, & + iatom, jatom, katom, latom, & + iset, jset, kset, lset, & + atomic_offset_ac, atomic_offset_ad, atomic_offset_bc, atomic_offset_bd, & + offset_ac_set, offset_ad_set, offset_bc_set, offset_bd_set, & + nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), & + la_min(iset), la_max(iset), lb_min(jset), lb_max(jset), & + lc_min(kset), lc_max(kset), ld_min(lset), ld_max(lset), & + nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d) - !! Compress the array for storage - spherical_estimate = 0.0_dp - DO i = 1, nints - spherical_estimate = MAX(spherical_estimate, ABS(primitive_integrals(i))) - END DO + end if + end if + ! Run using libint as integral engine + if_use_libint: if (use_libint) then + !! store current number of integrals, update total number and number of integrals in buffer + current_counter = nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset) + IF (buffer_overflow) THEN + neris_onthefly = neris_onthefly + current_counter + END IF - IF (spherical_estimate == 0.0_dp) spherical_estimate = TINY(spherical_estimate) - estimate_to_store_int = EXPONENT(spherical_estimate) - estimate_to_store_int = MAX(estimate_to_store_int, -15_int_8) - - IF (.NOT. buffer_overflow .AND. my_geo_change) THEN + !! Get integrals from buffer and update Kohn-Sham matrix + IF (.NOT. buffer_overflow .AND. .NOT. my_geo_change) THEN + nints = current_counter IF (.NOT. use_disk_storage) THEN - CALL hfx_add_single_cache_element( & + CALL hfx_get_single_cache_element( & estimate_to_store_int, 6, & maxval_cache, maxval_container, memory_parameter%actual_memory_usage, & - use_disk_storage, max_val_memory) + use_disk_storage) ELSE - CALL hfx_add_single_cache_element( & + CALL hfx_get_single_cache_element( & estimate_to_store_int, 6, & maxval_cache_disk, maxval_container_disk, memory_parameter%actual_memory_usage_disk, & use_disk_storage) END IF - END IF - spherical_estimate = SET_EXPONENT(1.0_dp, estimate_to_store_int + 1) - IF (spherical_estimate*pmax_entry < eps_schwarz) CYCLE - IF (.NOT. buffer_overflow) THEN + spherical_estimate = SET_EXPONENT(1.0_dp, estimate_to_store_int + 1) + IF (spherical_estimate*pmax_entry < eps_schwarz) CYCLE nbits = EXPONENT(ANINT(spherical_estimate*pmax_entry/eps_storage)) + 1 - - ! In case of a tight eps_storage threshold the number of significant - ! bits in the integer number NINT(value*pmax_entry/eps_storage) may - ! exceed the width of the storage element. As the compression algorithm - ! is designed for IEEE 754 double precision numbers, a 64-bit signed - ! integer variable which is used to store the result of this float-to- - ! integer conversion (we have no wish to use more memory for storing - ! compressed ERIs than it is needed for uncompressed ERIs) may overflow. - ! Abort with a meaningful message when it happens. - ! - ! The magic number 63 stands for the number of magnitude bits - ! (64 bits minus one sign bit). - IF (nbits > 63) THEN - WRITE (eps_schwarz_min_str, '(ES10.3E2)') & - spherical_estimate*pmax_entry/ & - (SET_EXPONENT(1.0_dp, 63)*memory_parameter%eps_storage_scaling) - - WRITE (eps_scaling_str, '(ES10.3E2)') & - spherical_estimate*pmax_entry/(SET_EXPONENT(1.0_dp, 63)*eps_schwarz) - - CALL cp_abort(__LOCATION__, & - "Overflow during ERI's compression. Please use a larger "// & - "EPS_SCHWARZ threshold (above "//TRIM(ADJUSTL(eps_schwarz_min_str))// & - ") or increase the EPS_STORAGE_SCALING factor above "// & - TRIM(ADJUSTL(eps_scaling_str))//".") - END IF - buffer_left = nints buffer_start = 1 IF (.NOT. use_disk_storage) THEN neris_incore = neris_incore + INT(nints, int_8) -! neris_incore = neris_incore+nints ELSE neris_disk = neris_disk + INT(nints, int_8) -! neris_disk = neris_disk+nints END IF DO WHILE (buffer_left > 0) - buffer_size = MIN(buffer_left, CACHE_SIZE) + buffer_size = MIN(buffer_left, cache_size) IF (.NOT. use_disk_storage) THEN - CALL hfx_add_mult_cache_elements(primitive_integrals(buffer_start), & + CALL hfx_get_mult_cache_elements(primitive_integrals(buffer_start), & buffer_size, nbits, & integral_caches(nbits), & integral_containers(nbits), & @@ -1646,7 +1610,7 @@ CONTAINS memory_parameter%actual_memory_usage, & use_disk_storage) ELSE - CALL hfx_add_mult_cache_elements(primitive_integrals(buffer_start), & + CALL hfx_get_mult_cache_elements(primitive_integrals(buffer_start), & buffer_size, nbits, & integral_caches_disk(nbits), & integral_containers_disk(nbits), & @@ -1657,65 +1621,193 @@ CONTAINS buffer_left = buffer_left - buffer_size buffer_start = buffer_start + buffer_size END DO - ELSE - !! In order to be consistent with in-core part, round all the eris wrt. eps_schwarz + END IF + !! Calculate integrals if we run out of buffer or the geometry did change + IF (my_geo_change .OR. buffer_overflow) THEN + + max_contraction_val = max_contraction(iset, iatom)* & + max_contraction(jset, jatom)* & + max_contraction(kset, katom)* & + max_contraction(lset, latom)*pmax_entry + tmp_R_1 => radii_pgf(:, :, jset, iset, jkind, ikind) + tmp_R_2 => radii_pgf(:, :, lset, kset, lkind, kkind) + tmp_screen_pgf1 => screen_coeffs_pgf(:, :, jset, iset, jkind, ikind) + tmp_screen_pgf2 => screen_coeffs_pgf(:, :, lset, kset, lkind, kkind) + + CALL coulomb4(private_lib, ra, rb, rc, rd, npgfa(iset), npgfb(jset), npgfc(kset), npgfd(lset), & + la_min(iset), la_max(iset), lb_min(jset), lb_max(jset), & + lc_min(kset), lc_max(kset), ld_min(lset), ld_max(lset), & + nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), & + sphi_a_u1, sphi_a_u2, sphi_a_u3, & + sphi_b_u1, sphi_b_u2, sphi_b_u3, & + sphi_c_u1, sphi_c_u2, sphi_c_u3, & + sphi_d_u1, sphi_d_u2, sphi_d_u3, & + zeta(1:npgfa(iset), iset), zetb(1:npgfb(jset), jset), & + zetc(1:npgfc(kset), kset), zetd(1:npgfd(lset), lset), & + primitive_integrals, & + potential_parameter, & + actual_x_data%neighbor_cells, screen_coeffs_set(jset, iset, jkind, ikind)%x, & + screen_coeffs_set(lset, kset, lkind, kkind)%x, eps_schwarz, & + max_contraction_val, cartesian_estimate, cell, neris_tmp, & + log10_pmax, log10_eps_schwarz, & + tmp_R_1, tmp_R_2, tmp_screen_pgf1, tmp_screen_pgf2, & + pgf_list_ij, pgf_list_kl, pgf_product_list, & + nsgfl_a(:, iset), nsgfl_b(:, jset), & + nsgfl_c(:, kset), nsgfl_d(:, lset), & + sphi_a_ext_set, & + sphi_b_ext_set, & + sphi_c_ext_set, & + sphi_d_ext_set, & + ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp, & + nimages, do_periodic, p_work) + + nints = nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset) + neris_total = neris_total + nints + nprim_ints = nprim_ints + neris_tmp + + !! Compress the array for storage + spherical_estimate = 0.0_dp DO i = 1, nints - primitive_integrals(i) = primitive_integrals(i)*pmax_entry - IF (ABS(primitive_integrals(i)) > eps_storage) THEN - primitive_integrals(i) = ANINT(primitive_integrals(i)/eps_storage, dp)*eps_storage/pmax_entry - ELSE - primitive_integrals(i) = 0.0_dp - END IF + spherical_estimate = MAX(spherical_estimate, ABS(primitive_integrals(i))) END DO - END IF - END IF - !!! DEBUG, print out primitive integrals and indices. Only works serial no OMP !!! - IF (.FALSE.) THEN - CALL print_integrals( & - iatom, jatom, katom, latom, shm_set_offset, shm_atomic_block_offset, & - iset, jset, kset, lset, nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), primitive_integrals) - END IF - IF (.NOT. is_anti_symmetric) THEN - !! Update Kohn-Sham matrix - CALL update_fock_matrix( & - nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), & - fac, symm_fac, full_density_alpha(:, 1), full_ks_alpha(:, 1), & - primitive_integrals, pbd_buf, pbc_buf, pad_buf, pac_buf, kbd_buf, & - kbc_buf, kad_buf, kac_buf, iatom, jatom, katom, latom, & - iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set, offset_ac_set, & - atomic_offset_bd, atomic_offset_bc, atomic_offset_ad, atomic_offset_ac) - IF (.NOT. treat_lsd_in_core) THEN - IF (my_nspins == 2) THEN - CALL update_fock_matrix( & - nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), & - fac, symm_fac, full_density_beta(:, 1), full_ks_beta(:, 1), & - primitive_integrals, pbd_buf, pbc_buf, pad_buf, pac_buf, kbd_buf, & - kbc_buf, kad_buf, kac_buf, iatom, jatom, katom, latom, & - iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set, offset_ac_set, & - atomic_offset_bd, atomic_offset_bc, atomic_offset_ad, atomic_offset_ac) + + IF (spherical_estimate == 0.0_dp) spherical_estimate = TINY(spherical_estimate) + estimate_to_store_int = EXPONENT(spherical_estimate) + estimate_to_store_int = MAX(estimate_to_store_int, -15_int_8) + + IF (.NOT. buffer_overflow .AND. my_geo_change) THEN + IF (.NOT. use_disk_storage) THEN + CALL hfx_add_single_cache_element( & + estimate_to_store_int, 6, & + maxval_cache, maxval_container, memory_parameter%actual_memory_usage, & + use_disk_storage, max_val_memory) + ELSE + CALL hfx_add_single_cache_element( & + estimate_to_store_int, 6, & + maxval_cache_disk, maxval_container_disk, memory_parameter%actual_memory_usage_disk, & + use_disk_storage) + END IF + END IF + spherical_estimate = SET_EXPONENT(1.0_dp, estimate_to_store_int + 1) + IF (spherical_estimate*pmax_entry < eps_schwarz) CYCLE + IF (.NOT. buffer_overflow) THEN + nbits = EXPONENT(ANINT(spherical_estimate*pmax_entry/eps_storage)) + 1 + + ! In case of a tight eps_storage threshold the number of significant + ! bits in the integer number NINT(value*pmax_entry/eps_storage) may + ! exceed the width of the storage element. As the compression algorithm + ! is designed for IEEE 754 double precision numbers, a 64-bit signed + ! integer variable which is used to store the result of this float-to- + ! integer conversion (we have no wish to use more memory for storing + ! compressed ERIs than it is needed for uncompressed ERIs) may overflow. + ! Abort with a meaningful message when it happens. + ! + ! The magic number 63 stands for the number of magnitude bits + ! (64 bits minus one sign bit). + IF (nbits > 63) THEN + WRITE (eps_schwarz_min_str, '(ES10.3E2)') & + spherical_estimate*pmax_entry/ & + (SET_EXPONENT(1.0_dp, 63)*memory_parameter%eps_storage_scaling) + + WRITE (eps_scaling_str, '(ES10.3E2)') & + spherical_estimate*pmax_entry/(SET_EXPONENT(1.0_dp, 63)*eps_schwarz) + + CALL cp_abort(__LOCATION__, & + "Overflow during ERI's compression. Please use a larger "// & + "EPS_SCHWARZ threshold (above "//TRIM(ADJUSTL(eps_schwarz_min_str))// & + ") or increase the EPS_STORAGE_SCALING factor above "// & + TRIM(ADJUSTL(eps_scaling_str))//".") + END IF + + buffer_left = nints + buffer_start = 1 + IF (.NOT. use_disk_storage) THEN + neris_incore = neris_incore + INT(nints, int_8) + ELSE + neris_disk = neris_disk + INT(nints, int_8) + END IF + DO WHILE (buffer_left > 0) + buffer_size = MIN(buffer_left, CACHE_SIZE) + IF (.NOT. use_disk_storage) THEN + CALL hfx_add_mult_cache_elements(primitive_integrals(buffer_start), & + buffer_size, nbits, & + integral_caches(nbits), & + integral_containers(nbits), & + eps_storage, pmax_entry, & + memory_parameter%actual_memory_usage, & + use_disk_storage) + ELSE + CALL hfx_add_mult_cache_elements(primitive_integrals(buffer_start), & + buffer_size, nbits, & + integral_caches_disk(nbits), & + integral_containers_disk(nbits), & + eps_storage, pmax_entry, & + memory_parameter%actual_memory_usage_disk, & + use_disk_storage) + END IF + buffer_left = buffer_left - buffer_size + buffer_start = buffer_start + buffer_size + END DO + ELSE + !! In order to be consistent with in-core part, round all the eris wrt. eps_schwarz + DO i = 1, nints + primitive_integrals(i) = primitive_integrals(i)*pmax_entry + IF (ABS(primitive_integrals(i)) > eps_storage) THEN + primitive_integrals(i) = ANINT(primitive_integrals(i)/eps_storage, dp)*eps_storage/pmax_entry + ELSE + primitive_integrals(i) = 0.0_dp + END IF + END DO END IF END IF - ELSE - !! Update Kohn-Sham matrix - CALL update_fock_matrix_as( & - nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), & - fac, symm_fac, full_density_alpha(:, 1), full_ks_alpha(:, 1), & - primitive_integrals, pbd_buf, pbc_buf, pad_buf, pac_buf, kbd_buf, & - kbc_buf, kad_buf, kac_buf, iatom, jatom, katom, latom, & - iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set, offset_ac_set, & - atomic_offset_bd, atomic_offset_bc, atomic_offset_ad, atomic_offset_ac) - IF (.NOT. treat_lsd_in_core) THEN - IF (my_nspins == 2) THEN - CALL update_fock_matrix_as( & - nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), & - fac, symm_fac, full_density_beta(:, 1), full_ks_beta(:, 1), & - primitive_integrals, pbd_buf, pbc_buf, pad_buf, pac_buf, kbd_buf, & - kbc_buf, kad_buf, kac_buf, iatom, jatom, katom, latom, & - iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set, offset_ac_set, & - atomic_offset_bd, atomic_offset_bc, atomic_offset_ad, atomic_offset_ac) + !!! DEBUG, print out primitive integrals and indices. Only works serial no OMP !!! + IF (.FALSE.) THEN + CALL print_integrals( & + iatom, jatom, katom, latom, shm_set_offset, shm_atomic_block_offset, & + iset, jset, kset, lset, nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), primitive_integrals) + END IF + IF (.NOT. is_anti_symmetric) THEN + !! Update Kohn-Sham matrix + CALL update_fock_matrix( & + nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), & + fac, symm_fac, full_density_alpha(:, 1), full_ks_alpha(:, 1), & + primitive_integrals, pbd_buf, pbc_buf, pad_buf, pac_buf, kbd_buf, & + kbc_buf, kad_buf, kac_buf, iatom, jatom, katom, latom, & + iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set, offset_ac_set, & + atomic_offset_bd, atomic_offset_bc, atomic_offset_ad, atomic_offset_ac) + IF (.NOT. treat_lsd_in_core) THEN + IF (my_nspins == 2) THEN + CALL update_fock_matrix( & + nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), & + fac, symm_fac, full_density_beta(:, 1), full_ks_beta(:, 1), & + primitive_integrals, pbd_buf, pbc_buf, pad_buf, pac_buf, kbd_buf, & + kbc_buf, kad_buf, kac_buf, iatom, jatom, katom, latom, & + iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set, offset_ac_set, & + atomic_offset_bd, atomic_offset_bc, atomic_offset_ad, atomic_offset_ac) + END IF + END IF + ELSE + !! Update Kohn-Sham matrix + CALL update_fock_matrix_as( & + nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), & + fac, symm_fac, full_density_alpha(:, 1), full_ks_alpha(:, 1), & + primitive_integrals, pbd_buf, pbc_buf, pad_buf, pac_buf, kbd_buf, & + kbc_buf, kad_buf, kac_buf, iatom, jatom, katom, latom, & + iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set, offset_ac_set, & + atomic_offset_bd, atomic_offset_bc, atomic_offset_ad, atomic_offset_ac) + IF (.NOT. treat_lsd_in_core) THEN + IF (my_nspins == 2) THEN + CALL update_fock_matrix_as( & + nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), & + fac, symm_fac, full_density_beta(:, 1), full_ks_beta(:, 1), & + primitive_integrals, pbd_buf, pbc_buf, pad_buf, pac_buf, kbd_buf, & + kbc_buf, kad_buf, kac_buf, iatom, jatom, katom, latom, & + iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set, offset_ac_set, & + atomic_offset_bd, atomic_offset_bc, atomic_offset_ad, atomic_offset_ac) + END IF END IF END IF - END IF + END IF if_use_libint END DO ! i_set_list_kl END DO ! i_set_list_ij IF (do_disk_storage) THEN @@ -1736,6 +1828,74 @@ CONTAINS !$OMP END MASTER END DO !bin + !!! + ! IF running both libint and libGint, uses the values from libint to check libGint values + if (use_libGint .and. use_libint) then + ! barrier needed after allocate full_ks_alpha(beta)_from_gpu before libGint get K + ! to feed libgint a valid position for full_ks_alpha(beta)_from_gpu + allocate (full_ks_alpha_from_gpu, Mold=full_ks_alpha) + !$OMP BARRIER + if (use_only_alpha_spin) then + call libGint_get_fock_matrix(full_ks_alpha_from_gpu) + else + allocate (full_ks_beta_from_gpu, Mold=full_ks_beta) + !$OMP BARRIER + call libGint_get_fock_matrix(full_ks_alpha_from_gpu, full_ks_beta_from_gpu) + end if + + !$OMP BARRIER + !$omp single + max_abs_delta_KS = 0.0_dp + do i = 1, size(full_ks_alpha) + max_abs_delta_KS = max(max_abs_delta_ks, abs(full_ks_alpha(i, 1) - full_ks_alpha_from_gpu(i, 1))) + end do + write (UNIT=iw, FMT=*) " ------------------------------------ " + write (UNIT=iw, FMT=*) " | rank ", para_env%mepos + write (UNIT=iw, FMT=*) " | max abs dif: ", max_abs_delta_KS + write (UNIT=iw, FMT=*) " ------------------------------------ " + + ! Print all F if delta is too large or (suspiciously) too small + if (max_abs_delta_KS > 100*eps_schwarz .or. max_abs_delta_KS < 1.E-20) then + write (UNIT=iw, FMT=*) " K from CPU | K from GPU | Delta | Ratio " + do i = 1, size(full_ks_alpha) + write (UNIT=iw, FMT=*) i, full_ks_alpha(i, 1), full_ks_alpha_from_gpu(i, 1), & + full_ks_alpha(i, 1) - full_ks_alpha_from_gpu(i, 1), & + full_ks_alpha(i, 1)/full_ks_alpha_from_gpu(i, 1) + end do + end if + + if (.not. use_only_alpha_spin) then + max_abs_delta_KS = 0.0_dp + do i = 1, size(full_ks_beta) + max_abs_delta_KS = max(max_abs_delta_ks, abs(full_ks_beta(i, 1) - full_ks_beta_from_gpu(i, 1))) + end do + write (UNIT=iw, FMT=*) " ------------------------------------ " + write (UNIT=iw, FMT=*) " | rank ", para_env%mepos + write (UNIT=iw, FMT=*) " | max abs dif B: ", max_abs_delta_KS + write (UNIT=iw, FMT=*) " ------------------------------------ " + + ! Print all F if delta is too large or (suspiciously) too small + if (max_abs_delta_KS > 100*eps_schwarz .or. max_abs_delta_KS < 1.E-18) then + write (UNIT=iw, FMT=*) " BETA K from CPU | K from GPU | Delta | Ratio " + do i = 1, size(full_ks_beta) + write (UNIT=iw, FMT=*) i, full_ks_beta(i, 1), full_ks_beta_from_gpu(i, 1), & + full_ks_beta(i, 1) - full_ks_beta_from_gpu(i, 1), & + full_ks_beta(i, 1)/full_ks_beta_from_gpu(i, 1) + end do + end if + end if + !$omp end single + end if + + ! Copy the KS matrix from device to host. + if (use_libGint .and. .not. use_libint) then + if (use_only_alpha_spin) then + call libGint_get_fock_matrix(full_ks_alpha) + else + call libGint_get_fock_matrix(full_ks_alpha, full_ks_beta) + end if + end if + !$OMP MASTER logger => cp_get_default_logger() do_print_load_balance_info = .FALSE. @@ -2318,6 +2478,7 @@ CONTAINS neris_tmp = 0 primitive_integrals = 0.0_dp max_l = la_max + lb_max + lc_max + ld_max + DO list_ij = 1, nelements_ij ZetaInv = pgf_list_ij(list_ij)%ZetaInv ipgf = pgf_list_ij(list_ij)%ipgf diff --git a/src/hfx_pair_list_methods.F b/src/hfx_pair_list_methods.F index c37159eb32..d0ea682b91 100644 --- a/src/hfx_pair_list_methods.F +++ b/src/hfx_pair_list_methods.F @@ -40,15 +40,142 @@ MODULE hfx_pair_list_methods build_pair_list_pgf, & build_pgf_product_list, & build_atomic_pair_list, & - pgf_product_list_size + pgf_product_list_size, & + build_pair_list_pbc_pgf, & + bra_t, allocate_bra ! an initial estimate for the size of the product list INTEGER, SAVE :: pgf_product_list_size = 128 -!*** + ! Store screened primitive gaussian function pairs between two shells (A and B) + TYPE :: Bra_t + REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pgf_scr + INTEGER, ALLOCATABLE, DIMENSION(:, :) :: pgf_idx, cell_idx + INTEGER :: cell_cnt = 0, pgf_cnt = 0 + END TYPE Bra_t CONTAINS +! ************************************************************************************************** +!> \brief Allocates all arrays within a Bra_t structure. +!> +!> This subroutine ensures that a given `Bra_t` pointer has allocated memory +!> for storing primitive Gaussian pair indices, screening data, and cell indices. +!> Existing allocations are safely deallocated first. +!> +!> \param[in,out] bra Pointer to the Bra_t structure to allocate. +!> \param[in] max_npgf Maximum number of primitives per shell. +!> \param[in] max_ncell Maximum number of neighbor cells. +!> +!> \note The arrays are allocated as follows: +!> - `pgf_idx(3, max_npgf^2 * max_ncell)` +!> - `pgf_scr(5, max_npgf^2 * max_ncell)` +!> - `cell_idx(3, max_ncell)` +!> +! ************************************************************************************************** + SUBROUTINE allocate_Bra(bra, max_npgf, max_ncell) + TYPE(bra_t), POINTER :: bra + INTEGER :: max_npgf, max_ncell + + IF (ALLOCATED(bra%pgf_idx)) DEALLOCATE (bra%pgf_idx) ! TODO add if n_prm changed + IF (ALLOCATED(bra%pgf_scr)) DEALLOCATE (bra%pgf_scr) ! TODO add if n_prm changed + IF (ALLOCATED(bra%cell_idx)) DEALLOCATE (bra%cell_idx) ! TODO add if n_prm changed + ALLOCATE (bra%pgf_idx(3, max_npgf*max_npgf*max_ncell)) + ALLOCATE (bra%pgf_scr(5, max_npgf*max_npgf*max_ncell)) + ALLOCATE (bra%cell_idx(3, max_ncell)) + END SUBROUTINE allocate_Bra + +! ************************************************************************************************** +!> \brief Builds a screened list of primitives from centers A and B, intersecting with another shell +!> +!> This subroutine populates the `Bra_t` structure with indices and screening data +!> for primitive Gaussian pairs (ipgf,jpgf,n3) that pass the Schwarz screening criterion +!> +!> \param[in] npgfa Number of primitives in shell A +!> \param[in] npgfb Number of primitives in shell B +!> \param[out] bra Pointer to Bra_t structure to populate +!> \param[in] screen1 Screening coeffiecients for the AB pair based on most diffuse pgf +!> \param[in] screen2 Screening coeffiecients for the CD pair based on most diffuse pgf +!> \param[in] pgf Screening coeeficents for primitive gaussians +!> \param[in] log10_pmax Log10 of maximum integral prefactor +!> \param[in] log10_eps_schwarz Log10 of Schwarz screening threshold +!> \param[in] ra Cartesian coordinates of center A +!> \param[in] rb Cartesian coordinates of center B +!> \param[out] nelements Total number of valid primitives found +!> \param[in] neighbor_cells Array of lattice vectors +!> \param[in] do_periodic Logical flag for periodic boundary conditions on B +!> +!> \note Loops over all cells, shifting the B shell position accordingly +!> Only pairs satisfying the screening condition are stored. +!> +!> \note Each valid cell contributes a block of pairs stored contiguously in +!> bra%pgf_idx, with cell metadata stored in bra%cell_idx +!> \note screening primitive pairs (ipgf, jpgf) using both coarse and fine screening thresholds +! ************************************************************************************************** + SUBROUTINE build_pair_list_pbc_pgf(npgfa, npgfb, bra, screen1, screen2, pgf, & + log10_pmax, log10_eps_schwarz, ra, rb, nelements, & + neighbor_cells, do_periodic) + + INTEGER, INTENT(IN) :: npgfa, npgfb + TYPE(bra_t), POINTER :: bra + REAL(dp), INTENT(IN) :: screen1(2), screen2(2) + TYPE(hfx_screen_coeff_type), DIMENSION(:, :), & + POINTER :: pgf + REAL(dp), INTENT(IN) :: log10_pmax, log10_eps_schwarz, ra(3), & + rb(3) + INTEGER, INTENT(OUT) :: nelements + TYPE(hfx_cell_type), DIMENSION(:), POINTER :: neighbor_cells + LOGICAL, INTENT(IN) :: do_periodic + + INTEGER :: cell_cnt, cell_idx, element_cnt, & + element_idx, element_off, ipgf, jpgf + REAL(dp) :: AB(3), im_B(3), pgf_max, rab2 + + element_idx = 0 + element_off = 0 + cell_cnt = 0 + DO cell_idx = 1, SIZE(neighbor_cells) + ! move B to this cell while keeping A fixed + ! NOTE rb has already been moved by build_pair_list + ! so that when cell_r is (000) AB is in the pbc cell + IF (do_periodic) THEN + im_B = rb + neighbor_cells(cell_idx)%cell_r(:) + ELSE + im_B = rb + END IF + AB = ra - im_B + rab2 = AB(1)**2 + AB(2)**2 + AB(3)**2 + ! First screening on the most diffuse AB gaussians along with the most diffuse CD gaussian + IF (screen1(1)*rab2 + screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) CYCLE + + element_cnt = 0 + DO ipgf = 1, npgfa + DO jpgf = 1, npgfb + + ! Second screening on this actual AB pair and the most diffuse CD from before + pgf_max = pgf(jpgf, ipgf)%x(1)*rab2 + pgf(jpgf, ipgf)%x(2) + IF (pgf_max + screen2(2) + log10_pmax < log10_eps_schwarz) CYCLE + ! This pair passed screening. Add it to the bra. + element_idx = element_idx + 1 + bra%pgf_idx(:, element_idx) = [ipgf, jpgf, cell_idx] + bra%pgf_scr(1, element_idx) = pgf_max + element_cnt = element_cnt + 1 + END DO + END DO + + ! If this cell produced any pair, add this cell to the bra + IF (element_cnt == 0) CYCLE + cell_cnt = cell_cnt + 1 + bra%cell_idx(:, cell_cnt) = [cell_idx, element_cnt, element_off] + element_off = element_off + element_cnt + + END DO ! cell + bra%cell_cnt = cell_cnt + nelements = element_idx + bra%pgf_cnt = nelements + + END SUBROUTINE build_pair_list_pbc_pgf + ! ************************************************************************************************** !> \brief ... !> \param list1 ... diff --git a/src/hfx_types.F b/src/hfx_types.F index 11fa339dc6..1c11cc9d77 100644 --- a/src/hfx_types.F +++ b/src/hfx_types.F @@ -10,6 +10,7 @@ !> \par History !> 04.2008 created [Manuel Guidon] !> 05.2019 Moved erfc_cutoff to common/mathlib (A. Bussy) +!> 10.2025 Added gcc from basis_parameter and hfx_library option !> \author Manuel Guidon ! ************************************************************************************************** MODULE hfx_types @@ -199,6 +200,7 @@ MODULE hfx_types ! ************************************************************************************************** TYPE hfx_general_type REAL(dp) :: fraction = 0.0_dp !! for hybrids + INTEGER :: hfx_library = 0 LOGICAL :: treat_lsd_in_core = .FALSE. END TYPE hfx_general_type @@ -277,6 +279,7 @@ MODULE hfx_types INTEGER, DIMENSION(:), POINTER :: nshell => NULL() REAL(dp), DIMENSION(:, :, :, :), POINTER & :: sphi_ext => NULL() + REAL(dp), DIMENSION(:, :, :), POINTER :: gcc => NULL() REAL(dp), DIMENSION(:), POINTER :: set_radius => NULL() REAL(dp), DIMENSION(:, :), POINTER :: pgf_radius => NULL() REAL(dp) :: kind_radius = 0.0_dp @@ -669,6 +672,9 @@ CONTAINS CALL section_vals_val_get(hfx_section, "TREAT_LSD_IN_CORE", l_val=logic_val, i_rep_section=irep) actual_x_data%general_parameter%treat_lsd_in_core = logic_val + CALL section_vals_val_get(hfx_section, "HFX_LIBRARY", i_val=int_val, i_rep_section=irep) + actual_x_data%general_parameter%hfx_library = int_val + hfx_ri_section => section_vals_get_subs_vals(hfx_section, "RI") CALL section_vals_val_get(hfx_ri_section, "_SECTION_PARAMETERS_", l_val=actual_x_data%do_hfx_ri) @@ -1747,6 +1753,7 @@ CONTAINS nsgf_set=basis_parameter(ikind)%nsgf, & first_sgf=basis_parameter(ikind)%first_sgf, & sphi=basis_parameter(ikind)%sphi, & + gcc=basis_parameter(ikind)%gcc, & nsgf=basis_parameter(ikind)%nsgf_total, & l=basis_parameter(ikind)%nl, & nshell=basis_parameter(ikind)%nshell, & diff --git a/src/input_constants.F b/src/input_constants.F index fdceaa510e..093c401959 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -13,6 +13,7 @@ !> 10.2008 Teodoro Laino [tlaino] - University of Zurich !> Cleaned this file of all arrays of strings. Here must be kept !> only integer/real global constants +!> 102025 M Puligheddu : added an enum for hfx_library ! ************************************************************************************************** MODULE input_constants @@ -860,6 +861,10 @@ MODULE input_constants hfx_ri_do_2c_diag = 2, & hfx_ri_do_2c_cholesky = 3 + INTEGER, PARAMETER, PUBLIC :: hfx_library_is_libint = 1, & + hfx_library_is_libGint = 2, & + hfx_library_is_both = 12 + ! mode selctive vibrational analysis INTEGER, PARAMETER, PUBLIC :: ms_guess_bfgs = 1, & ms_guess_atomic = 2, & diff --git a/src/input_cp2k_hfx.F b/src/input_cp2k_hfx.F index f7398abb39..d2b7d2a4b8 100644 --- a/src/input_cp2k_hfx.F +++ b/src/input_cp2k_hfx.F @@ -23,7 +23,8 @@ MODULE input_cp2k_hfx do_potential_coulomb, do_potential_gaussian, do_potential_id, do_potential_long, & do_potential_mix_cl, do_potential_mix_cl_trunc, do_potential_mix_lg, do_potential_short, & do_potential_truncated, ehrenfest, gaussian, hfx_ri_do_2c_cholesky, hfx_ri_do_2c_diag, & - hfx_ri_do_2c_iter + hfx_ri_do_2c_iter, hfx_library_is_both, hfx_library_is_libGint, & + hfx_library_is_libint USE input_keyword_types, ONLY: keyword_create, & keyword_release, & keyword_type @@ -35,6 +36,7 @@ MODULE input_cp2k_hfx USE input_val_types, ONLY: real_t USE kinds, ONLY: dp USE string_utilities, ONLY: s2a + #include "./base/base_uses.f90" IMPLICIT NONE @@ -85,6 +87,23 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create( & + keyword, __LOCATION__, & + name="HFX_LIBRARY", & + description="Which library should be used in the calculation of the HF exchange "// & + "(Libint (cpu, default), libGint(gpu, cuda), both(debug,temporary)", & + usage="HFX_LIBRARY libGint", & + enum_c_vals=s2a("libint", "libGint", "both"), & + enum_i_vals=[hfx_library_is_libint, hfx_library_is_libGint, hfx_library_is_both], & + enum_desc=s2a("libint: use the libint library to compute the 2 electron integrals for HFX/r", & + "libGint: use the libGint library to accelerate the calculation of the HF "// & + "exchange on (cuda) GPUs /r", & + "both: temporary debug option, will run both libint and libGint and "// & + "check if the fock matrix is within tolerance"), & + default_i_val=hfx_library_is_libint) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="PW_HFX", & description="Compute the Hartree-Fock energy also in the plane wave basis. "// & "The value is ignored, and intended for debugging only.", & diff --git a/src/libgint_wrapper.F b/src/libgint_wrapper.F new file mode 100644 index 0000000000..b82732ec63 --- /dev/null +++ b/src/libgint_wrapper.F @@ -0,0 +1,696 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2026 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Interface to the LibGint-Library. +!> \par History +!> 10.2024 Created +!> \author Marcello Puligheddu +! ************************************************************************************************** +MODULE libGint_wrapper + + USE kinds, ONLY: dp +#if(__LIBGINT) + USE input_constants, ONLY: do_potential_coulomb, do_potential_truncated + USE libGint, ONLY: libgint_init, libgint_set_Potential_Truncated, libgint_set_hf_fac, libgint_set_max_mem, & + libgint_set_P, libgint_set_P_polarized, libgint_set_K, libgint_set_K_polarized, & + libgint_get_K, libgint_get_K_polarized, libgint_set_Atom, libgint_set_Atom_L, & + libgint_set_cell, libgint_set_neighs, & + libgint_add_prm, libgint_add_shell, libgint_add_cell, libgint_add_qrt, & + libgint_add_qrtt, libgint_add_set + USE t_c_g0, ONLY: C0 +#endif + USE hfx_types, ONLY: hfx_type, hfx_memory_type, hfx_potential_type, & + hfx_screen_coeff_type, hfx_cell_type, hfx_basis_type + + USE cell_types, ONLY: cell_type + USE hfx_pair_list_methods, ONLY: build_pair_list_pbc_pgf, bra_t, allocate_bra + USE particle_types, ONLY: particle_type + + USE iso_C_binding, ONLY: c_ptr + +#include "./base/base_uses.f90" + IMPLICIT NONE + PRIVATE + INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: first_set_of_atom + TYPE(bra_t), TARGET, SAVE :: bra, ket + LOGICAL, SAVE :: first_call = .TRUE. + TYPE(c_ptr), SAVE :: libGint_handle + !$OMP THREADPRIVATE( first_set_of_atom, first_call ) + !$OMP THREADPRIVATE( bra,ket ) + !$OMP THREADPRIVATE( libGint_handle ) + + PUBLIC :: cp_libGint_init, libGint_update_env, libGint_set_density, libGint_coulomb4, & + libGint_update_fock_matrix, libGint_get_fock_matrix + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'libGint_wrapper' + +! Comunicates the current density to the libGint engine. +! Sets the Fock matrix on the device to zero + INTERFACE libGint_set_density + MODULE PROCEDURE libGint_set_density_A + MODULE PROCEDURE libGint_set_density_AB + END INTERFACE + + INTERFACE libGint_get_fock_matrix + MODULE PROCEDURE libGint_get_fock_matrix_A + MODULE PROCEDURE libGint_get_fock_matrix_AB + END INTERFACE + +CONTAINS +! ************************************************************************************************** +!> \brief Sets libGint internal enviroment, must be called at least once before libGint can be used +!> \param[in] actual_x_data pointer to hfx_data +! ************************************************************************************************** + SUBROUTINE cp_libGint_init(actual_x_data) + + TYPE(hfx_type), POINTER, INTENT(in) :: actual_x_data +#if(__LIBGINT) + ! Init the offload library, creates an handle unique to this omp thread + CALL libgint_init(libGint_handle) + ! Comunicate the chosen potential and its parameters to libGint + IF (actual_x_data%potential_parameter%potential_type == do_potential_truncated) THEN + ! truncated coulomb needs the C0 coefficients. We do not read or compute them, + ! they must be already saved in C0 from the t[runcated]_c[oulomb]_g0 module + CALL libgint_set_Potential_Truncated(libGint_handle, & + actual_x_data%potential_parameter%cutoff_radius, & + C0(:, :)) + ELSE + CPABORT("The selected interaction potential type is not available with libGint") + END IF + first_call = .FALSE. +#else + MARK_USED(actual_x_data) + CPABORT("This CP2K executable has not been linked against the required library libGint.") +#endif + END SUBROUTINE cp_libGint_init +! ************************************************************************************************** +!> \brief Initialize and update the libGint computational environment. Must be called at least once +!> after geo_change and before libGint can be used. +!> +!> This routine sets up data required by the libGint integral engine, including +!> Hartree–Fock scaling factors, memory limits, periodic cell information, and +!> per-atom/basis-set data. It also allocates CPU-side buffers (`bra` and `ket`) +!> used for storing screened Gaussian primitive pairs. +!> +!> \param[in] fac Fraction of exact exchange +!> \param[in] memory_parameter Pointer to memory configuration +!> \param[in] do_periodic Logical flag: whether to consider pbc +!> \param[in] cell primitive simulation cell +!> \param[in] actual_x_data HF exchange data +!> \param[in] nneighbors Lattice points +!> \param[in] max_pgf Maximum number of primitive Gaussians per shel +!> \param[in] natom Number of atoms in the system +!> \param[in] kind_of Array mapping atom indices to atomic kinds +!> \param[in] particle_set particle set, we extract the positions +!> \param[in] basis_parameter gaussian basis set parameters +!> +! ************************************************************************************************** + SUBROUTINE libGint_update_env(fac, memory_parameter, do_periodic, cell, actual_x_data, nneighbors, max_pgf, & + natom, kind_of, particle_set, basis_parameter) + + REAL(dp) :: fac + TYPE(hfx_memory_type), POINTER :: memory_parameter + LOGICAL :: do_periodic + TYPE(cell_type), POINTER :: cell + TYPE(hfx_type), POINTER :: actual_x_data + INTEGER :: nneighbors, max_pgf, natom + INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter +#if(__LIBGINT) + LOGICAL(1) :: do_pbc + REAL(dp), DIMENSION(:, :), ALLOCATABLE :: cell_r + REAL(dp), DIMENSION(:), ALLOCATABLE :: flat_gcc + INTEGER, DIMENSION(:), POINTER :: la_min, la_max, npgfa, nsgfa + INTEGER, DIMENSION(:, :), POINTER :: nsgfl_a + INTEGER :: i, l, iset, jset, iatom, ikind, nseta, inla, nla + REAL(dp) :: ra(3) + REAL(dp), DIMENSION(:, :), POINTER :: zeta + REAL(dp), DIMENSION(:, :, :), POINTER :: gcc + TYPE(bra_t), POINTER :: bra_p, ket_p + + ! Set the multiplicative factor fac (the fraction of exact exchange times the spin factor) for libGint + CALL libGint_set_hf_fac(libGint_handle, fac) + ! Comunicate max gpu mem per mpi + CALL libGint_set_max_mem(libGint_handle, memory_parameter%max_memory) + + ! Info about periodic cells and neighbouring cells + do_pbc = do_periodic + ALLOCATE (cell_r(3, nneighbors)) + DO i = 1, nneighbors + cell_r(:, i) = actual_x_data%neighbor_cells(i)%cell_r(:) + END DO + CALL libgint_set_cell(libGint_handle, do_pbc, cell%hmat, cell%h_inv) + CALL libgint_set_neighs(libGint_handle, cell_r, nneighbors) + + ! CPU side temporary arrays for info about the pairs + bra_p => bra + ket_p => ket + CALL allocate_bra(bra_p, max_pgf, nneighbors) + CALL allocate_bra(ket_p, max_pgf, nneighbors) + + ! Comunicate atomset and atom info to the engine + jset = 1 + IF (ALLOCATED(first_set_of_atom)) DEALLOCATE (first_set_of_atom) + ALLOCATE (first_set_of_atom(natom)) + DO iatom = 1, natom + ikind = kind_of(iatom) + ra = particle_set(iatom)%r(:) + nseta = basis_parameter(ikind)%nset + npgfa => basis_parameter(ikind)%npgf + la_min => basis_parameter(ikind)%lmin + la_max => basis_parameter(ikind)%lmax + zeta => basis_parameter(ikind)%zet + nsgfa => basis_parameter(ikind)%nsgf + nsgfl_a => basis_parameter(ikind)%nsgfl + gcc => basis_parameter(ikind)%gcc + first_set_of_atom(iatom) = jset + DO iset = 1, nseta + CALL libgint_set_Atom(libGint_handle, jset - 1, ra, zeta(:, iset), npgfa(iset)) + inla = 1 + DO l = la_min(iset), la_max(iset) + nla = nsgfl_a(l, iset) + IF (ALLOCATED(flat_gcc)) DEALLOCATE (flat_gcc) + ALLOCATE (flat_gcc(npgfa(iset)*nla)) + flat_gcc = PACK(gcc(1:npgfa(iset), inla:inla + nla - 1, iset), .TRUE.) + CALL libgint_set_Atom_L(libGint_handle, jset - 1, l, nla, flat_gcc) + inla = inla + nla + END DO + jset = jset + 1 + END DO + END DO +#else + MARK_USED(fac) + MARK_USED(memory_parameter) + MARK_USED(do_periodic) + MARK_USED(cell) + MARK_USED(actual_x_data) + MARK_USED(nneighbors) + MARK_USED(max_pgf) + MARK_USED(natom) + MARK_USED(kind_of) + MARK_USED(particle_set) + MARK_USED(basis_parameter) + CPABORT("This CP2K executable has not been linked against the required library libGint.") +#endif + END SUBROUTINE libGint_update_env + +! ************************************************************************************************** +!> \brief communicates the density to libGint, no spin case +!> +!> \param[in] full_density_alpha full density matrix array +! ************************************************************************************************** + SUBROUTINE libGint_set_density_A(full_density_alpha) + + REAL(dp), DIMENSION(:, :), POINTER :: full_density_alpha +#if(__LIBGINT) + CALL libgint_set_P(libGint_handle, full_density_alpha(:, 1)) +#else + MARK_USED(full_density_alpha) + CPABORT("This CP2K executable has not been linked against the required library libGint.") +#endif + END SUBROUTINE libGint_set_density_A + +! ************************************************************************************************** +!> \brief communicates the density to libGint, spin case +!> +!> \param[in] full_density_alpha full density matrix array on the alpha channel +!> \param[in] full_density_beta full density matrix array on the beta channel +! ************************************************************************************************** + SUBROUTINE libGint_set_density_AB(full_density_alpha, full_density_beta) + + REAL(dp), DIMENSION(:, :), POINTER :: full_density_alpha, full_density_beta +#if(__LIBGINT) + CALL libgint_set_P_polarized(libGint_handle, full_density_alpha(:, 1), full_density_beta(:, 1)) +#else + MARK_USED(full_density_alpha) + MARK_USED(full_density_beta) + CPABORT("This CP2K executable has not been linked against the required library libGint.") +#endif + END SUBROUTINE libGint_set_density_AB + +! ************************************************************************************************** +!> \brief requests the Fock matrix ( hf_fraction * D @@ I ) from libGint. When this function +!> returns, full_ks_alpha_from_gpu will contain the fock matrix for this MPI rank. +!> +!> \param[out] full_ks_alpha_from_gpu location, already allocated, for the Fock matrix +!> \note assumes full_ks_alpha_from_gpu is already allocated with enough memory +! ************************************************************************************************** + SUBROUTINE libGint_get_fock_matrix_A(full_ks_alpha_from_gpu) + + REAL(dp), DIMENSION(:, :), POINTER :: full_ks_alpha_from_gpu +#if(__LIBGINT) + CALL libGint_get_K(libGint_handle, full_ks_alpha_from_gpu(:, 1)) +#else + MARK_USED(full_ks_alpha_from_gpu) + CPABORT("This CP2K executable has not been linked against the required library libGint.") +#endif + END SUBROUTINE libGint_get_fock_matrix_A + +! ************************************************************************************************** +!> \brief requests the Fock matrices from libGint for both channels. When this function +!> returns, full_ks_alpha_from_gpu and beta will contain the fock matrix for this MPI rank. +!> +!> \param[out] full_ks_alpha_from_gpu location, already allocated, for the Fock matrix alpha +!> \param[out] full_ks_beta_from_gpu location, already allocated, for the Fock matrix beta +!> \note assumes full_ks_alpha_from_gpu and full_ks_beta_from_gpu are already allocated +! ************************************************************************************************** + SUBROUTINE libGint_get_fock_matrix_AB(full_ks_alpha_from_gpu, full_ks_beta_from_gpu) + + REAL(dp), DIMENSION(:, :), POINTER :: full_ks_alpha_from_gpu, full_ks_beta_from_gpu +#if(__LIBGINT) + CALL libgint_get_K_polarized(libGint_handle, full_ks_alpha_from_gpu(:, 1), full_ks_beta_from_gpu(:, 1)) +#else + MARK_USED(full_ks_alpha_from_gpu) + MARK_USED(full_ks_beta_from_gpu) + CPABORT("This CP2K executable has not been linked against the required library libGint.") +#endif + END SUBROUTINE libGint_get_fock_matrix_AB + +! ************************************************************************************************** +!> \brief Assign two-electron integrals of a quartet/shell to libGint +!> +!> First we build a list of ab primitives and cd primitives using build_pair_list_pbc_pgf +!> These lists will contain (before screening) n_cell * npgf * npgf. +!> We loop over both and send ab(g) | cd(h) to libgint, which will loop internally over n +!> +!> \param[in] iatom index of atom I / A +!> \param[in] jatom index of atom J / B +!> \param[in] katom index of atom K / C +!> \param[in] latom index of atom L / D +!> \param[in] iset index of the set for atom iatom we are adding +!> \param[in] jset index of the set for atom jatom we are adding +!> \param[in] kset index of the set for atom katom we are adding +!> \param[in] lset index of the set for atom latom we are adding +!> \param[in] ra position of atom iatom +!> \param[in] rb position of atom jatom +!> \param[in] rc position of atom katom +!> \param[in] rd position of atom latom +!> \param[in] npgfa number of gaussians in set iset of atom iatom +!> \param[in] npgfb number of gaussians in set jset of atom jatom +!> \param[in] npgfc number of gaussians in set kset of atom katom +!> \param[in] npgfd number of gaussians in set lset of atom latom +!> \param[in] potential_parameter information about the potential +!> \param[in] screen1 screening information for AB pair +!> \param[in] screen2 screening information for CD pair +!> \param[in] log10_pmax density screening factor +!> \param[in] log10_eps_schwarz screening tolerance +!> \param[in] pgf1 screening information for each AB pair primitive +!> \param[in] pgf2 screening information for each CD pair primitive +!> \param[in] neighbor_cells array with lattice vectors +!> \param[in] cell information about the simulation box +!> \param[in] do_periodic flag for pbc +!> \param[out] screened whether the whole quartet was screened out +! ************************************************************************************************** + + SUBROUTINE libGint_coulomb4(iatom, jatom, katom, latom, iset, jset, kset, lset, & + ra, rb, rc, rd, npgfa, npgfb, npgfc, npgfd, & + potential_parameter, & + screen1, screen2, log10_pmax, log10_eps_schwarz, & + pgf1, pgf2, & + neighbor_cells, cell, do_periodic, screened) + + INTEGER, INTENT(in) :: iatom, jatom, katom, latom, iset, jset, kset, lset + REAL(dp), INTENT(IN) :: ra(3), rb(3), rc(3), rd(3) + INTEGER, INTENT(IN) :: npgfa, npgfb, npgfc, npgfd + TYPE(hfx_potential_type) :: potential_parameter + REAL(dp), INTENT(IN) :: screen1(2), screen2(2) + REAL(dp), INTENT(IN) :: log10_pmax, log10_eps_schwarz + TYPE(hfx_screen_coeff_type), DIMENSION(:, :), & + POINTER :: pgf1, pgf2 + TYPE(hfx_cell_type), DIMENSION(:), POINTER :: neighbor_cells + TYPE(cell_type), POINTER :: cell + LOGICAL, INTENT(IN) :: do_periodic + LOGICAL, INTENT(out) :: screened + +#if(__LIBGINT) + TYPE(bra_t), POINTER :: bra_p, ket_p + LOGICAL :: cell_was_screened + INTEGER :: idx_n1, idx_n2, n1, n2, idx_ij, idx_kl, o_ij, n_ij, o_kl, n_kl !, n3 + INTEGER :: ipgf, jpgf, kpgf, lpgf, iatom_set, jatom_set, katom_set, latom_set + REAL(dp) :: pgf_max_1, pgf_max_2 ! , R1, R2, rpq2 + INTEGER :: nelements_ij, nelements_kl + + cell = cell + potential_parameter = potential_parameter + bra_p => bra + ket_p => ket + + screened = .TRUE. + iatom_set = first_set_of_atom(iatom) + iset - 2 + jatom_set = first_set_of_atom(jatom) + jset - 2 + katom_set = first_set_of_atom(katom) + kset - 2 + latom_set = first_set_of_atom(latom) + lset - 2 + + CALL build_pair_list_pbc_pgf(npgfa, npgfb, bra_p, screen1, screen2, & + pgf1, log10_pmax, log10_eps_schwarz, ra, rb, & + nelements_ij, neighbor_cells, do_periodic) + + CALL build_pair_list_pbc_pgf(npgfc, npgfd, ket_p, screen2, screen1, & + pgf2, log10_pmax, log10_eps_schwarz, rc, rd, & + nelements_kl, neighbor_cells, do_periodic) + + ! Note: we use 3 numbers n1 n2 and n3 as indices for the lattice traslantion vectors + ! n1 for the AB pair, n2 for the CD pair and n3 for the PQ pair + ! so that e.g. B = B0 + ua(n1) means B.x = B0.x + ua(n1).x (and same for y and z) + ! the ua, saved in neighbor_cells(:)%cell_r(:), are already computed + ! lattice translation vectors T = i a1 + j a2 + k a3 + ! where a1,a2 and a3 are lattice vectors and i,j and k integers. + ! So, B.x = B0.x + n1.i * a1.x + n1.j * a2.x + n1.k * a3.x (and same for y and z) + ! + + DO idx_n1 = 1, bra%cell_cnt + + n1 = bra%cell_idx(1, idx_n1) + n_ij = bra%cell_idx(2, idx_n1) + o_ij = bra%cell_idx(3, idx_n1) + + DO idx_n2 = 1, ket%cell_cnt + + n2 = ket%cell_idx(1, idx_n2) + n_kl = ket%cell_idx(2, idx_n2) + o_kl = ket%cell_idx(3, idx_n2) + + cell_was_screened = .TRUE. + DO idx_ij = o_ij + 1, o_ij + n_ij + + ipgf = bra%pgf_idx(1, idx_ij) + jpgf = bra%pgf_idx(2, idx_ij) + + pgf_max_1 = bra%pgf_scr(1, idx_ij) + + DO idx_kl = o_kl + 1, o_kl + n_kl + kpgf = ket%pgf_idx(1, idx_kl) + lpgf = ket%pgf_idx(2, idx_kl) + + pgf_max_2 = ket%pgf_scr(1, idx_kl) + + IF (pgf_max_1 + pgf_max_2 + log10_pmax < log10_eps_schwarz) CYCLE + + CALL libGint_add_prm(libGint_handle, ipgf - 1, jpgf - 1, kpgf - 1, lpgf - 1) + cell_was_screened = .FALSE. + + END DO ! ket pgf + END DO ! bra pgf + + IF (.NOT. cell_was_screened) THEN + CALL libgint_add_shell(libGint_handle, iatom_set, jatom_set, katom_set, latom_set, n1 - 1, n2 - 1) + screened = .FALSE. + END IF + + END DO ! ket n2 + END DO ! bra n1 +#else + MARK_USED(iatom) + MARK_USED(jatom) + MARK_USED(katom) + MARK_USED(latom) + MARK_USED(iset) + MARK_USED(jset) + MARK_USED(kset) + MARK_USED(lset) + MARK_USED(ra) + MARK_USED(rb) + MARK_USED(rc) + MARK_USED(rd) + MARK_USED(npgfa) + MARK_USED(npgfb) + MARK_USED(npgfc) + MARK_USED(npgfd) + MARK_USED(potential_parameter) + MARK_USED(screen1) + MARK_USED(screen2) + MARK_USED(log10_pmax) + MARK_USED(log10_eps_schwarz) + MARK_USED(pgf1) + MARK_USED(pgf2) + MARK_USED(neighbor_cells) + MARK_USED(cell) + MARK_USED(do_periodic) + MARK_USED(screened) + CPABORT("This CP2K executable has not been linked against the required library libGint.") +#endif + END SUBROUTINE libGint_coulomb4 + +! ************************************************************************************************** +!> \brief The previous coulomb_4 function assigned an integral between primitive gaussian. +!> Now we assign the nsgfa(iset)*b*c*d gcc integrals form this set, +!> along with the Pac Pad Pbc Pbd density to the Kbd Kbc Kad Kac Fock matrix +!> +!> \param[in] symm_fac simmetry factor from iatom=jatom and the like +!> \param[in] iatom index of atom A +!> \param[in] jatom index of atom B +!> \param[in] katom index of atom C +!> \param[in] latom index of atom D +!> \param[in] iset index of set for atom A +!> \param[in] jset index of set for atom B +!> \param[in] kset index of set for atom C +!> \param[in] lset index of set for atom D +!> \param[in] atomic_offset_ac global offset for the pair of A and C atom in the density ( and Fock) matrix +!> \param[in] atomic_offset_ad global offset for the pair of A and D atom in the density ( and Fock) matrix +!> \param[in] atomic_offset_bc global offset for the pair of B and C atom in the density ( and Fock) matrix +!> \param[in] atomic_offset_bd global offset for the pair of B and D atom in the density ( and Fock) matrix +!> \param[in] offset_ac_set matrix of sub_offset for sets in atomic_offset_ac +!> \param[in] offset_ad_set matrix of sub_offset for sets in atomic_offset_ad +!> \param[in] offset_bc_set matrix of sub_offset for sets in atomic_offset_bc +!> \param[in] offset_bd_set matrix of sub_offset for sets in atomic_offset_bd +!> \param[in] nsgfa total number of (spherical, contracted) integrals for iset. +!> Used as leading dimension of the AC, AD sublock, depending on transposition +!> \param[in] nsgfb total number of (spherical, contracted) integrals for jset. +!> Used as leading dimension of the BC, BD sublock, depending on transposition +!> \param[in] nsgfc total number of (spherical, contracted) integrals for kset. +!> Used as leading dimension of the AC, BC sublock, depending on transposition +!> \param[in] nsgfd total number of (spherical, contracted) integrals for lset. +!> Used as leading dimension of the AD, BD sublock, depending on transposition +!> \param[in] la_min minumum angular moment in set iset of atom A +!> \param[in] la_max maximum angular moment in set iset of atom A +!> \param[in] lb_min minumum angular moment in set iset of atom B +!> \param[in] lb_max maximum angular moment in set iset of atom B +!> \param[in] lc_min minumum angular moment in set iset of atom C +!> \param[in] lc_max maximum angular moment in set iset of atom C +!> \param[in] ld_min minumum angular moment in set iset of atom D +!> \param[in] ld_max maximum angular moment in set iset of atom D +!> \param[in] nsgfl_a matrix with the number of linear combinations of primitive gaussians +!> for each angular moment for each set in atom A. We only read iset +!> \param[in] nsgfl_b matrix with the number of linear combinations of primitive gaussians +!> for each angular moment for each set in atom B. We only read jset +!> \param[in] nsgfl_c matrix with the number of linear combinations of primitive gaussians +!> for each angular moment for each set in atom C. We only read kset +!> \param[in] nsgfl_d matrix with the number of linear combinations of primitive gaussians +!> for each angular moment for each set in atom D. We only read lset +!> \note +!> The atomic_offset_xy, offset_xy_set matrices provide set offsets which are combined with per-L and +!> per linear combination offsets to produce the final indices into the dense (but block-sparse) +!> density and Fock matrices +!> - The routine assumes Fortran column major order and contiguous storage for the per set +!> density subblocks as described in the code comments. +!> - The code computes transposition flags to only use the lower part of P and K +!> +!> +! ************************************************************************************************** + + SUBROUTINE libGint_update_fock_matrix( & + symm_fac, & + iatom, jatom, katom, latom, & + iset, jset, kset, lset, & + atomic_offset_ac, atomic_offset_ad, atomic_offset_bc, atomic_offset_bd, & + offset_ac_set, offset_ad_set, offset_bc_set, offset_bd_set, & + nsgfa, nsgfb, nsgfc, nsgfd, & + la_min, la_max, lb_min, lb_max, & + lc_min, lc_max, ld_min, ld_max, & + nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d) + + REAL(dp) :: symm_fac + INTEGER :: iatom, jatom, katom, latom + INTEGER :: iset, jset, kset, lset + INTEGER :: atomic_offset_ac, atomic_offset_ad, atomic_offset_bc, atomic_offset_bd + INTEGER, DIMENSION(:, :), POINTER :: offset_ac_set, offset_ad_set + INTEGER, DIMENSION(:, :), POINTER :: offset_bc_set, offset_bd_set + INTEGER :: nsgfa, nsgfb, nsgfc, nsgfd + INTEGER :: la_min, la_max, lb_min, lb_max + INTEGER :: lc_min, lc_max, ld_min, ld_max + INTEGER, DIMENSION(:, :), POINTER :: nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d + +#if(__LIBGINT) + !! (Hyp) + ! Let a be a set composed of 2 s and 1 p function. + ! Let c be a set composed of 1 s and 2 p function. + ! (1) The density matrix for the ac pair is a 5 x 7 matrix organized as + ! + ! / -------------------------------------------------------------------------------------------------------------------\ + ! | a_s1_0@c_s1_0 || a_s1_0@c_p1_0 | a_s1_0@c_p1_1 | a_s1_0@c_p1_2 || a_s1_0@c_p2_0 | a_s1_0@c_p2_1 | a_s1_0@c_p2_2 | + ! | a_s2_0@c_s1_0 || a_s2_0@c_p1_0 | a_s2_0@c_p1_1 | a_s2_0@c_p1_2 || a_s2_0@c_p2_0 | a_s2_0@c_p2_1 | a_s2_0@c_p2_2 | + ! | a_p1_0@c_s1_0 || a_p1_0@c_p1_0 | a_p1_0@c_p1_1 | a_p1_0@c_p1_2 || a_p1_0@c_p2_0 | a_p1_0@c_p2_1 | a_p1_0@c_p2_2 | + ! | a_p1_1@c_s1_0 || a_p1_1@c_p1_0 | a_p1_1@c_p1_1 | a_p1_1@c_p1_2 || a_p1_1@c_p2_0 | a_p1_1@c_p2_1 | a_p1_1@c_p2_2 | + ! | a_p1_2@c_s1_0 || a_p1_2@c_p1_0 | a_p1_2@c_p1_1 | a_p1_2@c_p1_2 || a_p1_2@c_p2_0 | a_p1_2@c_p2_1 | a_p1_2@c_p2_2 | + ! \ -------------------------------------------------------------------------------------------------------------------/ + ! + ! where A_LX_Y means the (Y+1) component of the Xth linear combination of the L angular moment for atom A + ! + ! (2) This matrix is dense, rectangular and contigous in memory, in fortran column major order. + ! (3) The big matrix with all pairs is block sparse triangular, only the lower part is valid. + + LOGICAL(1) :: Tac, Tad, Tbc, Tbd + INTEGER :: offset_ac_L_set, offset_ad_L_set, offset_bc_L_set, offset_bd_L_set + INTEGER :: s_offset_a, s_offset_b, s_offset_c, s_offset_d + INTEGER :: s_offset_a_l, s_offset_b_l, s_offset_c_l, s_offset_d_l + INTEGER :: s_offset_ac, s_offset_ad, s_offset_bc, s_offset_bd + INTEGER :: ld_ac_set, ld_ad_set, ld_bc_set, ld_bd_set + INTEGER :: la, lb, lc, ld, nla, nlb, nlc, nld, inla, inlb, inlc, inld + ! TODO rewrite as update_fock_matrix_gpu(libGint_handle, iatomset,jatomset,katomset,latomset ) + ! AFTER TODO communicate (the pointer to) atomic_offset to libGint AND + ! AFTER TODO communicate (the pointer to) set_offset to libGint AND + ! AFTER TODO check if this idea makes sense in general for other codes + ! + ! Note: this would not change the need to compute sub offsets and + ! the 8 loops, it would just transfer them to libGint + ! Except, if libGint can be sure every set has 1 l, it can collapse the l loops + ! and/or, if libGint can be sure every l has 1 nl, it can collapse the n loops + IF (jatom >= latom) THEN + offset_bd_L_set = offset_bd_set(jset, lset) + atomic_offset_bd - 2 + ld_bd_set = nsgfb + Tbd = .FALSE. + ELSE + offset_bd_L_set = offset_bd_set(lset, jset) + atomic_offset_bd - 2 + ld_bd_set = nsgfd + Tbd = .TRUE. + END IF + IF (jatom >= katom) THEN + offset_bc_L_set = offset_bc_set(jset, kset) + atomic_offset_bc - 2 + ld_bc_set = nsgfb + Tbc = .FALSE. + ELSE + offset_bc_L_set = offset_bc_set(kset, jset) + atomic_offset_bc - 2 + ld_bc_set = nsgfc + Tbc = .TRUE. + END IF + + IF (iatom >= latom) THEN + offset_ad_L_set = offset_ad_set(iset, lset) + atomic_offset_ad - 2 + ld_ad_set = nsgfa + Tad = .FALSE. + ELSE + offset_ad_L_set = offset_ad_set(lset, iset) + atomic_offset_ad - 2 + ld_ad_set = nsgfd + Tad = .TRUE. + END IF + + IF (iatom >= katom) THEN + offset_ac_L_set = offset_ac_set(iset, kset) + atomic_offset_ac - 2 + ld_ac_set = nsgfa + Tac = .FALSE. + ELSE + offset_ac_L_set = offset_ac_set(kset, iset) + atomic_offset_ac - 2 + ld_ac_set = nsgfc + Tac = .TRUE. + END IF + + s_offset_a_l = 0 + DO la = la_min, la_max + nla = nsgfl_a(la, iset) + s_offset_b_l = 0 + DO lb = lb_min, lb_max + nlb = nsgfl_b(lb, jset) + s_offset_c_l = 0 + DO lc = lc_min, lc_max + nlc = nsgfl_c(lc, kset) + s_offset_d_l = 0 + ld_loop: DO ld = ld_min, ld_max + nld = nsgfl_d(ld, lset) + CALL libgint_add_qrt(libGint_handle, la, lb, lc, ld, nla, nlb, nlc, nld) + DO inla = 1, nla + s_offset_a = s_offset_a_l + (inla - 1)*(2*la + 1) + DO inlb = 1, nlb + s_offset_b = s_offset_b_l + (inlb - 1)*(2*lb + 1) + DO inlc = 1, nlc + s_offset_c = s_offset_c_l + (inlc - 1)*(2*lc + 1) + DO inld = 1, nld + s_offset_d = s_offset_d_l + (inld - 1)*(2*ld + 1) + IF (.NOT. Tac) THEN + s_offset_ac = offset_ac_L_set + s_offset_c*ld_ac_set + s_offset_a + ELSE + s_offset_ac = offset_ac_L_set + s_offset_a*ld_ac_set + s_offset_c + END IF + + IF (.NOT. Tad) THEN + s_offset_ad = offset_ad_L_set + s_offset_d*ld_ad_set + s_offset_a + ELSE + s_offset_ad = offset_ad_L_set + s_offset_a*ld_ad_set + s_offset_d + END IF + + IF (.NOT. Tbc) THEN + s_offset_bc = offset_bc_L_set + s_offset_c*ld_bc_set + s_offset_b + ELSE + s_offset_bc = offset_bc_L_set + s_offset_b*ld_bc_set + s_offset_c + END IF + + IF (.NOT. Tbd) THEN + s_offset_bd = offset_bd_L_set + s_offset_d*ld_bd_set + s_offset_b + ELSE + s_offset_bd = offset_bd_L_set + s_offset_b*ld_bd_set + s_offset_d + END IF + + CALL libgint_add_qrtt(libGint_handle, symm_fac, & + la, lb, lc, ld, inla - 1, inlb - 1, inlc - 1, inld - 1, & + ld_ac_set, ld_ad_set, ld_bc_set, ld_bd_set, & + s_offset_ac, s_offset_ad, s_offset_bc, s_offset_bd, & + Tac, Tad, Tbc, Tbd) + + END DO + END DO + END DO + END DO + s_offset_d_l = s_offset_d_l + nld*(2*ld + 1) + END DO ld_loop + s_offset_c_l = s_offset_c_l + nlc*(2*lc + 1) + END DO + s_offset_b_l = s_offset_b_l + nlb*(2*lb + 1) + END DO + s_offset_a_l = s_offset_a_l + nla*(2*la + 1) + END DO + + CALL libgint_add_set(libGint_handle) + +#else + MARK_USED(symm_fac) + MARK_USED(iatom) + MARK_USED(jatom) + MARK_USED(katom) + MARK_USED(latom) + MARK_USED(iset) + MARK_USED(jset) + MARK_USED(kset) + MARK_USED(lset) + MARK_USED(atomic_offset_ac) + MARK_USED(atomic_offset_ad) + MARK_USED(atomic_offset_bc) + MARK_USED(atomic_offset_bd) + MARK_USED(offset_ac_set) + MARK_USED(offset_ad_set) + MARK_USED(offset_bc_set) + MARK_USED(offset_bd_set) + MARK_USED(nsgfa) + MARK_USED(nsgfb) + MARK_USED(nsgfc) + MARK_USED(nsgfd) + MARK_USED(la_min) + MARK_USED(la_max) + MARK_USED(lb_min) + MARK_USED(lb_max) + MARK_USED(lc_min) + MARK_USED(lc_max) + MARK_USED(ld_min) + MARK_USED(ld_max) + MARK_USED(nsgfl_a) + MARK_USED(nsgfl_b) + MARK_USED(nsgfl_c) + MARK_USED(nsgfl_d) + CPABORT("This CP2K executable has not been linked against the required library libGint.") +#endif + END SUBROUTINE libGint_update_fock_matrix + +END MODULE libGint_wrapper diff --git a/src/offload/offload_library.c b/src/offload/offload_library.c index e8cb6c9317..7b2de2c16b 100644 --- a/src/offload/offload_library.c +++ b/src/offload/offload_library.c @@ -23,7 +23,7 @@ #if defined(__OFFLOAD_PROFILING) #if defined(__OFFLOAD_CUDA) -#include +#include #elif defined(__OFFLOAD_HIP) && defined(__HIP_PLATFORM_AMD__) #include #endif diff --git a/tests/QS/regtest-libgint/H1.inp b/tests/QS/regtest-libgint/H1.inp new file mode 100644 index 0000000000..76bc888c6e --- /dev/null +++ b/tests/QS/regtest-libgint/H1.inp @@ -0,0 +1,67 @@ +&GLOBAL + EPS_CHECK_DIAG 1.0E-12 + PRINT_LEVEL Medium + PROJECT H2O-32 + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + METHOD QS + &DFT + BASIS_SET_FILE_NAME ./basis.dat + POTENTIAL_FILE_NAME POTENTIAL + &MGRID + CUTOFF 140 + REL_CUTOFF 10 + &END MGRID + &QS + EPS_DEFAULT 1.0E-12 + EXTRAPOLATION_ORDER 3 + WF_INTERPOLATION PS + &END QS + &SCF + EPS_SCF 1.0E-12 + MAX_SCF 400 + SCF_GUESS ATOMIC + &OT ON + MINIMIZER DIIS + &END OT + &PRINT + &RESTART OFF + &END RESTART + &END PRINT + &END SCF + &XC + &HF + HFX_LIBRARY both + &INTERACTION_POTENTIAL + CUTOFF_RADIUS 4 + POTENTIAL_TYPE TRUNCATED + T_C_G_DATA ./t_c_g.dat + &END INTERACTION_POTENTIAL + &MEMORY + EPS_STORAGE_SCALING 0.1 + MAX_MEMORY 1000 + &END MEMORY + &SCREENING + EPS_SCHWARZ 1.0E-12 + &END SCREENING + &END HF + &XC_FUNCTIONAL NONE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 10. 10. 10. + &END CELL + &COORD + H 0. 0. 0. + H 1. 0. 0. + &END COORD + &KIND H + BASIS_SET H_1 + POTENTIAL GTH-PADE-q1 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-libgint/H2O-HFX-6.inp b/tests/QS/regtest-libgint/H2O-HFX-6.inp new file mode 100644 index 0000000000..9e4bb52af2 --- /dev/null +++ b/tests/QS/regtest-libgint/H2O-HFX-6.inp @@ -0,0 +1,89 @@ +&GLOBAL + EPS_CHECK_DIAG 1.0E-8 + PRINT_LEVEL Medium + PROJECT H2O-6_low_prec + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + METHOD QS + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME POTENTIAL + &MGRID + CUTOFF 80 + REL_CUTOFF 40 + &END MGRID + &QS + EPS_DEFAULT 1.0E-8 + EXTRAPOLATION_ORDER 3 + WF_INTERPOLATION PS + &END QS + &SCF + EPS_SCF 1.0E-8 + IGNORE_CONVERGENCE_FAILURE + MAX_SCF 3 + SCF_GUESS ATOMIC + &OT ON + MINIMIZER DIIS + &END OT + &PRINT + &RESTART OFF + &END RESTART + &END PRINT + &END SCF + &XC + &HF + HFX_LIBRARY libgint + &INTERACTION_POTENTIAL + CUTOFF_RADIUS 4.5 + POTENTIAL_TYPE TRUNCATED + T_C_G_DATA ./t_c_g.dat + &END INTERACTION_POTENTIAL + &MEMORY + EPS_STORAGE_SCALING 0.1 + MAX_MEMORY 1000 + &END MEMORY + &SCREENING + EPS_SCHWARZ 1.0E-6 + &END SCREENING + &END HF + &XC_FUNCTIONAL NONE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 9.8528 9.8528 9.8528 + &END CELL + # 3 H2O cut from 32 H2O (TIP5P,1bar,300K) a = 9.8528 + &COORD + O 2.280398 9.146539 5.088696 + H 1.762019 9.820429 5.528454 + H 3.095987 9.107088 5.588186 + O 1.251703 2.406261 7.769908 + H 0.554129 2.982634 8.082024 + H 1.771257 2.954779 7.182181 + O 1.596302 6.920128 0.656695 + H 2.112148 6.126321 0.798136 + H 1.776389 7.463264 1.424030 + O 3.280398 10.146539 6.088696 + H 2.762019 10.820429 6.528454 + H 4.095987 10.107088 6.588186 + O 2.251703 3.406261 8.769908 + H 1.554129 3.982634 9.082024 + H 2.771257 3.954779 8.182181 + O 2.596302 7.920128 1.656695 + H 3.112148 7.126321 1.798136 + H 2.776389 8.463264 2.424030 + &END COORD + &KIND H + BASIS_SET DZVP-GTH + POTENTIAL GTH-PADE-q1 + &END KIND + &KIND O + BASIS_SET DZVP-GTH + POTENTIAL GTH-PADE-q6 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-libgint/TEST_FILES.toml b/tests/QS/regtest-libgint/TEST_FILES.toml new file mode 100644 index 0000000000..c0c455af4f --- /dev/null +++ b/tests/QS/regtest-libgint/TEST_FILES.toml @@ -0,0 +1,4 @@ +# The regtests of libGint +"H1.inp" = [{matcher="E_total", tol=3e-10, ref=-0.370379995347038}] +"H2O-HFX-6.inp" = [{matcher="E_total", tol=1e-1, ref=-98.48394326}] +#EOF diff --git a/tests/QS/regtest-libgint/basis.dat b/tests/QS/regtest-libgint/basis.dat new file mode 100644 index 0000000000..211fc65dfc --- /dev/null +++ b/tests/QS/regtest-libgint/basis.dat @@ -0,0 +1,10 @@ +# Hydrogen s + g single gaussian each +H H_1 + 2 + 1 0 0 1 1 + 1.0 1.0 + 1 1 1 1 1 + 1.0 1.0 + + +# diff --git a/tests/QS/regtest-libgint/tidy b/tests/QS/regtest-libgint/tidy new file mode 100644 index 0000000000..611f1ade63 --- /dev/null +++ b/tests/QS/regtest-libgint/tidy @@ -0,0 +1,18 @@ +#!/bin/tcsh +# Script to clean files after execution of a test +echo "This script will clean the directory" $PWD "from scratch files.." +echo +foreach file (`cvs update -A -P|grep '? '|awk '{print $2}'`) + if ( ${file:e} == 'inp' ) then + echo "This files could be a new input file you don't want to be deleted.. Do you really want it deleted?" + if ( $1 == "-int" ) then + echo "You requested to run the cleaning script interactively! Forget about complains!" + rm -f $file + else + rm -i $file + endif + else + rm -f $file + endif +end +rm -f *.bak diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 1a0e56c8f3..cf51c92537 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -10,6 +10,7 @@ QS/regtest-dft-d3-auto-ref libint QS/regtest-dft-d3-bj-auto-ref libint QS/regtest-dft-d3-lib-ref libint s_dftd3 QS/regtest-dft-d3bj-lib-ref libint s_dftd3 +QS/regtest-libgint libGint libint offload_cuda QS/regtest-pcc QS/regtest-new-pulay-mixing QS/regtest-gcscf diff --git a/tools/docker/Dockerfile.test_cuda_A100 b/tools/docker/Dockerfile.test_cuda_A100 index f2054502d0..99d11ab1e0 100644 --- a/tools/docker/Dockerfile.test_cuda_A100 +++ b/tools/docker/Dockerfile.test_cuda_A100 @@ -36,6 +36,7 @@ RUN ./install_cp2k_toolchain.sh \ --with-mpich=install \ --mpi-mode=mpich \ --enable-cuda=yes \ + --with-libgint=install \ --with-sirius=install \ --gpu-ver=A100 \ --dry-run diff --git a/tools/docker/Dockerfile.test_cuda_P100 b/tools/docker/Dockerfile.test_cuda_P100 index 5d41c5e35d..c977932697 100644 --- a/tools/docker/Dockerfile.test_cuda_P100 +++ b/tools/docker/Dockerfile.test_cuda_P100 @@ -36,6 +36,7 @@ RUN ./install_cp2k_toolchain.sh \ --with-mpich=install \ --mpi-mode=mpich \ --enable-cuda=yes \ + --with-libgint=install \ --with-sirius=install \ --gpu-ver=P100 \ --dry-run diff --git a/tools/docker/Dockerfile.test_cuda_V100 b/tools/docker/Dockerfile.test_cuda_V100 index 2c5b1c4421..6c63984131 100644 --- a/tools/docker/Dockerfile.test_cuda_V100 +++ b/tools/docker/Dockerfile.test_cuda_V100 @@ -36,6 +36,7 @@ RUN ./install_cp2k_toolchain.sh \ --with-mpich=install \ --mpi-mode=mpich \ --enable-cuda=yes \ + --with-libgint=install \ --with-sirius=install \ --gpu-ver=V100 \ --dry-run diff --git a/tools/docker/Dockerfile.test_performance_cuda_A100 b/tools/docker/Dockerfile.test_performance_cuda_A100 index 45e34422e5..d23a77bc63 100644 --- a/tools/docker/Dockerfile.test_performance_cuda_A100 +++ b/tools/docker/Dockerfile.test_performance_cuda_A100 @@ -36,6 +36,7 @@ RUN ./install_cp2k_toolchain.sh \ --with-mpich=install \ --mpi-mode=mpich \ --enable-cuda=yes \ + --with-libgint=install \ --with-sirius=install \ --gpu-ver=A100 \ --dry-run diff --git a/tools/docker/Dockerfile.test_performance_cuda_P100 b/tools/docker/Dockerfile.test_performance_cuda_P100 index 8f2d583a96..c62690acef 100644 --- a/tools/docker/Dockerfile.test_performance_cuda_P100 +++ b/tools/docker/Dockerfile.test_performance_cuda_P100 @@ -36,6 +36,7 @@ RUN ./install_cp2k_toolchain.sh \ --with-mpich=install \ --mpi-mode=mpich \ --enable-cuda=yes \ + --with-libgint=install \ --with-sirius=install \ --gpu-ver=P100 \ --dry-run diff --git a/tools/docker/Dockerfile.test_performance_cuda_V100 b/tools/docker/Dockerfile.test_performance_cuda_V100 index 7a1ea52a5b..6c3d8c59b3 100644 --- a/tools/docker/Dockerfile.test_performance_cuda_V100 +++ b/tools/docker/Dockerfile.test_performance_cuda_V100 @@ -36,6 +36,7 @@ RUN ./install_cp2k_toolchain.sh \ --with-mpich=install \ --mpi-mode=mpich \ --enable-cuda=yes \ + --with-libgint=install \ --with-sirius=install \ --gpu-ver=V100 \ --dry-run diff --git a/tools/docker/generate_dockerfiles.py b/tools/docker/generate_dockerfiles.py index 6d0d149b54..bfffdc4507 100755 --- a/tools/docker/generate_dockerfiles.py +++ b/tools/docker/generate_dockerfiles.py @@ -577,6 +577,7 @@ RUN apt-get update -qq && apt-get install -qq --no-install-recommends \ with_mpich="install", mpi_mode="mpich", enable_cuda="yes", + with_libgint="install", with_sirius="install", gpu_ver=gpu_ver, **kwargs, diff --git a/tools/toolchain/install_cp2k_toolchain.sh b/tools/toolchain/install_cp2k_toolchain.sh index a0e5e03840..7be7f09cb7 100755 --- a/tools/toolchain/install_cp2k_toolchain.sh +++ b/tools/toolchain/install_cp2k_toolchain.sh @@ -361,6 +361,8 @@ Specific options of --with-PKG: --with-cusolvermp NVIDIA cusolverMp: CUDA library for distributed dense linear algebra. Default = no + --with-libgint Enable the use of libGint for the calculation of the Hartree-Fock exchange on (nvidia) GPUs + Default = no FURTHER INSTRUCTIONS @@ -479,7 +481,7 @@ math_list="mkl acml openblas" lib_list="fftw eigen libint libxc gauxc libxsmm libxs libxstream cosma scalapack elpa dbcsr cusolvermp plumed spfft spla gsl spglib hdf5 libvdwxc sirius libvori libtorch deepmd ace dftd4 tblite pugixml libsmeagol fmt trexio - libfci greenx gmp mcl" + libfci greenx gmp mcl libgint" package_list="${tool_list} ${mpi_list} ${math_list} ${lib_list}" # ------------------------------------------------------------------------ @@ -526,7 +528,7 @@ with_dftd4="__DONTUSE__" with_tblite="__INSTALL__" with_libsmeagol="__DONTUSE__" with_mcl="__DONTUSE__" - +with_libgint="__DONTUSE__" # the math and mpi libraries are chosen by their respective modes. # default math library settings, MATH_MODE picks the math library # to use, and with_* defines the default method of installation if it @@ -898,6 +900,9 @@ Otherwise use option no." --with-gsl*) with_gsl=$(read_with "${1}") ;; + --with-libgint*) + with_libgint=$(read_with "${1}") + ;; --with-fmt*) with_fmt=$(read_with "${1}") ;; @@ -1165,6 +1170,12 @@ CMake is not enabled, so a new copy of CMake will be installed first." with_cmake="__INSTALL__" fi +#libGint installation requires cuda enabled +if [ "${with_libgint}" != "__DONTUSE__" ] && [ "${enable_cuda}" != "__TRUE__" ]; then + report_warning ${LINENO} "libGint requires the use of cuda. Disabling libGint" + with_libgint="__DONTUSE__" +fi + # SIRIUS dependencies if [ "${with_sirius}" = "__INSTALL__" ]; then [ "${with_spfft}" = "__DONTUSE__" ] && with_spfft="__INSTALL__" diff --git a/tools/toolchain/scripts/stage4/install_libGint.sh b/tools/toolchain/scripts/stage4/install_libGint.sh new file mode 100755 index 0000000000..98ce62560a --- /dev/null +++ b/tools/toolchain/scripts/stage4/install_libGint.sh @@ -0,0 +1,77 @@ +#!/bin/bash -e + +# TODO: Review and if possible fix shellcheck errors. +# shellcheck disable=all + +[ "${BASH_SOURCE[0]}" ] && SCRIPT_NAME="${BASH_SOURCE[0]}" || SCRIPT_NAME=$0 +SCRIPT_DIR="$(cd "$(dirname "$SCRIPT_NAME")/.." && pwd -P)" +libgint_ver="v1" +libgint_sha256=cc0dfeb6022ebfe0c3028a131045ff49b4c1005ad9cb77c5736fb3dac045b192 + +source "${SCRIPT_DIR}"/common_vars.sh +source "${SCRIPT_DIR}"/tool_kit.sh +source "${SCRIPT_DIR}"/signal_trap.sh +source "${INSTALLDIR}"/toolchain.conf +source "${INSTALLDIR}"/toolchain.env + +[ -f "${BUILDDIR}/setup_libGint" ] && rm "${BUILDDIR}/setup_libGint" + +! [ -d "${BUILDDIR}" ] && mkdir -p "${BUILDDIR}" +cd "${BUILDDIR}" + +case "$with_libgint" in + __INSTALL__) + echo "==================== Installing libGint ====================" + pkg_install_dir="${INSTALLDIR}/libGint-${libgint_ver}" + install_lock_file="$pkg_install_dir/install_successful" + if verify_checksums "${install_lock_file}"; then + echo "libGint-${libgint_ver} is already installed, skipping it." + else + retrieve_package "${libgint_sha256}" "libGint-${libgint_ver}.tar.gz" + # The tar called libGint-v1 expands into libGint-release_v1 + [ -d libGint-release_${libgint_ver} ] && rm -rf libGint-release_${libgint_ver} + tar -xzf libGint-${libgint_ver}.tar.gz + echo "Installing from scratch into ${pkg_install_dir}" + cd libGint-release_${libgint_ver} + + make -j $(get_nprocs) > make.log 2>&1 || tail -n ${LOG_LINES} make.log + make install PREFIX="${pkg_install_dir}" > install.log 2>&1 || tail -n ${LOG_LINES} install.log + write_checksums "${install_lock_file}" "${SCRIPT_DIR}/stage4/$(basename ${SCRIPT_NAME})" + fi + ;; + __SYSTEM__) + echo "==================== Finding libGint from system paths ====================" + check_lib -lGint "libGint" + pkg_install_dir="$(dirname $(dirname $(find_in_paths "libxs.*" $LIB_PATHS)))" + ;; + __DONTUSE__) ;; + *) + echo "==================== Linking libGint to user paths ====================" + echo "$with_libgint" + pkg_install_dir="$with_libgint" + check_dir "$pkg_install_dir/lib" + check_dir "$pkg_install_dir/include" + ;; +esac + +if [ "$with_libgint" != "__DONTUSE__" ]; then + cat << EOF > "${BUILDDIR}/setup_libGint" +export LIBGINT_VER="${libgint_ver}" +EOF + if [ "$with_libgint" != "__SYSTEM__" ]; then + cat << EOF >> "${BUILDDIR}/setup_libGint" +prepend_path LD_LIBRARY_PATH "$pkg_install_dir/lib" +prepend_path LD_RUN_PATH "$pkg_install_dir/lib" +prepend_path LIBRARY_PATH "$pkg_install_dir/lib" +prepend_path PKG_CONFIG_PATH "$pkg_install_dir/lib/pkgconfig" +prepend_path CMAKE_PREFIX_PATH "$pkg_install_dir" +EOF + fi + filter_setup "${BUILDDIR}/setup_libGint" "${SETUPFILE}" +fi +cd "${ROOTDIR}" + +load "${BUILDDIR}/setup_libGint" +write_toolchain_env "${INSTALLDIR}" + +report_timing "libGint" diff --git a/tools/toolchain/scripts/stage4/install_stage4.sh b/tools/toolchain/scripts/stage4/install_stage4.sh index 772b13a4b8..37bdcc7410 100755 --- a/tools/toolchain/scripts/stage4/install_stage4.sh +++ b/tools/toolchain/scripts/stage4/install_stage4.sh @@ -6,6 +6,7 @@ ./scripts/stage4/install_libxsmm.sh ./scripts/stage4/install_libxs.sh ./scripts/stage4/install_libxstream.sh +./scripts/stage4/install_libGint.sh ./scripts/stage4/install_scalapack.sh ./scripts/stage4/install_cusolvermp.sh ./scripts/stage4/install_cosma.sh