From 418b37c462ef1a73dd6b7a7a67b9e9215cacbae3 Mon Sep 17 00:00:00 2001 From: Dynamics of Condensed Matter <30792324+DCM-Uni-Paderborn@users.noreply.github.com> Date: Fri, 22 May 2026 09:41:49 +0200 Subject: [PATCH] GauXC: use finite-difference GAPW XC gradients for SKALA and high-l bases (#5263) Co-authored-by: Thomas D. Kuehne --- src/input_cp2k_xc.F | 3 +- src/qs_dispersion_utils.F | 3 +- src/xc/xc_gauxc_functional.F | 143 ++++++++++++++++-- src/xc/xc_gauxc_interface.F | 4 + .../QS/regtest-dft-vdw-corr-3/TEST_FILES.toml | 1 + .../regtest-dft-vdw-corr-3/dftd3bj_skala.inp | 66 ++++++++ .../regtest-gauxc/CH4_DIMER_GAUXC_PBE_D3.inp | 69 +++++++++ tests/QS/regtest-gauxc/TEST_FILES.toml | 1 + 8 files changed, 274 insertions(+), 16 deletions(-) create mode 100644 tests/QS/regtest-dft-vdw-corr-3/dftd3bj_skala.inp create mode 100644 tests/QS/regtest-gauxc/CH4_DIMER_GAUXC_PBE_D3.inp diff --git a/src/input_cp2k_xc.F b/src/input_cp2k_xc.F index 69f0bb8e1a..d1f8d59cc4 100644 --- a/src/input_cp2k_xc.F +++ b/src/input_cp2k_xc.F @@ -1002,7 +1002,8 @@ CONTAINS "For available DFT-D4 functionals and parameters see: "// & ". "// & "For available D3 and D3(BJ) parameters see: "// & - ".", & + ". "// & + "For Skala with D3(BJ), SKALA selects the B3LYP D3(BJ) parameters.", & usage="REFERENCE_FUNCTIONAL ", & type_of_var=char_t) CALL section_add_keyword(subsection, keyword) diff --git a/src/qs_dispersion_utils.F b/src/qs_dispersion_utils.F index 24a1fa41fc..d2ff6787c2 100644 --- a/src/qs_dispersion_utils.F +++ b/src/qs_dispersion_utils.F @@ -962,7 +962,7 @@ CONTAINS a1 = 0.4724_dp s8 = 3.5681_dp a2 = 4.9858_dp - CASE ("B3LYP") + CASE ("B3LYP", "SKALA", "SKALA-1.1", "SKALA1.1") s6 = 1.0000_dp a1 = 0.3981_dp s8 = 1.9889_dp @@ -1233,4 +1233,3 @@ CONTAINS ! ************************************************************************************************** END MODULE qs_dispersion_utils - diff --git a/src/xc/xc_gauxc_functional.F b/src/xc/xc_gauxc_functional.F index 8ff01e4c74..01faa3b034 100644 --- a/src/xc/xc_gauxc_functional.F +++ b/src/xc/xc_gauxc_functional.F @@ -20,7 +20,7 @@ MODULE xc_gauxc_functional USE cp_log_handling, ONLY: cp_logger_get_default_io_unit USE external_potential_types, ONLY: gth_potential_type,& sgp_potential_type - USE input_constants, ONLY: xc_vdw_fun_none + USE input_constants, ONLY: xc_vdw_fun_nonloc USE input_section_types, ONLY: section_vals_get_subs_vals,& section_vals_get_subs_vals2,& section_vals_type,& @@ -456,6 +456,93 @@ CONTAINS END SUBROUTINE gauxc_xc_energy_for_particles +! ************************************************************************************************** +!> \brief compute a finite-difference GauXC XC nuclear gradient at fixed density +!> \param particle_set ... +!> \param qs_kind_set ... +!> \param density_scalar ... +!> \param nspins ... +!> \param model_name ... +!> \param xc_fun_name ... +!> \param grid_type ... +!> \param radial_quadrature ... +!> \param pruning_scheme ... +!> \param lb_exec_space ... +!> \param int_exec_space ... +!> \param batch_size ... +!> \param dx ... +!> \param para_env ... +!> \param exc_grad ... +!> \param density_zeta ... +! ************************************************************************************************** + SUBROUTINE 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, batch_size, dx, para_env, exc_grad, density_zeta) + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set + REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: density_scalar + INTEGER, INTENT(IN) :: nspins + CHARACTER(len=*), INTENT(IN) :: model_name, xc_fun_name, grid_type, & + radial_quadrature, pruning_scheme, & + lb_exec_space, int_exec_space + INTEGER, INTENT(IN) :: batch_size + REAL(KIND=dp), INTENT(IN) :: dx + TYPE(mp_para_env_type), POINTER :: para_env + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & + INTENT(OUT) :: exc_grad + REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), & + OPTIONAL :: density_zeta + + INTEGER :: iatom, idir + REAL(KIND=dp) :: xc_minus, xc_plus + TYPE(particle_type), ALLOCATABLE, DIMENSION(:) :: particle_set_minus, particle_set_plus + + CPASSERT(ASSOCIATED(particle_set)) + CPASSERT(dx > 0.0_dp) + + ALLOCATE (exc_grad(3*SIZE(particle_set))) + exc_grad = 0.0_dp + + IF (para_env%mepos == 0) THEN + ALLOCATE (particle_set_minus(SIZE(particle_set)), particle_set_plus(SIZE(particle_set))) + + DO iatom = 1, SIZE(particle_set) + DO idir = 1, 3 + particle_set_minus = particle_set + particle_set_plus = particle_set + particle_set_minus(iatom)%r(idir) = particle_set_minus(iatom)%r(idir) - dx + particle_set_plus(iatom)%r(idir) = particle_set_plus(iatom)%r(idir) + dx + IF (PRESENT(density_zeta)) THEN + CALL gauxc_xc_energy_for_particles( & + particle_set_plus, qs_kind_set, density_scalar, nspins, model_name, & + xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, & + int_exec_space, batch_size, xc_plus, density_zeta=density_zeta) + CALL gauxc_xc_energy_for_particles( & + particle_set_minus, qs_kind_set, density_scalar, nspins, model_name, & + xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, & + int_exec_space, batch_size, xc_minus, density_zeta=density_zeta) + ELSE + CALL gauxc_xc_energy_for_particles( & + particle_set_plus, qs_kind_set, density_scalar, nspins, model_name, & + xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, & + int_exec_space, batch_size, xc_plus) + CALL gauxc_xc_energy_for_particles( & + particle_set_minus, qs_kind_set, density_scalar, nspins, model_name, & + xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, & + int_exec_space, batch_size, xc_minus) + END IF + exc_grad(3*iatom - 3 + idir) = (xc_plus - xc_minus)/(2.0_dp*dx) + END DO + END DO + + DEALLOCATE (particle_set_minus, particle_set_plus) + END IF + + CALL para_env%bcast(exc_grad, 0) + + END SUBROUTINE gauxc_xc_gradient_fd + ! ************************************************************************************************** !> \brief finite-difference check of the molecular GauXC XC virial diagnostic !> \param exc_grad ... @@ -574,7 +661,7 @@ CONTAINS CHARACTER(len=1), DIMENSION(3), PARAMETER :: label = ["x", "y", "z"] INTEGER :: i, iatom, iw, j - REAL(KIND=dp), DIMENSION(3) :: grad, grad_sum + REAL(KIND=dp), DIMENSION(3) :: center, displacement, grad, grad_sum REAL(KIND=dp), DIMENSION(3, 3) :: molecular_virial CPASSERT(ASSOCIATED(particle_set)) @@ -582,14 +669,21 @@ CONTAINS IF (para_env%mepos /= 0) RETURN + center = 0.0_dp + DO iatom = 1, SIZE(particle_set) + center = center + particle_set(iatom)%r + END DO + center = center/REAL(SIZE(particle_set), dp) + grad_sum = 0.0_dp molecular_virial = 0.0_dp DO iatom = 1, SIZE(particle_set) grad = exc_grad(3*iatom - 2:3*iatom) + displacement = particle_set(iatom)%r - center grad_sum = grad_sum + grad DO i = 1, 3 DO j = 1, 3 - molecular_virial(i, j) = molecular_virial(i, j) + grad(i)*particle_set(iatom)%r(j) + molecular_virial(i, j) = molecular_virial(i, j) + grad(i)*displacement(j) END DO END DO END DO @@ -678,11 +772,13 @@ CONTAINS TYPE(section_vals_type), INTENT(in), POINTER :: xc_section LOGICAL, INTENT(IN) :: calculate_forces - CHARACTER(len=*), PARAMETER :: abort_message = & - "GauXC does not support dispersion corrections. "// & - "Please disable the VDW_POTENTIAL section or remove the XC_FUNCTIONAL%GAUXC subsection.", & - gapw_xc_abort_message = "GauXC with METHOD GAPW_XC is not supported yet. "// & - "The GAPW_XC one-center XC correction needs a dedicated GauXC design." + CHARACTER(len=*), PARAMETER :: gapw_xc_abort_message = & + "GauXC with METHOD GAPW_XC is not supported yet. "// & + "The GAPW_XC one-center XC correction needs a dedicated GauXC design.", & + nonlocal_vdw_abort_message = & + "GauXC does not support non-local VDW_POTENTIAL corrections. "// & + "Use an additive PAIR_POTENTIAL dispersion correction or disable GauXC." + REAL(KIND=dp), PARAMETER :: gapw_fd_gradient_dx = 1.0E-4_dp CHARACTER(len=default_path_length) :: model_key, model_name CHARACTER(len=default_string_length) :: grid_key, grid_type, int_exec_space, lb_exec_space, & @@ -690,7 +786,7 @@ CONTAINS INTEGER :: batch_size, img, ispin, natom, nimages, & nspins LOGICAL :: grid_explicit, molecular_virial, molecular_virial_debug, pruning_explicit, & - use_gradient_self_runtime, use_onedft, use_self_runtime, use_skala_model + use_fd_gradient, use_gradient_self_runtime, use_onedft, use_self_runtime, use_skala_model 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 @@ -761,8 +857,8 @@ CONTAINS nspins = dft_control%nspins IF (ASSOCIATED(qs_env%dispersion_env)) THEN - IF (qs_env%dispersion_env%type /= xc_vdw_fun_none) THEN - CPABORT(abort_message) + IF (qs_env%dispersion_env%type == xc_vdw_fun_nonloc) THEN + CPABORT(nonlocal_vdw_abort_message) END IF END IF NULLIFY (vxc_zeta_tmp%matrix) @@ -860,6 +956,15 @@ CONTAINS particle_set, & gauxc_status) CALL gauxc_check_status(gauxc_status) + use_fd_gradient = dft_control%qs_control%gapw .AND. calculate_forces .AND. & + (use_skala_model .OR. gauxc_basis%max_l > 3) + IF (use_fd_gradient .AND. para_env%mepos == 0) THEN + CALL cp_warn( & + __LOCATION__, & + "Using finite-difference GauXC XC gradients for METHOD GAPW with SKALA or "// & + "all-electron 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. @@ -946,7 +1051,13 @@ CONTAINS model=TRIM(model_name)) CALL gauxc_check_status(gauxc_status) IF (calculate_forces) THEN - IF (use_gradient_self_runtime) THEN + IF (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, batch_size, gapw_fd_gradient_dx, & + para_env, exc_grad%exc_grad) + ELSE IF (use_gradient_self_runtime) THEN exc_grad = gauxc_compute_xc_gradient( & gauxc_gradient_integrator_result, & density_scalar, & @@ -1003,7 +1114,13 @@ CONTAINS model=TRIM(model_name)) CALL gauxc_check_status(gauxc_status) IF (calculate_forces) THEN - IF (use_gradient_self_runtime) THEN + IF (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, batch_size, gapw_fd_gradient_dx, & + para_env, exc_grad%exc_grad, density_zeta=density_zeta) + ELSE IF (use_gradient_self_runtime) THEN exc_grad = gauxc_compute_xc_gradient( & gauxc_gradient_integrator_result, & density_scalar, & diff --git a/src/xc/xc_gauxc_interface.F b/src/xc/xc_gauxc_interface.F index 6dcc2c3fa7..d29593bc43 100644 --- a/src/xc/xc_gauxc_interface.F +++ b/src/xc/xc_gauxc_interface.F @@ -134,6 +134,7 @@ MODULE xc_gauxc_interface END TYPE cp_gauxc_molecule_type TYPE cp_gauxc_basisset_type + INTEGER :: max_l = -1 END TYPE cp_gauxc_basisset_type TYPE cp_gauxc_grid_type @@ -155,6 +156,7 @@ MODULE xc_gauxc_interface TYPE cp_gauxc_basisset_type TYPE(gauxc_basisset_type) :: basis + INTEGER :: max_l = -1 END TYPE cp_gauxc_basisset_type TYPE cp_gauxc_grid_type @@ -382,6 +384,7 @@ CONTAINS ALLOCATE (shells(nshell_total)) shell_index = 0 + res%max_l = -1 DO iatom = 1, natoms ! for each atom atomic_kind => particle_set(iatom)%atomic_kind CALL get_atomic_kind(atomic_kind, kind_number=ikind) @@ -401,6 +404,7 @@ CONTAINS DO ishell = 1, gto_basis%nshell(iset) ! for each shell within the shell group shell_index = shell_index + 1 ! global shell index, flattened over atoms and groups lval = gto_basis%l(ishell, iset) + res%max_l = MAX(res%max_l, lval) shells(shell_index)%l = INT(lval, c_int32_t) ! FIXME hardcoded true param ! pure=1: spherical Gaussians; pure=0: cartesian Gaussians diff --git a/tests/QS/regtest-dft-vdw-corr-3/TEST_FILES.toml b/tests/QS/regtest-dft-vdw-corr-3/TEST_FILES.toml index 15bab70889..042e2361da 100644 --- a/tests/QS/regtest-dft-vdw-corr-3/TEST_FILES.toml +++ b/tests/QS/regtest-dft-vdw-corr-3/TEST_FILES.toml @@ -12,6 +12,7 @@ "argon-beef.inp" = [{matcher="E_total", tol=4e-13, ref=-42.46387625350550}] # BJ "dftd3bj_t1.inp" = [{matcher="M033", tol=1.0E-14, ref=-0.00355123783846}] +"dftd3bj_skala.inp" = [{matcher="M033", tol=1.0E-14, ref=-0.00559494624729}] "dftd3bj_t2.inp" = [{matcher="M033", tol=1.0E-14, ref=-0.05897356220363}] "dftd3bj_t3.inp" = [{matcher="M033", tol=1.0E-14, ref=-0.00112424003807}] "dftd3bj_t4.inp" = [{matcher="M007", tol=2.0E-12, ref=-84.2983390350}] diff --git a/tests/QS/regtest-dft-vdw-corr-3/dftd3bj_skala.inp b/tests/QS/regtest-dft-vdw-corr-3/dftd3bj_skala.inp new file mode 100644 index 0000000000..d9212ad92e --- /dev/null +++ b/tests/QS/regtest-dft-vdw-corr-3/dftd3bj_skala.inp @@ -0,0 +1,66 @@ +&GLOBAL + PRINT_LEVEL MEDIUM + PROJECT dftd3bj_skala + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + METHOD QS + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME POTENTIAL + &MGRID + CUTOFF 100 + NGRIDS 1 + &END MGRID + &QS + METHOD GPW + &END QS + &SCF + EPS_SCF 1.0e-0 + IGNORE_CONVERGENCE_FAILURE + MAX_SCF 1 + SCF_GUESS ATOMIC + &END SCF + &XC + &VDW_POTENTIAL + DISPERSION_FUNCTIONAL PAIR_POTENTIAL + &PAIR_POTENTIAL + PARAMETER_FILE_NAME dftd3.dat + REFERENCE_FUNCTIONAL SKALA + TYPE DFTD3(BJ) + &PRINT_DFTD + &END PRINT_DFTD + &END PAIR_POTENTIAL + &END VDW_POTENTIAL + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 6.0 6.0 6.0 + PERIODIC NONE + &END CELL + &COORD + C 7.499969 9.250001 7.500000 + C 7.499969 5.750001 7.500000 + H 8.205969 8.887002 6.753000 + H 6.499969 8.887002 7.262000 + H 7.793968 8.887002 8.485001 + H 7.499969 10.340001 7.500000 + H 7.199968 6.113001 6.517000 + H 8.501968 6.113001 7.731000 + H 6.798969 6.113001 8.252001 + H 7.499969 4.660001 7.500000 + &END COORD + &KIND H + BASIS_SET DZV-GTH + POTENTIAL GTH-PBE-q1 + &END KIND + &KIND C + BASIS_SET DZVP-GTH + POTENTIAL GTH-PBE-q4 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-gauxc/CH4_DIMER_GAUXC_PBE_D3.inp b/tests/QS/regtest-gauxc/CH4_DIMER_GAUXC_PBE_D3.inp new file mode 100644 index 0000000000..af481333b5 --- /dev/null +++ b/tests/QS/regtest-gauxc/CH4_DIMER_GAUXC_PBE_D3.inp @@ -0,0 +1,69 @@ +&GLOBAL + PRINT_LEVEL MEDIUM + PROJECT CH4_DIMER_GAUXC_PBE_D3 + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + METHOD QS + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 100 + NGRIDS 1 + &END MGRID + &QS + METHOD GPW + &END QS + &SCF + EPS_SCF 1.0e-0 + IGNORE_CONVERGENCE_FAILURE + MAX_SCF 1 + SCF_GUESS ATOMIC + &END SCF + &XC + &VDW_POTENTIAL + DISPERSION_FUNCTIONAL PAIR_POTENTIAL + &PAIR_POTENTIAL + PARAMETER_FILE_NAME dftd3.dat + REFERENCE_FUNCTIONAL PBE + TYPE DFTD3(BJ) + &PRINT_DFTD + &END PRINT_DFTD + &END PAIR_POTENTIAL + &END VDW_POTENTIAL + &XC_FUNCTIONAL + &GAUXC + FUNCTIONAL PBE + &END GAUXC + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 6.0 6.0 6.0 + PERIODIC NONE + &END CELL + &COORD + C 7.499969 9.250001 7.500000 + C 7.499969 5.750001 7.500000 + H 8.205969 8.887002 6.753000 + H 6.499969 8.887002 7.262000 + H 7.793968 8.887002 8.485001 + H 7.499969 10.340001 7.500000 + H 7.199968 6.113001 6.517000 + H 8.501968 6.113001 7.731000 + H 6.798969 6.113001 8.252001 + H 7.499969 4.660001 7.500000 + &END COORD + &KIND H + BASIS_SET DZV-GTH + POTENTIAL GTH-PBE-q1 + &END KIND + &KIND C + BASIS_SET DZVP-GTH + POTENTIAL GTH-PBE-q4 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-gauxc/TEST_FILES.toml b/tests/QS/regtest-gauxc/TEST_FILES.toml index 8bde4237c9..145b19af84 100644 --- a/tests/QS/regtest-gauxc/TEST_FILES.toml +++ b/tests/QS/regtest-gauxc/TEST_FILES.toml @@ -19,3 +19,4 @@ "NH3_ONEDFT_PBE_REFERENCE.inp" = [{matcher="E_total", tol=1e-9, ref=-11.722432805445091}] "NH3_ONEDFT_PBE.inp" = [{matcher="E_total", tol=1e-9, ref=-11.722558568119791}] "NH3_ONEDFT_PBE_FORCE_DEBUG.inp" = [{matcher="DEBUG_force_sum", tol=5e-5, ref=0.0}] +"CH4_DIMER_GAUXC_PBE_D3.inp" = [{matcher="M033", tol=1e-14, ref=-0.00355123783846}]