Skip to content

Commit

Permalink
2nd/3rd functional derivatives for TDDFT forces: (#3122)
Browse files Browse the repository at this point in the history
  • Loading branch information
@juerghutter
juerghutter committed Nov 21, 2023
1 parent 39fc2fe commit c3608e6
Show file tree
Hide file tree
Showing 9 changed files with 602 additions and 20 deletions.
6 changes: 6 additions & 0 deletions src/exstates_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,8 @@ MODULE exstates_types
INTEGER :: state
REAL(KIND=dp) :: evalue
INTEGER :: xc_kernel_method
REAL(KIND=dp) :: eps_delta_rho = 1.E-02_dp
INTEGER :: diff_order
LOGICAL :: debug_forces = .FALSE.
TYPE(cp_fm_type), POINTER, DIMENSION(:) :: evect => NULL()
TYPE(cp_fm_type), POINTER, DIMENSION(:) :: cpmos => NULL()
Expand Down Expand Up @@ -200,6 +202,10 @@ SUBROUTINE exstate_create(ex_env, excited_state, dft_section)
i_val=ex_env%xc_kernel_method)
CALL section_vals_val_get(dft_section, "EXCITED_STATES%DEBUG_FORCES", &
l_val=ex_env%debug_forces)
CALL section_vals_val_get(dft_section, "EXCITED_STATES%EPS_DELTA_RHO", &
r_val=ex_env%eps_delta_rho)
CALL section_vals_val_get(dft_section, "EXCITED_STATES%DIFF_ORDER", &
i_val=ex_env%diff_order)
ELSE
ex_env%state = 0
ex_env%xc_kernel_method = xc_kernel_method_best
Expand Down
15 changes: 15 additions & 0 deletions src/input_cp2k_exstate.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@ MODULE input_cp2k_exstate
USE input_section_types, ONLY: section_add_keyword,&
section_create,&
section_type
USE kinds, ONLY: dp
USE string_utilities, ONLY: s2a
#include "./base/base_uses.f90"

Expand Down Expand Up @@ -74,6 +75,20 @@ SUBROUTINE create_exstate_section(section)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="EPS_DELTA_RHO", &
description="Step size for finite difference calculation of functional derivatives.", &
usage="EPS_DELTA_RHO 1.E-02", &
default_r_val=1.E-03_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="DIFF_ORDER", &
description="Order of finite differentiation formula used for functional derivatives.", &
usage="DIFF_ORDER 4", &
default_i_val=6)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="OVERLAP_DELTAT", &
description="Keyword for the computation of the overlap matrix between two consecutive time steps.", &
usage="OVERLAP_DELTAT", &
Expand Down
136 changes: 135 additions & 1 deletion src/qs_fxc.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ MODULE qs_fxc
PRIVATE

! *** Public subroutines ***
PUBLIC :: qs_fxc_fdiff, qs_fxc_analytic, qs_fgxc_create, qs_fgxc_release
PUBLIC :: qs_fxc_fdiff, qs_fxc_analytic, qs_fgxc_gdiff, qs_fgxc_create, qs_fgxc_release

CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_fxc'

Expand Down Expand Up @@ -267,6 +267,140 @@ SUBROUTINE qs_fxc_fdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy,

END SUBROUTINE qs_fxc_fdiff

! **************************************************************************************************
!> \brief ...
!> \param ks_env ...
!> \param rho0_struct ...
!> \param rho1_struct ...
!> \param xc_section ...
!> \param accuracy ...
!> \param epsrho ...
!> \param is_triplet ...
!> \param fxc_rho ...
!> \param fxc_tau ...
!> \param gxc_rho ...
!> \param gxc_tau ...
! **************************************************************************************************
SUBROUTINE qs_fgxc_gdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, epsrho, &
is_triplet, fxc_rho, fxc_tau, gxc_rho, gxc_tau)

TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
TYPE(section_vals_type), POINTER :: xc_section
INTEGER, INTENT(IN) :: accuracy
REAL(KIND=dp), INTENT(IN) :: epsrho
LOGICAL, INTENT(IN) :: is_triplet
TYPE(pw_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau

CHARACTER(len=*), PARAMETER :: routineN = 'qs_fgxc_gdiff'

INTEGER :: handle, ispin, istep, nspins, nstep
REAL(KIND=dp) :: alpha, beta, exc, oeps1
REAL(KIND=dp), DIMENSION(-4:4) :: ak
TYPE(dft_control_type), POINTER :: dft_control
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
TYPE(pw_type), DIMENSION(:), POINTER :: v_tau_rspace, vxc00
TYPE(qs_rho_type), POINTER :: rhoin

CALL timeset(routineN, handle)

CPASSERT(.NOT. ASSOCIATED(fxc_rho))
CPASSERT(.NOT. ASSOCIATED(fxc_tau))
CPASSERT(.NOT. ASSOCIATED(gxc_rho))
CPASSERT(.NOT. ASSOCIATED(gxc_tau))
CPASSERT(ASSOCIATED(rho0_struct))
CPASSERT(ASSOCIATED(rho1_struct))

ak = 0.0_dp
SELECT CASE (accuracy)
CASE (:4)
nstep = 2
ak(-2:2) = (/1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp/)/12.0_dp
CASE (5:7)
nstep = 3
ak(-3:3) = (/-1.0_dp, 9.0_dp, -45.0_dp, 0.0_dp, 45.0_dp, -9.0_dp, 1.0_dp/)/60.0_dp
CASE (8:)
nstep = 4
ak(-4:4) = (/1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp/)/280.0_dp
END SELECT

CALL get_ks_env(ks_env, dft_control=dft_control, pw_env=pw_env)
CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)

nspins = dft_control%nspins
exc = 0.0_dp

CALL qs_fxc_fdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, &
fxc_rho, fxc_tau)

DO istep = -nstep, nstep

IF (ak(istep) /= 0.0_dp) THEN
alpha = 1.0_dp
beta = REAL(istep, KIND=dp)*epsrho
NULLIFY (rhoin)
ALLOCATE (rhoin)
CALL qs_rho_create(rhoin)
NULLIFY (vxc00, v_tau_rspace)
CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
CALL qs_fxc_fdiff(ks_env=ks_env, rho0_struct=rhoin, rho1_struct=rho1_struct, &
xc_section=xc_section, accuracy=accuracy, is_triplet=is_triplet, &
fxc_rho=vxc00, fxc_tau=v_tau_rspace)
CALL qs_rho_release(rhoin)
DEALLOCATE (rhoin)
IF (.NOT. ASSOCIATED(gxc_rho)) THEN
ALLOCATE (gxc_rho(nspins))
DO ispin = 1, nspins
CALL auxbas_pw_pool%create_pw(gxc_rho(ispin), &
in_space=REALSPACE, use_data=REALDATA3D)
CALL pw_zero(gxc_rho(ispin))
END DO
END IF
DO ispin = 1, nspins
CALL pw_axpy(vxc00(ispin), gxc_rho(ispin), ak(istep))
END DO
DO ispin = 1, SIZE(vxc00)
CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
END DO
DEALLOCATE (vxc00)
IF (ASSOCIATED(v_tau_rspace)) THEN
IF (.NOT. ASSOCIATED(gxc_tau)) THEN
ALLOCATE (gxc_tau(nspins))
DO ispin = 1, nspins
CALL auxbas_pw_pool%create_pw(gxc_tau(ispin), &
in_space=REALSPACE, use_data=REALDATA3D)
CALL pw_zero(gxc_tau(ispin))
END DO
END IF
DO ispin = 1, nspins
CALL pw_axpy(v_tau_rspace(ispin), gxc_tau(ispin), ak(istep))
END DO
DO ispin = 1, SIZE(v_tau_rspace)
CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
END DO
DEALLOCATE (v_tau_rspace)
END IF
END IF

END DO

oeps1 = 1.0_dp/epsrho
DO ispin = 1, nspins
CALL pw_scale(gxc_rho(ispin), oeps1)
END DO
IF (ASSOCIATED(gxc_tau)) THEN
DO ispin = 1, nspins
CALL pw_scale(gxc_tau(ispin), oeps1)
END DO
END IF

CALL timestop(handle)

END SUBROUTINE qs_fgxc_gdiff

! **************************************************************************************************
!> \brief ...
!> \param ks_env ...
Expand Down
74 changes: 60 additions & 14 deletions src/qs_tddfpt2_fhxc_forces.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,7 @@ MODULE qs_tddfpt2_fhxc_forces
set_qs_env
USE qs_force_types, ONLY: qs_force_type
USE qs_fxc, ONLY: qs_fgxc_create,&
qs_fgxc_gdiff,&
qs_fgxc_release
USE qs_gapw_densities, ONLY: prepare_gapw_den
USE qs_integrate_potential, ONLY: integrate_v_rspace
Expand Down Expand Up @@ -130,7 +131,8 @@ MODULE qs_tddfpt2_fhxc_forces
USE qs_tddfpt2_subgroups, ONLY: tddfpt_subgroup_env_type
USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos,&
tddfpt_work_matrices
USE qs_vxc_atom, ONLY: calculate_gfxc_atom
USE qs_vxc_atom, ONLY: calculate_gfxc_atom,&
gfxc_atom_diff
USE task_list_types, ONLY: task_list_type
USE util, ONLY: get_limit
USE virial_types, ONLY: virial_type
Expand Down Expand Up @@ -172,11 +174,11 @@ SUBROUTINE fhxc_force(qs_env, ex_env, gs_mos, full_kernel, debug_forces)
CHARACTER(LEN=default_string_length) :: basis_type
INTEGER :: handle, iounit, ispin, mspin, myfun, &
n_rep_hf, nao, nao_aux, natom, nkind, &
norb, nspins
norb, nspins, order
LOGICAL :: distribute_fock_matrix, do_admm, do_analytic, do_hfx, do_numeric, gapw, gapw_xc, &
hfx_treat_lsd_in_core, is_rks_triplets, s_mstruct_changed, use_virial
REAL(KIND=dp) :: eh1, eh1c, eps_fit, focc, fscal, fval, &
kval, xehartree
REAL(KIND=dp) :: eh1, eh1c, eps_delta, eps_fit, focc, &
fscal, fval, kval, xehartree
REAL(KIND=dp), DIMENSION(3) :: fodeb
TYPE(admm_type), POINTER :: admm_env
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
Expand Down Expand Up @@ -510,6 +512,8 @@ SUBROUTINE fhxc_force(qs_env, ex_env, gs_mos, full_kernel, debug_forces)
CASE DEFAULT
CPABORT("invalid xc_kernel_method")
END SELECT
order = ex_env%diff_order
eps_delta = ex_env%eps_delta_rho

IF (gapw_xc) THEN
CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, rho_xc=rho)
Expand All @@ -529,9 +533,16 @@ SUBROUTINE fhxc_force(qs_env, ex_env, gs_mos, full_kernel, debug_forces)
END IF
IF (do_analytic .AND. .NOT. do_numeric) THEN
CPABORT("Analytic 3rd EXC derivatives not available")
ELSEIF (do_numeric) THEN
IF (do_analytic) THEN
CALL qs_fgxc_gdiff(ks_env, rho, rhox, xc_section, order, eps_delta, is_rks_triplets, &
fxc_rho, fxc_tau, gxc_rho, gxc_tau)
ELSE
CALL qs_fgxc_create(ks_env, rho, rhox, xc_section, order, is_rks_triplets, &
fxc_rho, fxc_tau, gxc_rho, gxc_tau)
END IF
ELSE
CALL qs_fgxc_create(ks_env, rho, rhox, xc_section, 6, is_rks_triplets, &
fxc_rho, fxc_tau, gxc_rho, gxc_tau)
CPABORT("FHXC forces analytic/numeric")
END IF

! Well, this is a hack :-(
Expand All @@ -547,9 +558,21 @@ SUBROUTINE fhxc_force(qs_env, ex_env, gs_mos, full_kernel, debug_forces)
END DO
END IF
IF (gapw .OR. gapw_xc) THEN
CALL calculate_gfxc_atom(qs_env, ex_env%local_rho_set%rho_atom_set, &
local_rho_set_f%rho_atom_set, local_rho_set_g%rho_atom_set, &
qs_kind_set, xc_section, is_rks_triplets, 6)
IF (do_analytic .AND. .NOT. do_numeric) THEN
CPABORT("Analytic 3rd EXC derivatives not available")
ELSEIF (do_numeric) THEN
IF (do_analytic) THEN
CALL gfxc_atom_diff(qs_env, ex_env%local_rho_set%rho_atom_set, &
local_rho_set_f%rho_atom_set, local_rho_set_g%rho_atom_set, &
qs_kind_set, xc_section, is_rks_triplets, order, eps_delta)
ELSE
CALL calculate_gfxc_atom(qs_env, ex_env%local_rho_set%rho_atom_set, &
local_rho_set_f%rho_atom_set, local_rho_set_g%rho_atom_set, &
qs_kind_set, xc_section, is_rks_triplets, order)
END IF
ELSE
CPABORT("FHXC forces analytic/numeric")
END IF
END IF
IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
Expand Down Expand Up @@ -679,9 +702,16 @@ SUBROUTINE fhxc_force(qs_env, ex_env, gs_mos, full_kernel, debug_forces)
rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
IF (do_analytic .AND. .NOT. do_numeric) THEN
CPABORT("Analytic 3rd derivatives of EXC not available")
ELSEIF (do_numeric) THEN
IF (do_analytic) THEN
CALL qs_fgxc_gdiff(ks_env, rho_aux_fit, rhox_aux, xc_section, order, eps_delta, &
is_rks_triplets, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
ELSE
CALL qs_fgxc_create(ks_env, rho_aux_fit, rhox_aux, xc_section, &
order, is_rks_triplets, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
END IF
ELSE
CALL qs_fgxc_create(ks_env, rho_aux_fit, rhox_aux, xc_section, &
6, is_rks_triplets, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
CPABORT("FHXC forces analytic/numeric")
END IF
! Well, this is a hack :-(
Expand Down Expand Up @@ -740,9 +770,25 @@ SUBROUTINE fhxc_force(qs_env, ex_env, gs_mos, full_kernel, debug_forces)
rho_atom_set => admm_env%admm_gapw_env%local_rho_set%rho_atom_set
rho_atom_set_f => local_rho_set_f_admm%rho_atom_set
rho_atom_set_g => local_rho_set_g_admm%rho_atom_set
CALL calculate_gfxc_atom(qs_env, rho_atom_set, &
rho_atom_set_f, rho_atom_set_g, &
admm_env%admm_gapw_env%admm_kind_set, xc_section, is_rks_triplets, 6)

IF (do_analytic .AND. .NOT. do_numeric) THEN
CPABORT("Analytic 3rd EXC derivatives not available")
ELSEIF (do_numeric) THEN
IF (do_analytic) THEN
CALL gfxc_atom_diff(qs_env, rho_atom_set, &
rho_atom_set_f, rho_atom_set_g, &
admm_env%admm_gapw_env%admm_kind_set, xc_section, &
is_rks_triplets, order, eps_delta)
ELSE
CALL calculate_gfxc_atom(qs_env, rho_atom_set, &
rho_atom_set_f, rho_atom_set_g, &
admm_env%admm_gapw_env%admm_kind_set, xc_section, &
is_rks_triplets, order)
END IF
ELSE
CPABORT("FHXC forces analytic/numeric")
END IF

IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
IF (nspins == 1) THEN
CALL update_ks_atom(qs_env, mfx, matrix_px1_admm, forces=.TRUE., tddft=.TRUE., &
Expand Down

0 comments on commit c3608e6

Please sign in to comment.