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

Refactoring of KG code:

Browse source 
- TNADD correction separated from XC functionals, Input not compatible!
- Enforce minimum image convention for molecular neighbor lists
  This changes reference energies for most regtests!
- Add basic code for a new method (GPW-KG + LRI TNADD)


svn-origin-rev: 18233
This commit is contained in:
Jürg Hutter 2018-01-25 14:18:01 +00:00
parent 368684fa46
commit 28fde2b103
42 changed files with 942 additions and 517 deletions
Download patch file Download diff file Expand all files Collapse all files

1
src/input_constants.F
View file

@ -974,6 +974,7 @@ MODULE input_constants
kg_color_greedy = 1
! kg kinetic energy embedding potential
INTEGER, PARAMETER, PUBLIC :: kg_tnadd_embed = 100, &
kg_tnadd_embed_ri = 101, &
kg_tnadd_atomic = 200, &
kg_tnadd_none = 300

79
src/input_cp2k_dft.F
View file

@ -59,38 +59,38 @@ MODULE input_cp2k_dft
general_roks, hf_model, high_spin_roks, history_guess, jacobian_fd1, &
jacobian_fd1_backward, jacobian_fd1_central, jacobian_fd2, jacobian_fd2_backward, &
kg_cholesky, kg_color_dsatur, kg_color_greedy, kg_ec_diagonalization, &
kg_ec_functional_harris, kg_tnadd_atomic, kg_tnadd_embed, kg_tnadd_none, ls_2pnt, ls_3pnt, &
ls_gold, ls_none, mao_basis_ext, mao_basis_orb, mao_basis_prim, mao_projection, &
mopac_guess, no_excitations, no_guess, numerical, oe_gllb, oe_lb, oe_none, oe_saop, &
oe_sic, op_loc_berry, op_loc_boys, op_loc_pipek, orb_dx2, orb_dxy, orb_dy2, orb_dyz, &
orb_dz2, orb_dzx, orb_px, orb_py, orb_pz, orb_s, ot_algo_irac, ot_algo_taylor_or_diag, &
ot_chol_irac, ot_lwdn_irac, ot_mini_broyden, ot_mini_cg, ot_mini_diis, ot_mini_sd, &
ot_poly_irac, ot_precond_full_all, ot_precond_full_kinetic, ot_precond_full_single, &
ot_precond_full_single_inverse, ot_precond_none, ot_precond_s_inverse, &
ot_precond_solver_default, ot_precond_solver_direct, ot_precond_solver_inv_chol, &
ot_precond_solver_update, outer_scf_basis_center_opt, outer_scf_becke_constraint, &
outer_scf_cdft_constraint, outer_scf_ddapc_constraint, outer_scf_hirshfeld_constraint, &
outer_scf_none, outer_scf_optimizer_bisect, outer_scf_optimizer_broyden, &
outer_scf_optimizer_diis, outer_scf_optimizer_newton, outer_scf_optimizer_newton_ls, &
outer_scf_optimizer_none, outer_scf_optimizer_sd, outer_scf_optimizer_secant, &
outer_scf_s2_constraint, plus_u_lowdin, plus_u_mulliken, plus_u_mulliken_charges, &
pw_interp, radius_covalent, radius_default, radius_single, radius_user, radius_vdw, &
random_guess, real_time_propagation, ref_charge_atomic, ref_charge_mulliken, rel_dkh, &
rel_none, rel_pot_erfc, rel_pot_full, rel_sczora_mp, rel_trans_atom, rel_trans_full, &
rel_trans_molecule, rel_zora, rel_zora_full, rel_zora_mp, restart_guess, rsdft_model, &
sccs_andreussi, sccs_derivative_cd3, sccs_derivative_cd5, sccs_derivative_cd7, &
sccs_derivative_fft, sccs_fattebert_gygi, shape_function_density, shape_function_gaussian, &
sic_ad, sic_eo, sic_list_all, sic_list_unpaired, sic_mauri_spz, sic_mauri_us, sic_none, &
slater, smear_energy_window, smear_fermi_dirac, smear_list, sparse_guess, &
spline3_nopbc_interp, spline3_pbc_interp, tddfpt_davidson, tddfpt_excitations, &
tddfpt_lanczos, tddfpt_singlet, tddfpt_triplet, use_coulomb, use_diff, use_no, &
use_restart_wfn, use_rt_restart, use_scf_wfn, wannier_projection, weight_type_mass, &
weight_type_unit, wfi_aspc_nr, wfi_frozen_method_nr, wfi_linear_p_method_nr, &
wfi_linear_ps_method_nr, wfi_linear_wf_method_nr, wfi_ps_method_nr, &
wfi_use_guess_method_nr, wfi_use_prev_p_method_nr, wfi_use_prev_rho_r_method_nr, &
wfi_use_prev_wf_method_nr, xas_1s_type, xas_2p_type, xas_2s_type, xas_dip_len, &
xas_dip_vel, xas_dscf, xas_none, xas_tp_fh, xas_tp_flex, xas_tp_hh, xas_tp_xfh, &
xas_tp_xhh, xes_tp_val
kg_ec_functional_harris, kg_tnadd_atomic, kg_tnadd_embed, kg_tnadd_embed_ri, &
kg_tnadd_none, ls_2pnt, ls_3pnt, ls_gold, ls_none, mao_basis_ext, mao_basis_orb, &
mao_basis_prim, mao_projection, mopac_guess, no_excitations, no_guess, numerical, oe_gllb, &
oe_lb, oe_none, oe_saop, oe_sic, op_loc_berry, op_loc_boys, op_loc_pipek, orb_dx2, &
orb_dxy, orb_dy2, orb_dyz, orb_dz2, orb_dzx, orb_px, orb_py, orb_pz, orb_s, ot_algo_irac, &
ot_algo_taylor_or_diag, ot_chol_irac, ot_lwdn_irac, ot_mini_broyden, ot_mini_cg, &
ot_mini_diis, ot_mini_sd, ot_poly_irac, ot_precond_full_all, ot_precond_full_kinetic, &
ot_precond_full_single, ot_precond_full_single_inverse, ot_precond_none, &
ot_precond_s_inverse, ot_precond_solver_default, ot_precond_solver_direct, &
ot_precond_solver_inv_chol, ot_precond_solver_update, outer_scf_basis_center_opt, &
outer_scf_becke_constraint, outer_scf_cdft_constraint, outer_scf_ddapc_constraint, &
outer_scf_hirshfeld_constraint, outer_scf_none, outer_scf_optimizer_bisect, &
outer_scf_optimizer_broyden, outer_scf_optimizer_diis, outer_scf_optimizer_newton, &
outer_scf_optimizer_newton_ls, outer_scf_optimizer_none, outer_scf_optimizer_sd, &
outer_scf_optimizer_secant, outer_scf_s2_constraint, plus_u_lowdin, plus_u_mulliken, &
plus_u_mulliken_charges, pw_interp, radius_covalent, radius_default, radius_single, &
radius_user, radius_vdw, random_guess, real_time_propagation, ref_charge_atomic, &
ref_charge_mulliken, rel_dkh, rel_none, rel_pot_erfc, rel_pot_full, rel_sczora_mp, &
rel_trans_atom, rel_trans_full, rel_trans_molecule, rel_zora, rel_zora_full, rel_zora_mp, &
restart_guess, rsdft_model, sccs_andreussi, sccs_derivative_cd3, sccs_derivative_cd5, &
sccs_derivative_cd7, sccs_derivative_fft, sccs_fattebert_gygi, shape_function_density, &
shape_function_gaussian, sic_ad, sic_eo, sic_list_all, sic_list_unpaired, sic_mauri_spz, &
sic_mauri_us, sic_none, slater, smear_energy_window, smear_fermi_dirac, smear_list, &
sparse_guess, spline3_nopbc_interp, spline3_pbc_interp, tddfpt_davidson, &
tddfpt_excitations, tddfpt_lanczos, tddfpt_singlet, tddfpt_triplet, use_coulomb, use_diff, &
use_no, use_restart_wfn, use_rt_restart, use_scf_wfn, wannier_projection, &
weight_type_mass, weight_type_unit, wfi_aspc_nr, wfi_frozen_method_nr, &
wfi_linear_p_method_nr, wfi_linear_ps_method_nr, wfi_linear_wf_method_nr, &
wfi_ps_method_nr, wfi_use_guess_method_nr, wfi_use_prev_p_method_nr, &
wfi_use_prev_rho_r_method_nr, wfi_use_prev_wf_method_nr, xas_1s_type, xas_2p_type, &
xas_2s_type, xas_dip_len, xas_dip_vel, xas_dscf, xas_none, xas_tp_fh, xas_tp_flex, &
xas_tp_hh, xas_tp_xfh, xas_tp_xhh, xes_tp_val
USE input_cp2k_almo, ONLY: create_almo_scf_section
USE input_cp2k_distribution, ONLY: create_distribution_section
USE input_cp2k_field, ONLY: create_efield_section,&
@ -5177,6 +5177,16 @@ CONTAINS
NULLIFY (keyword, subsection, print_key)
! add a XC section
CALL create_xc_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
! add LRI section
CALL create_lrigpw_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
CALL keyword_create(keyword, name="COLORING_METHOD", &
description="Which algorithm to use for coloring.", &
usage="COLORING_METHOD GREEDY", &
@ -5192,11 +5202,12 @@ CONTAINS
description="Algorithm to use for the calculation of the nonadditive kinetic energy.", &
usage="TNADD_METHOD ATOMIC", &
default_i_val=kg_tnadd_embed, &
enum_c_vals=s2a("EMBEDDING", "ATOMIC", "NONE"), &
enum_c_vals=s2a("EMBEDDING", "RI_EMBEDDING", "ATOMIC", "NONE"), &
enum_desc=s2a("Use full embedding potential (see Iannuzzi et al)", &
"Use full embedding potential with RI density fitting", &
"Use sum of atomic model potentials", &
"Do not use kinetic energy embedding"), &
enum_i_vals=(/kg_tnadd_embed, kg_tnadd_atomic, kg_tnadd_none/))
enum_i_vals=(/kg_tnadd_embed, kg_tnadd_embed_ri, kg_tnadd_atomic, kg_tnadd_none/))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

476
src/kg_correction.F
View file

@ -18,10 +18,12 @@ MODULE kg_correction
dbcsr_trace
USE input_constants, ONLY: kg_tnadd_atomic,&
kg_tnadd_embed,&
kg_tnadd_embed_ri,&
kg_tnadd_none
USE kg_environment_types, ONLY: kg_environment_type
USE kinds, ONLY: dp
USE lri_environment_methods, ONLY: lri_kg_rho_update
USE lri_environment_methods, ONLY: calculate_lri_densities,&
lri_kg_rho_update
USE lri_environment_types, ONLY: lri_density_type,&
lri_environment_type,&
lri_kind_type
@ -63,8 +65,6 @@ CONTAINS
!> \brief Calculates the subsystem Hohenberg-Kohn kinetic energy and the forces
!> \param qs_env ...
!> \param ks_matrix ...
!> \param gapw ...
!> \param gapw_xc ...
!> \param ekin_mol ...
!> \param calc_force ...
!> \par History
@ -72,32 +72,35 @@ CONTAINS
!> 2014.01 added atomic potential option [JGH]
!> \author Martin Haeufel and Florian Schiffmann
! **************************************************************************************************
SUBROUTINE kg_ekin_subset(qs_env, ks_matrix, gapw, gapw_xc, ekin_mol, calc_force)
SUBROUTINE kg_ekin_subset(qs_env, ks_matrix, ekin_mol, calc_force)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix
LOGICAL, INTENT(in), OPTIONAL :: gapw, gapw_xc
REAL(KIND=dp), INTENT(out) :: ekin_mol
LOGICAL :: calc_force
CHARACTER(LEN=*), PARAMETER :: routineN = 'kg_ekin_subset', routineP = moduleN//':'//routineN
LOGICAL :: lrigpw
TYPE(dft_control_type), POINTER :: dft_control
TYPE(kg_environment_type), POINTER :: kg_env
kg_env => qs_env%kg_env
CALL get_qs_env(qs_env, kg_env=kg_env, dft_control=dft_control)
lrigpw = dft_control%qs_control%lrigpw
IF (kg_env%tnadd_method == kg_tnadd_embed) THEN
CALL kg_ekin_embed(qs_env, ks_matrix, gapw, gapw_xc, &
ekin_mol, calc_force)
IF (lrigpw) THEN
CALL kg_ekin_embed_lri(qs_env, kg_env, ks_matrix, ekin_mol, calc_force)
ELSE
CALL kg_ekin_embed(qs_env, kg_env, ks_matrix, ekin_mol, calc_force)
END IF
ELSE IF (kg_env%tnadd_method == kg_tnadd_embed_ri) THEN
CALL kg_ekin_ri_embed(qs_env, kg_env, ks_matrix, ekin_mol, calc_force)
ELSE IF (kg_env%tnadd_method == kg_tnadd_atomic) THEN
CALL kg_ekin_atomic(qs_env, ks_matrix, ekin_mol)
ELSE IF (kg_env%tnadd_method == kg_tnadd_none) THEN
ekin_mol = 0.0_dp
ELSE
CPABORT("")
CPABORT("Unknown KG embedding method")
END IF
END SUBROUTINE kg_ekin_subset
@ -105,25 +108,147 @@ CONTAINS
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param kg_env ...
!> \param ks_matrix ...
!> \param gapw ...
!> \param gapw_xc ...
!> \param ekin_mol ...
!> \param calc_force ...
! **************************************************************************************************
SUBROUTINE kg_ekin_embed(qs_env, ks_matrix, gapw, gapw_xc, ekin_mol, calc_force)
SUBROUTINE kg_ekin_embed(qs_env, kg_env, ks_matrix, ekin_mol, calc_force)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(kg_environment_type), POINTER :: kg_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix
LOGICAL, INTENT(in), OPTIONAL :: gapw, gapw_xc
REAL(KIND=dp), INTENT(out) :: ekin_mol
LOGICAL :: calc_force
CHARACTER(LEN=*), PARAMETER :: routineN = 'kg_ekin_embed', routineP = moduleN//':'//routineN
INTEGER :: handle, ispin, isub, natom, nspins
LOGICAL :: use_virial
REAL(KIND=dp) :: ekin_imol
REAL(KIND=dp), DIMENSION(3, 3) :: xcvirial
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: density_matrix
TYPE(dft_control_type), POINTER :: dft_control
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_p_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: old_rho, rho_struct
TYPE(virial_type), POINTER :: virial
CALL timeset(routineN, handle)
NULLIFY (vxc_rho, vxc_tau, old_rho, rho_struct, ks_env)
! get set of molecules, natom, dft_control, pw_env
CALL get_qs_env(qs_env, &
ks_env=ks_env, &
rho=old_rho, &
natom=natom, &
dft_control=dft_control, &
virial=virial, &
para_env=para_env, &
pw_env=pw_env)
nspins = dft_control%nspins
use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
use_virial = use_virial .AND. calc_force
CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
! get the density matrix
CALL qs_rho_get(old_rho, rho_ao=density_matrix)
! allocate and initialize the density
CALL qs_rho_create(rho_struct)
! set the density matrix to the blocked matrix
CALL qs_rho_set(rho_struct, rho_ao=density_matrix) ! blocked_matrix
CALL qs_rho_rebuild(rho_struct, qs_env, rebuild_ao=.FALSE., rebuild_grids=.TRUE.)
! full density kinetic energy term
CALL qs_rho_update_rho(rho_struct, qs_env)
ekin_imol = 0.0_dp
CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=kg_env%xc_section_kg, &
vxc_rho=vxc_rho, vxc_tau=vxc_tau, exc=ekin_imol)
IF (ASSOCIATED(vxc_tau)) THEN
CPABORT(" KG with meta-kinetic energy functionals not implemented")
END IF
DO ispin = 1, nspins
vxc_rho(ispin)%pw%cr3d = vxc_rho(ispin)%pw%cr3d*vxc_rho(ispin)%pw%pw_grid%dvol
CALL integrate_v_rspace(v_rspace=vxc_rho(ispin), &
pmat=density_matrix(ispin), hmat=ks_matrix(ispin), &
qs_env=qs_env, calculate_forces=calc_force)
CALL pw_pool_give_back_pw(auxbas_pw_pool, vxc_rho(ispin)%pw)
END DO
DEALLOCATE (vxc_rho)
ekin_mol = -ekin_imol
xcvirial(1:3, 1:3) = 0.0_dp
IF (use_virial) xcvirial(1:3, 1:3) = xcvirial(1:3, 1:3)+virial%pv_xc(1:3, 1:3)
! loop over all subsets
DO isub = 1, kg_env%nsubsets
! calculate the densities for the given blocked density matrix - pass the subset task_list
CALL qs_rho_update_rho(rho_struct, qs_env, &
task_list_external=kg_env%subset(isub)%task_list)
ekin_imol = 0.0_dp
! calc Hohenberg-Kohn kin. energy of the density corresp. to the remaining molecular block(s)
CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=kg_env%xc_section_kg, &
vxc_rho=vxc_rho, vxc_tau=vxc_tau, exc=ekin_imol)
ekin_mol = ekin_mol+ekin_imol
DO ispin = 1, nspins
vxc_rho(ispin)%pw%cr3d = -vxc_rho(ispin)%pw%cr3d*vxc_rho(ispin)%pw%pw_grid%dvol
CALL integrate_v_rspace(v_rspace=vxc_rho(ispin), &
pmat=density_matrix(ispin), &
hmat=ks_matrix(ispin), &
qs_env=qs_env, &
calculate_forces=calc_force, &
task_list_external=kg_env%subset(isub)%task_list)
! clean up vxc_rho
CALL pw_pool_give_back_pw(auxbas_pw_pool, vxc_rho(ispin)%pw)
END DO
DEALLOCATE (vxc_rho)
IF (use_virial) THEN
xcvirial(1:3, 1:3) = xcvirial(1:3, 1:3)+virial%pv_xc(1:3, 1:3)
END IF
END DO
IF (use_virial) THEN
virial%pv_xc(1:3, 1:3) = xcvirial(1:3, 1:3)
END IF
! clean up rho_struct
CALL qs_rho_set(rho_struct, rho_ao=Null())
CALL qs_rho_release(rho_struct)
CALL timestop(handle)
END SUBROUTINE kg_ekin_embed
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param kg_env ...
!> \param ks_matrix ...
!> \param ekin_mol ...
!> \param calc_force ...
! **************************************************************************************************
SUBROUTINE kg_ekin_embed_lri(qs_env, kg_env, ks_matrix, ekin_mol, calc_force)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(kg_environment_type), POINTER :: kg_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix
REAL(KIND=dp), INTENT(out) :: ekin_mol
LOGICAL :: calc_force
CHARACTER(LEN=*), PARAMETER :: routineN = 'kg_ekin_embed_lri', &
routineP = moduleN//':'//routineN
INTEGER :: color, handle, iatom, ikind, imol, &
ispin, isub, natom, nkind, nspins
INTEGER, ALLOCATABLE, DIMENSION(:) :: atomlist
LOGICAL :: lrigpw, my_gapw, my_gapw_xc, use_virial
LOGICAL :: use_virial
REAL(KIND=dp) :: ekin_imol
REAL(KIND=dp), DIMENSION(3, 3) :: xcvirial
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
@ -131,7 +256,6 @@ CONTAINS
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: density_matrix, ksmat
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: pmat
TYPE(dft_control_type), POINTER :: dft_control
TYPE(kg_environment_type), POINTER :: kg_env
TYPE(lri_density_type), POINTER :: lri_density
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int
@ -146,22 +270,14 @@ CONTAINS
NULLIFY (vxc_rho, vxc_tau, old_rho, rho_struct, ks_env)
CALL get_qs_env(qs_env, kg_env=kg_env, dft_control=dft_control)
lrigpw = dft_control%qs_control%lrigpw
my_gapw = .FALSE.
IF (PRESENT(gapw)) my_gapw = gapw
my_gapw_xc = .FALSE.
IF (PRESENT(gapw_xc)) my_gapw_xc = gapw_xc
ekin_mol = 0.0_dp
CALL get_qs_env(qs_env, dft_control=dft_control)
! get set of molecules, natom, dft_control, pw_env
CALL get_qs_env(qs_env, &
ks_env=ks_env, &
rho=old_rho, &
natom=natom, &
dft_control=dft_control, &
virial=virial, &
para_env=para_env, &
pw_env=pw_env)
@ -174,84 +290,72 @@ CONTAINS
! get the density matrix
CALL qs_rho_get(old_rho, rho_ao=density_matrix)
! allocate and initialize the density
CALL qs_rho_create(rho_struct)
! set the density matrix to the blocked matrix
CALL qs_rho_set(rho_struct, rho_ao=density_matrix) ! blocked_matrix
CALL qs_rho_rebuild(rho_struct, qs_env, rebuild_ao=.FALSE., rebuild_grids=.TRUE.)
xcvirial(1:3, 1:3) = 0.0_dp
IF (lrigpw) THEN
CALL get_qs_env(qs_env, lri_env=lri_env, lri_density=lri_density, nkind=nkind)
IF (lri_env%exact_1c_terms) THEN
CPABORT(" KG with LRI and exact one-center terms not implemented")
END IF
ALLOCATE (atomlist(natom))
DO ispin = 1, nspins
lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
DO ikind = 1, nkind
lri_v_int(ikind)%v_int = 0.0_dp
IF (calc_force) THEN
lri_v_int(ikind)%v_dadr = 0.0_dp
lri_v_int(ikind)%v_dfdr = 0.0_dp
END IF
END DO
END DO
CALL get_qs_env(qs_env, lri_env=lri_env, lri_density=lri_density, nkind=nkind)
IF (lri_env%exact_1c_terms) THEN
CPABORT(" KG with LRI and exact one-center terms not implemented")
END IF
ALLOCATE (atomlist(natom))
DO ispin = 1, nspins
lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
DO ikind = 1, nkind
lri_v_int(ikind)%v_int = 0.0_dp
IF (calc_force) THEN
lri_v_int(ikind)%v_dadr = 0.0_dp
lri_v_int(ikind)%v_dfdr = 0.0_dp
END IF
END DO
END DO
! full density kinetic energy term
atomlist = 1
CALL lri_kg_rho_update(rho_struct, qs_env, lri_env, lri_density, atomlist)
ekin_imol = 0.0_dp
CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=kg_env%xc_section_kg, &
vxc_rho=vxc_rho, vxc_tau=vxc_tau, exc=ekin_imol)
IF (ASSOCIATED(vxc_tau)) THEN
CPABORT(" KG with meta-kinetic energy functionals not implemented")
END IF
DO ispin = 1, nspins
vxc_rho(ispin)%pw%cr3d = vxc_rho(ispin)%pw%cr3d*vxc_rho(ispin)%pw%pw_grid%dvol
lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
CALL integrate_v_rspace_one_center(vxc_rho(ispin), qs_env, lri_v_int, calc_force, "LRI_AUX")
CALL pw_pool_give_back_pw(auxbas_pw_pool, vxc_rho(ispin)%pw)
END DO
DEALLOCATE (vxc_rho)
ekin_mol = -ekin_imol
xcvirial(1:3, 1:3) = 0.0_dp
IF (use_virial) xcvirial(1:3, 1:3) = xcvirial(1:3, 1:3)+virial%pv_xc(1:3, 1:3)
! loop over all subsets
DO isub = 1, kg_env%nsubsets
! for gapw the task_list_soft would be needed
IF (my_gapw .OR. my_gapw_xc) &
CPABORT(" KG is currently not compatible with all electron calculations.")
IF (lrigpw) THEN
atomlist = 0
DO iatom = 1, natom
imol = kg_env%atom_to_molecule(iatom)
color = kg_env%subset_of_mol(imol)
IF (color == isub) atomlist(iatom) = 1
END DO
CALL lri_kg_rho_update(rho_struct, qs_env, lri_env, atomlist)
ELSE
! calculate the densities for the given blocked density matrix - pass the subset task_list
CALL qs_rho_update_rho(rho_struct, qs_env, &
task_list_external=kg_env%subset(isub)%task_list)
END IF
atomlist = 0
DO iatom = 1, natom
imol = kg_env%atom_to_molecule(iatom)
color = kg_env%subset_of_mol(imol)
IF (color == isub) atomlist(iatom) = 1
END DO
CALL lri_kg_rho_update(rho_struct, qs_env, lri_env, lri_density, atomlist)
ekin_imol = 0.0_dp
! calc Hohenberg-Kohn kin. energy of the density corresp. to the remaining molecular block(s)
CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=kg_env%xc_section_kg, &
vxc_rho=vxc_rho, vxc_tau=vxc_tau, exc=ekin_imol)
ekin_mol = ekin_mol+ekin_imol
DO ispin = 1, nspins
vxc_rho(ispin)%pw%cr3d = -vxc_rho(ispin)%pw%cr3d*vxc_rho(ispin)%pw%pw_grid%dvol
IF (lrigpw) THEN
lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
CALL integrate_v_rspace_one_center(vxc_rho(ispin), qs_env, &
lri_v_int, calc_force, &
"LRI_AUX", atomlist=atomlist)
ELSE
CALL integrate_v_rspace(v_rspace=vxc_rho(ispin), &
pmat=density_matrix(ispin), &
hmat=ks_matrix(ispin), &
qs_env=qs_env, &
calculate_forces=calc_force, &
gapw=gapw, &
task_list_external=kg_env%subset(isub)%task_list)
END IF
lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
CALL integrate_v_rspace_one_center(vxc_rho(ispin), qs_env, &
lri_v_int, calc_force, &
"LRI_AUX", atomlist=atomlist)
! clean up vxc_rho
CALL pw_pool_give_back_pw(auxbas_pw_pool, vxc_rho(ispin)%pw)
END DO
DEALLOCATE (vxc_rho)
@ -259,41 +363,27 @@ CONTAINS
xcvirial(1:3, 1:3) = xcvirial(1:3, 1:3)+virial%pv_xc(1:3, 1:3)
END IF
IF (ASSOCIATED(vxc_tau)) THEN
! this case has NOT been tested yet...
CALL cp_warn(__LOCATION__, &
" KG is not tested with TPSS, Becke-Roussel, or other"// &
" Meta-GGA functionals.")
DO ispin = 1, nspins
! clean up vxc_tau
CALL pw_pool_give_back_pw(auxbas_pw_pool, vxc_tau(ispin)%pw)
END DO
DEALLOCATE (vxc_tau)
END IF
END DO
IF (use_virial) THEN
virial%pv_xc(1:3, 1:3) = xcvirial(1:3, 1:3)
END IF
IF (lrigpw) THEN
CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
ALLOCATE (ksmat(1))
DO ispin = 1, nspins
lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
DO ikind = 1, nkind
CALL mp_sum(lri_v_int(ikind)%v_int, para_env%group)
END DO
ksmat(1)%matrix => ks_matrix(ispin)%matrix
CALL calculate_lri_ks_matrix(lri_env, lri_v_int, ksmat, atomic_kind_set)
CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
ALLOCATE (ksmat(1))
DO ispin = 1, nspins
lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
DO ikind = 1, nkind
CALL mp_sum(lri_v_int(ikind)%v_int, para_env%group)
END DO
IF (calc_force) THEN
pmat(1:nspins, 1:1) => density_matrix(1:nspins)
CALL calculate_lri_forces(lri_env, lri_density, qs_env, pmat, atomic_kind_set)
ENDIF
DEALLOCATE (atomlist, ksmat)
END IF
ksmat(1)%matrix => ks_matrix(ispin)%matrix
CALL calculate_lri_ks_matrix(lri_env, lri_v_int, ksmat, atomic_kind_set)
END DO
IF (calc_force) THEN
pmat(1:nspins, 1:1) => density_matrix(1:nspins)
CALL calculate_lri_forces(lri_env, lri_density, qs_env, pmat, atomic_kind_set)
ENDIF
DEALLOCATE (atomlist, ksmat)
! clean up rho_struct
CALL qs_rho_set(rho_struct, rho_ao=Null())
@ -301,7 +391,163 @@ CONTAINS
CALL timestop(handle)
END SUBROUTINE kg_ekin_embed
END SUBROUTINE kg_ekin_embed_lri
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param kg_env ...
!> \param ks_matrix ...
!> \param ekin_mol ...
!> \param calc_force ...
! **************************************************************************************************
SUBROUTINE kg_ekin_ri_embed(qs_env, kg_env, ks_matrix, ekin_mol, calc_force)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(kg_environment_type), POINTER :: kg_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix
REAL(KIND=dp), INTENT(out) :: ekin_mol
LOGICAL :: calc_force
CHARACTER(LEN=*), PARAMETER :: routineN = 'kg_ekin_ri_embed', &
routineP = moduleN//':'//routineN
INTEGER :: color, handle, iatom, ikind, imol, &
ispin, isub, natom, nkind, nspins
INTEGER, ALLOCATABLE, DIMENSION(:) :: atomlist
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: use_virial
REAL(KIND=dp) :: ekin_imol
REAL(KIND=dp), DIMENSION(3, 3) :: xcvirial
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: density_matrix
TYPE(dft_control_type), POINTER :: dft_control
TYPE(lri_density_type), POINTER :: lri_density
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_p_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho, rho_struct
TYPE(virial_type), POINTER :: virial
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env, &
ks_env=ks_env, &
rho=rho, &
natom=natom, &
nkind=nkind, &
dft_control=dft_control, &
virial=virial, &
para_env=para_env, &
pw_env=pw_env)
nspins = dft_control%nspins
use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
use_virial = use_virial .AND. calc_force
CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
! get the density matrix
CALL qs_rho_get(rho, rho_ao_kp=density_matrix)
! allocate and initialize the density
NULLIFY (rho_struct)
CALL qs_rho_create(rho_struct)
! set the density matrix to the blocked matrix
CALL qs_rho_set(rho_struct, rho_ao_kp=density_matrix)
CALL qs_rho_rebuild(rho_struct, qs_env, rebuild_ao=.FALSE., rebuild_grids=.TRUE.)
CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
ALLOCATE (cell_to_index(1, 1, 1))
cell_to_index(1, 1, 1) = 1
lri_env => kg_env%lri_env
lri_density => kg_env%lri_density
CALL calculate_lri_densities(lri_env, lri_density, qs_env, density_matrix, cell_to_index, &
rho_struct, atomic_kind_set, para_env)
kg_env%lri_density => lri_density
! full density kinetic energy term
ekin_imol = 0.0_dp
NULLIFY (vxc_rho, vxc_tau)
CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=kg_env%xc_section_kg, &
vxc_rho=vxc_rho, vxc_tau=vxc_tau, exc=ekin_imol)
IF (ASSOCIATED(vxc_tau)) THEN
CPABORT(" KG with meta-kinetic energy functionals not implemented")
END IF
DO ispin = 1, nspins
vxc_rho(ispin)%pw%cr3d = vxc_rho(ispin)%pw%cr3d*vxc_rho(ispin)%pw%pw_grid%dvol
lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
CALL integrate_v_rspace_one_center(vxc_rho(ispin), qs_env, lri_v_int, calc_force, "LRI_AUX")
CALL pw_pool_give_back_pw(auxbas_pw_pool, vxc_rho(ispin)%pw)
END DO
DEALLOCATE (vxc_rho)
ekin_mol = -ekin_imol
xcvirial(1:3, 1:3) = 0.0_dp
IF (use_virial) xcvirial(1:3, 1:3) = xcvirial(1:3, 1:3)+virial%pv_xc(1:3, 1:3)
! loop over all subsets
ALLOCATE (atomlist(natom))
DO isub = 1, kg_env%nsubsets
atomlist = 0
DO iatom = 1, natom
imol = kg_env%atom_to_molecule(iatom)
color = kg_env%subset_of_mol(imol)
IF (color == isub) atomlist(iatom) = 1
END DO
CALL lri_kg_rho_update(rho_struct, qs_env, lri_env, lri_density, atomlist)
ekin_imol = 0.0_dp
! calc Hohenberg-Kohn kin. energy of the density corresp. to the remaining molecular block(s)
NULLIFY (vxc_rho, vxc_tau)
CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=kg_env%xc_section_kg, &
vxc_rho=vxc_rho, vxc_tau=vxc_tau, exc=ekin_imol)
ekin_mol = ekin_mol+ekin_imol
DO ispin = 1, nspins
vxc_rho(ispin)%pw%cr3d = -vxc_rho(ispin)%pw%cr3d*vxc_rho(ispin)%pw%pw_grid%dvol
lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
CALL integrate_v_rspace_one_center(vxc_rho(ispin), qs_env, &
lri_v_int, calc_force, &
"LRI_AUX", atomlist=atomlist)
! clean up vxc_rho
CALL pw_pool_give_back_pw(auxbas_pw_pool, vxc_rho(ispin)%pw)
END DO
DEALLOCATE (vxc_rho)
IF (use_virial) THEN
xcvirial(1:3, 1:3) = xcvirial(1:3, 1:3)+virial%pv_xc(1:3, 1:3)
END IF
END DO
IF (use_virial) THEN
virial%pv_xc(1:3, 1:3) = xcvirial(1:3, 1:3)
END IF
ALLOCATE (ksmat(1))
DO ispin = 1, nspins
lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
DO ikind = 1, nkind
CALL mp_sum(lri_v_int(ikind)%v_int, para_env%group)
END DO
ksmat(1)%matrix => ks_matrix(ispin)%matrix
CALL calculate_lri_ks_matrix(lri_env, lri_v_int, ksmat, atomic_kind_set)
END DO
IF (calc_force) THEN
CALL calculate_lri_forces(lri_env, lri_density, qs_env, density_matrix, atomic_kind_set)
ENDIF
DEALLOCATE (atomlist, ksmat)
! clean up rho_struct
CALL qs_rho_set(rho_struct, rho_ao=Null())
CALL qs_rho_release(rho_struct)
DEALLOCATE (cell_to_index)
CALL timestop(handle)
END SUBROUTINE kg_ekin_ri_embed
! **************************************************************************************************
!> \brief ...

6
src/kg_energy_corrections.F
View file

@ -912,19 +912,19 @@ CONTAINS
! overlap
CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
CALL build_neighbor_lists(ec_env%sab_orb, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, name="sab_orb")
subcells=subcells, nlname="sab_orb")
! pseudopotential
IF (gth_potential_present .OR. sgp_potential_present) THEN
IF (ANY(ppl_present)) THEN
CALL pair_radius_setup(orb_present, ppl_present, orb_radius, ppl_radius, pair_radius)
CALL build_neighbor_lists(ec_env%sac_ppl, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="ABC", name="sac_ppl")
subcells=subcells, operator_type="ABC", nlname="sac_ppl")
END IF
IF (ANY(ppnl_present)) THEN
CALL pair_radius_setup(orb_present, ppnl_present, orb_radius, ppnl_radius, pair_radius)
CALL build_neighbor_lists(ec_env%sap_ppnl, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="ABBA", name="sap_ppnl")
subcells=subcells, operator_type="ABBA", nlname="sap_ppnl")
END IF
END IF

169
src/kg_environment.F
View file

@ -32,20 +32,21 @@ MODULE kg_environment
USE input_constants, ONLY: kg_ec_functional_harris,&
kg_tnadd_atomic,&
kg_tnadd_embed,&
kg_tnadd_embed_ri,&
kg_tnadd_none,&
xc_funct_no_shortcut,&
xc_vdw_fun_nonloc,&
xc_vdw_fun_pairpot
USE input_section_types, ONLY: &
section_vals_create, section_vals_duplicate, section_vals_get, section_vals_get_subs_vals, &
section_vals_release, section_vals_set_subs_vals, section_vals_type, section_vals_val_get, &
section_vals_val_set
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
USE kg_environment_types, ONLY: kg_environment_type
USE kg_vertex_coloring_methods, ONLY: kg_vertex_coloring
USE kinds, ONLY: dp,&
int_4,&
int_4_size,&
int_8
USE lri_environment_init, ONLY: lri_env_basis,&
lri_env_init
USE message_passing, ONLY: mp_bcast,&
mp_gather,&
mp_max
@ -118,6 +119,7 @@ CONTAINS
!> \param input ...
!> \par History
!> 2012.07 created [Martin Haeufel]
!> 2018.01 TNADD correction {JGH]
!> \author Martin Haeufel
! **************************************************************************************************
SUBROUTINE init_kg_env(qs_env, kg_env, qs_kind_set, input)
@ -128,17 +130,15 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'init_kg_env', routineP = moduleN//':'//routineN
INTEGER :: handle, ikind, ke_sections, nkind
LOGICAL :: is_set
INTEGER :: handle, ikind, nkind
REAL(KIND=dp) :: eps_pgf_orb
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dft_control_type), POINTER :: dft_control
TYPE(gto_basis_set_type), POINTER :: basis_set, harris_basis
TYPE(gto_basis_set_type), POINTER :: basis_set, harris_basis, lri_aux_basis
TYPE(qs_dispersion_type), POINTER :: dispersion_env
TYPE(qs_kind_type), POINTER :: qs_kind
TYPE(section_vals_type), POINTER :: kef_section, nl_section, pp_section, &
xc_fun_section, xc_fun_section_kg, &
TYPE(section_vals_type), POINTER :: lri_section, nl_section, pp_section, &
xc_section, xc_section_kg
CALL timeset(routineN, handle)
@ -151,6 +151,8 @@ CONTAINS
NULLIFY (kg_env%subset_of_mol)
NULLIFY (kg_env%subset)
NULLIFY (kg_env%tnadd_mat)
NULLIFY (kg_env%lri_env)
NULLIFY (kg_env%lri_density)
NULLIFY (kg_env%ec_env%sab_orb, kg_env%ec_env%sac_ppl, kg_env%ec_env%sap_ppnl)
NULLIFY (kg_env%ec_env%matrix_ks, kg_env%ec_env%matrix_h, kg_env%ec_env%matrix_s)
NULLIFY (kg_env%ec_env%matrix_t, kg_env%ec_env%matrix_p)
@ -167,6 +169,7 @@ CONTAINS
CALL section_vals_val_get(input, "DFT%KG_METHOD%COLORING_METHOD", i_val=kg_env%coloring_method)
! get method for nonadditive kinetic energy embedding potential
CALL section_vals_val_get(input, "DFT%KG_METHOD%TNADD_METHOD", i_val=kg_env%tnadd_method)
!
CALL section_vals_val_get(input, "DFT%KG_METHOD%ENERGY_CORRECTION%_SECTION_PARAMETERS_", &
l_val=kg_env%energy_correction)
IF (kg_env%energy_correction) THEN
@ -250,17 +253,8 @@ CONTAINS
xc_section_kg => section_vals_get_subs_vals(input, "DFT%KG_METHOD%ENERGY_CORRECTION%XC")
IF (ASSOCIATED(xc_section_kg)) THEN
kg_env%ec_env%xc_section => xc_section_kg
! section1 => section_vals_get_subs_vals(section, "XC")
! section2 => section_vals_get_subs_vals(section, "XC%XC_FUNCTIONAL")
! CALL xc_functionals_expand(section2, section1)
ELSE
kg_env%ec_env%xc_section => xc_section
CPABORT("KG energy correction needs XC section, with no KE functional")
! CALL cp_assert(kg_env%tnadd_method /= kg_tnadd_embed, cp_fatal_level, cp_assertion_failed, routineP, &
! "KG energy correction with full embedding requires a special XC section "// &
! "in the ENERGY_CORRECTION section. There should be no KE functional in this section.")
END IF
! dispersion
ALLOCATE (dispersion_env)
@ -280,96 +274,33 @@ CONTAINS
kg_env%ec_env%dispersion_env => dispersion_env
END IF
IF (kg_env%tnadd_method == kg_tnadd_embed) THEN
! generate a new XC section with only the KE functional
NULLIFY (xc_fun_section_kg, xc_section_kg)
xc_section => section_vals_get_subs_vals(input, "DFT%XC")
xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
CALL section_vals_duplicate(xc_section, xc_section_kg)
! test for TF, TFW, KE_GGA, and KE_LIBXC sections in turn
ke_sections = 0
NULLIFY (kef_section)
kef_section => section_vals_get_subs_vals(xc_fun_section, "TF")
CALL section_vals_get(kef_section, explicit=is_set)
IF (is_set) THEN
NULLIFY (xc_fun_section_kg)
CALL section_vals_create(xc_fun_section_kg, xc_fun_section%section)
CALL section_vals_val_set(xc_fun_section_kg, "_SECTION_PARAMETERS_", &
i_val=xc_funct_no_shortcut)
CALL section_vals_set_subs_vals(xc_fun_section_kg, "TF", &
kef_section)
CALL section_vals_set_subs_vals(xc_section_kg, "XC_FUNCTIONAL", &
xc_fun_section_kg)
ke_sections = ke_sections+1
END IF
NULLIFY (kef_section)
kef_section => section_vals_get_subs_vals(xc_fun_section, "TFW")
CALL section_vals_get(kef_section, explicit=is_set)
IF (is_set) THEN
NULLIFY (xc_fun_section_kg)
CALL section_vals_create(xc_fun_section_kg, xc_fun_section%section)
CALL section_vals_val_set(xc_fun_section_kg, "_SECTION_PARAMETERS_", &
i_val=xc_funct_no_shortcut)
CALL section_vals_set_subs_vals(xc_fun_section_kg, "TFW", &
kef_section)
CALL section_vals_set_subs_vals(xc_section_kg, "XC_FUNCTIONAL", &
xc_fun_section_kg)
ke_sections = ke_sections+1
END IF
NULLIFY (kef_section)
kef_section => section_vals_get_subs_vals(xc_fun_section, "KE_GGA")
CALL section_vals_get(kef_section, explicit=is_set)
IF (is_set) THEN
NULLIFY (xc_fun_section_kg)
CALL section_vals_create(xc_fun_section_kg, xc_fun_section%section)
CALL section_vals_val_set(xc_fun_section_kg, "_SECTION_PARAMETERS_", &
i_val=xc_funct_no_shortcut)
CALL section_vals_set_subs_vals(xc_fun_section_kg, "KE_GGA", &
kef_section)
CALL section_vals_set_subs_vals(xc_section_kg, "XC_FUNCTIONAL", &
xc_fun_section_kg)
ke_sections = ke_sections+1
END IF
NULLIFY (kef_section)
kef_section => section_vals_get_subs_vals(xc_fun_section, "KE_LIBXC")
CALL section_vals_get(kef_section, explicit=is_set)
IF (is_set) THEN
NULLIFY (xc_fun_section_kg)
CALL section_vals_create(xc_fun_section_kg, xc_fun_section%section)
CALL section_vals_val_set(xc_fun_section_kg, "_SECTION_PARAMETERS_", &
i_val=xc_funct_no_shortcut)
CALL section_vals_set_subs_vals(xc_fun_section_kg, "KE_LIBXC", &
kef_section)
CALL section_vals_set_subs_vals(xc_section_kg, "XC_FUNCTIONAL", &
xc_fun_section_kg)
ke_sections = ke_sections+1
END IF
! stop if there is more than one kinetic energy functional section present
!CPPostcondition(ke_sections==1,cp_failure_level,routineP,failure)
IF (ke_sections /= 1) &
SELECT CASE (kg_env%tnadd_method)
CASE (kg_tnadd_embed, kg_tnadd_embed_ri)
! kinetic energy functional
kg_env%xc_section_kg => section_vals_get_subs_vals(input, "DFT%KG_METHOD%XC")
IF (.NOT. ASSOCIATED(kg_env%xc_section_kg)) THEN
CALL cp_abort(__LOCATION__, &
"KG runs require a kinetic energy functional. Exactly one of the following "// &
"subsections of XC_FUNCTIONAL needs to be present: KE_GGA, TF, TFW, or KE_LIBXC.")
kg_env%xc_section_kg => xc_section_kg
NULLIFY (kef_section, xc_fun_section)
CALL section_vals_release(xc_fun_section_kg)
ELSEIF (kg_env%tnadd_method == kg_tnadd_atomic) THEN
"KG runs require a kinetic energy functional set in &KG_METHOD")
END IF
CASE (kg_tnadd_atomic, kg_tnadd_none)
NULLIFY (kg_env%xc_section_kg)
ELSEIF (kg_env%tnadd_method == kg_tnadd_none) THEN
NULLIFY (kg_env%xc_section_kg)
ELSE
CPABORT("")
CASE DEFAULT
CPABORT("KG:TNADD METHOD")
END SELECT
IF (kg_env%tnadd_method == kg_tnadd_embed_ri) THEN
! initialize the LRI environment
! Check if LRI_AUX basis is available, auto-generate if needed
nkind = SIZE(qs_kind_set)
DO ikind = 1, nkind
qs_kind => qs_kind_set(ikind)
NULLIFY (lri_aux_basis)
CALL get_qs_kind(qs_kind, basis_set=lri_aux_basis, basis_type="LRI_AUX")
CPASSERT(ASSOCIATED(lri_aux_basis))
END DO
lri_section => section_vals_get_subs_vals(input, "DFT%KG_METHOD%LRIGPW")
CALL lri_env_init(kg_env%lri_env, lri_section)
CALL lri_env_basis("LRI", qs_env, kg_env%lri_env, qs_kind_set)
END IF
CALL timestop(handle)
@ -402,7 +333,7 @@ CONTAINS
routineP = moduleN//':'//routineN
INTEGER :: handle, ikind, nkind
LOGICAL :: molecule_only
LOGICAL :: mic, molecule_only
LOGICAL, ALLOCATABLE, DIMENSION(:) :: orb_present, tpot_present
REAL(dp) :: subcells
REAL(dp), ALLOCATABLE, DIMENSION(:) :: orb_radius, tpot_radius
@ -423,8 +354,11 @@ CONTAINS
CALL timeset(routineN, handle)
NULLIFY (para_env)
! restrict lists to molecular subgroups
molecule_only = .FALSE.
IF (PRESENT(molecular)) molecule_only = molecular
! enforce minimum image convention if we use molecules
mic = molecule_only
CALL get_qs_env(qs_env=qs_env, &
atomic_kind_set=atomic_kind_set, &
@ -467,22 +401,18 @@ CONTAINS
END DO
IF (PRESENT(sab_orb)) THEN
! Build the orbital-orbital overlap neighbor list
CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
IF (PRESENT(subset_of_mol)) THEN
CALL build_neighbor_lists(sab_orb, particle_set, atom2d, cell, pair_radius, &
mic=.FALSE., subcells=subcells, molecular=molecule_only, subset_of_mol=subset_of_mol, &
current_subset=current_subset, name="sab_orb")
mic=mic, subcells=subcells, molecular=molecule_only, subset_of_mol=subset_of_mol, &
current_subset=current_subset, nlname="sab_orb")
ELSE
CALL build_neighbor_lists(sab_orb, particle_set, atom2d, cell, pair_radius, &
mic=.FALSE., subcells=subcells, molecular=molecule_only, name="sab_orb")
mic=mic, subcells=subcells, molecular=molecule_only, nlname="sab_orb")
END IF
! Print out the neighborlist
neighbor_list_section => section_vals_get_subs_vals(qs_env%input, "DFT%KG_METHOD%PRINT%NEIGHBOR_LISTS")
IF (molecule_only) THEN
CALL write_neighbor_lists(sab_orb, particle_set, cell, para_env, neighbor_list_section, &
"/SAB_ORB_MOLECULAR", "sab_orb", "MOLECULAR SUBSET NEIGHBORLIST")
@ -490,7 +420,6 @@ CONTAINS
CALL write_neighbor_lists(sab_orb, particle_set, cell, para_env, neighbor_list_section, &
"/SAB_ORB_FULL", "sab_orb", "FULL NEIGHBORLIST")
END IF
END IF
IF (PRESENT(sac_kin)) THEN
@ -504,7 +433,7 @@ CONTAINS
END DO
CALL pair_radius_setup(orb_present, tpot_present, orb_radius, tpot_radius, pair_radius)
CALL build_neighbor_lists(sac_kin, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="ABC", name="sac_kin")
subcells=subcells, operator_type="ABC", nlname="sac_kin")
neighbor_list_section => section_vals_get_subs_vals(qs_env%input, &
"DFT%KG_METHOD%PRINT%NEIGHBOR_LISTS")
CALL write_neighbor_lists(sac_kin, particle_set, cell, para_env, neighbor_list_section, &
@ -570,7 +499,7 @@ CONTAINS
CALL timestop(handle)
END SUBROUTINE
END SUBROUTINE kg_remove_duplicates
! **************************************************************************************************
!> \brief writes the graph to file using the DIMACS standard format
@ -641,7 +570,7 @@ CONTAINS
CALL timestop(handle)
END SUBROUTINE
END SUBROUTINE write_to_file
! **************************************************************************************************
!> \brief ...
@ -845,6 +774,6 @@ CONTAINS
CALL timestop(handle)
END SUBROUTINE
END SUBROUTINE kg_build_subsets
END MODULE kg_environment

20
src/kg_environment_types.F
View file

@ -13,9 +13,12 @@
MODULE kg_environment_types
USE dbcsr_api, ONLY: dbcsr_deallocate_matrix_set,&
dbcsr_p_type
USE input_section_types, ONLY: section_vals_release,&
section_vals_type
USE input_section_types, ONLY: section_vals_type
USE kinds, ONLY: dp
USE lri_environment_types, ONLY: lri_density_release,&
lri_density_type,&
lri_env_release,&
lri_environment_type
USE molecule_types, ONLY: molecule_type
USE qs_dispersion_types, ONLY: qs_dispersion_release,&
qs_dispersion_type
@ -108,6 +111,9 @@ MODULE kg_environment_types
INTEGER :: tnadd_method
LOGICAL :: energy_correction
TYPE(energy_correction_type) :: ec_env
! LRI
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(lri_density_type), POINTER :: lri_density
END TYPE kg_environment_type
CONTAINS
@ -127,8 +133,6 @@ CONTAINS
CPASSERT(ASSOCIATED(kg_env))
CALL section_vals_release(kg_env%xc_section_kg)
IF (ASSOCIATED(kg_env%sab_orb_full)) THEN
DO iab = 1, SIZE(kg_env%sab_orb_full)
CALL deallocate_neighbor_list_set(kg_env%sab_orb_full(iab)%neighbor_list_set)
@ -159,6 +163,14 @@ CONTAINS
IF (ALLOCATED(kg_env%atom_to_molecule)) DEALLOCATE (kg_env%atom_to_molecule)
! LRI
IF (ASSOCIATED(kg_env%lri_env)) THEN
CALL lri_env_release(kg_env%lri_env)
END IF
IF (ASSOCIATED(kg_env%lri_density)) THEN
CALL lri_density_release(kg_env%lri_density)
END IF
! energy correction
IF (kg_env%energy_correction) THEN
! neighbor lists

19
src/lri_environment_init.F
View file

@ -34,8 +34,7 @@ MODULE lri_environment_init
pi
USE mathlib, ONLY: invert_matrix
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type,&
set_qs_env
qs_environment_type
USE qs_kind_types, ONLY: get_qs_kind,&
qs_kind_type
#include "./base/base_uses.f90"
@ -56,18 +55,17 @@ CONTAINS
! **************************************************************************************************
!> \brief initializes the lri env
!> \param qs_env ...
!> \param lri_env ...
!> \param lri_section ...
! **************************************************************************************************
SUBROUTINE lri_env_init(qs_env, lri_section)
SUBROUTINE lri_env_init(lri_env, lri_section)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(section_vals_type), POINTER :: lri_section
CHARACTER(len=*), PARAMETER :: routineN = 'lri_env_init', routineP = moduleN//':'//routineN
REAL(KIND=dp), DIMENSION(:), POINTER :: radii
TYPE(lri_environment_type), POINTER :: lri_env
NULLIFY (lri_env)
CALL lri_env_create(lri_env)
@ -111,8 +109,6 @@ CONTAINS
lri_env%r_in = radii(1)
lri_env%r_out = radii(2)
CALL set_qs_env(qs_env, lri_env=lri_env)
CALL cite_reference(Golze2017b)
IF (lri_env%use_shg_integrals) CALL cite_reference(Golze2017a)
@ -121,12 +117,14 @@ CONTAINS
!> \brief initializes the lri env
!> \param ri_type ...
!> \param qs_env ...
!> \param lri_env ...
!> \param qs_kind_set ...
! **************************************************************************************************
SUBROUTINE lri_env_basis(ri_type, qs_env, qs_kind_set)
SUBROUTINE lri_env_basis(ri_type, qs_env, lri_env, qs_kind_set)
CHARACTER(len=*), INTENT(IN) :: ri_type
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
CHARACTER(len=*), PARAMETER :: routineN = 'lri_env_basis', routineP = moduleN//':'//routineN
@ -139,9 +137,6 @@ CONTAINS
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(dft_control_type), POINTER :: dft_control
TYPE(gto_basis_set_type), POINTER :: orb_basis_set, ri_basis_set
TYPE(lri_environment_type), POINTER :: lri_env
CALL get_qs_env(qs_env=qs_env, lri_env=lri_env)
! initialize the basic basis sets (orb and ri)
CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)

20
src/lri_environment_methods.F
View file

@ -62,8 +62,7 @@ MODULE lri_environment_methods
USE pw_types, ONLY: pw_p_type
USE qs_collocate_density, ONLY: calculate_lri_rho_elec
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type,&
set_qs_env
qs_environment_type
USE qs_force_types, ONLY: qs_force_type
USE qs_kind_types, ONLY: qs_kind_type
USE qs_neighbor_list_types, ONLY: get_iterator_info,&
@ -337,13 +336,15 @@ CONTAINS
!> \param rho_struct ...
!> \param qs_env ...
!> \param lri_env ...
!> \param lri_density ...
!> \param atomlist ...
! **************************************************************************************************
SUBROUTINE lri_kg_rho_update(rho_struct, qs_env, lri_env, atomlist)
SUBROUTINE lri_kg_rho_update(rho_struct, qs_env, lri_env, lri_density, atomlist)
TYPE(qs_rho_type), POINTER :: rho_struct
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(lri_density_type), POINTER :: lri_density
INTEGER, DIMENSION(:), INTENT(IN) :: atomlist
CHARACTER(LEN=*), PARAMETER :: routineN = 'lri_kg_rho_update', &
@ -354,7 +355,6 @@ CONTAINS
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_type) :: pmat_diag
TYPE(dft_control_type), POINTER :: dft_control
TYPE(lri_density_type), POINTER :: lri_density
TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g, rho_r
CALL timeset(routineN, handle)
@ -362,7 +362,6 @@ CONTAINS
CPASSERT(ASSOCIATED(rho_struct))
CALL get_qs_env(qs_env, dft_control=dft_control, para_env=para_env)
CALL get_qs_env(qs_env, lri_density=lri_density)
CALL qs_rho_get(rho_struct, rho_r=rho_r, rho_g=rho_g, tot_rho_r=tot_rho_r)
@ -561,6 +560,7 @@ CONTAINS
!> \brief performs the fitting of the density and distributes the fitted
!> density on the grid
!> \param lri_env the lri environment
!> \param lri_density ...
!> \param qs_env ...
!> \param pmatrix ...
!> \param cell_to_index ...
@ -568,10 +568,11 @@ CONTAINS
!> \param atomic_kind_set ...
!> \param para_env ...
! **************************************************************************************************
SUBROUTINE calculate_lri_densities(lri_env, qs_env, pmatrix, cell_to_index, &
SUBROUTINE calculate_lri_densities(lri_env, lri_density, qs_env, pmatrix, cell_to_index, &
lri_rho_struct, atomic_kind_set, para_env)
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(lri_density_type), POINTER :: lri_density
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: pmatrix
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
@ -582,12 +583,7 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_lri_densities', &
routineP = moduleN//':'//routineN
TYPE(lri_density_type), POINTER :: lri_density
CALL get_qs_env(qs_env, lri_density=lri_density)
CALL calculate_avec_lri(lri_env, lri_density, pmatrix, cell_to_index)
CALL set_qs_env(qs_env, lri_density=lri_density)
CALL distribute_lri_density_on_the_grid(lri_env, lri_density, qs_env, &
lri_rho_struct, atomic_kind_set, para_env)
@ -963,8 +959,6 @@ CONTAINS
IF (lri_env%exact_1c_terms) CALL dbcsr_release(pmat_diag)
ENDDO
CALL set_qs_env(qs_env, lri_density=lri_density)
DEALLOCATE (atom_of_kind)
END IF

1
src/lri_environment_types.F
View file

@ -1401,6 +1401,7 @@ CONTAINS
DEALLOCATE (lri_rhos)
END IF
NULLIFY (lri_rhos)
END SUBROUTINE deallocate_lri_rhos

4
src/mp2_gpw.F
View file

@ -2295,7 +2295,7 @@ CONTAINS
CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
CALL build_neighbor_lists(sab_orb_sub, particle_set, atom2d, cell, pair_radius, &
mic=.FALSE., subcells=subcells, molecular=.FALSE., name="sab_orb_sub")
mic=.FALSE., subcells=subcells, molecular=.FALSE., nlname="sab_orb_sub")
CALL atom2d_cleanup(atom2d)
DEALLOCATE (atom2d)
DEALLOCATE (orb_present, orb_radius, pair_radius)
@ -2435,7 +2435,7 @@ CONTAINS
CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
CALL build_neighbor_lists(sab_orb_all, particle_set, atom2d, cell, pair_radius, &
mic=.FALSE., subcells=subcells, molecular=.FALSE., name="sab_orb_sub")
mic=.FALSE., subcells=subcells, molecular=.FALSE., nlname="sab_orb_sub")
CALL atom2d_cleanup(atom2d)
DEALLOCATE (atom2d)
DEALLOCATE (orb_present, orb_radius, pair_radius)

112
src/qs_energy_init.F
View file

@ -20,15 +20,18 @@ MODULE qs_energy_init
USE efield_utils, ONLY: calculate_ecore_efield
USE input_constants, ONLY: kg_tnadd_atomic,&
kg_tnadd_embed,&
kg_tnadd_embed_ri,&
kg_tnadd_none,&
outer_scf_hirshfeld_constraint
USE input_section_types, ONLY: section_vals_type
USE kg_environment, ONLY: kg_build_neighborlist,&
kg_build_subsets
USE kg_environment_types, ONLY: kg_environment_type
USE kinds, ONLY: dp
USE kpoint_methods, ONLY: kpoint_init_cell_index
USE kpoint_types, ONLY: kpoint_type
USE lri_environment_methods, ONLY: build_lri_matrices
USE lri_environment_methods, ONLY: build_lri_matrices,&
calculate_lri_integrals
USE lri_environment_types, ONLY: lri_environment_type
USE molecule_types, ONLY: molecule_of_atom,&
molecule_type
@ -94,17 +97,22 @@ CONTAINS
NULLIFY (ks_env, matrix_w, matrix_w_mp2, matrix_s, dft_control)
CALL qs_energies_init_kg(qs_env, molecule_only)
CALL get_qs_env(qs_env, dft_control=dft_control)
IF (dft_control%qs_control%do_kg) THEN
molecule_only = .TRUE.
CALL qs_energies_init_kg(qs_env)
ELSE
molecule_only = .FALSE.
END IF
CALL qs_energies_init_hamiltonians(qs_env, calc_forces, molecule_only)
CALL qs_ks_allocate_basics(qs_env)
CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
! if need forces allocate energy weighted density matrices
CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
IF (calc_forces .AND. .NOT. has_unit_metric) THEN
CALL get_qs_env(qs_env, &
ks_env=ks_env, &
matrix_s_kp=matrix_s, &
dft_control=dft_control)
matrix_s_kp=matrix_s)
nspin = dft_control%nspins
nimg = dft_control%nimages
matrix => matrix_s(1, 1)%matrix
@ -136,14 +144,12 @@ CONTAINS
!> \brief Refactoring of qs_energies_scf. Puts initialization of the Kim-Gordon
!> settings into separate subroutine
!> \param qs_env ...
!> \param molecule_only ...
!> \par History
!> 05.2013 created [Florian Schiffmann]
! **************************************************************************************************
SUBROUTINE qs_energies_init_kg(qs_env, molecule_only)
SUBROUTINE qs_energies_init_kg(qs_env)
TYPE(qs_environment_type), POINTER :: qs_env
LOGICAL :: molecule_only
CHARACTER(len=*), PARAMETER :: routineN = 'qs_energies_init_kg', &
routineP = moduleN//':'//routineN
@ -151,57 +157,63 @@ CONTAINS
INTEGER :: handle, isubset, natom
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dft_control_type), POINTER :: dft_control
TYPE(kg_environment_type), POINTER :: kg_env
TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: soo_list
CALL timeset(routineN, handle)
NULLIFY (dft_control, para_env)
molecule_only = .FALSE.
CALL get_qs_env(qs_env, dft_control=dft_control, para_env=para_env)
IF (dft_control%qs_control%do_kg) THEN
CPASSERT(dft_control%qs_control%do_kg)
! create neighbor lists with molecular blocks
molecule_only = .TRUE.
kg_env => qs_env%kg_env
! get the set of molecules
CALL get_qs_env(qs_env=qs_env, molecule_set=molecule_set, natom=natom)
qs_env%kg_env%natom = natom
! store set of molecules in kg_env
qs_env%kg_env%molecule_set => molecule_set
! build the (new) full neighborlist
CALL kg_build_neighborlist(qs_env, sab_orb=qs_env%kg_env%sab_orb_full)
IF (.NOT. ALLOCATED(qs_env%kg_env%atom_to_molecule)) THEN
ALLOCATE (qs_env%kg_env%atom_to_molecule(natom))
! get the mapping from atoms to molecules
CALL molecule_of_atom(molecule_set, atom_to_mol=qs_env%kg_env%atom_to_molecule)
END IF
IF (qs_env%kg_env%tnadd_method == kg_tnadd_embed) THEN
! allocate the subset list
IF (.NOT. ASSOCIATED(qs_env%kg_env%subset_of_mol)) ALLOCATE (qs_env%kg_env%subset_of_mol(SIZE(molecule_set)))
CALL kg_build_subsets(qs_env%kg_env, para_env)
DO isubset = 1, qs_env%kg_env%nsubsets
! build the (new) molecular neighborlist of the current subset
CALL kg_build_neighborlist(qs_env, sab_orb=qs_env%kg_env%subset(isubset)%sab_orb, molecular=.TRUE., &
subset_of_mol=qs_env%kg_env%subset_of_mol, current_subset=isubset)
END DO
ELSE IF (qs_env%kg_env%tnadd_method == kg_tnadd_atomic) THEN
! build the A-C list for the nonadditive kinetic energy potential
CALL kg_build_neighborlist(qs_env, sac_kin=qs_env%kg_env%sac_kin)
ELSE IF (qs_env%kg_env%tnadd_method == kg_tnadd_none) THEN
! nothing to do
ELSE
CPABORT("KG:TNADD METHOD")
END IF
! get the set of molecules
CALL get_qs_env(qs_env=qs_env, molecule_set=molecule_set, natom=natom)
kg_env%natom = natom
! store set of molecules in kg_env
kg_env%molecule_set => molecule_set
! build the (new) full neighborlist
CALL kg_build_neighborlist(qs_env, sab_orb=kg_env%sab_orb_full)
IF (.NOT. ALLOCATED(kg_env%atom_to_molecule)) THEN
ALLOCATE (kg_env%atom_to_molecule(natom))
! get the mapping from atoms to molecules
CALL molecule_of_atom(molecule_set, atom_to_mol=kg_env%atom_to_molecule)
END IF
SELECT CASE (kg_env%tnadd_method)
CASE (kg_tnadd_embed, kg_tnadd_embed_ri)
! allocate the subset list
IF (.NOT. ASSOCIATED(kg_env%subset_of_mol)) THEN
ALLOCATE (kg_env%subset_of_mol(SIZE(molecule_set)))
END IF
!
CALL kg_build_subsets(kg_env, para_env)
!
DO isubset = 1, kg_env%nsubsets
! build the (new) molecular neighborlist of the current subset
CALL kg_build_neighborlist(qs_env, sab_orb=kg_env%subset(isubset)%sab_orb, molecular=.TRUE., &
subset_of_mol=kg_env%subset_of_mol, current_subset=isubset)
END DO
CASE (kg_tnadd_atomic)
! build the A-C list for the nonadditive kinetic energy potential
CALL kg_build_neighborlist(qs_env, sac_kin=kg_env%sac_kin)
CASE (kg_tnadd_none)
! nothing to do
CASE DEFAULT
CPABORT("KG:TNADD METHOD")
END SELECT
! LRI neighborlist
IF (kg_env%tnadd_method == kg_tnadd_embed_ri) THEN
NULLIFY (soo_list)
CALL kg_build_neighborlist(qs_env, sab_orb=soo_list, molecular=.TRUE.)
kg_env%lri_env%soo_list => soo_list
CALL calculate_lri_integrals(kg_env%lri_env, qs_env)
END IF
CALL timestop(handle)
END SUBROUTINE qs_energies_init_kg

24
src/qs_environment.F
View file

@ -86,7 +86,7 @@ MODULE qs_environment
do_method_gapw_xc, do_method_gpw, do_method_lrigpw, do_method_mndo, do_method_mndod, &
do_method_ofgpw, do_method_pdg, do_method_pm3, do_method_pm6, do_method_pm6fm, &
do_method_pnnl, do_method_rigpw, do_method_rm1, do_qmmm_gauss, do_qmmm_swave, &
general_roks, rel_none, rel_trans_atom, vdw_pairpot_dftd3, wfi_aspc_nr, &
general_roks, kg_tnadd_embed_ri, rel_none, rel_trans_atom, vdw_pairpot_dftd3, wfi_aspc_nr, &
wfi_linear_ps_method_nr, wfi_linear_wf_method_nr, wfi_ps_method_nr, &
wfi_use_guess_method_nr, xc_vdw_fun_none, xc_vdw_fun_nonloc, xc_vdw_fun_pairpot
USE input_section_types, ONLY: section_vals_get,&
@ -242,6 +242,7 @@ CONTAINS
TYPE(cp_blacs_env_type), POINTER :: blacs_env
TYPE(dft_control_type), POINTER :: dft_control
TYPE(kpoint_type), POINTER :: kpoints
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_ks_env_type), POINTER :: ks_env
@ -341,11 +342,13 @@ CONTAINS
CALL get_qs_env(qs_env, dft_control=dft_control)
IF (method_id == do_method_lrigpw .OR. dft_control%qs_control%lri_optbas) THEN
CALL lri_env_basis("LRI", qs_env, qs_kind_set)
CALL get_qs_env(qs_env=qs_env, lri_env=lri_env)
CALL lri_env_basis("LRI", qs_env, lri_env, qs_kind_set)
ELSE IF (method_id == do_method_rigpw) THEN
CALL cp_warn(__LOCATION__, "Experimental code: "// &
"RIGPW should only be used for testing.")
CALL lri_env_basis("RI", qs_env, qs_kind_set)
CALL get_qs_env(qs_env=qs_env, lri_env=lri_env)
CALL lri_env_basis("RI", qs_env, lri_env, qs_kind_set)
END IF
! more kpoint stuff
@ -465,10 +468,10 @@ CONTAINS
INTEGER :: handle, ikind, ispin, iw, lmax_sphere, maxl, maxlgto, maxlgto_lri, maxlppl, &
maxlppnl, method_id, multiplicity, my_ival, n_ao, n_ao_aux_fit, n_mo_add, natom, &
nelectron, nkind, output_unit
nelectron, nkind, output_unit, tnadd_method
INTEGER, DIMENSION(2) :: n_mo, nelectron_spin
LOGICAL :: all_potential_present, be_silent, &
do_kpoints, has_unit_metric, &
do_kpoints, has_unit_metric, lribas, &
was_present
REAL(dp) :: ewald_rcut, maxocc, verlet_skin
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
@ -630,7 +633,8 @@ CONTAINS
IF (dft_control%qs_control%method_id == do_method_lrigpw .OR. &
dft_control%qs_control%lri_optbas .OR. &
dft_control%qs_control%method_id == do_method_rigpw) THEN
CALL lri_env_init(qs_env, lri_section)
CALL lri_env_init(lri_env, lri_section)
CALL set_qs_env(qs_env, lri_env=lri_env)
END IF
! *** Check basis and fill in missing parts ***
@ -667,7 +671,13 @@ CONTAINS
END DO
END IF
IF (dft_control%qs_control%method_id == do_method_lrigpw .OR. dft_control%qs_control%lri_optbas) THEN
lribas = .FALSE.
IF (dft_control%qs_control%method_id == do_method_lrigpw) lribas = .TRUE.
IF (dft_control%qs_control%do_kg) THEN
CALL section_vals_val_get(dft_section, "KG_METHOD%TNADD_METHOD", i_val=tnadd_method)
IF (tnadd_method == kg_tnadd_embed_ri) lribas = .TRUE.
END IF
IF (lribas) THEN
! Check if LRI_AUX basis is available, auto-generate if needed
CALL get_qs_env(qs_env, nkind=nkind, lri_env=lri_env)
DO ikind = 1, nkind

8
src/qs_ks_methods.F
View file

@ -648,10 +648,10 @@ CONTAINS
! calculate KG correction
IF (dft_control%qs_control%do_kg .AND. just_energy) THEN
CPASSERT(.NOT. (gapw .OR. gapw_xc))
CPASSERT(nimages == 1)
ksmat => ks_matrix(:, 1)
CALL kg_ekin_subset(qs_env, ksmat, gapw, gapw_xc, ekin_mol, &
calculate_forces)
CALL kg_ekin_subset(qs_env, ksmat, ekin_mol, calculate_forces)
! substract kg corr from the total energy
energy%exc = energy%exc-ekin_mol
@ -694,10 +694,10 @@ CONTAINS
becke, cdft_control, calculate_forces)
IF (dft_control%qs_control%do_kg) THEN
CPASSERT(.NOT. (gapw .OR. gapw_xc))
CPASSERT(nimages == 1)
ksmat => ks_matrix(:, 1)
CALL kg_ekin_subset(qs_env, ksmat, gapw, gapw_xc, ekin_mol, &
calculate_forces)
CALL kg_ekin_subset(qs_env, ksmat, ekin_mol, calculate_forces)
! substract kg corr from the total energy
energy%exc = energy%exc-ekin_mol
! virial corrections

107
src/qs_neighbor_lists.F
View file

@ -377,12 +377,8 @@ CONTAINS
neighbor_list_section => section_vals_get_subs_vals(force_env_section, "DFT%PRINT%NEIGHBOR_LISTS")
molecule_only = .FALSE.
IF (PRESENT(molecular)) molecule_only = molecular
! This sets the id number of the qs neighbor lists, new lists, means new version
! new version implies new sparsity of the matrices
last_qs_neighbor_list_id_nr = last_qs_neighbor_list_id_nr+1
CALL set_ks_env(ks_env=ks_env, neighbor_list_id=last_qs_neighbor_list_id_nr)
@ -414,6 +410,19 @@ CONTAINS
rigpw = (dft_control%qs_control%method_id == do_method_rigpw)
lri_optbas = dft_control%qs_control%lri_optbas
! molecular lists
molecule_only = .FALSE.
IF (PRESENT(molecular)) molecule_only = molecular
! minimum image convention (MIC)
mic = molecule_only
IF (dokp) THEN
! no MIC for kpoints
mic = .FALSE.
ELSEIF (nddo .OR. dftb) THEN
! enforce MIC for interaction lists in SE and DFTB
mic = .TRUE.
END IF
hfx_sections => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
CALL section_vals_get(hfx_sections, explicit=do_hfx)
@ -426,12 +435,9 @@ CONTAINS
CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
! Allocate work storage
nkind = SIZE(atomic_kind_set)
ALLOCATE (orb_present(nkind), aux_fit_present(nkind), aux_present(nkind), &
default_present(nkind), core_present(nkind))
ALLOCATE (orb_radius(nkind), aux_fit_radius(nkind), c_radius(nkind), &
core_radius(nkind), calpha(nkind), zeff(nkind))
orb_radius(:) = 0.0_dp
@ -442,19 +448,16 @@ CONTAINS
zeff(:) = 0.0_dp
ALLOCATE (pair_radius(nkind, nkind))
IF (gth_potential_present .OR. sgp_potential_present) THEN
ALLOCATE (ppl_present(nkind), ppl_radius(nkind))
ppl_radius = 0.0_dp
ALLOCATE (ppnl_present(nkind), ppnl_radius(nkind))
ppnl_radius = 0.0_dp
END IF
IF (paw_atom_present) THEN
ALLOCATE (oce_present(nkind), oce_radius(nkind))
oce_radius = 0.0_dp
END IF
IF (all_potential_present .OR. sgp_potential_present) THEN
ALLOCATE (all_present(nkind), all_pot_rad(nkind))
all_pot_rad = 0.0_dp
@ -462,7 +465,6 @@ CONTAINS
! Initialize the local data structures
ALLOCATE (atom2d(nkind))
CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
molecule_set, molecule_only, particle_set=particle_set)
@ -568,22 +570,10 @@ CONTAINS
END DO
IF (dokp) THEN
! no MIC for kpoints
mic = .FALSE.
ELSE
! enforce MIC for interaction lists in SE and DFTB
IF (nddo .OR. dftb) THEN
mic = .TRUE.
ELSE
mic = .FALSE.
END IF
END IF
! Build the orbital-orbital overlap neighbor lists
CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
CALL build_neighbor_lists(sab_orb, particle_set, atom2d, cell, pair_radius, &
mic=mic, subcells=subcells, molecular=molecule_only, name="sab_orb")
mic=mic, subcells=subcells, molecular=molecule_only, nlname="sab_orb")
CALL set_ks_env(ks_env=ks_env, sab_orb=sab_orb)
CALL write_neighbor_lists(sab_orb, particle_set, cell, para_env, neighbor_list_section, &
"/SAB_ORB", "sab_orb", "ORBITAL ORBITAL")
@ -593,15 +583,15 @@ CONTAINS
! might be not optimal. It should be verified for each operator.
IF (.NOT. (nddo .OR. dftb)) THEN
CALL build_neighbor_lists(sab_all, particle_set, atom2d, cell, pair_radius, &
mic=mic, symmetric=.FALSE., subcells=subcells, molecular=molecule_only, name="sab_all")
mic=mic, symmetric=.FALSE., subcells=subcells, molecular=molecule_only, nlname="sab_all")
CALL set_ks_env(ks_env=ks_env, sab_all=sab_all)
ENDIF
! Build the core-core overlap neighbor lists
IF (.NOT. (nddo .OR. dftb)) THEN
CALL pair_radius_setup(core_present, core_present, core_radius, core_radius, pair_radius)
CALL build_neighbor_lists(sab_core, particle_set, atom2d, cell, pair_radius, mic=mic, subcells=subcells, &
operator_type="PP", name="sab_core")
CALL build_neighbor_lists(sab_core, particle_set, atom2d, cell, pair_radius, subcells=subcells, &
operator_type="PP", nlname="sab_core")
CALL set_ks_env(ks_env=ks_env, sab_core=sab_core)
CALL write_neighbor_lists(sab_core, particle_set, cell, para_env, neighbor_list_section, &
"/SAB_CORE", "sab_core", "CORE CORE")
@ -610,14 +600,14 @@ CONTAINS
IF (dft_control%do_admm) THEN
CALL pair_radius_setup(aux_fit_present, aux_fit_present, aux_fit_radius, aux_fit_radius, pair_radius)
CALL build_neighbor_lists(sab_aux_fit, particle_set, atom2d, cell, pair_radius, &
mic=mic, molecular=molecule_only, subcells=subcells, name="sab_aux_fit")
mic=mic, molecular=molecule_only, subcells=subcells, nlname="sab_aux_fit")
CALL build_neighbor_lists(sab_aux_fit_asymm, particle_set, atom2d, cell, pair_radius, &
mic=mic, symmetric=.FALSE., molecular=molecule_only, subcells=subcells, &
name="sab_aux_fit_asymm")
nlname="sab_aux_fit_asymm")
CALL pair_radius_setup(aux_fit_present, orb_present, aux_fit_radius, orb_radius, pair_radius)
CALL build_neighbor_lists(sab_aux_fit_vs_orb, particle_set, atom2d, cell, pair_radius, &
mic=mic, symmetric=.FALSE., molecular=molecule_only, subcells=subcells, &
name="sab_aux_fit_vs_orb")
nlname="sab_aux_fit_vs_orb")
CALL set_ks_env(ks_env=ks_env, sab_aux_fit=sab_aux_fit)
CALL set_ks_env(ks_env=ks_env, sab_aux_fit_vs_orb=sab_aux_fit_vs_orb)
@ -660,7 +650,7 @@ CONTAINS
CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
END IF
CALL build_neighbor_lists(sab_kp, particle_set, atom2d, cell, pair_radius, &
mic=mic, subcells=subcells, molecular=molecule_only, name="sab_kp")
subcells=subcells, nlname="sab_kp")
CALL set_ks_env(ks_env=ks_env, sab_kp=sab_kp)
END IF
@ -669,14 +659,14 @@ CONTAINS
IF (ANY(ppl_present)) THEN
CALL pair_radius_setup(orb_present, ppl_present, orb_radius, ppl_radius, pair_radius)
CALL build_neighbor_lists(sac_ppl, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="ABC", name="sac_ppl")
subcells=subcells, operator_type="ABC", nlname="sac_ppl")
CALL set_ks_env(ks_env=ks_env, sac_ppl=sac_ppl)
CALL write_neighbor_lists(sac_ppl, particle_set, cell, para_env, neighbor_list_section, &
"/SAC_PPL", "sac_ppl", "ORBITAL GTH-PPL")
IF (lrigpw) THEN
IF (qs_env%lri_env%ppl_ri) THEN
CALL build_neighbor_lists(sac_lri, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, symmetric=.FALSE., operator_type="PP", name="sac_lri")
subcells=subcells, symmetric=.FALSE., operator_type="PP", nlname="sac_lri")
CALL set_ks_env(ks_env=ks_env, sac_lri=sac_lri)
END IF
END IF
@ -685,7 +675,7 @@ CONTAINS
IF (ANY(ppnl_present)) THEN
CALL pair_radius_setup(orb_present, ppnl_present, orb_radius, ppnl_radius, pair_radius)
CALL build_neighbor_lists(sap_ppnl, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="ABBA", name="sap_ppnl")
subcells=subcells, operator_type="ABBA", nlname="sap_ppnl")
CALL set_ks_env(ks_env=ks_env, sap_ppnl=sap_ppnl)
CALL write_neighbor_lists(sap_ppnl, particle_set, cell, para_env, neighbor_list_section, &
"/SAP_PPNL", "sap_ppnl", "ORBITAL GTH-PPNL")
@ -697,7 +687,7 @@ CONTAINS
IF (ANY(oce_present)) THEN
CALL pair_radius_setup(orb_present, oce_present, orb_radius, oce_radius, pair_radius)
CALL build_neighbor_lists(sap_oce, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="ABBA", name="sap_oce")
subcells=subcells, operator_type="ABBA", nlname="sap_oce")
CALL set_ks_env(ks_env=ks_env, sap_oce=sap_oce)
CALL write_neighbor_lists(sap_oce, particle_set, cell, para_env, neighbor_list_section, &
"/SAP_OCE", "sap_oce", "ORBITAL(A) PAW-PRJ")
@ -709,7 +699,7 @@ CONTAINS
IF (all_potential_present .OR. sgp_potential_present) THEN
CALL pair_radius_setup(orb_present, all_present, orb_radius, all_pot_rad, pair_radius)
CALL build_neighbor_lists(sac_ae, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="ABC", name="sac_ae")
subcells=subcells, operator_type="ABC", nlname="sac_ae")
CALL set_ks_env(ks_env=ks_env, sac_ae=sac_ae)
CALL write_neighbor_lists(sac_ae, particle_set, cell, para_env, neighbor_list_section, &
"/SAC_AE", "sac_ae", "ORBITAL ERFC POTENTIAL")
@ -725,10 +715,10 @@ CONTAINS
IF (dft_control%qs_control%se_control%do_ewald_gks) THEN
! Use MIC for the periodic code of GKS
CALL build_neighbor_lists(sab_se, particle_set, atom2d, cell, pair_radius, mic=mic, &
subcells=subcells, name="sab_se")
subcells=subcells, nlname="sab_se")
ELSE
CALL build_neighbor_lists(sab_se, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, name="sab_se")
subcells=subcells, nlname="sab_se")
END IF
CALL set_ks_env(ks_env=ks_env, sab_se=sab_se)
CALL write_neighbor_lists(sab_se, particle_set, cell, para_env, neighbor_list_section, &
@ -740,7 +730,7 @@ CONTAINS
c_radius = dft_control%qs_control%se_control%cutoff_lrc
CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
CALL build_neighbor_lists(sab_lrc, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, name="sab_lrc")
subcells=subcells, nlname="sab_lrc")
CALL set_ks_env(ks_env=ks_env, sab_lrc=sab_lrc)
CALL write_neighbor_lists(sab_lrc, particle_set, cell, para_env, neighbor_list_section, &
"/SAB_LRC", "sab_lrc", "SE LONG-RANGE CORRECTION")
@ -755,7 +745,7 @@ CONTAINS
c_radius = 0.5_dp*SQRT(-LOG(3.5_dp*alpha**3*1.e-12_dp))/alpha
CALL pair_radius_setup(orb_present, orb_present, c_radius, c_radius, pair_radius)
CALL build_neighbor_lists(sab_tbe, particle_set, atom2d, cell, pair_radius, mic=mic, &
subcells=subcells, name="sab_tbe")
subcells=subcells, nlname="sab_tbe")
CALL set_ks_env(ks_env=ks_env, sab_tbe=sab_tbe)
END IF
@ -764,13 +754,12 @@ CONTAINS
IF (dft_control%qs_control%dftb_control%dispersion_type == dispersion_uff) THEN
DO ikind = 1, nkind
CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom)
CALL get_dftb_atom_param(dftb_parameter=dftb_atom, &
rcdisp=c_radius(ikind))
CALL get_dftb_atom_param(dftb_parameter=dftb_atom, rcdisp=c_radius(ikind))
END DO
default_present = .TRUE.
CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
CALL build_neighbor_lists(sab_vdw, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, name="sab_vdw")
subcells=subcells, nlname="sab_vdw")
CALL set_ks_env(ks_env=ks_env, sab_vdw=sab_vdw)
END IF
END IF
@ -785,7 +774,7 @@ CONTAINS
default_present = .TRUE. !include all atoms in vdW (even without basis)
CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
CALL build_neighbor_lists(sab_vdw, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="PP", name="sab_vdw")
subcells=subcells, operator_type="PP", nlname="sab_vdw")
dispersion_env%sab_vdw => sab_vdw
IF (dispersion_env%pp_type == vdw_pairpot_dftd3 .OR. &
@ -797,7 +786,7 @@ CONTAINS
END DO
CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
CALL build_neighbor_lists(sab_cn, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="PP", name="sab_cn")
subcells=subcells, operator_type="PP", nlname="sab_cn")
dispersion_env%sab_cn => sab_cn
END IF
END IF
@ -807,7 +796,7 @@ CONTAINS
CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
soo_list => qs_env%lri_env%soo_list
CALL build_neighbor_lists(soo_list, particle_set, atom2d, cell, pair_radius, &
mic=mic, molecular=molecule_only, subcells=subcells, name="soo_list")
mic=mic, molecular=molecule_only, subcells=subcells, nlname="soo_list")
qs_env%lri_env%soo_list => soo_list
CALL write_neighbor_lists(soo_list, particle_set, cell, para_env, neighbor_list_section, &
"/SOO_LIST", "soo_list", "ORBITAL ORBITAL (RI)")
@ -828,20 +817,20 @@ CONTAINS
CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
soo_list => qs_env%lri_env%soo_list
CALL build_neighbor_lists(soo_list, particle_set, atom2d, cell, pair_radius, &
mic=mic, molecular=molecule_only, subcells=subcells, name="soo_list")
mic=mic, molecular=molecule_only, subcells=subcells, nlname="soo_list")
qs_env%lri_env%soo_list => soo_list
!
CALL pair_radius_setup(ri_present, ri_present, ri_radius, ri_radius, pair_radius)
saa_list => qs_env%lri_env%saa_list
CALL build_neighbor_lists(saa_list, particle_set, atom2d, cell, pair_radius, &
mic=mic, molecular=molecule_only, subcells=subcells, name="saa_list")
mic=mic, molecular=molecule_only, subcells=subcells, nlname="saa_list")
qs_env%lri_env%saa_list => saa_list
!
CALL pair_radius_setup(ri_present, orb_present, ri_radius, orb_radius, pair_radius)
soa_list => qs_env%lri_env%soa_list
CALL build_neighbor_lists(soa_list, particle_set, atom2d, cell, pair_radius, &
mic=mic, symmetric=.FALSE., molecular=molecule_only, &
subcells=subcells, operator_type="ABC", name="saa_list")
subcells=subcells, operator_type="ABC", nlname="saa_list")
qs_env%lri_env%soa_list => soa_list
END IF
@ -856,7 +845,7 @@ CONTAINS
default_present = .TRUE. !include all atoms (even without basis)
CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
CALL build_neighbor_lists(sab_almo, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="PP", name="sab_almo")
subcells=subcells, operator_type="PP", nlname="sab_almo")
CALL set_ks_env(ks_env=ks_env, sab_almo=sab_almo)
ENDIF
@ -913,7 +902,7 @@ CONTAINS
!> \param subset_of_mol ...
!> \param current_subset ...
!> \param operator_type ...
!> \param name ...
!> \param nlname ...
!> \date 20.03.2002
!> \par History
!> - Major refactoring (25.07.2010,jhu)
@ -922,7 +911,7 @@ CONTAINS
! **************************************************************************************************
SUBROUTINE build_neighbor_lists(ab_list, particle_set, atom, cell, pair_radius, subcells, &
mic, symmetric, molecular, subset_of_mol, current_subset, &
operator_type, name)
operator_type, nlname)
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: ab_list
@ -935,7 +924,7 @@ CONTAINS
INTEGER, DIMENSION(:), OPTIONAL, POINTER :: subset_of_mol
INTEGER, OPTIONAL :: current_subset
CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: operator_type
CHARACTER(LEN=*), INTENT(IN) :: name
CHARACTER(LEN=*), INTENT(IN) :: nlname
CHARACTER(len=*), PARAMETER :: routineN = 'build_neighbor_lists', &
routineP = moduleN//':'//routineN
@ -964,7 +953,7 @@ CONTAINS
POINTER :: subcell
REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: r_pbc
CALL timeset(routineN//"_"//TRIM(name), handle)
CALL timeset(routineN//"_"//TRIM(nlname), handle)
! input options
my_mic = .FALSE.
@ -1408,7 +1397,7 @@ CONTAINS
molecule_set, my_molecular, particle_set=particle_set)
CALL build_neighbor_lists(ab_list, particle_set, atom2d, cell, pair_radius, &
mic=my_mic, symmetric=my_symmetric, molecular=my_molecular, &
subcells=2.0_dp, name="AUX_NL")
subcells=2.0_dp, nlname="AUX_NL")
CALL atom2d_cleanup(atom2d)
@ -1617,7 +1606,7 @@ CONTAINS
!> \param neighbor_list_section ...
!> \param nl_type ...
!> \param middle_name ...
!> \param name ...
!> \param nlname ...
!> \date 04.03.2002
!> \par History
!> - Adapted to the new parallelized neighbor list version
@ -1626,7 +1615,7 @@ CONTAINS
!> \version 1.0
! **************************************************************************************************
SUBROUTINE write_neighbor_lists(ab, particle_set, cell, para_env, neighbor_list_section, &
nl_type, middle_name, name)
nl_type, middle_name, nlname)
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: ab
@ -1634,7 +1623,7 @@ CONTAINS
TYPE(cell_type), POINTER :: cell
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(section_vals_type), POINTER :: neighbor_list_section
CHARACTER(LEN=*), INTENT(IN) :: nl_type, middle_name, name
CHARACTER(LEN=*), INTENT(IN) :: nl_type, middle_name, nlname
CHARACTER(LEN=default_string_length) :: string, unit_str
INTEGER :: iatom, inode, iw, jatom, mype, &
@ -1666,7 +1655,7 @@ CONTAINS
! Print headline
string = ""
WRITE (UNIT=string, FMT="(A,I5,A)") &
TRIM(name)//" IN "//TRIM(unit_str)//" (PROCESS", mype, ")"
TRIM(nlname)//" IN "//TRIM(unit_str)//" (PROCESS", mype, ")"
CALL compress(string)
IF (iw > 0) WRITE (UNIT=iw, FMT="(/,/,T2,A)") TRIM(string)

11
src/qs_rho_methods.F
View file

@ -29,7 +29,8 @@ MODULE qs_rho_methods
USE kpoint_types, ONLY: get_kpoint_info,&
kpoint_type
USE lri_environment_methods, ONLY: calculate_lri_densities
USE lri_environment_types, ONLY: lri_environment_type
USE lri_environment_types, ONLY: lri_density_type,&
lri_environment_type
USE pw_env_types, ONLY: pw_env_get,&
pw_env_type
USE pw_methods, ONLY: pw_zero
@ -44,7 +45,8 @@ MODULE qs_rho_methods
USE qs_collocate_density, ONLY: calculate_drho_elec,&
calculate_rho_elec
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
qs_environment_type,&
set_qs_env
USE qs_ks_types, ONLY: get_ks_env,&
qs_ks_env_type
USE qs_local_rho_types, ONLY: local_rho_type
@ -349,6 +351,7 @@ CONTAINS
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp, rho_xc_ao
TYPE(dft_control_type), POINTER :: dft_control
TYPE(kpoint_type), POINTER :: kpoints
TYPE(lri_density_type), POINTER :: lri_density
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_p_type), DIMENSION(:), POINTER :: drho_g, drho_r, drho_xc_g, rho_g, rho_r, &
@ -405,10 +408,12 @@ CONTAINS
CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
CALL calculate_lri_densities(lri_env, qs_env, rho_ao_kp, cell_to_index, &
CALL get_qs_env(qs_env, lri_density=lri_density)
CALL calculate_lri_densities(lri_env, lri_density, qs_env, rho_ao_kp, cell_to_index, &
lri_rho_struct=rho_struct, &
atomic_kind_set=atomic_kind_set, &
para_env=para_env)
CALL set_qs_env(qs_env, lri_density=lri_density)
CALL qs_rho_set(rho_struct, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
ELSEIF (dft_control%qs_control%rigpw) THEN
CPASSERT(.NOT. dft_control%use_kinetic_energy_density)

2
src/qs_tddfpt2_subgroups.F
View file

@ -604,7 +604,7 @@ CONTAINS
NULLIFY (sab)
CALL section_vals_val_get(input, "DFT%SUBCELLS", r_val=subcells)
CALL build_neighbor_lists(sab, particle_set, atom2d, cell, pair_radius, &
mic=.FALSE., subcells=subcells, molecular=.FALSE., name="sab_orb")
mic=.FALSE., subcells=subcells, molecular=.FALSE., nlname="sab_orb")
CALL atom2d_cleanup(atom2d)
DEALLOCATE (atom2d, orb_present, orb_radius, pair_radius)

9
src/qs_update_s_mstruct.F
View file

@ -16,7 +16,8 @@ MODULE qs_update_s_mstruct
USE cp_ddapc_util, ONLY: cp_ddapc_init
USE input_constants, ONLY: do_ppl_analytic,&
do_ppl_grid,&
kg_tnadd_embed
kg_tnadd_embed,&
kg_tnadd_embed_ri
USE pw_methods, ONLY: pw_transfer
USE pw_types, ONLY: pw_p_type,&
pw_release
@ -211,7 +212,8 @@ CONTAINS
IF (dft_control%qs_control%do_kg) THEN
soft_valid = (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc)
IF (qs_env%kg_env%tnadd_method == kg_tnadd_embed) THEN
IF (qs_env%kg_env%tnadd_method == kg_tnadd_embed .OR. &
qs_env%kg_env%tnadd_method == kg_tnadd_embed_ri) THEN
IF (ASSOCIATED(qs_env%kg_env%subset)) THEN
DO isub = 1, qs_env%kg_env%nsubsets
@ -223,15 +225,12 @@ CONTAINS
END IF
DO isub = 1, qs_env%kg_env%nsubsets
CALL allocate_task_list(qs_env%kg_env%subset(isub)%task_list)
! generate the subset task list from the neighborlist
CALL generate_qs_task_list(ks_env, qs_env%kg_env%subset(isub)%task_list, &
reorder_rs_grid_ranks=.FALSE., soft_valid=soft_valid, &
skip_load_balance_distributed=skip_load_balance_distributed, &
sab_orb_external=qs_env%kg_env%subset(isub)%sab_orb)
END DO
END IF

8
tests/QS/regtest-kg/2H2O_ecmao.inp
View file

@ -20,6 +20,12 @@
PURIFICATION_METHOD TRS4
&END
&KG_METHOD
&XC
&XC_FUNCTIONAL
&TF
&END
&END
&END XC
&ENERGY_CORRECTION
HARRIS_BASIS HARRIS
MAO
@ -39,8 +45,6 @@
&XC_FUNCTIONAL
&PADE
&END
&TF
&END
&END
&END XC
&END DFT

8
tests/QS/regtest-kg/2H2O_ecmao2.inp
View file

@ -20,6 +20,12 @@
PURIFICATION_METHOD TRS4
&END
&KG_METHOD
&XC
&XC_FUNCTIONAL
&TF
&END
&END
&END XC
&ENERGY_CORRECTION
HARRIS_BASIS HARRIS
MAO
@ -41,8 +47,6 @@
&XC_FUNCTIONAL
&PADE
&END
&TF
&END
&END
&END XC
&END DFT

15
tests/QS/regtest-kg/H2-libxc-diag.inp
View file

@ -22,13 +22,22 @@
MU -0.1
PURIFICATION_METHOD TRS4
&END
&KG_METHOD
&XC
&XC_FUNCTIONAL
&KE_LIBXC
FUNCTIONAL XC_GGA_K_LLP
&END
&END
&XC_GRID
XC_DERIV SPLINE2
&END
&END XC
&END KG_METHOD
&XC
&XC_FUNCTIONAL
&PADE
&END
&KE_LIBXC
FUNCTIONAL XC_GGA_K_LLP
&END
&END
&XC_GRID
XC_DERIV SPLINE2

15
tests/QS/regtest-kg/H2-libxc-ot.inp
View file

@ -26,13 +26,22 @@
MU -0.1
PURIFICATION_METHOD TRS4
&END
&KG_METHOD
&XC
&XC_FUNCTIONAL
&KE_LIBXC
FUNCTIONAL XC_GGA_K_LLP
&END
&END
&XC_GRID
XC_DERIV SPLINE2
&END
&END XC
&END KG_METHOD
&XC
&XC_FUNCTIONAL
&PADE
&END
&KE_LIBXC
FUNCTIONAL XC_GGA_K_LLP
&END
&END
&XC_GRID
XC_DERIV SPLINE2

15
tests/QS/regtest-kg/H2-libxc.inp
View file

@ -22,13 +22,22 @@
MU -0.1
PURIFICATION_METHOD TRS4
&END
&KG_METHOD
&XC
&XC_FUNCTIONAL
&KE_LIBXC
FUNCTIONAL XC_GGA_K_LLP
&END
&END
&XC_GRID
XC_DERIV SPLINE2
&END
&END XC
&END KG_METHOD
&XC
&XC_FUNCTIONAL
&PADE
&END
&KE_LIBXC
FUNCTIONAL XC_GGA_K_LLP
&END
&END
&XC_GRID
XC_DERIV SPLINE2

76
tests/QS/regtest-kg/H2_H2O-kglri.inp Normal file
View file

@ -0,0 +1,76 @@
&FORCE_EVAL
METHOD Quickstep
&PRINT
&FORCES
&END
&END
&DFT
BASIS_SET_FILE_NAME BASIS_MOLOPT
POTENTIAL_FILE_NAME GTH_POTENTIALS
&MGRID
CUTOFF 400
&END MGRID
&QS
METHOD GPW
EPS_DEFAULT 5.0E-14
EPS_KG_ORB 1.0E-6
LS_SCF
KG_METHOD
&END QS
&LS_SCF
MAX_SCF 40
EPS_FILTER 1.0E-6
EPS_SCF 1.0E-7
MU -0.1
PURIFICATION_METHOD TRS4
&END
&SCF
EPS_SCF 1.0E-5
MAX_SCF 50
&END
&KG_METHOD
&XC
&XC_FUNCTIONAL
&KE_GGA
FUNCTIONAL T92
&END
&END
&END XC
&LRIGPW
RI_STATISTIC T
EPS_O3_INT 1.e-14
&END LRIGPW
&END KG_METHOD
&XC
&XC_FUNCTIONAL
&PADE
&END
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC 4.0 6.0 3.0
&END CELL
&COORD
O 0.000000 0.000000 -0.065587 H2O
H 0.000000 -0.757136 0.520545 H2O
H 0.000000 0.757136 0.520545 H2O
H 2.000000 2.500000 0.000000 H2
H 2.000000 -2.500000 0.000000 H2
&END COORD
&KIND H
BASIS_SET DZVP-MOLOPT-GTH
POTENTIAL GTH-PADE-q1
&END KIND
&KIND O
BASIS_SET DZVP-MOLOPT-GTH
POTENTIAL GTH-PADE-q6
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT H2
RUN_TYPE ENERGY_FORCE
PRINT_LEVEL MEDIUM
&END GLOBAL

14
tests/QS/regtest-kg/H2_H2O-lri.inp
View file

@ -7,8 +7,6 @@
&DFT
BASIS_SET_FILE_NAME BASIS_MOLOPT
POTENTIAL_FILE_NAME GTH_POTENTIALS
&KG_METHOD
&END
&MGRID
CUTOFF 400
&END MGRID
@ -39,13 +37,19 @@
EPS_SCF 1.0E-5
MAX_SCF 50
&END
&KG_METHOD
&XC
&XC_FUNCTIONAL
&KE_GGA
FUNCTIONAL T92
&END
&END
&END XC
&END KG_METHOD
&XC
&XC_FUNCTIONAL
&PADE
&END
&KE_GGA
FUNCTIONAL T92
&END
&END
&END XC
&END DFT

15
tests/QS/regtest-kg/H2_H2O-vdW.inp
View file

@ -27,13 +27,22 @@
EPS_SCF 1.0E-8
SCF_GUESS ATOMIC
&END
&KG_METHOD
&XC
&XC_FUNCTIONAL
&KE_GGA
FUNCTIONAL LLP
&END
&END
&XC_GRID
XC_DERIV SPLINE2
&END
&END XC
&END KG_METHOD
&XC
&XC_FUNCTIONAL
&PADE
&END
&KE_GGA
FUNCTIONAL LLP
&END
&END
&XC_GRID
XC_DERIV SPLINE2

15
tests/QS/regtest-kg/H2_H2O-xcLC.inp
View file

@ -27,13 +27,22 @@
EPS_SCF 1.0E-8
SCF_GUESS ATOMIC
&END
&KG_METHOD
&XC
&XC_FUNCTIONAL
&KE_GGA
FUNCTIONAL LC
&END
&END
&XC_GRID
XC_DERIV SPLINE2
&END
&END XC
&END KG_METHOD
&XC
&XC_FUNCTIONAL
&PADE
&END
&KE_GGA
FUNCTIONAL LC
&END
&END
&XC_GRID
XC_DERIV SPLINE2

15
tests/QS/regtest-kg/H2_H2O-xcLLP.inp
View file

@ -27,13 +27,22 @@
EPS_SCF 1.0E-8
SCF_GUESS ATOMIC
&END
&KG_METHOD
&XC
&XC_FUNCTIONAL
&KE_GGA
FUNCTIONAL LLP
&END
&END
&XC_GRID
XC_DERIV SPLINE2
&END
&END XC
&END KG_METHOD
&XC
&XC_FUNCTIONAL
&PADE
&END
&KE_GGA
FUNCTIONAL LLP
&END
&END
&XC_GRID
XC_DERIV SPLINE2

25
tests/QS/regtest-kg/H2_H2O-xcLLP_ec.inp
View file

@ -24,13 +24,23 @@
PURIFICATION_METHOD TRS4
&END
&KG_METHOD
&XC
&XC_FUNCTIONAL
&KE_GGA
FUNCTIONAL LLP
&END
&END
&XC_GRID
XC_DERIV SPLINE2
&END
&END XC
&ENERGY_CORRECTION
&XC
&XC_FUNCTIONAL
&PADE
&END
&END
&END XC
&XC
&XC_FUNCTIONAL
&PADE
&END
&END
&END XC
&END ENERGY_CORRECTION
&END KG_METHOD
&SCF
@ -41,9 +51,6 @@
&XC_FUNCTIONAL
&PADE
&END
&KE_GGA
FUNCTIONAL LLP
&END
&END
&XC_GRID
XC_DERIV SPLINE2

15
tests/QS/regtest-kg/H2_H2O-xcPW86.inp
View file

@ -27,13 +27,22 @@
EPS_SCF 1.0E-8
SCF_GUESS ATOMIC
&END
&KG_METHOD
&XC
&XC_FUNCTIONAL
&KE_GGA
FUNCTIONAL PW86
&END
&END
&XC_GRID
XC_DERIV SPLINE2
&END
&END XC
&END KG_METHOD
&XC
&XC_FUNCTIONAL
&PADE
&END
&KE_GGA
FUNCTIONAL PW86
&END
&END
&XC_GRID
XC_DERIV SPLINE2

15
tests/QS/regtest-kg/H2_H2O-xcPW91.inp
View file

@ -27,13 +27,22 @@
EPS_SCF 1.0E-8
SCF_GUESS ATOMIC
&END
&KG_METHOD
&XC
&XC_FUNCTIONAL
&KE_GGA
FUNCTIONAL PW91
&END
&END
&XC_GRID
XC_DERIV SPLINE2
&END
&END XC
&END KG_METHOD
&XC
&XC_FUNCTIONAL
&PADE
&END
&KE_GGA
FUNCTIONAL PW91
&END
&END
&XC_GRID
XC_DERIV SPLINE2

15
tests/QS/regtest-kg/H2_H2O-xcT92.inp
View file

@ -27,13 +27,22 @@
EPS_SCF 1.0E-8
SCF_GUESS ATOMIC
&END
&KG_METHOD
&XC
&XC_FUNCTIONAL
&KE_GGA
FUNCTIONAL T92
&END
&END
&XC_GRID
XC_DERIV SPLINE2
&END
&END XC
&END KG_METHOD
&XC
&XC_FUNCTIONAL
&PADE
&END
&KE_GGA
FUNCTIONAL T92
&END
&END
&XC_GRID
XC_DERIV SPLINE2

8
tests/QS/regtest-kg/H2_H2O_ec.inp
View file

@ -20,6 +20,12 @@
PURIFICATION_METHOD TRS4
&END
&KG_METHOD
&XC
&XC_FUNCTIONAL
&TF
&END
&END
&END XC
&ENERGY_CORRECTION
HARRIS_BASIS HARRIS
&XC
@ -38,8 +44,6 @@
&XC_FUNCTIONAL
&PADE
&END
&TF
&END
&END
&END XC
&END DFT

8
tests/QS/regtest-kg/H2_H2O_ecprim.inp
View file

@ -20,6 +20,12 @@
PURIFICATION_METHOD TRS4
&END
&KG_METHOD
&XC
&XC_FUNCTIONAL
&TF
&END
&END
&END XC
&ENERGY_CORRECTION
HARRIS_BASIS PRIMITIVE
&XC
@ -38,8 +44,6 @@
&XC_FUNCTIONAL
&PADE
&END
&TF
&END
&END
&END XC
&END DFT

3
tests/QS/regtest-kg/H2_H2O_lsks.inp
View file

@ -34,9 +34,6 @@
&XC_FUNCTIONAL
&PADE
&END
# &KE_GGA
# FUNCTIONAL LLP
# &END
&END
&XC_GRID
XC_DERIV SPLINE2

15
tests/QS/regtest-kg/H2_MD-2.inp
View file

@ -27,13 +27,22 @@
EPS_SCF 1.0E-8
SCF_GUESS ATOMIC
&END
&KG_METHOD
&XC
&XC_FUNCTIONAL
&KE_GGA
FUNCTIONAL LLP
&END
&END
&XC_GRID
XC_DERIV SPLINE2
&END
&END XC
&END KG_METHOD
&XC
&XC_FUNCTIONAL
&PADE
&END
&KE_GGA
FUNCTIONAL LLP
&END
&END
&XC_GRID
XC_DERIV SPLINE2

13
tests/QS/regtest-kg/H2_MD-3.inp
View file

@ -18,6 +18,16 @@
&END QS
&KG_METHOD
COLORING_METHOD GREEDY
&XC
&XC_FUNCTIONAL
&KE_GGA
FUNCTIONAL LLP
&END
&END
&XC_GRID
XC_DERIV SPLINE2
&END
&END XC
&END KG_METHOD
&LS_SCF
MAX_SCF 40
@ -34,9 +44,6 @@
&XC_FUNCTIONAL
&PADE
&END
&KE_GGA
FUNCTIONAL LLP
&END
&END
&XC_GRID
XC_DERIV SPLINE2

15
tests/QS/regtest-kg/H2_MD.inp
View file

@ -26,13 +26,22 @@
EPS_SCF 1.0E-8
SCF_GUESS ATOMIC
&END
&KG_METHOD
&XC
&XC_FUNCTIONAL
&KE_GGA
FUNCTIONAL LLP
&END
&END
&XC_GRID
XC_DERIV SPLINE2
&END
&END XC
&END KG_METHOD
&XC
&XC_FUNCTIONAL
&PADE
&END
&KE_GGA
FUNCTIONAL LLP
&END
&END
&XC_GRID
XC_DERIV SPLINE2

29
tests/QS/regtest-kg/TEST_FILES
View file

@ -4,26 +4,27 @@
H2_MD.inp 11 7e-07 -3.154884771619098
H2_MD-2.inp 11 7e-07 -3.1548873963582911
H2_MD-3.inp 11 7e-07 -2.1282196550651968
H2-libxc.inp 11 4e-07 -3.0639455259826192
H2-libxc-ot.inp 11 4e-07 -3.0639455258035309
H2-libxc-diag.inp 11 4e-07 -3.0639455259860058
geo-phase-1.inp 11 5e-07 -2.229647716115335
H2-libxc.inp 11 4e-07 -3.0453523116483912
H2-libxc-ot.inp 11 4e-07 -3.0453523114395449
H2-libxc-diag.inp 11 4e-07 -3.0453523116522888
geo-phase-1.inp 11 5e-07 -2.229629661035424
geo-phase-2.inp 0
H2_KG-1.inp 11 7e-07 -3.0857150864890102
H2_KG-2.inp 11 1e-06 -3.2269982118329681
H2_H2O-xcLC.inp 11 1e-10 -18.128250540244913
H2_H2O-xcLLP.inp 11 1e-10 -18.127917783808059
H2_H2O-xcPW86.inp 11 1e-10 -18.127360510720912
H2_H2O-xcPW91.inp 11 1e-10 -18.126985946101961
H2_H2O-xcT92.inp 11 1e-10 -18.128137231215568
H2_H2O-xcLLP_ec.inp 66 1e-08 -18.1302536445
H2_H2O-xcLC.inp 11 1e-10 -18.143825547885832
H2_H2O-xcLLP.inp 11 1e-10 -18.143534964747040
H2_H2O-xcPW86.inp 11 1e-10 -18.143043340247019
H2_H2O-xcPW91.inp 11 1e-10 -18.142668216945395
H2_H2O-xcT92.inp 11 1e-10 -18.143724634599963
H2_H2O-xcLLP_ec.inp 66 1e-08 -18.1323810074
H2_H2O_ks.inp 11 1e-10 -18.130229622911152
H2_H2O_lsks.inp 11 1e-10 -18.130225854904928
H2_H2O_ec.inp 66 1e-10 -18.4642032438
H2_H2O_ecprim.inp 66 1e-10 -18.4070965525
H2_H2O_ec.inp 66 1e-10 -18.4641086282
H2_H2O_ecprim.inp 66 1e-10 -18.4071025433
2H2O_ecmao.inp 66 1e-10 -34.0838901711
2H2O_ecmao2.inp 66 1e-08 -34.5024811487
H2-none.inp 11 7e-07 -3.359680469888914
H2_H2O-vdW.inp 11 1e-10 -18.133386209975370
H2_H2O-lri.inp 72 3e-02 0.00018044
H2_H2O-vdW.inp 11 1e-10 -18.149003390914348
H2_H2O-lri.inp 72 3e-02 0.00011458
H2_H2O-kglri.inp 11 3e-02 -18.315900381303461
#EOF

10
tests/QS/regtest-kg/geo-phase-1.inp
View file

@ -63,14 +63,18 @@
! satur 27 colors (~250s) , greedy 37 (~1s)
COLORING_METHOD DSATUR
! COLORING_METHOD GREEDY
&XC
&XC_FUNCTIONAL
&KE_GGA
FUNCTIONAL LLP
&END
&END XC_FUNCTIONAL
&END XC
&END
&XC
&XC_FUNCTIONAL
&PBE
&END
&KE_GGA
FUNCTIONAL LLP
&END
&END XC_FUNCTIONAL
&END XC
&END DFT

10
tests/QS/regtest-kg/geo-phase-2.inp
View file

@ -63,14 +63,18 @@
! satur 27 colors (~250s) , greedy 37 (~1s)
COLORING_METHOD DSATUR
! COLORING_METHOD GREEDY
&XC
&XC_FUNCTIONAL
&KE_GGA
FUNCTIONAL LLP
&END
&END XC_FUNCTIONAL
&END XC
&END
&XC
&XC_FUNCTIONAL
&PBE
&END
&KE_GGA
FUNCTIONAL LLP
&END
&END XC_FUNCTIONAL
&END XC
&END DFT

2
tests/QS/regtest-ls-rtp/TEST_FILES
View file

@ -13,7 +13,7 @@ H2O-ls-rtp-scf-restart.inp 1 6e-14 -17.08907732640686
H2O-ls-emd-scf-restart.inp 1 1e-11 -17.0890686606
H2O-ls-emd-mcweeny.inp 1 1e-11 -17.0890686606
H2O-ls-emd-mcweeny-rt-restart.inp 1 1e-11 -17.0890598734
emd-KG.inp 1 1e-14 -18.10682510355022
emd-KG.inp 1 1e-14 -18.11820168793887
H2plus-ls-rtp.inp 1 2e-12 -0.52809706445746
H2plus-ls-rtp2.inp 1 3e-13 -0.53624161700164
H2O-ls-emd-mixing.inp 1 2e-11 -17.0889404255

12
tests/QS/regtest-ls-rtp/emd-KG.inp
View file

@ -40,13 +40,19 @@
MAX_ITER_LANCZOS 30
SIGN_SQRT_ORDER 3
&END LS_SCF
&KG_METHOD
&XC
&XC_FUNCTIONAL
&KE_GGA
FUNCTIONAL LLP
&END
&END XC_FUNCTIONAL
&END XC
&END KG_METHOD
&XC
&XC_FUNCTIONAL
&PBE
&END
&KE_GGA
FUNCTIONAL LLP
&END
&END XC_FUNCTIONAL
&END XC
&END DFT