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

RI-HFX forces

Browse source
This commit is contained in:
abussy 2021-09-09 18:34:19 +02:00 • committed by Augustin Bussy
parent c346f4425c
commit 51cc574b6b
37 changed files with 5058 additions and 253 deletions
Download patch file Download diff file Expand all files Collapse all files

66
src/hfx_admm_utils.F
View file

@ -38,13 +38,15 @@ MODULE hfx_admm_utils
USE hfx_derivatives, ONLY: derivatives_four_center
USE hfx_energy_potential, ONLY: integrate_four_center
USE hfx_pw_methods, ONLY: pw_hfx
USE hfx_ri, ONLY: hfx_ri_update_ks
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_bee, do_admm_aux_exch_func_default, do_admm_aux_exch_func_none, &
do_admm_aux_exch_func_opt, do_admm_aux_exch_func_pbex, do_potential_coulomb, &
do_potential_long, do_potential_mix_cl, do_potential_mix_cl_trunc, do_potential_short, &
do_potential_truncated, xc_funct_no_shortcut, xc_none
USE input_cp2k_hfx, ONLY: ri_mo
USE input_section_types, ONLY: section_vals_duplicate,&
section_vals_get,&
section_vals_get_subs_vals,&
@ -421,10 +423,6 @@ CONTAINS
! Finally the real hfx calulation
ehfx = 0.0_dp
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, &
@ -432,11 +430,12 @@ CONTAINS
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
DO ispin = 1, nspins
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
@ -461,8 +460,22 @@ CONTAINS
ELSE
NULLIFY (rho_ao_resp)
END IF
CALL derivatives_four_center(qs_env, rho_ao_orb, rho_ao_resp, hfx_sections, &
para_env, irep, 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_orb, mos=mo_array, &
rho_ao_resp=rho_ao_resp, &
use_virial=use_virial)
ELSE
CALL derivatives_four_center(qs_env, rho_ao_orb, rho_ao_resp, hfx_sections, &
para_env, irep, use_virial)
END IF
!Scale auxiliary density matrix for ADMMP back with 1/gsi(ispin)
IF (dft_control%do_admm) THEN
CALL scale_dm(qs_env, rho_ao_orb, scale_back=.TRUE.)
@ -506,7 +519,7 @@ CONTAINS
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
DO ispin = 1, nspins
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
@ -523,8 +536,18 @@ CONTAINS
IF (calculate_forces .AND. .NOT. do_adiabatic_rescaling) THEN
NULLIFY (rho_ao_resp)
CALL derivatives_four_center(qs_env, rho_ao_orb, rho_ao_resp, hfx_sections, &
para_env, irep, 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_orb, mos=mo_array, &
use_virial=use_virial)
ELSE
CALL derivatives_four_center(qs_env, rho_ao_orb, rho_ao_resp, hfx_sections, &
para_env, irep, use_virial)
END IF
END IF
!! If required, the calculation of the forces will be done later with adiabatic rescaling
@ -1211,11 +1234,22 @@ CONTAINS
DO irep = 1, n_rep_hf
! the real hfx calulation
ehfx = 0.0_dp
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
IF (x_data(irep, 1)%do_hfx_ri) THEN
IF (x_data(irep, 1)%ri_data%flavor == ri_mo) THEN
CPABORT("NYI with RI_FLAVOR MO")
END IF
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)
ELSE
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
IF (my_update_energy) energy%ex = ehfx

1919
src/hfx_ri.F
View file

File diff suppressed because it is too large Load diff

72
src/hfx_types.F
View file

@ -36,14 +36,14 @@ MODULE hfx_types
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_type
USE cp_output_handling, ONLY: cp_print_key_finished_output,&
cp_print_key_unit_nr
cp_print_key_unit_nr,&
debug_print_level
USE cp_para_types, ONLY: cp_para_env_type
USE cp_units, ONLY: cp_unit_from_cp2k
USE dbcsr_tensor_api, ONLY: &
dbcsr_t_batched_contract_finalize, dbcsr_t_default_distvec, dbcsr_t_destroy, &
dbcsr_t_distribution_destroy, dbcsr_t_distribution_new, dbcsr_t_distribution_type, &
dbcsr_t_mp_dims_create, dbcsr_t_pgrid_create, dbcsr_t_pgrid_destroy, dbcsr_t_pgrid_type, &
dbcsr_t_type
dbcsr_t_create, dbcsr_t_default_distvec, dbcsr_t_destroy, dbcsr_t_distribution_destroy, &
dbcsr_t_distribution_new, dbcsr_t_distribution_type, dbcsr_t_mp_dims_create, &
dbcsr_t_pgrid_create, dbcsr_t_pgrid_destroy, dbcsr_t_pgrid_type, dbcsr_t_type
USE hfx_helpers, ONLY: count_cells_perd,&
next_image_cell_perd
USE input_constants, ONLY: &
@ -55,7 +55,8 @@ MODULE hfx_types
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
section_vals_val_get,&
section_vals_val_set
USE kinds, ONLY: default_path_length,&
default_string_length,&
dp,&
@ -375,7 +376,7 @@ MODULE hfx_types
REAL(KIND=dp) :: filter_eps, filter_eps_2c, filter_eps_storage, filter_eps_mo, &
eps_lanczos, eps_pgf_orb, eps_eigval
INTEGER :: t2c_sqrt_order, max_iter_lanczos, flavor, unit_nr_dbcsr, unit_nr, &
min_bsize, min_bsize_MO, t2c_method, nelectron_total
min_bsize, max_bsize_MO, t2c_method, nelectron_total
LOGICAL :: check_2c_inv, calc_condnum
TYPE(libint_potential_type) :: ri_metric
@ -383,6 +384,7 @@ MODULE hfx_types
! input parameters from hfx
TYPE(libint_potential_type) :: hfx_pot ! interaction potential
REAL(KIND=dp) :: eps_schwarz ! integral screening threshold
REAL(KIND=dp) :: eps_schwarz_forces ! integral derivatives screening threshold
LOGICAL :: same_op ! whether RI operator is same as HF potential
@ -411,12 +413,13 @@ MODULE hfx_types
! Note: changed static DIMENSION(1,1) of dbcsr_t_type to allocatables as workaround for gfortran 8.3.0,
! with static dimension gfortran gets stuck during compilation
! 2c tensors in (RI | RI) format for forces
TYPE(dbcsr_t_type), DIMENSION(:, :), ALLOCATABLE :: t_2c_inv
TYPE(dbcsr_t_type), DIMENSION(:, :), ALLOCATABLE :: t_2c_pot
! 2c tensor in (RI | RI) format for contraction
TYPE(dbcsr_t_type), DIMENSION(:, :), ALLOCATABLE :: t_2c_int
! 3c integral tensor in default format
TYPE(dbcsr_t_type), DIMENSION(:, :), ALLOCATABLE :: t_3c_int
! 3c integral tensor in (AO RI | AO) format for contraction
TYPE(dbcsr_t_type), DIMENSION(:, :), ALLOCATABLE :: t_3c_int_ctr_1
TYPE(block_ind_type), DIMENSION(:, :), ALLOCATABLE :: blk_indices
@ -974,7 +977,7 @@ CONTAINS
CALL hfx_ri_init_read_input_from_hfx(actual_x_data%ri_data, actual_x_data, hfx_section, &
hf_sub_section, qs_kind_set, &
particle_set, atomic_kind_set, dft_control, para_env, irep, do_ot, &
nelectron_total)
nelectron_total, my_do_exx)
ENDIF
END DO
@ -1001,9 +1004,11 @@ CONTAINS
!> \param irep ...
!> \param do_ot ...
!> \param nelectron_total ...
!> \param do_exx ...
! **************************************************************************************************
SUBROUTINE hfx_ri_init_read_input_from_hfx(ri_data, x_data, hfx_section, ri_section, qs_kind_set, &
particle_set, atomic_kind_set, dft_control, para_env, irep, do_ot, nelectron_total)
particle_set, atomic_kind_set, dft_control, para_env, irep, &
do_ot, nelectron_total, do_exx)
TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
TYPE(hfx_type), INTENT(INOUT) :: x_data
TYPE(section_vals_type), POINTER :: hfx_section, ri_section
@ -1015,6 +1020,7 @@ CONTAINS
INTEGER, INTENT(IN) :: irep
LOGICAL, INTENT(IN) :: do_ot
INTEGER, INTENT(IN) :: nelectron_total
LOGICAL, INTENT(IN) :: do_exx
CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_ri_init_read_input_from_hfx'
@ -1037,6 +1043,7 @@ CONTAINS
ri_data%ri_section => ri_section
ri_data%hfx_section => hfx_section
ri_data%eps_schwarz = x_data%screening_parameter%eps_schwarz
ri_data%eps_schwarz_forces = x_data%screening_parameter%eps_schwarz_forces
logger => cp_get_default_logger()
unit_nr_dbcsr = cp_print_key_unit_nr(logger, ri_data%ri_section, "RI_INFO", &
@ -1057,7 +1064,7 @@ CONTAINS
t_c_filename = char_val
END IF
IF (dft_control%do_admm) THEN
IF (dft_control%do_admm .AND. (.NOT. do_exx)) THEN
orb_basis_type = "AUX_FIT"
ELSE
orb_basis_type = "ORB"
@ -1112,6 +1119,7 @@ CONTAINS
INTEGER :: handle
LOGICAL :: explicit
REAL(dp) :: eps_storage_scaling
TYPE(section_vals_type), POINTER :: prog_run_info
CALL timeset(routineN, handle)
@ -1148,7 +1156,7 @@ CONTAINS
CALL section_vals_val_get(ri_section, "RI_FLAVOR", i_val=ri_data%flavor)
CALL section_vals_val_get(ri_section, "EPS_PGF_ORB", r_val=ri_data%eps_pgf_orb)
CALL section_vals_val_get(ri_section, "MIN_BLOCK_SIZE", i_val=ri_data%min_bsize)
CALL section_vals_val_get(ri_section, "MIN_BLOCK_SIZE_MO", i_val=ri_data%min_bsize_MO)
CALL section_vals_val_get(ri_section, "MAX_BLOCK_SIZE_MO", i_val=ri_data%max_bsize_MO)
CALL section_vals_val_get(ri_section, "MEMORY_CUT", i_val=ri_data%n_mem)
IF (ri_data%flavor == ri_pmat) THEN
@ -1169,6 +1177,9 @@ CONTAINS
ri_data%loc_subsection => section_vals_get_subs_vals(ri_section, "LOCALIZE")
ri_data%print_loc_subsection => section_vals_get_subs_vals(ri_data%loc_subsection, "PRINT")
prog_run_info => section_vals_get_subs_vals(ri_data%print_loc_subsection, "PROGRAM_RUN_INFO")
CALL section_vals_val_set(prog_run_info, keyword_name="_SECTION_PARAMETERS_", &
i_val=debug_print_level) !keep output clean
CALL section_vals_val_get(ri_data%loc_subsection, "_SECTION_PARAMETERS_", l_val=ri_data%do_loc)
@ -1347,12 +1358,12 @@ CONTAINS
DEALLOCATE (dist1, dist2)
ELSEIF (ri_data%flavor == ri_mo) THEN
ALLOCATE (ri_data%t_3c_int(1, 1))
ALLOCATE (ri_data%t_2c_int(2, 1))
CALL create_2c_tensor(ri_data%t_2c_int(1, 1), dist1, dist2, ri_data%pgrid_2d, &
ri_data%bsizes_RI_fit, ri_data%bsizes_RI_fit, &
name="(RI | RI)")
CALL dbcsr_t_create(ri_data%t_2c_int(1, 1), ri_data%t_2c_int(2, 1))
DEALLOCATE (dist1, dist2)
@ -1370,7 +1381,7 @@ CONTAINS
! is larger than this (it is however not a problem for load balancing if actual MO dimension
! is slightly smaller)
MO_dim = MAX((ri_data%nelectron_total/2 - 1)/ri_data%n_mem + 1, 1)
MO_dim = (MO_dim - 1)/ri_data%min_bsize_MO + 1
MO_dim = (MO_dim - 1)/ri_data%max_bsize_MO + 1
pdims = 0
CALL dbcsr_t_mp_dims_create(nproc, pdims, [SIZE(ri_data%bsizes_AO_split), SIZE(ri_data%bsizes_RI_split), MO_dim])
@ -1391,6 +1402,19 @@ CONTAINS
ENDIF
!For forces
ALLOCATE (ri_data%t_2c_inv(1, 1))
CALL create_2c_tensor(ri_data%t_2c_inv(1, 1), dist1, dist2, ri_data%pgrid_2d, &
ri_data%bsizes_RI_split, ri_data%bsizes_RI_split, &
name="(RI | RI)")
DEALLOCATE (dist1, dist2)
ALLOCATE (ri_data%t_2c_pot(1, 1))
CALL create_2c_tensor(ri_data%t_2c_pot(1, 1), dist1, dist2, ri_data%pgrid_2d, &
ri_data%bsizes_RI_split, ri_data%bsizes_RI_split, &
name="(RI | RI)")
DEALLOCATE (dist1, dist2)
ri_data%dbcsr_nflop = 0
ri_data%dbcsr_time = 0.0_dp
ri_data%num_pe = para_env%num_pe
@ -1499,15 +1523,14 @@ CONTAINS
DEALLOCATE (ri_data%t_3c_int_ctr_2)
DO ispin = 1, SIZE(ri_data%t_3c_int_mo, 1)
CALL dbcsr_t_batched_contract_finalize(ri_data%t_2c_int(ispin, 1))
CALL dbcsr_t_batched_contract_finalize(ri_data%t_3c_int_mo(ispin, 1, 1))
CALL dbcsr_t_batched_contract_finalize(ri_data%t_3c_ctr_RI(ispin, 1, 1))
CALL dbcsr_t_destroy(ri_data%t_2c_int(ispin, 1))
CALL dbcsr_t_destroy(ri_data%t_3c_int_mo(ispin, 1, 1))
CALL dbcsr_t_destroy(ri_data%t_3c_ctr_RI(ispin, 1, 1))
CALL dbcsr_t_destroy(ri_data%t_3c_ctr_KS(ispin, 1, 1))
CALL dbcsr_t_destroy(ri_data%t_3c_ctr_KS_copy(ispin, 1, 1))
ENDDO
DO ispin = 1, 2
CALL dbcsr_t_destroy(ri_data%t_2c_int(ispin, 1))
ENDDO
DEALLOCATE (ri_data%t_2c_int)
DEALLOCATE (ri_data%t_3c_int_mo)
DEALLOCATE (ri_data%t_3c_ctr_RI)
@ -1515,6 +1538,11 @@ CONTAINS
DEALLOCATE (ri_data%t_3c_ctr_KS_copy)
ENDIF
CALL dbcsr_t_destroy(ri_data%t_2c_inv(1, 1))
DEALLOCATE (ri_data%t_2c_inv)
CALL dbcsr_t_destroy(ri_data%t_2c_pot(1, 1))
DEALLOCATE (ri_data%t_2c_pot)
CALL timestop(handle)
END SUBROUTINE
@ -2633,10 +2661,12 @@ CONTAINS
"HFX_RI_INFO| EPS_PGF_ORB", ri_data%eps_pgf_orb
WRITE (UNIT=iw, FMT="((T3, A, T73, ES8.1))") &
"HFX_RI_INFO| EPS_SCHWARZ: ", ri_data%eps_schwarz
WRITE (UNIT=iw, FMT="((T3, A, T73, ES8.1))") &
"HFX_RI_INFO| EPS_SCHWARZ_FORCES: ", ri_data%eps_schwarz_forces
WRITE (UNIT=iw, FMT="(T3, A, T78, I3)") &
"HFX_RI_INFO| Minimum block size", ri_data%min_bsize
WRITE (UNIT=iw, FMT="(T3, A, T78, I3)") &
"HFX_RI_INFO| MO block size", ri_data%min_bsize_MO
"HFX_RI_INFO| MO block size", ri_data%max_bsize_MO
SELECT CASE (ri_data%flavor)
CASE (ri_mo)
WRITE (UNIT=iw, FMT="(T3, A, T79, I2)") &

4
src/input_cp2k_hfx.F
View file

@ -633,8 +633,8 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="MIN_BLOCK_SIZE_MO", &
description="Minimum tensor block size for MOs.", &
CALL keyword_create(keyword, __LOCATION__, name="MAX_BLOCK_SIZE_MO", &
description="Maximum tensor block size for MOs.", &
default_i_val=64)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

594
src/libint_2c_3c.F
View file

@ -20,9 +20,12 @@ MODULE libint_2c_3c
do_potential_truncated
USE kinds, ONLY: default_path_length,&
dp
USE libint_wrapper, ONLY: cp_libint_get_2eris,&
USE libint_wrapper, ONLY: cp_libint_get_2eri_derivs,&
cp_libint_get_2eris,&
cp_libint_get_3eri_derivs,&
cp_libint_get_3eris,&
cp_libint_set_params_eri,&
cp_libint_set_params_eri_deriv,&
cp_libint_t,&
prim_data_f_size
USE mathconstants, ONLY: pi
@ -37,7 +40,8 @@ MODULE libint_2c_3c
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'libint_2c_3c'
PUBLIC :: eri_3center, eri_2center, cutoff_screen_factor, libint_potential_type, compare_potential_types
PUBLIC :: eri_3center, eri_2center, cutoff_screen_factor, libint_potential_type, &
eri_3center_derivs, eri_2center_derivs, compare_potential_types
! For screening of integrals with a truncated potential, it is important to use a slightly larger
! cutoff radius due to the discontinuity of the truncated Coulomb potential at the cutoff radius.
@ -119,7 +123,7 @@ CONTAINS
REAL(KIND=dp), INTENT(IN) :: dab, dac, dbc
TYPE(cp_libint_t), INTENT(INOUT) :: lib
TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
REAL(dp), INTENT(OUT), OPTIONAL :: int_abc_ext
REAL(dp), INTENT(INOUT), OPTIONAL :: int_abc_ext
INTEGER :: a_mysize(1), a_offset, a_start, b_offset, b_start, c_offset, c_start, i, ipgf, j, &
jpgf, k, kpgf, li, lj, lk, ncoa, ncob, ncoc, op, p1, p2, p3
@ -175,7 +179,7 @@ CONTAINS
c_start = (kpgf - 1)*ncoset(lc_max)
!start with all the (c|ba) integrals (standard order) and keep to lb >= la
CALL set_params_3c(lib, ra, rb, rc, la_max, lb_max, lc_max, zeti, zetj, zetk, &
CALL set_params_3c(lib, ra, rb, rc, zeti, zetj, zetk, la_max, lb_max, lc_max, &
potential_parameter=potential_parameter, params_out=params)
DO li = la_min, la_max
@ -279,12 +283,12 @@ CONTAINS
!> \param ri ...
!> \param rj ...
!> \param rk ...
!> \param li_max ...
!> \param lj_max ...
!> \param lk_max ...
!> \param zeti ...
!> \param zetj ...
!> \param zetk ...
!> \param li_max ...
!> \param lj_max ...
!> \param lk_max ...
!> \param potential_parameter ...
!> \param params_in external parameters to use for libint
!> \param params_out returns the libint parameters computed based on the other arguments
@ -292,13 +296,13 @@ CONTAINS
!> centers 3 and 4 because of angular momenta and pretty much all the parameters of libint
!> remain the same upon such a change => might avoid recomputing things over and over again
! **************************************************************************************************
SUBROUTINE set_params_3c(lib, ri, rj, rk, li_max, lj_max, lk_max, zeti, zetj, zetk, &
SUBROUTINE set_params_3c(lib, ri, rj, rk, zeti, zetj, zetk, li_max, lj_max, lk_max, &
potential_parameter, params_in, params_out)
TYPE(cp_libint_t), INTENT(INOUT) :: lib
REAL(dp), DIMENSION(3), INTENT(IN) :: ri, rj, rk
INTEGER, INTENT(IN), OPTIONAL :: li_max, lj_max, lk_max
REAL(dp), INTENT(IN), OPTIONAL :: zeti, zetj, zetk
INTEGER, INTENT(IN), OPTIONAL :: li_max, lj_max, lk_max
TYPE(libint_potential_type), INTENT(IN), OPTIONAL :: potential_parameter
TYPE(params_3c), OPTIONAL, POINTER :: params_in, params_out
@ -393,6 +397,373 @@ CONTAINS
END SUBROUTINE set_params_3c
! **************************************************************************************************
!> \brief Computes the derivatives of the 3-center electron repulsion integrals (ab|c) for a given
!> set of cartesian gaussian orbitals. Returns x,y,z derivatives for 1st and 2nd center
!> \param der_abc_1 the derivatives for the 1st center (allocated before hand)
!> \param der_abc_2 the derivatives for the 2nd center (allocated before hand)
!> \param la_min ...
!> \param la_max ...
!> \param npgfa ...
!> \param zeta ...
!> \param rpgfa ...
!> \param ra ...
!> \param lb_min ...
!> \param lb_max ...
!> \param npgfb ...
!> \param zetb ...
!> \param rpgfb ...
!> \param rb ...
!> \param lc_min ...
!> \param lc_max ...
!> \param npgfc ...
!> \param zetc ...
!> \param rpgfc ...
!> \param rc ...
!> \param dab ...
!> \param dac ...
!> \param dbc ...
!> \param lib the libint_t object for evaluation (assume that it is initialized outside)
!> \param potential_parameter the info about the potential
!> \param der_abc_1_ext the extremal value of der_abc_1, i.e., MAXVAL(ABS(der_abc_1))
!> \param der_abc_2_ext ...
!> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
!> the libint library must be static initialized, and in case of truncated Coulomb operator,
!> the latter must be initialized too. Note that the derivative wrt to the third center
!> can be obtained via translational invariance
! **************************************************************************************************
SUBROUTINE eri_3center_derivs(der_abc_1, der_abc_2, &
la_min, la_max, npgfa, zeta, rpgfa, ra, &
lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
lc_min, lc_max, npgfc, zetc, rpgfc, rc, &
dab, dac, dbc, lib, potential_parameter, &
der_abc_1_ext, der_abc_2_ext)
REAL(dp), DIMENSION(:, :, :, :), INTENT(INOUT) :: der_abc_1, der_abc_2
INTEGER, INTENT(IN) :: la_min, la_max, npgfa
REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
REAL(dp), DIMENSION(3), INTENT(IN) :: ra
INTEGER, INTENT(IN) :: lb_min, lb_max, npgfb
REAL(dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
REAL(dp), DIMENSION(3), INTENT(IN) :: rb
INTEGER, INTENT(IN) :: lc_min, lc_max, npgfc
REAL(dp), DIMENSION(:), INTENT(IN) :: zetc, rpgfc
REAL(dp), DIMENSION(3), INTENT(IN) :: rc
REAL(KIND=dp), INTENT(IN) :: dab, dac, dbc
TYPE(cp_libint_t), INTENT(INOUT) :: lib
TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
REAL(dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: der_abc_1_ext, der_abc_2_ext
INTEGER :: a_mysize(1), a_offset, a_start, b_offset, b_start, c_offset, c_start, i, i_deriv, &
ipgf, j, jpgf, k, kpgf, li, lj, lk, ncoa, ncob, ncoc, op, p1, p2, p3
INTEGER, DIMENSION(3) :: permute_1, permute_2
LOGICAL :: do_ext
REAL(dp) :: dr_ab, dr_ac, dr_bc, zeti, zetj, zetk
REAL(dp), DIMENSION(3) :: der_abc_1_ext_prv, der_abc_2_ext_prv
REAL(dp), DIMENSION(:, :), POINTER :: p_deriv
TYPE(params_3c), POINTER :: params
NULLIFY (params, p_deriv)
ALLOCATE (params)
permute_1 = [4, 5, 6]
permute_2 = [7, 8, 9]
dr_ab = 0.0_dp
dr_bc = 0.0_dp
dr_ac = 0.0_dp
op = potential_parameter%potential_type
IF (op == do_potential_truncated .OR. op == do_potential_short) THEN
dr_bc = potential_parameter%cutoff_radius*cutoff_screen_factor
dr_ac = potential_parameter%cutoff_radius*cutoff_screen_factor
ELSEIF (op == do_potential_coulomb) THEN
dr_bc = 1000000.0_dp
dr_ac = 1000000.0_dp
ENDIF
do_ext = .FALSE.
IF (PRESENT(der_abc_1_ext) .OR. PRESENT(der_abc_2_ext)) do_ext = .TRUE.
der_abc_1_ext_prv = 0.0_dp
der_abc_2_ext_prv = 0.0_dp
!Note: we want to compute all possible integrals based on the 3-centers (ab|c) before
! having to switch to (ba|c) (or the other way around) due to angular momenta in libint
! For a triplet of centers (k|ji), we can only compute integrals for which lj >= li
!Looping over the pgfs
DO ipgf = 1, npgfa
zeti = zeta(ipgf)
a_start = (ipgf - 1)*ncoset(la_max)
DO jpgf = 1, npgfb
! screening
IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
zetj = zetb(jpgf)
b_start = (jpgf - 1)*ncoset(lb_max)
DO kpgf = 1, npgfc
! screening
IF (rpgfb(jpgf) + rpgfc(kpgf) + dr_bc < dbc) CYCLE
IF (rpgfa(ipgf) + rpgfc(kpgf) + dr_ac < dac) CYCLE
zetk = zetc(kpgf)
c_start = (kpgf - 1)*ncoset(lc_max)
!start with all the (c|ba) integrals (standard order) and keep to lb >= la
CALL set_params_3c_deriv(lib, ra, rb, rc, zeti, zetj, zetk, la_max, lb_max, lc_max, &
potential_parameter=potential_parameter, params_out=params)
DO li = la_min, la_max
a_offset = a_start + ncoset(li - 1)
ncoa = nco(li)
DO lj = MAX(li, lb_min), lb_max
b_offset = b_start + ncoset(lj - 1)
ncob = nco(lj)
DO lk = lc_min, lc_max
c_offset = c_start + ncoset(lk - 1)
ncoc = nco(lk)
a_mysize(1) = ncoa*ncob*ncoc
CALL cp_libint_get_3eri_derivs(li, lj, lk, lib, p_deriv, a_mysize)
IF (do_ext) THEN
DO k = 1, ncoc
p1 = (k - 1)*ncob
DO j = 1, ncob
p2 = (p1 + j - 1)*ncoa
DO i = 1, ncoa
p3 = p2 + i
DO i_deriv = 1, 3
der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
p_deriv(p3, permute_2(i_deriv))
der_abc_1_ext_prv(i_deriv) = MAX(der_abc_1_ext_prv(i_deriv), &
ABS(p_deriv(p3, permute_2(i_deriv))))
der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
p_deriv(p3, permute_1(i_deriv))
der_abc_2_ext_prv(i_deriv) = MAX(der_abc_2_ext_prv(i_deriv), &
ABS(p_deriv(p3, permute_1(i_deriv))))
END DO
END DO
END DO
END DO
ELSE
DO k = 1, ncoc
p1 = (k - 1)*ncob
DO j = 1, ncob
p2 = (p1 + j - 1)*ncoa
DO i = 1, ncoa
p3 = p2 + i
DO i_deriv = 1, 3
der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
p_deriv(p3, permute_2(i_deriv))
der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
p_deriv(p3, permute_1(i_deriv))
END DO
END DO
END DO
END DO
ENDIF
DEALLOCATE (p_deriv)
END DO !lk
END DO !lj
END DO !li
!swap centers 3 and 4 to compute (c|ab) with lb < la
CALL set_params_3c_deriv(lib, rb, ra, rc, zetj, zeti, zetk, params_in=params)
DO lj = lb_min, lb_max
b_offset = b_start + ncoset(lj - 1)
ncob = nco(lj)
DO li = MAX(lj + 1, la_min), la_max
a_offset = a_start + ncoset(li - 1)
ncoa = nco(li)
DO lk = lc_min, lc_max
c_offset = c_start + ncoset(lk - 1)
ncoc = nco(lk)
a_mysize(1) = ncoa*ncob*ncoc
CALL cp_libint_get_3eri_derivs(lj, li, lk, lib, p_deriv, a_mysize)
IF (do_ext) THEN
DO k = 1, ncoc
p1 = (k - 1)*ncoa
DO i = 1, ncoa
p2 = (p1 + i - 1)*ncob
DO j = 1, ncob
p3 = p2 + j
DO i_deriv = 1, 3
der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
p_deriv(p3, permute_1(i_deriv))
der_abc_1_ext_prv(i_deriv) = MAX(der_abc_1_ext_prv(i_deriv), &
ABS(p_deriv(p3, permute_1(i_deriv))))
der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
p_deriv(p3, permute_2(i_deriv))
der_abc_2_ext_prv(i_deriv) = MAX(der_abc_2_ext_prv(i_deriv), &
ABS(p_deriv(p3, permute_2(i_deriv))))
END DO
END DO
END DO
END DO
ELSE
DO k = 1, ncoc
p1 = (k - 1)*ncoa
DO i = 1, ncoa
p2 = (p1 + i - 1)*ncob
DO j = 1, ncob
p3 = p2 + j
DO i_deriv = 1, 3
der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
p_deriv(p3, permute_1(i_deriv))
der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
p_deriv(p3, permute_2(i_deriv))
END DO
END DO
END DO
END DO
ENDIF
DEALLOCATE (p_deriv)
END DO !lk
END DO !li
END DO !lj
END DO !kpgf
END DO !jpgf
END DO !ipgf
IF (PRESENT(der_abc_1_ext)) der_abc_1_ext = der_abc_1_ext_prv
IF (PRESENT(der_abc_2_ext)) der_abc_2_ext = der_abc_2_ext_prv
DEALLOCATE (params)
END SUBROUTINE eri_3center_derivs
! **************************************************************************************************
!> \brief Sets the internals of the cp_libint_t object for derivatives of integrals of type (k|ji)
!> \param lib ..
!> \param ri ...
!> \param rj ...
!> \param rk ...
!> \param zeti ...
!> \param zetj ...
!> \param zetk ...
!> \param li_max ...
!> \param lj_max ...
!> \param lk_max ...
!> \param potential_parameter ...
!> \param params_in ...
!> \param params_out ...
!> \note The use of params_in and params_out comes from the fact that one might have to swap
!> centers 3 and 4 because of angular momenta and pretty much all the parameters of libint
!> remain the same upon such a change => might avoid recomputing things over and over again
! **************************************************************************************************
SUBROUTINE set_params_3c_deriv(lib, ri, rj, rk, zeti, zetj, zetk, li_max, lj_max, lk_max, &
potential_parameter, params_in, params_out)
TYPE(cp_libint_t), INTENT(INOUT) :: lib
REAL(dp), DIMENSION(3), INTENT(IN) :: ri, rj, rk
REAL(dp), INTENT(IN) :: zeti, zetj, zetk
INTEGER, INTENT(IN), OPTIONAL :: li_max, lj_max, lk_max
TYPE(libint_potential_type), INTENT(IN), OPTIONAL :: potential_parameter
TYPE(params_3c), OPTIONAL, POINTER :: params_in, params_out
INTEGER :: l
LOGICAL :: use_gamma
REAL(dp) :: gammaq, omega2, omega_corr, omega_corr2, &
prefac, R, S1234, T, tmp
REAL(dp), ALLOCATABLE, DIMENSION(:) :: Fm
TYPE(params_3c), POINTER :: params
IF (PRESENT(params_in)) THEN
params => params_in
ELSE
params => params_out
params%m_max = li_max + lj_max + lk_max + 1
gammaq = zeti + zetj
params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/gammaq
params%ZetapEtaInv = 1._dp/(zetk + gammaq)
params%Q = (zeti*ri + zetj*rj)*params%EtaInv
params%W = (zetk*rk + gammaq*params%Q)*params%ZetapEtaInv
params%Rho = zetk*gammaq/(zetk + gammaq)
params%Fm = 0.0_dp
SELECT CASE (potential_parameter%potential_type)
CASE (do_potential_coulomb)
T = params%Rho*SUM((params%Q - rk)**2)
S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
CALL fgamma(params%m_max, T, params%Fm)
params%Fm = prefac*params%Fm
CASE (do_potential_truncated)
R = potential_parameter%cutoff_radius*SQRT(params%Rho)
T = params%Rho*SUM((params%Q - rk)**2)
S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
params%Fm = prefac*params%Fm
CASE (do_potential_short)
T = params%Rho*SUM((params%Q - rk)**2)
S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
CALL fgamma(params%m_max, T, params%Fm)
omega2 = potential_parameter%omega**2
omega_corr2 = omega2/(omega2 + params%Rho)
omega_corr = SQRT(omega_corr2)
T = T*omega_corr2
ALLOCATE (Fm(prim_data_f_size))
CALL fgamma(params%m_max, T, Fm)
tmp = -omega_corr
DO l = 1, params%m_max + 1
params%Fm(l) = params%Fm(l) + Fm(l)*tmp
tmp = tmp*omega_corr2
END DO
params%Fm = prefac*params%Fm
CASE (do_potential_id)
S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2) &
- gammaq*zetk*params%ZetapEtaInv*SUM((params%Q - rk)**2))
prefac = SQRT((pi*params%ZetapEtaInv)**3)*S1234
params%Fm(:) = prefac
CASE DEFAULT
CPABORT("Requested operator NYI")
END SELECT
END IF
CALL cp_libint_set_params_eri_deriv(lib, rk, rk, rj, ri, rk, &
params%Q, params%W, zetk, 0.0_dp, zetj, zeti, params%ZetaInv, &
params%EtaInv, params%ZetapEtaInv, params%Rho, params%m_max, params%Fm)
END SUBROUTINE set_params_3c_deriv
! **************************************************************************************************
!> \brief Computes the 2-center electron repulsion integrals (a|b) for a given set of cartesian
!> gaussian orbitals
@ -460,7 +831,7 @@ CONTAINS
!screening
IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
CALL set_params_2c(lib, ra, rb, la_max, lb_max, zeti, zetj, potential_parameter)
CALL set_params_2c(lib, ra, rb, zeti, zetj, la_max, lb_max, potential_parameter)
DO li = la_min, la_max
a_offset = a_start + ncoset(li - 1)
@ -486,25 +857,25 @@ CONTAINS
END DO
END DO
END SUBROUTINE
END SUBROUTINE eri_2center
! **************************************************************************************************
!> \brief Sets the internals of the cp_libint_t object for integrals of type (k|j)
!> \param lib ..
!> \param rj ...
!> \param rk ...
!> \param lj_max ...
!> \param lk_max ...
!> \param zetj ...
!> \param zetk ...
!> \param lj_max ...
!> \param lk_max ...
!> \param potential_parameter ...
! **************************************************************************************************
SUBROUTINE set_params_2c(lib, rj, rk, lj_max, lk_max, zetj, zetk, potential_parameter)
SUBROUTINE set_params_2c(lib, rj, rk, zetj, zetk, lj_max, lk_max, potential_parameter)
TYPE(cp_libint_t), INTENT(INOUT) :: lib
REAL(dp), DIMENSION(3), INTENT(IN) :: rj, rk
INTEGER, INTENT(IN) :: lj_max, lk_max
REAL(dp), INTENT(IN) :: zetj, zetk
INTEGER, INTENT(IN) :: lj_max, lk_max
TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
INTEGER :: l, op
@ -604,5 +975,198 @@ CONTAINS
END FUNCTION compare_potential_types
!> \brief Computes the 2-center derivatives of the electron repulsion integrals (a|b) for a given
!> set of cartesian gaussian orbitals. Returns the derivatives wrt to the first center
!> \param der_ab the derivatives as array of cartesian orbitals (allocated before hand)
!> \param la_min ...
!> \param la_max ...
!> \param npgfa ...
!> \param zeta ...
!> \param rpgfa ...
!> \param ra ...
!> \param lb_min ...
!> \param lb_max ...
!> \param npgfb ...
!> \param zetb ...
!> \param rpgfb ...
!> \param rb ...
!> \param dab ...
!> \param lib the libint_t object for evaluation (assume that it is initialized outside)
!> \param potential_parameter the info about the potential
!> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
!> the libint library must be static initialized, and in case of truncated Coulomb operator,
!> the latter must be initialized too
! **************************************************************************************************
SUBROUTINE eri_2center_derivs(der_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
dab, lib, potential_parameter)
REAL(dp), DIMENSION(:, :, :), INTENT(INOUT) :: der_ab
INTEGER, INTENT(IN) :: la_min, la_max, npgfa
REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
REAL(dp), DIMENSION(3), INTENT(IN) :: ra
INTEGER, INTENT(IN) :: lb_min, lb_max, npgfb
REAL(dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
REAL(dp), DIMENSION(3), INTENT(IN) :: rb
REAL(dp), INTENT(IN) :: dab
TYPE(cp_libint_t), INTENT(INOUT) :: lib
TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
INTEGER :: a_mysize(1), a_offset, a_start, &
b_offset, b_start, i, i_deriv, ipgf, &
j, jpgf, li, lj, ncoa, ncob, p1, p2
INTEGER, DIMENSION(3) :: permute
REAL(dp) :: dr_ab, zeti, zetj
REAL(dp), DIMENSION(:, :), POINTER :: p_deriv
NULLIFY (p_deriv)
permute = [4, 5, 6]
dr_ab = 0.0_dp
IF (potential_parameter%potential_type == do_potential_truncated .OR. &
potential_parameter%potential_type == do_potential_short) THEN
dr_ab = potential_parameter%cutoff_radius*cutoff_screen_factor
ELSEIF (potential_parameter%potential_type == do_potential_coulomb) THEN
dr_ab = 1000000.0_dp
ENDIF
!Looping over the pgfs
DO ipgf = 1, npgfa
zeti = zeta(ipgf)
a_start = (ipgf - 1)*ncoset(la_max)
DO jpgf = 1, npgfb
zetj = zetb(jpgf)
b_start = (jpgf - 1)*ncoset(lb_max)
!screening
IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
CALL set_params_2c_deriv(lib, ra, rb, zeti, zetj, la_max, lb_max, potential_parameter)
DO li = la_min, la_max
a_offset = a_start + ncoset(li - 1)
ncoa = nco(li)
DO lj = lb_min, lb_max
b_offset = b_start + ncoset(lj - 1)
ncob = nco(lj)
a_mysize(1) = ncoa*ncob
CALL cp_libint_get_2eri_derivs(li, lj, lib, p_deriv, a_mysize)
DO j = 1, ncob
p1 = (j - 1)*ncoa
DO i = 1, ncoa
p2 = p1 + i
DO i_deriv = 1, 3
der_ab(a_offset + i, b_offset + j, i_deriv) = p_deriv(p2, permute(i_deriv))
END DO
END DO
END DO
DEALLOCATE (p_deriv)
END DO
END DO
END DO
END DO
END SUBROUTINE eri_2center_derivs
! **************************************************************************************************
!> \brief Sets the internals of the cp_libint_t object for derivatives of integrals of type (k|j)
!> \param lib ..
!> \param rj ...
!> \param rk ...
!> \param zetj ...
!> \param zetk ...
!> \param lj_max ...
!> \param lk_max ...
!> \param potential_parameter ...
! **************************************************************************************************
SUBROUTINE set_params_2c_deriv(lib, rj, rk, zetj, zetk, lj_max, lk_max, potential_parameter)
TYPE(cp_libint_t), INTENT(INOUT) :: lib
REAL(dp), DIMENSION(3), INTENT(IN) :: rj, rk
REAL(dp), INTENT(IN) :: zetj, zetk
INTEGER, INTENT(IN) :: lj_max, lk_max
TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
INTEGER :: l, op
LOGICAL :: use_gamma
REAL(dp) :: omega2, omega_corr, omega_corr2, prefac, &
R, T, tmp
REAL(dp), ALLOCATABLE, DIMENSION(:) :: Fm
TYPE(params_2c) :: params
!The internal structure of libint2 is based on 4-center integrals
!For 2-center, two of those are dummy centers
!The integral is assumed to be (k|j) where the centers are ordered as:
!k -> 1, j -> 3 and (the centers #2 & #4 are dummy centers)
!Note: some variable of 4-center integrals simplify due to dummy centers:
! P -> rk, gammap -> zetk
! Q -> rj, gammaq -> zetj
op = potential_parameter%potential_type
params%m_max = lj_max + lk_max + 1
params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/zetj
params%ZetapEtaInv = 1._dp/(zetk + zetj)
params%W = (zetk*rk + zetj*rj)*params%ZetapEtaInv
params%Rho = zetk*zetj/(zetk + zetj)
params%Fm = 0.0_dp
SELECT CASE (op)
CASE (do_potential_coulomb)
T = params%Rho*SUM((rj - rk)**2)
prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
CALL fgamma(params%m_max, T, params%Fm)
params%Fm = prefac*params%Fm
CASE (do_potential_truncated)
R = potential_parameter%cutoff_radius*SQRT(params%Rho)
T = params%Rho*SUM((rj - rk)**2)
prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
params%Fm = prefac*params%Fm
CASE (do_potential_short)
T = params%Rho*SUM((rj - rk)**2)
prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
CALL fgamma(params%m_max, T, params%Fm)
omega2 = potential_parameter%omega**2
omega_corr2 = omega2/(omega2 + params%Rho)
omega_corr = SQRT(omega_corr2)
T = T*omega_corr2
ALLOCATE (Fm(prim_data_f_size))
CALL fgamma(params%m_max, T, Fm)
tmp = -omega_corr
DO l = 1, params%m_max + 1
params%Fm(l) = params%Fm(l) + Fm(l)*tmp
tmp = tmp*omega_corr2
END DO
params%Fm = prefac*params%Fm
CASE (do_potential_id)
prefac = SQRT((pi*params%ZetapEtaInv)**3)*EXP(-zetj*zetk*params%ZetapEtaInv*SUM((rk - rj)**2))
params%Fm(:) = prefac
CASE DEFAULT
CPABORT("Requested operator NYI")
END SELECT
CALL cp_libint_set_params_eri_deriv(lib, rk, rk, rj, rj, rk, rj, params%W, zetk, 0.0_dp, &
zetj, 0.0_dp, params%ZetaInv, params%EtaInv, &
params%ZetapEtaInv, params%Rho, &
params%m_max, params%Fm)
END SUBROUTINE set_params_2c_deriv
END MODULE libint_2c_3c

136
src/libint_wrapper.F
View file

@ -33,7 +33,8 @@ MODULE libint_wrapper
libint2_cleanup_eri1, libint2_init_eri, libint2_init_eri1, libint2_static_cleanup, &
libint2_static_init, libint_t, libint2_max_am_eri, libint2_init_3eri, libint2_cleanup_3eri, &
libint2_init_2eri, libint2_cleanup_2eri, &
libint2_build_2eri, libint2_build_3eri
libint2_build_2eri, libint2_build_3eri, libint2_build_3eri1, libint2_cleanup_3eri1, libint2_init_3eri1, &
libint2_build_2eri1, libint2_cleanup_2eri1, libint2_init_2eri1
#endif
USE orbital_pointers, ONLY: nco
#include "./base/base_uses.f90"
@ -47,7 +48,9 @@ MODULE libint_wrapper
get_ssss_f_val, cp_libint_set_contrdepth, cp_libint_set_params_eri_screen, &
cp_libint_set_params_eri, cp_libint_set_params_eri_deriv, &
cp_libint_init_3eri, cp_libint_cleanup_3eri, cp_libint_get_3eris, &
cp_libint_init_2eri, cp_libint_cleanup_2eri, cp_libint_get_2eris
cp_libint_init_2eri, cp_libint_cleanup_2eri, cp_libint_get_2eris, &
cp_libint_get_3eri_derivs, cp_libint_init_3eri1, cp_libint_cleanup_3eri1, &
cp_libint_get_2eri_derivs, cp_libint_init_2eri1, cp_libint_cleanup_2eri1
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'libint_wrapper'
@ -377,6 +380,91 @@ CONTAINS
END SUBROUTINE cp_libint_get_3eris
! **************************************************************************************************
!> \brief ...
!> \param n_c ...
!> \param n_b ...
!> \param n_a ...
!> \param lib ...
!> \param p_work ...
!> \param a_mysize ...
! **************************************************************************************************
SUBROUTINE cp_libint_get_3eri_derivs(n_c, n_b, n_a, lib, p_work, a_mysize)
INTEGER, INTENT(IN) :: n_c, n_b, n_a
TYPE(cp_libint_t) :: lib
INTEGER :: a_mysize(1)
REAL(dp), DIMENSION(:, :), POINTER :: p_work
REAL(dp), DIMENSION(:), POINTER :: p_work_tmp
#if(__LIBINT)
PROCEDURE(libint2_build), POINTER :: pbuild
INTEGER :: i
CALL C_F_PROCPOINTER(libint2_build_3eri1(n_c, n_b, n_a), pbuild)
CALL pbuild(lib%prv)
ALLOCATE (p_work(a_mysize(1), 9))
!Derivatives 1-3 can be obtained using translational invariance
DO i = 4, 9
NULLIFY (p_work_tmp)
CALL C_F_POINTER(lib%prv(1)%targets(i), p_work_tmp, SHAPE=a_mysize)
p_work(:, i) = p_work_tmp
ENDDO
#else
MARK_USED(n_c)
MARK_USED(n_b)
MARK_USED(n_a)
MARK_USED(lib)
MARK_USED(p_work)
MARK_USED(a_mysize)
CPABORT("This CP2K executable has not been linked against the required library libint.")
#endif
END SUBROUTINE cp_libint_get_3eri_derivs
! **************************************************************************************************
!> \brief ...
!> \param n_c ...
!> \param n_b ...
!> \param n_a ...
!> \param lib ...
!> \param p_work ...
!> \param a_mysize ...
! **************************************************************************************************
SUBROUTINE cp_libint_get_2eri_derivs(n_b, n_a, lib, p_work, a_mysize)
INTEGER, INTENT(IN) :: n_b, n_a
TYPE(cp_libint_t) :: lib
INTEGER :: a_mysize(1)
REAL(dp), DIMENSION(:, :), POINTER :: p_work
REAL(dp), DIMENSION(:), POINTER :: p_work_tmp
#if(__LIBINT)
PROCEDURE(libint2_build), POINTER :: pbuild
INTEGER :: i
CALL C_F_PROCPOINTER(libint2_build_2eri1(n_b, n_a), pbuild)
CALL pbuild(lib%prv)
ALLOCATE (p_work(a_mysize(1), 6))
!Derivatives 1-3 can be obtained using translational invariance
DO i = 4, 6
NULLIFY (p_work_tmp)
CALL C_F_POINTER(lib%prv(1)%targets(i), p_work_tmp, SHAPE=a_mysize)
p_work(:, i) = p_work_tmp
ENDDO
#else
MARK_USED(n_b)
MARK_USED(n_a)
MARK_USED(lib)
MARK_USED(p_work)
MARK_USED(a_mysize)
CPABORT("This CP2K executable has not been linked against the required library libint.")
#endif
END SUBROUTINE cp_libint_get_2eri_derivs
! **************************************************************************************************
!> \brief ...
!> \param n_c ...
@ -528,6 +616,30 @@ CONTAINS
#endif
END SUBROUTINE
SUBROUTINE cp_libint_init_3eri1(lib, max_am)
TYPE(cp_libint_t) :: lib
INTEGER :: max_am
#if(__LIBINT)
CALL libint2_init_3eri1(lib%prv, max_am, C_NULL_PTR)
#else
MARK_USED(lib)
MARK_USED(max_am)
CPABORT("This CP2K executable has not been linked against the required library libint.")
#endif
END SUBROUTINE
SUBROUTINE cp_libint_init_2eri1(lib, max_am)
TYPE(cp_libint_t) :: lib
INTEGER :: max_am
#if(__LIBINT)
CALL libint2_init_2eri1(lib%prv, max_am, C_NULL_PTR)
#else
MARK_USED(lib)
MARK_USED(max_am)
CPABORT("This CP2K executable has not been linked against the required library libint.")
#endif
END SUBROUTINE
SUBROUTINE cp_libint_init_2eri(lib, max_am)
TYPE(cp_libint_t) :: lib
INTEGER :: max_am
@ -570,6 +682,26 @@ CONTAINS
#endif
END SUBROUTINE
SUBROUTINE cp_libint_cleanup_3eri1(lib)
TYPE(cp_libint_t) :: lib
#if(__LIBINT)
CALL libint2_cleanup_3eri1(lib%prv)
#else
MARK_USED(lib)
CPABORT("This CP2K executable has not been linked against the required library libint.")
#endif
END SUBROUTINE
SUBROUTINE cp_libint_cleanup_2eri1(lib)
TYPE(cp_libint_t) :: lib
#if(__LIBINT)
CALL libint2_cleanup_2eri1(lib%prv)
#else
MARK_USED(lib)
CPABORT("This CP2K executable has not been linked against the required library libint.")
#endif
END SUBROUTINE
SUBROUTINE cp_libint_cleanup_2eri(lib)
TYPE(cp_libint_t) :: lib
#if(__LIBINT)

2
src/mp2.F
View file

@ -450,6 +450,7 @@ CONTAINS
DO irep = 1, n_rep_hf
DO i_thread = 0, n_threads - 1
actual_x_data => qs_env%x_data(irep, i_thread + 1)
IF (actual_x_data%do_hfx_ri) CYCLE
do_dynamic_load_balancing = .TRUE.
IF (n_threads == 1 .OR. actual_x_data%memory_parameter%do_disk_storage) do_dynamic_load_balancing = .FALSE.
@ -696,6 +697,7 @@ CONTAINS
DO irep = 1, n_rep_hf
DO i_thread = 0, n_threads - 1
actual_x_data => qs_env%x_data(irep, i_thread + 1)
IF (actual_x_data%do_hfx_ri) CYCLE
do_dynamic_load_balancing = .TRUE.
IF (n_threads == 1 .OR. actual_x_data%memory_parameter%do_disk_storage) do_dynamic_load_balancing = .FALSE.

25
src/mp2_cphf.F
View file

@ -48,12 +48,14 @@ MODULE mp2_cphf
dbcsr_set
USE hfx_admm_utils, ONLY: tddft_hfx_matrix
USE hfx_derivatives, ONLY: derivatives_four_center
USE hfx_ri, ONLY: hfx_ri_update_forces
USE hfx_types, ONLY: alloc_containers,&
hfx_container_type,&
hfx_init_container,&
hfx_type
USE input_constants, ONLY: do_admm_aux_exch_func_none,&
ot_precond_full_all
USE input_cp2k_hfx, ONLY: ri_mo
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type
@ -282,6 +284,7 @@ CONTAINS
DO irep = 1, n_rep_hf
DO i_thread = 0, n_threads - 1
actual_x_data => qs_env%x_data(irep, i_thread + 1)
IF (actual_x_data%do_hfx_ri) CYCLE
do_dynamic_load_balancing = .TRUE.
IF (n_threads == 1 .OR. actual_x_data%memory_parameter%do_disk_storage) do_dynamic_load_balancing = .FALSE.
@ -970,6 +973,7 @@ CONTAINS
rho_ao, rho_ao_aux, scrm
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp, scrm_kp
TYPE(dft_control_type), POINTER :: dft_control
TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
TYPE(linres_control_type), POINTER :: linres_control
TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
@ -1019,7 +1023,8 @@ CONTAINS
virial=virial, &
sab_orb=sab_orb, &
energy=energy, &
rho_core=rho_core)
rho_core=rho_core, &
x_data=x_data)
p_env => qs_env%mp2_env%ri_grad%p_env
@ -1364,8 +1369,22 @@ CONTAINS
rho1 => p_env%p1
END IF
CALL derivatives_four_center(qs_env, rho_ao_kp, rho1, hfx_sections, para_env, &
1, use_virial)
IF (x_data(1, 1)%do_hfx_ri) THEN
IF (x_data(1, 1)%ri_data%flavor == ri_mo) THEN
CPABORT("NYI with RI_FLAVOR MO")
END IF
CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
x_data(1, 1)%general_parameter%fraction, &
rho_ao=rho_ao_kp, rho_ao_resp=rho1, &
use_virial=use_virial)
ELSE
CALL derivatives_four_center(qs_env, rho_ao_kp, rho1, hfx_sections, para_env, &
1, use_virial)
END IF
IF (use_virial) THEN
virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
virial%pv_virial = virial%pv_virial - virial%pv_fock_4c

24
src/qs_linres_kernel.F
View file

@ -35,12 +35,14 @@ MODULE qs_linres_kernel
dbcsr_set
USE hartree_local_methods, ONLY: Vh_1c_gg_integrals
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_none,&
do_admm_basis_projection,&
do_admm_exch_scaling_none,&
do_admm_purify_none,&
kg_tnadd_embed
USE input_cp2k_hfx, ONLY: ri_mo
USE input_section_types, ONLY: section_get_ival,&
section_get_lval,&
section_get_rval,&
@ -851,10 +853,24 @@ CONTAINS
rho_ao_kp(1:nspins, 1:1) => rho_ao(1:nspins)
DO irep = 1, n_rep_hf
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)
END DO
eh1 = 0.0_dp
IF (x_data(irep, 1)%do_hfx_ri) THEN
IF (x_data(irep, 1)%ri_data%flavor == ri_mo) THEN
CPABORT("NYI with RI_FLAVOR MO")
END IF
CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_kp, eh1, &
rho_ao=rho_ao_kp, geometry_did_change=s_mstruct_changed, &
nspins=nspins, hf_fraction=x_data(irep, 1)%general_parameter%fraction)
ELSE
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)
END DO
END IF
END DO
CALL timestop(handle)

2
src/qs_tddfpt2_fhxc_forces.F
View file

@ -591,6 +591,7 @@ CONTAINS
mhe(ispin, 1)%matrix => matrix_hfx_admm(ispin)%matrix
mpe(ispin, 1)%matrix => matrix_px1_admm(ispin)%matrix
END DO
IF (x_data(1, 1)%do_hfx_ri) CPABORT("RI-HFX with TDDFPT NYI")
DO ispin = 1, mspin
eh1 = 0.0
CALL integrate_four_center(qs_env, x_data, mhe, eh1, mpe, hfx_section, &
@ -648,6 +649,7 @@ CONTAINS
mhe(ispin, 1)%matrix => matrix_hfx(ispin)%matrix
mpe(ispin, 1)%matrix => matrix_px1(ispin)%matrix
END DO
IF (x_data(1, 1)%do_hfx_ri) CPABORT("RI-HFX with TDDFPT NYI")
DO ispin = 1, mspin
eh1 = 0.0
CALL integrate_four_center(qs_env, x_data, mhe, eh1, mpe, hfx_section, &

2
src/qs_tddfpt2_forces.F
View file

@ -751,6 +751,7 @@ CONTAINS
CALL dbcsr_set(mhz(ispin, 1)%matrix, 0.0_dp)
mpe(ispin, 1)%matrix => matrix_pe_admm(ispin)%matrix
END DO
IF (x_data(1, 1)%do_hfx_ri) CPABORT("RI-HFX with TDDFPT NYI")
DO ispin = 1, mspin
eh1 = 0.0
CALL integrate_four_center(qs_env, x_data, mhz, eh1, mpe, hfx_section, &
@ -786,6 +787,7 @@ CONTAINS
mhz(ispin, 1)%matrix => matrix_hz(ispin)%matrix
mpe(ispin, 1)%matrix => matrix_pe(ispin)%matrix
END DO
IF (x_data(1, 1)%do_hfx_ri) CPABORT("RI-HFX with TDDFPT NYI")
DO ispin = 1, mspin
eh1 = 0.0
CALL integrate_four_center(qs_env, x_data, mhz, eh1, mpe, hfx_section, &

1395
src/qs_tensors.F
View file

File diff suppressed because it is too large Load diff

102
src/response_solver.F
View file

@ -53,12 +53,15 @@ MODULE response_solver
USE exstates_types, ONLY: excited_energy_type
USE hfx_derivatives, ONLY: derivatives_four_center
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, ec_ls_solver, ec_mo_solver, kg_tnadd_atomic, kg_tnadd_embed, &
kg_tnadd_embed_ri, ls_s_sqrt_ns, ls_s_sqrt_proot, ot_precond_full_all, &
ot_precond_full_kinetic, ot_precond_full_single, ot_precond_full_single_inverse, &
ot_precond_none, ot_precond_s_inverse, precond_mlp
USE input_cp2k_hfx, ONLY: ri_mo
USE input_section_types, ONLY: section_get_lval,&
section_vals_get,&
section_vals_get_subs_vals,&
@ -1609,18 +1612,36 @@ CONTAINS
mpd(ispin, 1)%matrix => matrix_p(ispin, 1)%matrix
END DO
!
DO ispin = 1, mspin
eh1 = 0.0
CALL integrate_four_center(qs_env, x_data, mhz, eh1, mpz, hfx_section, &
para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
ispin=ispin)
END DO
DO ispin = 1, mspin
eh1 = 0.0
CALL integrate_four_center(qs_env, x_data, mhd, eh1, mpd, hfx_section, &
para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
ispin=ispin)
END DO
IF (x_data(1, 1)%do_hfx_ri) THEN
IF (x_data(1, 1)%ri_data%flavor == ri_mo) THEN
CPABORT("NYI with RI_FLAVOR MO")
END IF
eh1 = 0.0_dp
CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhz, eh1, rho_ao=mpz, &
geometry_did_change=s_mstruct_changed, nspins=nspins, &
hf_fraction=x_data(1, 1)%general_parameter%fraction)
eh1 = 0.0_dp
CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhd, eh1, rho_ao=mpd, &
geometry_did_change=s_mstruct_changed, nspins=nspins, &
hf_fraction=x_data(1, 1)%general_parameter%fraction)
ELSE
DO ispin = 1, mspin
eh1 = 0.0
CALL integrate_four_center(qs_env, x_data, mhz, eh1, mpz, hfx_section, &
para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
ispin=ispin)
END DO
DO ispin = 1, mspin
eh1 = 0.0
CALL integrate_four_center(qs_env, x_data, mhd, eh1, mpd, hfx_section, &
para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
ispin=ispin)
END DO
END IF
!
CALL get_qs_env(qs_env, admm_env=admm_env)
CPASSERT(ASSOCIATED(admm_env%work_aux_orb))
@ -1671,12 +1692,25 @@ CONTAINS
mhz(ispin, 1)%matrix => matrix_hz(ispin)%matrix
mpz(ispin, 1)%matrix => mpa(ispin)%matrix
END DO
DO ispin = 1, mspin
eh1 = 0.0
CALL integrate_four_center(qs_env, x_data, mhz, eh1, mpz, hfx_section, &
para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
ispin=ispin)
END DO
IF (x_data(1, 1)%do_hfx_ri) THEN
IF (x_data(1, 1)%ri_data%flavor == ri_mo) THEN
CPABORT("NYI with RI_FLAVOR MO")
END IF
eh1 = 0.0_dp
CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhz, eh1, rho_ao=mpz, &
geometry_did_change=s_mstruct_changed, nspins=nspins, &
hf_fraction=x_data(1, 1)%general_parameter%fraction)
ELSE
DO ispin = 1, mspin
eh1 = 0.0
CALL integrate_four_center(qs_env, x_data, mhz, eh1, mpz, hfx_section, &
para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
ispin=ispin)
END DO
END IF
DEALLOCATE (mhz, mpz)
END IF
@ -1698,13 +1732,37 @@ CONTAINS
CALL dbcsr_copy(matrix_pza(ispin)%matrix, matrix_pz_admm(ispin)%matrix)
END IF
END DO
CALL derivatives_four_center(qs_env, matrix_p, matrix_pza, hfx_section, para_env, &
1, use_virial, resp_only=resp_only)
IF (x_data(1, 1)%do_hfx_ri) THEN
IF (x_data(1, 1)%ri_data%flavor == ri_mo) THEN
CPABORT("NYI with RI_FLAVOR MO")
END IF
CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
x_data(1, 1)%general_parameter%fraction, &
rho_ao=matrix_p, rho_ao_resp=matrix_pza, &
use_virial=use_virial, resp_only=resp_only)
ELSE
CALL derivatives_four_center(qs_env, matrix_p, matrix_pza, hfx_section, para_env, &
1, use_virial, resp_only=resp_only)
END IF
CALL dbcsr_deallocate_matrix_set(matrix_pza)
ELSE
CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
CALL derivatives_four_center(qs_env, matrix_p, mpa, hfx_section, para_env, &
1, use_virial, resp_only=resp_only)
IF (x_data(1, 1)%do_hfx_ri) THEN
IF (x_data(1, 1)%ri_data%flavor == ri_mo) THEN
CPABORT("NYI with RI_FLAVOR MO")
END IF
CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
x_data(1, 1)%general_parameter%fraction, &
rho_ao=matrix_p, rho_ao_resp=mpa, &
use_virial=use_virial, resp_only=resp_only)
ELSE
CALL derivatives_four_center(qs_env, matrix_p, mpa, hfx_section, para_env, &
1, use_virial, resp_only=resp_only)
END IF
END IF
IF (debug_forces) THEN
fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3)

32
src/rpa_axk.F
View file

@ -36,9 +36,11 @@ MODULE rpa_axk
dbcsr_copy, dbcsr_create, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, dbcsr_release, &
dbcsr_set, dbcsr_trace, dbcsr_type, dbcsr_type_no_symmetry
USE hfx_energy_potential, ONLY: integrate_four_center
USE hfx_ri, ONLY: hfx_ri_update_ks
USE hfx_types, ONLY: hfx_create,&
hfx_release,&
hfx_type
USE input_cp2k_hfx, ONLY: ri_mo
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type
@ -52,6 +54,7 @@ MODULE rpa_axk
USE qs_subsys_types, ONLY: qs_subsys_get,&
qs_subsys_type
USE rpa_communication, ONLY: gamma_fm_to_dbcsr
USE scf_control_types, ONLY: scf_control_type
USE util, ONLY: get_limit
#include "./base/base_uses.f90"
@ -434,9 +437,20 @@ CONTAINS
rho_ao_2d(1:ns, 1:1) => rho_work_ao(1:ns)
CALL dbcsr_set(mat_2d(1, 1)%matrix, 0.0_dp)
CALL integrate_four_center(qs_env, x_data, mat_2d, ehfx, rho_ao_2d, hfx_sections, &
para_env_sub, my_recalc_hfx_integrals, irep, .TRUE., &
ispin=1)
IF (x_data(irep, 1)%do_hfx_ri) THEN
IF (x_data(irep, 1)%ri_data%flavor == ri_mo) THEN
CPABORT("NYI with RI_FLAVOR MO")
END IF
CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, mat_2d, ehfx, &
rho_ao=rho_ao_2d, geometry_did_change=my_recalc_hfx_integrals, &
nspins=ns, hf_fraction=x_data(irep, 1)%general_parameter%fraction)
ELSE
CALL integrate_four_center(qs_env, x_data, mat_2d, ehfx, rho_ao_2d, hfx_sections, &
para_env_sub, my_recalc_hfx_integrals, irep, .TRUE., &
ispin=1)
END IF
END DO
my_recalc_hfx_integrals = .FALSE.
@ -479,7 +493,7 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_create_subgroup'
INTEGER :: handle
INTEGER :: handle, nelectron_total
LOGICAL :: do_hfx
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: my_cell
@ -487,15 +501,18 @@ CONTAINS
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_subsys_type), POINTER :: subsys
TYPE(scf_control_type), POINTER :: scf_control
TYPE(section_vals_type), POINTER :: input
CALL timeset(routineN, handle)
NULLIFY (my_cell, atomic_kind_set, particle_set, dft_control, x_data, qs_kind_set)
NULLIFY (my_cell, atomic_kind_set, particle_set, dft_control, x_data, qs_kind_set, scf_control)
CALL get_qs_env(qs_env, &
subsys=subsys, &
input=input)
input=input, &
scf_control=scf_control, &
nelectron_total=nelectron_total)
CALL qs_subsys_get(subsys, &
cell=my_cell, &
@ -512,7 +529,8 @@ CONTAINS
IF (do_hfx) THEN
! Retrieve particle_set and atomic_kind_set
CALL hfx_create(x_data, para_env_sub, hfx_section, atomic_kind_set, &
qs_kind_set, particle_set, dft_control, my_cell, do_exx=.TRUE.)
qs_kind_set, particle_set, dft_control, my_cell, do_exx=.TRUE., &
do_ot=scf_control%use_ot, nelectron_total=nelectron_total)
END IF
CALL timestop(handle)

37
src/rpa_gw_sigma.F
View file

@ -40,11 +40,13 @@ MODULE rpa_gw_sigma
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_ri, ONLY: hfx_ri_update_ks
USE input_constants, ONLY: do_admm_basis_projection,&
do_admm_purify_none,&
gw_print_exx,&
gw_read_exx,&
xc_none
USE input_cp2k_hfx, ONLY: ri_mo
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
@ -125,7 +127,8 @@ CONTAINS
TYPE(admm_type), POINTER :: admm_env
TYPE(cp_fm_type), POINTER :: mo_coeff
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, rho_ao
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, &
matrix_ks_aux_fit_hfx, rho_ao
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_2d, matrix_ks_kp_im, &
matrix_ks_kp_re, matrix_ks_transl, matrix_sigma_x_minus_vxc, matrix_sigma_x_minus_vxc_im, &
rho_ao_2d
@ -141,7 +144,7 @@ CONTAINS
NULLIFY (admm_env, matrix_ks, matrix_ks_aux_fit, rho_ao, matrix_sigma_x_minus_vxc, input, &
xc_section, xc_section_admm_aux, xc_section_admm_prim, hfx_sections, rho, &
dft_control, para_env, ks_env, mo_coeff, matrix_sigma_x_minus_vxc_im)
dft_control, para_env, ks_env, mo_coeff, matrix_sigma_x_minus_vxc_im, matrix_ks_aux_fit_hfx)
CALL timeset(routineN, handle)
@ -167,7 +170,8 @@ CONTAINS
dft_control=dft_control, &
para_env=para_env, &
ks_env=ks_env, &
energy=energy)
energy=energy, &
matrix_ks_aux_fit_hfx=matrix_ks_aux_fit_hfx)
! RPA/GW with ADMM for EXX or the exchange self-energy only implemented for
! ADMM_PURIFICATION_METHOD NONE
@ -430,11 +434,28 @@ CONTAINS
matrix_ks_2d(1:ns, 1:1) => matrix_ks(1:ns)
END IF
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, .TRUE., irep, .TRUE., &
ispin=1)
ehfx = ehfx + eh1
IF (qs_env%mp2_env%ri_rpa%x_data(irep, 1)%do_hfx_ri) THEN
IF (qs_env%mp2_env%ri_rpa%x_data(irep, 1)%ri_data%flavor == ri_mo) THEN
CPABORT("NYI with RI_FLAVOR MO")
END IF
CALL hfx_ri_update_ks(qs_env, qs_env%mp2_env%ri_rpa%x_data(irep, 1)%ri_data, matrix_ks_2d, ehfx, &
rho_ao=rho_ao_2d, geometry_did_change=.TRUE., nspins=nspins, &
hf_fraction=qs_env%mp2_env%ri_rpa%x_data(irep, 1)%general_parameter%fraction)
IF (do_admm_rpa) THEN
!for ADMMS, we need the exchange matrix k(d) for both spins
DO ispin = 1, nspins
CALL dbcsr_copy(matrix_ks_aux_fit_hfx(ispin)%matrix, matrix_ks_2d(ispin, 1)%matrix, &
name="HF exch. part of matrix_ks_aux_fit for ADMMS")
END DO
END IF
ELSE
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, .TRUE., irep, .TRUE., &
ispin=1)
ehfx = ehfx + eh1
END IF
END DO
END IF
energy_ex = ehfx

23
src/rpa_hfx.F
View file

@ -16,6 +16,8 @@ MODULE rpa_hfx
USE dbcsr_api, ONLY: dbcsr_p_type,&
dbcsr_set
USE hfx_energy_potential, ONLY: integrate_four_center
USE hfx_ri, ONLY: hfx_ri_update_ks
USE input_cp2k_hfx, ONLY: ri_mo
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type
@ -133,11 +135,22 @@ CONTAINS
ELSE
matrix_ks_2d(1:ns, 1:1) => matrix_ks(1:ns)
END IF
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, .TRUE., irep, .TRUE., &
ispin=1)
ehfx = ehfx + eh1
IF (qs_env%mp2_env%ri_rpa%x_data(irep, 1)%do_hfx_ri) THEN
IF (qs_env%mp2_env%ri_rpa%x_data(irep, 1)%ri_data%flavor == ri_mo) THEN
CPABORT("NYI with RI_FLAVOR MO")
END IF
CALL hfx_ri_update_ks(qs_env, qs_env%mp2_env%ri_rpa%x_data(irep, 1)%ri_data, matrix_ks_2d, ehfx, &
rho_ao=rho_ao_2d, geometry_did_change=.TRUE., nspins=ns, &
hf_fraction=qs_env%mp2_env%ri_rpa%x_data(irep, 1)%general_parameter%fraction)
ELSE
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, .TRUE., irep, .TRUE., &
ispin=1)
ehfx = ehfx + eh1
END IF
END DO
! include the EXX contribution to the total energy

38
src/rpa_rse.F
View file

@ -39,9 +39,12 @@ MODULE rpa_rse
dbcsr_init_p,&
dbcsr_p_type,&
dbcsr_release,&
dbcsr_scale,&
dbcsr_set,&
dbcsr_type_symmetric
USE hfx_energy_potential, ONLY: integrate_four_center
USE hfx_ri, ONLY: hfx_ri_update_ks
USE input_cp2k_hfx, ONLY: ri_mo
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
@ -353,15 +356,34 @@ CONTAINS
CALL dbcsr_set(mat_mu_nu(1)%matrix, 0.0_dp)
DO irep = 1, n_rep_hf
rho_ao_2d(1:ns, 1:1) => P_mu_nu(1:ns)
mat_2d(1:ns, 1:1) => mat_mu_nu(1:ns)
CALL integrate_four_center(qs_env, qs_env%mp2_env%ri_rpa%x_data, mat_2d, ehfx, rho_ao_2d, hfx_sections, &
para_env, my_recalc_hfx_integrals, irep, .TRUE., &
ispin=1)
IF (qs_env%mp2_env%ri_rpa%x_data(1, 1)%do_hfx_ri) THEN
IF (qs_env%mp2_env%ri_rpa%x_data(1, 1)%ri_data%flavor == ri_mo) THEN
CPABORT("NYI with RI_FLAVOR MO")
END IF
my_recalc_hfx_integrals = .FALSE.
END DO
DO irep = 1, n_rep_hf
rho_ao_2d(1:ns, 1:1) => P_mu_nu(1:ns)
mat_2d(1:ns, 1:1) => mat_mu_nu(1:ns)
CALL hfx_ri_update_ks(qs_env, qs_env%mp2_env%ri_rpa%x_data(irep, 1)%ri_data, mat_2d, ehfx, &
rho_ao=rho_ao_2d, geometry_did_change=my_recalc_hfx_integrals, nspins=1, &
hf_fraction=qs_env%mp2_env%ri_rpa%x_data(irep, 1)%general_parameter%fraction)
IF (ns == 2) CALL dbcsr_scale(mat_mu_nu(1)%matrix, 2.0_dp)
my_recalc_hfx_integrals = .FALSE.
END DO
ELSE
DO irep = 1, n_rep_hf
rho_ao_2d(1:ns, 1:1) => P_mu_nu(1:ns)
mat_2d(1:ns, 1:1) => mat_mu_nu(1:ns)
CALL integrate_four_center(qs_env, qs_env%mp2_env%ri_rpa%x_data, mat_2d, ehfx, rho_ao_2d, hfx_sections, &
para_env, my_recalc_hfx_integrals, irep, .TRUE., &
ispin=1)
my_recalc_hfx_integrals = .FALSE.
END DO
END IF
! copy back to fm
CALL cp_fm_set_all(fm_X_ao, 0.0_dp)

8
src/xc_adiabatic_utils.F
View file

@ -20,6 +20,7 @@ MODULE xc_adiabatic_utils
USE dbcsr_api, ONLY: dbcsr_p_type
USE hfx_communication, ONLY: scale_and_add_fock_to_ks_matrix
USE hfx_derivatives, ONLY: derivatives_four_center
USE hfx_types, ONLY: hfx_type
USE input_constants, ONLY: do_adiabatic_hybrid_mcy3,&
do_adiabatic_model_pade
USE input_section_types, ONLY: section_vals_get,&
@ -88,6 +89,7 @@ CONTAINS
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_resp
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
TYPE(dft_control_type), POINTER :: dft_control
TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho_xc
TYPE(section_vals_type), POINTER :: adiabatic_rescaling_section, &
@ -95,15 +97,17 @@ CONTAINS
CALL timeset(routineN, handle)
NULLIFY (para_env, dft_control, adiabatic_rescaling_section, hfx_sections, &
input, xc_section, rho_xc, ks_env, rho_ao, rho_ao_resp)
input, xc_section, rho_xc, ks_env, rho_ao, rho_ao_resp, x_data)
CALL get_qs_env(qs_env, &
dft_control=dft_control, &
para_env=para_env, &
input=input, &
rho_xc=rho_xc, &
ks_env=ks_env)
ks_env=ks_env, &
x_data=x_data)
IF (x_data(1, 1)%do_hfx_ri) CPABORT("RI-HFX not compatible with this kinf of functionals")
nimages = dft_control%nimages
CPASSERT(nimages == 1)

61
tests/QS/regtest-hfx-ri-2/CH-hfx-ri-mo.inp Normal file
View file

@ -0,0 +1,61 @@
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME EMSL_BASIS_SETS
POTENTIAL_FILE_NAME POTENTIAL
AUTO_BASIS RI_HFX SMALL
LSD
&MGRID
CUTOFF 300
REL_CUTOFF 50
&END MGRID
&QS
METHOD GAPW
&END QS
&SCF
SCF_GUESS ATOMIC
MAX_SCF 20
EPS_SCF 1.0E-07
&OT
PRECONDITIONER FULL_ALL
&END
&END SCF
&XC
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&HF
&RI
RI_FLAVOR MO
&END
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC 5.0 5.0 5.0
PERIODIC NONE
&END CELL
&COORD
C 0.000000 0.000000 0.2581
H 0.000000 0.000000 -0.9487
&END COORD
&KIND C
BASIS_SET 6-31Gx
POTENTIAL ALL
&END KIND
&KIND H
BASIS_SET 6-31Gx
POTENTIAL ALL
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT CH-hfx-ri-mo
RUN_TYPE GEO_OPT
PRINT_LEVEL MEDIUM
&END GLOBAL
&MOTION
&GEO_OPT
MAX_ITER 1
&END
&END MOTION

58
tests/QS/regtest-hfx-ri-2/CH-hfx-ri-rho.inp Normal file
View file

@ -0,0 +1,58 @@
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME EMSL_BASIS_SETS
POTENTIAL_FILE_NAME POTENTIAL
AUTO_BASIS RI_HFX SMALL
LSD
&MGRID
CUTOFF 300
REL_CUTOFF 50
&END MGRID
&QS
METHOD GAPW
&END QS
&SCF
SCF_GUESS ATOMIC
MAX_SCF 20
EPS_SCF 1.0E-08
&END SCF
&XC
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&HF
&RI
RI_FLAVOR RHO
&END
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC 5.0 5.0 5.0
PERIODIC NONE
&END CELL
&COORD
C 0.000000 0.000000 0.2581
H 0.000000 0.000000 -0.9487
&END COORD
&KIND C
BASIS_SET 6-31Gx
POTENTIAL ALL
&END KIND
&KIND H
BASIS_SET 6-31Gx
POTENTIAL ALL
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT CH-hfx-ri-rho
RUN_TYPE GEO_OPT
PRINT_LEVEL MEDIUM
&END GLOBAL
&MOTION
&GEO_OPT
MAX_ITER 1
&END
&END MOTION

75
tests/QS/regtest-hfx-ri-2/CH3-b3lyp-ADMM.inp Normal file
View file

@ -0,0 +1,75 @@
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS_MOLOPT
BASIS_SET_FILE_NAME BASIS_ADMM
POTENTIAL_FILE_NAME POTENTIAL_UZH
AUTO_BASIS RI_HFX SMALL
LSD
&MGRID
CUTOFF 200
REL_CUTOFF 30
&END MGRID
&QS
METHOD GPW
&END QS
&AUXILIARY_DENSITY_MATRIX_METHOD
&END
&POISSON
PERIODIC NONE
PSOLVER MT
&END
&SCF
EPS_SCF 1.0E-6
SCF_GUESS ATOMIC
MAX_SCF 5
&OT
PRECONDITIONER FULL_ALL
&END
&END SCF
&XC
&XC_FUNCTIONAL
&LIBXC
FUNCTIONAL HYB_GGA_XC_B3LYP
&END
&END XC_FUNCTIONAL
&HF
FRACTION 0.2
&RI
&END
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC 8.0 8.0 8.0
PERIODIC NONE
&END CELL
&COORD
C 0.0000 0.0000 0.5000
H 0.0000 1.0728 0.0000
H 0.9291 -0.5364 0.0000
H -0.9291 -0.5364 0.0000
&END COORD
&KIND H
BASIS_SET DZVP-MOLOPT-GTH
BASIS_SET AUX_FIT FIT3
POTENTIAL GTH-HYB-q1
&END KIND
&KIND C
BASIS_SET DZVP-MOLOPT-GTH
BASIS_SET AUX_FIT FIT3
POTENTIAL GTH-HYB-q4
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT CH3-b3lyp-ADMM
PRINT_LEVEL MEDIUM
RUN_TYPE GEO_OPT
&END GLOBAL
&MOTION
&GEO_OPT
MAX_ITER 1
&END GEO_OPT
&END MOTION

59
tests/QS/regtest-hfx-ri-2/H2O-hfx-stress-identity.inp Normal file
View file

@ -0,0 +1,59 @@
&FORCE_EVAL
METHOD Quickstep
STRESS_TENSOR ANALYTICAL
&PRINT
&STRESS_TENSOR
&END
&END
&DFT
BASIS_SET_FILE_NAME EMSL_BASIS_SETS
POTENTIAL_FILE_NAME POTENTIAL
AUTO_BASIS RI_HFX SMALL
&MGRID
CUTOFF 300
REL_CUTOFF 50
&END MGRID
&QS
METHOD GAPW
&END QS
&SCF
SCF_GUESS ATOMIC
MAX_SCF 20
EPS_SCF 1.0E-07
&END SCF
&XC
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&HF
&RI
&END
&INTERACTION_POTENTIAL
POTENTIAL_TYPE IDENTITY
&END
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC 6.0 6.0 6.0
&END CELL
&COORD
O 0.000000 0.000000 -0.065587
H 0.000000 -0.757136 0.520545
H 0.000000 0.757136 0.520545
&END COORD
&KIND O
BASIS_SET Ahlrichs-def2-SVP
POTENTIAL ALL
&END KIND
&KIND H
BASIS_SET Ahlrichs-def2-SVP
POTENTIAL ALL
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT H2O-hfx-stress-identity
RUN_TYPE ENERGY_FORCE
PRINT_LEVEL MEDIUM
&END GLOBAL

65
tests/QS/regtest-hfx-ri-2/H2O-pbe0-stress-truncated.inp Normal file
View file

@ -0,0 +1,65 @@
&FORCE_EVAL
METHOD Quickstep
STRESS_TENSOR ANALYTICAL
&PRINT
&STRESS_TENSOR
&END
&END
&DFT
BASIS_SET_FILE_NAME BASIS_MOLOPT
POTENTIAL_FILE_NAME POTENTIAL_UZH
AUTO_BASIS RI_HFX SMALL
&MGRID
CUTOFF 200
REL_CUTOFF 30
&END MGRID
&QS
METHOD GPW
&END QS
&SCF
SCF_GUESS ATOMIC
MAX_SCF 20
EPS_SCF 1.0E-07
&END SCF
&XC
&XC_FUNCTIONAL PBE
&PBE
SCALE_C 1.0
SCALE_X 0.75
&END
&END XC_FUNCTIONAL
&HF
FRACTION 0.25
&RI
&END
&INTERACTION_POTENTIAL
POTENTIAL_TYPE TRUNCATED
CUTOFF_RADIUS 2.0
&END
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC 6.0 6.0 6.0
&END CELL
&COORD
O 0.000000 0.000000 -0.065587
H 0.000000 -0.757136 0.520545
H 0.000000 0.757136 0.520545
&END COORD
&KIND O
BASIS_SET SZV-MOLOPT-GTH
POTENTIAL GTH-PBE0-q6
&END KIND
&KIND H
BASIS_SET SZV-MOLOPT-GTH
POTENTIAL GTH-PBE0-q1
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT H2O-pbe0-stress-truncated
RUN_TYPE ENERGY_FORCE
PRINT_LEVEL MEDIUM
&END GLOBAL

63
tests/QS/regtest-hfx-ri-2/Ne-hfx-pbc-metric-mo.inp Normal file
View file

@ -0,0 +1,63 @@
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME EMSL_BASIS_SETS
POTENTIAL_FILE_NAME POTENTIAL
AUTO_BASIS RI_HFX SMALL
&MGRID
CUTOFF 300
REL_CUTOFF 50
&END MGRID
&QS
METHOD GAPW
&END QS
&SCF
EPS_SCF 1.0E-7
SCF_GUESS ATOMIC
MAX_SCF 5
&OT ON
PRECONDITIONER FULL_ALL
&END
&END SCF
&XC
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&HF
&RI
RI_FLAVOR MO
RI_METRIC IDENTITY
&END
&INTERACTION_POTENTIAL
POTENTIAL_TYPE SHORTRANGE
OMEGA 0.11
&END
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC 5.0 5.0 10.0
&END CELL
&COORD
Ne 0.000000 0.000000 0.000000
Ne 0.000000 0.000000 2.800000
Ne 0.000000 0.000000 4.000000
Ne 0.000000 0.000000 6.100000
Ne 0.000000 0.000000 8.900000
&END COORD
&KIND Ne
BASIS_SET 3-21Gx
POTENTIAL ALL
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT Ne-hfx-pbc-metric-mo
PRINT_LEVEL MEDIUM
RUN_TYPE MD
&END GLOBAL
&MOTION
&MD
MAX_STEPS 1
&END
&END MOTION

63
tests/QS/regtest-hfx-ri-2/Ne-hfx-pbc-metric-rho.inp Normal file
View file

@ -0,0 +1,63 @@
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME EMSL_BASIS_SETS
POTENTIAL_FILE_NAME POTENTIAL
AUTO_BASIS RI_HFX SMALL
&MGRID
CUTOFF 300
REL_CUTOFF 50
&END MGRID
&QS
METHOD GAPW
&END QS
&SCF
EPS_SCF 1.0E-7
SCF_GUESS ATOMIC
MAX_SCF 5
&OT ON
PRECONDITIONER FULL_ALL
&END
&END SCF
&XC
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&HF
&RI
RI_FLAVOR RHO
RI_METRIC IDENTITY
&END
&INTERACTION_POTENTIAL
POTENTIAL_TYPE TRUNCATED
CUTOFF_RADIUS 2.0
&END
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC 5.0 5.0 10.0
&END CELL
&COORD
Ne 0.000000 0.000000 0.000000
Ne 0.000000 0.000000 2.800000
Ne 0.000000 0.000000 4.000000
Ne 0.000000 0.000000 6.100000
Ne 0.000000 0.000000 8.900000
&END COORD
&KIND Ne
BASIS_SET 3-21Gx
POTENTIAL ALL
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT Ne-hfx-pbc-metric-rho
PRINT_LEVEL MEDIUM
RUN_TYPE MD
&END GLOBAL
&MOTION
&MD
MAX_STEPS 1
&END
&END MOTION

9
tests/QS/regtest-hfx-ri-2/TEST_FILES Normal file
View file

@ -0,0 +1,9 @@
#Testing the forces:
CH-hfx-ri-rho.inp 11 1.0E-8 -38.259646458859827
CH-hfx-ri-mo.inp 11 1.0E-9 -38.262303893579102
Ne-hfx-pbc-metric-rho.inp 11 1.0E-9 -633.525916555354343
Ne-hfx-pbc-metric-mo.inp 11 1.0E-9 -632.285120513072002
CH3-b3lyp-ADMM.inp 11 1.0E-9 -7.413742783916254
H2O-pbe0-stress-truncated.inp 31 1.0E-9 2.01366131588E-01
H2O-hfx-stress-identity.inp 31 1.0E-9 6.44733085054E-02
#EOF

2
tests/QS/regtest-hfx-ri/CH3-ADMM.inp
View file

@ -78,7 +78,7 @@
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT CH3-BP-MO_DIAG
PROJECT CH3-ADMM
PRINT_LEVEL MEDIUM
RUN_TYPE ENERGY
&TIMINGS

2
tests/QS/regtest-hfx-ri/CH3-hfx-converged.inp
View file

@ -61,7 +61,7 @@
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT CH3-TZV2P-converged
PROJECT CH3-hfx-converged
PRINT_LEVEL MEDIUM
RUN_TYPE ENERGY
&END GLOBAL

2
tests/QS/regtest-hfx-ri/H2O-hfx-periodic-ri-truncated.inp
View file

@ -65,7 +65,7 @@
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT H2O-ri-trunc
PROJECT H2O-hfx-periodic-ri-truncated
PRINT_LEVEL MEDIUM
RUN_TYPE ENERGY
&END GLOBAL

2
tests/QS/regtest-hfx-ri/Ne-hybrid-periodic-shortrange.inp
View file

@ -69,7 +69,7 @@
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT NE-hybrid-HSE06-lda
PROJECT Ne-hybrid-periodic-shortrange
PRINT_LEVEL MEDIUM
RUN_TYPE ENERGY
&END GLOBAL

103
tests/QS/regtest-mp2-grad/H2O_grad_ri-hfx.inp Normal file
View file

@ -0,0 +1,103 @@
&GLOBAL
PROJECT GRAD_H2O_gpw
PRINT_LEVEL LOW
RUN_TYPE GEO_OPT
&TIMINGS
THRESHOLD 0.001
&END
&END GLOBAL
&MOTION
&GEO_OPT
MAX_ITER 1
&END
&END MOTION
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME GTH_BASIS_SETS
BASIS_SET_FILE_NAME HFX_BASIS
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 150
REL_CUTOFF 30
&END MGRID
&QS
METHOD GPW
EPS_DEFAULT 1.0E-12
&END QS
&SCF
SCF_GUESS ATOMIC
EPS_SCF 1.0E-6
MAX_SCF 100
&END SCF
&XC
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&HF
FRACTION 1.0000000
&SCREENING
SCREEN_ON_INITIAL_P .FALSE.
EPS_SCHWARZ 1.0E-9
EPS_SCHWARZ_FORCES 1.0E-9
&END SCREENING
&INTERACTION_POTENTIAL
POTENTIAL_TYPE TRUNCATED
CUTOFF_RADIUS 2.0
T_C_G_DATA t_c_g.dat
&END
&RI
&END RI
&END HF
&WF_CORRELATION
&RI_MP2
BLOCK_SIZE 1
EPS_CANONICAL 0.0001
FREE_HFX_BUFFER .TRUE.
&CPHF
EPS_CONV 1.0E-4
MAX_ITER 10
&END
&END
&INTEGRALS
&WFC_GPW
CUTOFF 100
REL_CUTOFF 30
EPS_FILTER 1.0E-6
EPS_GRID 1.0E-6
&END WFC_GPW
ERI_METHOD MME
&END INTEGRALS
MEMORY 1.00
NUMBER_PROC 1
&END
&END XC
&END DFT
&PRINT
&FORCES
&END
&END
&SUBSYS
&CELL
ABC [angstrom] 5.0 5.0 5.0
&END CELL
&KIND H
BASIS_SET SZV-GTH
BASIS_SET RI_AUX RI_DZVP-GTH
POTENTIAL GTH-HF-q1
&END KIND
&KIND O
BASIS_SET SZV-GTH
BASIS_SET RI_AUX RI_DZVP-GTH
POTENTIAL GTH-HF-q6
&END KIND
&COORD
O 0.000000 0.000000 -0.211000
H 0.000000 -0.844000 0.495000
H 0.000000 0.744000 0.495000
&END
&TOPOLOGY
&CENTER_COORDINATES
&END
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL

94
tests/QS/regtest-mp2-grad/O2_dyn_ri-hfx.inp Normal file
View file

@ -0,0 +1,94 @@
&GLOBAL
PROJECT O2_dyn_ri-hfx
PRINT_LEVEL LOW
&TIMINGS
THRESHOLD 0.01
&END
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVE
STEPS 1
&END
&END MOTION
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME HFX_BASIS
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 100
REL_CUTOFF 20
&END MGRID
&QS
METHOD GPW
EPS_DEFAULT 1.0E-10
EPS_PGF_ORB 1.0E-20
&END QS
&SCF
SCF_GUESS ATOMIC
EPS_SCF 1.0E-5
MAX_SCF 100
&END SCF
&XC
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&HF
FRACTION 1.0000000
&SCREENING
EPS_SCHWARZ 1.0E-6
SCREEN_ON_INITIAL_P FALSE
&END SCREENING
&INTERACTION_POTENTIAL
POTENTIAL_TYPE TRUNCATED
CUTOFF_RADIUS 1.5
T_C_G_DATA t_c_g.dat
&END
&RI
&END
&END HF
&WF_CORRELATION
&RI_MP2
BLOCK_SIZE 1
EPS_CANONICAL 0.0001
FREE_HFX_BUFFER .TRUE.
&END RI_MP2
&INTEGRALS
&WFC_GPW
CUTOFF 50
REL_CUTOFF 25
EPS_FILTER 1.0E-5
EPS_GRID 1.0E-4
&END WFC_GPW
&END INTEGRALS
MEMORY 500.0
NUMBER_PROC 1
&END
&END XC
UKS
MULTIPLICITY 3
&END DFT
&SUBSYS
&VELOCITY
0.0 0.0 0.0
0.0 0.0 0.0
&END VELOCITY
&CELL
ABC [angstrom] 6.000 6.000 6.000
!PERIODIC NONE
&END CELL
&KIND O
BASIS_SET DZVP-GTH
BASIS_SET RI_AUX RI_DZVP-GTH
POTENTIAL GTH-HF-q6
&END KIND
&COORD
O 4.0000000084 4.0000000084 4.6623718822
O 3.9999999905 3.9999999905 3.3376281178
&END
&TOPOLOGY
&CENTER_COORDINATES
&END
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL

2
tests/QS/regtest-mp2-grad/TEST_FILES
View file

@ -7,4 +7,6 @@ CH_dyn_screen.inp 11 2e-04
MOM_MP2_geoopt.inp 11 2e-06 -13.969926262647746
H2O_MD_mme.inp 11 1e-10 -17.056807598425291
H2_MP2_debug.inp 11 1e-10 -1.146241031776471
H2O_grad_ri-hfx.inp 11 1e-10 -16.763446635143474
O2_dyn_ri-hfx.inp 11 1e-10 -31.518462685641584
#EOF

85
tests/QS/regtest-rpa-cubic-scaling/Cubic_RPA_CH3_ri-hfx.inp Normal file
View file

@ -0,0 +1,85 @@
&GLOBAL
PROJECT RI_RPA_H2O_minimax
PRINT_LEVEL MEDIUM
RUN_TYPE ENERGY
&TIMINGS
THRESHOLD 0.01
&END
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME HFX_BASIS
SORT_BASIS EXP
POTENTIAL_FILE_NAME GTH_POTENTIALS
UKS
MULTIPLICITY 2
&MGRID
CUTOFF 100
REL_CUTOFF 20
&END MGRID
&POISSON
PERIODIC NONE
POISSON_SOLVER MT
&END POISSON
&QS
METHOD GPW
EPS_DEFAULT 1.0E-15
EPS_PGF_ORB 1.0E-30
&END QS
&SCF
SCF_GUESS ATOMIC
EPS_SCF 1.0E-7
MAX_SCF 100
&PRINT
&RESTART OFF
&END
&END
&END SCF
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&WF_CORRELATION
&LOW_SCALING
MEMORY_CUT 1
&END
&RI_RPA
RPA_NUM_QUAD_POINTS 6
&HF
FRACTION 1.0000000
&SCREENING
EPS_SCHWARZ 1.0E-8
SCREEN_ON_INITIAL_P FALSE
&END SCREENING
&RI
&END
&END HF
&END RI_RPA
MEMORY 200.
NUMBER_PROC 1
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC [angstrom] 9.000 9.000 9.000
PERIODIC NONE
&END CELL
&KIND H
BASIS_SET DZVP-GTH
BASIS_SET RI_AUX RI_DZVP-GTH
POTENTIAL GTH-PBE-q1
&END KIND
&KIND C
BASIS_SET DZVP-GTH
BASIS_SET RI_AUX RI_DZVP-GTH
POTENTIAL GTH-PBE-q4
&END KIND
&TOPOLOGY
COORD_FILE_NAME CH3.xyz
COORD_FILE_FORMAT xyz
&CENTER_COORDINATES
&END
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL

82
tests/QS/regtest-rpa-cubic-scaling/Cubic_RPA_H2O_ri-hfx.inp Normal file
View file

@ -0,0 +1,82 @@
&GLOBAL
PROJECT RI_RPA_H2O_minimax
PRINT_LEVEL MEDIUM
RUN_TYPE ENERGY
&TIMINGS
THRESHOLD 0.01
&END
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME HFX_BASIS
SORT_BASIS EXP
POTENTIAL_FILE_NAME GTH_POTENTIALS
&MGRID
CUTOFF 100
REL_CUTOFF 20
&END MGRID
&POISSON
PERIODIC NONE
POISSON_SOLVER WAVELET
&END POISSON
&QS
METHOD GPW
EPS_DEFAULT 1.0E-15
EPS_PGF_ORB 1.0E-30
&END QS
&SCF
SCF_GUESS ATOMIC
EPS_SCF 1.0E-7
MAX_SCF 100
&PRINT
&RESTART OFF
&END
&END
&END SCF
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&WF_CORRELATION
&LOW_SCALING
&END
&RI_RPA
RPA_NUM_QUAD_POINTS 6
&HF
FRACTION 1.0000000
&SCREENING
EPS_SCHWARZ 1.0E-8
SCREEN_ON_INITIAL_P FALSE
&END SCREENING
&RI
&END RI
&END HF
&END RI_RPA
MEMORY 200.
NUMBER_PROC 1
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC [angstrom] 8.000 8.000 8.000
PERIODIC NONE
&END CELL
&KIND H
BASIS_SET DZVP-GTH
BASIS_SET RI_AUX RI_DZVP-GTH
POTENTIAL GTH-PBE-q1
&END KIND
&KIND O
BASIS_SET DZVP-GTH
BASIS_SET RI_AUX RI_DZVP-GTH
POTENTIAL GTH-PBE-q6
&END KIND
&TOPOLOGY
COORD_FILE_NAME H2O_gas.xyz
COORD_FILE_FORMAT xyz
&CENTER_COORDINATES
&END
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL

2
tests/QS/regtest-rpa-cubic-scaling/TEST_FILES
View file

@ -6,4 +6,6 @@ Cubic_RPA_H2O_standard_svd.inp 11 1e-08 -
RPA_kpoints_H2O.inp 11 1e-08 -17.384612021187685
RPA_kpoints_H2O_batched.inp 11 1e-08 -17.384612021187685
RPA_kpoints_from_Gamma_H2O.inp 11 1e-08 -17.384615079520682
Cubic_RPA_CH3_ri-hfx.inp 11 1e-08 -7.414799528509679
Cubic_RPA_H2O_ri-hfx.inp 11 1e-08 -17.159865017693086
#EOF

1
tests/TEST_DIRS
View file

@ -15,6 +15,7 @@ QS/regtest-mp2-block libint parallel mpir
QS/regtest-mp2-admm-grad libint
QS/regtest-double-hybrid-stress libint
QS/regtest-double-hybrid-grad libint
QS/regtest-hfx-ri-2 libint libxc
QS/regtest-rma libint mpi3 mpiranks==4
LIBTEST/libvori libvori
LIBTEST/libbqb libbqb