mirror of
https://github.com/cp2k/cp2k.git
synced 2026-07-28 14:15:19 -04:00
2nd/3rd functional derivatives for TDDFT forces: (#3122)
This commit is contained in:
Juerg Hutter
GitHub
parent
39fc2feded
commit
c3608e6d9b
9 changed files with 602 additions and 20 deletions
|
|
@ -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()
|
||||
|
|
@ -200,6 +202,10 @@ CONTAINS
|
|||
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
|
||||
|
|
|
|||
|
|
@ -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"
|
||||
|
||||
|
|
@ -74,6 +75,20 @@ CONTAINS
|
|||
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", &
|
||||
|
|
|
|||
136
src/qs_fxc.F
136
src/qs_fxc.F
|
|
@ -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'
|
||||
|
||||
|
|
@ -267,6 +267,140 @@ CONTAINS
|
|||
|
||||
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 ...
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
@ -172,11 +174,11 @@ CONTAINS
|
|||
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
|
||||
|
|
@ -510,6 +512,8 @@ CONTAINS
|
|||
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)
|
||||
|
|
@ -529,9 +533,16 @@ CONTAINS
|
|||
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 :-(
|
||||
|
|
@ -547,9 +558,21 @@ CONTAINS
|
|||
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)
|
||||
|
|
@ -679,9 +702,16 @@ CONTAINS
|
|||
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 :-(
|
||||
|
|
@ -740,9 +770,25 @@ CONTAINS
|
|||
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., &
|
||||
|
|
|
|||
|
|
@ -65,6 +65,7 @@ MODULE qs_vxc_atom
|
|||
calculate_xc_2nd_deriv_atom, &
|
||||
calc_rho_angular, &
|
||||
calculate_gfxc_atom, &
|
||||
gfxc_atom_diff, &
|
||||
gaVxcgb_noGC
|
||||
|
||||
CONTAINS
|
||||
|
|
@ -1137,6 +1138,386 @@ CONTAINS
|
|||
|
||||
END SUBROUTINE calculate_gfxc_atom
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param qs_env ...
|
||||
!> \param rho0_atom_set ...
|
||||
!> \param rho1_atom_set ...
|
||||
!> \param rho2_atom_set ...
|
||||
!> \param kind_set ...
|
||||
!> \param xc_section ...
|
||||
!> \param is_triplet ...
|
||||
!> \param accuracy ...
|
||||
!> \param epsrho ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE gfxc_atom_diff(qs_env, rho0_atom_set, rho1_atom_set, rho2_atom_set, &
|
||||
kind_set, xc_section, is_triplet, accuracy, epsrho)
|
||||
|
||||
TYPE(qs_environment_type), POINTER :: qs_env
|
||||
TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho1_atom_set, &
|
||||
rho2_atom_set
|
||||
TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set
|
||||
TYPE(section_vals_type), OPTIONAL, POINTER :: xc_section
|
||||
LOGICAL, INTENT(IN) :: is_triplet
|
||||
INTEGER, INTENT(IN) :: accuracy
|
||||
REAL(KIND=dp), INTENT(IN) :: epsrho
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'gfxc_atom_diff'
|
||||
|
||||
INTEGER :: bo(2), handle, iat, iatom, ikind, ir, &
|
||||
istep, mspins, myfun, na, natom, nf, &
|
||||
nr, ns, nspins, nstep, num_pe
|
||||
INTEGER, DIMENSION(2, 3) :: bounds
|
||||
INTEGER, DIMENSION(:), POINTER :: atom_list
|
||||
LOGICAL :: donlcc, gradient_f, lsd, nlcc, paw_atom, &
|
||||
tau_f
|
||||
REAL(dp) :: beta, density_cut, gradient_cut, oeps1, &
|
||||
tau_cut
|
||||
REAL(dp), CONTIGUOUS, DIMENSION(:, :), POINTER :: rho_nlcc, weight
|
||||
REAL(dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: rho0_h, rho0_s, rho1_h, rho1_s, rho_h, &
|
||||
rho_s, tau0_h, tau0_s, tau1_h, tau1_s, &
|
||||
tau_h, tau_s, vtau_h, vtau_s, vxc_h, &
|
||||
vxc_s
|
||||
REAL(dp), CONTIGUOUS, DIMENSION(:, :, :, :), &
|
||||
POINTER :: drho0_h, drho0_s, drho1_h, drho1_s, &
|
||||
drho_h, drho_s, vxg_h, vxg_s
|
||||
REAL(dp), DIMENSION(1, 1, 1) :: tau_d
|
||||
REAL(dp), DIMENSION(1, 1, 1, 1) :: rho_d
|
||||
REAL(KIND=dp), DIMENSION(-4:4) :: ak
|
||||
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
|
||||
TYPE(dft_control_type), POINTER :: dft_control
|
||||
TYPE(grid_atom_type), POINTER :: grid_atom
|
||||
TYPE(gto_basis_set_type), POINTER :: basis_1c
|
||||
TYPE(harmonics_atom_type), POINTER :: harmonics
|
||||
TYPE(mp_para_env_type), POINTER :: para_env
|
||||
TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, fint_hh, fint_ss, int_hh, &
|
||||
int_ss, r_h, r_s
|
||||
TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
|
||||
TYPE(rho_atom_type), POINTER :: rho0_atom, rho1_atom, rho2_atom
|
||||
TYPE(section_vals_type), POINTER :: xc_fun_section
|
||||
TYPE(xc_derivative_set_type) :: deriv_set
|
||||
TYPE(xc_rho_cflags_type) :: needs
|
||||
TYPE(xc_rho_set_type) :: rho1_set_h, rho1_set_s, rho_set_h, &
|
||||
rho_set_s
|
||||
|
||||
CALL timeset(routineN, handle)
|
||||
|
||||
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
|
||||
oeps1 = 1.0_dp/epsrho
|
||||
|
||||
CALL get_qs_env(qs_env=qs_env, &
|
||||
dft_control=dft_control, &
|
||||
para_env=para_env, &
|
||||
atomic_kind_set=atomic_kind_set)
|
||||
|
||||
xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
|
||||
CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
|
||||
|
||||
IF (myfun == xc_none) THEN
|
||||
! no action needed?
|
||||
ELSE
|
||||
! calculate fxc
|
||||
CALL calculate_xc_2nd_deriv_atom(rho0_atom_set, rho1_atom_set, qs_env, xc_section, para_env, &
|
||||
do_triplet=is_triplet, kind_set_external=kind_set)
|
||||
|
||||
CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
|
||||
CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
|
||||
CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
|
||||
|
||||
nlcc = has_nlcc(kind_set)
|
||||
lsd = dft_control%lsd
|
||||
nspins = dft_control%nspins
|
||||
mspins = nspins
|
||||
IF (is_triplet) THEN
|
||||
CPASSERT(nspins == 1)
|
||||
lsd = .TRUE.
|
||||
mspins = 2
|
||||
END IF
|
||||
needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, calc_potential=.TRUE.)
|
||||
gradient_f = (needs%drho .OR. needs%drho_spin)
|
||||
tau_f = (needs%tau .OR. needs%tau_spin)
|
||||
|
||||
! Here starts the loop over all the atoms
|
||||
DO ikind = 1, SIZE(atomic_kind_set)
|
||||
CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
|
||||
CALL get_qs_kind(kind_set(ikind), paw_atom=paw_atom, &
|
||||
harmonics=harmonics, grid_atom=grid_atom)
|
||||
CALL get_qs_kind(kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
|
||||
|
||||
IF (.NOT. paw_atom) CYCLE
|
||||
|
||||
nr = grid_atom%nr
|
||||
na = grid_atom%ng_sphere
|
||||
|
||||
! Prepare the structures needed to calculate and store the xc derivatives
|
||||
|
||||
! Array dimension: here anly one dimensional arrays are used,
|
||||
! i.e. only the first column of deriv_data is read.
|
||||
! The other to dimensions are set to size equal 1
|
||||
bounds(1:2, 1:3) = 1
|
||||
bounds(2, 1) = na
|
||||
bounds(2, 2) = nr
|
||||
|
||||
! create a place where to put the derivatives
|
||||
CALL xc_dset_create(deriv_set, local_bounds=bounds)
|
||||
! create the place where to store the argument for the functionals
|
||||
CALL xc_rho_set_create(rho_set_h, bounds, rho_cutoff=density_cut, &
|
||||
drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
|
||||
CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
|
||||
drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
|
||||
CALL xc_rho_set_create(rho1_set_h, bounds, rho_cutoff=density_cut, &
|
||||
drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
|
||||
CALL xc_rho_set_create(rho1_set_s, bounds, rho_cutoff=density_cut, &
|
||||
drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
|
||||
|
||||
! allocate the required 3d arrays where to store rho and drho
|
||||
CALL xc_rho_set_atom_update(rho_set_h, needs, mspins, bounds)
|
||||
CALL xc_rho_set_atom_update(rho_set_s, needs, mspins, bounds)
|
||||
CALL xc_rho_set_atom_update(rho1_set_h, needs, mspins, bounds)
|
||||
CALL xc_rho_set_atom_update(rho1_set_s, needs, mspins, bounds)
|
||||
|
||||
weight => grid_atom%weight
|
||||
|
||||
ALLOCATE (rho_h(na, nr, nspins), rho_s(na, nr, nspins), &
|
||||
rho0_h(na, nr, nspins), rho0_s(na, nr, nspins), &
|
||||
rho1_h(na, nr, nspins), rho1_s(na, nr, nspins))
|
||||
ALLOCATE (vxc_h(na, nr, nspins), vxc_s(na, nr, nspins))
|
||||
IF (gradient_f) THEN
|
||||
ALLOCATE (drho_h(4, na, nr, nspins), drho_s(4, na, nr, nspins), &
|
||||
drho0_h(4, na, nr, nspins), drho0_s(4, na, nr, nspins), &
|
||||
drho1_h(4, na, nr, nspins), drho1_s(4, na, nr, nspins))
|
||||
ALLOCATE (vxg_h(3, na, nr, nspins), vxg_s(3, na, nr, nspins))
|
||||
END IF
|
||||
IF (tau_f) THEN
|
||||
ALLOCATE (tau_h(na, nr, nspins), tau_s(na, nr, nspins), &
|
||||
tau0_h(na, nr, nspins), tau0_s(na, nr, nspins), &
|
||||
tau1_h(na, nr, nspins), tau1_s(na, nr, nspins))
|
||||
ALLOCATE (vtau_h(na, nr, nspins), vtau_s(na, nr, nspins))
|
||||
END IF
|
||||
!
|
||||
! NLCC: prepare rho and drho of the core charge for this KIND
|
||||
donlcc = .FALSE.
|
||||
IF (nlcc) THEN
|
||||
NULLIFY (rho_nlcc)
|
||||
rho_nlcc => kind_set(ikind)%nlcc_pot
|
||||
IF (ASSOCIATED(rho_nlcc)) donlcc = .TRUE.
|
||||
END IF
|
||||
|
||||
! Distribute the atoms of this kind
|
||||
num_pe = para_env%num_pe
|
||||
bo = get_limit(natom, num_pe, para_env%mepos)
|
||||
|
||||
DO iat = bo(1), bo(2)
|
||||
iatom = atom_list(iat)
|
||||
!
|
||||
NULLIFY (int_hh, int_ss)
|
||||
rho0_atom => rho0_atom_set(iatom)
|
||||
CALL get_rho_atom(rho_atom=rho0_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
|
||||
ALLOCATE (fint_ss(nspins), fint_hh(nspins))
|
||||
DO ns = 1, nspins
|
||||
nf = SIZE(int_ss(ns)%r_coef, 1)
|
||||
ALLOCATE (fint_ss(ns)%r_coef(nf, nf))
|
||||
nf = SIZE(int_hh(ns)%r_coef, 1)
|
||||
ALLOCATE (fint_hh(ns)%r_coef(nf, nf))
|
||||
END DO
|
||||
|
||||
! RHO0
|
||||
rho0_h = 0.0_dp
|
||||
rho0_s = 0.0_dp
|
||||
rho0_atom => rho0_atom_set(iatom)
|
||||
IF (gradient_f) THEN
|
||||
NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
|
||||
CALL get_rho_atom(rho_atom=rho0_atom, rho_rad_h=r_h, rho_rad_s=r_s, drho_rad_h=dr_h, &
|
||||
drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
|
||||
drho0_h = 0.0_dp
|
||||
drho0_s = 0.0_dp
|
||||
ELSE
|
||||
NULLIFY (r_h, r_s)
|
||||
CALL get_rho_atom(rho_atom=rho0_atom, rho_rad_h=r_h, rho_rad_s=r_s)
|
||||
rho_d = 0.0_dp
|
||||
END IF
|
||||
DO ir = 1, nr
|
||||
CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
|
||||
ir, r_h, r_s, rho0_h, rho0_s, dr_h, dr_s, &
|
||||
r_h_d, r_s_d, drho0_h, drho0_s)
|
||||
IF (donlcc) THEN
|
||||
CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
|
||||
ir, rho_nlcc(:, 1), rho0_h, rho0_s, rho_nlcc(:, 2), drho0_h, drho0_s)
|
||||
END IF
|
||||
END DO
|
||||
IF (tau_f) THEN
|
||||
!compute tau on the grid all at once
|
||||
CALL calc_tau_atom(tau0_h, tau0_s, rho0_atom, kind_set(ikind), nspins)
|
||||
ELSE
|
||||
tau_d = 0.0_dp
|
||||
END IF
|
||||
! RHO1
|
||||
rho1_h = 0.0_dp
|
||||
rho1_s = 0.0_dp
|
||||
rho1_atom => rho1_atom_set(iatom)
|
||||
IF (gradient_f) THEN
|
||||
NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
|
||||
CALL get_rho_atom(rho_atom=rho1_atom, rho_rad_h=r_h, rho_rad_s=r_s, drho_rad_h=dr_h, &
|
||||
drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
|
||||
drho1_h = 0.0_dp
|
||||
drho1_s = 0.0_dp
|
||||
ELSE
|
||||
NULLIFY (r_h, r_s)
|
||||
CALL get_rho_atom(rho_atom=rho1_atom, rho_rad_h=r_h, rho_rad_s=r_s)
|
||||
END IF
|
||||
DO ir = 1, nr
|
||||
CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
|
||||
ir, r_h, r_s, rho1_h, rho1_s, dr_h, dr_s, &
|
||||
r_h_d, r_s_d, drho1_h, drho1_s)
|
||||
END DO
|
||||
IF (tau_f) THEN
|
||||
!compute tau on the grid all at once
|
||||
CALL calc_tau_atom(tau1_h, tau1_s, rho1_atom, kind_set(ikind), nspins)
|
||||
END IF
|
||||
|
||||
DO ir = 1, nr
|
||||
IF (tau_f) THEN
|
||||
CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, tau1_h, na, ir)
|
||||
CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, tau1_s, na, ir)
|
||||
ELSE IF (gradient_f) THEN
|
||||
CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, tau_d, na, ir)
|
||||
CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, tau_d, na, ir)
|
||||
ELSE
|
||||
CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, rho_d, tau_d, na, ir)
|
||||
CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, rho_d, tau_d, na, ir)
|
||||
END IF
|
||||
END DO
|
||||
|
||||
! RHO2
|
||||
rho2_atom => rho2_atom_set(iatom)
|
||||
|
||||
DO istep = -nstep, nstep
|
||||
|
||||
beta = REAL(istep, KIND=dp)*epsrho
|
||||
|
||||
rho_h = rho0_h + beta*rho1_h
|
||||
rho_s = rho0_s + beta*rho1_s
|
||||
IF (gradient_f) THEN
|
||||
drho_h = drho0_h + beta*drho1_h
|
||||
drho_s = drho0_s + beta*drho1_s
|
||||
END IF
|
||||
IF (tau_f) THEN
|
||||
tau_h = tau0_h + beta*tau1_h
|
||||
tau_s = tau0_s + beta*tau1_s
|
||||
END IF
|
||||
!
|
||||
IF (gradient_f) THEN
|
||||
drho_h(4, :, :, :) = SQRT( &
|
||||
drho_h(1, :, :, :)*drho_h(1, :, :, :) + &
|
||||
drho_h(2, :, :, :)*drho_h(2, :, :, :) + &
|
||||
drho_h(3, :, :, :)*drho_h(3, :, :, :))
|
||||
|
||||
drho_s(4, :, :, :) = SQRT( &
|
||||
drho_s(1, :, :, :)*drho_s(1, :, :, :) + &
|
||||
drho_s(2, :, :, :)*drho_s(2, :, :, :) + &
|
||||
drho_s(3, :, :, :)*drho_s(3, :, :, :))
|
||||
END IF
|
||||
|
||||
DO ir = 1, nr
|
||||
IF (tau_f) THEN
|
||||
CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
|
||||
CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
|
||||
ELSE IF (gradient_f) THEN
|
||||
CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_d, na, ir)
|
||||
CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_d, na, ir)
|
||||
ELSE
|
||||
CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rho_d, tau_d, na, ir)
|
||||
CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, rho_d, tau_d, na, ir)
|
||||
END IF
|
||||
END DO
|
||||
|
||||
! hard atom density !
|
||||
CALL xc_dset_zero_all(deriv_set)
|
||||
!deb CALL vxc_of_r_new(xc_fun_section, rho_set_h, deriv_set, 1, needs, weight, &
|
||||
!deb lsd, na, nr, exc_h, vxc_h, vxg_h, vtau_h)
|
||||
CALL xc_2nd_deriv_of_r(xc_section=xc_section, &
|
||||
rho_set=rho_set_h, rho1_set=rho1_set_h, &
|
||||
deriv_set=deriv_set, &
|
||||
w=weight, vxc=vxc_h, vxg=vxg_h, do_triplet=is_triplet)
|
||||
! soft atom density !
|
||||
CALL xc_dset_zero_all(deriv_set)
|
||||
!deb CALL vxc_of_r_new(xc_fun_section, rho_set_s, deriv_set, 1, needs, weight, &
|
||||
!deb lsd, na, nr, exc_s, vxc_s, vxg_s, vtau_s)
|
||||
CALL xc_2nd_deriv_of_r(xc_section=xc_section, &
|
||||
rho_set=rho_set_s, rho1_set=rho1_set_s, &
|
||||
deriv_set=deriv_set, &
|
||||
w=weight, vxc=vxc_s, vxg=vxg_s, do_triplet=is_triplet)
|
||||
! potentials
|
||||
DO ns = 1, nspins
|
||||
fint_hh(ns)%r_coef(:, :) = 0.0_dp
|
||||
fint_ss(ns)%r_coef(:, :) = 0.0_dp
|
||||
END DO
|
||||
IF (gradient_f) THEN
|
||||
CALL gaVxcgb_GC(vxc_h, vxc_s, vxg_h, vxg_s, fint_hh, fint_ss, &
|
||||
grid_atom, basis_1c, harmonics, nspins)
|
||||
ELSE
|
||||
CALL gaVxcgb_noGC(vxc_h, vxc_s, fint_hh, fint_ss, &
|
||||
grid_atom, basis_1c, harmonics, nspins)
|
||||
END IF
|
||||
IF (tau_f) THEN
|
||||
CALL dgaVtaudgb(vtau_h, vtau_s, fint_hh, fint_ss, &
|
||||
grid_atom, basis_1c, harmonics, nspins)
|
||||
END IF
|
||||
! second derivative gxc
|
||||
NULLIFY (int_hh, int_ss)
|
||||
CALL get_rho_atom(rho_atom=rho2_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
|
||||
DO ns = 1, nspins
|
||||
int_ss(ns)%r_coef(:, :) = int_ss(ns)%r_coef(:, :) + oeps1*ak(istep)*fint_ss(ns)%r_coef(:, :)
|
||||
int_hh(ns)%r_coef(:, :) = int_hh(ns)%r_coef(:, :) + oeps1*ak(istep)*fint_hh(ns)%r_coef(:, :)
|
||||
END DO
|
||||
END DO
|
||||
!
|
||||
DO ns = 1, nspins
|
||||
DEALLOCATE (fint_ss(ns)%r_coef)
|
||||
DEALLOCATE (fint_hh(ns)%r_coef)
|
||||
END DO
|
||||
DEALLOCATE (fint_ss, fint_hh)
|
||||
|
||||
END DO ! iat
|
||||
|
||||
! Release the xc structure used to store the xc derivatives
|
||||
CALL xc_dset_release(deriv_set)
|
||||
CALL xc_rho_set_release(rho_set_h)
|
||||
CALL xc_rho_set_release(rho_set_s)
|
||||
CALL xc_rho_set_release(rho1_set_h)
|
||||
CALL xc_rho_set_release(rho1_set_s)
|
||||
|
||||
DEALLOCATE (rho_h, rho_s, rho0_h, rho0_s, rho1_h, rho1_s)
|
||||
DEALLOCATE (vxc_h, vxc_s)
|
||||
IF (gradient_f) THEN
|
||||
DEALLOCATE (drho_h, drho_s, drho0_h, drho0_s, drho1_h, drho1_s)
|
||||
DEALLOCATE (vxg_h, vxg_s)
|
||||
END IF
|
||||
IF (tau_f) THEN
|
||||
DEALLOCATE (tau_h, tau_s, tau0_h, tau0_s, tau1_h, tau1_s)
|
||||
DEALLOCATE (vtau_h, vtau_s)
|
||||
END IF
|
||||
|
||||
END DO ! ikind
|
||||
|
||||
END IF !xc_none
|
||||
|
||||
CALL timestop(handle)
|
||||
|
||||
END SUBROUTINE gfxc_atom_diff
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param grid_atom ...
|
||||
|
|
|
|||
|
|
@ -16,5 +16,5 @@ h2o_dip09.inp 86 1e-09
|
|||
h2o_dip10.inp 86 1e-09 0.150644266924E+00
|
||||
h2o_dip11.inp 86 1e-09 0.804659694314E-01
|
||||
h2o_dip11b.inp 86 1e-09 0.104678363910E+00
|
||||
h2o_dip12.inp 86 5e-09 0.380069249337E+00
|
||||
h2o_dip12.inp 86 5e-09 0.379667325157E+00
|
||||
#EOF
|
||||
|
|
|
|||
|
|
@ -8,7 +8,7 @@ h2o_dip13.inp 86 1e-08
|
|||
h2o_dip13b.inp 86 1e-08 0.151671831678E+00
|
||||
h2o_dip13c.inp 86 1e-08 0.168627194788E+00
|
||||
h2o_dip14.inp 86 1e-08 0.815745822946E+00
|
||||
h2o_dip15.inp 86 1e-08 0.176997002028E+00
|
||||
h2o_dip15.inp 86 1e-08 0.176996979844E+00
|
||||
h2o_dip15b.inp 86 1e-08 0.153260708546E+00
|
||||
h2o_dip16.inp 86 1e-08 0.157128093302E+00
|
||||
h2o_dip16b.inp 86 1e-08 0.131279651123E+00
|
||||
|
|
|
|||
|
|
@ -4,8 +4,8 @@
|
|||
# 1 compares the last total energy in the file
|
||||
# for details see cp2k/tools/do_regtest
|
||||
#
|
||||
h2o_a01.inp 1 6.0E-11 -17.24005864917290
|
||||
h2o_a01.inp 1 6.0E-11 -17.24005860897760
|
||||
h2o_a02.inp 1 6.0E-11 -17.13768744679261
|
||||
h2o_a03.inp 1 6.0E-11 -17.20804254307528
|
||||
h2o_a04.inp 1 6.0E-11 -17.04272708960406
|
||||
h2o_a04.inp 1 6.0E-11 -17.04346339096423
|
||||
#EOF
|
||||
|
|
|
|||
|
|
@ -19,5 +19,5 @@ h2o_f32.inp 1 1.0E-11
|
|||
h2o_f33.inp 1 5.0E-10 -16.75717077859229
|
||||
h2o_f34.inp 1 5.0E-10 -16.75710415920780
|
||||
h2o_f35.inp 1 5.0E-10 -16.75717077859468
|
||||
h2o_f36.inp 1 5.0E-10 -16.68896977690735
|
||||
h2o_f36.inp 1 5.0E-10 -16.68272790093562
|
||||
#EOF
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue