LCOV - code coverage report
Current view: top level - src - libgrpp_integrals.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 39.9 % 278 111
Test Date: 2025-07-25 12:55:17 Functions: 50.0 % 4 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              : ! **************************************************************************************************
       9              : !> \brief Local and semi-local ECP integrals using the libgrpp library
      10              : ! **************************************************************************************************
      11              : 
      12              : MODULE libgrpp_integrals
      13              :    USE kinds, ONLY: dp
      14              :    USE mathconstants, ONLY: pi
      15              :    USE ai_derivatives, ONLY: dabdr_noscreen, adbdr, dabdr
      16              :    USE orbital_pointers, ONLY: nco, &
      17              :                                ncoset
      18              : #if defined(__LIBGRPP)
      19              :    USE libgrpp, ONLY: libgrpp_init, libgrpp_type1_integrals, libgrpp_type2_integrals, &
      20              :                       libgrpp_type1_integrals_gradient, libgrpp_type2_integrals_gradient
      21              : #endif
      22              : #include "./base/base_uses.f90"
      23              : 
      24              :    IMPLICIT NONE
      25              :    PRIVATE
      26              : 
      27              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'libgrpp_integrals'
      28              : 
      29              :    PUBLIC :: libgrpp_semilocal_integrals, libgrpp_local_integrals, &
      30              :              libgrpp_local_forces_ref, libgrpp_semilocal_forces_ref
      31              : 
      32              : CONTAINS
      33              : 
      34              : ! **************************************************************************************************
      35              : !> \brief Local ECP integrals using libgrpp
      36              : !> \param la_max_set ...
      37              : !> \param la_min_set ...
      38              : !> \param npgfa ...
      39              : !> \param rpgfa ...
      40              : !> \param zeta ...
      41              : !> \param lb_max_set ...
      42              : !> \param lb_min_set ...
      43              : !> \param npgfb ...
      44              : !> \param rpgfb ...
      45              : !> \param zetb ...
      46              : !> \param npot_ecp ...
      47              : !> \param alpha_ecp ...
      48              : !> \param coeffs_ecp ...
      49              : !> \param nrpot_ecp ...
      50              : !> \param rpgfc ...
      51              : !> \param rab ...
      52              : !> \param dab ...
      53              : !> \param rac ...
      54              : !> \param dac ...
      55              : !> \param dbc ...
      56              : !> \param vab ...
      57              : !> \param pab ...
      58              : !> \param force_a ...
      59              : !> \param force_b ...
      60              : ! **************************************************************************************************
      61            0 :    SUBROUTINE libgrpp_local_integrals(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
      62            0 :                                       lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
      63            0 :                                       npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
      64            0 :                                       rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
      65              : 
      66              :       INTEGER, INTENT(IN)                                :: la_max_set, la_min_set, npgfa
      67              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rpgfa, zeta
      68              :       INTEGER, INTENT(IN)                                :: lb_max_set, lb_min_set, npgfb
      69              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rpgfb, zetb
      70              :       INTEGER, INTENT(IN)                                :: npot_ecp
      71              :       REAL(KIND=dp), DIMENSION(1:npot_ecp), INTENT(IN)   :: alpha_ecp, coeffs_ecp
      72              :       INTEGER, DIMENSION(1:npot_ecp), INTENT(IN)         :: nrpot_ecp
      73              :       REAL(KIND=dp), INTENT(IN)                          :: rpgfc
      74              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
      75              :       REAL(KIND=dp), INTENT(IN)                          :: dab
      76              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rac
      77              :       REAL(KIND=dp), INTENT(IN)                          :: dac
      78              :       REAL(KIND=dp), INTENT(IN)                          :: dbc
      79              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: vab
      80              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
      81              :          OPTIONAL                                        :: pab
      82              :       REAL(KIND=dp), DIMENSION(3), INTENT(INOUT), &
      83              :          OPTIONAL                                        :: force_a, force_b
      84              : 
      85              : #if defined(__LIBGRPP)
      86              :       INTEGER                                            :: a_offset, a_start, b_offset, b_start, i, &
      87              :                                                             ipgf, j, jpgf, li, lj, ncoa, ncob
      88              :       LOGICAL                                            :: calc_forces
      89              :       REAL(dp)                                           :: expi, expj, normi, normj, prefi, prefj, &
      90              :                                                             zeti, zetj, mindist, fac_a, fac_b
      91            0 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: tmp, tmpx, tmpy, tmpz
      92              :       REAL(dp), DIMENSION(3)                             :: ra, rb, rc
      93              : 
      94            0 :       CALL libgrpp_init()
      95              : 
      96            0 :       calc_forces = .FALSE.
      97            0 :       IF (PRESENT(pab) .AND. PRESENT(force_a) .AND. PRESENT(force_b)) calc_forces = .TRUE.
      98              : 
      99              :       IF (calc_forces) THEN
     100              : 
     101              :          !Note: warning against numerical stability of libgrpp gradients. The day the library becomes
     102              :          !      stable, this routine can be used immediatly as is, and the warning removed.
     103              :          CALL cp_warn(__LOCATION__, &
     104              :                       "ECP gradients calculated with the libgrpp library are, to this date, not numerically stable. "// &
     105            0 :                       "Please use the reference routine 'libgrpp_local_forces_ref' instead.")
     106              : 
     107              :          !there is a weird feature of libgrpp gradients, which is such that the gradient is calculated
     108              :          !for a point in space, and not with respect to an atomic center. For example, if atoms A and
     109              :          !B are the same (and C is different), then d<A | U_C | B>/dPx = d<A | U_C | B>/dAx + d<A | U_C | B>/dBx
     110              :          !Because we want the forces on centers A and B seprately, we need a case study on atomic positions
     111              :          !We always calculate the gradient wrt to atomic position of A and B, and we scale accordingly
     112              : 
     113            0 :          mindist = 1.0E-6_dp
     114              :          !If ra != rb != rc
     115            0 :          IF (dab >= mindist .AND. dbc >= mindist .AND. dac >= mindist) THEN
     116              :             fac_a = 1.0_dp
     117              :             fac_b = 1.0_dp
     118              : 
     119              :             !If ra = rb, but ra != rc
     120            0 :          ELSE IF (dab < mindist .AND. dac >= mindist) THEN
     121              :             fac_a = 0.5_dp
     122              :             fac_b = 0.5_dp
     123              : 
     124              :             !IF ra != rb but ra = rc
     125            0 :          ELSE IF (dab >= mindist .AND. dac < mindist) THEN
     126              :             fac_a = 0.5_dp
     127              :             fac_b = 1.0_dp
     128              : 
     129              :             !IF ra != rb but rb = rc
     130            0 :          ELSE IF (dab >= mindist .AND. dbc < mindist) THEN
     131              :             fac_a = 1.0_dp
     132              :             fac_b = 0.5_dp
     133              : 
     134              :             !If all atoms the same --> no force
     135              :          ELSE
     136            0 :             calc_forces = .FALSE.
     137              :          END IF
     138              :       END IF
     139              : 
     140              :       !libgrpp requires absolute positions, not relative ones
     141            0 :       ra(:) = 0.0_dp
     142            0 :       rb(:) = rab(:)
     143            0 :       rc(:) = rac(:)
     144              : 
     145            0 :       ALLOCATE (tmp(nco(la_max_set)*nco(lb_max_set)))
     146            0 :       IF (calc_forces) THEN
     147            0 :          ALLOCATE (tmpx(nco(la_max_set)*nco(lb_max_set)))
     148            0 :          ALLOCATE (tmpy(nco(la_max_set)*nco(lb_max_set)))
     149            0 :          ALLOCATE (tmpz(nco(la_max_set)*nco(lb_max_set)))
     150              :       END IF
     151              : 
     152            0 :       DO ipgf = 1, npgfa
     153            0 :          IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
     154            0 :          zeti = zeta(ipgf)
     155            0 :          a_start = (ipgf - 1)*ncoset(la_max_set)
     156              : 
     157            0 :          DO jpgf = 1, npgfb
     158            0 :             IF (rpgfb(jpgf) + rpgfc < dbc) CYCLE
     159            0 :             IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) CYCLE
     160            0 :             zetj = zetb(jpgf)
     161            0 :             b_start = (jpgf - 1)*ncoset(lb_max_set)
     162              : 
     163            0 :             DO li = la_min_set, la_max_set
     164            0 :                a_offset = a_start + ncoset(li - 1)
     165            0 :                ncoa = nco(li)
     166            0 :                prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
     167            0 :                expi = 0.25_dp*REAL(2*li + 3, dp)
     168            0 :                normi = 1.0_dp/(prefi*zeti**expi)
     169              : 
     170            0 :                DO lj = lb_min_set, lb_max_set
     171            0 :                   b_offset = b_start + ncoset(lj - 1)
     172            0 :                   ncob = nco(lj)
     173            0 :                   prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
     174            0 :                   expj = 0.25_dp*REAL(2*lj + 3, dp)
     175            0 :                   normj = 1.0_dp/(prefj*zetj**expj)
     176              : 
     177            0 :                   tmp(1:ncoa*ncob) = 0.0_dp
     178              :                   !libgrpp implicitely normalizes cartesian Gaussian. In CP2K, we do not, hence
     179              :                   !the 1/norm coefficients for PGFi and PGFj
     180              :                   CALL libgrpp_type1_integrals(ra, li, 1, [normi], [zeti], &
     181              :                                                rb, lj, 1, [normj], [zetj], &
     182              :                                                rc, [npot_ecp], nrpot_ecp, &
     183            0 :                                                coeffs_ecp, alpha_ecp, tmp)
     184              : 
     185              :                   !note: tmp array is in C row-major ordering
     186            0 :                   DO j = 1, ncob
     187            0 :                      DO i = 1, ncoa
     188            0 :                         vab(a_offset + i, b_offset + j) = vab(a_offset + i, b_offset + j) + tmp((i - 1)*ncob + j)
     189              :                      END DO
     190              :                   END DO
     191              : 
     192            0 :                   IF (calc_forces) THEN
     193            0 :                      tmpx(1:ncoa*ncob) = 0.0_dp
     194            0 :                      tmpy(1:ncoa*ncob) = 0.0_dp
     195            0 :                      tmpz(1:ncoa*ncob) = 0.0_dp
     196              : 
     197              :                      !force wrt to atomic position A
     198              :                      CALL libgrpp_type1_integrals_gradient(ra, li, 1, [normi], [zeti], &
     199              :                                                            rb, lj, 1, [normj], [zetj], &
     200              :                                                            rc, [npot_ecp], nrpot_ecp, &
     201              :                                                            coeffs_ecp, alpha_ecp, ra, &
     202            0 :                                                            tmpx, tmpy, tmpz)
     203              : 
     204              :                      !note: tmp array is in C row-major ordering
     205              :                      !note: zero-gradients sometime comes out as NaN, hence tampval==tmpval check
     206            0 :                      DO j = 1, ncob
     207            0 :                         DO i = 1, ncoa
     208            0 :                            force_a(1) = force_a(1) + fac_a*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
     209            0 :                            force_a(2) = force_a(2) + fac_a*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
     210            0 :                            force_a(3) = force_a(3) + fac_a*pab(a_offset + i, b_offset + j)*tmpz((i - 1)*ncob + j)
     211              :                         END DO
     212              :                      END DO
     213              : 
     214            0 :                      tmpx(1:ncoa*ncob) = 0.0_dp
     215            0 :                      tmpy(1:ncoa*ncob) = 0.0_dp
     216            0 :                      tmpz(1:ncoa*ncob) = 0.0_dp
     217              : 
     218              :                      !force wrt to atomic position B
     219              :                      CALL libgrpp_type1_integrals_gradient(ra, li, 1, [normi], [zeti], &
     220              :                                                            rb, lj, 1, [normj], [zetj], &
     221              :                                                            rc, [npot_ecp], nrpot_ecp, &
     222              :                                                            coeffs_ecp, alpha_ecp, rb, &
     223            0 :                                                            tmpx, tmpy, tmpz)
     224              : 
     225              :                      !note: tmp array is in C row-major ordering
     226              :                      !note: zero-gradients sometime comes out as NaN, hence tampval==tmpval check
     227            0 :                      DO j = 1, ncob
     228            0 :                         DO i = 1, ncoa
     229            0 :                            force_b(1) = force_b(1) + fac_b*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
     230            0 :                            force_b(2) = force_b(2) + fac_b*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
     231            0 :                            force_b(3) = force_b(3) + fac_b*pab(a_offset + i, b_offset + j)*tmpz((i - 1)*ncob + j)
     232              :                         END DO
     233              :                      END DO
     234              :                   END IF
     235              : 
     236              :                END DO !lj
     237              :             END DO !li
     238              : 
     239              :          END DO !jpgf
     240              :       END DO !ipgf
     241              : #else
     242              : 
     243              :       MARK_USED(la_max_set)
     244              :       MARK_USED(la_min_set)
     245              :       MARK_USED(npgfa)
     246              :       MARK_USED(rpgfa)
     247              :       MARK_USED(zeta)
     248              :       MARK_USED(lb_max_set)
     249              :       MARK_USED(lb_min_set)
     250              :       MARK_USED(npgfb)
     251              :       MARK_USED(rpgfb)
     252              :       MARK_USED(zetb)
     253              :       MARK_USED(npot_ecp)
     254              :       MARK_USED(alpha_ecp)
     255              :       MARK_USED(coeffs_ecp)
     256              :       MARK_USED(nrpot_ecp)
     257              :       MARK_USED(rpgfc)
     258              :       MARK_USED(rab)
     259              :       MARK_USED(dab)
     260              :       MARK_USED(rac)
     261              :       MARK_USED(dac)
     262              :       MARK_USED(dbc)
     263              :       MARK_USED(vab)
     264              :       MARK_USED(pab)
     265              :       MARK_USED(force_a)
     266              :       MARK_USED(force_b)
     267              : 
     268              :       CPABORT("Please compile CP2K with libgrpp support for calculations with ECPs")
     269              : #endif
     270              : 
     271            0 :    END SUBROUTINE libgrpp_local_integrals
     272              : 
     273              : ! **************************************************************************************************
     274              : !> \brief Semi-local ECP integrals using libgrpp.
     275              : !> \param la_max_set ...
     276              : !> \param la_min_set ...
     277              : !> \param npgfa ...
     278              : !> \param rpgfa ...
     279              : !> \param zeta ...
     280              : !> \param lb_max_set ...
     281              : !> \param lb_min_set ...
     282              : !> \param npgfb ...
     283              : !> \param rpgfb ...
     284              : !> \param zetb ...
     285              : !> \param lmax_ecp ...
     286              : !> \param npot_ecp ...
     287              : !> \param alpha_ecp ...
     288              : !> \param coeffs_ecp ...
     289              : !> \param nrpot_ecp ...
     290              : !> \param rpgfc ...
     291              : !> \param rab ...
     292              : !> \param dab ...
     293              : !> \param rac ...
     294              : !> \param dac ...
     295              : !> \param dbc ...
     296              : !> \param vab ...
     297              : !> \param pab ...
     298              : !> \param force_a ...
     299              : !> \param force_b ...
     300              : ! **************************************************************************************************
     301         9712 :    SUBROUTINE libgrpp_semilocal_integrals(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
     302         9712 :                                           lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
     303              :                                           lmax_ecp, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
     304        19424 :                                           rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
     305              : 
     306              :       INTEGER, INTENT(IN)                                :: la_max_set, la_min_set, npgfa
     307              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rpgfa, zeta
     308              :       INTEGER, INTENT(IN)                                :: lb_max_set, lb_min_set, npgfb
     309              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rpgfb, zetb
     310              :       INTEGER, INTENT(IN)                                :: lmax_ecp
     311              :       INTEGER, DIMENSION(0:10), INTENT(IN)               :: npot_ecp
     312              :       REAL(KIND=dp), DIMENSION(1:15, 0:10), INTENT(IN)   :: alpha_ecp, coeffs_ecp
     313              :       INTEGER, DIMENSION(1:15, 0:10), INTENT(IN)         :: nrpot_ecp
     314              :       REAL(KIND=dp), INTENT(IN)                          :: rpgfc
     315              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     316              :       REAL(KIND=dp), INTENT(IN)                          :: dab
     317              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rac
     318              :       REAL(KIND=dp), INTENT(IN)                          :: dac
     319              :       REAL(KIND=dp), INTENT(IN)                          :: dbc
     320              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: vab
     321              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
     322              :          OPTIONAL                                        :: pab
     323              :       REAL(KIND=dp), DIMENSION(3), INTENT(INOUT), &
     324              :          OPTIONAL                                        :: force_a, force_b
     325              : 
     326              : #if defined(__LIBGRPP)
     327              :       INTEGER                                            :: a_offset, a_start, b_offset, b_start, i, &
     328              :                                                             ipgf, j, jpgf, li, lj, lk, ncoa, ncob
     329              :       LOGICAL                                            :: calc_forces
     330              :       REAL(dp)                                           :: expi, expj, normi, normj, prefi, prefj, &
     331              :                                                             zeti, zetj, mindist, fac_a, fac_b
     332         9712 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: tmp, tmpx, tmpz, tmpy
     333              :       REAL(dp), DIMENSION(3)                             :: ra, rb, rc
     334              : 
     335         9712 :       CALL libgrpp_init()
     336              : 
     337         9712 :       calc_forces = .FALSE.
     338         9712 :       IF (PRESENT(pab) .AND. PRESENT(force_a) .AND. PRESENT(force_b)) calc_forces = .TRUE.
     339              : 
     340              :       IF (calc_forces) THEN
     341              : 
     342              :          !Note: warning against numerical stability of libgrpp gradients. The day the library becomes
     343              :          !      stable, this routine can be used immediatly as is, and the warning removed.
     344              :          CALL cp_warn(__LOCATION__, &
     345              :                       "ECP gradients calculated with the libgrpp library are, to this date, not numerically stable. "// &
     346            0 :                       "Please use the reference routine 'libgrpp_semilocal_forces_ref' instead.")
     347              : 
     348              :          !there is a weird feature of libgrpp gradients, which is such that the gradient is calculated
     349              :          !for a point in space, and not with respect to an atomic center. For example, if atoms A and
     350              :          !B are the same (and C is different), then d<A | U_C | B>/dPx = d<A | U_C | B>/dAx + d<A | U_C | B>/dBx
     351              :          !Because we want the forces on centers A and B seprately, we need a case study on atomic positions
     352              :          !We always calculate the gradient wrt to atomic position of A and B, and we scale accordingly
     353              : 
     354            0 :          mindist = 1.0E-6_dp
     355              :          !If ra != rb != rc
     356            0 :          IF (dab >= mindist .AND. dbc >= mindist .AND. dac >= mindist) THEN
     357              :             fac_a = 1.0_dp
     358              :             fac_b = 1.0_dp
     359              : 
     360              :             !If ra = rb, but ra != rc
     361            0 :          ELSE IF (dab < mindist .AND. dac >= mindist) THEN
     362              :             fac_a = 0.5_dp
     363              :             fac_b = 0.5_dp
     364              : 
     365              :             !IF ra != rb but ra = rc
     366            0 :          ELSE IF (dab >= mindist .AND. dac < mindist) THEN
     367              :             fac_a = 0.5_dp
     368              :             fac_b = 1.0_dp
     369              : 
     370              :             !IF ra != rb but rb = rc
     371            0 :          ELSE IF (dab >= mindist .AND. dbc < mindist) THEN
     372              :             fac_a = 1.0_dp
     373              :             fac_b = 0.5_dp
     374              : 
     375              :             !If all atoms the same --> no force
     376              :          ELSE
     377            0 :             calc_forces = .FALSE.
     378              :          END IF
     379              :       END IF
     380              : 
     381              :       !libgrpp requires absolute positions, not relative ones
     382         9712 :       ra(:) = 0.0_dp
     383         9712 :       rb(:) = rab(:)
     384         9712 :       rc(:) = rac(:)
     385              : 
     386        29136 :       ALLOCATE (tmp(nco(la_max_set)*nco(lb_max_set)))
     387         9712 :       IF (calc_forces) THEN
     388            0 :          ALLOCATE (tmpx(nco(la_max_set)*nco(lb_max_set)))
     389            0 :          ALLOCATE (tmpy(nco(la_max_set)*nco(lb_max_set)))
     390            0 :          ALLOCATE (tmpz(nco(la_max_set)*nco(lb_max_set)))
     391              :       END IF
     392              : 
     393        23817 :       DO ipgf = 1, npgfa
     394        14105 :          IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
     395        11683 :          zeti = zeta(ipgf)
     396        11683 :          a_start = (ipgf - 1)*ncoset(la_max_set)
     397              : 
     398        39365 :          DO jpgf = 1, npgfb
     399        17970 :             IF (rpgfb(jpgf) + rpgfc < dbc) CYCLE
     400        15236 :             IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) CYCLE
     401        14328 :             zetj = zetb(jpgf)
     402        14328 :             b_start = (jpgf - 1)*ncoset(lb_max_set)
     403              : 
     404        42921 :             DO li = la_min_set, la_max_set
     405        14488 :                a_offset = a_start + ncoset(li - 1)
     406        14488 :                ncoa = nco(li)
     407        14488 :                prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
     408        14488 :                expi = 0.25_dp*REAL(2*li + 3, dp)
     409        14488 :                normi = 1.0_dp/(prefi*zeti**expi)
     410              : 
     411        47362 :                DO lj = lb_min_set, lb_max_set
     412        14904 :                   b_offset = b_start + ncoset(lj - 1)
     413        14904 :                   ncob = nco(lj)
     414        14904 :                   prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
     415        14904 :                   expj = 0.25_dp*REAL(2*lj + 3, dp)
     416        14904 :                   normj = 1.0_dp/(prefj*zetj**expj)
     417              : 
     418              :                   !Loop over ECP angular momentum
     419        90836 :                   DO lk = 0, lmax_ecp
     420       529916 :                      tmp(1:ncoa*ncob) = 0.0_dp
     421              :                      !libgrpp implicitely normalizes cartesian Gaussian. In CP2K, we do not, hence
     422              :                      !the 1/norm coefficients for PGFi and PGFj
     423              :                      CALL libgrpp_type2_integrals(ra, li, 1, [normi], [zeti], &
     424              :                                                   rb, lj, 1, [normj], [zetj], &
     425              :                                                   rc, lk, [npot_ecp(lk)], nrpot_ecp(:, lk), &
     426       368664 :                                                   coeffs_ecp(:, lk), alpha_ecp(:, lk), tmp)
     427              : 
     428              :                      !note: tmp array is in C row-major ordering
     429       231435 :                      DO j = 1, ncob
     430       699907 :                         DO i = 1, ncoa
     431       638463 :                            vab(a_offset + i, b_offset + j) = vab(a_offset + i, b_offset + j) + tmp((i - 1)*ncob + j)
     432              :                         END DO
     433              :                      END DO
     434              : 
     435        76348 :                      IF (calc_forces) THEN
     436              : 
     437            0 :                         tmpx(1:ncoa*ncob) = 0.0_dp
     438            0 :                         tmpy(1:ncoa*ncob) = 0.0_dp
     439            0 :                         tmpz(1:ncoa*ncob) = 0.0_dp
     440              : 
     441              :                         !force wrt to atomic position A
     442              :                         CALL libgrpp_type2_integrals_gradient(ra, li, 1, [normi], [zeti], &
     443              :                                                               rb, lj, 1, [normj], [zetj], &
     444              :                                                               rc, lk, [npot_ecp(lk)], nrpot_ecp(:, lk), &
     445              :                                                               coeffs_ecp(:, lk), alpha_ecp(:, lk), ra, &
     446            0 :                                                               tmpx, tmpy, tmpz)
     447              : 
     448              :                         !note: tmp array is in C row-major ordering
     449              :                         !note: zero-gradients sometime comes out as NaN, hence tampval==tmpval check
     450            0 :                         DO j = 1, ncob
     451            0 :                            DO i = 1, ncoa
     452            0 :                               force_a(1) = force_a(1) + fac_a*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
     453            0 :                               force_a(2) = force_a(2) + fac_a*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
     454            0 :                               force_a(3) = force_a(3) + fac_a*pab(a_offset + i, b_offset + j)*tmpz((i - 1)*ncob + j)
     455              :                            END DO
     456              :                         END DO
     457              : 
     458            0 :                         tmpx(1:ncoa*ncob) = 0.0_dp
     459            0 :                         tmpy(1:ncoa*ncob) = 0.0_dp
     460            0 :                         tmpz(1:ncoa*ncob) = 0.0_dp
     461              : 
     462              :                         !force wrt to atomic position B
     463              :                         CALL libgrpp_type2_integrals_gradient(ra, li, 1, [normi], [zeti], &
     464              :                                                               rb, lj, 1, [normj], [zetj], &
     465              :                                                               rc, lk, [npot_ecp(lk)], nrpot_ecp(:, lk), &
     466              :                                                               coeffs_ecp(:, lk), alpha_ecp(:, lk), rb, &
     467            0 :                                                               tmpx, tmpy, tmpz)
     468              :                         !note: tmp array is in C row-major ordering
     469              :                         !note: zero-gradients sometime comes out as NaN, hence tampval==tmpval check
     470            0 :                         DO j = 1, ncob
     471            0 :                            DO i = 1, ncoa
     472            0 :                               force_b(1) = force_b(1) + fac_b*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
     473            0 :                               force_b(2) = force_b(2) + fac_b*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
     474            0 :                               force_b(3) = force_b(3) + fac_b*pab(a_offset + i, b_offset + j)*tmpz((i - 1)*ncob + j)
     475              :                            END DO
     476              :                         END DO
     477              : 
     478              :                      END IF !calc_forces
     479              : 
     480              :                   END DO !lk
     481              : 
     482              :                END DO !lj
     483              :             END DO !li
     484              : 
     485              :          END DO !jpgf
     486              :       END DO !ipgf
     487              : 
     488              : #else
     489              : 
     490              :       MARK_USED(la_max_set)
     491              :       MARK_USED(la_min_set)
     492              :       MARK_USED(npgfa)
     493              :       MARK_USED(rpgfa)
     494              :       MARK_USED(zeta)
     495              :       MARK_USED(lb_max_set)
     496              :       MARK_USED(lb_min_set)
     497              :       MARK_USED(npgfb)
     498              :       MARK_USED(rpgfb)
     499              :       MARK_USED(zetb)
     500              :       MARK_USED(lmax_ecp)
     501              :       MARK_USED(npot_ecp)
     502              :       MARK_USED(alpha_ecp)
     503              :       MARK_USED(coeffs_ecp)
     504              :       MARK_USED(nrpot_ecp)
     505              :       MARK_USED(rpgfc)
     506              :       MARK_USED(rab)
     507              :       MARK_USED(dab)
     508              :       MARK_USED(rac)
     509              :       MARK_USED(dac)
     510              :       MARK_USED(dbc)
     511              :       MARK_USED(vab)
     512              :       MARK_USED(pab)
     513              :       MARK_USED(force_a)
     514              :       MARK_USED(force_b)
     515              : 
     516              :       CPABORT("Please compile CP2K with libgrpp support for calculations with ECPs")
     517              : #endif
     518              : 
     519         9712 :    END SUBROUTINE libgrpp_semilocal_integrals
     520              : 
     521              : ! **************************************************************************************************
     522              : !> \brief Reference local ECP force routine using l+-1 integrals. No call is made to the numerically
     523              : !>        unstable gradient routine of libgrpp. Calculates both the integrals and the forces.
     524              : !> \param la_max_set ...
     525              : !> \param la_min_set ...
     526              : !> \param npgfa ...
     527              : !> \param rpgfa ...
     528              : !> \param zeta ...
     529              : !> \param lb_max_set ...
     530              : !> \param lb_min_set ...
     531              : !> \param npgfb ...
     532              : !> \param rpgfb ...
     533              : !> \param zetb ...
     534              : !> \param npot_ecp ...
     535              : !> \param alpha_ecp ...
     536              : !> \param coeffs_ecp ...
     537              : !> \param nrpot_ecp ...
     538              : !> \param rpgfc ...
     539              : !> \param rab ...
     540              : !> \param dab ...
     541              : !> \param rac ...
     542              : !> \param dac ...
     543              : !> \param dbc ...
     544              : !> \param vab ...
     545              : !> \param pab ...
     546              : !> \param force_a ...
     547              : !> \param force_b ...
     548              : !> \note: this is a reference routine, which has no reason to be used once the libgrpp gradients
     549              : !>        become numerically stable
     550              : ! **************************************************************************************************
     551            0 :    SUBROUTINE libgrpp_local_forces_ref(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
     552            0 :                                        lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
     553            0 :                                        npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
     554            0 :                                        rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
     555              : 
     556              :       INTEGER, INTENT(IN)                                :: la_max_set, la_min_set, npgfa
     557              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rpgfa, zeta
     558              :       INTEGER, INTENT(IN)                                :: lb_max_set, lb_min_set, npgfb
     559              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rpgfb, zetb
     560              :       INTEGER, INTENT(IN)                                :: npot_ecp
     561              :       REAL(KIND=dp), DIMENSION(1:npot_ecp), INTENT(IN)   :: alpha_ecp, coeffs_ecp
     562              :       INTEGER, DIMENSION(1:npot_ecp), INTENT(IN)         :: nrpot_ecp
     563              :       REAL(KIND=dp), INTENT(IN)                          :: rpgfc
     564              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     565              :       REAL(KIND=dp), INTENT(IN)                          :: dab
     566              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rac
     567              :       REAL(KIND=dp), INTENT(IN)                          :: dac
     568              :       REAL(KIND=dp), INTENT(IN)                          :: dbc
     569              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: vab
     570              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: pab
     571              :       REAL(KIND=dp), DIMENSION(3), INTENT(INOUT)         :: force_a, force_b
     572              : 
     573              : #if defined(__LIBGRPP)
     574              :       INTEGER                                            :: a_offset, a_start, b_offset, b_start, i, &
     575              :                                                             ipgf, j, jpgf, li, lj, ncoa, ncob, a_offset_f, &
     576              :                                                             b_offset_f, a_start_f, b_start_f
     577              :       REAL(dp)                                           :: expi, expj, normi, normj, prefi, prefj, &
     578              :                                                             zeti, zetj
     579            0 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: tmp
     580            0 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: vab_f, tmpx, tmpy, tmpz
     581              :       REAL(dp), DIMENSION(3)                             :: ra, rb, rc
     582              : 
     583            0 :       CALL libgrpp_init()
     584              : 
     585              :       !Contains the integrals necessary for the forces, with angular momenta from lmin-1 to lmax+1
     586            0 :       ALLOCATE (vab_f(npgfa*ncoset(la_max_set + 1), npgfb*ncoset(lb_max_set + 1)))
     587            0 :       vab_f(:, :) = 0.0_dp
     588              : 
     589              :       !libgrpp requires absolute positions, not relative ones
     590            0 :       ra(:) = 0.0_dp
     591            0 :       rb(:) = rab(:)
     592            0 :       rc(:) = rac(:)
     593              : 
     594            0 :       ALLOCATE (tmp(nco(la_max_set + 1)*nco(lb_max_set + 1)))
     595              : 
     596            0 :       DO ipgf = 1, npgfa
     597            0 :          IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
     598            0 :          zeti = zeta(ipgf)
     599            0 :          a_start = (ipgf - 1)*ncoset(la_max_set)
     600            0 :          a_start_f = (ipgf - 1)*ncoset(la_max_set + 1)
     601              : 
     602            0 :          DO jpgf = 1, npgfb
     603            0 :             IF (rpgfb(jpgf) + rpgfc < dbc) CYCLE
     604            0 :             IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) CYCLE
     605            0 :             zetj = zetb(jpgf)
     606            0 :             b_start = (jpgf - 1)*ncoset(lb_max_set)
     607            0 :             b_start_f = (jpgf - 1)*ncoset(lb_max_set + 1)
     608              : 
     609            0 :             DO li = MAX(0, la_min_set - 1), la_max_set + 1
     610            0 :                a_offset = a_start + ncoset(li - 1)
     611            0 :                a_offset_f = a_start_f + ncoset(li - 1)
     612            0 :                ncoa = nco(li)
     613            0 :                prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
     614            0 :                expi = 0.25_dp*REAL(2*li + 3, dp)
     615            0 :                normi = 1.0_dp/(prefi*zeti**expi)
     616              : 
     617            0 :                DO lj = MAX(0, lb_min_set - 1), lb_max_set + 1
     618            0 :                   b_offset = b_start + ncoset(lj - 1)
     619            0 :                   b_offset_f = b_start_f + ncoset(lj - 1)
     620            0 :                   ncob = nco(lj)
     621            0 :                   prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
     622            0 :                   expj = 0.25_dp*REAL(2*lj + 3, dp)
     623            0 :                   normj = 1.0_dp/(prefj*zetj**expj)
     624              : 
     625            0 :                   tmp(1:ncoa*ncob) = 0.0_dp
     626              :                   !libgrpp implicitely normalizes cartesian Gaussian. In CP2K, we do not, hence
     627              :                   !the 1/norm coefficients for PGFi and PGFj
     628              :                   CALL libgrpp_type1_integrals(ra, li, 1, [normi], [zeti], &
     629              :                                                rb, lj, 1, [normj], [zetj], &
     630              :                                                rc, [npot_ecp], nrpot_ecp, &
     631            0 :                                                coeffs_ecp, alpha_ecp, tmp)
     632              : 
     633              :                   !the l+-1 integrals for gradient calculation
     634            0 :                   DO j = 1, ncob
     635            0 :                      DO i = 1, ncoa
     636              :                         vab_f(a_offset_f + i, b_offset_f + j) = &
     637            0 :                            vab_f(a_offset_f + i, b_offset_f + j) + tmp((i - 1)*ncob + j)
     638              :                      END DO
     639              :                   END DO
     640              : 
     641              :                   !the actual integrals
     642            0 :                   IF (li >= la_min_set .AND. li <= la_max_set .AND. lj >= lb_min_set .AND. lj <= lb_max_set) THEN
     643            0 :                      DO j = 1, ncob
     644            0 :                         DO i = 1, ncoa
     645            0 :                            vab(a_offset + i, b_offset + j) = vab(a_offset + i, b_offset + j) + tmp((i - 1)*ncob + j)
     646              :                         END DO
     647              :                      END DO
     648              :                   END IF
     649              : 
     650              :                END DO !lj
     651              :             END DO !li
     652              : 
     653              :          END DO !jpgf
     654              :       END DO !ipgf
     655              : 
     656            0 :       ALLOCATE (tmpx(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
     657            0 :       ALLOCATE (tmpy(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
     658            0 :       ALLOCATE (tmpz(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
     659              : 
     660              :       !Derivative wrt to center A
     661            0 :       tmpx(:, :) = 0.0_dp
     662            0 :       tmpy(:, :) = 0.0_dp
     663            0 :       tmpz(:, :) = 0.0_dp
     664              :       CALL dabdr(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max_set, npgfb, rpgfb, lb_min_set, &
     665            0 :                  dab, vab_f, tmpx, tmpy, tmpz)
     666            0 :       DO j = 1, npgfb*ncoset(lb_max_set)
     667            0 :          DO i = 1, npgfa*ncoset(la_max_set)
     668            0 :             force_a(1) = force_a(1) + tmpx(i, j)*pab(i, j)
     669            0 :             force_a(2) = force_a(2) + tmpy(i, j)*pab(i, j)
     670            0 :             force_a(3) = force_a(3) + tmpz(i, j)*pab(i, j)
     671              :          END DO
     672              :       END DO
     673              : 
     674              :       !Derivative wrt to center B
     675            0 :       tmpx(:, :) = 0.0_dp
     676            0 :       tmpy(:, :) = 0.0_dp
     677            0 :       tmpz(:, :) = 0.0_dp
     678              :       CALL adbdr(la_max_set, npgfa, rpgfa, la_min_set, lb_max_set, npgfb, zetb, rpgfb, lb_min_set, &
     679            0 :                  dab, vab_f, tmpx, tmpy, tmpz)
     680            0 :       DO j = 1, npgfb*ncoset(lb_max_set)
     681            0 :          DO i = 1, npgfa*ncoset(la_max_set)
     682            0 :             force_b(1) = force_b(1) + tmpx(i, j)*pab(i, j)
     683            0 :             force_b(2) = force_b(2) + tmpy(i, j)*pab(i, j)
     684            0 :             force_b(3) = force_b(3) + tmpz(i, j)*pab(i, j)
     685              :          END DO
     686              :       END DO
     687            0 :       DEALLOCATE (tmpx, tmpy, tmpz)
     688              : 
     689              : #else
     690              : 
     691              :       MARK_USED(la_max_set)
     692              :       MARK_USED(la_min_set)
     693              :       MARK_USED(npgfa)
     694              :       MARK_USED(rpgfa)
     695              :       MARK_USED(zeta)
     696              :       MARK_USED(lb_max_set)
     697              :       MARK_USED(lb_min_set)
     698              :       MARK_USED(npgfb)
     699              :       MARK_USED(rpgfb)
     700              :       MARK_USED(zetb)
     701              :       MARK_USED(npot_ecp)
     702              :       MARK_USED(alpha_ecp)
     703              :       MARK_USED(coeffs_ecp)
     704              :       MARK_USED(nrpot_ecp)
     705              :       MARK_USED(rpgfc)
     706              :       MARK_USED(rab)
     707              :       MARK_USED(dab)
     708              :       MARK_USED(rac)
     709              :       MARK_USED(dac)
     710              :       MARK_USED(dbc)
     711              :       MARK_USED(pab)
     712              :       MARK_USED(vab)
     713              :       MARK_USED(force_a)
     714              :       MARK_USED(force_b)
     715              : 
     716              :       CPABORT("Please compile CP2K with libgrpp support for calculations with ECPs")
     717              : #endif
     718              : 
     719            0 :    END SUBROUTINE libgrpp_local_forces_ref
     720              : 
     721              : ! **************************************************************************************************
     722              : !> \brief Reference semi-local ECP forces using l+-1 integrals. No call is made to the numerically
     723              : !>        unstable gradient routine of libgrpp. Calculates both the integrals and the forces.
     724              : !> \param la_max_set ...
     725              : !> \param la_min_set ...
     726              : !> \param npgfa ...
     727              : !> \param rpgfa ...
     728              : !> \param zeta ...
     729              : !> \param lb_max_set ...
     730              : !> \param lb_min_set ...
     731              : !> \param npgfb ...
     732              : !> \param rpgfb ...
     733              : !> \param zetb ...
     734              : !> \param lmax_ecp ...
     735              : !> \param npot_ecp ...
     736              : !> \param alpha_ecp ...
     737              : !> \param coeffs_ecp ...
     738              : !> \param nrpot_ecp ...
     739              : !> \param rpgfc ...
     740              : !> \param rab ...
     741              : !> \param dab ...
     742              : !> \param rac ...
     743              : !> \param dac ...
     744              : !> \param dbc ...
     745              : !> \param vab ...
     746              : !> \param pab ...
     747              : !> \param force_a ...
     748              : !> \param force_b ...
     749              : !> \note: this is a reference routine, which has no reason to be used once the libgrpp gradients
     750              : !>        become numerically stable
     751              : ! **************************************************************************************************
     752         2888 :    SUBROUTINE libgrpp_semilocal_forces_ref(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
     753         2888 :                                            lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
     754              :                                            lmax_ecp, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
     755         2888 :                                            rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
     756              : 
     757              :       INTEGER, INTENT(IN)                                :: la_max_set, la_min_set, npgfa
     758              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rpgfa, zeta
     759              :       INTEGER, INTENT(IN)                                :: lb_max_set, lb_min_set, npgfb
     760              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rpgfb, zetb
     761              :       INTEGER, INTENT(IN)                                :: lmax_ecp
     762              :       INTEGER, DIMENSION(0:10), INTENT(IN)               :: npot_ecp
     763              :       REAL(KIND=dp), DIMENSION(1:15, 0:10), INTENT(IN)   :: alpha_ecp, coeffs_ecp
     764              :       INTEGER, DIMENSION(1:15, 0:10), INTENT(IN)         :: nrpot_ecp
     765              :       REAL(KIND=dp), INTENT(IN)                          :: rpgfc
     766              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     767              :       REAL(KIND=dp), INTENT(IN)                          :: dab
     768              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rac
     769              :       REAL(KIND=dp), INTENT(IN)                          :: dac
     770              :       REAL(KIND=dp), INTENT(IN)                          :: dbc
     771              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: vab
     772              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: pab
     773              :       REAL(KIND=dp), DIMENSION(3), INTENT(INOUT)         :: force_a, force_b
     774              : 
     775              : #if defined(__LIBGRPP)
     776              :       INTEGER                                            :: a_offset, a_start, b_offset, b_start, i, &
     777              :                                                             ipgf, j, jpgf, li, lj, lk, ncoa, ncob, &
     778              :                                                             a_start_f, b_start_f, a_offset_f, b_offset_f
     779              :       REAL(dp)                                           :: expi, expj, normi, normj, prefi, prefj, &
     780              :                                                             zeti, zetj
     781         2888 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: tmp
     782         2888 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: vab_f, tmpx, tmpy, tmpz
     783              :       REAL(dp), DIMENSION(3)                             :: ra, rb, rc
     784              : 
     785         2888 :       CALL libgrpp_init()
     786              : 
     787              :       !Contains the integrals necessary for the forces, with angular momenta from lmin-1 to lmax+1
     788        11552 :       ALLOCATE (vab_f(npgfa*ncoset(la_max_set + 1), npgfb*ncoset(lb_max_set + 1)))
     789       586688 :       vab_f(:, :) = 0.0_dp
     790              : 
     791              :       !libgrpp requires absolute positions, not relative ones
     792         2888 :       ra(:) = 0.0_dp
     793         2888 :       rb(:) = rab(:)
     794         2888 :       rc(:) = rac(:)
     795              : 
     796         8664 :       ALLOCATE (tmp(nco(la_max_set + 1)*nco(lb_max_set + 1)))
     797              : 
     798         7442 :       DO ipgf = 1, npgfa
     799         4554 :          IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
     800         3772 :          zeti = zeta(ipgf)
     801         3772 :          a_start = (ipgf - 1)*ncoset(la_max_set)
     802         3772 :          a_start_f = (ipgf - 1)*ncoset(la_max_set + 1)
     803              : 
     804        13006 :          DO jpgf = 1, npgfb
     805         6346 :             IF (rpgfb(jpgf) + rpgfc < dbc) CYCLE
     806         5412 :             IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) CYCLE
     807         4943 :             zetj = zetb(jpgf)
     808         4943 :             b_start = (jpgf - 1)*ncoset(lb_max_set)
     809         4943 :             b_start_f = (jpgf - 1)*ncoset(lb_max_set + 1)
     810              : 
     811        22112 :             DO li = MAX(0, la_min_set - 1), la_max_set + 1
     812        12615 :                a_offset = a_start + ncoset(li - 1)
     813        12615 :                a_offset_f = a_start_f + ncoset(li - 1)
     814        12615 :                ncoa = nco(li)
     815        12615 :                prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
     816        12615 :                expi = 0.25_dp*REAL(2*li + 3, dp)
     817        12615 :                normi = 1.0_dp/(prefi*zeti**expi)
     818              : 
     819        51515 :                DO lj = MAX(0, lb_min_set - 1), lb_max_set + 1
     820        32554 :                   b_offset = b_start + ncoset(lj - 1)
     821        32554 :                   b_offset_f = b_start_f + ncoset(lj - 1)
     822        32554 :                   ncob = nco(lj)
     823        32554 :                   prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
     824        32554 :                   expj = 0.25_dp*REAL(2*lj + 3, dp)
     825        32554 :                   normj = 1.0_dp/(prefj*zetj**expj)
     826              : 
     827              :                   !Loop over ECP angular momentum
     828       177901 :                   DO lk = 0, lmax_ecp
     829      1667965 :                      tmp(1:ncoa*ncob) = 0.0_dp
     830              :                      !libgrpp implicitely normalizes cartesian Gaussian. In CP2K, we do not, hence
     831              :                      !the 1/norm coefficients for PGFi and PGFj
     832              :                      CALL libgrpp_type2_integrals(ra, li, 1, [normi], [zeti], &
     833              :                                                   rb, lj, 1, [normj], [zetj], &
     834              :                                                   rc, lk, [npot_ecp(lk)], nrpot_ecp(:, lk), &
     835       796392 :                                                   coeffs_ecp(:, lk), alpha_ecp(:, lk), tmp)
     836              : 
     837              :                      !the l+-1 integrals for gradient calculation
     838       586239 :                      DO j = 1, ncob
     839      2121472 :                         DO i = 1, ncoa
     840              :                            vab_f(a_offset_f + i, b_offset_f + j) = &
     841      1988740 :                               vab_f(a_offset_f + i, b_offset_f + j) + tmp((i - 1)*ncob + j)
     842              :                         END DO
     843              :                      END DO
     844              : 
     845              :                      !the actual integrals
     846       165286 :                      IF (li >= la_min_set .AND. li <= la_max_set .AND. lj >= lb_min_set .AND. lj <= lb_max_set) THEN
     847        72268 :                         DO j = 1, ncob
     848       206965 :                            DO i = 1, ncoa
     849       186927 :                               vab(a_offset + i, b_offset + j) = vab(a_offset + i, b_offset + j) + tmp((i - 1)*ncob + j)
     850              :                            END DO
     851              :                         END DO
     852              :                      END IF
     853              : 
     854              :                   END DO !lk
     855              : 
     856              :                END DO !lj
     857              :             END DO !li
     858              : 
     859              :          END DO !jpgf
     860              :       END DO !ipgf
     861              : 
     862        11552 :       ALLOCATE (tmpx(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
     863         8664 :       ALLOCATE (tmpy(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
     864         8664 :       ALLOCATE (tmpz(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
     865              : 
     866              :       !Derivative wrt to center A
     867       104379 :       tmpx(:, :) = 0.0_dp
     868       104379 :       tmpy(:, :) = 0.0_dp
     869       104379 :       tmpz(:, :) = 0.0_dp
     870              :       CALL dabdr(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max_set, npgfb, rpgfb, lb_min_set, &
     871         2888 :                  0.0_dp, vab_f, tmpx, tmpy, tmpz)
     872        18655 :       DO j = 1, npgfb*ncoset(lb_max_set)
     873       104379 :          DO i = 1, npgfa*ncoset(la_max_set)
     874        85724 :             force_a(1) = force_a(1) + tmpx(i, j)*pab(i, j)
     875        85724 :             force_a(2) = force_a(2) + tmpy(i, j)*pab(i, j)
     876       101491 :             force_a(3) = force_a(3) + tmpz(i, j)*pab(i, j)
     877              :          END DO
     878              :       END DO
     879              : 
     880              :       !Derivative wrt to center B
     881       104379 :       tmpx(:, :) = 0.0_dp
     882       104379 :       tmpy(:, :) = 0.0_dp
     883       104379 :       tmpz(:, :) = 0.0_dp
     884              :       CALL adbdr(la_max_set, npgfa, rpgfa, la_min_set, lb_max_set, npgfb, zetb, rpgfb, lb_min_set, &
     885         2888 :                  0.0_dp, vab_f, tmpx, tmpy, tmpz)
     886        18655 :       DO j = 1, npgfb*ncoset(lb_max_set)
     887       104379 :          DO i = 1, npgfa*ncoset(la_max_set)
     888        85724 :             force_b(1) = force_b(1) + tmpx(i, j)*pab(i, j)
     889        85724 :             force_b(2) = force_b(2) + tmpy(i, j)*pab(i, j)
     890       101491 :             force_b(3) = force_b(3) + tmpz(i, j)*pab(i, j)
     891              :          END DO
     892              :       END DO
     893         2888 :       DEALLOCATE (tmpx, tmpy, tmpz)
     894              : 
     895              : #else
     896              : 
     897              :       MARK_USED(la_max_set)
     898              :       MARK_USED(la_min_set)
     899              :       MARK_USED(npgfa)
     900              :       MARK_USED(rpgfa)
     901              :       MARK_USED(zeta)
     902              :       MARK_USED(lb_max_set)
     903              :       MARK_USED(lb_min_set)
     904              :       MARK_USED(npgfb)
     905              :       MARK_USED(rpgfb)
     906              :       MARK_USED(zetb)
     907              :       MARK_USED(lmax_ecp)
     908              :       MARK_USED(npot_ecp)
     909              :       MARK_USED(alpha_ecp)
     910              :       MARK_USED(coeffs_ecp)
     911              :       MARK_USED(nrpot_ecp)
     912              :       MARK_USED(rpgfc)
     913              :       MARK_USED(rab)
     914              :       MARK_USED(dab)
     915              :       MARK_USED(rac)
     916              :       MARK_USED(dac)
     917              :       MARK_USED(dbc)
     918              :       MARK_USED(pab)
     919              :       MARK_USED(vab)
     920              :       MARK_USED(force_a)
     921              :       MARK_USED(force_b)
     922              : 
     923              :       CPABORT("Please compile CP2K with libgrpp support for calculations with ECPs")
     924              : #endif
     925              : 
     926         2888 :    END SUBROUTINE libgrpp_semilocal_forces_ref
     927              : 
     928              : END MODULE libgrpp_integrals
        

Generated by: LCOV version 2.0-1