XAS_TDP| Added EPS_PGF_XAS keyword to define basis function radii

independently from ground state calculation.

src/auto_basis.F

@@ -23,6 +23,7 @@ MODULE auto_basis
                                               gamma1,&
                                               pi,&
                                               rootpi
+   USE orbital_pointers,                ONLY: init_orbital_pointers
    USE periodic_table,                  ONLY: get_ptable_info
    USE powell,                          ONLY: opt_state_type,&
                                               powell_optimize
@@ -84,6 +85,9 @@ SUBROUTINE create_ri_aux_basis_set(ri_aux_basis_set, qs_kind, basis_cntrl)
       CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, basis_type="ORB")
       IF (ASSOCIATED(orb_basis_set)) THEN
          CALL get_basis_keyfigures(orb_basis_set, lmax, zmin, zmax, zeff)
+         !Note: RI basis coud require lmax up to 2*orb_lmax. This ensures that all orbital pointers
+         !      are properly initialized before building the basis
+         CALL init_orbital_pointers(2*lmax)
          CALL get_basis_products(lmax, zmin, zmax, zeff, pmin, pmax, peff)
          CALL get_qs_kind(qs_kind, zeff=zval, elec_conf=econf, element_symbol=element_symbol)
          CALL get_ptable_info(element_symbol, ielement=z)

src/input_cp2k_dft.F

@@ -7888,7 +7888,7 @@ SUBROUTINE create_xas_tdp_section(section)
                           description="XAS simulations using linear-response TDDFT. 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=12, n_subsections=4, repeats=.FALSE.)
+                          n_keywords=13, n_subsections=4, repeats=.FALSE.)
 
       NULLIFY (keyword, subsection, print_key)
 
@@ -7945,6 +7945,19 @@ SUBROUTINE create_xas_tdp_section(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create(keyword, __LOCATION__, name="EPS_PGF_XAS", &
+                          variants=s2a("EPS_PGF", "EPS_PGF_XAS_TDP"), &
+                          description="The threshold used to determine the spacial extent of all "// &
+                          "primitive Gaussian functions used for the construction "// &
+                          "of neighbor lists in the XAS_TDP method. "// &
+                          "By default, takes the value of QS%EPS_PGF_ORB. Useful if "// &
+                          "the former value is tiny due to possible ground state HFX "// &
+                          "contributions.", &
+                          usage="EPS_PGS_XAS {real}", &
+                          type_of_var=real_t)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
       CALL keyword_create(keyword, __LOCATION__, name="DIPOLE_FORM", &
                           variants=(/"DIP_FORM"/), &
                           description="Type of integral to get the oscillator strengths "// &
@@ -8272,7 +8285,7 @@ SUBROUTINE create_xas_tdp_section(section)
                           "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=0.0_dp, &
+                          default_r_val=1.0E-10_dp, &
                           repeats=.FALSE.)
       CALL section_add_keyword(subsubsection, keyword)
       CALL keyword_release(keyword)

src/qs_environment.F

@@ -803,9 +803,7 @@ 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)
+               ! Generate a default basis
                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

src/qs_o3c_methods.F

@@ -344,50 +344,63 @@ SUBROUTINE calculate_o3c_libint_integrals(o3c, op, t_c_filename, para_env, r_cut
       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_o3c_libint_integrals', &
          routineP = moduleN//':'//routineN
 
-      INTEGER :: egfa, egfb, egfc, handle, ibasis, ikind, imax, iset, jkind, jset, kkind, kset, &
-         m_max, maxli, maxlj, maxlk, mepos, ncoa, ncob, ncoc, ni, nj, nk, nseta, nsetb, nsetc, &
-         nspin, nthread, sgfa, sgfb, sgfc, unit_id
+      INTEGER :: egfa, egfb, egfc, handle, i, ibasis, ikind, ilist, imax, iset, jkind, jset, &
+         kkind, kset, m_max, max_nset, maxli, maxlj, maxlk, mepos, nbasis, ncoa, ncob, ncoc, ni, &
+         nj, nk, nseta, nsetb, nsetc, nspin, nthread, sgfa, sgfb, sgfc, unit_id
       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, lc_max, &
                                                             lc_min, npgfa, npgfb, npgfc, nsgfa, &
                                                             nsgfb, nsgfc
       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb, first_sgfc
       LOGICAL                                            :: do_screen
-      REAL(dp)                                           :: my_eps_screen, my_omega, my_r_cutoff
+      REAL(dp)                                           :: max_val, my_eps_screen, my_omega, &
+                                                            my_r_cutoff
+      REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: max_contr, max_contra, max_contrb, &
+                                                            max_contrc
       REAL(dp), ALLOCATABLE, DIMENSION(:, :, :)          :: sabc
       REAL(dp), DIMENSION(3)                             :: ri, rij, rik, rj, rk
       REAL(dp), DIMENSION(:, :), POINTER                 :: sphi_a, sphi_b, sphi_c, tvec, zeta, &
                                                             zetb, zetc
       REAL(dp), DIMENSION(:, :, :), POINTER              :: iabc
       TYPE(cp_libint_t)                                  :: lib
-      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
+      TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list, basis_set_list_a, &
+                                                            basis_set_list_b, basis_set_list_c
+      TYPE(gto_basis_set_type), POINTER                  :: basis_set, basis_set_a, basis_set_b, &
+                                                            basis_set_c
       TYPE(o3c_iterator_type)                            :: o3c_iterator
 
       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)
+      NULLIFY (basis_set, basis_set_list)
 
       CALL timeset(routineN, handle)
 
       CALL get_o3c_container(o3c, nspin=nspin, basis_set_list_a=basis_set_list_a, &
                              basis_set_list_b=basis_set_list_b, basis_set_list_c=basis_set_list_c)
 
-      !Need the max l for each basis for libint
+      !Need the max l for each basis for libint (and overall max #of sets for screening)
+      nbasis = SIZE(basis_set_list_a)
+      max_nset = 0
       maxli = 0
-      DO ibasis = 1, SIZE(basis_set_list_a)
-         CALL get_gto_basis_set(gto_basis_set=basis_set_list_a(ibasis)%gto_basis_set, maxl=imax)
+      DO ibasis = 1, nbasis
+         CALL get_gto_basis_set(gto_basis_set=basis_set_list_a(ibasis)%gto_basis_set, &
+                                maxl=imax, nset=iset)
          maxli = MAX(maxli, imax)
+         max_nset = MAX(max_nset, iset)
       END DO
       maxlj = 0
-      DO ibasis = 1, SIZE(basis_set_list_b)
-         CALL get_gto_basis_set(gto_basis_set=basis_set_list_b(ibasis)%gto_basis_set, maxl=imax)
+      DO ibasis = 1, nbasis
+         CALL get_gto_basis_set(gto_basis_set=basis_set_list_b(ibasis)%gto_basis_set, &
+                                maxl=imax, nset=iset)
          maxlj = MAX(maxlj, imax)
+         max_nset = MAX(max_nset, iset)
       END DO
       maxlk = 0
-      DO ibasis = 1, SIZE(basis_set_list_c)
-         CALL get_gto_basis_set(gto_basis_set=basis_set_list_c(ibasis)%gto_basis_set, maxl=imax)
+      DO ibasis = 1, nbasis
+         CALL get_gto_basis_set(gto_basis_set=basis_set_list_c(ibasis)%gto_basis_set, &
+                                maxl=imax, nset=iset)
          maxlk = MAX(maxlk, imax)
+         max_nset = MAX(max_nset, iset)
       END DO
       m_max = maxli+maxlj+maxlk
 
@@ -429,13 +442,14 @@ SUBROUTINE calculate_o3c_libint_integrals(o3c, op, t_c_filename, para_env, r_cut
       CALL o3c_iterator_create(o3c, o3c_iterator, nthread=nthread)
 
 !$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,do_screen,my_eps_screen) &
+!$OMP SHARED (nthread,o3c_iterator,nspin,basis_set_list_a, basis_set_list_b,basis_set_list_c,nbasis,&
+!$OMP         maxli,maxlj,maxlk,my_r_cutoff,my_omega,ncoset,op,do_screen,my_eps_screen,max_nset) &
 !$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,&
 !$OMP          first_sgfc,lc_max,lc_min,npgfc,nsetc,nsgfc,zetc,kset,nk,rk,ncoc,sgfc,egfc,sphi_c,&
-!$OMP          sabc,lib)
+!$OMP          sabc,lib,max_contra,max_contrb,max_contrc,ibasis,i,basis_set,basis_set_list,ilist,&
+!$OMP          max_contr,max_val)
 
       mepos = 0
 !$    mepos = omp_get_thread_num()
@@ -444,6 +458,45 @@ SUBROUTINE calculate_o3c_libint_integrals(o3c, op, t_c_filename, para_env, r_cut
       CALL cp_libint_init_3eri(lib, MAX(maxli, maxlj, maxlk))
       CALL cp_libint_set_contrdepth(lib, 1)
 
+      !get the max_contraction values before we loop, on each thread => least amount of computation
+      !and false sharing
+      IF (do_screen) THEN
+
+         !Allocate max_contraction arrays such that we have a specific value for each set/kind
+         ALLOCATE (max_contr(max_nset, nbasis), max_contra(max_nset, nbasis), &
+                   max_contrb(max_nset, nbasis), max_contrc(max_nset, nbasis))
+
+         !Not the most elegent, but better than copying 3 times the same
+         DO ilist = 1, 3
+
+            IF (ilist == 1) basis_set_list => basis_set_list_a
+            IF (ilist == 2) basis_set_list => basis_set_list_b
+            IF (ilist == 3) basis_set_list => basis_set_list_c
+
+            max_contr = 0.0_dp
+
+            DO ibasis = 1, nbasis
+               basis_set => basis_set_list(ibasis)%gto_basis_set
+
+               DO iset = 1, basis_set%nset
+
+                  ncoa = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
+                  sgfa = basis_set%first_sgf(1, iset)
+                  egfa = sgfa+basis_set%nsgf_set(iset)-1
+
+                  max_contr(iset, ibasis) = &
+                     MAXVAL((/(SUM(ABS(basis_set%sphi(1:ncoa, i))), i=sgfa, egfa)/))
+
+               END DO !iset
+            END DO !ibasis
+
+            IF (ilist == 1) max_contra(:, :) = max_contr(:, :)
+            IF (ilist == 2) max_contrb(:, :) = max_contr(:, :)
+            IF (ilist == 3) max_contrc(:, :) = max_contr(:, :)
+         END DO !ilist
+         DEALLOCATE (max_contr)
+      END IF !do_screen
+
       DO WHILE (o3c_iterate(o3c_iterator, mepos=mepos) == 0)
 
          CALL get_o3c_iterator_info(o3c_iterator, mepos=mepos, ikind=ikind, jkind=jkind, &
@@ -520,7 +573,9 @@ SUBROUTINE calculate_o3c_libint_integrals(o3c, op, t_c_filename, para_env, r_cut
                                    rk, lib, op, omega=my_omega, r_cutoff=my_r_cutoff)
 
                   IF (do_screen) THEN
-                     IF (MAXVAL(ABS(sabc)) < my_eps_screen) THEN
+                     max_val = MAXVAL(ABS(sabc))*max_contra(iset, ikind)*max_contrb(jset, jkind) &
+                               *max_contrc(kset, kkind)
+                     IF (max_val < my_eps_screen) THEN
                         DEALLOCATE (sabc)
                         CYCLE
                      END IF

src/xas_tdp_atom.F

@@ -874,15 +874,15 @@ SUBROUTINE calculate_density_coeffs(xas_atom_env, qs_env)
       TYPE(distribution_2d_type), POINTER                :: opt_dist2d
       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_orb, basis_set_ri
       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
-         POINTER                                         :: ab_list, ac_list, sab_orb, sab_ri
+         POINTER                                         :: ab_list, ac_list, sab_ri
       TYPE(o3c_container_type), POINTER                  :: o3c
       TYPE(o3c_iterator_type)                            :: o3c_iterator
       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
       TYPE(qs_rho_type), POINTER                         :: rho
 
       NULLIFY (qs_kind_set, basis_set_ri, basis_set_orb, ac_list, rho, rho_ao, sab_ri, ri_mats)
-      NULLIFY (o3c, tvec, particle_set, para_env, all_ri_atoms, sab_orb, overlap)
+      NULLIFY (o3c, tvec, particle_set, para_env, all_ri_atoms, overlap)
       NULLIFY (blacs_env, pblock, pblockt, ab_list, opt_dist2d, rho_redist)
 
       !Idea: We don't do a full RI here as it would be too expensive in many ways (inversion of a
@@ -897,7 +897,7 @@ SUBROUTINE calculate_density_coeffs(xas_atom_env, qs_env)
 
       CALL timeset(routineN, handle)
 
-      CALL get_qs_env(qs_env, nkind=nkind, sab_orb=sab_orb, qs_kind_set=qs_kind_set, rho=rho, &
+      CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, rho=rho, &
                       natom=natom, particle_set=particle_set, &
                       para_env=para_env, blacs_env=blacs_env)
       nspins = xas_atom_env%nspins
@@ -950,10 +950,12 @@ 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_ovlp_nl(ab_list, basis_set_orb, basis_set_orb, qs_env)
       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, &
+      CALL get_opt_3c_dist2d(opt_dist2d, ab_list, ac_list, basis_set_orb, basis_set_orb, &
                              basis_set_ri, qs_env)
+      CALL release_sab(ab_list)
       CALL release_sab(ac_list)
 
       !get the ab and ac neighbor lists based on the optimized distribution

src/xas_tdp_kernel.F

@@ -618,8 +618,8 @@ SUBROUTINE kernel_exchange(ex_ker, donor_state, xas_tdp_env, xas_tdp_control, qs
       blk_size => donor_state%blk_size
 
       !compute the off-diag spin-conserving only if not TDA and anything that is spin-conserving
-      do_off_sc = (.NOT. xas_tdp_control%tamm_dancoff) .AND. (xas_tdp_control%do_spin_cons &
-                                                              .OR. xas_tdp_control%do_singlet .OR. xas_tdp_control%do_triplet)
+      do_off_sc = (.NOT. xas_tdp_control%tamm_dancoff) .AND. &
+                  (xas_tdp_control%do_spin_cons .OR. xas_tdp_control%do_singlet .OR. xas_tdp_control%do_triplet)
 
       ! Need the once contracted integrals (aI|P)
       CALL contract_o3c_int(contr1_int, "EXCHANGE", donor_state, xas_tdp_env, xas_tdp_control, qs_env)
@@ -1072,9 +1072,9 @@ SUBROUTINE change_dist(mat_in, mat_out, para_env)
       num_pe = para_env%num_pe
 
       !Allocate a pointer array which will point on the block => we won't need to send more than
-      !nblk**2/(2*num_pe) message per processor (only look at upper diagonal, blocks well distributed)
-      ALLOCATE (send_buff(nblk**2/(2*num_pe)+1), recv_buff(nblk**2/(2*num_pe)+1))
-      ALLOCATE (send_req(nblk**2/(2*num_pe)+1), recv_req(nblk**2/(2*num_pe)+1))
+      !nblk**2/num_pe message per processor (assumes blocks well distributed)
+      ALLOCATE (send_buff(nblk**2/num_pe+1), recv_buff(nblk**2/num_pe+1))
+      ALLOCATE (send_req(nblk**2/num_pe+1), recv_req(nblk**2/num_pe+1))
       is = 0; ir = 0
 
       !Iterate over input matrix