LCOV - code coverage report
Current view: top level - src - core_ppl.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:15a58fb) Lines: 112 116 96.6 %
Date: 2025-02-18 08:24:35 Functions: 2 2 100.0 %

          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       17367 :    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       17367 :                              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       17367 :       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       17367 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, &
     130       17367 :                                                             nct_lpot, npgfa, npgfb, nsgfa, nsgfb
     131       17367 :       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       17367 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: work
     139       17367 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: hab2_w, ppl_fwork, ppl_work
     140       17367 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: hab, pab
     141             :       REAL(KIND=dp), ALLOCATABLE, &
     142       17367 :          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       17367 :          DIMENSION(:), POINTER                           :: ap_iterator
     149             :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
     150       17367 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
     151             :       TYPE(gth_potential_type), POINTER                  :: gth_potential
     152       34734 :       REAL(KIND=dp), DIMENSION(SIZE(particle_set))       :: at_thread
     153             :       REAL(KIND=dp), DIMENSION(nexp_max)                 :: alpha_ppl
     154       17367 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: cval_lpot, h1_1block, h1_2block, &
     155       17367 :                                                             h1_3block, h_block, p_block, rpgfa, &
     156       17367 :                                                             rpgfb, sphi_a, sphi_b, zeta, zetb
     157       17367 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: a_local, alpha_lpot, c_local, cexp_ppl, &
     158       17367 :                                                             set_radius_a, set_radius_b
     159             :       REAL(KIND=dp), DIMENSION(4, nexp_max)              :: cval_ppl
     160       34734 :       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       17367 : !$       ALLOCATABLE, DIMENSION(:) :: locks
     165             : !$    INTEGER                                            :: lock_num, hash, hash1, hash2
     166             : !$    INTEGER(KIND=int_8)                                :: iatom8
     167             : !$    INTEGER, PARAMETER                                 :: nlock = 501
     168             : 
     169       17367 :       do_dR = PRESENT(deltaR)
     170       17367 :       doat = PRESENT(atcore)
     171             :       MARK_USED(int_8)
     172             : 
     173             :       ! Use internal integral routine for local ECP terms or use libgrrp
     174       17367 :       libgrpp_local = .FALSE.
     175             : 
     176       34734 :       IF (calculate_forces) THEN
     177        7139 :          CALL timeset(routineN//"_forces", handle)
     178             :       ELSE
     179       10228 :          CALL timeset(routineN, handle)
     180             :       END IF
     181             : 
     182       17367 :       nkind = SIZE(atomic_kind_set)
     183       17367 :       natom = SIZE(particle_set)
     184             : 
     185       17367 :       dokp = (nimages > 1)
     186             : 
     187       17367 :       IF (dokp) THEN
     188         234 :          CPASSERT(PRESENT(cell_to_index) .AND. ASSOCIATED(cell_to_index))
     189             :       END IF
     190             : 
     191       17367 :       IF (calculate_forces .OR. doat) THEN
     192        7201 :          IF (SIZE(matrix_p, 1) == 2) THEN
     193        2414 :             DO img = 1, nimages
     194             :                CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
     195        1490 :                               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        2414 :                               alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
     198             :             END DO
     199             :          END IF
     200             :       END IF
     201      276647 :       force_thread = 0.0_dp
     202       82187 :       at_thread = 0.0_dp
     203             : 
     204       17367 :       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       17367 :                            basis_type=basis_type)
     209             : 
     210       17367 :       maxl = MAX(maxlgto, maxlppl)
     211       17367 :       CALL init_orbital_pointers(2*maxl + 2*nder + 1)
     212             : 
     213       17367 :       ldsab = MAX(maxco, ncoset(maxlppl), maxsgf, maxlppl)
     214       17367 :       ldai = ncoset(maxl + nder + 1)
     215             : 
     216       82922 :       ALLOCATE (basis_set_list(nkind))
     217       48188 :       DO ikind = 1, nkind
     218       30821 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a, basis_type=basis_type)
     219       48188 :          IF (ASSOCIATED(basis_set_a)) THEN
     220       30821 :             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       17367 :       pv_thread = 0.0_dp
     227             : 
     228             :       nthread = 1
     229       17367 : !$    nthread = omp_get_max_threads()
     230             : 
     231             :       ! iterator for basis/potential list
     232       17367 :       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, &
     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       17367 : !$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       17367 :       CALL neighbor_list_iterator_release(ap_iterator)
     727             : 
     728       17367 :       DEALLOCATE (basis_set_list)
     729             : 
     730       17367 :       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        7201 :          IF (SIZE(matrix_p, 1) == 2) THEN
     734        2414 :             DO img = 1, nimages
     735             :                CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
     736        1490 :                               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        2414 :                               alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
     739             :             END DO
     740             :          END IF
     741             :       END IF
     742             : 
     743       17367 :       IF (calculate_forces) THEN
     744        7139 :          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       25333 :             atom_a = atom_of_kind(iatom)
     748       25333 :             ikind = kind_of(iatom)
     749      101332 :             force(ikind)%gth_ppl(:, atom_a) = force(ikind)%gth_ppl(:, atom_a) + force_thread(:, iatom)
     750             :          END DO
     751             : !$OMP END DO
     752        7139 :          DEALLOCATE (atom_of_kind, kind_of)
     753             :       END IF
     754       17367 :       IF (doat) THEN
     755         280 :          atcore(1:natom) = atcore(1:natom) + at_thread(1:natom)
     756             :       END IF
     757             : 
     758       17367 :       IF (calculate_forces .AND. use_virial) THEN
     759       10634 :          virial%pv_ppl = virial%pv_ppl + pv_thread
     760       10634 :          virial%pv_virial = virial%pv_virial + pv_thread
     761             :       END IF
     762             : 
     763       17367 :       CALL timestop(handle)
     764             : 
     765       52101 :    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,&
     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 1.15