LCOV - code coverage report
Current view: top level - src - qs_integrate_potential_single.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 83.4 % 571 476
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 7 7

            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 Build up the plane wave density by collocating the primitive Gaussian
      10              : !>      functions (pgf).
      11              : !> \par History
      12              : !>      Joost VandeVondele (02.2002)
      13              : !>            1) rewrote collocate_pgf for increased accuracy and speed
      14              : !>            2) collocate_core hack for PGI compiler
      15              : !>            3) added multiple grid feature
      16              : !>            4) new way to go over the grid
      17              : !>      Joost VandeVondele (05.2002)
      18              : !>            1) prelim. introduction of the real space grid type
      19              : !>      JGH [30.08.02] multigrid arrays independent from potential
      20              : !>      JGH [17.07.03] distributed real space code
      21              : !>      JGH [23.11.03] refactoring and new loop ordering
      22              : !>      JGH [04.12.03] OpneMP parallelization of main loops
      23              : !>      Joost VandeVondele (12.2003)
      24              : !>           1) modified to compute tau
      25              : !>      Joost removed incremental build feature
      26              : !>      Joost introduced map consistent
      27              : !>      Rewrote grid integration/collocation routines, [Joost VandeVondele,03.2007]
      28              : !> \author Matthias Krack (03.04.2001)
      29              : ! **************************************************************************************************
      30              : MODULE qs_integrate_potential_single
      31              :    USE ao_util,                         ONLY: exp_radius_very_extended
      32              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      33              :                                               get_atomic_kind,&
      34              :                                               get_atomic_kind_set
      35              :    USE atprop_types,                    ONLY: atprop_array_init,&
      36              :                                               atprop_type
      37              :    USE basis_set_types,                 ONLY: get_gto_basis_set,&
      38              :                                               gto_basis_set_type
      39              :    USE cell_types,                      ONLY: cell_type,&
      40              :                                               pbc
      41              :    USE cp_control_types,                ONLY: dft_control_type
      42              :    USE cp_dbcsr_api,                    ONLY: dbcsr_get_block_p,&
      43              :                                               dbcsr_type
      44              :    USE external_potential_types,        ONLY: get_potential,&
      45              :                                               gth_potential_type
      46              :    USE gaussian_gridlevels,             ONLY: gaussian_gridlevel,&
      47              :                                               gridlevel_info_type
      48              :    USE grid_api,                        ONLY: integrate_pgf_product
      49              :    USE kinds,                           ONLY: dp
      50              :    USE lri_environment_types,           ONLY: lri_kind_type
      51              :    USE memory_utilities,                ONLY: reallocate
      52              :    USE message_passing,                 ONLY: mp_para_env_type
      53              :    USE orbital_pointers,                ONLY: coset,&
      54              :                                               ncoset
      55              :    USE particle_types,                  ONLY: particle_type
      56              :    USE pw_env_types,                    ONLY: pw_env_get,&
      57              :                                               pw_env_type
      58              :    USE pw_types,                        ONLY: pw_r3d_rs_type
      59              :    USE qs_environment_types,            ONLY: get_qs_env,&
      60              :                                               qs_environment_type
      61              :    USE qs_force_types,                  ONLY: qs_force_type
      62              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      63              :                                               get_qs_kind_set,&
      64              :                                               qs_kind_type
      65              :    USE realspace_grid_types,            ONLY: map_gaussian_here,&
      66              :                                               realspace_grid_type,&
      67              :                                               rs_grid_zero,&
      68              :                                               transfer_pw2rs
      69              :    USE rs_pw_interface,                 ONLY: potential_pw2rs
      70              :    USE virial_types,                    ONLY: virial_type
      71              : #include "./base/base_uses.f90"
      72              : 
      73              :    IMPLICIT NONE
      74              : 
      75              :    PRIVATE
      76              : 
      77              :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
      78              : 
      79              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_integrate_potential_single'
      80              : 
      81              : ! *** Public subroutines ***
      82              : ! *** Don't include this routines directly, use the interface to
      83              : ! *** qs_integrate_potential
      84              : 
      85              :    PUBLIC :: integrate_v_rspace_one_center, &
      86              :              integrate_v_rspace_diagonal, &
      87              :              integrate_v_core_rspace, &
      88              :              integrate_v_gaussian_rspace, &
      89              :              integrate_function, &
      90              :              integrate_ppl_rspace, &
      91              :              integrate_rho_nlcc
      92              : 
      93              : CONTAINS
      94              : 
      95              : ! **************************************************************************************************
      96              : !> \brief computes the forces/virial due to the local pseudopotential
      97              : !> \param rho_rspace ...
      98              : !> \param qs_env ...
      99              : ! **************************************************************************************************
     100           14 :    SUBROUTINE integrate_ppl_rspace(rho_rspace, qs_env)
     101              :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: rho_rspace
     102              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     103              : 
     104              :       CHARACTER(len=*), PARAMETER :: routineN = 'integrate_ppl_rspace'
     105              : 
     106              :       INTEGER                                            :: atom_a, handle, iatom, ikind, j, lppl, &
     107              :                                                             n, natom_of_kind, ni, npme
     108           14 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list, cores
     109              :       LOGICAL                                            :: use_virial
     110              :       REAL(KIND=dp)                                      :: alpha, eps_rho_rspace, radius
     111              :       REAL(KIND=dp), DIMENSION(3)                        :: force_a, force_b, ra
     112              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: my_virial_a, my_virial_b
     113           14 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: cexp_ppl
     114           14 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: hab, pab
     115           14 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     116              :       TYPE(cell_type), POINTER                           :: cell
     117              :       TYPE(dft_control_type), POINTER                    :: dft_control
     118              :       TYPE(gth_potential_type), POINTER                  :: gth_potential
     119           14 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     120              :       TYPE(pw_env_type), POINTER                         :: pw_env
     121           14 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     122           14 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     123              :       TYPE(realspace_grid_type), POINTER                 :: rs_v
     124              :       TYPE(virial_type), POINTER                         :: virial
     125              : 
     126           14 :       CALL timeset(routineN, handle)
     127              : 
     128           14 :       NULLIFY (pw_env, cores)
     129              : 
     130           14 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
     131           14 :       CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
     132              : 
     133           14 :       CALL transfer_pw2rs(rs_v, rho_rspace)
     134              : 
     135              :       CALL get_qs_env(qs_env=qs_env, &
     136              :                       atomic_kind_set=atomic_kind_set, &
     137              :                       qs_kind_set=qs_kind_set, &
     138              :                       cell=cell, &
     139              :                       dft_control=dft_control, &
     140              :                       particle_set=particle_set, &
     141              :                       pw_env=pw_env, &
     142           14 :                       force=force, virial=virial)
     143              : 
     144           14 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     145              : 
     146           14 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
     147              : 
     148           34 :       DO ikind = 1, SIZE(atomic_kind_set)
     149              : 
     150           20 :          CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
     151           20 :          CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential)
     152              : 
     153           20 :          IF (.NOT. ASSOCIATED(gth_potential)) CYCLE
     154           20 :          CALL get_potential(potential=gth_potential, alpha_ppl=alpha, nexp_ppl=lppl, cexp_ppl=cexp_ppl)
     155              : 
     156           20 :          IF (lppl <= 0) CYCLE
     157              : 
     158           20 :          ni = ncoset(2*lppl - 2)
     159          100 :          ALLOCATE (hab(ni, 1), pab(ni, 1))
     160          240 :          pab = 0._dp
     161              : 
     162           20 :          CALL reallocate(cores, 1, natom_of_kind)
     163           20 :          npme = 0
     164           70 :          cores = 0
     165              : 
     166              :          ! prepare core function
     167           60 :          DO j = 1, lppl
     168           20 :             SELECT CASE (j)
     169              :             CASE (1)
     170           20 :                pab(1, 1) = cexp_ppl(1)
     171              :             CASE (2)
     172           20 :                n = coset(2, 0, 0)
     173           20 :                pab(n, 1) = cexp_ppl(2)
     174           20 :                n = coset(0, 2, 0)
     175           20 :                pab(n, 1) = cexp_ppl(2)
     176           20 :                n = coset(0, 0, 2)
     177           20 :                pab(n, 1) = cexp_ppl(2)
     178              :             CASE (3)
     179            0 :                n = coset(4, 0, 0)
     180            0 :                pab(n, 1) = cexp_ppl(3)
     181            0 :                n = coset(0, 4, 0)
     182            0 :                pab(n, 1) = cexp_ppl(3)
     183            0 :                n = coset(0, 0, 4)
     184            0 :                pab(n, 1) = cexp_ppl(3)
     185            0 :                n = coset(2, 2, 0)
     186            0 :                pab(n, 1) = 2._dp*cexp_ppl(3)
     187            0 :                n = coset(2, 0, 2)
     188            0 :                pab(n, 1) = 2._dp*cexp_ppl(3)
     189            0 :                n = coset(0, 2, 2)
     190            0 :                pab(n, 1) = 2._dp*cexp_ppl(3)
     191              :             CASE (4)
     192            0 :                n = coset(6, 0, 0)
     193            0 :                pab(n, 1) = cexp_ppl(4)
     194            0 :                n = coset(0, 6, 0)
     195            0 :                pab(n, 1) = cexp_ppl(4)
     196            0 :                n = coset(0, 0, 6)
     197            0 :                pab(n, 1) = cexp_ppl(4)
     198            0 :                n = coset(4, 2, 0)
     199            0 :                pab(n, 1) = 3._dp*cexp_ppl(4)
     200            0 :                n = coset(4, 0, 2)
     201            0 :                pab(n, 1) = 3._dp*cexp_ppl(4)
     202            0 :                n = coset(2, 4, 0)
     203            0 :                pab(n, 1) = 3._dp*cexp_ppl(4)
     204            0 :                n = coset(2, 0, 4)
     205            0 :                pab(n, 1) = 3._dp*cexp_ppl(4)
     206            0 :                n = coset(0, 4, 2)
     207            0 :                pab(n, 1) = 3._dp*cexp_ppl(4)
     208            0 :                n = coset(0, 2, 4)
     209            0 :                pab(n, 1) = 3._dp*cexp_ppl(4)
     210            0 :                n = coset(2, 2, 2)
     211            0 :                pab(n, 1) = 6._dp*cexp_ppl(4)
     212              :             CASE DEFAULT
     213           40 :                CPABORT("")
     214              :             END SELECT
     215              :          END DO
     216              : 
     217           70 :          DO iatom = 1, natom_of_kind
     218           50 :             atom_a = atom_list(iatom)
     219           50 :             ra(:) = pbc(particle_set(atom_a)%r, cell)
     220           70 :             IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
     221              :                ! replicated realspace grid, split the atoms up between procs
     222           50 :                IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
     223           25 :                   npme = npme + 1
     224           25 :                   cores(npme) = iatom
     225              :                END IF
     226              :             ELSE
     227            0 :                npme = npme + 1
     228            0 :                cores(npme) = iatom
     229              :             END IF
     230              :          END DO
     231              : 
     232           45 :          DO j = 1, npme
     233              : 
     234           25 :             iatom = cores(j)
     235           25 :             atom_a = atom_list(iatom)
     236           25 :             ra(:) = pbc(particle_set(atom_a)%r, cell)
     237          275 :             hab(:, 1) = 0.0_dp
     238           25 :             force_a(:) = 0.0_dp
     239           25 :             force_b(:) = 0.0_dp
     240           25 :             IF (use_virial) THEN
     241            0 :                my_virial_a = 0.0_dp
     242            0 :                my_virial_b = 0.0_dp
     243              :             END IF
     244           25 :             ni = 2*lppl - 2
     245              : 
     246              :             radius = exp_radius_very_extended(la_min=0, la_max=ni, lb_min=0, lb_max=0, &
     247              :                                               ra=ra, rb=ra, rp=ra, &
     248              :                                               zetp=alpha, eps=eps_rho_rspace, &
     249              :                                               pab=pab, o1=0, o2=0, &  ! without map_consistent
     250           25 :                                               prefactor=1.0_dp, cutoff=1.0_dp)
     251              : 
     252              :             CALL integrate_pgf_product(ni, alpha, 0, &
     253              :                                        0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
     254              :                                        rs_v, hab, pab=pab, o1=0, o2=0, &
     255              :                                        radius=radius, &
     256              :                                        calculate_forces=.TRUE., force_a=force_a, &
     257              :                                        force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, &
     258           25 :                                        my_virial_b=my_virial_b, use_subpatch=.TRUE., subpatch_pattern=0)
     259              : 
     260              :             force(ikind)%gth_ppl(:, iatom) = &
     261          100 :                force(ikind)%gth_ppl(:, iatom) + force_a(:)*rho_rspace%pw_grid%dvol
     262              : 
     263           45 :             IF (use_virial) THEN
     264            0 :                virial%pv_ppl = virial%pv_ppl + my_virial_a*rho_rspace%pw_grid%dvol
     265            0 :                virial%pv_virial = virial%pv_virial + my_virial_a*rho_rspace%pw_grid%dvol
     266            0 :                CPABORT("Virial not debuged for CORE_PPL")
     267              :             END IF
     268              :          END DO
     269              : 
     270           74 :          DEALLOCATE (hab, pab)
     271              : 
     272              :       END DO
     273              : 
     274           14 :       DEALLOCATE (cores)
     275              : 
     276           14 :       CALL timestop(handle)
     277              : 
     278           14 :    END SUBROUTINE integrate_ppl_rspace
     279              : 
     280              : ! **************************************************************************************************
     281              : !> \brief computes the forces/virial due to the nlcc pseudopotential
     282              : !> \param rho_rspace ...
     283              : !> \param qs_env ...
     284              : ! **************************************************************************************************
     285           30 :    SUBROUTINE integrate_rho_nlcc(rho_rspace, qs_env)
     286              :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: rho_rspace
     287              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     288              : 
     289              :       CHARACTER(len=*), PARAMETER :: routineN = 'integrate_rho_nlcc'
     290              : 
     291              :       INTEGER                                            :: atom_a, handle, iatom, iexp_nlcc, ikind, &
     292              :                                                             ithread, j, n, natom, nc, nexp_nlcc, &
     293              :                                                             ni, npme, nthread
     294           30 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list, cores, nct_nlcc
     295              :       LOGICAL                                            :: nlcc, use_virial
     296              :       REAL(KIND=dp)                                      :: alpha, eps_rho_rspace, radius
     297              :       REAL(KIND=dp), DIMENSION(3)                        :: force_a, force_b, ra
     298              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: my_virial_a, my_virial_b
     299           30 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: alpha_nlcc
     300           30 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: cval_nlcc, hab, pab
     301           30 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     302              :       TYPE(cell_type), POINTER                           :: cell
     303              :       TYPE(dft_control_type), POINTER                    :: dft_control
     304              :       TYPE(gth_potential_type), POINTER                  :: gth_potential
     305           30 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     306              :       TYPE(pw_env_type), POINTER                         :: pw_env
     307           30 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     308           30 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     309              :       TYPE(realspace_grid_type), POINTER                 :: rs_v
     310              :       TYPE(virial_type), POINTER                         :: virial
     311              : 
     312           30 :       CALL timeset(routineN, handle)
     313              : 
     314           30 :       NULLIFY (pw_env, cores)
     315              : 
     316           30 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
     317           30 :       CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
     318              : 
     319           30 :       CALL transfer_pw2rs(rs_v, rho_rspace)
     320              : 
     321              :       CALL get_qs_env(qs_env=qs_env, &
     322              :                       atomic_kind_set=atomic_kind_set, &
     323              :                       qs_kind_set=qs_kind_set, &
     324              :                       cell=cell, &
     325              :                       dft_control=dft_control, &
     326              :                       particle_set=particle_set, &
     327              :                       pw_env=pw_env, &
     328           30 :                       force=force, virial=virial)
     329              : 
     330           30 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     331              : 
     332           30 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
     333              : 
     334           74 :       DO ikind = 1, SIZE(atomic_kind_set)
     335              : 
     336           44 :          CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
     337           44 :          CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential)
     338              : 
     339           44 :          IF (.NOT. ASSOCIATED(gth_potential)) CYCLE
     340              :          CALL get_potential(potential=gth_potential, nlcc_present=nlcc, nexp_nlcc=nexp_nlcc, &
     341           44 :                             alpha_nlcc=alpha_nlcc, nct_nlcc=nct_nlcc, cval_nlcc=cval_nlcc)
     342              : 
     343           44 :          IF (.NOT. nlcc) CYCLE
     344              : 
     345          202 :          DO iexp_nlcc = 1, nexp_nlcc
     346              : 
     347           42 :             alpha = alpha_nlcc(iexp_nlcc)
     348           42 :             nc = nct_nlcc(iexp_nlcc)
     349              : 
     350           42 :             ni = ncoset(2*nc - 2)
     351              : 
     352           42 :             nthread = 1
     353           42 :             ithread = 0
     354              : 
     355          168 :             ALLOCATE (hab(ni, 1), pab(ni, 1))
     356          270 :             pab = 0._dp
     357              : 
     358           42 :             CALL reallocate(cores, 1, natom)
     359           42 :             npme = 0
     360          172 :             cores = 0
     361              : 
     362              :             ! prepare core function
     363          100 :             DO j = 1, nc
     364           42 :                SELECT CASE (j)
     365              :                CASE (1)
     366           42 :                   pab(1, 1) = cval_nlcc(1, iexp_nlcc)
     367              :                CASE (2)
     368           16 :                   n = coset(2, 0, 0)
     369           16 :                   pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
     370           16 :                   n = coset(0, 2, 0)
     371           16 :                   pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
     372           16 :                   n = coset(0, 0, 2)
     373           16 :                   pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
     374              :                CASE (3)
     375            0 :                   n = coset(4, 0, 0)
     376            0 :                   pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
     377            0 :                   n = coset(0, 4, 0)
     378            0 :                   pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
     379            0 :                   n = coset(0, 0, 4)
     380            0 :                   pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
     381            0 :                   n = coset(2, 2, 0)
     382            0 :                   pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
     383            0 :                   n = coset(2, 0, 2)
     384            0 :                   pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
     385            0 :                   n = coset(0, 2, 2)
     386            0 :                   pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
     387              :                CASE (4)
     388            0 :                   n = coset(6, 0, 0)
     389            0 :                   pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
     390            0 :                   n = coset(0, 6, 0)
     391            0 :                   pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
     392            0 :                   n = coset(0, 0, 6)
     393            0 :                   pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
     394            0 :                   n = coset(4, 2, 0)
     395            0 :                   pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
     396            0 :                   n = coset(4, 0, 2)
     397            0 :                   pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
     398            0 :                   n = coset(2, 4, 0)
     399            0 :                   pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
     400            0 :                   n = coset(2, 0, 4)
     401            0 :                   pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
     402            0 :                   n = coset(0, 4, 2)
     403            0 :                   pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
     404            0 :                   n = coset(0, 2, 4)
     405            0 :                   pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
     406            0 :                   n = coset(2, 2, 2)
     407            0 :                   pab(n, 1) = 6._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
     408              :                CASE DEFAULT
     409           58 :                   CPABORT("")
     410              :                END SELECT
     411              :             END DO
     412           42 :             IF (dft_control%nspins == 2) pab = pab*0.5_dp
     413              : 
     414          172 :             DO iatom = 1, natom
     415          130 :                atom_a = atom_list(iatom)
     416          130 :                ra(:) = pbc(particle_set(atom_a)%r, cell)
     417          172 :                IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
     418              :                   ! replicated realspace grid, split the atoms up between procs
     419          130 :                   IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
     420           65 :                      npme = npme + 1
     421           65 :                      cores(npme) = iatom
     422              :                   END IF
     423              :                ELSE
     424            0 :                   npme = npme + 1
     425            0 :                   cores(npme) = iatom
     426              :                END IF
     427              :             END DO
     428              : 
     429          107 :             DO j = 1, npme
     430              : 
     431           65 :                iatom = cores(j)
     432           65 :                atom_a = atom_list(iatom)
     433           65 :                ra(:) = pbc(particle_set(atom_a)%r, cell)
     434          274 :                hab(:, 1) = 0.0_dp
     435           65 :                force_a(:) = 0.0_dp
     436           65 :                force_b(:) = 0.0_dp
     437           65 :                IF (use_virial) THEN
     438           48 :                   my_virial_a = 0.0_dp
     439           48 :                   my_virial_b = 0.0_dp
     440              :                END IF
     441           65 :                ni = 2*nc - 2
     442              : 
     443              :                radius = exp_radius_very_extended(la_min=0, la_max=ni, lb_min=0, lb_max=0, &
     444              :                                                  ra=ra, rb=ra, rp=ra, &
     445              :                                                  zetp=1/(2*alpha**2), eps=eps_rho_rspace, &
     446              :                                                  pab=pab, o1=0, o2=0, &  ! without map_consistent
     447           65 :                                                  prefactor=1.0_dp, cutoff=1.0_dp)
     448              : 
     449              :                CALL integrate_pgf_product(ni, 1/(2*alpha**2), 0, &
     450              :                                           0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
     451              :                                           rs_v, hab, pab=pab, o1=0, o2=0, &
     452              :                                           radius=radius, &
     453              :                                           calculate_forces=.TRUE., force_a=force_a, &
     454              :                                           force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, &
     455           65 :                                           my_virial_b=my_virial_b, use_subpatch=.TRUE., subpatch_pattern=0)
     456              : 
     457              :                force(ikind)%gth_nlcc(:, iatom) = &
     458          260 :                   force(ikind)%gth_nlcc(:, iatom) + force_a(:)*rho_rspace%pw_grid%dvol
     459              : 
     460          107 :                IF (use_virial) THEN
     461          624 :                   virial%pv_nlcc = virial%pv_nlcc + my_virial_a*rho_rspace%pw_grid%dvol
     462          624 :                   virial%pv_virial = virial%pv_virial + my_virial_a*rho_rspace%pw_grid%dvol
     463              :                END IF
     464              :             END DO
     465              : 
     466           86 :             DEALLOCATE (hab, pab)
     467              : 
     468              :          END DO
     469              : 
     470              :       END DO
     471              : 
     472           30 :       DEALLOCATE (cores)
     473              : 
     474           30 :       CALL timestop(handle)
     475              : 
     476           30 :    END SUBROUTINE integrate_rho_nlcc
     477              : 
     478              : ! **************************************************************************************************
     479              : !> \brief computes the forces/virial due to the ionic cores with a potential on
     480              : !>      grid
     481              : !> \param v_rspace ...
     482              : !> \param qs_env ...
     483              : !> \param atecc ...
     484              : ! **************************************************************************************************
     485         7443 :    SUBROUTINE integrate_v_core_rspace(v_rspace, qs_env, atecc)
     486              :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: v_rspace
     487              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     488              :       REAL(KIND=dp), DIMENSION(:), OPTIONAL              :: atecc
     489              : 
     490              :       CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_core_rspace'
     491              : 
     492              :       INTEGER                                            :: atom_a, handle, iatom, ikind, j, natom, &
     493              :                                                             natom_of_kind, npme
     494         7443 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list, cores
     495              :       LOGICAL                                            :: paw_atom, skip_fcore, use_virial
     496              :       REAL(KIND=dp)                                      :: alpha_core_charge, ccore_charge, &
     497              :                                                             eps_rho_rspace, radius
     498              :       REAL(KIND=dp), DIMENSION(3)                        :: force_a, force_b, ra
     499              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: my_virial_a, my_virial_b
     500         7443 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: hab, pab
     501         7443 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     502              :       TYPE(atprop_type), POINTER                         :: atprop
     503              :       TYPE(cell_type), POINTER                           :: cell
     504              :       TYPE(dft_control_type), POINTER                    :: dft_control
     505         7443 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     506              :       TYPE(pw_env_type), POINTER                         :: pw_env
     507         7443 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     508         7443 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     509              :       TYPE(realspace_grid_type), POINTER                 :: rs_v
     510              :       TYPE(virial_type), POINTER                         :: virial
     511              : 
     512         7443 :       CALL timeset(routineN, handle)
     513         7443 :       NULLIFY (virial, force, atprop, dft_control)
     514              : 
     515         7443 :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
     516              : 
     517              :       !If gapw, check for gpw kinds
     518         7443 :       skip_fcore = .FALSE.
     519         7443 :       IF (dft_control%qs_control%gapw) THEN
     520          510 :          IF (.NOT. dft_control%qs_control%gapw_control%nopaw_as_gpw) skip_fcore = .TRUE.
     521              :       END IF
     522              : 
     523              :       IF (.NOT. skip_fcore) THEN
     524         6989 :          NULLIFY (pw_env)
     525         6989 :          ALLOCATE (cores(1))
     526         6989 :          ALLOCATE (hab(1, 1))
     527         6989 :          ALLOCATE (pab(1, 1))
     528              : 
     529         6989 :          CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
     530         6989 :          CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
     531              : 
     532         6989 :          CALL transfer_pw2rs(rs_v, v_rspace)
     533              : 
     534              :          CALL get_qs_env(qs_env=qs_env, &
     535              :                          atomic_kind_set=atomic_kind_set, &
     536              :                          qs_kind_set=qs_kind_set, &
     537              :                          cell=cell, &
     538              :                          dft_control=dft_control, &
     539              :                          particle_set=particle_set, &
     540              :                          pw_env=pw_env, &
     541              :                          force=force, &
     542              :                          virial=virial, &
     543         6989 :                          atprop=atprop)
     544              : 
     545              :          ! atomic energy contributions
     546         6989 :          natom = SIZE(particle_set)
     547         6989 :          IF (ASSOCIATED(atprop)) THEN
     548         6989 :             CALL atprop_array_init(atprop%ateb, natom)
     549              :          END IF
     550              : 
     551         6989 :          use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     552              : 
     553         6989 :          eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
     554              : 
     555        19273 :          DO ikind = 1, SIZE(atomic_kind_set)
     556              : 
     557        12284 :             CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
     558              :             CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom, &
     559              :                              alpha_core_charge=alpha_core_charge, &
     560        12284 :                              ccore_charge=ccore_charge)
     561              : 
     562        12284 :             IF (dft_control%qs_control%gapw .AND. paw_atom) CYCLE
     563              : 
     564        12222 :             pab(1, 1) = -ccore_charge
     565        12222 :             IF (alpha_core_charge == 0.0_dp .OR. pab(1, 1) == 0.0_dp) CYCLE
     566              : 
     567        12186 :             CALL reallocate(cores, 1, natom_of_kind)
     568        12186 :             npme = 0
     569        36889 :             cores = 0
     570              : 
     571        36889 :             DO iatom = 1, natom_of_kind
     572        24703 :                atom_a = atom_list(iatom)
     573        24703 :                ra(:) = pbc(particle_set(atom_a)%r, cell)
     574        36889 :                IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
     575              :                   ! replicated realspace grid, split the atoms up between procs
     576        23954 :                   IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
     577        11977 :                      npme = npme + 1
     578        11977 :                      cores(npme) = iatom
     579              :                   END IF
     580              :                ELSE
     581          749 :                   npme = npme + 1
     582          749 :                   cores(npme) = iatom
     583              :                END IF
     584              :             END DO
     585              : 
     586        44185 :             DO j = 1, npme
     587              : 
     588        12726 :                iatom = cores(j)
     589        12726 :                atom_a = atom_list(iatom)
     590        12726 :                ra(:) = pbc(particle_set(atom_a)%r, cell)
     591        12726 :                hab(1, 1) = 0.0_dp
     592        12726 :                force_a(:) = 0.0_dp
     593        12726 :                force_b(:) = 0.0_dp
     594        12726 :                IF (use_virial) THEN
     595         1600 :                   my_virial_a = 0.0_dp
     596         1600 :                   my_virial_b = 0.0_dp
     597              :                END IF
     598              : 
     599              :                radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
     600              :                                                  ra=ra, rb=ra, rp=ra, &
     601              :                                                  zetp=alpha_core_charge, eps=eps_rho_rspace, &
     602              :                                                  pab=pab, o1=0, o2=0, &  ! without map_consistent
     603        12726 :                                                  prefactor=1.0_dp, cutoff=1.0_dp)
     604              : 
     605              :                CALL integrate_pgf_product(0, alpha_core_charge, 0, &
     606              :                                           0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
     607              :                                           rs_v, hab, pab=pab, o1=0, o2=0, &
     608              :                                           radius=radius, &
     609              :                                           calculate_forces=.TRUE., force_a=force_a, &
     610              :                                           force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, &
     611        12726 :                                           my_virial_b=my_virial_b, use_subpatch=.TRUE., subpatch_pattern=0)
     612              : 
     613        12726 :                IF (ASSOCIATED(force)) THEN
     614        50516 :                   force(ikind)%rho_core(:, iatom) = force(ikind)%rho_core(:, iatom) + force_a(:)
     615              :                END IF
     616        12726 :                IF (use_virial) THEN
     617        20800 :                   virial%pv_ehartree = virial%pv_ehartree + my_virial_a
     618        20800 :                   virial%pv_virial = virial%pv_virial + my_virial_a
     619              :                END IF
     620        12726 :                IF (ASSOCIATED(atprop)) THEN
     621        12726 :                   atprop%ateb(atom_a) = atprop%ateb(atom_a) + 0.5_dp*hab(1, 1)*pab(1, 1)
     622              :                END IF
     623        25010 :                IF (PRESENT(atecc)) THEN
     624           47 :                   atecc(atom_a) = atecc(atom_a) + 0.5_dp*hab(1, 1)*pab(1, 1)
     625              :                END IF
     626              : 
     627              :             END DO
     628              : 
     629              :          END DO
     630              : 
     631         6989 :          DEALLOCATE (hab, pab, cores)
     632              : 
     633              :       END IF
     634              : 
     635         7443 :       CALL timestop(handle)
     636              : 
     637         7443 :    END SUBROUTINE integrate_v_core_rspace
     638              : 
     639              : ! **************************************************************************************************
     640              : !> \brief computes the overlap of a set of Gaussians with a potential on grid
     641              : !> \param v_rspace ...
     642              : !> \param qs_env ...
     643              : !> \param alpha ...
     644              : !> \param ccore ...
     645              : !> \param atecc ...
     646              : ! **************************************************************************************************
     647            2 :    SUBROUTINE integrate_v_gaussian_rspace(v_rspace, qs_env, alpha, ccore, atecc)
     648              :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: v_rspace
     649              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     650              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: alpha, ccore
     651              :       REAL(KIND=dp), DIMENSION(:)                        :: atecc
     652              : 
     653              :       CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_gaussian_rspace'
     654              : 
     655              :       INTEGER                                            :: atom_a, handle, iatom, ikind, j, natom, &
     656              :                                                             natom_of_kind, npme
     657            2 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list, cores
     658              :       REAL(KIND=dp)                                      :: alpha_core_charge, eps_rho_rspace, radius
     659              :       REAL(KIND=dp), DIMENSION(3)                        :: ra
     660            2 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: hab, pab
     661            2 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     662              :       TYPE(cell_type), POINTER                           :: cell
     663              :       TYPE(dft_control_type), POINTER                    :: dft_control
     664            2 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     665              :       TYPE(pw_env_type), POINTER                         :: pw_env
     666            2 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     667              :       TYPE(realspace_grid_type), POINTER                 :: rs_v
     668              : 
     669            2 :       CALL timeset(routineN, handle)
     670              : 
     671            2 :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
     672              : 
     673              :       !If gapw, check for gpw kinds
     674            2 :       CPASSERT(.NOT. dft_control%qs_control%gapw)
     675              : 
     676            2 :       NULLIFY (pw_env)
     677            2 :       ALLOCATE (cores(1))
     678            2 :       ALLOCATE (hab(1, 1))
     679            2 :       ALLOCATE (pab(1, 1))
     680              : 
     681            2 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
     682            2 :       CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
     683              : 
     684            2 :       CALL transfer_pw2rs(rs_v, v_rspace)
     685              : 
     686              :       CALL get_qs_env(qs_env=qs_env, &
     687              :                       atomic_kind_set=atomic_kind_set, &
     688              :                       qs_kind_set=qs_kind_set, &
     689              :                       cell=cell, &
     690              :                       dft_control=dft_control, &
     691              :                       particle_set=particle_set, &
     692            2 :                       pw_env=pw_env)
     693              : 
     694              :       ! atomic energy contributions
     695            2 :       natom = SIZE(particle_set)
     696            2 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
     697              : 
     698            6 :       DO ikind = 1, SIZE(atomic_kind_set)
     699              : 
     700            4 :          CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
     701            4 :          pab(1, 1) = -ccore(ikind)
     702            4 :          alpha_core_charge = alpha(ikind)
     703            4 :          IF (alpha_core_charge == 0.0_dp .OR. pab(1, 1) == 0.0_dp) CYCLE
     704              : 
     705            4 :          CALL reallocate(cores, 1, natom_of_kind)
     706            4 :          npme = 0
     707           10 :          cores = 0
     708              : 
     709           10 :          DO iatom = 1, natom_of_kind
     710            6 :             atom_a = atom_list(iatom)
     711            6 :             ra(:) = pbc(particle_set(atom_a)%r, cell)
     712           10 :             IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
     713              :                ! replicated realspace grid, split the atoms up between procs
     714            6 :                IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
     715            3 :                   npme = npme + 1
     716            3 :                   cores(npme) = iatom
     717              :                END IF
     718              :             ELSE
     719            0 :                npme = npme + 1
     720            0 :                cores(npme) = iatom
     721              :             END IF
     722              :          END DO
     723              : 
     724           13 :          DO j = 1, npme
     725              : 
     726            3 :             iatom = cores(j)
     727            3 :             atom_a = atom_list(iatom)
     728            3 :             ra(:) = pbc(particle_set(atom_a)%r, cell)
     729            3 :             hab(1, 1) = 0.0_dp
     730              : 
     731              :             radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
     732              :                                               ra=ra, rb=ra, rp=ra, &
     733              :                                               zetp=alpha_core_charge, eps=eps_rho_rspace, &
     734              :                                               pab=pab, o1=0, o2=0, &  ! without map_consistent
     735            3 :                                               prefactor=1.0_dp, cutoff=1.0_dp)
     736              : 
     737              :             CALL integrate_pgf_product(0, alpha_core_charge, 0, &
     738              :                                        0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
     739              :                                        rs_v, hab, pab=pab, o1=0, o2=0, &
     740              :                                        radius=radius, calculate_forces=.FALSE., &
     741            3 :                                        use_subpatch=.TRUE., subpatch_pattern=0)
     742            7 :             atecc(atom_a) = atecc(atom_a) + 0.5_dp*hab(1, 1)*pab(1, 1)
     743              : 
     744              :          END DO
     745              : 
     746              :       END DO
     747              : 
     748            2 :       DEALLOCATE (hab, pab, cores)
     749              : 
     750            2 :       CALL timestop(handle)
     751              : 
     752            2 :    END SUBROUTINE integrate_v_gaussian_rspace
     753              : ! **************************************************************************************************
     754              : !> \brief computes integrals of product of v_rspace times a one-center function
     755              : !>        required for LRIGPW
     756              : !> \param v_rspace ...
     757              : !> \param qs_env ...
     758              : !> \param int_res ...
     759              : !> \param calculate_forces ...
     760              : !> \param basis_type ...
     761              : !> \param atomlist ...
     762              : !> \author Dorothea Golze
     763              : ! **************************************************************************************************
     764          822 :    SUBROUTINE integrate_v_rspace_one_center(v_rspace, qs_env, int_res, &
     765          822 :                                             calculate_forces, basis_type, atomlist)
     766              :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: v_rspace
     767              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     768              :       TYPE(lri_kind_type), DIMENSION(:), POINTER         :: int_res
     769              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     770              :       CHARACTER(len=*), INTENT(IN)                       :: basis_type
     771              :       INTEGER, DIMENSION(:), OPTIONAL                    :: atomlist
     772              : 
     773              :       CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_rspace_one_center'
     774              : 
     775              :       INTEGER :: atom_a, group_size, handle, i, iatom, igrid_level, ikind, ipgf, iset, m1, maxco, &
     776              :          maxsgf_set, my_pos, na1, natom_of_kind, ncoa, nkind, nseta, offset, sgfa
     777          822 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list, la_max, la_min, npgfa, &
     778          822 :                                                             nsgf_seta
     779          822 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa
     780              :       LOGICAL                                            :: use_virial
     781          822 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: map_it
     782              :       REAL(KIND=dp)                                      :: eps_rho_rspace, radius
     783              :       REAL(KIND=dp), DIMENSION(3)                        :: force_a, force_b, ra
     784              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: my_virial_a, my_virial_b
     785          822 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a
     786          822 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: hab, pab, rpgfa, sphi_a, work_f, work_i, &
     787          822 :                                                             zeta
     788          822 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     789              :       TYPE(cell_type), POINTER                           :: cell
     790              :       TYPE(dft_control_type), POINTER                    :: dft_control
     791              :       TYPE(gridlevel_info_type), POINTER                 :: gridlevel_info
     792              :       TYPE(gto_basis_set_type), POINTER                  :: lri_basis_set
     793          822 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     794              :       TYPE(pw_env_type), POINTER                         :: pw_env
     795          822 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     796          822 :       TYPE(realspace_grid_type), DIMENSION(:), POINTER   :: rs_v
     797              :       TYPE(realspace_grid_type), POINTER                 :: rs_grid
     798              :       TYPE(virial_type), POINTER                         :: virial
     799              : 
     800          822 :       CALL timeset(routineN, handle)
     801              : 
     802          822 :       NULLIFY (atomic_kind_set, qs_kind_set, atom_list, cell, dft_control, &
     803          822 :                first_sgfa, gridlevel_info, hab, la_max, la_min, lri_basis_set, &
     804          822 :                npgfa, nsgf_seta, pab, particle_set, pw_env, rpgfa, &
     805          822 :                rs_grid, rs_v, virial, set_radius_a, sphi_a, work_f, &
     806          822 :                work_i, zeta)
     807              : 
     808          822 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
     809              : 
     810          822 :       CALL pw_env_get(pw_env, rs_grids=rs_v)
     811         4074 :       DO i = 1, SIZE(rs_v)
     812         4074 :          CALL rs_grid_zero(rs_v(i))
     813              :       END DO
     814              : 
     815          822 :       gridlevel_info => pw_env%gridlevel_info
     816              : 
     817          822 :       CALL potential_pw2rs(rs_v, v_rspace, pw_env)
     818              : 
     819              :       CALL get_qs_env(qs_env=qs_env, &
     820              :                       atomic_kind_set=atomic_kind_set, &
     821              :                       qs_kind_set=qs_kind_set, &
     822              :                       cell=cell, &
     823              :                       dft_control=dft_control, &
     824              :                       nkind=nkind, &
     825              :                       particle_set=particle_set, &
     826              :                       pw_env=pw_env, &
     827          822 :                       virial=virial)
     828              : 
     829          822 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     830              : 
     831          822 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
     832              : 
     833          822 :       offset = 0
     834          822 :       my_pos = v_rspace%pw_grid%para%group%mepos
     835          822 :       group_size = v_rspace%pw_grid%para%group%num_pe
     836              : 
     837         2448 :       DO ikind = 1, nkind
     838              : 
     839         1626 :          CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
     840         1626 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis_set, basis_type=basis_type)
     841              :          CALL get_gto_basis_set(gto_basis_set=lri_basis_set, &
     842              :                                 first_sgf=first_sgfa, &
     843              :                                 lmax=la_max, &
     844              :                                 lmin=la_min, &
     845              :                                 maxco=maxco, &
     846              :                                 maxsgf_set=maxsgf_set, &
     847              :                                 npgf=npgfa, &
     848              :                                 nset=nseta, &
     849              :                                 nsgf_set=nsgf_seta, &
     850              :                                 pgf_radius=rpgfa, &
     851              :                                 set_radius=set_radius_a, &
     852              :                                 sphi=sphi_a, &
     853         1626 :                                 zet=zeta)
     854              : 
     855         6504 :          ALLOCATE (hab(maxco, 1), pab(maxco, 1))
     856        43850 :          hab = 0._dp
     857        42224 :          pab(:, 1) = 0._dp
     858              : 
     859         4746 :          DO iatom = 1, natom_of_kind
     860              : 
     861         3120 :             atom_a = atom_list(iatom)
     862         3120 :             IF (PRESENT(atomlist)) THEN
     863          400 :                IF (atomlist(atom_a) == 0) CYCLE
     864              :             END IF
     865         2920 :             ra(:) = pbc(particle_set(atom_a)%r, cell)
     866         2920 :             force_a(:) = 0._dp
     867         2920 :             force_b(:) = 0._dp
     868         2920 :             my_virial_a(:, :) = 0._dp
     869         2920 :             my_virial_b(:, :) = 0._dp
     870              : 
     871        42868 :             m1 = MAXVAL(npgfa(1:nseta))
     872         8760 :             ALLOCATE (map_it(m1))
     873              : 
     874        42868 :             DO iset = 1, nseta
     875              :                !
     876        84320 :                map_it = .FALSE.
     877        84128 :                DO ipgf = 1, npgfa(iset)
     878        44180 :                   igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
     879        44180 :                   rs_grid => rs_v(igrid_level)
     880        84128 :                   map_it(ipgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
     881              :                END DO
     882        39948 :                offset = offset + 1
     883              :                !
     884        64958 :                IF (ANY(map_it(1:npgfa(iset)))) THEN
     885        19974 :                   sgfa = first_sgfa(1, iset)
     886        19974 :                   ncoa = npgfa(iset)*ncoset(la_max(iset))
     887       424846 :                   hab(:, 1) = 0._dp
     888        59922 :                   ALLOCATE (work_i(nsgf_seta(iset), 1))
     889       192924 :                   work_i = 0.0_dp
     890              : 
     891              :                   ! get fit coefficients for forces
     892        19974 :                   IF (calculate_forces) THEN
     893          895 :                      m1 = sgfa + nsgf_seta(iset) - 1
     894         1790 :                      ALLOCATE (work_f(nsgf_seta(iset), 1))
     895         6465 :                      work_f(1:nsgf_seta(iset), 1) = int_res(ikind)%acoef(iatom, sgfa:m1)
     896              :                      CALL dgemm("N", "N", ncoa, 1, nsgf_seta(iset), 1.0_dp, sphi_a(1, sgfa), &
     897              :                                 SIZE(sphi_a, 1), work_f(1, 1), SIZE(work_f, 1), 0.0_dp, pab(1, 1), &
     898          895 :                                 SIZE(pab, 1))
     899          895 :                      DEALLOCATE (work_f)
     900              :                   END IF
     901              : 
     902        42064 :                   DO ipgf = 1, npgfa(iset)
     903        22090 :                      na1 = (ipgf - 1)*ncoset(la_max(iset))
     904        22090 :                      igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
     905        22090 :                      rs_grid => rs_v(igrid_level)
     906              : 
     907              :                      radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
     908              :                                                        lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
     909              :                                                        zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
     910        22090 :                                                        prefactor=1.0_dp, cutoff=1.0_dp)
     911              : 
     912        42064 :                      IF (map_it(ipgf)) THEN
     913        22090 :                         IF (.NOT. calculate_forces) THEN
     914              :                            CALL integrate_pgf_product(la_max=la_max(iset), &
     915              :                                                       zeta=zeta(ipgf, iset), la_min=la_min(iset), &
     916              :                                                       lb_max=0, zetb=0.0_dp, lb_min=0, &
     917              :                                                       ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
     918              :                                                       rsgrid=rs_grid, &
     919              :                                                       hab=hab, o1=na1, o2=0, radius=radius, &
     920        20831 :                                                       calculate_forces=calculate_forces)
     921              :                         ELSE
     922              :                            CALL integrate_pgf_product(la_max=la_max(iset), &
     923              :                                                       zeta=zeta(ipgf, iset), la_min=la_min(iset), &
     924              :                                                       lb_max=0, zetb=0.0_dp, lb_min=0, &
     925              :                                                       ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
     926              :                                                       rsgrid=rs_grid, &
     927              :                                                       hab=hab, pab=pab, o1=na1, o2=0, radius=radius, &
     928              :                                                       calculate_forces=calculate_forces, &
     929              :                                                       force_a=force_a, force_b=force_b, &
     930              :                                                       use_virial=use_virial, &
     931         1259 :                                                       my_virial_a=my_virial_a, my_virial_b=my_virial_b)
     932              :                         END IF
     933              :                      END IF
     934              :                   END DO
     935              :                   ! contract hab
     936              :                   CALL dgemm("T", "N", nsgf_seta(iset), 1, ncoa, 1.0_dp, sphi_a(1, sgfa), &
     937        19974 :                              SIZE(sphi_a, 1), hab(1, 1), SIZE(hab, 1), 0.0_dp, work_i(1, 1), SIZE(work_i, 1))
     938              : 
     939              :                   int_res(ikind)%v_int(iatom, sgfa:sgfa - 1 + nsgf_seta(iset)) = &
     940       325926 :                      int_res(ikind)%v_int(iatom, sgfa:sgfa - 1 + nsgf_seta(iset)) + work_i(1:nsgf_seta(iset), 1)
     941        19974 :                   DEALLOCATE (work_i)
     942              :                END IF
     943              :             END DO
     944              :             !
     945         2920 :             IF (calculate_forces) THEN
     946          568 :                int_res(ikind)%v_dfdr(iatom, :) = int_res(ikind)%v_dfdr(iatom, :) + force_a(:)
     947          142 :                IF (use_virial) THEN
     948          156 :                   virial%pv_lrigpw = virial%pv_lrigpw + my_virial_a
     949          156 :                   virial%pv_virial = virial%pv_virial + my_virial_a
     950              :                END IF
     951              :             END IF
     952              : 
     953         4746 :             DEALLOCATE (map_it)
     954              : 
     955              :          END DO
     956              : 
     957         5700 :          DEALLOCATE (hab, pab)
     958              :       END DO
     959              : 
     960          822 :       CALL timestop(handle)
     961              : 
     962         1644 :    END SUBROUTINE integrate_v_rspace_one_center
     963              : 
     964              : ! **************************************************************************************************
     965              : !> \brief computes integrals of product of v_rspace times the diagonal block basis functions
     966              : !>        required for LRIGPW with exact 1c terms
     967              : !> \param v_rspace ...
     968              : !> \param ksmat ...
     969              : !> \param pmat ...
     970              : !> \param qs_env ...
     971              : !> \param calculate_forces ...
     972              : !> \param basis_type ...
     973              : !> \author JGH
     974              : ! **************************************************************************************************
     975           36 :    SUBROUTINE integrate_v_rspace_diagonal(v_rspace, ksmat, pmat, qs_env, calculate_forces, basis_type)
     976              :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: v_rspace
     977              :       TYPE(dbcsr_type), INTENT(INOUT)                    :: ksmat, pmat
     978              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     979              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     980              :       CHARACTER(len=*), INTENT(IN)                       :: basis_type
     981              : 
     982              :       CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_rspace_diagonal'
     983              : 
     984              :       INTEGER :: atom_a, group_size, handle, iatom, igrid_level, ikind, ipgf, iset, jpgf, jset, &
     985              :          m1, maxco, maxsgf_set, my_pos, na1, na2, natom_of_kind, nb1, nb2, ncoa, ncob, nkind, &
     986              :          nseta, nsgfa, offset, sgfa, sgfb
     987           36 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list, la_max, la_min, npgfa, &
     988           36 :                                                             nsgf_seta
     989           36 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa
     990              :       LOGICAL                                            :: found, use_virial
     991           36 :       LOGICAL, ALLOCATABLE, DIMENSION(:, :)              :: map_it2
     992              :       REAL(KIND=dp)                                      :: eps_rho_rspace, radius, zetp
     993              :       REAL(KIND=dp), DIMENSION(3)                        :: force_a, force_b, ra
     994              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: my_virial_a, my_virial_b
     995           36 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a
     996           36 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: h_block, hab, hmat, p_block, pab, pblk, &
     997           36 :                                                             rpgfa, sphi_a, work, zeta
     998           36 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     999              :       TYPE(cell_type), POINTER                           :: cell
    1000              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1001              :       TYPE(gridlevel_info_type), POINTER                 :: gridlevel_info
    1002              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
    1003              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1004           36 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1005              :       TYPE(pw_env_type), POINTER                         :: pw_env
    1006           36 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
    1007           36 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1008           36 :       TYPE(realspace_grid_type), DIMENSION(:), POINTER   :: rs_v
    1009              :       TYPE(realspace_grid_type), POINTER                 :: rs_grid
    1010              :       TYPE(virial_type), POINTER                         :: virial
    1011              : 
    1012           36 :       CALL timeset(routineN, handle)
    1013              : 
    1014           36 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
    1015           36 :       CALL pw_env_get(pw_env, rs_grids=rs_v)
    1016           36 :       CALL potential_pw2rs(rs_v, v_rspace, pw_env)
    1017              : 
    1018           36 :       gridlevel_info => pw_env%gridlevel_info
    1019              : 
    1020              :       CALL get_qs_env(qs_env=qs_env, &
    1021              :                       atomic_kind_set=atomic_kind_set, &
    1022              :                       qs_kind_set=qs_kind_set, &
    1023              :                       cell=cell, &
    1024              :                       dft_control=dft_control, &
    1025              :                       nkind=nkind, &
    1026              :                       particle_set=particle_set, &
    1027              :                       force=force, &
    1028              :                       virial=virial, &
    1029           36 :                       para_env=para_env)
    1030              : 
    1031           36 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
    1032           36 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
    1033              : 
    1034           36 :       offset = 0
    1035           36 :       my_pos = v_rspace%pw_grid%para%group%mepos
    1036           36 :       group_size = v_rspace%pw_grid%para%group%num_pe
    1037              : 
    1038          108 :       DO ikind = 1, nkind
    1039              : 
    1040           72 :          CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
    1041           72 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=basis_type)
    1042              :          CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
    1043              :                                 lmax=la_max, lmin=la_min, maxco=maxco, maxsgf_set=maxsgf_set, &
    1044              :                                 npgf=npgfa, nset=nseta, nsgf_set=nsgf_seta, nsgf=nsgfa, &
    1045              :                                 first_sgf=first_sgfa, pgf_radius=rpgfa, set_radius=set_radius_a, &
    1046           72 :                                 sphi=sphi_a, zet=zeta)
    1047              : 
    1048          720 :          ALLOCATE (hab(maxco, maxco), work(maxco, maxsgf_set), hmat(nsgfa, nsgfa))
    1049           72 :          IF (calculate_forces) ALLOCATE (pab(maxco, maxco), pblk(nsgfa, nsgfa))
    1050              : 
    1051          288 :          DO iatom = 1, natom_of_kind
    1052          216 :             atom_a = atom_list(iatom)
    1053          216 :             ra(:) = pbc(particle_set(atom_a)%r, cell)
    1054        17640 :             hmat = 0.0_dp
    1055          216 :             IF (calculate_forces) THEN
    1056            0 :                CALL dbcsr_get_block_p(matrix=pmat, row=atom_a, col=atom_a, BLOCK=p_block, found=found)
    1057            0 :                IF (found) THEN
    1058            0 :                   pblk(1:nsgfa, 1:nsgfa) = p_block(1:nsgfa, 1:nsgfa)
    1059              :                ELSE
    1060            0 :                   pblk = 0.0_dp
    1061              :                END IF
    1062            0 :                CALL para_env%sum(pblk)
    1063            0 :                force_a(:) = 0._dp
    1064            0 :                force_b(:) = 0._dp
    1065            0 :                IF (use_virial) THEN
    1066            0 :                   my_virial_a = 0.0_dp
    1067            0 :                   my_virial_b = 0.0_dp
    1068              :                END IF
    1069              :             END IF
    1070          432 :             m1 = MAXVAL(npgfa(1:nseta))
    1071          864 :             ALLOCATE (map_it2(m1, m1))
    1072          432 :             DO iset = 1, nseta
    1073          216 :                sgfa = first_sgfa(1, iset)
    1074          216 :                ncoa = npgfa(iset)*ncoset(la_max(iset))
    1075          648 :                DO jset = 1, nseta
    1076          216 :                   sgfb = first_sgfa(1, jset)
    1077          216 :                   ncob = npgfa(jset)*ncoset(la_max(jset))
    1078              :                   !
    1079        12312 :                   map_it2 = .FALSE.
    1080         1728 :                   DO ipgf = 1, npgfa(iset)
    1081        12312 :                      DO jpgf = 1, npgfa(jset)
    1082        10584 :                         zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
    1083        10584 :                         igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
    1084        10584 :                         rs_grid => rs_v(igrid_level)
    1085        12096 :                         map_it2(ipgf, jpgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
    1086              :                      END DO
    1087              :                   END DO
    1088          216 :                   offset = offset + 1
    1089              :                   !
    1090         6480 :                   IF (ANY(map_it2(1:npgfa(iset), 1:npgfa(jset)))) THEN
    1091       237492 :                      hab(:, :) = 0._dp
    1092          108 :                      IF (calculate_forces) THEN
    1093              :                         CALL dgemm("N", "N", ncoa, nsgf_seta(jset), nsgf_seta(iset), &
    1094              :                                    1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
    1095              :                                    pblk(sgfa, sgfb), SIZE(pblk, 1), &
    1096            0 :                                    0.0_dp, work(1, 1), SIZE(work, 1))
    1097              :                         CALL dgemm("N", "T", ncoa, ncob, nsgf_seta(jset), &
    1098              :                                    1.0_dp, work(1, 1), SIZE(work, 1), &
    1099              :                                    sphi_a(1, sgfb), SIZE(sphi_a, 1), &
    1100            0 :                                    0.0_dp, pab(1, 1), SIZE(pab, 1))
    1101              :                      END IF
    1102              : 
    1103          864 :                      DO ipgf = 1, npgfa(iset)
    1104          756 :                         na1 = (ipgf - 1)*ncoset(la_max(iset))
    1105          756 :                         na2 = ipgf*ncoset(la_max(iset))
    1106         6156 :                         DO jpgf = 1, npgfa(jset)
    1107         5292 :                            nb1 = (jpgf - 1)*ncoset(la_max(jset))
    1108         5292 :                            nb2 = jpgf*ncoset(la_max(jset))
    1109         5292 :                            zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
    1110         5292 :                            igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
    1111         5292 :                            rs_grid => rs_v(igrid_level)
    1112              : 
    1113              :                            radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
    1114              :                                                              lb_min=la_min(jset), lb_max=la_max(jset), &
    1115              :                                                              ra=ra, rb=ra, rp=ra, &
    1116              :                                                              zetp=zetp, eps=eps_rho_rspace, &
    1117         5292 :                                                              prefactor=1.0_dp, cutoff=1.0_dp)
    1118              : 
    1119         6048 :                            IF (map_it2(ipgf, jpgf)) THEN
    1120         5292 :                               IF (calculate_forces) THEN
    1121              :                                  CALL integrate_pgf_product( &
    1122              :                                     la_max(iset), zeta(ipgf, iset), la_min(iset), &
    1123              :                                     la_max(jset), zeta(jpgf, jset), la_min(jset), &
    1124              :                                     ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
    1125              :                                     rsgrid=rs_v(igrid_level), &
    1126              :                                     hab=hab, pab=pab, o1=na1, o2=nb1, &
    1127              :                                     radius=radius, &
    1128              :                                     calculate_forces=.TRUE., &
    1129              :                                     force_a=force_a, force_b=force_b, &
    1130            0 :                                     use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b)
    1131              :                               ELSE
    1132              :                                  CALL integrate_pgf_product( &
    1133              :                                     la_max(iset), zeta(ipgf, iset), la_min(iset), &
    1134              :                                     la_max(jset), zeta(jpgf, jset), la_min(jset), &
    1135              :                                     ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
    1136              :                                     rsgrid=rs_v(igrid_level), &
    1137              :                                     hab=hab, o1=na1, o2=nb1, &
    1138              :                                     radius=radius, &
    1139         5292 :                                     calculate_forces=.FALSE.)
    1140              :                               END IF
    1141              :                            END IF
    1142              :                         END DO
    1143              :                      END DO
    1144              :                      ! contract hab
    1145              :                      CALL dgemm("N", "N", ncoa, nsgf_seta(jset), ncob, &
    1146              :                                 1.0_dp, hab(1, 1), SIZE(hab, 1), &
    1147              :                                 sphi_a(1, sgfb), SIZE(sphi_a, 1), &
    1148          108 :                                 0.0_dp, work(1, 1), SIZE(work, 1))
    1149              :                      CALL dgemm("T", "N", nsgf_seta(iset), nsgf_seta(jset), ncoa, &
    1150              :                                 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
    1151              :                                 work(1, 1), SIZE(work, 1), &
    1152          108 :                                 1.0_dp, hmat(sgfa, sgfb), SIZE(hmat, 1))
    1153              :                   END IF
    1154              :                END DO
    1155              :             END DO
    1156          216 :             DEALLOCATE (map_it2)
    1157              :             ! update KS matrix
    1158        35064 :             CALL para_env%sum(hmat)
    1159          216 :             CALL dbcsr_get_block_p(matrix=ksmat, row=atom_a, col=atom_a, BLOCK=h_block, found=found)
    1160          216 :             IF (found) THEN
    1161        17532 :                h_block(1:nsgfa, 1:nsgfa) = h_block(1:nsgfa, 1:nsgfa) + hmat(1:nsgfa, 1:nsgfa)
    1162              :             END IF
    1163          504 :             IF (calculate_forces) THEN
    1164            0 :                force(ikind)%rho_elec(:, iatom) = force(ikind)%rho_elec(:, iatom) + 2.0_dp*force_a(:)
    1165            0 :                IF (use_virial) THEN
    1166              :                   IF (use_virial .AND. calculate_forces) THEN
    1167            0 :                      virial%pv_lrigpw = virial%pv_lrigpw + 2.0_dp*my_virial_a
    1168            0 :                      virial%pv_virial = virial%pv_virial + 2.0_dp*my_virial_a
    1169              :                   END IF
    1170              :                END IF
    1171              :             END IF
    1172              :          END DO
    1173           72 :          DEALLOCATE (hab, work, hmat)
    1174          252 :          IF (calculate_forces) DEALLOCATE (pab, pblk)
    1175              :       END DO
    1176              : 
    1177           36 :       CALL timestop(handle)
    1178              : 
    1179           72 :    END SUBROUTINE integrate_v_rspace_diagonal
    1180              : 
    1181              : ! **************************************************************************************************
    1182              : !> \brief computes integrals of product of v_rspace times a basis function (vector function)
    1183              : !>        and possible forces
    1184              : !> \param qs_env ...
    1185              : !> \param v_rspace ...
    1186              : !> \param f_coef ...
    1187              : !> \param f_integral ...
    1188              : !> \param calculate_forces ...
    1189              : !> \param basis_type ...
    1190              : !> \author JGH [8.2024]
    1191              : ! **************************************************************************************************
    1192            4 :    SUBROUTINE integrate_function(qs_env, v_rspace, f_coef, f_integral, calculate_forces, basis_type)
    1193              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1194              :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: v_rspace
    1195              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: f_coef
    1196              :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: f_integral
    1197              :       LOGICAL, INTENT(IN)                                :: calculate_forces
    1198              :       CHARACTER(len=*), INTENT(IN)                       :: basis_type
    1199              : 
    1200              :       CHARACTER(len=*), PARAMETER :: routineN = 'integrate_function'
    1201              : 
    1202              :       INTEGER :: atom_a, group_size, handle, i, iatom, igrid_level, ikind, ipgf, iset, maxco, &
    1203              :          maxsgf_set, my_pos, na1, natom, ncoa, nkind, nseta, offset, sgfa
    1204            4 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind
    1205            4 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, npgfa, nsgfa
    1206            4 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa
    1207              :       LOGICAL                                            :: use_virial
    1208              :       REAL(KIND=dp)                                      :: eps_rho_rspace, radius
    1209              :       REAL(KIND=dp), DIMENSION(3)                        :: force_a, force_b, ra
    1210              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: my_virial_a, my_virial_b
    1211            4 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: hab, pab, sphi_a, work, zeta
    1212            4 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1213              :       TYPE(cell_type), POINTER                           :: cell
    1214              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1215              :       TYPE(gridlevel_info_type), POINTER                 :: gridlevel_info
    1216              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
    1217            4 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1218              :       TYPE(pw_env_type), POINTER                         :: pw_env
    1219            4 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
    1220            4 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1221            4 :       TYPE(realspace_grid_type), DIMENSION(:), POINTER   :: rs_v
    1222              :       TYPE(realspace_grid_type), POINTER                 :: rs_grid
    1223              :       TYPE(virial_type), POINTER                         :: virial
    1224              : 
    1225            4 :       CALL timeset(routineN, handle)
    1226              : 
    1227            4 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
    1228            4 :       gridlevel_info => pw_env%gridlevel_info
    1229              : 
    1230            4 :       CALL pw_env_get(pw_env, rs_grids=rs_v)
    1231            4 :       CALL potential_pw2rs(rs_v, v_rspace, pw_env)
    1232              : 
    1233              :       CALL get_qs_env(qs_env=qs_env, &
    1234              :                       atomic_kind_set=atomic_kind_set, &
    1235              :                       qs_kind_set=qs_kind_set, &
    1236              :                       force=force, &
    1237              :                       cell=cell, &
    1238              :                       dft_control=dft_control, &
    1239              :                       nkind=nkind, &
    1240              :                       natom=natom, &
    1241              :                       particle_set=particle_set, &
    1242              :                       pw_env=pw_env, &
    1243            4 :                       virial=virial)
    1244            4 :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
    1245              : 
    1246            4 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
    1247            4 :       IF (use_virial) THEN
    1248            0 :          CPABORT("Virial NYA")
    1249              :       END IF
    1250              : 
    1251            4 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
    1252              : 
    1253              :       CALL get_qs_kind_set(qs_kind_set, &
    1254            4 :                            maxco=maxco, maxsgf_set=maxsgf_set, basis_type=basis_type)
    1255           20 :       ALLOCATE (hab(maxco, 1), pab(maxco, 1), work(maxco, 1))
    1256              : 
    1257            4 :       offset = 0
    1258            4 :       my_pos = v_rspace%pw_grid%para%group%mepos
    1259            4 :       group_size = v_rspace%pw_grid%para%group%num_pe
    1260              : 
    1261           20 :       DO iatom = 1, natom
    1262           16 :          ikind = particle_set(iatom)%atomic_kind%kind_number
    1263           16 :          atom_a = atom_of_kind(iatom)
    1264           16 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=basis_type)
    1265              :          CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
    1266              :                                 first_sgf=first_sgfa, &
    1267              :                                 lmax=la_max, &
    1268              :                                 lmin=la_min, &
    1269              :                                 npgf=npgfa, &
    1270              :                                 nset=nseta, &
    1271              :                                 nsgf_set=nsgfa, &
    1272              :                                 sphi=sphi_a, &
    1273           16 :                                 zet=zeta)
    1274           16 :          ra(:) = pbc(particle_set(iatom)%r, cell)
    1275              : 
    1276           16 :          force_a(:) = 0._dp
    1277           16 :          force_b(:) = 0._dp
    1278           16 :          my_virial_a(:, :) = 0._dp
    1279           16 :          my_virial_b(:, :) = 0._dp
    1280              : 
    1281           32 :          DO iset = 1, nseta
    1282              : 
    1283           16 :             ncoa = npgfa(iset)*ncoset(la_max(iset))
    1284           16 :             sgfa = first_sgfa(1, iset)
    1285              : 
    1286          144 :             hab = 0._dp
    1287          144 :             pab = 0._dp
    1288              : 
    1289          120 :             DO i = 1, nsgfa(iset)
    1290          120 :                work(i, 1) = f_coef(offset + i)
    1291              :             END DO
    1292              : 
    1293              :             CALL dgemm("N", "N", ncoa, 1, nsgfa(iset), &
    1294              :                        1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
    1295              :                        work(1, 1), SIZE(work, 1), &
    1296           16 :                        0.0_dp, pab(1, 1), SIZE(pab, 1))
    1297              : 
    1298          120 :             DO ipgf = 1, npgfa(iset)
    1299              : 
    1300          104 :                na1 = (ipgf - 1)*ncoset(la_max(iset))
    1301              : 
    1302          104 :                igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
    1303          104 :                rs_grid => rs_v(igrid_level)
    1304              : 
    1305          120 :                IF (map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)) THEN
    1306              :                   radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
    1307              :                                                     lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
    1308              :                                                     zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
    1309           52 :                                                     prefactor=1.0_dp, cutoff=1.0_dp)
    1310              : 
    1311           52 :                   IF (calculate_forces) THEN
    1312              :                      CALL integrate_pgf_product(la_max=la_max(iset), &
    1313              :                                                 zeta=zeta(ipgf, iset), la_min=la_min(iset), &
    1314              :                                                 lb_max=0, zetb=0.0_dp, lb_min=0, &
    1315              :                                                 ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
    1316              :                                                 rsgrid=rs_grid, &
    1317              :                                                 hab=hab, pab=pab, o1=na1, o2=0, radius=radius, &
    1318              :                                                 calculate_forces=.TRUE., &
    1319              :                                                 force_a=force_a, force_b=force_b, &
    1320              :                                                 use_virial=use_virial, &
    1321           52 :                                                 my_virial_a=my_virial_a, my_virial_b=my_virial_b)
    1322              :                   ELSE
    1323              :                      CALL integrate_pgf_product(la_max=la_max(iset), &
    1324              :                                                 zeta=zeta(ipgf, iset), la_min=la_min(iset), &
    1325              :                                                 lb_max=0, zetb=0.0_dp, lb_min=0, &
    1326              :                                                 ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
    1327              :                                                 rsgrid=rs_grid, &
    1328              :                                                 hab=hab, o1=na1, o2=0, radius=radius, &
    1329            0 :                                                 calculate_forces=.FALSE.)
    1330              :                   END IF
    1331              : 
    1332              :                END IF
    1333              : 
    1334              :             END DO
    1335              :             !
    1336              :             CALL dgemm("T", "N", nsgfa(iset), 1, ncoa, 1.0_dp, sphi_a(1, sgfa), &
    1337           16 :                        SIZE(sphi_a, 1), hab(1, 1), SIZE(hab, 1), 0.0_dp, work(1, 1), SIZE(work, 1))
    1338          120 :             DO i = 1, nsgfa(iset)
    1339          120 :                f_integral(offset + i) = work(i, 1)
    1340              :             END DO
    1341              : 
    1342           32 :             offset = offset + nsgfa(iset)
    1343              : 
    1344              :          END DO
    1345              : 
    1346           36 :          IF (calculate_forces) THEN
    1347           64 :             force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + force_a(:)
    1348           16 :             IF (use_virial) THEN
    1349            0 :                virial%pv_virial = virial%pv_virial + my_virial_a
    1350              :             END IF
    1351              :          END IF
    1352              : 
    1353              :       END DO
    1354              : 
    1355            4 :       DEALLOCATE (hab, pab, work)
    1356              : 
    1357            4 :       CALL timestop(handle)
    1358              : 
    1359           12 :    END SUBROUTINE integrate_function
    1360              : 
    1361              : END MODULE qs_integrate_potential_single
        

Generated by: LCOV version 2.0-1