LCOV - code coverage report
Current view: top level - src - core_ppl.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 96.6 % 116 112
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 2 2

            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              : !> \brief Calculation of the local pseudopotential contribution to the core Hamiltonian
       9              : !>         <a|V(local)|b> = <a|Sum e^a*rc**2|b>
      10              : !> \par History
      11              : !>      - core_ppnl refactored from qs_core_hamiltonian [Joost VandeVondele, 2008-11-01]
      12              : !>      - adapted for PPL [jhu, 2009-02-23]
      13              : !>      - OpenMP added [Iain Bethune, Fiona Reid, 2013-11-13]
      14              : !>      - Bug fix: correct orbital pointer range [07.2014,JGH]
      15              : !>      - k-point aware [07.2015,JGH]
      16              : !>      - Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
      17              : ! **************************************************************************************************
      18              : MODULE core_ppl
      19              : 
      20              :    USE ai_overlap_ppl,                  ONLY: ecploc_integral,&
      21              :                                               ppl_integral,&
      22              :                                               ppl_integral_ri
      23              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      24              :                                               get_atomic_kind_set
      25              :    USE basis_set_types,                 ONLY: gto_basis_set_p_type,&
      26              :                                               gto_basis_set_type
      27              :    USE cp_dbcsr_api,                    ONLY: dbcsr_add,&
      28              :                                               dbcsr_get_block_p,&
      29              :                                               dbcsr_p_type
      30              :    USE external_potential_types,        ONLY: get_potential,&
      31              :                                               gth_potential_type,&
      32              :                                               sgp_potential_type
      33              :    USE kinds,                           ONLY: dp,&
      34              :                                               int_8
      35              :    USE libgrpp_integrals,               ONLY: libgrpp_local_forces_ref,&
      36              :                                               libgrpp_local_integrals,&
      37              :                                               libgrpp_semilocal_forces_ref,&
      38              :                                               libgrpp_semilocal_integrals
      39              :    USE lri_environment_types,           ONLY: lri_kind_type
      40              :    USE orbital_pointers,                ONLY: init_orbital_pointers,&
      41              :                                               ncoset
      42              :    USE particle_types,                  ONLY: particle_type
      43              :    USE qs_force_types,                  ONLY: qs_force_type
      44              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      45              :                                               get_qs_kind_set,&
      46              :                                               qs_kind_type
      47              :    USE qs_neighbor_list_types,          ONLY: get_iterator_info,&
      48              :                                               neighbor_list_iterator_create,&
      49              :                                               neighbor_list_iterator_p_type,&
      50              :                                               neighbor_list_iterator_release,&
      51              :                                               neighbor_list_set_p_type,&
      52              :                                               nl_set_sub_iterator,&
      53              :                                               nl_sub_iterate
      54              :    USE virial_methods,                  ONLY: virial_pair_force
      55              :    USE virial_types,                    ONLY: virial_type
      56              : 
      57              : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
      58              : !$ USE OMP_LIB, ONLY: omp_lock_kind, &
      59              : !$                    omp_init_lock, omp_set_lock, &
      60              : !$                    omp_unset_lock, omp_destroy_lock
      61              : 
      62              : #include "./base/base_uses.f90"
      63              : 
      64              :    IMPLICIT NONE
      65              : 
      66              :    PRIVATE
      67              : 
      68              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'core_ppl'
      69              : 
      70              :    PUBLIC :: build_core_ppl, build_core_ppl_ri
      71              : 
      72              : CONTAINS
      73              : 
      74              : ! **************************************************************************************************
      75              : !> \brief ...
      76              : !> \param matrix_h ...
      77              : !> \param matrix_p ...
      78              : !> \param force ...
      79              : !> \param virial ...
      80              : !> \param calculate_forces ...
      81              : !> \param use_virial ...
      82              : !> \param nder ...
      83              : !> \param qs_kind_set ...
      84              : !> \param atomic_kind_set ...
      85              : !> \param particle_set ...
      86              : !> \param sab_orb ...
      87              : !> \param sac_ppl ...
      88              : !> \param nimages ...
      89              : !> \param cell_to_index ...
      90              : !> \param basis_type ...
      91              : !> \param deltaR Weighting factors of the derivatives wrt. nuclear positions
      92              : !> \param atcore ...
      93              : ! **************************************************************************************************
      94        17695 :    SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
      95              :                              qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
      96        17695 :                              nimages, cell_to_index, basis_type, deltaR, atcore)
      97              : 
      98              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_h, matrix_p
      99              :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     100              :       TYPE(virial_type), POINTER                         :: virial
     101              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     102              :       LOGICAL                                            :: use_virial
     103              :       INTEGER                                            :: nder
     104              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     105              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     106              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     107              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     108              :          POINTER                                         :: sab_orb, sac_ppl
     109              :       INTEGER, INTENT(IN)                                :: nimages
     110              :       INTEGER, DIMENSION(:, :, :), OPTIONAL, POINTER     :: cell_to_index
     111              :       CHARACTER(LEN=*), INTENT(IN)                       :: basis_type
     112              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
     113              :          OPTIONAL                                        :: deltaR
     114              :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
     115              :          OPTIONAL                                        :: atcore
     116              : 
     117              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'build_core_ppl'
     118              :       INTEGER, PARAMETER                                 :: nexp_max = 30
     119              : 
     120              :       INTEGER :: atom_a, handle, i, iatom, icol, ikind, img, irow, iset, jatom, jkind, jset, &
     121              :          katom, kkind, ldai, ldsab, maxco, maxder, maxl, maxlgto, maxlppl, maxnset, maxsgf, mepos, &
     122              :          n_local, natom, ncoa, ncob, nexp_lpot, nexp_ppl, nkind, nloc, nseta, nsetb, nthread, &
     123              :          sgfa, sgfb, slmax, slot
     124        17695 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, kind_of
     125              :       INTEGER, DIMENSION(0:10)                           :: npot
     126              :       INTEGER, DIMENSION(1:10)                           :: nrloc
     127              :       INTEGER, DIMENSION(1:15, 0:10)                     :: nrpot
     128              :       INTEGER, DIMENSION(3)                              :: cellind
     129        17695 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, &
     130        17695 :                                                             nct_lpot, npgfa, npgfb, nsgfa, nsgfb
     131        17695 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
     132              :       INTEGER, DIMENSION(nexp_max)                       :: nct_ppl
     133              :       LOGICAL                                            :: do_dR, doat, dokp, ecp_local, &
     134              :                                                             ecp_semi_local, found, libgrpp_local, &
     135              :                                                             lpotextended, only_gaussians
     136              :       REAL(KIND=dp)                                      :: alpha, atk0, atk1, dab, dac, dbc, f0, &
     137              :                                                             ppl_radius
     138        17695 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: work
     139        17695 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: hab2_w, ppl_fwork, ppl_work
     140        17695 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: hab, pab
     141              :       REAL(KIND=dp), ALLOCATABLE, &
     142        17695 :          DIMENSION(:, :, :, :, :)                        :: hab2
     143              :       REAL(KIND=dp), DIMENSION(1:10)                     :: aloc, bloc
     144              :       REAL(KIND=dp), DIMENSION(1:15, 0:10)               :: apot, bpot
     145              :       REAL(KIND=dp), DIMENSION(3)                        :: force_a, force_b, rab, rac, rbc
     146              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_thread
     147              :       TYPE(neighbor_list_iterator_p_type), &
     148        17695 :          DIMENSION(:), POINTER                           :: ap_iterator
     149              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
     150        17695 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
     151              :       TYPE(gth_potential_type), POINTER                  :: gth_potential
     152        35390 :       REAL(KIND=dp), DIMENSION(SIZE(particle_set))       :: at_thread
     153              :       REAL(KIND=dp), DIMENSION(nexp_max)                 :: alpha_ppl
     154        17695 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: cval_lpot, h1_1block, h1_2block, &
     155        17695 :                                                             h1_3block, h_block, p_block, rpgfa, &
     156        17695 :                                                             rpgfb, sphi_a, sphi_b, zeta, zetb
     157        17695 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: a_local, alpha_lpot, c_local, cexp_ppl, &
     158        17695 :                                                             set_radius_a, set_radius_b
     159              :       REAL(KIND=dp), DIMENSION(4, nexp_max)              :: cval_ppl
     160        35390 :       REAL(KIND=dp), DIMENSION(3, SIZE(particle_set))    :: force_thread
     161              :       TYPE(sgp_potential_type), POINTER                  :: sgp_potential
     162              : 
     163              : !$    INTEGER(kind=omp_lock_kind), &
     164        17695 : !$       ALLOCATABLE, DIMENSION(:) :: locks
     165              : !$    INTEGER                                            :: lock_num, hash, hash1, hash2
     166              : !$    INTEGER(KIND=int_8)                                :: iatom8
     167              : !$    INTEGER, PARAMETER                                 :: nlock = 501
     168              : 
     169        17695 :       do_dR = PRESENT(deltaR)
     170        17695 :       doat = PRESENT(atcore)
     171              :       MARK_USED(int_8)
     172              : 
     173              :       ! Use internal integral routine for local ECP terms or use libgrrp
     174        17695 :       libgrpp_local = .FALSE.
     175              : 
     176        35390 :       IF (calculate_forces) THEN
     177         7327 :          CALL timeset(routineN//"_forces", handle)
     178              :       ELSE
     179        10368 :          CALL timeset(routineN, handle)
     180              :       END IF
     181              : 
     182        17695 :       nkind = SIZE(atomic_kind_set)
     183        17695 :       natom = SIZE(particle_set)
     184              : 
     185        17695 :       dokp = (nimages > 1)
     186              : 
     187        17695 :       IF (dokp) THEN
     188          238 :          CPASSERT(PRESENT(cell_to_index) .AND. ASSOCIATED(cell_to_index))
     189              :       END IF
     190              : 
     191        17695 :       IF (calculate_forces .OR. doat) THEN
     192         7389 :          IF (SIZE(matrix_p, 1) == 2) THEN
     193         2454 :             DO img = 1, nimages
     194              :                CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
     195         1528 :                               alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
     196              :                CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
     197         2454 :                               alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
     198              :             END DO
     199              :          END IF
     200              :       END IF
     201       281647 :       force_thread = 0.0_dp
     202        83683 :       at_thread = 0.0_dp
     203              : 
     204        17695 :       maxder = ncoset(nder)
     205              : 
     206              :       CALL get_qs_kind_set(qs_kind_set, maxco=maxco, maxlgto=maxlgto, &
     207              :                            maxsgf=maxsgf, maxnset=maxnset, maxlppl=maxlppl, &
     208        17695 :                            basis_type=basis_type)
     209              : 
     210        17695 :       maxl = MAX(maxlgto, maxlppl)
     211        17695 :       CALL init_orbital_pointers(2*maxl + 2*nder + 1)
     212              : 
     213        17695 :       ldsab = MAX(maxco, ncoset(maxlppl), maxsgf, maxlppl)
     214        17695 :       ldai = ncoset(maxl + nder + 1)
     215              : 
     216        84534 :       ALLOCATE (basis_set_list(nkind))
     217        49144 :       DO ikind = 1, nkind
     218        31449 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a, basis_type=basis_type)
     219        49144 :          IF (ASSOCIATED(basis_set_a)) THEN
     220        31449 :             basis_set_list(ikind)%gto_basis_set => basis_set_a
     221              :          ELSE
     222            0 :             NULLIFY (basis_set_list(ikind)%gto_basis_set)
     223              :          END IF
     224              :       END DO
     225              : 
     226        17695 :       pv_thread = 0.0_dp
     227              : 
     228              :       nthread = 1
     229        17695 : !$    nthread = omp_get_max_threads()
     230              : 
     231              :       ! iterator for basis/potential list
     232        17695 :       CALL neighbor_list_iterator_create(ap_iterator, sac_ppl, search=.TRUE., nthread=nthread)
     233              : 
     234              : !$OMP PARALLEL &
     235              : !$OMP DEFAULT (NONE) &
     236              : !$OMP SHARED  (ap_iterator, basis_set_list, calculate_forces, use_virial, &
     237              : !$OMP          matrix_h, matrix_p, atomic_kind_set, qs_kind_set, particle_set, &
     238              : !$OMP          sab_orb, sac_ppl, nthread, ncoset, nkind, cell_to_index, &
     239              : !$OMP          ldsab,  maxnset, maxder, do_dR, deltaR, doat, libgrpp_local, &
     240              : !$OMP          maxlgto, nder, maxco, dokp, locks, natom) &
     241              : !$OMP PRIVATE (ikind, jkind, iatom, jatom, rab, basis_set_a, basis_set_b, &
     242              : !$OMP          first_sgfa, la_max, la_min, npgfa, nsgfa, sphi_a, &
     243              : !$OMP          zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, sphi_b, &
     244              : !$OMP          zetb, dab, irow, icol, h_block, found, iset, ncoa, lock_num, &
     245              : !$OMP          sgfa, jset, ncob, sgfb, nsgfb, p_block, work, pab, hab, hab2, hab2_w, &
     246              : !$OMP          atk0, atk1, h1_1block, h1_2block, h1_3block, kkind, nseta, &
     247              : !$OMP          gth_potential, sgp_potential, alpha, cexp_ppl, lpotextended, &
     248              : !$OMP          ppl_radius, nexp_lpot, nexp_ppl, alpha_ppl, alpha_lpot, nct_ppl, &
     249              : !$OMP          nct_lpot, cval_ppl, cval_lpot, rac, dac, rbc, dbc, &
     250              : !$OMP          set_radius_a,  rpgfa, force_a, force_b, ppl_fwork, mepos, &
     251              : !$OMP          slot, f0, katom, ppl_work, cellind, img, ecp_local, ecp_semi_local, &
     252              : !$OMP          nloc, nrloc, aloc, bloc, n_local, a_local, c_local, &
     253              : !$OMP          slmax, npot, nrpot, apot, bpot, only_gaussians, &
     254              : !$OMP          ldai, hash, hash1, hash2, iatom8) &
     255        17695 : !$OMP REDUCTION (+ : pv_thread, force_thread, at_thread )
     256              : 
     257              : !$OMP SINGLE
     258              : !$    ALLOCATE (locks(nlock))
     259              : !$OMP END SINGLE
     260              : 
     261              : !$OMP DO
     262              : !$    DO lock_num = 1, nlock
     263              : !$       call omp_init_lock(locks(lock_num))
     264              : !$    END DO
     265              : !$OMP END DO
     266              : 
     267              :       mepos = 0
     268              : !$    mepos = omp_get_thread_num()
     269              : 
     270              :       ALLOCATE (hab(ldsab, ldsab, maxnset, maxnset), work(ldsab, ldsab*maxder))
     271              :       ldai = ncoset(2*maxlgto + 2*nder)
     272              :       ALLOCATE (ppl_work(ldai, ldai, MAX(maxder, 2*maxlgto + 2*nder + 1)))
     273              :       IF (calculate_forces .OR. doat) THEN
     274              :          ALLOCATE (pab(maxco, maxco, maxnset, maxnset))
     275              :          ldai = ncoset(maxlgto)
     276              :          ALLOCATE (ppl_fwork(ldai, ldai, maxder))
     277              :       END IF
     278              : 
     279              : !$OMP DO SCHEDULE(GUIDED)
     280              :       DO slot = 1, sab_orb(1)%nl_size
     281              :          !SL
     282              :          IF (do_dR) THEN
     283              :             ALLOCATE (hab2(ldsab, ldsab, 4, maxnset, maxnset))
     284              :             ALLOCATE (hab2_w(ldsab, ldsab, 6))
     285              :             ALLOCATE (ppl_fwork(ldai, ldai, maxder))
     286              :          END IF
     287              : 
     288              :          ikind = sab_orb(1)%nlist_task(slot)%ikind
     289              :          jkind = sab_orb(1)%nlist_task(slot)%jkind
     290              :          iatom = sab_orb(1)%nlist_task(slot)%iatom
     291              :          jatom = sab_orb(1)%nlist_task(slot)%jatom
     292              :          cellind(:) = sab_orb(1)%nlist_task(slot)%cell(:)
     293              :          rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
     294              : 
     295              :          basis_set_a => basis_set_list(ikind)%gto_basis_set
     296              :          IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
     297              :          basis_set_b => basis_set_list(jkind)%gto_basis_set
     298              :          IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
     299              : 
     300              : !$       iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
     301              : !$       hash1 = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
     302              : 
     303              :          ! basis ikind
     304              :          first_sgfa => basis_set_a%first_sgf
     305              :          la_max => basis_set_a%lmax
     306              :          la_min => basis_set_a%lmin
     307              :          npgfa => basis_set_a%npgf
     308              :          nseta = basis_set_a%nset
     309              :          nsgfa => basis_set_a%nsgf_set
     310              :          rpgfa => basis_set_a%pgf_radius
     311              :          set_radius_a => basis_set_a%set_radius
     312              :          sphi_a => basis_set_a%sphi
     313              :          zeta => basis_set_a%zet
     314              :          ! basis jkind
     315              :          first_sgfb => basis_set_b%first_sgf
     316              :          lb_max => basis_set_b%lmax
     317              :          lb_min => basis_set_b%lmin
     318              :          npgfb => basis_set_b%npgf
     319              :          nsetb = basis_set_b%nset
     320              :          nsgfb => basis_set_b%nsgf_set
     321              :          rpgfb => basis_set_b%pgf_radius
     322              :          set_radius_b => basis_set_b%set_radius
     323              :          sphi_b => basis_set_b%sphi
     324              :          zetb => basis_set_b%zet
     325              : 
     326              :          dab = SQRT(SUM(rab*rab))
     327              : 
     328              :          IF (dokp) THEN
     329              :             img = cell_to_index(cellind(1), cellind(2), cellind(3))
     330              :          ELSE
     331              :             img = 1
     332              :          END IF
     333              : 
     334              :          ! *** Use the symmetry of the first derivatives ***
     335              :          IF (iatom == jatom) THEN
     336              :             f0 = 1.0_dp
     337              :          ELSE
     338              :             f0 = 2.0_dp
     339              :          END IF
     340              : 
     341              :          ! *** Create matrix blocks for a new matrix block column ***
     342              :          IF (iatom <= jatom) THEN
     343              :             irow = iatom
     344              :             icol = jatom
     345              :          ELSE
     346              :             irow = jatom
     347              :             icol = iatom
     348              :          END IF
     349              :          NULLIFY (h_block)
     350              : 
     351              :          IF (do_dR) THEN
     352              :             NULLIFY (h1_1block, h1_2block, h1_3block)
     353              : 
     354              :             CALL dbcsr_get_block_p(matrix=matrix_h(1, img)%matrix, &
     355              :                                    row=irow, col=icol, BLOCK=h1_1block, found=found)
     356              :             CALL dbcsr_get_block_p(matrix=matrix_h(2, img)%matrix, &
     357              :                                    row=irow, col=icol, BLOCK=h1_2block, found=found)
     358              :             CALL dbcsr_get_block_p(matrix=matrix_h(3, img)%matrix, &
     359              :                                    row=irow, col=icol, BLOCK=h1_3block, found=found)
     360              :          END IF
     361              : 
     362              :          CALL dbcsr_get_block_p(matrix_h(1, img)%matrix, irow, icol, h_block, found)
     363              :          CPASSERT(found)
     364              :          IF (calculate_forces .OR. doat) THEN
     365              :             NULLIFY (p_block)
     366              :             CALL dbcsr_get_block_p(matrix_p(1, img)%matrix, irow, icol, p_block, found)
     367              :             IF (ASSOCIATED(p_block)) THEN
     368              :                DO iset = 1, nseta
     369              :                   ncoa = npgfa(iset)*ncoset(la_max(iset))
     370              :                   sgfa = first_sgfa(1, iset)
     371              :                   DO jset = 1, nsetb
     372              :                      ncob = npgfb(jset)*ncoset(lb_max(jset))
     373              :                      sgfb = first_sgfb(1, jset)
     374              : 
     375              :                      ! *** Decontract density matrix block ***
     376              :                      IF (iatom <= jatom) THEN
     377              :                         work(1:ncoa, 1:nsgfb(jset)) = MATMUL(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1), &
     378              :                                                              p_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1))
     379              :                      ELSE
     380              :                         work(1:ncoa, 1:nsgfb(jset)) = MATMUL(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1), &
     381              :                                                        TRANSPOSE(p_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1)))
     382              :                      END IF
     383              : 
     384              :                      pab(1:ncoa, 1:ncob, iset, jset) = MATMUL(work(1:ncoa, 1:nsgfb(jset)), &
     385              :                                                               TRANSPOSE(sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1)))
     386              :                   END DO
     387              :                END DO
     388              :             END IF
     389              :          END IF
     390              : 
     391              :          hab = 0._dp
     392              :          IF (do_dr) hab2 = 0._dp
     393              : 
     394              :          ! loop over all kinds for pseudopotential atoms
     395              :          DO kkind = 1, nkind
     396              : 
     397              :             CALL get_qs_kind(qs_kind_set(kkind), gth_potential=gth_potential, &
     398              :                              sgp_potential=sgp_potential)
     399              :             ecp_semi_local = .FALSE.
     400              :             only_gaussians = .TRUE.
     401              :             IF (ASSOCIATED(gth_potential)) THEN
     402              :                CALL get_potential(potential=gth_potential, &
     403              :                                   alpha_ppl=alpha, cexp_ppl=cexp_ppl, &
     404              :                                   lpot_present=lpotextended, ppl_radius=ppl_radius)
     405              :                nexp_ppl = 1
     406              :                alpha_ppl(1) = alpha
     407              :                nct_ppl(1) = SIZE(cexp_ppl)
     408              :                cval_ppl(1:nct_ppl(1), 1) = cexp_ppl(1:nct_ppl(1))
     409              :                IF (lpotextended) THEN
     410              :                   CALL get_potential(potential=gth_potential, &
     411              :                                      nexp_lpot=nexp_lpot, alpha_lpot=alpha_lpot, nct_lpot=nct_lpot, &
     412              :                                      cval_lpot=cval_lpot)
     413              :                   CPASSERT(nexp_lpot < nexp_max)
     414              :                   nexp_ppl = nexp_lpot + 1
     415              :                   alpha_ppl(2:nexp_lpot + 1) = alpha_lpot(1:nexp_lpot)
     416              :                   nct_ppl(2:nexp_lpot + 1) = nct_lpot(1:nexp_lpot)
     417              :                   DO i = 1, nexp_lpot
     418              :                      cval_ppl(1:nct_lpot(i), i + 1) = cval_lpot(1:nct_lpot(i), i)
     419              :                   END DO
     420              :                END IF
     421              :             ELSE IF (ASSOCIATED(sgp_potential)) THEN
     422              :                CALL get_potential(potential=sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local, &
     423              :                                   ppl_radius=ppl_radius)
     424              :                IF (ecp_local) THEN
     425              :                   CALL get_potential(potential=sgp_potential, nloc=nloc, nrloc=nrloc, aloc=aloc, bloc=bloc)
     426              :                   nexp_ppl = nloc
     427              :                   CPASSERT(nexp_ppl <= nexp_max)
     428              :                   nct_ppl(1:nloc) = nrloc(1:nloc)
     429              :                   alpha_ppl(1:nloc) = bloc(1:nloc)
     430              :                   cval_ppl(1, 1:nloc) = aloc(1:nloc)
     431              :                   only_gaussians = .FALSE.
     432              :                ELSE
     433              :                   CALL get_potential(potential=sgp_potential, n_local=n_local, a_local=a_local, c_local=c_local)
     434              :                   nexp_ppl = n_local
     435              :                   CPASSERT(nexp_ppl <= nexp_max)
     436              :                   nct_ppl(1:n_local) = 1
     437              :                   alpha_ppl(1:n_local) = a_local(1:n_local)
     438              :                   cval_ppl(1, 1:n_local) = c_local(1:n_local)
     439              :                END IF
     440              :                IF (ecp_semi_local) THEN
     441              :                   CALL get_potential(potential=sgp_potential, sl_lmax=slmax, &
     442              :                                      npot=npot, nrpot=nrpot, apot=apot, bpot=bpot)
     443              :                ELSEIF (ecp_local) THEN
     444              :                   IF (SUM(ABS(aloc(1:nloc))) < 1.0e-12_dp) CYCLE
     445              :                END IF
     446              :             ELSE
     447              :                CYCLE
     448              :             END IF
     449              : 
     450              :             CALL nl_set_sub_iterator(ap_iterator, ikind, kkind, iatom, mepos=mepos)
     451              : 
     452              :             DO WHILE (nl_sub_iterate(ap_iterator, mepos=mepos) == 0)
     453              : 
     454              :                CALL get_iterator_info(ap_iterator, mepos=mepos, jatom=katom, r=rac)
     455              : 
     456              :                dac = SQRT(SUM(rac*rac))
     457              :                rbc(:) = rac(:) - rab(:)
     458              :                dbc = SQRT(SUM(rbc*rbc))
     459              :                IF ((MAXVAL(set_radius_a(:)) + ppl_radius < dac) .OR. &
     460              :                    (MAXVAL(set_radius_b(:)) + ppl_radius < dbc)) THEN
     461              :                   CYCLE
     462              :                END IF
     463              : 
     464              :                DO iset = 1, nseta
     465              :                   IF (set_radius_a(iset) + ppl_radius < dac) CYCLE
     466              :                   ncoa = npgfa(iset)*ncoset(la_max(iset))
     467              :                   sgfa = first_sgfa(1, iset)
     468              :                   DO jset = 1, nsetb
     469              :                      IF (set_radius_b(jset) + ppl_radius < dbc) CYCLE
     470              :                      ncob = npgfb(jset)*ncoset(lb_max(jset))
     471              :                      sgfb = first_sgfb(1, jset)
     472              :                      IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
     473              :                      ! *** Calculate the GTH pseudo potential forces ***
     474              :                      IF (doat) THEN
     475              :                         atk0 = f0*SUM(hab(1:ncoa, 1:ncob, iset, jset)* &
     476              :                                       pab(1:ncoa, 1:ncob, iset, jset))
     477              :                      END IF
     478              :                      IF (calculate_forces) THEN
     479              : 
     480              :                         force_a(:) = 0.0_dp
     481              :                         force_b(:) = 0.0_dp
     482              : 
     483              :                         IF (only_gaussians) THEN
     484              :                            CALL ppl_integral( &
     485              :                               la_max(iset), la_min(iset), npgfa(iset), &
     486              :                               rpgfa(:, iset), zeta(:, iset), &
     487              :                               lb_max(jset), lb_min(jset), npgfb(jset), &
     488              :                               rpgfb(:, jset), zetb(:, jset), &
     489              :                               nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
     490              :                               rab, dab, rac, dac, rbc, dbc, &
     491              :                               hab(:, :, iset, jset), ppl_work, pab(:, :, iset, jset), &
     492              :                               force_a, force_b, ppl_fwork)
     493              :                         ELSEIF (libgrpp_local) THEN
     494              : !$OMP CRITICAL(type1)
     495              :                            CALL libgrpp_local_forces_ref(la_max(iset), la_min(iset), npgfa(iset), &
     496              :                                                          rpgfa(:, iset), zeta(:, iset), &
     497              :                                                          lb_max(jset), lb_min(jset), npgfb(jset), &
     498              :                                                          rpgfb(:, jset), zetb(:, jset), &
     499              :                                                          nexp_ppl, alpha_ppl, cval_ppl(1, :), nct_ppl, &
     500              :                                                          ppl_radius, rab, dab, rac, dac, dbc, &
     501              :                                                          hab(:, :, iset, jset), pab(:, :, iset, jset), &
     502              :                                                          force_a, force_b)
     503              : !$OMP END CRITICAL(type1)
     504              :                         ELSE
     505              :                            CALL ecploc_integral( &
     506              :                               la_max(iset), la_min(iset), npgfa(iset), &
     507              :                               rpgfa(:, iset), zeta(:, iset), &
     508              :                               lb_max(jset), lb_min(jset), npgfb(jset), &
     509              :                               rpgfb(:, jset), zetb(:, jset), &
     510              :                               nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
     511              :                               rab, dab, rac, dac, rbc, dbc, &
     512              :                               hab(:, :, iset, jset), ppl_work, pab(:, :, iset, jset), &
     513              :                               force_a, force_b, ppl_fwork)
     514              :                         END IF
     515              : 
     516              :                         IF (ecp_semi_local) THEN
     517              : 
     518              : !$OMP CRITICAL(type2)
     519              :                            CALL libgrpp_semilocal_forces_ref(la_max(iset), la_min(iset), npgfa(iset), &
     520              :                                                              rpgfa(:, iset), zeta(:, iset), &
     521              :                                                              lb_max(jset), lb_min(jset), npgfb(jset), &
     522              :                                                              rpgfb(:, jset), zetb(:, jset), &
     523              :                                                              slmax, npot, bpot, apot, nrpot, &
     524              :                                                              ppl_radius, rab, dab, rac, dac, dbc, &
     525              :                                                              hab(:, :, iset, jset), pab(:, :, iset, jset), &
     526              :                                                              force_a, force_b)
     527              : !$OMP END CRITICAL(type2)
     528              :                         END IF
     529              :                         ! *** The derivatives w.r.t. atomic center c are    ***
     530              :                         ! *** calculated using the translational invariance ***
     531              :                         ! *** of the first derivatives                      ***
     532              : 
     533              :                         force_thread(1, iatom) = force_thread(1, iatom) + f0*force_a(1)
     534              :                         force_thread(2, iatom) = force_thread(2, iatom) + f0*force_a(2)
     535              :                         force_thread(3, iatom) = force_thread(3, iatom) + f0*force_a(3)
     536              :                         force_thread(1, katom) = force_thread(1, katom) - f0*force_a(1)
     537              :                         force_thread(2, katom) = force_thread(2, katom) - f0*force_a(2)
     538              :                         force_thread(3, katom) = force_thread(3, katom) - f0*force_a(3)
     539              : 
     540              :                         force_thread(1, jatom) = force_thread(1, jatom) + f0*force_b(1)
     541              :                         force_thread(2, jatom) = force_thread(2, jatom) + f0*force_b(2)
     542              :                         force_thread(3, jatom) = force_thread(3, jatom) + f0*force_b(3)
     543              :                         force_thread(1, katom) = force_thread(1, katom) - f0*force_b(1)
     544              :                         force_thread(2, katom) = force_thread(2, katom) - f0*force_b(2)
     545              :                         force_thread(3, katom) = force_thread(3, katom) - f0*force_b(3)
     546              : 
     547              :                         IF (use_virial) THEN
     548              :                            CALL virial_pair_force(pv_thread, f0, force_a, rac)
     549              :                            CALL virial_pair_force(pv_thread, f0, force_b, rbc)
     550              :                         END IF
     551              :                      ELSEIF (do_dR) THEN
     552              :                         hab2_w = 0._dp
     553              :                         CALL ppl_integral( &
     554              :                            la_max(iset), la_min(iset), npgfa(iset), &
     555              :                            rpgfa(:, iset), zeta(:, iset), &
     556              :                            lb_max(jset), lb_min(jset), npgfb(jset), &
     557              :                            rpgfb(:, jset), zetb(:, jset), &
     558              :                            nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
     559              :                            rab, dab, rac, dac, rbc, dbc, &
     560              :                            vab=hab(:, :, iset, jset), s=ppl_work, &
     561              :                            hab2=hab2(:, :, :, iset, jset), hab2_work=hab2_w, fs=ppl_fwork, &
     562              :                            deltaR=deltaR, iatom=iatom, jatom=jatom, katom=katom)
     563              :                         IF (ecp_semi_local) THEN
     564              :                            ! semi local ECP part
     565              :                            CPABORT("Option not implemented")
     566              :                         END IF
     567              :                      ELSE
     568              :                         IF (only_gaussians) THEN
     569              :                            !If the local part of the pseudo-potential only has Gaussian functions
     570              :                            !we can use CP2K native code, that can run without libgrpp installation
     571              :                            CALL ppl_integral( &
     572              :                               la_max(iset), la_min(iset), npgfa(iset), &
     573              :                               rpgfa(:, iset), zeta(:, iset), &
     574              :                               lb_max(jset), lb_min(jset), npgfb(jset), &
     575              :                               rpgfb(:, jset), zetb(:, jset), &
     576              :                               nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
     577              :                               rab, dab, rac, dac, rbc, dbc, hab(:, :, iset, jset), ppl_work)
     578              : 
     579              :                         ELSEIF (libgrpp_local) THEN
     580              :                            !If the local part of the potential is more complex, we need libgrpp
     581              : !$OMP CRITICAL(type1)
     582              :                            CALL libgrpp_local_integrals(la_max(iset), la_min(iset), npgfa(iset), &
     583              :                                                         rpgfa(:, iset), zeta(:, iset), &
     584              :                                                         lb_max(jset), lb_min(jset), npgfb(jset), &
     585              :                                                         rpgfb(:, jset), zetb(:, jset), &
     586              :                                                         nexp_ppl, alpha_ppl, cval_ppl(1, :), nct_ppl, &
     587              :                                                         ppl_radius, rab, dab, rac, dac, dbc, &
     588              :                                                         hab(:, :, iset, jset))
     589              : !$OMP END CRITICAL(type1)
     590              :                         ELSE
     591              :                            CALL ecploc_integral( &
     592              :                               la_max(iset), la_min(iset), npgfa(iset), &
     593              :                               rpgfa(:, iset), zeta(:, iset), &
     594              :                               lb_max(jset), lb_min(jset), npgfb(jset), &
     595              :                               rpgfb(:, jset), zetb(:, jset), &
     596              :                               nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
     597              :                               rab, dab, rac, dac, rbc, dbc, hab(:, :, iset, jset), ppl_work)
     598              :                         END IF
     599              : 
     600              :                         IF (ecp_semi_local) THEN
     601              :                            ! semi local ECP part
     602              : !$OMP CRITICAL(type2)
     603              :                            CALL libgrpp_semilocal_integrals(la_max(iset), la_min(iset), npgfa(iset), &
     604              :                                                             rpgfa(:, iset), zeta(:, iset), &
     605              :                                                             lb_max(jset), lb_min(jset), npgfb(jset), &
     606              :                                                             rpgfb(:, jset), zetb(:, jset), &
     607              :                                                             slmax, npot, bpot, apot, nrpot, &
     608              :                                                             ppl_radius, rab, dab, rac, dac, dbc, &
     609              :                                                             hab(:, :, iset, jset))
     610              : !$OMP END CRITICAL(type2)
     611              :                         END IF
     612              :                      END IF
     613              :                      ! calculate atomic contributions
     614              :                      IF (doat) THEN
     615              :                         atk1 = f0*SUM(hab(1:ncoa, 1:ncob, iset, jset)* &
     616              :                                       pab(1:ncoa, 1:ncob, iset, jset))
     617              :                         at_thread(katom) = at_thread(katom) + (atk1 - atk0)
     618              :                      END IF
     619              :                   END DO
     620              :                END DO
     621              :             END DO
     622              :          END DO
     623              : 
     624              :          ! *** Contract PPL integrals
     625              :          IF (.NOT. do_dR) THEN
     626              :          DO iset = 1, nseta
     627              :             ncoa = npgfa(iset)*ncoset(la_max(iset))
     628              :             sgfa = first_sgfa(1, iset)
     629              :             DO jset = 1, nsetb
     630              :                ncob = npgfb(jset)*ncoset(lb_max(jset))
     631              :                sgfb = first_sgfb(1, jset)
     632              : 
     633              : !$             hash2 = MOD((iset - 1)*nsetb + jset, nlock) + 1
     634              : !$             hash = MOD(hash1 + hash2, nlock) + 1
     635              : 
     636              :                work(1:ncoa, 1:nsgfb(jset)) = MATMUL(hab(1:ncoa, 1:ncob, iset, jset), &
     637              :                                                     sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
     638              : !$             CALL omp_set_lock(locks(hash))
     639              :                IF (iatom <= jatom) THEN
     640              :                   h_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
     641              :                      h_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
     642              :                      MATMUL(TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
     643              :                ELSE
     644              :                   h_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
     645              :                      h_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
     646              :                      MATMUL(TRANSPOSE(work(1:ncoa, 1:nsgfb(jset))), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
     647              :                END IF
     648              : !$             CALL omp_unset_lock(locks(hash))
     649              : 
     650              :             END DO
     651              :          END DO
     652              :          ELSE  ! do_dr == .true.
     653              :          DO iset = 1, nseta
     654              :             ncoa = npgfa(iset)*ncoset(la_max(iset))
     655              :             sgfa = first_sgfa(1, iset)
     656              :             DO jset = 1, nsetb
     657              :                ncob = npgfb(jset)*ncoset(lb_max(jset))
     658              :                sgfb = first_sgfb(1, jset)
     659              :                work(1:ncoa, 1:nsgfb(jset)) = MATMUL(hab2(1:ncoa, 1:ncob, 1, iset, jset), &
     660              :                                                     sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
     661              : 
     662              : !$OMP CRITICAL(h1_1block_critical)
     663              :                IF (iatom <= jatom) THEN
     664              :                   h1_1block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
     665              :                      h1_1block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
     666              :                      MATMUL(TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
     667              : 
     668              :                ELSE
     669              :                   h1_1block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
     670              :                      h1_1block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
     671              :                      MATMUL(TRANSPOSE(work(1:ncoa, 1:nsgfb(jset))), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
     672              :                END IF
     673              : !$OMP END CRITICAL(h1_1block_critical)
     674              :                work(1:ncoa, 1:nsgfb(jset)) = MATMUL(hab2(1:ncoa, 1:ncob, 2, iset, jset), &
     675              :                                                     sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
     676              : 
     677              : !$OMP CRITICAL(h1_2block_critical)
     678              :                IF (iatom <= jatom) THEN
     679              :                   h1_2block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
     680              :                      h1_2block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
     681              :                      MATMUL(TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
     682              : 
     683              :                ELSE
     684              :                   h1_2block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
     685              :                      h1_2block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
     686              :                      MATMUL(TRANSPOSE(work(1:ncoa, 1:nsgfb(jset))), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
     687              :                END IF
     688              : !$OMP END CRITICAL(h1_2block_critical)
     689              :                work(1:ncoa, 1:nsgfb(jset)) = MATMUL(hab2(1:ncoa, 1:ncob, 3, iset, jset), &
     690              :                                                     sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
     691              : !$OMP CRITICAL(h1_3block_critical)
     692              :                IF (iatom <= jatom) THEN
     693              :                   h1_3block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
     694              :                      h1_3block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
     695              :                      MATMUL(TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
     696              : 
     697              :                ELSE
     698              :                   h1_3block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
     699              :                      h1_3block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
     700              :                      MATMUL(TRANSPOSE(work(1:ncoa, 1:nsgfb(jset))), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
     701              :                END IF
     702              : !$OMP END CRITICAL(h1_3block_critical)
     703              :             END DO
     704              :          END DO
     705              :          END IF
     706              :          IF (do_dR) DEALLOCATE (hab2, ppl_fwork, hab2_w)
     707              :       END DO ! slot
     708              : 
     709              :       DEALLOCATE (hab, work, ppl_work)
     710              :       IF (calculate_forces .OR. doat) THEN
     711              :          DEALLOCATE (pab, ppl_fwork)
     712              :       END IF
     713              : 
     714              : !$OMP DO
     715              : !$    DO lock_num = 1, nlock
     716              : !$       call omp_destroy_lock(locks(lock_num))
     717              : !$    END DO
     718              : !$OMP END DO
     719              : 
     720              : !$OMP SINGLE
     721              : !$    DEALLOCATE (locks)
     722              : !$OMP END SINGLE NOWAIT
     723              : 
     724              : !$OMP END PARALLEL
     725              : 
     726        17695 :       CALL neighbor_list_iterator_release(ap_iterator)
     727              : 
     728        17695 :       DEALLOCATE (basis_set_list)
     729              : 
     730        17695 :       IF (calculate_forces .OR. doat) THEN
     731              :          ! *** If LSD, then recover alpha density and beta density     ***
     732              :          ! *** from the total density (1) and the spin density (2)     ***
     733         7389 :          IF (SIZE(matrix_p, 1) == 2) THEN
     734         2454 :             DO img = 1, nimages
     735              :                CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
     736         1528 :                               alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
     737              :                CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
     738         2454 :                               alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
     739              :             END DO
     740              :          END IF
     741              :       END IF
     742              : 
     743        17695 :       IF (calculate_forces) THEN
     744         7327 :          CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
     745              : !$OMP DO
     746              :          DO iatom = 1, natom
     747        25947 :             atom_a = atom_of_kind(iatom)
     748        25947 :             ikind = kind_of(iatom)
     749       103788 :             force(ikind)%gth_ppl(:, atom_a) = force(ikind)%gth_ppl(:, atom_a) + force_thread(:, iatom)
     750              :          END DO
     751              : !$OMP END DO
     752         7327 :          DEALLOCATE (atom_of_kind, kind_of)
     753              :       END IF
     754        17695 :       IF (doat) THEN
     755          280 :          atcore(1:natom) = atcore(1:natom) + at_thread(1:natom)
     756              :       END IF
     757              : 
     758        17695 :       IF (calculate_forces .AND. use_virial) THEN
     759        10920 :          virial%pv_ppl = virial%pv_ppl + pv_thread
     760        10920 :          virial%pv_virial = virial%pv_virial + pv_thread
     761              :       END IF
     762              : 
     763        17695 :       CALL timestop(handle)
     764              : 
     765        53085 :    END SUBROUTINE build_core_ppl
     766              : 
     767              : ! **************************************************************************************************
     768              : !> \brief ...
     769              : !> \param lri_ppl_coef ...
     770              : !> \param force ...
     771              : !> \param virial ...
     772              : !> \param calculate_forces ...
     773              : !> \param use_virial ...
     774              : !> \param qs_kind_set ...
     775              : !> \param atomic_kind_set ...
     776              : !> \param particle_set ...
     777              : !> \param sac_ppl ...
     778              : !> \param basis_type ...
     779              : ! **************************************************************************************************
     780            4 :    SUBROUTINE build_core_ppl_ri(lri_ppl_coef, force, virial, calculate_forces, use_virial, &
     781              :                                 qs_kind_set, atomic_kind_set, particle_set, sac_ppl, &
     782              :                                 basis_type)
     783              : 
     784              :       TYPE(lri_kind_type), DIMENSION(:), POINTER         :: lri_ppl_coef
     785              :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     786              :       TYPE(virial_type), POINTER                         :: virial
     787              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     788              :       LOGICAL                                            :: use_virial
     789              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     790              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     791              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     792              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     793              :          POINTER                                         :: sac_ppl
     794              :       CHARACTER(LEN=*), INTENT(IN)                       :: basis_type
     795              : 
     796              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'build_core_ppl_ri'
     797              :       INTEGER, PARAMETER                                 :: nexp_max = 30
     798              : 
     799              :       INTEGER :: atom_a, handle, i, iatom, ikind, iset, katom, kkind, maxco, maxsgf, n_local, &
     800              :          natom, ncoa, nexp_lpot, nexp_ppl, nfun, nkind, nloc, nseta, sgfa, sgfb, slot
     801            4 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, kind_of
     802              :       INTEGER, DIMENSION(1:10)                           :: nrloc
     803            4 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, nct_lpot, npgfa, nsgfa
     804            4 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa
     805              :       INTEGER, DIMENSION(nexp_max)                       :: nct_ppl
     806              :       LOGICAL                                            :: ecp_local, ecp_semi_local, lpotextended
     807              :       REAL(KIND=dp)                                      :: alpha, dac, ppl_radius
     808            4 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: va, work
     809            4 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: dva, dvas
     810              :       REAL(KIND=dp), DIMENSION(1:10)                     :: aloc, bloc
     811              :       REAL(KIND=dp), DIMENSION(3)                        :: force_a, rac
     812              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_thread
     813              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set
     814            4 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
     815              :       TYPE(gth_potential_type), POINTER                  :: gth_potential
     816              :       REAL(KIND=dp), DIMENSION(nexp_max)                 :: alpha_ppl
     817            4 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: bcon, cval_lpot, rpgfa, sphi_a, zeta
     818            4 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: a_local, alpha_lpot, c_local, cexp_ppl, &
     819            4 :                                                             set_radius_a
     820              :       REAL(KIND=dp), DIMENSION(4, nexp_max)              :: cval_ppl
     821            8 :       REAL(KIND=dp), DIMENSION(3, SIZE(particle_set))    :: force_thread
     822              :       TYPE(sgp_potential_type), POINTER                  :: sgp_potential
     823              : 
     824              : !$    INTEGER(kind=omp_lock_kind), &
     825            4 : !$       ALLOCATABLE, DIMENSION(:) :: locks
     826              : !$    INTEGER                                            :: lock_num, hash
     827              : !$    INTEGER, PARAMETER                                 :: nlock = 501
     828              : 
     829            8 :       IF (calculate_forces) THEN
     830            2 :          CALL timeset(routineN//"_forces", handle)
     831              :       ELSE
     832            2 :          CALL timeset(routineN, handle)
     833              :       END IF
     834              : 
     835            4 :       nkind = SIZE(atomic_kind_set)
     836            4 :       natom = SIZE(particle_set)
     837              : 
     838           52 :       force_thread = 0.0_dp
     839            4 :       pv_thread = 0.0_dp
     840            4 :       CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
     841              : 
     842           20 :       ALLOCATE (basis_set_list(nkind))
     843           12 :       DO ikind = 1, nkind
     844            8 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, basis_type=basis_type)
     845           12 :          IF (ASSOCIATED(basis_set)) THEN
     846            8 :             basis_set_list(ikind)%gto_basis_set => basis_set
     847              :          ELSE
     848            0 :             NULLIFY (basis_set_list(ikind)%gto_basis_set)
     849              :          END IF
     850              :       END DO
     851              : 
     852            4 :       CALL get_qs_kind_set(qs_kind_set, maxco=maxco, maxsgf=maxsgf, basis_type=basis_type)
     853              : 
     854              : !$OMP PARALLEL &
     855              : !$OMP DEFAULT (NONE) &
     856              : !$OMP SHARED  (maxco,maxsgf,basis_set_list,calculate_forces,lri_ppl_coef,qs_kind_set,&
     857              : !$OMP          locks,natom,use_virial,virial,ncoset,atom_of_kind,sac_ppl) &
     858              : !$OMP PRIVATE (ikind,kkind,iatom,katom,atom_a,rac,va,dva,dvas,basis_set,slot,&
     859              : !$OMP          first_sgfa,la_max,la_min,npgfa,nseta,nsgfa,rpgfa,set_radius_a,lock_num,&
     860              : !$OMP          sphi_a,zeta,gth_potential,sgp_potential,alpha,cexp_ppl,lpotextended,ppl_radius,&
     861              : !$OMP          nexp_ppl,alpha_ppl,nct_ppl,cval_ppl,nloc,n_local,nrloc,a_local,aloc,bloc,c_local,nfun,work,&
     862              : !$OMP          hash,dac,force_a,iset,sgfa,sgfb,ncoa,bcon,cval_lpot,nct_lpot,alpha_lpot,nexp_lpot,&
     863              : !$OMP          ecp_local,ecp_semi_local) &
     864            4 : !$OMP REDUCTION (+ : pv_thread, force_thread )
     865              : 
     866              : !$OMP SINGLE
     867              : !$    ALLOCATE (locks(nlock))
     868              : !$OMP END SINGLE
     869              : 
     870              : !$OMP DO
     871              : !$    DO lock_num = 1, nlock
     872              : !$       call omp_init_lock(locks(lock_num))
     873              : !$    END DO
     874              : !$OMP END DO
     875              : 
     876              :       ALLOCATE (va(maxco), work(maxsgf))
     877              :       IF (calculate_forces) THEN
     878              :          ALLOCATE (dva(maxco, 3), dvas(maxco, 3))
     879              :       END IF
     880              : 
     881              : !$OMP DO SCHEDULE(GUIDED)
     882              :       DO slot = 1, sac_ppl(1)%nl_size
     883              : 
     884              :          ikind = sac_ppl(1)%nlist_task(slot)%ikind
     885              :          kkind = sac_ppl(1)%nlist_task(slot)%jkind
     886              :          iatom = sac_ppl(1)%nlist_task(slot)%iatom
     887              :          katom = sac_ppl(1)%nlist_task(slot)%jatom
     888              :          rac(1:3) = sac_ppl(1)%nlist_task(slot)%r(1:3)
     889              :          atom_a = atom_of_kind(iatom)
     890              : 
     891              :          basis_set => basis_set_list(ikind)%gto_basis_set
     892              :          IF (.NOT. ASSOCIATED(basis_set)) CYCLE
     893              : 
     894              :          ! basis ikind
     895              :          first_sgfa => basis_set%first_sgf
     896              :          la_max => basis_set%lmax
     897              :          la_min => basis_set%lmin
     898              :          npgfa => basis_set%npgf
     899              :          nseta = basis_set%nset
     900              :          nsgfa => basis_set%nsgf_set
     901              :          nfun = basis_set%nsgf
     902              :          rpgfa => basis_set%pgf_radius
     903              :          set_radius_a => basis_set%set_radius
     904              :          sphi_a => basis_set%sphi
     905              :          zeta => basis_set%zet
     906              : 
     907              :          CALL get_qs_kind(qs_kind_set(kkind), gth_potential=gth_potential, &
     908              :                           sgp_potential=sgp_potential)
     909              :          ecp_semi_local = .FALSE.
     910              :          IF (ASSOCIATED(gth_potential)) THEN
     911              :             CALL get_potential(potential=gth_potential, &
     912              :                                alpha_ppl=alpha, cexp_ppl=cexp_ppl, &
     913              :                                lpot_present=lpotextended, ppl_radius=ppl_radius)
     914              :             nexp_ppl = 1
     915              :             alpha_ppl(1) = alpha
     916              :             nct_ppl(1) = SIZE(cexp_ppl)
     917              :             cval_ppl(1:nct_ppl(1), 1) = cexp_ppl(1:nct_ppl(1))
     918              :             IF (lpotextended) THEN
     919              :                CALL get_potential(potential=gth_potential, &
     920              :                                   nexp_lpot=nexp_lpot, alpha_lpot=alpha_lpot, nct_lpot=nct_lpot, cval_lpot=cval_lpot)
     921              :                CPASSERT(nexp_lpot < nexp_max)
     922              :                nexp_ppl = nexp_lpot + 1
     923              :                alpha_ppl(2:nexp_lpot + 1) = alpha_lpot(1:nexp_lpot)
     924              :                nct_ppl(2:nexp_lpot + 1) = nct_lpot(1:nexp_lpot)
     925              :                DO i = 1, nexp_lpot
     926              :                   cval_ppl(1:nct_lpot(i), i + 1) = cval_lpot(1:nct_lpot(i), i)
     927              :                END DO
     928              :             END IF
     929              :          ELSE IF (ASSOCIATED(sgp_potential)) THEN
     930              :             CALL get_potential(potential=sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local, &
     931              :                                ppl_radius=ppl_radius)
     932              :             CPASSERT(.NOT. ecp_semi_local)
     933              :             IF (ecp_local) THEN
     934              :                CALL get_potential(potential=sgp_potential, nloc=nloc, nrloc=nrloc, aloc=aloc, bloc=bloc)
     935              :                IF (SUM(ABS(aloc(1:nloc))) < 1.0e-12_dp) CYCLE
     936              :                nexp_ppl = nloc
     937              :                CPASSERT(nexp_ppl <= nexp_max)
     938              :                nct_ppl(1:nloc) = nrloc(1:nloc)
     939              :                alpha_ppl(1:nloc) = bloc(1:nloc)
     940              :                cval_ppl(1, 1:nloc) = aloc(1:nloc)
     941              :             ELSE
     942              :                CALL get_potential(potential=sgp_potential, n_local=n_local, a_local=a_local, c_local=c_local)
     943              :                nexp_ppl = n_local
     944              :                CPASSERT(nexp_ppl <= nexp_max)
     945              :                nct_ppl(1:n_local) = 1
     946              :                alpha_ppl(1:n_local) = a_local(1:n_local)
     947              :                cval_ppl(1, 1:n_local) = c_local(1:n_local)
     948              :             END IF
     949              :          ELSE
     950              :             CYCLE
     951              :          END IF
     952              : 
     953              :          dac = SQRT(SUM(rac*rac))
     954              :          IF ((MAXVAL(set_radius_a(:)) + ppl_radius < dac)) CYCLE
     955              :          IF (calculate_forces) force_a = 0.0_dp
     956              :          work(1:nfun) = 0.0_dp
     957              : 
     958              :          DO iset = 1, nseta
     959              :             IF (set_radius_a(iset) + ppl_radius < dac) CYCLE
     960              :             ! integrals
     961              :             IF (calculate_forces) THEN
     962              :                va = 0.0_dp
     963              :                dva = 0.0_dp
     964              :                CALL ppl_integral_ri( &
     965              :                   la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
     966              :                   nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
     967              :                   -rac, dac, va, dva)
     968              :             ELSE
     969              :                va = 0.0_dp
     970              :                CALL ppl_integral_ri( &
     971              :                   la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
     972              :                   nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
     973              :                   -rac, dac, va)
     974              :             END IF
     975              :             ! contraction
     976              :             sgfa = first_sgfa(1, iset)
     977              :             sgfb = sgfa + nsgfa(iset) - 1
     978              :             ncoa = npgfa(iset)*ncoset(la_max(iset))
     979              :             bcon => sphi_a(1:ncoa, sgfa:sgfb)
     980              :             work(sgfa:sgfb) = MATMUL(TRANSPOSE(bcon), va(1:ncoa))
     981              :             IF (calculate_forces) THEN
     982              :                dvas(1:nsgfa(iset), 1:3) = MATMUL(TRANSPOSE(bcon), dva(1:ncoa, 1:3))
     983              :                force_a(1) = force_a(1) + SUM(lri_ppl_coef(ikind)%acoef(atom_a, sgfa:sgfb)*dvas(1:nsgfa(iset), 1))
     984              :                force_a(2) = force_a(2) + SUM(lri_ppl_coef(ikind)%acoef(atom_a, sgfa:sgfb)*dvas(1:nsgfa(iset), 2))
     985              :                force_a(3) = force_a(3) + SUM(lri_ppl_coef(ikind)%acoef(atom_a, sgfa:sgfb)*dvas(1:nsgfa(iset), 3))
     986              :             END IF
     987              :          END DO
     988              : !$       hash = MOD(iatom, nlock) + 1
     989              : !$       CALL omp_set_lock(locks(hash))
     990              :          lri_ppl_coef(ikind)%v_int(atom_a, 1:nfun) = lri_ppl_coef(ikind)%v_int(atom_a, 1:nfun) + work(1:nfun)
     991              : !$       CALL omp_unset_lock(locks(hash))
     992              :          IF (calculate_forces) THEN
     993              :             force_thread(1, iatom) = force_thread(1, iatom) + force_a(1)
     994              :             force_thread(2, iatom) = force_thread(2, iatom) + force_a(2)
     995              :             force_thread(3, iatom) = force_thread(3, iatom) + force_a(3)
     996              :             force_thread(1, katom) = force_thread(1, katom) - force_a(1)
     997              :             force_thread(2, katom) = force_thread(2, katom) - force_a(2)
     998              :             force_thread(3, katom) = force_thread(3, katom) - force_a(3)
     999              :             IF (use_virial) THEN
    1000              :                CALL virial_pair_force(pv_thread, 1.0_dp, force_a, rac)
    1001              :             END IF
    1002              :          END IF
    1003              :       END DO
    1004              : 
    1005              :       DEALLOCATE (va, work)
    1006              :       IF (calculate_forces) THEN
    1007              :          DEALLOCATE (dva, dvas)
    1008              :       END IF
    1009              : 
    1010              : !$OMP END PARALLEL
    1011              : 
    1012            4 :       IF (calculate_forces) THEN
    1013            8 :          DO iatom = 1, natom
    1014            6 :             atom_a = atom_of_kind(iatom)
    1015            6 :             ikind = kind_of(iatom)
    1016            6 :             force(ikind)%gth_ppl(1, atom_a) = force(ikind)%gth_ppl(1, atom_a) + force_thread(1, iatom)
    1017            6 :             force(ikind)%gth_ppl(2, atom_a) = force(ikind)%gth_ppl(2, atom_a) + force_thread(2, iatom)
    1018            8 :             force(ikind)%gth_ppl(3, atom_a) = force(ikind)%gth_ppl(3, atom_a) + force_thread(3, iatom)
    1019              :          END DO
    1020              :       END IF
    1021            4 :       DEALLOCATE (atom_of_kind, kind_of)
    1022              : 
    1023            4 :       IF (calculate_forces .AND. use_virial) THEN
    1024            0 :          virial%pv_ppl = virial%pv_ppl + pv_thread
    1025            0 :          virial%pv_virial = virial%pv_virial + pv_thread
    1026              :       END IF
    1027              : 
    1028            4 :       DEALLOCATE (basis_set_list)
    1029              : 
    1030            4 :       CALL timestop(handle)
    1031              : 
    1032           12 :    END SUBROUTINE build_core_ppl_ri
    1033              : 
    1034              : ! **************************************************************************************************
    1035              : 
    1036              : END MODULE core_ppl
        

Generated by: LCOV version 2.0-1