LCOV - code coverage report
Current view: top level - src - kpoint_coulomb_2c.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 99.3 % 299 297
Test Date: 2025-07-25 12:55:17 Functions: 87.5 % 16 14

            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 Routines to compute the Coulomb integral V_(alpha beta)(k) for a k-point k using lattice
      10              : !>        summation in real space. These integrals are e.g. needed in periodic RI for RPA, GW
      11              : !> \par History
      12              : !>       2018.05 created [Jan Wilhelm]
      13              : !> \author Jan Wilhelm
      14              : ! **************************************************************************************************
      15              : MODULE kpoint_coulomb_2c
      16              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      17              :                                               get_atomic_kind_set
      18              :    USE basis_set_types,                 ONLY: gto_basis_set_type
      19              :    USE cell_types,                      ONLY: cell_type,&
      20              :                                               get_cell,&
      21              :                                               pbc
      22              :    USE constants_operator,              ONLY: operator_coulomb
      23              :    USE cp_dbcsr_api,                    ONLY: &
      24              :         dbcsr_create, dbcsr_init_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
      25              :         dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
      26              :         dbcsr_release_p, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
      27              :    USE cp_dbcsr_contrib,                ONLY: dbcsr_reserve_all_blocks
      28              :    USE generic_shg_integrals,           ONLY: int_operators_r12_ab_shg
      29              :    USE generic_shg_integrals_init,      ONLY: contraction_matrix_shg
      30              :    USE kinds,                           ONLY: dp
      31              :    USE kpoint_types,                    ONLY: get_kpoint_info,&
      32              :                                               kpoint_type
      33              :    USE mathconstants,                   ONLY: gaussi,&
      34              :                                               twopi,&
      35              :                                               z_one
      36              :    USE message_passing,                 ONLY: mp_para_env_type
      37              :    USE particle_types,                  ONLY: particle_type
      38              :    USE qs_environment_types,            ONLY: get_qs_env,&
      39              :                                               qs_environment_type
      40              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      41              :                                               qs_kind_type
      42              : #include "./base/base_uses.f90"
      43              : 
      44              :    IMPLICIT NONE
      45              : 
      46              :    PRIVATE
      47              : 
      48              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kpoint_coulomb_2c'
      49              : 
      50              :    PUBLIC :: build_2c_coulomb_matrix_kp, build_2c_coulomb_matrix_kp_small_cell
      51              : 
      52              : ! **************************************************************************************************
      53              : 
      54              :    TYPE two_d_util_type
      55              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)  :: block
      56              :    END TYPE
      57              : 
      58              : CONTAINS
      59              : 
      60              : ! **************************************************************************************************
      61              : !> \brief ...
      62              : !> \param matrix_v_kp ...
      63              : !> \param kpoints ...
      64              : !> \param basis_type ...
      65              : !> \param cell ...
      66              : !> \param particle_set ...
      67              : !> \param qs_kind_set ...
      68              : !> \param atomic_kind_set ...
      69              : !> \param size_lattice_sum ...
      70              : !> \param operator_type ...
      71              : !> \param ikp_start ...
      72              : !> \param ikp_end ...
      73              : ! **************************************************************************************************
      74          100 :    SUBROUTINE build_2c_coulomb_matrix_kp(matrix_v_kp, kpoints, basis_type, cell, particle_set, qs_kind_set, &
      75              :                                          atomic_kind_set, size_lattice_sum, operator_type, ikp_start, ikp_end)
      76              : 
      77              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_v_kp
      78              :       TYPE(kpoint_type), POINTER                         :: kpoints
      79              :       CHARACTER(LEN=*), INTENT(IN)                       :: basis_type
      80              :       TYPE(cell_type), POINTER                           :: cell
      81              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      82              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      83              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      84              :       INTEGER                                            :: size_lattice_sum, operator_type, &
      85              :                                                             ikp_start, ikp_end
      86              : 
      87              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'build_2c_coulomb_matrix_kp'
      88              : 
      89              :       INTEGER                                            :: handle
      90              :       TYPE(dbcsr_type), POINTER                          :: matrix_v_L_tmp
      91              : 
      92          100 :       CALL timeset(routineN, handle)
      93              : 
      94          100 :       CALL allocate_tmp(matrix_v_L_tmp, matrix_v_kp, ikp_start)
      95              : 
      96              :       CALL lattice_sum(matrix_v_kp, kpoints, basis_type, cell, particle_set, &
      97              :                        qs_kind_set, atomic_kind_set, size_lattice_sum, matrix_v_L_tmp, &
      98          100 :                        operator_type, ikp_start, ikp_end)
      99              : 
     100          100 :       CALL deallocate_tmp(matrix_v_L_tmp)
     101              : 
     102          100 :       CALL timestop(handle)
     103              : 
     104          100 :    END SUBROUTINE build_2c_coulomb_matrix_kp
     105              : 
     106              : ! **************************************************************************************************
     107              : !> \brief ...
     108              : !> \param matrix_v_kp ...
     109              : !> \param kpoints ...
     110              : !> \param basis_type ...
     111              : !> \param cell ...
     112              : !> \param particle_set ...
     113              : !> \param qs_kind_set ...
     114              : !> \param atomic_kind_set ...
     115              : !> \param size_lattice_sum ...
     116              : !> \param matrix_v_L_tmp ...
     117              : !> \param operator_type ...
     118              : !> \param ikp_start ...
     119              : !> \param ikp_end ...
     120              : ! **************************************************************************************************
     121          100 :    SUBROUTINE lattice_sum(matrix_v_kp, kpoints, basis_type, cell, particle_set, &
     122              :                           qs_kind_set, atomic_kind_set, size_lattice_sum, matrix_v_L_tmp, &
     123              :                           operator_type, ikp_start, ikp_end)
     124              : 
     125              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_v_kp
     126              :       TYPE(kpoint_type), POINTER                         :: kpoints
     127              :       CHARACTER(LEN=*), INTENT(IN)                       :: basis_type
     128              :       TYPE(cell_type), POINTER                           :: cell
     129              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     130              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     131              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     132              :       INTEGER                                            :: size_lattice_sum
     133              :       TYPE(dbcsr_type), POINTER                          :: matrix_v_L_tmp
     134              :       INTEGER                                            :: operator_type, ikp_start, ikp_end
     135              : 
     136              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'lattice_sum'
     137              : 
     138              :       INTEGER :: factor, handle, handle2, i_block, i_x, i_x_inner, i_x_outer, ik, j_y, j_y_inner, &
     139              :          j_y_outer, k_z, k_z_inner, k_z_outer, x_max, x_min, y_max, y_min, z_max, z_min
     140              :       INTEGER, DIMENSION(3)                              :: nkp_grid
     141              :       REAL(KIND=dp)                                      :: coskl, sinkl
     142              :       REAL(KIND=dp), DIMENSION(3)                        :: vec_L, vec_s
     143              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: hmat
     144          100 :       TYPE(two_d_util_type), ALLOCATABLE, DIMENSION(:)   :: blocks_v_L, blocks_v_L_store
     145              :       TYPE(two_d_util_type), ALLOCATABLE, &
     146          100 :          DIMENSION(:, :, :)                              :: blocks_v_kp
     147              : 
     148          100 :       CALL timeset(routineN, handle)
     149              : 
     150              :       CALL get_factor_and_xyz_min_max(cell, kpoints, size_lattice_sum, factor, hmat, &
     151          100 :                                       x_min, x_max, y_min, y_max, z_min, z_max, nkp_grid)
     152              : 
     153          100 :       CALL allocate_blocks_v_kp(blocks_v_kp, matrix_v_kp, ikp_start, ikp_end)
     154          100 :       CALL allocate_blocks_v_L(blocks_v_L, matrix_v_L_tmp)
     155          100 :       CALL allocate_blocks_v_L(blocks_v_L_store, matrix_v_L_tmp)
     156              : 
     157          488 :       DO i_x_inner = 0, 2*nkp_grid(1) - 1
     158         3784 :          DO j_y_inner = 0, 2*nkp_grid(2) - 1
     159        24100 :             DO k_z_inner = 0, 2*nkp_grid(3) - 1
     160              : 
     161        50944 :                DO i_x_outer = x_min, x_max + nkp_grid(1), 2*nkp_grid(1)
     162       160000 :                   DO j_y_outer = y_min, y_max + nkp_grid(2), 2*nkp_grid(2)
     163       547520 :                      DO k_z_outer = z_min, z_max + nkp_grid(3), 2*nkp_grid(3)
     164              : 
     165       407936 :                         i_x = i_x_inner + i_x_outer
     166       407936 :                         j_y = j_y_inner + j_y_outer
     167       407936 :                         k_z = k_z_inner + k_z_outer
     168              : 
     169              :                         IF (i_x > x_max .OR. i_x < x_min .OR. &
     170              :                             j_y > y_max .OR. j_y < y_min .OR. &
     171       407936 :                             k_z > z_max .OR. k_z < z_min) CYCLE
     172              : 
     173       625408 :                         vec_s = [REAL(i_x, dp), REAL(j_y, dp), REAL(k_z, dp)]
     174              : 
     175      2032576 :                         vec_L = MATMUL(hmat, vec_s)
     176              : 
     177              :                         ! Compute (P 0 | Q vec_L) and store it in matrix_v_L_tmp
     178              :                         CALL compute_v_transl(matrix_v_L_tmp, blocks_v_L, vec_L, particle_set, &
     179              :                                               qs_kind_set, atomic_kind_set, basis_type, cell, &
     180       156352 :                                               operator_type)
     181              : 
     182       845092 :                         DO i_block = 1, SIZE(blocks_v_L)
     183              :                            blocks_v_L_store(i_block)%block(:, :) = blocks_v_L_store(i_block)%block(:, :) &
     184    139785176 :                                                                    + blocks_v_L(i_block)%block(:, :)
     185              :                         END DO
     186              : 
     187              :                      END DO
     188              :                   END DO
     189              :                END DO
     190              : 
     191        20416 :                CALL timeset(routineN//"_R_to_k", handle2)
     192              : 
     193              :                ! add exp(iq*vec_L) * (P 0 | Q vec_L) to V_PQ(q)
     194       169728 :                DO ik = ikp_start, ikp_end
     195              : 
     196              :                   ! coskl and sinkl are identical for all i_x_outer, j_y_outer, k_z_outer
     197       597248 :                   coskl = COS(twopi*DOT_PRODUCT(vec_s(1:3), kpoints%xkp(1:3, ik)))
     198       597248 :                   sinkl = SIN(twopi*DOT_PRODUCT(vec_s(1:3), kpoints%xkp(1:3, ik)))
     199              : 
     200       748352 :                   DO i_block = 1, SIZE(blocks_v_L)
     201              : 
     202              :                      blocks_v_kp(ik, 1, i_block)%block(:, :) = blocks_v_kp(ik, 1, i_block)%block(:, :) &
     203    282880576 :                                                                + coskl*blocks_v_L_store(i_block)%block(:, :)
     204              :                      blocks_v_kp(ik, 2, i_block)%block(:, :) = blocks_v_kp(ik, 2, i_block)%block(:, :) &
     205    283029888 :                                                                + sinkl*blocks_v_L_store(i_block)%block(:, :)
     206              : 
     207              :                   END DO
     208              : 
     209              :                END DO
     210              : 
     211        93088 :                DO i_block = 1, SIZE(blocks_v_L)
     212              : 
     213     19736480 :                   blocks_v_L_store(i_block)%block(:, :) = 0.0_dp
     214              : 
     215              :                END DO
     216              : 
     217        44128 :                CALL timestop(handle2)
     218              : 
     219              :             END DO
     220              :          END DO
     221              :       END DO
     222              : 
     223          100 :       CALL set_blocks_to_matrix_v_kp(matrix_v_kp, blocks_v_kp, ikp_start, ikp_end)
     224              : 
     225          100 :       CALL deallocate_blocks_v_kp(blocks_v_kp)
     226          100 :       CALL deallocate_blocks_v_L(blocks_v_L)
     227          100 :       CALL deallocate_blocks_v_L(blocks_v_L_store)
     228              : 
     229          100 :       CALL timestop(handle)
     230              : 
     231          100 :    END SUBROUTINE lattice_sum
     232              : 
     233              : ! **************************************************************************************************
     234              : !> \brief ...
     235              : !> \param matrix_v_kp ...
     236              : !> \param blocks_v_kp ...
     237              : !> \param ikp_start ...
     238              : !> \param ikp_end ...
     239              : ! **************************************************************************************************
     240          100 :    SUBROUTINE set_blocks_to_matrix_v_kp(matrix_v_kp, blocks_v_kp, ikp_start, ikp_end)
     241              : 
     242              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_v_kp
     243              :       TYPE(two_d_util_type), ALLOCATABLE, &
     244              :          DIMENSION(:, :, :)                              :: blocks_v_kp
     245              :       INTEGER                                            :: ikp_start, ikp_end
     246              : 
     247              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'set_blocks_to_matrix_v_kp'
     248              : 
     249              :       INTEGER                                            :: col, handle, i_block, i_real_im, ik, row
     250          100 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: data_block
     251              :       TYPE(dbcsr_iterator_type)                          :: iter
     252              : 
     253          100 :       CALL timeset(routineN, handle)
     254              : 
     255          740 :       DO ik = ikp_start, ikp_end
     256              : 
     257         2020 :          DO i_real_im = 1, 2
     258              : 
     259         1280 :             i_block = 1
     260              : 
     261         1280 :             CALL dbcsr_iterator_start(iter, matrix_v_kp(ik, i_real_im)%matrix)
     262              : 
     263         6220 :             DO WHILE (dbcsr_iterator_blocks_left(iter))
     264              : 
     265         4940 :                CALL dbcsr_iterator_next_block(iter, row, col, data_block)
     266              : 
     267      2423116 :                data_block(:, :) = blocks_v_kp(ik, i_real_im, i_block)%block(:, :)
     268              : 
     269         4940 :                i_block = i_block + 1
     270              : 
     271              :             END DO
     272              : 
     273         3200 :             CALL dbcsr_iterator_stop(iter)
     274              : 
     275              :          END DO
     276              : 
     277              :       END DO
     278              : 
     279          100 :       CALL timestop(handle)
     280              : 
     281          100 :    END SUBROUTINE set_blocks_to_matrix_v_kp
     282              : 
     283              : ! **************************************************************************************************
     284              : !> \brief ...
     285              : !> \param matrix_v_L_tmp ...
     286              : !> \param blocks_v_L ...
     287              : !> \param vec_L ...
     288              : !> \param particle_set ...
     289              : !> \param qs_kind_set ...
     290              : !> \param atomic_kind_set ...
     291              : !> \param basis_type ...
     292              : !> \param cell ...
     293              : !> \param operator_type ...
     294              : ! **************************************************************************************************
     295       156352 :    SUBROUTINE compute_v_transl(matrix_v_L_tmp, blocks_v_L, vec_L, particle_set, &
     296              :                                qs_kind_set, atomic_kind_set, basis_type, cell, operator_type)
     297              : 
     298              :       TYPE(dbcsr_type), POINTER                          :: matrix_v_L_tmp
     299              :       TYPE(two_d_util_type), ALLOCATABLE, DIMENSION(:)   :: blocks_v_L
     300              :       REAL(KIND=dp), DIMENSION(3)                        :: vec_L
     301              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     302              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     303              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     304              :       CHARACTER(LEN=*), INTENT(IN)                       :: basis_type
     305              :       TYPE(cell_type), POINTER                           :: cell
     306              :       INTEGER                                            :: operator_type
     307              : 
     308              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'compute_v_transl'
     309              : 
     310              :       INTEGER                                            :: col, handle, i_block, kind_a, kind_b, row
     311       156352 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: kind_of
     312              :       REAL(dp), DIMENSION(3)                             :: ra, rab_L, rb
     313       156352 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: data_block
     314       156352 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: contr_a, contr_b
     315              :       TYPE(dbcsr_iterator_type)                          :: iter
     316              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
     317              : 
     318       156352 :       CALL timeset(routineN, handle)
     319              : 
     320       156352 :       NULLIFY (basis_set_a, basis_set_b, data_block)
     321              : 
     322       156352 :       CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
     323              : 
     324       156352 :       CALL dbcsr_set(matrix_v_L_tmp, 0.0_dp)
     325              : 
     326       156352 :       i_block = 1
     327              : 
     328       156352 :       CALL dbcsr_iterator_start(iter, matrix_v_L_tmp)
     329              : 
     330       736036 :       DO WHILE (dbcsr_iterator_blocks_left(iter))
     331              : 
     332       579684 :          CALL dbcsr_iterator_next_block(iter, row, col, data_block)
     333              : 
     334       579684 :          kind_a = kind_of(row)
     335       579684 :          kind_b = kind_of(col)
     336              : 
     337       579684 :          CALL get_qs_kind(qs_kind=qs_kind_set(kind_a), basis_set=basis_set_a, basis_type=basis_type)
     338       579684 :          CALL get_qs_kind(qs_kind=qs_kind_set(kind_b), basis_set=basis_set_b, basis_type=basis_type)
     339              : 
     340       579684 :          ra(1:3) = pbc(particle_set(row)%r(1:3), cell)
     341       579684 :          rb(1:3) = pbc(particle_set(col)%r(1:3), cell)
     342              : 
     343      2318736 :          rab_L(1:3) = rb(1:3) - ra(1:3) + vec_L(1:3)
     344              : 
     345       579684 :          CALL contraction_matrix_shg(basis_set_a, contr_a)
     346       579684 :          CALL contraction_matrix_shg(basis_set_b, contr_b)
     347              : 
     348    139377240 :          blocks_v_L(i_block)%block = 0.0_dp
     349              : 
     350              :          CALL int_operators_r12_ab_shg(operator_type, blocks_v_L(i_block)%block, rab=rab_L, &
     351              :                                        fba=basis_set_a, fbb=basis_set_b, scona_shg=contr_a, sconb_shg=contr_b, &
     352       579684 :                                        calculate_forces=.FALSE.)
     353              : 
     354       579684 :          i_block = i_block + 1
     355              : 
     356       579684 :          DEALLOCATE (contr_a, contr_b)
     357              : 
     358              :       END DO
     359              : 
     360       156352 :       CALL dbcsr_iterator_stop(iter)
     361              : 
     362       156352 :       DEALLOCATE (kind_of)
     363              : 
     364       156352 :       CALL timestop(handle)
     365              : 
     366       156352 :    END SUBROUTINE compute_v_transl
     367              : 
     368              : ! **************************************************************************************************
     369              : !> \brief ...
     370              : !> \param blocks_v_kp ...
     371              : ! **************************************************************************************************
     372          100 :    SUBROUTINE deallocate_blocks_v_kp(blocks_v_kp)
     373              :       TYPE(two_d_util_type), ALLOCATABLE, &
     374              :          DIMENSION(:, :, :)                              :: blocks_v_kp
     375              : 
     376              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'deallocate_blocks_v_kp'
     377              : 
     378              :       INTEGER                                            :: handle, i_block, i_real_img, ik
     379              : 
     380          100 :       CALL timeset(routineN, handle)
     381              : 
     382          940 :       DO ik = LBOUND(blocks_v_kp, 1), UBOUND(blocks_v_kp, 1)
     383         2020 :          DO i_real_img = 1, SIZE(blocks_v_kp, 2)
     384         6860 :             DO i_block = 1, SIZE(blocks_v_kp, 3)
     385         6220 :                DEALLOCATE (blocks_v_kp(ik, i_real_img, i_block)%block)
     386              :             END DO
     387              :          END DO
     388              :       END DO
     389              : 
     390         5040 :       DEALLOCATE (blocks_v_kp)
     391              : 
     392          100 :       CALL timestop(handle)
     393              : 
     394          100 :    END SUBROUTINE
     395              : 
     396              : ! **************************************************************************************************
     397              : !> \brief ...
     398              : !> \param blocks_v_L ...
     399              : ! **************************************************************************************************
     400          200 :    SUBROUTINE deallocate_blocks_v_L(blocks_v_L)
     401              :       TYPE(two_d_util_type), ALLOCATABLE, DIMENSION(:)   :: blocks_v_L
     402              : 
     403              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'deallocate_blocks_v_L'
     404              : 
     405              :       INTEGER                                            :: handle, i_block
     406              : 
     407          200 :       CALL timeset(routineN, handle)
     408              : 
     409          900 :       DO i_block = 1, SIZE(blocks_v_L, 1)
     410          900 :          DEALLOCATE (blocks_v_L(i_block)%block)
     411              :       END DO
     412              : 
     413          900 :       DEALLOCATE (blocks_v_L)
     414              : 
     415          200 :       CALL timestop(handle)
     416              : 
     417          200 :    END SUBROUTINE
     418              : 
     419              : ! **************************************************************************************************
     420              : !> \brief ...
     421              : !> \param blocks_v_L ...
     422              : !> \param matrix_v_L_tmp ...
     423              : ! **************************************************************************************************
     424          200 :    SUBROUTINE allocate_blocks_v_L(blocks_v_L, matrix_v_L_tmp)
     425              :       TYPE(two_d_util_type), ALLOCATABLE, DIMENSION(:)   :: blocks_v_L
     426              :       TYPE(dbcsr_type), POINTER                          :: matrix_v_L_tmp
     427              : 
     428              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_blocks_v_L'
     429              : 
     430              :       INTEGER                                            :: col, handle, i_block, nblocks, row
     431          200 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: data_block
     432              :       TYPE(dbcsr_iterator_type)                          :: iter
     433              : 
     434          200 :       CALL timeset(routineN, handle)
     435              : 
     436          200 :       nblocks = 0
     437              : 
     438          200 :       CALL dbcsr_iterator_start(iter, matrix_v_L_tmp)
     439              : 
     440          900 :       DO WHILE (dbcsr_iterator_blocks_left(iter))
     441              : 
     442          700 :          CALL dbcsr_iterator_next_block(iter, row, col, data_block)
     443              : 
     444          700 :          nblocks = nblocks + 1
     445              : 
     446              :       END DO
     447              : 
     448          200 :       CALL dbcsr_iterator_stop(iter)
     449              : 
     450         1300 :       ALLOCATE (blocks_v_L(nblocks))
     451              : 
     452          200 :       i_block = 1
     453              : 
     454          200 :       CALL dbcsr_iterator_start(iter, matrix_v_L_tmp)
     455              : 
     456          900 :       DO WHILE (dbcsr_iterator_blocks_left(iter))
     457              : 
     458          700 :          CALL dbcsr_iterator_next_block(iter, row, col, data_block)
     459              : 
     460         2800 :          ALLOCATE (blocks_v_L(i_block)%block(SIZE(data_block, 1), SIZE(data_block, 2)))
     461       217232 :          blocks_v_L(i_block)%block = 0.0_dp
     462              : 
     463         1600 :          i_block = i_block + 1
     464              : 
     465              :       END DO
     466              : 
     467          200 :       CALL dbcsr_iterator_stop(iter)
     468              : 
     469          200 :       CALL timestop(handle)
     470              : 
     471          400 :    END SUBROUTINE
     472              : 
     473              : ! **************************************************************************************************
     474              : !> \brief ...
     475              : !> \param blocks_v_kp ...
     476              : !> \param matrix_v_kp ...
     477              : !> \param ikp_start ...
     478              : !> \param ikp_end ...
     479              : ! **************************************************************************************************
     480          100 :    SUBROUTINE allocate_blocks_v_kp(blocks_v_kp, matrix_v_kp, ikp_start, ikp_end)
     481              :       TYPE(two_d_util_type), ALLOCATABLE, &
     482              :          DIMENSION(:, :, :)                              :: blocks_v_kp
     483              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_v_kp
     484              :       INTEGER                                            :: ikp_start, ikp_end
     485              : 
     486              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_blocks_v_kp'
     487              : 
     488              :       INTEGER                                            :: col, handle, i_block, i_real_img, ik, &
     489              :                                                             nblocks, row
     490          100 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: data_block
     491              :       TYPE(dbcsr_iterator_type)                          :: iter
     492              : 
     493          100 :       CALL timeset(routineN, handle)
     494              : 
     495          100 :       nblocks = 0
     496              : 
     497          100 :       CALL dbcsr_iterator_start(iter, matrix_v_kp(ikp_start, 1)%matrix)
     498              : 
     499          450 :       DO WHILE (dbcsr_iterator_blocks_left(iter))
     500              : 
     501          350 :          CALL dbcsr_iterator_next_block(iter, row, col, data_block)
     502              : 
     503          350 :          nblocks = nblocks + 1
     504              : 
     505              :       END DO
     506              : 
     507          100 :       CALL dbcsr_iterator_stop(iter)
     508              : 
     509         6490 :       ALLOCATE (blocks_v_kp(ikp_start:ikp_end, SIZE(matrix_v_kp, 2), nblocks))
     510              : 
     511          740 :       DO ik = ikp_start, ikp_end
     512              : 
     513         2020 :          DO i_real_img = 1, SIZE(matrix_v_kp, 2)
     514              : 
     515         1280 :             i_block = 1
     516              : 
     517         1280 :             CALL dbcsr_iterator_start(iter, matrix_v_kp(ik, i_real_img)%matrix)
     518              : 
     519         6220 :             DO WHILE (dbcsr_iterator_blocks_left(iter))
     520              : 
     521         4940 :                CALL dbcsr_iterator_next_block(iter, row, col, data_block)
     522              : 
     523            0 :                ALLOCATE (blocks_v_kp(ik, i_real_img, i_block)%block(SIZE(data_block, 1), &
     524        19760 :                                                                     SIZE(data_block, 2)))
     525      2423116 :                blocks_v_kp(ik, i_real_img, i_block)%block = 0.0_dp
     526              : 
     527        11160 :                i_block = i_block + 1
     528              : 
     529              :             END DO
     530              : 
     531         3200 :             CALL dbcsr_iterator_stop(iter)
     532              : 
     533              :          END DO
     534              : 
     535              :       END DO
     536              : 
     537          100 :       CALL timestop(handle)
     538              : 
     539          100 :    END SUBROUTINE
     540              : 
     541              : ! **************************************************************************************************
     542              : !> \brief ...
     543              : !> \param cell ...
     544              : !> \param kpoints ...
     545              : !> \param size_lattice_sum ...
     546              : !> \param factor ...
     547              : !> \param hmat ...
     548              : !> \param x_min ...
     549              : !> \param x_max ...
     550              : !> \param y_min ...
     551              : !> \param y_max ...
     552              : !> \param z_min ...
     553              : !> \param z_max ...
     554              : !> \param nkp_grid ...
     555              : ! **************************************************************************************************
     556          112 :    SUBROUTINE get_factor_and_xyz_min_max(cell, kpoints, size_lattice_sum, factor, hmat, &
     557              :                                          x_min, x_max, y_min, y_max, z_min, z_max, nkp_grid)
     558              : 
     559              :       TYPE(cell_type), POINTER                           :: cell
     560              :       TYPE(kpoint_type), POINTER                         :: kpoints
     561              :       INTEGER                                            :: size_lattice_sum, factor
     562              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: hmat
     563              :       INTEGER                                            :: x_min, x_max, y_min, y_max, z_min, z_max
     564              :       INTEGER, DIMENSION(3)                              :: nkp_grid
     565              : 
     566              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'get_factor_and_xyz_min_max'
     567              : 
     568              :       INTEGER                                            :: handle, nkp
     569              :       INTEGER, DIMENSION(3)                              :: periodic
     570              : 
     571          112 :       CALL timeset(routineN, handle)
     572              : 
     573          112 :       CALL get_kpoint_info(kpoints, nkp_grid=nkp_grid, nkp=nkp)
     574          112 :       CALL get_cell(cell=cell, h=hmat, periodic=periodic)
     575              : 
     576          112 :       IF (periodic(1) == 0) THEN
     577           86 :          CPASSERT(nkp_grid(1) == 1)
     578              :       END IF
     579          112 :       IF (periodic(2) == 0) THEN
     580           20 :          CPASSERT(nkp_grid(2) == 1)
     581              :       END IF
     582          112 :       IF (periodic(3) == 0) THEN
     583           26 :          CPASSERT(nkp_grid(3) == 1)
     584              :       END IF
     585              : 
     586          112 :       IF (MODULO(nkp_grid(1), 2) == 1) THEN
     587           86 :          factor = 3**(size_lattice_sum - 1)
     588           26 :       ELSE IF (MODULO(nkp_grid(1), 2) == 0) THEN
     589           26 :          factor = 2**(size_lattice_sum - 1)
     590              :       END IF
     591              : 
     592          112 :       IF (MODULO(nkp_grid(1), 2) == 1) THEN
     593           86 :          x_min = -(factor*nkp_grid(1) - 1)/2
     594           86 :          x_max = (factor*nkp_grid(1) - 1)/2
     595           26 :       ELSE IF (MODULO(nkp_grid(1), 2) == 0) THEN
     596           26 :          x_min = -factor*nkp_grid(1)/2
     597           26 :          x_max = factor*nkp_grid(1)/2 - 1
     598              :       END IF
     599          112 :       IF (periodic(1) == 0) THEN
     600           86 :          x_min = 0
     601           86 :          x_max = 0
     602              :       END IF
     603              : 
     604          112 :       IF (MODULO(nkp_grid(2), 2) == 1) THEN
     605           20 :          y_min = -(factor*nkp_grid(2) - 1)/2
     606           20 :          y_max = (factor*nkp_grid(2) - 1)/2
     607           92 :       ELSE IF (MODULO(nkp_grid(2), 2) == 0) THEN
     608           92 :          y_min = -factor*nkp_grid(2)/2
     609           92 :          y_max = factor*nkp_grid(2)/2 - 1
     610              :       END IF
     611          112 :       IF (periodic(2) == 0) THEN
     612           20 :          y_min = 0
     613           20 :          y_max = 0
     614              :       END IF
     615              : 
     616          112 :       IF (MODULO(nkp_grid(3), 2) == 1) THEN
     617           26 :          z_min = -(factor*nkp_grid(3) - 1)/2
     618           26 :          z_max = (factor*nkp_grid(3) - 1)/2
     619           86 :       ELSE IF (MODULO(nkp_grid(3), 2) == 0) THEN
     620           86 :          z_min = -factor*nkp_grid(3)/2
     621           86 :          z_max = factor*nkp_grid(3)/2 - 1
     622              :       END IF
     623          112 :       IF (periodic(3) == 0) THEN
     624           26 :          z_min = 0
     625           26 :          z_max = 0
     626              :       END IF
     627              : 
     628          112 :       CALL timestop(handle)
     629              : 
     630          112 :    END SUBROUTINE
     631              : 
     632              : ! **************************************************************************************************
     633              : !> \brief ...
     634              : !> \param matrix_v_L_tmp ...
     635              : !> \param matrix_v_kp ...
     636              : !> \param ikp_start ...
     637              : ! **************************************************************************************************
     638          100 :    SUBROUTINE allocate_tmp(matrix_v_L_tmp, matrix_v_kp, ikp_start)
     639              : 
     640              :       TYPE(dbcsr_type), POINTER                          :: matrix_v_L_tmp
     641              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_v_kp
     642              :       INTEGER                                            :: ikp_start
     643              : 
     644              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'allocate_tmp'
     645              : 
     646              :       INTEGER                                            :: handle
     647              : 
     648          100 :       CALL timeset(routineN, handle)
     649              : 
     650          100 :       NULLIFY (matrix_v_L_tmp)
     651          100 :       CALL dbcsr_init_p(matrix_v_L_tmp)
     652              :       CALL dbcsr_create(matrix=matrix_v_L_tmp, &
     653              :                         template=matrix_v_kp(ikp_start, 1)%matrix, &
     654          100 :                         matrix_type=dbcsr_type_no_symmetry)
     655          100 :       CALL dbcsr_reserve_all_blocks(matrix_v_L_tmp)
     656          100 :       CALL dbcsr_set(matrix_v_L_tmp, 0.0_dp)
     657              : 
     658          100 :       CALL timestop(handle)
     659              : 
     660          100 :    END SUBROUTINE
     661              : 
     662              : ! **************************************************************************************************
     663              : !> \brief ...
     664              : !> \param matrix_v_L_tmp ...
     665              : ! **************************************************************************************************
     666          100 :    SUBROUTINE deallocate_tmp(matrix_v_L_tmp)
     667              : 
     668              :       TYPE(dbcsr_type), POINTER                          :: matrix_v_L_tmp
     669              : 
     670              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'deallocate_tmp'
     671              : 
     672              :       INTEGER                                            :: handle
     673              : 
     674          100 :       CALL timeset(routineN, handle)
     675              : 
     676          100 :       CALL dbcsr_release_p(matrix_v_L_tmp)
     677              : 
     678          100 :       CALL timestop(handle)
     679              : 
     680          100 :    END SUBROUTINE
     681              : 
     682              : ! **************************************************************************************************
     683              : !> \brief ...
     684              : !> \param V_k ...
     685              : !> \param qs_env ...
     686              : !> \param kpoints ...
     687              : !> \param size_lattice_sum ...
     688              : !> \param basis_type ...
     689              : !> \param ikp_start ...
     690              : !> \param ikp_end ...
     691              : ! **************************************************************************************************
     692           12 :    SUBROUTINE build_2c_coulomb_matrix_kp_small_cell(V_k, qs_env, kpoints, size_lattice_sum, &
     693              :                                                     basis_type, ikp_start, ikp_end)
     694              :       COMPLEX(KIND=dp), DIMENSION(:, :, :)               :: V_k
     695              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     696              :       TYPE(kpoint_type), POINTER                         :: kpoints
     697              :       INTEGER                                            :: size_lattice_sum
     698              :       CHARACTER(LEN=*), INTENT(IN)                       :: basis_type
     699              :       INTEGER                                            :: ikp_start, ikp_end
     700              : 
     701              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'build_2c_coulomb_matrix_kp_small_cell'
     702              : 
     703              :       INTEGER :: factor, handle, handle2, i_cell, i_x, i_x_inner, i_x_outer, ik, ikp_local, j_y, &
     704              :          j_y_inner, j_y_outer, k_z, k_z_inner, k_z_outer, n_atom, n_bf, x_max, x_min, y_max, &
     705              :          y_min, z_max, z_min
     706           12 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: bf_end_from_atom, bf_start_from_atom
     707              :       INTEGER, DIMENSION(3)                              :: nkp_grid
     708              :       REAL(KIND=dp)                                      :: coskl, sinkl
     709              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: V_L
     710              :       REAL(KIND=dp), DIMENSION(3)                        :: vec_L, vec_s
     711              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: hmat
     712           12 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     713              :       TYPE(cell_type), POINTER                           :: cell
     714              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     715           12 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     716           12 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     717              : 
     718           12 :       CALL timeset(routineN, handle)
     719              : 
     720              :       CALL get_qs_env(qs_env=qs_env, &
     721              :                       para_env=para_env, &
     722              :                       particle_set=particle_set, &
     723              :                       cell=cell, &
     724              :                       qs_kind_set=qs_kind_set, &
     725           12 :                       atomic_kind_set=atomic_kind_set)
     726              : 
     727              :       CALL get_factor_and_xyz_min_max(cell, kpoints, size_lattice_sum, factor, hmat, &
     728           12 :                                       x_min, x_max, y_min, y_max, z_min, z_max, nkp_grid)
     729              : 
     730           12 :       CALL get_basis_sizes(qs_env, n_atom, basis_type, bf_start_from_atom, bf_end_from_atom, n_bf)
     731              : 
     732           48 :       ALLOCATE (V_L(n_bf, n_bf))
     733              : 
     734          220 :       DO i_x_inner = 0, 2*nkp_grid(1) - 1
     735        11228 :          DO j_y_inner = 0, 2*nkp_grid(2) - 1
     736        72656 :             DO k_z_inner = 0, 2*nkp_grid(3) - 1
     737              : 
     738      2396160 :                V_L(:, :) = 0.0_dp
     739        61440 :                i_cell = 0
     740              : 
     741       163840 :                DO i_x_outer = x_min, x_max + nkp_grid(1), 2*nkp_grid(1)
     742       552960 :                   DO j_y_outer = y_min, y_max + nkp_grid(2), 2*nkp_grid(2)
     743      1699840 :                      DO k_z_outer = z_min, z_max + nkp_grid(3), 2*nkp_grid(3)
     744              : 
     745      1208320 :                         i_x = i_x_inner + i_x_outer
     746      1208320 :                         j_y = j_y_inner + j_y_outer
     747      1208320 :                         k_z = k_z_inner + k_z_outer
     748              : 
     749              :                         IF (i_x > x_max .OR. i_x < x_min .OR. &
     750              :                             j_y > y_max .OR. j_y < y_min .OR. &
     751      1208320 :                             k_z > z_max .OR. k_z < z_min) CYCLE
     752              : 
     753       455680 :                         i_cell = i_cell + 1
     754              : 
     755      1822720 :                         vec_s = [REAL(i_x, dp), REAL(j_y, dp), REAL(k_z, dp)]
     756              : 
     757       455680 :                         IF (MODULO(i_cell, para_env%num_pe) .NE. para_env%mepos) CYCLE
     758              : 
     759      2961920 :                         vec_L = MATMUL(hmat, vec_s)
     760              : 
     761              :                         ! Compute (P 0 | Q vec_L) and add it to V_R
     762              :                         CALL add_V_L(V_L, vec_L, n_atom, bf_start_from_atom, bf_end_from_atom, &
     763      1597440 :                                      particle_set, qs_kind_set, atomic_kind_set, basis_type, cell)
     764              : 
     765              :                      END DO
     766              :                   END DO
     767              :                END DO
     768              : 
     769        61440 :                CALL para_env%sync()
     770        61440 :                CALL para_env%sum(V_L)
     771              : 
     772        61440 :                CALL timeset(routineN//"_R_to_k", handle2)
     773              : 
     774        61440 :                ikp_local = 0
     775              : 
     776              :                ! add exp(iq*vec_L) * (P 0 | Q vec_L) to V_PQ(q)
     777     33091584 :                DO ik = 1, ikp_end
     778              : 
     779     33030144 :                   IF (MODULO(ik, para_env%num_pe) .NE. para_env%mepos) CYCLE
     780              : 
     781     16515072 :                   ikp_local = ikp_local + 1
     782              : 
     783     16515072 :                   IF (ik < ikp_start) CYCLE
     784              : 
     785              :                   ! coskl and sinkl are identical for all i_x_outer, j_y_outer, k_z_outer
     786     53477376 :                   coskl = COS(twopi*DOT_PRODUCT(vec_s(1:3), kpoints%xkp(1:3, ik)))
     787     53477376 :                   sinkl = SIN(twopi*DOT_PRODUCT(vec_s(1:3), kpoints%xkp(1:3, ik)))
     788              : 
     789              :                   V_k(:, :, ikp_local) = V_k(:, :, ikp_local) + z_one*coskl*V_L(:, :) + &
     790    524611584 :                                          gaussi*sinkl*V_L(:, :)
     791              : 
     792              :                END DO
     793              : 
     794       133888 :                CALL timestop(handle2)
     795              : 
     796              :             END DO
     797              :          END DO
     798              :       END DO
     799              : 
     800           12 :       CALL timestop(handle)
     801              : 
     802           24 :    END SUBROUTINE build_2c_coulomb_matrix_kp_small_cell
     803              : 
     804              : ! **************************************************************************************************
     805              : !> \brief ...
     806              : !> \param qs_env ...
     807              : !> \param n_atom ...
     808              : !> \param basis_type ...
     809              : !> \param bf_start_from_atom ...
     810              : !> \param bf_end_from_atom ...
     811              : !> \param n_bf ...
     812              : ! **************************************************************************************************
     813           12 :    SUBROUTINE get_basis_sizes(qs_env, n_atom, basis_type, bf_start_from_atom, bf_end_from_atom, n_bf)
     814              : 
     815              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     816              :       INTEGER                                            :: n_atom
     817              :       CHARACTER(LEN=*), INTENT(IN)                       :: basis_type
     818              :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: bf_start_from_atom, bf_end_from_atom
     819              :       INTEGER                                            :: n_bf
     820              : 
     821              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'get_basis_sizes'
     822              : 
     823              :       INTEGER                                            :: handle, iatom, ikind, n_kind, nsgf
     824           12 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: kind_of
     825           12 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     826              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a
     827           12 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     828           12 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     829              : 
     830           12 :       CALL timeset(routineN, handle)
     831              : 
     832              :       CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, &
     833           12 :                       qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set)
     834           12 :       CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
     835              : 
     836           12 :       n_atom = SIZE(particle_set)
     837           12 :       n_kind = SIZE(qs_kind_set)
     838              : 
     839           36 :       DO ikind = 1, n_kind
     840              :          CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, &
     841           24 :                           basis_type=basis_type)
     842           36 :          CPASSERT(ASSOCIATED(basis_set_a))
     843              :       END DO
     844              : 
     845           60 :       ALLOCATE (bf_start_from_atom(n_atom), bf_end_from_atom(n_atom))
     846              : 
     847           12 :       n_bf = 0
     848           44 :       DO iatom = 1, n_atom
     849           32 :          bf_start_from_atom(iatom) = n_bf + 1
     850           32 :          ikind = kind_of(iatom)
     851           32 :          CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type=basis_type)
     852           32 :          n_bf = n_bf + nsgf
     853           44 :          bf_end_from_atom(iatom) = n_bf
     854              :       END DO
     855              : 
     856           12 :       CALL timestop(handle)
     857              : 
     858           24 :    END SUBROUTINE get_basis_sizes
     859              : 
     860              : ! **************************************************************************************************
     861              : !> \brief ...
     862              : !> \param V_L ...
     863              : !> \param vec_L ...
     864              : !> \param n_atom ...
     865              : !> \param bf_start_from_atom ...
     866              : !> \param bf_end_from_atom ...
     867              : !> \param particle_set ...
     868              : !> \param qs_kind_set ...
     869              : !> \param atomic_kind_set ...
     870              : !> \param basis_type ...
     871              : !> \param cell ...
     872              : ! **************************************************************************************************
     873       227840 :    SUBROUTINE add_V_L(V_L, vec_L, n_atom, bf_start_from_atom, bf_end_from_atom, &
     874              :                       particle_set, qs_kind_set, atomic_kind_set, basis_type, cell)
     875              : 
     876              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: V_L
     877              :       REAL(KIND=dp), DIMENSION(3)                        :: vec_L
     878              :       INTEGER                                            :: n_atom
     879              :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: bf_start_from_atom, bf_end_from_atom
     880              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     881              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     882              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     883              :       CHARACTER(LEN=*), INTENT(IN)                       :: basis_type
     884              :       TYPE(cell_type), POINTER                           :: cell
     885              : 
     886              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'add_V_L'
     887              : 
     888              :       INTEGER                                            :: a_1, a_2, atom_a, atom_b, b_1, b_2, &
     889              :                                                             handle, kind_a, kind_b
     890       227840 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: kind_of
     891              :       REAL(dp), DIMENSION(3)                             :: ra, rab_L, rb
     892       227840 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: V_L_ab
     893       227840 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: contr_a, contr_b
     894              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
     895              : 
     896       227840 :       CALL timeset(routineN, handle)
     897              : 
     898       227840 :       NULLIFY (basis_set_a, basis_set_b)
     899              : 
     900       227840 :       CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
     901              : 
     902       890880 :       DO atom_a = 1, n_atom
     903              : 
     904      2839040 :          DO atom_b = 1, n_atom
     905              : 
     906      1948160 :             kind_a = kind_of(atom_a)
     907      1948160 :             kind_b = kind_of(atom_b)
     908              : 
     909              :             CALL get_qs_kind(qs_kind=qs_kind_set(kind_a), basis_set=basis_set_a, &
     910      1948160 :                              basis_type=basis_type)
     911              :             CALL get_qs_kind(qs_kind=qs_kind_set(kind_b), basis_set=basis_set_b, &
     912      1948160 :                              basis_type=basis_type)
     913              : 
     914      1948160 :             ra(1:3) = pbc(particle_set(atom_a)%r(1:3), cell)
     915      1948160 :             rb(1:3) = pbc(particle_set(atom_b)%r(1:3), cell)
     916              : 
     917      7792640 :             rab_L(1:3) = rb(1:3) - ra(1:3) + vec_L(1:3)
     918              : 
     919      1948160 :             CALL contraction_matrix_shg(basis_set_a, contr_a)
     920      1948160 :             CALL contraction_matrix_shg(basis_set_b, contr_b)
     921              : 
     922      1948160 :             a_1 = bf_start_from_atom(atom_a)
     923      1948160 :             a_2 = bf_end_from_atom(atom_a)
     924      1948160 :             b_1 = bf_start_from_atom(atom_b)
     925      1948160 :             b_2 = bf_end_from_atom(atom_b)
     926              : 
     927      7792640 :             ALLOCATE (V_L_ab(a_2 - a_1 + 1, b_2 - b_1 + 1))
     928              : 
     929              :             CALL int_operators_r12_ab_shg(operator_coulomb, V_L_ab, rab=rab_L, &
     930              :                                           fba=basis_set_a, fbb=basis_set_b, &
     931              :                                           scona_shg=contr_a, sconb_shg=contr_b, &
     932      1948160 :                                           calculate_forces=.FALSE.)
     933              : 
     934     13862400 :             V_L(a_1:a_2, b_1:b_2) = V_L(a_1:a_2, b_1:b_2) + V_L_ab(:, :)
     935              : 
     936      2611200 :             DEALLOCATE (contr_a, contr_b, V_L_ab)
     937              : 
     938              :          END DO
     939              : 
     940              :       END DO
     941              : 
     942       227840 :       DEALLOCATE (kind_of)
     943              : 
     944       227840 :       CALL timestop(handle)
     945              : 
     946       227840 :    END SUBROUTINE add_V_L
     947              : 
     948            0 : END MODULE kpoint_coulomb_2c
        

Generated by: LCOV version 2.0-1