diff --git a/src/exstates_types.F b/src/exstates_types.F index 5f18de00cc..2fdaa83c24 100644 --- a/src/exstates_types.F +++ b/src/exstates_types.F @@ -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 diff --git a/src/input_cp2k_exstate.F b/src/input_cp2k_exstate.F index c14c8f525d..c2b2e6ecc7 100644 --- a/src/input_cp2k_exstate.F +++ b/src/input_cp2k_exstate.F @@ -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", & diff --git a/src/qs_fxc.F b/src/qs_fxc.F index cdd5bce42c..2bdac8cc72 100644 --- a/src/qs_fxc.F +++ b/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 ... diff --git a/src/qs_tddfpt2_fhxc_forces.F b/src/qs_tddfpt2_fhxc_forces.F index 4f52a6de16..d6a084aa6d 100644 --- a/src/qs_tddfpt2_fhxc_forces.F +++ b/src/qs_tddfpt2_fhxc_forces.F @@ -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., & diff --git a/src/qs_vxc_atom.F b/src/qs_vxc_atom.F index 99c5a0fd21..f4dfbba9a8 100644 --- a/src/qs_vxc_atom.F +++ b/src/qs_vxc_atom.F @@ -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 ... diff --git a/tests/QS/regtest-debug-4/TEST_FILES b/tests/QS/regtest-debug-4/TEST_FILES index 1a0836e131..c4efefb33f 100644 --- a/tests/QS/regtest-debug-4/TEST_FILES +++ b/tests/QS/regtest-debug-4/TEST_FILES @@ -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 diff --git a/tests/QS/regtest-debug-5/TEST_FILES b/tests/QS/regtest-debug-5/TEST_FILES index ce075aa6c8..e88e492be4 100644 --- a/tests/QS/regtest-debug-5/TEST_FILES +++ b/tests/QS/regtest-debug-5/TEST_FILES @@ -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 diff --git a/tests/QS/regtest-tddfpt-admm/TEST_FILES b/tests/QS/regtest-tddfpt-admm/TEST_FILES index 7604c96d5f..4b4ffcb02b 100644 --- a/tests/QS/regtest-tddfpt-admm/TEST_FILES +++ b/tests/QS/regtest-tddfpt-admm/TEST_FILES @@ -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 diff --git a/tests/QS/regtest-tddfpt-force-3/TEST_FILES b/tests/QS/regtest-tddfpt-force-3/TEST_FILES index 90a5c7b7fc..b81f450938 100644 --- a/tests/QS/regtest-tddfpt-force-3/TEST_FILES +++ b/tests/QS/regtest-tddfpt-force-3/TEST_FILES @@ -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