Skip to content

Commit

Permalink
Add tau-dependent 2nd derivatives
Browse files Browse the repository at this point in the history
  • Loading branch information
@fstein93
Frederick Stein authored and fstein93 committed Mar 4, 2022
1 parent 5d3faee commit 839adb6
Show file tree
Hide file tree
Showing 37 changed files with 1,911 additions and 451 deletions.
13 changes: 5 additions & 8 deletions src/cp_control_utils.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,9 +62,9 @@ MODULE cp_control_utils
USE qs_cdft_utils, ONLY: read_cdft_control_section
USE string_utilities, ONLY: uppercase
USE util, ONLY: sort
USE xc_derivatives, ONLY: xc_functionals_get_needs
USE xc, ONLY: xc_uses_kinetic_energy_density,&
xc_uses_norm_drho
USE xc_input_constants, ONLY: xc_deriv_collocate
USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
USE xc_write_output, ONLY: xc_write
#include "./base/base_uses.f90"

Expand Down Expand Up @@ -106,7 +106,6 @@ SUBROUTINE read_dft_control(dft_control, dft_section)
TYPE(section_vals_type), POINTER :: maxwell_section, sccs_section, &
scf_section, tmp_section, &
xc_fun_section, xc_section
TYPE(xc_rho_cflags_type) :: needs

was_present = .FALSE.

Expand Down Expand Up @@ -168,13 +167,11 @@ SUBROUTINE read_dft_control(dft_control, dft_section)
END IF

dft_control%lsd = (dft_control%nspins > 1)
needs = xc_functionals_get_needs(xc_fun_section, &
lsd=dft_control%lsd, &
add_basic_components=.TRUE.)
dft_control%use_kinetic_energy_density = (needs%tau_spin .OR. needs%tau)
dft_control%use_kinetic_energy_density = xc_uses_kinetic_energy_density(xc_fun_section, dft_control%lsd)

xc_deriv_method_id = section_get_ival(xc_section, "XC_GRID%XC_DERIV")
dft_control%drho_by_collocation = (needs%norm_drho .AND. (xc_deriv_method_id == xc_deriv_collocate))
dft_control%drho_by_collocation = (xc_uses_norm_drho(xc_fun_section, dft_control%lsd) &
.AND. (xc_deriv_method_id == xc_deriv_collocate))
IF (dft_control%drho_by_collocation) THEN
CPABORT("derivatives by collocation not implemented")
END IF
Expand Down
2 changes: 1 addition & 1 deletion src/ec_env_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,7 @@ MODULE ec_env_types
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: cpmos
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz, matrix_z, matrix_wz, z_admm
! Harris (rhoout), and response density (rhoz) on grid
TYPE(pw_p_type), DIMENSION(:), POINTER :: rhoout_r, rhoz_r
TYPE(pw_p_type), DIMENSION(:), POINTER :: rhoout_r, rhoz_r, tauout_r
! potentials from input density
TYPE(pw_p_type), POINTER :: vh_rspace
TYPE(pw_p_type), DIMENSION(:), POINTER :: vxc_rspace, vtau_rspace, vadmm_rspace
Expand Down
19 changes: 17 additions & 2 deletions src/ec_environment.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,8 @@ MODULE ec_environment
ls_s_preconditioner_atomic, ls_s_preconditioner_molecular, ls_s_preconditioner_none, &
ls_s_sqrt_ns, ls_s_sqrt_proot, xc_vdw_fun_nonloc, xc_vdw_fun_pairpot
USE input_cp2k_check, ONLY: xc_functionals_expand
USE input_section_types, ONLY: section_vals_get,&
USE input_section_types, ONLY: section_get_ival,&
section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
Expand All @@ -55,7 +56,11 @@ MODULE ec_environment
USE qs_kind_types, ONLY: get_qs_kind,&
get_qs_kind_set,&
qs_kind_type
USE qs_rho_types, ONLY: qs_rho_type
USE string_utilities, ONLY: uppercase
USE xc, ONLY: xc_uses_kinetic_energy_density,&
xc_uses_norm_drho
USE xc_input_constants, ONLY: xc_deriv_collocate
#include "./base/base_uses.f90"

IMPLICIT NONE
Expand Down Expand Up @@ -120,8 +125,9 @@ SUBROUTINE init_ec_env(qs_env, ec_env, dft_section, ec_section)
TYPE(qs_dispersion_type), POINTER :: dispersion_env
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_kind_type), POINTER :: qs_kind
TYPE(qs_rho_type), POINTER :: rho
TYPE(section_vals_type), POINTER :: nl_section, pp_section, section1, &
section2, xc_section
section2, xc_fun_section, xc_section

CALL timeset(routineN, handle)

Expand Down Expand Up @@ -259,6 +265,15 @@ SUBROUTINE init_ec_env(qs_env, ec_env, dft_section, ec_section)
ELSE
ec_env%xc_section => xc_section
END IF
! Check whether energy correction requires the kinetic energy density and rebuild rho if necessary
CALL get_qs_env(qs_env, dft_control=dft_control, rho=rho)
xc_fun_section => section_vals_get_subs_vals(ec_env%xc_section, "XC_FUNCTIONAL")
dft_control%use_kinetic_energy_density = dft_control%use_kinetic_energy_density .OR. &
xc_uses_kinetic_energy_density(xc_fun_section, dft_control%lsd)
! Same for density gradient
dft_control%drho_by_collocation = dft_control%drho_by_collocation .OR. &
(xc_uses_norm_drho(xc_fun_section, dft_control%lsd) .AND. &
(section_get_ival(xc_section, "XC_GRID%XC_DERIV") == xc_deriv_collocate))
! dispersion
ALLOCATE (dispersion_env)
NULLIFY (xc_section)
Expand Down
20 changes: 12 additions & 8 deletions src/ec_methods.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,22 +71,23 @@ MODULE ec_methods
!> \param vxc will contain the partially integrated second derivative
!> taken with respect to rho, evaluated in rho and folded with rho1
!> vxc is allocated here and needs to be deallocated by the caller.
!> \param vxc_tau ...
!> \param rho density at which derivatives were calculated
!> \param rho1_r density in r-space with which to fold
!> \param rho1_g density in g-space with which to fold
!> \param tau1_r ...
!> \param xc_section XC parameters of this potential
!> \param compute_virial Enable stress-tensor calculation
!> \param virial_xc Will contain GGA contribution of XC-kernel to stress-tensor
!> \date 11.2019
!> \author fbelle
! **************************************************************************************************
SUBROUTINE create_kernel(qs_env, vxc, rho, rho1_r, rho1_g, xc_section, compute_virial, virial_xc)
SUBROUTINE create_kernel(qs_env, vxc, vxc_tau, rho, rho1_r, rho1_g, tau1_r, xc_section, compute_virial, virial_xc)

TYPE(qs_environment_type), POINTER :: qs_env
TYPE(pw_p_type), DIMENSION(:), INTENT(OUT), &
POINTER :: vxc
TYPE(pw_p_type), DIMENSION(:), POINTER :: vxc, vxc_tau
TYPE(qs_rho_type), INTENT(IN), POINTER :: rho
TYPE(pw_p_type), DIMENSION(:), INTENT(IN), POINTER :: rho1_r, rho1_g
TYPE(pw_p_type), DIMENSION(:), INTENT(IN), POINTER :: rho1_r, rho1_g, tau1_r
TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
LOGICAL, INTENT(IN), OPTIONAL :: compute_virial
REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
Expand All @@ -96,32 +97,35 @@ SUBROUTINE create_kernel(qs_env, vxc, rho, rho1_r, rho1_g, xc_section, compute_v

INTEGER :: handle
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r
TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r, tau_r
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
TYPE(xc_derivative_set_type), POINTER :: deriv_set
TYPE(xc_rho_set_type), POINTER :: rho_set

CALL timeset(routineN, handle)

NULLIFY (auxbas_pw_pool, deriv_set, pw_env, rho_r, rho_set)
NULLIFY (auxbas_pw_pool, deriv_set, pw_env, rho_r, rho_set, vxc, vxc_tau)

CALL get_qs_env(qs_env, pw_env=pw_env)
CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool)
! Get density
CALL qs_rho_get(rho, rho_r=rho_r)
CALL qs_rho_get(rho, rho_r=rho_r, tau_r=tau_r)

CALL xc_prep_2nd_deriv(deriv_set=deriv_set, & ! containing potentials
rho_set=rho_set, & ! density at which derivs are calculated
rho_r=rho_r, & ! place where derivative is evaluated
tau_r=tau_r, &
pw_pool=auxbas_pw_pool, & ! pool for grids
xc_section=xc_section)

! folding of second deriv with density in rho1_set
CALL xc_calc_2nd_deriv(v_xc=vxc, & ! XC-potential
CALL xc_calc_2nd_deriv(v_xc=vxc, & ! XC-potential, rho-dependent
v_xc_tau=vxc_tau, & ! XC-potential, tau-dependent
deriv_set=deriv_set, & ! deriv of xc-potential
rho_set=rho_set, & ! density at which deriv are calculated
rho1_r=rho1_r, & ! density with which to fold
rho1_g=rho1_g, & ! density with which to fold
tau1_r=tau1_r, &
pw_pool=auxbas_pw_pool, & ! pool for grids
xc_section=xc_section, &
gapw=.FALSE., &
Expand Down
29 changes: 25 additions & 4 deletions src/ec_orth_solver.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1133,7 +1133,7 @@ SUBROUTINE hessian_op2(qs_env, matrix_Ax, matrix_p, matrix_s_sqrt_inv, rho1, eps
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_p_type) :: rho_tot_gspace, v_hartree_gspace, &
v_hartree_rspace
TYPE(pw_p_type), DIMENSION(:), POINTER :: rho1_g, rho1_r, v_xc
TYPE(pw_p_type), DIMENSION(:), POINTER :: rho1_g, rho1_r, tau1_r, v_xc, v_xc_tau
TYPE(pw_poisson_type), POINTER :: poisson_env
TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
Expand All @@ -1155,7 +1155,7 @@ SUBROUTINE hessian_op2(qs_env, matrix_Ax, matrix_p, matrix_s_sqrt_inv, rho1, eps
! Get non-ortho input density matrix on grid
CALL qs_rho_get(rho, rho_ao=rho_ao)
! Get non-ortho trial density
CALL qs_rho_get(rho1, rho_g=rho1_g, rho_r=rho1_r)
CALL qs_rho_get(rho1, rho_g=rho1_g, rho_r=rho1_r, tau_r=tau1_r)

NULLIFY (pw_env)
CALL get_qs_env(qs_env, pw_env=pw_env)
Expand Down Expand Up @@ -1183,7 +1183,7 @@ SUBROUTINE hessian_op2(qs_env, matrix_Ax, matrix_p, matrix_s_sqrt_inv, rho1, eps
in_space=REALSPACE)

! XC-Kernel
NULLIFY (v_xc, xc_section)
NULLIFY (v_xc, v_xc_tau, xc_section)
xc_section => section_vals_get_subs_vals(input, "DFT%XC")

! No HFX allowed
Expand All @@ -1198,9 +1198,11 @@ SUBROUTINE hessian_op2(qs_env, matrix_Ax, matrix_p, matrix_s_sqrt_inv, rho1, eps
! add xc-kernel
CALL create_kernel(qs_env, &
vxc=v_xc, &
vxc_tau=v_xc_tau, &
rho=rho, &
rho1_r=rho1_r, &
rho1_g=rho1_g, &
tau1_r=tau1_r, &
xc_section=xc_section)

DO ispin = 1, nspins
Expand All @@ -1223,7 +1225,10 @@ SUBROUTINE hessian_op2(qs_env, matrix_Ax, matrix_p, matrix_s_sqrt_inv, rho1, eps
DO ispin = 1, nspins
CALL pw_axpy(v_hartree_rspace%pw, v_xc(ispin)%pw)
END DO
IF (nspins == 1) CALL pw_scale(v_xc(1)%pw, 2.0_dp)
IF (nspins == 1) THEN
CALL pw_scale(v_xc(1)%pw, 2.0_dp)
IF (ASSOCIATED(v_xc_tau)) CALL pw_scale(v_xc_tau(1)%pw, 2.0_dp)
END IF

! Init response kernel matrix
! matrix G(B)
Expand All @@ -1241,6 +1246,16 @@ SUBROUTINE hessian_op2(qs_env, matrix_Ax, matrix_p, matrix_s_sqrt_inv, rho1, eps
qs_env=qs_env, &
calculate_forces=.FALSE., &
basis_type="ORB")

IF (ASSOCIATED(v_xc_tau)) THEN
CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
pmat=rho_ao(ispin), &
hmat=matrix_G(ispin), &
qs_env=qs_env, &
compute_tau=.TRUE., &
calculate_forces=.FALSE., &
basis_type="ORB")
END IF
END DO

! Calculate KG correction to kernel
Expand Down Expand Up @@ -1278,6 +1293,12 @@ SUBROUTINE hessian_op2(qs_env, matrix_Ax, matrix_p, matrix_s_sqrt_inv, rho1, eps
CALL pw_pool_give_back_pw(auxbas_pw_pool, v_xc(ispin)%pw)
END DO
DEALLOCATE (v_xc)
IF (ASSOCIATED(v_xc_tau)) THEN
DO ispin = 1, nspins
CALL pw_pool_give_back_pw(auxbas_pw_pool, v_xc_tau(ispin)%pw)
END DO
DEALLOCATE (v_xc_tau)
END IF

CALL dbcsr_deallocate_matrix_set(matrix_G)

Expand Down

0 comments on commit 839adb6

Please sign in to comment.