diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 15ee5b58c9..a6c351768f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -207,6 +207,7 @@ list( hfx_compression_methods.F hfx_derivatives.F hfx_energy_potential.F + hfx_exx.F hfx_helpers.F hfx_libint_interface.F hfx_load_balance_methods.F @@ -687,7 +688,6 @@ list( rpa_gw_im_time_util.F rpa_gw_kpoints_util.F rpa_gw_sigma_x.F - rpa_hfx.F rpa_im_time.F rpa_main.F rpa_rse.F diff --git a/src/ec_env_types.F b/src/ec_env_types.F index 187225907f..031d86640b 100644 --- a/src/ec_env_types.F +++ b/src/ec_env_types.F @@ -16,7 +16,10 @@ MODULE ec_env_types USE dbcsr_api, ONLY: dbcsr_p_type USE dm_ls_scf_types, ONLY: ls_scf_env_type,& ls_scf_release - USE input_section_types, ONLY: section_vals_type + USE hfx_types, ONLY: hfx_release,& + hfx_type + USE input_section_types, ONLY: section_vals_release,& + section_vals_type USE kinds, ONLY: dp USE pw_types, ONLY: pw_release,& pw_type @@ -55,8 +58,12 @@ MODULE ec_env_types INTEGER :: factorization INTEGER :: ec_initial_guess REAL(KIND=dp) :: eps_default + LOGICAL :: do_ec_admm + LOGICAL :: do_ec_hfx LOGICAL :: should_update LOGICAL :: use_ls_solver + LOGICAL :: reuse_hfx + LOGICAL :: basis_inconsistent ! debug LOGICAL :: debug_forces = .FALSE. LOGICAL :: debug_stress = .FALSE. @@ -70,10 +77,10 @@ MODULE ec_env_types REAL(KIND=dp) :: mao_eps1 INTEGER :: mao_iolevel ! energy components - REAL(KIND=dp) :: etotal + REAL(KIND=dp) :: etotal, old_etotal REAL(KIND=dp) :: eband, ecore, exc, ehartree, vhxc REAL(KIND=dp) :: edispersion, efield_elec, & - efield_nuclear, exc_aux_fit + efield_nuclear, ex, exc_aux_fit ! forces TYPE(qs_force_type), DIMENSION(:), POINTER :: force => Null() ! full neighbor lists and corresponding task list @@ -96,16 +103,21 @@ MODULE ec_env_types TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef ! CP equations TYPE(qs_p_env_type), POINTER :: p_env - TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz, matrix_z, matrix_wz, z_admm + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz, matrix_wz, matrix_z, z_admm ! Harris (rhoout), and response density (rhoz) on grid - TYPE(pw_type), DIMENSION(:), POINTER :: rhoout_r, rhoz_r + TYPE(pw_type), DIMENSION(:), POINTER :: rhoout_r, rhoz_r ! potentials from input density TYPE(pw_type) :: vh_rspace - TYPE(pw_type), DIMENSION(:), POINTER :: vxc_rspace, vtau_rspace, vadmm_rspace + TYPE(pw_type), DIMENSION(:), POINTER :: vxc_rspace, vtau_rspace, vadmm_rspace ! efield TYPE(efield_berry_type), POINTER :: efield => NULL() ! LS matrices and types TYPE(ls_scf_env_type), POINTER :: ls_env + ! Environment for Hartree-Fock exchange + TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data + ! ADMM XC environments + TYPE(section_vals_type), POINTER :: xc_section_primary => Null(), & + xc_section_aux => Null() END TYPE energy_correction_type CONTAINS @@ -147,9 +159,9 @@ CONTAINS CALL qs_dispersion_release(ec_env%dispersion_env) END IF - IF (ASSOCIATED(ec_env%matrix_z)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_z) IF (ASSOCIATED(ec_env%matrix_hz)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_hz) IF (ASSOCIATED(ec_env%matrix_wz)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_wz) + IF (ASSOCIATED(ec_env%matrix_z)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_z) IF (ASSOCIATED(ec_env%z_admm)) CALL dbcsr_deallocate_matrix_set(ec_env%z_admm) NULLIFY (ec_env%matrix_z, ec_env%matrix_hz, ec_env%matrix_wz) NULLIFY (ec_env%z_admm) @@ -186,6 +198,13 @@ CONTAINS CALL ls_scf_release(ec_env%ls_env) END IF + IF (.NOT. ec_env%reuse_hfx) THEN + IF (ASSOCIATED(ec_env%x_data)) CALL hfx_release(ec_env%x_data) + END IF + + IF (ASSOCIATED(ec_env%xc_section_aux)) CALL section_vals_release(ec_env%xc_section_aux) + IF (ASSOCIATED(ec_env%xc_section_primary)) CALL section_vals_release(ec_env%xc_section_primary) + DEALLOCATE (ec_env) END IF diff --git a/src/ec_environment.F b/src/ec_environment.F index 5a005e0095..619881516a 100644 --- a/src/ec_environment.F +++ b/src/ec_environment.F @@ -126,8 +126,9 @@ CONTAINS 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_fun_section, xc_section + TYPE(section_vals_type), POINTER :: ec_hfx_section, nl_section, pp_section, & + section1, section2, xc_fun_section, & + xc_section CALL timeset(routineN, handle) @@ -149,8 +150,12 @@ CONTAINS NULLIFY (ec_env%vtau_rspace) NULLIFY (ec_env%vadmm_rspace) NULLIFY (ec_env%rhoout_r, ec_env%rhoz_r) + NULLIFY (ec_env%x_data) ec_env%should_update = .TRUE. ec_env%mao = .FALSE. + ec_env%do_ec_admm = .FALSE. + ec_env%do_ec_hfx = .FALSE. + ec_env%reuse_hfx = .FALSE. IF (qs_env%energy_correction) THEN @@ -193,6 +198,8 @@ CONTAINS l_val=ec_env%debug_forces) CALL section_vals_val_get(ec_section, "DEBUG_STRESS", & l_val=ec_env%debug_stress) + ! ADMM + CALL section_vals_val_get(ec_section, "ADMM", l_val=ec_env%do_ec_admm) ec_env%do_skip = .FALSE. @@ -250,10 +257,11 @@ CONTAINS CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, basis_type="HARRIS") CALL init_orbital_pointers(maxlgto + 1) ! - ! Density-corrected DFT must be performed with the same basis as ground-state CALL uppercase(ec_env%basis) + ! Basis may only differ from ground-state if explicitly added - IF (ec_env%energy_functional == ec_functional_dc .AND. ec_env%basis == "HARRIS") THEN + ec_env%basis_inconsistent = .FALSE. + IF (ec_env%basis == "HARRIS") THEN DO ikind = 1, nkind qs_kind => qs_kind_set(ikind) ! Basis sets of ground-state @@ -262,10 +270,17 @@ CONTAINS CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS") IF (basis_set%name .NE. harris_basis%name) THEN - CPABORT("DC-DFT: Correction and ground state need to use the same basis") + ec_env%basis_inconsistent = .TRUE. END IF END DO END IF + + !Density-corrected DFT must be performed with the same basis as ground-state + IF (ec_env%energy_functional == ec_functional_dc .AND. ec_env%basis_inconsistent) THEN + CALL cp_abort(__LOCATION__, & + "DC-DFT: Correction and ground state need to use the same basis."// & + "Checked by comparing basis set names only.") + END IF ! ! set functional SELECT CASE (ec_env%energy_functional) @@ -314,6 +329,10 @@ CONTAINS END IF ec_env%dispersion_env => dispersion_env + ! Check if hybrid functional are used + ec_hfx_section => section_vals_get_subs_vals(ec_section, "XC%HF") + CALL section_vals_get(ec_hfx_section, explicit=ec_env%do_ec_hfx) + ! Initialize Harris LS solver environment ec_env%use_ls_solver = .FALSE. ec_env%use_ls_solver = (ec_env%ks_solver .EQ. ec_matrix_sign) & @@ -324,11 +343,6 @@ CONTAINS CALL ec_ls_create(qs_env, ec_env) END IF - ! Write input - IF (unit_nr > 0) THEN - CALL ec_write_input(ec_env, unit_nr) - END IF - END IF CALL timestop(handle) @@ -455,158 +469,177 @@ CONTAINS !> \brief Print out the energy correction input section !> !> \param ec_env ... -!> \param unit_nr ... !> \par History !> 2020.10 created [Fabian Belleflamme] !> \author Fabian Belleflamme ! ************************************************************************************************** - SUBROUTINE ec_write_input(ec_env, unit_nr) + SUBROUTINE ec_write_input(ec_env) TYPE(energy_correction_type), POINTER :: ec_env - INTEGER, INTENT(IN) :: unit_nr CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_write_input' - INTEGER :: handle + INTEGER :: handle, unit_nr + TYPE(cp_logger_type), POINTER :: logger TYPE(ls_scf_env_type), POINTER :: ls_env CALL timeset(routineN, handle) - WRITE (unit_nr, '(T2,A)') & - "!"//REPEAT("-", 29)//" Energy Correction "//REPEAT("-", 29)//"!" + logger => cp_get_default_logger() + IF (logger%para_env%is_source()) THEN + unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) + ELSE + unit_nr = -1 + END IF - ! Type of energy correction - SELECT CASE (ec_env%energy_functional) - CASE (ec_functional_harris) - WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Correction: ", "HARRIS FUNCTIONAL" - CASE (ec_functional_dc) - WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Correction: ", "DC-DFT" - END SELECT - WRITE (unit_nr, '()') + IF (unit_nr > 0) THEN - ! Energy correction parameters - WRITE (unit_nr, '(T2,A,T61,E20.3)') "eps_default:", ec_env%eps_default + WRITE (unit_nr, '(T2,A)') & + "!"//REPEAT("-", 29)//" Energy Correction "//REPEAT("-", 29)//"!" - CALL uppercase(ec_env%basis) - SELECT CASE (ec_env%basis) - CASE ("ORBITAL") - WRITE (unit_nr, '(T2,A,T61,A20)') "EC basis: ", "ORBITAL" - CASE ("PRIMITIVE") - WRITE (unit_nr, '(T2,A,T61,A20)') "EC basis: ", "PRIMITIVE" - CASE ("HARRIS") - WRITE (unit_nr, '(T2,A,T61,A20)') "EC Basis: ", "HARRIS" - END SELECT - - ! Parameters for Harris functional solver - IF (ec_env%energy_functional == ec_functional_harris) THEN - - ! Algorithm - SELECT CASE (ec_env%ks_solver) - CASE (ec_diagonalization) - WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "DIAGONALIZATION" - CASE (ec_ot_diag) - WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "OT DIAGONALIZATION" - CASE (ec_matrix_sign) - WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "MATRIX_SIGN" - CASE (ec_matrix_trs4) - WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "TRS4" - CALL cite_reference(Niklasson2003) - CASE (ec_matrix_tc2) - WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "TC2" - CALL cite_reference(Niklasson2014) + ! Type of energy correction + SELECT CASE (ec_env%energy_functional) + CASE (ec_functional_harris) + WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Correction: ", "HARRIS FUNCTIONAL" + CASE (ec_functional_dc) + WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Correction: ", "DC-DFT" END SELECT WRITE (unit_nr, '()') - ! MAO - IF (ec_env%mao) THEN - WRITE (unit_nr, '(T2,A,T61,L20)') "MAO:", ec_env%mao - WRITE (unit_nr, '(T2,A,T61,L20)') "MAO_IOLEVEL:", ec_env%mao_iolevel - WRITE (unit_nr, '(T2,A,T61,I20)') "MAO_MAX_ITER:", ec_env%mao_max_iter - WRITE (unit_nr, '(T2,A,T61,E20.3)') "MAO_EPS_GRAD:", ec_env%mao_eps_grad - WRITE (unit_nr, '(T2,A,T61,E20.3)') "MAO_EPS1:", ec_env%mao_eps1 - WRITE (unit_nr, '()') - END IF + ! Energy correction parameters + WRITE (unit_nr, '(T2,A,T61,E20.3)') "eps_default:", ec_env%eps_default - ! Parameters for linear response solver - IF (.NOT. ec_env%use_ls_solver) THEN + CALL uppercase(ec_env%basis) + SELECT CASE (ec_env%basis) + CASE ("ORBITAL") + WRITE (unit_nr, '(T2,A,T61,A20)') "EC basis: ", "ORBITAL" + CASE ("PRIMITIVE") + WRITE (unit_nr, '(T2,A,T61,A20)') "EC basis: ", "PRIMITIVE" + CASE ("HARRIS") + WRITE (unit_nr, '(T2,A,T61,A20)') "EC Basis: ", "HARRIS" + END SELECT - WRITE (unit_nr, '(T2,A)') "MO Solver" - WRITE (unit_nr, '()') + ! Info how HFX in energy correction is treated + IF (ec_env%do_ec_hfx) THEN + WRITE (unit_nr, '(T2,A,T61,L20)') "DC-DFT with HFX", ec_env%do_ec_hfx + WRITE (unit_nr, '(T2,A,T61,L20)') "Reuse HFX integrals", ec_env%reuse_hfx + WRITE (unit_nr, '(T2,A,T61,L20)') "DC-DFT HFX with ADMM", ec_env%do_ec_admm + + END IF ! ec_env%do_ec_hfx + + ! Parameters for Harris functional solver + IF (ec_env%energy_functional == ec_functional_harris) THEN + + ! Algorithm SELECT CASE (ec_env%ks_solver) CASE (ec_diagonalization) - - SELECT CASE (ec_env%factorization) - CASE (kg_cholesky) - WRITE (unit_nr, '(T2,A,T61,A20)') "Factorization: ", "CHOLESKY" - END SELECT - + WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "DIAGONALIZATION" CASE (ec_ot_diag) - - ! OT Diagonalization - ! Initial guess : 1) block diagonal initial guess - ! 2) GS-density matrix (might require trafo if basis diff) - - SELECT CASE (ec_env%ec_initial_guess) - CASE (ec_ot_atomic) - WRITE (unit_nr, '(T2,A,T61,A20)') "OT Diag initial guess: ", "ATOMIC" - CASE (ec_ot_gs) - WRITE (unit_nr, '(T2,A,T61,A20)') "OT Diag initial guess: ", "GROUND STATE DM" - END SELECT - - CASE DEFAULT - CPABORT("Unknown Diagonalization algorithm for Harris functional") + WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "OT DIAGONALIZATION" + CASE (ec_matrix_sign) + WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "MATRIX_SIGN" + CASE (ec_matrix_trs4) + WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "TRS4" + CALL cite_reference(Niklasson2003) + CASE (ec_matrix_tc2) + WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "TC2" + CALL cite_reference(Niklasson2014) END SELECT - - ELSE - - WRITE (unit_nr, '(T2,A)') "AO Solver" WRITE (unit_nr, '()') - ls_env => ec_env%ls_env - WRITE (unit_nr, '(T2,A,T61,E20.3)') "eps_filter:", ls_env%eps_filter - WRITE (unit_nr, '(T2,A,T61,L20)') "fixed chemical potential (mu)", ls_env%fixed_mu - WRITE (unit_nr, '(T2,A,T61,L20)') "Computing inv(S):", ls_env%needs_s_inv - WRITE (unit_nr, '(T2,A,T61,L20)') "Computing sqrt(S):", ls_env%use_s_sqrt - WRITE (unit_nr, '(T2,A,T61,L20)') "Computing S preconditioner ", ls_env%has_s_preconditioner - WRITE (unit_nr, '(T2,A,T61,L20)') "Use single precision matrices", & - ls_env%ls_mstruct%single_precision - - IF (ls_env%use_s_sqrt) THEN - SELECT CASE (ls_env%s_sqrt_method) - CASE (ls_s_sqrt_ns) - WRITE (unit_nr, '(T2,A,T61,A20)') "S sqrt method:", "NEWTONSCHULZ" - CASE (ls_s_sqrt_proot) - WRITE (unit_nr, '(T2,A,T61,A20)') "S sqrt method:", "PROOT" - CASE DEFAULT - CPABORT("Unknown sqrt method.") - END SELECT - WRITE (unit_nr, '(T2,A,T61,I20)') "S sqrt order:", ls_env%s_sqrt_order + ! MAO + IF (ec_env%mao) THEN + WRITE (unit_nr, '(T2,A,T61,L20)') "MAO:", ec_env%mao + WRITE (unit_nr, '(T2,A,T61,L20)') "MAO_IOLEVEL:", ec_env%mao_iolevel + WRITE (unit_nr, '(T2,A,T61,I20)') "MAO_MAX_ITER:", ec_env%mao_max_iter + WRITE (unit_nr, '(T2,A,T61,E20.3)') "MAO_EPS_GRAD:", ec_env%mao_eps_grad + WRITE (unit_nr, '(T2,A,T61,E20.3)') "MAO_EPS1:", ec_env%mao_eps1 + WRITE (unit_nr, '()') END IF - SELECT CASE (ls_env%s_preconditioner_type) - CASE (ls_s_preconditioner_none) - WRITE (unit_nr, '(T2,A,T61,A20)') "S preconditioner type ", "NONE" - CASE (ls_s_preconditioner_atomic) - WRITE (unit_nr, '(T2,A,T61,A20)') "S preconditioner type ", "ATOMIC" - CASE (ls_s_preconditioner_molecular) - WRITE (unit_nr, '(T2,A,T61,A20)') "S preconditioner type ", "MOLECULAR" - END SELECT + ! Parameters for linear response solver + IF (.NOT. ec_env%use_ls_solver) THEN - SELECT CASE (ls_env%ls_mstruct%cluster_type) - CASE (ls_cluster_atomic) - WRITE (unit_nr, '(T2,A,T61,A20)') "Cluster type", ADJUSTR("ATOMIC") - CASE (ls_cluster_molecular) - WRITE (unit_nr, '(T2,A,T61,A20)') "Cluster type", ADJUSTR("MOLECULAR") - CASE DEFAULT - CPABORT("Unknown cluster type") - END SELECT + WRITE (unit_nr, '(T2,A)') "MO Solver" + WRITE (unit_nr, '()') - END IF + SELECT CASE (ec_env%ks_solver) + CASE (ec_diagonalization) - END IF + SELECT CASE (ec_env%factorization) + CASE (kg_cholesky) + WRITE (unit_nr, '(T2,A,T61,A20)') "Factorization: ", "CHOLESKY" + END SELECT - WRITE (unit_nr, '(T2,A)') REPEAT("-", 79) - WRITE (unit_nr, '()') + CASE (ec_ot_diag) + + ! OT Diagonalization + ! Initial guess : 1) block diagonal initial guess + ! 2) GS-density matrix (might require trafo if basis diff) + + SELECT CASE (ec_env%ec_initial_guess) + CASE (ec_ot_atomic) + WRITE (unit_nr, '(T2,A,T61,A20)') "OT Diag initial guess: ", "ATOMIC" + CASE (ec_ot_gs) + WRITE (unit_nr, '(T2,A,T61,A20)') "OT Diag initial guess: ", "GROUND STATE DM" + END SELECT + + CASE DEFAULT + CPABORT("Unknown Diagonalization algorithm for Harris functional") + END SELECT + + ELSE + + WRITE (unit_nr, '(T2,A)') "AO Solver" + WRITE (unit_nr, '()') + + ls_env => ec_env%ls_env + WRITE (unit_nr, '(T2,A,T61,E20.3)') "eps_filter:", ls_env%eps_filter + WRITE (unit_nr, '(T2,A,T61,L20)') "fixed chemical potential (mu)", ls_env%fixed_mu + WRITE (unit_nr, '(T2,A,T61,L20)') "Computing inv(S):", ls_env%needs_s_inv + WRITE (unit_nr, '(T2,A,T61,L20)') "Computing sqrt(S):", ls_env%use_s_sqrt + WRITE (unit_nr, '(T2,A,T61,L20)') "Computing S preconditioner ", ls_env%has_s_preconditioner + WRITE (unit_nr, '(T2,A,T61,L20)') "Use single precision matrices", & + ls_env%ls_mstruct%single_precision + + IF (ls_env%use_s_sqrt) THEN + SELECT CASE (ls_env%s_sqrt_method) + CASE (ls_s_sqrt_ns) + WRITE (unit_nr, '(T2,A,T61,A20)') "S sqrt method:", "NEWTONSCHULZ" + CASE (ls_s_sqrt_proot) + WRITE (unit_nr, '(T2,A,T61,A20)') "S sqrt method:", "PROOT" + CASE DEFAULT + CPABORT("Unknown sqrt method.") + END SELECT + WRITE (unit_nr, '(T2,A,T61,I20)') "S sqrt order:", ls_env%s_sqrt_order + END IF + + SELECT CASE (ls_env%s_preconditioner_type) + CASE (ls_s_preconditioner_none) + WRITE (unit_nr, '(T2,A,T61,A20)') "S preconditioner type ", "NONE" + CASE (ls_s_preconditioner_atomic) + WRITE (unit_nr, '(T2,A,T61,A20)') "S preconditioner type ", "ATOMIC" + CASE (ls_s_preconditioner_molecular) + WRITE (unit_nr, '(T2,A,T61,A20)') "S preconditioner type ", "MOLECULAR" + END SELECT + + SELECT CASE (ls_env%ls_mstruct%cluster_type) + CASE (ls_cluster_atomic) + WRITE (unit_nr, '(T2,A,T61,A20)') "Cluster type", ADJUSTR("ATOMIC") + CASE (ls_cluster_molecular) + WRITE (unit_nr, '(T2,A,T61,A20)') "Cluster type", ADJUSTR("MOLECULAR") + CASE DEFAULT + CPABORT("Unknown cluster type") + END SELECT + + END IF + + END IF ! if ec_functional_harris + + WRITE (unit_nr, '(T2,A)') REPEAT("-", 79) + WRITE (unit_nr, '()') + + END IF ! unit_nr CALL timestop(handle) diff --git a/src/energy_corrections.F b/src/energy_corrections.F index 520e34a880..9fcafe30f0 100644 --- a/src/energy_corrections.F +++ b/src/energy_corrections.F @@ -11,10 +11,12 @@ !> 03.2014 created !> 09.2019 Moved from KG to Kohn-Sham !> 08.2022 Add Density-Corrected DFT methods +!> 04.2023 Add hybrid functionals for DC-DFT !> \author JGH ! ************************************************************************************************** MODULE energy_corrections - USE admm_methods, ONLY: admm_mo_merge_ks_matrix + USE admm_dm_methods, ONLY: admm_dm_calc_rho_aux + USE admm_methods, ONLY: admm_mo_calc_rho_aux USE atomic_kind_types, ONLY: atomic_kind_type,& get_atomic_kind USE basis_set_types, ONLY: get_gto_basis_set,& @@ -82,8 +84,9 @@ MODULE energy_corrections USE external_potential_types, ONLY: get_potential,& gth_potential_type,& sgp_potential_type + USE hfx_exx, ONLY: add_exx_to_rhs,& + calculate_exx USE input_constants, ONLY: & - do_admm_basis_projection, do_admm_exch_scaling_none, do_admm_purify_none, & ec_diagonalization, ec_functional_dc, ec_functional_harris, ec_matrix_sign, ec_matrix_tc2, & ec_matrix_trs4, ec_ot_atomic, ec_ot_diag, ec_ot_gs, ot_precond_full_single_inverse, & ot_precond_solver_default, vdw_pairpot_dftd3, vdw_pairpot_dftd3bj, xc_vdw_fun_pairpot @@ -154,7 +157,6 @@ MODULE energy_corrections USE qs_ks_methods, ONLY: calc_rho_tot_gspace USE qs_ks_reference, ONLY: ks_ref_potential USE qs_ks_types, ONLY: qs_ks_env_type - USE qs_linres_kernel, ONLY: hfx_matrix USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues,& make_basis_sm USE qs_mo_types, ONLY: deallocate_mo_set,& @@ -174,6 +176,7 @@ MODULE energy_corrections USE qs_vxc, ONLY: qs_vxc_create USE response_solver, ONLY: response_calculation,& response_force + USE rtp_admm_methods, ONLY: rtp_admm_calc_rho_aux USE string_utilities, ONLY: uppercase USE task_list_methods, ONLY: generate_qs_task_list USE task_list_types, ONLY: allocate_task_list,& @@ -247,13 +250,20 @@ CONTAINS IF (PRESENT(calculate_forces)) my_calc_forces = calculate_forces IF (ec_env%should_update) THEN + ec_env%old_etotal = 0.0_dp ec_env%etotal = 0.0_dp ec_env%eband = 0.0_dp ec_env%ehartree = 0.0_dp + ec_env%ex = 0.0_dp ec_env%exc = 0.0_dp ec_env%vhxc = 0.0_dp ec_env%edispersion = 0.0_dp ec_env%exc_aux_fit = 0.0_dp + + ! Save total energy of reference calculation + CALL get_qs_env(qs_env, energy=energy) + ec_env%old_etotal = energy%total + END IF IF (my_calc_forces) THEN @@ -273,9 +283,9 @@ CONTAINS ! Perform the energy correction CALL energy_correction_low(qs_env, ec_env, my_calc_forces, unit_nr) - CALL get_qs_env(qs_env, energy=energy) + ! Update total energy in qs environment and amount fo correction IF (ec_env%should_update) THEN - energy%nonscf_correction = ec_env%etotal - energy%total + energy%nonscf_correction = ec_env%etotal - ec_env%old_etotal energy%total = ec_env%etotal END IF @@ -603,18 +613,18 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_dc_build_ks_matrix_force' INTEGER :: handle, i, iounit, ispin, natom, nspins - LOGICAL :: debug_forces, debug_stress, & - do_adiabatic_rescaling, do_hfx, & + LOGICAL :: debug_forces, debug_stress, do_ec_hfx, & use_virial - REAL(dp) :: ehartree, eovrl, exc, fconv + REAL(dp) :: dummy_real, dummy_real2(2), ehartree, & + eovrl, exc, fconv REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ftot - REAL(dp), DIMENSION(3) :: fodeb + REAL(dp), DIMENSION(3) :: fodeb, fodeb2 REAL(KIND=dp), DIMENSION(3, 3) :: h_stress, pv_loc, stdeb, sttot TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(cell_type), POINTER :: cell TYPE(cp_logger_type), POINTER :: logger - TYPE(dbcsr_p_type) :: scrm - TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dbcsr_work, matrix_ks, matrix_p, matrix_s + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, scrm + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p TYPE(dft_control_type), POINTER :: dft_control TYPE(mp_para_env_type), POINTER :: para_env TYPE(neighbor_list_set_p_type), DIMENSION(:), & @@ -630,8 +640,7 @@ CONTAINS TYPE(qs_force_type), DIMENSION(:), POINTER :: force TYPE(qs_ks_env_type), POINTER :: ks_env TYPE(qs_rho_type), POINTER :: rho - TYPE(section_vals_type), POINTER :: adiabatic_rescaling_section, & - hfx_sections, input + TYPE(section_vals_type), POINTER :: ec_hfx_sections TYPE(virial_type), POINTER :: virial CALL timeset(routineN, handle) @@ -652,7 +661,6 @@ CONTAINS cell=cell, & dft_control=dft_control, & force=force, & - input=input, & ks_env=ks_env, & matrix_ks=matrix_ks, & matrix_s=matrix_s, & @@ -671,6 +679,9 @@ CONTAINS sttot = virial%pv_virial END IF + ! Get density matrix of reference calculation + CALL qs_rho_get(rho, rho_ao_kp=matrix_p) + NULLIFY (auxbas_pw_pool, poisson_env) ! gets the tmp grids CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, & @@ -716,7 +727,7 @@ CONTAINS CALL pw_transfer(v_hartree_gspace, v_hartree_rspace) CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol) - ! Save Harris on real space grid for use in properties + ! Save density on real space grid for use in properties CALL qs_rho_get(rho, rho_r=rho_r) ALLOCATE (ec_env%rhoout_r(nspins)) DO ispin = 1, nspins @@ -749,6 +760,7 @@ CONTAINS ! only activate stress calculation if IF (use_virial) virial%pv_calculate = .TRUE. + ! Exchange-correlation potential CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, & vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.) @@ -767,6 +779,16 @@ CONTAINS ! virial%pv_xc will be zeroed in the xc routines END IF + ! initialize srcm matrix + NULLIFY (scrm) + CALL dbcsr_allocate_matrix_set(scrm, nspins) + DO ispin = 1, nspins + ALLOCATE (scrm(ispin)%matrix) + CALL dbcsr_create(scrm(ispin)%matrix, template=ec_env%matrix_ks(ispin, 1)%matrix) + CALL dbcsr_copy(scrm(ispin)%matrix, ec_env%matrix_ks(ispin, 1)%matrix) + CALL dbcsr_set(scrm(ispin)%matrix, 0.0_dp) + END DO + pw_grid => v_hartree_rspace%pw_grid ALLOCATE (v_rspace_in(nspins)) DO ispin = 1, nspins @@ -782,11 +804,52 @@ CONTAINS CALL pw_axpy(ec_env%vh_rspace, v_rspace_in(ispin)) END DO - ! initialize src matrix - NULLIFY (scrm%matrix) - ALLOCATE (scrm%matrix) - CALL dbcsr_create(scrm%matrix, template=matrix_s(1)%matrix) - CALL cp_dbcsr_alloc_block_from_nbl(scrm%matrix, ec_env%sab_orb) +!------------------------------------------------ + + ! If hybrid functional in DC-DFT + ec_hfx_sections => section_vals_get_subs_vals(qs_env%input, "DFT%ENERGY_CORRECTION%XC%HF") + CALL section_vals_get(ec_hfx_sections, explicit=do_ec_hfx) + + IF (do_ec_hfx) THEN + + IF (debug_forces) fodeb(1:3) = force(1)%fock_4c(1:3, 1) + IF (debug_forces) fodeb2(1:3) = force(1)%overlap_admm(1:3, 1) + + ! Calculate direct HFX forces here + ! Virial contribution (fock_4c) done inside calculate_exx + dummy_real = 0.0_dp + CALL calculate_exx(qs_env=qs_env, & + unit_nr=iounit, & + hfx_sections=ec_hfx_sections, & + x_data=ec_env%x_data, & + do_gw=.FALSE., & + do_admm=ec_env%do_ec_admm, & + calc_forces=.TRUE., & + reuse_hfx=ec_env%reuse_hfx, & + do_im_time=.FALSE., & + E_ex_from_GW=dummy_real, & + E_admm_from_GW=dummy_real2, & + t3=dummy_real) + + IF (debug_forces) THEN + fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3) + CALL para_env%sum(fodeb) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*hfx_DC ", fodeb + + fodeb2(1:3) = force(1)%overlap_admm(1:3, 1) - fodeb2(1:3) + CALL para_env%sum(fodeb2) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*hfx_DC*S ", fodeb2 + END IF + IF (debug_stress .AND. use_virial) THEN + stdeb = -1.0_dp*fconv*virial%pv_fock_4c + CALL para_env%sum(stdeb) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| P*hfx_DC ', one_third_sum_diag(stdeb), det_3x3(stdeb) + END IF + + END IF + +!------------------------------------------------ ! Stress-tensor contribution derivative of integrand ! int v_Hxc[n^în]*n^out @@ -794,7 +857,6 @@ CONTAINS pv_loc = virial%pv_virial END IF - CALL qs_rho_get(rho, rho_ao=matrix_p) IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1) IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial @@ -804,8 +866,8 @@ CONTAINS CALL pw_axpy(v_hartree_rspace, v_rspace(ispin)) ! integrate over potential CALL integrate_v_rspace(v_rspace=v_rspace(ispin), & - hmat=scrm, & - pmat=matrix_p(ispin), & + hmat=scrm(ispin), & + pmat=matrix_p(ispin, 1), & qs_env=qs_env, & calculate_forces=.TRUE., & basis_type="HARRIS", & @@ -831,8 +893,8 @@ CONTAINS CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol) ! integrate over Tau-potential CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), & - hmat=scrm, & - pmat=matrix_p(ispin), & + hmat=scrm(ispin), & + pmat=matrix_p(ispin, 1), & qs_env=qs_env, & calculate_forces=.TRUE., & compute_tau=.TRUE., & @@ -858,14 +920,22 @@ CONTAINS virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc) END IF - ! delete scr matrix - CALL dbcsr_release(scrm%matrix) - DEALLOCATE (scrm%matrix) + ! delete scrm matrix + CALL dbcsr_deallocate_matrix_set(scrm) !---------------------------------------------------- ! Right-hand-side matrix B for linear response equations AX = B !---------------------------------------------------- + ! RHS = int v_Hxc[n]_DC - v_Hxc[n]_GS dr + alpha_DC * E_X[n] - alpha_gs * E_X[n] + ! = int v_Hxc[n]_DC - v_Hxc[n]_GS dr + alpha_DC / alpha_GS * E_X[n]_GS - E_X[n]_GS + ! + ! with v_Hxc[n] = v_H[n] + v_xc[n] + ! + ! Actually v_H[n_in] same for DC and GS, just there for convenience + ! v_xc[n_in]_GS = 0 if GS is HF BUT =/0 if hybrid + ! so, we keep this general form + NULLIFY (ec_env%matrix_hz) CALL dbcsr_allocate_matrix_set(ec_env%matrix_hz, nspins) DO ispin = 1, nspins @@ -876,14 +946,15 @@ CONTAINS END DO DO ispin = 1, nspins - ! v_rspace = v_Hxc[n_in]_DC - v_Hxc[n_in]_GS + ! v_rspace = v_rspace - v_rspace_in + ! = v_Hxc[n_in]_DC - v_Hxc[n_in]_GS CALL pw_axpy(v_rspace_in(ispin), v_rspace(ispin), -1.0_dp) END DO DO ispin = 1, nspins CALL integrate_v_rspace(v_rspace=v_rspace(ispin), & hmat=ec_env%matrix_hz(ispin), & - pmat=matrix_p(ispin), & + pmat=matrix_p(ispin, 1), & qs_env=qs_env, & calculate_forces=.FALSE., & basis_type="HARRIS", & @@ -911,7 +982,7 @@ CONTAINS ! integrate over Tau-potential CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), & hmat=ec_env%matrix_hz(ispin), & - pmat=matrix_p(ispin), & + pmat=matrix_p(ispin, 1), & qs_env=qs_env, & calculate_forces=.FALSE., compute_tau=.TRUE., & basis_type="HARRIS", & @@ -919,82 +990,17 @@ CONTAINS END DO END IF - ! Need to also subtract HFX contribution from ec_env%matrix_hz - hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF") - CALL section_vals_get(hfx_sections, explicit=do_hfx) - adiabatic_rescaling_section => section_vals_get_subs_vals(input, "DFT%XC%ADIABATIC_RESCALING") - CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling) - - IF (do_hfx .AND. .NOT. do_adiabatic_rescaling) THEN - - ! HFX-ADMM - IF (dft_control%do_admm) THEN - - IF (dft_control%admm_control%purification_method /= do_admm_purify_none) THEN - CPABORT("ADMM: Linear Response needs purification_method=none") - END IF - IF (dft_control%admm_control%scaling_model /= do_admm_exch_scaling_none) THEN - CPABORT("ADMM: Linear Response needs scaling_model=none") - END IF - IF (dft_control%admm_control%method /= do_admm_basis_projection) THEN - CPABORT("ADMM: Linear Response needs admm_method=basis_projection") - END IF - - ! HFX matrix in ADMM environment already calculated during ground-state - ! Here, need to subtract it from RHS matrix for the response equations - - ! make a work matrix to isolate the HFX matrix - ! dbcsr_work = old KS matrix - ALLOCATE (dbcsr_work(nspins)) - DO ispin = 1, nspins - ALLOCATE (dbcsr_work(ispin)%matrix) - CALL dbcsr_copy(dbcsr_work(ispin)%matrix, matrix_ks(ispin)%matrix) - END DO - - ! Can actually only call subroutine merge_ks_matrix_none - ! Adds to matrix_ks of qs_env - CALL admm_mo_merge_ks_matrix(qs_env) - - ! dbcsr_work = HFX matrix = (KS+HFX) - KS - DO ispin = 1, nspins - CALL dbcsr_add(dbcsr_work(ispin)%matrix, matrix_ks(ispin)%matrix, -1.0_dp, 1.0_dp) - END DO - - ! Subtract from RHS matrix - DO ispin = 1, nspins - CALL dbcsr_add(ec_env%matrix_hz(ispin)%matrix, dbcsr_work(ispin)%matrix, & - 1.0_dp, -1.0_dp) - END DO - - ! Restore Kohn-Sham matrix, to make it consistent for response calculation - DO ispin = 1, nspins - CALL dbcsr_copy(matrix_ks(ispin)%matrix, ec_env%matrix_ks(ispin, 1)%matrix) - END DO - - CALL dbcsr_deallocate_matrix_set(dbcsr_work) - - ! conventional HFX - ELSE - ALLOCATE (dbcsr_work(nspins)) - DO ispin = 1, nspins - ALLOCATE (dbcsr_work(ispin)%matrix) - CALL dbcsr_copy(dbcsr_work(ispin)%matrix, matrix_s(ispin)%matrix) - CALL dbcsr_set(dbcsr_work(ispin)%matrix, 0.0_dp) - END DO - - ! Get the HFX contribution to the Kohn-Sham matrix - CALL hfx_matrix(dbcsr_work, matrix_p, qs_env, hfx_sections) - - ! Subtract from RHS matrix - DO ispin = 1, nspins - CALL dbcsr_add(ec_env%matrix_hz(ispin)%matrix, dbcsr_work(ispin)%matrix, & - 1.0_dp, -1.0_dp) - END DO - CALL dbcsr_deallocate_matrix_set(dbcsr_work) - - END IF ! do_admm - - END IF ! do_fhx + ! Need to also subtract HFX contribution of reference calculation from ec_env%matrix_hz + ! and/or add HFX contribution if DC-DFT ueses hybrid XC-functional + CALL add_exx_to_rhs(rhs=ec_env%matrix_hz, & + qs_env=qs_env, & + ext_hfx_section=ec_hfx_sections, & + x_data=ec_env%x_data, & + recalc_integrals=.FALSE., & + do_admm=ec_env%do_ec_admm, & + do_ec=.TRUE., & + do_exx=.FALSE., & + reuse_hfx=ec_env%reuse_hfx) ! Core overlap IF (debug_forces) fodeb(1:3) = force(1)%core_overlap(1:3, 1) @@ -1249,20 +1255,35 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_build_ks_matrix' CHARACTER(LEN=default_string_length) :: headline - INTEGER :: handle, ispin, nspins - LOGICAL :: calculate_forces, use_virial - REAL(dp) :: eexc, evhxc + INTEGER :: handle, iounit, ispin, nspins + LOGICAL :: calculate_forces, & + do_adiabatic_rescaling, do_ec_hfx, & + hfx_treat_lsd_in_core, use_virial + REAL(dp) :: dummy_real, dummy_real2(2), eexc, evhxc, & + t3 + TYPE(cp_logger_type), POINTER :: logger + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_mat 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 :: rho_r, tau_r, v_rspace, v_tau_rspace + TYPE(qs_energy_type), POINTER :: energy TYPE(qs_ks_env_type), POINTER :: ks_env TYPE(qs_rho_type), POINTER :: rho + TYPE(section_vals_type), POINTER :: adiabatic_rescaling_section, & + ec_hfx_sections, ec_section CALL timeset(routineN, handle) + logger => cp_get_default_logger() + IF (logger%para_env%is_source()) THEN + iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.) + ELSE + iounit = -1 + END IF + ! get all information on the electronic density - NULLIFY (auxbas_pw_pool, dft_control, ks_env, rho, rho_r, tau_r) + NULLIFY (auxbas_pw_pool, dft_control, energy, ks_env, rho, rho_r, tau_r) CALL get_qs_env(qs_env=qs_env, & dft_control=dft_control, & ks_env=ks_env, & @@ -1287,6 +1308,78 @@ CONTAINS CALL get_qs_env(qs_env=qs_env, pw_env=pw_env) CPASSERT(ASSOCIATED(pw_env)) + ! Exact exchange contribution (hybrid functionals) + ec_section => section_vals_get_subs_vals(qs_env%input, "DFT%ENERGY_CORRECTION") + ec_hfx_sections => section_vals_get_subs_vals(ec_section, "XC%HF") + CALL section_vals_get(ec_hfx_sections, explicit=do_ec_hfx) + + IF (do_ec_hfx) THEN + + ! Check what works + adiabatic_rescaling_section => section_vals_get_subs_vals(ec_section, "XC%ADIABATIC_RESCALING") + CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling) + IF (do_adiabatic_rescaling) THEN + CALL cp_abort(__LOCATION__, "Adiabatic rescaling NYI for energy correction") + END IF + CALL section_vals_val_get(ec_hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core) + IF (hfx_treat_lsd_in_core) THEN + CALL cp_abort(__LOCATION__, "HFX_TREAT_LSD_IN_CORE NYI for energy correction") + END IF + + ! calculate the density matrix for the fitted mo_coeffs + IF (dft_control%do_admm) THEN + + IF (dft_control%do_admm_mo) THEN + IF (qs_env%run_rtp) THEN + CALL rtp_admm_calc_rho_aux(qs_env) + ELSE + CALL admm_mo_calc_rho_aux(qs_env) + END IF + ELSEIF (dft_control%do_admm_dm) THEN + CALL admm_dm_calc_rho_aux(qs_env) + END IF + END IF + + ! Get exact exchange energy + dummy_real = 0.0_dp + t3 = 0.0_dp + CALL get_qs_env(qs_env, energy=energy) + CALL calculate_exx(qs_env=qs_env, & + unit_nr=iounit, & + hfx_sections=ec_hfx_sections, & + x_data=ec_env%x_data, & + do_gw=.FALSE., & + do_admm=ec_env%do_ec_admm, & + calc_forces=.FALSE., & + reuse_hfx=ec_env%reuse_hfx, & + do_im_time=.FALSE., & + E_ex_from_GW=dummy_real, & + E_admm_from_GW=dummy_real2, & + t3=dummy_real) + + ! Save exchange energy + ec_env%ex = energy%ex + ! Save EXX ADMM XC correction + IF (ec_env%do_ec_admm) THEN + ec_env%exc_aux_fit = energy%exc_aux_fit + energy%exc + END IF + + ! Add exact echange contribution of EC to EC Hamiltonian + ! do_ec = .FALSE prevents subtraction of HFX contribution of reference calculation + ! do_exx = .FALSE. prevents subtraction of reference XC contribution + ks_mat => ec_env%matrix_ks(:, 1) + CALL add_exx_to_rhs(rhs=ks_mat, & + qs_env=qs_env, & + ext_hfx_section=ec_hfx_sections, & + x_data=ec_env%x_data, & + recalc_integrals=.FALSE., & + do_admm=ec_env%do_ec_admm, & + do_ec=.FALSE., & + do_exx=.FALSE., & + reuse_hfx=ec_env%reuse_hfx) + + END IF + ! v_rspace and v_tau_rspace are generated from the auxbas pool CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) NULLIFY (v_rspace, v_tau_rspace) @@ -1599,16 +1692,17 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_build_ks_matrix_force' INTEGER :: handle, i, iounit, ispin, natom, nspins - LOGICAL :: debug_forces, debug_stress, use_virial - REAL(dp) :: dehartree, eexc, ehartree, eovrl, exc, & - fconv + LOGICAL :: debug_forces, debug_stress, do_ec_hfx, & + use_virial + REAL(dp) :: dehartree, dummy_real, dummy_real2(2), & + eexc, ehartree, eovrl, exc, fconv REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ftot REAL(dp), DIMENSION(3) :: fodeb REAL(KIND=dp), DIMENSION(3, 3) :: h_stress, pv_loc, stdeb, sttot TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(cell_type), POINTER :: cell TYPE(cp_logger_type), POINTER :: logger - TYPE(dbcsr_p_type) :: scrm + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, scrm TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s TYPE(dft_control_type), POINTER :: dft_control TYPE(mp_para_env_type), POINTER :: para_env @@ -1627,7 +1721,7 @@ CONTAINS TYPE(qs_force_type), DIMENSION(:), POINTER :: force TYPE(qs_ks_env_type), POINTER :: ks_env TYPE(qs_rho_type), POINTER :: rho - TYPE(section_vals_type), POINTER :: xc_section + TYPE(section_vals_type), POINTER :: ec_hfx_sections, xc_section TYPE(virial_type), POINTER :: virial CALL timeset(routineN, handle) @@ -1644,13 +1738,14 @@ CONTAINS ! get all information on the electronic density NULLIFY (atomic_kind_set, cell, dft_control, force, ks_env, & - matrix_p, matrix_s, para_env, rho, rho_core, rho_g, rho_r, sab_orb, & - tau_r, virial) + matrix_ks, matrix_p, matrix_s, para_env, rho, rho_core, & + rho_g, rho_r, sab_orb, tau_r, virial) CALL get_qs_env(qs_env=qs_env, & cell=cell, & dft_control=dft_control, & force=force, & ks_env=ks_env, & + matrix_ks=matrix_ks, & para_env=para_env, & rho=rho, & sab_orb=sab_orb, & @@ -1972,11 +2067,15 @@ CONTAINS virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc) END IF - ! initialize src matrix - NULLIFY (scrm%matrix) - ALLOCATE (scrm%matrix) - CALL dbcsr_create(scrm%matrix, template=ec_env%matrix_s(1, 1)%matrix) - CALL cp_dbcsr_alloc_block_from_nbl(scrm%matrix, ec_env%sab_orb) + ! initialize srcm matrix + NULLIFY (scrm) + CALL dbcsr_allocate_matrix_set(scrm, nspins) + DO ispin = 1, nspins + ALLOCATE (scrm(ispin)%matrix) + CALL dbcsr_create(scrm(ispin)%matrix, template=ec_env%matrix_ks(ispin, 1)%matrix) + CALL dbcsr_copy(scrm(ispin)%matrix, ec_env%matrix_ks(ispin, 1)%matrix) + CALL dbcsr_set(scrm(ispin)%matrix, 0.0_dp) + END DO ! v_rspace and v_tau_rspace are generated from the auxbas pool NULLIFY (v_rspace, v_tau_rspace) @@ -2023,7 +2122,7 @@ CONTAINS CALL pw_axpy(v_hartree_rspace, v_rspace(ispin)) ! integrate over potential CALL integrate_v_rspace(v_rspace=v_rspace(ispin), & - hmat=scrm, & + hmat=scrm(ispin), & pmat=ec_env%matrix_p(ispin, 1), & qs_env=qs_env, & calculate_forces=.TRUE., & @@ -2054,7 +2153,7 @@ CONTAINS ! integrate over Tau-potential CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol) CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), & - hmat=scrm, & + hmat=scrm(ispin), & pmat=ec_env%matrix_p(ispin, 1), & qs_env=qs_env, & calculate_forces=.TRUE., & @@ -2069,9 +2168,55 @@ CONTAINS END IF END IF - ! delete scr matrix - CALL dbcsr_release(scrm%matrix) - DEALLOCATE (scrm%matrix) +!------------------------------------------------------------------------------ +! HFX direct force +!------------------------------------------------------------------------------ + + ! If hybrid functional + ec_hfx_sections => section_vals_get_subs_vals(qs_env%input, "DFT%ENERGY_CORRECTION%XC%HF") + CALL section_vals_get(ec_hfx_sections, explicit=do_ec_hfx) + + IF (do_ec_hfx) THEN + + IF (debug_forces) fodeb(1:3) = force(1)%fock_4c(1:3, 1) + IF (use_virial) virial%pv_fock_4c = 0.0_dp + + CALL calculate_exx(qs_env=qs_env, & + unit_nr=iounit, & + hfx_sections=ec_hfx_sections, & + x_data=ec_env%x_data, & + do_gw=.FALSE., & + do_admm=ec_env%do_ec_admm, & + calc_forces=.TRUE., & + reuse_hfx=ec_env%reuse_hfx, & + do_im_time=.FALSE., & + E_ex_from_GW=dummy_real, & + E_admm_from_GW=dummy_real2, & + t3=dummy_real) + + IF (use_virial) THEN + virial%pv_exx = virial%pv_exx - virial%pv_fock_4c + virial%pv_virial = virial%pv_virial - virial%pv_fock_4c + virial%pv_calculate = .FALSE. + END IF + IF (debug_forces) THEN + fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3) + CALL para_env%sum(fodeb) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*hfx ", fodeb + END IF + IF (debug_stress .AND. use_virial) THEN + stdeb = -1.0_dp*fconv*virial%pv_fock_4c + CALL para_env%sum(stdeb) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| Pout*hfx ', one_third_sum_diag(stdeb), det_3x3(stdeb) + END IF + + END IF + +!------------------------------------------------------------------------------ + + ! delete scrm matrix + CALL dbcsr_deallocate_matrix_set(scrm) ! return pw grids CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_rspace) @@ -2282,6 +2427,57 @@ CONTAINS CPASSERT(.FALSE.) END SELECT + ! OUtput density available now + +! HFX contribution to Harris functional and energy +! Can't calculate this earlier, cause ec_env%matrix_p doesnt exist yet +!------------------------------------------------------------------------------ + +! ! Exact exchange contribution (hybrid functionals) +! ec_section => section_vals_get_subs_vals(qs_env%input, "DFT%ENERGY_CORRECTION") +! ec_hfx_sections => section_vals_get_subs_vals(ec_section, "XC%HF") +! CALL section_vals_get(ec_hfx_sections, explicit=do_ec_hfx) +! +! IF (do_ec_hfx) THEN +! +! ! Check what works +! IF (dft_control%do_admm) THEN +! CALL cp_warn(__LOCATION__, "Energy correction with hybrid functional does not use ADMM.") +! END IF +! +! adiabatic_rescaling_section => section_vals_get_subs_vals(ec_section, "XC%ADIABATIC_RESCALING") +! CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling) +! IF (do_adiabatic_rescaling) THEN +! CALL cp_abort(__LOCATION__, "Adiabatic rescaling NYI for energy correction") +! END IF +! CALL section_vals_val_get(ec_hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core) +! IF (hfx_treat_lsd_in_core) THEN +! CALL cp_abort(__LOCATION__, "HFX_TREAT_LSD_IN_CORE NYI for energy correction") +! END IF +! +! ! Exchange matrix +! IF (ASSOCIATED(ec_env%matrix_x)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_x) +! CALL dbcsr_allocate_matrix_set(ec_env%matrix_x, nspins) +! DO ispin = 1, nspins +! headline = "EXCHANGE MATRIX" +! ALLOCATE (ec_env%matrix_x(ispin)%matrix) +! CALL dbcsr_create(ec_env%matrix_x(ispin)%matrix, name=TRIM(headline), & +! template=ec_env%matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric) +! CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_x(ispin)%matrix, ec_env%sab_orb) +! CALL dbcsr_set(ec_env%matrix_x(ispin)%matrix, 0.0_dp) +! END DO +! +! ! Get exact exchange energy (fraction) and its contribution to the EC hamiltonian +! should_update=.TRUE. +! ks_mat => ec_env%matrix_ks(:,1) +! CALL ec_hfx_contributions(qs_env, ks_mat, matrix_p, & +! ec_hfx_sections, ec_env%x_data, use_virial, & +! should_update, calculate_forces, matrix_x = ec_env%matrix_x, ex = ec_env%ex) +! +! END IF + +!------------------------------------------------------------------------------ + IF (ec_env%mao) THEN CALL mao_release_matrices(ec_env, ksmat, smat, pmat, wmat) END IF @@ -2572,11 +2768,12 @@ CONTAINS ec_env%eband = eband + ec_env%efield_nuclear ! Add Harris functional "correction" terms - ec_env%etotal = ec_env%eband + ec_env%ehartree + ec_env%exc - ec_env%vhxc + ec_env%edispersion + ec_env%etotal = ec_env%eband + ec_env%ehartree + ec_env%exc - ec_env%vhxc + ec_env%edispersion - ec_env%ex IF (unit_nr > 0) THEN WRITE (unit_nr, '(T3,A,T56,F25.15)') "Eband ", ec_env%eband WRITE (unit_nr, '(T3,A,T56,F25.15)') "Ehartree ", ec_env%ehartree WRITE (unit_nr, '(T3,A,T56,F25.15)') "Exc ", ec_env%exc + WRITE (unit_nr, '(T3,A,T56,F25.15)') "Ex ", ec_env%ex WRITE (unit_nr, '(T3,A,T56,F25.15)') "Evhxc ", ec_env%vhxc WRITE (unit_nr, '(T3,A,T56,F25.15)') "Edisp ", ec_env%edispersion WRITE (unit_nr, '(T3,A,T56,F25.15)') "Etotal Harris Functional ", ec_env%etotal @@ -2588,11 +2785,15 @@ CONTAINS CALL calculate_ptrace(ec_env%matrix_h, ec_env%matrix_p, ec_env%ecore, SIZE(ec_env%matrix_p, 1)) ec_env%ecore = ec_env%ecore + ec_env%efield_nuclear - ec_env%etotal = ec_env%ecore + ec_env%ehartree + ec_env%exc + ec_env%edispersion + ec_env%etotal = ec_env%ecore + ec_env%ehartree + ec_env%exc + ec_env%edispersion & + + ec_env%ex + ec_env%exc_aux_fit + IF (unit_nr > 0) THEN WRITE (unit_nr, '(T3,A,T56,F25.15)') "Ecore ", ec_env%ecore WRITE (unit_nr, '(T3,A,T56,F25.15)') "Ehartree ", ec_env%ehartree WRITE (unit_nr, '(T3,A,T56,F25.15)') "Exc ", ec_env%exc + WRITE (unit_nr, '(T3,A,T56,F25.15)') "Ex ", ec_env%ex + WRITE (unit_nr, '(T3,A,T56,F25.15)') "Exc_aux_fit", ec_env%exc_aux_fit WRITE (unit_nr, '(T3,A,T56,F25.15)') "Edisp ", ec_env%edispersion WRITE (unit_nr, '(T3,A,T56,F25.15)') "Etotal Energy Functional ", ec_env%etotal END IF diff --git a/src/rpa_hfx.F b/src/hfx_exx.F similarity index 70% rename from src/rpa_hfx.F rename to src/hfx_exx.F index fc321db596..0b24131411 100644 --- a/src/rpa_hfx.F +++ b/src/hfx_exx.F @@ -6,13 +6,14 @@ !--------------------------------------------------------------------------------------------------! ! ************************************************************************************************** -!> \brief Routines to calculate EXX in RPA +!> \brief Routines to calculate EXX in RPA and energy correction methods !> \par History !> 07.2020 separated from mp2.F [F. Stein, code by Jan Wilhelm] !> 06.2022 EXX contribution to the forces [A. Bussy] -!> \author Jan Wilhelm, Frederick Stein, Augustin Bussy +!> 03.2023 Generalized for energy correction methods +!> \author Jan Wilhelm, Frederick Stein, Augustin Bussy, Fabian Belleflamme ! ************************************************************************************************** -MODULE rpa_hfx +MODULE hfx_exx USE admm_methods, ONLY: admm_projection_derivative USE admm_types, ONLY: admm_env_create,& admm_env_release,& @@ -23,6 +24,9 @@ MODULE rpa_hfx copy_fm_to_dbcsr,& dbcsr_allocate_matrix_set,& dbcsr_deallocate_matrix_set + USE cp_log_handling, ONLY: cp_get_default_logger,& + cp_logger_get_default_unit_nr,& + cp_logger_type USE dbcsr_api, ONLY: dbcsr_add,& dbcsr_copy,& dbcsr_create,& @@ -37,6 +41,7 @@ MODULE rpa_hfx USE hfx_energy_potential, ONLY: integrate_four_center USE hfx_ri, ONLY: hfx_ri_update_forces,& hfx_ri_update_ks + USE hfx_types, ONLY: hfx_type USE input_constants, ONLY: do_admm_aux_exch_func_none USE input_section_types, ONLY: section_vals_create,& section_vals_duplicate,& @@ -75,9 +80,9 @@ MODULE rpa_hfx PRIVATE - CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_hfx' + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hfx_exx' - PUBLIC :: calculate_exx, add_exx_to_rhs, calc_ri_rpa_admm_xc_contributions, ri_rpa_pre_hfx, ri_rpa_post_hfx + PUBLIC :: calculate_exx, add_exx_to_rhs, calc_exx_admm_xc_contributions, exx_pre_hfx, exx_post_hfx CONTAINS @@ -85,17 +90,26 @@ CONTAINS !> \brief ... !> \param qs_env ... !> \param unit_nr ... +!> \param hfx_sections ... +!> \param x_data ... !> \param do_gw ... !> \param do_admm ... !> \param calc_forces ... +!> \param reuse_hfx ... +!> \param do_im_time ... !> \param E_ex_from_GW ... !> \param E_admm_from_GW ... !> \param t3 ... ! ************************************************************************************************** - SUBROUTINE calculate_exx(qs_env, unit_nr, do_gw, do_admm, calc_forces, E_ex_from_GW, E_admm_from_GW, t3) + SUBROUTINE calculate_exx(qs_env, unit_nr, hfx_sections, x_data, & + do_gw, do_admm, calc_forces, reuse_hfx, do_im_time, & + E_ex_from_GW, E_admm_from_GW, t3) TYPE(qs_environment_type), POINTER :: qs_env INTEGER, INTENT(IN) :: unit_nr - LOGICAL, INTENT(IN) :: do_gw, do_admm, calc_forces + TYPE(section_vals_type), POINTER :: hfx_sections + TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data + LOGICAL, INTENT(IN) :: do_gw, do_admm, calc_forces, reuse_hfx, & + do_im_time REAL(KIND=dp), INTENT(IN) :: E_ex_from_GW, E_admm_from_GW(2), t3 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_exx' @@ -104,7 +118,7 @@ CONTAINS nspins LOGICAL :: calc_ints, hfx_treat_lsd_in_core, & use_virial - REAL(KIND=dp) :: eh1, ehfx, t1, t2 + REAL(KIND=dp) :: eh1, ehfx, t1, t2, tf1, tf2 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, rho_ao, & rho_ao_aux_fit, rho_ao_resp TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_2d, rho_ao_2d @@ -112,17 +126,17 @@ CONTAINS TYPE(mp_para_env_type), POINTER :: para_env TYPE(qs_energy_type), POINTER :: energy TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit - TYPE(section_vals_type), POINTER :: hfx_sections, input + TYPE(section_vals_type), POINTER :: input TYPE(virial_type), POINTER :: virial CALL timeset(routineN, handle) t1 = m_walltime() - NULLIFY (hfx_sections, input, para_env, matrix_ks, matrix_ks_aux_fit, rho, rho_ao, virial, & + NULLIFY (input, para_env, matrix_ks, matrix_ks_aux_fit, rho, rho_ao, virial, & dft_control, rho_aux_fit, rho_ao_aux_fit) - CALL ri_rpa_pre_hfx(qs_env) + CALL exx_pre_hfx(hfx_sections, x_data, reuse_hfx) CALL get_qs_env(qs_env=qs_env, & input=input, & @@ -135,16 +149,16 @@ CONTAINS CALL qs_rho_get(rho, rho_ao=rho_ao) IF (do_admm) THEN - CALL get_admm_env(qs_env%admm_env, matrix_ks_aux_fit=matrix_ks_aux_fit, rho_aux_fit=rho_aux_fit) + CALL get_admm_env(qs_env%admm_env, & + matrix_ks_aux_fit=matrix_ks_aux_fit, & + rho_aux_fit=rho_aux_fit) CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux_fit) IF (qs_env%admm_env%do_gapw) THEN - CPABORT("RPA with ADMM EXX only implmented with GPW") + CPABORT("ADMM EXX only implmented with GPW") END IF END IF - hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%WF_CORRELATION%RI_RPA%HF") - CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf) CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, & i_rep_section=1) @@ -188,8 +202,8 @@ CONTAINS IF (hfx_treat_lsd_in_core) mspin = nspins calc_ints = .TRUE. - IF (qs_env%mp2_env%ri_rpa%reuse_hfx) calc_ints = .FALSE. - IF (calc_forces .AND. qs_env%mp2_env%do_im_time) calc_ints = .FALSE. + IF (reuse_hfx) calc_ints = .FALSE. + IF (calc_forces .AND. do_im_time) calc_ints = .FALSE. ehfx = 0.0_dp IF (do_admm) THEN @@ -202,15 +216,17 @@ CONTAINS DO irep = 1, n_rep_hf - IF (qs_env%mp2_env%ri_rpa%x_data(irep, 1)%do_hfx_ri) THEN - CALL hfx_ri_update_ks(qs_env, qs_env%mp2_env%ri_rpa%x_data(irep, 1)%ri_data, matrix_ks_2d, ehfx, & + IF (x_data(irep, 1)%do_hfx_ri) THEN + CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_2d, ehfx, & rho_ao=rho_ao_2d, geometry_did_change=calc_ints, nspins=nspins, & - hf_fraction=qs_env%mp2_env%ri_rpa%x_data(irep, 1)%general_parameter%fraction) + hf_fraction=x_data(irep, 1)%general_parameter%fraction) ELSE DO ispin = 1, mspin - CALL integrate_four_center(qs_env, qs_env%mp2_env%ri_rpa%x_data, matrix_ks_2d, eh1, & - rho_ao_2d, hfx_sections, para_env, calc_ints, irep, .TRUE., ispin=ispin) + + CALL integrate_four_center(qs_env, x_data, matrix_ks_2d, eh1, & + rho_ao_2d, hfx_sections, para_env, & + calc_ints, irep, .TRUE., ispin=ispin) ehfx = ehfx + eh1 END DO END IF @@ -227,25 +243,31 @@ CONTAINS END IF IF (calc_forces) THEN + tf1 = m_walltime() + !Note: no need to remove xc forces: they are not even calculated in the first place NULLIFY (rho_ao_resp) DO irep = 1, n_rep_hf - IF (qs_env%mp2_env%ri_rpa%x_data(irep, 1)%do_hfx_ri) THEN - - CALL hfx_ri_update_forces(qs_env, qs_env%mp2_env%ri_rpa%x_data(irep, 1)%ri_data, nspins, & - qs_env%mp2_env%ri_rpa%x_data(irep, 1)%general_parameter%fraction, & - rho_ao=rho_ao_2d, rho_ao_resp=rho_ao_resp, use_virial=use_virial) + IF (x_data(irep, 1)%do_hfx_ri) THEN + CALL hfx_ri_update_forces(qs_env, x_data(irep, 1)%ri_data, nspins, & + x_data(irep, 1)%general_parameter%fraction, & + rho_ao=rho_ao_2d, rho_ao_resp=rho_ao_resp, & + use_virial=use_virial) ELSE - CALL derivatives_four_center(qs_env, rho_ao_2d, rho_ao_resp, hfx_sections, para_env, irep, & - use_virial, external_x_data=qs_env%mp2_env%ri_rpa%x_data) + CALL derivatives_four_center(qs_env, rho_ao_2d, rho_ao_resp, & + hfx_sections, para_env, irep, & + use_virial, external_x_data=x_data) END IF END DO !irep + + tf2 = m_walltime() + IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.6)') 'Total EXX Force Time=', tf2 - tf1 END IF IF (use_virial) THEN @@ -256,8 +278,9 @@ CONTAINS ! ADMM XC correction IF (do_admm) THEN - CALL calc_ri_rpa_admm_xc_contributions(qs_env, matrix_ks, matrix_ks_aux_fit, energy%exc, & - energy%exc_aux_fit, calc_forces, use_virial) + CALL calc_exx_admm_xc_contributions(qs_env, matrix_ks, matrix_ks_aux_fit, x_data, & + energy%exc, energy%exc_aux_fit, calc_forces, & + use_virial) ! ADMM overlap forces IF (calc_forces) CALL admm_projection_derivative(qs_env, matrix_ks_aux_fit, rho_ao) @@ -280,7 +303,7 @@ CONTAINS END IF END IF - CALL ri_rpa_post_hfx(qs_env) + CALL exx_post_hfx(qs_env, x_data, reuse_hfx) CALL timestop(handle) @@ -290,20 +313,33 @@ CONTAINS !> \brief Add the EXX contribution to the RHS of the Z-vector equation, namely the HF Hamiltonian !> \param rhs ... !> \param qs_env ... +!> \param ext_hfx_section ... +!> \param x_data ... !> \param recalc_integrals ... +!> \param do_admm ... +!> \param do_ec ... +!> \param do_exx ... +!> \param reuse_hfx ... ! ************************************************************************************************** - SUBROUTINE add_exx_to_rhs(rhs, qs_env, recalc_integrals) + SUBROUTINE add_exx_to_rhs(rhs, qs_env, ext_hfx_section, x_data, & + recalc_integrals, do_admm, do_ec, do_exx, reuse_hfx) TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: rhs TYPE(qs_environment_type), POINTER :: qs_env - LOGICAL, INTENT(IN), OPTIONAL :: recalc_integrals + TYPE(section_vals_type), POINTER :: ext_hfx_section + TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data + LOGICAL, INTENT(IN), OPTIONAL :: recalc_integrals, do_admm, do_ec, & + do_exx, reuse_hfx CHARACTER(LEN=*), PARAMETER :: routineN = 'add_exx_to_rhs' - INTEGER :: handle, ispin, nao, nao_aux, nspins - LOGICAL :: calc_ints, do_hfx, my_recalc_integrals + INTEGER :: handle, ispin, nao, nao_aux, nspins, & + unit_nr + LOGICAL :: calc_ints, do_hfx, my_do_ec, my_do_exx, & + my_recalc_integrals REAL(dp) :: dummy_real1, dummy_real2 TYPE(admm_type), POINTER :: admm_env + TYPE(cp_logger_type), POINTER :: logger TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dbcsr_work, matrix_ks, matrix_s_aux, & rho_ao, rho_ao_aux, work_admm TYPE(dbcsr_type) :: dbcsr_tmp @@ -320,20 +356,36 @@ CONTAINS auxbas_pw_pool, dft_control, rho_ao, rho_aux_fit, rho_ao_aux, work_admm, & matrix_s_aux, admm_env, task_list_aux_fit) + logger => cp_get_default_logger() + IF (logger%para_env%is_source()) THEN + unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) + ELSE + unit_nr = -1 + END IF + CALL timeset(routineN, handle) my_recalc_integrals = .FALSE. IF (PRESENT(recalc_integrals)) my_recalc_integrals = recalc_integrals + my_do_ec = .FALSE. + IF (PRESENT(do_ec)) my_do_ec = do_ec + my_do_exx = .TRUE. + IF (PRESENT(do_exx)) my_do_exx = do_exx ! Strategy: we take the ks_matrix, remove the current xc contribution, and then add the RPA HF one CALL get_qs_env(qs_env, matrix_ks=matrix_ks, rho=rho, pw_env=pw_env, dft_control=dft_control) nspins = dft_control%nspins + ! do_exx ; subtract XC and EX + ! do_dcdft: subtract EX + CALL dbcsr_allocate_matrix_set(dbcsr_work, nspins) DO ispin = 1, nspins ALLOCATE (dbcsr_work(ispin)%matrix) CALL dbcsr_copy(dbcsr_work(ispin)%matrix, matrix_ks(ispin)%matrix) + CALL dbcsr_set(dbcsr_work(ispin)%matrix, 0.0_dp) END DO + IF (dft_control%do_admm) THEN CALL get_qs_env(qs_env, admm_env=admm_env) CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux, task_list_aux_fit=task_list_aux_fit, & @@ -353,55 +405,74 @@ CONTAINS CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux) END IF - !Remove the standard XC + HFX contribution CALL qs_rho_get(rho, rho_ao=rho_ao) + + ! Remove the standard XC + HFX contribution CALL ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, dummy_real1, dummy_real2) - DO ispin = 1, nspins - CALL pw_scale(vxc_rspace(ispin), -1.0_dp) - CALL integrate_v_rspace(v_rspace=vxc_rspace(ispin), hmat=dbcsr_work(ispin), qs_env=qs_env, & - calculate_forces=.FALSE.) - IF (ASSOCIATED(vtau_rspace)) THEN - CALL pw_scale(vtau_rspace(ispin), -1.0_dp) - CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), hmat=dbcsr_work(ispin), qs_env=qs_env, & - calculate_forces=.FALSE., compute_tau=.TRUE.) - END IF - IF (dft_control%do_admm) THEN - !note: factor -1.0 taken care of later, after HFX ADMM contribution is taken - IF (.NOT. qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN - CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), hmat=work_admm(ispin), & - qs_env=qs_env, calculate_forces=.FALSE., basis_type="AUX_FIT", & - task_list_external=task_list_aux_fit) + + ! Only remove standard XC and/or HFX + ! (due to different requirements for RHS) + IF (my_do_exx) THEN + + DO ispin = 1, nspins + CALL dbcsr_copy(dbcsr_work(ispin)%matrix, matrix_ks(ispin)%matrix) + END DO + + DO ispin = 1, nspins + CALL pw_scale(vxc_rspace(ispin), -1.0_dp) + CALL integrate_v_rspace(v_rspace=vxc_rspace(ispin), hmat=dbcsr_work(ispin), qs_env=qs_env, & + calculate_forces=.FALSE.) + IF (ASSOCIATED(vtau_rspace)) THEN + CALL pw_scale(vtau_rspace(ispin), -1.0_dp) + CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), hmat=dbcsr_work(ispin), qs_env=qs_env, & + calculate_forces=.FALSE., compute_tau=.TRUE.) END IF - END IF - END DO + END DO - hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF") - CALL section_vals_get(hfx_section, explicit=do_hfx) - IF (do_hfx) THEN - IF (dft_control%do_admm) THEN + END IF ! do_exx - CALL tddft_hfx_matrix(work_admm, rho_ao_aux, qs_env, .FALSE., my_recalc_integrals) + ! Remove standard HFX + IF (my_do_exx .OR. my_do_ec) THEN - DO ispin = 1, nspins - CALL copy_dbcsr_to_fm(work_admm(ispin)%matrix, admm_env%work_aux_aux) - CALL parallel_gemm('N', 'N', nao_aux, nao, nao_aux, 1.0_dp, admm_env%work_aux_aux, admm_env%A, & - 0.0_dp, admm_env%work_aux_orb) - CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, admm_env%work_aux_orb, & - 0.0_dp, admm_env%work_orb_orb) - CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbcsr_tmp, keep_sparsity=.TRUE.) - CALL dbcsr_add(dbcsr_work(ispin)%matrix, dbcsr_tmp, 1.0_dp, -1.0_dp) - END DO - ELSE - DO ispin = 1, nspins - CALL dbcsr_scale(rho_ao(ispin)%matrix, -1.0_dp) - END DO - CALL tddft_hfx_matrix(dbcsr_work, rho_ao, qs_env, .FALSE., my_recalc_integrals) - DO ispin = 1, nspins - CALL dbcsr_scale(rho_ao(ispin)%matrix, -1.0_dp) - END DO - END IF - END IF !do_hfx + hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF") + CALL section_vals_get(hfx_section, explicit=do_hfx) + IF (do_hfx) THEN + IF (dft_control%do_admm) THEN + !note: factor -1.0 taken care of later, after HFX ADMM contribution is taken + IF (.NOT. qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN + DO ispin = 1, nspins + CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), hmat=work_admm(ispin), & + qs_env=qs_env, calculate_forces=.FALSE., basis_type="AUX_FIT", & + task_list_external=task_list_aux_fit) + END DO + END IF + + CALL tddft_hfx_matrix(work_admm, rho_ao_aux, qs_env, .FALSE., my_recalc_integrals) + + DO ispin = 1, nspins + CALL copy_dbcsr_to_fm(work_admm(ispin)%matrix, admm_env%work_aux_aux) + CALL parallel_gemm('N', 'N', nao_aux, nao, nao_aux, 1.0_dp, admm_env%work_aux_aux, admm_env%A, & + 0.0_dp, admm_env%work_aux_orb) + CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, admm_env%work_aux_orb, & + 0.0_dp, admm_env%work_orb_orb) + CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbcsr_tmp, keep_sparsity=.TRUE.) + CALL dbcsr_add(dbcsr_work(ispin)%matrix, dbcsr_tmp, 1.0_dp, -1.0_dp) + END DO + ELSE + DO ispin = 1, nspins + CALL dbcsr_scale(rho_ao(ispin)%matrix, -1.0_dp) + END DO + CALL tddft_hfx_matrix(dbcsr_work, rho_ao, qs_env, .FALSE., my_recalc_integrals) + DO ispin = 1, nspins + CALL dbcsr_scale(rho_ao(ispin)%matrix, -1.0_dp) + END DO + END IF + END IF !do_hfx + + END IF ! do_exx + + ! Clean CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) CALL pw_pool_give_back_pw(auxbas_pw_pool, vh_rspace) DO ispin = 1, nspins @@ -409,30 +480,35 @@ CONTAINS IF (ASSOCIATED(vtau_rspace)) THEN CALL pw_pool_give_back_pw(auxbas_pw_pool, vtau_rspace(ispin)) END IF - IF (ASSOCIATED(vadmm_rspace)) THEN - CALL pw_pool_give_back_pw(auxbas_pw_pool, vadmm_rspace(ispin)) - END IF END DO DEALLOCATE (vxc_rspace) IF (ASSOCIATED(vtau_rspace)) DEALLOCATE (vtau_rspace) - IF (ASSOCIATED(vadmm_rspace)) DEALLOCATE (vadmm_rspace) - !Add the HF contribution from RI_RPA to the ks matrix - CALL ri_rpa_pre_hfx(qs_env) + IF (dft_control%do_admm) THEN + DO ispin = 1, nspins + IF (ASSOCIATED(vadmm_rspace)) THEN + CALL pw_pool_give_back_pw(auxbas_pw_pool, vadmm_rspace(ispin)) + END IF + END DO + IF (ASSOCIATED(vadmm_rspace)) DEALLOCATE (vadmm_rspace) + END IF + + !Add the HF contribution from RI_RPA/EC_ENV to the ks matrix + CALL exx_pre_hfx(ext_hfx_section, x_data, reuse_hfx) + calc_ints = .TRUE. - IF (qs_env%mp2_env%ri_rpa%reuse_hfx) calc_ints = .FALSE. - hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF") - IF (qs_env%mp2_env%ri_rpa%do_admm) THEN + IF (reuse_hfx) calc_ints = .FALSE. + IF (do_admm) THEN DO ispin = 1, nspins CALL dbcsr_set(work_admm(ispin)%matrix, 0.0_dp) END DO - CALL tddft_hfx_matrix(work_admm, rho_ao_aux, qs_env, .FALSE., calc_ints, hfx_section, & - qs_env%mp2_env%ri_rpa%x_data) + CALL tddft_hfx_matrix(work_admm, rho_ao_aux, qs_env, .FALSE., calc_ints, ext_hfx_section, & + x_data) !ADMM XC correction - CALL calc_ri_rpa_admm_xc_contributions(qs_env, dbcsr_work, work_admm, dummy_real1, & - dummy_real2, .FALSE., .FALSE.) + CALL calc_exx_admm_xc_contributions(qs_env, dbcsr_work, work_admm, x_data, dummy_real1, & + dummy_real2, .FALSE., .FALSE.) DO ispin = 1, nspins CALL copy_dbcsr_to_fm(work_admm(ispin)%matrix, admm_env%work_aux_aux) @@ -445,9 +521,9 @@ CONTAINS END DO ELSE - CALL tddft_hfx_matrix(dbcsr_work, rho_ao, qs_env, .FALSE., calc_ints, hfx_section, qs_env%mp2_env%ri_rpa%x_data) + CALL tddft_hfx_matrix(dbcsr_work, rho_ao, qs_env, .FALSE., calc_ints, ext_hfx_section, x_data) END IF - CALL ri_rpa_post_hfx(qs_env) + CALL exx_post_hfx(qs_env, x_data, reuse_hfx) !Update the RHS DO ispin = 1, nspins @@ -465,13 +541,15 @@ CONTAINS END SUBROUTINE add_exx_to_rhs ! ************************************************************************************************** -!> \brief get the ADMM XC section from the ri_rpa type if available, create and store them otherwise +!> \brief get the ADMM XC section from the ri_rpa/ec_env type if available, create and store them otherwise !> \param qs_env ... +!> \param x_data ... !> \param xc_section_aux ... !> \param xc_section_primary ... ! ************************************************************************************************** - SUBROUTINE get_ri_rpa_admm_xc_sections(qs_env, xc_section_aux, xc_section_primary) + SUBROUTINE get_exx_admm_xc_sections(qs_env, x_data, xc_section_aux, xc_section_primary) TYPE(qs_environment_type), POINTER :: qs_env + TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data TYPE(section_vals_type), POINTER :: xc_section_aux, xc_section_primary INTEGER :: natom @@ -484,11 +562,21 @@ CONTAINS NULLIFY (qs_admm_env, tmp_admm_env, para_env, xc_section, xc_section_empty, xc_fun_empty, & xc_fun, dft_control) - IF (ASSOCIATED(qs_env%mp2_env%ri_rpa%xc_section_aux) .AND. & - ASSOCIATED(qs_env%mp2_env%ri_rpa%xc_section_primary)) THEN - xc_section_aux => qs_env%mp2_env%ri_rpa%xc_section_aux - xc_section_primary => qs_env%mp2_env%ri_rpa%xc_section_primary - ELSE + IF (ASSOCIATED(qs_env%mp2_env)) THEN + IF (ASSOCIATED(qs_env%mp2_env%ri_rpa%xc_section_aux) .AND. & + ASSOCIATED(qs_env%mp2_env%ri_rpa%xc_section_primary)) THEN + xc_section_aux => qs_env%mp2_env%ri_rpa%xc_section_aux + xc_section_primary => qs_env%mp2_env%ri_rpa%xc_section_primary + END IF + ELSEIF (qs_env%energy_correction) THEN + IF (ASSOCIATED(qs_env%ec_env%xc_section_aux) .AND. & + ASSOCIATED(qs_env%ec_env%xc_section_primary)) THEN + xc_section_aux => qs_env%ec_env%xc_section_aux + xc_section_primary => qs_env%ec_env%xc_section_primary + END IF + END IF + + IF (.NOT. ASSOCIATED(xc_section_aux) .OR. .NOT. ASSOCIATED(xc_section_primary)) THEN CALL get_qs_env(qs_env, admm_env=qs_admm_env, natom=natom, para_env=para_env, dft_control=dft_control) CPASSERT(ASSOCIATED(qs_admm_env)) @@ -503,14 +591,20 @@ CONTAINS CALL admm_env_create(tmp_admm_env, dft_control%admm_control, qs_admm_env%mos_aux_fit, & para_env, natom, qs_admm_env%nao_aux_fit) - CALL create_admm_xc_section(x_data=qs_env%mp2_env%ri_rpa%x_data, xc_section=xc_section_empty, & + CALL create_admm_xc_section(x_data=x_data, xc_section=xc_section_empty, & admm_env=tmp_admm_env) CALL section_vals_duplicate(tmp_admm_env%xc_section_aux, xc_section_aux) CALL section_vals_duplicate(tmp_admm_env%xc_section_primary, xc_section_primary) - qs_env%mp2_env%ri_rpa%xc_section_aux => xc_section_aux - qs_env%mp2_env%ri_rpa%xc_section_primary => xc_section_primary + IF (ASSOCIATED(qs_env%mp2_env)) THEN + qs_env%mp2_env%ri_rpa%xc_section_aux => xc_section_aux + qs_env%mp2_env%ri_rpa%xc_section_primary => xc_section_primary + END IF + IF (qs_env%energy_correction) THEN + qs_env%ec_env%xc_section_aux => xc_section_aux + qs_env%ec_env%xc_section_primary => xc_section_primary + END IF CALL section_vals_release(xc_section_empty) CALL section_vals_release(xc_fun_empty) @@ -518,22 +612,24 @@ CONTAINS END IF - END SUBROUTINE get_ri_rpa_admm_xc_sections + END SUBROUTINE get_exx_admm_xc_sections ! ************************************************************************************************** -!> \brief Calculate the RI_RPA%HF ADMM XC contributions to the KS matrices and the respective energies +!> \brief Calculate the RI_RPA%HF / EC_ENV%HF ADMM XC contributions to the KS matrices and the respective energies !> \param qs_env ... !> \param matrix_prim ... !> \param matrix_aux ... +!> \param x_data ... !> \param exc ... !> \param exc_aux_fit ... !> \param calc_forces ... !> \param use_virial ... ! ************************************************************************************************** - SUBROUTINE calc_ri_rpa_admm_xc_contributions(qs_env, matrix_prim, matrix_aux, exc, exc_aux_fit, & - calc_forces, use_virial) + SUBROUTINE calc_exx_admm_xc_contributions(qs_env, matrix_prim, matrix_aux, x_data, exc, exc_aux_fit, & + calc_forces, use_virial) TYPE(qs_environment_type), POINTER :: qs_env TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_prim, matrix_aux + TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data REAL(dp), INTENT(INOUT) :: exc, exc_aux_fit LOGICAL, INTENT(IN) :: calc_forces, use_virial @@ -560,7 +656,7 @@ CONTAINS CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux_fit) CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) - CALL get_ri_rpa_admm_xc_sections(qs_env, xc_section_aux, xc_section_primary) + CALL get_exx_admm_xc_sections(qs_env, x_data, xc_section_aux, xc_section_primary) IF (use_virial) virial%pv_xc = 0.0_dp CALL qs_vxc_create(qs_env%ks_env, rho_struct=rho_aux_fit, xc_section=xc_section_aux, & @@ -627,45 +723,50 @@ CONTAINS DEALLOCATE (v_dummy) END IF - END SUBROUTINE calc_ri_rpa_admm_xc_contributions + END SUBROUTINE calc_exx_admm_xc_contributions ! ************************************************************************************************** -!> \brief Prepare the ri_rpa%x_data for integration. Simply change the HFX fraction in case the +!> \brief Prepare the external x_data for integration. Simply change the HFX fraction in case the !> qs_env%x_data is reused -!> \param qs_env ... +!> \param ext_hfx_section ... +!> \param x_data ... +!> \param reuse_hfx ... ! ************************************************************************************************** - SUBROUTINE ri_rpa_pre_hfx(qs_env) - TYPE(qs_environment_type), POINTER :: qs_env + SUBROUTINE exx_pre_hfx(ext_hfx_section, x_data, reuse_hfx) + TYPE(section_vals_type), POINTER :: ext_hfx_section + TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data + LOGICAL :: reuse_hfx INTEGER :: irep, n_rep_hf REAL(dp) :: frac - TYPE(section_vals_type), POINTER :: input, rpa_hfx_section - IF (.NOT. qs_env%mp2_env%ri_rpa%reuse_hfx) RETURN + IF (.NOT. reuse_hfx) RETURN - CALL get_qs_env(qs_env, input=input) - rpa_hfx_section => section_vals_get_subs_vals(input, "DFT%XC%WF_CORRELATION%RI_RPA%HF") - CALL section_vals_get(rpa_hfx_section, n_repetition=n_rep_hf) + CALL section_vals_get(ext_hfx_section, n_repetition=n_rep_hf) DO irep = 1, n_rep_hf - CALL section_vals_val_get(rpa_hfx_section, "FRACTION", r_val=frac, i_rep_section=irep) - qs_env%mp2_env%ri_rpa%x_data(irep, :)%general_parameter%fraction = frac + CALL section_vals_val_get(ext_hfx_section, "FRACTION", r_val=frac, i_rep_section=irep) + x_data(irep, :)%general_parameter%fraction = frac END DO - END SUBROUTINE ri_rpa_pre_hfx + END SUBROUTINE exx_pre_hfx ! ************************************************************************************************** !> \brief Revert back to the proper HFX fraction in case qs_env%x_data is reused !> \param qs_env ... +!> \param x_data ... +!> \param reuse_hfx ... ! ************************************************************************************************** - SUBROUTINE ri_rpa_post_hfx(qs_env) + SUBROUTINE exx_post_hfx(qs_env, x_data, reuse_hfx) TYPE(qs_environment_type), POINTER :: qs_env + TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data + LOGICAL :: reuse_hfx INTEGER :: irep, n_rep_hf REAL(dp) :: frac TYPE(section_vals_type), POINTER :: input, qs_hfx_section - IF (.NOT. qs_env%mp2_env%ri_rpa%reuse_hfx) RETURN + IF (.NOT. reuse_hfx) RETURN CALL get_qs_env(qs_env, input=input) qs_hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF") @@ -673,10 +774,10 @@ CONTAINS DO irep = 1, n_rep_hf CALL section_vals_val_get(qs_hfx_section, "FRACTION", r_val=frac, i_rep_section=irep) - qs_env%mp2_env%ri_rpa%x_data(irep, :)%general_parameter%fraction = frac + x_data(irep, :)%general_parameter%fraction = frac END DO - END SUBROUTINE ri_rpa_post_hfx + END SUBROUTINE exx_post_hfx ! ************************************************************************************************** !> \brief ... @@ -697,5 +798,5 @@ CONTAINS END SUBROUTINE -END MODULE rpa_hfx +END MODULE hfx_exx diff --git a/src/input_cp2k_ec.F b/src/input_cp2k_ec.F index 77e082cca7..191195abe8 100644 --- a/src/input_cp2k_ec.F +++ b/src/input_cp2k_ec.F @@ -347,6 +347,19 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create( & + keyword, __LOCATION__, & + name="ADMM", & + description="Decide whether to perform ADMM in the exact exchange calc. for DC-DFT. "// & + "The ADMM XC correction is governed by the AUXILIARY_DENSITY_MATRIX_METHOD section in &DFT. "// & + "In most cases, the Hartree-Fock exchange is not too expensive and there is no need for ADMM, "// & + "ADMM can however provide significant speedup and memory savings in case of diffuse basis sets. ", & + usage="ADMM", & + default_l_val=.FALSE., & + lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + END SUBROUTINE create_ec_section ! ************************************************************************************************** diff --git a/src/mp2.F b/src/mp2.F index 70ac971f0b..a0d07b8285 100644 --- a/src/mp2.F +++ b/src/mp2.F @@ -42,6 +42,7 @@ MODULE mp2 cp_print_key_unit_nr USE dbcsr_api, ONLY: dbcsr_get_info,& dbcsr_p_type + USE hfx_exx, ONLY: calculate_exx USE hfx_types, ONLY: & alloc_containers, dealloc_containers, hfx_basis_info_type, hfx_basis_type, & hfx_container_type, hfx_create_basis_types, hfx_init_container, hfx_release_basis_types, & @@ -85,7 +86,6 @@ MODULE mp2 eigensolver_symm USE qs_scf_types, ONLY: qs_scf_env_type USE rpa_gw_sigma_x, ONLY: compute_vec_Sigma_x_minus_vxc_gw - USE rpa_hfx, ONLY: calculate_exx USE scf_control_types, ONLY: scf_control_type USE virial_types, ONLY: virial_type @@ -121,7 +121,7 @@ CONTAINS INTEGER(KIND=int_8) :: mem INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, nelec LOGICAL :: calc_ex, do_admm, do_admm_rpa_exx, do_dynamic_load_balancing, do_exx, do_gw, & - do_kpoints_cubic_RPA, free_hfx_buffer, update_xc_energy + do_im_time, do_kpoints_cubic_RPA, free_hfx_buffer, reuse_hfx, update_xc_energy REAL(KIND=dp) :: E_admm_from_GW(2), E_ex_from_GW, Emp2, Emp2_AA, Emp2_AA_Cou, Emp2_AA_ex, & Emp2_AB, Emp2_AB_Cou, Emp2_AB_ex, Emp2_BB, Emp2_BB_Cou, Emp2_BB_ex, Emp2_Cou, Emp2_ex, & Emp2_S, Emp2_T, maxocc, mem_real, t1, t2, t3 @@ -726,7 +726,22 @@ IF (.NOT. (mp2_env%method == ri_mp2_method_gpw .OR. mp2_env%method == ri_rpa_met IF (do_exx) THEN do_gw = mp2_env%ri_rpa%do_ri_g0w0 do_admm = mp2_env%ri_rpa%do_admm - CALL calculate_exx(qs_env, unit_nr, do_gw, do_admm, .FALSE., E_ex_from_GW, E_admm_from_GW, t3) + reuse_hfx = qs_env%mp2_env%ri_rpa%reuse_hfx + do_im_time = qs_env%mp2_env%do_im_time + + CALL calculate_exx(qs_env=qs_env, & + unit_nr=unit_nr, & + hfx_sections=hfx_sections, & + x_data=qs_env%mp2_env%ri_rpa%x_data, & + do_gw=do_gw, & + do_admm=do_admm, & + calc_forces=.FALSE., & + reuse_hfx=reuse_hfx, & + do_im_time=do_im_time, & + E_ex_from_GW=E_ex_from_GW, & + E_admm_from_GW=E_admm_from_GW, & + t3=t3) + END IF END IF diff --git a/src/mp2_cphf.F b/src/mp2_cphf.F index 0f5242d0a8..d788c05e12 100644 --- a/src/mp2_cphf.F +++ b/src/mp2_cphf.F @@ -48,6 +48,7 @@ MODULE mp2_cphf dbcsr_set USE hfx_admm_utils, ONLY: tddft_hfx_matrix USE hfx_derivatives, ONLY: derivatives_four_center + USE hfx_exx, ONLY: add_exx_to_rhs USE hfx_ri, ONLY: hfx_ri_update_forces USE hfx_types, ONLY: alloc_containers,& hfx_container_type,& @@ -124,7 +125,6 @@ MODULE mp2_cphf qs_p_env_type USE qs_rho_types, ONLY: qs_rho_get,& qs_rho_type - USE rpa_hfx, ONLY: add_exx_to_rhs USE task_list_types, ONLY: task_list_type USE virial_types, ONLY: virial_type,& zero_virial @@ -329,7 +329,14 @@ CONTAINS CALL section_vals_get(hfx_section, explicit=do_exx) END IF IF (do_exx) THEN - CALL add_exx_to_rhs(P_mu_nu, qs_env, mp2_env%ri_grad%free_hfx_buffer) + CALL add_exx_to_rhs(rhs=P_mu_nu, & + qs_env=qs_env, & + ext_hfx_section=hfx_section, & + x_data=mp2_env%ri_rpa%x_data, & + recalc_integrals=.FALSE., & + do_admm=mp2_env%ri_rpa%do_admm, & + do_exx=do_exx, & + reuse_hfx=mp2_env%ri_rpa%reuse_hfx) focc = 1.0_dp IF (nspins == 1) focc = 2.0_dp diff --git a/src/qs_environment.F b/src/qs_environment.F index 013475a60b..a155ac2f06 100644 --- a/src/qs_environment.F +++ b/src/qs_environment.F @@ -64,7 +64,8 @@ MODULE qs_environment USE distribution_methods, ONLY: distribute_molecules_1d USE dm_ls_scf_create, ONLY: ls_scf_create USE ec_env_types, ONLY: energy_correction_type - USE ec_environment, ONLY: ec_env_create + USE ec_environment, ONLY: ec_env_create,& + ec_write_input USE et_coupling_types, ONLY: et_coupling_create USE ewald_environment_types, ONLY: ewald_env_create,& ewald_env_get,& @@ -250,7 +251,7 @@ CONTAINS LOGICAL, INTENT(IN) :: use_motion_section INTEGER :: ikind, method_id, nelectron_total, nkind - LOGICAL :: do_admm_rpa, do_et, do_exx, do_hfx, do_kpoints, is_identical, is_semi, & + LOGICAL :: do_admm_rpa, do_ec_hfx, do_et, do_exx, do_hfx, do_kpoints, is_identical, is_semi, & mp2_present, my_qmmm, qmmm_decoupl, same_except_frac, silent, use_ref_cell REAL(KIND=dp), DIMENSION(:, :), POINTER :: rtmat TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set @@ -268,8 +269,9 @@ CONTAINS TYPE(qs_wf_history_type), POINTER :: wf_history TYPE(rel_control_type), POINTER :: rel_control TYPE(scf_control_type), POINTER :: scf_control - TYPE(section_vals_type), POINTER :: dft_section, ec_section, et_coupling_section, & - hfx_section, kpoint_section, mp2_section, rpa_hfx_section, transport_section + TYPE(section_vals_type), POINTER :: dft_section, ec_hfx_section, ec_section, & + et_coupling_section, hfx_section, kpoint_section, mp2_section, rpa_hfx_section, & + transport_section NULLIFY (my_cell, my_cell_ref, atomic_kind_set, particle_set, & qs_kind_set, kpoint_section, dft_section, ec_section, & @@ -439,9 +441,6 @@ CONTAINS CALL kg_env_create(qs_env, qs_env%kg_env, qs_kind_set, qs_env%input) END IF - dft_section => section_vals_get_subs_vals(qs_env%input, "DFT") - CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%_SECTION_PARAMETERS_", & - l_val=qs_env%energy_correction) dft_section => section_vals_get_subs_vals(qs_env%input, "DFT") CALL section_vals_val_get(dft_section, "EXCITED_STATES%_SECTION_PARAMETERS_", & l_val=qs_env%excited_state) @@ -472,6 +471,46 @@ CONTAINS CALL ec_env_create(qs_env, ec_env, dft_section, ec_section) CALL set_qs_env(qs_env, ec_env=ec_env) + IF (qs_env%energy_correction) THEN + ! Energy correction with Hartree-Fock exchange + ec_hfx_section => section_vals_get_subs_vals(ec_section, "XC%HF") + CALL section_vals_get(ec_hfx_section, explicit=do_ec_hfx) + + IF (ec_env%do_ec_hfx) THEN + + ! Hybrid functionals require same basis + IF (ec_env%basis_inconsistent) THEN + CALL cp_abort(__LOCATION__, & + "Energy correction methods with hybrid functionals:"// & + "correction and ground state need to use the same basis."// & + "Checked by comparing basis set names only.") + END IF + + ! Similar to RPA_HFX we can check if HFX integrals from the qs_env can be reused + IF (ec_env%do_ec_admm .AND. .NOT. dft_control%do_admm) THEN + CALL cp_abort(__LOCATION__, "Need an ADMM input section for ADMM EC to work") + END IF + + ec_env%reuse_hfx = .TRUE. + IF (.NOT. do_hfx) ec_env%reuse_hfx = .FALSE. + CALL compare_hfx_sections(hfx_section, ec_hfx_section, is_identical, same_except_frac) + IF (.NOT. (is_identical .OR. same_except_frac)) ec_env%reuse_hfx = .FALSE. + IF (dft_control%do_admm .AND. .NOT. ec_env%do_ec_admm) ec_env%reuse_hfx = .FALSE. + + IF (.NOT. ec_env%reuse_hfx) THEN + CALL hfx_create(ec_env%x_data, para_env, ec_hfx_section, atomic_kind_set, & + qs_kind_set, particle_set, dft_control, my_cell, do_exx=(.NOT. ec_env%do_ec_admm), & + nelectron_total=nelectron_total) + ELSE + ec_env%x_data => qs_env%x_data + END IF + END IF + + ! Print information of the EC section + CALL ec_write_input(ec_env) + + END IF + IF (dft_control%qs_control%do_almo_scf) THEN CALL almo_scf_env_create(qs_env) END IF diff --git a/src/qs_force.F b/src/qs_force.F index 1593dcb5b5..4cbb4edf8c 100644 --- a/src/qs_force.F +++ b/src/qs_force.F @@ -31,6 +31,7 @@ MODULE qs_force USE efield_utils, ONLY: calculate_ecore_efield USE energy_corrections, ONLY: energy_correction USE excited_states, ONLY: excited_state_energy + USE hfx_exx, ONLY: calculate_exx USE input_constants, ONLY: ri_mp2_laplace,& ri_mp2_method_gpw,& ri_rpa_method_gpw @@ -69,7 +70,6 @@ MODULE qs_force USE qs_subsys_types, ONLY: qs_subsys_set,& qs_subsys_type USE ri_environment_methods, ONLY: build_ri_matrices - USE rpa_hfx, ONLY: calculate_exx USE rt_propagation_forces, ONLY: calc_c_mat_force,& rt_admm_force USE se_core_core, ONLY: se_core_core_interaction @@ -131,8 +131,9 @@ CONTAINS ispin, iw, natom, nkind, nspin, & output_unit INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, natom_of_kind - LOGICAL :: do_admm, do_exx, do_gw, has_unit_metric, & - omit_headers, perform_ec + LOGICAL :: do_admm, do_exx, do_gw, do_im_time, & + has_unit_metric, omit_headers, & + perform_ec, reuse_hfx REAL(dp) :: dummy_real, dummy_real2(2) TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(cp_logger_type), POINTER :: logger @@ -327,9 +328,23 @@ CONTAINS IF (do_exx) THEN do_gw = qs_env%mp2_env%ri_rpa%do_ri_g0w0 do_admm = qs_env%mp2_env%ri_rpa%do_admm + reuse_hfx = qs_env%mp2_env%ri_rpa%reuse_hfx + do_im_time = qs_env%mp2_env%do_im_time output_unit = cp_logger_get_default_io_unit() dummy_real = 0.0_dp - CALL calculate_exx(qs_env, output_unit, do_gw, do_admm, .TRUE., dummy_real, dummy_real2, dummy_real) + + CALL calculate_exx(qs_env=qs_env, & + unit_nr=output_unit, & + hfx_sections=hfx_sections, & + x_data=qs_env%mp2_env%ri_rpa%x_data, & + do_gw=do_gw, & + do_admm=do_admm, & + calc_forces=.TRUE., & + reuse_hfx=reuse_hfx, & + do_im_time=do_im_time, & + E_ex_from_GW=dummy_real, & + E_admm_from_GW=dummy_real2, & + t3=dummy_real) END IF END IF ELSEIF (.NOT. perform_ec) THEN diff --git a/src/qs_linres_kernel.F b/src/qs_linres_kernel.F index b8767787ea..1b47affac6 100644 --- a/src/qs_linres_kernel.F +++ b/src/qs_linres_kernel.F @@ -747,13 +747,17 @@ CONTAINS !> \param rho_ao ... !> \param qs_env ... !> \param hfx_sections ... +!> \param external_x_data ... +!> \param ex ... !> \note !> Simplified version of subroutine hfx_ks_matrix() ! ************************************************************************************************** - SUBROUTINE hfx_matrix(matrix_ks, rho_ao, qs_env, hfx_sections) + SUBROUTINE hfx_matrix(matrix_ks, rho_ao, qs_env, hfx_sections, external_x_data, ex) TYPE(dbcsr_p_type), DIMENSION(:), TARGET :: matrix_ks, rho_ao TYPE(qs_environment_type), POINTER :: qs_env TYPE(section_vals_type), POINTER :: hfx_sections + TYPE(hfx_type), DIMENSION(:, :), OPTIONAL, TARGET :: external_x_data + REAL(KIND=dp), OPTIONAL :: ex CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_matrix' @@ -762,7 +766,7 @@ CONTAINS LOGICAL :: distribute_fock_matrix, & hfx_treat_lsd_in_core, & s_mstruct_changed - REAL(KIND=dp) :: eh1 + REAL(KIND=dp) :: eh1, ehfx TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, rho_ao_kp TYPE(dft_control_type), POINTER :: dft_control TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data @@ -770,7 +774,7 @@ CONTAINS CALL timeset(routineN, handle) - NULLIFY (dft_control, para_env, matrix_ks_kp, rho_ao_kp) + NULLIFY (dft_control, para_env, matrix_ks_kp, rho_ao_kp, x_data) CALL get_qs_env(qs_env=qs_env, & dft_control=dft_control, & @@ -778,6 +782,8 @@ CONTAINS s_mstruct_changed=s_mstruct_changed, & x_data=x_data) + IF (PRESENT(external_x_data)) x_data => external_x_data + CPASSERT(dft_control%nimages == 1) nspins = dft_control%nspins @@ -795,10 +801,10 @@ CONTAINS rho_ao_kp(1:nspins, 1:1) => rho_ao(1:nspins) DO irep = 1, n_rep_hf - eh1 = 0.0_dp + ehfx = 0.0_dp IF (x_data(irep, 1)%do_hfx_ri) THEN - CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_kp, eh1, & + CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_kp, ehfx, & rho_ao=rho_ao_kp, geometry_did_change=s_mstruct_changed, & nspins=nspins, hf_fraction=x_data(irep, 1)%general_parameter%fraction) @@ -807,11 +813,15 @@ CONTAINS DO ispin = 1, mspin CALL integrate_four_center(qs_env, x_data, matrix_ks_kp, eh1, rho_ao_kp, hfx_sections, para_env, & s_mstruct_changed, irep, distribute_fock_matrix, ispin=ispin) + ehfx = ehfx + eh1 END DO END IF END DO + ! Export energy + IF (PRESENT(ex)) ex = ehfx + CALL timestop(handle) END SUBROUTINE hfx_matrix diff --git a/src/rpa_gw_sigma_x.F b/src/rpa_gw_sigma_x.F index f9cca17133..f7353a6bfd 100644 --- a/src/rpa_gw_sigma_x.F +++ b/src/rpa_gw_sigma_x.F @@ -38,6 +38,9 @@ MODULE rpa_gw_sigma_x dbcsr_p_type, dbcsr_release, dbcsr_release_p, dbcsr_set, dbcsr_type, & dbcsr_type_antisymmetric, dbcsr_type_symmetric USE hfx_energy_potential, ONLY: integrate_four_center + USE hfx_exx, ONLY: calc_exx_admm_xc_contributions,& + exx_post_hfx,& + exx_pre_hfx USE hfx_ri, ONLY: hfx_ri_update_ks USE input_constants, ONLY: do_admm_basis_projection,& do_admm_purify_none,& @@ -77,9 +80,6 @@ MODULE rpa_gw_sigma_x USE rpa_gw, ONLY: compute_minus_vxc_kpoints,& trafo_to_mo_and_kpoints USE rpa_gw_kpoints_util, ONLY: get_bandstruc_and_k_dependent_MOs - USE rpa_hfx, ONLY: calc_ri_rpa_admm_xc_contributions,& - ri_rpa_post_hfx,& - ri_rpa_pre_hfx !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads @@ -424,7 +424,7 @@ CONTAINS ehfx = 0.0_dp IF (.NOT. do_ri_Sigma_x) THEN - CALL ri_rpa_pre_hfx(qs_env) + CALL exx_pre_hfx(hfx_sections, qs_env%mp2_env%ri_rpa%x_data, qs_env%mp2_env%ri_rpa%reuse_hfx) calc_ints = .NOT. qs_env%mp2_env%ri_rpa%reuse_hfx ! add here HFX (=Sigma_exchange) to matrix_sigma_x_minus_vxc @@ -460,9 +460,14 @@ CONTAINS !ADMM XC correction IF (do_admm_rpa) THEN - CALL calc_ri_rpa_admm_xc_contributions(qs_env, matrix_ks, matrix_ks_aux_fit, & - energy_xc_admm(1), energy_xc_admm(2), & - .FALSE., .FALSE.) + CALL calc_exx_admm_xc_contributions(qs_env=qs_env, & + matrix_prim=matrix_ks, & + matrix_aux=matrix_ks_aux_fit, & + x_data=qs_env%mp2_env%ri_rpa%x_data, & + exc=energy_xc_admm(1), & + exc_aux_fit=energy_xc_admm(2), & + calc_forces=.FALSE., & + use_virial=.FALSE.) END IF IF (do_kpoints_from_Gamma .AND. print_exx == gw_print_exx) THEN @@ -478,7 +483,7 @@ CONTAINS END IF - CALL ri_rpa_post_hfx(qs_env) + CALL exx_post_hfx(qs_env, qs_env%mp2_env%ri_rpa%x_data, qs_env%mp2_env%ri_rpa%reuse_hfx) END IF energy_ex = ehfx diff --git a/src/rpa_im_time_force_methods.F b/src/rpa_im_time_force_methods.F index 7945c27e93..72bc71292c 100644 --- a/src/rpa_im_time_force_methods.F +++ b/src/rpa_im_time_force_methods.F @@ -59,6 +59,7 @@ MODULE rpa_im_time_force_methods USE gaussian_gridlevels, ONLY: gaussian_gridlevel USE hfx_admm_utils, ONLY: tddft_hfx_matrix USE hfx_derivatives, ONLY: derivatives_four_center + USE hfx_exx, ONLY: add_exx_to_rhs USE hfx_ri, ONLY: get_2c_der_force,& get_force_from_3c_trace,& get_idx_to_atom,& @@ -158,7 +159,6 @@ MODULE rpa_im_time_force_methods USE realspace_grid_types, ONLY: map_gaussian_here,& realspace_grid_type USE response_solver, ONLY: response_equation_new - USE rpa_hfx, ONLY: add_exx_to_rhs USE rpa_im_time, ONLY: compute_mat_dm_global USE rpa_im_time_force_types, ONLY: im_time_force_type USE rs_pw_interface, ONLY: potential_pw2rs @@ -2322,7 +2322,17 @@ CONTAINS hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF") CALL section_vals_get(hfx_section, explicit=do_exx) END IF - IF (do_exx) CALL add_exx_to_rhs(force_data%sum_O_tau, qs_env) + + IF (do_exx) THEN + CALL add_exx_to_rhs(rhs=force_data%sum_O_tau, & + qs_env=qs_env, & + ext_hfx_section=hfx_section, & + x_data=qs_env%mp2_env%ri_rpa%x_data, & + recalc_integrals=.FALSE., & + do_admm=qs_env%mp2_env%ri_rpa%do_admm, & + do_exx=do_exx, & + reuse_hfx=qs_env%mp2_env%ri_rpa%reuse_hfx) + END IF focc = 2.0_dp IF (nspins == 1) focc = 4.0_dp diff --git a/src/rpa_rse.F b/src/rpa_rse.F index ea5ed7229b..d9786b109f 100644 --- a/src/rpa_rse.F +++ b/src/rpa_rse.F @@ -40,6 +40,8 @@ MODULE rpa_rse dbcsr_set,& dbcsr_type_symmetric USE hfx_energy_potential, ONLY: integrate_four_center + USE hfx_exx, ONLY: exx_post_hfx,& + exx_pre_hfx USE hfx_ri, ONLY: hfx_ri_update_ks USE input_section_types, ONLY: section_vals_get,& section_vals_get_subs_vals,& @@ -60,8 +62,6 @@ MODULE rpa_rse USE qs_rho_types, ONLY: qs_rho_get,& qs_rho_type USE qs_vxc, ONLY: qs_vxc_create - USE rpa_hfx, ONLY: ri_rpa_post_hfx,& - ri_rpa_pre_hfx !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads @@ -346,7 +346,7 @@ CONTAINS CALL dbcsr_set(P_mu_nu(is)%matrix, 0.0_dp) END DO - CALL ri_rpa_pre_hfx(qs_env) + CALL exx_pre_hfx(hfx_sections, qs_env%mp2_env%ri_rpa%x_data, qs_env%mp2_env%ri_rpa%reuse_hfx) DO is = 1, ns CALL copy_fm_to_dbcsr(fm_P_mu_nu(is), P_mu_nu(1)%matrix, keep_sparsity=.TRUE.) @@ -392,7 +392,7 @@ CONTAINS fm_X_ao_mo, mo_coeff(is), 1.0_dp, fm_X_mo(is)) END DO - CALL ri_rpa_post_hfx(qs_env) + CALL exx_post_hfx(qs_env, qs_env%mp2_env%ri_rpa%x_data, qs_env%mp2_env%ri_rpa%reuse_hfx) ! Release dbcsr objects DO is = 1, SIZE(P_mu_nu) diff --git a/tests/QS/regtest-dcdft-hfx/N2_t01.inp b/tests/QS/regtest-dcdft-hfx/N2_t01.inp new file mode 100644 index 0000000000..ae4f553157 --- /dev/null +++ b/tests/QS/regtest-dcdft-hfx/N2_t01.inp @@ -0,0 +1,73 @@ +&FORCE_EVAL + METHOD Quickstep + + &PRINT + &FORCES ON + &END FORCES + &END PRINT + + &DFT + BASIS_SET_FILE_NAME HFX_BASIS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 100 + REL_CUTOFF 20 + &END MGRID + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END POISSON + &QS + EPS_DEFAULT 1.E-12 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS HARRIS + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + EPS 1.0E-6 + &END RESPONSE_SOLVER + &XC + &XC_FUNCTIONAL B3LYP + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PBE0 + &END + &END XC + &END DFT + + &SUBSYS + &CELL + ABC 5.0 5.0 5.0 + PERIODIC NONE + &END CELL + &COORD + N 0.000000 0.000000 0.650000 + N 0.000000 0.000000 -0.650000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + BASIS_SET HARRIS DZVP-GTH + POTENTIAL GTH-PBE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW + &REFERENCES OFF + &END REFERENCES +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-hfx/N2_t02.inp b/tests/QS/regtest-dcdft-hfx/N2_t02.inp new file mode 100644 index 0000000000..9a00b12a38 --- /dev/null +++ b/tests/QS/regtest-dcdft-hfx/N2_t02.inp @@ -0,0 +1,73 @@ +&FORCE_EVAL + METHOD Quickstep + + &PRINT + &FORCES ON + &END FORCES + &END PRINT + + &DFT + BASIS_SET_FILE_NAME HFX_BASIS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + REL_CUTOFF 30 + &END MGRID + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END POISSON + &QS + EPS_DEFAULT 1.E-12 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS HARRIS + &RESPONSE_SOLVER + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + EPS 1.0E-6 + &END RESPONSE_SOLVER + &XC + &XC_FUNCTIONAL B3LYP + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PBE0 + &END + &END XC + &END DFT + + &SUBSYS + &CELL + ABC 5.0 5.0 5.0 + PERIODIC NONE + &END CELL + &COORD + N 0.000000 0.000000 0.650000 + N 0.000000 0.000000 -0.650000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + BASIS_SET HARRIS DZVP-GTH + POTENTIAL GTH-PBE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW + &REFERENCES OFF + &END REFERENCES +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-hfx/N2_t03.inp b/tests/QS/regtest-dcdft-hfx/N2_t03.inp new file mode 100644 index 0000000000..83eb4fcc7c --- /dev/null +++ b/tests/QS/regtest-dcdft-hfx/N2_t03.inp @@ -0,0 +1,112 @@ +&FORCE_EVAL + METHOD Quickstep + + &PRINT + &FORCES ON + &END FORCES + &END PRINT + + &DFT + BASIS_SET_FILE_NAME HFX_BASIS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 100 + REL_CUTOFF 30 + &END MGRID + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END POISSON + &QS + EPS_DEFAULT 1.E-12 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS HARRIS + &RESPONSE_SOLVER + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + EPS 1.0E-6 + &END RESPONSE_SOLVER + &XC + ! use a PBE0 functional + ! By using a different EPS_SCHWARZ + ! the HFX section of the ground-state calculation is not reused + &XC_FUNCTIONAL + &PBE + ! 75% GGA exchange + SCALE_X 0.75 + ! 100% GGA correlation + SCALE_C 1.0 + &END PBE + &END XC_FUNCTIONAL + &HF + ! 25 % HFX exchange + FRACTION 0.25 + &SCREENING + EPS_SCHWARZ 1.0E-4 + &END + &INTERACTION_POTENTIAL + POTENTIAL_TYPE TRUNCATED + CUTOFF_RADIUS 2.4 + T_C_G_DATA t_c_g.dat + &END INTERACTION_POTENTIAL + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL + &PBE + ! 75% GGA exchange + SCALE_X 0.75 + ! 100% GGA correlation + SCALE_C 1.0 + &END PBE + &END XC_FUNCTIONAL + &HF + ! 25 % HFX exchange + FRACTION 0.25 + &SCREENING + EPS_SCHWARZ 1.0E-3 + &END + &INTERACTION_POTENTIAL + POTENTIAL_TYPE TRUNCATED + CUTOFF_RADIUS 2.4 + T_C_G_DATA t_c_g.dat + &END INTERACTION_POTENTIAL + &END + &END XC + &END DFT + + &SUBSYS + &CELL + ABC 5.0 5.0 5.0 + PERIODIC NONE + &END CELL + &COORD + N 0.000000 0.000000 0.650000 + N 0.000000 0.000000 -0.650000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + BASIS_SET HARRIS DZVP-GTH + POTENTIAL GTH-PBE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW + &REFERENCES OFF + &END REFERENCES +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-hfx/N2_t04.inp b/tests/QS/regtest-dcdft-hfx/N2_t04.inp new file mode 100644 index 0000000000..4844c6cc48 --- /dev/null +++ b/tests/QS/regtest-dcdft-hfx/N2_t04.inp @@ -0,0 +1,107 @@ +&FORCE_EVAL + METHOD Quickstep + + &PRINT + &FORCES ON + &END FORCES + &END PRINT + + &DFT + BASIS_SET_FILE_NAME HFX_BASIS + POTENTIAL_FILE_NAME GTH_POTENTIALS + BASIS_SET_FILE_NAME BASIS_ADMM + &AUXILIARY_DENSITY_MATRIX_METHOD + ADMM_PURIFICATION_METHOD NONE + EXCH_CORRECTION_FUNC NONE + EXCH_SCALING_MODEL NONE + METHOD BASIS_PROJECTION + &END + &MGRID + CUTOFF 100 + REL_CUTOFF 30 + &END MGRID + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END POISSON + &QS + EPS_DEFAULT 1.E-12 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS HARRIS + &RESPONSE_SOLVER + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + EPS 1.0E-6 + &END RESPONSE_SOLVER + &XC + ! use a PBE0 functional + ! By using an interaction potential and different EPS_SCHWARZ + ! we do not reuse the HFX section of the ground-state calculation + &XC_FUNCTIONAL + &PBE + ! 75% GGA exchange + SCALE_X 0.75 + ! 100% GGA correlation + SCALE_C 1.0 + &END PBE + &END XC_FUNCTIONAL + &HF + ! 25 % HFX exchange + FRACTION 0.25 + &SCREENING + EPS_SCHWARZ 1.0E-5 + &END + &INTERACTION_POTENTIAL + POTENTIAL_TYPE TRUNCATED + CUTOFF_RADIUS 2.4 + T_C_G_DATA t_c_g.dat + &END INTERACTION_POTENTIAL + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL NONE + &END + &HF + &SCREENING + EPS_SCHWARZ 1.0E-4 + &END + &END HF + &END XC + &END DFT + + &SUBSYS + &CELL + ABC 5.0 5.0 5.0 + PERIODIC NONE + &END CELL + &COORD + N 0.000000 0.000000 0.650000 + N 0.000000 0.000000 -0.650000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + BASIS_SET HARRIS DZVP-GTH + BASIS_SET AUX_FIT cFIT3 + POTENTIAL GTH-PBE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW + &REFERENCES OFF + &END REFERENCES +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-hfx/N2_t05.inp b/tests/QS/regtest-dcdft-hfx/N2_t05.inp new file mode 100644 index 0000000000..8efa76eaa8 --- /dev/null +++ b/tests/QS/regtest-dcdft-hfx/N2_t05.inp @@ -0,0 +1,107 @@ +&FORCE_EVAL + METHOD Quickstep + + &PRINT + &FORCES ON + &END FORCES + &END PRINT + + &DFT + BASIS_SET_FILE_NAME HFX_BASIS + POTENTIAL_FILE_NAME GTH_POTENTIALS + BASIS_SET_FILE_NAME BASIS_ADMM + &AUXILIARY_DENSITY_MATRIX_METHOD + ADMM_PURIFICATION_METHOD NONE + EXCH_CORRECTION_FUNC PBEX + EXCH_SCALING_MODEL NONE + METHOD BASIS_PROJECTION + &END + &MGRID + CUTOFF 100 + REL_CUTOFF 30 + &END MGRID + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END POISSON + &QS + EPS_DEFAULT 1.E-12 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS HARRIS + &RESPONSE_SOLVER + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + EPS 1.0E-6 + &END RESPONSE_SOLVER + &XC + ! use a PBE0 functional + ! By using an interaction potential and different EPS_SCHWARZ + ! we do not reuse the HFX section of the ground-state calculation + &XC_FUNCTIONAL + &PBE + ! 75% GGA exchange + SCALE_X 0.75 + ! 100% GGA correlation + SCALE_C 1.0 + &END PBE + &END XC_FUNCTIONAL + &HF + ! 25 % HFX exchange + FRACTION 0.25 + &SCREENING + EPS_SCHWARZ 1.0E-5 + &END + &INTERACTION_POTENTIAL + POTENTIAL_TYPE TRUNCATED + CUTOFF_RADIUS 2.4 + T_C_G_DATA t_c_g.dat + &END INTERACTION_POTENTIAL + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL NONE + &END + &HF + &SCREENING + EPS_SCHWARZ 1.0E-4 + &END + &END HF + &END XC + &END DFT + + &SUBSYS + &CELL + ABC 5.0 5.0 5.0 + PERIODIC NONE + &END CELL + &COORD + N 0.000000 0.000000 0.650000 + N 0.000000 0.000000 -0.650000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + BASIS_SET HARRIS DZVP-GTH + BASIS_SET AUX_FIT cFIT3 + POTENTIAL GTH-PBE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW + &REFERENCES OFF + &END REFERENCES +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-hfx/N2_t06.inp b/tests/QS/regtest-dcdft-hfx/N2_t06.inp new file mode 100644 index 0000000000..404e234be7 --- /dev/null +++ b/tests/QS/regtest-dcdft-hfx/N2_t06.inp @@ -0,0 +1,108 @@ +&FORCE_EVAL + METHOD Quickstep + + &PRINT + &FORCES ON + &END FORCES + &END PRINT + + &DFT + BASIS_SET_FILE_NAME HFX_BASIS + POTENTIAL_FILE_NAME GTH_POTENTIALS + BASIS_SET_FILE_NAME BASIS_ADMM + &AUXILIARY_DENSITY_MATRIX_METHOD + ADMM_PURIFICATION_METHOD NONE + EXCH_CORRECTION_FUNC PBEX + EXCH_SCALING_MODEL NONE + METHOD BASIS_PROJECTION + &END + &MGRID + CUTOFF 100 + REL_CUTOFF 30 + &END MGRID + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END POISSON + &QS + EPS_DEFAULT 1.E-12 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS HARRIS + ADMM + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + EPS 1.0E-6 + &END RESPONSE_SOLVER + &XC + ! use a PBE0 functional + ! By using an interaction potential and different EPS_SCHWARZ + ! we do not reuse the HFX section of the ground-state calculation + &XC_FUNCTIONAL + &PBE + ! 75% GGA exchange + SCALE_X 0.75 + ! 100% GGA correlation + SCALE_C 1.0 + &END PBE + &END XC_FUNCTIONAL + &HF + ! 25 % HFX exchange + FRACTION 0.25 + &SCREENING + EPS_SCHWARZ 1.0E-5 + &END + &INTERACTION_POTENTIAL + POTENTIAL_TYPE TRUNCATED + CUTOFF_RADIUS 2.4 + T_C_G_DATA t_c_g.dat + &END INTERACTION_POTENTIAL + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL NONE + &END + &HF + &SCREENING + EPS_SCHWARZ 1.0E-4 + &END + &END HF + &END XC + &END DFT + + &SUBSYS + &CELL + ABC 5.0 5.0 5.0 + PERIODIC NONE + &END CELL + &COORD + N 0.000000 0.000000 0.650000 + N 0.000000 0.000000 -0.650000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + BASIS_SET HARRIS DZVP-GTH + BASIS_SET AUX_FIT cFIT3 + POTENTIAL GTH-PBE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW + &REFERENCES OFF + &END REFERENCES +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-hfx/N2_t07.inp b/tests/QS/regtest-dcdft-hfx/N2_t07.inp new file mode 100644 index 0000000000..e4f07982ae --- /dev/null +++ b/tests/QS/regtest-dcdft-hfx/N2_t07.inp @@ -0,0 +1,103 @@ +&FORCE_EVAL + METHOD Quickstep + + &PRINT + &FORCES ON + &END FORCES + &END PRINT + + &DFT + BASIS_SET_FILE_NAME HFX_BASIS + POTENTIAL_FILE_NAME GTH_POTENTIALS + BASIS_SET_FILE_NAME BASIS_ADMM + &AUXILIARY_DENSITY_MATRIX_METHOD + ADMM_PURIFICATION_METHOD NONE + EXCH_CORRECTION_FUNC PBEX + EXCH_SCALING_MODEL NONE + METHOD BASIS_PROJECTION + &END + &MGRID + CUTOFF 100 + REL_CUTOFF 30 + &END MGRID + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END POISSON + &QS + EPS_DEFAULT 1.E-12 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS HARRIS + ADMM + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + EPS 1.0E-6 + &END RESPONSE_SOLVER + &XC + ! use a PBE0 functional + ! By using an interaction potential and different EPS_SCHWARZ + ! we do not reuse the HFX section of the ground-state calculation + &XC_FUNCTIONAL + &PBE + ! 75% GGA exchange + SCALE_X 0.75 + ! 100% GGA correlation + SCALE_C 1.0 + &END PBE + &END XC_FUNCTIONAL + &HF + ! 25 % HFX exchange + FRACTION 0.25 + &SCREENING + EPS_SCHWARZ 1.0E-4 + &END + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL NONE + &END + &HF + &SCREENING + EPS_SCHWARZ 1.0E-4 + &END + &END HF + &END XC + &END DFT + + &SUBSYS + &CELL + ABC 5.0 5.0 5.0 + PERIODIC NONE + &END CELL + &COORD + N 0.000000 0.000000 0.650000 + N 0.000000 0.000000 -0.650000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + BASIS_SET HARRIS DZVP-GTH + BASIS_SET AUX_FIT cFIT3 + POTENTIAL GTH-PBE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW + &REFERENCES OFF + &END REFERENCES +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-hfx/N2_t08.inp b/tests/QS/regtest-dcdft-hfx/N2_t08.inp new file mode 100644 index 0000000000..f94da12c28 --- /dev/null +++ b/tests/QS/regtest-dcdft-hfx/N2_t08.inp @@ -0,0 +1,103 @@ +&FORCE_EVAL + METHOD Quickstep + + &PRINT + &FORCES ON + &END FORCES + &END PRINT + + &DFT + BASIS_SET_FILE_NAME HFX_BASIS + POTENTIAL_FILE_NAME GTH_POTENTIALS + BASIS_SET_FILE_NAME BASIS_ADMM + &AUXILIARY_DENSITY_MATRIX_METHOD + ADMM_PURIFICATION_METHOD NONE + EXCH_CORRECTION_FUNC PBEX + EXCH_SCALING_MODEL NONE + METHOD BASIS_PROJECTION + &END + &MGRID + CUTOFF 100 + REL_CUTOFF 30 + &END MGRID + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END POISSON + &QS + EPS_DEFAULT 1.E-12 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS HARRIS + ADMM + &RESPONSE_SOLVER + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + EPS 1.0E-6 + &END RESPONSE_SOLVER + &XC + ! use a PBE0 functional + ! By using an interaction potential and different EPS_SCHWARZ + ! we do not reuse the HFX section of the ground-state calculation + &XC_FUNCTIONAL + &PBE + ! 75% GGA exchange + SCALE_X 0.75 + ! 100% GGA correlation + SCALE_C 1.0 + &END PBE + &END XC_FUNCTIONAL + &HF + ! 25 % HFX exchange + FRACTION 0.25 + &SCREENING + EPS_SCHWARZ 1.0E-4 + &END + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL NONE + &END + &HF + &SCREENING + EPS_SCHWARZ 1.0E-4 + &END + &END HF + &END XC + &END DFT + + &SUBSYS + &CELL + ABC 5.0 5.0 5.0 + PERIODIC NONE + &END CELL + &COORD + N 0.000000 0.000000 0.650000 + N 0.000000 0.000000 -0.650000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + BASIS_SET HARRIS DZVP-GTH + BASIS_SET AUX_FIT cFIT3 + POTENTIAL GTH-PBE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW + &REFERENCES OFF + &END REFERENCES +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-hfx/TEST_FILES b/tests/QS/regtest-dcdft-hfx/TEST_FILES new file mode 100644 index 0000000000..42e3930c8b --- /dev/null +++ b/tests/QS/regtest-dcdft-hfx/TEST_FILES @@ -0,0 +1,21 @@ +# runs are executed in the same order as in this file +# the second field tells which test should be run in order to compare with the last available output +# see regtest/TEST_FILES +# Density-corrected DFT with exact exchange contribution +# PBE0 - B3LYP - MO solver +N2_t01.inp 11 2e-11 -19.832131213062343 +# PBE0 - B3LYP - AO solver +N2_t02.inp 11 2e-11 -19.833869176317489 +# PBE0 - PBE0 (no hfx_reuse) +N2_t03.inp 11 2e-11 -19.844567389563924 +# HFX(ADMM-NONE) - PBE0 (no hfx_reuse) +N2_t04.inp 11 2e-11 -19.828808345190410 +# HFX(ADMM-PBEX) - PBE0 (no hfx_reuse) +N2_t05.inp 11 2e-11 -19.834993600346927 +# HFX(ADMM-PBEX) - PBE0(ADMM) - MO solver (no hfx_reuse) +N2_t06.inp 11 2e-11 -19.826274785794372 +# HFX(ADMM-PBEX) - PBE0(ADMM) - MO solver (hfx_reuse) +N2_t07.inp 11 2e-11 -19.830296445771687 +# HFX(ADMM-PBEX) - B3LYP(ADMM) - AO solver (hfx_reuse) +N2_t08.inp 11 2e-11 -19.830322314126654 +#EOF diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index ad34d63c0b..04368b3a9a 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -5,6 +5,7 @@ # the order will be regularly checked and modified... QS/regtest-wfn-restart QS/regtest-dcdft-force libint libxc +QS/regtest-dcdft-hfx libint libxc QS/regtest-dcdft-stress libint libxc QS/regtest-ec-stress libint libxc QS/regtest-ri-rpa-grad libint