diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 21439a11e8..a31daadb31 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -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 diff --git a/src/qs_environment_types.F b/src/qs_environment_types.F index ac6c8fa2a6..a7391d5956 100644 --- a/src/qs_environment_types.F +++ b/src/qs_environment_types.F @@ -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 diff --git a/src/qs_ks_methods.F b/src/qs_ks_methods.F index 652655dae1..9f6a420799 100644 --- a/src/qs_ks_methods.F +++ b/src/qs_ks_methods.F @@ -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 diff --git a/src/qs_vxc_atom.F b/src/qs_vxc_atom.F index 77ca466195..84b2bd4516 100644 --- a/src/qs_vxc_atom.F +++ b/src/qs_vxc_atom.F @@ -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 diff --git a/src/skala_gpw_functional.F b/src/skala_gpw_functional.F index ea71c9ead5..212747e875 100644 --- a/src/skala_gpw_functional.F +++ b/src/skala_gpw_functional.F @@ -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 diff --git a/src/torch_c_api.cpp b/src/torch_c_api.cpp index 9349653ebf..9293fd6425 100644 --- a/src/torch_c_api.cpp +++ b/src/torch_c_api.cpp @@ -7,6 +7,7 @@ #if defined(__LIBTORCH) +#include #include #include #include @@ -16,6 +17,7 @@ #include #include +#include #include #include #include @@ -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(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)); diff --git a/src/xc/xc_derivatives.F b/src/xc/xc_derivatives.F index 431bb25b6e..846f254921 100644 --- a/src/xc/xc_derivatives.F +++ b/src/xc/xc_derivatives.F @@ -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 diff --git a/src/xc/xc_gauxc_cache.F b/src/xc/xc_gauxc_cache.F new file mode 100644 index 0000000000..5da9d0aee3 --- /dev/null +++ b/src/xc/xc_gauxc_cache.F @@ -0,0 +1,430 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2026 CP2K developers group ! +! ! +! 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 diff --git a/src/xc/xc_gauxc_functional.F b/src/xc/xc_gauxc_functional.F index 3c12631d2f..50adf633b9 100644 --- a/src/xc/xc_gauxc_functional.F +++ b/src/xc/xc_gauxc_functional.F @@ -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 diff --git a/src/xc/xc_gauxc_interface.F b/src/xc/xc_gauxc_interface.F index ed1356ebde..dadfbc0c99 100644 --- a/src/xc/xc_gauxc_interface.F +++ b/src/xc/xc_gauxc_interface.F @@ -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, & diff --git a/src/xc/xc_libxc.F b/src/xc/xc_libxc.F index 3faed929f3..4c58f994a0 100644 --- a/src/xc/xc_libxc.F +++ b/src/xc/xc_libxc.F @@ -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, "// & diff --git a/tests/QS/GAPW_ACCURATE_XCINT_COVERAGE.md b/tests/QS/GAPW_ACCURATE_XCINT_COVERAGE.md index c7e9b8747a..b369135c48 100644 --- a/tests/QS/GAPW_ACCURATE_XCINT_COVERAGE.md +++ b/tests/QS/GAPW_ACCURATE_XCINT_COVERAGE.md @@ -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`. diff --git a/tests/QS/regtest-acc-1/TEST_FILES.toml b/tests/QS/regtest-acc-1/TEST_FILES.toml index 41e4817c71..d26801d538 100644 --- a/tests/QS/regtest-acc-1/TEST_FILES.toml +++ b/tests/QS/regtest-acc-1/TEST_FILES.toml @@ -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 diff --git a/tests/QS/regtest-acc-1/h2o-uzh-gapw-force-1.inp b/tests/QS/regtest-acc-1/h2o-uzh-gapw-force-1.inp new file mode 100644 index 0000000000..9e6684b46a --- /dev/null +++ b/tests/QS/regtest-acc-1/h2o-uzh-gapw-force-1.inp @@ -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 diff --git a/tests/QS/regtest-acc-1/h2o-uzh-gapw-stress-debug-1.inp b/tests/QS/regtest-acc-1/h2o-uzh-gapw-stress-debug-1.inp new file mode 100644 index 0000000000..920bbbf7a1 --- /dev/null +++ b/tests/QS/regtest-acc-1/h2o-uzh-gapw-stress-debug-1.inp @@ -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 diff --git a/tests/QS/regtest-acc-1/h2o-uzh-gth-gapw-force-1.inp b/tests/QS/regtest-acc-1/h2o-uzh-gth-gapw-force-1.inp new file mode 100644 index 0000000000..924f9670d2 --- /dev/null +++ b/tests/QS/regtest-acc-1/h2o-uzh-gth-gapw-force-1.inp @@ -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 diff --git a/tests/QS/regtest-acc-1/h2o-uzh-gth-gapw-stress-debug-1.inp b/tests/QS/regtest-acc-1/h2o-uzh-gth-gapw-stress-debug-1.inp new file mode 100644 index 0000000000..3c80d1546e --- /dev/null +++ b/tests/QS/regtest-acc-1/h2o-uzh-gth-gapw-stress-debug-1.inp @@ -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 diff --git a/tests/QS/regtest-acc-1/h2o-uzh-gth-gapw_xc-force-1.inp b/tests/QS/regtest-acc-1/h2o-uzh-gth-gapw_xc-force-1.inp new file mode 100644 index 0000000000..412457eafe --- /dev/null +++ b/tests/QS/regtest-acc-1/h2o-uzh-gth-gapw_xc-force-1.inp @@ -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 diff --git a/tests/QS/regtest-acc-1/h2o-uzh-gth-gapw_xc-stress-debug-1.inp b/tests/QS/regtest-acc-1/h2o-uzh-gth-gapw_xc-stress-debug-1.inp new file mode 100644 index 0000000000..5bceffa537 --- /dev/null +++ b/tests/QS/regtest-acc-1/h2o-uzh-gth-gapw_xc-stress-debug-1.inp @@ -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 diff --git a/tests/QS/regtest-ecp/ICl_lanl2dz_gapw_force.inp b/tests/QS/regtest-ecp/ICl_lanl2dz_gapw_force.inp new file mode 100644 index 0000000000..35a8618cc7 --- /dev/null +++ b/tests/QS/regtest-ecp/ICl_lanl2dz_gapw_force.inp @@ -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 diff --git a/tests/QS/regtest-ecp/ICl_lanl2dz_gapw_stress.inp b/tests/QS/regtest-ecp/ICl_lanl2dz_gapw_stress.inp new file mode 100644 index 0000000000..aa73463cd3 --- /dev/null +++ b/tests/QS/regtest-ecp/ICl_lanl2dz_gapw_stress.inp @@ -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 diff --git a/tests/QS/regtest-ecp/ICl_lanl2dz_gapw_xc_force.inp b/tests/QS/regtest-ecp/ICl_lanl2dz_gapw_xc_force.inp new file mode 100644 index 0000000000..a1fd01ce1e --- /dev/null +++ b/tests/QS/regtest-ecp/ICl_lanl2dz_gapw_xc_force.inp @@ -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 diff --git a/tests/QS/regtest-ecp/ICl_lanl2dz_gapw_xc_stress.inp b/tests/QS/regtest-ecp/ICl_lanl2dz_gapw_xc_stress.inp new file mode 100644 index 0000000000..b94b1a47a1 --- /dev/null +++ b/tests/QS/regtest-ecp/ICl_lanl2dz_gapw_xc_stress.inp @@ -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 diff --git a/tests/QS/regtest-ecp/TEST_FILES.toml b/tests/QS/regtest-ecp/TEST_FILES.toml index 77e7b87d8e..914bf20fb4 100644 --- a/tests/QS/regtest-ecp/TEST_FILES.toml +++ b/tests/QS/regtest-ecp/TEST_FILES.toml @@ -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" = [] diff --git a/tests/QS/regtest-gauxc-api/TEST_FILES.toml b/tests/QS/regtest-gauxc-api/TEST_FILES.toml index c03ead6223..e9ee7d1cee 100644 --- a/tests/QS/regtest-gauxc-api/TEST_FILES.toml +++ b/tests/QS/regtest-gauxc-api/TEST_FILES.toml @@ -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}] diff --git a/tests/QS/regtest-gauxc-cdft/TEST_FILES.toml b/tests/QS/regtest-gauxc-cdft/TEST_FILES.toml index 3bd87ffe71..db816e3a27 100644 --- a/tests/QS/regtest-gauxc-cdft/TEST_FILES.toml +++ b/tests/QS/regtest-gauxc-cdft/TEST_FILES.toml @@ -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}] diff --git a/tests/QS/regtest-gauxc-gapw-ecp/TEST_FILES.toml b/tests/QS/regtest-gauxc-gapw-ecp/TEST_FILES.toml index 274511571a..babec12d78 100644 --- a/tests/QS/regtest-gauxc-gapw-ecp/TEST_FILES.toml +++ b/tests/QS/regtest-gauxc-gapw-ecp/TEST_FILES.toml @@ -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}] diff --git a/tests/QS/regtest-gauxc-gapw-gth-kp/TEST_FILES.toml b/tests/QS/regtest-gauxc-gapw-gth-kp/TEST_FILES.toml index 80fac5377e..dea40ed44e 100644 --- a/tests/QS/regtest-gauxc-gapw-gth-kp/TEST_FILES.toml +++ b/tests/QS/regtest-gauxc-gapw-gth-kp/TEST_FILES.toml @@ -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}] diff --git a/tests/QS/regtest-gauxc-gapw-gth/TEST_FILES.toml b/tests/QS/regtest-gauxc-gapw-gth/TEST_FILES.toml index 300b8bb9d2..fbeaab6fca 100644 --- a/tests/QS/regtest-gauxc-gapw-gth/TEST_FILES.toml +++ b/tests/QS/regtest-gauxc-gapw-gth/TEST_FILES.toml @@ -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}] diff --git a/tools/toolchain/scripts/stage6/gauxc-1.1-skala-cp2k-fixes.patch b/tools/toolchain/scripts/stage6/gauxc-1.1-skala-cp2k-fixes.patch index c129f3c11c..62a7ddf812 100644 --- a/tools/toolchain/scripts/stage6/gauxc-1.1-skala-cp2k-fixes.patch +++ b/tools/toolchain/scripts/stage6/gauxc-1.1-skala-cp2k-fixes.patch @@ -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:: + 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:: + + } + } +- +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:: + IntegratorSettingsEXC_GRAD exc_grad_settings; + if( auto* tmp = dynamic_cast(&settings) ) { + exc_grad_settings = *tmp; ++ } else if( auto* ks_tmp = dynamic_cast(&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:: + 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:: + // 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:: + 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:: + 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:: + + if( func.is_mgga() and is_gks ) GAUXC_GENERIC_EXCEPTION("GKS mGGAs NYI!"); + ++ IntegratorSettingsKS ks_settings; ++ if( auto* tmp = dynamic_cast(&settings) ) { ++ ks_settings = *tmp; ++ } ++ + // Get basis map + BasisSetMap basis_map(basis,mol); + +@@ -312,7 +318,8 @@ void IncoreReplicatedXCDeviceIntegrator:: + 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:: + 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:: + + } + } +- +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(&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:: + IntegratorSettingsEXC_GRAD exc_grad_settings; + if( auto* tmp = dynamic_cast(&settings) ) { + exc_grad_settings = *tmp; ++ } else if( auto* ks_tmp = dynamic_cast(&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:: + + // 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:: + + + // 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:: + + + // 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 +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 ShellBatchedReplicatedXCIntegratorload_balancer_->basis(); +@@ -84,7 +84,7 @@ void ShellBatchedReplicatedXCIntegratortimer_.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 +@@ -31,7 +31,7 @@ void ShellBatchedReplicatedXCIntegrator +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 ShellBatchedReplicatedXCIntegratorload_balancer_->basis(); +@@ -98,7 +98,7 @@ void ShellBatchedReplicatedXCIntegratortimer_.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