comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 22:25:32 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

PR for adding libGint to CP2K (#5446)

Browse source
This commit is contained in:
Marcello Puligheddu 2026-07-02 13:30:11 +01:00 • committed by GitHub
parent 3a6a80ebfc
commit de7340f8be
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
29 changed files with 1608 additions and 232 deletions
Download patch file Download diff file Expand all files Collapse all files

18
CMakeLists.txt
View file

@ -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()

39
cmake/modules/FindLibGint.cmake Normal file
View file

@ -0,0 +1,39 @@
#!-------------------------------------------------------------------------------------------------!
#! CP2K: A general program to perform molecular dynamics simulations !
#! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
#! !
#! 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)

9
docs/technologies/libraries.md
View file

@ -346,3 +346,12 @@ MiMiC - Multiscale simulation framework
<https://https://mimic-project.org>
- For more information about the framework and supported programs see <https://mimic-project.org>
- 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.

3
src/CMakeLists.txt
View file

@ -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(
$<$<BOOL:${CP2K_USE_MIMIC}>:CP2K::MIMIC::mcl>
$<$<BOOL:${CP2K_USE_LIBINT2}>:Libint2::int2>
$<$<BOOL:${CP2K_USE_LIBFCI}>:libfci::fci>
$<$<BOOL:${CP2K_USE_LIBGINT}>:cp2k::LibGint::libGint>
$<$<BOOL:${CP2K_USE_COSMA}>:cp2k::cosma>
$<$<BOOL:${CP2K_USE_DLAF}>:${CP2K_DLAF_FORTRAN_TARGET}>
$<$<BOOL:${CP2K_USE_GREENX}>:cp2k::greenx>
@ -1861,6 +1863,7 @@ target_compile_definitions(
$<$<BOOL:${CP2K_USE_GREENX}>:__GREENX>
$<$<OR:$<BOOL:${CP2K_USE_FFTW3_}>,$<BOOL:${CP2K_USE_FFTW3_MKL_}>>:__FFTW3>
$<$<BOOL:${CP2K_USE_LIBINT2}>:__LIBINT>
$<$<BOOL:${CP2K_USE_LIBGINT}>:__LIBGINT>
$<$<BOOL:${CP2K_USE_LIBFCI}>:__LIBFCI>
$<$<BOOL:${CP2K_USE_LIBTORCH}>:__LIBTORCH>
$<$<BOOL:${CP2K_USE_COSMA}>:__COSMA>

6
src/common/t_c_g0.F
View file

@ -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

3
src/cp2k_info.F
View file

@ -335,6 +335,9 @@ CONTAINS
flags = TRIM(flags)//" openpmd"
#endif
#if defined(__LIBGINT)
flags = TRIM(flags)//" libGint"
#endif
END FUNCTION cp2k_flags
! **************************************************************************************************

611
src/hfx_energy_potential.F
View file

@ -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

131
src/hfx_pair_list_methods.F
View file

@ -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 ...

7
src/hfx_types.F
View file

@ -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, &

5
src/input_constants.F
View file

@ -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, &

21
src/input_cp2k_hfx.F
View file

@ -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.", &

696
src/libgint_wrapper.F Normal file
View file

@ -0,0 +1,696 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
! !
! 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 <AB(g) and CD(n)> 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

2
src/offload/offload_library.c
View file

@ -23,7 +23,7 @@
#if defined(__OFFLOAD_PROFILING)
#if defined(__OFFLOAD_CUDA)
#include <nvToolsExt.h>
#include <nvtx3/nvToolsExt.h>
#elif defined(__OFFLOAD_HIP) && defined(__HIP_PLATFORM_AMD__)
#include <roctracer/roctx.h>
#endif

67
tests/QS/regtest-libgint/H1.inp Normal file
View file

@ -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

89
tests/QS/regtest-libgint/H2O-HFX-6.inp Normal file
View file

@ -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

4
tests/QS/regtest-libgint/TEST_FILES.toml Normal file
View file

@ -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

10
tests/QS/regtest-libgint/basis.dat Normal file
View file

@ -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
#

18
tests/QS/regtest-libgint/tidy Normal file
View file

@ -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

1
tests/TEST_DIRS
View file

@ -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

1
tools/docker/Dockerfile.test_cuda_A100
View file

@ -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

1
tools/docker/Dockerfile.test_cuda_P100
View file

@ -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

1
tools/docker/Dockerfile.test_cuda_V100
View file

@ -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

1
tools/docker/Dockerfile.test_performance_cuda_A100
View file

@ -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

1
tools/docker/Dockerfile.test_performance_cuda_P100
View file

@ -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

1
tools/docker/Dockerfile.test_performance_cuda_V100
View file

@ -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

1
tools/docker/generate_dockerfiles.py
View file

@ -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,

15
tools/toolchain/install_cp2k_toolchain.sh
View file

@ -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__"

77
tools/toolchain/scripts/stage4/install_libGint.sh Executable file
View file

@ -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"

1
tools/toolchain/scripts/stage4/install_stage4.sh
View file

@ -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