Skip to content

Commit

Permalink
fxc potential (#2626)
Browse files Browse the repository at this point in the history
  • Loading branch information
@juerghutter
juerghutter committed Feb 24, 2023
1 parent be1541d commit 90e3ec8
Show file tree
Hide file tree
Showing 16 changed files with 618 additions and 64 deletions.
4 changes: 4 additions & 0 deletions src/input_constants.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -561,6 +561,10 @@ MODULE input_constants
xc_funct_b3lyp = 10, &
xc_funct_pbe0 = 11, &
xc_funct_beefvdw = 12
INTEGER, PARAMETER, PUBLIC :: fxc_none = 0, &
fxc_funct_pade = 1, &
fxc_funct_lda = 2, &
fxc_funct_gga = 3
INTEGER, PARAMETER, PUBLIC :: sic_none = 0, &
sic_mauri_us = 1, &
sic_mauri_spz = 2, &
Expand Down
49 changes: 27 additions & 22 deletions src/input_cp2k_xc.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,12 +20,13 @@ MODULE input_cp2k_xc
cp_print_key_section_create,&
high_print_level
USE input_constants, ONLY: &
do_adiabatic_hybrid_mcy3, do_adiabatic_model_pade, gaussian, slater, vdw_nl_drsll, &
vdw_nl_lmkll, vdw_nl_rvv10, vdw_pairpot_dftd2, vdw_pairpot_dftd3, vdw_pairpot_dftd3bj, &
xc_funct_b3lyp, xc_funct_beefvdw, xc_funct_blyp, xc_funct_bp, xc_funct_hcth120, &
xc_funct_no_shortcut, xc_funct_olyp, xc_funct_pade, xc_funct_pbe, xc_funct_pbe0, &
xc_funct_tpss, xc_none, xc_pot_energy_none, xc_pot_energy_sum_eigenvalues, &
xc_pot_energy_xc_functional, xc_vdw_fun_none, xc_vdw_fun_nonloc, xc_vdw_fun_pairpot
do_adiabatic_hybrid_mcy3, do_adiabatic_model_pade, fxc_funct_gga, fxc_funct_lda, &
fxc_funct_pade, fxc_none, gaussian, slater, vdw_nl_drsll, vdw_nl_lmkll, vdw_nl_rvv10, &
vdw_pairpot_dftd2, vdw_pairpot_dftd3, vdw_pairpot_dftd3bj, xc_funct_b3lyp, &
xc_funct_beefvdw, xc_funct_blyp, xc_funct_bp, xc_funct_hcth120, xc_funct_no_shortcut, &
xc_funct_olyp, xc_funct_pade, xc_funct_pbe, xc_funct_pbe0, xc_funct_tpss, xc_none, &
xc_pot_energy_none, xc_pot_energy_sum_eigenvalues, xc_pot_energy_xc_functional, &
xc_vdw_fun_none, xc_vdw_fun_nonloc, xc_vdw_fun_pairpot
USE input_cp2k_hfx, ONLY: create_hfx_section
USE input_cp2k_mp2, ONLY: create_mp2_section
USE input_keyword_types, ONLY: keyword_create,&
Expand Down Expand Up @@ -783,7 +784,6 @@ SUBROUTINE create_xc_kernel_section(section)
TYPE(section_type), POINTER :: section

TYPE(keyword_type), POINTER :: keyword
TYPE(section_type), POINTER :: subsection

CPASSERT(.NOT. ASSOCIATED(section))
CALL section_create(section, __LOCATION__, name="XC_KERNEL", &
Expand All @@ -793,23 +793,28 @@ SUBROUTINE create_xc_kernel_section(section)
"Cannot be combined with XC_FUNCTIONAL or XC_POTENTIAL.", &
n_keywords=1, n_subsections=1, repeats=.FALSE.)

NULLIFY (subsection, keyword)
CALL section_create(subsection, __LOCATION__, name="UZH2022", &
description="Uses the UZH2022 xc kernel", &
n_keywords=3, n_subsections=0, repeats=.TRUE.)
CALL keyword_create(keyword, __LOCATION__, name="ALPHA", &
description="Value of the alpha parameter (default = 1.19).", &
usage="ALPHA 1.19", default_r_val=1.19_dp)
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="BETA", &
description="Value of the beta parameter (default = 0.01).", &
usage="BETA 0.01", default_r_val=0.01_dp)
CALL section_add_keyword(subsection, keyword)
NULLIFY (keyword)
CALL keyword_create( &
keyword, __LOCATION__, name="_SECTION_PARAMETERS_", &
description="Selection of kernel functionals.", &
usage="&XC_KERNEL LDAfxc", &
enum_c_vals=s2a("PADEfxc", "LDAfxc", "GGAfxc", "NONE"), &
enum_i_vals=(/fxc_funct_pade, fxc_funct_lda, fxc_funct_gga, fxc_none/), &
enum_desc=s2a("Fxc based on LDA PADE approximation", &
"Fxc based on LDA functionals", &
"Fxc model from fit to PBE functional", &
"NONE"), &
default_i_val=fxc_none, &
lone_keyword_i_val=fxc_none)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)

CALL keyword_create(keyword, __LOCATION__, name="PARAMETER", &
description="List of parameters specific to the kernel function", &
usage="PARAMETER <REAL> .. <REAL>", &
type_of_var=real_t, n_var=-1)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="SCALE_X", &
description="Scaling parameter for exchange kernel.", &
usage="SCALE_X 0.2", default_r_val=1.0_dp)
Expand Down
2 changes: 1 addition & 1 deletion src/pw/pw_methods.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1205,7 +1205,7 @@ SUBROUTINE pw_multiply(pw_out, pw1, pw2, alpha)
pw_out%in_use == REALDATA1D) THEN
IF (my_alpha == 1.0_dp) THEN
!$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw_out, pw1, pw2)
pw_out%cr = pw_out%cr + pw1%cr + pw2%cr
pw_out%cr = pw_out%cr + pw1%cr*pw2%cr
!$OMP END PARALLEL WORKSHARE
flop = REAL(2*SIZE(pw2%cr), KIND=dp)
ELSE
Expand Down
85 changes: 84 additions & 1 deletion src/qs_kernel_methods.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,8 @@ MODULE qs_kernel_methods
USE kinds, ONLY: dp
USE pw_env_types, ONLY: pw_env_get,&
pw_env_type
USE pw_methods, ONLY: pw_axpy,&
pw_zero
USE pw_pool_types, ONLY: pw_pool_type
USE pw_types, ONLY: pw_type
USE qs_environment_types, ONLY: get_qs_env,&
Expand All @@ -34,6 +36,7 @@ MODULE qs_kernel_methods
tddfpt_subgroup_env_type
USE xc, ONLY: xc_prep_2nd_deriv
USE xc_derivatives, ONLY: xc_functionals_get_needs
USE xc_fxc_kernel, ONLY: calc_fxc_kernel
USE xc_rho_set_types, ONLY: xc_rho_set_create
#include "./base/base_uses.f90"

Expand All @@ -45,7 +48,7 @@ MODULE qs_kernel_methods

LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.

PUBLIC :: create_kernel_env
PUBLIC :: create_kernel_env, create_fxc_kernel

CONTAINS

Expand Down Expand Up @@ -212,4 +215,84 @@ SUBROUTINE create_kernel_env(kernel_env, xc_section, is_rks_triplets, rho_struct

END SUBROUTINE create_kernel_env

! **************************************************************************************************
!> \brief Create the xc kernel potential for the approximate Fxc kernel model
!> \param rho_struct ...
!> \param fxc_rspace ...
!> \param xc_section ...
!> \param is_rks_triplets ...
!> \param sub_env ...
!> \param qs_env ...
! **************************************************************************************************
SUBROUTINE create_fxc_kernel(rho_struct, fxc_rspace, xc_section, is_rks_triplets, sub_env, qs_env)
TYPE(qs_rho_type), POINTER :: rho_struct
TYPE(pw_type), DIMENSION(:), POINTER :: fxc_rspace
TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
LOGICAL, INTENT(IN) :: is_rks_triplets
TYPE(tddfpt_subgroup_env_type), INTENT(IN) :: sub_env
TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env

CHARACTER(LEN=*), PARAMETER :: routineN = 'create_fxc_kernel'

INTEGER :: handle, ispin, nspins
LOGICAL :: rho_g_valid, tau_r_valid
REAL(KIND=dp) :: factor
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
TYPE(pw_type), DIMENSION(:), POINTER :: rho_g, rho_r, tau_r
TYPE(pw_type), POINTER :: rho_nlcc, rho_nlcc_g
TYPE(section_vals_type), POINTER :: xc_kernel

CALL timeset(routineN, handle)

pw_env => sub_env%pw_env
CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)

NULLIFY (rho_r, rho_g, tau_r)
CALL qs_rho_get(rho_struct, &
tau_r_valid=tau_r_valid, &
rho_g_valid=rho_g_valid, &
rho_r=rho_r, &
rho_g=rho_g, &
tau_r=tau_r)

IF (.NOT. tau_r_valid) NULLIFY (tau_r)
IF (.NOT. rho_g_valid) NULLIFY (rho_g)

nspins = SIZE(rho_r)

NULLIFY (rho_nlcc, rho_nlcc_g)!deb
CALL get_qs_env(qs_env, &
rho_nlcc=rho_nlcc, &
rho_nlcc_g=rho_nlcc_g)
! add the nlcc densities
IF (ASSOCIATED(rho_nlcc)) THEN
factor = 1.0_dp
DO ispin = 1, nspins
CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
END DO
END IF

DO ispin = 1, SIZE(fxc_rspace)
CALL pw_zero(fxc_rspace(ispin))
END DO

xc_kernel => section_vals_get_subs_vals(xc_section, "XC_KERNEL")
CALL calc_fxc_kernel(fxc_rspace, rho_r, rho_g, tau_r, &
xc_kernel, is_rks_triplets, auxbas_pw_pool)

! remove the nlcc densities (keep stuff in original state)
IF (ASSOCIATED(rho_nlcc)) THEN
factor = -1.0_dp
DO ispin = 1, nspins
CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
END DO
END IF

CALL timestop(handle)

END SUBROUTINE create_fxc_kernel

END MODULE qs_kernel_methods
1 change: 1 addition & 0 deletions src/qs_kernel_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,7 @@ MODULE qs_kernel_types
!> first and second derivatives of exchange-correlation functional
TYPE(xc_derivative_set_type) :: xc_deriv_set
!> XC input section
LOGICAL :: do_exck
TYPE(section_vals_type), POINTER :: xc_section => Null()
!> flags which indicate required components of the exchange-correlation functional
!> (density, gradient, etc)
Expand Down
32 changes: 21 additions & 11 deletions src/qs_tddfpt2_fhxc.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,8 @@ MODULE qs_tddfpt2_fhxc
USE qs_tddfpt2_densities, ONLY: tddfpt_construct_aux_fit_density
USE qs_tddfpt2_lri_utils, ONLY: tddfpt2_lri_Amat
USE qs_tddfpt2_operators, ONLY: tddfpt_apply_coulomb,&
tddfpt_apply_xc
tddfpt_apply_xc,&
tddfpt_apply_xc_potential
USE qs_tddfpt2_stda_types, ONLY: stda_env_type
USE qs_tddfpt2_stda_utils, ONLY: stda_calculate_kernel
USE qs_tddfpt2_subgroups, ONLY: tddfpt_subgroup_env_type
Expand Down Expand Up @@ -229,11 +230,15 @@ SUBROUTINE fhxc_kernel(Aop_evects, evects, is_rks_triplets, &
! C_x d^{2}E_{x}^{DFT}[\rho] / d\rho^2
! + C_{HF} d^{2}E_{x, ADMM}^{DFT}[\rho] / d\rho^2 in case of ADMM calculation
IF (gapw_xc) THEN
CALL tddfpt_apply_xc(A_ia_rspace=work_matrices%A_ia_rspace_sub, kernel_env=kernel_env, &
rho_ia_struct=work_matrices%rho_xc_struct_sub, &
is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, &
work_v_xc=work_matrices%wpw_rspace_sub, &
work_v_xc_tau=work_matrices%wpw_tau_rspace_sub)
IF (kernel_env%do_exck) THEN
CPABORT("NYA")
ELSE
CALL tddfpt_apply_xc(A_ia_rspace=work_matrices%A_ia_rspace_sub, kernel_env=kernel_env, &
rho_ia_struct=work_matrices%rho_xc_struct_sub, &
is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, &
work_v_xc=work_matrices%wpw_rspace_sub, &
work_v_xc_tau=work_matrices%wpw_tau_rspace_sub)
END IF
DO ispin = 1, nspins
CALL pw_scale(work_matrices%A_ia_rspace_sub(ispin), &
work_matrices%A_ia_rspace_sub(ispin)%pw_grid%dvol)
Expand All @@ -245,11 +250,16 @@ SUBROUTINE fhxc_kernel(Aop_evects, evects, is_rks_triplets, &
CALL pw_zero(work_matrices%A_ia_rspace_sub(ispin))
END DO
ELSE
CALL tddfpt_apply_xc(A_ia_rspace=work_matrices%A_ia_rspace_sub, kernel_env=kernel_env, &
rho_ia_struct=work_matrices%rho_orb_struct_sub, &
is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, &
work_v_xc=work_matrices%wpw_rspace_sub, &
work_v_xc_tau=work_matrices%wpw_tau_rspace_sub)
IF (kernel_env%do_exck) THEN
CALL tddfpt_apply_xc_potential(work_matrices%A_ia_rspace_sub, work_matrices%fxc_rspace_sub, &
work_matrices%rho_orb_struct_sub, is_rks_triplets)
ELSE
CALL tddfpt_apply_xc(A_ia_rspace=work_matrices%A_ia_rspace_sub, kernel_env=kernel_env, &
rho_ia_struct=work_matrices%rho_orb_struct_sub, &
is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, &
work_v_xc=work_matrices%wpw_rspace_sub, &
work_v_xc_tau=work_matrices%wpw_tau_rspace_sub)
END IF
END IF
IF (gapw .OR. gapw_xc) THEN
rho_atom_set => sub_env%local_rho_set%rho_atom_set
Expand Down
3 changes: 3 additions & 0 deletions src/qs_tddfpt2_fhxc_forces.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -445,6 +445,9 @@ SUBROUTINE fhxc_force(qs_env, ex_env, gs_mos, full_kernel, debug_forces)
END IF

! XC
IF (full_kernel%do_exck) THEN
CPABORT("NYA")
END IF
NULLIFY (fxc_rho, fxc_tau, gxc_rho, gxc_tau)
xc_section => full_kernel%xc_section
CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
Expand Down
17 changes: 13 additions & 4 deletions src/qs_tddfpt2_methods.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,8 @@ MODULE qs_tddfpt2_methods
USE qs_2nd_kernel_ao, ONLY: build_dm_response
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_kernel_methods, ONLY: create_kernel_env
USE qs_kernel_methods, ONLY: create_fxc_kernel,&
create_kernel_env
USE qs_kernel_types, ONLY: full_kernel_env_type,&
kernel_env_type,&
release_kernel_env
Expand Down Expand Up @@ -247,7 +248,8 @@ SUBROUTINE tddfpt(qs_env, calc_forces)

! allocate pools and work matrices
nstates = tddfpt_control%nstates
CALL tddfpt_create_work_matrices(work_matrices, gs_mos, nstates, do_hfx, do_admm, qs_env, sub_env)
CALL tddfpt_create_work_matrices(work_matrices, gs_mos, nstates, do_hfx, do_admm, do_exck, &
qs_env, sub_env)

CALL tddfpt_construct_ground_state_orb_density(rho_orb_struct=work_matrices%rho_orb_struct_sub, &
rho_xc_struct=work_matrices%rho_xc_struct_sub, &
Expand Down Expand Up @@ -294,6 +296,13 @@ SUBROUTINE tddfpt(qs_env, calc_forces)
kernel_env%full_kernel => full_kernel_env
NULLIFY (kernel_env%admm_kernel)
END IF
! Fxc from kernel definition
kernel_env%full_kernel%do_exck = do_exck
! initilize xc kernel
IF (do_exck) THEN
CALL create_fxc_kernel(work_matrices%rho_orb_struct_sub, work_matrices%fxc_rspace_sub, &
xc_section, tddfpt_control%rks_triplets, sub_env, qs_env)
END IF
ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
! sTDA kernel
nactive = 0
Expand Down Expand Up @@ -715,12 +724,12 @@ SUBROUTINE tddfpt_input(qs_env, do_hfx, do_admm, do_exck, do_hfxlr, &
"the one used for the ground-state calculation. A ground-state 'admm_env' cannot be reused.")
END IF
! SET HFX_KERNEL and/or XC_KERNEL
IF (exhfxk) THEN
IF (exk) THEN
do_exck = .TRUE.
ELSE
do_exck = .FALSE.
END IF
IF (exk) THEN
IF (exhfxk) THEN
do_hfxlr = .TRUE.
ELSE
do_hfxlr = .FALSE.
Expand Down

0 comments on commit 90e3ec8

Please sign in to comment.