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

            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 Calculation of the integrals over solid harmonic Gaussian(SHG) functions.
      10              : !>        Routines for (a|O(r12)|b) and overlap integrals (ab), (aba) and (abb).
      11              : !> \par Literature (partly)
      12              : !>      T.J. Giese and D. M. York, J. Chem. Phys, 128, 064104 (2008)
      13              : !>      T. Helgaker, P Joergensen, J. Olsen, Molecular Electronic-Structure
      14              : !>                                           Theory, Wiley
      15              : !> \par History
      16              : !>      created [04.2015]
      17              : !> \author Dorothea Golze
      18              : ! **************************************************************************************************
      19              : MODULE construct_shg
      20              :    USE kinds,                           ONLY: dp
      21              :    USE mathconstants,                   ONLY: dfac,&
      22              :                                               fac
      23              :    USE orbital_pointers,                ONLY: indso,&
      24              :                                               indso_inv,&
      25              :                                               nsoset
      26              : #include "../base/base_uses.f90"
      27              : 
      28              :    IMPLICIT NONE
      29              : 
      30              :    PRIVATE
      31              : 
      32              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'construct_shg'
      33              : 
      34              : ! *** Public subroutines ***
      35              :    PUBLIC :: get_real_scaled_solid_harmonic, get_W_matrix, get_dW_matrix, &
      36              :              construct_int_shg_ab, construct_dev_shg_ab, construct_overlap_shg_aba, &
      37              :              dev_overlap_shg_aba, construct_overlap_shg_abb, dev_overlap_shg_abb
      38              : 
      39              : CONTAINS
      40              : 
      41              : ! **************************************************************************************************
      42              : !> \brief computes the real scaled solid harmonics Rlm up to a given l
      43              : !> \param Rlm_c cosine part of real scaled soldi harmonics
      44              : !> \param Rlm_s sine part of real scaled soldi harmonics
      45              : !> \param l maximal l quantum up to where Rlm is calculated
      46              : !> \param r distance vector between a and b
      47              : !> \param r2 square of distance vector
      48              : ! **************************************************************************************************
      49      2542961 :    SUBROUTINE get_real_scaled_solid_harmonic(Rlm_c, Rlm_s, l, r, r2)
      50              : 
      51              :       INTEGER, INTENT(IN)                                :: l
      52              :       REAL(KIND=dp), DIMENSION(0:l, -2*l:2*l), &
      53              :          INTENT(OUT)                                     :: Rlm_s, Rlm_c
      54              :       REAL(KIND=dp), DIMENSION(3)                        :: r
      55              :       REAL(KIND=dp)                                      :: r2
      56              : 
      57              :       INTEGER                                            :: li, mi, prefac
      58              :       REAL(KIND=dp)                                      :: Rc, Rc_00, Rlm, Rmlm, Rplm, Rs, Rs_00, &
      59              :                                                             temp_c
      60              : 
      61      2542961 :       Rc_00 = 1.0_dp
      62      2542961 :       Rs_00 = 0.0_dp
      63              : 
      64      2542961 :       Rlm_c(0, 0) = Rc_00
      65      2542961 :       Rlm_s(0, 0) = Rs_00
      66              : 
      67              :       ! generate elements Rmm
      68              :       ! start
      69      2542961 :       IF (l > 0) THEN
      70      1535101 :          Rc = -0.5_dp*r(1)*Rc_00
      71      1535101 :          Rs = -0.5_dp*r(2)*Rc_00
      72      1535101 :          Rlm_c(1, 1) = Rc
      73      1535101 :          Rlm_s(1, 1) = Rs
      74      1535101 :          Rlm_c(1, -1) = -Rc
      75      1535101 :          Rlm_s(1, -1) = Rs
      76              :       END IF
      77      3054060 :       DO li = 2, l
      78       511099 :          temp_c = (-r(1)*Rc + r(2)*Rs)/(REAL(2*(li - 1) + 2, dp))
      79       511099 :          Rs = (-r(2)*Rc - r(1)*Rs)/(REAL(2*(li - 1) + 2, dp))
      80       511099 :          Rc = temp_c
      81       511099 :          Rlm_c(li, li) = Rc
      82       511099 :          Rlm_s(li, li) = Rs
      83      3054060 :          IF (MODULO(li, 2) /= 0) THEN
      84       122986 :             Rlm_c(li, -li) = -Rc
      85       122986 :             Rlm_s(li, -li) = Rs
      86              :          ELSE
      87       388113 :             Rlm_c(li, -li) = Rc
      88       388113 :             Rlm_s(li, -li) = -Rs
      89              :          END IF
      90              :       END DO
      91              : 
      92      4589161 :       DO mi = 0, l - 1
      93      2046200 :          Rmlm = Rlm_c(mi, mi)
      94      2046200 :          Rlm = r(3)*Rlm_c(mi, mi)
      95      2046200 :          Rlm_c(mi + 1, mi) = Rlm
      96      2046200 :          IF (MODULO(mi, 2) /= 0) THEN
      97       388113 :             Rlm_c(mi + 1, -mi) = -Rlm
      98              :          ELSE
      99      1658087 :             Rlm_c(mi + 1, -mi) = Rlm
     100              :          END IF
     101      5296958 :          DO li = mi + 2, l
     102       707797 :             prefac = (li + mi)*(li - mi)
     103       707797 :             Rplm = (REAL(2*li - 1, dp)*r(3)*Rlm - r2*Rmlm)/REAL(prefac, dp)
     104       707797 :             Rmlm = Rlm
     105       707797 :             Rlm = Rplm
     106       707797 :             Rlm_c(li, mi) = Rlm
     107      2753997 :             IF (MODULO(mi, 2) /= 0) THEN
     108       159842 :                Rlm_c(li, -mi) = -Rlm
     109              :             ELSE
     110       547955 :                Rlm_c(li, -mi) = Rlm
     111              :             END IF
     112              :          END DO
     113              :       END DO
     114      3054060 :       DO mi = 1, l - 1
     115       511099 :          Rmlm = Rlm_s(mi, mi)
     116       511099 :          Rlm = r(3)*Rlm_s(mi, mi)
     117       511099 :          Rlm_s(mi + 1, mi) = Rlm
     118       511099 :          IF (MODULO(mi, 2) /= 0) THEN
     119       388113 :             Rlm_s(mi + 1, -mi) = Rlm
     120              :          ELSE
     121       122986 :             Rlm_s(mi + 1, -mi) = -Rlm
     122              :          END IF
     123      3250758 :          DO li = mi + 2, l
     124       196698 :             prefac = (li + mi)*(li - mi)
     125       196698 :             Rplm = (REAL(2*li - 1, dp)*r(3)*Rlm - r2*Rmlm)/REAL(prefac, dp)
     126       196698 :             Rmlm = Rlm
     127       196698 :             Rlm = Rplm
     128       196698 :             Rlm_s(li, mi) = Rlm
     129       707797 :             IF (MODULO(mi, 2) /= 0) THEN
     130       159842 :                Rlm_s(li, -mi) = Rlm
     131              :             ELSE
     132        36856 :                Rlm_s(li, -mi) = -Rlm
     133              :             END IF
     134              :          END DO
     135              :       END DO
     136              : 
     137      2542961 :    END SUBROUTINE get_real_scaled_solid_harmonic
     138              : 
     139              : ! **************************************************************************************************
     140              : !> \brief Calculate the prefactor A(l,m) = (-1)^m \sqrt[(2-delta(m,0))(l+m)!(l-m)!]
     141              : !> \param lmax maximal l quantum number
     142              : !> \param A matrix storing the prefactor for a given l and m
     143              : ! **************************************************************************************************
     144      2542961 :    SUBROUTINE get_Alm(lmax, A)
     145              : 
     146              :       INTEGER, INTENT(IN)                                :: lmax
     147              :       REAL(KIND=dp), DIMENSION(0:lmax, 0:lmax), &
     148              :          INTENT(INOUT)                                   :: A
     149              : 
     150              :       INTEGER                                            :: l, m
     151              :       REAL(KIND=dp)                                      :: temp
     152              : 
     153      7132122 :       DO l = 0, lmax
     154     14475280 :          DO m = 0, l
     155      7343158 :             temp = SQRT(fac(l + m)*fac(l - m))
     156      7343158 :             IF (MODULO(m, 2) /= 0) temp = -temp
     157      7343158 :             IF (m /= 0) temp = temp*SQRT(2.0_dp)
     158     11932319 :             A(l, m) = temp
     159              :          END DO
     160              :       END DO
     161              : 
     162      2542961 :    END SUBROUTINE get_Alm
     163              : 
     164              : ! **************************************************************************************************
     165              : !> \brief calculates the prefactors for the derivatives of the W matrix
     166              : !> \param lmax maximal l quantum number
     167              : !> \param dA_p = A(l,m)/A(l-1,m+1)
     168              : !> \param dA_m = A(l,m)/A(l-1,m-1)
     169              : !> \param dA = A(l,m)/A(l-1,m)
     170              : !> \note for m=0, W_l-1,-1 can't be read from Waux_mat, but we use
     171              : !>       W_l-1,-1 = -W_l-1,1 [cc(1), cs(2)] or W_l-1,-1 = W_l-1,1 [[sc(3), ss(4)], i.e.
     172              : !>       effectively we multiply dA_p by 2
     173              : ! **************************************************************************************************
     174         6396 :    SUBROUTINE get_dA_prefactors(lmax, dA_p, dA_m, dA)
     175              : 
     176              :       INTEGER, INTENT(IN)                                :: lmax
     177              :       REAL(KIND=dp), DIMENSION(0:lmax, 0:lmax), &
     178              :          INTENT(INOUT)                                   :: dA_p, dA_m, dA
     179              : 
     180              :       INTEGER                                            :: l, m
     181              :       REAL(KIND=dp)                                      :: bm, bm_m, bm_p
     182              : 
     183        45356 :       DO l = 0, lmax
     184       185688 :          DO m = 0, l
     185       140332 :             bm = 1.0_dp
     186       140332 :             bm_m = 1.0_dp
     187       140332 :             bm_p = 1.0_dp
     188       140332 :             IF (m /= 0) bm = SQRT(2.0_dp)
     189       101372 :             IF (m - 1 /= 0) bm_m = SQRT(2.0_dp)
     190       107768 :             IF (m + 1 /= 0) bm_p = SQRT(2.0_dp)
     191       140332 :             dA_p(l, m) = -bm/bm_p*SQRT(REAL((l - m)*(l - m - 1), dp))
     192       140332 :             dA_m(l, m) = -bm/bm_m*SQRT(REAL((l + m)*(l + m - 1), dp))
     193       140332 :             dA(l, m) = 2.0_dp*SQRT(REAL((l + m)*(l - m), dp))
     194       179292 :             IF (m == 0) dA_p(l, m) = 2.0_dp*dA_p(l, m)
     195              :          END DO
     196              :       END DO
     197         6396 :    END SUBROUTINE get_dA_prefactors
     198              : 
     199              : ! **************************************************************************************************
     200              : !> \brief calculates the angular dependent-part of the SHG integrals,
     201              : !>        transformation matrix W, see literature above
     202              : !> \param lamax array of maximal l quantum number on a;
     203              : !>        lamax(lb) with lb= 0..lbmax
     204              : !> \param lbmax maximal l quantum number on b
     205              : !> \param lmax maximal l quantum number
     206              : !> \param Rc cosine part of real scaled solid harmonics
     207              : !> \param Rs sine part of real scaled solid harmonics
     208              : !> \param Waux_mat stores the angular-dependent part of the SHG integrals
     209              : ! **************************************************************************************************
     210      2542961 :    SUBROUTINE get_W_matrix(lamax, lbmax, lmax, Rc, Rs, Waux_mat)
     211              : 
     212              :       INTEGER, DIMENSION(:), POINTER                     :: lamax
     213              :       INTEGER, INTENT(IN)                                :: lbmax, lmax
     214              :       REAL(KIND=dp), DIMENSION(0:lmax, -2*lmax:2*lmax), &
     215              :          INTENT(IN)                                      :: Rc, Rs
     216              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: Waux_mat
     217              : 
     218              :       INTEGER                                            :: j, k, la, labmin, laj, lb, lbj, ma, &
     219              :                                                             ma_m, ma_p, mb, mb_m, mb_p, nla, nlb
     220              :       REAL(KIND=dp) :: A_jk, A_lama, A_lbmb, Alm_fac, delta_k, prefac, Rca_m, Rca_p, Rcb_m, Rcb_p, &
     221              :          Rsa_m, Rsa_p, Rsb_m, Rsb_p, sign_fac, Wa(4), Wb(4), Wmat(4)
     222      2542961 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: A
     223              : 
     224      2542961 :       Wa(:) = 0.0_dp
     225      2542961 :       Wb(:) = 0.0_dp
     226      2542961 :       Wmat(:) = 0.0_dp
     227              : 
     228     10171844 :       ALLOCATE (A(0:lmax, 0:lmax))
     229      2542961 :       CALL get_Alm(lmax, A)
     230              : 
     231      6467253 :       DO lb = 0, lbmax
     232      3924292 :          nlb = nsoset(lb - 1)
     233     13366409 :          DO la = 0, lamax(lb)
     234      6899156 :             nla = nsoset(la - 1)
     235      6899156 :             labmin = MIN(la, lb)
     236     22064935 :             DO mb = 0, lb
     237     11241487 :                A_lbmb = A(lb, mb)
     238     11241487 :                IF (MODULO(lb, 2) /= 0) A_lbmb = -A_lbmb
     239     38632144 :                DO ma = 0, la
     240     20491501 :                   A_lama = A(la, ma)
     241     20491501 :                   Alm_fac = A_lama*A_lbmb
     242     69764431 :                   DO j = 0, labmin
     243     38031443 :                      laj = la - j
     244     38031443 :                      lbj = lb - j
     245     38031443 :                      prefac = Alm_fac*REAL(2**(la + lb - j), dp)*dfac(2*j - 1)
     246     38031443 :                      delta_k = 0.5_dp
     247     38031443 :                      Wmat = 0.0_dp
     248    101063762 :                      DO k = 0, j
     249     63032319 :                         ma_m = ma - k
     250     63032319 :                         ma_p = ma + k
     251     63032319 :                         IF (laj < ABS(ma_m) .AND. laj < ABS(ma_p)) CYCLE
     252     46707644 :                         mb_m = mb - k
     253     46707644 :                         mb_p = mb + k
     254     46707644 :                         IF (lbj < ABS(mb_m) .AND. lbj < ABS(mb_p)) CYCLE
     255     37066246 :                         IF (k /= 0) delta_k = 1.0_dp
     256     37066246 :                         A_jk = fac(j + k)*fac(j - k)
     257     37066246 :                         IF (k /= 0) A_jk = 2.0_dp*A_jk
     258     37066246 :                         IF (MODULO(k, 2) /= 0) THEN
     259              :                            sign_fac = -1.0_dp
     260              :                         ELSE
     261     27655873 :                            sign_fac = 1.0_dp
     262              :                         END IF
     263     37066246 :                         Rca_m = Rc(laj, ma_m)
     264     37066246 :                         Rsa_m = Rs(laj, ma_m)
     265     37066246 :                         Rca_p = Rc(laj, ma_p)
     266     37066246 :                         Rsa_p = Rs(laj, ma_p)
     267     37066246 :                         Rcb_m = Rc(lbj, mb_m)
     268     37066246 :                         Rsb_m = Rs(lbj, mb_m)
     269     37066246 :                         Rcb_p = Rc(lbj, mb_p)
     270     37066246 :                         Rsb_p = Rs(lbj, mb_p)
     271     37066246 :                         Wa(1) = delta_k*(Rca_m + sign_fac*Rca_p)
     272     37066246 :                         Wb(1) = delta_k*(Rcb_m + sign_fac*Rcb_p)
     273     37066246 :                         Wa(2) = -Rsa_m + sign_fac*Rsa_p
     274     37066246 :                         Wb(2) = -Rsb_m + sign_fac*Rsb_p
     275     37066246 :                         Wmat(1) = Wmat(1) + prefac/A_jk*(Wa(1)*Wb(1) + Wa(2)*Wb(2))
     276     37066246 :                         IF (mb > 0) THEN
     277     19855456 :                            Wb(3) = delta_k*(Rsb_m + sign_fac*Rsb_p)
     278     19855456 :                            Wb(4) = Rcb_m - sign_fac*Rcb_p
     279     19855456 :                            Wmat(2) = Wmat(2) + prefac/A_jk*(Wa(1)*Wb(3) + Wa(2)*Wb(4))
     280              :                         END IF
     281     37066246 :                         IF (ma > 0) THEN
     282     20524281 :                            Wa(3) = delta_k*(Rsa_m + sign_fac*Rsa_p)
     283     20524281 :                            Wa(4) = Rca_m - sign_fac*Rca_p
     284     20524281 :                            Wmat(3) = Wmat(3) + prefac/A_jk*(Wa(3)*Wb(1) + Wa(4)*Wb(2))
     285              :                         END IF
     286     75097689 :                         IF (ma > 0 .AND. mb > 0) THEN
     287     12767595 :                            Wmat(4) = Wmat(4) + prefac/A_jk*(Wa(3)*Wb(3) + Wa(4)*Wb(4))
     288              :                         END IF
     289              :                      END DO
     290     38031443 :                      Waux_mat(j + 1, nla + la + 1 + ma, nlb + lb + 1 + mb) = Wmat(1)
     291     38031443 :                      IF (mb > 0) Waux_mat(j + 1, nla + la + 1 + ma, nlb + lb + 1 - mb) = Wmat(2)
     292     38031443 :                      IF (ma > 0) Waux_mat(j + 1, nla + la + 1 - ma, nlb + lb + 1 + mb) = Wmat(3)
     293     58522944 :                      IF (ma > 0 .AND. mb > 0) Waux_mat(j + 1, nla + la + 1 - ma, nlb + lb + 1 - mb) = Wmat(4)
     294              :                   END DO
     295              :                END DO
     296              :             END DO
     297              :          END DO
     298              :       END DO
     299              : 
     300      2542961 :       DEALLOCATE (A)
     301              : 
     302      2542961 :    END SUBROUTINE get_W_matrix
     303              : 
     304              : ! **************************************************************************************************
     305              : !> \brief calculates derivatives of transformation matrix W,
     306              : !> \param lamax array of maximal l quantum number on a;
     307              : !>        lamax(lb) with lb= 0..lbmax
     308              : !> \param lbmax maximal l quantum number on b
     309              : !> \param Waux_mat stores the angular-dependent part of the SHG integrals
     310              : !> \param dWaux_mat stores the derivatives of the angular-dependent part of
     311              : !>        the SHG integrals
     312              : ! **************************************************************************************************
     313         6396 :    SUBROUTINE get_dW_matrix(lamax, lbmax, Waux_mat, dWaux_mat)
     314              : 
     315              :       INTEGER, DIMENSION(:), POINTER                     :: lamax
     316              :       INTEGER, INTENT(IN)                                :: lbmax
     317              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: Waux_mat
     318              :       REAL(KIND=dp), DIMENSION(:, :, :, :), &
     319              :          INTENT(INOUT)                                   :: dWaux_mat
     320              : 
     321              :       INTEGER                                            :: ima, imam, imb, imbm, ipa, ipam, ipb, &
     322              :                                                             ipbm, j, jmax, la, labm, labmin, lamb, &
     323              :                                                             lb, lmax, ma, mb, nla, nlam, nlb, nlbm
     324              :       REAL(KIND=dp)                                      :: dAa, dAa_m, dAa_p, dAb, dAb_m, dAb_p
     325         6396 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: dA, dA_m, dA_p, Wam, Wamm, Wamp, Wbm, &
     326         6396 :                                                             Wbmm, Wbmp
     327              : 
     328        32680 :       jmax = MIN(MAXVAL(lamax), lbmax)
     329        38376 :       ALLOCATE (Wam(0:jmax, 4), Wamm(0:jmax, 4), Wamp(0:jmax, 4))
     330        25584 :       ALLOCATE (Wbm(0:jmax, 4), Wbmm(0:jmax, 4), Wbmp(0:jmax, 4))
     331              : 
     332              :       !*** get dA_p=A(l,m)/A(l-1,m+1)
     333              :       !*** get dA_m=A(l,m)/A(l-1,m-1)
     334              :       !*** get dA=2*A(l,m)/A(l-1,m)
     335        32680 :       lmax = MAX(MAXVAL(lamax), lbmax)
     336        51168 :       ALLOCATE (dA_p(0:lmax, 0:lmax), dA_m(0:lmax, 0:lmax), dA(0:lmax, 0:lmax))
     337         6396 :       CALL get_dA_prefactors(lmax, dA_p, dA_m, dA)
     338              : 
     339        32680 :       DO lb = 0, lbmax
     340        26284 :          nlb = nsoset(lb - 1)
     341        26284 :          nlbm = 0
     342        26284 :          IF (lb > 0) nlbm = nsoset(lb - 2)
     343       179986 :          DO la = 0, lamax(lb)
     344       147306 :             nla = nsoset(la - 1)
     345       147306 :             nlam = 0
     346       147306 :             IF (la > 0) nlam = nsoset(la - 2)
     347       147306 :             labmin = MIN(la, lb)
     348       147306 :             lamb = MIN(la - 1, lb)
     349       147306 :             labm = MIN(la, lb - 1)
     350       537322 :             DO mb = 0, lb
     351       363732 :                dAb = dA(lb, mb)
     352       363732 :                dAb_p = dA_p(lb, mb)
     353       363732 :                dAb_m = dA_m(lb, mb)
     354       363732 :                ipb = nlb + lb + mb + 1
     355       363732 :                imb = nlb + lb - mb + 1
     356       363732 :                ipbm = nlbm + lb + mb
     357       363732 :                imbm = nlbm + lb - mb
     358      1730862 :                DO ma = 0, la
     359      1219824 :                   dAa = dA(la, ma)
     360      1219824 :                   dAa_p = dA_p(la, ma)
     361      1219824 :                   dAa_m = dA_m(la, ma)
     362      1219824 :                   ipa = nla + la + ma + 1
     363      1219824 :                   ima = nla + la - ma + 1
     364      1219824 :                   ipam = nlam + la + ma
     365      1219824 :                   imam = nlam + la - ma
     366     27537348 :                   Wam(:, :) = 0.0_dp
     367     27537348 :                   Wamm(:, :) = 0.0_dp
     368     27537348 :                   Wamp(:, :) = 0.0_dp
     369              :                   !*** Wam: la-1, ma
     370      1219824 :                   IF (ma <= la - 1) THEN
     371      2924270 :                      Wam(0:lamb, 1) = Waux_mat(1:lamb + 1, ipam, ipb)
     372      2183157 :                      IF (mb > 0) Wam(0:lamb, 2) = Waux_mat(1:lamb + 1, ipam, imb)
     373      2305787 :                      IF (ma > 0) Wam(0:lamb, 3) = Waux_mat(1:lamb + 1, imam, ipb)
     374      1787750 :                      IF (ma > 0 .AND. mb > 0) Wam(0:lamb, 4) = Waux_mat(1:lamb + 1, imam, imb)
     375              :                   END IF
     376              :                   !*** Wamm: la-1, ma-1
     377      1219824 :                   IF (ma - 1 >= 0) THEN
     378      2924270 :                      Wamm(0:lamb, 1) = Waux_mat(1:lamb + 1, ipam - 1, ipb)
     379      2183157 :                      IF (mb > 0) Wamm(0:lamb, 2) = Waux_mat(1:lamb + 1, ipam - 1, imb)
     380      2305787 :                      IF (ma - 1 > 0) Wamm(0:lamb, 3) = Waux_mat(1:lamb + 1, imam + 1, ipb) !order: e.g. -1 0 1, if < 0 |m|, -1 means  -m+1
     381      1787750 :                      IF (ma - 1 > 0 .AND. mb > 0) Wamm(0:lamb, 4) = Waux_mat(1:lamb + 1, imam + 1, imb)
     382              :                   END IF
     383              :                   !*** Wamp: la-1, ma+1
     384      1219824 :                   IF (ma + 1 <= la - 1) THEN
     385      2009941 :                      Wamp(0:lamb, 1) = Waux_mat(1:lamb + 1, ipam + 1, ipb)
     386      1491904 :                      IF (mb > 0) Wamp(0:lamb, 2) = Waux_mat(1:lamb + 1, ipam + 1, imb)
     387      2009941 :                      IF (ma + 1 > 0) Wamp(0:lamb, 3) = Waux_mat(1:lamb + 1, imam - 1, ipb)
     388      1491904 :                      IF (ma + 1 > 0 .AND. mb > 0) Wamp(0:lamb, 4) = Waux_mat(1:lamb + 1, imam - 1, imb)
     389              :                   END IF
     390     27537348 :                   Wbm(:, :) = 0.0_dp
     391     27537348 :                   Wbmm(:, :) = 0.0_dp
     392     27537348 :                   Wbmp(:, :) = 0.0_dp
     393              :                   !*** Wbm: lb-1, mb
     394      1219824 :                   IF (mb <= lb - 1) THEN
     395      2266824 :                      Wbm(0:labm, 1) = Waux_mat(1:labm + 1, ipa, ipbm)
     396      1597563 :                      IF (mb > 0) Wbm(0:labm, 2) = Waux_mat(1:labm + 1, ipa, imbm)
     397      1838473 :                      IF (ma > 0) Wbm(0:labm, 3) = Waux_mat(1:labm + 1, ima, ipbm)
     398      1354185 :                      IF (ma > 0 .AND. mb > 0) Wbm(0:labm, 4) = Waux_mat(1:labm + 1, ima, imbm)
     399              :                   END IF
     400              :                   !*** Wbmm: lb-1, mb-1
     401      1219824 :                   IF (mb - 1 >= 0) THEN
     402      2266824 :                      Wbmm(0:labm, 1) = Waux_mat(1:labm + 1, ipa, ipbm - 1)
     403      1597563 :                      IF (mb - 1 > 0) Wbmm(0:labm, 2) = Waux_mat(1:labm + 1, ipa, imbm + 1)
     404      1838473 :                      IF (ma > 0) Wbmm(0:labm, 3) = Waux_mat(1:labm + 1, ima, ipbm - 1)
     405      1354185 :                      IF (ma > 0 .AND. mb - 1 > 0) Wbmm(0:labm, 4) = Waux_mat(1:labm + 1, ima, imbm + 1)
     406              :                   END IF
     407              :                   !*** Wbmp: lb-1, mb+1
     408      1219824 :                   IF (mb + 1 <= lb - 1) THEN
     409      1229384 :                      Wbmp(0:labm, 1) = Waux_mat(1:labm + 1, ipa, ipbm + 1)
     410      1229384 :                      IF (mb + 1 > 0) Wbmp(0:labm, 2) = Waux_mat(1:labm + 1, ipa, imbm - 1)
     411       986006 :                      IF (ma > 0) Wbmp(0:labm, 3) = Waux_mat(1:labm + 1, ima, ipbm + 1)
     412       986006 :                      IF (ma > 0 .AND. mb + 1 > 0) Wbmp(0:labm, 4) = Waux_mat(1:labm + 1, ima, imbm - 1)
     413              :                   END IF
     414      4751406 :                   DO j = 0, labmin
     415              :                      !*** x component
     416              :                      dWaux_mat(1, j + 1, ipa, ipb) = dAa_p*Wamp(j, 1) - dAa_m*Wamm(j, 1) &
     417      3167850 :                                                      - dAb_p*Wbmp(j, 1) + dAb_m*Wbmm(j, 1)
     418      3167850 :                      IF (mb > 0) THEN
     419              :                         dWaux_mat(1, j + 1, ipa, imb) = dAa_p*Wamp(j, 2) - dAa_m*Wamm(j, 2) &
     420      2064253 :                                                         - dAb_p*Wbmp(j, 2) + dAb_m*Wbmm(j, 2)
     421              :                      END IF
     422      3167850 :                      IF (ma > 0) THEN
     423              :                         dWaux_mat(1, j + 1, ima, ipb) = dAa_p*Wamp(j, 3) - dAa_m*Wamm(j, 3) &
     424      2336904 :                                                         - dAb_p*Wbmp(j, 3) + dAb_m*Wbmm(j, 3)
     425              :                      END IF
     426      3167850 :                      IF (ma > 0 .AND. mb > 0) THEN
     427              :                         dWaux_mat(1, j + 1, ima, imb) = dAa_p*Wamp(j, 4) - dAa_m*Wamm(j, 4) &
     428      1523984 :                                                         - dAb_p*Wbmp(j, 4) + dAb_m*Wbmm(j, 4)
     429              :                      END IF
     430              : 
     431              :                      !**** y component
     432              :                      dWaux_mat(2, j + 1, ipa, ipb) = dAa_p*Wamp(j, 3) + dAa_m*Wamm(j, 3) &
     433      3167850 :                                                      - dAb_p*Wbmp(j, 2) - dAb_m*Wbmm(j, 2)
     434      3167850 :                      IF (mb > 0) THEN
     435              :                         dWaux_mat(2, j + 1, ipa, imb) = dAa_p*Wamp(j, 4) + dAa_m*Wamm(j, 4) &
     436      2064253 :                                                         + dAb_p*Wbmp(j, 1) + dAb_m*Wbmm(j, 1)
     437              :                      END IF
     438      3167850 :                      IF (ma > 0) THEN
     439              :                         dWaux_mat(2, j + 1, ima, ipb) = -dAa_p*Wamp(j, 1) - dAa_m*Wamm(j, 1) &
     440      2336904 :                                                         - dAb_p*Wbmp(j, 4) - dAb_m*Wbmm(j, 4)
     441              :                      END IF
     442      3167850 :                      IF (ma > 0 .AND. mb > 0) THEN
     443              :                         dWaux_mat(2, j + 1, ima, imb) = -dAa_p*Wamp(j, 2) - dAa_m*Wamm(j, 2) &
     444      1523984 :                                                         + dAb_p*Wbmp(j, 3) + dAb_m*Wbmm(j, 3)
     445              :                      END IF
     446              :                      !**** z compnent
     447      3167850 :                      dWaux_mat(3, j + 1, ipa, ipb) = dAa*Wam(j, 1) - dAb*Wbm(j, 1)
     448      3167850 :                      IF (mb > 0) THEN
     449      2064253 :                         dWaux_mat(3, j + 1, ipa, imb) = dAa*Wam(j, 2) - dAb*Wbm(j, 2)
     450              :                      END IF
     451      3167850 :                      IF (ma > 0) THEN
     452      2336904 :                         dWaux_mat(3, j + 1, ima, ipb) = dAa*Wam(j, 3) - dAb*Wbm(j, 3)
     453              :                      END IF
     454      4387674 :                      IF (ma > 0 .AND. mb > 0) THEN
     455      1523984 :                         dWaux_mat(3, j + 1, ima, imb) = dAa*Wam(j, 4) - dAb*Wbm(j, 4)
     456              :                      END IF
     457              : 
     458              :                   END DO
     459              :                END DO
     460              :             END DO
     461              :          END DO
     462              :       END DO
     463              : 
     464         6396 :       DEALLOCATE (Wam, Wamm, Wamp)
     465         6396 :       DEALLOCATE (Wbm, Wbmm, Wbmp)
     466         6396 :       DEALLOCATE (dA, dA_p, dA_m)
     467              : 
     468         6396 :    END SUBROUTINE get_dW_matrix
     469              : 
     470              : ! **************************************************************************************************
     471              : !> \brief calculates [ab] SHG overlap integrals using precomputed angular-
     472              : !>        dependent part
     473              : !> \param la set of l quantum number on a
     474              : !> \param first_sgfa indexing
     475              : !> \param nshella number of shells for a
     476              : !> \param lb set of l quantum number on b
     477              : !> \param first_sgfb indexing
     478              : !> \param nshellb number of shells for b
     479              : !> \param swork_cont contracted and normalized [s|s] integrals
     480              : !> \param Waux_mat precomputed angular-dependent part
     481              : !> \param sab contracted integral of spherical harmonic Gaussianslm
     482              : ! **************************************************************************************************
     483     20393168 :    SUBROUTINE construct_int_shg_ab(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
     484     20393168 :                                    swork_cont, Waux_mat, sab)
     485              : 
     486              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: la, first_sgfa
     487              :       INTEGER, INTENT(IN)                                :: nshella
     488              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: lb, first_sgfb
     489              :       INTEGER, INTENT(IN)                                :: nshellb
     490              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: swork_cont, Waux_mat
     491              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: sab
     492              : 
     493              :       INTEGER                                            :: fnla, fnlb, fsgfa, fsgfb, ishella, j, &
     494              :                                                             jshellb, labmin, lai, lbj, lnla, lnlb, &
     495              :                                                             lsgfa, lsgfb, mai, mbj
     496              :       REAL(KIND=dp)                                      :: prefac
     497              : 
     498     41010759 :       DO jshellb = 1, nshellb
     499     20617591 :          lbj = lb(jshellb)
     500     20617591 :          fnlb = nsoset(lbj - 1) + 1
     501     20617591 :          lnlb = nsoset(lbj)
     502     20617591 :          fsgfb = first_sgfb(jshellb)
     503     20617591 :          lsgfb = fsgfb + 2*lbj
     504     62348090 :          DO ishella = 1, nshella
     505     21337331 :             lai = la(ishella)
     506     21337331 :             fnla = nsoset(lai - 1) + 1
     507     21337331 :             lnla = nsoset(lai)
     508     21337331 :             fsgfa = first_sgfa(ishella)
     509     21337331 :             lsgfa = fsgfa + 2*lai
     510     21337331 :             labmin = MIN(lai, lbj)
     511     98271251 :             DO mbj = 0, 2*lbj
     512    232316947 :                DO mai = 0, 2*lai
     513    524733240 :                   DO j = 0, labmin
     514    313753624 :                      prefac = swork_cont(lai + lbj - j + 1, ishella, jshellb)
     515              :                      sab(fsgfa + mai, fsgfb + mbj) = sab(fsgfa + mai, fsgfb + mbj) &
     516    468416911 :                                                      + prefac*Waux_mat(j + 1, fnla + mai, fnlb + mbj)
     517              :                   END DO
     518              :                END DO
     519              :             END DO
     520              :          END DO
     521              :       END DO
     522              : 
     523     20393168 :    END SUBROUTINE construct_int_shg_ab
     524              : 
     525              : ! **************************************************************************************************
     526              : !> \brief calculates derivatives of [ab] SHG overlap integrals using precomputed
     527              : !>        angular-dependent part
     528              : !> \param la set of l quantum number on a
     529              : !> \param first_sgfa indexing
     530              : !> \param nshella number of shells for a
     531              : !> \param lb set of l quantum number on b
     532              : !> \param first_sgfb indexing
     533              : !> \param nshellb number of shells for b
     534              : !> \param rab distance vector Ra-Rb
     535              : !> \param swork_cont contracted and normalized [s|s] integrals
     536              : !> \param Waux_mat precomputed angular-dependent part
     537              : !> \param dWaux_mat ...
     538              : !> \param dsab derivative of contracted integral of spherical harmonic Gaussians
     539              : ! **************************************************************************************************
     540        77242 :    SUBROUTINE construct_dev_shg_ab(la, first_sgfa, nshella, lb, first_sgfb, nshellb, rab, &
     541        77242 :                                    swork_cont, Waux_mat, dWaux_mat, dsab)
     542              : 
     543              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: la, first_sgfa
     544              :       INTEGER, INTENT(IN)                                :: nshella
     545              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: lb, first_sgfb
     546              :       INTEGER, INTENT(IN)                                :: nshellb
     547              :       REAL(KIND=dp), INTENT(IN)                          :: rab(3)
     548              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: swork_cont, Waux_mat
     549              :       REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: dWaux_mat
     550              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: dsab
     551              : 
     552              :       INTEGER                                            :: fnla, fnlb, fsgfa, fsgfb, i, ishella, j, &
     553              :                                                             jshellb, labmin, lai, lbj, lnla, lnlb, &
     554              :                                                             lsgfa, lsgfb
     555              :       REAL(KIND=dp)                                      :: dprefac, prefac, rabx2(3)
     556              : 
     557       308968 :       rabx2(:) = 2.0_dp*rab
     558       324154 :       DO jshellb = 1, nshellb
     559       246912 :          lbj = lb(jshellb)
     560       246912 :          fnlb = nsoset(lbj - 1) + 1
     561       246912 :          lnlb = nsoset(lbj)
     562       246912 :          fsgfb = first_sgfb(jshellb)
     563       246912 :          lsgfb = fsgfb + 2*lbj
     564      1108796 :          DO ishella = 1, nshella
     565       784642 :             lai = la(ishella)
     566       784642 :             fnla = nsoset(lai - 1) + 1
     567       784642 :             lnla = nsoset(lai)
     568       784642 :             fsgfa = first_sgfa(ishella)
     569       784642 :             lsgfa = fsgfa + 2*lai
     570       784642 :             labmin = MIN(lai, lbj)
     571      2345125 :             DO j = 0, labmin
     572      1313571 :                prefac = swork_cont(lai + lbj - j + 1, ishella, jshellb)
     573      1313571 :                dprefac = swork_cont(lai + lbj - j + 2, ishella, jshellb) !j+1
     574      6038926 :                DO i = 1, 3
     575              :                   dsab(fsgfa:lsgfa, fsgfb:lsgfb, i) = dsab(fsgfa:lsgfa, fsgfb:lsgfb, i) &
     576              :                                                       + rabx2(i)*dprefac*Waux_mat(j + 1, fnla:lnla, fnlb:lnlb) &
     577     97498992 :                                                       + prefac*dWaux_mat(i, j + 1, fnla:lnla, fnlb:lnlb)
     578              :                END DO
     579              :             END DO
     580              :          END DO
     581              :       END DO
     582              : 
     583        77242 :    END SUBROUTINE construct_dev_shg_ab
     584              : 
     585              : ! **************************************************************************************************
     586              : !> \brief calculates [aba] SHG overlap integrals using precomputed angular-
     587              : !>        dependent part
     588              : !> \param la set of l quantum number on a, orbital basis
     589              : !> \param first_sgfa indexing
     590              : !> \param nshella number of shells for a, orbital basis
     591              : !> \param lb set of l quantum number on b. orbital basis
     592              : !> \param first_sgfb indexing
     593              : !> \param nshellb number of shells for b, orbital basis
     594              : !> \param lca of l quantum number on a, aux basis
     595              : !> \param first_sgfca indexing
     596              : !> \param nshellca  number of shells for a, aux basis
     597              : !> \param cg_coeff Clebsch-Gordon coefficients
     598              : !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
     599              : !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
     600              : !> \param swork_cont contracted and normalized [s|ra^n|s] integrals
     601              : !> \param Waux_mat precomputed angular-dependent part
     602              : !> \param saba contracted overlap [aba] of spherical harmonic Gaussians
     603              : ! **************************************************************************************************
     604        40040 :    SUBROUTINE construct_overlap_shg_aba(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
     605        40040 :                                         lca, first_sgfca, nshellca, cg_coeff, cg_none0_list, &
     606        40040 :                                         ncg_none0, swork_cont, Waux_mat, saba)
     607              : 
     608              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: la, first_sgfa
     609              :       INTEGER, INTENT(IN)                                :: nshella
     610              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: lb, first_sgfb
     611              :       INTEGER, INTENT(IN)                                :: nshellb
     612              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: lca, first_sgfca
     613              :       INTEGER, INTENT(IN)                                :: nshellca
     614              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: cg_coeff
     615              :       INTEGER, DIMENSION(:, :, :), INTENT(IN)            :: cg_none0_list
     616              :       INTEGER, DIMENSION(:, :), INTENT(IN)               :: ncg_none0
     617              :       REAL(KIND=dp), DIMENSION(:, 0:, :, :, :), &
     618              :          INTENT(IN)                                      :: swork_cont
     619              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: Waux_mat
     620              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: saba
     621              : 
     622              :       INTEGER :: ia, il, ilist, ishella, isoa1, isoa2, isoaa, j, jb, jshellb, ka, kshella, laa, &
     623              :          labmin, lai, lak, lbj, maa, mai, mak, mbj, nla, nlb, sgfa, sgfb, sgfca
     624              :       REAL(KIND=dp)                                      :: prefac, stemp
     625              : 
     626       143017 :       DO kshella = 1, nshellca
     627       102977 :          lak = lca(kshella)
     628       102977 :          sgfca = first_sgfca(kshella)
     629       102977 :          ka = sgfca + lak
     630       599463 :          DO jshellb = 1, nshellb
     631       456446 :             lbj = lb(jshellb)
     632       456446 :             nlb = nsoset(lbj - 1) + lbj + 1
     633       456446 :             sgfb = first_sgfb(jshellb)
     634       456446 :             jb = sgfb + lbj
     635      2711645 :             DO ishella = 1, nshella
     636      2152222 :                lai = la(ishella)
     637      2152222 :                sgfa = first_sgfa(ishella)
     638      2152222 :                ia = sgfa + lai
     639      8069462 :                DO mai = -lai, lai, 1
     640     26951258 :                   DO mak = -lak, lak, 1
     641     19338242 :                      isoa1 = indso_inv(lai, mai)
     642     19338242 :                      isoa2 = indso_inv(lak, mak)
     643     73280668 :                      DO mbj = -lbj, lbj, 1
     644    167368307 :                         DO ilist = 1, ncg_none0(isoa1, isoa2)
     645     99548433 :                            isoaa = cg_none0_list(isoa1, isoa2, ilist)
     646     99548433 :                            laa = indso(1, isoaa)
     647     99548433 :                            maa = indso(2, isoaa)
     648     99548433 :                            nla = nsoset(laa - 1) + laa + 1
     649     99548433 :                            labmin = MIN(laa, lbj)
     650     99548433 :                            il = INT((lai + lak - laa)/2)
     651     99548433 :                            stemp = 0.0_dp
     652    306983753 :                            DO j = 0, labmin
     653    207435320 :                               prefac = swork_cont(laa + lbj - j + 1, il, ishella, jshellb, kshella)
     654    306983753 :                               stemp = stemp + prefac*Waux_mat(j + 1, nla + maa, nlb + mbj)
     655              :                            END DO
     656    148030065 :                        saba(ia + mai, jb + mbj, ka + mak) = saba(ia + mai, jb + mbj, ka + mak) + cg_coeff(isoa1, isoa2, isoaa)*stemp
     657              :                         END DO
     658              :                      END DO
     659              :                   END DO
     660              :                END DO
     661              :             END DO
     662              :          END DO
     663              :       END DO
     664              : 
     665        40040 :    END SUBROUTINE construct_overlap_shg_aba
     666              : 
     667              : ! **************************************************************************************************
     668              : !> \brief calculates derivatives of [aba] SHG overlap integrals using
     669              : !>        precomputed angular-dependent part
     670              : !> \param la set of l quantum number on a, orbital basis
     671              : !> \param first_sgfa indexing
     672              : !> \param nshella number of shells for a, orbital basis
     673              : !> \param lb set of l quantum number on b. orbital basis
     674              : !> \param first_sgfb indexing
     675              : !> \param nshellb number of shells for b, orbital basis
     676              : !> \param lca of l quantum number on a, aux basis
     677              : !> \param first_sgfca indexing
     678              : !> \param nshellca  number of shells for a, aux basis
     679              : !> \param cg_coeff Clebsch-Gordon coefficients
     680              : !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
     681              : !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
     682              : !> \param rab distance vector Ra-Rb
     683              : !> \param swork_cont contracted and normalized [s|ra^n|s] integrals
     684              : !> \param Waux_mat precomputed angular-dependent part
     685              : !> \param dWaux_mat derivatives of precomputed angular-dependent part
     686              : !> \param dsaba derivative of contracted overlap [aba] of spherical harmonic
     687              : !>              Gaussians
     688              : ! **************************************************************************************************
     689        27388 :    SUBROUTINE dev_overlap_shg_aba(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
     690        54776 :                                   lca, first_sgfca, nshellca, cg_coeff, cg_none0_list, &
     691        27388 :                                   ncg_none0, rab, swork_cont, Waux_mat, dWaux_mat, dsaba)
     692              : 
     693              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: la, first_sgfa
     694              :       INTEGER, INTENT(IN)                                :: nshella
     695              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: lb, first_sgfb
     696              :       INTEGER, INTENT(IN)                                :: nshellb
     697              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: lca, first_sgfca
     698              :       INTEGER, INTENT(IN)                                :: nshellca
     699              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: cg_coeff
     700              :       INTEGER, DIMENSION(:, :, :), INTENT(IN)            :: cg_none0_list
     701              :       INTEGER, DIMENSION(:, :), INTENT(IN)               :: ncg_none0
     702              :       REAL(KIND=dp), INTENT(IN)                          :: rab(3)
     703              :       REAL(KIND=dp), DIMENSION(:, 0:, :, :, :), &
     704              :          INTENT(IN)                                      :: swork_cont
     705              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: Waux_mat
     706              :       REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: dWaux_mat
     707              :       REAL(KIND=dp), DIMENSION(:, :, :, :), &
     708              :          INTENT(INOUT)                                   :: dsaba
     709              : 
     710              :       INTEGER :: i, ia, il, ilist, ishella, isoa1, isoa2, isoaa, j, jb, jshellb, ka, kshella, laa, &
     711              :          labmin, lai, lak, lbj, maa, mai, mak, mbj, nla, nlb, sgfa, sgfb, sgfca
     712              :       REAL(KIND=dp)                                      :: dprefac, dtemp(3), prefac, rabx2(3)
     713              : 
     714       109552 :       rabx2(:) = 2.0_dp*rab
     715              : 
     716        95801 :       DO kshella = 1, nshellca
     717        68413 :          lak = lca(kshella)
     718        68413 :          sgfca = first_sgfca(kshella)
     719        68413 :          ka = sgfca + lak
     720       428276 :          DO jshellb = 1, nshellb
     721       332475 :             lbj = lb(jshellb)
     722       332475 :             nlb = nsoset(lbj - 1) + lbj + 1
     723       332475 :             sgfb = first_sgfb(jshellb)
     724       332475 :             jb = sgfb + lbj
     725      2038117 :             DO ishella = 1, nshella
     726      1637229 :                lai = la(ishella)
     727      1637229 :                sgfa = first_sgfa(ishella)
     728      1637229 :                ia = sgfa + lai
     729      6211191 :                DO mai = -lai, lai, 1
     730     20427013 :                   DO mak = -lak, lak, 1
     731     14548297 :                      isoa1 = indso_inv(lai, mai)
     732     14548297 :                      isoa2 = indso_inv(lak, mak)
     733     56467669 :                      DO mbj = -lbj, lbj, 1
     734    128860480 :                         DO ilist = 1, ncg_none0(isoa1, isoa2)
     735     76634298 :                            isoaa = cg_none0_list(isoa1, isoa2, ilist)
     736     76634298 :                            laa = indso(1, isoaa)
     737     76634298 :                            maa = indso(2, isoaa)
     738     76634298 :                            nla = nsoset(laa - 1) + laa + 1
     739     76634298 :                            labmin = MIN(laa, lbj)
     740     76634298 :                            il = (lai + lak - laa)/2 ! lai+lak-laa always even
     741     76634298 :                            dtemp = 0.0_dp
     742    238324798 :                            DO j = 0, labmin
     743    161690500 :                               prefac = swork_cont(laa + lbj - j + 1, il, ishella, jshellb, kshella)
     744    161690500 :                               dprefac = swork_cont(laa + lbj - j + 2, il, ishella, jshellb, kshella)
     745    723396298 :                               DO i = 1, 3
     746              :                                  dtemp(i) = dtemp(i) + rabx2(i)*dprefac*Waux_mat(j + 1, nla + maa, nlb + mbj) &
     747    646762000 :                                             + prefac*dWaux_mat(i, j + 1, nla + maa, nlb + mbj)
     748              :                               END DO
     749              :                            END DO
     750    344215077 :                            DO i = 1, 3
     751              :                               dsaba(ia + mai, jb + mbj, ka + mak, i) = dsaba(ia + mai, jb + mbj, ka + mak, i) &
     752    306537192 :                                                                        + cg_coeff(isoa1, isoa2, isoaa)*dtemp(i)
     753              :                            END DO
     754              :                         END DO
     755              :                      END DO
     756              :                   END DO
     757              :                END DO
     758              :             END DO
     759              :          END DO
     760              :       END DO
     761              : 
     762        27388 :    END SUBROUTINE dev_overlap_shg_aba
     763              : 
     764              : ! **************************************************************************************************
     765              : !> \brief calculates [abb] SHG overlap integrals using precomputed angular-
     766              : !>        dependent part
     767              : !> \param la set of l quantum number on a, orbital basis
     768              : !> \param first_sgfa indexing
     769              : !> \param nshella number of shells for a, orbital basis
     770              : !> \param lb set of l quantum number on b. orbital basis
     771              : !> \param first_sgfb indexing
     772              : !> \param nshellb number of shells for b, orbital basis
     773              : !> \param lcb l quantum number on b, aux basis
     774              : !> \param first_sgfcb indexing
     775              : !> \param nshellcb number of shells for b, aux basis
     776              : !> \param cg_coeff Clebsch-Gordon coefficients
     777              : !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
     778              : !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
     779              : !> \param swork_cont contracted and normalized [s|rb^n|s] integrals
     780              : !> \param Waux_mat precomputed angular-dependent part
     781              : !> \param sabb contracted overlap [abb] of spherical harmonic Gaussians
     782              : ! **************************************************************************************************
     783         3277 :    SUBROUTINE construct_overlap_shg_abb(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
     784         3277 :                                         lcb, first_sgfcb, nshellcb, cg_coeff, cg_none0_list, &
     785         3277 :                                         ncg_none0, swork_cont, Waux_mat, sabb)
     786              : 
     787              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: la, first_sgfa
     788              :       INTEGER, INTENT(IN)                                :: nshella
     789              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: lb, first_sgfb
     790              :       INTEGER, INTENT(IN)                                :: nshellb
     791              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: lcb, first_sgfcb
     792              :       INTEGER, INTENT(IN)                                :: nshellcb
     793              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: cg_coeff
     794              :       INTEGER, DIMENSION(:, :, :), INTENT(IN)            :: cg_none0_list
     795              :       INTEGER, DIMENSION(:, :), INTENT(IN)               :: ncg_none0
     796              :       REAL(KIND=dp), DIMENSION(:, 0:, :, :, :), &
     797              :          INTENT(IN)                                      :: swork_cont
     798              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: Waux_mat
     799              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: sabb
     800              : 
     801              :       INTEGER :: ia, il, ilist, ishella, isob1, isob2, isobb, j, jb, jshellb, kb, kshellb, labmin, &
     802              :          lai, lbb, lbj, lbk, mai, mbb, mbj, mbk, nla, nlb, sgfa, sgfb, sgfcb
     803              :       REAL(KIND=dp)                                      :: prefac, stemp, tsign
     804              : 
     805        14952 :       DO kshellb = 1, nshellcb
     806        11675 :          lbk = lcb(kshellb)
     807        11675 :          sgfcb = first_sgfcb(kshellb)
     808        11675 :          kb = sgfcb + lbk
     809        62545 :          DO jshellb = 1, nshellb
     810        47593 :             lbj = lb(jshellb)
     811        47593 :             sgfb = first_sgfb(jshellb)
     812        47593 :             jb = sgfb + lbj
     813       212112 :             DO ishella = 1, nshella
     814       152844 :                lai = la(ishella)
     815       152844 :                nla = nsoset(lai - 1) + lai + 1
     816       152844 :                sgfa = first_sgfa(ishella)
     817       152844 :                ia = sgfa + lai
     818       565795 :                DO mbj = -lbj, lbj, 1
     819      2165118 :                   DO mbk = -lbk, lbk, 1
     820      1646916 :                      isob1 = indso_inv(lbj, mbj)
     821      1646916 :                      isob2 = indso_inv(lbk, mbk)
     822      4963078 :                      DO mai = -lai, lai, 1
     823     11318995 :                         DO ilist = 1, ncg_none0(isob1, isob2)
     824      6721275 :                            isobb = cg_none0_list(isob1, isob2, ilist)
     825      6721275 :                            lbb = indso(1, isobb)
     826      6721275 :                            mbb = indso(2, isobb)
     827      6721275 :                            nlb = nsoset(lbb - 1) + lbb + 1
     828              :                            ! tsgin: because we take the transpose of auxmat (calculated for (la,lb))
     829      6721275 :                            tsign = 1.0_dp
     830      6721275 :                            IF (MODULO(lbb - lai, 2) /= 0) tsign = -1.0_dp
     831      6721275 :                            labmin = MIN(lai, lbb)
     832      6721275 :                            il = INT((lbj + lbk - lbb)/2)
     833      6721275 :                            stemp = 0.0_dp
     834     18145266 :                            DO j = 0, labmin
     835     11423991 :                               prefac = swork_cont(lai + lbb - j + 1, il, ishella, jshellb, kshellb)
     836     18145266 :                               stemp = stemp + prefac*Waux_mat(j + 1, nlb + mbb, nla + mai)
     837              :                            END DO
     838      9672079 :                  sabb(ia + mai, jb + mbj, kb + mbk) = sabb(ia + mai, jb + mbj, kb + mbk) + tsign*cg_coeff(isob1, isob2, isobb)*stemp
     839              :                         END DO
     840              :                      END DO
     841              :                   END DO
     842              :                END DO
     843              :             END DO
     844              :          END DO
     845              :       END DO
     846              : 
     847         3277 :    END SUBROUTINE construct_overlap_shg_abb
     848              : 
     849              : ! **************************************************************************************************
     850              : !> \brief calculates derivatives of [abb] SHG overlap integrals using
     851              : !>        precomputed angular-dependent part
     852              : !> \param la set of l quantum number on a, orbital basis
     853              : !> \param first_sgfa indexing
     854              : !> \param nshella number of shells for a, orbital basis
     855              : !> \param lb set of l quantum number on b. orbital basis
     856              : !> \param first_sgfb indexing
     857              : !> \param nshellb number of shells for b, orbital basis
     858              : !> \param lcb l quantum number on b, aux basis
     859              : !> \param first_sgfcb indexing
     860              : !> \param nshellcb number of shells for b, aux basis
     861              : !> \param cg_coeff Clebsch-Gordon coefficients
     862              : !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
     863              : !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
     864              : !> \param rab distance vector Ra-Rb
     865              : !> \param swork_cont contracted and normalized [s|rb^n|s] integrals
     866              : !> \param Waux_mat precomputed angular-dependent part
     867              : !> \param dWaux_mat derivatives of precomputed angular-dependent part
     868              : !> \param dsabb derivative of contracted overlap [abb] of spherical harmonic
     869              : !>        Gaussians
     870              : ! **************************************************************************************************
     871          967 :    SUBROUTINE dev_overlap_shg_abb(la, first_sgfa, nshella, lb, first_sgfb, nshellb, &
     872         1934 :                                   lcb, first_sgfcb, nshellcb, cg_coeff, cg_none0_list, &
     873          967 :                                   ncg_none0, rab, swork_cont, Waux_mat, dWaux_mat, dsabb)
     874              : 
     875              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: la, first_sgfa
     876              :       INTEGER, INTENT(IN)                                :: nshella
     877              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: lb, first_sgfb
     878              :       INTEGER, INTENT(IN)                                :: nshellb
     879              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: lcb, first_sgfcb
     880              :       INTEGER, INTENT(IN)                                :: nshellcb
     881              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: cg_coeff
     882              :       INTEGER, DIMENSION(:, :, :), INTENT(IN)            :: cg_none0_list
     883              :       INTEGER, DIMENSION(:, :), INTENT(IN)               :: ncg_none0
     884              :       REAL(KIND=dp), INTENT(IN)                          :: rab(3)
     885              :       REAL(KIND=dp), DIMENSION(:, 0:, :, :, :), &
     886              :          INTENT(IN)                                      :: swork_cont
     887              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: Waux_mat
     888              :       REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: dWaux_mat
     889              :       REAL(KIND=dp), DIMENSION(:, :, :, :), &
     890              :          INTENT(INOUT)                                   :: dsabb
     891              : 
     892              :       INTEGER :: i, ia, il, ilist, ishella, isob1, isob2, isobb, j, jb, jshellb, kb, kshellb, &
     893              :          labmin, lai, lbb, lbj, lbk, mai, mbb, mbj, mbk, nla, nlb, sgfa, sgfb, sgfcb
     894              :       REAL(KIND=dp)                                      :: dprefac, dtemp(3), prefac, rabx2(3), &
     895              :                                                             tsign
     896              : 
     897         3868 :       rabx2(:) = 2.0_dp*rab
     898              : 
     899         3794 :       DO kshellb = 1, nshellcb
     900         2827 :          lbk = lcb(kshellb)
     901         2827 :          sgfcb = first_sgfcb(kshellb)
     902         2827 :          kb = sgfcb + lbk
     903        12463 :          DO jshellb = 1, nshellb
     904         8669 :             lbj = lb(jshellb)
     905         8669 :             sgfb = first_sgfb(jshellb)
     906         8669 :             jb = sgfb + lbj
     907        34899 :             DO ishella = 1, nshella
     908        23403 :                lai = la(ishella)
     909        23403 :                nla = nsoset(lai - 1) + lai + 1
     910        23403 :                sgfa = first_sgfa(ishella)
     911        23403 :                ia = sgfa + lai
     912        92645 :                DO mbj = -lbj, lbj, 1
     913       359327 :                   DO mbk = -lbk, lbk, 1
     914       275351 :                      isob1 = indso_inv(lbj, mbj)
     915       275351 :                      isob2 = indso_inv(lbk, mbk)
     916       805479 :                      DO mai = -lai, lai, 1
     917      1982969 :                         DO ilist = 1, ncg_none0(isob1, isob2)
     918      1238063 :                            isobb = cg_none0_list(isob1, isob2, ilist)
     919      1238063 :                            lbb = indso(1, isobb)
     920      1238063 :                            mbb = indso(2, isobb)
     921      1238063 :                            nlb = nsoset(lbb - 1) + lbb + 1
     922              :                            ! tsgin: because we take the transpose of auxmat (calculated for (la,lb))
     923      1238063 :                            tsign = 1.0_dp
     924      1238063 :                            IF (MODULO(lbb - lai, 2) /= 0) tsign = -1.0_dp
     925      1238063 :                            labmin = MIN(lai, lbb)
     926      1238063 :                            il = (lbj + lbk - lbb)/2
     927      1238063 :                            dtemp = 0.0_dp
     928      3490537 :                            DO j = 0, labmin
     929      2252474 :                               prefac = swork_cont(lai + lbb - j + 1, il, ishella, jshellb, kshellb)
     930      2252474 :                               dprefac = swork_cont(lai + lbb - j + 2, il, ishella, jshellb, kshellb)
     931     10247959 :                               DO i = 1, 3
     932              :                                  dtemp(i) = dtemp(i) + rabx2(i)*dprefac*Waux_mat(j + 1, nlb + mbb, nla + mai) &
     933      9009896 :                                             + prefac*dWaux_mat(i, j + 1, nlb + mbb, nla + mai)
     934              :                               END DO
     935              :                            END DO
     936      5421807 :                            DO i = 1, 3
     937              :                               dsabb(ia + mai, jb + mbj, kb + mbk, i) = dsabb(ia + mai, jb + mbj, kb + mbk, i) &
     938      4952252 :                                                                        + tsign*cg_coeff(isob1, isob2, isobb)*dtemp(i)
     939              :                            END DO
     940              :                         END DO
     941              :                      END DO
     942              :                   END DO
     943              :                END DO
     944              :             END DO
     945              :          END DO
     946              :       END DO
     947              : 
     948          967 :    END SUBROUTINE dev_overlap_shg_abb
     949              : 
     950              : END MODULE construct_shg
     951              : 
        

Generated by: LCOV version 2.0-1