comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-21 14:35:15 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Gauxc object caching (#5340)

Browse source 
Co-authored-by: Thomas D. Kuehne <tkuehne@cp2k.org>
Co-authored-by: Thomas Kuehne <kuehne88@users.noreply.github.com>
This commit is contained in:
Franz Pöschel 2026-07-15 13:48:13 +02:00 • committed by GitHub
parent ea663e99cf
commit 39fe4c47c6
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
30 changed files with 2146 additions and 328 deletions
Download patch file Download diff file Expand all files Collapse all files

1
src/CMakeLists.txt
View file

@ -1378,6 +1378,7 @@ list(
xc/xc.F
xc/xc_functionals_utilities.F
xc/xc_fxc_kernel.F
xc/xc_gauxc_cache.F
xc/xc_gauxc_functional.F
xc/xc_gauxc_interface.F
xc/xc_hcth.F

7
src/qs_environment_types.F
View file

@ -163,6 +163,8 @@ MODULE qs_environment_types
USE wannier_states_types, ONLY: wannier_centres_type
USE xas_env_types, ONLY: xas_env_release,&
xas_environment_type
USE xc_gauxc_cache, ONLY: cp_gauxc_cache_type,&
gauxc_cache_release
#include "./base/base_uses.f90"
IMPLICIT NONE
@ -322,6 +324,7 @@ MODULE qs_environment_types
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_last_converged => NULL()
! tblite
TYPE(tblite_type), POINTER :: tb_tblite => Null()
TYPE(cp_gauxc_cache_type) :: gauxc_cache
END TYPE qs_environment_type
CONTAINS
@ -1680,6 +1683,8 @@ CONTAINS
CALL deallocate_tblite_type(qs_env%tb_tblite)
END IF
CALL gauxc_cache_release(qs_env%gauxc_cache)
END SUBROUTINE qs_env_release
! **************************************************************************************************
@ -1904,6 +1909,8 @@ CONTAINS
CALL deallocate_tblite_type(qs_env%tb_tblite)
END IF
CALL gauxc_cache_release(qs_env%gauxc_cache)
END SUBROUTINE qs_env_part_release
END MODULE qs_environment_types

27
src/qs_ks_methods.F
View file

@ -143,6 +143,7 @@ MODULE qs_ks_methods
get_gauxc_section,&
xc_section_uses_native_skala_grid
USE smeagol_interface, ONLY: smeagol_shift_v_hartree
USE string_utilities, ONLY: uppercase
USE surface_dipole, ONLY: calc_dipsurf_potential
USE tblite_ks_matrix, ONLY: build_tblite_ks_matrix
USE virial_types, ONLY: virial_type
@ -199,14 +200,15 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_ks_build_kohn_sham_matrix'
CHARACTER(len=default_string_length) :: name
CHARACTER(len=default_string_length) :: gauxc_model_name, name
INTEGER :: ace_rebuild_frequency, atom_a, handle, &
iatom, ikind, img, ispin, natom, &
nimages, nspins, output_unit
INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
LOGICAL :: ace_active, do_adiabatic_rescaling, do_ddapc, do_hfx, do_kpoints, do_ppl, dokp, &
gapw, gapw_xc, just_energy_xc, lrigpw, my_print, native_grid_diagnostics, &
native_grid_use_cuda, native_skala_restore_exc, rigpw, use_gauxc_matrix, use_virial
gapw, gapw_xc, gauxc_model_none, just_energy_xc, lrigpw, my_print, &
native_grid_diagnostics, native_grid_use_cuda, native_skala_restore_exc, rigpw, &
use_gauxc_matrix, use_virial
LOGICAL, SAVE :: native_grid_cpu_kpoints_warned = .FALSE.
REAL(KIND=dp) :: ecore_ppl, edisp, ee_ener, ekin_mol, &
mulliken_order_p, &
@ -922,12 +924,19 @@ CONTAINS
NULLIFY (rho1)
IF (dft_control%use_gauxc .AND. (gapw .OR. gapw_xc) .AND. &
.NOT. xc_section_uses_native_skala_grid(xc_section)) THEN
! Molecular GauXC evaluates the XC term outside xc_derivatives.
! The accurate-XCINT force correction would otherwise try to
! evaluate the GAUXC section through CP2K's local functional path.
! Native-grid SKALA can provide this correction through the
! CP2K grid path and needs it for GAPW_XC force consistency.
CONTINUE
gauxc_model_none = .FALSE.
gauxc_section => get_gauxc_section(xc_section)
IF (ASSOCIATED(gauxc_section)) THEN
CALL section_vals_val_get(gauxc_section, "MODEL", c_val=gauxc_model_name)
gauxc_model_name = ADJUSTL(gauxc_model_name)
CALL uppercase(gauxc_model_name)
gauxc_model_none = (TRIM(gauxc_model_name) == "" .OR. &
TRIM(gauxc_model_name) == "NONE")
END IF
IF (gauxc_model_none .AND. &
(gapw_xc .OR. gauxc_gapw_has_paw_pseudopotentials(qs_kind_set))) THEN
CALL accint_weight_force(qs_env, rho_struct, rho1, 0, xc_section)
END IF
ELSE
CALL accint_weight_force(qs_env, rho_struct, rho1, 0, xc_section)
END IF

20
src/qs_vxc_atom.F
View file

@ -15,6 +15,8 @@ MODULE qs_vxc_atom
USE basis_set_types, ONLY: get_gto_basis_set,&
gto_basis_set_type
USE cp_control_types, ONLY: dft_control_type
USE external_potential_types, ONLY: gth_potential_type,&
sgp_potential_type
USE input_constants, ONLY: xc_none
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type,&
@ -120,14 +122,16 @@ CONTAINS
INTEGER :: bo(2), gapw_density_partition, handle, &
iat, iatom, idir, ikind, ir, jdir, &
myfun, na, natom, nr, nspins, num_pe
myfun, na, natom, nr, nspins, num_pe, &
zatom
INTEGER, DIMENSION(2, 3) :: bounds
INTEGER, DIMENSION(:), POINTER :: atom_list
LOGICAL :: accint, donlcc, evaluate_hard, evaluate_soft, gradient_f, lsd, &
my_calculate_forces, nlcc, paw_atom, skala_atom_grid, tau_f, use_virial
REAL(dp) :: agr, alpha, density_cut, exc_h, exc_s, &
gradient_cut, &
my_adiabatic_rescale_factor, tau_cut
my_adiabatic_rescale_factor, tau_cut, &
zeff
REAL(dp), DIMENSION(1, 1, 1) :: tau_d
REAL(dp), DIMENSION(1, 1, 1, 1) :: rho_d
REAL(dp), DIMENSION(3) :: skala_atom_force_h, skala_atom_force_s
@ -140,6 +144,7 @@ CONTAINS
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(dft_control_type), POINTER :: dft_control
TYPE(grid_atom_type), POINTER :: grid_atom
TYPE(gth_potential_type), POINTER :: gth_potential
TYPE(gto_basis_set_type), POINTER :: basis_1c
TYPE(harmonics_atom_type), POINTER :: harmonics
TYPE(mp_para_env_type), POINTER :: para_env
@ -151,6 +156,7 @@ CONTAINS
TYPE(rho_atom_type), DIMENSION(:), POINTER :: my_rho_atom_set
TYPE(rho_atom_type), POINTER :: rho_atom
TYPE(section_vals_type), POINTER :: input, my_xc_section, xc_fun_section
TYPE(sgp_potential_type), POINTER :: sgp_potential
TYPE(tau_basis_cache_type) :: tau_basis_cache
TYPE(virial_type), POINTER :: virial
TYPE(xc_derivative_set_type) :: deriv_set
@ -165,6 +171,7 @@ CONTAINS
NULLIFY (my_kind_set)
NULLIFY (atomic_kind_set)
NULLIFY (grid_atom)
NULLIFY (gth_potential)
NULLIFY (force)
NULLIFY (harmonics)
NULLIFY (input)
@ -173,6 +180,7 @@ CONTAINS
NULLIFY (rho_atom)
NULLIFY (my_rho_atom_set)
NULLIFY (rho_nlcc)
NULLIFY (sgp_potential)
NULLIFY (virial)
my_calculate_forces = .FALSE.
IF (PRESENT(calculate_forces)) my_calculate_forces = calculate_forces
@ -250,11 +258,17 @@ CONTAINS
DO ikind = 1, SIZE(atomic_kind_set)
CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
NULLIFY (gth_potential, sgp_potential)
CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
harmonics=harmonics, grid_atom=grid_atom)
gth_potential=gth_potential, harmonics=harmonics, &
grid_atom=grid_atom, sgp_potential=sgp_potential, &
zatom=zatom, zeff=zeff)
CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
IF (.NOT. paw_atom) CYCLE
IF (skala_atom_grid .AND. &
(ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) .AND. &
ABS(zeff - REAL(zatom, dp)) <= 1.0E-10_dp) CYCLE
nr = grid_atom%nr
na = grid_atom%ng_sphere

10
src/skala_gpw_functional.F
View file

@ -109,7 +109,7 @@ CONTAINS
END FUNCTION xc_section_uses_native_skala_grid
! **************************************************************************************************
!> \brief Return true if the GAUXC subsection requests a Skala-style model.
!> \brief Return true if the GAUXC subsection requests a model evaluation.
!> \param xc_section ...
!> \return ...
! **************************************************************************************************
@ -117,16 +117,20 @@ CONTAINS
TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
LOGICAL :: uses_gauxc_model
CHARACTER(len=default_path_length) :: model_key, model_name
CHARACTER(len=default_path_length) :: model_key, model_name, xc_key, xc_name
TYPE(section_vals_type), POINTER :: gauxc_section
uses_gauxc_model = .FALSE.
gauxc_section => get_gauxc_section(xc_section)
IF (ASSOCIATED(gauxc_section)) THEN
CALL section_vals_val_get(gauxc_section, "MODEL", c_val=model_name)
CALL section_vals_val_get(gauxc_section, "FUNCTIONAL", c_val=xc_name)
model_key = ADJUSTL(model_name)
xc_key = ADJUSTL(xc_name)
CALL uppercase(model_key)
uses_gauxc_model = (TRIM(model_key) /= "" .AND. TRIM(model_key) /= "NONE")
CALL uppercase(xc_key)
uses_gauxc_model = (TRIM(model_key) /= "" .AND. TRIM(model_key) /= "NONE" .AND. &
TRIM(model_key) /= TRIM(xc_key))
END IF
END FUNCTION xc_section_uses_gauxc_model

34
src/torch_c_api.cpp
View file

@ -7,6 +7,7 @@
#if defined(__LIBTORCH)
#include <ATen/Parallel.h>
#include <c10/core/DeviceGuard.h>
#include <torch/csrc/api/include/torch/cuda.h>
#include <torch/script.h>
@ -16,6 +17,7 @@
#include <cassert>
#include <cfenv>
#include <climits>
#include <cstdlib>
#include <cstring>
#include <string>
@ -47,7 +49,38 @@ private:
******************************************************************************/
static bool use_cuda_if_available = true;
static bool get_positive_int_env(const char *name, int &value) {
const char *raw = std::getenv(name);
if (raw == nullptr || raw[0] == '\0') {
return false;
}
char *end = nullptr;
const long parsed = std::strtol(raw, &end, 10);
if (end == raw || *end != '\0' || parsed <= 0 || parsed > INT_MAX) {
return false;
}
value = static_cast<int>(parsed);
return true;
}
static void initialize_torch_threads_from_env() {
static bool initialized = false;
if (initialized) {
return;
}
initialized = true;
int num_threads = 0;
if (get_positive_int_env("CP2K_TORCH_NUM_THREADS", num_threads)) {
at::set_num_threads(num_threads);
}
if (get_positive_int_env("CP2K_TORCH_NUM_INTEROP_THREADS", num_threads)) {
at::set_num_interop_threads(num_threads);
}
}
static torch::Device get_device() {
initialize_torch_threads_from_env();
if (!use_cuda_if_available || !torch::cuda::is_available()) {
return torch::kCPU;
}
@ -107,6 +140,7 @@ static torch_c_tensor_t *tensor_from_array(const torch::Dtype dtype,
const bool req_grad, const int ndims,
const int64_t sizes[],
void *source) {
initialize_torch_threads_from_env();
const auto opts = torch::TensorOptions().dtype(dtype).requires_grad(req_grad);
const auto sizes_ref = c10::IntArrayRef(sizes, ndims);
return new torch_c_tensor_t(torch::from_blob(source, sizes_ref, opts));

169
src/xc/xc_derivatives.F
View file

@ -11,7 +11,9 @@ MODULE xc_derivatives
USE input_section_types, ONLY: section_vals_get_subs_vals2,&
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: dp
USE kinds, ONLY: default_string_length,&
dp
USE string_utilities, ONLY: uppercase
USE xc_b97, ONLY: b97_lda_eval,&
b97_lda_info,&
b97_lsd_eval,&
@ -281,9 +283,14 @@ CONTAINS
CALL pbe_lda_info(functional, reference, shortform, needs, max_deriv)
END IF
CASE ("GAUXC")
CALL skala_info(functional, lsd, reference, shortform, needs, max_deriv)
! Note: SKALA functional routes through apply_gauxc in qs_ks_methods.F
! when USE_GAUXC = .TRUE. (requires dft_control%use_gauxc to be set)
IF (gauxc_model_none_selected(functional)) THEN
CALL gauxc_model_none_xc_info(functional, lsd, reference, shortform, &
needs, max_deriv, print_warn)
ELSE
CALL skala_info(functional, lsd, reference, shortform, needs, max_deriv)
! Note: SKALA functional routes through apply_gauxc in qs_ks_methods.F
! when USE_GAUXC = .TRUE. (requires dft_control%use_gauxc to be set)
END IF
CASE ("XWPBE")
IF (lsd) THEN
CALL xwpbe_lsd_info(reference, shortform, needs, max_deriv)
@ -507,7 +514,11 @@ CONTAINS
CALL pbe_lda_eval(rho_set, deriv_set, deriv_order, functional)
END IF
CASE ("GAUXC")
CPABORT(abort_message_skala)
IF (gauxc_model_none_selected(functional)) THEN
CALL gauxc_model_none_xc_eval(functional, lsd, rho_set, deriv_set, deriv_order)
ELSE
CPABORT(abort_message_skala)
END IF
CASE ("XWPBE")
IF (lsd) THEN
CALL xwpbe_lsd_eval(rho_set, deriv_set, deriv_order, functional)
@ -552,6 +563,154 @@ CONTAINS
CALL timestop(handle)
END SUBROUTINE xc_functional_eval
! **************************************************************************************************
!> \brief true for GAUXC sections that wrap a conventional LibXC functional
!> \param functional the GAUXC section
!> \return whether MODEL NONE is active
! **************************************************************************************************
FUNCTION gauxc_model_none_selected(functional)
TYPE(section_vals_type), POINTER :: functional
LOGICAL :: gauxc_model_none_selected
CHARACTER(LEN=default_string_length) :: model_key, model_name, xc_fun_name, &
xc_key
CALL section_vals_val_get(functional, "MODEL", c_val=model_name)
CALL section_vals_val_get(functional, "FUNCTIONAL", c_val=xc_fun_name)
model_key = ADJUSTL(model_name)
xc_key = ADJUSTL(xc_fun_name)
CALL uppercase(model_key)
CALL uppercase(xc_key)
gauxc_model_none_selected = (TRIM(model_key) == "" .OR. TRIM(model_key) == "NONE" .OR. &
TRIM(model_key) == TRIM(xc_key))
END FUNCTION gauxc_model_none_selected
! **************************************************************************************************
!> \brief map GAUXC MODEL NONE shorthand names to LibXC exchange/correlation components
!> \param functional the GAUXC section
!> \param xc_fun_name the GAUXC FUNCTIONAL value
!> \param libxc_names LibXC section names to evaluate and add
!> \param nfunc number of LibXC components
! **************************************************************************************************
SUBROUTINE gauxc_model_none_libxc_names(functional, xc_fun_name, libxc_names, nfunc)
TYPE(section_vals_type), POINTER :: functional
CHARACTER(LEN=*), INTENT(OUT) :: xc_fun_name
CHARACTER(LEN=*), DIMENSION(:), INTENT(OUT) :: libxc_names
INTEGER, INTENT(OUT) :: nfunc
CHARACTER(LEN=default_string_length) :: xc_key
CALL section_vals_val_get(functional, "FUNCTIONAL", c_val=xc_fun_name)
xc_key = ADJUSTL(xc_fun_name)
CALL uppercase(xc_key)
libxc_names(:) = ""
SELECT CASE (TRIM(xc_key))
CASE ("LDA", "PADE")
nfunc = 2
libxc_names(1) = "LDA_X"
libxc_names(2) = "LDA_C_PW"
CASE ("VWN")
nfunc = 2
libxc_names(1) = "LDA_X"
libxc_names(2) = "LDA_C_VWN"
CASE ("PBE")
nfunc = 2
libxc_names(1) = "GGA_X_PBE"
libxc_names(2) = "GGA_C_PBE"
CASE ("BLYP")
nfunc = 2
libxc_names(1) = "GGA_X_B88"
libxc_names(2) = "GGA_C_LYP"
CASE ("BP")
nfunc = 2
libxc_names(1) = "GGA_X_B88"
libxc_names(2) = "GGA_C_P86"
CASE ("TPSS")
nfunc = 2
libxc_names(1) = "MGGA_X_TPSS"
libxc_names(2) = "MGGA_C_TPSS"
CASE ("R2SCAN")
nfunc = 2
libxc_names(1) = "MGGA_X_R2SCAN"
libxc_names(2) = "MGGA_C_R2SCAN"
CASE DEFAULT
nfunc = 1
libxc_names(1) = TRIM(xc_key)
END SELECT
END SUBROUTINE gauxc_model_none_libxc_names
! **************************************************************************************************
!> \brief needs information for GAUXC MODEL NONE one-center GAPW corrections
!> \param functional the GAUXC section
!> \param lsd whether spin-polarized derivatives are needed
!> \param reference reference string for the wrapped functional
!> \param shortform short name for printout
!> \param needs density ingredients needed by the wrapped functional
!> \param max_deriv maximum implemented derivative order
!> \param print_warn whether LibXC should print development warnings
! **************************************************************************************************
SUBROUTINE gauxc_model_none_xc_info(functional, lsd, reference, shortform, &
needs, max_deriv, print_warn)
TYPE(section_vals_type), POINTER :: functional
LOGICAL, INTENT(in) :: lsd
CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
INTEGER, INTENT(out), OPTIONAL :: max_deriv
LOGICAL, INTENT(IN), OPTIONAL :: print_warn
CHARACTER(LEN=default_string_length) :: libxc_names(2), xc_fun_name
INTEGER :: ifunc, max_deriv_i, max_deriv_min, nfunc
CALL gauxc_model_none_libxc_names(functional, xc_fun_name, libxc_names, nfunc)
max_deriv_min = HUGE(max_deriv_min)
DO ifunc = 1, nfunc
IF (lsd) THEN
CALL libxc_lsd_info(functional, needs=needs, max_deriv=max_deriv_i, &
print_warn=print_warn, &
func_name_override=TRIM(libxc_names(ifunc)))
ELSE
CALL libxc_lda_info(functional, needs=needs, max_deriv=max_deriv_i, &
print_warn=print_warn, &
func_name_override=TRIM(libxc_names(ifunc)))
END IF
max_deriv_min = MIN(max_deriv_min, max_deriv_i)
END DO
IF (PRESENT(max_deriv)) max_deriv = max_deriv_min
IF (PRESENT(reference)) &
reference = "Functional computed by GauXC (underlying: "//TRIM(xc_fun_name)//")"
IF (PRESENT(shortform)) shortform = "GAUXC ("//TRIM(xc_fun_name)//")"
END SUBROUTINE gauxc_model_none_xc_info
! **************************************************************************************************
!> \brief one-center GAPW correction evaluation for GAUXC MODEL NONE
!> \param functional the GAUXC section
!> \param lsd whether spin-polarized derivatives are evaluated
!> \param rho_set density ingredients
!> \param deriv_set derivative accumulator
!> \param deriv_order derivative order to evaluate
! **************************************************************************************************
SUBROUTINE gauxc_model_none_xc_eval(functional, lsd, rho_set, deriv_set, deriv_order)
TYPE(section_vals_type), POINTER :: functional
LOGICAL, INTENT(in) :: lsd
TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
INTEGER, INTENT(IN) :: deriv_order
CHARACTER(LEN=default_string_length) :: libxc_names(2), xc_fun_name
INTEGER :: ifunc, nfunc
CALL gauxc_model_none_libxc_names(functional, xc_fun_name, libxc_names, nfunc)
DO ifunc = 1, nfunc
IF (lsd) THEN
CALL libxc_lsd_eval(rho_set, deriv_set, deriv_order, functional, &
func_name_override=TRIM(libxc_names(ifunc)))
ELSE
CALL libxc_lda_eval(rho_set, deriv_set, deriv_order, functional, &
func_name_override=TRIM(libxc_names(ifunc)))
END IF
END DO
END SUBROUTINE gauxc_model_none_xc_eval
! **************************************************************************************************
!> \brief ...
!> \param functionals a section containing the functional combination to be

430
src/xc/xc_gauxc_cache.F Normal file
View file

@ -0,0 +1,430 @@
!--------------------------------------------------------------------------------------------------!
! 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 !
!--------------------------------------------------------------------------------------------------!
#ifdef __GAUXC
#include "gauxc/gauxc_config.f"
#endif
MODULE xc_gauxc_cache
USE iso_c_binding, ONLY: c_double
USE kinds, ONLY: default_path_length,&
default_string_length
USE message_passing, ONLY: mp_comm_self,&
mp_para_env_type
USE particle_types, ONLY: particle_type
USE qs_kind_types, ONLY: qs_kind_type
USE xc_gauxc_interface, ONLY: &
cp_gauxc_basisset_type, cp_gauxc_grid_type, cp_gauxc_integrator_type, &
cp_gauxc_molecule_type, cp_gauxc_status_type, gauxc_check_status, gauxc_create_basisset, &
gauxc_create_grid, gauxc_create_integrator, gauxc_create_molecule, gauxc_destroy_basisset, &
gauxc_destroy_grid, gauxc_destroy_integrator, gauxc_destroy_molecule
#include "../base/base_uses.f90"
IMPLICIT NONE
PRIVATE
TYPE cp_gauxc_cache_params
INTEGER :: natom = -1
INTEGER :: nspins = -1
INTEGER :: batch_size = -1
REAL(c_double) :: device_runtime_fill_fraction = 0.0_c_double
CHARACTER(LEN=default_string_length) :: xc_fun_name = ""
CHARACTER(LEN=default_string_length) :: grid_type = ""
CHARACTER(LEN=default_string_length) :: radial_quadrature = ""
CHARACTER(LEN=default_string_length) :: pruning_scheme = ""
CHARACTER(LEN=default_string_length) :: lb_exec_space = ""
CHARACTER(LEN=default_string_length) :: int_exec_space = ""
CHARACTER(LEN=default_path_length) :: model_eval_name = ""
CHARACTER(len=default_string_length) :: lwd_kernel = ""
LOGICAL :: use_mpi_runtime = .FALSE.
LOGICAL :: use_gradient_mpi_runtime = .FALSE.
LOGICAL :: use_self_runtime = .FALSE.
LOGICAL :: use_gradient_self_runtime = .FALSE.
LOGICAL :: use_fd_gradient = .FALSE.
LOGICAL :: use_gauxc_model = .FALSE.
END TYPE cp_gauxc_cache_params
TYPE, extends(cp_gauxc_cache_params) :: cp_gauxc_cache_type
TYPE(cp_gauxc_molecule_type) :: molecule
TYPE(cp_gauxc_basisset_type) :: basisset
TYPE(cp_gauxc_grid_type) :: grid
TYPE(cp_gauxc_grid_type) :: gradient_grid
TYPE(cp_gauxc_integrator_type) :: integrator
TYPE(cp_gauxc_integrator_type) :: gradient_integrator
INTEGER :: mpi_comm
LOGICAL :: is_init = .FALSE.
REAL(c_double), DIMENSION(:, :), ALLOCATABLE :: last_positions
CONTAINS
PROCEDURE, PUBLIC :: needs_gradient_grid => gauxc_cache_type_needs_gradient_grid
PROCEDURE, PUBLIC :: max_l => gauxc_cache_type_max_l
END TYPE cp_gauxc_cache_type
PUBLIC :: &
cp_gauxc_cache_type, &
gauxc_cache_init, &
gauxc_cache_is_valid, &
gauxc_cache_release, cp_gauxc_cache_params
CONTAINS
! **************************************************************************************************
!> \brief ...
!> \param this ...
!> \return ...
! **************************************************************************************************
FUNCTION gauxc_cache_type_needs_gradient_grid(this) RESULT(res)
class(cp_gauxc_cache_type) :: this
LOGICAL :: res
res = this%use_gradient_self_runtime
END FUNCTION gauxc_cache_type_needs_gradient_grid
! **************************************************************************************************
!> \brief ...
!> \param this ...
!> \return ...
! **************************************************************************************************
FUNCTION gauxc_cache_type_max_l(this) RESULT(res)
class(cp_gauxc_cache_type) :: this
INTEGER :: res
res = this%basisset%max_l
END FUNCTION gauxc_cache_type_max_l
! **************************************************************************************************
!> \brief ...
!> \param cache ...
!> \param params ...
!> \param para_env ...
!> \param particle_set ...
!> \param qs_kind_set ...
!> \param status ...
! **************************************************************************************************
SUBROUTINE gauxc_cache_init( &
cache, &
params, &
para_env, &
particle_set, &
qs_kind_set, &
status)
TYPE(cp_gauxc_cache_type), INTENT(inout) :: cache
TYPE(cp_gauxc_cache_params), INTENT(inout) :: params
TYPE(mp_para_env_type), INTENT(in), POINTER :: para_env
TYPE(particle_type), DIMENSION(:), INTENT(in), &
POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), INTENT(in), &
POINTER :: qs_kind_set
TYPE(cp_gauxc_status_type), INTENT(inout) :: status
#ifdef __GAUXC
IF (gauxc_cache_is_valid( &
cache, &
params, &
para_env, &
particle_set)) THEN
status%status%code = 0
RETURN
END IF
CALL gauxc_cache_release(cache)
cache%molecule = gauxc_create_molecule( &
particle_set, &
status)
IF (status%status%code /= 0) GO TO 100
cache%basisset = gauxc_create_basisset( &
qs_kind_set, &
particle_set, &
status)
IF (status%status%code /= 0) GO TO 200
params%use_fd_gradient = &
params%use_fd_gradient .AND. (cache%max_l() > 3)
params%use_gradient_mpi_runtime = &
params%use_gradient_mpi_runtime .AND. .NOT. params%use_fd_gradient
params%use_gradient_self_runtime = &
params%use_gradient_self_runtime .AND. .NOT. params%use_gradient_mpi_runtime
IF (params%use_fd_gradient .AND. para_env%mepos == 0) THEN
CALL cp_warn( &
__LOCATION__, &
"Using finite-difference GauXC XC gradients for METHOD GAPW/GAPW_XC with "// &
"basis functions beyond f shells. The upstream analytical GauXC gradient path is "// &
"not yet reliable for this case.")
END IF
IF (params%use_self_runtime .OR. .NOT. params%use_gauxc_model) THEN
! SKALA currently needs a replicated molecular runtime for reproducible
! open-shell densities across CP2K MPI ranks.
! Conventional GauXC also uses a replicated runtime here: CP2K has
! already allreduced the dense AO density matrix on every rank.
cache%grid = gauxc_create_grid( &
cache%molecule, &
cache%basisset, &
params%grid_type, &
params%radial_quadrature, &
params%pruning_scheme, &
params%lb_exec_space, &
params%batch_size, &
params%device_runtime_fill_fraction, &
status, &
mpi_comm=mp_comm_self%get_handle(), &
force_new_runtime=.TRUE.)
ELSE
! Use the QS force-evaluation communicator rather than GauXC's global
! runtime. In mixed CDFT the diabatic states can be built on disjoint
! MPI subgroups; using the global communicator would make GauXC wait
! for ranks that are working on another state.
cache%grid = gauxc_create_grid( &
cache%molecule, &
cache%basisset, &
params%grid_type, &
params%radial_quadrature, &
params%pruning_scheme, &
params%lb_exec_space, &
params%batch_size, &
params%device_runtime_fill_fraction, &
status, &
mpi_comm=para_env%get_handle())
END IF
IF (status%status%code /= 0) GO TO 300
cache%integrator = gauxc_create_integrator( &
TRIM(params%xc_fun_name), &
cache%grid, &
params%int_exec_space, &
params%lwd_kernel, &
params%nspins, &
status)
IF (status%status%code /= 0) GO TO 400
IF (params%use_gradient_self_runtime) THEN
! Upstream GauXC does not yet support OneDFT/SKALA nuclear gradients
! on an MPI runtime. Keep the energy/VXC path on the normal MPI
! runtime and use an isolated runtime only for the replicated gradient.
cache%gradient_grid = gauxc_create_grid( &
cache%molecule, &
cache%basisset, &
params%grid_type, &
params%radial_quadrature, &
params%pruning_scheme, &
params%lb_exec_space, &
params%batch_size, &
params%device_runtime_fill_fraction, &
status, &
mpi_comm=mp_comm_self%get_handle(), &
force_new_runtime=.TRUE.)
IF (status%status%code /= 0) GO TO 500
cache%gradient_integrator = gauxc_create_integrator( &
TRIM(params%xc_fun_name), &
cache%gradient_grid, &
params%int_exec_space, &
params%lwd_kernel, &
params%nspins, &
status)
IF (status%status%code /= 0) GO TO 600
END IF
cache%natom = params%natom
cache%nspins = params%nspins
cache%batch_size = params%batch_size
cache%device_runtime_fill_fraction = REAL(params%device_runtime_fill_fraction, c_double)
cache%xc_fun_name = TRIM(params%xc_fun_name)
cache%grid_type = TRIM(params%grid_type)
cache%radial_quadrature = TRIM(params%radial_quadrature)
cache%pruning_scheme = TRIM(params%pruning_scheme)
cache%lb_exec_space = TRIM(params%lb_exec_space)
cache%int_exec_space = TRIM(params%int_exec_space)
cache%model_eval_name = TRIM(params%model_eval_name)
cache%lwd_kernel = TRIM(params%lwd_kernel)
cache%use_gradient_mpi_runtime = params%use_gradient_mpi_runtime
cache%use_mpi_runtime = params%use_mpi_runtime
cache%use_gradient_self_runtime = params%use_gradient_self_runtime
cache%use_self_runtime = params%use_self_runtime
cache%use_fd_gradient = params%use_fd_gradient
IF (ALLOCATED(cache%last_positions)) &
DEALLOCATE (cache%last_positions)
ALLOCATE (cache%last_positions(3, params%natom))
BLOCK
INTEGER :: iatom
DO iatom = 1, params%natom
cache%last_positions(:, iatom) = REAL(particle_set(iatom)%r(:), c_double)
END DO
END BLOCK
cache%mpi_comm = para_env%get_handle()
cache%is_init = .TRUE.
RETURN
! Cleanup section in the error case.
100 CONTINUE
CALL gauxc_destroy_molecule(cache%molecule, status)
CALL gauxc_check_status(status)
GO TO 999
200 CONTINUE
CALL gauxc_destroy_basisset(cache%basisset, status)
CALL gauxc_check_status(status)
GO TO 100
300 CONTINUE
CALL gauxc_destroy_grid(cache%grid, status)
CALL gauxc_check_status(status)
GO TO 200
400 CONTINUE
CALL gauxc_destroy_integrator(cache%integrator, status)
CALL gauxc_check_status(status)
GO TO 300
500 CONTINUE
CALL gauxc_destroy_grid(cache%gradient_grid, status)
CALL gauxc_check_status(status)
GO TO 400
600 CONTINUE
CALL gauxc_destroy_integrator(cache%gradient_integrator, status)
CALL gauxc_check_status(status)
999 CONTINUE
IF (ALLOCATED(cache%last_positions)) DEALLOCATE (cache%last_positions)
cache%is_init = .FALSE.
#else
MARK_USED(params)
MARK_USED(para_env)
MARK_USED(particle_set)
MARK_USED(qs_kind_set)
MARK_USED(status)
cache%is_init = .FALSE.
#endif
END SUBROUTINE gauxc_cache_init
! **************************************************************************************************
!> \brief Release all GauXC objects in a cache and deallocate arrays
!> \param cache ...
! **************************************************************************************************
SUBROUTINE gauxc_cache_release(cache)
TYPE(cp_gauxc_cache_type), INTENT(INOUT) :: cache
#ifdef __GAUXC
TYPE(cp_gauxc_status_type) :: status
IF (cache%is_init) THEN
IF (cache%needs_gradient_grid()) THEN
CALL gauxc_destroy_integrator(cache%gradient_integrator, status)
CALL gauxc_check_status(status)
CALL gauxc_destroy_grid(cache%gradient_grid, status)
CALL gauxc_check_status(status)
END IF
CALL gauxc_destroy_integrator(cache%integrator, status)
CALL gauxc_check_status(status)
CALL gauxc_destroy_grid(cache%grid, status)
CALL gauxc_check_status(status)
CALL gauxc_destroy_basisset(cache%basisset, status)
CALL gauxc_check_status(status)
CALL gauxc_destroy_molecule(cache%molecule, status)
CALL gauxc_check_status(status)
END IF
IF (ALLOCATED(cache%last_positions)) DEALLOCATE (cache%last_positions)
cache%is_init = .FALSE.
#else
MARK_USED(cache)
IF (ALLOCATED(cache%last_positions)) DEALLOCATE (cache%last_positions)
cache%is_init = .FALSE.
#endif
END SUBROUTINE gauxc_cache_release
! **************************************************************************************************
!> \brief Check if the GauXC cache is valid for the given parameters
!> \param cache ...
!> \param params ...
!> \param para_env ...
!> \param particle_set ...
!> \return ...
!> \retval cache_valid ...
! **************************************************************************************************
FUNCTION gauxc_cache_is_valid( &
cache, &
params, &
para_env, &
particle_set) &
RESULT(cache_valid)
TYPE(cp_gauxc_cache_type), INTENT(IN) :: cache
TYPE(cp_gauxc_cache_params), INTENT(in) :: params
TYPE(mp_para_env_type), INTENT(in), POINTER :: para_env
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
LOGICAL :: use_gradient_self_runtime, use_gradient_mpi_runtime, use_fd_gradient
LOGICAL :: cache_valid
#ifdef __GAUXC
INTEGER :: iatom
REAL(c_double), PARAMETER :: pos_tol = 1.0E-8_c_double
cache_valid = .FALSE.
IF (.NOT. cache%is_init) RETURN
IF (cache%natom /= params%natom) RETURN
IF (cache%nspins /= params%nspins) RETURN
IF (cache%batch_size /= params%batch_size) RETURN
IF (ABS(cache%device_runtime_fill_fraction - REAL(params%device_runtime_fill_fraction, c_double)) > &
1.0E-10_c_double) RETURN
IF (TRIM(cache%xc_fun_name) /= TRIM(params%xc_fun_name)) RETURN
IF (TRIM(cache%grid_type) /= TRIM(params%grid_type)) RETURN
IF (TRIM(cache%radial_quadrature) /= TRIM(params%radial_quadrature)) RETURN
IF (TRIM(cache%pruning_scheme) /= TRIM(params%pruning_scheme)) RETURN
IF (TRIM(cache%lb_exec_space) /= TRIM(params%lb_exec_space)) RETURN
IF (TRIM(cache%int_exec_space) /= TRIM(params%int_exec_space)) RETURN
IF (TRIM(cache%model_eval_name) /= TRIM(params%model_eval_name)) RETURN
IF (TRIM(cache%lwd_kernel) /= TRIM(params%lwd_kernel)) RETURN
IF (cache%use_self_runtime .NEQV. params%use_self_runtime) RETURN
IF (cache%use_mpi_runtime .NEQV. params%use_mpi_runtime) RETURN
! These three are modified during cache init, so we must emulate these changes here
use_gradient_self_runtime = params%use_gradient_self_runtime
use_gradient_mpi_runtime = params%use_gradient_mpi_runtime
use_fd_gradient = params%use_fd_gradient
use_fd_gradient = &
use_fd_gradient .AND. (cache%max_l() > 3)
use_gradient_mpi_runtime = &
use_gradient_mpi_runtime .AND. .NOT. use_fd_gradient
use_gradient_self_runtime = &
use_gradient_self_runtime .AND. .NOT. use_gradient_mpi_runtime
IF (cache%use_gradient_mpi_runtime .NEQV. use_gradient_mpi_runtime) RETURN
IF (cache%use_fd_gradient .NEQV. use_fd_gradient) RETURN
IF (cache%use_gradient_self_runtime .NEQV. use_gradient_self_runtime) RETURN
IF (cache%mpi_comm /= para_env%get_handle()) RETURN
IF (.NOT. ALLOCATED(cache%last_positions)) RETURN
IF (SIZE(cache%last_positions, 2) /= params%natom) RETURN
DO iatom = 1, params%natom
IF (ANY(ABS(cache%last_positions(:, iatom) - &
REAL(particle_set(iatom)%r(:), c_double)) > pos_tol)) RETURN
END DO
cache_valid = .TRUE.
#else
MARK_USED(cache)
MARK_USED(particle_set)
MARK_USED(params)
MARK_USED(use_gradient_self_runtime)
MARK_USED(use_gradient_mpi_runtime)
MARK_USED(use_fd_gradient)
MARK_USED(para_env)
cache_valid = .FALSE.
#endif
END FUNCTION gauxc_cache_is_valid
END MODULE xc_gauxc_cache

397
src/xc/xc_gauxc_functional.F
View file

@ -49,6 +49,9 @@ MODULE xc_gauxc_functional
qs_rho_type
USE qs_scf_types, ONLY: qs_scf_env_type
USE string_utilities, ONLY: uppercase
USE xc_gauxc_cache, ONLY: cp_gauxc_cache_params,&
cp_gauxc_cache_type,&
gauxc_cache_init
USE xc_gauxc_interface, ONLY: &
cp_gauxc_basisset_type, cp_gauxc_grid_type, cp_gauxc_integrator_type, &
cp_gauxc_molecule_type, cp_gauxc_status_type, cp_gauxc_xc_gradient_type, cp_gauxc_xc_type, &
@ -117,24 +120,26 @@ CONTAINS
!> \brief ...
!> \param dbcsr_mat ...
!> \param dense_mat ...
!> \param para_env ...
! **************************************************************************************************
SUBROUTINE dbcsr_to_dense(dbcsr_mat, dense_mat)
USE cp_dbcsr_api, ONLY: dbcsr_get_info, dbcsr_get_matrix_type, dbcsr_iterator_type, &
dbcsr_iterator_start, dbcsr_iterator_next_block, dbcsr_iterator_blocks_left, &
dbcsr_iterator_stop, dbcsr_type_antisymmetric, dbcsr_type_symmetric
SUBROUTINE dbcsr_to_dense(dbcsr_mat, dense_mat, para_env)
USE cp_dbcsr_api, ONLY: dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_get_info, &
dbcsr_get_matrix_type, dbcsr_get_readonly_block_p, &
dbcsr_get_stored_coordinates, dbcsr_type_antisymmetric, &
dbcsr_type_symmetric
TYPE(dbcsr_p_type), INTENT(IN) :: dbcsr_mat
REAL(c_double), ALLOCATABLE, DIMENSION(:, :), &
INTENT(INOUT) :: dense_mat
TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
CHARACTER :: matrix_type
INTEGER :: col, col_end, col_start, icol, irow, &
nblkcols_total, nblkrows_total, ncols, &
nrows, row, row_end, row_start
INTEGER :: col, col_end, col_start, icol, irow, mynode, nblkcols_total, nblkrows_total, &
ncols, nrows, numnodes, owner, row, row_end, row_start
INTEGER, ALLOCATABLE, DIMENSION(:) :: c_offset, r_offset
INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
LOGICAL :: transposed
LOGICAL :: found
REAL(c_double), POINTER :: block(:, :)
TYPE(dbcsr_iterator_type) :: iter
TYPE(dbcsr_distribution_type) :: dist
CALL dbcsr_get_info(dbcsr_mat%matrix, &
row_blk_size=row_blk_size, &
@ -142,7 +147,9 @@ CONTAINS
nblkrows_total=nblkrows_total, &
nblkcols_total=nblkcols_total, &
nfullrows_total=nrows, &
nfullcols_total=ncols)
nfullcols_total=ncols, &
distribution=dist)
CALL dbcsr_distribution_get(dist, mynode=mynode, numnodes=numnodes)
matrix_type = dbcsr_get_matrix_type(dbcsr_mat%matrix)
IF (.NOT. ALLOCATED(dense_mat)) THEN
@ -166,23 +173,19 @@ CONTAINS
c_offset(col) = c_offset(col - 1) + col_blk_size(col - 1)
END DO
CALL dbcsr_iterator_start(iter, dbcsr_mat%matrix)
DO WHILE (dbcsr_iterator_blocks_left(iter))
CALL dbcsr_iterator_next_block(iter, irow, icol, block, transposed=transposed)
row_start = r_offset(irow)
row_end = row_start + row_blk_size(irow) - 1
col_start = c_offset(icol)
col_end = col_start + col_blk_size(icol) - 1
IF (transposed) THEN
dense_mat(row_start:row_end, col_start:col_end) = TRANSPOSE(block)
IF (irow /= icol) THEN
IF (matrix_type == dbcsr_type_symmetric) THEN
dense_mat(col_start:col_end, row_start:row_end) = block
ELSE IF (matrix_type == dbcsr_type_antisymmetric) THEN
dense_mat(col_start:col_end, row_start:row_end) = -block
END IF
END IF
ELSE
! Replicated DBCSR blocks must enter the following MPI sum exactly once.
DO irow = 1, nblkrows_total
DO icol = 1, nblkcols_total
IF (numnodes == 1 .AND. para_env%num_pe > 1 .AND. para_env%mepos /= 0) CYCLE
CALL dbcsr_get_stored_coordinates(dbcsr_mat%matrix, irow, icol, owner)
IF (owner /= mynode) CYCLE
CALL dbcsr_get_readonly_block_p(matrix=dbcsr_mat%matrix, row=irow, col=icol, &
block=block, found=found)
IF (.NOT. found) CYCLE
row_start = r_offset(irow)
row_end = row_start + row_blk_size(irow) - 1
col_start = c_offset(icol)
col_end = col_start + col_blk_size(icol) - 1
dense_mat(row_start:row_end, col_start:col_end) = block
IF (irow /= icol) THEN
IF (matrix_type == dbcsr_type_symmetric) THEN
@ -191,9 +194,8 @@ CONTAINS
dense_mat(col_start:col_end, row_start:row_end) = -TRANSPOSE(block)
END IF
END IF
END IF
END DO
END DO
CALL dbcsr_iterator_stop(iter)
DEALLOCATE (r_offset, c_offset)
@ -572,7 +574,8 @@ CONTAINS
batch_size, &
device_runtime_fill_fraction, &
gauxc_status, &
mpi_comm=mp_comm_self%get_handle())
mpi_comm=mp_comm_self%get_handle(), &
force_new_runtime=.TRUE.)
CALL gauxc_check_status(gauxc_status)
gauxc_integrator_fd = gauxc_create_integrator( &
TRIM(xc_fun_name), &
@ -895,31 +898,36 @@ CONTAINS
INTEGER, INTENT(out), OPTIONAL :: max_deriv
CHARACTER(len=default_path_length) :: model_key, model_name
CHARACTER(len=default_string_length) :: xc_fun_name
CHARACTER(len=default_string_length) :: xc_fun_key, xc_fun_name
LOGICAL :: native_grid
CALL section_vals_val_get(functional, "FUNCTIONAL", c_val=xc_fun_name)
CALL section_vals_val_get(functional, "MODEL", c_val=model_name)
CALL section_vals_val_get(functional, "NATIVE_GRID", l_val=native_grid)
model_key = ADJUSTL(model_name)
xc_fun_key = ADJUSTL(xc_fun_name)
CALL uppercase(model_key)
CALL uppercase(xc_fun_key)
IF (PRESENT(reference)) THEN
IF (TRIM(model_key) == "NONE" .OR. TRIM(model_key) == "") THEN
IF (TRIM(model_key) == "NONE" .OR. TRIM(model_key) == "" .OR. &
TRIM(model_key) == TRIM(xc_fun_key)) THEN
reference = "Functional computed by GauXC (underlying: "//TRIM(xc_fun_name)//")"
ELSE
reference = "Functional computed by GauXC Skala model "//TRIM(model_name)
END IF
END IF
IF (PRESENT(shortform)) THEN
IF (TRIM(model_key) == "NONE" .OR. TRIM(model_key) == "") THEN
IF (TRIM(model_key) == "NONE" .OR. TRIM(model_key) == "" .OR. &
TRIM(model_key) == TRIM(xc_fun_key)) THEN
shortform = "GAUXC ("//TRIM(xc_fun_name)//")"
ELSE
shortform = "GAUXC Skala"
END IF
END IF
IF (PRESENT(needs)) THEN
IF (native_grid .AND. TRIM(model_key) /= "NONE" .AND. TRIM(model_key) /= "") THEN
IF (native_grid .AND. TRIM(model_key) /= "NONE" .AND. TRIM(model_key) /= "" .AND. &
TRIM(model_key) /= TRIM(xc_fun_key)) THEN
IF (lsd) THEN
needs%rho_spin = .TRUE.
needs%drho_spin = .TRUE.
@ -961,26 +969,21 @@ CONTAINS
REAL(KIND=dp), PARAMETER :: gapw_fd_gradient_dx = 1.0E-4_dp
CHARACTER(len=default_path_length) :: model_key, model_name, output_path
CHARACTER(len=default_string_length) :: gradient_runtime, gradient_runtime_key, grid_key, &
grid_type, int_exec_space, lb_exec_space, lwd_kernel, pruning_key, pruning_scheme, &
radial_quadrature, skala_runtime, skala_runtime_key, xc_fun_name
INTEGER :: atom_chunk_size, batch_size, env_status, &
img, ispin, natom, nimages, nspins
CHARACTER(len=default_string_length) :: gradient_runtime, gradient_runtime_key, &
grid_key, pruning_key, skala_runtime, &
skala_runtime_key, xc_fun_key
INTEGER :: atom_chunk_size, env_status, img, ispin, &
nimages
LOGICAL :: atom_chunk_size_explicit, do_kpoints, gapw_method, gapw_paw_pseudopotentials, &
gapw_pseudopotentials, grid_explicit, hdf5_output, is_periodic, molecular_virial, &
molecular_virial_debug, need_xc_gradient, periodic_reference, pruning_explicit, &
use_fd_gradient, use_gauxc_model, use_gradient_mpi_runtime, use_gradient_self_runtime, &
use_self_runtime, use_skala_model, write_hdf5_output
REAL(KIND=dp) :: device_runtime_fill_fraction, &
molecular_virial_debug_dx
use_skala_model, write_hdf5_output
REAL(KIND=dp) :: molecular_virial_debug_dx
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: density_scalar, density_zeta
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cp_gauxc_basisset_type) :: gauxc_basis
TYPE(cp_gauxc_grid_type) :: gauxc_gradient_grid_result, &
gauxc_grid_result
TYPE(cp_gauxc_integrator_type) :: gauxc_gradient_integrator_result, gauxc_integrator_result
TYPE(cp_gauxc_molecule_type) :: gauxc_mol
TYPE(cp_gauxc_cache_params) :: params
TYPE(cp_gauxc_cache_type), POINTER :: cache
TYPE(cp_gauxc_status_type) :: gauxc_status
TYPE(cp_gauxc_xc_gradient_type) :: exc_grad
TYPE(cp_gauxc_xc_type) :: gauxc_xc_result
@ -1023,7 +1026,7 @@ CONTAINS
energy=energy, &
ks_env=ks_env, &
matrix_vxc=matrix_vxc, &
natom=natom, &
natom=params%natom, &
atomic_kind_set=atomic_kind_set, &
force=force, &
para_env=para_env, &
@ -1050,7 +1053,7 @@ CONTAINS
rho_ao_kp=rho_ao)
nimages = dft_control%nimages
nspins = dft_control%nspins
params%nspins = dft_control%nspins
is_periodic = .FALSE.
IF (ASSOCIATED(cell)) is_periodic = ANY(cell%perd /= 0)
@ -1065,7 +1068,7 @@ CONTAINS
CALL section_vals_val_get( &
gauxc_functional_section, &
"FUNCTIONAL", &
c_val=xc_fun_name)
c_val=params%xc_fun_name)
CALL section_vals_val_get( &
gauxc_functional_section, &
"MODEL", &
@ -1073,25 +1076,25 @@ CONTAINS
CALL section_vals_val_get( &
gauxc_functional_section, &
"GRID", &
c_val=grid_type, &
c_val=params%grid_type, &
explicit=grid_explicit)
CALL section_vals_val_get( &
gauxc_functional_section, &
"RADIAL_QUADRATURE", &
c_val=radial_quadrature)
c_val=params%radial_quadrature)
CALL section_vals_val_get( &
gauxc_functional_section, &
"PRUNING_SCHEME", &
c_val=pruning_scheme, &
c_val=params%pruning_scheme, &
explicit=pruning_explicit)
CALL section_vals_val_get( &
gauxc_functional_section, &
"BATCH_SIZE", &
i_val=batch_size)
i_val=params%batch_size)
CALL section_vals_val_get( &
gauxc_functional_section, &
"DEVICE_RUNTIME_FILL_FRACTION", &
r_val=device_runtime_fill_fraction)
r_val=params%device_runtime_fill_fraction)
CALL section_vals_val_get( &
gauxc_functional_section, &
"MODEL_ATOM_CHUNK_SIZE", &
@ -1116,15 +1119,15 @@ CONTAINS
CALL section_vals_val_get( &
gauxc_functional_section, &
"LB_EXECUTION_SPACE", &
c_val=lb_exec_space)
c_val=params%lb_exec_space)
CALL section_vals_val_get( &
gauxc_functional_section, &
"INT_EXECUTION_SPACE", &
c_val=int_exec_space)
c_val=params%int_exec_space)
CALL section_vals_val_get( &
gauxc_functional_section, &
"LWD_KERNEL", &
c_val=lwd_kernel)
c_val=params%lwd_kernel)
CALL section_vals_val_get( &
gauxc_functional_section, &
"SKALA_RUNTIME", &
@ -1140,20 +1143,21 @@ CONTAINS
model_key = ADJUSTL(model_name)
CALL uppercase(model_key)
xc_fun_key = ADJUSTL(params%xc_fun_name)
CALL uppercase(xc_fun_key)
skala_runtime_key = ADJUSTL(skala_runtime)
CALL uppercase(skala_runtime_key)
gradient_runtime_key = ADJUSTL(gradient_runtime)
CALL uppercase(gradient_runtime_key)
use_gauxc_model = (TRIM(model_key) /= "" .AND. TRIM(model_key) /= "NONE")
params%use_gauxc_model = (TRIM(model_key) /= "" .AND. TRIM(model_key) /= "NONE" .AND. &
TRIM(model_key) /= TRIM(xc_fun_key))
use_skala_model = (INDEX(TRIM(model_key), "SKALA") > 0)
IF (gapw_pseudopotentials .AND. .NOT. use_gauxc_model) THEN
CALL cp_abort(__LOCATION__, &
"GauXC with METHOD GAPW/GAPW_XC and pseudopotentials is supported only for "// &
"Skala-style models that replace the molecular XC term. "// &
"Use POTENTIAL ALL for local/semi-local GauXC GAPW validation or METHOD GPW "// &
"with pseudopotentials.")
params%model_eval_name = model_name
IF (.NOT. params%use_gauxc_model) THEN
! MODEL NONE and MODEL equal to FUNCTIONAL select conventional GauXC.
params%model_eval_name = "NONE"
END IF
IF (gapw_pseudopotentials .AND. use_gauxc_model .AND. .NOT. dft_control%qs_control%gapw_xc .AND. &
IF (gapw_pseudopotentials .AND. params%use_gauxc_model .AND. .NOT. dft_control%qs_control%gapw_xc .AND. &
.NOT. gapw_paw_pseudopotentials .AND. para_env%mepos == 0 .AND. ASSOCIATED(scf_env)) THEN
IF (scf_env%iter_count == 1) THEN
CALL cp_warn( &
@ -1163,7 +1167,7 @@ CONTAINS
"XC correction is used for those regular-grid kinds.")
END IF
END IF
IF (device_runtime_fill_fraction <= 0.0_dp .OR. device_runtime_fill_fraction > 1.0_dp) THEN
IF (params%device_runtime_fill_fraction <= 0.0_dp .OR. params%device_runtime_fill_fraction > 1.0_dp) THEN
CALL cp_abort(__LOCATION__, &
"GAUXC%DEVICE_RUNTIME_FILL_FRACTION must be > 0 and <= 1.")
END IF
@ -1191,21 +1195,21 @@ CONTAINS
END IF
END IF
END IF
IF (use_gauxc_model) THEN
IF (params%use_gauxc_model) THEN
IF (has_nlcc(qs_kind_set)) THEN
CALL cp_abort(__LOCATION__, &
"GauXC Skala with NLCC pseudopotentials is not implemented. "// &
"The frozen core density would need a SKALA-consistent feature definition.")
END IF
END IF
IF (use_gauxc_model) THEN
IF (params%use_gauxc_model) THEN
CALL set_gauxc_model_atom_chunk_env( &
atom_chunk_size, atom_chunk_size_explicit)
IF (.NOT. grid_explicit) grid_type = "SUPERFINE"
IF (.NOT. pruning_explicit) pruning_scheme = "UNPRUNED"
IF (.NOT. grid_explicit) params%grid_type = "SUPERFINE"
IF (.NOT. pruning_explicit) params%pruning_scheme = "UNPRUNED"
grid_key = ADJUSTL(grid_type)
pruning_key = ADJUSTL(pruning_scheme)
grid_key = ADJUSTL(params%grid_type)
pruning_key = ADJUSTL(params%pruning_scheme)
CALL uppercase(grid_key)
CALL uppercase(pruning_key)
IF (use_skala_model .AND. need_xc_gradient .AND. &
@ -1220,35 +1224,36 @@ CONTAINS
IF (env_status /= 0 .OR. LEN_TRIM(model_name) == 0) THEN
CPABORT("MODEL SKALA requires the GAUXC_SKALA_MODEL environment variable")
END IF
params%model_eval_name = model_name
END IF
END IF
SELECT CASE (TRIM(skala_runtime_key))
CASE ("AUTO")
use_self_runtime = use_skala_model .AND. para_env%num_pe > 1 .AND. nspins > 1
params%use_self_runtime = use_skala_model .AND. para_env%num_pe > 1 .AND. params%nspins > 1
CASE ("MPI")
use_self_runtime = .FALSE.
params%use_self_runtime = .FALSE.
CASE ("SELF")
use_self_runtime = use_skala_model .AND. para_env%num_pe > 1
params%use_self_runtime = use_skala_model .AND. para_env%num_pe > 1
CASE DEFAULT
CALL cp_abort(__LOCATION__, "Unknown GAUXC%SKALA_RUNTIME value.")
END SELECT
IF (.NOT. use_skala_model) use_self_runtime = .FALSE.
IF (.NOT. use_skala_model) params%use_self_runtime = .FALSE.
SELECT CASE (TRIM(gradient_runtime_key))
CASE ("AUTO", "SELF")
use_gradient_mpi_runtime = .FALSE.
use_gradient_self_runtime = need_xc_gradient .AND. use_gauxc_model .AND. &
para_env%num_pe > 1 .AND. .NOT. use_self_runtime
params%use_gradient_mpi_runtime = .FALSE.
params%use_gradient_self_runtime = need_xc_gradient .AND. params%use_gauxc_model .AND. &
para_env%num_pe > 1 .AND. .NOT. params%use_self_runtime
CASE ("MPI")
use_gradient_mpi_runtime = need_xc_gradient .AND. use_gauxc_model .AND. para_env%num_pe > 1
use_gradient_self_runtime = .FALSE.
params%use_gradient_mpi_runtime = need_xc_gradient .AND. params%use_gauxc_model .AND. para_env%num_pe > 1
params%use_gradient_self_runtime = .FALSE.
CASE DEFAULT
CALL cp_abort(__LOCATION__, "Unknown GAUXC%MODEL_GRADIENT_RUNTIME value.")
END SELECT
IF (.NOT. use_gauxc_model) THEN
use_gradient_mpi_runtime = .FALSE.
use_gradient_self_runtime = .FALSE.
IF (.NOT. params%use_gauxc_model) THEN
params%use_gradient_mpi_runtime = .FALSE.
params%use_gradient_self_runtime = .FALSE.
END IF
IF (use_skala_model .AND. para_env%num_pe > 1 .AND. .NOT. use_self_runtime .AND. &
IF (use_skala_model .AND. para_env%num_pe > 1 .AND. .NOT. params%use_self_runtime .AND. &
para_env%mepos == 0 .AND. ASSOCIATED(scf_env)) THEN
IF (scf_env%iter_count == 1) THEN
CALL cp_warn( &
@ -1259,15 +1264,18 @@ CONTAINS
END IF
END IF
gauxc_mol = gauxc_create_molecule( &
particle_set, &
gauxc_status)
CALL gauxc_check_status(gauxc_status)
gauxc_basis = gauxc_create_basisset( &
qs_kind_set, &
particle_set, &
gauxc_status)
CALL gauxc_check_status(gauxc_status)
! After creating the basisset, we will have to check max_l>3 as a further condition
params%use_fd_gradient = gapw_method .AND. need_xc_gradient
cache => qs_env%gauxc_cache
CALL gauxc_cache_init( &
cache, &
params, &
para_env, &
particle_set, &
qs_kind_set, &
gauxc_status)
hdf5_output = (TRIM(output_path) /= "")
write_hdf5_output = hdf5_output .AND. para_env%mepos == 0
IF (write_hdf5_output .AND. ASSOCIATED(scf_env)) THEN
@ -1275,98 +1283,20 @@ CONTAINS
END IF
IF (write_hdf5_output) THEN
CALL gauxc_write_molecule_hdf5( &
gauxc_mol, &
cache%molecule, &
output_path, &
"molecule.h5", &
"molecule", &
gauxc_status)
CALL gauxc_check_status(gauxc_status)
CALL gauxc_write_basisset_hdf5( &
gauxc_basis, &
cache%basisset, &
output_path, &
"basisset.h5", &
"basisset", &
gauxc_status)
CALL gauxc_check_status(gauxc_status)
END IF
use_fd_gradient = gapw_method .AND. need_xc_gradient .AND. (gauxc_basis%max_l > 3)
IF (use_fd_gradient) use_gradient_mpi_runtime = .FALSE.
IF (use_gradient_mpi_runtime) THEN
use_gradient_self_runtime = .FALSE.
END IF
IF (use_fd_gradient .AND. para_env%mepos == 0) THEN
CALL cp_warn( &
__LOCATION__, &
"Using finite-difference GauXC XC gradients for METHOD GAPW/GAPW_XC with "// &
"basis functions beyond f shells. The upstream analytical GauXC gradient path is "// &
"not yet reliable for this case.")
END IF
IF (use_self_runtime) THEN
! SKALA currently needs a replicated molecular runtime for reproducible
! open-shell densities across CP2K MPI ranks.
gauxc_grid_result = gauxc_create_grid( &
gauxc_mol, &
gauxc_basis, &
grid_type, &
radial_quadrature, &
pruning_scheme, &
lb_exec_space, &
batch_size, &
device_runtime_fill_fraction, &
gauxc_status, &
mpi_comm=mp_comm_self%get_handle())
ELSE
! Use the QS force-evaluation communicator rather than GauXC's global
! runtime. In mixed CDFT the diabatic states can be built on disjoint
! MPI subgroups; using the global communicator would make GauXC wait
! for ranks that are working on another state.
gauxc_grid_result = gauxc_create_grid( &
gauxc_mol, &
gauxc_basis, &
grid_type, &
radial_quadrature, &
pruning_scheme, &
lb_exec_space, &
batch_size, &
device_runtime_fill_fraction, &
gauxc_status, &
mpi_comm=para_env%get_handle())
END IF
CALL gauxc_check_status(gauxc_status)
gauxc_integrator_result = gauxc_create_integrator( &
TRIM(xc_fun_name), &
gauxc_grid_result, &
int_exec_space, &
lwd_kernel, &
nspins, &
gauxc_status)
CALL gauxc_check_status(gauxc_status)
IF (use_gradient_self_runtime) THEN
! Upstream GauXC does not yet support Skala nuclear gradients
! on an MPI runtime. Keep the energy/VXC path on the normal MPI
! runtime and use an isolated runtime only for the replicated gradient.
gauxc_gradient_grid_result = gauxc_create_grid( &
gauxc_mol, &
gauxc_basis, &
grid_type, &
radial_quadrature, &
pruning_scheme, &
lb_exec_space, &
batch_size, &
device_runtime_fill_fraction, &
gauxc_status, &
mpi_comm=mp_comm_self%get_handle())
CALL gauxc_check_status(gauxc_status)
gauxc_gradient_integrator_result = gauxc_create_integrator( &
TRIM(xc_fun_name), &
gauxc_gradient_grid_result, &
int_exec_space, &
lwd_kernel, &
nspins, &
gauxc_status)
CALL gauxc_check_status(gauxc_status)
END IF
IF (qs_env%run_rtp) THEN
CPABORT("GAUXC XC energy currently does not support real-time propagation")
@ -1375,46 +1305,46 @@ CONTAINS
energy%exc = 0
IF (ASSOCIATED(matrix_vxc)) CALL dbcsr_deallocate_matrix_set(matrix_vxc)
CALL dbcsr_allocate_matrix_set(matrix_vxc, nspins)
CALL dbcsr_allocate_matrix_set(matrix_vxc, params%nspins)
DO img = 1, nimages
IF (img > 1) THEN
CPABORT("UNIMPLEMENTED: Handling nimg>1 in k-point integration")
END IF
CALL dbcsr_to_dense(rho_ao(1, img), density_scalar)
CALL dbcsr_to_dense(rho_ao(1, img), density_scalar, para_env)
CALL para_env%sum(density_scalar)
IF (nspins == 1) THEN
IF (params%nspins == 1) THEN
gauxc_xc_result = gauxc_compute_xc( &
gauxc_integrator_result, &
cache%integrator, &
density_scalar, &
nspins=nspins, &
nspins=params%nspins, &
status=gauxc_status, &
model=TRIM(model_name))
model=TRIM(params%model_eval_name))
CALL gauxc_check_status(gauxc_status)
IF (need_xc_gradient) THEN
IF (use_fd_gradient) THEN
IF (params%use_fd_gradient) THEN
CALL gauxc_xc_gradient_fd( &
particle_set, qs_kind_set, density_scalar, nspins, model_name, &
xc_fun_name, grid_type, radial_quadrature, pruning_scheme, &
lb_exec_space, int_exec_space, lwd_kernel, batch_size, &
device_runtime_fill_fraction, gapw_fd_gradient_dx, para_env, &
particle_set, qs_kind_set, density_scalar, params%nspins, params%model_eval_name, &
params%xc_fun_name, params%grid_type, params%radial_quadrature, params%pruning_scheme, &
params%lb_exec_space, params%int_exec_space, params%lwd_kernel, params%batch_size, &
params%device_runtime_fill_fraction, gapw_fd_gradient_dx, para_env, &
exc_grad%exc_grad)
ELSE IF (use_gradient_self_runtime) THEN
ELSE IF (params%use_gradient_self_runtime) THEN
exc_grad = gauxc_compute_xc_gradient( &
gauxc_gradient_integrator_result, &
cache%gradient_integrator, &
density_scalar, &
nspins=nspins, &
natom=natom, &
nspins=params%nspins, &
natom=params%natom, &
status=gauxc_status, &
model=TRIM(model_name))
model=TRIM(params%model_eval_name))
ELSE
exc_grad = gauxc_compute_xc_gradient( &
gauxc_integrator_result, &
cache%integrator, &
density_scalar, &
nspins=nspins, &
natom=natom, &
nspins=params%nspins, &
natom=params%natom, &
status=gauxc_status, &
model=TRIM(model_name))
model=TRIM(params%model_eval_name))
END IF
CALL gauxc_check_status(gauxc_status)
IF (calculate_forces) THEN
@ -1429,19 +1359,19 @@ CONTAINS
END IF
IF (molecular_virial_debug) THEN
CALL debug_gauxc_molecular_virial( &
exc_grad%exc_grad, particle_set, qs_kind_set, density_scalar, nspins, &
model_name, xc_fun_name, grid_type, radial_quadrature, pruning_scheme, &
lb_exec_space, int_exec_space, lwd_kernel, batch_size, &
device_runtime_fill_fraction, molecular_virial_debug_dx, para_env)
exc_grad%exc_grad, particle_set, qs_kind_set, density_scalar, params%nspins, &
params%model_eval_name, params%xc_fun_name, params%grid_type, params%radial_quadrature, params%pruning_scheme, &
params%lb_exec_space, params%int_exec_space, params%lwd_kernel, params%batch_size, &
params%device_runtime_fill_fraction, molecular_virial_debug_dx, para_env)
END IF
DEALLOCATE (exc_grad%exc_grad)
END IF
ELSE
CPASSERT(nspins == 2)
CPASSERT(params%nspins == 2)
! In here:
! scalar <- rho_ao(1, :) + rho_ao(2, :)
! zeta <- rho_ao(1, :) - rho_ao(2, :)
CALL dbcsr_to_dense(rho_ao(2, img), density_zeta)
CALL dbcsr_to_dense(rho_ao(2, img), density_zeta, para_env)
CALL para_env%sum(density_zeta)
! Do NOT reorder the following lines!
density_scalar(:, :) = density_scalar(:, :) + density_zeta(:, :)
@ -1451,39 +1381,39 @@ CONTAINS
! This style lowers memory footprint.
density_zeta(:, :) = density_scalar(:, :) - 2.0_dp*density_zeta(:, :)
gauxc_xc_result = gauxc_compute_xc( &
gauxc_integrator_result, &
cache%integrator, &
density_scalar, &
density_zeta, &
nspins, &
params%nspins, &
gauxc_status, &
model=TRIM(model_name))
model=TRIM(params%model_eval_name))
CALL gauxc_check_status(gauxc_status)
IF (need_xc_gradient) THEN
IF (use_fd_gradient) THEN
IF (params%use_fd_gradient) THEN
CALL gauxc_xc_gradient_fd( &
particle_set, qs_kind_set, density_scalar, nspins, model_name, &
xc_fun_name, grid_type, radial_quadrature, pruning_scheme, &
lb_exec_space, int_exec_space, lwd_kernel, batch_size, &
device_runtime_fill_fraction, gapw_fd_gradient_dx, para_env, &
particle_set, qs_kind_set, density_scalar, params%nspins, params%model_eval_name, &
params%xc_fun_name, params%grid_type, params%radial_quadrature, params%pruning_scheme, &
params%lb_exec_space, params%int_exec_space, params%lwd_kernel, params%batch_size, &
params%device_runtime_fill_fraction, gapw_fd_gradient_dx, para_env, &
exc_grad%exc_grad, density_zeta=density_zeta)
ELSE IF (use_gradient_self_runtime) THEN
ELSE IF (params%use_gradient_self_runtime) THEN
exc_grad = gauxc_compute_xc_gradient( &
gauxc_gradient_integrator_result, &
cache%gradient_integrator, &
density_scalar, &
density_zeta, &
nspins, &
natom, &
params%nspins, &
params%natom, &
gauxc_status, &
model=TRIM(model_name))
model=TRIM(params%model_eval_name))
ELSE
exc_grad = gauxc_compute_xc_gradient( &
gauxc_integrator_result, &
cache%integrator, &
density_scalar, &
density_zeta, &
nspins, &
natom, &
params%nspins, &
params%natom, &
gauxc_status, &
model=TRIM(model_name))
model=TRIM(params%model_eval_name))
END IF
CALL gauxc_check_status(gauxc_status)
IF (calculate_forces) THEN
@ -1498,10 +1428,10 @@ CONTAINS
END IF
IF (molecular_virial_debug) THEN
CALL debug_gauxc_molecular_virial( &
exc_grad%exc_grad, particle_set, qs_kind_set, density_scalar, nspins, &
model_name, xc_fun_name, grid_type, radial_quadrature, pruning_scheme, &
lb_exec_space, int_exec_space, lwd_kernel, batch_size, &
device_runtime_fill_fraction, molecular_virial_debug_dx, para_env, &
exc_grad%exc_grad, particle_set, qs_kind_set, density_scalar, params%nspins, &
params%model_eval_name, params%xc_fun_name, params%grid_type, params%radial_quadrature, params%pruning_scheme, &
params%lb_exec_space, params%int_exec_space, params%lwd_kernel, params%batch_size, &
params%device_runtime_fill_fraction, molecular_virial_debug_dx, para_env, &
density_zeta=density_zeta)
END IF
DEALLOCATE (exc_grad%exc_grad)
@ -1510,14 +1440,14 @@ CONTAINS
energy%exc = energy%exc + gauxc_xc_result%exc
IF (nspins == 1) THEN
IF (params%nspins == 1) THEN
IF (img == 1) THEN
matrix_vxc(1) = dense_to_dbcsr(gauxc_xc_result%vxc_scalar, rho_ao(1, img))
ELSE
CPABORT("UNIMPLEMENTED: Handling multiple result matrices in k-point integration")
END IF
ELSE
CPASSERT(nspins == 2)
CPASSERT(params%nspins == 2)
! Transform derivatives from total/spin density back to alpha/beta channels.
vxc_zeta_tmp = dense_to_dbcsr(gauxc_xc_result%vxc_zeta, rho_ao(1, img))
IF (img == 1) THEN
@ -1544,25 +1474,10 @@ CONTAINS
IF (ALLOCATED(gauxc_xc_result%vxc_zeta)) DEALLOCATE (gauxc_xc_result%vxc_zeta)
CALL set_ks_env(ks_env, matrix_vxc=matrix_vxc)
DO ispin = 1, nspins
DO ispin = 1, params%nspins
CALL dbcsr_finalize(matrix_vxc(ispin)%matrix)
END DO
IF (use_gradient_self_runtime) THEN
CALL gauxc_destroy_integrator(gauxc_gradient_integrator_result, gauxc_status)
CALL gauxc_check_status(gauxc_status)
CALL gauxc_destroy_grid(gauxc_gradient_grid_result, gauxc_status)
CALL gauxc_check_status(gauxc_status)
END IF
CALL gauxc_destroy_integrator(gauxc_integrator_result, gauxc_status)
CALL gauxc_check_status(gauxc_status)
CALL gauxc_destroy_grid(gauxc_grid_result, gauxc_status)
CALL gauxc_check_status(gauxc_status)
CALL gauxc_destroy_basisset(gauxc_basis, gauxc_status)
CALL gauxc_check_status(gauxc_status)
CALL gauxc_destroy_molecule(gauxc_mol, gauxc_status)
CALL gauxc_check_status(gauxc_status)
END SUBROUTINE apply_gauxc
END MODULE xc_gauxc_functional

18
src/xc/xc_gauxc_interface.F
View file

@ -506,6 +506,7 @@ CONTAINS
!> \param device_runtime_fill_fraction ...
!> \param status ...
!> \param mpi_comm optional communicator for a grid-local GauXC runtime
!> \param force_new_runtime force creation of a grid-local GauXC runtime
!> \return ...
! **************************************************************************************************
FUNCTION gauxc_create_grid( &
@ -518,7 +519,8 @@ CONTAINS
batch_size, &
device_runtime_fill_fraction, &
status, &
mpi_comm) RESULT(res)
mpi_comm, &
force_new_runtime) RESULT(res)
TYPE(cp_gauxc_molecule_type), INTENT(IN) :: molecule
TYPE(cp_gauxc_basisset_type), INTENT(in) :: basis
@ -528,13 +530,14 @@ CONTAINS
REAL(c_double), INTENT(IN) :: device_runtime_fill_fraction
TYPE(cp_gauxc_status_type), INTENT(OUT) :: status
INTEGER, INTENT(IN), OPTIONAL :: mpi_comm
LOGICAL, INTENT(IN), OPTIONAL :: force_new_runtime
TYPE(cp_gauxc_grid_type) :: res
#ifdef __GAUXC
INTEGER(c_int) :: grid_type_local, int_exec_space_local, &
lb_exec_space_local, &
pruning_scheme_local, radial_quad_local
LOGICAL :: use_device_runtime
LOGICAL :: force_new_runtime_local, use_device_runtime
grid_type_local = read_atomic_grid_size(grid_type)
radial_quad_local = read_radial_quad(radial_quadrature)
@ -542,6 +545,8 @@ CONTAINS
lb_exec_space_local = read_execution_space(lb_exec_space)
int_exec_space_local = read_execution_space("host")
use_device_runtime = (lb_exec_space_local == gauxc_executionspace%device)
force_new_runtime_local = .FALSE.
IF (PRESENT(force_new_runtime)) force_new_runtime_local = force_new_runtime
res%owns_rt = .FALSE.
IF (use_device_runtime) THEN
@ -570,7 +575,8 @@ CONTAINS
IF (PRESENT(mpi_comm)) THEN
! Reuse the global runtime when the requested communicator matches
! the communicator used during gauxc_init.
IF (.NOT. rt_has_mpi_comm .OR. mpi_comm /= rt_mpi_comm) THEN
IF (force_new_runtime_local .OR. .NOT. rt_has_mpi_comm .OR. &
mpi_comm /= rt_mpi_comm) THEN
res%rt = gauxc_runtime_environment_new(status%status, mpi_comm)
GAUXC_RETURN_IF_ERROR(status)
res%owns_rt = .TRUE.
@ -578,6 +584,7 @@ CONTAINS
END IF
#else
MARK_USED(mpi_comm)
MARK_USED(force_new_runtime)
#endif
END IF
@ -637,6 +644,7 @@ CONTAINS
MARK_USED(grid_type)
MARK_USED(lb_exec_space)
MARK_USED(mpi_comm)
MARK_USED(force_new_runtime)
MARK_USED(molecule)
MARK_USED(pruning_scheme)
MARK_USED(radial_quadrature)
@ -809,10 +817,6 @@ CONTAINS
END IF
IF (nspins == 1) THEN
! xmat factor 2 is applied by both CP2K and GauXC
! "unapply" it here to even things back out.
! This is NOT necessary in the Skala branch.
density_scalar = 0.5_dp*density_scalar
CALL gauxc_integrator_eval_exc_vxc_rks( &
status%status, &
integrator%integrator, &

66
src/xc/xc_libxc.F
View file

@ -307,9 +307,11 @@ CONTAINS
!> true (does not set the unneeded components to false)
!> \param max_deriv maximum implemented derivative of the xc functional
!> \param print_warn whether to print warning about development status of a functional
!> \param func_name_override optional LibXC functional name overriding the section name
!> \author F. Tran
! **************************************************************************************************
SUBROUTINE libxc_lda_info(libxc_params, reference, shortform, needs, max_deriv, print_warn)
SUBROUTINE libxc_lda_info(libxc_params, reference, shortform, needs, max_deriv, print_warn, &
func_name_override)
TYPE(section_vals_type), POINTER :: libxc_params
CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
@ -317,6 +319,7 @@ CONTAINS
INTENT(inout), OPTIONAL :: needs
INTEGER, INTENT(out), OPTIONAL :: max_deriv
LOGICAL, INTENT(IN), OPTIONAL :: print_warn
CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: func_name_override
#if defined (__LIBXC)
CHARACTER(LEN=128) :: s1, s2
@ -326,8 +329,13 @@ CONTAINS
TYPE(xc_f03_func_t) :: xc_func
TYPE(xc_f03_func_info_t) :: xc_info
func_name = libxc_params%section%name
CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale)
IF (PRESENT(func_name_override)) THEN
func_name = func_name_override
func_scale = 1.0_dp
ELSE
func_name = libxc_params%section%name
CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale)
END IF
CALL cite_reference(Marques2012)
CALL cite_reference(Lehtola2018)
@ -398,6 +406,7 @@ CONTAINS
MARK_USED(needs)
MARK_USED(max_deriv)
MARK_USED(print_warn)
MARK_USED(func_name_override)
CALL cp_abort(__LOCATION__, "Unknown functional! If you are asking "// &
"for a functional of the LibXC library, "// &
@ -415,9 +424,11 @@ CONTAINS
!> true (does not set the unneeded components to false)
!> \param max_deriv maximum implemented derivative of the xc functional
!> \param print_warn whether to print warning about development status of a functional
!> \param func_name_override optional LibXC functional name overriding the section name
!> \author F. Tran
! **************************************************************************************************
SUBROUTINE libxc_lsd_info(libxc_params, reference, shortform, needs, max_deriv, print_warn)
SUBROUTINE libxc_lsd_info(libxc_params, reference, shortform, needs, max_deriv, print_warn, &
func_name_override)
TYPE(section_vals_type), POINTER :: libxc_params
CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
@ -425,6 +436,7 @@ CONTAINS
INTENT(inout), OPTIONAL :: needs
INTEGER, INTENT(out), OPTIONAL :: max_deriv
LOGICAL, INTENT(IN), OPTIONAL :: print_warn
CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: func_name_override
#if defined (__LIBXC)
CHARACTER(LEN=128) :: s1, s2
@ -434,8 +446,13 @@ CONTAINS
TYPE(xc_f03_func_t) :: xc_func
TYPE(xc_f03_func_info_t) :: xc_info
func_name = libxc_params%section%name
CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale)
IF (PRESENT(func_name_override)) THEN
func_name = func_name_override
func_scale = 1.0_dp
ELSE
func_name = libxc_params%section%name
CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale)
END IF
CALL cite_reference(Marques2012)
CALL cite_reference(Lehtola2018)
@ -508,6 +525,7 @@ CONTAINS
MARK_USED(needs)
MARK_USED(max_deriv)
MARK_USED(print_warn)
MARK_USED(func_name_override)
CALL cp_abort(__LOCATION__, "Unknown functional! If you are "// &
"asking for a functional of the LibXC library, "// &
@ -542,14 +560,16 @@ CONTAINS
!> if positive all the derivatives up to the given degree are evaluated,
!> if negative only the given degree is calculated
!> \param libxc_params input parameter (functional name, scaling and parameters)
!> \param func_name_override optional LibXC functional name overriding the section name
!> \author F. Tran
! **************************************************************************************************
SUBROUTINE libxc_lda_eval(rho_set, deriv_set, grad_deriv, libxc_params)
SUBROUTINE libxc_lda_eval(rho_set, deriv_set, grad_deriv, libxc_params, func_name_override)
TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
INTEGER, INTENT(in) :: grad_deriv
TYPE(section_vals_type), POINTER :: libxc_params
CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: func_name_override
#if defined (__LIBXC)
CHARACTER(len=*), PARAMETER :: routineN = 'libxc_lda_eval'
@ -574,15 +594,22 @@ CONTAINS
NULLIFY (dummy)
NULLIFY (rho, norm_drho, laplace_rho, tau)
func_name = libxc_params%section%name
CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale)
IF (PRESENT(func_name_override)) THEN
func_name = func_name_override
func_scale = 1.0_dp
ELSE
func_name = libxc_params%section%name
CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale)
END IF
IF (ABS(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp
func_id = xc_libxc_wrap_functional_get_number(func_name)
CALL xc_f03_func_init(xc_func, func_id, XC_UNPOLARIZED)
xc_info = xc_f03_func_get_info(xc_func)
CALL xc_libxc_wrap_functional_set_params(xc_func, xc_info, libxc_params, no_exc)
no_exc = .FALSE.
IF (.NOT. PRESENT(func_name_override)) &
CALL xc_libxc_wrap_functional_set_params(xc_func, xc_info, libxc_params, no_exc)
CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., &
rho=rho, norm_drho=norm_drho, laplace_rho=laplace_rho, &
@ -762,6 +789,7 @@ CONTAINS
MARK_USED(deriv_set)
MARK_USED(grad_deriv)
MARK_USED(libxc_params)
MARK_USED(func_name_override)
CALL cp_abort(__LOCATION__, "Unknown functional! If you are asking "// &
"for a functional of the LibXC library, "// &
"you have to download and install the library!")
@ -777,14 +805,16 @@ CONTAINS
!> if positive all the derivatives up to the given degree are evaluated,
!> if negative only the given degree is calculated
!> \param libxc_params input parameter (functional name, scaling and parameters)
!> \param func_name_override optional LibXC functional name overriding the section name
!> \author F. Tran
! **************************************************************************************************
SUBROUTINE libxc_lsd_eval(rho_set, deriv_set, grad_deriv, libxc_params)
SUBROUTINE libxc_lsd_eval(rho_set, deriv_set, grad_deriv, libxc_params, func_name_override)
TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
INTEGER, INTENT(in) :: grad_deriv
TYPE(section_vals_type), POINTER :: libxc_params
CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: func_name_override
#if defined (__LIBXC)
CHARACTER(len=*), PARAMETER :: routineN = 'libxc_lsd_eval'
@ -823,15 +853,22 @@ CONTAINS
NULLIFY (rhoa, rhob, norm_drho, norm_drhoa, norm_drhob, laplace_rhoa, &
laplace_rhob, tau_a, tau_b)
func_name = libxc_params%section%name
CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale)
IF (PRESENT(func_name_override)) THEN
func_name = func_name_override
func_scale = 1.0_dp
ELSE
func_name = libxc_params%section%name
CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale)
END IF
IF (ABS(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp
func_id = xc_libxc_wrap_functional_get_number(func_name)
CALL xc_f03_func_init(xc_func, func_id, XC_POLARIZED)
xc_info = xc_f03_func_get_info(xc_func)
CALL xc_libxc_wrap_functional_set_params(xc_func, xc_info, libxc_params, no_exc)
no_exc = .FALSE.
IF (.NOT. PRESENT(func_name_override)) &
CALL xc_libxc_wrap_functional_set_params(xc_func, xc_info, libxc_params, no_exc)
CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., &
rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, &
@ -1290,6 +1327,7 @@ CONTAINS
MARK_USED(deriv_set)
MARK_USED(grad_deriv)
MARK_USED(libxc_params)
MARK_USED(func_name_override)
CALL cp_abort(__LOCATION__, "Unknown functional! If you are asking "// &
"for a functional of the LibXC library, "// &

12
tests/QS/GAPW_ACCURATE_XCINT_COVERAGE.md
View file

@ -20,6 +20,8 @@ used to move towards making accurate XC integration the GAPW default in a future
| DC-DFT/Energy Correction | `regtest-acc-2/HF-ec1.inp`, `HF-ec2.inp`, `HF-ec3.inp` | `HF-ec4.inp`, `HF-ec5.inp`, `HF-ec6.inp`, and `HF-ec7.inp` now also use explicit accurate integration |
| TDDFPT forces | `regtest-acc-5/h2o_f01.inp` | `regtest-acc-5/h2o_f01_fine.inp` |
| ADMM-GAPW TDDFPT response | `regtest-acc-5/ft3.inp` | `regtest-acc-5/ft3_fine.inp` |
| UZH basis/potential data | all-electron UZH HF force-debug tests in `regtest-acc-1` | GAPW all-electron and GAPW/GAPW_XC GTH finite-difference force and stress checks with `BASIS_MOLOPT_UZH` and `POTENTIAL_UZH` |
| def2-ECP data | def2-SVP ECP GAPW energy checks in `regtest-ecp` and `regtest-ecp-2` | GAPW mixed all-electron/ECP and GAPW_XC ECP-only finite-difference force and stress checks |
| XAS/RT response | XAS_TDP and RTBSE coverage in `regtest-xastdp` and `regtest-rtbse` | open-shell GAPW XAS_TDP and GAPW RTBSE smoke tests in `regtest-acc-3` |
| KG embedding | existing `regtest-kg` GPW cases | GAPW/GAPW_XC energy and stress checks for libxc KG, meta-GGA/tau KG, and RI embedding |
| KG atomic potential | `regtest-kg/H2_KG-1.inp` | GAPW/GAPW_XC energy and stress checks for `TNADD_METHOD ATOMIC` |
@ -55,6 +57,16 @@ Finite-difference checks with `STOP_ON_MISMATCH` are included for:
is kept diagonal in the targeted debug tests above.
- GAPW_XC forces and diagonal stress: `regtest-acc-1/h2o-gapw_xc-force-1.inp`,
`regtest-acc-1/h2o-gapw_xc-stress-debug-1.inp`.
- UZH all-electron GAPW forces and analytical stress smoke:
`regtest-acc-1/h2o-uzh-gapw-force-1.inp`, `regtest-acc-1/h2o-uzh-gapw-stress-debug-1.inp`.
- UZH GTH GAPW/GAPW_XC forces and analytical stress smoke:
`regtest-acc-1/h2o-uzh-gth-gapw-force-1.inp`, `regtest-acc-1/h2o-uzh-gth-gapw-stress-debug-1.inp`,
`regtest-acc-1/h2o-uzh-gth-gapw_xc-force-1.inp`, and
`regtest-acc-1/h2o-uzh-gth-gapw_xc-stress-debug-1.inp`.
- ECP GAPW/GAPW_XC targeted force and analytical stress smoke checks:
`regtest-ecp/ICl_lanl2dz_gapw_force.inp`, `ICl_lanl2dz_gapw_stress.inp`,
`ICl_lanl2dz_gapw_xc_force.inp`, and `ICl_lanl2dz_gapw_xc_stress.inp`; the def2-ECP energy path
remains covered by `regtest-ecp/SbH3_def2_gapw.inp`.
- ADMM-GAPW forces: `regtest-acc-1/HF-d5.inp`, `regtest-acc-5/ft3_fine.inp`.
- ADMM-GAPW diagonal stress: `regtest-acc-2/h2o-admm-gapw-stress-debug-1.inp`,
`regtest-acc-2/h2o-admm-gapw-pbe-stress-debug-1.inp`.

6
tests/QS/regtest-acc-1/TEST_FILES.toml
View file

@ -36,4 +36,10 @@
"h2o-gapw_xc-stress-debug-1.inp" = []
"h2o-gapw_xc-full-stress-1.inp" = [{matcher="M011", tol=1e-08, ref=-17.256333793570239},
{matcher="M031", tol=1e-08, ref=3.91898639036E+03}]
"h2o-uzh-gapw-force-1.inp" = []
"h2o-uzh-gapw-stress-debug-1.inp" = []
"h2o-uzh-gth-gapw-force-1.inp" = []
"h2o-uzh-gth-gapw-stress-debug-1.inp" = []
"h2o-uzh-gth-gapw_xc-force-1.inp" = []
"h2o-uzh-gth-gapw_xc-stress-debug-1.inp" = []
#EOF

63
tests/QS/regtest-acc-1/h2o-uzh-gapw-force-1.inp Normal file
View file

@ -0,0 +1,63 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT h2o-uzh-gapw-force-1
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
CHECK_ATOM_FORCE 1 x
DEBUG_FORCES T
DEBUG_STRESS_TENSOR F
DX 0.0005
MAX_RELATIVE_ERROR 0.10
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS_MOLOPT_UZH
POTENTIAL_FILE_NAME POTENTIAL_UZH
&MGRID
CUTOFF 240
NGRIDS 5
REL_CUTOFF 40
&END MGRID
&QS
ALPHA_WEIGHTS 6.5
EPS_DEFAULT 1.0E-10
GAPW_ACCURATE_XCINT T
METHOD GAPW
&END QS
&SCF
EPS_SCF 1.0E-7
MAX_SCF 80
SCF_GUESS ATOMIC
&END SCF
&XC
DENSITY_CUTOFF 1.0E-11
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
A 5.0 0.0 0.0
B 0.25 5.1 0.0
C 0.15 0.35 5.2
&END CELL
&COORD
O 0.110000 0.170000 -0.045587
H 0.260000 -0.627136 0.570545
H -0.180000 0.897136 0.440545
&END COORD
&KIND H
BASIS_SET SVP-MOLOPT-PBE-ae
POTENTIAL ALL
&END KIND
&KIND O
BASIS_SET SVP-MOLOPT-PBE-ae
POTENTIAL ALL
&END KIND
&END SUBSYS
&END FORCE_EVAL

60
tests/QS/regtest-acc-1/h2o-uzh-gapw-stress-debug-1.inp Normal file
View file

@ -0,0 +1,60 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT h2o-uzh-gapw-stress-debug-1
RUN_TYPE ENERGY_FORCE
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
STRESS_TENSOR ANALYTICAL
&DFT
BASIS_SET_FILE_NAME BASIS_MOLOPT_UZH
POTENTIAL_FILE_NAME POTENTIAL_UZH
&MGRID
CUTOFF 240
NGRIDS 5
REL_CUTOFF 40
&END MGRID
&QS
ALPHA_WEIGHTS 6.5
EPS_DEFAULT 1.0E-10
GAPW_ACCURATE_XCINT T
METHOD GAPW
&END QS
&SCF
EPS_SCF 1.0E-7
MAX_SCF 80
SCF_GUESS ATOMIC
&END SCF
&XC
DENSITY_CUTOFF 1.0E-11
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PRINT
&STRESS_TENSOR
COMPONENTS
&END STRESS_TENSOR
&END PRINT
&SUBSYS
&CELL
A 5.0 0.0 0.0
B 0.25 5.1 0.0
C 0.15 0.35 5.2
&END CELL
&COORD
O 0.110000 0.170000 -0.045587
H 0.260000 -0.627136 0.570545
H -0.180000 0.897136 0.440545
&END COORD
&KIND H
BASIS_SET SVP-MOLOPT-PBE-ae
POTENTIAL ALL
&END KIND
&KIND O
BASIS_SET SVP-MOLOPT-PBE-ae
POTENTIAL ALL
&END KIND
&END SUBSYS
&END FORCE_EVAL

67
tests/QS/regtest-acc-1/h2o-uzh-gth-gapw-force-1.inp Normal file
View file

@ -0,0 +1,67 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT h2o-uzh-gth-gapw-force-1
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
CHECK_ATOM_FORCE 1 x
DEBUG_FORCES T
DEBUG_STRESS_TENSOR F
DX 0.0005
MAX_RELATIVE_ERROR 0.10
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS_MOLOPT_UZH
POTENTIAL_FILE_NAME POTENTIAL_UZH
&MGRID
CUTOFF 240
NGRIDS 5
REL_CUTOFF 40
&END MGRID
&QS
ALPHA_WEIGHTS 6.5
EPS_DEFAULT 1.0E-10
GAPW_ACCURATE_XCINT T
METHOD GAPW
&END QS
&SCF
EPS_SCF 1.0E-7
MAX_SCF 80
SCF_GUESS ATOMIC
&END SCF
&XC
DENSITY_CUTOFF 1.0E-11
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
A 5.0 0.0 0.0
B 0.25 5.1 0.0
C 0.15 0.35 5.2
&END CELL
&COORD
O 0.110000 0.170000 -0.045587
H 0.260000 -0.627136 0.570545
H -0.180000 0.897136 0.440545
&END COORD
&KIND H
BASIS_SET DZVP-MOLOPT-PBE-GTH-q1
LEBEDEV_GRID 50
POTENTIAL GTH-PBE-q1
RADIAL_GRID 50
&END KIND
&KIND O
BASIS_SET DZVP-MOLOPT-PBE-GTH-q6
LEBEDEV_GRID 50
POTENTIAL GTH-PBE-q6
RADIAL_GRID 50
&END KIND
&END SUBSYS
&END FORCE_EVAL

64
tests/QS/regtest-acc-1/h2o-uzh-gth-gapw-stress-debug-1.inp Normal file
View file

@ -0,0 +1,64 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT h2o-uzh-gth-gapw-stress-debug-1
RUN_TYPE ENERGY_FORCE
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
STRESS_TENSOR ANALYTICAL
&DFT
BASIS_SET_FILE_NAME BASIS_MOLOPT_UZH
POTENTIAL_FILE_NAME POTENTIAL_UZH
&MGRID
CUTOFF 240
NGRIDS 5
REL_CUTOFF 40
&END MGRID
&QS
ALPHA_WEIGHTS 6.5
EPS_DEFAULT 1.0E-10
GAPW_ACCURATE_XCINT T
METHOD GAPW
&END QS
&SCF
EPS_SCF 1.0E-7
MAX_SCF 80
SCF_GUESS ATOMIC
&END SCF
&XC
DENSITY_CUTOFF 1.0E-11
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PRINT
&STRESS_TENSOR
COMPONENTS
&END STRESS_TENSOR
&END PRINT
&SUBSYS
&CELL
A 5.0 0.0 0.0
B 0.25 5.1 0.0
C 0.15 0.35 5.2
&END CELL
&COORD
O 0.110000 0.170000 -0.045587
H 0.260000 -0.627136 0.570545
H -0.180000 0.897136 0.440545
&END COORD
&KIND H
BASIS_SET DZVP-MOLOPT-PBE-GTH-q1
LEBEDEV_GRID 50
POTENTIAL GTH-PBE-q1
RADIAL_GRID 50
&END KIND
&KIND O
BASIS_SET DZVP-MOLOPT-PBE-GTH-q6
LEBEDEV_GRID 50
POTENTIAL GTH-PBE-q6
RADIAL_GRID 50
&END KIND
&END SUBSYS
&END FORCE_EVAL

69
tests/QS/regtest-acc-1/h2o-uzh-gth-gapw_xc-force-1.inp Normal file
View file

@ -0,0 +1,69 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT h2o-uzh-gth-gapw_xc-force-1
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
CHECK_ATOM_FORCE 1 x
DEBUG_FORCES T
DEBUG_STRESS_TENSOR F
DX 0.0005
MAX_RELATIVE_ERROR 0.10
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS_MOLOPT_UZH
POTENTIAL_FILE_NAME POTENTIAL_UZH
&MGRID
CUTOFF 240
NGRIDS 5
REL_CUTOFF 40
&END MGRID
&QS
ALPHA_WEIGHTS 6.5
EPS_DEFAULT 1.0E-10
FORCE_PAW
GAPW_1C_BASIS EXT_SMALL
GAPW_ACCURATE_XCINT T
METHOD GAPW_XC
&END QS
&SCF
EPS_SCF 1.0E-7
MAX_SCF 80
SCF_GUESS ATOMIC
&END SCF
&XC
DENSITY_CUTOFF 1.0E-11
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
A 5.0 0.0 0.0
B 0.25 5.1 0.0
C 0.15 0.35 5.2
&END CELL
&COORD
O 0.110000 0.170000 -0.045587
H 0.260000 -0.627136 0.570545
H -0.180000 0.897136 0.440545
&END COORD
&KIND H
BASIS_SET DZVP-MOLOPT-PBE-GTH-q1
LEBEDEV_GRID 50
POTENTIAL GTH-PBE-q1
RADIAL_GRID 50
&END KIND
&KIND O
BASIS_SET DZVP-MOLOPT-PBE-GTH-q6
LEBEDEV_GRID 50
POTENTIAL GTH-PBE-q6
RADIAL_GRID 50
&END KIND
&END SUBSYS
&END FORCE_EVAL

66
tests/QS/regtest-acc-1/h2o-uzh-gth-gapw_xc-stress-debug-1.inp Normal file
View file

@ -0,0 +1,66 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT h2o-uzh-gth-gapw_xc-stress-debug-1
RUN_TYPE ENERGY_FORCE
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
STRESS_TENSOR ANALYTICAL
&DFT
BASIS_SET_FILE_NAME BASIS_MOLOPT_UZH
POTENTIAL_FILE_NAME POTENTIAL_UZH
&MGRID
CUTOFF 240
NGRIDS 5
REL_CUTOFF 40
&END MGRID
&QS
ALPHA_WEIGHTS 6.5
EPS_DEFAULT 1.0E-10
FORCE_PAW
GAPW_1C_BASIS EXT_SMALL
GAPW_ACCURATE_XCINT T
METHOD GAPW_XC
&END QS
&SCF
EPS_SCF 1.0E-7
MAX_SCF 80
SCF_GUESS ATOMIC
&END SCF
&XC
DENSITY_CUTOFF 1.0E-11
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PRINT
&STRESS_TENSOR
COMPONENTS
&END STRESS_TENSOR
&END PRINT
&SUBSYS
&CELL
A 5.0 0.0 0.0
B 0.25 5.1 0.0
C 0.15 0.35 5.2
&END CELL
&COORD
O 0.110000 0.170000 -0.045587
H 0.260000 -0.627136 0.570545
H -0.180000 0.897136 0.440545
&END COORD
&KIND H
BASIS_SET DZVP-MOLOPT-PBE-GTH-q1
LEBEDEV_GRID 50
POTENTIAL GTH-PBE-q1
RADIAL_GRID 50
&END KIND
&KIND O
BASIS_SET DZVP-MOLOPT-PBE-GTH-q6
LEBEDEV_GRID 50
POTENTIAL GTH-PBE-q6
RADIAL_GRID 50
&END KIND
&END SUBSYS
&END FORCE_EVAL

62
tests/QS/regtest-ecp/ICl_lanl2dz_gapw_force.inp Normal file
View file

@ -0,0 +1,62 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT ICl_lanl2dz_gapw_force
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
CHECK_ATOM_FORCE 1 x
DEBUG_FORCES T
DEBUG_STRESS_TENSOR F
DX 0.0005
MAX_RELATIVE_ERROR 0.10
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME ./ECP_BASIS_POT
POTENTIAL_FILE_NAME ./ECP_BASIS_POT
&MGRID
CUTOFF 160
NGRIDS 4
REL_CUTOFF 30
&END MGRID
&QS
ALPHA_WEIGHTS 6.5
EPS_DEFAULT 1.0E-9
GAPW_ACCURATE_XCINT T
METHOD GAPW
&END QS
&SCF
EPS_SCF 1.0E-7
MAX_SCF 80
SCF_GUESS ATOMIC
&END SCF
&XC
DENSITY_CUTOFF 1.0E-11
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
A 5.0 0.0 0.0
B 0.25 5.1 0.0
C 0.15 0.35 5.2
&END CELL
&COORD
Cl 0.000000 0.000000 0.000000
I 0.250000 0.100000 2.490000
&END COORD
&KIND Cl
BASIS_SET LANL2DZ
POTENTIAL ECP LANL2DZ_ECP
&END KIND
&KIND I
BASIS_SET LANL2DZ
POTENTIAL ECP LANL2DZ_ECP
&END KIND
&END SUBSYS
&END FORCE_EVAL

59
tests/QS/regtest-ecp/ICl_lanl2dz_gapw_stress.inp Normal file
View file

@ -0,0 +1,59 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT ICl_lanl2dz_gapw_stress
RUN_TYPE ENERGY_FORCE
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
STRESS_TENSOR ANALYTICAL
&DFT
BASIS_SET_FILE_NAME ./ECP_BASIS_POT
POTENTIAL_FILE_NAME ./ECP_BASIS_POT
&MGRID
CUTOFF 160
NGRIDS 4
REL_CUTOFF 30
&END MGRID
&QS
ALPHA_WEIGHTS 6.5
EPS_DEFAULT 1.0E-9
GAPW_ACCURATE_XCINT T
METHOD GAPW
&END QS
&SCF
EPS_SCF 1.0E-7
MAX_SCF 80
SCF_GUESS ATOMIC
&END SCF
&XC
DENSITY_CUTOFF 1.0E-11
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PRINT
&STRESS_TENSOR
COMPONENTS
&END STRESS_TENSOR
&END PRINT
&SUBSYS
&CELL
A 5.0 0.0 0.0
B 0.25 5.1 0.0
C 0.15 0.35 5.2
&END CELL
&COORD
Cl 0.000000 0.000000 0.000000
I 0.250000 0.100000 2.490000
&END COORD
&KIND Cl
BASIS_SET LANL2DZ
POTENTIAL ECP LANL2DZ_ECP
&END KIND
&KIND I
BASIS_SET LANL2DZ
POTENTIAL ECP LANL2DZ_ECP
&END KIND
&END SUBSYS
&END FORCE_EVAL

64
tests/QS/regtest-ecp/ICl_lanl2dz_gapw_xc_force.inp Normal file
View file

@ -0,0 +1,64 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT ICl_lanl2dz_gapw_xc_force
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
CHECK_ATOM_FORCE 1 x
DEBUG_FORCES T
DEBUG_STRESS_TENSOR F
DX 0.0005
MAX_RELATIVE_ERROR 0.10
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME ./ECP_BASIS_POT
POTENTIAL_FILE_NAME ./ECP_BASIS_POT
&MGRID
CUTOFF 160
NGRIDS 4
REL_CUTOFF 30
&END MGRID
&QS
ALPHA_WEIGHTS 6.5
EPS_DEFAULT 1.0E-9
FORCE_PAW
GAPW_1C_BASIS EXT_SMALL
GAPW_ACCURATE_XCINT T
METHOD GAPW_XC
&END QS
&SCF
EPS_SCF 1.0E-7
MAX_SCF 80
SCF_GUESS ATOMIC
&END SCF
&XC
DENSITY_CUTOFF 1.0E-11
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
A 6.0 0.0 0.0
B 0.25 6.1 0.0
C 0.15 0.35 6.2
&END CELL
&COORD
Cl 0.000000 0.000000 0.000000
I 0.250000 0.100000 2.490000
&END COORD
&KIND Cl
BASIS_SET LANL2DZ
POTENTIAL ECP LANL2DZ_ECP
&END KIND
&KIND I
BASIS_SET LANL2DZ
POTENTIAL ECP LANL2DZ_ECP
&END KIND
&END SUBSYS
&END FORCE_EVAL

61
tests/QS/regtest-ecp/ICl_lanl2dz_gapw_xc_stress.inp Normal file
View file

@ -0,0 +1,61 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT ICl_lanl2dz_gapw_xc_stress
RUN_TYPE ENERGY_FORCE
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
STRESS_TENSOR ANALYTICAL
&DFT
BASIS_SET_FILE_NAME ./ECP_BASIS_POT
POTENTIAL_FILE_NAME ./ECP_BASIS_POT
&MGRID
CUTOFF 160
NGRIDS 4
REL_CUTOFF 30
&END MGRID
&QS
ALPHA_WEIGHTS 6.5
EPS_DEFAULT 1.0E-9
FORCE_PAW
GAPW_1C_BASIS EXT_SMALL
GAPW_ACCURATE_XCINT T
METHOD GAPW_XC
&END QS
&SCF
EPS_SCF 1.0E-7
MAX_SCF 80
SCF_GUESS ATOMIC
&END SCF
&XC
DENSITY_CUTOFF 1.0E-11
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PRINT
&STRESS_TENSOR
COMPONENTS
&END STRESS_TENSOR
&END PRINT
&SUBSYS
&CELL
A 6.0 0.0 0.0
B 0.25 6.1 0.0
C 0.15 0.35 6.2
&END CELL
&COORD
Cl 0.000000 0.000000 0.000000
I 0.250000 0.100000 2.490000
&END COORD
&KIND Cl
BASIS_SET LANL2DZ
POTENTIAL ECP LANL2DZ_ECP
&END KIND
&KIND I
BASIS_SET LANL2DZ
POTENTIAL ECP LANL2DZ_ECP
&END KIND
&END SUBSYS
&END FORCE_EVAL

4
tests/QS/regtest-ecp/TEST_FILES.toml
View file

@ -3,3 +3,7 @@
"SbH3_def2_gapw.inp" = [{matcher="M011", tol=1.0E-11, ref=-241.414769303399083}]
"HCl_ccECP.inp" = [{matcher="M011", tol=1.0E-12, ref=-15.457709431781252}]
"HCl_force.inp" = []
"ICl_lanl2dz_gapw_force.inp" = []
"ICl_lanl2dz_gapw_stress.inp" = []
"ICl_lanl2dz_gapw_xc_force.inp" = []
"ICl_lanl2dz_gapw_xc_stress.inp" = []

4
tests/QS/regtest-gauxc-api/TEST_FILES.toml
View file

@ -1,7 +1,7 @@
"1H2_GAUXC_MODEL_PBE_REFERENCE.inp" = [{matcher="E_total", tol=1e-10, ref=-1.163125608332391}]
"1H2_GAUXC_PBE.inp" = [{matcher="E_total", tol=1e-10, ref=-1.163121772046629}]
"1H2_GAUXC_MODEL_PBE.inp" = [{matcher="E_total", tol=1e-10, ref=-1.163121899608445}]
"1H2_GAUXC_MODEL_PBE.inp" = [{matcher="E_total", tol=1e-10, ref=-1.163121772046629}]
"NH3_GAUXC_MODEL_PBE_REFERENCE.inp" = [{matcher="E_total", tol=1e-9, ref=-11.722432805445091}]
"NH3_GAUXC_MODEL_PBE.inp" = [{matcher="E_total", tol=1e-9, ref=-11.722558568119791}]
"NH3_GAUXC_MODEL_PBE.inp" = [{matcher="E_total", tol=1e-9, ref=-11.722557896081605}]
"OH_GAUXC_MODEL_PBE_UKS.inp" = [{matcher="E_total", tol=5e-6, ref=-16.541584062034670}]
"CH4_DIMER_GAUXC_PBE_D3.inp" = [{matcher="M033", tol=1e-14, ref=-0.00355123783846}]

6
tests/QS/regtest-gauxc-cdft/TEST_FILES.toml
View file

@ -1,11 +1,11 @@
"H2_PBE_CDFT_REFERENCE.inp" = [{matcher="E_total", tol=1e-9, ref=-1.157232743213941},
{matcher="M071", tol=1e-8, ref=0.197625731046}]
"H2_GAUXC_MODEL_PBE_CDFT.inp" = [{matcher="E_total", tol=1e-9, ref=-1.157229554329094},
"H2_GAUXC_MODEL_PBE_CDFT.inp" = [{matcher="E_total", tol=1e-9, ref=-1.157229340943540},
{matcher="M071", tol=1e-8, ref=0.197624070585}]
"H2_SKALA_CDFT.inp" = [{matcher="M071", tol=1e-3, ref=0.216557647288}]
"H2_GAPW_SKALA_CDFT.inp" = [{matcher="M071", tol=1e-3, ref=0.217417757774}]
"H2_PBE_CDFT_CI_REFERENCE.inp" = [{matcher="M073", tol=1e-8, ref=345.364329058819},
{matcher="M077", tol=1e-8, ref=-1.16295181842678}]
"H2_GAUXC_MODEL_PBE_CDFT_CI_NGROUPS.inp" = [{matcher="M073", tol=1e-8, ref=345.325437357859},
{matcher="M077", tol=1e-8, ref=-1.16294823026735}]
"H2_GAUXC_MODEL_PBE_CDFT_CI_NGROUPS.inp" = [{matcher="M073", tol=1e-8, ref=345.325380057211},
{matcher="M077", tol=1e-8, ref=-1.16294801624740}]
"H2_SKALA_CDFT_CI.inp" = [{matcher="M077", tol=5e-8, ref=-1.39020335005227}]

10
tests/QS/regtest-gauxc-gapw-ecp/TEST_FILES.toml
View file

@ -1,8 +1,8 @@
"HCl_GAPW_SKALA_ECP_ENERGY.inp" = [{matcher="E_total", tol=5e-6, ref=-15.464581508021762}]
"HCl_NATIVE_SKALA_GAPW_ECP_STRESS.inp" = [{matcher="E_total", tol=5e-6, ref=-15.45263701}]
"HCl_NATIVE_SKALA_GAPW_ECP_STRESS.inp" = [{matcher="E_total", tol=5e-6, ref=-15.44142596263930}]
"HCl_NATIVE_SKALA_GAPW_GPWTYPE_ECP_STRESS.inp" = [{matcher="E_total", tol=5e-6, ref=-15.36078871}]
"HCl_NATIVE_SKALA_GAPW_XC_ECP_STRESS.inp" = [{matcher="E_total", tol=5e-6, ref=-15.36291913}]
"HCl_NATIVE_SKALA_GAPW_XC_ECP_STRESS.inp" = [{matcher="E_total", tol=5e-6, ref=-15.35170807520943}]
"H2_NATIVE_SKALA_GAPW_ECP_STRESS_DEBUG.inp" = [{matcher="DEBUG_stress_sum", tol=6e-1, ref=0.0}]
"HCl_NATIVE_SKALA_GAPW_ECP_KP_INV_STRESS.inp" = [{matcher="E_total", tol=5e-6, ref=-15.44577629}]
"HCl_NATIVE_SKALA_GAPW_ECP_KP_SYM_STRESS.inp" = [{matcher="E_total", tol=5e-6, ref=-15.44577629}]
"HCl_NATIVE_SKALA_GAPW_ECP_KP_SPGLIB_STRESS.inp" = [{matcher="E_total", tol=5e-6, ref=-15.44577629}]
"HCl_NATIVE_SKALA_GAPW_ECP_KP_INV_STRESS.inp" = [{matcher="E_total", tol=5e-6, ref=-15.43466946670853}]
"HCl_NATIVE_SKALA_GAPW_ECP_KP_SYM_STRESS.inp" = [{matcher="E_total", tol=5e-6, ref=-15.43466946670853}]
"HCl_NATIVE_SKALA_GAPW_ECP_KP_SPGLIB_STRESS.inp" = [{matcher="E_total", tol=5e-6, ref=-15.43466946670853}]

6
tests/QS/regtest-gauxc-gapw-gth-kp/TEST_FILES.toml
View file

@ -1,3 +1,3 @@
"H2_NATIVE_SKALA_GAPW_GTH_KP_INV_STRESS.inp" = [{matcher="E_total", tol=1e-8, ref=-0.97063781133388}]
"H2_NATIVE_SKALA_GAPW_GTH_KP_SYM_STRESS.inp" = [{matcher="E_total", tol=1e-8, ref=-0.97063781133388}]
"H2_NATIVE_SKALA_GAPW_GTH_KP_SPGLIB_STRESS.inp" = [{matcher="E_total", tol=1e-8, ref=-0.97063781133388}]
"H2_NATIVE_SKALA_GAPW_GTH_KP_INV_STRESS.inp" = [{matcher="E_total", tol=1e-8, ref=-0.94951222851689}]
"H2_NATIVE_SKALA_GAPW_GTH_KP_SYM_STRESS.inp" = [{matcher="E_total", tol=1e-8, ref=-0.94951222851689}]
"H2_NATIVE_SKALA_GAPW_GTH_KP_SPGLIB_STRESS.inp" = [{matcher="E_total", tol=1e-8, ref=-0.94951222851689}]

6
tests/QS/regtest-gauxc-gapw-gth/TEST_FILES.toml
View file

@ -1,8 +1,8 @@
"H2_NATIVE_SKALA_GAPW_GTH_STRESS.inp" = [{matcher="E_total", tol=5e-8, ref=-0.9912833172}]
"H2_NATIVE_SKALA_GAPW_GTH_STRESS_CHUNK_REQUEST.inp" = [{matcher="E_total", tol=5e-8, ref=-0.9912833169}]
"H2_NATIVE_SKALA_GAPW_GTH_STRESS.inp" = [{matcher="E_total", tol=5e-8, ref=-0.97019845356640}]
"H2_NATIVE_SKALA_GAPW_GTH_STRESS_CHUNK_REQUEST.inp" = [{matcher="E_total", tol=5e-8, ref=-0.97019845356640}]
"H2_NATIVE_SKALA_GAPW_GPWTYPE_STRESS.inp" = [{matcher="E_total", tol=1e-8, ref=-0.96567258919119},
{matcher="M072", tol=1e-4, ref=5.37252718E-05},
{matcher="M031", tol=1e-5, ref=-8.58976831594E+04}]
"H2_NATIVE_SKALA_GAPW_XC_GTH_STRESS.inp" = [{matcher="E_total", tol=5e-8, ref=-0.9698629747}]
"H2_NATIVE_SKALA_GAPW_XC_GTH_STRESS.inp" = [{matcher="E_total", tol=5e-8, ref=-0.94877811114650}]
"H2_NATIVE_SKALA_GAPW_GTH_FORCE_DEBUG.inp" = [{matcher="DEBUG_force_sum", tol=2e-1, ref=0.0}]
"H2_NATIVE_SKALA_GAPW_GTH_STRESS_DEBUG.inp" = [{matcher="DEBUG_stress_sum", tol=2e-1, ref=0.0}]

606
tools/toolchain/scripts/stage6/gauxc-1.1-skala-cp2k-fixes.patch
View file

@ -1,8 +1,11 @@
diff --git a/cmake/gauxc-config.cmake.in b/cmake/gauxc-config.cmake.in
index c6c4c2d4..35e426af 100644
index c6c4c2d4..b517ccc5 100644
--- a/cmake/gauxc-config.cmake.in
+++ b/cmake/gauxc-config.cmake.in
@@ -8,0 +9,8 @@ include(CMakeFindDependencyMacro)
@@ -6,6 +6,14 @@ list(PREPEND CMAKE_MODULE_PATH ${GauXC_CMAKE_DIR} )
list(PREPEND CMAKE_MODULE_PATH ${GauXC_CMAKE_DIR}/linalg-cmake-modules )
include(CMakeFindDependencyMacro)
+if(POLICY CMP0144)
+ cmake_policy(PUSH)
+ cmake_policy(SET CMP0144 NEW)
@ -11,7 +14,13 @@ index c6c4c2d4..35e426af 100644
+ endif()
+ set(CMAKE_POLICY_DEFAULT_CMP0144 NEW)
+endif()
@@ -93,0 +102,9 @@ endif()
# Always Required Dependencies
find_dependency( ExchCXX )
find_dependency( IntegratorXX )
@@ -91,4 +99,13 @@ if(NOT TARGET gauxc::gauxc)
include("${GauXC_CMAKE_DIR}/gauxc-targets.cmake")
endif()
+if(POLICY CMP0144)
+ if(DEFINED _GAUXC_PREV_CMAKE_POLICY_DEFAULT_CMP0144)
+ set(CMAKE_POLICY_DEFAULT_CMP0144 "${_GAUXC_PREV_CMAKE_POLICY_DEFAULT_CMP0144}")
@ -21,33 +30,7 @@ index c6c4c2d4..35e426af 100644
+ unset(_GAUXC_PREV_CMAKE_POLICY_DEFAULT_CMP0144)
+ cmake_policy(POP)
+endif()
diff --git a/external/gau2grid/generated_source/gau2grid/gau2grid_pragma.h b/external/gau2grid/generated_source/gau2grid/gau2grid_pragma.h
--- a/external/gau2grid/generated_source/gau2grid/gau2grid_pragma.h
+++ b/external/gau2grid/generated_source/gau2grid/gau2grid_pragma.h
@@ -14,0 +15,8 @@
+#if (defined(__GNUC__) || defined(__GNUG__)) && \
+ !defined(__APPLE__) && !defined(__clang__)
+static inline void* gg_aligned_alloc(const size_t alignment, const size_t size) {
+ const size_t aligned_size =
+ ((size + alignment - 1) / alignment) * alignment;
+ return aligned_alloc(alignment, aligned_size);
+}
+#endif
@@ -15,0 +24 @@
+
@@ -70 +79 @@
- #define ALIGNED_MALLOC(alignment, size) aligned_alloc(alignment, size)
+ #define ALIGNED_MALLOC(alignment, size) gg_aligned_alloc(alignment, size)
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index a9e60fd8..29b0949d 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -122,1 +122,4 @@ if ( GAUXC_HAS_ONEDFT )
- target_link_libraries( gauxc PUBLIC "${TORCH_LIBRARIES}" nlohmann_json::nlohmann_json)
+ get_target_property(GAUXC_NLOHMANN_JSON_INCLUDE_DIRS
+ nlohmann_json::nlohmann_json INTERFACE_INCLUDE_DIRECTORIES)
+ target_include_directories(gauxc PRIVATE ${GAUXC_NLOHMANN_JSON_INCLUDE_DIRS})
+ target_link_libraries( gauxc PUBLIC "${TORCH_LIBRARIES}")
set(GauXC_LIBRARIES gauxc::gauxc)
diff --git a/cmake/gauxc-exchcxx.cmake b/cmake/gauxc-exchcxx.cmake
index 412df9b3..011ae844 100644
--- a/cmake/gauxc-exchcxx.cmake
@ -98,6 +81,129 @@ index b6bbbf0e..502067d2 100644
+ endif()
+endif()
diff --git a/external/gau2grid/generated_source/gau2grid/gau2grid_pragma.h b/external/gau2grid/generated_source/gau2grid/gau2grid_pragma.h
index f6033886..4bbbf13f 100644
--- a/external/gau2grid/generated_source/gau2grid/gau2grid_pragma.h
+++ b/external/gau2grid/generated_source/gau2grid/gau2grid_pragma.h
@@ -14,0 +15,9 @@
+#if (defined(__GNUC__) || defined(__GNUG__)) && \
+ !defined(__APPLE__) && !defined(__clang__)
+static inline void* gg_aligned_alloc(const size_t alignment, const size_t size) {
+ const size_t aligned_size =
+ ((size + alignment - 1) / alignment) * alignment;
+ return aligned_alloc(alignment, aligned_size);
+}
+#endif
+
@@ -67,7 +76,7 @@
#elif defined(__GNUC__) || defined(__GNUG__)
// pragmas for GCC
- #define ALIGNED_MALLOC(alignment, size) aligned_alloc(alignment, size)
+ #define ALIGNED_MALLOC(alignment, size) gg_aligned_alloc(alignment, size)
#define ALIGNED_FREE(ptr) free(ptr)
#define ASSUME_ALIGNED(ptr, width)
diff --git a/include/gauxc/xc_integrator_settings.hpp b/include/gauxc/xc_integrator_settings.hpp
index a63899e6..9031be78 100644
--- a/include/gauxc/xc_integrator_settings.hpp
+++ b/include/gauxc/xc_integrator_settings.hpp
@@ -23,6 +23,9 @@ struct IntegratorSettingsSNLinK : public IntegratorSettingsEXX {
struct IntegratorSettingsXC { virtual ~IntegratorSettingsXC() noexcept = default; };
struct IntegratorSettingsKS : public IntegratorSettingsXC {
double gks_dtol = 1e-12;
+ // RKS density matrices are interpreted as one-spin densities by default.
+ // Set this when the caller provides the spin-summed closed-shell density.
+ bool rks_density_matrix_is_spin_summed = false;
};
struct OneDFTSettings : public IntegratorSettingsXC {
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index a9e60fd8..29b0949d 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -119,7 +119,10 @@ if( GAUXC_HAS_MPI )
endif()
if ( GAUXC_HAS_ONEDFT )
- target_link_libraries( gauxc PUBLIC "${TORCH_LIBRARIES}" nlohmann_json::nlohmann_json)
+ get_target_property(GAUXC_NLOHMANN_JSON_INCLUDE_DIRS
+ nlohmann_json::nlohmann_json INTERFACE_INCLUDE_DIRECTORIES)
+ target_include_directories(gauxc PRIVATE ${GAUXC_NLOHMANN_JSON_INCLUDE_DIRS})
+ target_link_libraries( gauxc PUBLIC "${TORCH_LIBRARIES}")
endif()
add_subdirectory( runtime_environment )
diff --git a/src/c-api/c_xc_integrator.cxx b/src/c-api/c_xc_integrator.cxx
index 534eb3d0..f1efc430 100644
--- a/src/c-api/c_xc_integrator.cxx
+++ b/src/c-api/c_xc_integrator.cxx
@@ -200,12 +200,14 @@ void gauxc_integrator_eval_exc_rks(
detail::gauxc_status_handle(status, 1, "Exc output pointer cannot be null");
return;
}
+ IntegratorSettingsKS ks_settings{};
+ ks_settings.rks_density_matrix_is_spin_summed = true;
try {
detail::get_xc_integrator_ptr(integrator)->eval_exc(
m, n,
density_matrix, ldp,
exc,
- IntegratorSettingsXC{} );
+ ks_settings );
} catch (std::exception& e) {
detail::gauxc_status_handle(status, 1, e.what());
}
@@ -333,13 +335,15 @@ void gauxc_integrator_eval_exc_vxc_rks(
detail::gauxc_status_handle(status, 1, "VXC matrix pointer cannot be null");
return;
}
+ IntegratorSettingsKS ks_settings{};
+ ks_settings.rks_density_matrix_is_spin_summed = true;
try {
detail::get_xc_integrator_ptr(integrator)->eval_exc_vxc(
m, n,
density_matrix, ldp,
vxc_matrix, vxc_ld,
exc,
- IntegratorSettingsXC{} );
+ ks_settings );
} catch (std::exception& e) {
detail::gauxc_status_handle(status, 1, e.what());
}
@@ -565,12 +569,14 @@ void gauxc_integrator_eval_exc_grad_rks(
detail::gauxc_status_handle(status, 1, "Exc gradient output pointer cannot be null");
return;
}
+ IntegratorSettingsEXC_GRAD exc_grad_settings{};
+ exc_grad_settings.rks_density_matrix_is_spin_summed = true;
try {
detail::get_xc_integrator_ptr(integrator)->eval_exc_grad(
m, n,
density_matrix, ldp,
exc_grad,
- IntegratorSettingsXC{} );
+ exc_grad_settings );
} catch (std::exception& e) {
detail::gauxc_status_handle(status, 1, e.what());
}
@@ -726,13 +732,15 @@ void gauxc_integrator_eval_fxc_contraction_rks(
detail::gauxc_status_handle(status, 1, "FXC output pointer cannot be null");
return;
}
+ IntegratorSettingsKS ks_settings{};
+ ks_settings.rks_density_matrix_is_spin_summed = true;
try {
detail::get_xc_integrator_ptr(integrator)->eval_fxc_contraction(
m, n,
density_matrix, ldp,
t_density_matrix, ldtp,
fxc, ldfxc,
- IntegratorSettingsXC{} );
+ ks_settings );
} catch (std::exception& e) {
detail::gauxc_status_handle(status, 1, e.what());
}
diff --git a/src/xc_integrator/local_work_driver/factory.cxx b/src/xc_integrator/local_work_driver/factory.cxx
index fd6b86ad..3547e1d3 100644
--- a/src/xc_integrator/local_work_driver/factory.cxx
@ -135,3 +241,443 @@ index fd6b86ad..3547e1d3 100644
(void)(settings);
switch(ex) {
diff --git a/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator.hpp b/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator.hpp
index 1287835e..91027539 100644
--- a/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator.hpp
+++ b/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator.hpp
@@ -128,7 +128,8 @@ protected:
const value_type* Py, int64_t ldpy,
const value_type* Px, int64_t ldpx,
host_task_iterator task_begin, host_task_iterator task_end,
- XCDeviceData& device_data, bool do_vxc );
+ XCDeviceData& device_data, bool do_vxc,
+ const IntegratorSettingsXC& settings );
void exc_vxc_local_work_( const basis_type& basis, const value_type* Ps, int64_t ldps,
const value_type* Pz, int64_t ldpz,
@@ -139,7 +140,7 @@ protected:
value_type* VXCy, int64_t ldvxcy,
value_type* VXCx, int64_t ldvxcx, value_type* EXC, value_type *N_EL,
host_task_iterator task_begin, host_task_iterator task_end,
- XCDeviceData& device_data );
+ XCDeviceData& device_data, const IntegratorSettingsXC& settings );
void pre_onedft_local_work_( const basis_type& basis, const value_type* Ps, int64_t ldps,
const value_type* Pz, int64_t ldpz,
@@ -159,7 +160,7 @@ protected:
const value_type* tPs, int64_t ldtps,
const value_type* tPz, int64_t ldtpz,
host_task_iterator task_begin, host_task_iterator task_end,
- XCDeviceData& device_data);
+ XCDeviceData& device_data, const IntegratorSettingsXC& settings);
void fxc_contraction_local_work_( const basis_type& basis, const value_type* Ps, int64_t ldps,
const value_type* Pz, int64_t ldpz,
@@ -169,7 +170,7 @@ protected:
value_type* FXCs, int64_t ldfxcs,
value_type* FXCz, int64_t ldfxcz,
host_task_iterator task_begin, host_task_iterator task_end,
- XCDeviceData& device_data );
+ XCDeviceData& device_data, const IntegratorSettingsXC& settings );
void eval_exc_grad_local_work_( const basis_type& basis, const value_type* Ps, int64_t ldps,
const value_type* Pz, int64_t ldpz,
diff --git a/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_exc.hpp b/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_exc.hpp
index 9a2a7cf4..b6f4a3da 100644
--- a/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_exc.hpp
+++ b/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_exc.hpp
@@ -59,7 +59,7 @@ void IncoreReplicatedXCDeviceIntegrator<ValueType>::
exc_vxc_local_work_( basis, Ps, ldps, Pz, ldpz, Py, ldpy, Px, ldpx,
// Passing nullptr for VXCs disables VXC entirely
nullptr, 0, nullptr, 0, nullptr, 0, nullptr, 0, EXC, &N_EL,
- tasks.begin(), tasks.end(), *device_data_ptr);
+ tasks.begin(), tasks.end(), *device_data_ptr, settings);
});
GAUXC_MPI_CODE(
@@ -100,4 +100,3 @@ void IncoreReplicatedXCDeviceIntegrator<ValueType>::
}
}
-
diff --git a/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_exc_grad.hpp b/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_exc_grad.hpp
index 6c030bc2..15230d68 100644
--- a/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_exc_grad.hpp
+++ b/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_exc_grad.hpp
@@ -168,6 +168,10 @@ void IncoreReplicatedXCDeviceIntegrator<ValueType>::
IntegratorSettingsEXC_GRAD exc_grad_settings;
if( auto* tmp = dynamic_cast<const IntegratorSettingsEXC_GRAD*>(&settings) ) {
exc_grad_settings = *tmp;
+ } else if( auto* ks_tmp = dynamic_cast<const IntegratorSettingsKS*>(&settings) ) {
+ exc_grad_settings.gks_dtol = ks_tmp->gks_dtol;
+ exc_grad_settings.rks_density_matrix_is_spin_summed =
+ ks_tmp->rks_density_matrix_is_spin_summed;
}
// Check that Partition Weights have been calculated
@@ -221,7 +225,8 @@ void IncoreReplicatedXCDeviceIntegrator<ValueType>::
else lwd->eval_collocation_gradient( &device_data );
// Evaluate X matrix and V vars
- const auto xmat_fac = is_rks ? 2.0 : 1.0;
+ const auto xmat_fac =
+ (is_rks and not exc_grad_settings.rks_density_matrix_is_spin_summed) ? 2.0 : 1.0;
const auto need_lapl = func.needs_laplacian();
const auto need_xmat_grad = not func.is_lda();
auto do_xmat_vvar = [&](density_id den_id) {
diff --git a/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_exc_vxc.hpp b/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_exc_vxc.hpp
index 6a27521d..42723b94 100644
--- a/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_exc_vxc.hpp
+++ b/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_exc_vxc.hpp
@@ -110,8 +110,8 @@ void IncoreReplicatedXCDeviceIntegrator<ValueType>::
// If we can do reductions on the device (e.g. NCCL)
// Don't communicate data back to the host before reduction
this->timer_.time_op("XCIntegrator.LocalWork_EXC_VXC", [&](){
- exc_vxc_local_work_( basis, Ps, ldps, Pz, ldpz, Py, ldpy, Px, ldpx, tasks.begin(), tasks.end(),
- *device_data_ptr, true);
+ exc_vxc_local_work_( basis, Ps, ldps, Pz, ldpz, Py, ldpy, Px, ldpx, tasks.begin(), tasks.end(),
+ *device_data_ptr, true, settings);
});
GAUXC_MPI_CODE(
@@ -177,7 +177,7 @@ void IncoreReplicatedXCDeviceIntegrator<ValueType>::
this->timer_.time_op("XCIntegrator.LocalWork_EXC_VXC", [&](){
exc_vxc_local_work_( basis, Ps, ldps, Pz, ldpz, Py, ldpy, Px, ldpx,
VXCs, ldvxcs, VXCz, ldvxcz, VXCy, ldvxcy, VXCx, ldvxcx, EXC,
- &N_EL, tasks.begin(), tasks.end(), *device_data_ptr);
+ &N_EL, tasks.begin(), tasks.end(), *device_data_ptr, settings);
});
GAUXC_MPI_CODE(
@@ -225,7 +225,8 @@ void IncoreReplicatedXCDeviceIntegrator<ValueType>::
const value_type* Py, int64_t ldpy,
const value_type* Px, int64_t ldpx,
host_task_iterator task_begin, host_task_iterator task_end,
- XCDeviceData& device_data, bool do_vxc ) {
+ XCDeviceData& device_data, bool do_vxc,
+ const IntegratorSettingsXC& settings ) {
const bool is_gks = (Pz != nullptr) and (Py != nullptr) and (Px != nullptr);
const bool is_uks = (Pz != nullptr) and (Py == nullptr) and (Px == nullptr);
const bool is_rks = (Ps != nullptr) and (not is_uks and not is_gks);
@@ -243,6 +244,11 @@ void IncoreReplicatedXCDeviceIntegrator<ValueType>::
if( func.is_mgga() and is_gks ) GAUXC_GENERIC_EXCEPTION("GKS mGGAs NYI!");
+ IntegratorSettingsKS ks_settings;
+ if( auto* tmp = dynamic_cast<const IntegratorSettingsKS*>(&settings) ) {
+ ks_settings = *tmp;
+ }
+
// Get basis map
BasisSetMap basis_map(basis,mol);
@@ -312,7 +318,8 @@ void IncoreReplicatedXCDeviceIntegrator<ValueType>::
else if( func.is_gga() ) lwd->eval_collocation_gradient( &device_data );
else lwd->eval_collocation( &device_data );
- const double xmat_fac = is_rks ? 2.0 : 1.0;
+ const double xmat_fac =
+ (is_rks and not ks_settings.rks_density_matrix_is_spin_summed) ? 2.0 : 1.0;
const bool need_xmat_grad = func.is_mgga();
// Evaluate X matrix and V vars
@@ -396,11 +403,12 @@ void IncoreReplicatedXCDeviceIntegrator<ValueType>::
value_type* VXCy, int64_t ldvxcy,
value_type* VXCx, int64_t ldvxcx, value_type* EXC, value_type *N_EL,
host_task_iterator task_begin, host_task_iterator task_end,
- XCDeviceData& device_data ) {
+ XCDeviceData& device_data, const IntegratorSettingsXC& settings ) {
// Get integrate and keep data on device
const bool do_vxc = VXCs;
- exc_vxc_local_work_( basis, Ps, ldps, Pz, ldpz, Py, ldpy, Px, ldpx, task_begin, task_end, device_data, do_vxc );
+ exc_vxc_local_work_( basis, Ps, ldps, Pz, ldpz, Py, ldpy, Px, ldpx,
+ task_begin, task_end, device_data, do_vxc, settings );
auto rt = detail::as_device_runtime(this->load_balancer_->runtime());
rt.device_backend()->master_queue_synchronize();
@@ -414,4 +422,3 @@ void IncoreReplicatedXCDeviceIntegrator<ValueType>::
}
}
-
diff --git a/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_fxc_contraction.hpp b/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_fxc_contraction.hpp
index ffc0ca41..0b813924 100644
--- a/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_fxc_contraction.hpp
+++ b/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_fxc_contraction.hpp
@@ -85,7 +85,7 @@ namespace GauXC::detail {
// Don't communicate data back to the host before reduction
this->timer_.time_op("XCIntegrator.LocalWork_FXC", [&](){
fxc_contraction_local_work_( basis, Ps, ldps, Pz, ldpz, tPs, ldtps, tPz, ldtpz,
- tasks.begin(), tasks.end(), *device_data_ptr);
+ tasks.begin(), tasks.end(), *device_data_ptr, ks_settings);
});
GAUXC_MPI_CODE(
@@ -127,7 +127,8 @@ namespace GauXC::detail {
// data from device
this->timer_.time_op("XCIntegrator.LocalWork_FXC", [&](){
fxc_contraction_local_work_( basis, Ps, ldps, Pz, ldpz, tPs, ldtps, tPz, ldtpz, &N_EL,
- FXCs, ldfxcs, FXCz, ldfxcz, tasks.begin(), tasks.end(), *device_data_ptr);
+ FXCs, ldfxcs, FXCz, ldfxcz, tasks.begin(), tasks.end(), *device_data_ptr,
+ ks_settings);
});
GAUXC_MPI_CODE(
@@ -160,7 +161,7 @@ namespace GauXC::detail {
const value_type* tPs, int64_t ldtps,
const value_type* tPz, int64_t ldtpz,
host_task_iterator task_begin, host_task_iterator task_end,
- XCDeviceData& device_data) {
+ XCDeviceData& device_data, const IntegratorSettingsXC& settings) {
const bool is_uks = (Pz != nullptr);
const bool is_rks = !is_uks;
if (not is_rks and not is_uks) {
@@ -175,6 +176,11 @@ namespace GauXC::detail {
const auto& func = *this->func_;
const auto& mol = this->load_balancer_->molecule();
+ IntegratorSettingsKS ks_settings;
+ if( auto* tmp = dynamic_cast<const IntegratorSettingsKS*>(&settings) ) {
+ ks_settings = *tmp;
+ }
+
// Get basis map
BasisSetMap basis_map(basis,mol);
@@ -243,7 +249,8 @@ namespace GauXC::detail {
else if( func.is_gga() ) lwd->eval_collocation_gradient( &device_data );
else lwd->eval_collocation( &device_data );
- const double xmat_fac = is_rks ? 2.0 : 1.0;
+ const double xmat_fac =
+ (is_rks and not ks_settings.rks_density_matrix_is_spin_summed) ? 2.0 : 1.0;
const bool need_xmat_grad = func.is_mgga();
// Evaluate X matrix and V vars
@@ -327,11 +334,11 @@ namespace GauXC::detail {
value_type* FXCs, int64_t ldfxcs,
value_type* FXCz, int64_t ldfxcz,
host_task_iterator task_begin, host_task_iterator task_end,
- XCDeviceData& device_data ) {
+ XCDeviceData& device_data, const IntegratorSettingsXC& settings ) {
// Get integrate and keep data on device
fxc_contraction_local_work_( basis, Ps, ldps, Pz, ldpz, tPs, ldtps, tPz, ldtpz,
- task_begin, task_end, device_data);
+ task_begin, task_end, device_data, settings);
auto rt = detail::as_device_runtime(this->load_balancer_->runtime());
rt.device_backend()->master_queue_synchronize();
diff --git a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_exc_grad.hpp b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_exc_grad.hpp
index f04ae24b..b3c5db95 100644
--- a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_exc_grad.hpp
+++ b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_exc_grad.hpp
@@ -126,6 +126,10 @@ void ReferenceReplicatedXCHostIntegrator<ValueType>::
IntegratorSettingsEXC_GRAD exc_grad_settings;
if( auto* tmp = dynamic_cast<const IntegratorSettingsEXC_GRAD*>(&settings) ) {
exc_grad_settings = *tmp;
+ } else if( auto* ks_tmp = dynamic_cast<const IntegratorSettingsKS*>(&settings) ) {
+ exc_grad_settings.gks_dtol = ks_tmp->gks_dtol;
+ exc_grad_settings.rks_density_matrix_is_spin_summed =
+ ks_tmp->rks_density_matrix_is_spin_summed;
}
// Get basis map
@@ -330,7 +334,8 @@ void ReferenceReplicatedXCHostIntegrator<ValueType>::
// Evaluate X matrix (2 * P * B/Bx/By/Bz) -> store in Z
// XXX: This assumes that bfn + gradients are contiguous in memory
- const auto xmat_fac = is_rks ? 2.0 : 1.0;
+ const auto xmat_fac =
+ (is_rks and not exc_grad_settings.rks_density_matrix_is_spin_summed) ? 2.0 : 1.0;
const int xmat_len = func.is_lda() ? 1 : 4;
lwd->eval_xmat( xmat_len*npts, nbf, nbe, submat_map, xmat_fac, Ps, ldps, basis_eval, nbe,
xNmat, nbe, nbe_scr );
diff --git a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_exc_vxc.hpp b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_exc_vxc.hpp
index 29878c50..94774eb5 100644
--- a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_exc_vxc.hpp
+++ b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_exc_vxc.hpp
@@ -385,7 +385,8 @@ void ReferenceReplicatedXCHostIntegrator<ValueType>::
// Evaluate X matrix (fac * P * B) -> store in Z
- const auto xmat_fac = is_rks ? 2.0 : 1.0; // TODO Fix for spinor RKS input
+ const auto xmat_fac =
+ (is_rks and not ks_settings.rks_density_matrix_is_spin_summed) ? 2.0 : 1.0;
lwd->eval_xmat( mgga_dim_scal * npts, nbf, nbe, submat_map, xmat_fac, Ps, ldps, basis_eval, nbe,
zmat, nbe, nbe_scr );
diff --git a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_fxc_contraction.hpp b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_fxc_contraction.hpp
index 192fe0f8..473e56f7 100644
--- a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_fxc_contraction.hpp
+++ b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_fxc_contraction.hpp
@@ -387,7 +387,8 @@ void ReferenceReplicatedXCHostIntegrator<ValueType>::
// Evaluate X matrix (fac * P * B) -> store in Z
- const auto xmat_fac = is_rks ? 2.0 : 1.0; // TODO Fix for spinor RKS input
+ const auto xmat_fac =
+ (is_rks and not ks_settings.rks_density_matrix_is_spin_summed) ? 2.0 : 1.0;
lwd->eval_xmat( mgga_dim_scal * npts, nbf, nbe, submat_map, xmat_fac, Ps, ldps, basis_eval, nbe,
zmat, nbe, nbe_scr );
// X matrix for Pz
diff --git a/src/xc_integrator/shell_batched/shell_batched_replicated_xc_integrator.hpp b/src/xc_integrator/shell_batched/shell_batched_replicated_xc_integrator.hpp
index 40e9512f..ceaaee1a 100644
--- a/src/xc_integrator/shell_batched/shell_batched_replicated_xc_integrator.hpp
+++ b/src/xc_integrator/shell_batched/shell_batched_replicated_xc_integrator.hpp
@@ -139,7 +139,9 @@ protected:
value_type* VXCy, int64_t ldvxcy,
value_type* VXCx, int64_t ldvxcx,
value_type* EXC, value_type *N_EL,
- host_task_iterator task_begin, host_task_iterator task_end, incore_integrator_type& incore_integrator
+ host_task_iterator task_begin, host_task_iterator task_end,
+ incore_integrator_type& incore_integrator,
+ const IntegratorSettingsXC& ks_settings
);
@@ -152,7 +154,9 @@ protected:
value_type* VXCz, int64_t ldvxcz,
value_type* VXCy, int64_t ldvxcy,
value_type* VXCx, int64_t ldvxcx,
- value_type* EXC, value_type* N_EL, incore_integrator_type& incore_integrator);
+ value_type* EXC, value_type* N_EL,
+ incore_integrator_type& incore_integrator,
+ const IntegratorSettingsXC& ks_settings);
public:
template <typename... Args>
diff --git a/src/xc_integrator/shell_batched/shell_batched_replicated_xc_integrator_exc.hpp b/src/xc_integrator/shell_batched/shell_batched_replicated_xc_integrator_exc.hpp
index 2a5565c9..02126e12 100644
--- a/src/xc_integrator/shell_batched/shell_batched_replicated_xc_integrator_exc.hpp
+++ b/src/xc_integrator/shell_batched/shell_batched_replicated_xc_integrator_exc.hpp
@@ -38,7 +38,7 @@ void ShellBatchedReplicatedXCIntegrator<BaseIntegratorType, IncoreIntegratorType
const value_type* Pz, int64_t ldpz,
const value_type* Py, int64_t ldpy,
const value_type* Px, int64_t ldpx,
- value_type* EXC, const IntegratorSettingsXC& /*ks_settings*/) {
+ value_type* EXC, const IntegratorSettingsXC& ks_settings) {
const auto& basis = this->load_balancer_->basis();
@@ -84,7 +84,7 @@ void ShellBatchedReplicatedXCIntegrator<BaseIntegratorType, IncoreIntegratorType
this->timer_.time_op("XCIntegrator.LocalWork", [&](){
exc_vxc_local_work_( basis, Ps, ldps, Pz, ldpz, Py, ldpy, Px, ldpx,
nullptr, 0, nullptr, 0, nullptr, 0, nullptr, 0, EXC,
- &N_EL, tasks.begin(), tasks.end(), incore_integrator );
+ &N_EL, tasks.begin(), tasks.end(), incore_integrator, ks_settings );
});
// Release ownership of LWD back to this integrator instance
@@ -134,4 +134,3 @@ void ShellBatchedReplicatedXCIntegrator<BaseIntegratorType, IncoreIntegratorType
}
}
-
diff --git a/src/xc_integrator/shell_batched/shell_batched_replicated_xc_integrator_exc_grad.hpp b/src/xc_integrator/shell_batched/shell_batched_replicated_xc_integrator_exc_grad.hpp
index 1b594e31..caae56d1 100644
--- a/src/xc_integrator/shell_batched/shell_batched_replicated_xc_integrator_exc_grad.hpp
+++ b/src/xc_integrator/shell_batched/shell_batched_replicated_xc_integrator_exc_grad.hpp
@@ -22,7 +22,7 @@ void ShellBatchedReplicatedXCIntegrator<BaseIntegratorType, IncoreIntegratorType
eval_exc_grad_( int64_t m, int64_t n, const value_type* P, int64_t ldp, value_type* EXC_GRAD, const IntegratorSettingsXC& settings ) {
GAUXC_GENERIC_EXCEPTION("ShellBatched exc_grad NYI" );
- util::unused(m,n,P,ldp,EXC_GRAD);
+ util::unused(m,n,P,ldp,EXC_GRAD,settings);
}
template <typename BaseIntegratorType, typename IncoreIntegratorType>
@@ -31,7 +31,7 @@ void ShellBatchedReplicatedXCIntegrator<BaseIntegratorType, IncoreIntegratorType
const value_type* Pz, int64_t lpdz, value_type* EXC_GRAD, const IntegratorSettingsXC& settings ) {
GAUXC_GENERIC_EXCEPTION("ShellBatched exc_grad NYI" );
- util::unused(m,n,Ps,ldps,Pz,lpdz,EXC_GRAD);
+ util::unused(m,n,Ps,ldps,Pz,lpdz,EXC_GRAD,settings);
}
template <typename BaseIntegratorType, typename IncoreIntegratorType>
diff --git a/src/xc_integrator/shell_batched/shell_batched_replicated_xc_integrator_exc_vxc.hpp b/src/xc_integrator/shell_batched/shell_batched_replicated_xc_integrator_exc_vxc.hpp
index 3dd43f4d..961c3f46 100644
--- a/src/xc_integrator/shell_batched/shell_batched_replicated_xc_integrator_exc_vxc.hpp
+++ b/src/xc_integrator/shell_batched/shell_batched_replicated_xc_integrator_exc_vxc.hpp
@@ -42,7 +42,7 @@ void ShellBatchedReplicatedXCIntegrator<BaseIntegratorType, IncoreIntegratorType
value_type* VXCz, int64_t ldvxcz,
value_type* VXCy, int64_t ldvxcy,
value_type* VXCx, int64_t ldvxcx,
- value_type* EXC, const IntegratorSettingsXC& /*ks_settings*/) {
+ value_type* EXC, const IntegratorSettingsXC& ks_settings) {
const auto& basis = this->load_balancer_->basis();
@@ -98,7 +98,7 @@ void ShellBatchedReplicatedXCIntegrator<BaseIntegratorType, IncoreIntegratorType
this->timer_.time_op("XCIntegrator.LocalWork", [&](){
exc_vxc_local_work_( basis, Ps, ldps, Pz, ldpz, Py, ldpy, Px, ldpx,
VXCs, ldvxcs, VXCz, ldvxcz, VXCy, ldvxcy, VXCx, ldvxcx, EXC,
- &N_EL, tasks.begin(), tasks.end(), incore_integrator );
+ &N_EL, tasks.begin(), tasks.end(), incore_integrator, ks_settings );
});
// Release ownership of LWD back to this integrator instance
@@ -166,7 +166,8 @@ void ShellBatchedReplicatedXCIntegrator<BaseIntegratorType, IncoreIntegratorType
value_type* VXCx, int64_t ldvxcx,
value_type* EXC, value_type *N_EL,
host_task_iterator task_begin, host_task_iterator task_end,
- incore_integrator_type& incore_integrator ) {
+ incore_integrator_type& incore_integrator,
+ const IntegratorSettingsXC& ks_settings ) {
//incore_integrator.exc_vxc_local_work( basis, P, ldp, VXC, ldvxc, EXC, N_EL, task_begin, task_end, device_data );
//return;
@@ -227,7 +228,7 @@ void ShellBatchedReplicatedXCIntegrator<BaseIntegratorType, IncoreIntegratorType
// Execute task
execute_task_batch( next_task, basis, mol, Ps, ldps, Pz, ldpz,
Py, ldpy, Px, ldpx, VXCs, ldvxcs, VXCz, ldvxcz, VXCy, ldvxcy,
- VXCx, ldvxcx, EXC, N_EL, incore_integrator );
+ VXCx, ldvxcx, EXC, N_EL, incore_integrator, ks_settings );
};
@@ -287,7 +288,8 @@ void ShellBatchedReplicatedXCIntegrator<BaseIntegratorType, IncoreIntegratorType
value_type* VXCy, int64_t ldvxcy,
value_type* VXCx, int64_t ldvxcx,
value_type* EXC, value_type *N_EL,
- incore_integrator_type& incore_integrator ) {
+ incore_integrator_type& incore_integrator,
+ const IntegratorSettingsXC& ks_settings ) {
// Alias information
@@ -398,13 +400,13 @@ void ShellBatchedReplicatedXCIntegrator<BaseIntegratorType, IncoreIntegratorType
incore_integrator.exc_vxc_local_work( basis_subset, Ps_submat, nbe,
Pz_submat, nbe, Py_submat, nbe, Px_submat, nbe, VXCs_submat, nbe,
VXCz_submat, nbe, VXCy_submat, nbe, VXCx_submat, nbe,
- &EXC_tmp, &NEL_tmp, task_begin, task_end, *device_data_ptr_ );
+ &EXC_tmp, &NEL_tmp, task_begin, task_end, *device_data_ptr_, ks_settings );
} else if constexpr (not IncoreIntegratorType::is_device) {
#endif
incore_integrator.exc_vxc_local_work( basis_subset, Ps_submat, nbe,
Pz_submat, nbe, Py_submat, nbe, Px_submat, nbe, VXCs_submat, nbe,
VXCz_submat, nbe, VXCy_submat, nbe, VXCx_submat, nbe,
- &EXC_tmp, &NEL_tmp, IntegratorSettingsKS{}, task_begin, task_end );
+ &EXC_tmp, &NEL_tmp, ks_settings, task_begin, task_end );
#ifdef GAUXC_HAS_DEVICE
}
#endif
@@ -444,4 +446,3 @@ void ShellBatchedReplicatedXCIntegrator<BaseIntegratorType, IncoreIntegratorType
}
}
-
diff --git a/src/xc_integrator/shell_batched/shell_batched_replicated_xc_integrator_fxc_contraction.hpp b/src/xc_integrator/shell_batched/shell_batched_replicated_xc_integrator_fxc_contraction.hpp
index 289de960..6dc916cb 100644
--- a/src/xc_integrator/shell_batched/shell_batched_replicated_xc_integrator_fxc_contraction.hpp
+++ b/src/xc_integrator/shell_batched/shell_batched_replicated_xc_integrator_fxc_contraction.hpp
@@ -41,7 +41,7 @@ void ShellBatchedReplicatedXCIntegrator<BaseIntegratorType, IncoreIntegratorType
const IntegratorSettingsXC& ks_settings ) {
GAUXC_GENERIC_EXCEPTION("ShellBatched FXC contraction NYI");
util::unused(m,n,Ps,ldps,Pz,ldpz,tPs,ldtps,tPz,ldtpz,
- FXCs,ldfxcs,FXCz,ldfxcz);
+ FXCs,ldfxcs,FXCz,ldfxcz,ks_settings);
}