Skip to content

Commit

Permalink
RI implementation of Hartree-Fock exchange
Browse files Browse the repository at this point in the history
  • Loading branch information
@pseewald
pseewald committed Mar 7, 2020
1 parent 295a4e9 commit f874287
Show file tree
Hide file tree
Showing 29 changed files with 2,146 additions and 86 deletions.
5 changes: 4 additions & 1 deletion src/aobasis/basis_set_container_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,8 @@ MODULE basis_set_container_types
ri_hxc_basis = 111, &
ri_k_basis = 112, &
ri_xas_basis = 113, &
aux_fit_soft_basis = 114
aux_fit_soft_basis = 114, &
ri_hfx_basis = 115
! **************************************************************************************************
TYPE basis_set_container_type
PRIVATE
Expand Down Expand Up @@ -99,6 +100,8 @@ FUNCTION get_basis_type(basis_set_type) RESULT(basis_type_nr)
basis_type_nr = ri_aux_basis
CASE ("RI_HXC")
basis_type_nr = ri_hxc_basis
CASE ("RI_HFX")
basis_type_nr = ri_hfx_basis
CASE ("RI_K")
basis_type_nr = ri_k_basis
CASE ("LRI_AUX")
Expand Down
10 changes: 8 additions & 2 deletions src/auto_basis.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,13 +47,15 @@ MODULE auto_basis
!> \param qs_kind ...
!> \param basis_cntrl ...
!> \param basis_sort ...
!> \param basis_type ...
!> \date 01.11.2017
!> \author JGH
! **************************************************************************************************
SUBROUTINE create_ri_aux_basis_set(ri_aux_basis_set, qs_kind, basis_cntrl, basis_sort)
SUBROUTINE create_ri_aux_basis_set(ri_aux_basis_set, qs_kind, basis_cntrl, basis_type, basis_sort)
TYPE(gto_basis_set_type), POINTER :: ri_aux_basis_set
TYPE(qs_kind_type), INTENT(IN) :: qs_kind
INTEGER, INTENT(IN) :: basis_cntrl
CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: basis_type
INTEGER, INTENT(IN), OPTIONAL :: basis_sort

CHARACTER(len=*), PARAMETER :: routineN = 'create_ri_aux_basis_set', &
Expand Down Expand Up @@ -84,7 +86,11 @@ SUBROUTINE create_ri_aux_basis_set(ri_aux_basis_set, qs_kind, basis_cntrl, basis
!
CPASSERT(.NOT. ASSOCIATED(ri_aux_basis_set))
NULLIFY (orb_basis_set)
CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, basis_type="ORB")
IF (.NOT. PRESENT(basis_type)) THEN
CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, basis_type="ORB")
ELSE
CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, basis_type=basis_type)
ENDIF
IF (ASSOCIATED(orb_basis_set)) THEN
CALL get_basis_keyfigures(orb_basis_set, lmax, zmin, zmax, zeff)
!Note: RI basis coud require lmax up to 2*orb_lmax. This ensures that all orbital pointers
Expand Down
9 changes: 7 additions & 2 deletions src/basis_set_output.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,8 +61,8 @@ SUBROUTINE print_basis_set_file(qs_env, base_section)
INTEGER, SAVE :: ncalls = 0
TYPE(cp_logger_type), POINTER :: logger
TYPE(gto_basis_set_type), POINTER :: aux_fit_basis, lri_aux_basis, orb_basis, &
ri_aux_basis, ri_hxc_basis, &
ri_xas_basis
ri_aux_basis, ri_hfx_basis, &
ri_hxc_basis, ri_xas_basis
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_kind_type), POINTER :: qs_kind

Expand Down Expand Up @@ -93,6 +93,7 @@ SUBROUTINE print_basis_set_file(qs_env, base_section)
CALL get_qs_kind(qs_kind, basis_set=orb_basis, basis_type="ORB")
CALL get_qs_kind(qs_kind, basis_set=ri_aux_basis, basis_type="RI_AUX")
CALL get_qs_kind(qs_kind, basis_set=ri_hxc_basis, basis_type="RI_HXC")
CALL get_qs_kind(qs_kind, basis_set=ri_hfx_basis, basis_type="RI_HFX")
CALL get_qs_kind(qs_kind, basis_set=lri_aux_basis, basis_type="LRI_AUX")
CALL get_qs_kind(qs_kind, basis_set=aux_fit_basis, basis_type="AUX_FIT")
CALL get_qs_kind(qs_kind, basis_set=ri_xas_basis, basis_type="RI_XAS")
Expand Down Expand Up @@ -121,6 +122,10 @@ SUBROUTINE print_basis_set_file(qs_env, base_section)
bname = "local_ri_xas"
CALL basis_out(ri_xas_basis, element_symbol, bname, iunit)
END IF
IF (ASSOCIATED(ri_hfx_basis)) THEN
bname = "local_ri_hfx"
CALL basis_out(ri_hfx_basis, element_symbol, bname, iunit)
ENDIF
ENDIF
END DO

Expand Down
3 changes: 2 additions & 1 deletion src/cp_control_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -495,7 +495,8 @@ MODULE cp_control_types
auto_basis_aux_fit = 1, &
auto_basis_lri_aux = 1, &
auto_basis_ri_hxc = 1, &
auto_basis_ri_xas = 1
auto_basis_ri_xas = 1, &
auto_basis_ri_hfx
REAL(KIND=dp) :: relax_multiplicity, &
sic_scaling_a, &
sic_scaling_b, &
Expand Down
2 changes: 2 additions & 0 deletions src/cp_control_utils.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -207,6 +207,8 @@ SUBROUTINE read_dft_control(dft_control, dft_section)
dft_control%auto_basis_ri_hxc = isize
CASE ("RI_XAS")
dft_control%auto_basis_ri_xas = isize
CASE ("RI_HFX")
dft_control%auto_basis_ri_hfx = isize
CASE DEFAULT
CPWARN("Unknown basis type in AUTO_BASIS keyword:"//TRIM(tmpstringlist(1)))
END SELECT
Expand Down
77 changes: 62 additions & 15 deletions src/hfx_admm_utils.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,12 +35,14 @@ MODULE hfx_admm_utils
cp_print_key_unit_nr
USE cp_para_types, ONLY: cp_para_env_type
USE dbcsr_api, ONLY: dbcsr_add,&
dbcsr_copy,&
dbcsr_p_type,&
dbcsr_set,&
dbcsr_type
USE external_potential_types, ONLY: copy_potential
USE hfx_derivatives, ONLY: derivatives_four_center
USE hfx_energy_potential, ONLY: integrate_four_center
USE hfx_ri, ONLY: hfx_ri_update_ks
USE hfx_types, ONLY: hfx_type
USE input_constants, ONLY: &
do_admm_aux_exch_func_bee, do_admm_aux_exch_func_default, do_admm_aux_exch_func_none, &
Expand Down Expand Up @@ -345,6 +347,7 @@ SUBROUTINE hfx_ks_matrix(qs_env, matrix_ks, rho, energy, calculate_forces, &
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_ks_orb, rho_ao_orb
TYPE(dft_control_type), POINTER :: dft_control
TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mo_array
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_poisson_type), POINTER :: poisson_env
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
Expand All @@ -358,9 +361,9 @@ SUBROUTINE hfx_ks_matrix(qs_env, matrix_ks, rho, energy, calculate_forces, &
CALL timeset(routineN, handle)

NULLIFY (auxbas_pw_pool, dft_control, force, hfx_sections, input, &
para_env, poisson_env, pw_env, virial, &
matrix_ks_im, matrix_ks_orb, rho_ao_orb, &
matrix_h, matrix_ks_aux_fit, matrix_ks_aux_fit_im, matrix_ks_aux_fit_hfx)
para_env, poisson_env, pw_env, virial, matrix_ks_im, &
matrix_ks_orb, rho_ao_orb, matrix_h, matrix_ks_aux_fit, &
matrix_ks_aux_fit_im, matrix_ks_aux_fit_hfx)

CALL get_qs_env(qs_env=qs_env, &
dft_control=dft_control, &
Expand All @@ -376,6 +379,12 @@ SUBROUTINE hfx_ks_matrix(qs_env, matrix_ks, rho, energy, calculate_forces, &
s_mstruct_changed=s_mstruct_changed, &
x_data=x_data)

IF (dft_control%do_admm) THEN
CALL get_qs_env(qs_env=qs_env, mos_aux_fit=mo_array)
ELSE
CALL get_qs_env(qs_env=qs_env, mos=mo_array)
ENDIF

nspins = dft_control%nspins
nimages = dft_control%nimages

Expand Down Expand Up @@ -433,12 +442,33 @@ SUBROUTINE hfx_ks_matrix(qs_env, matrix_ks, rho, energy, calculate_forces, &
CALL qs_rho_get(rho_struct=rho_orb, rho_ao_kp=rho_ao_orb)
! Finally the real hfx calulation
ehfx = 0.0_dp
DO ispin = 1, mspin
CALL integrate_four_center(qs_env, x_data, matrix_ks_orb, eh1, rho_ao_orb, hfx_sections, &
para_env, s_mstruct_changed, irep, distribute_fock_matrix, &
ispin=ispin)
ehfx = ehfx + eh1
END DO

IF (x_data(irep, 1)%do_hfx_ri .AND. calculate_forces) THEN
CPABORT("RI forces not yet implemented in HFX")
ENDIF

IF (x_data(irep, 1)%do_hfx_ri) THEN
CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_orb, ehfx, &
mo_array, rho_ao_orb, &
s_mstruct_changed, nspins, &
x_data(irep, 1)%general_parameter%fraction)
IF (dft_control%do_admm) THEN
!for ADMMS, we need the exchange matrix k(d) for both spins

DO ispin = 1, mspin
CALL dbcsr_copy(matrix_ks_aux_fit_hfx(ispin)%matrix, matrix_ks_orb(ispin, 1)%matrix, &
name="HF exch. part of matrix_ks_aux_fit for ADMMS")
ENDDO
END IF
ELSE

DO ispin = 1, mspin
CALL integrate_four_center(qs_env, x_data, matrix_ks_orb, eh1, rho_ao_orb, hfx_sections, &
para_env, s_mstruct_changed, irep, distribute_fock_matrix, &
ispin=ispin)
ehfx = ehfx + eh1
END DO
ENDIF

IF (calculate_forces .AND. .NOT. do_adiabatic_rescaling) THEN
!Scale auxiliary density matrix for ADMMP (see Merlot2014) with gsi(ispin) to scale force
Expand Down Expand Up @@ -487,12 +517,29 @@ SUBROUTINE hfx_ks_matrix(qs_env, matrix_ks, rho, energy, calculate_forces, &
rho_ao_orb(1:ns, 1:1) => rho_ao_1d(1:ns)

ehfxrt = 0.0_dp
DO ispin = 1, mspin
CALL integrate_four_center(qs_env, x_data, matrix_ks_orb, eh1, rho_ao_orb, hfx_sections, &
para_env, .FALSE., irep, distribute_fock_matrix, &
ispin=ispin)
ehfxrt = ehfxrt + eh1
END DO

IF (x_data(irep, 1)%do_hfx_ri) THEN
CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_orb, ehfx, &
mo_array, rho_ao_orb, &
.FALSE., nspins, &
x_data(irep, 1)%general_parameter%fraction)
IF (dft_control%do_admm) THEN
!for ADMMS, we need the exchange matrix k(d) for both spins

DO ispin = 1, mspin
CALL dbcsr_copy(matrix_ks_aux_fit_hfx(ispin)%matrix, matrix_ks_orb(ispin, 1)%matrix, &
name="HF exch. part of matrix_ks_aux_fit for ADMMS")
ENDDO
END IF

ELSE
DO ispin = 1, mspin
CALL integrate_four_center(qs_env, x_data, matrix_ks_orb, eh1, rho_ao_orb, hfx_sections, &
para_env, .FALSE., irep, distribute_fock_matrix, &
ispin=ispin)
ehfxrt = ehfxrt + eh1
END DO
ENDIF

IF (calculate_forces .AND. .NOT. do_adiabatic_rescaling) THEN
NULLIFY (rho_ao_resp)
Expand Down

0 comments on commit f874287

Please sign in to comment.