LCOV - code coverage report
Current view: top level - src - qs_collocate_density.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 77.0 % 1071 825
Test Date: 2025-07-25 12:55:17 Functions: 94.1 % 17 16

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Calculate the plane wave density by collocating the primitive Gaussian
      10              : !>      functions (pgf).
      11              : !> \par History
      12              : !>      - rewrote collocate for increased accuracy and speed
      13              : !>      - introduced the PGI hack for increased speed with that compiler
      14              : !>        (22.02.02)
      15              : !>      - Added Multiple Grid feature
      16              : !>      - new way to get over the grid (01.03.02)
      17              : !>      - removed timing calls since they were getting expensive
      18              : !>      - Updated with the new QS data structures (09.04.02,MK)
      19              : !>      - introduction of the real space grid type ( prelim. version JVdV 05.02)
      20              : !>      - parallel FFT (JGH 22.05.02)
      21              : !>      - multigrid arrays independent from density (JGH 30.08.02)
      22              : !>      - old density stored in g space (JGH 30.08.02)
      23              : !>      - distributed real space code (JGH 17.07.03)
      24              : !>      - refactoring and new loop ordering (JGH 23.11.03)
      25              : !>      - OpenMP parallelization (JGH 03.12.03)
      26              : !>      - Modified to compute tau (Joost 12.03)
      27              : !>      - removed the incremental density rebuild (Joost 01.04)
      28              : !>      - introduced realspace multigridding (Joost 02.04)
      29              : !>      - introduced map_consistent (Joost 02.04)
      30              : !>      - Addition of the subroutine calculate_atomic_charge_density (TdK, 08.05)
      31              : !>      - rewrite of the collocate/integrate kernels (Joost VandeVondele, 03.07)
      32              : !>      - Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
      33              : !> \author Matthias Krack (03.04.2001)
      34              : !>      1) Joost VandeVondele (01.2002)
      35              : !>      Thomas D. Kuehne (04.08.2005)
      36              : !>      Ole Schuett (2020)
      37              : ! **************************************************************************************************
      38              : MODULE qs_collocate_density
      39              :    USE admm_types, ONLY: get_admm_env
      40              :    USE ao_util, ONLY: exp_radius_very_extended
      41              :    USE atomic_kind_types, ONLY: atomic_kind_type, &
      42              :                                 get_atomic_kind, &
      43              :                                 get_atomic_kind_set
      44              :    USE basis_set_types, ONLY: get_gto_basis_set, &
      45              :                               gto_basis_set_type
      46              :    USE cell_types, ONLY: cell_type, &
      47              :                          pbc
      48              :    USE cp_control_types, ONLY: dft_control_type
      49              :    USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set
      50              :    USE cp_fm_types, ONLY: cp_fm_get_element, &
      51              :                           cp_fm_get_info, &
      52              :                           cp_fm_type
      53              :    USE cp_dbcsr_api, ONLY: dbcsr_copy, &
      54              :                            dbcsr_get_block_p, &
      55              :                            dbcsr_p_type, &
      56              :                            dbcsr_type
      57              :    USE external_potential_types, ONLY: get_potential, &
      58              :                                        gth_potential_type
      59              :    USE gaussian_gridlevels, ONLY: gaussian_gridlevel, &
      60              :                                   gridlevel_info_type
      61              :    USE grid_api, ONLY: &
      62              :       GRID_FUNC_AB, GRID_FUNC_CORE_X, GRID_FUNC_CORE_Y, GRID_FUNC_CORE_Z, GRID_FUNC_DAB_X, &
      63              :       GRID_FUNC_DAB_Y, GRID_FUNC_DAB_Z, GRID_FUNC_DABpADB_X, GRID_FUNC_DABpADB_Y, &
      64              :       GRID_FUNC_DABpADB_Z, GRID_FUNC_DADB, GRID_FUNC_DX, GRID_FUNC_DXDX, GRID_FUNC_DXDY, &
      65              :       GRID_FUNC_DY, GRID_FUNC_DYDY, GRID_FUNC_DYDZ, GRID_FUNC_DZ, GRID_FUNC_DZDX, &
      66              :       GRID_FUNC_DZDZ, collocate_pgf_product, grid_collocate_task_list
      67              :    USE input_constants, ONLY: &
      68              :       orb_dx2, orb_dxy, orb_dy2, orb_dyz, orb_dz2, orb_dzx, orb_px, orb_py, orb_pz, orb_s
      69              :    USE kinds, ONLY: default_string_length, &
      70              :                     dp
      71              :    USE lri_environment_types, ONLY: lri_kind_type
      72              :    USE memory_utilities, ONLY: reallocate
      73              :    USE message_passing, ONLY: mp_comm_type
      74              :    USE orbital_pointers, ONLY: coset, &
      75              :                                ncoset
      76              :    USE particle_types, ONLY: particle_type
      77              :    USE pw_env_types, ONLY: pw_env_get, &
      78              :                            pw_env_type
      79              :    USE pw_methods, ONLY: pw_axpy, &
      80              :                          pw_integrate_function, &
      81              :                          pw_transfer, &
      82              :                          pw_zero
      83              :    USE pw_pool_types, ONLY: pw_pool_p_type, &
      84              :                             pw_pool_type, &
      85              :                             pw_pools_create_pws, &
      86              :                             pw_pools_give_back_pws
      87              :    USE pw_types, ONLY: pw_r3d_rs_type, &
      88              :                        pw_c1d_gs_type, &
      89              :                        pw_r3d_rs_type
      90              :    USE qs_environment_types, ONLY: get_qs_env, &
      91              :                                    qs_environment_type
      92              :    USE qs_kind_types, ONLY: get_qs_kind, &
      93              :                             get_qs_kind_set, &
      94              :                             qs_kind_type
      95              :    USE qs_ks_types, ONLY: get_ks_env, &
      96              :                           qs_ks_env_type
      97              :    USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
      98              :    USE realspace_grid_types, ONLY: map_gaussian_here, &
      99              :                                    realspace_grid_desc_p_type, &
     100              :                                    realspace_grid_type, &
     101              :                                    rs_grid_zero, &
     102              :                                    transfer_rs2pw
     103              :    USE rs_pw_interface, ONLY: density_rs2pw
     104              :    USE task_list_methods, ONLY: rs_copy_to_buffer, &
     105              :                                 rs_distribute_matrix, &
     106              :                                 rs_scatter_matrices
     107              :    USE task_list_types, ONLY: atom_pair_type, &
     108              :                               task_list_type, &
     109              :                               task_type
     110              : 
     111              : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
     112              : 
     113              : #include "./base/base_uses.f90"
     114              : 
     115              :    IMPLICIT NONE
     116              : 
     117              :    PRIVATE
     118              : 
     119              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_collocate_density'
     120              : ! *** Public subroutines ***
     121              : 
     122              :    PUBLIC :: calculate_ppl_grid, &
     123              :              calculate_rho_core, &
     124              :              calculate_lri_rho_elec, &
     125              :              calculate_rho_single_gaussian, &
     126              :              calculate_rho_metal, &
     127              :              calculate_rho_resp_single, &
     128              :              calculate_rho_resp_all, &
     129              :              calculate_rho_elec, &
     130              :              calculate_drho_elec, &
     131              :              calculate_wavefunction, &
     132              :              collocate_function, &
     133              :              calculate_rho_nlcc, &
     134              :              calculate_drho_elec_dR, &
     135              :              calculate_drho_core, &
     136              :              collocate_single_gaussian
     137              : 
     138              :    INTERFACE calculate_rho_core
     139              :       MODULE PROCEDURE calculate_rho_core_r3d_rs
     140              :       MODULE PROCEDURE calculate_rho_core_c1d_gs
     141              :    END INTERFACE
     142              : 
     143              :    INTERFACE calculate_rho_resp_all
     144              :       MODULE PROCEDURE calculate_rho_resp_all_r3d_rs, calculate_rho_resp_all_c1d_gs
     145              :    END INTERFACE
     146              : 
     147              : CONTAINS
     148              : 
     149              : ! **************************************************************************************************
     150              : !> \brief computes the density of the non-linear core correction on the grid
     151              : !> \param rho_nlcc ...
     152              : !> \param qs_env ...
     153              : ! **************************************************************************************************
     154           36 :    SUBROUTINE calculate_rho_nlcc(rho_nlcc, qs_env)
     155              : 
     156              :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                       :: rho_nlcc
     157              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     158              : 
     159              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_nlcc'
     160              : 
     161              :       INTEGER                                            :: atom_a, handle, iatom, iexp_nlcc, ikind, &
     162              :                                                             ithread, j, n, natom, nc, nexp_nlcc, &
     163              :                                                             ni, npme, nthread, subpatch_pattern
     164           36 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list, cores, nct_nlcc
     165              :       LOGICAL                                            :: nlcc
     166              :       REAL(KIND=dp)                                      :: alpha, eps_rho_rspace, radius
     167              :       REAL(KIND=dp), DIMENSION(3)                        :: ra
     168           36 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: alpha_nlcc
     169           36 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: cval_nlcc, pab
     170           36 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     171              :       TYPE(cell_type), POINTER                           :: cell
     172              :       TYPE(dft_control_type), POINTER                    :: dft_control
     173              :       TYPE(gth_potential_type), POINTER                  :: gth_potential
     174           36 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     175              :       TYPE(pw_env_type), POINTER                         :: pw_env
     176              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     177           36 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     178              :       TYPE(realspace_grid_type), POINTER                 :: rs_rho
     179              : 
     180           36 :       CALL timeset(routineN, handle)
     181              : 
     182           36 :       NULLIFY (cell, dft_control, pab, particle_set, atomic_kind_set, &
     183           36 :                qs_kind_set, atom_list, pw_env, rs_rho, auxbas_pw_pool, cores)
     184              : 
     185              :       CALL get_qs_env(qs_env=qs_env, &
     186              :                       atomic_kind_set=atomic_kind_set, &
     187              :                       qs_kind_set=qs_kind_set, &
     188              :                       cell=cell, &
     189              :                       dft_control=dft_control, &
     190              :                       particle_set=particle_set, &
     191           36 :                       pw_env=pw_env)
     192              :       CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
     193           36 :                       auxbas_pw_pool=auxbas_pw_pool)
     194              :       ! be careful in parallel nsmax is chosen with multigrid in mind!
     195           36 :       CALL rs_grid_zero(rs_rho)
     196              : 
     197           36 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
     198              : 
     199           92 :       DO ikind = 1, SIZE(atomic_kind_set)
     200           56 :          CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
     201           56 :          CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential)
     202              : 
     203           56 :          IF (.NOT. ASSOCIATED(gth_potential)) CYCLE
     204              :          CALL get_potential(potential=gth_potential, nlcc_present=nlcc, nexp_nlcc=nexp_nlcc, &
     205           56 :                             alpha_nlcc=alpha_nlcc, nct_nlcc=nct_nlcc, cval_nlcc=cval_nlcc)
     206              : 
     207           56 :          IF (.NOT. nlcc) CYCLE
     208              : 
     209          256 :          DO iexp_nlcc = 1, nexp_nlcc
     210              : 
     211           54 :             alpha = alpha_nlcc(iexp_nlcc)
     212           54 :             nc = nct_nlcc(iexp_nlcc)
     213              : 
     214           54 :             ni = ncoset(2*nc - 2)
     215          162 :             ALLOCATE (pab(ni, 1))
     216          306 :             pab = 0._dp
     217              : 
     218           54 :             nthread = 1
     219           54 :             ithread = 0
     220              : 
     221           54 :             CALL reallocate(cores, 1, natom)
     222           54 :             npme = 0
     223          232 :             cores = 0
     224              : 
     225              :             ! prepare core function
     226          124 :             DO j = 1, nc
     227           54 :                SELECT CASE (j)
     228              :                CASE (1)
     229           54 :                   pab(1, 1) = cval_nlcc(1, iexp_nlcc)
     230              :                CASE (2)
     231           16 :                   n = coset(2, 0, 0)
     232           16 :                   pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
     233           16 :                   n = coset(0, 2, 0)
     234           16 :                   pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
     235           16 :                   n = coset(0, 0, 2)
     236           16 :                   pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
     237              :                CASE (3)
     238            0 :                   n = coset(4, 0, 0)
     239            0 :                   pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
     240            0 :                   n = coset(0, 4, 0)
     241            0 :                   pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
     242            0 :                   n = coset(0, 0, 4)
     243            0 :                   pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
     244            0 :                   n = coset(2, 2, 0)
     245            0 :                   pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
     246            0 :                   n = coset(2, 0, 2)
     247            0 :                   pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
     248            0 :                   n = coset(0, 2, 2)
     249            0 :                   pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
     250              :                CASE (4)
     251            0 :                   n = coset(6, 0, 0)
     252            0 :                   pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
     253            0 :                   n = coset(0, 6, 0)
     254            0 :                   pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
     255            0 :                   n = coset(0, 0, 6)
     256            0 :                   pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
     257            0 :                   n = coset(4, 2, 0)
     258            0 :                   pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
     259            0 :                   n = coset(4, 0, 2)
     260            0 :                   pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
     261            0 :                   n = coset(2, 4, 0)
     262            0 :                   pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
     263            0 :                   n = coset(2, 0, 4)
     264            0 :                   pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
     265            0 :                   n = coset(0, 4, 2)
     266            0 :                   pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
     267            0 :                   n = coset(0, 2, 4)
     268            0 :                   pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
     269            0 :                   n = coset(2, 2, 2)
     270            0 :                   pab(n, 1) = 6._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
     271              :                CASE DEFAULT
     272           70 :                   CPABORT("")
     273              :                END SELECT
     274              :             END DO
     275           54 :             IF (dft_control%nspins == 2) pab = pab*0.5_dp
     276              : 
     277          232 :             DO iatom = 1, natom
     278          178 :                atom_a = atom_list(iatom)
     279          178 :                ra(:) = pbc(particle_set(atom_a)%r, cell)
     280          232 :                IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
     281              :                   ! replicated realspace grid, split the atoms up between procs
     282          178 :                   IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
     283           89 :                      npme = npme + 1
     284           89 :                      cores(npme) = iatom
     285              :                   END IF
     286              :                ELSE
     287            0 :                   npme = npme + 1
     288            0 :                   cores(npme) = iatom
     289              :                END IF
     290              :             END DO
     291              : 
     292          143 :             DO j = 1, npme
     293              : 
     294           89 :                iatom = cores(j)
     295           89 :                atom_a = atom_list(iatom)
     296           89 :                ra(:) = pbc(particle_set(atom_a)%r, cell)
     297           89 :                subpatch_pattern = 0
     298           89 :                ni = 2*nc - 2
     299              :                radius = exp_radius_very_extended(la_min=0, la_max=ni, lb_min=0, lb_max=0, &
     300              :                                                  ra=ra, rb=ra, rp=ra, &
     301              :                                                  zetp=1/(2*alpha**2), eps=eps_rho_rspace, &
     302              :                                                  pab=pab, o1=0, o2=0, &  ! without map_consistent
     303           89 :                                                  prefactor=1.0_dp, cutoff=0.0_dp)
     304              : 
     305              :                CALL collocate_pgf_product(ni, 1/(2*alpha**2), 0, 0, 0.0_dp, 0, ra, &
     306              :                                           (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, rs_rho, &
     307              :                                           ga_gb_function=GRID_FUNC_AB, radius=radius, &
     308          143 :                                           use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern)
     309              : 
     310              :             END DO
     311              : 
     312          110 :             DEALLOCATE (pab)
     313              : 
     314              :          END DO
     315              : 
     316              :       END DO
     317              : 
     318           36 :       IF (ASSOCIATED(cores)) THEN
     319           36 :          DEALLOCATE (cores)
     320              :       END IF
     321              : 
     322           36 :       CALL transfer_rs2pw(rs_rho, rho_nlcc)
     323              : 
     324           36 :       CALL timestop(handle)
     325              : 
     326           36 :    END SUBROUTINE calculate_rho_nlcc
     327              : 
     328              : ! **************************************************************************************************
     329              : !> \brief computes the local pseudopotential (without erf term) on the grid
     330              : !> \param vppl ...
     331              : !> \param qs_env ...
     332              : ! **************************************************************************************************
     333           12 :    SUBROUTINE calculate_ppl_grid(vppl, qs_env)
     334              : 
     335              :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                       :: vppl
     336              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     337              : 
     338              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_ppl_grid'
     339              : 
     340              :       INTEGER                                            :: atom_a, handle, iatom, ikind, ithread, &
     341              :                                                             j, lppl, n, natom, ni, npme, nthread, &
     342              :                                                             subpatch_pattern
     343           12 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list, cores
     344              :       REAL(KIND=dp)                                      :: alpha, eps_rho_rspace, radius
     345              :       REAL(KIND=dp), DIMENSION(3)                        :: ra
     346           12 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: cexp_ppl
     347           12 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: pab
     348           12 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     349              :       TYPE(cell_type), POINTER                           :: cell
     350              :       TYPE(dft_control_type), POINTER                    :: dft_control
     351              :       TYPE(gth_potential_type), POINTER                  :: gth_potential
     352           12 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     353              :       TYPE(pw_env_type), POINTER                         :: pw_env
     354              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     355           12 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     356              :       TYPE(realspace_grid_type), POINTER                 :: rs_rho
     357              : 
     358           12 :       CALL timeset(routineN, handle)
     359              : 
     360           12 :       NULLIFY (cell, dft_control, pab, atomic_kind_set, qs_kind_set, particle_set, &
     361           12 :                atom_list, pw_env, rs_rho, auxbas_pw_pool, cores)
     362              : 
     363              :       CALL get_qs_env(qs_env=qs_env, &
     364              :                       atomic_kind_set=atomic_kind_set, &
     365              :                       qs_kind_set=qs_kind_set, &
     366              :                       cell=cell, &
     367              :                       dft_control=dft_control, &
     368              :                       particle_set=particle_set, &
     369           12 :                       pw_env=pw_env)
     370              :       CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
     371           12 :                       auxbas_pw_pool=auxbas_pw_pool)
     372              :       ! be careful in parallel nsmax is chosen with multigrid in mind!
     373           12 :       CALL rs_grid_zero(rs_rho)
     374              : 
     375           12 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
     376              : 
     377           28 :       DO ikind = 1, SIZE(atomic_kind_set)
     378           16 :          CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
     379           16 :          CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential)
     380              : 
     381           16 :          IF (.NOT. ASSOCIATED(gth_potential)) CYCLE
     382           16 :          CALL get_potential(potential=gth_potential, alpha_ppl=alpha, nexp_ppl=lppl, cexp_ppl=cexp_ppl)
     383              : 
     384           16 :          IF (lppl <= 0) CYCLE
     385              : 
     386           16 :          ni = ncoset(2*lppl - 2)
     387           48 :          ALLOCATE (pab(ni, 1))
     388          192 :          pab = 0._dp
     389              : 
     390           16 :          nthread = 1
     391           16 :          ithread = 0
     392              : 
     393           16 :          CALL reallocate(cores, 1, natom)
     394           16 :          npme = 0
     395           60 :          cores = 0
     396              : 
     397              :          ! prepare core function
     398           48 :          DO j = 1, lppl
     399           16 :             SELECT CASE (j)
     400              :             CASE (1)
     401           16 :                pab(1, 1) = cexp_ppl(1)
     402              :             CASE (2)
     403           16 :                n = coset(2, 0, 0)
     404           16 :                pab(n, 1) = cexp_ppl(2)
     405           16 :                n = coset(0, 2, 0)
     406           16 :                pab(n, 1) = cexp_ppl(2)
     407           16 :                n = coset(0, 0, 2)
     408           16 :                pab(n, 1) = cexp_ppl(2)
     409              :             CASE (3)
     410            0 :                n = coset(4, 0, 0)
     411            0 :                pab(n, 1) = cexp_ppl(3)
     412            0 :                n = coset(0, 4, 0)
     413            0 :                pab(n, 1) = cexp_ppl(3)
     414            0 :                n = coset(0, 0, 4)
     415            0 :                pab(n, 1) = cexp_ppl(3)
     416            0 :                n = coset(2, 2, 0)
     417            0 :                pab(n, 1) = 2._dp*cexp_ppl(3)
     418            0 :                n = coset(2, 0, 2)
     419            0 :                pab(n, 1) = 2._dp*cexp_ppl(3)
     420            0 :                n = coset(0, 2, 2)
     421            0 :                pab(n, 1) = 2._dp*cexp_ppl(3)
     422              :             CASE (4)
     423            0 :                n = coset(6, 0, 0)
     424            0 :                pab(n, 1) = cexp_ppl(4)
     425            0 :                n = coset(0, 6, 0)
     426            0 :                pab(n, 1) = cexp_ppl(4)
     427            0 :                n = coset(0, 0, 6)
     428            0 :                pab(n, 1) = cexp_ppl(4)
     429            0 :                n = coset(4, 2, 0)
     430            0 :                pab(n, 1) = 3._dp*cexp_ppl(4)
     431            0 :                n = coset(4, 0, 2)
     432            0 :                pab(n, 1) = 3._dp*cexp_ppl(4)
     433            0 :                n = coset(2, 4, 0)
     434            0 :                pab(n, 1) = 3._dp*cexp_ppl(4)
     435            0 :                n = coset(2, 0, 4)
     436            0 :                pab(n, 1) = 3._dp*cexp_ppl(4)
     437            0 :                n = coset(0, 4, 2)
     438            0 :                pab(n, 1) = 3._dp*cexp_ppl(4)
     439            0 :                n = coset(0, 2, 4)
     440            0 :                pab(n, 1) = 3._dp*cexp_ppl(4)
     441            0 :                n = coset(2, 2, 2)
     442            0 :                pab(n, 1) = 6._dp*cexp_ppl(4)
     443              :             CASE DEFAULT
     444           32 :                CPABORT("")
     445              :             END SELECT
     446              :          END DO
     447              : 
     448           60 :          DO iatom = 1, natom
     449           44 :             atom_a = atom_list(iatom)
     450           44 :             ra(:) = pbc(particle_set(atom_a)%r, cell)
     451           60 :             IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
     452              :                ! replicated realspace grid, split the atoms up between procs
     453           44 :                IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
     454           22 :                   npme = npme + 1
     455           22 :                   cores(npme) = iatom
     456              :                END IF
     457              :             ELSE
     458            0 :                npme = npme + 1
     459            0 :                cores(npme) = iatom
     460              :             END IF
     461              :          END DO
     462              : 
     463           16 :          IF (npme .GT. 0) THEN
     464           36 :             DO j = 1, npme
     465              : 
     466           22 :                iatom = cores(j)
     467           22 :                atom_a = atom_list(iatom)
     468           22 :                ra(:) = pbc(particle_set(atom_a)%r, cell)
     469           22 :                subpatch_pattern = 0
     470           22 :                ni = 2*lppl - 2
     471              : 
     472              :                radius = exp_radius_very_extended(la_min=0, la_max=ni, &
     473              :                                                  lb_min=0, lb_max=0, &
     474              :                                                  ra=ra, rb=ra, rp=ra, &
     475              :                                                  zetp=alpha, eps=eps_rho_rspace, &
     476              :                                                  pab=pab, o1=0, o2=0, &  ! without map_consistent
     477           22 :                                                  prefactor=1.0_dp, cutoff=0.0_dp)
     478              : 
     479              :                CALL collocate_pgf_product(ni, alpha, 0, 0, 0.0_dp, 0, ra, &
     480              :                                           (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, rs_rho, &
     481              :                                           radius=radius, ga_gb_function=GRID_FUNC_AB, &
     482           36 :                                           use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern)
     483              : 
     484              :             END DO
     485              :          END IF
     486              : 
     487           60 :          DEALLOCATE (pab)
     488              : 
     489              :       END DO
     490              : 
     491           12 :       IF (ASSOCIATED(cores)) THEN
     492           12 :          DEALLOCATE (cores)
     493              :       END IF
     494              : 
     495           12 :       CALL transfer_rs2pw(rs_rho, vppl)
     496              : 
     497           12 :       CALL timestop(handle)
     498              : 
     499           12 :    END SUBROUTINE calculate_ppl_grid
     500              : 
     501              : ! **************************************************************************************************
     502              : !> \brief Collocates the fitted lri density on a grid.
     503              : !> \param lri_rho_g ...
     504              : !> \param lri_rho_r ...
     505              : !> \param qs_env ...
     506              : !> \param lri_coef ...
     507              : !> \param total_rho ...
     508              : !> \param basis_type ...
     509              : !> \param exact_1c_terms ...
     510              : !> \param pmat replicated block diagonal density matrix (optional)
     511              : !> \param atomlist list of atoms to be included (optional)
     512              : !> \par History
     513              : !>      04.2013
     514              : !> \author Dorothea Golze
     515              : ! **************************************************************************************************
     516          908 :    SUBROUTINE calculate_lri_rho_elec(lri_rho_g, lri_rho_r, qs_env, &
     517          908 :                                      lri_coef, total_rho, basis_type, exact_1c_terms, pmat, atomlist)
     518              : 
     519              :       TYPE(pw_c1d_gs_type), INTENT(INOUT) :: lri_rho_g
     520              :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                       ::  lri_rho_r
     521              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     522              :       TYPE(lri_kind_type), DIMENSION(:), POINTER         :: lri_coef
     523              :       REAL(KIND=dp), INTENT(OUT)                         :: total_rho
     524              :       CHARACTER(len=*), INTENT(IN)                       :: basis_type
     525              :       LOGICAL, INTENT(IN)                                :: exact_1c_terms
     526              :       TYPE(dbcsr_type), OPTIONAL                         :: pmat
     527              :       INTEGER, DIMENSION(:), OPTIONAL                    :: atomlist
     528              : 
     529              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_lri_rho_elec'
     530              : 
     531              :       INTEGER :: atom_a, group_size, handle, iatom, igrid_level, ikind, ipgf, iset, jpgf, jset, &
     532              :                  m1, maxco, maxsgf_set, my_pos, na1, natom, nb1, ncoa, ncob, nseta, offset, sgfa, sgfb
     533          908 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list, la_max, la_min, npgfa, nsgfa
     534          908 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa
     535              :       LOGICAL                                            :: found
     536          908 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: map_it
     537          908 :       LOGICAL, ALLOCATABLE, DIMENSION(:, :)              :: map_it2
     538              :       REAL(KIND=dp)                                      :: eps_rho_rspace, radius, zetp
     539              :       REAL(KIND=dp), DIMENSION(3)                        :: ra
     540          908 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: aci
     541          908 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: p_block, pab, sphi_a, work, zeta
     542          908 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     543              :       TYPE(cell_type), POINTER                           :: cell
     544              :       TYPE(dft_control_type), POINTER                    :: dft_control
     545              :       TYPE(gridlevel_info_type), POINTER                 :: gridlevel_info
     546              :       TYPE(gto_basis_set_type), POINTER                  :: lri_basis_set, orb_basis_set
     547          908 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     548              :       TYPE(pw_env_type), POINTER                         :: pw_env
     549          908 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
     550          908 :       TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:) :: mgrid_gspace
     551          908 :       TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:)           ::  mgrid_rspace
     552          908 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     553          908 :       TYPE(realspace_grid_type), DIMENSION(:), POINTER   :: rs_rho
     554              :       TYPE(realspace_grid_type), POINTER                 :: rs_grid
     555              : 
     556          908 :       NULLIFY (aci, atomic_kind_set, qs_kind_set, atom_list, cell, &
     557          908 :                dft_control, first_sgfa, gridlevel_info, la_max, &
     558          908 :                la_min, lri_basis_set, npgfa, nsgfa, &
     559          908 :                pab, particle_set, pw_env, pw_pools, rs_grid, rs_rho, sphi_a, &
     560          908 :                work, zeta)
     561              : 
     562          908 :       CALL timeset(routineN, handle)
     563              : 
     564          908 :       IF (exact_1c_terms) THEN
     565           48 :          CPASSERT(PRESENT(pmat))
     566              :       END IF
     567              : 
     568              :       CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
     569              :                       atomic_kind_set=atomic_kind_set, &
     570              :                       cell=cell, particle_set=particle_set, &
     571              :                       pw_env=pw_env, &
     572          908 :                       dft_control=dft_control)
     573              : 
     574          908 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
     575          908 :       gridlevel_info => pw_env%gridlevel_info
     576              : 
     577              :       ! *** set up the pw multi-grids *** !
     578          908 :       CPASSERT(ASSOCIATED(pw_env))
     579          908 :       CALL pw_env_get(pw_env=pw_env, rs_grids=rs_rho, pw_pools=pw_pools)
     580              : 
     581          908 :       CALL pw_pools_create_pws(pw_pools, mgrid_rspace)
     582              : 
     583          908 :       CALL pw_pools_create_pws(pw_pools, mgrid_gspace)
     584              : 
     585              :       ! *** set up the rs multi-grids *** !
     586         4480 :       DO igrid_level = 1, gridlevel_info%ngrid_levels
     587         4480 :          CALL rs_grid_zero(rs_rho(igrid_level))
     588              :       END DO
     589              : 
     590              :       !take maxco from the LRI basis set!
     591              :       CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
     592          908 :                            maxco=maxco, basis_type=basis_type)
     593              : 
     594         2724 :       ALLOCATE (pab(maxco, 1))
     595          908 :       offset = 0
     596          908 :       my_pos = mgrid_rspace(1)%pw_grid%para%group%mepos
     597          908 :       group_size = mgrid_rspace(1)%pw_grid%para%group%num_pe
     598              : 
     599         2718 :       DO ikind = 1, SIZE(atomic_kind_set)
     600              : 
     601         1810 :          CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
     602         1810 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis_set, basis_type=basis_type)
     603              : 
     604              :          !Take the lri basis set here!
     605              :          CALL get_gto_basis_set(gto_basis_set=lri_basis_set, lmax=la_max, &
     606              :                                 lmin=la_min, zet=zeta, nset=nseta, npgf=npgfa, &
     607         1810 :                                 sphi=sphi_a, first_sgf=first_sgfa, nsgf_set=nsgfa)
     608              : 
     609         7950 :          DO iatom = 1, natom
     610         3422 :             atom_a = atom_list(iatom)
     611         3422 :             IF (PRESENT(ATOMLIST)) THEN
     612          500 :                IF (atomlist(atom_a) == 0) CYCLE
     613              :             END IF
     614         3222 :             ra(:) = pbc(particle_set(atom_a)%r, cell)
     615         3222 :             aci => lri_coef(ikind)%acoef(iatom, :)
     616              : 
     617        48280 :             m1 = MAXVAL(npgfa(1:nseta))
     618         9666 :             ALLOCATE (map_it(m1))
     619        48280 :             DO iset = 1, nseta
     620              :                ! collocate this set locally?
     621        95420 :                map_it = .FALSE.
     622        95228 :                DO ipgf = 1, npgfa(iset)
     623        50170 :                   igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
     624        50170 :                   rs_grid => rs_rho(igrid_level)
     625        95228 :                   map_it(ipgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
     626              :                END DO
     627        45058 :                offset = offset + 1
     628              : 
     629        73365 :                IF (ANY(map_it(1:npgfa(iset)))) THEN
     630        22529 :                   sgfa = first_sgfa(1, iset)
     631        22529 :                   ncoa = npgfa(iset)*ncoset(la_max(iset))
     632        22529 :                   m1 = sgfa + nsgfa(iset) - 1
     633        67587 :                   ALLOCATE (work(nsgfa(iset), 1))
     634       373257 :                   work(1:nsgfa(iset), 1) = aci(sgfa:m1)
     635       668770 :                   pab = 0._dp
     636              : 
     637              :                   CALL dgemm("N", "N", ncoa, 1, nsgfa(iset), 1.0_dp, lri_basis_set%sphi(1, sgfa), &
     638              :                              SIZE(lri_basis_set%sphi, 1), work(1, 1), SIZE(work, 1), 0.0_dp, pab(1, 1), &
     639        22529 :                              SIZE(pab, 1))
     640              : 
     641        47614 :                   DO ipgf = 1, npgfa(iset)
     642        25085 :                      na1 = (ipgf - 1)*ncoset(la_max(iset))
     643        25085 :                      igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
     644        25085 :                      rs_grid => rs_rho(igrid_level)
     645        47614 :                      IF (map_it(ipgf)) THEN
     646              :                         radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
     647              :                                                           lb_min=0, lb_max=0, &
     648              :                                                           ra=ra, rb=ra, rp=ra, &
     649              :                                                           zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
     650        25085 :                                                           prefactor=1.0_dp, cutoff=1.0_dp)
     651              : 
     652              :                         CALL collocate_pgf_product(la_max=la_max(iset), &
     653              :                                                    zeta=zeta(ipgf, iset), &
     654              :                                                    la_min=la_min(iset), &
     655              :                                                    lb_max=0, zetb=0.0_dp, lb_min=0, &
     656              :                                                    ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
     657              :                                                    scale=1._dp, &
     658              :                                                    pab=pab, o1=na1, o2=0, &
     659              :                                                    rsgrid=rs_grid, &
     660              :                                                    radius=radius, &
     661        25085 :                                                    ga_gb_function=GRID_FUNC_AB)
     662              :                      END IF
     663              :                   END DO
     664        22529 :                   DEALLOCATE (work)
     665              :                END IF
     666              :             END DO
     667         5232 :             DEALLOCATE (map_it)
     668              :          END DO
     669              :       END DO
     670              : 
     671          908 :       DEALLOCATE (pab)
     672              : 
     673              :       ! process the one-center terms
     674          908 :       IF (exact_1c_terms) THEN
     675              :          ! find maximum numbers
     676           48 :          offset = 0
     677              :          CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
     678              :                               maxco=maxco, &
     679              :                               maxsgf_set=maxsgf_set, &
     680           48 :                               basis_type="ORB")
     681          336 :          ALLOCATE (pab(maxco, maxco), work(maxco, maxsgf_set))
     682              : 
     683          144 :          DO ikind = 1, SIZE(atomic_kind_set)
     684           96 :             CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
     685           96 :             CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB")
     686              :             CALL get_gto_basis_set(gto_basis_set=orb_basis_set, lmax=la_max, &
     687              :                                    lmin=la_min, zet=zeta, nset=nseta, npgf=npgfa, &
     688           96 :                                    sphi=sphi_a, first_sgf=first_sgfa, nsgf_set=nsgfa)
     689          528 :             DO iatom = 1, natom
     690          288 :                atom_a = atom_list(iatom)
     691          288 :                ra(:) = pbc(particle_set(atom_a)%r, cell)
     692          288 :                CALL dbcsr_get_block_p(matrix=pmat, row=atom_a, col=atom_a, BLOCK=p_block, found=found)
     693          576 :                m1 = MAXVAL(npgfa(1:nseta))
     694         1152 :                ALLOCATE (map_it2(m1, m1))
     695          576 :                DO iset = 1, nseta
     696          864 :                   DO jset = 1, nseta
     697              :                      ! processor mappint
     698        16416 :                      map_it2 = .FALSE.
     699         2304 :                      DO ipgf = 1, npgfa(iset)
     700        16416 :                         DO jpgf = 1, npgfa(jset)
     701        14112 :                            zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
     702        14112 :                            igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
     703        14112 :                            rs_grid => rs_rho(igrid_level)
     704        16128 :                            map_it2(ipgf, jpgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
     705              :                         END DO
     706              :                      END DO
     707          288 :                      offset = offset + 1
     708              :                      !
     709         8640 :                      IF (ANY(map_it2(1:npgfa(iset), 1:npgfa(jset)))) THEN
     710          144 :                         ncoa = npgfa(iset)*ncoset(la_max(iset))
     711          144 :                         sgfa = first_sgfa(1, iset)
     712          144 :                         ncob = npgfa(jset)*ncoset(la_max(jset))
     713          144 :                         sgfb = first_sgfa(1, jset)
     714              :                         ! decontract density block
     715              :                         CALL dgemm("N", "N", ncoa, nsgfa(jset), nsgfa(iset), &
     716              :                                    1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
     717              :                                    p_block(sgfa, sgfb), SIZE(p_block, 1), &
     718          144 :                                    0.0_dp, work(1, 1), maxco)
     719              :                         CALL dgemm("N", "T", ncoa, ncob, nsgfa(jset), &
     720              :                                    1.0_dp, work(1, 1), maxco, &
     721              :                                    sphi_a(1, sgfb), SIZE(sphi_a, 1), &
     722          144 :                                    0.0_dp, pab(1, 1), maxco)
     723         1152 :                         DO ipgf = 1, npgfa(iset)
     724         8208 :                            DO jpgf = 1, npgfa(jset)
     725         7056 :                               zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
     726         7056 :                               igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
     727         7056 :                               rs_grid => rs_rho(igrid_level)
     728              : 
     729         7056 :                               na1 = (ipgf - 1)*ncoset(la_max(iset))
     730         7056 :                               nb1 = (jpgf - 1)*ncoset(la_max(jset))
     731              : 
     732         8064 :                               IF (map_it2(ipgf, jpgf)) THEN
     733              :                                  radius = exp_radius_very_extended(la_min=la_min(iset), &
     734              :                                                                    la_max=la_max(iset), &
     735              :                                                                    lb_min=la_min(jset), &
     736              :                                                                    lb_max=la_max(jset), &
     737              :                                                                    ra=ra, rb=ra, rp=ra, &
     738              :                                                                    zetp=zetp, eps=eps_rho_rspace, &
     739         7056 :                                                                    prefactor=1.0_dp, cutoff=1.0_dp)
     740              : 
     741              :                                  CALL collocate_pgf_product( &
     742              :                                     la_max(iset), zeta(ipgf, iset), la_min(iset), &
     743              :                                     la_max(jset), zeta(jpgf, jset), la_min(jset), &
     744              :                                     ra, (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, na1, nb1, &
     745              :                                     rs_grid, &
     746         7056 :                                     radius=radius, ga_gb_function=GRID_FUNC_AB)
     747              :                               END IF
     748              :                            END DO
     749              :                         END DO
     750              :                      END IF
     751              :                   END DO
     752              :                END DO
     753          672 :                DEALLOCATE (map_it2)
     754              :                !
     755              :             END DO
     756              :          END DO
     757           96 :          DEALLOCATE (pab, work)
     758              :       END IF
     759              : 
     760          908 :       CALL pw_zero(lri_rho_g)
     761          908 :       CALL pw_zero(lri_rho_r)
     762              : 
     763         4480 :       DO igrid_level = 1, gridlevel_info%ngrid_levels
     764         3572 :          CALL pw_zero(mgrid_rspace(igrid_level))
     765              :          CALL transfer_rs2pw(rs=rs_rho(igrid_level), &
     766         4480 :                              pw=mgrid_rspace(igrid_level))
     767              :       END DO
     768              : 
     769         4480 :       DO igrid_level = 1, gridlevel_info%ngrid_levels
     770         3572 :          CALL pw_zero(mgrid_gspace(igrid_level))
     771              :          CALL pw_transfer(mgrid_rspace(igrid_level), &
     772         3572 :                           mgrid_gspace(igrid_level))
     773         4480 :          CALL pw_axpy(mgrid_gspace(igrid_level), lri_rho_g)
     774              :       END DO
     775          908 :       CALL pw_transfer(lri_rho_g, lri_rho_r)
     776          908 :       total_rho = pw_integrate_function(lri_rho_r, isign=-1)
     777              : 
     778              :       ! *** give back the multi-grids *** !
     779          908 :       CALL pw_pools_give_back_pws(pw_pools, mgrid_gspace)
     780          908 :       CALL pw_pools_give_back_pws(pw_pools, mgrid_rspace)
     781              : 
     782          908 :       CALL timestop(handle)
     783              : 
     784         4540 :    END SUBROUTINE calculate_lri_rho_elec
     785              : 
     786              :    #:for kind in ["r3d_rs", "c1d_gs"]
     787              : ! **************************************************************************************************
     788              : !> \brief computes the density of the core charges on the grid
     789              : !> \param rho_core ...
     790              : !> \param total_rho ...
     791              : !> \param qs_env ...
     792              : !> \param calpha ...
     793              : !> \param ccore ...
     794              : !> \param only_nopaw ...
     795              : ! **************************************************************************************************
     796         9454 :       SUBROUTINE calculate_rho_core_${kind}$ (rho_core, total_rho, qs_env, calpha, ccore, only_nopaw)
     797              : 
     798              :          TYPE(pw_${kind}$_type), INTENT(INOUT)                       :: rho_core
     799              :          REAL(KIND=dp), INTENT(OUT)                         :: total_rho
     800              :          TYPE(qs_environment_type), POINTER                 :: qs_env
     801              :          REAL(KIND=dp), DIMENSION(:), OPTIONAL              :: calpha, ccore
     802              :          LOGICAL, INTENT(IN), OPTIONAL                      :: only_nopaw
     803              : 
     804              :          CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_core'
     805              : 
     806              :          INTEGER                                            :: atom_a, handle, iatom, ikind, ithread, &
     807              :                                                                j, natom, npme, nthread, &
     808              :                                                                subpatch_pattern
     809         9454 :          INTEGER, DIMENSION(:), POINTER                     :: atom_list, cores
     810              :          LOGICAL                                            :: my_only_nopaw, paw_atom
     811              :          REAL(KIND=dp)                                      :: alpha, eps_rho_rspace, radius
     812              :          REAL(KIND=dp), DIMENSION(3)                        :: ra
     813         9454 :          REAL(KIND=dp), DIMENSION(:, :), POINTER            :: pab
     814         9454 :          TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     815              :          TYPE(cell_type), POINTER                           :: cell
     816              :          TYPE(dft_control_type), POINTER                    :: dft_control
     817         9454 :          TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     818              :          TYPE(pw_env_type), POINTER                         :: pw_env
     819              :          TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     820              :          TYPE(pw_r3d_rs_type)                                      :: rhoc_r
     821         9454 :          TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     822              :          TYPE(realspace_grid_type), POINTER                 :: rs_rho
     823              : 
     824         9454 :          CALL timeset(routineN, handle)
     825         9454 :          NULLIFY (cell, dft_control, pab, atomic_kind_set, qs_kind_set, particle_set, &
     826         9454 :                   atom_list, pw_env, rs_rho, auxbas_pw_pool, cores)
     827         9454 :          ALLOCATE (pab(1, 1))
     828              : 
     829         9454 :          my_only_nopaw = .FALSE.
     830         9454 :          IF (PRESENT(only_nopaw)) my_only_nopaw = only_nopaw
     831         9454 :          IF (PRESENT(calpha)) THEN
     832            2 :             CPASSERT(PRESENT(ccore))
     833              :          END IF
     834              : 
     835              :          CALL get_qs_env(qs_env=qs_env, &
     836              :                          atomic_kind_set=atomic_kind_set, &
     837              :                          qs_kind_set=qs_kind_set, &
     838              :                          cell=cell, &
     839              :                          dft_control=dft_control, &
     840              :                          particle_set=particle_set, &
     841         9454 :                          pw_env=pw_env)
     842              :          CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
     843         9454 :                          auxbas_pw_pool=auxbas_pw_pool)
     844              :          ! be careful in parallel nsmax is chosen with multigrid in mind!
     845         9454 :          CALL rs_grid_zero(rs_rho)
     846              : 
     847         9454 :          eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
     848              : 
     849        26393 :          DO ikind = 1, SIZE(atomic_kind_set)
     850        16939 :             CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
     851        16939 :             IF (PRESENT(calpha)) THEN
     852            4 :                alpha = calpha(ikind)
     853            4 :                pab(1, 1) = ccore(ikind)
     854              :             ELSE
     855              :                CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom, &
     856        16935 :                                 alpha_core_charge=alpha, ccore_charge=pab(1, 1))
     857              :             END IF
     858              : 
     859        16939 :             IF (my_only_nopaw .AND. paw_atom) CYCLE
     860        16781 :             IF (alpha == 0.0_dp .OR. pab(1, 1) == 0.0_dp) CYCLE
     861              : 
     862        16605 :             nthread = 1
     863        16605 :             ithread = 0
     864              : 
     865        16605 :             CALL reallocate(cores, 1, natom)
     866        16605 :             npme = 0
     867        53670 :             cores = 0
     868              : 
     869        53670 :             DO iatom = 1, natom
     870        53670 :                IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
     871              :                   ! replicated realspace grid, split the atoms up between procs
     872        36238 :                   IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
     873        18119 :                      npme = npme + 1
     874        18119 :                      cores(npme) = iatom
     875              :                   END IF
     876              :                ELSE
     877          827 :                   npme = npme + 1
     878          827 :                   cores(npme) = iatom
     879              :                END IF
     880              :             END DO
     881              : 
     882        42998 :             IF (npme .GT. 0) THEN
     883        32107 :                DO j = 1, npme
     884              : 
     885        18946 :                   iatom = cores(j)
     886        18946 :                   atom_a = atom_list(iatom)
     887        18946 :                   ra(:) = pbc(particle_set(atom_a)%r, cell)
     888        18946 :                   subpatch_pattern = 0
     889              :                   radius = exp_radius_very_extended(la_min=0, la_max=0, &
     890              :                                                     lb_min=0, lb_max=0, &
     891              :                                                     ra=ra, rb=ra, rp=ra, &
     892              :                                                     zetp=alpha, eps=eps_rho_rspace, &
     893              :                                                     pab=pab, o1=0, o2=0, &  ! without map_consistent
     894        18946 :                                                     prefactor=-1.0_dp, cutoff=0.0_dp)
     895              : 
     896              :                   CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
     897              :                                              (/0.0_dp, 0.0_dp, 0.0_dp/), -1.0_dp, pab, 0, 0, rs_rho, &
     898              :                                              radius=radius, ga_gb_function=GRID_FUNC_AB, &
     899        32107 :                                              use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern)
     900              : 
     901              :                END DO
     902              :             END IF
     903              : 
     904              :          END DO
     905              : 
     906         9454 :          IF (ASSOCIATED(cores)) THEN
     907         9446 :             DEALLOCATE (cores)
     908              :          END IF
     909         9454 :          DEALLOCATE (pab)
     910              : 
     911         9454 :          CALL auxbas_pw_pool%create_pw(rhoc_r)
     912              : 
     913         9454 :          CALL transfer_rs2pw(rs_rho, rhoc_r)
     914              : 
     915         9454 :          total_rho = pw_integrate_function(rhoc_r, isign=-1)
     916              : 
     917         9454 :          CALL pw_transfer(rhoc_r, rho_core)
     918              : 
     919         9454 :          CALL auxbas_pw_pool%give_back_pw(rhoc_r)
     920              : 
     921         9454 :          CALL timestop(handle)
     922              : 
     923         9454 :       END SUBROUTINE calculate_rho_core_${kind}$
     924              :    #:endfor
     925              : 
     926              : ! *****************************************************************************
     927              : !> \brief Computes the derivative of the density of the core charges with
     928              : !>        respect to the nuclear coordinates on the grid.
     929              : !> \param drho_core The resulting density derivative
     930              : !> \param qs_env ...
     931              : !> \param beta Derivative direction
     932              : !> \param lambda Atom index
     933              : !> \note SL November 2014, ED 2021
     934              : ! **************************************************************************************************
     935          216 :    SUBROUTINE calculate_drho_core(drho_core, qs_env, beta, lambda)
     936              : 
     937              :       TYPE(pw_c1d_gs_type), INTENT(INOUT)                       :: drho_core
     938              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     939              :       INTEGER, INTENT(IN)                                :: beta, lambda
     940              : 
     941              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_drho_core'
     942              : 
     943              :       INTEGER                                            :: atom_a, dabqadb_func, handle, iatom, &
     944              :                                                             ikind, ithread, j, natom, npme, &
     945              :                                                             nthread, subpatch_pattern
     946          216 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list, cores
     947              :       REAL(KIND=dp)                                      :: alpha, eps_rho_rspace, radius
     948              :       REAL(KIND=dp), DIMENSION(3)                        :: ra
     949          216 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: pab
     950          216 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     951              :       TYPE(cell_type), POINTER                           :: cell
     952              :       TYPE(dft_control_type), POINTER                    :: dft_control
     953          216 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     954              :       TYPE(pw_env_type), POINTER                         :: pw_env
     955              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     956              :       TYPE(pw_r3d_rs_type)                                      :: rhoc_r
     957          216 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     958              :       TYPE(realspace_grid_type), POINTER                 :: rs_rho
     959              : 
     960          216 :       CALL timeset(routineN, handle)
     961          216 :       NULLIFY (cell, dft_control, pab, atomic_kind_set, qs_kind_set, particle_set, &
     962          216 :                atom_list, pw_env, rs_rho, auxbas_pw_pool, cores)
     963          216 :       ALLOCATE (pab(1, 1))
     964              : 
     965              :       CALL get_qs_env(qs_env=qs_env, &
     966              :                       atomic_kind_set=atomic_kind_set, &
     967              :                       qs_kind_set=qs_kind_set, &
     968              :                       cell=cell, &
     969              :                       dft_control=dft_control, &
     970              :                       particle_set=particle_set, &
     971          216 :                       pw_env=pw_env)
     972              :       CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
     973          216 :                       auxbas_pw_pool=auxbas_pw_pool)
     974              :       ! be careful in parallel nsmax is chosen with multigrid in mind!
     975          216 :       CALL rs_grid_zero(rs_rho)
     976              : 
     977          216 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
     978              : 
     979          288 :       SELECT CASE (beta)
     980              :       CASE (1)
     981           72 :          dabqadb_func = GRID_FUNC_CORE_X
     982              :       CASE (2)
     983           72 :          dabqadb_func = GRID_FUNC_CORE_Y
     984              :       CASE (3)
     985           72 :          dabqadb_func = GRID_FUNC_CORE_Z
     986              :       CASE DEFAULT
     987          216 :          CPABORT("invalid beta")
     988              :       END SELECT
     989          648 :       DO ikind = 1, SIZE(atomic_kind_set)
     990          432 :          CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
     991              :          CALL get_qs_kind(qs_kind_set(ikind), &
     992          432 :                           alpha_core_charge=alpha, ccore_charge=pab(1, 1))
     993              : 
     994          432 :          IF (alpha == 0.0_dp .OR. pab(1, 1) == 0.0_dp) CYCLE
     995              : 
     996          432 :          nthread = 1
     997          432 :          ithread = 0
     998              : 
     999          432 :          CALL reallocate(cores, 1, natom)
    1000          432 :          npme = 0
    1001         1080 :          cores = 0
    1002              : 
    1003         1080 :          DO iatom = 1, natom
    1004         1080 :             IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
    1005              :                ! replicated realspace grid, split the atoms up between procs
    1006          648 :                IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
    1007          324 :                   npme = npme + 1
    1008          324 :                   cores(npme) = iatom
    1009              :                END IF
    1010              :             ELSE
    1011            0 :                npme = npme + 1
    1012            0 :                cores(npme) = iatom
    1013              :             END IF
    1014              :          END DO
    1015              : 
    1016         1080 :          IF (npme .GT. 0) THEN
    1017          648 :             DO j = 1, npme
    1018              : 
    1019          324 :                iatom = cores(j)
    1020          324 :                atom_a = atom_list(iatom)
    1021          324 :                IF (atom_a /= lambda) CYCLE
    1022          108 :                ra(:) = pbc(particle_set(atom_a)%r, cell)
    1023          108 :                subpatch_pattern = 0
    1024              :                radius = exp_radius_very_extended(la_min=0, la_max=0, &
    1025              :                                                  lb_min=0, lb_max=0, &
    1026              :                                                  ra=ra, rb=ra, rp=ra, &
    1027              :                                                  zetp=alpha, eps=eps_rho_rspace, &
    1028              :                                                  pab=pab, o1=0, o2=0, &  ! without map_consistent
    1029          108 :                                                  prefactor=-1.0_dp, cutoff=0.0_dp)
    1030              : 
    1031              :                CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
    1032              :                                           (/0.0_dp, 0.0_dp, 0.0_dp/), -1.0_dp, pab, 0, 0, rs_rho, &
    1033              :                                           radius=radius, ga_gb_function=dabqadb_func, &
    1034          648 :                                           use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern)
    1035              : 
    1036              :             END DO
    1037              :          END IF
    1038              : 
    1039              :       END DO
    1040              : 
    1041          216 :       IF (ASSOCIATED(cores)) THEN
    1042          216 :          DEALLOCATE (cores)
    1043              :       END IF
    1044          216 :       DEALLOCATE (pab)
    1045              : 
    1046          216 :       CALL auxbas_pw_pool%create_pw(rhoc_r)
    1047              : 
    1048          216 :       CALL transfer_rs2pw(rs_rho, rhoc_r)
    1049              : 
    1050          216 :       CALL pw_transfer(rhoc_r, drho_core)
    1051              : 
    1052          216 :       CALL auxbas_pw_pool%give_back_pw(rhoc_r)
    1053              : 
    1054          216 :       CALL timestop(handle)
    1055              : 
    1056          216 :    END SUBROUTINE calculate_drho_core
    1057              : 
    1058              : ! **************************************************************************************************
    1059              : !> \brief collocate a single Gaussian on the grid
    1060              : !> \param rho_gb charge density generated by a single gaussian
    1061              : !> \param qs_env qs environment
    1062              : !> \param iatom_in atom index
    1063              : !> \par History
    1064              : !>        12.2011 created
    1065              : !> \author Dorothea Golze
    1066              : ! **************************************************************************************************
    1067            4 :    SUBROUTINE calculate_rho_single_gaussian(rho_gb, qs_env, iatom_in)
    1068              : 
    1069              :       TYPE(pw_c1d_gs_type), INTENT(INOUT)                       :: rho_gb
    1070              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1071              :       INTEGER, INTENT(IN)                                :: iatom_in
    1072              : 
    1073              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_single_gaussian'
    1074              : 
    1075              :       INTEGER                                            :: atom_a, handle, iatom, npme, &
    1076              :                                                             subpatch_pattern
    1077              :       REAL(KIND=dp)                                      :: eps_rho_rspace, radius
    1078              :       REAL(KIND=dp), DIMENSION(3)                        :: ra
    1079            4 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: pab
    1080              :       TYPE(cell_type), POINTER                           :: cell
    1081              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1082              :       TYPE(pw_env_type), POINTER                         :: pw_env
    1083              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    1084              :       TYPE(pw_r3d_rs_type)                                      :: rhoc_r
    1085              :       TYPE(realspace_grid_type), POINTER                 :: rs_rho
    1086              : 
    1087            4 :       CALL timeset(routineN, handle)
    1088            4 :       NULLIFY (cell, dft_control, pab, pw_env, rs_rho, auxbas_pw_pool)
    1089              : 
    1090            4 :       ALLOCATE (pab(1, 1))
    1091              : 
    1092              :       CALL get_qs_env(qs_env=qs_env, &
    1093              :                       cell=cell, &
    1094              :                       dft_control=dft_control, &
    1095            4 :                       pw_env=pw_env)
    1096              :       CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
    1097            4 :                       auxbas_pw_pool=auxbas_pw_pool)
    1098            4 :       CALL rs_grid_zero(rs_rho)
    1099              : 
    1100            4 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
    1101            4 :       pab(1, 1) = 1.0_dp
    1102            4 :       iatom = iatom_in
    1103              : 
    1104            4 :       npme = 0
    1105              : 
    1106            4 :       IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
    1107            4 :          IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
    1108              :             npme = npme + 1
    1109              :          END IF
    1110              :       ELSE
    1111              :          npme = npme + 1
    1112              :       END IF
    1113              : 
    1114              :       IF (npme .GT. 0) THEN
    1115            2 :          atom_a = qs_env%qmmm_env_qm%image_charge_pot%image_mm_list(iatom)
    1116            2 :          ra(:) = pbc(qs_env%qmmm_env_qm%image_charge_pot%particles_all(atom_a)%r, cell)
    1117            2 :          subpatch_pattern = 0
    1118              :          radius = exp_radius_very_extended(la_min=0, la_max=0, &
    1119              :                                            lb_min=0, lb_max=0, &
    1120              :                                            ra=ra, rb=ra, rp=ra, &
    1121              :                                            zetp=qs_env%qmmm_env_qm%image_charge_pot%eta, &
    1122              :                                            eps=eps_rho_rspace, &
    1123              :                                            pab=pab, o1=0, o2=0, &  ! without map_consistent
    1124            2 :                                            prefactor=1.0_dp, cutoff=0.0_dp)
    1125              : 
    1126              :          CALL collocate_pgf_product(0, qs_env%qmmm_env_qm%image_charge_pot%eta, &
    1127              :                                     0, 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, rs_rho, &
    1128              :                                     radius=radius, ga_gb_function=GRID_FUNC_AB, &
    1129            2 :                                     use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern)
    1130              :       END IF
    1131              : 
    1132            4 :       DEALLOCATE (pab)
    1133              : 
    1134            4 :       CALL auxbas_pw_pool%create_pw(rhoc_r)
    1135              : 
    1136            4 :       CALL transfer_rs2pw(rs_rho, rhoc_r)
    1137              : 
    1138            4 :       CALL pw_transfer(rhoc_r, rho_gb)
    1139              : 
    1140            4 :       CALL auxbas_pw_pool%give_back_pw(rhoc_r)
    1141              : 
    1142            4 :       CALL timestop(handle)
    1143              : 
    1144            4 :    END SUBROUTINE calculate_rho_single_gaussian
    1145              : 
    1146              : ! **************************************************************************************************
    1147              : !> \brief computes the image charge density on the grid (including coeffcients)
    1148              : !> \param rho_metal image charge density
    1149              : !> \param coeff expansion coefficients of the image charge density, i.e.
    1150              : !>        rho_metal=sum_a c_a*g_a
    1151              : !> \param total_rho_metal total induced image charge density
    1152              : !> \param qs_env qs environment
    1153              : !> \par History
    1154              : !>        01.2012 created
    1155              : !> \author Dorothea Golze
    1156              : ! **************************************************************************************************
    1157           90 :    SUBROUTINE calculate_rho_metal(rho_metal, coeff, total_rho_metal, qs_env)
    1158              : 
    1159              :       TYPE(pw_c1d_gs_type), INTENT(INOUT)                       :: rho_metal
    1160              :       REAL(KIND=dp), DIMENSION(:), POINTER               :: coeff
    1161              :       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: total_rho_metal
    1162              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1163              : 
    1164              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_metal'
    1165              : 
    1166              :       INTEGER                                            :: atom_a, handle, iatom, j, natom, npme, &
    1167              :                                                             subpatch_pattern
    1168           90 :       INTEGER, DIMENSION(:), POINTER                     :: cores
    1169              :       REAL(KIND=dp)                                      :: eps_rho_rspace, radius
    1170              :       REAL(KIND=dp), DIMENSION(3)                        :: ra
    1171           90 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: pab
    1172              :       TYPE(cell_type), POINTER                           :: cell
    1173              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1174              :       TYPE(pw_env_type), POINTER                         :: pw_env
    1175              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    1176              :       TYPE(pw_r3d_rs_type)                                      :: rhoc_r
    1177              :       TYPE(realspace_grid_type), POINTER                 :: rs_rho
    1178              : 
    1179           90 :       CALL timeset(routineN, handle)
    1180              : 
    1181           90 :       NULLIFY (cell, dft_control, pab, pw_env, rs_rho, auxbas_pw_pool, cores)
    1182              : 
    1183           90 :       ALLOCATE (pab(1, 1))
    1184              : 
    1185              :       CALL get_qs_env(qs_env=qs_env, &
    1186              :                       cell=cell, &
    1187              :                       dft_control=dft_control, &
    1188           90 :                       pw_env=pw_env)
    1189              :       CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
    1190           90 :                       auxbas_pw_pool=auxbas_pw_pool)
    1191           90 :       CALL rs_grid_zero(rs_rho)
    1192              : 
    1193           90 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
    1194           90 :       pab(1, 1) = 1.0_dp
    1195              : 
    1196           90 :       natom = SIZE(qs_env%qmmm_env_qm%image_charge_pot%image_mm_list)
    1197              : 
    1198           90 :       CALL reallocate(cores, 1, natom)
    1199           90 :       npme = 0
    1200          270 :       cores = 0
    1201              : 
    1202          270 :       DO iatom = 1, natom
    1203          270 :          IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
    1204          180 :             IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
    1205           90 :                npme = npme + 1
    1206           90 :                cores(npme) = iatom
    1207              :             END IF
    1208              :          ELSE
    1209            0 :             npme = npme + 1
    1210            0 :             cores(npme) = iatom
    1211              :          END IF
    1212              :       END DO
    1213              : 
    1214           90 :       IF (npme .GT. 0) THEN
    1215          180 :          DO j = 1, npme
    1216           90 :             iatom = cores(j)
    1217           90 :             atom_a = qs_env%qmmm_env_qm%image_charge_pot%image_mm_list(iatom)
    1218           90 :             ra(:) = pbc(qs_env%qmmm_env_qm%image_charge_pot%particles_all(atom_a)%r, cell)
    1219           90 :             subpatch_pattern = 0
    1220              :             radius = exp_radius_very_extended(la_min=0, la_max=0, &
    1221              :                                               lb_min=0, lb_max=0, &
    1222              :                                               ra=ra, rb=ra, rp=ra, &
    1223              :                                               zetp=qs_env%qmmm_env_qm%image_charge_pot%eta, &
    1224              :                                               eps=eps_rho_rspace, &
    1225              :                                               pab=pab, o1=0, o2=0, &  ! without map_consistent
    1226           90 :                                               prefactor=coeff(iatom), cutoff=0.0_dp)
    1227              : 
    1228              :             CALL collocate_pgf_product( &
    1229              :                0, qs_env%qmmm_env_qm%image_charge_pot%eta, &
    1230              :                0, 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), coeff(iatom), pab, 0, 0, rs_rho, &
    1231              :                radius=radius, ga_gb_function=GRID_FUNC_AB, &
    1232          180 :                use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern)
    1233              :          END DO
    1234              :       END IF
    1235              : 
    1236           90 :       DEALLOCATE (pab, cores)
    1237              : 
    1238           90 :       CALL auxbas_pw_pool%create_pw(rhoc_r)
    1239              : 
    1240           90 :       CALL transfer_rs2pw(rs_rho, rhoc_r)
    1241              : 
    1242           90 :       IF (PRESENT(total_rho_metal)) &
    1243              :          !minus sign: account for the fact that rho_metal has opposite sign
    1244           90 :          total_rho_metal = pw_integrate_function(rhoc_r, isign=-1)
    1245              : 
    1246           90 :       CALL pw_transfer(rhoc_r, rho_metal)
    1247           90 :       CALL auxbas_pw_pool%give_back_pw(rhoc_r)
    1248              : 
    1249           90 :       CALL timestop(handle)
    1250              : 
    1251           90 :    END SUBROUTINE calculate_rho_metal
    1252              : 
    1253              : ! **************************************************************************************************
    1254              : !> \brief collocate a single Gaussian on the grid for periodic RESP fitting
    1255              : !> \param rho_gb charge density generated by a single gaussian
    1256              : !> \param qs_env qs environment
    1257              : !> \param eta width of single Gaussian
    1258              : !> \param iatom_in atom index
    1259              : !> \par History
    1260              : !>        06.2012 created
    1261              : !> \author Dorothea Golze
    1262              : ! **************************************************************************************************
    1263           66 :    SUBROUTINE calculate_rho_resp_single(rho_gb, qs_env, eta, iatom_in)
    1264              : 
    1265              :       TYPE(pw_c1d_gs_type), INTENT(INOUT)                       :: rho_gb
    1266              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1267              :       REAL(KIND=dp), INTENT(IN)                          :: eta
    1268              :       INTEGER, INTENT(IN)                                :: iatom_in
    1269              : 
    1270              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_resp_single'
    1271              : 
    1272              :       INTEGER                                            :: handle, iatom, npme, subpatch_pattern
    1273              :       REAL(KIND=dp)                                      :: eps_rho_rspace, radius
    1274              :       REAL(KIND=dp), DIMENSION(3)                        :: ra
    1275           66 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: pab
    1276              :       TYPE(cell_type), POINTER                           :: cell
    1277              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1278           66 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1279              :       TYPE(pw_env_type), POINTER                         :: pw_env
    1280              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    1281              :       TYPE(pw_r3d_rs_type)                                      :: rhoc_r
    1282              :       TYPE(realspace_grid_type), POINTER                 :: rs_rho
    1283              : 
    1284           66 :       CALL timeset(routineN, handle)
    1285           66 :       NULLIFY (cell, dft_control, pab, pw_env, rs_rho, auxbas_pw_pool, &
    1286           66 :                particle_set)
    1287              : 
    1288           66 :       ALLOCATE (pab(1, 1))
    1289              : 
    1290              :       CALL get_qs_env(qs_env=qs_env, &
    1291              :                       cell=cell, &
    1292              :                       dft_control=dft_control, &
    1293              :                       particle_set=particle_set, &
    1294           66 :                       pw_env=pw_env)
    1295              :       CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
    1296           66 :                       auxbas_pw_pool=auxbas_pw_pool)
    1297           66 :       CALL rs_grid_zero(rs_rho)
    1298              : 
    1299           66 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
    1300           66 :       pab(1, 1) = 1.0_dp
    1301           66 :       iatom = iatom_in
    1302              : 
    1303           66 :       npme = 0
    1304              : 
    1305           66 :       IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
    1306           66 :          IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
    1307              :             npme = npme + 1
    1308              :          END IF
    1309              :       ELSE
    1310              :          npme = npme + 1
    1311              :       END IF
    1312              : 
    1313              :       IF (npme .GT. 0) THEN
    1314           33 :          ra(:) = pbc(particle_set(iatom)%r, cell)
    1315           33 :          subpatch_pattern = 0
    1316              :          radius = exp_radius_very_extended(la_min=0, la_max=0, &
    1317              :                                            lb_min=0, lb_max=0, &
    1318              :                                            ra=ra, rb=ra, rp=ra, &
    1319              :                                            zetp=eta, eps=eps_rho_rspace, &
    1320              :                                            pab=pab, o1=0, o2=0, &  ! without map_consistent
    1321           33 :                                            prefactor=1.0_dp, cutoff=0.0_dp)
    1322              : 
    1323              :          CALL collocate_pgf_product(0, eta, 0, 0, 0.0_dp, 0, ra, &
    1324              :                                     (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, rs_rho, &
    1325              :                                     radius=radius, ga_gb_function=GRID_FUNC_AB, &
    1326           33 :                                     use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern)
    1327              :       END IF
    1328              : 
    1329           66 :       DEALLOCATE (pab)
    1330              : 
    1331           66 :       CALL auxbas_pw_pool%create_pw(rhoc_r)
    1332              : 
    1333           66 :       CALL transfer_rs2pw(rs_rho, rhoc_r)
    1334              : 
    1335           66 :       CALL pw_transfer(rhoc_r, rho_gb)
    1336              : 
    1337           66 :       CALL auxbas_pw_pool%give_back_pw(rhoc_r)
    1338              : 
    1339           66 :       CALL timestop(handle)
    1340              : 
    1341           66 :    END SUBROUTINE calculate_rho_resp_single
    1342              : 
    1343              :    #:for kind in ["r3d_rs", "c1d_gs"]
    1344              : ! **************************************************************************************************
    1345              : !> \brief computes the RESP charge density on a grid based on the RESP charges
    1346              : !> \param rho_resp RESP charge density
    1347              : !> \param coeff RESP charges, take care of normalization factor
    1348              : !>        (eta/pi)**1.5 later
    1349              : !> \param natom number of atoms
    1350              : !> \param eta width of single Gaussian
    1351              : !> \param qs_env qs environment
    1352              : !> \par History
    1353              : !>        01.2012 created
    1354              : !> \author Dorothea Golze
    1355              : ! **************************************************************************************************
    1356           24 :       SUBROUTINE calculate_rho_resp_all_${kind}$ (rho_resp, coeff, natom, eta, qs_env)
    1357              : 
    1358              :          TYPE(pw_${kind}$_type), INTENT(INOUT)                       :: rho_resp
    1359              :          REAL(KIND=dp), DIMENSION(:), POINTER               :: coeff
    1360              :          INTEGER, INTENT(IN)                                :: natom
    1361              :          REAL(KIND=dp), INTENT(IN)                          :: eta
    1362              :          TYPE(qs_environment_type), POINTER                 :: qs_env
    1363              : 
    1364              :          CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_resp_all'
    1365              : 
    1366              :          INTEGER                                            :: handle, iatom, j, npme, subpatch_pattern
    1367           24 :          INTEGER, DIMENSION(:), POINTER                     :: cores
    1368              :          REAL(KIND=dp)                                      :: eps_rho_rspace, radius
    1369              :          REAL(KIND=dp), DIMENSION(3)                        :: ra
    1370           24 :          REAL(KIND=dp), DIMENSION(:, :), POINTER            :: pab
    1371              :          TYPE(cell_type), POINTER                           :: cell
    1372              :          TYPE(dft_control_type), POINTER                    :: dft_control
    1373           24 :          TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1374              :          TYPE(pw_env_type), POINTER                         :: pw_env
    1375              :          TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    1376              :          TYPE(pw_r3d_rs_type)                                      :: rhoc_r
    1377              :          TYPE(realspace_grid_type), POINTER                 :: rs_rho
    1378              : 
    1379           24 :          CALL timeset(routineN, handle)
    1380              : 
    1381           24 :          NULLIFY (cell, cores, dft_control, pab, pw_env, rs_rho, auxbas_pw_pool, &
    1382           24 :                   particle_set)
    1383              : 
    1384           24 :          ALLOCATE (pab(1, 1))
    1385              : 
    1386              :          CALL get_qs_env(qs_env=qs_env, &
    1387              :                          cell=cell, &
    1388              :                          dft_control=dft_control, &
    1389              :                          particle_set=particle_set, &
    1390           24 :                          pw_env=pw_env)
    1391              :          CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
    1392           24 :                          auxbas_pw_pool=auxbas_pw_pool)
    1393           24 :          CALL rs_grid_zero(rs_rho)
    1394              : 
    1395           24 :          eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
    1396           24 :          pab(1, 1) = 1.0_dp
    1397              : 
    1398           24 :          CALL reallocate(cores, 1, natom)
    1399           24 :          npme = 0
    1400          142 :          cores = 0
    1401              : 
    1402          142 :          DO iatom = 1, natom
    1403          142 :             IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
    1404          118 :                IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
    1405           59 :                   npme = npme + 1
    1406           59 :                   cores(npme) = iatom
    1407              :                END IF
    1408              :             ELSE
    1409            0 :                npme = npme + 1
    1410            0 :                cores(npme) = iatom
    1411              :             END IF
    1412              :          END DO
    1413              : 
    1414           24 :          IF (npme .GT. 0) THEN
    1415           83 :             DO j = 1, npme
    1416           59 :                iatom = cores(j)
    1417           59 :                ra(:) = pbc(particle_set(iatom)%r, cell)
    1418           59 :                subpatch_pattern = 0
    1419              :                radius = exp_radius_very_extended(la_min=0, la_max=0, &
    1420              :                                                  lb_min=0, lb_max=0, &
    1421              :                                                  ra=ra, rb=ra, rp=ra, &
    1422              :                                                  zetp=eta, eps=eps_rho_rspace, &
    1423              :                                                  pab=pab, o1=0, o2=0, &  ! without map_consistent
    1424           59 :                                                  prefactor=coeff(iatom), cutoff=0.0_dp)
    1425              : 
    1426              :                CALL collocate_pgf_product( &
    1427              :                   0, eta, &
    1428              :                   0, 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), coeff(iatom), pab, 0, 0, rs_rho, &
    1429              :                   radius=radius, ga_gb_function=GRID_FUNC_AB, &
    1430           83 :                   use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern)
    1431              :             END DO
    1432              :          END IF
    1433              : 
    1434           24 :          DEALLOCATE (pab, cores)
    1435              : 
    1436           24 :          CALL auxbas_pw_pool%create_pw(rhoc_r)
    1437              : 
    1438           24 :          CALL transfer_rs2pw(rs_rho, rhoc_r)
    1439              : 
    1440           24 :          CALL pw_transfer(rhoc_r, rho_resp)
    1441           24 :          CALL auxbas_pw_pool%give_back_pw(rhoc_r)
    1442              : 
    1443           24 :          CALL timestop(handle)
    1444              : 
    1445           24 :       END SUBROUTINE calculate_rho_resp_all_${kind}$
    1446              :    #:endfor
    1447              : 
    1448              : ! **************************************************************************************************
    1449              : !> \brief computes the density corresponding to a given density matrix on the grid
    1450              : !> \param matrix_p ...
    1451              : !> \param matrix_p_kp ...
    1452              : !> \param rho ...
    1453              : !> \param rho_gspace ...
    1454              : !> \param total_rho ...
    1455              : !> \param ks_env ...
    1456              : !> \param soft_valid ...
    1457              : !> \param compute_tau ...
    1458              : !> \param compute_grad ...
    1459              : !> \param basis_type ...
    1460              : !> \param der_type ...
    1461              : !> \param idir ...
    1462              : !> \param task_list_external ...
    1463              : !> \param pw_env_external ...
    1464              : !> \par History
    1465              : !>      IAB (15-Feb-2010): Added OpenMP parallelisation to task loop
    1466              : !>                         (c) The Numerical Algorithms Group (NAG) Ltd, 2010 on behalf of the HECToR project
    1467              : !>      Anything that is not the default ORB basis_type requires an external_task_list 12.2019, (A.Bussy)
    1468              : !>      Ole Schuett (2020): Migrated to C, see grid_api.F
    1469              : !> \note
    1470              : !>      both rho and rho_gspace contain the new rho
    1471              : !>      (in real and g-space respectively)
    1472              : ! **************************************************************************************************
    1473       201618 :    SUBROUTINE calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, &
    1474              :                                  ks_env, soft_valid, compute_tau, compute_grad, &
    1475              :                                  basis_type, der_type, idir, task_list_external, pw_env_external)
    1476              : 
    1477              :       TYPE(dbcsr_type), OPTIONAL, TARGET                 :: matrix_p
    1478              :       TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
    1479              :          POINTER                                         :: matrix_p_kp
    1480              :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                       :: rho
    1481              :       TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_gspace
    1482              :       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: total_rho
    1483              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
    1484              :       LOGICAL, INTENT(IN), OPTIONAL                      :: soft_valid, compute_tau, compute_grad
    1485              :       CHARACTER(LEN=*), INTENT(IN), OPTIONAL             :: basis_type
    1486              :       INTEGER, INTENT(IN), OPTIONAL                      :: der_type, idir
    1487              :       TYPE(task_list_type), OPTIONAL, POINTER            :: task_list_external
    1488              :       TYPE(pw_env_type), OPTIONAL, POINTER               :: pw_env_external
    1489              : 
    1490              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_elec'
    1491              : 
    1492              :       CHARACTER(LEN=default_string_length)               :: my_basis_type
    1493              :       INTEGER                                            :: ga_gb_function, handle, ilevel, img, &
    1494              :                                                             nimages, nlevels
    1495              :       LOGICAL                                            :: any_distributed, my_compute_grad, &
    1496              :                                                             my_compute_tau, my_soft_valid
    1497       201618 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_images
    1498              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1499              :       TYPE(mp_comm_type)                                 :: group
    1500              :       TYPE(pw_env_type), POINTER                         :: pw_env
    1501       201618 :       TYPE(realspace_grid_type), DIMENSION(:), POINTER   :: rs_rho
    1502              :       TYPE(task_list_type), POINTER                      :: task_list
    1503              : 
    1504       201618 :       CALL timeset(routineN, handle)
    1505              : 
    1506       201618 :       NULLIFY (matrix_images, dft_control, pw_env, rs_rho, task_list)
    1507              : 
    1508              :       ! Figure out which function to collocate.
    1509       201618 :       my_compute_tau = .FALSE.
    1510       201618 :       IF (PRESENT(compute_tau)) my_compute_tau = compute_tau
    1511       201618 :       my_compute_grad = .FALSE.
    1512       201618 :       IF (PRESENT(compute_grad)) my_compute_grad = compute_grad
    1513       201618 :       IF (PRESENT(der_type)) THEN
    1514           84 :          SELECT CASE (der_type)
    1515              :          CASE (orb_s)
    1516           36 :             ga_gb_function = GRID_FUNC_AB
    1517              :          CASE (orb_px)
    1518            0 :             ga_gb_function = GRID_FUNC_DX
    1519              :          CASE (orb_py)
    1520            0 :             ga_gb_function = GRID_FUNC_DY
    1521              :          CASE (orb_pz)
    1522           12 :             ga_gb_function = GRID_FUNC_DZ
    1523              :          CASE (orb_dxy)
    1524            0 :             ga_gb_function = GRID_FUNC_DXDY
    1525              :          CASE (orb_dyz)
    1526            0 :             ga_gb_function = GRID_FUNC_DYDZ
    1527              :          CASE (orb_dzx)
    1528            0 :             ga_gb_function = GRID_FUNC_DZDX
    1529              :          CASE (orb_dx2)
    1530            0 :             ga_gb_function = GRID_FUNC_DXDX
    1531              :          CASE (orb_dy2)
    1532            0 :             ga_gb_function = GRID_FUNC_DYDY
    1533              :          CASE (orb_dz2)
    1534            0 :             ga_gb_function = GRID_FUNC_DZDZ
    1535              :          CASE DEFAULT
    1536           48 :             CPABORT("Unknown der_type")
    1537              :          END SELECT
    1538       201570 :       ELSE IF (my_compute_tau) THEN
    1539         4536 :          ga_gb_function = GRID_FUNC_DADB
    1540       197034 :       ELSE IF (my_compute_grad) THEN
    1541          318 :          CPASSERT(PRESENT(idir))
    1542          424 :          SELECT CASE (idir)
    1543              :          CASE (1)
    1544          106 :             ga_gb_function = GRID_FUNC_DABpADB_X
    1545              :          CASE (2)
    1546          106 :             ga_gb_function = GRID_FUNC_DABpADB_Y
    1547              :          CASE (3)
    1548          106 :             ga_gb_function = GRID_FUNC_DABpADB_Z
    1549              :          CASE DEFAULT
    1550          318 :             CPABORT("invalid idir")
    1551              :          END SELECT
    1552              :       ELSE
    1553       196716 :          ga_gb_function = GRID_FUNC_AB
    1554              :       END IF
    1555              : 
    1556              :       ! Figure out which basis_type to use.
    1557       201618 :       my_basis_type = "ORB"  ! by default, the full density is calculated
    1558       201618 :       IF (PRESENT(basis_type)) my_basis_type = basis_type
    1559       201618 :       CPASSERT(my_basis_type == "ORB" .OR. PRESENT(task_list_external))
    1560              : 
    1561              :       ! Figure out which task_list to use.
    1562       201618 :       my_soft_valid = .FALSE.
    1563       201618 :       IF (PRESENT(soft_valid)) my_soft_valid = soft_valid
    1564       201618 :       IF (PRESENT(task_list_external)) THEN
    1565        38524 :          task_list => task_list_external
    1566       163094 :       ELSEIF (my_soft_valid) THEN
    1567        24058 :          CALL get_ks_env(ks_env, task_list_soft=task_list)
    1568              :       ELSE
    1569       139036 :          CALL get_ks_env(ks_env, task_list=task_list)
    1570              :       END IF
    1571       201618 :       CPASSERT(ASSOCIATED(task_list))
    1572              : 
    1573              :       ! Figure out which pw_env to use.
    1574       201618 :       IF (PRESENT(pw_env_external)) THEN
    1575        20792 :          pw_env => pw_env_external
    1576              :       ELSE
    1577       180826 :          CALL get_ks_env(ks_env, pw_env=pw_env)
    1578              :       END IF
    1579       201618 :       CPASSERT(ASSOCIATED(pw_env))
    1580              : 
    1581              :       ! Get grids.
    1582       201618 :       CALL pw_env_get(pw_env, rs_grids=rs_rho)
    1583       201618 :       nlevels = SIZE(rs_rho)
    1584       201618 :       group = rs_rho(1)%desc%group
    1585              : 
    1586              :       ! Check if any of the grids is distributed.
    1587       201618 :       any_distributed = .FALSE.
    1588       999350 :       DO ilevel = 1, nlevels
    1589      1796174 :          any_distributed = any_distributed .OR. rs_rho(ilevel)%desc%distributed
    1590              :       END DO
    1591              : 
    1592              :       ! Gather all matrix images in a single array.
    1593       201618 :       CALL get_ks_env(ks_env, dft_control=dft_control)
    1594       201618 :       nimages = dft_control%nimages
    1595       932852 :       ALLOCATE (matrix_images(nimages))
    1596       201618 :       IF (PRESENT(matrix_p_kp)) THEN
    1597       169974 :          CPASSERT(.NOT. PRESENT(matrix_p))
    1598       466328 :          DO img = 1, nimages
    1599       466328 :             matrix_images(img)%matrix => matrix_p_kp(img)%matrix
    1600              :          END DO
    1601              :       ELSE
    1602        31644 :          CPASSERT(PRESENT(matrix_p) .AND. nimages == 1)
    1603        31644 :          matrix_images(1)%matrix => matrix_p
    1604              :       END IF
    1605              : 
    1606              :       ! Distribute matrix blocks.
    1607       201618 :       IF (any_distributed) THEN
    1608          230 :          CALL rs_scatter_matrices(matrix_images, task_list%pab_buffer, task_list, group)
    1609              :       ELSE
    1610       201388 :          CALL rs_copy_to_buffer(matrix_images, task_list%pab_buffer, task_list)
    1611              :       END IF
    1612       201618 :       DEALLOCATE (matrix_images)
    1613              : 
    1614              :       ! Map all tasks onto the grids
    1615              :       CALL grid_collocate_task_list(task_list=task_list%grid_task_list, &
    1616              :                                     ga_gb_function=ga_gb_function, &
    1617              :                                     pab_blocks=task_list%pab_buffer, &
    1618       201618 :                                     rs_grids=rs_rho)
    1619              : 
    1620              :       ! Merge realspace multi-grids into single planewave grid.
    1621       201618 :       CALL density_rs2pw(pw_env, rs_rho, rho, rho_gspace)
    1622       201618 :       IF (PRESENT(total_rho)) total_rho = pw_integrate_function(rho, isign=-1)
    1623              : 
    1624       201618 :       CALL timestop(handle)
    1625              : 
    1626       201618 :    END SUBROUTINE calculate_rho_elec
    1627              : 
    1628              : ! **************************************************************************************************
    1629              : !> \brief computes the gradient of the density corresponding to a given
    1630              : !>        density matrix on the grid
    1631              : !> \param matrix_p ...
    1632              : !> \param matrix_p_kp ...
    1633              : !> \param drho ...
    1634              : !> \param drho_gspace ...
    1635              : !> \param qs_env ...
    1636              : !> \param soft_valid ...
    1637              : !> \param basis_type ...
    1638              : !> \note  this is an alternative to calculate the gradient through FFTs
    1639              : ! **************************************************************************************************
    1640            0 :    SUBROUTINE calculate_drho_elec(matrix_p, matrix_p_kp, drho, drho_gspace, qs_env, &
    1641              :                                   soft_valid, basis_type)
    1642              : 
    1643              :       TYPE(dbcsr_type), OPTIONAL, TARGET                 :: matrix_p
    1644              :       TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
    1645              :          POINTER                                         :: matrix_p_kp
    1646              :       TYPE(pw_r3d_rs_type), DIMENSION(3), INTENT(INOUT)         :: drho
    1647              :       TYPE(pw_c1d_gs_type), DIMENSION(3), INTENT(INOUT) :: drho_gspace
    1648              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1649              :       LOGICAL, INTENT(IN), OPTIONAL                      :: soft_valid
    1650              :       CHARACTER(LEN=*), INTENT(IN), OPTIONAL             :: basis_type
    1651              : 
    1652              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_drho_elec'
    1653              : 
    1654              :       CHARACTER(LEN=default_string_length)               :: my_basis_type
    1655              :       INTEGER :: bcol, brow, dabqadb_func, handle, iatom, iatom_old, idir, igrid_level, ikind, &
    1656              :                  ikind_old, img, img_old, ipgf, iset, iset_old, itask, ithread, jatom, jatom_old, jkind, &
    1657              :                  jkind_old, jpgf, jset, jset_old, maxco, maxsgf_set, na1, na2, natoms, nb1, nb2, ncoa, &
    1658              :                  ncob, nimages, nseta, nsetb, ntasks, nthread, sgfa, sgfb
    1659            0 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
    1660            0 :                                                             npgfb, nsgfa, nsgfb
    1661            0 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
    1662              :       LOGICAL                                            :: atom_pair_changed, distributed_rs_grids, &
    1663              :                                                             do_kp, found, my_soft, use_subpatch
    1664              :       REAL(KIND=dp)                                      :: eps_rho_rspace, f, prefactor, radius, &
    1665              :                                                             scale, zetp
    1666              :       REAL(KIND=dp), DIMENSION(3)                        :: ra, rab, rab_inv, rb, rp
    1667            0 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: p_block, pab, sphi_a, sphi_b, work, &
    1668            0 :                                                             zeta, zetb
    1669            0 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: pabt, workt
    1670            0 :       TYPE(atom_pair_type), DIMENSION(:), POINTER        :: atom_pair_recv, atom_pair_send
    1671              :       TYPE(cell_type), POINTER                           :: cell
    1672            0 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: deltap
    1673              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1674              :       TYPE(gridlevel_info_type), POINTER                 :: gridlevel_info
    1675              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
    1676              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    1677            0 :          POINTER                                         :: sab_orb
    1678            0 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1679              :       TYPE(pw_env_type), POINTER                         :: pw_env
    1680            0 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1681              :       TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
    1682            0 :          POINTER                                         :: rs_descs
    1683            0 :       TYPE(realspace_grid_type), DIMENSION(:), POINTER   :: rs_rho
    1684              :       TYPE(task_list_type), POINTER                      :: task_list, task_list_soft
    1685            0 :       TYPE(task_type), DIMENSION(:), POINTER             :: tasks
    1686              : 
    1687            0 :       CALL timeset(routineN, handle)
    1688              : 
    1689            0 :       CPASSERT(PRESENT(matrix_p) .OR. PRESENT(matrix_p_kp))
    1690            0 :       do_kp = PRESENT(matrix_p_kp)
    1691              : 
    1692            0 :       NULLIFY (cell, dft_control, orb_basis_set, deltap, qs_kind_set, &
    1693            0 :                sab_orb, particle_set, rs_rho, pw_env, rs_descs, la_max, la_min, &
    1694            0 :                lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, p_block, sphi_a, &
    1695            0 :                sphi_b, zeta, zetb, first_sgfa, first_sgfb, tasks, pabt, workt)
    1696              : 
    1697              :       ! by default, the full density is calculated
    1698            0 :       my_soft = .FALSE.
    1699            0 :       IF (PRESENT(soft_valid)) my_soft = soft_valid
    1700              : 
    1701            0 :       IF (PRESENT(basis_type)) THEN
    1702            0 :          my_basis_type = basis_type
    1703              :       ELSE
    1704            0 :          my_basis_type = "ORB"
    1705              :       END IF
    1706              : 
    1707              :       CALL get_qs_env(qs_env=qs_env, &
    1708              :                       qs_kind_set=qs_kind_set, &
    1709              :                       cell=cell, &
    1710              :                       dft_control=dft_control, &
    1711              :                       particle_set=particle_set, &
    1712              :                       sab_orb=sab_orb, &
    1713            0 :                       pw_env=pw_env)
    1714              : 
    1715            0 :       SELECT CASE (my_basis_type)
    1716              :       CASE ("ORB")
    1717              :          CALL get_qs_env(qs_env=qs_env, &
    1718              :                          task_list=task_list, &
    1719            0 :                          task_list_soft=task_list_soft)
    1720              :       CASE ("AUX_FIT")
    1721              :          CALL get_qs_env(qs_env=qs_env, &
    1722            0 :                          task_list_soft=task_list_soft)
    1723            0 :          CALL get_admm_env(qs_env%admm_env, task_list_aux_fit=task_list)
    1724              :       END SELECT
    1725              : 
    1726              :       ! *** assign from pw_env
    1727            0 :       gridlevel_info => pw_env%gridlevel_info
    1728              : 
    1729              :       !   *** Allocate work storage ***
    1730            0 :       nthread = 1
    1731              :       CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
    1732              :                            maxco=maxco, &
    1733              :                            maxsgf_set=maxsgf_set, &
    1734            0 :                            basis_type=my_basis_type)
    1735            0 :       CALL reallocate(pabt, 1, maxco, 1, maxco, 0, nthread - 1)
    1736            0 :       CALL reallocate(workt, 1, maxco, 1, maxsgf_set, 0, nthread - 1)
    1737              : 
    1738              :       ! find maximum numbers
    1739            0 :       nimages = dft_control%nimages
    1740            0 :       CPASSERT(nimages == 1 .OR. do_kp)
    1741              : 
    1742            0 :       natoms = SIZE(particle_set)
    1743              : 
    1744              :       ! get the task lists
    1745            0 :       IF (my_soft) task_list => task_list_soft
    1746            0 :       CPASSERT(ASSOCIATED(task_list))
    1747            0 :       tasks => task_list%tasks
    1748            0 :       atom_pair_send => task_list%atom_pair_send
    1749            0 :       atom_pair_recv => task_list%atom_pair_recv
    1750            0 :       ntasks = task_list%ntasks
    1751              : 
    1752              :       ! *** set up the rs multi-grids
    1753            0 :       CPASSERT(ASSOCIATED(pw_env))
    1754            0 :       CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_rho)
    1755            0 :       DO igrid_level = 1, gridlevel_info%ngrid_levels
    1756            0 :          distributed_rs_grids = rs_rho(igrid_level)%desc%distributed
    1757              :       END DO
    1758              : 
    1759            0 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
    1760              : 
    1761              :       !   *** Initialize working density matrix ***
    1762              :       ! distributed rs grids require a matrix that will be changed
    1763              :       ! whereas this is not the case for replicated grids
    1764            0 :       ALLOCATE (deltap(nimages))
    1765            0 :       IF (distributed_rs_grids) THEN
    1766            0 :          DO img = 1, nimages
    1767              :          END DO
    1768              :          ! this matrix has no strict sparsity pattern in parallel
    1769              :          ! deltap%sparsity_id=-1
    1770            0 :          IF (do_kp) THEN
    1771            0 :             DO img = 1, nimages
    1772              :                CALL dbcsr_copy(deltap(img)%matrix, matrix_p_kp(img)%matrix, &
    1773            0 :                                name="DeltaP")
    1774              :             END DO
    1775              :          ELSE
    1776            0 :             CALL dbcsr_copy(deltap(1)%matrix, matrix_p, name="DeltaP")
    1777              :          END IF
    1778              :       ELSE
    1779            0 :          IF (do_kp) THEN
    1780            0 :             DO img = 1, nimages
    1781            0 :                deltap(img)%matrix => matrix_p_kp(img)%matrix
    1782              :             END DO
    1783              :          ELSE
    1784            0 :             deltap(1)%matrix => matrix_p
    1785              :          END IF
    1786              :       END IF
    1787              : 
    1788              :       ! distribute the matrix
    1789            0 :       IF (distributed_rs_grids) THEN
    1790              :          CALL rs_distribute_matrix(rs_descs=rs_descs, pmats=deltap, &
    1791              :                                    atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
    1792            0 :                                    nimages=nimages, scatter=.TRUE.)
    1793              :       END IF
    1794              : 
    1795              :       ! map all tasks on the grids
    1796              : 
    1797            0 :       ithread = 0
    1798            0 :       pab => pabt(:, :, ithread)
    1799            0 :       work => workt(:, :, ithread)
    1800              : 
    1801            0 :       loop_xyz: DO idir = 1, 3
    1802              : 
    1803            0 :          DO igrid_level = 1, gridlevel_info%ngrid_levels
    1804            0 :             CALL rs_grid_zero(rs_rho(igrid_level))
    1805              :          END DO
    1806              : 
    1807              :          iatom_old = -1; jatom_old = -1; iset_old = -1; jset_old = -1
    1808              :          ikind_old = -1; jkind_old = -1; img_old = -1
    1809            0 :          loop_tasks: DO itask = 1, ntasks
    1810              : 
    1811              :             !decode the atom pair and basis info
    1812            0 :             igrid_level = tasks(itask)%grid_level
    1813            0 :             img = tasks(itask)%image
    1814            0 :             iatom = tasks(itask)%iatom
    1815            0 :             jatom = tasks(itask)%jatom
    1816            0 :             iset = tasks(itask)%iset
    1817            0 :             jset = tasks(itask)%jset
    1818            0 :             ipgf = tasks(itask)%ipgf
    1819            0 :             jpgf = tasks(itask)%jpgf
    1820              : 
    1821            0 :             ikind = particle_set(iatom)%atomic_kind%kind_number
    1822            0 :             jkind = particle_set(jatom)%atomic_kind%kind_number
    1823              : 
    1824            0 :             IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old .OR. img .NE. img_old) THEN
    1825              : 
    1826            0 :                IF (iatom .NE. iatom_old) ra(:) = pbc(particle_set(iatom)%r, cell)
    1827              : 
    1828            0 :                IF (iatom <= jatom) THEN
    1829            0 :                   brow = iatom
    1830            0 :                   bcol = jatom
    1831              :                ELSE
    1832            0 :                   brow = jatom
    1833            0 :                   bcol = iatom
    1834              :                END IF
    1835              : 
    1836            0 :                IF (ikind .NE. ikind_old) THEN
    1837            0 :                   IF (my_soft) THEN
    1838              :                      CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
    1839            0 :                                       basis_type="ORB_SOFT")
    1840              :                   ELSE
    1841              :                      CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
    1842            0 :                                       basis_type=my_basis_type)
    1843              :                   END IF
    1844              :                   CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
    1845              :                                          first_sgf=first_sgfa, &
    1846              :                                          lmax=la_max, &
    1847              :                                          lmin=la_min, &
    1848              :                                          npgf=npgfa, &
    1849              :                                          nset=nseta, &
    1850              :                                          nsgf_set=nsgfa, &
    1851              :                                          sphi=sphi_a, &
    1852            0 :                                          zet=zeta)
    1853              :                END IF
    1854              : 
    1855            0 :                IF (jkind .NE. jkind_old) THEN
    1856            0 :                   IF (my_soft) THEN
    1857              :                      CALL get_qs_kind(qs_kind_set(jkind), basis_set=orb_basis_set, &
    1858            0 :                                       basis_type="ORB_SOFT")
    1859              :                   ELSE
    1860              :                      CALL get_qs_kind(qs_kind_set(jkind), basis_set=orb_basis_set, &
    1861            0 :                                       basis_type=my_basis_type)
    1862              :                   END IF
    1863              :                   CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
    1864              :                                          first_sgf=first_sgfb, &
    1865              :                                          lmax=lb_max, &
    1866              :                                          lmin=lb_min, &
    1867              :                                          npgf=npgfb, &
    1868              :                                          nset=nsetb, &
    1869              :                                          nsgf_set=nsgfb, &
    1870              :                                          sphi=sphi_b, &
    1871            0 :                                          zet=zetb)
    1872              :                END IF
    1873              : 
    1874              :                CALL dbcsr_get_block_p(matrix=deltap(img)%matrix, &
    1875            0 :                                       row=brow, col=bcol, BLOCK=p_block, found=found)
    1876            0 :                CPASSERT(found)
    1877              : 
    1878              :                iatom_old = iatom
    1879              :                jatom_old = jatom
    1880              :                ikind_old = ikind
    1881              :                jkind_old = jkind
    1882              :                img_old = img
    1883              :                atom_pair_changed = .TRUE.
    1884              : 
    1885              :             ELSE
    1886              : 
    1887              :                atom_pair_changed = .FALSE.
    1888              : 
    1889              :             END IF
    1890              : 
    1891            0 :             IF (atom_pair_changed .OR. iset_old .NE. iset .OR. jset_old .NE. jset) THEN
    1892              : 
    1893            0 :                ncoa = npgfa(iset)*ncoset(la_max(iset))
    1894            0 :                sgfa = first_sgfa(1, iset)
    1895            0 :                ncob = npgfb(jset)*ncoset(lb_max(jset))
    1896            0 :                sgfb = first_sgfb(1, jset)
    1897              : 
    1898            0 :                IF (iatom <= jatom) THEN
    1899              :                   CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
    1900              :                              1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
    1901              :                              p_block(sgfa, sgfb), SIZE(p_block, 1), &
    1902            0 :                              0.0_dp, work(1, 1), maxco)
    1903              :                   CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
    1904              :                              1.0_dp, work(1, 1), maxco, &
    1905              :                              sphi_b(1, sgfb), SIZE(sphi_b, 1), &
    1906            0 :                              0.0_dp, pab(1, 1), maxco)
    1907              :                ELSE
    1908              :                   CALL dgemm("N", "N", ncob, nsgfa(iset), nsgfb(jset), &
    1909              :                              1.0_dp, sphi_b(1, sgfb), SIZE(sphi_b, 1), &
    1910              :                              p_block(sgfb, sgfa), SIZE(p_block, 1), &
    1911            0 :                              0.0_dp, work(1, 1), maxco)
    1912              :                   CALL dgemm("N", "T", ncob, ncoa, nsgfa(iset), &
    1913              :                              1.0_dp, work(1, 1), maxco, &
    1914              :                              sphi_a(1, sgfa), SIZE(sphi_a, 1), &
    1915            0 :                              0.0_dp, pab(1, 1), maxco)
    1916              :                END IF
    1917              : 
    1918              :                iset_old = iset
    1919              :                jset_old = jset
    1920              : 
    1921              :             END IF
    1922              : 
    1923            0 :             rab(:) = tasks(itask)%rab
    1924            0 :             rb(:) = ra(:) + rab(:)
    1925            0 :             zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
    1926              : 
    1927            0 :             f = zetb(jpgf, jset)/zetp
    1928            0 :             rp(:) = ra(:) + f*rab(:)
    1929            0 :             prefactor = EXP(-zeta(ipgf, iset)*f*DOT_PRODUCT(rab, rab))
    1930              :             radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
    1931              :                                               lb_min=lb_min(jset), lb_max=lb_max(jset), &
    1932              :                                               ra=ra, rb=rb, rp=rp, &
    1933              :                                               zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
    1934            0 :                                               prefactor=prefactor, cutoff=1.0_dp)
    1935              : 
    1936            0 :             na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
    1937            0 :             na2 = ipgf*ncoset(la_max(iset))
    1938            0 :             nb1 = (jpgf - 1)*ncoset(lb_max(jset)) + 1
    1939            0 :             nb2 = jpgf*ncoset(lb_max(jset))
    1940              : 
    1941              :             ! takes the density matrix symmetry in account, i.e. off-diagonal blocks need to be mapped 'twice'
    1942            0 :             IF (iatom == jatom .AND. img == 1) THEN
    1943            0 :                scale = 1.0_dp
    1944              :             ELSE
    1945            0 :                scale = 2.0_dp
    1946              :             END IF
    1947              : 
    1948              :             ! check whether we need to use fawzi's generalised collocation scheme
    1949            0 :             IF (rs_rho(igrid_level)%desc%distributed) THEN
    1950              :                !tasks(4,:) is 0 for replicated, 1 for distributed 2 for exceptional distributed tasks
    1951            0 :                IF (tasks(itask)%dist_type .EQ. 2) THEN
    1952            0 :                   use_subpatch = .TRUE.
    1953              :                ELSE
    1954            0 :                   use_subpatch = .FALSE.
    1955              :                END IF
    1956              :             ELSE
    1957            0 :                use_subpatch = .FALSE.
    1958              :             END IF
    1959              : 
    1960            0 :             SELECT CASE (idir)
    1961              :             CASE (1)
    1962            0 :                dabqadb_func = GRID_FUNC_DABpADB_X
    1963              :             CASE (2)
    1964            0 :                dabqadb_func = GRID_FUNC_DABpADB_Y
    1965              :             CASE (3)
    1966            0 :                dabqadb_func = GRID_FUNC_DABpADB_Z
    1967              :             CASE DEFAULT
    1968            0 :                CPABORT("invalid idir")
    1969              :             END SELECT
    1970              : 
    1971            0 :             IF (iatom <= jatom) THEN
    1972              :                CALL collocate_pgf_product( &
    1973              :                   la_max(iset), zeta(ipgf, iset), la_min(iset), &
    1974              :                   lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
    1975              :                   ra, rab, scale, pab, na1 - 1, nb1 - 1, &
    1976              :                   rs_rho(igrid_level), &
    1977              :                   radius=radius, ga_gb_function=dabqadb_func, &
    1978            0 :                   use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern)
    1979              :             ELSE
    1980            0 :                rab_inv = -rab
    1981              :                CALL collocate_pgf_product( &
    1982              :                   lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
    1983              :                   la_max(iset), zeta(ipgf, iset), la_min(iset), &
    1984              :                   rb, rab_inv, scale, pab, nb1 - 1, na1 - 1, &
    1985              :                   rs_rho(igrid_level), &
    1986              :                   radius=radius, ga_gb_function=dabqadb_func, &
    1987            0 :                   use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern)
    1988              :             END IF
    1989              : 
    1990              :          END DO loop_tasks
    1991              : 
    1992            0 :          CALL density_rs2pw(pw_env, rs_rho, drho(idir), drho_gspace(idir))
    1993              : 
    1994              :       END DO loop_xyz
    1995              : 
    1996              :       !   *** Release work storage ***
    1997            0 :       IF (distributed_rs_grids) THEN
    1998            0 :          CALL dbcsr_deallocate_matrix_set(deltap)
    1999              :       ELSE
    2000            0 :          DO img = 1, nimages
    2001            0 :             NULLIFY (deltap(img)%matrix)
    2002              :          END DO
    2003            0 :          DEALLOCATE (deltap)
    2004              :       END IF
    2005              : 
    2006            0 :       DEALLOCATE (pabt, workt)
    2007              : 
    2008            0 :       CALL timestop(handle)
    2009              : 
    2010            0 :    END SUBROUTINE calculate_drho_elec
    2011              : 
    2012              : ! **************************************************************************************************
    2013              : !> \brief Computes the gradient wrt. nuclear coordinates of a density on the grid
    2014              : !>        The density is given in terms of the density matrix_p
    2015              : !> \param matrix_p Density matrix
    2016              : !> \param matrix_p_kp ...
    2017              : !> \param drho Density gradient on the grid
    2018              : !> \param drho_gspace Density gradient on the reciprocal grid
    2019              : !> \param qs_env ...
    2020              : !> \param soft_valid ...
    2021              : !> \param basis_type ...
    2022              : !> \param beta Derivative direction
    2023              : !> \param lambda Atom index
    2024              : !> \note SL, ED 2021
    2025              : !>       Adapted from calculate_drho_elec
    2026              : ! **************************************************************************************************
    2027          252 :    SUBROUTINE calculate_drho_elec_dR(matrix_p, matrix_p_kp, drho, drho_gspace, qs_env, &
    2028              :                                      soft_valid, basis_type, beta, lambda)
    2029              : 
    2030              :       TYPE(dbcsr_type), OPTIONAL, TARGET                 :: matrix_p
    2031              :       TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
    2032              :          POINTER                                         :: matrix_p_kp
    2033              :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                       :: drho
    2034              :       TYPE(pw_c1d_gs_type), INTENT(INOUT) :: drho_gspace
    2035              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2036              :       LOGICAL, INTENT(IN), OPTIONAL                      :: soft_valid
    2037              :       CHARACTER(LEN=*), INTENT(IN), OPTIONAL             :: basis_type
    2038              :       INTEGER, INTENT(IN)                                :: beta, lambda
    2039              : 
    2040              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_drho_elec_dR'
    2041              : 
    2042              :       CHARACTER(LEN=default_string_length)               :: my_basis_type
    2043              :       INTEGER :: bcol, brow, dabqadb_func, handle, iatom, iatom_old, igrid_level, ikind, &
    2044              :                  ikind_old, img, img_old, ipgf, iset, iset_old, itask, ithread, jatom, jatom_old, jkind, &
    2045              :                  jkind_old, jpgf, jset, jset_old, maxco, maxsgf_set, na1, na2, natoms, nb1, nb2, ncoa, &
    2046              :                  ncob, nimages, nseta, nsetb, ntasks, nthread, sgfa, sgfb
    2047          252 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
    2048          252 :                                                             npgfb, nsgfa, nsgfb
    2049          252 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
    2050              :       LOGICAL                                            :: atom_pair_changed, distributed_rs_grids, &
    2051              :                                                             do_kp, found, my_soft, use_subpatch
    2052              :       REAL(KIND=dp)                                      :: eps_rho_rspace, f, prefactor, radius, &
    2053              :                                                             scale, zetp
    2054              :       REAL(KIND=dp), DIMENSION(3)                        :: ra, rab, rab_inv, rb, rp
    2055          252 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: p_block, pab, sphi_a, sphi_b, work, &
    2056          252 :                                                             zeta, zetb
    2057          252 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: pabt, workt
    2058          252 :       TYPE(atom_pair_type), DIMENSION(:), POINTER        :: atom_pair_recv, atom_pair_send
    2059              :       TYPE(cell_type), POINTER                           :: cell
    2060          252 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: deltap
    2061              :       TYPE(dft_control_type), POINTER                    :: dft_control
    2062              :       TYPE(gridlevel_info_type), POINTER                 :: gridlevel_info
    2063              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
    2064          252 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    2065              :       TYPE(pw_env_type), POINTER                         :: pw_env
    2066          252 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    2067              :       TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
    2068          252 :          POINTER                                         :: rs_descs
    2069          252 :       TYPE(realspace_grid_type), DIMENSION(:), POINTER   :: rs_rho
    2070              :       TYPE(task_list_type), POINTER                      :: task_list, task_list_soft
    2071          252 :       TYPE(task_type), DIMENSION(:), POINTER             :: tasks
    2072              : 
    2073          252 :       CALL timeset(routineN, handle)
    2074              : 
    2075          252 :       CPASSERT(PRESENT(matrix_p) .OR. PRESENT(matrix_p_kp))
    2076          252 :       do_kp = PRESENT(matrix_p_kp)
    2077              : 
    2078          252 :       NULLIFY (cell, dft_control, orb_basis_set, deltap, qs_kind_set, &
    2079          252 :                particle_set, rs_rho, pw_env, rs_descs, la_max, la_min, lb_max, &
    2080          252 :                lb_min, npgfa, npgfb, nsgfa, nsgfb, p_block, sphi_a, sphi_b, &
    2081          252 :                zeta, zetb, first_sgfa, first_sgfb, tasks, pabt, workt)
    2082              : 
    2083              :       ! by default, the full density is calculated
    2084          252 :       my_soft = .FALSE.
    2085          252 :       IF (PRESENT(soft_valid)) my_soft = soft_valid
    2086              : 
    2087          252 :       IF (PRESENT(basis_type)) THEN
    2088            0 :          my_basis_type = basis_type
    2089              :       ELSE
    2090          252 :          my_basis_type = "ORB"
    2091              :       END IF
    2092              : 
    2093              :       CALL get_qs_env(qs_env=qs_env, &
    2094              :                       qs_kind_set=qs_kind_set, &
    2095              :                       cell=cell, &
    2096              :                       dft_control=dft_control, &
    2097              :                       particle_set=particle_set, &
    2098          252 :                       pw_env=pw_env)
    2099              : 
    2100          252 :       SELECT CASE (my_basis_type)
    2101              :       CASE ("ORB")
    2102              :          CALL get_qs_env(qs_env=qs_env, &
    2103              :                          task_list=task_list, &
    2104          252 :                          task_list_soft=task_list_soft)
    2105              :       CASE ("AUX_FIT")
    2106              :          CALL get_qs_env(qs_env=qs_env, &
    2107            0 :                          task_list_soft=task_list_soft)
    2108          252 :          CALL get_admm_env(qs_env%admm_env, task_list_aux_fit=task_list)
    2109              :       END SELECT
    2110              : 
    2111              :       ! *** assign from pw_env
    2112          252 :       gridlevel_info => pw_env%gridlevel_info
    2113              : 
    2114              :       !   *** Allocate work storage ***
    2115          252 :       nthread = 1
    2116              :       CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
    2117              :                            maxco=maxco, &
    2118              :                            maxsgf_set=maxsgf_set, &
    2119          252 :                            basis_type=my_basis_type)
    2120          252 :       CALL reallocate(pabt, 1, maxco, 1, maxco, 0, nthread - 1)
    2121          252 :       CALL reallocate(workt, 1, maxco, 1, maxsgf_set, 0, nthread - 1)
    2122              : 
    2123              :       ! find maximum numbers
    2124          252 :       nimages = dft_control%nimages
    2125          252 :       CPASSERT(nimages == 1 .OR. do_kp)
    2126              : 
    2127          252 :       natoms = SIZE(particle_set)
    2128              : 
    2129              :       ! get the task lists
    2130          252 :       IF (my_soft) task_list => task_list_soft
    2131          252 :       CPASSERT(ASSOCIATED(task_list))
    2132          252 :       tasks => task_list%tasks
    2133          252 :       atom_pair_send => task_list%atom_pair_send
    2134          252 :       atom_pair_recv => task_list%atom_pair_recv
    2135          252 :       ntasks = task_list%ntasks
    2136              : 
    2137              :       ! *** set up the rs multi-grids
    2138          252 :       CPASSERT(ASSOCIATED(pw_env))
    2139          252 :       CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_rho)
    2140          774 :       DO igrid_level = 1, gridlevel_info%ngrid_levels
    2141          774 :          distributed_rs_grids = rs_rho(igrid_level)%desc%distributed
    2142              :       END DO
    2143              : 
    2144          252 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
    2145              : 
    2146              :       !   *** Initialize working density matrix ***
    2147              :       ! distributed rs grids require a matrix that will be changed
    2148              :       ! whereas this is not the case for replicated grids
    2149         1008 :       ALLOCATE (deltap(nimages))
    2150          252 :       IF (distributed_rs_grids) THEN
    2151            0 :          DO img = 1, nimages
    2152              :          END DO
    2153              :          ! this matrix has no strict sparsity pattern in parallel
    2154              :          ! deltap%sparsity_id=-1
    2155            0 :          IF (do_kp) THEN
    2156            0 :             DO img = 1, nimages
    2157              :                CALL dbcsr_copy(deltap(img)%matrix, matrix_p_kp(img)%matrix, &
    2158            0 :                                name="DeltaP")
    2159              :             END DO
    2160              :          ELSE
    2161            0 :             CALL dbcsr_copy(deltap(1)%matrix, matrix_p, name="DeltaP")
    2162              :          END IF
    2163              :       ELSE
    2164          252 :          IF (do_kp) THEN
    2165            0 :             DO img = 1, nimages
    2166            0 :                deltap(img)%matrix => matrix_p_kp(img)%matrix
    2167              :             END DO
    2168              :          ELSE
    2169          252 :             deltap(1)%matrix => matrix_p
    2170              :          END IF
    2171              :       END IF
    2172              : 
    2173              :       ! distribute the matrix
    2174          252 :       IF (distributed_rs_grids) THEN
    2175              :          CALL rs_distribute_matrix(rs_descs=rs_descs, pmats=deltap, &
    2176              :                                    atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
    2177            0 :                                    nimages=nimages, scatter=.TRUE.)
    2178              :       END IF
    2179              : 
    2180              :       ! map all tasks on the grids
    2181              : 
    2182          252 :       ithread = 0
    2183          252 :       pab => pabt(:, :, ithread)
    2184          252 :       work => workt(:, :, ithread)
    2185              : 
    2186          774 :       DO igrid_level = 1, gridlevel_info%ngrid_levels
    2187          774 :          CALL rs_grid_zero(rs_rho(igrid_level))
    2188              :       END DO
    2189              : 
    2190              :       iatom_old = -1; jatom_old = -1; iset_old = -1; jset_old = -1
    2191              :       ikind_old = -1; jkind_old = -1; img_old = -1
    2192        16506 :       loop_tasks: DO itask = 1, ntasks
    2193              : 
    2194              :          !decode the atom pair and basis info
    2195        16254 :          igrid_level = tasks(itask)%grid_level
    2196        16254 :          img = tasks(itask)%image
    2197        16254 :          iatom = tasks(itask)%iatom
    2198        16254 :          jatom = tasks(itask)%jatom
    2199        16254 :          iset = tasks(itask)%iset
    2200        16254 :          jset = tasks(itask)%jset
    2201        16254 :          ipgf = tasks(itask)%ipgf
    2202        16254 :          jpgf = tasks(itask)%jpgf
    2203              : 
    2204        16254 :          ikind = particle_set(iatom)%atomic_kind%kind_number
    2205        16254 :          jkind = particle_set(jatom)%atomic_kind%kind_number
    2206              : 
    2207        16254 :          IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old .OR. img .NE. img_old) THEN
    2208              : 
    2209         1296 :             IF (iatom .NE. iatom_old) ra(:) = pbc(particle_set(iatom)%r, cell)
    2210              : 
    2211         1296 :             IF (iatom <= jatom) THEN
    2212          864 :                brow = iatom
    2213          864 :                bcol = jatom
    2214              :             ELSE
    2215          432 :                brow = jatom
    2216          432 :                bcol = iatom
    2217              :             END IF
    2218              : 
    2219         1296 :             IF (ikind .NE. ikind_old) THEN
    2220          252 :                IF (my_soft) THEN
    2221              :                   CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
    2222            0 :                                    basis_type="ORB_SOFT")
    2223              :                ELSE
    2224              :                   CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
    2225          252 :                                    basis_type=my_basis_type)
    2226              :                END IF
    2227              :                CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
    2228              :                                       first_sgf=first_sgfa, &
    2229              :                                       lmax=la_max, &
    2230              :                                       lmin=la_min, &
    2231              :                                       npgf=npgfa, &
    2232              :                                       nset=nseta, &
    2233              :                                       nsgf_set=nsgfa, &
    2234              :                                       sphi=sphi_a, &
    2235          252 :                                       zet=zeta)
    2236              :             END IF
    2237              : 
    2238         1296 :             IF (jkind .NE. jkind_old) THEN
    2239          864 :                IF (my_soft) THEN
    2240              :                   CALL get_qs_kind(qs_kind_set(jkind), basis_set=orb_basis_set, &
    2241            0 :                                    basis_type="ORB_SOFT")
    2242              :                ELSE
    2243              :                   CALL get_qs_kind(qs_kind_set(jkind), basis_set=orb_basis_set, &
    2244          864 :                                    basis_type=my_basis_type)
    2245              :                END IF
    2246              :                CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
    2247              :                                       first_sgf=first_sgfb, &
    2248              :                                       lmax=lb_max, &
    2249              :                                       lmin=lb_min, &
    2250              :                                       npgf=npgfb, &
    2251              :                                       nset=nsetb, &
    2252              :                                       nsgf_set=nsgfb, &
    2253              :                                       sphi=sphi_b, &
    2254          864 :                                       zet=zetb)
    2255              :             END IF
    2256              : 
    2257              :             CALL dbcsr_get_block_p(matrix=deltap(img)%matrix, &
    2258         1296 :                                    row=brow, col=bcol, BLOCK=p_block, found=found)
    2259         1296 :             CPASSERT(found)
    2260              : 
    2261              :             iatom_old = iatom
    2262              :             jatom_old = jatom
    2263              :             ikind_old = ikind
    2264              :             jkind_old = jkind
    2265              :             img_old = img
    2266              :             atom_pair_changed = .TRUE.
    2267              : 
    2268              :          ELSE
    2269              : 
    2270              :             atom_pair_changed = .FALSE.
    2271              : 
    2272              :          END IF
    2273              : 
    2274        16254 :          IF (atom_pair_changed .OR. iset_old .NE. iset .OR. jset_old .NE. jset) THEN
    2275              : 
    2276         1296 :             ncoa = npgfa(iset)*ncoset(la_max(iset))
    2277         1296 :             sgfa = first_sgfa(1, iset)
    2278         1296 :             ncob = npgfb(jset)*ncoset(lb_max(jset))
    2279         1296 :             sgfb = first_sgfb(1, jset)
    2280              : 
    2281         1296 :             IF (iatom <= jatom) THEN
    2282              :                CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
    2283              :                           1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
    2284              :                           p_block(sgfa, sgfb), SIZE(p_block, 1), &
    2285          864 :                           0.0_dp, work(1, 1), maxco)
    2286              :                CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
    2287              :                           1.0_dp, work(1, 1), maxco, &
    2288              :                           sphi_b(1, sgfb), SIZE(sphi_b, 1), &
    2289          864 :                           0.0_dp, pab(1, 1), maxco)
    2290              :             ELSE
    2291              :                CALL dgemm("N", "N", ncob, nsgfa(iset), nsgfb(jset), &
    2292              :                           1.0_dp, sphi_b(1, sgfb), SIZE(sphi_b, 1), &
    2293              :                           p_block(sgfb, sgfa), SIZE(p_block, 1), &
    2294          432 :                           0.0_dp, work(1, 1), maxco)
    2295              :                CALL dgemm("N", "T", ncob, ncoa, nsgfa(iset), &
    2296              :                           1.0_dp, work(1, 1), maxco, &
    2297              :                           sphi_a(1, sgfa), SIZE(sphi_a, 1), &
    2298          432 :                           0.0_dp, pab(1, 1), maxco)
    2299              :             END IF
    2300              : 
    2301              :             iset_old = iset
    2302              :             jset_old = jset
    2303              : 
    2304              :          END IF
    2305              : 
    2306        65016 :          rab(:) = tasks(itask)%rab
    2307        65016 :          rb(:) = ra(:) + rab(:)
    2308        16254 :          zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
    2309              : 
    2310        16254 :          f = zetb(jpgf, jset)/zetp
    2311        65016 :          rp(:) = ra(:) + f*rab(:)
    2312        65016 :          prefactor = EXP(-zeta(ipgf, iset)*f*DOT_PRODUCT(rab, rab))
    2313              :          radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
    2314              :                                            lb_min=lb_min(jset), lb_max=lb_max(jset), &
    2315              :                                            ra=ra, rb=rb, rp=rp, &
    2316              :                                            zetp=zetp, eps=eps_rho_rspace, &
    2317        16254 :                                            prefactor=prefactor, cutoff=1.0_dp)
    2318              : 
    2319        16254 :          na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
    2320        16254 :          na2 = ipgf*ncoset(la_max(iset))
    2321        16254 :          nb1 = (jpgf - 1)*ncoset(lb_max(jset)) + 1
    2322        16254 :          nb2 = jpgf*ncoset(lb_max(jset))
    2323              : 
    2324              :          ! takes the density matrix symmetry in account, i.e. off-diagonal blocks need to be mapped 'twice'
    2325        16254 :          IF (iatom == jatom .AND. img == 1) THEN
    2326         8100 :             scale = 1.0_dp
    2327              :          ELSE
    2328         8154 :             scale = 2.0_dp
    2329              :          END IF
    2330              : 
    2331              :          ! check whether we need to use fawzi's generalised collocation scheme
    2332        16254 :          IF (rs_rho(igrid_level)%desc%distributed) THEN
    2333              :             !tasks(4,:) is 0 for replicated, 1 for distributed 2 for exceptional distributed tasks
    2334            0 :             IF (tasks(itask)%dist_type .EQ. 2) THEN
    2335            0 :                use_subpatch = .TRUE.
    2336              :             ELSE
    2337            0 :                use_subpatch = .FALSE.
    2338              :             END IF
    2339              :          ELSE
    2340        16254 :             use_subpatch = .FALSE.
    2341              :          END IF
    2342              : 
    2343        21672 :          SELECT CASE (beta)
    2344              :          CASE (1)
    2345         5418 :             dabqadb_func = GRID_FUNC_DAB_X
    2346              :          CASE (2)
    2347         5418 :             dabqadb_func = GRID_FUNC_DAB_Y
    2348              :          CASE (3)
    2349         5418 :             dabqadb_func = GRID_FUNC_DAB_Z
    2350              :          CASE DEFAULT
    2351        16254 :             CPABORT("invalid beta")
    2352              :          END SELECT
    2353              : 
    2354        16506 :          IF (iatom <= jatom) THEN
    2355        10854 :             IF (iatom == lambda) &
    2356              :                CALL collocate_pgf_product( &
    2357              :                la_max(iset), zeta(ipgf, iset), la_min(iset), &
    2358              :                lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
    2359              :                ra, rab, scale, pab, na1 - 1, nb1 - 1, &
    2360              :                rsgrid=rs_rho(igrid_level), &
    2361              :                ga_gb_function=dabqadb_func, radius=radius, &
    2362              :                use_subpatch=use_subpatch, &
    2363         3618 :                subpatch_pattern=tasks(itask)%subpatch_pattern)
    2364        10854 :             IF (jatom == lambda) &
    2365              :                CALL collocate_pgf_product( &
    2366              :                la_max(iset), zeta(ipgf, iset), la_min(iset), &
    2367              :                lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
    2368              :                ra, rab, scale, pab, na1 - 1, nb1 - 1, &
    2369              :                rsgrid=rs_rho(igrid_level), &
    2370              :                ga_gb_function=dabqadb_func + 3, radius=radius, &
    2371              :                use_subpatch=use_subpatch, &
    2372         3618 :                subpatch_pattern=tasks(itask)%subpatch_pattern)
    2373              :          ELSE
    2374        21600 :             rab_inv = -rab
    2375         5400 :             IF (jatom == lambda) &
    2376              :                CALL collocate_pgf_product( &
    2377              :                lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
    2378              :                la_max(iset), zeta(ipgf, iset), la_min(iset), &
    2379              :                rb, rab_inv, scale, pab, nb1 - 1, na1 - 1, &
    2380              :                rs_rho(igrid_level), &
    2381              :                ga_gb_function=dabqadb_func, radius=radius, &
    2382              :                use_subpatch=use_subpatch, &
    2383         1800 :                subpatch_pattern=tasks(itask)%subpatch_pattern)
    2384         5400 :             IF (iatom == lambda) &
    2385              :                CALL collocate_pgf_product( &
    2386              :                lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
    2387              :                la_max(iset), zeta(ipgf, iset), la_min(iset), &
    2388              :                rb, rab_inv, scale, pab, nb1 - 1, na1 - 1, &
    2389              :                rs_rho(igrid_level), &
    2390              :                ga_gb_function=dabqadb_func + 3, radius=radius, &
    2391              :                use_subpatch=use_subpatch, &
    2392         1800 :                subpatch_pattern=tasks(itask)%subpatch_pattern)
    2393              :          END IF
    2394              : 
    2395              :       END DO loop_tasks
    2396              : 
    2397          252 :       CALL density_rs2pw(pw_env, rs_rho, drho, drho_gspace)
    2398              : 
    2399              :       !   *** Release work storage ***
    2400          252 :       IF (distributed_rs_grids) THEN
    2401            0 :          CALL dbcsr_deallocate_matrix_set(deltap)
    2402              :       ELSE
    2403          504 :          DO img = 1, nimages
    2404          504 :             NULLIFY (deltap(img)%matrix)
    2405              :          END DO
    2406          252 :          DEALLOCATE (deltap)
    2407              :       END IF
    2408              : 
    2409          252 :       DEALLOCATE (pabt, workt)
    2410              : 
    2411          252 :       CALL timestop(handle)
    2412              : 
    2413          504 :    END SUBROUTINE calculate_drho_elec_dR
    2414              : 
    2415              : ! **************************************************************************************************
    2416              : !> \brief maps a single gaussian on the grid
    2417              : !> \param rho ...
    2418              : !> \param rho_gspace ...
    2419              : !> \param atomic_kind_set ...
    2420              : !> \param qs_kind_set ...
    2421              : !> \param cell ...
    2422              : !> \param dft_control ...
    2423              : !> \param particle_set ...
    2424              : !> \param pw_env ...
    2425              : !> \param required_function ...
    2426              : !> \param basis_type ...
    2427              : !> \par History
    2428              : !>      08.2022 created from calculate_wavefunction
    2429              : !> \note
    2430              : !>      modified calculate_wave function assuming that the collocation of only a single Gaussian is required.
    2431              : !>      chooses a basis function (in contrast to calculate_rho_core or calculate_rho_single_gaussian)
    2432              : ! **************************************************************************************************
    2433        28573 :    SUBROUTINE collocate_single_gaussian(rho, rho_gspace, &
    2434              :                                         atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, &
    2435              :                                         pw_env, required_function, basis_type)
    2436              : 
    2437              :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                       :: rho
    2438              :       TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_gspace
    2439              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    2440              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    2441              :       TYPE(cell_type), POINTER                           :: cell
    2442              :       TYPE(dft_control_type), POINTER                    :: dft_control
    2443              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    2444              :       TYPE(pw_env_type), POINTER                         :: pw_env
    2445              :       INTEGER, INTENT(IN)                                :: required_function
    2446              :       CHARACTER(LEN=*), INTENT(IN), OPTIONAL             :: basis_type
    2447              : 
    2448              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'collocate_single_gaussian'
    2449              : 
    2450              :       CHARACTER(LEN=default_string_length)               :: my_basis_type
    2451              :       INTEGER :: group_size, handle, i, iatom, igrid_level, ikind, ipgf, iset, maxco, maxsgf_set, &
    2452              :                  my_index, my_pos, na1, na2, natom, ncoa, nseta, offset, sgfa
    2453        28573 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: where_is_the_point
    2454        28573 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, npgfa, nsgfa
    2455        28573 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa
    2456              :       LOGICAL                                            :: found
    2457              :       REAL(KIND=dp)                                      :: dab, eps_rho_rspace, radius, scale
    2458              :       REAL(KIND=dp), DIMENSION(3)                        :: ra
    2459        28573 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: pab, sphi_a, zeta
    2460              :       TYPE(gridlevel_info_type), POINTER                 :: gridlevel_info
    2461              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
    2462              :       TYPE(mp_comm_type)                                 :: group
    2463        28573 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
    2464        28573 :       TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:) :: mgrid_gspace
    2465        28573 :       TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:)           ::  mgrid_rspace
    2466        28573 :       TYPE(realspace_grid_type), DIMENSION(:), POINTER   :: rs_rho
    2467              : 
    2468        28573 :       IF (PRESENT(basis_type)) THEN
    2469        28573 :          my_basis_type = basis_type
    2470              :       ELSE
    2471            0 :          my_basis_type = "ORB"
    2472              :       END IF
    2473              : 
    2474        28573 :       CALL timeset(routineN, handle)
    2475              : 
    2476        28573 :       NULLIFY (orb_basis_set, pab, la_max, la_min, npgfa, nsgfa, sphi_a, &
    2477        28573 :                zeta, first_sgfa, rs_rho, pw_pools)
    2478              : 
    2479              :       ! *** set up the pw multi-grids
    2480        28573 :       CPASSERT(ASSOCIATED(pw_env))
    2481              :       CALL pw_env_get(pw_env, rs_grids=rs_rho, pw_pools=pw_pools, &
    2482        28573 :                       gridlevel_info=gridlevel_info)
    2483              : 
    2484        28573 :       CALL pw_pools_create_pws(pw_pools, mgrid_gspace)
    2485        28573 :       CALL pw_pools_create_pws(pw_pools, mgrid_rspace)
    2486              : 
    2487              :       ! *** set up rs multi-grids
    2488       142865 :       DO igrid_level = 1, gridlevel_info%ngrid_levels
    2489       142865 :          CALL rs_grid_zero(rs_rho(igrid_level))
    2490              :       END DO
    2491              : 
    2492        28573 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
    2493              : !   *** Allocate work storage ***
    2494        28573 :       CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
    2495              :       CALL get_qs_kind_set(qs_kind_set, &
    2496              :                            maxco=maxco, &
    2497              :                            maxsgf_set=maxsgf_set, &
    2498        28573 :                            basis_type=my_basis_type)
    2499              : 
    2500        85719 :       ALLOCATE (pab(maxco, 1))
    2501              : 
    2502        28573 :       offset = 0
    2503        28573 :       group = mgrid_rspace(1)%pw_grid%para%group
    2504        28573 :       my_pos = mgrid_rspace(1)%pw_grid%para%group%mepos
    2505        28573 :       group_size = mgrid_rspace(1)%pw_grid%para%group%num_pe
    2506        85719 :       ALLOCATE (where_is_the_point(0:group_size - 1))
    2507              : 
    2508       117551 :       DO iatom = 1, natom
    2509        88978 :          ikind = particle_set(iatom)%atomic_kind%kind_number
    2510        88978 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=my_basis_type)
    2511              :          CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
    2512              :                                 first_sgf=first_sgfa, &
    2513              :                                 lmax=la_max, &
    2514              :                                 lmin=la_min, &
    2515              :                                 npgf=npgfa, &
    2516              :                                 nset=nseta, &
    2517              :                                 nsgf_set=nsgfa, &
    2518              :                                 sphi=sphi_a, &
    2519        88978 :                                 zet=zeta)
    2520        88978 :          ra(:) = pbc(particle_set(iatom)%r, cell)
    2521        88978 :          dab = 0.0_dp
    2522              : 
    2523      1047571 :          DO iset = 1, nseta
    2524              : 
    2525       841042 :             ncoa = npgfa(iset)*ncoset(la_max(iset))
    2526       841042 :             sgfa = first_sgfa(1, iset)
    2527              : 
    2528       841042 :             found = .FALSE.
    2529       841042 :             my_index = 0
    2530      3176013 :             DO i = 1, nsgfa(iset)
    2531      3176013 :                IF (offset + i == required_function) THEN
    2532              :                   my_index = i
    2533              :                   found = .TRUE.
    2534              :                   EXIT
    2535              :                END IF
    2536              :             END DO
    2537              : 
    2538       841042 :             IF (found) THEN
    2539              : 
    2540       523769 :                pab(1:ncoa, 1) = sphi_a(1:ncoa, sgfa + my_index - 1)
    2541              : 
    2542        58202 :                DO ipgf = 1, npgfa(iset)
    2543              : 
    2544        29629 :                   na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
    2545        29629 :                   na2 = ipgf*ncoset(la_max(iset))
    2546              : 
    2547        29629 :                   scale = 1.0_dp
    2548        29629 :                   igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
    2549              : 
    2550        58202 :                   IF (map_gaussian_here(rs_rho(igrid_level), cell%h_inv, ra, offset, group_size, my_pos)) THEN
    2551              :                      radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
    2552              :                                                        lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
    2553              :                                                        zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
    2554        27872 :                                                        prefactor=1.0_dp, cutoff=1.0_dp)
    2555              : 
    2556              :                      CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), la_min(iset), &
    2557              :                                                 0, 0.0_dp, 0, &
    2558              :                                                 ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
    2559              :                                                 scale, pab, na1 - 1, 0, rs_rho(igrid_level), &
    2560        27872 :                                                 radius=radius, ga_gb_function=GRID_FUNC_AB)
    2561              :                   END IF
    2562              : 
    2563              :                END DO
    2564              : 
    2565              :             END IF
    2566              : 
    2567       930020 :             offset = offset + nsgfa(iset)
    2568              : 
    2569              :          END DO
    2570              : 
    2571              :       END DO
    2572              : 
    2573       142865 :       DO igrid_level = 1, gridlevel_info%ngrid_levels
    2574              :          CALL transfer_rs2pw(rs_rho(igrid_level), &
    2575       142865 :                              mgrid_rspace(igrid_level))
    2576              :       END DO
    2577              : 
    2578        28573 :       CALL pw_zero(rho_gspace)
    2579       142865 :       DO igrid_level = 1, gridlevel_info%ngrid_levels
    2580              :          CALL pw_transfer(mgrid_rspace(igrid_level), &
    2581       114292 :                           mgrid_gspace(igrid_level))
    2582       142865 :          CALL pw_axpy(mgrid_gspace(igrid_level), rho_gspace)
    2583              :       END DO
    2584              : 
    2585        28573 :       CALL pw_transfer(rho_gspace, rho)
    2586              : 
    2587              :       ! Release work storage
    2588        28573 :       DEALLOCATE (pab)
    2589              : 
    2590              :       ! give back the pw multi-grids
    2591        28573 :       CALL pw_pools_give_back_pws(pw_pools, mgrid_gspace)
    2592        28573 :       CALL pw_pools_give_back_pws(pw_pools, mgrid_rspace)
    2593              : 
    2594        28573 :       CALL timestop(handle)
    2595              : 
    2596       171438 :    END SUBROUTINE collocate_single_gaussian
    2597              : 
    2598              : ! **************************************************************************************************
    2599              : !> \brief maps a given wavefunction on the grid
    2600              : !> \param mo_vectors ...
    2601              : !> \param ivector ...
    2602              : !> \param rho ...
    2603              : !> \param rho_gspace ...
    2604              : !> \param atomic_kind_set ...
    2605              : !> \param qs_kind_set ...
    2606              : !> \param cell ...
    2607              : !> \param dft_control ...
    2608              : !> \param particle_set ...
    2609              : !> \param pw_env ...
    2610              : !> \param basis_type ...
    2611              : !> \par History
    2612              : !>      08.2002 created [Joost VandeVondele]
    2613              : !>      03.2006 made independent of qs_env [Joost VandeVondele]
    2614              : !>      08.2024 call collocate_function [JGH]
    2615              : ! **************************************************************************************************
    2616         1354 :    SUBROUTINE calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, &
    2617              :                                      atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, &
    2618              :                                      pw_env, basis_type)
    2619              :       TYPE(cp_fm_type), INTENT(IN)                       :: mo_vectors
    2620              :       INTEGER, INTENT(IN)                                :: ivector
    2621              :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                :: rho
    2622              :       TYPE(pw_c1d_gs_type), INTENT(INOUT)                :: rho_gspace
    2623              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    2624              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    2625              :       TYPE(cell_type), POINTER                           :: cell
    2626              :       TYPE(dft_control_type), POINTER                    :: dft_control
    2627              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    2628              :       TYPE(pw_env_type), POINTER                         :: pw_env
    2629              :       CHARACTER(LEN=*), INTENT(IN), OPTIONAL             :: basis_type
    2630              : 
    2631              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_wavefunction'
    2632              : 
    2633              :       INTEGER                                            :: handle, i, nao
    2634              :       LOGICAL                                            :: local
    2635              :       REAL(KIND=dp)                                      :: eps_rho_rspace
    2636              :       REAL(KIND=dp), DIMENSION(:), POINTER               :: eigenvector
    2637              : 
    2638         1354 :       CALL timeset(routineN, handle)
    2639              : 
    2640         1354 :       CALL cp_fm_get_info(matrix=mo_vectors, nrow_global=nao)
    2641         4062 :       ALLOCATE (eigenvector(nao))
    2642        25180 :       DO i = 1, nao
    2643        25180 :          CALL cp_fm_get_element(mo_vectors, i, ivector, eigenvector(i), local)
    2644              :       END DO
    2645              : 
    2646         1354 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
    2647              : 
    2648              :       CALL collocate_function(eigenvector, rho, rho_gspace, &
    2649              :                               atomic_kind_set, qs_kind_set, cell, particle_set, pw_env, &
    2650         2708 :                               eps_rho_rspace, basis_type)
    2651              : 
    2652         1354 :       DEALLOCATE (eigenvector)
    2653              : 
    2654         1354 :       CALL timestop(handle)
    2655              : 
    2656         1354 :    END SUBROUTINE calculate_wavefunction
    2657              : 
    2658              : ! **************************************************************************************************
    2659              : !> \brief maps a given function on the grid
    2660              : !> \param vector ...
    2661              : !> \param rho ...
    2662              : !> \param rho_gspace ...
    2663              : !> \param atomic_kind_set ...
    2664              : !> \param qs_kind_set ...
    2665              : !> \param cell ...
    2666              : !> \param particle_set ...
    2667              : !> \param pw_env ...
    2668              : !> \param eps_rho_rspace ...
    2669              : !> \param basis_type ...
    2670              : !> \par History
    2671              : !>      08.2002 created [Joost VandeVondele]
    2672              : !>      03.2006 made independent of qs_env [Joost VandeVondele]
    2673              : !>      08.2024 specialized version from calculate_wavefunction [JGH]
    2674              : !> \notes
    2675              : !>      modified calculate_rho_elec, should write the wavefunction represented by vector
    2676              : !>      it's presumably dominated by the FFT and the rs->pw and back routines
    2677              : ! **************************************************************************************************
    2678        39624 :    SUBROUTINE collocate_function(vector, rho, rho_gspace, &
    2679              :                                  atomic_kind_set, qs_kind_set, cell, particle_set, pw_env, &
    2680              :                                  eps_rho_rspace, basis_type)
    2681              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: vector
    2682              :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                :: rho
    2683              :       TYPE(pw_c1d_gs_type), INTENT(INOUT)                :: rho_gspace
    2684              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    2685              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    2686              :       TYPE(cell_type), POINTER                           :: cell
    2687              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    2688              :       TYPE(pw_env_type), POINTER                         :: pw_env
    2689              :       REAL(KIND=dp), INTENT(IN)                          :: eps_rho_rspace
    2690              :       CHARACTER(LEN=*), INTENT(IN), OPTIONAL             :: basis_type
    2691              : 
    2692              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'collocate_function'
    2693              : 
    2694              :       CHARACTER(LEN=default_string_length)               :: my_basis_type
    2695              :       INTEGER :: group_size, handle, i, iatom, igrid_level, ikind, ipgf, iset, maxco, maxsgf_set, &
    2696              :                  my_pos, na1, na2, natom, ncoa, nseta, offset, sgfa
    2697        19812 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: where_is_the_point
    2698        19812 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, npgfa, nsgfa
    2699        19812 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa
    2700              :       REAL(KIND=dp)                                      :: dab, radius, scale
    2701              :       REAL(KIND=dp), DIMENSION(3)                        :: ra
    2702        19812 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: pab, sphi_a, work, zeta
    2703              :       TYPE(gridlevel_info_type), POINTER                 :: gridlevel_info
    2704              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
    2705              :       TYPE(mp_comm_type)                                 :: group
    2706        19812 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
    2707        19812 :       TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:)    :: mgrid_gspace
    2708        19812 :       TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:)    :: mgrid_rspace
    2709        19812 :       TYPE(realspace_grid_type), DIMENSION(:), POINTER   :: rs_rho
    2710              : 
    2711        19812 :       CALL timeset(routineN, handle)
    2712              : 
    2713        19812 :       IF (PRESENT(basis_type)) THEN
    2714        18154 :          my_basis_type = basis_type
    2715              :       ELSE
    2716         1658 :          my_basis_type = "ORB"
    2717              :       END IF
    2718              : 
    2719        19812 :       NULLIFY (orb_basis_set, pab, work, la_max, la_min, &
    2720        19812 :                npgfa, nsgfa, sphi_a, zeta, first_sgfa, rs_rho, pw_pools)
    2721              : 
    2722              :       ! *** set up the pw multi-grids
    2723        19812 :       CPASSERT(ASSOCIATED(pw_env))
    2724              :       CALL pw_env_get(pw_env, rs_grids=rs_rho, pw_pools=pw_pools, &
    2725        19812 :                       gridlevel_info=gridlevel_info)
    2726              : 
    2727        19812 :       CALL pw_pools_create_pws(pw_pools, mgrid_gspace)
    2728        19812 :       CALL pw_pools_create_pws(pw_pools, mgrid_rspace)
    2729              : 
    2730              :       ! *** set up rs multi-grids
    2731        98784 :       DO igrid_level = 1, gridlevel_info%ngrid_levels
    2732        98784 :          CALL rs_grid_zero(rs_rho(igrid_level))
    2733              :       END DO
    2734              : 
    2735              :       !   *** Allocate work storage ***
    2736        19812 :       CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
    2737              :       CALL get_qs_kind_set(qs_kind_set, &
    2738              :                            maxco=maxco, &
    2739              :                            maxsgf_set=maxsgf_set, &
    2740        19812 :                            basis_type=my_basis_type)
    2741              : 
    2742        59436 :       ALLOCATE (pab(maxco, 1))
    2743        39624 :       ALLOCATE (work(maxco, 1))
    2744              : 
    2745        19812 :       offset = 0
    2746        19812 :       group = mgrid_rspace(1)%pw_grid%para%group
    2747        19812 :       my_pos = mgrid_rspace(1)%pw_grid%para%group%mepos
    2748        19812 :       group_size = mgrid_rspace(1)%pw_grid%para%group%num_pe
    2749        59436 :       ALLOCATE (where_is_the_point(0:group_size - 1))
    2750              : 
    2751        82112 :       DO iatom = 1, natom
    2752        62300 :          ikind = particle_set(iatom)%atomic_kind%kind_number
    2753        62300 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=my_basis_type)
    2754              :          CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
    2755              :                                 first_sgf=first_sgfa, &
    2756              :                                 lmax=la_max, &
    2757              :                                 lmin=la_min, &
    2758              :                                 npgf=npgfa, &
    2759              :                                 nset=nseta, &
    2760              :                                 nsgf_set=nsgfa, &
    2761              :                                 sphi=sphi_a, &
    2762        62300 :                                 zet=zeta)
    2763        62300 :          ra(:) = pbc(particle_set(iatom)%r, cell)
    2764        62300 :          dab = 0.0_dp
    2765              : 
    2766       689795 :          DO iset = 1, nseta
    2767              : 
    2768       545383 :             ncoa = npgfa(iset)*ncoset(la_max(iset))
    2769       545383 :             sgfa = first_sgfa(1, iset)
    2770              : 
    2771      2108432 :             DO i = 1, nsgfa(iset)
    2772      2108432 :                work(i, 1) = vector(offset + i)
    2773              :             END DO
    2774              : 
    2775              :             CALL dgemm("N", "N", ncoa, 1, nsgfa(iset), &
    2776              :                        1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
    2777              :                        work(1, 1), SIZE(work, 1), &
    2778       545383 :                        0.0_dp, pab(1, 1), SIZE(pab, 1))
    2779              : 
    2780      1116759 :             DO ipgf = 1, npgfa(iset)
    2781              : 
    2782       571376 :                na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
    2783       571376 :                na2 = ipgf*ncoset(la_max(iset))
    2784              : 
    2785       571376 :                scale = 1.0_dp
    2786       571376 :                igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
    2787              : 
    2788      1116759 :                IF (map_gaussian_here(rs_rho(igrid_level), cell%h_inv, ra, offset, group_size, my_pos)) THEN
    2789              :                   radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
    2790              :                                                     lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
    2791              :                                                     zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
    2792       522063 :                                                     prefactor=1.0_dp, cutoff=1.0_dp)
    2793              : 
    2794              :                   CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), la_min(iset), &
    2795              :                                              0, 0.0_dp, 0, &
    2796              :                                              ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
    2797              :                                              scale, pab, na1 - 1, 0, rs_rho(igrid_level), &
    2798       522063 :                                              radius=radius, ga_gb_function=GRID_FUNC_AB)
    2799              :                END IF
    2800              : 
    2801              :             END DO
    2802              : 
    2803       607683 :             offset = offset + nsgfa(iset)
    2804              : 
    2805              :          END DO
    2806              : 
    2807              :       END DO
    2808              : 
    2809        98784 :       DO igrid_level = 1, gridlevel_info%ngrid_levels
    2810              :          CALL transfer_rs2pw(rs_rho(igrid_level), &
    2811        98784 :                              mgrid_rspace(igrid_level))
    2812              :       END DO
    2813              : 
    2814        19812 :       CALL pw_zero(rho_gspace)
    2815        98784 :       DO igrid_level = 1, gridlevel_info%ngrid_levels
    2816              :          CALL pw_transfer(mgrid_rspace(igrid_level), &
    2817        78972 :                           mgrid_gspace(igrid_level))
    2818        98784 :          CALL pw_axpy(mgrid_gspace(igrid_level), rho_gspace)
    2819              :       END DO
    2820              : 
    2821        19812 :       CALL pw_transfer(rho_gspace, rho)
    2822              : 
    2823              :       ! Release work storage
    2824        19812 :       DEALLOCATE (pab)
    2825        19812 :       DEALLOCATE (work)
    2826              : 
    2827              :       ! give back the pw multi-grids
    2828        19812 :       CALL pw_pools_give_back_pws(pw_pools, mgrid_gspace)
    2829        19812 :       CALL pw_pools_give_back_pws(pw_pools, mgrid_rspace)
    2830              : 
    2831        19812 :       CALL timestop(handle)
    2832              : 
    2833        99060 :    END SUBROUTINE collocate_function
    2834              : 
    2835              : END MODULE qs_collocate_density
        

Generated by: LCOV version 2.0-1