LCOV - code coverage report
Current view: top level - src - semi_empirical_int_gks.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 54.0 % 681 368
Test Date: 2025-07-25 12:55:17 Functions: 66.7 % 12 8

            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 Integral GKS scheme: The order of the integrals in makeCoul reflects
      10              : !>        the standard order by MOPAC
      11              : !> \par History
      12              : !>      Teodoro Laino [tlaino] - 04.2009 : Adapted size arrays to d-orbitals and
      13              : !>                               get rid of the alternative ordering. Using the
      14              : !>                               CP2K one.
      15              : !>      Teodoro Laino [tlaino] - 04.2009 : Skip nullification (speed-up)
      16              : !>      Teodoro Laino [tlaino] - 04.2009 : Speed-up due to fortran arrays order
      17              : !>                               optimization and collection of common pieces of
      18              : !>                               code
      19              : ! **************************************************************************************************
      20              : MODULE semi_empirical_int_gks
      21              : 
      22              :    USE dg_rho0_types,                   ONLY: dg_rho0_type
      23              :    USE dg_types,                        ONLY: dg_get,&
      24              :                                               dg_type
      25              :    USE kinds,                           ONLY: dp
      26              :    USE mathconstants,                   ONLY: fourpi,&
      27              :                                               oorootpi
      28              :    USE multipole_types,                 ONLY: do_multipole_none
      29              :    USE pw_grid_types,                   ONLY: pw_grid_type
      30              :    USE pw_pool_types,                   ONLY: pw_pool_type
      31              :    USE semi_empirical_int_arrays,       ONLY: indexb,&
      32              :                                               rij_threshold
      33              :    USE semi_empirical_mpole_types,      ONLY: semi_empirical_mpole_type
      34              :    USE semi_empirical_types,            ONLY: se_int_control_type,&
      35              :                                               semi_empirical_type,&
      36              :                                               setup_se_int_control_type
      37              : #include "./base/base_uses.f90"
      38              : 
      39              :    IMPLICIT NONE
      40              :    PRIVATE
      41              : 
      42              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_int_gks'
      43              :    LOGICAL, PARAMETER, PRIVATE          :: debug_this_module = .FALSE.
      44              : 
      45              :    PUBLIC :: corecore_gks, rotnuc_gks, drotnuc_gks, rotint_gks, drotint_gks
      46              : 
      47              : CONTAINS
      48              : 
      49              : ! **************************************************************************************************
      50              : !> \brief Computes the electron-nuclei integrals
      51              : !> \param sepi ...
      52              : !> \param sepj ...
      53              : !> \param rij ...
      54              : !> \param e1b ...
      55              : !> \param e2a ...
      56              : !> \param se_int_control ...
      57              : ! **************************************************************************************************
      58         3813 :    SUBROUTINE rotnuc_gks(sepi, sepj, rij, e1b, e2a, se_int_control)
      59              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
      60              :       REAL(dp), DIMENSION(3), INTENT(IN)                 :: rij
      61              :       REAL(dp), DIMENSION(45), INTENT(OUT), OPTIONAL     :: e1b, e2a
      62              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
      63              : 
      64              :       INTEGER                                            :: i, mu, nu
      65              :       REAL(KIND=dp), DIMENSION(3)                        :: rab
      66              :       REAL(kind=dp), DIMENSION(45, 45)                   :: Coul
      67              : 
      68        15252 :       rab = -rij
      69              : 
      70         3813 :       IF (se_int_control%do_ewald_gks) THEN
      71           24 :          IF (DOT_PRODUCT(rij, rij) > rij_threshold) THEN
      72            3 :             CALL makeCoulE(rab, sepi, sepj, Coul, se_int_control)
      73              :          ELSE
      74            3 :             CALL makeCoulE0(sepi, Coul, se_int_control)
      75              :          END IF
      76              :       ELSE
      77         3807 :          CALL makeCoul(rab, sepi, sepj, Coul, se_int_control)
      78              :       END IF
      79              : 
      80         3813 :       i = 0
      81        19065 :       DO mu = 1, sepi%natorb
      82        57195 :          DO nu = 1, mu
      83        38130 :             i = i + 1
      84        53382 :             e1b(i) = -Coul(i, 1)*sepj%zeff
      85              :          END DO
      86              :       END DO
      87              : 
      88         3813 :       i = 0
      89        19065 :       DO mu = 1, sepj%natorb
      90        57195 :          DO nu = 1, mu
      91        38130 :             i = i + 1
      92        53382 :             e2a(i) = -Coul(1, i)*sepi%zeff
      93              :          END DO
      94              :       END DO
      95              : 
      96         3813 :    END SUBROUTINE rotnuc_gks
      97              : 
      98              : ! **************************************************************************************************
      99              : !> \brief Computes the electron-electron integrals
     100              : !> \param sepi ...
     101              : !> \param sepj ...
     102              : !> \param rij ...
     103              : !> \param w ...
     104              : !> \param se_int_control ...
     105              : ! **************************************************************************************************
     106         4754 :    SUBROUTINE rotint_gks(sepi, sepj, rij, w, se_int_control)
     107              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
     108              :       REAL(dp), DIMENSION(3), INTENT(IN)                 :: rij
     109              :       REAL(dp), DIMENSION(2025), INTENT(OUT), OPTIONAL   :: w
     110              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
     111              : 
     112              :       INTEGER                                            :: i, ind1, ind2, lam, mu, nu, sig
     113              :       REAL(KIND=dp), DIMENSION(3)                        :: rab
     114              :       REAL(kind=dp), DIMENSION(45, 45)                   :: Coul
     115              : 
     116        19016 :       rab = -rij
     117              : 
     118         4754 :       IF (se_int_control%do_ewald_gks) THEN
     119           24 :          IF (DOT_PRODUCT(rij, rij) > rij_threshold) THEN
     120            3 :             CALL makeCoulE(rab, sepi, sepj, Coul, se_int_control)
     121              :          ELSE
     122            3 :             CALL makeCoulE0(sepi, Coul, se_int_control)
     123              :          END IF
     124              :       ELSE
     125         4748 :          CALL makeCoul(rab, sepi, sepj, Coul, se_int_control)
     126              :       END IF
     127              : 
     128         4754 :       i = 0
     129         4754 :       ind1 = 0
     130        23770 :       DO mu = 1, sepi%natorb
     131        71310 :          DO nu = 1, mu
     132        47540 :             ind1 = ind1 + 1
     133        47540 :             ind2 = 0
     134       256716 :             DO lam = 1, sepj%natorb
     135       713100 :                DO sig = 1, lam
     136       475400 :                   i = i + 1
     137       475400 :                   ind2 = ind2 + 1
     138       665560 :                   w(i) = Coul(ind1, ind2)
     139              :                END DO
     140              :             END DO
     141              :          END DO
     142              :       END DO
     143              : 
     144         4754 :    END SUBROUTINE rotint_gks
     145              : 
     146              : ! **************************************************************************************************
     147              : !> \brief Computes the derivatives of the electron-nuclei integrals
     148              : !> \param sepi ...
     149              : !> \param sepj ...
     150              : !> \param rij ...
     151              : !> \param de1b ...
     152              : !> \param de2a ...
     153              : !> \param se_int_control ...
     154              : ! **************************************************************************************************
     155            0 :    SUBROUTINE drotnuc_gks(sepi, sepj, rij, de1b, de2a, se_int_control)
     156              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
     157              :       REAL(dp), DIMENSION(3), INTENT(IN)                 :: rij
     158              :       REAL(dp), DIMENSION(3, 45), INTENT(OUT), OPTIONAL  :: de1b, de2a
     159              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
     160              : 
     161              :       INTEGER                                            :: i, mu, nu
     162              :       REAL(KIND=dp), DIMENSION(3)                        :: rab
     163              :       REAL(kind=dp), DIMENSION(3, 45, 45)                :: dCoul
     164              : 
     165            0 :       rab = -rij
     166              : 
     167            0 :       IF (se_int_control%do_ewald_gks) THEN
     168            0 :          CALL makedCoulE(rab, sepi, sepj, dCoul, se_int_control)
     169              :       ELSE
     170            0 :          CALL makedCoul(rab, sepi, sepj, dCoul, se_int_control)
     171              :       END IF
     172              : 
     173            0 :       i = 0
     174            0 :       DO mu = 1, sepi%natorb
     175            0 :          DO nu = 1, mu
     176            0 :             i = i + 1
     177            0 :             de1b(1, i) = dCoul(1, i, 1)*sepj%zeff
     178            0 :             de1b(2, i) = dCoul(2, i, 1)*sepj%zeff
     179            0 :             de1b(3, i) = dCoul(3, i, 1)*sepj%zeff
     180              :          END DO
     181              :       END DO
     182              : 
     183            0 :       i = 0
     184            0 :       DO mu = 1, sepj%natorb
     185            0 :          DO nu = 1, mu
     186            0 :             i = i + 1
     187            0 :             de2a(1, i) = dCoul(1, 1, i)*sepi%zeff
     188            0 :             de2a(2, i) = dCoul(2, 1, i)*sepi%zeff
     189            0 :             de2a(3, i) = dCoul(3, 1, i)*sepi%zeff
     190              :          END DO
     191              :       END DO
     192              : 
     193            0 :    END SUBROUTINE drotnuc_gks
     194              : 
     195              : ! **************************************************************************************************
     196              : !> \brief Computes the derivatives of the electron-electron integrals
     197              : !> \param sepi ...
     198              : !> \param sepj ...
     199              : !> \param rij ...
     200              : !> \param dw ...
     201              : !> \param se_int_control ...
     202              : ! **************************************************************************************************
     203            0 :    SUBROUTINE drotint_gks(sepi, sepj, rij, dw, se_int_control)
     204              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
     205              :       REAL(dp), DIMENSION(3), INTENT(IN)                 :: rij
     206              :       REAL(dp), DIMENSION(3, 2025), INTENT(OUT)          :: dw
     207              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
     208              : 
     209              :       INTEGER                                            :: i, ind1, ind2, lam, mu, nu, sig
     210              :       REAL(KIND=dp), DIMENSION(3)                        :: rab
     211              :       REAL(kind=dp), DIMENSION(3, 45, 45)                :: dCoul
     212              : 
     213            0 :       rab = -rij
     214              : 
     215            0 :       IF (se_int_control%do_ewald_gks) THEN
     216            0 :          CALL makedCoulE(rab, sepi, sepj, dCoul, se_int_control)
     217              :       ELSE
     218            0 :          CALL makedCoul(rab, sepi, sepj, dCoul, se_int_control)
     219              :       END IF
     220              : 
     221            0 :       i = 0
     222            0 :       ind1 = 0
     223            0 :       DO mu = 1, sepi%natorb
     224            0 :          DO nu = 1, mu
     225            0 :             ind1 = ind1 + 1
     226            0 :             ind2 = 0
     227            0 :             DO lam = 1, sepj%natorb
     228            0 :                DO sig = 1, lam
     229            0 :                   i = i + 1
     230            0 :                   ind2 = ind2 + 1
     231            0 :                   dw(1, i) = -dCoul(1, ind1, ind2)
     232            0 :                   dw(2, i) = -dCoul(2, ind1, ind2)
     233            0 :                   dw(3, i) = -dCoul(3, ind1, ind2)
     234              :                END DO
     235              :             END DO
     236              :          END DO
     237              :       END DO
     238              : 
     239            0 :    END SUBROUTINE drotint_gks
     240              : 
     241              : ! **************************************************************************************************
     242              : !> \brief Computes the primitives of the integrals (non-periodic case)
     243              : !> \param RAB ...
     244              : !> \param sepi ...
     245              : !> \param sepj ...
     246              : !> \param Coul ...
     247              : !> \param se_int_control ...
     248              : ! **************************************************************************************************
     249        32338 :    SUBROUTINE makeCoul(RAB, sepi, sepj, Coul, se_int_control)
     250              :       REAL(kind=dp), DIMENSION(3)                        :: RAB
     251              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
     252              :       REAL(kind=dp), DIMENSION(45, 45), INTENT(OUT)      :: Coul
     253              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
     254              : 
     255              :       INTEGER                                            :: iA, iB, imA, imB, jA, jB, k1, k2, k3, k4
     256              :       LOGICAL                                            :: shortrange
     257              :       REAL(kind=dp)                                      :: a2, ACOULA, ACOULB, d1, d1f(3), d2, &
     258              :                                                             d2f(3, 3), d3, d3f(3, 3, 3), d4, &
     259              :                                                             d4f(3, 3, 3, 3), f, rr, w, w0, w1, w2, &
     260              :                                                             w3, w4, w5
     261              :       REAL(kind=dp), DIMENSION(3)                        :: v
     262              :       REAL(kind=dp), DIMENSION(3, 3, 45)                 :: M2A, M2B
     263              :       REAL(kind=dp), DIMENSION(3, 45)                    :: M1A, M1B
     264              :       REAL(kind=dp), DIMENSION(45)                       :: M0A, M0B
     265              : 
     266        16169 :       shortrange = se_int_control%shortrange
     267        16169 :       CALL get_se_slater_multipole(sepi, M0A, M1A, M2A, ACOULA)
     268        16169 :       CALL get_se_slater_multipole(sepj, M0B, M1B, M2B, ACOULB)
     269              : 
     270        16169 :       v(1) = RAB(1)
     271        16169 :       v(2) = RAB(2)
     272        16169 :       v(3) = RAB(3)
     273        64676 :       rr = SQRT(DOT_PRODUCT(v, v))
     274              : 
     275        16169 :       a2 = 0.5_dp*(1.0_dp/ACOULA + 1.0_dp/ACOULB)
     276        16169 :       w0 = a2*rr
     277        16169 :       w = EXP(-w0)
     278        16169 :       w1 = (1.0_dp + 0.5_dp*w0)
     279        16169 :       w2 = (w1 + 0.5_dp*w0 + 0.5_dp*w0**2)
     280        16169 :       w3 = (w2 + w0**3/6.0_dp)
     281        16169 :       w4 = (w3 + w0**4/30.0_dp)
     282        16169 :       w5 = (w3 + 8.0_dp*w0**4/210.0_dp + w0**5/210.0_dp)
     283              : 
     284        16169 :       IF (shortrange) THEN
     285        11412 :          f = (-w*w1)/rr
     286        11412 :          d1 = -1.0_dp*(-w*w2)/rr**3
     287        11412 :          d2 = 3.0_dp*(-w*w3)/rr**5
     288        11412 :          d3 = -15.0_dp*(-w*w4)/rr**7
     289        11412 :          d4 = 105.0_dp*(-w*w5)/rr**9
     290              :       ELSE
     291         4757 :          f = (1.0_dp - w*w1)/rr
     292         4757 :          d1 = -1.0_dp*(1.0_dp - w*w2)/rr**3
     293         4757 :          d2 = 3.0_dp*(1.0_dp - w*w3)/rr**5
     294         4757 :          d3 = -15.0_dp*(1.0_dp - w*w4)/rr**7
     295         4757 :          d4 = 105.0_dp*(1.0_dp - w*w5)/rr**9
     296              :       END IF
     297              : 
     298        16169 :       CALL build_d_tensor_gks(d1f, d2f, d3f, d4f, v=v, d1=d1, d2=d2, d3=d3, d4=d4)
     299              : 
     300        16169 :       imA = 0
     301        80845 :       DO iA = 1, sepi%natorb
     302       242535 :          DO jA = 1, iA
     303       161690 :             imA = imA + 1
     304              : 
     305       161690 :             imB = 0
     306       873126 :             DO iB = 1, sepj%natorb
     307      2425350 :                DO jB = 1, iB
     308      1616900 :                   imB = imB + 1
     309              : 
     310      1616900 :                   w = M0A(imA)*M0B(imB)*f
     311      6467600 :                   DO k1 = 1, 3
     312      6467600 :                      w = w + (M1A(k1, imA)*M0B(imB) - M0A(imA)*M1B(k1, imB))*d1f(k1)
     313              :                   END DO
     314      6467600 :                   DO k2 = 1, 3
     315     21019700 :                      DO k1 = 1, 3
     316     19402800 :                         w = w + (M2A(k1, k2, imA)*M0B(imB) - M1A(k1, imA)*M1B(k2, imB) + M0A(imA)*M2B(k1, k2, imB))*d2f(k1, k2)
     317              :                      END DO
     318              :                   END DO
     319      6467600 :                   DO k3 = 1, 3
     320     21019700 :                      DO k2 = 1, 3
     321     63059100 :                         DO k1 = 1, 3
     322     58208400 :                            w = w + (-M2A(k1, k2, imA)*M1B(k3, imB) + M1A(k1, imA)*M2B(k2, k3, imB))*d3f(k1, k2, k3)
     323              :                         END DO
     324              :                      END DO
     325              :                   END DO
     326              : 
     327      6467600 :                   DO k4 = 1, 3
     328     21019700 :                      DO k3 = 1, 3
     329     63059100 :                         DO k2 = 1, 3
     330    189177300 :                            DO k1 = 1, 3
     331    174625200 :                               w = w + M2A(k1, k2, imA)*M2B(k3, k4, imB)*d4f(k1, k2, k3, k4)
     332              :                            END DO
     333              :                         END DO
     334              :                      END DO
     335              :                   END DO
     336              : 
     337      2263660 :                   Coul(imA, imB) = w
     338              :                END DO
     339              :             END DO
     340              :          END DO
     341              :       END DO
     342              : 
     343        16169 :    END SUBROUTINE makeCoul
     344              : 
     345              : ! **************************************************************************************************
     346              : !> \brief Computes the derivatives of the primitives of the integrals
     347              : !>        (non-periodic case)
     348              : !> \param RAB ...
     349              : !> \param sepi ...
     350              : !> \param sepj ...
     351              : !> \param dCoul ...
     352              : !> \param se_int_control ...
     353              : ! **************************************************************************************************
     354            0 :    SUBROUTINE makedCoul(RAB, sepi, sepj, dCoul, se_int_control)
     355              :       REAL(kind=dp), DIMENSION(3)                        :: RAB
     356              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
     357              :       REAL(kind=dp), DIMENSION(3, 45, 45), INTENT(OUT)   :: dCoul
     358              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
     359              : 
     360              :       INTEGER                                            :: iA, iB, imA, imB, jA, jB, k1, k2, k3, k4
     361              :       LOGICAL                                            :: shortrange
     362              :       REAL(kind=dp) :: a2, ACOULA, ACOULB, d1, d1f(3), d2, d2f(3, 3), d3, d3f(3, 3, 3), d4, &
     363              :          d4f(3, 3, 3, 3), d5, d5f(3, 3, 3, 3, 3), f, rr, tmp, w, w0, w1, w2, w3, w4, w5, w6
     364              :       REAL(kind=dp), DIMENSION(3)                        :: v, wv
     365              :       REAL(kind=dp), DIMENSION(3, 3, 45)                 :: M2A, M2B
     366              :       REAL(kind=dp), DIMENSION(3, 45)                    :: M1A, M1B
     367              :       REAL(kind=dp), DIMENSION(45)                       :: M0A, M0B
     368              : 
     369            0 :       shortrange = se_int_control%shortrange
     370            0 :       CALL get_se_slater_multipole(sepi, M0A, M1A, M2A, ACOULA)
     371            0 :       CALL get_se_slater_multipole(sepj, M0B, M1B, M2B, ACOULB)
     372              : 
     373            0 :       v(1) = RAB(1)
     374            0 :       v(2) = RAB(2)
     375            0 :       v(3) = RAB(3)
     376            0 :       rr = SQRT(DOT_PRODUCT(v, v))
     377              : 
     378            0 :       a2 = 0.5_dp*(1.0_dp/ACOULA + 1.0_dp/ACOULB)
     379            0 :       w0 = a2*rr
     380            0 :       w = EXP(-w0)
     381            0 :       w1 = (1.0_dp + 0.5_dp*w0)
     382            0 :       w2 = (w1 + 0.5_dp*w0 + 0.5_dp*w0**2)
     383            0 :       w3 = (w2 + w0**3/6.0_dp)
     384            0 :       w4 = (w3 + w0**4/30.0_dp)
     385            0 :       w5 = (w3 + 4.0_dp*w0**4/105.0_dp + w0**5/210.0_dp)
     386            0 :       w6 = (w3 + 15.0_dp*w0**4/378.0_dp + 2.0_dp*w0**5/315.0_dp + w0**6/1890.0_dp)
     387              : 
     388            0 :       IF (shortrange) THEN
     389            0 :          f = (-w*w1)/rr
     390            0 :          d1 = -1.0_dp*(-w*w2)/rr**3
     391            0 :          d2 = 3.0_dp*(-w*w3)/rr**5
     392            0 :          d3 = -15.0_dp*(-w*w4)/rr**7
     393            0 :          d4 = 105.0_dp*(-w*w5)/rr**9
     394            0 :          d5 = -945.0_dp*(-w*w6)/rr**11
     395              :       ELSE
     396            0 :          f = (1.0_dp - w*w1)/rr
     397            0 :          d1 = -1.0_dp*(1.0_dp - w*w2)/rr**3
     398            0 :          d2 = 3.0_dp*(1.0_dp - w*w3)/rr**5
     399            0 :          d3 = -15.0_dp*(1.0_dp - w*w4)/rr**7
     400            0 :          d4 = 105.0_dp*(1.0_dp - w*w5)/rr**9
     401            0 :          d5 = -945.0_dp*(1.0_dp - w*w6)/rr**11
     402              :       END IF
     403              : 
     404            0 :       CALL build_d_tensor_gks(d1f, d2f, d3f, d4f, d5f, v, d1, d2, d3, d4, d5)
     405              : 
     406            0 :       imA = 0
     407            0 :       DO iA = 1, sepi%natorb
     408            0 :          DO jA = 1, iA
     409            0 :             imA = imA + 1
     410              : 
     411            0 :             imB = 0
     412            0 :             DO iB = 1, sepj%natorb
     413            0 :                DO jB = 1, iB
     414            0 :                   imB = imB + 1
     415              : 
     416            0 :                   tmp = M0A(imA)*M0B(imB)
     417            0 :                   wv(1) = tmp*d1f(1)
     418            0 :                   wv(2) = tmp*d1f(2)
     419            0 :                   wv(3) = tmp*d1f(3)
     420            0 :                   DO k1 = 1, 3
     421            0 :                      tmp = M1A(k1, imA)*M0B(imB) - M0A(imA)*M1B(k1, imB)
     422            0 :                      wv(1) = wv(1) + tmp*d2f(1, k1)
     423            0 :                      wv(2) = wv(2) + tmp*d2f(2, k1)
     424            0 :                      wv(3) = wv(3) + tmp*d2f(3, k1)
     425              :                   END DO
     426            0 :                   DO k2 = 1, 3
     427            0 :                      DO k1 = 1, 3
     428            0 :                         tmp = M2A(k1, k2, imA)*M0B(imB) - M1A(k1, imA)*M1B(k2, imB) + M0A(imA)*M2B(k1, k2, imB)
     429            0 :                         wv(1) = wv(1) + tmp*d3f(1, k1, k2)
     430            0 :                         wv(2) = wv(2) + tmp*d3f(2, k1, k2)
     431            0 :                         wv(3) = wv(3) + tmp*d3f(3, k1, k2)
     432              :                      END DO
     433              :                   END DO
     434            0 :                   DO k3 = 1, 3
     435            0 :                      DO k2 = 1, 3
     436            0 :                         DO k1 = 1, 3
     437            0 :                            tmp = -M2A(k1, k2, imA)*M1B(k3, imB) + M1A(k1, imA)*M2B(k2, k3, imB)
     438            0 :                            wv(1) = wv(1) + tmp*d4f(1, k1, k2, k3)
     439            0 :                            wv(2) = wv(2) + tmp*d4f(2, k1, k2, k3)
     440            0 :                            wv(3) = wv(3) + tmp*d4f(3, k1, k2, k3)
     441              :                         END DO
     442              :                      END DO
     443              :                   END DO
     444              : 
     445            0 :                   DO k4 = 1, 3
     446            0 :                      DO k3 = 1, 3
     447            0 :                         DO k2 = 1, 3
     448            0 :                            DO k1 = 1, 3
     449            0 :                               tmp = M2A(k1, k2, imA)*M2B(k3, k4, imB)
     450            0 :                               wv(1) = wv(1) + tmp*d5f(1, k1, k2, k3, k4)
     451            0 :                               wv(2) = wv(2) + tmp*d5f(2, k1, k2, k3, k4)
     452            0 :                               wv(3) = wv(3) + tmp*d5f(3, k1, k2, k3, k4)
     453              :                            END DO
     454              :                         END DO
     455              :                      END DO
     456              :                   END DO
     457              : 
     458            0 :                   dCoul(1, imA, imB) = wv(1)
     459            0 :                   dCoul(2, imA, imB) = wv(2)
     460            0 :                   dCoul(3, imA, imB) = wv(3)
     461              :                END DO
     462              :             END DO
     463              :          END DO
     464              :       END DO
     465              : 
     466            0 :    END SUBROUTINE makedCoul
     467              : 
     468              : ! **************************************************************************************************
     469              : !> \brief Computes nuclei-nuclei interactions
     470              : !> \param sepi ...
     471              : !> \param sepj ...
     472              : !> \param rijv ...
     473              : !> \param enuc ...
     474              : !> \param denuc ...
     475              : !> \param se_int_control ...
     476              : ! **************************************************************************************************
     477         3813 :    SUBROUTINE corecore_gks(sepi, sepj, rijv, enuc, denuc, se_int_control)
     478              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
     479              :       REAL(dp), DIMENSION(3), INTENT(IN)                 :: rijv
     480              :       REAL(dp), INTENT(OUT), OPTIONAL                    :: enuc
     481              :       REAL(dp), DIMENSION(3), INTENT(OUT), OPTIONAL      :: denuc
     482              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
     483              : 
     484              :       LOGICAL                                            :: l_denuc, l_enuc
     485              :       REAL(dp)                                           :: alpi, alpj, dscale, rij, scale, zz
     486              :       REAL(kind=dp), DIMENSION(3, 45, 45)                :: dCoul, dCoulE
     487              :       REAL(kind=dp), DIMENSION(45, 45)                   :: Coul, CoulE
     488              :       TYPE(se_int_control_type)                          :: se_int_control_off
     489              : 
     490        15252 :       rij = DOT_PRODUCT(rijv, rijv)
     491              : 
     492         3813 :       l_enuc = PRESENT(enuc)
     493         3813 :       l_denuc = PRESENT(denuc)
     494         3813 :       IF ((rij > rij_threshold) .AND. (l_enuc .OR. l_denuc)) THEN
     495              : 
     496         3810 :          rij = SQRT(rij)
     497              : 
     498         3810 :          IF (se_int_control%shortrange) THEN
     499              :             CALL setup_se_int_control_type(se_int_control_off, shortrange=.FALSE., do_ewald_r3=.FALSE., &
     500              :                                            do_ewald_gks=.FALSE., integral_screening=se_int_control%integral_screening, &
     501         3807 :                                            max_multipole=do_multipole_none, pc_coulomb_int=.FALSE.)
     502         3807 :             CALL makeCoul(rijv, sepi, sepj, Coul, se_int_control_off)
     503         3807 :             IF (l_denuc) CALL makedCoul(rijv, sepi, sepj, dCoul, se_int_control_off)
     504         3807 :             IF (se_int_control%do_ewald_gks) THEN
     505            3 :                CALL makeCoulE(rijv, sepi, sepj, CoulE, se_int_control)
     506            3 :                IF (l_denuc) CALL makedCoulE(rijv, sepi, sepj, dCoulE, se_int_control)
     507              :             ELSE
     508         3804 :                CALL makeCoul(rijv, sepi, sepj, CoulE, se_int_control)
     509         3804 :                IF (l_denuc) CALL makedCoul(rijv, sepi, sepj, dCoulE, se_int_control)
     510              :             END IF
     511              :          ELSE
     512            3 :             CALL makeCoul(rijv, sepi, sepj, Coul, se_int_control)
     513            3 :             CoulE = Coul
     514            3 :             IF (l_denuc) CALL makedCoul(rijv, sepi, sepj, dCoul, se_int_control)
     515            0 :             IF (l_denuc) dCoulE = dCoul
     516              :          END IF
     517              : 
     518         3810 :          scale = 0.0_dp
     519         3810 :          dscale = 0.0_dp
     520         3810 :          zz = sepi%zeff*sepj%zeff
     521         3810 :          alpi = sepi%alp
     522         3810 :          alpj = sepj%alp
     523         3810 :          scale = EXP(-alpi*rij) + EXP(-alpj*rij)
     524         3810 :          IF (l_enuc) THEN
     525         3810 :             enuc = zz*CoulE(1, 1) + scale*zz*Coul(1, 1)
     526              :          END IF
     527         3810 :          IF (l_denuc) THEN
     528            0 :             dscale = -alpi*EXP(-alpi*rij) - alpj*EXP(-alpj*rij)
     529            0 :             denuc(1) = zz*dCoulE(1, 1, 1) + dscale*(rijv(1)/rij)*zz*Coul(1, 1) + scale*zz*dCoul(1, 1, 1)
     530            0 :             denuc(2) = zz*dCoulE(2, 1, 1) + dscale*(rijv(2)/rij)*zz*Coul(1, 1) + scale*zz*dCoul(2, 1, 1)
     531            0 :             denuc(3) = zz*dCoulE(3, 1, 1) + dscale*(rijv(3)/rij)*zz*Coul(1, 1) + scale*zz*dCoul(3, 1, 1)
     532              :          END IF
     533              : 
     534              :       ELSE
     535              : 
     536            3 :          IF (se_int_control%do_ewald_gks) THEN
     537            3 :             zz = sepi%zeff*sepi%zeff
     538            3 :             CALL makeCoulE0(sepi, CoulE, se_int_control)
     539            3 :             IF (l_enuc) THEN
     540            3 :                enuc = zz*CoulE(1, 1)
     541              :             END IF
     542            3 :             IF (l_denuc) THEN
     543            0 :                denuc = 0._dp
     544              :             END IF
     545              :          END IF
     546              : 
     547              :       END IF
     548         3813 :    END SUBROUTINE corecore_gks
     549              : 
     550              : ! **************************************************************************************************
     551              : !> \brief Computes the primitives of the integrals (periodic case)
     552              : !> \param RAB ...
     553              : !> \param sepi ...
     554              : !> \param sepj ...
     555              : !> \param Coul ...
     556              : !> \param se_int_control ...
     557              : ! **************************************************************************************************
     558           27 :    SUBROUTINE makeCoulE(RAB, sepi, sepj, Coul, se_int_control)
     559              :       REAL(KIND=dp), DIMENSION(3)                        :: RAB
     560              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
     561              :       REAL(KIND=dp), DIMENSION(45, 45), INTENT(OUT)      :: Coul
     562              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
     563              : 
     564              :       INTEGER                                            :: gpt, imA, imB, k1, k2, k3, k4, lp, mp, np
     565            9 :       INTEGER, DIMENSION(:, :), POINTER                  :: bds
     566              :       REAL(KIND=dp) :: a2, ACOULA, ACOULB, alpha, cc, d1, d1f(3), d2, d2f(3, 3), d3, d3f(3, 3, 3), &
     567              :          d4, d4f(3, 3, 3, 3), f, ff, kr, kr2, r1, r2, r3, r5, r7, r9, rr, ss, w, w0, w1, w2, w3, &
     568              :          w4, w5
     569              :       REAL(KIND=dp), DIMENSION(3)                        :: kk, v
     570              :       REAL(KIND=dp), DIMENSION(3, 3, 45)                 :: M2A, M2B
     571              :       REAL(KIND=dp), DIMENSION(3, 45)                    :: M1A, M1B
     572              :       REAL(KIND=dp), DIMENSION(45)                       :: M0A, M0B
     573            9 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: rho0
     574              :       TYPE(dg_rho0_type), POINTER                        :: dg_rho0
     575              :       TYPE(dg_type), POINTER                             :: dg
     576              :       TYPE(pw_grid_type), POINTER                        :: pw_grid
     577              :       TYPE(pw_pool_type), POINTER                        :: pw_pool
     578              : 
     579            9 :       alpha = se_int_control%ewald_gks%alpha
     580            9 :       pw_pool => se_int_control%ewald_gks%pw_pool
     581            9 :       dg => se_int_control%ewald_gks%dg
     582            9 :       CALL dg_get(dg, dg_rho0=dg_rho0)
     583            9 :       rho0 => dg_rho0%density%array
     584            9 :       pw_grid => pw_pool%pw_grid
     585            9 :       bds => pw_grid%bounds
     586              : 
     587            9 :       CALL get_se_slater_multipole(sepi, M0A, M1A, M2A, ACOULA)
     588            9 :       CALL get_se_slater_multipole(sepj, M0B, M1B, M2B, ACOULB)
     589              : 
     590            9 :       v(1) = RAB(1)
     591            9 :       v(2) = RAB(2)
     592            9 :       v(3) = RAB(3)
     593           36 :       rr = SQRT(DOT_PRODUCT(v, v))
     594              : 
     595            9 :       r1 = 1.0_dp/rr
     596            9 :       r2 = r1*r1
     597            9 :       r3 = r2*r1
     598            9 :       r5 = r3*r2
     599            9 :       r7 = r5*r2
     600            9 :       r9 = r7*r2
     601              : 
     602            9 :       a2 = 0.5_dp*(1.0_dp/ACOULA + 1.0_dp/ACOULB)
     603              : 
     604            9 :       w0 = a2*rr
     605            9 :       w = EXP(-w0)
     606            9 :       w1 = (1.0_dp + 0.5_dp*w0)
     607            9 :       w2 = (w1 + 0.5_dp*w0 + 0.5_dp*w0**2)
     608            9 :       w3 = (w2 + w0**3/6.0_dp)
     609            9 :       w4 = (w3 + w0**4/30.0_dp)
     610            9 :       w5 = (w3 + 8.0_dp*w0**4/210.0_dp + w0**5/210.0_dp)
     611              : 
     612            9 :       f = (1.0_dp - w*w1)*r1
     613            9 :       d1 = -1.0_dp*(1.0_dp - w*w2)*r3
     614            9 :       d2 = 3.0_dp*(1.0_dp - w*w3)*r5
     615            9 :       d3 = -15.0_dp*(1.0_dp - w*w4)*r7
     616            9 :       d4 = 105.0_dp*(1.0_dp - w*w5)*r9
     617              : 
     618            9 :       kr = alpha*rr
     619            9 :       kr2 = kr*kr
     620            9 :       w0 = 1.0_dp - erfc(kr)
     621            9 :       w1 = 2.0_dp*oorootpi*EXP(-kr2)
     622            9 :       w2 = w1*kr
     623              : 
     624            9 :       f = f - w0*r1
     625            9 :       d1 = d1 + (-w2 + w0)*r3
     626            9 :       d2 = d2 + (w2*(3.0_dp + kr2*2.0_dp) - 3.0_dp*w0)*r5
     627            9 :       d3 = d3 + (-w2*(15.0_dp + kr2*(10.0_dp + kr2*4.0_dp)) + 15.0_dp*w0)*r7
     628            9 :       d4 = d4 + (w2*(105.0_dp + kr2*(70.0_dp + kr2*(28.0_dp + kr2*8.0_dp))) - 105.0_dp*w0)*r9
     629              : 
     630            9 :       CALL build_d_tensor_gks(d1f, d2f, d3f, d4f, v=v, d1=d1, d2=d2, d3=d3, d4=d4)
     631              : 
     632           99 :       DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2
     633          999 :          DO imB = 1, (sepj%natorb*(sepj%natorb + 1))/2
     634              : 
     635          900 :             w = M0A(imA)*M0B(imB)*f
     636         3600 :             DO k1 = 1, 3
     637         3600 :                w = w + (M1A(k1, imA)*M0B(imB) - M0A(imA)*M1B(k1, imB))*d1f(k1)
     638              :             END DO
     639         3600 :             DO k2 = 1, 3
     640        11700 :                DO k1 = 1, 3
     641        10800 :                   w = w + (M2A(k1, k2, imA)*M0B(imB) - M1A(k1, imA)*M1B(k2, imB) + M0A(imA)*M2B(k1, k2, imB))*d2f(k1, k2)
     642              :                END DO
     643              :             END DO
     644         3600 :             DO k3 = 1, 3
     645        11700 :                DO k2 = 1, 3
     646        35100 :                   DO k1 = 1, 3
     647        32400 :                      w = w + (-M2A(k1, k2, imA)*M1B(k3, imB) + M1A(k1, imA)*M2B(k2, k3, imB))*d3f(k1, k2, k3)
     648              :                   END DO
     649              :                END DO
     650              :             END DO
     651              : 
     652         3600 :             DO k4 = 1, 3
     653        11700 :                DO k3 = 1, 3
     654        35100 :                   DO k2 = 1, 3
     655       105300 :                      DO k1 = 1, 3
     656        97200 :                         w = w + M2A(k1, k2, imA)*M2B(k3, k4, imB)*d4f(k1, k2, k3, k4)
     657              :                      END DO
     658              :                   END DO
     659              :                END DO
     660              :             END DO
     661              : 
     662          990 :             Coul(imA, imB) = w
     663              :          END DO
     664              :       END DO
     665              : 
     666              :       v(1) = RAB(1)
     667              :       v(2) = RAB(2)
     668              :       v(3) = RAB(3)
     669              : 
     670            9 :       f = 0.0_dp
     671            9 :       d1f = 0.0_dp
     672            9 :       d2f = 0.0_dp
     673            9 :       d3f = 0.0_dp
     674            9 :       d4f = 0.0_dp
     675              : 
     676       150183 :       DO gpt = 1, pw_grid%ngpts_cut_local
     677       150174 :          lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
     678       150174 :          mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
     679       150174 :          np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
     680              : 
     681       150174 :          lp = lp + bds(1, 1)
     682       150174 :          mp = mp + bds(1, 2)
     683       150174 :          np = np + bds(1, 3)
     684              : 
     685       150174 :          IF (pw_grid%gsq(gpt) == 0.0_dp) CYCLE
     686       600660 :          kk(:) = pw_grid%g(:, gpt)
     687       150165 :          ff = 2.0_dp*fourpi*rho0(lp, mp, np)**2*pw_grid%vol/pw_grid%gsq(gpt)
     688              : 
     689       600660 :          kr = DOT_PRODUCT(kk, v)
     690       150165 :          cc = COS(kr)
     691       150165 :          ss = SIN(kr)
     692              : 
     693       150165 :          f = f + cc*ff
     694       600660 :          DO k1 = 1, 3
     695       600660 :             d1f(k1) = d1f(k1) - kk(k1)*ss*ff
     696              :          END DO
     697       600660 :          DO k2 = 1, 3
     698      1952145 :             DO k1 = 1, 3
     699      1801980 :                d2f(k1, k2) = d2f(k1, k2) - kk(k1)*kk(k2)*cc*ff
     700              :             END DO
     701              :          END DO
     702       600660 :          DO k3 = 1, 3
     703      1952145 :             DO k2 = 1, 3
     704      5856435 :                DO k1 = 1, 3
     705      5405940 :                   d3f(k1, k2, k3) = d3f(k1, k2, k3) + kk(k1)*kk(k2)*kk(k3)*ss*ff
     706              :                END DO
     707              :             END DO
     708              :          END DO
     709       600669 :          DO k4 = 1, 3
     710      1952154 :             DO k3 = 1, 3
     711      5856435 :                DO k2 = 1, 3
     712     17569305 :                   DO k1 = 1, 3
     713     16217820 :                      d4f(k1, k2, k3, k4) = d4f(k1, k2, k3, k4) + kk(k1)*kk(k2)*kk(k3)*kk(k4)*cc*ff
     714              :                   END DO
     715              :                END DO
     716              :             END DO
     717              :          END DO
     718              : 
     719              :       END DO
     720              : 
     721           99 :       DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2
     722          999 :          DO imB = 1, (sepj%natorb*(sepj%natorb + 1))/2
     723              : 
     724          900 :             w = M0A(imA)*M0B(imB)*f
     725         3600 :             DO k1 = 1, 3
     726         3600 :                w = w + (M1A(k1, imA)*M0B(imB) - M0A(imA)*M1B(k1, imB))*d1f(k1)
     727              :             END DO
     728         3600 :             DO k2 = 1, 3
     729        11700 :                DO k1 = 1, 3
     730        10800 :                   w = w + (M2A(k1, k2, imA)*M0B(imB) - M1A(k1, imA)*M1B(k2, imB) + M0A(imA)*M2B(k1, k2, imB))*d2f(k1, k2)
     731              :                END DO
     732              :             END DO
     733         3600 :             DO k3 = 1, 3
     734        11700 :                DO k2 = 1, 3
     735        35100 :                   DO k1 = 1, 3
     736        32400 :                      w = w + (-M2A(k1, k2, imA)*M1B(k3, imB) + M1A(k1, imA)*M2B(k2, k3, imB))*d3f(k1, k2, k3)
     737              :                   END DO
     738              :                END DO
     739              :             END DO
     740              : 
     741         3600 :             DO k4 = 1, 3
     742        11700 :                DO k3 = 1, 3
     743        35100 :                   DO k2 = 1, 3
     744       105300 :                      DO k1 = 1, 3
     745        97200 :                         w = w + M2A(k1, k2, imA)*M2B(k3, k4, imB)*d4f(k1, k2, k3, k4)
     746              :                      END DO
     747              :                   END DO
     748              :                END DO
     749              :             END DO
     750              : 
     751          990 :             Coul(imA, imB) = Coul(imA, imB) + w
     752              : 
     753              :          END DO
     754              :       END DO
     755              : 
     756           99 :       DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2
     757          999 :          DO imB = 1, (sepj%natorb*(sepj%natorb + 1))/2
     758          900 :             w = -M0A(imA)*M0B(imB)*0.25_dp*fourpi/(pw_grid%vol*alpha**2)
     759          990 :             Coul(imA, imB) = Coul(imA, imB) + w
     760              :          END DO
     761              :       END DO
     762              : 
     763            9 :    END SUBROUTINE makeCoulE
     764              : 
     765              : ! **************************************************************************************************
     766              : !> \brief Computes the derivatives of the primitives of the integrals
     767              : !>        (periodic case)
     768              : !> \param RAB ...
     769              : !> \param sepi ...
     770              : !> \param sepj ...
     771              : !> \param dCoul ...
     772              : !> \param se_int_control ...
     773              : ! **************************************************************************************************
     774            0 :    SUBROUTINE makedCoulE(RAB, sepi, sepj, dCoul, se_int_control)
     775              :       REAL(KIND=dp), DIMENSION(3)                        :: RAB
     776              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
     777              :       REAL(KIND=dp), DIMENSION(3, 45, 45), INTENT(OUT)   :: dCoul
     778              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
     779              : 
     780              :       INTEGER                                            :: gpt, imA, imB, k1, k2, k3, k4, k5, lp, &
     781              :                                                             mp, np
     782            0 :       INTEGER, DIMENSION(:, :), POINTER                  :: bds
     783              :       REAL(KIND=dp) :: a2, ACOULA, ACOULB, alpha, cc, d1, d1f(3), d2, d2f(3, 3), d3, d3f(3, 3, 3), &
     784              :          d4, d4f(3, 3, 3, 3), d5, d5f(3, 3, 3, 3, 3), f, ff, kr, kr2, r1, r11, r2, r3, r5, r7, r9, &
     785              :          rr, ss, tmp, w, w0, w1, w2, w3, w4, w5, w6
     786              :       REAL(KIND=dp), DIMENSION(3)                        :: kk, v, wv
     787              :       REAL(kind=dp), DIMENSION(3, 3, 45)                 :: M2A, M2B
     788              :       REAL(kind=dp), DIMENSION(3, 45)                    :: M1A, M1B
     789              :       REAL(kind=dp), DIMENSION(45)                       :: M0A, M0B
     790            0 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: rho0
     791              :       TYPE(dg_rho0_type), POINTER                        :: dg_rho0
     792              :       TYPE(dg_type), POINTER                             :: dg
     793              :       TYPE(pw_grid_type), POINTER                        :: pw_grid
     794              :       TYPE(pw_pool_type), POINTER                        :: pw_pool
     795              : 
     796            0 :       alpha = se_int_control%ewald_gks%alpha
     797            0 :       pw_pool => se_int_control%ewald_gks%pw_pool
     798            0 :       dg => se_int_control%ewald_gks%dg
     799            0 :       CALL dg_get(dg, dg_rho0=dg_rho0)
     800            0 :       rho0 => dg_rho0%density%array
     801            0 :       pw_grid => pw_pool%pw_grid
     802            0 :       bds => pw_grid%bounds
     803              : 
     804            0 :       CALL get_se_slater_multipole(sepi, M0A, M1A, M2A, ACOULA)
     805            0 :       CALL get_se_slater_multipole(sepj, M0B, M1B, M2B, ACOULB)
     806              : 
     807            0 :       v(1) = RAB(1)
     808            0 :       v(2) = RAB(2)
     809            0 :       v(3) = RAB(3)
     810            0 :       rr = SQRT(DOT_PRODUCT(v, v))
     811              : 
     812            0 :       a2 = 0.5_dp*(1.0_dp/ACOULA + 1.0_dp/ACOULB)
     813              : 
     814            0 :       r1 = 1.0_dp/rr
     815            0 :       r2 = r1*r1
     816            0 :       r3 = r2*r1
     817            0 :       r5 = r3*r2
     818            0 :       r7 = r5*r2
     819            0 :       r9 = r7*r2
     820            0 :       r11 = r9*r2
     821              : 
     822            0 :       w0 = a2*rr
     823            0 :       w = EXP(-w0)
     824            0 :       w1 = (1.0_dp + 0.5_dp*w0)
     825            0 :       w2 = (w1 + 0.5_dp*w0 + 0.5_dp*w0**2)
     826            0 :       w3 = (w2 + w0**3/6.0_dp)
     827            0 :       w4 = (w3 + w0**4/30.0_dp)
     828            0 :       w5 = (w3 + 8.0_dp*w0**4/210.0_dp + w0**5/210.0_dp)
     829            0 :       w6 = (w3 + 5.0_dp*w0**4/126.0_dp + 2.0_dp*w0**5/315.0_dp + w0**6/1890.0_dp)
     830              : 
     831            0 :       f = (1.0_dp - w*w1)*r1
     832            0 :       d1 = -1.0_dp*(1.0_dp - w*w2)*r3
     833            0 :       d2 = 3.0_dp*(1.0_dp - w*w3)*r5
     834            0 :       d3 = -15.0_dp*(1.0_dp - w*w4)*r7
     835            0 :       d4 = 105.0_dp*(1.0_dp - w*w5)*r9
     836            0 :       d5 = -945.0_dp*(1.0_dp - w*w6)*r11
     837              : 
     838            0 :       kr = alpha*rr
     839            0 :       kr2 = kr*kr
     840            0 :       w0 = 1.0_dp - erfc(kr)
     841            0 :       w1 = 2.0_dp*oorootpi*EXP(-kr2)
     842            0 :       w2 = w1*kr
     843              : 
     844            0 :       f = f - w0*r1
     845            0 :       d1 = d1 + (-w2 + w0)*r3
     846            0 :       d2 = d2 + (w2*(3.0_dp + kr2*2.0_dp) - 3.0_dp*w0)*r5
     847            0 :       d3 = d3 + (-w2*(15.0_dp + kr2*(10.0_dp + kr2*4.0_dp)) + 15.0_dp*w0)*r7
     848            0 :       d4 = d4 + (w2*(105.0_dp + kr2*(70.0_dp + kr2*(28.0_dp + kr2*8.0_dp))) - 105.0_dp*w0)*r9
     849            0 :       d5 = d5 + (-w2*(945.0_dp + kr2*(630.0_dp + kr2*(252.0_dp + kr2*(72.0_dp + kr2*16.0_dp)))) + 945.0_dp*w0)*r11
     850              : 
     851            0 :       CALL build_d_tensor_gks(d1f, d2f, d3f, d4f, d5f, v, d1, d2, d3, d4, d5)
     852              : 
     853            0 :       DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2
     854            0 :          DO imB = 1, (sepj%natorb*(sepj%natorb + 1))/2
     855              : 
     856            0 :             tmp = M0A(imA)*M0B(imB)
     857            0 :             wv(1) = tmp*d1f(1)
     858            0 :             wv(2) = tmp*d1f(2)
     859            0 :             wv(3) = tmp*d1f(3)
     860              : 
     861            0 :             DO k1 = 1, 3
     862            0 :                tmp = M1A(k1, imA)*M0B(imB) - M0A(imA)*M1B(k1, imB)
     863            0 :                wv(1) = wv(1) + tmp*d2f(1, k1)
     864            0 :                wv(2) = wv(2) + tmp*d2f(2, k1)
     865            0 :                wv(3) = wv(3) + tmp*d2f(3, k1)
     866              :             END DO
     867            0 :             DO k2 = 1, 3
     868            0 :                DO k1 = 1, 3
     869            0 :                   tmp = M2A(k1, k2, imA)*M0B(imB) - M1A(k1, imA)*M1B(k2, imB) + M0A(imA)*M2B(k1, k2, imB)
     870            0 :                   wv(1) = wv(1) + tmp*d3f(1, k1, k2)
     871            0 :                   wv(2) = wv(2) + tmp*d3f(2, k1, k2)
     872            0 :                   wv(3) = wv(3) + tmp*d3f(3, k1, k2)
     873              :                END DO
     874              :             END DO
     875            0 :             DO k3 = 1, 3
     876            0 :                DO k2 = 1, 3
     877            0 :                   DO k1 = 1, 3
     878            0 :                      tmp = -M2A(k1, k2, imA)*M1B(k3, imB) + M1A(k1, imA)*M2B(k2, k3, imB)
     879            0 :                      wv(1) = wv(1) + tmp*d4f(1, k1, k2, k3)
     880            0 :                      wv(2) = wv(2) + tmp*d4f(2, k1, k2, k3)
     881            0 :                      wv(3) = wv(3) + tmp*d4f(3, k1, k2, k3)
     882              :                   END DO
     883              :                END DO
     884              :             END DO
     885              : 
     886            0 :             DO k4 = 1, 3
     887            0 :                DO k3 = 1, 3
     888            0 :                   DO k2 = 1, 3
     889            0 :                      DO k1 = 1, 3
     890            0 :                         tmp = M2A(k1, k2, imA)*M2B(k3, k4, imB)
     891            0 :                         wv(1) = wv(1) + tmp*d5f(1, k1, k2, k3, k4)
     892            0 :                         wv(2) = wv(2) + tmp*d5f(2, k1, k2, k3, k4)
     893            0 :                         wv(3) = wv(3) + tmp*d5f(3, k1, k2, k3, k4)
     894              :                      END DO
     895              :                   END DO
     896              :                END DO
     897              :             END DO
     898              : 
     899            0 :             dCoul(1, imA, imB) = wv(1)
     900            0 :             dCoul(2, imA, imB) = wv(2)
     901            0 :             dCoul(3, imA, imB) = wv(3)
     902              :          END DO
     903              :       END DO
     904              : 
     905              :       v(1) = RAB(1)
     906              :       v(2) = RAB(2)
     907              :       v(3) = RAB(3)
     908              : 
     909            0 :       f = 0.0_dp
     910            0 :       d1f = 0.0_dp
     911            0 :       d2f = 0.0_dp
     912            0 :       d3f = 0.0_dp
     913            0 :       d4f = 0.0_dp
     914            0 :       d5f = 0.0_dp
     915              : 
     916            0 :       DO gpt = 1, pw_grid%ngpts_cut_local
     917            0 :          lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
     918            0 :          mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
     919            0 :          np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
     920              : 
     921            0 :          lp = lp + bds(1, 1)
     922            0 :          mp = mp + bds(1, 2)
     923            0 :          np = np + bds(1, 3)
     924              : 
     925            0 :          IF (pw_grid%gsq(gpt) == 0.0_dp) CYCLE
     926            0 :          kk(:) = pw_grid%g(:, gpt)
     927            0 :          ff = 2.0_dp*fourpi*rho0(lp, mp, np)**2*pw_grid%vol/pw_grid%gsq(gpt)
     928              : 
     929            0 :          kr = DOT_PRODUCT(kk, v)
     930            0 :          cc = COS(kr)
     931            0 :          ss = SIN(kr)
     932              : 
     933            0 :          f = f + cc*ff
     934            0 :          DO k1 = 1, 3
     935            0 :             d1f(k1) = d1f(k1) - kk(k1)*ss*ff
     936              :          END DO
     937            0 :          DO k2 = 1, 3
     938            0 :             DO k1 = 1, 3
     939            0 :                d2f(k1, k2) = d2f(k1, k2) - kk(k1)*kk(k2)*cc*ff
     940              :             END DO
     941              :          END DO
     942            0 :          DO k3 = 1, 3
     943            0 :             DO k2 = 1, 3
     944            0 :                DO k1 = 1, 3
     945            0 :                   d3f(k1, k2, k3) = d3f(k1, k2, k3) + kk(k1)*kk(k2)*kk(k3)*ss*ff
     946              :                END DO
     947              :             END DO
     948              :          END DO
     949            0 :          DO k4 = 1, 3
     950            0 :             DO k3 = 1, 3
     951            0 :                DO k2 = 1, 3
     952            0 :                   DO k1 = 1, 3
     953            0 :                      d4f(k1, k2, k3, k4) = d4f(k1, k2, k3, k4) + kk(k1)*kk(k2)*kk(k3)*kk(k4)*cc*ff
     954              :                   END DO
     955              :                END DO
     956              :             END DO
     957              :          END DO
     958            0 :          DO k5 = 1, 3
     959            0 :             DO k4 = 1, 3
     960            0 :                DO k3 = 1, 3
     961            0 :                   DO k2 = 1, 3
     962            0 :                      DO k1 = 1, 3
     963            0 :                         d5f(k1, k2, k3, k4, k5) = d5f(k1, k2, k3, k4, k5) - kk(k1)*kk(k2)*kk(k3)*kk(k4)*kk(k5)*ss*ff
     964              :                      END DO
     965              :                   END DO
     966              :                END DO
     967              :             END DO
     968              :          END DO
     969              :       END DO
     970              : 
     971            0 :       DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2
     972            0 :          DO imB = 1, (sepj%natorb*(sepj%natorb + 1))/2
     973            0 :             tmp = M0A(imA)*M0B(imB)
     974            0 :             wv(1) = tmp*d1f(1)
     975            0 :             wv(2) = tmp*d1f(2)
     976            0 :             wv(3) = tmp*d1f(3)
     977            0 :             DO k1 = 1, 3
     978            0 :                tmp = M1A(k1, imA)*M0B(imB) - M0A(imA)*M1B(k1, imB)
     979            0 :                wv(1) = wv(1) + tmp*d2f(1, k1)
     980            0 :                wv(2) = wv(2) + tmp*d2f(2, k1)
     981            0 :                wv(3) = wv(3) + tmp*d2f(3, k1)
     982              :             END DO
     983            0 :             DO k2 = 1, 3
     984            0 :                DO k1 = 1, 3
     985            0 :                   tmp = M2A(k1, k2, imA)*M0B(imB) - M1A(k1, imA)*M1B(k2, imB) + M0A(imA)*M2B(k1, k2, imB)
     986            0 :                   wv(1) = wv(1) + tmp*d3f(1, k1, k2)
     987            0 :                   wv(2) = wv(2) + tmp*d3f(2, k1, k2)
     988            0 :                   wv(3) = wv(3) + tmp*d3f(3, k1, k2)
     989              :                END DO
     990              :             END DO
     991            0 :             DO k3 = 1, 3
     992            0 :                DO k2 = 1, 3
     993            0 :                   DO k1 = 1, 3
     994            0 :                      tmp = -M2A(k1, k2, imA)*M1B(k3, imB) + M1A(k1, imA)*M2B(k2, k3, imB)
     995            0 :                      wv(1) = wv(1) + tmp*d4f(1, k1, k2, k3)
     996            0 :                      wv(2) = wv(2) + tmp*d4f(2, k1, k2, k3)
     997            0 :                      wv(3) = wv(3) + tmp*d4f(3, k1, k2, k3)
     998              :                   END DO
     999              :                END DO
    1000              :             END DO
    1001              : 
    1002            0 :             DO k4 = 1, 3
    1003            0 :                DO k3 = 1, 3
    1004            0 :                   DO k2 = 1, 3
    1005            0 :                      DO k1 = 1, 3
    1006            0 :                         tmp = M2A(k1, k2, imA)*M2B(k3, k4, imB)
    1007            0 :                         wv(1) = wv(1) + tmp*d5f(1, k1, k2, k3, k4)
    1008            0 :                         wv(2) = wv(2) + tmp*d5f(2, k1, k2, k3, k4)
    1009            0 :                         wv(3) = wv(3) + tmp*d5f(3, k1, k2, k3, k4)
    1010              :                      END DO
    1011              :                   END DO
    1012              :                END DO
    1013              :             END DO
    1014              : 
    1015            0 :             dCoul(1, imA, imB) = dCoul(1, imA, imB) + wv(1)
    1016            0 :             dCoul(2, imA, imB) = dCoul(2, imA, imB) + wv(2)
    1017            0 :             dCoul(3, imA, imB) = dCoul(3, imA, imB) + wv(3)
    1018              :          END DO
    1019              :       END DO
    1020              : 
    1021            0 :    END SUBROUTINE makedCoulE
    1022              : 
    1023              : ! **************************************************************************************************
    1024              : !> \brief Builds the tensor for the evaluation of the integrals with the
    1025              : !>        cartesian multipoles
    1026              : !> \param d1f ...
    1027              : !> \param d2f ...
    1028              : !> \param d3f ...
    1029              : !> \param d4f ...
    1030              : !> \param d5f ...
    1031              : !> \param v ...
    1032              : !> \param d1 ...
    1033              : !> \param d2 ...
    1034              : !> \param d3 ...
    1035              : !> \param d4 ...
    1036              : !> \param d5 ...
    1037              : ! **************************************************************************************************
    1038        16178 :    SUBROUTINE build_d_tensor_gks(d1f, d2f, d3f, d4f, d5f, v, d1, d2, d3, d4, d5)
    1039              :       REAL(KIND=dp), DIMENSION(3), INTENT(OUT)           :: d1f
    1040              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT)        :: d2f
    1041              :       REAL(KIND=dp), DIMENSION(3, 3, 3), INTENT(OUT)     :: d3f
    1042              :       REAL(KIND=dp), DIMENSION(3, 3, 3, 3), INTENT(OUT)  :: d4f
    1043              :       REAL(KIND=dp), DIMENSION(3, 3, 3, 3, 3), &
    1044              :          INTENT(OUT), OPTIONAL                           :: d5f
    1045              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: v
    1046              :       REAL(KIND=dp), INTENT(IN)                          :: d1, d2, d3, d4
    1047              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: d5
    1048              : 
    1049              :       INTEGER                                            :: k1, k2, k3, k4, k5
    1050              :       REAL(KIND=dp)                                      :: w
    1051              : 
    1052        16178 :       d1f = 0.0_dp
    1053        16178 :       d2f = 0.0_dp
    1054        16178 :       d3f = 0.0_dp
    1055        16178 :       d4f = 0.0_dp
    1056        64712 :       DO k1 = 1, 3
    1057        64712 :          d1f(k1) = d1f(k1) + v(k1)*d1
    1058              :       END DO
    1059        64712 :       DO k1 = 1, 3
    1060       194136 :          DO k2 = 1, 3
    1061       194136 :             d2f(k2, k1) = d2f(k2, k1) + v(k1)*v(k2)*d2
    1062              :          END DO
    1063        64712 :          d2f(k1, k1) = d2f(k1, k1) + d1
    1064              :       END DO
    1065        64712 :       DO k1 = 1, 3
    1066       210314 :          DO k2 = 1, 3
    1067       582408 :             DO k3 = 1, 3
    1068       582408 :                d3f(k3, k2, k1) = d3f(k3, k2, k1) + v(k1)*v(k2)*v(k3)*d3
    1069              :             END DO
    1070       145602 :             w = v(k1)*d2
    1071       145602 :             d3f(k1, k2, k2) = d3f(k1, k2, k2) + w
    1072       145602 :             d3f(k2, k1, k2) = d3f(k2, k1, k2) + w
    1073       194136 :             d3f(k2, k2, k1) = d3f(k2, k2, k1) + w
    1074              :          END DO
    1075              :       END DO
    1076        64712 :       DO k1 = 1, 3
    1077       210314 :          DO k2 = 1, 3
    1078       582408 :             DO k3 = 1, 3
    1079      1747224 :                DO k4 = 1, 3
    1080      1747224 :                   d4f(k4, k3, k2, k1) = d4f(k4, k3, k2, k1) + v(k1)*v(k2)*v(k3)*v(k4)*d4
    1081              :                END DO
    1082       436806 :                w = v(k1)*v(k2)*d3
    1083       436806 :                d4f(k1, k2, k3, k3) = d4f(k1, k2, k3, k3) + w
    1084       436806 :                d4f(k1, k3, k2, k3) = d4f(k1, k3, k2, k3) + w
    1085       436806 :                d4f(k3, k1, k2, k3) = d4f(k3, k1, k2, k3) + w
    1086       436806 :                d4f(k1, k3, k3, k2) = d4f(k1, k3, k3, k2) + w
    1087       436806 :                d4f(k3, k1, k3, k2) = d4f(k3, k1, k3, k2) + w
    1088       582408 :                d4f(k3, k3, k1, k2) = d4f(k3, k3, k1, k2) + w
    1089              :             END DO
    1090       145602 :             d4f(k1, k1, k2, k2) = d4f(k1, k1, k2, k2) + d2
    1091       145602 :             d4f(k1, k2, k1, k2) = d4f(k1, k2, k1, k2) + d2
    1092       194136 :             d4f(k1, k2, k2, k1) = d4f(k1, k2, k2, k1) + d2
    1093              :          END DO
    1094              :       END DO
    1095        16178 :       IF (PRESENT(d5f) .AND. PRESENT(d5)) THEN
    1096            0 :          d5f = 0.0_dp
    1097              : 
    1098            0 :          DO k1 = 1, 3
    1099            0 :             DO k2 = 1, 3
    1100            0 :                DO k3 = 1, 3
    1101            0 :                   DO k4 = 1, 3
    1102            0 :                      DO k5 = 1, 3
    1103            0 :                         d5f(k5, k4, k3, k2, k1) = d5f(k5, k4, k3, k2, k1) + v(k1)*v(k2)*v(k3)*v(k4)*v(k5)*d5
    1104              :                      END DO
    1105            0 :                      w = v(k1)*v(k2)*v(k3)*d4
    1106            0 :                      d5f(k1, k2, k3, k4, k4) = d5f(k1, k2, k3, k4, k4) + w
    1107            0 :                      d5f(k1, k2, k4, k3, k4) = d5f(k1, k2, k4, k3, k4) + w
    1108            0 :                      d5f(k1, k4, k2, k3, k4) = d5f(k1, k4, k2, k3, k4) + w
    1109            0 :                      d5f(k4, k1, k2, k3, k4) = d5f(k4, k1, k2, k3, k4) + w
    1110            0 :                      d5f(k1, k2, k4, k4, k3) = d5f(k1, k2, k4, k4, k3) + w
    1111            0 :                      d5f(k1, k4, k2, k4, k3) = d5f(k1, k4, k2, k4, k3) + w
    1112            0 :                      d5f(k4, k1, k2, k4, k3) = d5f(k4, k1, k2, k4, k3) + w
    1113            0 :                      d5f(k1, k4, k4, k2, k3) = d5f(k1, k4, k4, k2, k3) + w
    1114            0 :                      d5f(k4, k1, k4, k2, k3) = d5f(k4, k1, k4, k2, k3) + w
    1115            0 :                      d5f(k4, k4, k1, k2, k3) = d5f(k4, k4, k1, k2, k3) + w
    1116              :                   END DO
    1117            0 :                   w = v(k1)*d3
    1118            0 :                   d5f(k1, k2, k2, k3, k3) = d5f(k1, k2, k2, k3, k3) + w
    1119            0 :                   d5f(k1, k2, k3, k2, k3) = d5f(k1, k2, k3, k2, k3) + w
    1120            0 :                   d5f(k1, k2, k3, k3, k2) = d5f(k1, k2, k3, k3, k2) + w
    1121            0 :                   d5f(k2, k1, k2, k3, k3) = d5f(k2, k1, k2, k3, k3) + w
    1122            0 :                   d5f(k2, k1, k3, k2, k3) = d5f(k2, k1, k3, k2, k3) + w
    1123            0 :                   d5f(k2, k1, k3, k3, k2) = d5f(k2, k1, k3, k3, k2) + w
    1124            0 :                   d5f(k2, k2, k1, k3, k3) = d5f(k2, k2, k1, k3, k3) + w
    1125            0 :                   d5f(k2, k3, k1, k2, k3) = d5f(k2, k3, k1, k2, k3) + w
    1126            0 :                   d5f(k2, k3, k1, k3, k2) = d5f(k2, k3, k1, k3, k2) + w
    1127            0 :                   d5f(k2, k2, k3, k1, k3) = d5f(k2, k2, k3, k1, k3) + w
    1128            0 :                   d5f(k2, k3, k2, k1, k3) = d5f(k2, k3, k2, k1, k3) + w
    1129            0 :                   d5f(k2, k3, k3, k1, k2) = d5f(k2, k3, k3, k1, k2) + w
    1130            0 :                   d5f(k2, k2, k3, k3, k1) = d5f(k2, k2, k3, k3, k1) + w
    1131            0 :                   d5f(k2, k3, k2, k3, k1) = d5f(k2, k3, k2, k3, k1) + w
    1132            0 :                   d5f(k2, k3, k3, k2, k1) = d5f(k2, k3, k3, k2, k1) + w
    1133              :                END DO
    1134              :             END DO
    1135              :          END DO
    1136              :       END IF
    1137        16178 :    END SUBROUTINE build_d_tensor_gks
    1138              : 
    1139              : ! **************************************************************************************************
    1140              : !> \brief ...
    1141              : !> \param sepi ...
    1142              : !> \param Coul ...
    1143              : !> \param se_int_control ...
    1144              : ! **************************************************************************************************
    1145            9 :    SUBROUTINE makeCoulE0(sepi, Coul, se_int_control)
    1146              :       TYPE(semi_empirical_type), POINTER                 :: sepi
    1147              :       REAL(KIND=dp), DIMENSION(45, 45), INTENT(OUT)      :: Coul
    1148              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
    1149              : 
    1150              :       INTEGER                                            :: gpt, imA, imB, k1, k2, k3, k4, lp, mp, np
    1151            9 :       INTEGER, DIMENSION(:, :), POINTER                  :: bds
    1152              :       REAL(KIND=dp)                                      :: alpha, d2f(3, 3), d4f(3, 3, 3, 3), f, &
    1153              :                                                             ff, w
    1154              :       REAL(KIND=dp), DIMENSION(3)                        :: kk
    1155              :       REAL(KIND=dp), DIMENSION(3, 3, 45)                 :: M2A
    1156              :       REAL(KIND=dp), DIMENSION(3, 45)                    :: M1A
    1157              :       REAL(KIND=dp), DIMENSION(45)                       :: M0A
    1158            9 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: rho0
    1159              :       TYPE(dg_rho0_type), POINTER                        :: dg_rho0
    1160              :       TYPE(dg_type), POINTER                             :: dg
    1161              :       TYPE(pw_grid_type), POINTER                        :: pw_grid
    1162              :       TYPE(pw_pool_type), POINTER                        :: pw_pool
    1163              : 
    1164            9 :       alpha = se_int_control%ewald_gks%alpha
    1165            9 :       pw_pool => se_int_control%ewald_gks%pw_pool
    1166            9 :       dg => se_int_control%ewald_gks%dg
    1167            9 :       CALL dg_get(dg, dg_rho0=dg_rho0)
    1168            9 :       rho0 => dg_rho0%density%array
    1169            9 :       pw_grid => pw_pool%pw_grid
    1170            9 :       bds => pw_grid%bounds
    1171              : 
    1172            9 :       CALL get_se_slater_multipole(sepi, M0A, M1A, M2A)
    1173              : 
    1174            9 :       f = 0.0_dp
    1175            9 :       d2f = 0.0_dp
    1176            9 :       d4f = 0.0_dp
    1177              : 
    1178       150183 :       DO gpt = 1, pw_grid%ngpts_cut_local
    1179       150174 :          lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
    1180       150174 :          mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
    1181       150174 :          np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
    1182              : 
    1183       150174 :          lp = lp + bds(1, 1)
    1184       150174 :          mp = mp + bds(1, 2)
    1185       150174 :          np = np + bds(1, 3)
    1186              : 
    1187       150174 :          IF (pw_grid%gsq(gpt) == 0.0_dp) CYCLE
    1188       600660 :          kk(:) = pw_grid%g(:, gpt)
    1189       150165 :          ff = 2.0_dp*fourpi*rho0(lp, mp, np)**2*pw_grid%vol/pw_grid%gsq(gpt)
    1190              : 
    1191       150165 :          f = f + ff
    1192       600660 :          DO k2 = 1, 3
    1193      1952145 :             DO k1 = 1, 3
    1194      1801980 :                d2f(k1, k2) = d2f(k1, k2) - kk(k1)*kk(k2)*ff
    1195              :             END DO
    1196              :          END DO
    1197       600669 :          DO k4 = 1, 3
    1198      1952154 :             DO k3 = 1, 3
    1199      5856435 :                DO k2 = 1, 3
    1200     17569305 :                   DO k1 = 1, 3
    1201     16217820 :                      d4f(k1, k2, k3, k4) = d4f(k1, k2, k3, k4) + kk(k1)*kk(k2)*kk(k3)*kk(k4)*ff
    1202              :                   END DO
    1203              :                END DO
    1204              :             END DO
    1205              :          END DO
    1206              : 
    1207              :       END DO
    1208              : 
    1209           99 :       DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2
    1210          999 :          DO imB = 1, (sepi%natorb*(sepi%natorb + 1))/2
    1211              : 
    1212          900 :             w = M0A(imA)*M0A(imB)*f
    1213         3600 :             DO k2 = 1, 3
    1214        11700 :                DO k1 = 1, 3
    1215        10800 :                   w = w + (M2A(k1, k2, imA)*M0A(imB) - M1A(k1, imA)*M1A(k2, imB) + M0A(imA)*M2A(k1, k2, imB))*d2f(k1, k2)
    1216              :                END DO
    1217              :             END DO
    1218              : 
    1219         3600 :             DO k4 = 1, 3
    1220        11700 :                DO k3 = 1, 3
    1221        35100 :                   DO k2 = 1, 3
    1222       105300 :                      DO k1 = 1, 3
    1223        97200 :                         w = w + M2A(k1, k2, imA)*M2A(k3, k4, imB)*d4f(k1, k2, k3, k4)
    1224              :                      END DO
    1225              :                   END DO
    1226              :                END DO
    1227              :             END DO
    1228              : 
    1229          990 :             Coul(imA, imB) = w
    1230              : 
    1231              :          END DO
    1232              :       END DO
    1233              : 
    1234           99 :       DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2
    1235          999 :          DO imB = 1, (sepi%natorb*(sepi%natorb + 1))/2
    1236          900 :             w = -M0A(imA)*M0A(imB)*0.25_dp*fourpi/(pw_grid%vol*alpha**2)
    1237          990 :             Coul(imA, imB) = Coul(imA, imB) + w
    1238              :          END DO
    1239              :       END DO
    1240              : 
    1241           99 :       DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2
    1242          999 :          DO imB = 1, (sepi%natorb*(sepi%natorb + 1))/2
    1243              : 
    1244          900 :             w = M0A(imA)*M0A(imB)
    1245          900 :             Coul(imA, imB) = Coul(imA, imB) - 2.0_dp*alpha*oorootpi*w
    1246              : 
    1247          900 :             w = 0.0_dp
    1248         3600 :             DO k1 = 1, 3
    1249         2700 :                w = w + M1A(k1, imA)*M1A(k1, imB)
    1250         2700 :                w = w - M0A(imA)*M2A(k1, k1, imB)
    1251         3600 :                w = w - M2A(k1, k1, imA)*M0A(imB)
    1252              :             END DO
    1253          900 :             Coul(imA, imB) = Coul(imA, imB) - 4.0_dp*alpha**3*oorootpi*w/3.0_dp
    1254              : 
    1255          900 :             w = 0.0_dp
    1256         3600 :             DO k2 = 1, 3
    1257        11700 :                DO k1 = 1, 3
    1258         8100 :                   w = w + 2.0_dp*M2A(k1, k2, imA)*M2A(k1, k2, imB)
    1259        10800 :                   w = w + M2A(k1, k1, imA)*M2A(k2, k2, imB)
    1260              :                END DO
    1261              :             END DO
    1262          990 :             Coul(imA, imB) = Coul(imA, imB) - 8.0_dp*alpha**5*oorootpi*w/5.0_dp
    1263              :          END DO
    1264              :       END DO
    1265            9 :    END SUBROUTINE makeCoulE0
    1266              : 
    1267              : ! **************************************************************************************************
    1268              : !> \brief Retrieves the multipole for the Slater integral evaluation
    1269              : !> \param sepi ...
    1270              : !> \param M0 ...
    1271              : !> \param M1 ...
    1272              : !> \param M2 ...
    1273              : !> \param ACOUL ...
    1274              : ! **************************************************************************************************
    1275        32365 :    SUBROUTINE get_se_slater_multipole(sepi, M0, M1, M2, ACOUL)
    1276              :       TYPE(semi_empirical_type), POINTER                 :: sepi
    1277              :       REAL(kind=dp), DIMENSION(45), INTENT(OUT)          :: M0
    1278              :       REAL(kind=dp), DIMENSION(3, 45), INTENT(OUT)       :: M1
    1279              :       REAL(kind=dp), DIMENSION(3, 3, 45), INTENT(OUT)    :: M2
    1280              :       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: ACOUL
    1281              : 
    1282              :       INTEGER                                            :: i, j, jint, size_1c_int
    1283              :       TYPE(semi_empirical_mpole_type), POINTER           :: mpole
    1284              : 
    1285        32365 :       NULLIFY (mpole)
    1286        32365 :       size_1c_int = SIZE(sepi%w_mpole)
    1287       356015 :       DO jint = 1, size_1c_int
    1288       323650 :          mpole => sepi%w_mpole(jint)%mpole
    1289       323650 :          i = mpole%indi
    1290       323650 :          j = mpole%indj
    1291       323650 :          M0(indexb(i, j)) = -mpole%cs
    1292              : 
    1293       323650 :          M1(1, indexb(i, j)) = -mpole%ds(1)
    1294       323650 :          M1(2, indexb(i, j)) = -mpole%ds(2)
    1295       323650 :          M1(3, indexb(i, j)) = -mpole%ds(3)
    1296              : 
    1297       323650 :          M2(1, 1, indexb(i, j)) = -mpole%qq(1, 1)/3.0_dp
    1298       323650 :          M2(2, 1, indexb(i, j)) = -mpole%qq(2, 1)/3.0_dp
    1299       323650 :          M2(3, 1, indexb(i, j)) = -mpole%qq(3, 1)/3.0_dp
    1300              : 
    1301       323650 :          M2(1, 2, indexb(i, j)) = -mpole%qq(1, 2)/3.0_dp
    1302       323650 :          M2(2, 2, indexb(i, j)) = -mpole%qq(2, 2)/3.0_dp
    1303       323650 :          M2(3, 2, indexb(i, j)) = -mpole%qq(3, 2)/3.0_dp
    1304              : 
    1305       323650 :          M2(1, 3, indexb(i, j)) = -mpole%qq(1, 3)/3.0_dp
    1306       323650 :          M2(2, 3, indexb(i, j)) = -mpole%qq(2, 3)/3.0_dp
    1307       356015 :          M2(3, 3, indexb(i, j)) = -mpole%qq(3, 3)/3.0_dp
    1308              :       END DO
    1309        32365 :       IF (PRESENT(ACOUL)) ACOUL = sepi%acoul
    1310        32365 :    END SUBROUTINE get_se_slater_multipole
    1311              : 
    1312              : END MODULE semi_empirical_int_gks
        

Generated by: LCOV version 2.0-1