comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 14:15:19 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

SHG integrals: speed-up of [aba] and [abb] integrals + inclusion in libtest

Browse source 
svn-origin-rev: 17197
This commit is contained in:
Dorothea Golze 2016-08-05 21:04:08 +00:00
parent ca28d9b427
commit 889ab9ec78
10 changed files with 1551 additions and 843 deletions
Download patch file Download diff file Expand all files Collapse all files

466
src/aobasis/ai_s_contract_shg.F
View file

@ -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

434
src/generic_shg_integrals.F
View file

@ -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

351
src/generic_shg_integrals_init.F Normal file
View file

@ -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

210
src/lri_environment_init.F
View file

@ -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

40
src/lri_environment_methods.F
View file

@ -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

7
src/lri_environment_types.F
View file

@ -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

802
src/shg_integrals_test.F
View file

@ -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

1
tests/LIBTEST/TEST_FILES
View file

@ -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

11
tests/LIBTEST/test_shg_integrals_01.inp
View file

@ -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

72
tests/LIBTEST/test_shg_integrals_02.inp Normal file
View file

@ -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