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

GauXC: use finite-difference GAPW XC gradients for SKALA and high-l bases (#5263)

Browse source 
Co-authored-by: Thomas D. Kuehne <tkuehne@cp2k.org>
This commit is contained in:
Dynamics of Condensed Matter 2026-05-22 09:41:49 +02:00 • committed by GitHub
parent 919a8b9f61
commit 418b37c462
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
8 changed files with 274 additions and 16 deletions
Download patch file Download diff file Expand all files Collapse all files

3
src/input_cp2k_xc.F
View file

@ -1002,7 +1002,8 @@ CONTAINS
"For available DFT-D4 functionals and parameters see: "// &
"<https://github.com/dftd4/dftd4/blob/main/src/dftd4/param.f90>. "// &
"For available D3 and D3(BJ) parameters see: "// &
"<https://www.chemie.uni-bonn.de/grimme/de/software/dft-d3>.", &
"<https://www.chemie.uni-bonn.de/grimme/de/software/dft-d3>. "// &
"For Skala with D3(BJ), SKALA selects the B3LYP D3(BJ) parameters.", &
usage="REFERENCE_FUNCTIONAL <functional>", &
type_of_var=char_t)
CALL section_add_keyword(subsection, keyword)

3
src/qs_dispersion_utils.F
View file

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

143
src/xc/xc_gauxc_functional.F
View file

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

4
src/xc/xc_gauxc_interface.F
View file

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

1
tests/QS/regtest-dft-vdw-corr-3/TEST_FILES.toml
View file

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

66
tests/QS/regtest-dft-vdw-corr-3/dftd3bj_skala.inp Normal file
View file

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

69
tests/QS/regtest-gauxc/CH4_DIMER_GAUXC_PBE_D3.inp Normal file
View file

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

1
tests/QS/regtest-gauxc/TEST_FILES.toml
View file

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