From 889ab9ec78c7cafaa8381ac5b2fa7c17264aaae3 Mon Sep 17 00:00:00 2001 From: Dorothea Golze Date: Fri, 5 Aug 2016 21:04:08 +0000 Subject: [PATCH] SHG integrals: speed-up of [aba] and [abb] integrals + inclusion in libtest svn-origin-rev: 17197 --- src/aobasis/ai_s_contract_shg.F | 466 +++++++------- src/generic_shg_integrals.F | 434 +++++++------ src/generic_shg_integrals_init.F | 351 +++++++++++ src/lri_environment_init.F | 210 +------ src/lri_environment_methods.F | 40 +- src/lri_environment_types.F | 7 + src/shg_integrals_test.F | 802 ++++++++++++++++++------ tests/LIBTEST/TEST_FILES | 1 + tests/LIBTEST/test_shg_integrals_01.inp | 11 +- tests/LIBTEST/test_shg_integrals_02.inp | 72 +++ 10 files changed, 1551 insertions(+), 843 deletions(-) create mode 100644 src/generic_shg_integrals_init.F create mode 100644 tests/LIBTEST/test_shg_integrals_02.inp diff --git a/src/aobasis/ai_s_contract_shg.F b/src/aobasis/ai_s_contract_shg.F index 68f068ef93..0f9b5b0f92 100644 --- a/src/aobasis/ai_s_contract_shg.F +++ b/src/aobasis/ai_s_contract_shg.F @@ -20,8 +20,8 @@ MODULE ai_s_contract_shg USE gamma, ONLY: fgamma=>fgamma_0 USE kinds, ONLY: dp - USE mathconstants, ONLY: fac,& - ifac,& + USE mathconstants, ONLY: dfac,& + fac,& pi #include "../base/base_uses.f90" @@ -32,7 +32,7 @@ MODULE ai_s_contract_shg CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_s_contract_shg' ! *** Public subroutines *** - PUBLIC :: s_overlap_ab, s_overlap_abx, s_coulomb_ab, s_verf_ab, s_verfc_ab, & + PUBLIC :: s_overlap_ab, s_overlap_abb, s_coulomb_ab, s_verf_ab, s_verfc_ab, & s_vgauss_ab, s_gauss_ab PUBLIC :: contract_sint_ab_clow, contract_sint_ab_chigh, contract_s_overlap_aba, & @@ -102,44 +102,43 @@ CONTAINS ! ************************************************************************************************** !> \brief calculates [s|ra^n|s] integrals for [aba] and the [s|rb^n|s] !> integrals for [abb] -!> abx ... x stands either for a and b respectively !> \param la_max maximal l quantum number on a, orbital basis !> \param npgfa number of primitive Gaussian on a, orbital basis !> \param zeta set of exponents on a, orbital basis !> \param lb_max maximal l quantum number on b, orbital basis !> \param npgfb number of primitive Gaussian on a, orbital basis !> \param zetb set of exponents on b, orbital basis -!> \param lx_max maximal l quantum number of aux basis on x (a or b) -!> \param npgfx number of primitive Gaussian on x, aux basis -!> \param zetx set of exponents on x, aux basis +!> \param lcb_max maximal l quantum number of aux basis on b +!> \param npgfcb number of primitive Gaussian on b. aux basis +!> \param zetcb set of exponents on b, aux basis !> \param rab distance vector between a and b !> \param s uncontracted [s|r^n|s] integrals !> \param calculate_forces ... -!> \param calc_aba flag if [aba] or [abb] is calculated ! ************************************************************************************************** - SUBROUTINE s_overlap_abx(la_max, npgfa, zeta, lb_max, npgfb, zetb, lx_max, npgfx, zetx, & - rab, s, calculate_forces, calc_aba) + SUBROUTINE s_overlap_abb(la_max, npgfa, zeta, lb_max, npgfb, zetb, lcb_max, npgfcb, zetcb, & + rab, s, calculate_forces) INTEGER, INTENT(IN) :: la_max, npgfa REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta INTEGER, INTENT(IN) :: lb_max, npgfb REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb - INTEGER, INTENT(IN) :: lx_max, npgfx - REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetx + INTEGER, INTENT(IN) :: lcb_max, npgfcb + REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetcb REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab - REAL(KIND=dp), DIMENSION(:, :, :, 0:, :), & + REAL(KIND=dp), DIMENSION(:, :, :, :, :), & INTENT(INOUT) :: s - LOGICAL, INTENT(IN) :: calculate_forces, calc_aba + LOGICAL, INTENT(IN) :: calculate_forces - CHARACTER(len=*), PARAMETER :: routineN = 's_overlap_abx', routineP = moduleN//':'//routineN + CHARACTER(len=*), PARAMETER :: routineN = 's_overlap_abb', routineP = moduleN//':'//routineN - INTEGER :: handle, ids, il, ipgfa, jpgfb, kpgfx, & - laa_max, lbb_max, lmax, n, ndev, nds, & - nl - REAL(KIND=dp) :: a, b, exp_rab2, k, pfac, prefac, rab2, & - sqrt_pi3, sqrt_zet, sr_int, temp, x, & - xhi, zet + INTEGER :: handle, i, ids, il, ipgfa, j, jpgfb, & + kpgfb, lbb_max, lmax, n, ndev, nds, & + nfac, nl + REAL(KIND=dp) :: a, b, dfsr_int, exp_rab2, k, pfac, & + prefac, rab2, sqrt_pi3, sqrt_zet, & + sr_int, temp, xhi, zet REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dsr_int, dtemp + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: coeff_srs CALL timeset(routineN, handle) @@ -148,43 +147,31 @@ CONTAINS ndev = 0 IF (calculate_forces) ndev = 1 - IF (calc_aba) THEN - laa_max = la_max+lx_max - nl = INT(laa_max/2) - IF (la_max == 0 .OR. lx_max == 0) nl = 0 - lmax = laa_max+lb_max - ELSE - lbb_max = lb_max+lx_max - nl = INT(lbb_max/2) - IF (lb_max == 0 .OR. lx_max == 0) nl = 0 - lmax = la_max+lbb_max - ENDIF + lbb_max = lb_max+lcb_max + nl = INT(lbb_max/2) + IF (lb_max == 0 .OR. lcb_max == 0) nl = 0 + lmax = la_max+lbb_max ALLOCATE (dtemp(nl+1), dsr_int(nl+1)) + ALLOCATE (coeff_srs(nl+1, nl+1, nl+1)) + IF (nl > 5) CALL get_prefac_sabb(nl, coeff_srs) sqrt_pi3 = SQRT(pi**3) ! Loops over all pairs of primitive Gaussian-type functions - DO ipgfa = 1, npgfa + DO kpgfb = 1, npgfcb DO jpgfb = 1, npgfb - DO kpgfx = 1, npgfx + DO ipgfa = 1, npgfa - ! Calculate some prefactors - IF (calc_aba) THEN - a = zeta(ipgfa)+zetx(kpgfx) - b = zetb(jpgfb) - x = b - ELSE - a = zeta(ipgfa) - b = zetb(jpgfb)+zetx(kpgfx) - x = a - ENDIF + !Calculate some prefactors + a = zeta(ipgfa) + b = zetb(jpgfb)+zetcb(kpgfb) zet = a+b xhi = a*b/zet exp_rab2 = EXP(-xhi*rab2) - pfac = x**2/zet + pfac = a**2/zet sqrt_zet = SQRT(zet) DO il = 0, nl @@ -192,36 +179,36 @@ CONTAINS SELECT CASE (il) CASE (0) ! [s|s] integral - s(ipgfa, jpgfb, kpgfx, il, 1) = (pi/zet)**(1.5_dp)*exp_rab2 + s(ipgfa, jpgfb, kpgfb, il+1, 1) = (pi/zet)**(1.5_dp)*exp_rab2 DO ids = 2, nds n = ids-1 - s(ipgfa, jpgfb, kpgfx, il, ids) = (-xhi)**n*s(ipgfa, jpgfb, kpgfx, il, 1) + s(ipgfa, jpgfb, kpgfb, il+1, ids) = (-xhi)**n*s(ipgfa, jpgfb, kpgfb, il+1, 1) ENDDO CASE (1) ![s|r^2|s] integral sr_int = sqrt_pi3/sqrt_zet**5*(3.0_dp+2.0_dp*pfac*rab2)/2.0_dp - s(ipgfa, jpgfb, kpgfx, il, 1) = exp_rab2*sr_int - k = sqrt_pi3*x**2/sqrt_zet**7 + s(ipgfa, jpgfb, kpgfb, il+1, 1) = exp_rab2*sr_int + k = sqrt_pi3*a**2/sqrt_zet**7 DO ids = 2, nds n = ids-1 - s(ipgfa, jpgfb, kpgfx, il, ids) = (-xhi)**n*exp_rab2*sr_int & - +n*(-xhi)**(n-1)*k*exp_rab2 + s(ipgfa, jpgfb, kpgfb, il+1, ids) = (-xhi)**n*exp_rab2*sr_int & + +n*(-xhi)**(n-1)*k*exp_rab2 ENDDO CASE (2) ![s|r^4|s] integral prefac = sqrt_pi3/4.0_dp/sqrt_zet**7 temp = 15.0_dp+20.0_dp*pfac*rab2+4.0_dp*(pfac*rab2)**2 sr_int = prefac*temp - s(ipgfa, jpgfb, kpgfx, il, 1) = exp_rab2*sr_int + s(ipgfa, jpgfb, kpgfb, il+1, 1) = exp_rab2*sr_int !** derivatives - k = sqrt_pi3*x**4/sqrt_zet**11 + k = sqrt_pi3*a**4/sqrt_zet**11 dsr_int(1) = prefac*(20.0_dp*pfac+8.0_dp*pfac**2*rab2) DO ids = 2, nds n = ids-1 dtemp(1) = (-xhi)**n*exp_rab2*sr_int dtemp(2) = n*(-xhi)**(n-1)*exp_rab2*dsr_int(1) dtemp(3) = (n**2-n)*(-xhi)**(n-2)*k*exp_rab2 - s(ipgfa, jpgfb, kpgfx, il, ids) = dtemp(1)+dtemp(2)+dtemp(3) + s(ipgfa, jpgfb, kpgfb, il+1, ids) = dtemp(1)+dtemp(2)+dtemp(3) ENDDO CASE (3) ![s|r^6|s] integral @@ -229,9 +216,9 @@ CONTAINS temp = 105.0_dp+210.0_dp*pfac*rab2 temp = temp+84.0_dp*(pfac*rab2)**2+8.0_dp*(pfac*rab2)**3 sr_int = prefac*temp - s(ipgfa, jpgfb, kpgfx, il, 1) = exp_rab2*sr_int + s(ipgfa, jpgfb, kpgfb, il+1, 1) = exp_rab2*sr_int !** derivatives - k = sqrt_pi3*x**6/sqrt_zet**15 + k = sqrt_pi3*a**6/sqrt_zet**15 dsr_int(1) = prefac*(210.0_dp*pfac+168.0_dp*pfac**2*rab2 & +24.0_dp*pfac**3*rab2**2) dsr_int(2) = prefac*(168.0_dp*pfac**2+48.0_dp*pfac**3*rab2) @@ -242,8 +229,8 @@ CONTAINS dtemp(3) = REAL(n**2-n, dp)/2.0_dp*(-xhi)**(n-2) & *exp_rab2*dsr_int(2) dtemp(4) = REAL(n*(n-1)*(n-2), dp)*(-xhi)**(n-3)*k*exp_rab2 - s(ipgfa, jpgfb, kpgfx, il, ids) = dtemp(1)+dtemp(2) & - +dtemp(3)+dtemp(4) + s(ipgfa, jpgfb, kpgfb, il+1, ids) = dtemp(1)+dtemp(2) & + +dtemp(3)+dtemp(4) ENDDO CASE (4) ![s|r^8|s] integral @@ -251,9 +238,9 @@ CONTAINS temp = 945.0_dp+2520.0_dp*pfac*rab2+1512.0_dp*(pfac*rab2)**2 temp = temp+288.0_dp*(pfac*rab2)**3+16.0_dp*(pfac*rab2)**4 sr_int = prefac*temp - s(ipgfa, jpgfb, kpgfx, il, 1) = exp_rab2*sr_int + s(ipgfa, jpgfb, kpgfb, il+1, 1) = exp_rab2*sr_int !** derivatives - k = sqrt_pi3*x**8/sqrt_zet**19 + k = sqrt_pi3*a**8/sqrt_zet**19 dsr_int(1) = 2520.0_dp*pfac+3024.0_dp*pfac**2*rab2 dsr_int(1) = dsr_int(1)+864.0_dp*pfac**3*rab2**2 & +64.0_dp*pfac**4*rab2**3 @@ -273,8 +260,8 @@ CONTAINS *exp_rab2*dsr_int(3) dtemp(5) = REAL(n*(n-1)*(n-2)*(n-3), dp)*(-xhi)**(n-4) & *k*exp_rab2 - s(ipgfa, jpgfb, kpgfx, il, ids) = dtemp(1)+dtemp(2)+dtemp(3) & - +dtemp(4)+dtemp(5) + s(ipgfa, jpgfb, kpgfb, il+1, ids) = dtemp(1)+dtemp(2)+dtemp(3) & + +dtemp(4)+dtemp(5) ENDDO CASE (5) ![s|r^10|s] integral @@ -283,9 +270,9 @@ CONTAINS temp = temp+27720.0_dp*(pfac*rab2)**2+7920.0_dp*(pfac*rab2)**3 temp = temp+880.0_dp*(pfac*rab2)**4+32.0_dp*(pfac*rab2)**5 sr_int = prefac*temp - s(ipgfa, jpgfb, kpgfx, il, 1) = exp_rab2*sr_int + s(ipgfa, jpgfb, kpgfb, il+1, 1) = exp_rab2*sr_int !** derivatives - k = sqrt_pi3*x**10/sqrt_zet**23 + k = sqrt_pi3*a**10/sqrt_zet**23 dsr_int(1) = 34650.0_dp*pfac+55440.0_dp*pfac**2*rab2 dsr_int(1) = dsr_int(1)+23760.0_dp*pfac**3*rab2**2 dsr_int(1) = dsr_int(1)+3520.0_dp*pfac**4*rab2**3 @@ -312,13 +299,33 @@ CONTAINS *exp_rab2*dsr_int(4) dtemp(6) = REAL(n*(n-1)*(n-2)*(n-3)*(n-4), dp)*(-xhi)**(n-5) & *k*exp_rab2 - s(ipgfa, jpgfb, kpgfx, il, ids) = dtemp(1)+dtemp(2)+dtemp(3) & - +dtemp(4)+dtemp(5)+dtemp(6) + s(ipgfa, jpgfb, kpgfb, il+1, ids) = dtemp(1)+dtemp(2)+dtemp(3) & + +dtemp(4)+dtemp(5)+dtemp(6) ENDDO CASE DEFAULT - CALL cp_abort(__LOCATION__, & - "SHG integrals not implemented when sum of l quantum"// & - " number of orbital and ri basis larger than 11") + !*** general formula; factor 1.5-2 slower than explicit expressions + prefac = exp_rab2/sqrt_zet**(2*il+3) + sr_int = 0.0_dp + DO i = 0, il + sr_int = sr_int+coeff_srs(i+1, 1, il+1)*(pfac)**i*rab2**i + ENDDO + s(ipgfa, jpgfb, kpgfb, il+1, 1) = prefac*sr_int + nds = lmax-2*il+ndev+1 + DO ids = 2, nds + n = ids-1 + nfac = 1 + dfsr_int = (-xhi)**n*sr_int + DO j = 1, il + temp = 0.0_dp + DO i = j, il + temp = temp+coeff_srs(i+1, j+1, il+1)*(pfac)**i*rab2**(i-j) + ENDDO + nfac = nfac*(n-j+1) + dfsr_int = dfsr_int+temp*REAL(nfac, dp)/fac(j)*(-xhi)**(n-j) + ENDDO + s(ipgfa, jpgfb, kpgfb, il+1, ids) = prefac*dfsr_int + ENDDO + END SELECT ENDDO @@ -328,10 +335,39 @@ CONTAINS END DO DEALLOCATE (dtemp, dsr_int) + DEALLOCATE (coeff_srs) CALL timestop(handle) - END SUBROUTINE s_overlap_abx + END SUBROUTINE s_overlap_abb + +! ************************************************************************************************** +!> \brief prefactor for the general formula to calculate the (0a|0b|0b) overlap integrals +!> \param nl ... +!> \param prefac ... +! ************************************************************************************************** + SUBROUTINE get_prefac_sabb(nl, prefac) + INTEGER, INTENT(IN) :: nl + REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: prefac + + CHARACTER(len=*), PARAMETER :: routineN = 'get_prefac_sabb', & + routineP = moduleN//':'//routineN + + INTEGER :: il, j, k + REAL(KIND=dp) :: sqrt_pi3, temp + + sqrt_pi3 = SQRT(pi**3) + + DO il = 0, nl + temp = dfac(2*il+1)*sqrt_pi3*fac(il)/2.0_dp**il + DO j = 0, il + DO k = j, il + prefac(k+1, j+1, il+1) = temp*2.0_dp**k/dfac(2*k+1)/fac(il-k)/fac(k-j) + ENDDO + ENDDO + ENDDO + + END SUBROUTINE get_prefac_sabb ! ************************************************************************************************** !> \brief calculates the uncontracted, not normalized [s|1/r12|s] two-center coulomb integral @@ -685,28 +721,23 @@ CONTAINS !> \param la set of l quantum numbers for a !> \param npgfa number of primitive Gaussian on a !> \param nshella number of shells for a -!> \param zeta set of exponents on a !> \param scona_shg SHG contraction/normalization matrix for a !> \param lb set of l quantum numbers for b !> \param npgfb number of primitive Gaussian on b !> \param nshellb number of shells for b -!> \param zetb set of exponents on b !> \param sconb_shg SHG contraction/normalization matrix for b !> \param swork matrix storing the uncontracted and unnormalized s-integrals and derivatives !> \param swork_cont matrix storing the contracted and normalized s-integrals and derivatives !> \param calculate_forces ... ! ************************************************************************************************** - SUBROUTINE contract_sint_ab_clow(la, npgfa, nshella, zeta, scona_shg, & - lb, npgfb, nshellb, zetb, sconb_shg, & + SUBROUTINE contract_sint_ab_clow(la, npgfa, nshella, scona_shg, lb, npgfb, nshellb, sconb_shg, & swork, swork_cont, calculate_forces) INTEGER, DIMENSION(:), INTENT(IN) :: la INTEGER, INTENT(IN) :: npgfa, nshella - REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: scona_shg INTEGER, DIMENSION(:), INTENT(IN) :: lb INTEGER, INTENT(IN) :: npgfb, nshellb - REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sconb_shg REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: swork, swork_cont LOGICAL, INTENT(IN) :: calculate_forces @@ -716,7 +747,6 @@ CONTAINS INTEGER :: handle, ids, ids_start, ipgfa, ishella, & jpgfb, jshellb, lai, lbj, ndev, nds - REAL(KIND=dp) :: aif, bjf CALL timeset(routineN, handle) @@ -731,14 +761,12 @@ CONTAINS nds = lai+lbj+1 ids_start = nds-MIN(lai, lbj) DO ipgfa = 1, npgfa - aif = 1.0_dp/((2._dp*zeta(ipgfa))**lai) DO jpgfb = 1, npgfb - bjf = 1.0_dp/((2._dp*zetb(jpgfb))**lbj) DO ids = ids_start, nds+ndev swork_cont(ids, ishella, jshellb) = swork_cont(ids, ishella, jshellb) & +scona_shg(ipgfa, ishella) & *sconb_shg(jpgfb, jshellb) & - *aif*bjf*swork(ipgfa, jpgfb, ids) + *swork(ipgfa, jpgfb, ids) ENDDO ENDDO ENDDO @@ -795,154 +823,42 @@ CONTAINS END SUBROUTINE contract_sint_ab_chigh -! ************************************************************************************************** -!> \brief Contraction and normalization of the (aba) overlap -!> \param la set of l quantum numbers on a; orbital basis -!> \param npgfa number of primitive Gaussians functions on a; orbital basis -!> \param nshella number of shells for a; orbital basis -!> \param zeta set of exponents on a; orbital basis -!> \param scona_shg SHG contraction/normalization matrix for a; orbital basis -!> \param lb set of l quantum numbers on b; orbital basis -!> \param npgfb number of primitive Gaussians functions on b; orbital basis -!> \param nshellb number of shells for b; orbital basis -!> \param zetb set of exponents on b; orbital basis -!> \param sconb_shg SHG contraction/normalization matrix for b; orbital basis -!> \param lca set of l quantum numbers on a; ri basis -!> \param npgfca number of primitive Gaussians functions on a; ri basis -!> \param nshellca number of shells for a; ri basis -!> \param zetca set of exponents on a; ri basis -!> \param sconca_shg SHG contraction/normalization matrix for a; ri basis -!> \param swork matrix of storing the uncontracted and unnormalized s-integrals and derivatives -!> \param swork_cont matrix storing the contracted and normalized s-integrals and derivatives -!> \param calculate_forces ... -! ************************************************************************************************** - SUBROUTINE contract_s_overlap_aba(la, npgfa, nshella, zeta, scona_shg, & - lb, npgfb, nshellb, zetb, sconb_shg, & - lca, npgfca, nshellca, zetca, sconca_shg, & - swork, swork_cont, calculate_forces) - - INTEGER, DIMENSION(:), INTENT(IN) :: la - INTEGER, INTENT(IN) :: npgfa, nshella - REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta - REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: scona_shg - INTEGER, DIMENSION(:), INTENT(IN) :: lb - INTEGER, INTENT(IN) :: npgfb, nshellb - REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb - REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sconb_shg - INTEGER, DIMENSION(:), INTENT(IN) :: lca - INTEGER, INTENT(IN) :: npgfca, nshellca - REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetca - REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sconca_shg - REAL(KIND=dp), DIMENSION(:, :, :, 0:, :), & - INTENT(IN) :: swork - REAL(KIND=dp), DIMENSION(:, 0:, :, :, :), & - INTENT(INOUT) :: swork_cont - LOGICAL, INTENT(IN) :: calculate_forces - - CHARACTER(len=*), PARAMETER :: routineN = 'contract_s_overlap_aba', & - routineP = moduleN//':'//routineN - - INTEGER :: handle, ids, ids_start, iil, il, ipgfa, & - ishella, jpgfb, jshellb, kpgfa, & - kshella, laa, lai, lak, lbj, ndev, & - nds, nl - REAL(KIND=dp) :: a, aif, bjf, const, consta, constca, & - prefac, shg_fac(0:5) - - CALL timeset(routineN, handle) - - ndev = 0 - IF (calculate_forces) ndev = 1 - shg_fac(0) = 1.0_dp - - swork_cont = 0.0_dp - DO ishella = 1, nshella - lai = la(ishella) - consta = SQRT(1.0_dp/REAL(2*lai+1, dp)) - DO jshellb = 1, nshellb - lbj = lb(jshellb) - DO kshella = 1, nshellca - lak = lca(kshella) - constca = SQRT(1.0_dp/REAL(2*lak+1, dp)) - nl = INT((lai+lak)/2) - IF (lai == 0 .OR. lak == 0) nl = 0 - DO il = 0, nl - laa = lai+lak-2*il - nds = laa+lbj+1 - ids_start = nds-MIN(laa, lbj) - const = consta*constca*2.0_dp*SQRT(pi*REAL(2*laa+1, dp)) - DO iil = 1, il - shg_fac(iil) = fac(laa+iil-1)*ifac(laa)*REAL(laa, dp) & - *fac(il)/fac(il-iil)/fac(iil) - ENDDO - DO ipgfa = 1, npgfa - DO jpgfb = 1, npgfb - bjf = 1.0_dp/((2._dp*zetb(jpgfb))**lbj) - DO kpgfa = 1, npgfca - a = zeta(ipgfa)+zetca(kpgfa) - aif = 1.0_dp/((2._dp*a)**laa) - prefac = const*aif*bjf & - *scona_shg(ipgfa, ishella) & - *sconb_shg(jpgfb, jshellb) & - *sconca_shg(kpgfa, kshella) - DO ids = ids_start, nds+ndev - DO iil = 0, il - swork_cont(ids, il, ishella, jshellb, kshella) = & - swork_cont(ids, il, ishella, jshellb, kshella) & - +prefac*shg_fac(iil)*swork(ipgfa, jpgfb, kpgfa, il-iil, ids)/a**iil - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - ENDDO - - CALL timestop(handle) - - END SUBROUTINE contract_s_overlap_aba - ! ************************************************************************************************** !> \brief Contraction and normalization of the (abb) overlap !> \param la set of l quantum numbers on a !> \param npgfa number of primitive Gaussians functions on a; orbital basis !> \param nshella number of shells for a; orbital basis -!> \param zeta set of exponents on a; orbital basis !> \param scona_shg SHG contraction/normalization matrix for a; orbital basis !> \param lb set of l quantum numbers on b; orbital basis !> \param npgfb number of primitive Gaussians functions on b; orbital basis !> \param nshellb number of shells for b; orbital basis -!> \param zetb set of exponents on b; orbital basis -!> \param sconb_shg SHG contraction/normalization matrix for b; orbital basis !> \param lcb set of l quantum numbers on b; ri basis !> \param npgfcb number of primitive Gaussians functions on b; ri basis !> \param nshellcb number of shells for b; ri basis -!> \param zetcb set of exponents on b; ri basis -!> \param sconcb_shg SHG contraction/normalization matrix for b; ri basis +!> \param orbb_index index for orbital basis at B for sconb_mix +!> \param rib_index index for ri basis at B for sconb_mix +!> \param sconb_mix mixed contraction matrix for orbital and ri basis at B +!> \param nl_max related to the parameter m in (a|rb^(2m)|b) +!> \param nds_max derivative with respect to rab**2 !> \param swork matrix of storing the uncontracted and unnormalized s-integrals and derivatives !> \param swork_cont matrix storing the contracted and normalized s-integrals and derivatives !> \param calculate_forces ... ! ************************************************************************************************** - SUBROUTINE contract_s_overlap_abb(la, npgfa, nshella, zeta, scona_shg, & - lb, npgfb, nshellb, zetb, sconb_shg, & - lcb, npgfcb, nshellcb, zetcb, sconcb_shg, & - swork, swork_cont, calculate_forces) + SUBROUTINE contract_s_overlap_abb(la, npgfa, nshella, scona_shg, lb, npgfb, nshellb, & + lcb, npgfcb, nshellcb, orbb_index, rib_index, sconb_mix, & + nl_max, nds_max, swork, swork_cont, calculate_forces) INTEGER, DIMENSION(:), INTENT(IN) :: la INTEGER, INTENT(IN) :: npgfa, nshella - REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: scona_shg INTEGER, DIMENSION(:), INTENT(IN) :: lb INTEGER, INTENT(IN) :: npgfb, nshellb - REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb - REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sconb_shg INTEGER, DIMENSION(:), INTENT(IN) :: lcb INTEGER, INTENT(IN) :: npgfcb, nshellcb - REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetcb - REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sconcb_shg - REAL(KIND=dp), DIMENSION(:, :, :, 0:, :), & + INTEGER, DIMENSION(:, :), INTENT(IN) :: orbb_index, rib_index + REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: sconb_mix + INTEGER, INTENT(IN) :: nl_max, nds_max + REAL(KIND=dp), DIMENSION(:, :, :, :, :), & INTENT(IN) :: swork REAL(KIND=dp), DIMENSION(:, 0:, :, :, :), & INTENT(INOUT) :: swork_cont @@ -951,26 +867,27 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'contract_s_overlap_abb', & routineP = moduleN//':'//routineN - INTEGER :: handle, ids, ids_start, iil, il, ipgfa, & - ishella, jpgfb, jshellb, kpgfb, & + INTEGER :: forb, fri, handle, ids, ids_start, iil, & + il, ishella, jpgfb, jshellb, kpgfb, & kshellb, lai, lbb, lbj, lbk, ndev, & nds, nl - REAL(KIND=dp) :: aif, b, bjf, const, constb, constcb, & - prefac, shg_fac(0:5) + REAL(KIND=dp), ALLOCATABLE, & + DIMENSION(:, :, :, :, :) :: work_ppc CALL timeset(routineN, handle) ndev = 0 IF (calculate_forces) ndev = 1 - shg_fac(0) = 1.0_dp + ALLOCATE (work_ppc(npgfb, npgfcb, nl_max, nds_max, nshella)) + work_ppc = 0.0_dp + CALL dgemm("T", "N", npgfb*npgfcb*nl_max*nds_max, nshella, npgfa, 1._dp, swork, npgfa, & + scona_shg, npgfa, 0.0_dp, work_ppc, npgfb*npgfcb*nl_max*nds_max) swork_cont = 0.0_dp DO kshellb = 1, nshellcb lbk = lcb(kshellb) - constcb = SQRT(1.0_dp/REAL(2*lbk+1, dp)) DO jshellb = 1, nshellb lbj = lb(jshellb) - constb = SQRT(1.0_dp/REAL(2*lbj+1, dp)) nl = INT((lbj+lbk)/2) IF (lbj == 0 .OR. lbk == 0) nl = 0 DO ishella = 1, nshella @@ -979,27 +896,15 @@ CONTAINS lbb = lbj+lbk-2*il nds = lai+lbb+1 ids_start = nds-MIN(lai, lbb) - const = constb*constcb*2.0_dp*SQRT(pi*REAL(2*lbb+1, dp)) - DO iil = 1, il - shg_fac(iil) = fac(lbb+iil-1)*ifac(lbb)*REAL(lbb, dp) & - *fac(il)/fac(il-iil)/fac(iil) - ENDDO - DO ipgfa = 1, npgfa - aif = 1.0_dp/((2._dp*zeta(ipgfa))**lai) - DO jpgfb = 1, npgfb - DO kpgfb = 1, npgfcb - b = zetb(jpgfb)+zetcb(kpgfb) - bjf = 1.0_dp/((2._dp*b)**lbb) - prefac = const*aif*bjf & - *scona_shg(ipgfa, ishella) & - *sconb_shg(jpgfb, jshellb) & - *sconcb_shg(kpgfb, kshellb) - DO ids = ids_start, nds+ndev - DO iil = 0, il - swork_cont(ids, il, ishella, jshellb, kshellb) = & - swork_cont(ids, il, ishella, jshellb, kshellb) & - +prefac*shg_fac(iil)*swork(ipgfa, jpgfb, kpgfb, il-iil, ids)/b**iil - ENDDO + DO jpgfb = 1, npgfb + forb = orbb_index(jpgfb, jshellb) + DO kpgfb = 1, npgfcb + fri = rib_index(kpgfb, kshellb) + DO ids = ids_start, nds+ndev + DO iil = 0, il + swork_cont(ids, il, ishella, jshellb, kshellb) = & + swork_cont(ids, il, ishella, jshellb, kshellb) & + +sconb_mix(iil+1, fri, forb, il+1)*work_ppc(jpgfb, kpgfb, il-iil+1, ids, ishella) ENDDO ENDDO ENDDO @@ -1008,8 +913,105 @@ CONTAINS ENDDO ENDDO ENDDO + DEALLOCATE (work_ppc) CALL timestop(handle) END SUBROUTINE contract_s_overlap_abb +! ************************************************************************************************** +!> \brief Contraction and normalization of the (aba) overlap +!> \param la set of l quantum numbers on a; orbital basis +!> \param npgfa number of primitive Gaussians functions on a; orbital basis +!> \param nshella number of shells for a; orbital basis +!> \param lb set of l quantum numbers on b; orbital basis +!> \param npgfb number of primitive Gaussians functions on b; orbital basis +!> \param nshellb number of shells for b; orbital basis +!> \param sconb_shg SHG contraction/normalization matrix for b; orbital basis +!> \param lca set of l quantum numbers on a; ri basis +!> \param npgfca number of primitive Gaussians functions on a; ri basis +!> \param nshellca number of shells for a; ri basis +!> \param orba_index index for orbital basis at A for scona_mix +!> \param ria_index index for ri basis at A for scona_mix +!> \param scona_mix mixed contraction matrix for orbital and ri basis at A +!> \param nl_max related to the parameter m in (a|ra^(2m)|b) +!> \param nds_max maximal derivative with respect to rab**2 +!> \param swork matrix of storing the uncontracted and unnormalized s-integrals and derivatives +!> \param swork_cont matrix storing the contracted and normalized s-integrals and derivatives +!> \param calculate_forces ... +! ************************************************************************************************** + SUBROUTINE contract_s_overlap_aba(la, npgfa, nshella, lb, npgfb, nshellb, sconb_shg, & + lca, npgfca, nshellca, orba_index, ria_index, scona_mix, & + nl_max, nds_max, swork, swork_cont, calculate_forces) + + INTEGER, DIMENSION(:), INTENT(IN) :: la + INTEGER, INTENT(IN) :: npgfa, nshella + INTEGER, DIMENSION(:), INTENT(IN) :: lb + INTEGER, INTENT(IN) :: npgfb, nshellb + REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sconb_shg + INTEGER, DIMENSION(:), INTENT(IN) :: lca + INTEGER, INTENT(IN) :: npgfca, nshellca + INTEGER, DIMENSION(:, :), INTENT(IN) :: orba_index, ria_index + REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: scona_mix + INTEGER, INTENT(IN) :: nl_max, nds_max + REAL(KIND=dp), DIMENSION(:, :, :, :, :), & + INTENT(IN) :: swork + REAL(KIND=dp), DIMENSION(:, 0:, :, :, :), & + INTENT(INOUT) :: swork_cont + LOGICAL, INTENT(IN) :: calculate_forces + + CHARACTER(len=*), PARAMETER :: routineN = 'contract_s_overlap_aba', & + routineP = moduleN//':'//routineN + + INTEGER :: forb, fri, handle, ids, ids_start, iil, & + il, ipgfa, ishella, jshellb, kpgfa, & + kshella, laa, lai, lak, lbj, ndev, & + nds, nl + REAL(KIND=dp), ALLOCATABLE, & + DIMENSION(:, :, :, :, :) :: work_ppc + + CALL timeset(routineN, handle) + + ndev = 0 + IF (calculate_forces) ndev = 1 + + ALLOCATE (work_ppc(npgfa, npgfca, nl_max, nds_max, nshellb)) + work_ppc = 0.0_dp + CALL dgemm("T", "N", npgfa*npgfca*nl_max*nds_max, nshellb, npgfb, 1._dp, swork, npgfb, & + sconb_shg, npgfb, 0.0_dp, work_ppc, npgfa*npgfca*nl_max*nds_max) + swork_cont = 0.0_dp + DO kshella = 1, nshellca + lak = lca(kshella) + DO jshellb = 1, nshellb + lbj = lb(jshellb) + DO ishella = 1, nshella + lai = la(ishella) + nl = INT((lai+lak)/2) + IF (lai == 0 .OR. lak == 0) nl = 0 + DO il = 0, nl + laa = lai+lak-2*il + nds = laa+lbj+1 + ids_start = nds-MIN(laa, lbj) + DO kpgfa = 1, npgfca + fri = ria_index(kpgfa, kshella) + DO ipgfa = 1, npgfa + forb = orba_index(ipgfa, ishella) + DO ids = ids_start, nds+ndev + DO iil = 0, il + swork_cont(ids, il, ishella, jshellb, kshella) = & + swork_cont(ids, il, ishella, jshellb, kshella) & + +scona_mix(iil+1, fri, forb, il+1)*work_ppc(ipgfa, kpgfa, il-iil+1, ids, jshellb) + ENDDO + ENDDO + ENDDO + ENDDO + ENDDO + ENDDO + ENDDO + ENDDO + + DEALLOCATE (work_ppc) + CALL timestop(handle) + + END SUBROUTINE contract_s_overlap_aba + END MODULE ai_s_contract_shg diff --git a/src/generic_shg_integrals.F b/src/generic_shg_integrals.F index 0670834c02..0c932bf3cf 100644 --- a/src/generic_shg_integrals.F +++ b/src/generic_shg_integrals.F @@ -23,7 +23,7 @@ MODULE generic_shg_integrals get_dW_matrix, get_real_scaled_solid_harmonic USE ai_s_contract_shg, ONLY: & contract_s_overlap_aba, contract_s_overlap_abb, contract_sint_ab_chigh, & - contract_sint_ab_clow, s_coulomb_ab, s_gauss_ab, s_overlap_ab, s_overlap_abx, s_verf_ab, & + contract_sint_ab_clow, s_coulomb_ab, s_gauss_ab, s_overlap_ab, s_overlap_abb, s_verf_ab, & s_verfc_ab, s_vgauss_ab USE basis_set_types, ONLY: gto_basis_set_type USE input_constants, ONLY: operator_coulomb,& @@ -42,7 +42,8 @@ MODULE generic_shg_integrals CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'generic_shg_integrals' PUBLIC :: two_center_shg_integrals, lri_precalc_angular_shg_part, & - int_overlap_ab_shg, int_overlap_aba_shg, int_overlap_abb_shg + int_overlap_ab_shg, int_overlap_aba_shg, int_overlap_abb_shg, & + precalc_angular_shg_part ABSTRACT INTERFACE ! ************************************************************************************************** @@ -105,7 +106,7 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'two_center_shg_integrals', & routineP = moduleN//':'//routineN - INTEGER :: handle + INTEGER :: handle, la_max, lb_max REAL(KIND=dp) :: my_omega REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Waux_mat REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dWaux_mat @@ -115,7 +116,10 @@ CONTAINS NULLIFY (s_operator_ab) CALL timeset(routineN, handle) - CALL precalc_angular_shg_part(fba, fbb, rab, Waux_mat, dWaux_mat, calculate_forces) + la_max = MAXVAL(fba%lmax) + lb_max = MAXVAL(fbb%lmax) + + CALL precalc_angular_shg_part(la_max, lb_max, rab, Waux_mat, dWaux_mat, calculate_forces) my_omega = 1.0_dp SELECT CASE (r12_operator) @@ -148,16 +152,16 @@ CONTAINS ! ************************************************************************************************** !> \brief precalculates the angular part of the SHG integrals -!> \param fba basis at center A -!> \param fbb basis at center B +!> \param la_max ... +!> \param lb_max ... !> \param rab distance vector between a and b !> \param Waux_mat W matrix that contains the angular-dependent part !> \param dWaux_mat derivative of the W matrix !> \param calculate_forces ... ! ************************************************************************************************** - SUBROUTINE precalc_angular_shg_part(fba, fbb, rab, Waux_mat, dWaux_mat, calculate_forces) + SUBROUTINE precalc_angular_shg_part(la_max, lb_max, rab, Waux_mat, dWaux_mat, calculate_forces) - TYPE(gto_basis_set_type), POINTER :: fba, fbb + INTEGER, INTENT(IN) :: la_max, lb_max REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), & INTENT(OUT) :: Waux_mat @@ -168,7 +172,7 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'precalc_angular_shg_part', & routineP = moduleN//':'//routineN - INTEGER :: handle, la_max, lb_max, lmax, mdim(3) + INTEGER :: handle, lmax, mdim(3) INTEGER, DIMENSION(:), POINTER :: la_max_all REAL(KIND=dp) :: rab2 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Rc, Rs @@ -177,9 +181,6 @@ CONTAINS NULLIFY (la_max_all) rab2 = rab(1)*rab(1)+rab(2)*rab(2)+rab(3)*rab(3) - - la_max = MAXVAL(fba%lmax) - lb_max = MAXVAL(fbb%lmax) lmax = MAX(la_max, lb_max) ALLOCATE (la_max_all(0:lb_max)) @@ -239,7 +240,7 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'int_operator_ab_shg', routineP = moduleN//':'//routineN INTEGER :: handle, iset, jset, la_max_set, lb_max_set, ndev, nds, nds_max, npgfa_set, & - npgfb_set, nseta, nsetb, nsgfa_set, nsgfb_set, nshella_set, nshellb_set, sgfa, sgfb + npgfb_set, nseta, nsetb, nsgfa_set, nsgfb_set, nshella_set, nshellb_set INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, npgfa, npgfb, nsgfa, & nsgfb, nshella, nshellb INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, la, lb @@ -295,17 +296,17 @@ CONTAINS dvab = 0.0_dp DO iset = 1, nseta - sgfa = first_sgfa(1, iset) DO jset = 1, nsetb - sgfb = first_sgfb(1, jset) CALL s_operator_ab(la_max(iset), npgfa(iset), zeta(:, iset), & lb_max(jset), npgfb(jset), zetb(:, jset), & omega, rab, swork, calculate_forces) nds = la_max(iset)+lb_max(jset)+ndev+1 - CALL contract_sint_ab_chigh(npgfa(iset), nshella(iset), scona_shg(1:npgfa(iset), 1:nshella(iset), iset), & - npgfb(jset), nshellb(jset), sconb_shg(1:npgfb(jset), 1:nshellb(jset), jset), & + CALL contract_sint_ab_chigh(npgfa(iset), nshella(iset), & + scona_shg(1:npgfa(iset), 1:nshella(iset), iset), & + npgfb(jset), nshellb(jset), & + sconb_shg(1:npgfb(jset), 1:nshellb(jset), jset), & nds, swork(1:npgfa(iset), 1:npgfb(jset), 1:nds), & swork_cont(1:nds, 1:nshella(iset), 1:nshellb(jset))) CALL construct_int_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), & @@ -337,9 +338,10 @@ CONTAINS !> \param Waux_mat W matrix that contains the angular-dependent part !> \param dWaux_mat derivative of the W matrix !> \param calculate_forces ... +!> \param contraction_high ... ! ************************************************************************************************** SUBROUTINE int_overlap_ab_shg(sab, dsab, rab, fba, fbb, scona_shg, sconb_shg, Waux_mat, dWaux_mat, & - calculate_forces) + calculate_forces, contraction_high) REAL(KIND=dp), DIMENSION(:, :), POINTER :: sab REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dsab @@ -349,17 +351,17 @@ CONTAINS REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Waux_mat REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dWaux_mat LOGICAL, INTENT(IN) :: calculate_forces + LOGICAL, INTENT(IN), OPTIONAL :: contraction_high CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_ab_shg', & routineP = moduleN//':'//routineN - INTEGER :: handle, iset, jset, la_max_set, & - lb_max_set, ndev, nds_max, npgfa_set, & - npgfb_set, nseta, nsetb, nshella_set, & - nshellb_set + INTEGER :: handle, iset, jset, la_max_set, lb_max_set, ndev, nds, nds_max, npgfa_set, & + npgfb_set, nseta, nsetb, nshella_set, nshellb_set INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, npgfa, npgfb, nshella, & nshellb INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, la, lb + LOGICAL :: my_contraction_high REAL(KIND=dp) :: dab REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: swork, swork_cont REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b @@ -404,6 +406,12 @@ CONTAINS ALLOCATE (swork(npgfa_set, npgfb_set, nds_max)) ALLOCATE (swork_cont(nds_max, nshella_set, nshellb_set)) + sab = 0.0_dp + dsab = 0.0_dp + + my_contraction_high = .TRUE. + IF (PRESENT(contraction_high)) my_contraction_high = contraction_high + DO iset = 1, nseta DO jset = 1, nsetb @@ -413,11 +421,21 @@ CONTAINS CALL s_overlap_ab(la_max(iset), npgfa(iset), zeta(:, iset), & lb_max(jset), npgfb(jset), zetb(:, jset), & rab, swork, calculate_forces) - CALL contract_sint_ab_clow(la(:, iset), npgfa(iset), nshella(iset), & - zeta(:, iset), scona_shg(:, :, iset), & - lb(:, jset), npgfb(jset), nshellb(jset), & - zetb(:, jset), sconb_shg(:, :, jset), & - swork, swork_cont, calculate_forces) + IF (my_contraction_high) THEN + nds = la_max(iset)+lb_max(jset)+ndev+1 + CALL contract_sint_ab_chigh(npgfa(iset), nshella(iset), & + scona_shg(1:npgfa(iset), 1:nshella(iset), iset), & + npgfb(jset), nshellb(jset), & + sconb_shg(1:npgfb(jset), 1:nshellb(jset), jset), & + nds, swork(1:npgfa(iset), 1:npgfb(jset), 1:nds), & + swork_cont(1:nds, 1:nshella(iset), 1:nshellb(jset))) + ELSE + CALL contract_sint_ab_clow(la(:, iset), npgfa(iset), nshella(iset), & + scona_shg(:, :, iset), & + lb(:, jset), npgfb(jset), nshellb(jset), & + sconb_shg(:, :, jset), & + swork, swork_cont, calculate_forces) + ENDIF CALL construct_int_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), & lb(:, jset), first_sgfb(:, jset), nshellb(jset), & swork_cont, Waux_mat, sab) @@ -436,159 +454,6 @@ CONTAINS END SUBROUTINE int_overlap_ab_shg -! ************************************************************************************************** -!> \brief calculate integrals (a,b,fa) -!> \param abaint integral (a,b,fa) -!> \param dabdaint ... -!> \param rab distance vector between A and B -!> \param oba orbital basis at center A -!> \param obb orbital basis at center B -!> \param fba auxiliary basis set at center A -!> \param scona_shg contraction matrix for orb bas on A -!> \param sconb_shg contraction matrix for orb bas on B -!> \param sconca_shg contraction matrix for aux bas on A -!> \param cg_coeff Clebsch-Gordon coefficients -!> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients -!> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients -!> \param Waux_mat W matrix that contains the angular-dependent part -!> \param dWaux_mat derivative of the W matrix -!> \param calculate_forces ... -! ************************************************************************************************** - SUBROUTINE int_overlap_aba_shg(abaint, dabdaint, rab, oba, obb, fba, scona_shg, & - sconb_shg, sconca_shg, cg_coeff, cg_none0_list, ncg_none0, & - Waux_mat, dWaux_mat, calculate_forces) - - REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: abaint - REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: dabdaint - REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab - TYPE(gto_basis_set_type), POINTER :: oba, obb, fba - REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: scona_shg, sconb_shg, sconca_shg, & - cg_coeff - INTEGER, DIMENSION(:, :, :), INTENT(IN) :: cg_none0_list - INTEGER, DIMENSION(:, :), INTENT(IN) :: ncg_none0 - REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Waux_mat - REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dWaux_mat - LOGICAL, INTENT(IN) :: calculate_forces - - CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_aba_shg', & - routineP = moduleN//':'//routineN - - INTEGER :: handle, iset, jset, kset, la_max_set, laa_max_set, lb_max_set, lca_max_set, ndev, & - nds_max, nl_set, npgfa_set, npgfb_set, npgfca_set, nseta, nsetb, nsetca, nshella_set, & - nshellb_set, nshellca_set - INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, lca_max, npgfa, npgfb, & - npgfca, nshella, nshellb, nshellca - INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, first_sgfca, la, & - lb, lca - REAL(KIND=dp) :: dab, rab2 - REAL(KIND=dp), ALLOCATABLE, & - DIMENSION(:, :, :, :, :) :: swork, swork_cont - REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b, set_radius_ca - REAL(KIND=dp), DIMENSION(:, :), POINTER :: zeta, zetb, zetca - - CALL timeset(routineN, handle) - NULLIFY (la_max, lb_max, lca_max, npgfa, npgfb, npgfca) - NULLIFY (first_sgfa, first_sgfb, first_sgfca, set_radius_a, set_radius_b, & - set_radius_ca, zeta, zetb, zetca) - - ! basis ikind - first_sgfa => oba%first_sgf - la_max => oba%lmax - la => oba%l - npgfa => oba%npgf - nshella => oba%nshell - nseta = oba%nset - set_radius_a => oba%set_radius - zeta => oba%zet - ! basis jkind - first_sgfb => obb%first_sgf - lb_max => obb%lmax - lb => obb%l - npgfb => obb%npgf - nshellb => obb%nshell - nsetb = obb%nset - set_radius_b => obb%set_radius - zetb => obb%zet - - ! basis RI A - first_sgfca => fba%first_sgf - lca_max => fba%lmax - lca => fba%l - npgfca => fba%npgf - nshellca => fba%nshell - nsetca = fba%nset - set_radius_ca => fba%set_radius - zetca => fba%zet - - dab = SQRT(SUM(rab**2)) - rab2 = rab(1)*rab(1)+rab(2)*rab(2)+rab(3)*rab(3) - - la_max_set = MAXVAL(la_max) - lb_max_set = MAXVAL(lb_max) - lca_max_set = MAXVAL(lca_max) - npgfa_set = MAXVAL(npgfa) - npgfb_set = MAXVAL(npgfb) - npgfca_set = MAXVAL(npgfca) - nshella_set = MAXVAL(nshella) - nshellb_set = MAXVAL(nshellb) - nshellca_set = MAXVAL(nshellca) - !*** for forces: derivative+1 in auxiliary vector required - ndev = 0 - IF (calculate_forces) ndev = 1 - - laa_max_set = la_max_set+lca_max_set - - ! allocate some work storage.... - nds_max = laa_max_set+lb_max_set+ndev+1 - nl_set = INT((laa_max_set)/2) - ALLOCATE (swork(npgfa_set, npgfb_set, npgfca_set, 0:nl_set, nds_max)) - ALLOCATE (swork_cont(nds_max, 0:nl_set, nshella_set, nshellb_set, nshellca_set)) - - DO iset = 1, nseta - - DO jset = 1, nsetb - - IF (set_radius_a(iset)+set_radius_b(jset) < dab) CYCLE - - DO kset = 1, nsetca - - IF (set_radius_b(jset)+set_radius_ca(kset) < dab) CYCLE - - CALL s_overlap_abx(la_max(iset), npgfa(iset), zeta(:, iset), & - lb_max(jset), npgfb(jset), zetb(:, jset), & - lca_max(kset), npgfca(kset), zetca(:, kset), & - rab, swork, calculate_forces, calc_aba=.TRUE.) - - CALL contract_s_overlap_aba(la(:, iset), npgfa(iset), nshella(iset), zeta(:, iset), scona_shg(:, :, iset), & - lb(:, jset), npgfb(jset), nshellb(jset), zetb(:, jset), sconb_shg(:, :, jset), & - lca(:, kset), npgfca(kset), nshellca(kset), zetca(:, kset), sconca_shg(:, :, kset), & - swork, swork_cont, calculate_forces) - - CALL construct_overlap_shg_aba(la(:, iset), first_sgfa(:, iset), nshella(iset), & - lb(:, jset), first_sgfb(:, jset), nshellb(jset), & - lca(:, kset), first_sgfca(:, kset), nshellca(kset), & - cg_coeff, cg_none0_list, ncg_none0, & - swork_cont, Waux_mat, abaint) - IF (calculate_forces) THEN - !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra - CALL dev_overlap_shg_aba(la(:, iset), first_sgfa(:, iset), nshella(iset), & - lb(:, jset), first_sgfb(:, jset), nshellb(jset), & - lca(:, kset), first_sgfca(:, kset), nshellca(kset), & - cg_coeff, cg_none0_list, ncg_none0, & - -rab, swork_cont, Waux_mat, dWaux_mat, dabdaint) - ENDIF - - END DO - END DO - END DO - - DEALLOCATE (swork_cont) - DEALLOCATE (swork) - - CALL timestop(handle) - - END SUBROUTINE int_overlap_aba_shg - ! ************************************************************************************************** !> \brief calculate integrals (a,b,fb) !> \param abbint integral (a,b,fb) @@ -597,9 +462,10 @@ CONTAINS !> \param oba orbital basis at center A !> \param obb orbital basis at center B !> \param fbb auxiliary basis set at center B -!> \param scona_shg contraction matrix for orb bas on A -!> \param sconb_shg contraction matrix for orb bas on B -!> \param sconcb_shg contraction matrix for aux bas on B +!> \param scon_oba contraction matrix for orb bas on A +!> \param sconb_mix mixed contraction matrix orb + ri basis on B +!> \param obb_index orbital basis index for sconb_mix +!> \param fbb_index ri basis index for sconb_mix !> \param cg_coeff Clebsch-Gordon coefficients !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients @@ -607,16 +473,18 @@ CONTAINS !> \param dWaux_mat derivative of the W matrix !> \param calculate_forces ... ! ************************************************************************************************** - SUBROUTINE int_overlap_abb_shg(abbint, dabbint, rab, oba, obb, fbb, scona_shg, & - sconb_shg, sconcb_shg, cg_coeff, cg_none0_list, ncg_none0, & - Waux_mat, dWaux_mat, calculate_forces) + SUBROUTINE int_overlap_abb_shg(abbint, dabbint, rab, oba, obb, fbb, scon_oba, sconb_mix, & + obb_index, fbb_index, cg_coeff, cg_none0_list, & + ncg_none0, Waux_mat, dWaux_mat, calculate_forces) REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: abbint REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: dabbint REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab TYPE(gto_basis_set_type), POINTER :: oba, obb, fbb - REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: scona_shg, sconb_shg, sconcb_shg, & - cg_coeff + REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: scon_oba + REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: sconb_mix + INTEGER, DIMENSION(:, :, :), INTENT(IN) :: obb_index, fbb_index + REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: cg_coeff INTEGER, DIMENSION(:, :, :), INTENT(IN) :: cg_none0_list INTEGER, DIMENSION(:, :), INTENT(IN) :: ncg_none0 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Waux_mat @@ -626,9 +494,9 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_abb_shg', & routineP = moduleN//':'//routineN - INTEGER :: handle, iset, jset, kset, la_max_set, lb_max_set, lbb_max_set, lcb_max_set, ndev, & - nds_max, nl_set, npgfa_set, npgfb_set, npgfcb_set, nseta, nsetb, nsetcb, nshella_set, & - nshellb_set, nshellcb_set + INTEGER :: handle, iset, jset, kset, la_max_set, lb_max_set, lbb_max, lbb_max_set, & + lcb_max_set, ndev, nds, nds_max, nl, nl_set, npgfa_set, npgfb_set, npgfcb_set, nseta, & + nsetb, nsetcb, nshella_set, nshellb_set, nshellcb_set INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, lcb_max, npgfa, npgfb, & npgfcb, nshella, nshellb, nshellcb INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, first_sgfcb, la, & @@ -694,7 +562,7 @@ CONTAINS ! allocate some work storage.... nds_max = la_max_set+lbb_max_set+ndev+1 nl_set = INT((lbb_max_set)/2) - ALLOCATE (swork(1:npgfa_set, 1:npgfb_set, 1:npgfcb_set, 0:nl_set, 1:nds_max)) + ALLOCATE (swork(npgfa_set, npgfb_set, npgfcb_set, nl_set+1, nds_max)) ALLOCATE (swork_cont(nds_max, 0:nl_set, nshella_set, nshellb_set, nshellcb_set)) DO iset = 1, nseta @@ -707,15 +575,22 @@ CONTAINS IF (set_radius_a(iset)+set_radius_cb(kset) < dab) CYCLE - CALL s_overlap_abx(la_max(iset), npgfa(iset), zeta(:, iset), & + CALL s_overlap_abb(la_max(iset), npgfa(iset), zeta(:, iset), & lb_max(jset), npgfb(jset), zetb(:, jset), & lcb_max(kset), npgfcb(kset), zetcb(:, kset), & - rab, swork, calculate_forces, calc_aba=.FALSE.) + rab, swork, calculate_forces) - CALL contract_s_overlap_abb(la(:, iset), npgfa(iset), nshella(iset), zeta(:, iset), scona_shg(:, :, iset), & - lb(:, jset), npgfb(jset), nshellb(jset), zetb(:, jset), sconb_shg(:, :, jset), & - lcb(:, kset), npgfcb(kset), nshellcb(kset), zetcb(:, kset), sconcb_shg(:, :, kset), & - swork, swork_cont, calculate_forces) + lbb_max = lb_max(jset)+lcb_max(kset) + nds = la_max(iset)+lbb_max+ndev+1 + nl = INT((lbb_max)/2)+1 + + CALL contract_s_overlap_abb(la(:, iset), npgfa(iset), nshella(iset), & + scon_oba(1:npgfa(iset), 1:nshella(iset), iset), & + lb(:, jset), npgfb(jset), nshellb(jset), & + lcb(:, kset), npgfcb(kset), nshellcb(kset), & + obb_index(:, :, jset), fbb_index(:, :, kset), sconb_mix, nl, nds, & + swork(1:npgfa(iset), 1:npgfb(jset), 1:npgfcb(kset), 1:nl, 1:nds), & + swork_cont, calculate_forces) CALL construct_overlap_shg_abb(la(:, iset), first_sgfa(:, iset), nshella(iset), & lb(:, jset), first_sgfb(:, jset), nshellb(jset), & @@ -742,6 +617,169 @@ CONTAINS END SUBROUTINE int_overlap_abb_shg +! ************************************************************************************************** +!> \brief calculate integrals (a,b,fa) +!> \param abaint integral (a,b,fa) +!> \param dabdaint ... +!> \param rab distance vector between A and B +!> \param oba orbital basis at center A +!> \param obb orbital basis at center B +!> \param fba auxiliary basis set at center A +!> \param scon_obb contraction matrix for orb bas on B +!> \param scona_mix mixed contraction matrix orb + ri basis on A +!> \param oba_index orbital basis index for scona_mix +!> \param fba_index ri basis index for scona_mix +!> \param cg_coeff Clebsch-Gordon coefficients +!> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients +!> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients +!> \param Waux_mat W matrix that contains the angular-dependent part +!> \param dWaux_mat derivative of the W matrix +!> \param calculate_forces ... +! ************************************************************************************************** + SUBROUTINE int_overlap_aba_shg(abaint, dabdaint, rab, oba, obb, fba, scon_obb, scona_mix, & + oba_index, fba_index, cg_coeff, cg_none0_list, & + ncg_none0, Waux_mat, dWaux_mat, calculate_forces) + + REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: abaint + REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: dabdaint + REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab + TYPE(gto_basis_set_type), POINTER :: oba, obb, fba + REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: scon_obb + REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: scona_mix + INTEGER, DIMENSION(:, :, :), INTENT(IN) :: oba_index, fba_index + REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: cg_coeff + INTEGER, DIMENSION(:, :, :), INTENT(IN) :: cg_none0_list + INTEGER, DIMENSION(:, :), INTENT(IN) :: ncg_none0 + REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Waux_mat + REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dWaux_mat + LOGICAL, INTENT(IN) :: calculate_forces + + CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_aba_shg', & + routineP = moduleN//':'//routineN + + INTEGER :: handle, iset, jset, kset, la_max_set, laa_max, laa_max_set, lb_max_set, & + lca_max_set, ndev, nds, nds_max, nl, nl_set, npgfa_set, npgfb_set, npgfca_set, nseta, & + nsetb, nsetca, nshella_set, nshellb_set, nshellca_set + INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, lca_max, npgfa, npgfb, & + npgfca, nshella, nshellb, nshellca + INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, first_sgfca, la, & + lb, lca + REAL(KIND=dp) :: dab, rab2 + REAL(KIND=dp), ALLOCATABLE, & + DIMENSION(:, :, :, :, :) :: swork, swork_cont + REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b, set_radius_ca + REAL(KIND=dp), DIMENSION(:, :), POINTER :: zeta, zetb, zetca + + CALL timeset(routineN, handle) + NULLIFY (la_max, lb_max, lca_max, npgfa, npgfb, npgfca) + NULLIFY (first_sgfa, first_sgfb, first_sgfca, set_radius_a, set_radius_b, & + set_radius_ca, zeta, zetb, zetca) + + ! basis ikind + first_sgfa => oba%first_sgf + la_max => oba%lmax + la => oba%l + npgfa => oba%npgf + nshella => oba%nshell + nseta = oba%nset + set_radius_a => oba%set_radius + zeta => oba%zet + ! basis jkind + first_sgfb => obb%first_sgf + lb_max => obb%lmax + lb => obb%l + npgfb => obb%npgf + nshellb => obb%nshell + nsetb = obb%nset + set_radius_b => obb%set_radius + zetb => obb%zet + + ! basis RI A + first_sgfca => fba%first_sgf + lca_max => fba%lmax + lca => fba%l + npgfca => fba%npgf + nshellca => fba%nshell + nsetca = fba%nset + set_radius_ca => fba%set_radius + zetca => fba%zet + + dab = SQRT(SUM(rab**2)) + rab2 = rab(1)*rab(1)+rab(2)*rab(2)+rab(3)*rab(3) + + la_max_set = MAXVAL(la_max) + lb_max_set = MAXVAL(lb_max) + lca_max_set = MAXVAL(lca_max) + npgfa_set = MAXVAL(npgfa) + npgfb_set = MAXVAL(npgfb) + npgfca_set = MAXVAL(npgfca) + nshella_set = MAXVAL(nshella) + nshellb_set = MAXVAL(nshellb) + nshellca_set = MAXVAL(nshellca) + !*** for forces: derivative+1 in auxiliary vector required + ndev = 0 + IF (calculate_forces) ndev = 1 + + laa_max_set = la_max_set+lca_max_set + + ! allocate some work storage.... + nds_max = laa_max_set+lb_max_set+ndev+1 + nl_set = INT((laa_max_set)/2) + ALLOCATE (swork(npgfb_set, npgfa_set, npgfca_set, nl_set+1, nds_max)) + ALLOCATE (swork_cont(nds_max, 0:nl_set, nshella_set, nshellb_set, nshellca_set)) + + DO iset = 1, nseta + + DO jset = 1, nsetb + + IF (set_radius_a(iset)+set_radius_b(jset) < dab) CYCLE + + DO kset = 1, nsetca + + IF (set_radius_b(jset)+set_radius_ca(kset) < dab) CYCLE + + !*** calculate s_baa here + CALL s_overlap_abb(lb_max(jset), npgfb(jset), zetb(:, jset), & + la_max(iset), npgfa(iset), zeta(:, iset), & + lca_max(kset), npgfca(kset), zetca(:, kset), & + rab, swork, calculate_forces) + + laa_max = la_max(iset)+lca_max(kset) + nds = laa_max+lb_max(jset)+ndev+1 + nl = INT(laa_max/2)+1 + + CALL contract_s_overlap_aba(la(:, iset), npgfa(iset), nshella(iset), & + lb(:, jset), npgfb(jset), nshellb(jset), & + scon_obb(1:npgfb(jset), 1:nshellb(jset), jset), & + lca(:, kset), npgfca(kset), nshellca(kset), & + oba_index(:, :, iset), fba_index(:, :, kset), scona_mix, nl, nds, & + swork(1:npgfb(jset), 1:npgfa(iset), 1:npgfca(kset), 1:nl, 1:nds), & + swork_cont, calculate_forces) + CALL construct_overlap_shg_aba(la(:, iset), first_sgfa(:, iset), nshella(iset), & + lb(:, jset), first_sgfb(:, jset), nshellb(jset), & + lca(:, kset), first_sgfca(:, kset), nshellca(kset), & + cg_coeff, cg_none0_list, ncg_none0, & + swork_cont, Waux_mat, abaint) + IF (calculate_forces) THEN + !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra + CALL dev_overlap_shg_aba(la(:, iset), first_sgfa(:, iset), nshella(iset), & + lb(:, jset), first_sgfb(:, jset), nshellb(jset), & + lca(:, kset), first_sgfca(:, kset), nshellca(kset), & + cg_coeff, cg_none0_list, ncg_none0, & + -rab, swork_cont, Waux_mat, dWaux_mat, dabdaint) + ENDIF + + END DO + END DO + END DO + + DEALLOCATE (swork_cont) + DEALLOCATE (swork) + + CALL timestop(handle) + + END SUBROUTINE int_overlap_aba_shg + ! ************************************************************************************************** !> \brief precalculates the angular part of the SHG integrals for the matrices !> (fa,fb), (a,b), (a,b,fa) and (b,fb,a); the same Waux_mat can then be used for all diff --git a/src/generic_shg_integrals_init.F b/src/generic_shg_integrals_init.F new file mode 100644 index 0000000000..26ec7e20c3 --- /dev/null +++ b/src/generic_shg_integrals_init.F @@ -0,0 +1,351 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright (C) 2000 - 2016 CP2K developers group ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Initialization for solid harmonic Gaussian (SHG) integral scheme. Scheme for calculation +!> of contracted, spherical Gaussian integrals using the solid harmonics. Initialization of +!> the contraction matrices +!> \par History +!> created [08.2016] +!> \author Dorothea Golze +! ************************************************************************************************** +MODULE generic_shg_integrals_init + USE basis_set_types, ONLY: gto_basis_set_type + USE kinds, ONLY: dp + USE mathconstants, ONLY: fac,& + ifac,& + pi + USE memory_utilities, ONLY: reallocate + USE orbital_pointers, ONLY: indso,& + nsoset + USE spherical_harmonics, ONLY: clebsch_gordon,& + clebsch_gordon_deallocate,& + clebsch_gordon_init +#include "base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'generic_shg_integrals_init' + + PUBLIC :: contraction_matrix_shg, contraction_matrix_shg_mix, & + get_clebsch_gordon_coefficients + +! ************************************************************************************************** + +CONTAINS + +! ************************************************************************************************** +!> \brief contraction matrix for SHG integrals +!> \param basis ... +!> \param scon_shg contraction matrix +! ************************************************************************************************** + SUBROUTINE contraction_matrix_shg(basis, scon_shg) + + TYPE(gto_basis_set_type), POINTER :: basis + REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scon_shg + + CHARACTER(len=*), PARAMETER :: routineN = 'contraction_matrix_shg', & + routineP = moduleN//':'//routineN + + INTEGER :: ipgf, iset, ishell, l, maxpgf, maxshell, & + nset + INTEGER, DIMENSION(:), POINTER :: npgf, nshell + REAL(KIND=dp) :: aif, gcc + REAL(KIND=dp), DIMENSION(:, :), POINTER :: norm, zet + + nset = basis%nset + npgf => basis%npgf + nshell => basis%nshell + zet => basis%zet + + maxpgf = SIZE(basis%gcc, 1) + maxshell = SIZE(basis%gcc, 2) + ALLOCATE (norm(basis%nset, maxshell)) + ALLOCATE (scon_shg(maxpgf, maxshell, nset)) + scon_shg = 0.0_dp + + CALL basis_norm_shg(basis, norm) + + DO iset = 1, nset + DO ishell = 1, nshell(iset) + l = basis%l(ishell, iset) + DO ipgf = 1, npgf(iset) + aif = 1.0_dp/((2._dp*zet(ipgf, iset))**l) + gcc = basis%gcc(ipgf, ishell, iset) + scon_shg(ipgf, ishell, iset) = norm(iset, ishell)*gcc*aif + END DO + END DO + END DO + + DEALLOCATE (norm) + + END SUBROUTINE contraction_matrix_shg + +!*************************************************************************************************** +!> \brief normalization solid harmonic Gaussians (SHG) +!> \param basis ... +!> \param norm ... +! ************************************************************************************************** + SUBROUTINE basis_norm_shg(basis, norm) + + TYPE(gto_basis_set_type), POINTER :: basis + REAL(dp), DIMENSION(:, :), POINTER :: norm + + CHARACTER(len=*), PARAMETER :: routineN = 'basis_norm_shg', routineP = moduleN//':'//routineN + + INTEGER :: ipgf, iset, ishell, jpgf, l + REAL(KIND=dp) :: aai, aaj, cci, ccj, expa, ppl + + norm = 0._dp + + DO iset = 1, basis%nset + DO ishell = 1, basis%nshell(iset) + l = basis%l(ishell, iset) + expa = 0.5_dp*REAL(2*l+3, dp) + ppl = fac(2*l+2)*pi**(1.5_dp)/fac(l+1) + ppl = ppl/(2._dp**REAL(2*l+1, dp)) + ppl = ppl/REAL(2*l+1, dp) + DO ipgf = 1, basis%npgf(iset) + cci = basis%gcc(ipgf, ishell, iset) + aai = basis%zet(ipgf, iset) + DO jpgf = 1, basis%npgf(iset) + ccj = basis%gcc(jpgf, ishell, iset) + aaj = basis%zet(jpgf, iset) + norm(iset, ishell) = norm(iset, ishell)+cci*ccj*ppl/(aai+aaj)**expa + END DO + END DO + norm(iset, ishell) = 1.0_dp/SQRT(norm(iset, ishell)) + END DO + END DO + + END SUBROUTINE basis_norm_shg + +! ************************************************************************************************** +!> \brief mixed contraction matrix for SHG integrals [aba] and [abb] for orbital and ri basis +!> at the same atom +!> \param orb_basis orbital basis +!> \param ri_basis ... +!> \param orb_index index for orbital basis +!> \param ri_index index for ri basis +!> \param scon_mix mixed contraction matrix +! ************************************************************************************************** + SUBROUTINE contraction_matrix_shg_mix(orb_basis, ri_basis, orb_index, ri_index, scon_mix) + + TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis + INTEGER, DIMENSION(:, :, :), POINTER :: orb_index, ri_index + REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: scon_mix + + CHARACTER(len=*), PARAMETER :: routineN = 'contraction_matrix_shg_mix', & + routineP = moduleN//':'//routineN + + INTEGER :: forb, fri, iil, il, ipgf, iset, ishell, jpgf, jset, jshell, l, l1, l2, lmax_orb, & + lmax_ri, maxpgf_orb, maxpgf_ri, maxshell_orb, maxshell_ri, nf_orb, nf_ri, nl, nl_max, & + nset_orb, nset_ri + INTEGER, DIMENSION(:), POINTER :: npgf_orb, npgf_ri, nshell_orb, nshell_ri + REAL(KIND=dp) :: cjf, const, const1, const2, gcc_orb, & + gcc_ri, prefac, scon1, scon2, & + shg_fac(0:5), zet + REAL(KIND=dp), DIMENSION(:, :), POINTER :: norm_orb, norm_ri, zet_orb, zet_ri + + nset_orb = orb_basis%nset + npgf_orb => orb_basis%npgf + nshell_orb => orb_basis%nshell + zet_orb => orb_basis%zet + nset_ri = ri_basis%nset + npgf_ri => ri_basis%npgf + nshell_ri => ri_basis%nshell + zet_ri => ri_basis%zet + + maxpgf_orb = SIZE(orb_basis%gcc, 1) + maxshell_orb = SIZE(orb_basis%gcc, 2) + ALLOCATE (norm_orb(nset_orb, maxshell_orb)) + maxpgf_ri = SIZE(ri_basis%gcc, 1) + maxshell_ri = SIZE(ri_basis%gcc, 2) + ALLOCATE (norm_ri(nset_ri, maxshell_ri)) + + CALL basis_norm_shg(orb_basis, norm_orb) + CALL basis_norm_shg(ri_basis, norm_ri) + + ALLOCATE (orb_index(maxpgf_orb, maxshell_orb, nset_orb)) + ALLOCATE (ri_index(maxpgf_ri, maxshell_ri, nset_ri)) + + !** index orbital basis set + nf_orb = 0 + DO iset = 1, nset_orb + DO ishell = 1, nshell_orb(iset) + DO ipgf = 1, npgf_orb(iset) + nf_orb = nf_orb+1 + orb_index(ipgf, ishell, iset) = nf_orb + END DO + END DO + END DO + + !** index ri basis set + nf_ri = 0 + DO iset = 1, nset_ri + DO ishell = 1, nshell_ri(iset) + DO ipgf = 1, npgf_ri(iset) + nf_ri = nf_ri+1 + ri_index(ipgf, ishell, iset) = nf_ri + END DO + END DO + END DO + + lmax_orb = MAXVAL(orb_basis%lmax) + lmax_ri = MAXVAL(ri_basis%lmax) + nl_max = INT((lmax_orb+lmax_ri)/2)+1 + ALLOCATE (scon_mix(nl_max, nf_ri, nf_orb, nl_max)) + + shg_fac(0) = 1.0_dp + + DO iset = 1, nset_orb + DO ishell = 1, nshell_orb(iset) + l1 = orb_basis%l(ishell, iset) + const1 = SQRT(1.0_dp/REAL(2*l1+1, dp)) + DO jset = 1, nset_ri + DO jshell = 1, nshell_ri(jset) + l2 = ri_basis%l(jshell, jset) + const2 = SQRT(1.0_dp/REAL(2*l2+1, dp)) + nl = INT((l1+l2)/2) + IF (l1 == 0 .OR. l2 == 0) nl = 0 + DO il = 0, nl + l = l1+l2-2*il + const = const1*const2*2.0_dp*SQRT(pi*REAL(2*l+1, dp)) + DO iil = 1, il + shg_fac(iil) = fac(l+iil-1)*ifac(l)*REAL(l, dp) & + *fac(il)/fac(il-iil)/fac(iil) + ENDDO + DO ipgf = 1, npgf_orb(iset) + forb = orb_index(ipgf, ishell, iset) + gcc_orb = orb_basis%gcc(ipgf, ishell, iset) + scon1 = norm_orb(iset, ishell)*gcc_orb + DO jpgf = 1, npgf_ri(jset) + fri = ri_index(jpgf, jshell, jset) + gcc_ri = ri_basis%gcc(jpgf, jshell, jset) + scon2 = norm_ri(jset, jshell)*gcc_ri + zet = zet_orb(ipgf, iset)+zet_ri(jpgf, jset) + cjf = 1.0_dp/((2._dp*zet)**l) + prefac = const*cjf*scon1*scon2 + DO iil = 0, il + scon_mix(iil+1, fri, forb, il+1) = prefac*shg_fac(iil)/zet**iil + ENDDO + ENDDO + ENDDO + ENDDO + ENDDO + ENDDO + ENDDO + ENDDO + + DEALLOCATE (norm_orb, norm_ri) + + END SUBROUTINE contraction_matrix_shg_mix + +! ************************************************************************************************** +!> \brief calculate the Clebsch-Gordon (CG) coefficients for expansion of the +!> product of two spherical harmonic Gaussians +!> \param my_cg matrix storing CG coefficients +!> \param cg_none0_list list of none-zero CG coefficients +!> \param ncg_none0 number of none-zero CG coefficients +!> \param maxl1 maximal l quantum number of 1st spherical function +!> \param maxl2 maximal l quantum number of 2nd spherical function +! ************************************************************************************************** + SUBROUTINE get_clebsch_gordon_coefficients(my_cg, cg_none0_list, ncg_none0, & + maxl1, maxl2) + + REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: my_cg + INTEGER, DIMENSION(:, :, :), POINTER :: cg_none0_list + INTEGER, DIMENSION(:, :), POINTER :: ncg_none0 + INTEGER, INTENT(IN) :: maxl1, maxl2 + + CHARACTER(len=*), PARAMETER :: routineN = 'get_clebsch_gordon_coefficients', & + routineP = moduleN//':'//routineN + + INTEGER :: il, ilist, iso, iso1, iso2, l1, l1l2, & + l2, lc1, lc2, lp, m1, m2, maxl, mm, & + mp, nlist, nlist_max, nsfunc, nsfunc1, & + nsfunc2 + REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rga + + nlist_max = 6 + nsfunc1 = nsoset(maxl1) + nsfunc2 = nsoset(maxl2) + maxl = maxl1+maxl2 + nsfunc = nsoset(maxl) + + CALL clebsch_gordon_init(maxl) + + ALLOCATE (my_cg(nsfunc1, nsfunc2, nsfunc)) + my_cg = 0.0_dp + ALLOCATE (ncg_none0(nsfunc1, nsfunc2)) + ncg_none0 = 0 + ALLOCATE (cg_none0_list(nsfunc1, nsfunc2, nlist_max)) + cg_none0_list = 0 + + ALLOCATE (rga(maxl, 2)) + rga = 0.0_dp + DO lc1 = 0, maxl1 + DO iso1 = nsoset(lc1-1)+1, nsoset(lc1) + l1 = indso(1, iso1) + m1 = indso(2, iso1) + DO lc2 = 0, maxl2 + DO iso2 = nsoset(lc2-1)+1, nsoset(lc2) + l2 = indso(1, iso2) + m2 = indso(2, iso2) + CALL clebsch_gordon(l1, m1, l2, m2, rga) + l1l2 = l1+l2 + mp = m1+m2 + mm = m1-m2 + IF (m1*m2 < 0 .OR. (m1*m2 == 0 .AND. (m1 < 0 .OR. m2 < 0))) THEN + mp = -ABS(mp) + mm = -ABS(mm) + ELSE + mp = ABS(mp) + mm = ABS(mm) + END IF + DO lp = MOD(l1+l2, 2), l1l2, 2 + il = lp/2+1 + IF (ABS(mp) <= lp) THEN + IF (mp >= 0) THEN + iso = nsoset(lp-1)+lp+1+mp + ELSE + iso = nsoset(lp-1)+lp+1-ABS(mp) + END IF + my_cg(iso1, iso2, iso) = rga(il, 1) + ENDIF + IF (mp /= mm .AND. ABS(mm) <= lp) THEN + IF (mm >= 0) THEN + iso = nsoset(lp-1)+lp+1+mm + ELSE + iso = nsoset(lp-1)+lp+1-ABS(mm) + END IF + my_cg(iso1, iso2, iso) = rga(il, 2) + ENDIF + END DO + nlist = 0 + DO ilist = 1, nsfunc + IF (ABS(my_cg(iso1, iso2, ilist)) > 1.E-8_dp) THEN + nlist = nlist+1 + IF (nlist > nlist_max) THEN + CALL reallocate(cg_none0_list, 1, nsfunc1, 1, nsfunc2, 1, nlist) + nlist_max = nlist + ENDIF + cg_none0_list(iso1, iso2, nlist) = ilist + ENDIF + ENDDO + ncg_none0(iso1, iso2) = nlist + ENDDO ! iso2 + ENDDO ! lc2 + ENDDO ! iso1 + ENDDO ! lc1 + + DEALLOCATE (rga) + CALL clebsch_gordon_deallocate() + + END SUBROUTINE get_clebsch_gordon_coefficients + +END MODULE generic_shg_integrals_init diff --git a/src/lri_environment_init.F b/src/lri_environment_init.F index 7fff0e2a49..5ffe3d9884 100644 --- a/src/lri_environment_init.F +++ b/src/lri_environment_init.F @@ -14,6 +14,9 @@ MODULE lri_environment_init USE atomic_kind_types, ONLY: atomic_kind_type USE basis_set_types, ONLY: copy_gto_basis_set,& gto_basis_set_type + USE generic_shg_integrals_init, ONLY: contraction_matrix_shg,& + contraction_matrix_shg_mix,& + get_clebsch_gordon_coefficients USE input_section_types, ONLY: section_vals_type,& section_vals_val_get USE kinds, ONLY: dp @@ -22,17 +25,11 @@ MODULE lri_environment_init lri_environment_type USE mathconstants, ONLY: fac,& pi - USE memory_utilities, ONLY: reallocate - USE orbital_pointers, ONLY: indso,& - nsoset USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type,& set_qs_env USE qs_kind_types, ONLY: get_qs_kind,& qs_kind_type - USE spherical_harmonics, ONLY: clebsch_gordon,& - clebsch_gordon_deallocate,& - clebsch_gordon_init #include "./base/base_uses.f90" IMPLICIT NONE @@ -131,10 +128,12 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'lri_basis_init', routineP = moduleN//':'//routineN INTEGER :: ikind, nkind + INTEGER, DIMENSION(:, :, :), POINTER :: orb_index, ri_index REAL(KIND=dp), DIMENSION(:), POINTER :: orb_norm_r, ri_int_fbas, ri_norm_r, & ri_norm_s REAL(KIND=dp), DIMENSION(:, :), POINTER :: orb_ovlp, ri_ovlp REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scon_orb, scon_ri + REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: scon_mix TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis IF (ASSOCIATED(lri_env)) THEN @@ -170,6 +169,12 @@ CONTAINS NULLIFY (scon_orb) CALL contraction_matrix_shg(orb_basis, scon_orb) lri_env%bas_prop(ikind)%scon_orb => scon_orb + NULLIFY (scon_mix) + CALL contraction_matrix_shg_mix(orb_basis, ri_basis, & + orb_index, ri_index, scon_mix) + lri_env%bas_prop(ikind)%scon_mix => scon_mix + lri_env%bas_prop(ikind)%orb_index => orb_index + lri_env%bas_prop(ikind)%ri_index => ri_index DEALLOCATE (orb_norm_r, ri_norm_r, ri_norm_s) END IF END DO @@ -178,7 +183,7 @@ CONTAINS END SUBROUTINE lri_basis_init -!***************************************************************************** +! ************************************************************************************************** !> \brief normalization for a contracted Gaussian s-function, !> spherical = cartesian Gaussian for s-functions !> \param basis ... @@ -224,7 +229,7 @@ CONTAINS END SUBROUTINE basis_norm_s_func -!***************************************************************************** +! ************************************************************************************************** !> \brief normalization for radial part of contracted spherical Gaussian !> functions !> \param basis ... @@ -269,45 +274,6 @@ CONTAINS END SUBROUTINE basis_norm_radial -!***************************************************************************** -!> \brief normalization solid harmonic Gaussians (SHG) -!> \param basis ... -!> \param norm ... -! ************************************************************************************************** - SUBROUTINE basis_norm_shg(basis, norm) - - TYPE(gto_basis_set_type), POINTER :: basis - REAL(dp), DIMENSION(:, :), POINTER :: norm - - CHARACTER(len=*), PARAMETER :: routineN = 'basis_norm_shg', routineP = moduleN//':'//routineN - - INTEGER :: ipgf, iset, ishell, jpgf, l - REAL(KIND=dp) :: aai, aaj, cci, ccj, expa, ppl - - norm = 0._dp - - DO iset = 1, basis%nset - DO ishell = 1, basis%nshell(iset) - l = basis%l(ishell, iset) - expa = 0.5_dp*REAL(2*l+3, dp) - ppl = fac(2*l+2)*pi**(1.5_dp)/fac(l+1) - ppl = ppl/(2._dp**REAL(2*l+1, dp)) - ppl = ppl/REAL(2*l+1, dp) - DO ipgf = 1, basis%npgf(iset) - cci = basis%gcc(ipgf, ishell, iset) - aai = basis%zet(ipgf, iset) - DO jpgf = 1, basis%npgf(iset) - ccj = basis%gcc(jpgf, ishell, iset) - aaj = basis%zet(jpgf, iset) - norm(iset, ishell) = norm(iset, ishell)+cci*ccj*ppl/(aai+aaj)**expa - END DO - END DO - norm(iset, ishell) = 1.0_dp/SQRT(norm(iset, ishell)) - END DO - END DO - - END SUBROUTINE basis_norm_shg - !***************************************************************************** !> \brief integral over a single (contracted) lri auxiliary basis function, !> integral is zero for all but s-functions @@ -410,154 +376,4 @@ CONTAINS END SUBROUTINE basis_ovlp -! ************************************************************************************************** -!> \brief contraction matrix for SHG integrals -!> \param basis ... -!> \param scon_shg contraction matrix -! ************************************************************************************************** - SUBROUTINE contraction_matrix_shg(basis, scon_shg) - - TYPE(gto_basis_set_type), POINTER :: basis - REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scon_shg - - CHARACTER(len=*), PARAMETER :: routineN = 'contraction_matrix_shg', & - routineP = moduleN//':'//routineN - - INTEGER :: ipgf, iset, ishell, l, maxpgf, maxshell, & - nset - INTEGER, DIMENSION(:), POINTER :: npgf, nshell - REAL(KIND=dp) :: gcc, norm_w - REAL(KIND=dp), DIMENSION(:, :), POINTER :: norm, zet - - nset = basis%nset - npgf => basis%npgf - nshell => basis%nshell - zet => basis%zet - - maxpgf = SIZE(basis%gcc, 1) - maxshell = SIZE(basis%gcc, 2) - ALLOCATE (norm(basis%nset, maxshell)) - ALLOCATE (scon_shg(maxpgf, maxshell, nset)) - scon_shg = 0.0_dp - - CALL basis_norm_shg(basis, norm) - - DO iset = 1, nset - DO ishell = 1, nshell(iset) - l = basis%l(ishell, iset) - norm_w = norm(iset, ishell) - DO ipgf = 1, npgf(iset) - gcc = basis%gcc(ipgf, ishell, iset) - scon_shg(ipgf, ishell, iset) = norm(iset, ishell)*gcc - END DO - END DO - END DO - - DEALLOCATE (norm) - - END SUBROUTINE contraction_matrix_shg - -!***************************************************************************** -!> \brief calculate the Clebsch-Gordon (CG) coefficients for expansion of the -!> product of two spherical harmonic Gaussians -!> \param my_cg matrix storing CG coefficients -!> \param cg_none0_list list of none-zero CG coefficients -!> \param ncg_none0 number of none-zero CG coefficients -!> \param maxl1 maximal l quantum number of 1st spherical function -!> \param maxl2 maximal l quantum number of 2nd spherical function -! ************************************************************************************************** - SUBROUTINE get_clebsch_gordon_coefficients(my_cg, cg_none0_list, ncg_none0, & - maxl1, maxl2) - - REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: my_cg - INTEGER, DIMENSION(:, :, :), POINTER :: cg_none0_list - INTEGER, DIMENSION(:, :), POINTER :: ncg_none0 - INTEGER, INTENT(IN) :: maxl1, maxl2 - - CHARACTER(len=*), PARAMETER :: routineN = 'get_clebsch_gordon_coefficients', & - routineP = moduleN//':'//routineN - - INTEGER :: il, ilist, iso, iso1, iso2, l1, l1l2, & - l2, lc1, lc2, lp, m1, m2, maxl, mm, & - mp, nlist, nlist_max, nsfunc, nsfunc1, & - nsfunc2 - REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rga - - nlist_max = 6 - nsfunc1 = nsoset(maxl1) - nsfunc2 = nsoset(maxl2) - maxl = maxl1+maxl2 - nsfunc = nsoset(maxl) - - CALL clebsch_gordon_init(maxl) - - ALLOCATE (my_cg(nsfunc1, nsfunc2, nsfunc)) - my_cg = 0.0_dp - ALLOCATE (ncg_none0(nsfunc1, nsfunc2)) - ncg_none0 = 0 - ALLOCATE (cg_none0_list(nsfunc1, nsfunc2, nlist_max)) - cg_none0_list = 0 - - ALLOCATE (rga(maxl, 2)) - rga = 0.0_dp - DO lc1 = 0, maxl1 - DO iso1 = nsoset(lc1-1)+1, nsoset(lc1) - l1 = indso(1, iso1) - m1 = indso(2, iso1) - DO lc2 = 0, maxl2 - DO iso2 = nsoset(lc2-1)+1, nsoset(lc2) - l2 = indso(1, iso2) - m2 = indso(2, iso2) - CALL clebsch_gordon(l1, m1, l2, m2, rga) - l1l2 = l1+l2 - mp = m1+m2 - mm = m1-m2 - IF (m1*m2 < 0 .OR. (m1*m2 == 0 .AND. (m1 < 0 .OR. m2 < 0))) THEN - mp = -ABS(mp) - mm = -ABS(mm) - ELSE - mp = ABS(mp) - mm = ABS(mm) - END IF - DO lp = MOD(l1+l2, 2), l1l2, 2 - il = lp/2+1 - IF (ABS(mp) <= lp) THEN - IF (mp >= 0) THEN - iso = nsoset(lp-1)+lp+1+mp - ELSE - iso = nsoset(lp-1)+lp+1-ABS(mp) - END IF - my_cg(iso1, iso2, iso) = rga(il, 1) - ENDIF - IF (mp /= mm .AND. ABS(mm) <= lp) THEN - IF (mm >= 0) THEN - iso = nsoset(lp-1)+lp+1+mm - ELSE - iso = nsoset(lp-1)+lp+1-ABS(mm) - END IF - my_cg(iso1, iso2, iso) = rga(il, 2) - ENDIF - END DO - nlist = 0 - DO ilist = 1, nsfunc - IF (ABS(my_cg(iso1, iso2, ilist)) > 1.E-8_dp) THEN - nlist = nlist+1 - IF (nlist > nlist_max) THEN - CALL reallocate(cg_none0_list, 1, nsfunc1, 1, nsfunc2, 1, nlist) - nlist_max = nlist - ENDIF - cg_none0_list(iso1, iso2, nlist) = ilist - ENDIF - ENDDO - ncg_none0(iso1, iso2) = nlist - ENDDO ! iso2 - ENDDO ! lc2 - ENDDO ! iso1 - ENDDO ! lc1 - - DEALLOCATE (rga) - CALL clebsch_gordon_deallocate() - - END SUBROUTINE get_clebsch_gordon_coefficients - END MODULE lri_environment_init diff --git a/src/lri_environment_methods.F b/src/lri_environment_methods.F index 081f21cd23..5c44c34b55 100644 --- a/src/lri_environment_methods.F +++ b/src/lri_environment_methods.F @@ -730,13 +730,17 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'lri_int_shg', routineP = moduleN//':'//routineN INTEGER :: handle, nba, nbb, nfa, nfb + INTEGER, DIMENSION(:, :, :), POINTER :: fba_index, fbb_index, oba_index, & + obb_index REAL(KIND=dp) :: dab REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Waux_mat REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dWaux_mat - REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scon_orba, scon_orbb, scon_ria, scon_rib + REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scon_fba, scon_fbb, scon_oba, scon_obb + REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: scona_mix, sconb_mix CALL timeset(routineN, handle) - NULLIFY (scon_orba, scon_orbb, scon_ria, scon_rib) + NULLIFY (scon_oba, scon_obb, scon_fba, scon_fbb, scona_mix, sconb_mix, & + oba_index, obb_index, fba_index, fbb_index) dab = SQRT(SUM(rab*rab)) nba = obasa%nsgf nbb = obasb%nsgf @@ -758,27 +762,35 @@ CONTAINS lrii%dabdaint = 0.0_dp lrii%dabbint = 0.0_dp ELSE - scon_orba => lri_env%bas_prop(ikind)%scon_orb - scon_orbb => lri_env%bas_prop(jkind)%scon_orb - scon_ria => lri_env%bas_prop(ikind)%scon_ri - scon_rib => lri_env%bas_prop(jkind)%scon_ri + scon_oba => lri_env%bas_prop(ikind)%scon_orb + scon_obb => lri_env%bas_prop(jkind)%scon_orb + scon_fba => lri_env%bas_prop(ikind)%scon_ri + scon_fbb => lri_env%bas_prop(jkind)%scon_ri + scona_mix => lri_env%bas_prop(ikind)%scon_mix + sconb_mix => lri_env%bas_prop(jkind)%scon_mix + oba_index => lri_env%bas_prop(ikind)%orb_index + fba_index => lri_env%bas_prop(ikind)%ri_index + obb_index => lri_env%bas_prop(jkind)%orb_index + fbb_index => lri_env%bas_prop(jkind)%ri_index CALL lri_precalc_angular_shg_part(obasa, obasb, fbasa, fbasb, rab, Waux_mat, dWaux_mat, lrii%calc_force_pair) !*** (fa,fb) - CALL int_overlap_ab_shg(lrii%sab, lrii%dsab, rab, fbasa, fbasb, scon_ria, scon_rib, & - Waux_mat, dWaux_mat, lrii%calc_force_pair) + CALL int_overlap_ab_shg(lrii%sab, lrii%dsab, rab, fbasa, fbasb, scon_fba, scon_fbb, & + Waux_mat, dWaux_mat, lrii%calc_force_pair, contraction_high=.FALSE.) !*** (a,b) - CALL int_overlap_ab_shg(lrii%soo, lrii%dsoo, rab, obasa, obasb, scon_orba, scon_orbb, & - Waux_mat, dWaux_mat, lrii%calc_force_pair) + CALL int_overlap_ab_shg(lrii%soo, lrii%dsoo, rab, obasa, obasb, scon_oba, scon_obb, & + Waux_mat, dWaux_mat, lrii%calc_force_pair, contraction_high=.TRUE.) !*** (a,b,fa) CALL int_overlap_aba_shg(lrii%abaint, lrii%dabdaint, rab, obasa, obasb, fbasa, & - scon_orba, scon_orbb, scon_ria, & + scon_obb, scona_mix, oba_index, fba_index, & lri_env%cg_shg%cg_coeff, lri_env%cg_shg%cg_none0_list, & - lri_env%cg_shg%ncg_none0, Waux_mat, dWaux_mat, lrii%calc_force_pair) + lri_env%cg_shg%ncg_none0, & + Waux_mat, dWaux_mat, lrii%calc_force_pair) !*** (a,b,fb) CALL int_overlap_abb_shg(lrii%abbint, lrii%dabbint, rab, obasa, obasb, fbasb, & - scon_orba, scon_orbb, scon_rib, & + scon_oba, sconb_mix, obb_index, fbb_index, & lri_env%cg_shg%cg_coeff, lri_env%cg_shg%cg_none0_list, & - lri_env%cg_shg%ncg_none0, Waux_mat, dWaux_mat, lrii%calc_force_pair) + lri_env%cg_shg%ncg_none0, & + Waux_mat, dWaux_mat, lrii%calc_force_pair) DEALLOCATE (Waux_mat, dWaux_mat) ENDIF diff --git a/src/lri_environment_types.F b/src/lri_environment_types.F index f0bb79abcf..413d9c5fb1 100644 --- a/src/lri_environment_types.F +++ b/src/lri_environment_types.F @@ -139,6 +139,8 @@ MODULE lri_environment_types ! ************************************************************************************************** TYPE lri_bas_type + INTEGER, DIMENSION(:, :, :), POINTER :: orb_index + INTEGER, DIMENSION(:, :, :), POINTER :: ri_index ! integral of ri basis fbas REAL(KIND=dp), DIMENSION(:), POINTER :: int_fbas ! self overlap ri basis @@ -149,6 +151,8 @@ MODULE lri_environment_types REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scon_ri ! contraction matrix for SHG integrals orb basis REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scon_orb + ! contraction matrix for SHG integrals aba/abb + REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: scon_mix END TYPE lri_bas_type ! ************************************************************************************************** @@ -837,6 +841,9 @@ CONTAINS DEALLOCATE (lri_env%bas_prop(i)%orb_ovlp) DEALLOCATE (lri_env%bas_prop(i)%scon_ri) DEALLOCATE (lri_env%bas_prop(i)%scon_orb) + DEALLOCATE (lri_env%bas_prop(i)%orb_index) + DEALLOCATE (lri_env%bas_prop(i)%ri_index) + DEALLOCATE (lri_env%bas_prop(i)%scon_mix) END DO DEALLOCATE (lri_env%bas_prop) END IF diff --git a/src/shg_integrals_test.F b/src/shg_integrals_test.F index 721cd6f256..3d84d43b16 100644 --- a/src/shg_integrals_test.F +++ b/src/shg_integrals_test.F @@ -15,8 +15,18 @@ MODULE shg_integrals_test gto_basis_set_type,& init_orb_basis_set,& read_gto_basis_set - USE generic_os_integrals, ONLY: two_center_os_integrals - USE generic_shg_integrals, ONLY: two_center_shg_integrals + USE generic_os_integrals, ONLY: int_overlap_ab_os,& + int_overlap_aba_os,& + int_overlap_abb_os,& + two_center_os_integrals + USE generic_shg_integrals, ONLY: int_overlap_ab_shg,& + int_overlap_aba_shg,& + int_overlap_abb_shg,& + precalc_angular_shg_part,& + two_center_shg_integrals + USE generic_shg_integrals_init, ONLY: contraction_matrix_shg,& + contraction_matrix_shg_mix,& + get_clebsch_gordon_coefficients USE input_constants, ONLY: operator_coulomb,& operator_gauss,& operator_verf,& @@ -33,8 +43,6 @@ MODULE shg_integrals_test USE input_val_types, ONLY: real_t USE kinds, ONLY: default_string_length,& dp - USE mathconstants, ONLY: fac,& - pi USE orbital_pointers, ONLY: init_orbital_pointers USE orbital_transformation_matrices, ONLY: init_spherical_harmonics #include "./base/base_uses.f90" @@ -105,11 +113,73 @@ CONTAINS CALL keyword_create(keyword, name="CHECK_ACCURACY", & description="Causes abortion when SHG and OS integrals differ "// & - "more than 1.0E-08.", & + "more what's given by ACCURACY_LEVEL.", & default_l_val=.TRUE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, name="ACCURACY_LEVEL", & + description="Level of accuracy for comparison of SHG and OS "// & + "integrals.", & + default_r_val=1.0E-8_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, name="CALCULATE_DERIVATIVES", & + description="Calculates also the derivatives of the integrals.", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, name="TEST_OVERLAP", & + description="Calculates the integrals (a|b).", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_release(keyword) + CALL keyword_create(keyword, name="TEST_COULOMB", & + description="Calculates the integrals (a|1/r12|b).", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, name="TEST_VERF", & + description="Calculates the integrals (a|erf(omega*r12)/r12|b).", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, name="TEST_VERFC", & + description="Calculates the integrals (a|erfc(omega*r12)/r12|b).", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, name="TEST_VGAUSS", & + description="Calculates the integrals (a|exp(omega*r12^2)/r12|b).", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, name="TEST_GAUSS", & + description="Calculates the integrals (a|exp(omega*r12^2)|b).", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, name="TEST_OVERLAP_ABA", & + description="Calculates the integrals (a|b|b).", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, name="TEST_OVERLAP_ABB", & + description="Calculates the integrals (a|b|b).", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + END SUBROUTINE create_shg_integrals_test_section ! ************************************************************************************************** @@ -125,49 +195,77 @@ CONTAINS routineP = moduleN//':'//routineN CHARACTER(LEN=default_string_length) :: basis_type - INTEGER :: count_ab, iab, irep, jab, kab, lamax, & - lbmax, lmax, nab, nab_min, nab_xyz, & - nfa, nfb, nrep, nrep_bas - LOGICAL :: acc_check - REAL(KIND=dp) :: acc_param, ddmax_coulomb, ddmax_gauss, ddmax_verf, ddmax_verfc, & - ddmax_vgauss, ddtemp, dmax_coulomb, dmax_gauss, dmax_verf, dmax_verfc, dmax_vgauss, & - dtemp, omega + INTEGER :: count_ab, handle, iab, jab, kab, lamax, & + lbmax, lcamax, lcbmax, lmax, nab, & + nab_min, nab_xyz, nrep, nrep_bas + LOGICAL :: acc_check, calc_derivatives, & + test_overlap_aba, test_overlap_abb + REAL(KIND=dp) :: acc_param REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rab REAL(KIND=dp), DIMENSION(:), POINTER :: cell_par - REAL(KIND=dp), DIMENSION(:, :), POINTER :: vab_os, vab_shg - REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dvab_os, dvab_shg, scona_shg, sconb_shg - TYPE(gto_basis_set_type), POINTER :: fba, fbb + REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scona_shg, sconb_shg + TYPE(gto_basis_set_type), POINTER :: fba, fbb, oba, obb TYPE(section_vals_type), POINTER :: basis_section - NULLIFY (fba, fbb, basis_section, cell_par) + CALL timeset(routineN, handle) + NULLIFY (oba, obb, fba, fbb, basis_section, cell_par) CALL section_vals_val_get(shg_integrals_test_section, "ABC", r_vals=cell_par) CALL section_vals_val_get(shg_integrals_test_section, "NAB_MIN", i_val=nab_min) CALL section_vals_val_get(shg_integrals_test_section, "NREP", i_val=nrep) CALL section_vals_val_get(shg_integrals_test_section, "CHECK_ACCURACY", l_val=acc_check) + CALL section_vals_val_get(shg_integrals_test_section, "ACCURACY_LEVEL", r_val=acc_param) + CALL section_vals_val_get(shg_integrals_test_section, "CALCULATE_DERIVATIVES", l_val=calc_derivatives) + CALL section_vals_val_get(shg_integrals_test_section, "TEST_OVERLAP_ABA", l_val=test_overlap_aba) + CALL section_vals_val_get(shg_integrals_test_section, "TEST_OVERLAP_ABB", l_val=test_overlap_abb) !*** Read the basis set information basis_section => section_vals_get_subs_vals(shg_integrals_test_section, "BASIS") CALL section_vals_get(basis_section, n_repetition=nrep_bas) - IF (nrep_bas /= 2) RETURN - CALL allocate_gto_basis_set(fba) - CALL read_gto_basis_set(TRIM("A"), basis_type, fba, basis_section, irep=1) - lamax = MAXVAL(fba%lmax) - CALL allocate_gto_basis_set(fbb) - CALL read_gto_basis_set(TRIM("B"), basis_type, fbb, basis_section, irep=2) - lbmax = MAXVAL(fbb%lmax) + IF (.NOT. (nrep_bas == 2 .OR. nrep_bas == 3)) THEN + CALL cp_abort(__LOCATION__, & + "Provide basis sets") + END IF + CALL allocate_gto_basis_set(oba) + CALL read_gto_basis_set(TRIM("A"), basis_type, oba, basis_section, irep=1) + lamax = MAXVAL(oba%lmax) + CALL allocate_gto_basis_set(obb) + CALL read_gto_basis_set(TRIM("B"), basis_type, obb, basis_section, irep=2) + lbmax = MAXVAL(obb%lmax) lmax = MAX(lamax, lbmax) - + IF (test_overlap_aba) THEN + CALL allocate_gto_basis_set(fba) + CALL read_gto_basis_set(TRIM("CA"), basis_type, fba, basis_section, irep=3) + lcamax = MAXVAL(fba%lmax) + lmax = MAX(lamax+lcamax, lbmax) + ENDIF + IF (test_overlap_abb) THEN + CALL allocate_gto_basis_set(fbb) + CALL read_gto_basis_set(TRIM("CB"), basis_type, fbb, basis_section, irep=3) + lcbmax = MAXVAL(fbb%lmax) + lmax = MAX(lamax, lbmax+lcbmax) + ENDIF + IF (test_overlap_aba .AND. test_overlap_abb) THEN + lmax = MAX(MAX(lamax+lcamax, lbmax), MAX(lamax, lbmax+lcbmax)) + ENDIF !*** Initialize basis set information CALL init_orbital_pointers(lmax+1) CALL init_spherical_harmonics(lmax, output_unit=-100) - fba%norm_type = 2 - CALL init_orb_basis_set(fba) - fbb%norm_type = 2 - CALL init_orb_basis_set(fbb) + oba%norm_type = 2 + CALL init_orb_basis_set(oba) + obb%norm_type = 2 + CALL init_orb_basis_set(obb) + IF (test_overlap_aba) THEN + fba%norm_type = 2 + CALL init_orb_basis_set(fba) + ENDIF + IF (test_overlap_abb) THEN + fbb%norm_type = 2 + CALL init_orb_basis_set(fbb) + ENDIF ! if shg integrals are later actually used in the code, contraction_matrix_shg should be ! moved to init_orb_basis_set and scon_shg should become an element of gto_basis_set_type - CALL contraction_matrix_shg(fba, scona_shg) - CALL contraction_matrix_shg(fbb, sconb_shg) + CALL contraction_matrix_shg(oba, scona_shg) + CALL contraction_matrix_shg(obb, sconb_shg) !*** Create range of rab (atomic distances) to be tested nab_xyz = CEILING(REAL(nab_min, KIND=dp)**(1.0_dp/3.0_dp)-1.0E-06) @@ -185,197 +283,453 @@ CONTAINS ENDDO !*** Calculate the SHG integrals - nfa = fba%nsgf - nfb = fbb%nsgf - ALLOCATE (vab_shg(nfa, nfb), dvab_shg(nfa, nfb, 3)) - ALLOCATE (vab_os(nfa, nfb), dvab_os(nfa, nfb, 3)) - omega = 2.3_dp - dmax_coulomb = 0.0_dp - ddmax_coulomb = 0.0_dp - dmax_verf = 0.0_dp - ddmax_verf = 0.0_dp - dmax_verfc = 0.0_dp - ddmax_verfc = 0.0_dp - dmax_vgauss = 0.0_dp - ddmax_vgauss = 0.0_dp - dmax_gauss = 0.0_dp - ddmax_gauss = 0.0_dp - DO irep = 1, nrep - DO iab = 1, nab - !*** Coulomb: (a|1/r12|b) - CALL two_center_shg_integrals(operator_coulomb, vab_shg, dvab_shg, rab(:, iab), & - fba, fbb, scona_shg, sconb_shg, calculate_forces=.TRUE.) - CALL two_center_os_integrals(operator_coulomb, vab_os, dvab_os, rab(:, iab), fba, fbb, & - calculate_forces=.TRUE.) - CALL calculate_deviation(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp) - dmax_coulomb = MAX(dmax_coulomb, dtemp) - ddmax_coulomb = MAX(ddmax_coulomb, ddtemp) - !*** verf: (a|erf(omega*r12)/r12|b) - CALL two_center_shg_integrals(operator_verf, vab_shg, dvab_shg, rab(:, iab), & - fba, fbb, scona_shg, sconb_shg, omega, & - calculate_forces=.TRUE.) - CALL two_center_os_integrals(operator_verf, vab_os, dvab_os, rab(:, iab), fba, fbb, & - omega, calculate_forces=.TRUE.) - CALL calculate_deviation(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp) - dmax_verf = MAX(dmax_verf, dtemp) - ddmax_verf = MAX(ddmax_verf, ddtemp) - !*** verfc: (a|erfc(omega*r12)/r12|b) - CALL two_center_shg_integrals(operator_verfc, vab_shg, dvab_shg, rab(:, iab), & - fba, fbb, scona_shg, sconb_shg, omega, & - calculate_forces=.TRUE.) - CALL two_center_os_integrals(operator_verfc, vab_os, dvab_os, rab(:, iab), fba, fbb, & - omega, calculate_forces=.TRUE.) - CALL calculate_deviation(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp) - dmax_verfc = MAX(dmax_verfc, dtemp) - ddmax_verfc = MAX(ddmax_verfc, ddtemp) - !*** vgauss: (a|exp(omega*r12^2)/r12|b) - CALL two_center_shg_integrals(operator_vgauss, vab_shg, dvab_shg, rab(:, iab), & - fba, fbb, scona_shg, sconb_shg, omega, & - calculate_forces=.TRUE.) - CALL two_center_os_integrals(operator_vgauss, vab_os, dvab_os, rab(:, iab), fba, fbb, & - omega, calculate_forces=.TRUE.) - CALL calculate_deviation(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp) - dmax_vgauss = MAX(dmax_vgauss, dtemp) - ddmax_vgauss = MAX(ddmax_vgauss, ddtemp) - !*** gauss: (a|exp(omega*r12^2)|b) - CALL two_center_shg_integrals(operator_gauss, vab_shg, dvab_shg, rab(:, iab), & - fba, fbb, scona_shg, sconb_shg, omega, & - calculate_forces=.TRUE.) - CALL two_center_os_integrals(operator_gauss, vab_os, dvab_os, rab(:, iab), fba, fbb, & - omega, calculate_forces=.TRUE.) - CALL calculate_deviation(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp) - dmax_gauss = MAX(dmax_gauss, dtemp) - ddmax_gauss = MAX(ddmax_gauss, ddtemp) - ENDDO - ENDDO - IF (iw > 0) THEN - WRITE (iw, FMT="(/,T2,A)") "TEST INFO FOR 2-CENTER SHG and OS INTEGRALS:" - WRITE (iw, FMT="(T2,A)") "Maximal deviation between SHG and OS integrals and their derivatives" - WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|1/r12|b]", & - dmax_coulomb, ddmax_coulomb - WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|erf(omega*r12)/r12|b]", & - dmax_verf, ddmax_verf - WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|erfc(omega*r12)/r12|b]", & - dmax_verfc, ddmax_verfc - WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|exp(-omega*r12^2)/r12|b]", & - dmax_vgauss, ddmax_vgauss - WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|exp(-omega*r12^2)|b]", & - dmax_gauss, ddmax_gauss + CALL test_shg_operator12_integrals(oba, obb, rab, nrep, scona_shg, sconb_shg, & + shg_integrals_test_section, acc_check, & + acc_param, calc_derivatives, iw) - IF (acc_check) THEN - acc_param = 1.0E-08_dp - IF ((dmax_coulomb >= acc_param) .OR. (ddmax_coulomb >= acc_param)) THEN - CPABORT("[a|1/r12|b]: Deviation between SHG and OS larger than 1.0E-08") - ENDIF - IF ((dmax_verf >= acc_param) .OR. (ddmax_verf >= acc_param)) THEN - CPABORT("[a|erf(omega*r12)/r12|b]: Deviation between SHG and OS larger than 1.0E-08") - ENDIF - IF ((dmax_verfc >= acc_param) .OR. (ddmax_verfc >= acc_param)) THEN - CPABORT("[a|erfc(omega*r12)/r12|b]: Deviation between SHG and OS larger than 1.0E-08") - ENDIF - IF ((dmax_vgauss >= acc_param) .OR. (ddmax_vgauss >= acc_param)) THEN - CPABORT("[a|exp(-omega*r12^2)/r12|b]: Dev. between SHG and OS larger than 1.0E-08") - ENDIF - IF ((dmax_gauss >= acc_param) .OR. (ddmax_gauss >= acc_param)) THEN - CPABORT("[a|exp(-omega*r12^2)|b]: Deviation between SHG and OS larger than 1.0E-08") - ENDIF - ENDIF - ENDIF - DEALLOCATE (scona_shg, sconb_shg, rab, vab_shg, dvab_shg, vab_os, dvab_os) - CALL deallocate_gto_basis_set(fba) - CALL deallocate_gto_basis_set(fbb) + CALL test_shg_overlap_integrals(oba, obb, rab, nrep, scona_shg, sconb_shg, & + shg_integrals_test_section, acc_check, & + acc_param, calc_derivatives, iw) + + CALL test_shg_overlap_aba_integrals(oba, obb, fba, fbb, rab, nrep, scona_shg, sconb_shg, & + shg_integrals_test_section, acc_check, & + acc_param, calc_derivatives, iw) + + DEALLOCATE (scona_shg, sconb_shg, rab) + CALL deallocate_gto_basis_set(oba) + CALL deallocate_gto_basis_set(obb) + IF (test_overlap_aba) CALL deallocate_gto_basis_set(fba) + IF (test_overlap_abb) CALL deallocate_gto_basis_set(fbb) + + CALL timestop(handle) END SUBROUTINE shg_integrals_perf_acc_test ! ************************************************************************************************** -!> \brief contraction matrix for the SHG integrals. If these integrals are actually used, -!> this routine should be moved, e.g. in basis_set_types.F, where the basis set -!> information are initialized -!> \param basis ... -!> \param scon_shg contraction matrix +!> \brief ... +!> \param oba ... +!> \param obb ... +!> \param rab ... +!> \param nrep ... +!> \param scona_shg ... +!> \param sconb_shg ... +!> \param shg_integrals_test_section ... +!> \param acc_check ... +!> \param acc_param ... +!> \param calc_derivatives ... +!> \param iw ... ! ************************************************************************************************** - SUBROUTINE contraction_matrix_shg(basis, scon_shg) + SUBROUTINE test_shg_operator12_integrals(oba, obb, rab, nrep, scona_shg, sconb_shg, & + shg_integrals_test_section, acc_check, & + acc_param, calc_derivatives, iw) + TYPE(gto_basis_set_type), POINTER :: oba, obb + REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rab + INTEGER, INTENT(IN) :: nrep + REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scona_shg, sconb_shg + TYPE(section_vals_type), INTENT(IN), POINTER :: shg_integrals_test_section + LOGICAL, INTENT(IN) :: acc_check + REAL(KIND=dp), INTENT(IN) :: acc_param + LOGICAL, INTENT(IN) :: calc_derivatives + INTEGER, INTENT(IN) :: iw - TYPE(gto_basis_set_type), POINTER :: basis - REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scon_shg - - CHARACTER(len=*), PARAMETER :: routineN = 'contraction_matrix_shg', & + CHARACTER(len=*), PARAMETER :: routineN = 'test_shg_operator12_integrals', & routineP = moduleN//':'//routineN - INTEGER :: ipgf, iset, ishell, l, maxpgf, maxshell, & - nset - INTEGER, DIMENSION(:), POINTER :: npgf, nshell - REAL(KIND=dp) :: aif, gcc, norm_w - REAL(KIND=dp), DIMENSION(:, :), POINTER :: norm, zet + INTEGER :: iab, irep, nab, nfa, nfb + LOGICAL :: test_any, test_coulomb, test_gauss, & + test_verf, test_verfc, test_vgauss + REAL(KIND=dp) :: ddmax_coulomb, ddmax_gauss, ddmax_verf, ddmax_verfc, ddmax_vgauss, ddtemp, & + dmax_coulomb, dmax_gauss, dmax_verf, dmax_verfc, dmax_vgauss, dtemp, omega + REAL(KIND=dp), DIMENSION(:, :), POINTER :: vab_os, vab_shg + REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dvab_os, dvab_shg - nset = basis%nset - npgf => basis%npgf - nshell => basis%nshell - zet => basis%zet + CALL section_vals_val_get(shg_integrals_test_section, "TEST_COULOMB", l_val=test_coulomb) + CALL section_vals_val_get(shg_integrals_test_section, "TEST_VERF", l_val=test_verf) + CALL section_vals_val_get(shg_integrals_test_section, "TEST_VERFC", l_val=test_verfc) + CALL section_vals_val_get(shg_integrals_test_section, "TEST_VGAUSS", l_val=test_vgauss) + CALL section_vals_val_get(shg_integrals_test_section, "TEST_GAUSS", l_val=test_gauss) - maxpgf = SIZE(basis%gcc, 1) - maxshell = SIZE(basis%gcc, 2) - ALLOCATE (norm(basis%nset, maxshell)) - ALLOCATE (scon_shg(maxpgf, maxshell, nset)) - scon_shg = 0.0_dp + test_any = (test_coulomb .OR. test_verf .OR. test_verfc .OR. test_vgauss .OR. test_gauss) - CALL basis_norm_shg(basis, norm) + IF (test_any) THEN + nfa = oba%nsgf + nfb = obb%nsgf + ALLOCATE (vab_shg(nfa, nfb), dvab_shg(nfa, nfb, 3)) + ALLOCATE (vab_os(nfa, nfb), dvab_os(nfa, nfb, 3)) + omega = 2.3_dp + dmax_coulomb = 0.0_dp + ddmax_coulomb = 0.0_dp + dmax_verf = 0.0_dp + ddmax_verf = 0.0_dp + dmax_verfc = 0.0_dp + ddmax_verfc = 0.0_dp + dmax_vgauss = 0.0_dp + ddmax_vgauss = 0.0_dp + dmax_gauss = 0.0_dp + ddmax_gauss = 0.0_dp - DO iset = 1, nset - DO ishell = 1, nshell(iset) - l = basis%l(ishell, iset) - norm_w = norm(iset, ishell) - DO ipgf = 1, npgf(iset) - aif = 1.0_dp/((2._dp*zet(ipgf, iset))**l) - gcc = basis%gcc(ipgf, ishell, iset) - scon_shg(ipgf, ishell, iset) = norm(iset, ishell)*gcc*aif - END DO - END DO - END DO + nab = SIZE(rab, 2) + DO irep = 1, nrep + DO iab = 1, nab + !*** Coulomb: (a|1/r12|b) + IF (test_coulomb) THEN + CALL two_center_shg_integrals(operator_coulomb, vab_shg, dvab_shg, rab(:, iab), & + oba, obb, scona_shg, sconb_shg, & + calculate_forces=calc_derivatives) + CALL two_center_os_integrals(operator_coulomb, vab_os, dvab_os, rab(:, iab), oba, obb, & + calculate_forces=calc_derivatives) + CALL calculate_deviation_ab(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp) + dmax_coulomb = MAX(dmax_coulomb, dtemp) + ddmax_coulomb = MAX(ddmax_coulomb, ddtemp) + ENDIF + !*** verf: (a|erf(omega*r12)/r12|b) + IF (test_verf) THEN + CALL two_center_shg_integrals(operator_verf, vab_shg, dvab_shg, rab(:, iab), & + oba, obb, scona_shg, sconb_shg, omega, & + calc_derivatives) + CALL two_center_os_integrals(operator_verf, vab_os, dvab_os, rab(:, iab), oba, obb, & + omega, calc_derivatives) + CALL calculate_deviation_ab(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp) + dmax_verf = MAX(dmax_verf, dtemp) + ddmax_verf = MAX(ddmax_verf, ddtemp) + ENDIF + !*** verfc: (a|erfc(omega*r12)/r12|b) + IF (test_verfc) THEN + CALL two_center_shg_integrals(operator_verfc, vab_shg, dvab_shg, rab(:, iab), & + oba, obb, scona_shg, sconb_shg, omega, & + calc_derivatives) + CALL two_center_os_integrals(operator_verfc, vab_os, dvab_os, rab(:, iab), oba, obb, & + omega, calc_derivatives) + CALL calculate_deviation_ab(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp) + dmax_verfc = MAX(dmax_verfc, dtemp) + ddmax_verfc = MAX(ddmax_verfc, ddtemp) + ENDIF + !*** vgauss: (a|exp(omega*r12^2)/r12|b) + IF (test_vgauss) THEN + CALL two_center_shg_integrals(operator_vgauss, vab_shg, dvab_shg, rab(:, iab), & + oba, obb, scona_shg, sconb_shg, omega, & + calc_derivatives) + CALL two_center_os_integrals(operator_vgauss, vab_os, dvab_os, rab(:, iab), oba, obb, & + omega, calc_derivatives) + CALL calculate_deviation_ab(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp) + dmax_vgauss = MAX(dmax_vgauss, dtemp) + ddmax_vgauss = MAX(ddmax_vgauss, ddtemp) + ENDIF + !*** gauss: (a|exp(omega*r12^2)|b) + IF (test_gauss) THEN + CALL two_center_shg_integrals(operator_gauss, vab_shg, dvab_shg, rab(:, iab), & + oba, obb, scona_shg, sconb_shg, omega, & + calc_derivatives) + CALL two_center_os_integrals(operator_gauss, vab_os, dvab_os, rab(:, iab), oba, obb, & + omega, calc_derivatives) + CALL calculate_deviation_ab(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp) + dmax_gauss = MAX(dmax_gauss, dtemp) + ddmax_gauss = MAX(ddmax_gauss, ddtemp) + ENDIF + ENDDO + ENDDO - DEALLOCATE (norm) + IF (iw > 0) THEN + WRITE (iw, FMT="(/,T2,A)") "TEST INFO FOR 2-CENTER SHG and OS INTEGRALS:" + WRITE (iw, FMT="(T2,A)") "Maximal deviation between SHG and OS integrals and their derivatives" + IF (test_coulomb) THEN + WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|1/r12|b]", & + dmax_coulomb, ddmax_coulomb + ENDIF + IF (test_verf) THEN + WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|erf(omega*r12)/r12|b]", & + dmax_verf, ddmax_verf + ENDIF + IF (test_verfc) THEN + WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|erfc(omega*r12)/r12|b]", & + dmax_verfc, ddmax_verfc + ENDIF + IF (test_vgauss) THEN + WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|exp(-omega*r12^2)/r12|b]", & + dmax_vgauss, ddmax_vgauss + ENDIF + IF (test_gauss) THEN + WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|exp(-omega*r12^2)|b]", & + dmax_gauss, ddmax_gauss + ENDIF - END SUBROUTINE contraction_matrix_shg + IF (acc_check) THEN + IF ((dmax_coulomb >= acc_param) .OR. (ddmax_coulomb >= acc_param)) THEN + CPABORT("[a|1/r12|b]: Deviation between SHG and OS larger than 1.0E-08") + ENDIF + IF ((dmax_verf >= acc_param) .OR. (ddmax_verf >= acc_param)) THEN + CPABORT("[a|erf(omega*r12)/r12|b]: Dev. between SHG and OS larger than 1.0E-08") + ENDIF + IF ((dmax_verfc >= acc_param) .OR. (ddmax_verfc >= acc_param)) THEN + CPABORT("[a|erfc(omega*r12)/r12|b]: Dev. between SHG and OS larger than 1.0E-08") + ENDIF + IF ((dmax_vgauss >= acc_param) .OR. (ddmax_vgauss >= acc_param)) THEN + CPABORT("[a|exp(-omega*r12^2)/r12|b]: Dev. between SHG and OS larger than 1.0E-08") + ENDIF + IF ((dmax_gauss >= acc_param) .OR. (ddmax_gauss >= acc_param)) THEN + CPABORT("[a|exp(-omega*r12^2)|b]: Dev. between SHG and OS larger than 1.0E-08") + ENDIF + ENDIF + ENDIF + DEALLOCATE (vab_shg, vab_os, dvab_shg, dvab_os) + ENDIF + + END SUBROUTINE test_shg_operator12_integrals -!***************************************************************************** -!> \brief normalization solid harmonic Gaussians (SHG) -!> \param basis ... -!> \param norm ... ! ************************************************************************************************** - SUBROUTINE basis_norm_shg(basis, norm) +!> \brief ... +!> \param oba ... +!> \param obb ... +!> \param rab ... +!> \param nrep ... +!> \param scona_shg ... +!> \param sconb_shg ... +!> \param shg_integrals_test_section ... +!> \param acc_check ... +!> \param acc_param ... +!> \param calc_derivatives ... +!> \param iw ... +! ************************************************************************************************** + SUBROUTINE test_shg_overlap_integrals(oba, obb, rab, nrep, scona_shg, sconb_shg, & + shg_integrals_test_section, acc_check, & + acc_param, calc_derivatives, iw) + TYPE(gto_basis_set_type), POINTER :: oba, obb + REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rab + INTEGER, INTENT(IN) :: nrep + REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scona_shg, sconb_shg + TYPE(section_vals_type), INTENT(IN), POINTER :: shg_integrals_test_section + LOGICAL, INTENT(IN) :: acc_check + REAL(KIND=dp), INTENT(IN) :: acc_param + LOGICAL, INTENT(IN) :: calc_derivatives + INTEGER, INTENT(IN) :: iw - TYPE(gto_basis_set_type), POINTER :: basis - REAL(dp), DIMENSION(:, :), POINTER :: norm + CHARACTER(len=*), PARAMETER :: routineN = 'test_shg_overlap_integrals', & + routineP = moduleN//':'//routineN - CHARACTER(len=*), PARAMETER :: routineN = 'basis_norm_shg', routineP = moduleN//':'//routineN + INTEGER :: iab, irep, la_max, lb_max, nab, nfa, nfb + LOGICAL :: test_overlap + REAL(KIND=dp) :: ddmax_overlap, ddtemp, dmax_overlap, & + dtemp, dummy, ra(3), rb(3) + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Waux_mat + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dWaux_mat + REAL(KIND=dp), DIMENSION(:, :), POINTER :: sab_os, sab_shg + REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dsab_os, dsab_shg - INTEGER :: ipgf, iset, ishell, jpgf, l - REAL(KIND=dp) :: aai, aaj, cci, ccj, expa, ppl + CALL section_vals_val_get(shg_integrals_test_section, "TEST_OVERLAP", & + l_val=test_overlap) + IF (test_overlap) THEN + !effectively switch off screening; makes no sense for the tests + oba%set_radius(:) = HUGE(0.0_dp) + obb%set_radius(:) = HUGE(0.0_dp) + oba%pgf_radius(:, :) = HUGE(0.0_dp) + obb%pgf_radius(:, :) = HUGE(0.0_dp) + nfa = oba%nsgf + nfb = obb%nsgf + la_max = MAXVAL(oba%lmax) + lb_max = MAXVAL(obb%lmax) + ra(:) = 0.0_dp + rb(:) = 0.0_dp + dummy = 0.0_dp + dmax_overlap = 0.0_dp + ddmax_overlap = 0.0_dp + ALLOCATE (sab_shg(nfa, nfb), dsab_shg(nfa, nfb, 3)) + ALLOCATE (sab_os(nfa, nfb), dsab_os(nfa, nfb, 3)) + nab = SIZE(rab, 2) + DO irep = 1, nrep + DO iab = 1, nab + CALL precalc_angular_shg_part(la_max, lb_max, rab(:, iab), Waux_mat, dWaux_mat, & + calc_derivatives) + CALL int_overlap_ab_shg(sab_shg, dsab_shg, rab(:, iab), oba, obb, scona_shg, & + sconb_shg, Waux_mat, dWaux_mat, calc_derivatives, & + contraction_high=.TRUE.) + CALL int_overlap_ab_os(sab_os, dsab_os, ra, rb, rab(:, iab), oba, obb, & + calc_derivatives, debug=.FALSE., dmax=dummy) + CALL calculate_deviation_ab(sab_shg, sab_os, dsab_shg, dsab_os, dtemp, ddtemp) + dmax_overlap = MAX(dmax_overlap, dtemp) + ddmax_overlap = MAX(ddmax_overlap, ddtemp) + ENDDO + ENDDO - norm = 0._dp + IF (iw > 0) THEN + WRITE (iw, FMT="(/,T2,A)") "TEST INFO FOR 2-CENTER OVERLAP SHG and OS INTEGRALS:" + WRITE (iw, FMT="(T2,A)") "Maximal deviation between SHG and OS integrals and their derivatives" + WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|b]", & + dmax_overlap, ddmax_overlap + ENDIF + IF (acc_check) THEN + IF ((dmax_overlap >= acc_param) .OR. (ddmax_overlap >= acc_param)) THEN + CPABORT("[a|b]: Deviation between SHG and OS larger than 1.0E-08") + ENDIF + ENDIF + DEALLOCATE (Waux_mat, dWaux_mat) + DEALLOCATE (sab_shg, sab_os, dsab_shg, dsab_os) + ENDIF - DO iset = 1, basis%nset - DO ishell = 1, basis%nshell(iset) - l = basis%l(ishell, iset) - expa = 0.5_dp*REAL(2*l+3, dp) - ppl = fac(2*l+2)*pi**(1.5_dp)/fac(l+1) - ppl = ppl/(2._dp**REAL(2*l+1, dp)) - ppl = ppl/REAL(2*l+1, dp) - DO ipgf = 1, basis%npgf(iset) - cci = basis%gcc(ipgf, ishell, iset) - aai = basis%zet(ipgf, iset) - DO jpgf = 1, basis%npgf(iset) - ccj = basis%gcc(jpgf, ishell, iset) - aaj = basis%zet(jpgf, iset) - norm(iset, ishell) = norm(iset, ishell)+cci*ccj*ppl/(aai+aaj)**expa - END DO - END DO - norm(iset, ishell) = 1.0_dp/SQRT(norm(iset, ishell)) - END DO - END DO + END SUBROUTINE test_shg_overlap_integrals - END SUBROUTINE basis_norm_shg +! ************************************************************************************************** +!> \brief ... +!> \param oba ... +!> \param obb ... +!> \param fba ... +!> \param fbb ... +!> \param rab ... +!> \param nrep ... +!> \param scon_oba ... +!> \param scon_obb ... +!> \param shg_integrals_test_section ... +!> \param acc_check ... +!> \param acc_param ... +!> \param calc_derivatives ... +!> \param iw ... +! ************************************************************************************************** + SUBROUTINE test_shg_overlap_aba_integrals(oba, obb, fba, fbb, rab, nrep, scon_oba, scon_obb, & + shg_integrals_test_section, acc_check, & + acc_param, calc_derivatives, iw) + TYPE(gto_basis_set_type), POINTER :: oba, obb, fba, fbb + REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rab + INTEGER, INTENT(IN) :: nrep + REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scon_oba, scon_obb + TYPE(section_vals_type), INTENT(IN), POINTER :: shg_integrals_test_section + LOGICAL, INTENT(IN) :: acc_check + REAL(KIND=dp), INTENT(IN) :: acc_param + LOGICAL, INTENT(IN) :: calc_derivatives + INTEGER, INTENT(IN) :: iw + + CHARACTER(len=*), PARAMETER :: routineN = 'test_shg_overlap_aba_integrals', & + routineP = moduleN//':'//routineN + + INTEGER :: iab, irep, la_max, laa_max, lb_max, & + lbb_max, maxl_orb, maxl_ri, nab, nba, & + nbb, nfa, nfb + INTEGER, DIMENSION(:, :), POINTER :: ncg_none0 + INTEGER, DIMENSION(:, :, :), POINTER :: cg_none0_list, fba_index, fbb_index, & + oba_index, obb_index + LOGICAL :: test_overlap_aba, test_overlap_abb + REAL(KIND=dp) :: ddmax_overlap_aba, ddmax_overlap_abb, & + ddtemp, dmax_overlap_aba, & + dmax_overlap_abb, dtemp, dummy, ra(3), & + rb(3) + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Waux_mat + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dWaux_mat + REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: cg_coeff, saba_os, saba_shg, sabb_os, & + sabb_shg + REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: dsaba_os, dsaba_shg, dsabb_os, & + dsabb_shg, scona_mix, sconb_mix + + CALL section_vals_val_get(shg_integrals_test_section, "TEST_OVERLAP_ABA", & + l_val=test_overlap_aba) + CALL section_vals_val_get(shg_integrals_test_section, "TEST_OVERLAP_ABB", & + l_val=test_overlap_abb) + IF (test_overlap_aba .OR. test_overlap_abb) THEN + !effectively switch off screening; makes no sense for the tests + oba%set_radius(:) = HUGE(0.0_dp) + obb%set_radius(:) = HUGE(0.0_dp) + oba%pgf_radius(:, :) = HUGE(0.0_dp) + obb%pgf_radius(:, :) = HUGE(0.0_dp) + nba = oba%nsgf + nbb = obb%nsgf + maxl_orb = MAX(MAXVAL(oba%lmax), MAXVAL(obb%lmax)) + la_max = MAXVAL(oba%lmax) + lb_max = MAXVAL(obb%lmax) + IF (test_overlap_aba) THEN + fba%set_radius(:) = HUGE(0.0_dp) + fba%pgf_radius(:, :) = HUGE(0.0_dp) + nfa = fba%nsgf + maxl_ri = MAXVAL(fba%lmax) + laa_max = MAXVAL(oba%lmax)+MAXVAL(fba%lmax) + ALLOCATE (saba_shg(nba, nbb, nfa), dsaba_shg(nba, nbb, nfa, 3)) + ALLOCATE (saba_os(nba, nbb, nfa), dsaba_os(nba, nbb, nfa, 3)) + CALL contraction_matrix_shg_mix(oba, fba, oba_index, fba_index, scona_mix) + ENDIF + IF (test_overlap_abb) THEN + fbb%set_radius(:) = HUGE(0.0_dp) + fbb%pgf_radius(:, :) = HUGE(0.0_dp) + nfb = fbb%nsgf + maxl_ri = MAXVAL(fbb%lmax) + lbb_max = MAXVAL(obb%lmax)+MAXVAL(fbb%lmax) + ALLOCATE (sabb_shg(nba, nbb, nfb), dsabb_shg(nba, nbb, nfb, 3)) + ALLOCATE (sabb_os(nba, nbb, nfb), dsabb_os(nba, nbb, nfb, 3)) + CALL contraction_matrix_shg_mix(obb, fbb, obb_index, fbb_index, sconb_mix) + ENDIF + ra(:) = 0.0_dp + rb(:) = 0.0_dp + dummy = 0.0_dp + dmax_overlap_aba = 0.0_dp + ddmax_overlap_aba = 0.0_dp + dmax_overlap_abb = 0.0_dp + ddmax_overlap_abb = 0.0_dp + CALL get_clebsch_gordon_coefficients(cg_coeff, cg_none0_list, ncg_none0, maxl_orb, maxl_ri) + nab = SIZE(rab, 2) + IF (test_overlap_aba) THEN + DO irep = 1, nrep + DO iab = 1, nab + saba_shg = 0.0_dp + dsaba_shg = 0.0_dp + CALL precalc_angular_shg_part(laa_max, lb_max, rab(:, iab), Waux_mat, dWaux_mat, & + calc_derivatives) + CALL int_overlap_aba_shg(saba_shg, dsaba_shg, rab(:, iab), oba, obb, fba, & + scon_obb, scona_mix, oba_index, fba_index, & + cg_coeff, cg_none0_list, ncg_none0, & + Waux_mat, dWaux_mat, calc_derivatives) + CALL int_overlap_aba_os(saba_os, dsaba_os, ra, rb, rab(:, iab), oba, obb, fba, & + calc_derivatives, debug=.FALSE., dmax=dummy) + CALL calculate_deviation_abx(saba_shg, saba_os, dsaba_shg, dsaba_os, dtemp, ddtemp) + dmax_overlap_aba = MAX(dmax_overlap_aba, dtemp) + ddmax_overlap_aba = MAX(ddmax_overlap_aba, ddtemp) + DEALLOCATE (Waux_mat, dWaux_mat) + ENDDO + ENDDO + DEALLOCATE (oba_index, fba_index, scona_mix) + DEALLOCATE (saba_shg, saba_os, dsaba_shg, dsaba_os) + ENDIF + IF (test_overlap_abb) THEN + DO irep = 1, nrep + DO iab = 1, nab + sabb_shg = 0.0_dp + dsabb_shg = 0.0_dp + CALL precalc_angular_shg_part(lbb_max, la_max, rab(:, iab), Waux_mat, dWaux_mat, & + calc_derivatives) + CALL int_overlap_abb_shg(sabb_shg, dsabb_shg, rab(:, iab), oba, obb, fbb, & + scon_oba, sconb_mix, obb_index, fbb_index, & + cg_coeff, cg_none0_list, ncg_none0, & + Waux_mat, dWaux_mat, calc_derivatives) + CALL int_overlap_abb_os(sabb_os, dsabb_os, ra, rb, rab(:, iab), oba, obb, fbb, & + calc_derivatives, debug=.FALSE., dmax=dummy) + CALL calculate_deviation_abx(sabb_shg, sabb_os, dsabb_shg, dsabb_os, dtemp, ddtemp) + dmax_overlap_abb = MAX(dmax_overlap_abb, dtemp) + ddmax_overlap_abb = MAX(ddmax_overlap_abb, ddtemp) + DEALLOCATE (Waux_mat, dWaux_mat) + ENDDO + ENDDO + DEALLOCATE (obb_index, fbb_index, sconb_mix) + DEALLOCATE (sabb_shg, sabb_os, dsabb_shg, dsabb_os) + ENDIF + IF (iw > 0) THEN + WRITE (iw, FMT="(/,T2,A)") "TEST INFO [a|b|x] OVERLAP SHG and OS INTEGRALS:" + WRITE (iw, FMT="(T2,A)") "Maximal deviation between SHG and OS integrals and their derivatives" + IF (test_overlap_aba) THEN + WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|b|a]", & + dmax_overlap_aba, ddmax_overlap_aba + ENDIF + IF (test_overlap_abb) THEN + WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|b|b]", & + dmax_overlap_abb, ddmax_overlap_abb + ENDIF + ENDIF + IF (acc_check) THEN + IF ((dmax_overlap_aba >= acc_param) .OR. (ddmax_overlap_aba >= acc_param)) THEN + CPABORT("[a|b|a]: Deviation between SHG and OS larger than 1.0E-08") + ENDIF + IF ((dmax_overlap_abb >= acc_param) .OR. (ddmax_overlap_abb >= acc_param)) THEN + CPABORT("[a|b|b]: Deviation between SHG and OS larger than 1.0E-08") + ENDIF + ENDIF + DEALLOCATE (cg_coeff, cg_none0_list, ncg_none0) + ENDIF + + END SUBROUTINE test_shg_overlap_aba_integrals ! ************************************************************************************************** !> \brief Calculation of the deviation between SHG and OS integrals @@ -386,13 +740,13 @@ CONTAINS !> \param dmax maximal deviation of vab matrices !> \param ddmax maximal deviation of dvab matrices ! ************************************************************************************************** - SUBROUTINE calculate_deviation(vab_shg, vab_os, dvab_shg, dvab_os, dmax, ddmax) + SUBROUTINE calculate_deviation_ab(vab_shg, vab_os, dvab_shg, dvab_os, dmax, ddmax) REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: vab_shg, vab_os REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: dvab_shg, dvab_os REAL(KIND=dp), INTENT(OUT) :: dmax, ddmax - CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_deviation', & + CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_deviation_ab', & routineP = moduleN//':'//routineN INTEGER :: i, j, k @@ -419,5 +773,53 @@ CONTAINS ENDDO ENDDO - END SUBROUTINE + END SUBROUTINE calculate_deviation_ab + +! ************************************************************************************************** +!> \brief Calculation of the deviation between SHG and OS integrals +!> \param vab_shg integral matrix obtained from the SHG scheme +!> \param vab_os integral matrix obtained from the OS scheme +!> \param dvab_shg derivative of the integrals, SHG +!> \param dvab_os derivative of the integrals, OS +!> \param dmax maximal deviation of vab matrices +!> \param ddmax maximal deviation of dvab matrices +! ************************************************************************************************** + SUBROUTINE calculate_deviation_abx(vab_shg, vab_os, dvab_shg, dvab_os, dmax, ddmax) + + REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: vab_shg, vab_os + REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dvab_shg, dvab_os + REAL(KIND=dp), INTENT(OUT) :: dmax, ddmax + + CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_deviation_abx', & + routineP = moduleN//':'//routineN + + INTEGER :: i, j, k, l + REAL(KIND=dp) :: diff + + dmax = 0.0_dp + ddmax = 0.0_dp + + ! integrals vab + DO k = 1, SIZE(vab_shg, 3) + DO j = 1, SIZE(vab_shg, 2) + DO i = 1, SIZE(vab_shg, 1) + diff = ABS(vab_shg(i, j, k)-vab_os(i, j, k)) + dmax = MAX(dmax, diff) + ENDDO + ENDDO + ENDDO + + ! derivatives dvab + DO l = 1, 3 + DO k = 1, SIZE(dvab_shg, 3) + DO j = 1, SIZE(dvab_shg, 2) + DO i = 1, SIZE(dvab_shg, 1) + diff = ABS(dvab_shg(i, j, k, l)-dvab_os(i, j, k, l)) + ddmax = MAX(ddmax, diff) + ENDDO + ENDDO + ENDDO + ENDDO + + END SUBROUTINE calculate_deviation_abx END MODULE shg_integrals_test diff --git a/tests/LIBTEST/TEST_FILES b/tests/LIBTEST/TEST_FILES index aaeb762fa9..01d11496ce 100644 --- a/tests/LIBTEST/TEST_FILES +++ b/tests/LIBTEST/TEST_FILES @@ -30,4 +30,5 @@ test_eri_mme_accuracy.inp 0 test_eri_mme_performance.inp 0 test_minimax.inp 0 test_shg_integrals_01.inp 0 +test_shg_integrals_02.inp 0 #EOF diff --git a/tests/LIBTEST/test_shg_integrals_01.inp b/tests/LIBTEST/test_shg_integrals_01.inp index 2cfc805a75..11db85edd9 100644 --- a/tests/LIBTEST/test_shg_integrals_01.inp +++ b/tests/LIBTEST/test_shg_integrals_01.inp @@ -1,5 +1,5 @@ &GLOBAL - PROJECT test_shg_integrals + PROJECT test_shg_integrals_ab PRINT_LEVEL MEDIUM PROGRAM_NAME TEST RUN_TYPE NONE @@ -9,7 +9,14 @@ ABC 10.0 10.0 10.0 NAB_MIN 8 NREP 2 - CHECK_ACCURACY + CHECK_ACCURACY + CALCULATE_DERIVATIVES + TEST_OVERLAP + TEST_COULOMB + TEST_VERF + TEST_VERFC + TEST_VGAUSS + TEST_GAUSS &BASIS 10 2 0 0 1 1 diff --git a/tests/LIBTEST/test_shg_integrals_02.inp b/tests/LIBTEST/test_shg_integrals_02.inp new file mode 100644 index 0000000000..8eecde5e44 --- /dev/null +++ b/tests/LIBTEST/test_shg_integrals_02.inp @@ -0,0 +1,72 @@ +&GLOBAL + PROJECT test_shg_integrals_abx + PRINT_LEVEL MEDIUM + PROGRAM_NAME TEST + RUN_TYPE NONE +&END GLOBAL +&TEST + &SHG_INTEGRALS_TEST + ABC 5.0 5.0 5.0 + NAB_MIN 4 + NREP 1 + CHECK_ACCURACY + CALCULATE_DERIVATIVES + TEST_OVERLAP_ABA + TEST_OVERLAP_ABB + &BASIS + 1 + 2 0 1 7 2 1 + 11.478000339908 0.024916243200 -0.012512421400 0.024510918200 + 3.700758562763 0.079825490000 -0.056449071100 0.058140794100 + 1.446884268432 0.128862675300 0.011242684700 0.444709498500 + 0.716814589696 0.379448894600 -0.418587548300 0.646207973100 + 0.247918564176 0.324552432600 0.590363216700 0.803385018200 + 0.066918004004 0.037148121400 0.438703133000 0.892971208700 + 0.021708243634 -0.001125195500 -0.059693171300 0.120101316500 + &END + &BASIS + 1 + 2 0 2 7 2 2 1 + 12.015954705512 -0.060190841200 0.065738617900 0.036543638800 -0.034210557400 0.014807054400 + 5.108150287385 -0.129597923300 0.110885902200 0.120927648700 -0.120619770900 0.068186159300 + 2.048398039874 0.118175889400 -0.053732406400 0.251093670300 -0.213719464600 0.290576499200 + 0.832381575582 0.462964485000 -0.572670666200 0.352639910300 -0.473674858400 1.063344189500 + 0.352316246455 0.450353782600 0.186760006700 0.294708645200 0.484848376400 0.307656114200 + 0.142977330880 0.092715833600 0.387201458600 0.173039869300 0.717465919700 0.318346834400 + 0.046760918300 -0.000255945800 0.003825849600 0.009726110600 0.032498979400 -0.005771736600 + &END + &BASIS + 15 + 2 0 0 1 1 + 24.031909411024 1.0 + 2 0 0 1 1 + 16.167926705922 1.0 + 2 0 1 1 1 1 + 10.877281929506 1.0 1.0 + 2 0 1 1 1 1 + 7.317899463920 1.0 1.0 + 2 0 2 1 1 1 1 + 4.923256831173 1.0 1.0 1.0 + 2 0 3 1 1 1 1 1 + 3.312215198528 1.0 1.0 1.0 1.0 + 2 0 3 1 1 1 1 1 + 2.228356126354 1.0 1.0 1.0 1.0 + 2 0 4 1 1 1 1 1 1 + 1.499169204968 1.0 1.0 1.0 1.0 1.0 + 2 0 4 1 1 1 1 1 1 + 1.008594756710 1.0 1.0 1.0 1.0 1.0 + 2 0 4 1 1 1 1 1 1 + 0.678551413605 1.0 1.0 1.0 1.0 1.0 + 2 0 4 1 1 1 1 1 1 + 0.456508441910 1.0 1.0 1.0 1.0 1.0 + 2 0 4 1 1 1 1 1 1 + 0.307124785767 1.0 1.0 1.0 1.0 1.0 + 2 0 4 1 1 1 1 1 1 + 0.206624073890 1.0 1.0 1.0 1.0 1.0 + 2 0 4 1 1 1 1 1 1 + 0.139010297734 1.0 1.0 1.0 1.0 1.0 + 2 0 4 1 1 1 1 1 1 + 0.093521836600 1.0 1.0 1.0 1.0 1.0 + &END + &END SHG_INTEGRALS_TEST +&END TEST