Skip to content
Permalink
Browse files

XAS_TDP| Added option of doing stendard DFT (i.e. without kernel)

  • Loading branch information...
abussy authored and pseewald committed Jul 1, 2019
1 parent e2f7194 commit 9319f59c14301744b95f7708028e4ae7fe107059
Showing with 3,015 additions and 2,624 deletions.
  1. +16 −5 src/aobasis/ai_contraction_sphi.F
  2. +8 −40 src/distribution_methods.F
  3. +13 −2 src/input_cp2k_dft.F
  4. +1 −1 src/qs_ks_types.F
  5. +21 −8 src/qs_o3c_methods.F
  6. +194 −229 src/xas_tdp_atom.F
  7. +424 −419 src/xas_tdp_kernel.F
  8. +723 −727 src/xas_tdp_methods.F
  9. +155 −170 src/xas_tdp_types.F
  10. +1,460 −1,023 src/xas_tdp_utils.F
@@ -7,8 +7,7 @@
!> \brief Contraction of integrals over primitive Cartesian Gaussians based on the contraction
!> matrix sphi which is part of the gto_basis_set_type
!> \par History
!> -Removed timings (routines might be called 100,000s => cost of timeset might actually become
!> high, also not much info to be gained by these timings + cheap anyways) A. Bussy (05.2019)
!> none
!> \author Dorothea Golze (05.2016)
! **************************************************************************************************
MODULE ai_contraction_sphi
@@ -45,9 +44,11 @@ SUBROUTINE ab_contract(abint, sab, sphi_a, sphi_b, ncoa, ncob, nsgfa, nsgfb)

CHARACTER(LEN=*), PARAMETER :: routineN = 'ab_contract', routineP = moduleN//':'//routineN

INTEGER :: m1, m2, msphia, msphib, nn
INTEGER :: handle, m1, m2, msphia, msphib, nn
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cpp

CALL timeset(routineN, handle)

msphia = SIZE(sphi_a, 1)
msphib = SIZE(sphi_b, 1)

@@ -64,6 +65,8 @@ SUBROUTINE ab_contract(abint, sab, sphi_a, sphi_b, ncoa, ncob, nsgfa, nsgfb)

DEALLOCATE (cpp)

CALL timestop(handle)

END SUBROUTINE ab_contract

! **************************************************************************************************
@@ -90,10 +93,12 @@ SUBROUTINE abc_contract(abcint, sabc, sphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc,

CHARACTER(LEN=*), PARAMETER :: routineN = 'abc_contract', routineP = moduleN//':'//routineN

INTEGER :: i, m1, m2, m3, msphia, msphib, &
INTEGER :: handle, i, m1, m2, m3, msphia, msphib, &
msphic
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: cpc, cpp

CALL timeset(routineN, handle)

CPASSERT(SIZE(abcint, 1) == nsgfa)
CPASSERT(SIZE(abcint, 2) == nsgfb)

@@ -117,6 +122,8 @@ SUBROUTINE abc_contract(abcint, sabc, sphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc,

DEALLOCATE (cpp, cpc)

CALL timestop(handle)

END SUBROUTINE abc_contract

! **************************************************************************************************
@@ -149,12 +156,14 @@ SUBROUTINE abcd_contract(abcdint, sabcd, sphi_a, sphi_b, sphi_c, sphi_d, ncoa, n

CHARACTER(LEN=*), PARAMETER :: routineN = 'abcd_contract', routineP = moduleN//':'//routineN

INTEGER :: isgfc, isgfd, m1, m2, m3, m4, &
INTEGER :: handle, isgfc, isgfd, m1, m2, m3, m4, &
msphia, msphib, msphic, msphid
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: temp_cccc, work_cpcc
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: temp_cpcc, work_cppc
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: cpcc, cppc, cppp

CALL timeset(routineN, handle)

msphia = SIZE(sphi_a, 1)
msphib = SIZE(sphi_b, 1)
msphic = SIZE(sphi_c, 1)
@@ -197,6 +206,8 @@ SUBROUTINE abcd_contract(abcdint, sabcd, sphi_a, sphi_b, sphi_c, sphi_d, ncoa, n
DEALLOCATE (cpcc, cppc, cppp)
DEALLOCATE (work_cpcc, work_cppc, temp_cpcc, temp_cccc)

CALL timestop(handle)

END SUBROUTINE abcd_contract

END MODULE ai_contraction_sphi
@@ -81,8 +81,7 @@ MODULE distribution_methods
! *** Public subroutines ***

PUBLIC :: distribute_molecules_1d, &
distribute_molecules_2d, &
make_basic_spatial_distribution
distribute_molecules_2d

CONTAINS

@@ -906,27 +905,23 @@ END SUBROUTINE make_basic_distribution
!> \param row_distribution ...
!> \param npcols ...
!> \param col_distribution ...
!> \param maximize_spread tries to put atoms that are close together on different procs
!> \par History
!> - Created 2010-11-11 Joost VandeVondele
!> - Added the maximize_spread argument 05.2019 A. Bussy
! **************************************************************************************************
SUBROUTINE make_basic_spatial_distribution(pbc_scaled_coords, costs, nprows, row_distribution, &
npcols, col_distribution, maximize_spread)
SUBROUTINE make_basic_spatial_distribution(pbc_scaled_coords, costs, &
nprows, row_distribution, npcols, col_distribution)
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: pbc_scaled_coords
INTEGER, DIMENSION(:), INTENT(IN) :: costs
INTEGER, INTENT(IN) :: nprows
INTEGER, DIMENSION(:), INTENT(OUT) :: row_distribution
INTEGER, INTENT(IN) :: npcols
INTEGER, DIMENSION(:), INTENT(OUT) :: col_distribution
LOGICAL, INTENT(IN), OPTIONAL :: maximize_spread

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

INTEGER :: handle, iatom, natoms, nbins, pgrid_gcd
INTEGER, ALLOCATABLE, DIMENSION(:) :: bin_costs, distribution
LOGICAL :: my_max_dist

CALL timeset(routineN, handle)

@@ -936,13 +931,7 @@ SUBROUTINE make_basic_spatial_distribution(pbc_scaled_coords, costs, nprows, row
ALLOCATE (bin_costs(nbins), distribution(natoms))
bin_costs = 0

my_max_dist = .FALSE.
IF (PRESENT(maximize_spread)) THEN
my_max_dist = maximize_spread
END IF

CALL spatial_recurse(pbc_scaled_coords, costs, (/(iatom, iatom=1, natoms)/), bin_costs, &
distribution, 0, my_max_dist)
CALL spatial_recurse(pbc_scaled_coords, costs, (/(iatom, iatom=1, natoms)/), bin_costs, distribution, 0)

! WRITE(*, *) "Final bin costs: ", bin_costs

@@ -966,19 +955,14 @@ END SUBROUTINE make_basic_spatial_distribution
!> \param bin_costs ...
!> \param distribution ...
!> \param level ...
!> \param maximize_spread ...
! **************************************************************************************************
RECURSIVE SUBROUTINE spatial_recurse(pbc_scaled_coords, costs, indices, bin_costs, distribution,&
level, maximize_spread)
RECURSIVE SUBROUTINE spatial_recurse(pbc_scaled_coords, costs, indices, bin_costs, distribution, level)
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: pbc_scaled_coords
INTEGER, DIMENSION(:), INTENT(IN) :: costs, indices
INTEGER, DIMENSION(:), INTENT(INOUT) :: bin_costs, distribution
INTEGER, INTENT(IN) :: level
LOGICAL, INTENT(IN) :: maximize_spread

INTEGER :: iatom, ibin, natoms, nbins, nhalf, i, &
nquarter, itmp
REAL(KIND=dp) :: rtmp
INTEGER :: iatom, ibin, natoms, nbins, nhalf
INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_costs_sorted, atom_permutation, &
bin_costs_sorted, permutation
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coord
@@ -1012,26 +996,10 @@ RECURSIVE SUBROUTINE spatial_recurse(pbc_scaled_coords, costs, indices, bin_cost
ALLOCATE (coord(natoms), permutation(natoms))
coord(:) = pbc_scaled_coords(MOD(level, 3)+1, :)
CALL sort(coord, natoms, permutation)

IF (maximize_spread) THEN
! reorder the coord (and permutation) such that atoms that are far away are put together
nquarter = (natoms+1)/2
DO i = 1, nquarter
rtmp = coord(i)
coord(i) = coord(natoms+1-i)
coord(natoms+1-i) = rtmp
itmp = permutation(i)
permutation(i) = permutation(natoms+1-i)
permutation(natoms+1-i) = itmp
END DO
END IF

CALL spatial_recurse(pbc_scaled_coords(:, permutation(1:nhalf)), costs(permutation(1:nhalf)), &
indices(permutation(1:nhalf)), bin_costs, distribution, level+1, &
maximize_spread)
indices(permutation(1:nhalf)), bin_costs, distribution, level+1)
CALL spatial_recurse(pbc_scaled_coords(:, permutation(nhalf+1:)), costs(permutation(nhalf+1:)), &
indices(permutation(nhalf+1:)), bin_costs, distribution, level+1, &
maximize_spread)
indices(permutation(nhalf+1:)), bin_costs, distribution, level+1)
DEALLOCATE (coord, permutation)
ENDIF

@@ -8196,7 +8196,7 @@ SUBROUTINE create_xas_tdp_section(section)
CALL section_create(subsubsection, __LOCATION__, name="EXACT_EXCHANGE", &
description="Whther exact-exchange should be added to the kernel and " //&
"if so, with which scaling.", &
n_keywords=6, &
n_keywords=7, &
n_subsections=0, &
repeats=.FALSE.)

@@ -8250,12 +8250,23 @@ SUBROUTINE create_xas_tdp_section(section)

CALL keyword_create(keyword, __LOCATION__, name="EPS_RANGE", &
description="The threshold to determine the effective range of the short range "//&
"operator: erfc(omega*eff_range)/eff_range = EPS_EANGE", &
"operator: erfc(omega*eff_range)/eff_range = EPS_RANGE", &
usage="EPS_RANGE = {double}", &
default_r_val=1.0E-6_dp, &
repeats=.FALSE.)
CALL section_add_keyword(subsubsection, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="EPS_SCREENING", &
variants=s2a("EPS_SCREEN"), &
description="A threshold to determine which primitive 3-center integrals"//&
"are kept for contraction, as the latter operation can be "//&
"expensive (especially for large basis sets )." //&
"If |(ab|c)| < EPS_SCREENNING, it is discarded.", &
default_r_val=1.0E-10_dp, &
repeats=.FALSE.)
CALL section_add_keyword(subsubsection, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="SCALE", &
description="Scaling of the exact exchange contribution.", &
@@ -76,7 +76,7 @@ MODULE qs_ks_types
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_types'

PUBLIC :: qs_ks_env_type, qs_ks_env_create, qs_ks_did_change
PUBLIC :: qs_ks_release, qs_ks_retain, get_ks_env, set_ks_env
PUBLIC :: qs_ks_release, qs_ks_retain, get_ks_env, set_ks_env, release_sab

! **************************************************************************************************
!> \brief calculation environment to calculate the ks matrix,
@@ -323,17 +323,19 @@ END SUBROUTINE calculate_o3c_integrals
!> \param t_c_filename the name of the file with the truncated Coulomb data
!> \param r_cutoff the cutoff radius for the truncated coulomb operator
!> \param omega the range parameter for the erfc short range operator
!> \param eps_screen the screening threshold for sicarding integrals before contraction
!> \note The static initialization of the libint library needs to be done before hand
!> In case the truncated coulomb operator is used, the potential parameter file must be read
!> in advance too
! **************************************************************************************************
SUBROUTINE calculate_o3c_libint_integrals(o3c, op, t_c_filename, para_env, r_cutoff, omega)
SUBROUTINE calculate_o3c_libint_integrals(o3c, op, t_c_filename, para_env, r_cutoff, omega, &
eps_screen)

TYPE(o3c_container_type), POINTER :: o3c
INTEGER, INTENT(IN) :: op
CHARACTER(len=*), INTENT(IN), OPTIONAL :: t_c_filename
TYPE(cp_para_env_type), POINTER, OPTIONAL :: para_env
REAL(dp), INTENT(IN), OPTIONAL :: r_cutoff, omega
REAL(dp), INTENT(IN), OPTIONAL :: r_cutoff, omega, eps_screen

CHARACTER(len=*), PARAMETER :: routineN = "calculate_o3c_libint_integrals", &
routineP = moduleN//":"//routineN
@@ -344,6 +346,7 @@ SUBROUTINE calculate_o3c_libint_integrals(o3c, op, t_c_filename, para_env, r_cut
ncob, ncoc, sgfa, sgfb, sgfc, egfa, &
egfb, egfc, maxli, maxlj, maxlk, imax, &
ibasis, unit_id, m_max
LOGICAL :: do_screen
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b, &
basis_set_list_c
TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b, basis_set_c
@@ -358,13 +361,12 @@ SUBROUTINE calculate_o3c_libint_integrals(o3c, op, t_c_filename, para_env, r_cut
nsgfb, nsgfc
INTEGER, DIMENSION(:,:), POINTER :: first_sgfa, first_sgfb, first_sgfc
TYPE(cp_libint_t) :: lib
REAL(dp) :: my_omega, my_r_cutoff
REAL(dp) :: my_omega, my_r_cutoff, my_eps_screen

NULLIFY(basis_set_list_a, basis_set_list_b, basis_set_list_c, basis_set_a, basis_set_b)
NULLIFY(basis_set_c, iabc, tvec, first_sgfa, first_sgfb, first_sgfc, la_max, la_min, lb_max)
NULLIFY(lb_min, lc_max, lc_min, npgfa, npgfb, npgfc, nsgfa, nsgfb, nsgfc)


CALL timeset(routineN, handle)

CALL get_o3c_container(o3c, nspin=nspin, basis_set_list_a=basis_set_list_a, &
@@ -388,6 +390,13 @@ SUBROUTINE calculate_o3c_libint_integrals(o3c, op, t_c_filename, para_env, r_cut
END DO
m_max = maxli+maxlj+maxlk

!Screening
do_screen = .FALSE.
IF (PRESENT(eps_screen)) THEN
do_screen = .TRUE.
my_eps_screen = eps_screen
END IF

!Short range parameters
my_omega = 0.0_dp
IF (PRESENT(omega)) my_omega = omega
@@ -420,7 +429,7 @@ SUBROUTINE calculate_o3c_libint_integrals(o3c, op, t_c_filename, para_env, r_cut

!$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED (nthread,o3c_iterator,nspin,basis_set_list_a, basis_set_list_b,basis_set_list_c, &
!$OMP maxli,maxlj,maxlk,my_r_cutoff,my_omega,ncoset,op) &
!$OMP maxli,maxlj,maxlk,my_r_cutoff,my_omega,ncoset,op,do_screen,my_eps_screen) &
!$OMP PRIVATE (basis_set_a,basis_set_b,basis_set_c,mepos,ikind,jkind,kkind,iabc,tvec,rij,rik, &
!$OMP first_sgfa,la_max,la_min,npgfa,nseta,nsgfa,zeta,iset,ni,ri,ncoa,sgfa,egfa,sphi_a,&
!$OMP first_sgfb,lb_max,lb_min,npgfb,nsetb,nsgfb,zetb,jset,nj,rj,ncob,sgfb,egfb,sphi_b,&
@@ -439,9 +448,6 @@ SUBROUTINE calculate_o3c_libint_integrals(o3c, op, t_c_filename, para_env, r_cut
CALL get_o3c_iterator_info(o3c_iterator, mepos=mepos, ikind=ikind, jkind=jkind, &
kkind=kkind, rij=rij, rik=rik, integral=iabc, tvec=tvec)

CPASSERT(.NOT. ASSOCIATED(iabc))
CPASSERT(.NOT. ASSOCIATED(tvec))

!basis
basis_set_a => basis_set_list_a(ikind)%gto_basis_set
basis_set_b => basis_set_list_b(jkind)%gto_basis_set
@@ -512,6 +518,13 @@ SUBROUTINE calculate_o3c_libint_integrals(o3c, op, t_c_filename, para_env, r_cut
rj, lc_min(kset), lc_max(kset), npgfc(kset), zetc(:,kset), &
rk, lib, op, omega=my_omega, r_cutoff=my_r_cutoff)

IF (do_screen) THEN
IF (MAXVAL(ABS(sabc)) < my_eps_screen) THEN
DEALLOCATE(sabc)
CYCLE
END IF
END IF

CALL abc_contract(iabc(sgfa:egfa, sgfb:egfb, sgfc:egfc), sabc, &
sphi_a(:, sgfa:), sphi_b(:, sgfb:), sphi_c(:, sgfc:), &
ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfc(kset))

0 comments on commit 9319f59

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