low-scaling RPA/GW: sort basis sets with respect to exponent

This increases sparsity since low-scaling RPA/GW uses blocks smaller
than atomic.

This optimization can not be enabled automatically for RPA since this
would lead to inconsistencies (e.g. when restarting a RPA calculation
from a DFT wavefunction). Therefore a new keyword SORT_BASIS in DFT
input section was introduced to activate this optimization.

src/aobasis/basis_set_types.F

@@ -59,6 +59,11 @@ MODULE basis_set_types
 
    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'basis_set_types'
 
+   ! basis set sort criteria
+   INTEGER, PARAMETER, PUBLIC               :: basis_sort_default = 0, &
+                                               basis_sort_radius = 1, &
+                                               basis_sort_zet = 2
+
 ! **************************************************************************************************
    ! Define the Gaussian-type orbital basis set type
 
@@ -79,6 +84,7 @@ MODULE basis_set_types
       INTEGER, DIMENSION(:, :), POINTER           :: first_cgf, first_sgf, l, &
                                                      last_cgf, last_sgf, n
       REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc
+      INTEGER                                    :: sort = basis_sort_default
    END TYPE gto_basis_set_type
 
    TYPE gto_basis_set_p_type
@@ -119,6 +125,7 @@ MODULE basis_set_types
              create_primitive_basis_set, &
              combine_basis_sets, &
              set_gto_basis_set, &
+             sort_gto_basis_set, &
              write_gto_basis_set, &
              write_orb_basis_set
 
@@ -263,6 +270,7 @@ SUBROUTINE copy_gto_basis_set(basis_set_in, basis_set_out)
       bout%aliases = bin%aliases
       bout%kind_radius = bin%kind_radius
       bout%norm_type = bin%norm_type
+      bout%sort = bin%sort
       bout%nset = bin%nset
       bout%ncgf = bin%ncgf
       bout%nsgf = bin%nsgf
@@ -379,6 +387,7 @@ SUBROUTINE create_primitive_basis_set(basis_set, pbasis)
       pbasis%kind_radius = basis_set%kind_radius
       pbasis%short_kind_radius = basis_set%short_kind_radius
       pbasis%norm_type = basis_set%norm_type
+      pbasis%sort = basis_set%sort
       nset = lm + 1
       pbasis%nset = nset
       ALLOCATE (pbasis%lmax(nset), pbasis%lmin(nset), pbasis%npgf(nset), pbasis%nshell(nset))
@@ -2714,6 +2723,198 @@ FUNCTION srules(z, ne, n, l)
       srules = (REAL(z, dp) - s)/xns(nn)
    END FUNCTION srules
 
+! **************************************************************************************************
+!> \brief sort basis sets w.r.t. radius
+!> \param basis_set ...
+! **************************************************************************************************
+   SUBROUTINE sort_gto_basis_set(basis_set)
+      TYPE(gto_basis_set_type), POINTER                  :: basis_set
+
+      CHARACTER(LEN=12), DIMENSION(:), POINTER           :: cgf_symbol
+      CHARACTER(LEN=6), DIMENSION(:), POINTER            :: sgf_symbol
+      INTEGER :: ic, ic_max, icgf, icgf_new, icgf_old, ico, is, is_max, iset, isgf, isgf_new, &
+         isgf_old, ishell, lshell, maxco, maxpgf, maxshell, mm, nc, ncgf, ns, nset
+      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: sort_index
+      INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: icgf_set, isgf_set
+      INTEGER, DIMENSION(:), POINTER                     :: lx, ly, lz, m, npgf
+      REAL(dp), ALLOCATABLE, DIMENSION(:)                :: tmp
+      REAL(dp), DIMENSION(:), POINTER                    :: set_radius
+      REAL(dp), DIMENSION(:, :), POINTER                 :: zet
+      REAL(KIND=dp), DIMENSION(:), POINTER               :: norm_cgf
+      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: cphi, scon, sphi
+      TYPE(gto_basis_set_type), POINTER                  :: b
+
+      CPASSERT(ASSOCIATED(basis_set))
+
+      IF (basis_set%sort == basis_sort_default) RETURN
+
+      NULLIFY (set_radius, zet)
+
+      b => basis_set
+      CALL get_gto_basis_set(gto_basis_set=b, &
+                             nset=nset, &
+                             maxshell=maxshell, &
+                             maxpgf=maxpgf, &
+                             maxco=maxco, &
+                             ncgf=ncgf, &
+                             npgf=npgf, &
+                             set_radius=set_radius, &
+                             zet=zet)
+
+      ALLOCATE (sort_index(nset))
+      ALLOCATE (tmp(nset))
+      SELECT CASE (b%sort)
+      CASE (basis_sort_radius)
+         tmp(:) = set_radius(UBOUND(set_radius, 1):LBOUND(set_radius, 1):-1)
+      CASE (basis_sort_zet)
+         DO iset = 1, nset
+            tmp(iset) = MINVAL(basis_set%zet(:npgf(iset), iset))
+         ENDDO
+      CASE DEFAULT
+         CPABORT("Request basis sort criterion not implemented.")
+      END SELECT
+
+      CALL sort(tmp(1:nset), nset, sort_index)
+
+      ic_max = 0
+      is_max = 0
+      DO iset = 1, nset
+         ic = 0
+         is = 0
+         DO ishell = 1, b%nshell(iset)
+            DO ico = 1, nco(b%l(ishell, iset))
+               ic = ic + 1
+               IF (ic > ic_max) ic_max = ic
+            ENDDO
+            lshell = b%l(ishell, iset)
+            DO mm = -lshell, lshell
+               is = is + 1
+               IF (is > is_max) is_max = is
+            ENDDO
+         ENDDO
+      ENDDO
+
+      icgf = 0
+      isgf = 0
+      ALLOCATE (icgf_set(nset, ic_max))
+      icgf_set(:, :) = 0
+      ALLOCATE (isgf_set(nset, is_max))
+      isgf_set(:, :) = 0
+
+      DO iset = 1, nset
+         ic = 0
+         is = 0
+         DO ishell = 1, b%nshell(iset)
+            DO ico = 1, nco(b%l(ishell, iset))
+               icgf = icgf + 1
+               ic = ic + 1
+               icgf_set(iset, ic) = icgf
+            ENDDO
+            lshell = b%l(ishell, iset)
+            DO mm = -lshell, lshell
+               isgf = isgf + 1
+               is = is + 1
+               isgf_set(iset, is) = isgf
+            ENDDO
+         ENDDO
+      ENDDO
+
+      ALLOCATE (cgf_symbol(SIZE(b%cgf_symbol)))
+      ALLOCATE (norm_cgf(SIZE(b%norm_cgf)))
+      ALLOCATE (lx(SIZE(b%lx)))
+      ALLOCATE (ly(SIZE(b%ly)))
+      ALLOCATE (lz(SIZE(b%lz)))
+      ALLOCATE (cphi(SIZE(b%cphi, 1), SIZE(b%cphi, 2)))
+      cphi = 0.0_dp
+      ALLOCATE (sphi(SIZE(b%sphi, 1), SIZE(b%sphi, 2)))
+      sphi = 0.0_dp
+      ALLOCATE (scon(SIZE(b%scon, 1), SIZE(b%scon, 2)))
+      scon = 0.0_dp
+
+      ALLOCATE (sgf_symbol(SIZE(b%sgf_symbol)))
+      ALLOCATE (m(SIZE(b%m)))
+
+      icgf_new = 0
+      isgf_new = 0
+      DO iset = 1, nset
+         DO ic = 1, ic_max
+            icgf_old = icgf_set(sort_index(iset), ic)
+            IF (icgf_old == 0) CYCLE
+            icgf_new = icgf_new + 1
+            norm_cgf(icgf_new) = b%norm_cgf(icgf_old)
+            lx(icgf_new) = b%lx(icgf_old)
+            ly(icgf_new) = b%ly(icgf_old)
+            lz(icgf_new) = b%lz(icgf_old)
+            cphi(:, icgf_new) = b%cphi(:, icgf_old)
+            cgf_symbol(icgf_new) = b%cgf_symbol(icgf_old)
+         ENDDO
+         DO is = 1, is_max
+            isgf_old = isgf_set(sort_index(iset), is)
+            IF (isgf_old == 0) CYCLE
+            isgf_new = isgf_new + 1
+            m(isgf_new) = b%m(isgf_old)
+            sphi(:, isgf_new) = b%sphi(:, isgf_old)
+            scon(:, isgf_new) = b%scon(:, isgf_old)
+            sgf_symbol(isgf_new) = b%sgf_symbol(isgf_old)
+         ENDDO
+      ENDDO
+
+      DEALLOCATE (b%cgf_symbol)
+      b%cgf_symbol => cgf_symbol
+      DEALLOCATE (b%norm_cgf)
+      b%norm_cgf => norm_cgf
+      DEALLOCATE (b%lx)
+      b%lx => lx
+      DEALLOCATE (b%ly)
+      b%ly => ly
+      DEALLOCATE (b%lz)
+      b%lz => lz
+      DEALLOCATE (b%cphi)
+      b%cphi => cphi
+      DEALLOCATE (b%sphi)
+      b%sphi => sphi
+      DEALLOCATE (b%scon)
+      b%scon => scon
+
+      DEALLOCATE (b%m)
+      b%m => m
+      DEALLOCATE (b%sgf_symbol)
+      b%sgf_symbol => sgf_symbol
+
+      b%lmax = b%lmax(sort_index)
+      b%lmin = b%lmin(sort_index)
+      b%npgf = b%npgf(sort_index)
+      b%nshell = b%nshell(sort_index)
+      b%ncgf_set = b%ncgf_set(sort_index)
+      b%nsgf_set = b%nsgf_set(sort_index)
+
+      b%n(:, :) = b%n(:, sort_index)
+      b%l(:, :) = b%l(:, sort_index)
+      b%zet(:, :) = b%zet(:, sort_index)
+
+      b%gcc(:, :, :) = b%gcc(:, :, sort_index)
+      b%set_radius(:) = b%set_radius(sort_index)
+      b%pgf_radius(:, :) = b%pgf_radius(:, sort_index)
+
+      nc = 0
+      ns = 0
+      DO iset = 1, nset
+         DO ishell = 1, b%nshell(iset)
+            lshell = b%l(ishell, iset)
+            b%first_cgf(ishell, iset) = nc + 1
+            nc = nc + nco(lshell)
+            b%last_cgf(ishell, iset) = nc
+            b%first_sgf(ishell, iset) = ns + 1
+            ns = ns + nso(lshell)
+            b%last_sgf(ishell, iset) = ns
+         END DO
+      END DO
+
+      ! make sure that basis sets are sorted only once
+      b%sort = basis_sort_default
+
+   END SUBROUTINE
+
 ! **************************************************************************************************
 
 END MODULE basis_set_types

src/auto_basis.F

@@ -46,13 +46,15 @@ MODULE auto_basis
 !> \param ri_aux_basis_set ...
 !> \param qs_kind ...
 !> \param basis_cntrl ...
+!> \param basis_sort ...
 !> \date    01.11.2017
 !> \author  JGH
 ! **************************************************************************************************
-   SUBROUTINE create_ri_aux_basis_set(ri_aux_basis_set, qs_kind, basis_cntrl)
+   SUBROUTINE create_ri_aux_basis_set(ri_aux_basis_set, qs_kind, basis_cntrl, basis_sort)
       TYPE(gto_basis_set_type), POINTER                  :: ri_aux_basis_set
       TYPE(qs_kind_type), INTENT(IN)                     :: qs_kind
       INTEGER, INTENT(IN)                                :: basis_cntrl
+      INTEGER, INTENT(IN), OPTIONAL                      :: basis_sort
 
       CHARACTER(len=*), PARAMETER :: routineN = 'create_ri_aux_basis_set', &
          routineP = moduleN//':'//routineN
@@ -192,7 +194,9 @@ SUBROUTINE create_ri_aux_basis_set(ri_aux_basis_set, qs_kind, basis_cntrl)
          bsname = TRIM(element_symbol)//"-RI-AUX-"//TRIM(orb_basis_set%name)
          !
          CALL create_aux_basis(ri_aux_basis_set, bsname, nsets, ls1, ls2, nl, npgf, zet)
-         !
+
+         IF (PRESENT(basis_sort)) ri_aux_basis_set%sort = basis_sort
+
          DEALLOCATE (zet)
       END IF
 

src/input_cp2k_dft.F

@@ -10,6 +10,9 @@
 !> \author fawzi
 ! **************************************************************************************************
 MODULE input_cp2k_dft
+   USE basis_set_types,                 ONLY: basis_sort_default,&
+                                              basis_sort_radius,&
+                                              basis_sort_zet
    USE bibliography,                    ONLY: &
         Andermatt2016, Andreussi2012, Avezac2005, BaniHashemian2016, Becke1988b, Bengtsson1999, &
         Blochl1995, Brelaz1979, Dewar1977, Dewar1985, Dudarev1997, Dudarev1998, Ehrhardt1985, &
@@ -362,6 +365,18 @@ SUBROUTINE create_dft_section(section)
                           lone_keyword_l_val=.TRUE.)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
+
+      CALL keyword_create(keyword, __LOCATION__, &
+                          name="SORT_BASIS", &
+                          description="Sort basis sets according to a certain criterion.", &
+                          enum_c_vals=s2a("DEFAULT", "RADIUS", "EXP"), &
+                          enum_i_vals=(/basis_sort_default, basis_sort_radius, basis_sort_zet/), &
+                          enum_desc=s2a("don't sort", "sort w.r.t. set radius", "sort w.r.t. exponent"), &
+                          default_i_val=basis_sort_default, &
+                          usage="SORT_BASIS EXP")
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
       NULLIFY (subsection)
       CALL create_scf_section(subsection)
       CALL section_add_subsection(section, subsection)

src/mp2_integrals.F

@@ -946,7 +946,6 @@ SUBROUTINE create_contraction_batches(sizes, nbatches, starts_array, ends_array,
          bsum = bsum + SUM(sizes(starts_array_block(imem):ends_array_block(imem)))
          ends_array(imem) = bsum
       ENDDO
-
    END SUBROUTINE
 
 ! **************************************************************************************************

src/qs_environment.F

@@ -22,7 +22,8 @@ MODULE qs_environment
                                               create_lri_aux_basis_set,&
                                               create_ri_aux_basis_set
    USE basis_set_container_types,       ONLY: add_basis_set_to_container
-   USE basis_set_types,                 ONLY: get_gto_basis_set,&
+   USE basis_set_types,                 ONLY: basis_sort_zet,&
+                                              get_gto_basis_set,&
                                               gto_basis_set_type
    USE bibliography,                    ONLY: Iannuzzi2006,&
                                               Iannuzzi2007,&
@@ -503,10 +504,10 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
       CHARACTER(len=2)                                   :: element_symbol
       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, ngauss, nkind, output_unit, tnadd_method
+         nelectron, ngauss, nkind, output_unit, sort_basis, tnadd_method
       INTEGER, DIMENSION(2)                              :: n_mo, nelectron_spin
-      LOGICAL :: all_potential_present, be_silent, do_kpoints, e1terms, has_unit_metric, lribas, &
-         mp2_present, orb_gradient, was_present
+      LOGICAL :: all_potential_present, be_silent, do_kpoints, do_wfc_im_time, e1terms, &
+         has_unit_metric, lribas, mp2_present, orb_gradient
       REAL(dp)                                           :: alpha, ccore, ewald_rcut, maxocc, &
                                                             verlet_skin, zeff_correction
       REAL(KIND=dp)                                      :: total_zeff_corr
@@ -558,7 +559,6 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
       NULLIFY (logger)
       logger => cp_get_default_logger()
       output_unit = cp_logger_get_default_io_unit(logger)
-      was_present = .FALSE.
 
       be_silent = .FALSE.
       IF (PRESENT(silent)) be_silent = silent
@@ -812,10 +812,21 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
          END DO
       END IF
 
-      ! Check if RI_AUX basis (for MP2/RPA) is given, auto-generate if not
       mp2_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION")
       CALL section_vals_get(mp2_section, explicit=mp2_present)
       IF (mp2_present) THEN
+
+         ! basis should be sorted for imaginary time RPA/GW
+         CALL section_vals_val_get(qs_env%input, "DFT%SORT_BASIS", i_val=sort_basis)
+         CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%LOW_SCALING%_SECTION_PARAMETERS_", &
+                                   l_val=do_wfc_im_time)
+
+         IF (do_wfc_im_time .AND. sort_basis /= basis_sort_zet) THEN
+            CALL cp_warn(__LOCATION__, &
+                         "Low-scaling RPA requires SORT_BASIS EXP keyword (in DFT input section) for good performance")
+         ENDIF
+
+         ! Check if RI_AUX basis (for MP2/RPA) is given, auto-generate if not
          CALL mp2_env_create(qs_env%mp2_env)
          CALL get_qs_env(qs_env, mp2_env=mp2_env, nkind=nkind)
          CALL section_vals_val_get(mp2_section, "METHOD", i_val=mp2_env%method)
@@ -832,7 +843,7 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
                   CALL cp_warn(__LOCATION__, "Automatic Generation of RI_AUX basis. "// &
                                "This is experimental code.")
                   ! Generate a default basis
-                  CALL create_ri_aux_basis_set(ri_aux_basis_set, qs_kind, dft_control%auto_basis_ri_aux)
+                  CALL create_ri_aux_basis_set(ri_aux_basis_set, qs_kind, dft_control%auto_basis_ri_aux, basis_sort=sort_basis)
                   CALL add_basis_set_to_container(qs_kind%basis_sets, ri_aux_basis_set, "RI_AUX")
                END IF
             END DO

src/qs_interactions.F

@@ -19,7 +19,8 @@ MODULE qs_interactions
                                               get_atomic_kind
    USE basis_set_types,                 ONLY: get_gto_basis_set,&
                                               gto_basis_set_type,&
-                                              set_gto_basis_set
+                                              set_gto_basis_set,&
+                                              sort_gto_basis_set
    USE cp_control_types,                ONLY: qs_control_type,&
                                               semi_empirical_control_type
    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
@@ -431,6 +432,8 @@ SUBROUTINE init_interaction_radii_orb_basis(orb_basis_set, eps_pgf_orb, eps_pgf_
                                 kind_radius=kind_radius, &
                                 short_kind_radius=short_radius)
 
+         CALL sort_gto_basis_set(orb_basis_set)
+
       END IF
 
    END SUBROUTINE init_interaction_radii_orb_basis