comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-29 06:35:28 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

EC| Enable DC-DFT with HFX-ADMM for reference and DC calculation (#2780)

Browse source
This commit is contained in:
fbelle 2023-05-22 11:24:52 +02:00 • committed by GitHub
parent d6d429ae51
commit d4cb79545b
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
24 changed files with 1704 additions and 428 deletions
Download patch file Download diff file Expand all files Collapse all files

2
src/CMakeLists.txt
View file

@ -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

33
src/ec_env_types.F
View file

@ -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

291
src/ec_environment.F
View file

@ -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)

461
src/energy_corrections.F
View file

@ -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 <a|V|b>
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 <nabla.a|V|nabla.b>
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 <nabla.a|V|nabla.b>
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 <a|V|b>
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 <nabla.a|V|nabla.b>
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

353
src/rpa_hfx.F → src/hfx_exx.F
View file

@ -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

13
src/input_cp2k_ec.F
View file

@ -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
! **************************************************************************************************

21
src/mp2.F
View file

@ -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

11
src/mp2_cphf.F
View file

@ -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

53
src/qs_environment.F
View file

@ -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

23
src/qs_force.F
View file

@ -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

20
src/qs_linres_kernel.F
View file

@ -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

21
src/rpa_gw_sigma_x.F
View file

@ -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

14
src/rpa_im_time_force_methods.F
View file

@ -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

8
src/rpa_rse.F
View file

@ -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)

73
tests/QS/regtest-dcdft-hfx/N2_t01.inp Normal file
View file

@ -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

73
tests/QS/regtest-dcdft-hfx/N2_t02.inp Normal file
View file

@ -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

112
tests/QS/regtest-dcdft-hfx/N2_t03.inp Normal file
View file

@ -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

107
tests/QS/regtest-dcdft-hfx/N2_t04.inp Normal file
View file

@ -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

107
tests/QS/regtest-dcdft-hfx/N2_t05.inp Normal file
View file

@ -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

108
tests/QS/regtest-dcdft-hfx/N2_t06.inp Normal file
View file

@ -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

103
tests/QS/regtest-dcdft-hfx/N2_t07.inp Normal file
View file

@ -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

103
tests/QS/regtest-dcdft-hfx/N2_t08.inp Normal file
View file

@ -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

21
tests/QS/regtest-dcdft-hfx/TEST_FILES Normal file
View file

@ -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

1
tests/TEST_DIRS
View file

@ -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