Skip to content
Permalink
Browse files

XAS_TDP| Bug fixes

  • Loading branch information...
abussy authored and pseewald committed Jul 24, 2019
1 parent f224986 commit f1cd4a37765e969984604a9e58dba34a8dd9097f
@@ -69,8 +69,9 @@ MODULE hfx_libint_wrapper
#endif
END TYPE

!This is a hack because 3-center ERIs are bugged on libint 2.5
TYPE(C_FUNPTR), DIMENSION(0:5,0:5,0:7), BIND(C) :: libint2_build_3eri
!This is a safeguard/hack because 3-center ERIs are bugged on libint 2.5
TYPE(C_FUNPTR), DIMENSION(0:LIBINT2_MAX_AM,0:LIBINT2_MAX_AM,0:LIBINT2_MAX_AM_3eri), &
BIND(C) :: libint2_build_3eri

CONTAINS

@@ -96,8 +96,8 @@ MODULE input_cp2k_dft
xas_2s_type, xas_3d_type, xas_3p_type, xas_3s_type, xas_4d_type, xas_4f_type, xas_4p_type, &
xas_4s_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, xas_not_excited, xas_tdp_by_kind, &
xas_tdp_by_index, tddfpt_spin_cons, tddfpt_spin_flip, do_hfx_potential_coulomb, &
do_hfx_potential_truncated, do_hfx_potential_short, xas_tdp_diag_std, xas_tdp_diag_syevr
xas_tdp_by_index, tddfpt_spin_cons, tddfpt_spin_flip, do_potential_coulomb, &
do_potential_truncated, do_potential_short, xas_tdp_diag_std, xas_tdp_diag_syevr
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,&
@@ -7885,7 +7885,7 @@ SUBROUTINE create_xas_tdp_section(section)
description="XAS simulations using TDDFPT. Excitation from specified "//&
"core orbitals are considered one at a time. In case of high symmetry "//&
"structures, donor core orbitals should be localized.", &
n_keywords=11, n_subsections=5, repeats=.FALSE.)
n_keywords=11, n_subsections=4, repeats=.FALSE.)

NULLIFY (keyword, subsection, print_key)

@@ -7907,6 +7907,7 @@ SUBROUTINE create_xas_tdp_section(section)
"option allows for cheap verification of the input parameters", &
usage="CHECK_ONLY {logical}", &
default_l_val=.FALSE., &
lone_keyword_l_val=.TRUE., &
repeats=.FALSE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
@@ -7946,7 +7947,8 @@ SUBROUTINE create_xas_tdp_section(section)
"oscillator strenghts (in the length representation with "// &
"the origin set on the relevant excited atom)", &
usage="QUADRUPOLE {logical}", &
default_l_val=.FALSE.)
default_l_val=.FALSE.,&
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

@@ -7955,7 +7957,7 @@ SUBROUTINE create_xas_tdp_section(section)
CALL section_create(subsection, __LOCATION__, name="DONOR_STATES", &
description="Specifications for the donor states from which core "// &
"electrons are excited", &
n_keywords=5, &
n_keywords=6, &
n_subsections=0, &
repeats=.FALSE.)

@@ -8008,13 +8010,23 @@ SUBROUTINE create_xas_tdp_section(section)
CALL keyword_create(keyword, __LOCATION__, name="N_SEARCH", &
description="Number of MOs (per spin) to search to find specified donor core "//&
"orbitals, starting from the lowest in energy and upward. By default, "//&
"all HOMOs are searched. If the LOCALIZE subsection is used, "//&
"all HOMOs are searched. If the LOCALIZE keyword is used, "//&
"then all searched states are first localized.",&
usage="N_SEARCH {integer}",&
default_i_val=-1)
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="LOCALIZE", &
variants=s2a("LOC", "DO_LOC"), &
description="Whether the N_SEARCH potential donor states are localized."//&
"Necessary in case of symmetry equivalent excited atoms.", &
usage="LOCALIZE {logical}", &
default_l_val=.FALSE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)

CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
! End of the donor states subsection
@@ -8034,7 +8046,8 @@ SUBROUTINE create_xas_tdp_section(section)
variants=(/"TDA"/), &
description="Whether the Tamm-Dancoff approximation should be used.", &
usage="TAMM_DANCOFF {logical}", &
default_l_val=.FALSE.)
default_l_val=.TRUE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

@@ -8177,13 +8190,13 @@ SUBROUTINE create_xas_tdp_section(section)
"Coulomb operator cannot be used in periodic systems.", &
usage="OPERATOR {string}", &
repeats=.FALSE., &
default_i_val=do_hfx_potential_coulomb, &
default_i_val=do_potential_coulomb, &
enum_c_vals=s2a("COULOMB", "TRUNCATED", "SHORTRANGE"), &
enum_desc=s2a("Stendard Coulomb operator: 1/r", &
"Truncated Coulomb operator: 1/r if (r<R_c), 0 otherwise ",&
"Short range: erfc(omega*r)/r"),&
enum_i_vals=(/do_hfx_potential_coulomb, do_hfx_potential_truncated, &
do_hfx_potential_short/))
enum_i_vals=(/do_potential_coulomb, do_potential_truncated, &
do_potential_short/))
CALL section_add_keyword(subsubsection, keyword)
CALL keyword_release(keyword)

@@ -8245,10 +8258,6 @@ SUBROUTINE create_xas_tdp_section(section)
CALL section_release(subsection)
! End of Kernel subsection

CALL create_localize_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)

CALL section_create(subsection, __LOCATION__, "PRINT", "Controls the printing of information during "//&
"XAS TDP calculations", repeats=.FALSE.)
CALL cp_print_key_section_create(print_key, __LOCATION__, name="PROGRAM_RUN_INFO",&
@@ -13,8 +13,8 @@ MODULE libint_3center
USE gamma, ONLY: fgamma => fgamma_0
USE hfx_libint_wrapper, ONLY: cp_libint_set_params_eri, cp_libint_get_3eris, &
cp_libint_t, prim_data_f_size
USE input_constants, ONLY: do_hfx_potential_coulomb, do_hfx_potential_short, &
do_hfx_potential_truncated
USE input_constants, ONLY: do_potential_coulomb, do_potential_short, &
do_potential_truncated
USE kinds, ONLY: dp
USE mathconstants, ONLY: pi
USE orbital_pointers, ONLY: nco, ncoset
@@ -252,15 +252,15 @@ SUBROUTINE set_params(lib, ri, rj, rk, li_max, lj_max, lk_max, zeti, zetj, zetk,
params%Rho = zetk*gammaq/(zetk+gammaq)

SELECT CASE(op)
CASE (do_hfx_potential_coulomb)
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

params%Fm = 0.0_dp
CALL fgamma(params%m_max, T, params%Fm)
params%Fm = prefac*params%Fm
CASE (do_hfx_potential_truncated)
CASE (do_potential_truncated)
R = r_cutoff*SQRT(params%Rho)
T = params%Rho*SUM((params%Q-rk)**2)
S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj-ri)**2))
@@ -270,7 +270,7 @@ SUBROUTINE set_params(lib, ri, rj, rk, li_max, lj_max, lk_max, zeti, zetj, zetk,
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_hfx_potential_short)
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
@@ -571,7 +571,7 @@ SUBROUTINE build_atomic_relmat(matrix_h, atomic_kind_set, qs_kind_set, particle_
IF (iatom == jatom) THEN
ikind = kind_of(iatom)
CALL get_qs_kind(qs_kind_set(ikind), reltmat=reltmat)
hblock = hblock+reltmat
IF(ASSOCIATED(reltmat)) hblock = hblock+reltmat
END IF
END DO
CALL dbcsr_iterator_stop(iter)
@@ -802,6 +802,9 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
IF (.NOT. ASSOCIATED(ri_xas_basis)) THEN
CALL cp_warn(__LOCATION__, "Automatic Generation of RI_XAS basis. "// &
"This is experimental code.")
!Because might need nco(l) with high l<=max_orbl in create_ri_aux_basis_set
CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto)
CALL init_orbital_pointers(2*maxlgto)
CALL create_ri_aux_basis_set(ri_xas_basis, qs_kind, dft_control%auto_basis_ri_xas)
CALL add_basis_set_to_container(qs_kind%basis_sets, ri_xas_basis, "RI_XAS")
END IF
@@ -24,7 +24,7 @@ MODULE qs_o3c_methods
USE gamma, ONLY: init_md_ftable
USE hfx_libint_wrapper, ONLY: cp_libint_t, cp_libint_cleanup_3eri, &
cp_libint_init_3eri, cp_libint_set_contrdepth
USE input_constants, ONLY: do_hfx_potential_truncated
USE input_constants, ONLY: do_potential_truncated
USE kinds, ONLY: dp
USE libint_3center, ONLY: eri_3center
USE orbital_pointers, ONLY: ncoset
@@ -404,7 +404,7 @@ SUBROUTINE calculate_o3c_libint_integrals(o3c, op, t_c_filename, para_env, r_cut
IF (PRESENT(r_cutoff)) my_r_cutoff = r_cutoff

!Init the truncated Coulomb operator
IF (op == do_hfx_potential_truncated) THEN
IF (op == do_potential_truncated) THEN
CPASSERT(PRESENT(t_c_filename))
CPASSERT(PRESENT(para_env))

@@ -35,7 +35,7 @@ MODULE xas_tdp_atom
dbcsr_distribution_new, dbcsr_distribution_release, &
dbcsr_complete_redistribute
USE distribution_2d_types, ONLY: distribution_2d_type, distribution_2d_release
USE input_constants, ONLY: do_hfx_potential_id
USE input_constants, ONLY: do_potential_id
USE input_section_types, ONLY: section_vals_type, &
section_vals_get_subs_vals
USE kinds, ONLY: dp, default_string_length
@@ -193,27 +193,30 @@ SUBROUTINE init_xas_atom_grid_harmo(xas_atom_env, grid_info, do_xc, qs_env)

TYPE(xas_atom_env_type), POINTER :: xas_atom_env
CHARACTER(len=default_string_length), &
DIMENSION(:, :), POINTER :: grid_info
LOGICAL, INTENT(IN) :: do_xc
TYPE(qs_environment_type), POINTER :: qs_env
DIMENSION(:,:), POINTER :: grid_info
TYPE(qs_environment_type), POINTER :: qs_env
LOGICAL, INTENT(IN) :: do_xc

CHARACTER(len=*), PARAMETER :: routineN = 'init_xas_atom_grid_harmo', &
routineP = moduleN//':'//routineN
CHARACTER(len=*), PARAMETER :: routineN = "init_xas_atom_grid_harmo", &
routineP = moduleN//":"//routineN

CHARACTER(len=2) :: symbol
INTEGER :: igrid, ikind, il, iso, iso1, iso2, l1, l1l2, l2, la, lc1, lc2, lcleb, ll, llmax, &
lp, m1, m2, max_s_harm, max_s_set, maxl, maxlgto, maxs, mm, mp, na, nr, quadrature, stat
REAL(dp) :: kind_radius
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rga
REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
TYPE(dft_control_type), POINTER :: dft_control
TYPE(grid_atom_type), POINTER :: grid_atom
TYPE(gto_basis_set_type), POINTER :: tmp_basis
TYPE(harmonics_atom_type), POINTER :: harmonics
TYPE(qs_control_type), POINTER :: qs_control
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
REAL(dp) :: kind_radius
REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rga
INTEGER :: llmax, maxlgto, max_s_harm, max_s_set, &
lcleb, lc1, iso1, l1, m1, lc2, iso2, &
l2, m2, l1l2, mp, mm, lp, il, iso, &
ikind, quadrature, ll, na, la, nr, &
maxs, maxl, igrid, stat
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(gto_basis_set_type), POINTER :: tmp_basis
TYPE(grid_atom_type), POINTER :: grid_atom
TYPE(harmonics_atom_type), POINTER :: harmonics
TYPE(qs_control_type), POINTER :: qs_control
TYPE(dft_control_type), POINTER :: dft_control
CHARACTER(LEN=default_string_length) :: kind_name

NULLIFY (my_CG, qs_kind_set, tmp_basis, grid_atom, harmonics, qs_control, dft_control)
NULLIFY(my_CG, qs_kind_set, tmp_basis, grid_atom, harmonics, qs_control, dft_control)

! Initialization of some integer for the CG coeff generation
CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, dft_control=dft_control)
@@ -301,12 +304,12 @@ SUBROUTINE init_xas_atom_grid_harmo(xas_atom_env, grid_info, do_xc, qs_env)
harmonics => xas_atom_env%harmonics_atom_set(ikind)%harmonics_atom

! Initialize some integers
CALL get_qs_kind(qs_kind_set(ikind), ngrid_rad=nr, ngrid_ang=na, element_symbol=symbol)
CALL get_qs_kind(qs_kind_set(ikind), ngrid_rad=nr, ngrid_ang=na, name=kind_name)

!take the grid dimension given as input, if none, take the GAPW ones above
IF (ANY(grid_info == symbol)) THEN
DO igrid = 1, SIZE(grid_info, 1)
IF (grid_info(igrid, 1) == symbol) THEN
IF (ANY(grid_info == kind_name)) THEN
DO igrid = 1, SIZE(grid_info,1)
IF (grid_info(igrid,1) == kind_name) THEN

!hack to convert string into integer
READ (grid_info(igrid, 2), *, iostat=stat) na
@@ -316,6 +319,9 @@ SUBROUTINE init_xas_atom_grid_harmo(xas_atom_env, grid_info, do_xc, qs_env)
EXIT
END IF
END DO
ELSE IF (do_xc .AND. ANY(xas_atom_env%excited_kinds == ikind)) THEN
!need good integration grids for the xc kernel, but taking the default GAPW grid
CPWARN("The default (and possibly small) GAPW grid is beeing used for one excited KIND")
END IF

ll = get_number_of_lebedev_grid(n=na)
@@ -925,15 +931,15 @@ SUBROUTINE calculate_density_coeffs(xas_atom_env, qs_env)

!Get the optimal distribution_2d for the overlap integrals (abc), centers c on excited atoms
!and their neighbors defined by RI_RADIUS
CALL build_xas_tdp_3c_nl(ac_list, basis_set_orb, basis_set_ri, do_hfx_potential_id, &
CALL build_xas_tdp_3c_nl(ac_list, basis_set_orb, basis_set_ri, do_potential_id, &
qs_env, excited_atoms=all_ri_atoms)
CALL get_opt_3c_dist2d(opt_dist2d, sab_orb, ac_list, basis_set_orb, basis_set_orb, &
basis_set_ri, qs_env)
CALL release_sab(ac_list)

!get the ab and ac neighbor lists based on the optimized distribution
CALL build_xas_tdp_ovlp_nl(ab_list, basis_set_orb, basis_set_orb, qs_env, ext_dist2d=opt_dist2d)
CALL build_xas_tdp_3c_nl(ac_list, basis_set_orb, basis_set_ri, do_hfx_potential_id, &
CALL build_xas_tdp_3c_nl(ac_list, basis_set_orb, basis_set_ri, do_potential_id, &
qs_env, excited_atoms=all_ri_atoms, ext_dist2d=opt_dist2d)

!get the AO density matrix in the optimized distribution for compatibility
@@ -984,8 +990,8 @@ SUBROUTINE calculate_density_coeffs(xas_atom_env, qs_env)
CALL o3c_iterator_create(o3c, o3c_iterator, nthread=nthread)

!$OMP PARALLEL DEFAULT(NONE)&
!$OMP SHARED (o3c_iterator,xas_atom_env,nsgf,nspins,factor,nnb,neighbors,iex,idx_to_nb,ri_sinv,inb)&
!$OMP PRIVATE (mepos,katom,tvec,work,jnb,ri_at,found,foundt,pblock,pblockt)
!$OMP SHARED (o3c_iterator,xas_atom_env,nsgf,nspins,factor,nnb,neighbors,iex,idx_to_nb,ri_sinv)&
!$OMP PRIVATE (mepos,katom,tvec,work,inb,jnb,ri_at,found,foundt,pblock,pblockt)

mepos = 0
!$ mepos = omp_get_thread_num()

0 comments on commit f1cd4a3

Please sign in to comment.
You can’t perform that action at this time.