(git:58e3e09)
s_contract_shg.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Routines for calculating the s-integrals and their scalar derivatives with respect to rab2
10 !> over solid harmonic Gaussian (SHG) functions + contraction routines for these integrals
11 !> i) (s|O(r12)|s) where O(r12) is the overlap, coulomb operator etc.
12 !> ii) (aba) and (abb) s-overlaps
13 !> \par Literature (partly)
14 !> T.J. Giese and D. M. York, J. Chem. Phys, 128, 064104 (2008)
15 !> T. Helgaker, P Joergensen, J. Olsen, Molecular Electronic-Structure
16 !> Theory, Wiley
17 !> R. Ahlrichs, PCCP, 8, 3072 (2006)
18 !> \par History
19 !> created [05.2016]
20 !> \author Dorothea Golze
21 ! **************************************************************************************************
23  USE gamma, ONLY: fgamma => fgamma_0
24  USE kinds, ONLY: dp
25  USE mathconstants, ONLY: dfac,&
26  fac,&
27  pi
28 #include "../base/base_uses.f90"
29 
30  IMPLICIT NONE
31 
32  PRIVATE
33 
34  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 's_contract_shg'
35 
36 ! *** Public subroutines ***
39 
42 
43 CONTAINS
44 
45 ! **************************************************************************************************
46 !> \brief calculates the uncontracted, not normalized [s|s] overlap
47 !> \param la_max maximal l quantum number on a
48 !> \param npgfa number of primitive Gaussian on a
49 !> \param zeta set of exponents on a
50 !> \param lb_max maximal l quantum number on b
51 !> \param npgfb number of primitive Gaussian on a
52 !> \param zetb set of exponents on a
53 !> \param rab distance vector between a and b
54 !> \param s uncontracted overlap of s functions
55 !> \param calculate_forces ...
56 ! **************************************************************************************************
57  SUBROUTINE s_overlap_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, rab, s, calculate_forces)
58 
59  INTEGER, INTENT(IN) :: la_max, npgfa
60  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta
61  INTEGER, INTENT(IN) :: lb_max, npgfb
62  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb
63  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
64  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: s
65  LOGICAL, INTENT(IN) :: calculate_forces
66 
67  INTEGER :: ids, ipgfa, jpgfb, ndev
68  REAL(kind=dp) :: a, b, rab2, xhi, zet
69 
70  ! Distance of the centers a and b
71  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
72  ndev = 0
73  IF (calculate_forces) ndev = 1
74  ! Loops over all pairs of primitive Gaussian-type functions
75  DO ipgfa = 1, npgfa
76  DO jpgfb = 1, npgfb
77 
78  ! Distance Screening !maybe later
79 
80  ! Calculate some prefactors
81  a = zeta(ipgfa)
82  b = zetb(jpgfb)
83  zet = a + b
84  xhi = a*b/zet
85 
86  ! [s|s] integral
87  s(ipgfa, jpgfb, 1) = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
88 
89  DO ids = 2, la_max + lb_max + ndev + 1
90  s(ipgfa, jpgfb, ids) = -xhi*s(ipgfa, jpgfb, ids - 1)
91  END DO
92 
93  END DO
94  END DO
95 
96  END SUBROUTINE s_overlap_ab
97 
98 ! **************************************************************************************************
99 !> \brief calculates [s|ra^n|s] integrals for [aba] and the [s|rb^n|s]
100 !> integrals for [abb]
101 !> \param la_max maximal l quantum number on a, orbital basis
102 !> \param npgfa number of primitive Gaussian on a, orbital basis
103 !> \param zeta set of exponents on a, orbital basis
104 !> \param lb_max maximal l quantum number on b, orbital basis
105 !> \param npgfb number of primitive Gaussian on a, orbital basis
106 !> \param zetb set of exponents on b, orbital basis
107 !> \param lcb_max maximal l quantum number of aux basis on b
108 !> \param npgfcb number of primitive Gaussian on b. aux basis
109 !> \param zetcb set of exponents on b, aux basis
110 !> \param rab distance vector between a and b
111 !> \param s uncontracted [s|r^n|s] integrals
112 !> \param calculate_forces ...
113 ! **************************************************************************************************
114  SUBROUTINE s_overlap_abb(la_max, npgfa, zeta, lb_max, npgfb, zetb, lcb_max, npgfcb, zetcb, &
115  rab, s, calculate_forces)
116 
117  INTEGER, INTENT(IN) :: la_max, npgfa
118  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta
119  INTEGER, INTENT(IN) :: lb_max, npgfb
120  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb
121  INTEGER, INTENT(IN) :: lcb_max, npgfcb
122  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetcb
123  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
124  REAL(kind=dp), DIMENSION(:, :, :, :, :), &
125  INTENT(INOUT) :: s
126  LOGICAL, INTENT(IN) :: calculate_forces
127 
128  INTEGER :: i, ids, il, ipgfa, j, jpgfb, kpgfb, &
129  lbb_max, lmax, n, ndev, nds, nfac, nl
130  REAL(kind=dp) :: a, b, dfsr_int, exp_rab2, k, pfac, &
131  prefac, rab2, sqrt_pi3, sqrt_zet, &
132  sr_int, temp, xhi, zet
133  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dsr_int, dtemp
134  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: coeff_srs
135 
136  ! Distance of the centers a and b
137  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
138  ndev = 0
139  IF (calculate_forces) ndev = 1
140 
141  lbb_max = lb_max + lcb_max
142  nl = int(lbb_max/2)
143  IF (lb_max == 0 .OR. lcb_max == 0) nl = 0
144  lmax = la_max + lbb_max
145 
146  ALLOCATE (dtemp(nl + 1), dsr_int(nl + 1))
147  ALLOCATE (coeff_srs(nl + 1, nl + 1, nl + 1))
148  IF (nl > 5) CALL get_prefac_sabb(nl, coeff_srs)
149 
150  sqrt_pi3 = sqrt(pi**3)
151 
152  ! Loops over all pairs of primitive Gaussian-type functions
153  DO kpgfb = 1, npgfcb
154  DO jpgfb = 1, npgfb
155  DO ipgfa = 1, npgfa
156 
157  !Calculate some prefactors
158  a = zeta(ipgfa)
159  b = zetb(jpgfb) + zetcb(kpgfb)
160 
161  zet = a + b
162  xhi = a*b/zet
163  exp_rab2 = exp(-xhi*rab2)
164 
165  pfac = a**2/zet
166  sqrt_zet = sqrt(zet)
167 
168  DO il = 0, nl
169  nds = lmax - 2*il + ndev + 1
170  SELECT CASE (il)
171  CASE (0)
172  ! [s|s] integral
173  s(ipgfa, jpgfb, kpgfb, il + 1, 1) = (pi/zet)**(1.5_dp)*exp_rab2
174  DO ids = 2, nds
175  n = ids - 1
176  s(ipgfa, jpgfb, kpgfb, il + 1, ids) = (-xhi)**n*s(ipgfa, jpgfb, kpgfb, il + 1, 1)
177  END DO
178  CASE (1)
179  ![s|r^2|s] integral
180  sr_int = sqrt_pi3/sqrt_zet**5*(3.0_dp + 2.0_dp*pfac*rab2)/2.0_dp
181  s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int
182  k = sqrt_pi3*a**2/sqrt_zet**7
183  DO ids = 2, nds
184  n = ids - 1
185  s(ipgfa, jpgfb, kpgfb, il + 1, ids) = (-xhi)**n*exp_rab2*sr_int &
186  + n*(-xhi)**(n - 1)*k*exp_rab2
187  END DO
188  CASE (2)
189  ![s|r^4|s] integral
190  prefac = sqrt_pi3/4.0_dp/sqrt_zet**7
191  temp = 15.0_dp + 20.0_dp*pfac*rab2 + 4.0_dp*(pfac*rab2)**2
192  sr_int = prefac*temp
193  s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int
194  !** derivatives
195  k = sqrt_pi3*a**4/sqrt_zet**11
196  dsr_int(1) = prefac*(20.0_dp*pfac + 8.0_dp*pfac**2*rab2)
197  DO ids = 2, nds
198  n = ids - 1
199  dtemp(1) = (-xhi)**n*exp_rab2*sr_int
200  dtemp(2) = n*(-xhi)**(n - 1)*exp_rab2*dsr_int(1)
201  dtemp(3) = (n**2 - n)*(-xhi)**(n - 2)*k*exp_rab2
202  s(ipgfa, jpgfb, kpgfb, il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3)
203  END DO
204  CASE (3)
205  ![s|r^6|s] integral
206  prefac = sqrt_pi3/8.0_dp/sqrt_zet**9
207  temp = 105.0_dp + 210.0_dp*pfac*rab2
208  temp = temp + 84.0_dp*(pfac*rab2)**2 + 8.0_dp*(pfac*rab2)**3
209  sr_int = prefac*temp
210  s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int
211  !** derivatives
212  k = sqrt_pi3*a**6/sqrt_zet**15
213  dsr_int(1) = prefac*(210.0_dp*pfac + 168.0_dp*pfac**2*rab2 &
214  + 24.0_dp*pfac**3*rab2**2)
215  dsr_int(2) = prefac*(168.0_dp*pfac**2 + 48.0_dp*pfac**3*rab2)
216  DO ids = 2, nds
217  n = ids - 1
218  dtemp(1) = (-xhi)**n*exp_rab2*sr_int
219  dtemp(2) = real(n, dp)*(-xhi)**(n - 1)*exp_rab2*dsr_int(1)
220  dtemp(3) = real(n**2 - n, dp)/2.0_dp*(-xhi)**(n - 2) &
221  *exp_rab2*dsr_int(2)
222  dtemp(4) = real(n*(n - 1)*(n - 2), dp)*(-xhi)**(n - 3)*k*exp_rab2
223  s(ipgfa, jpgfb, kpgfb, il + 1, ids) = dtemp(1) + dtemp(2) &
224  + dtemp(3) + dtemp(4)
225  END DO
226  CASE (4)
227  ![s|r^8|s] integral
228  prefac = sqrt_pi3/16.0_dp/sqrt_zet**11
229  temp = 945.0_dp + 2520.0_dp*pfac*rab2 + 1512.0_dp*(pfac*rab2)**2
230  temp = temp + 288.0_dp*(pfac*rab2)**3 + 16.0_dp*(pfac*rab2)**4
231  sr_int = prefac*temp
232  s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int
233  !** derivatives
234  k = sqrt_pi3*a**8/sqrt_zet**19
235  dsr_int(1) = 2520.0_dp*pfac + 3024.0_dp*pfac**2*rab2
236  dsr_int(1) = dsr_int(1) + 864.0_dp*pfac**3*rab2**2 &
237  + 64.0_dp*pfac**4*rab2**3
238  dsr_int(1) = prefac*dsr_int(1)
239  dsr_int(2) = 3024.0_dp*pfac**2 + 1728.0_dp*pfac**3*rab2
240  dsr_int(2) = dsr_int(2) + 192.0_dp*pfac**4*rab2**2
241  dsr_int(2) = prefac*dsr_int(2)
242  dsr_int(3) = 1728.0_dp*pfac**3 + 384.0_dp*pfac**4*rab2
243  dsr_int(3) = prefac*dsr_int(3)
244  DO ids = 2, nds
245  n = ids - 1
246  dtemp(1) = (-xhi)**n*exp_rab2*sr_int
247  dtemp(2) = real(n, dp)*(-xhi)**(n - 1)*exp_rab2*dsr_int(1)
248  dtemp(3) = real(n**2 - n, dp)/2.0_dp*(-xhi)**(n - 2) &
249  *exp_rab2*dsr_int(2)
250  dtemp(4) = real(n*(n - 1)*(n - 2), dp)/6.0_dp*(-xhi)**(n - 3) &
251  *exp_rab2*dsr_int(3)
252  dtemp(5) = real(n*(n - 1)*(n - 2)*(n - 3), dp)*(-xhi)**(n - 4) &
253  *k*exp_rab2
254  s(ipgfa, jpgfb, kpgfb, il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) &
255  + dtemp(4) + dtemp(5)
256  END DO
257  CASE (5)
258  ![s|r^10|s] integral
259  prefac = sqrt_pi3/32.0_dp/sqrt_zet**13
260  temp = 10395.0_dp + 34650.0_dp*pfac*rab2
261  temp = temp + 27720.0_dp*(pfac*rab2)**2 + 7920.0_dp*(pfac*rab2)**3
262  temp = temp + 880.0_dp*(pfac*rab2)**4 + 32.0_dp*(pfac*rab2)**5
263  sr_int = prefac*temp
264  s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int
265  !** derivatives
266  k = sqrt_pi3*a**10/sqrt_zet**23
267  dsr_int(1) = 34650.0_dp*pfac + 55440.0_dp*pfac**2*rab2
268  dsr_int(1) = dsr_int(1) + 23760.0_dp*pfac**3*rab2**2
269  dsr_int(1) = dsr_int(1) + 3520.0_dp*pfac**4*rab2**3
270  dsr_int(1) = dsr_int(1) + 160.0_dp*pfac**5*rab2**4
271  dsr_int(1) = prefac*dsr_int(1)
272  dsr_int(2) = 55440.0_dp*pfac**2 + 47520.0_dp*pfac**3*rab2
273  dsr_int(2) = dsr_int(2) + 10560.0_dp*pfac**4*rab2**2
274  dsr_int(2) = dsr_int(2) + 640.0_dp*pfac**5*rab2**3
275  dsr_int(2) = prefac*dsr_int(2)
276  dsr_int(3) = 47520.0_dp*pfac**3 + 21120.0_dp*pfac**4*rab2
277  dsr_int(3) = dsr_int(3) + 1920.0_dp*pfac**5*rab2**2
278  dsr_int(3) = prefac*dsr_int(3)
279  dsr_int(4) = 21120.0_dp*pfac**4 + 3840.0_dp*pfac**5*rab2
280  dsr_int(4) = prefac*dsr_int(4)
281  DO ids = 2, nds
282  n = ids - 1
283  dtemp(1) = (-xhi)**n*exp_rab2*sr_int
284  dtemp(2) = real(n, dp)*(-xhi)**(n - 1)*exp_rab2*dsr_int(1)
285  dtemp(3) = real(n**2 - n, dp)/2.0_dp*(-xhi)**(n - 2) &
286  *exp_rab2*dsr_int(2)
287  dtemp(4) = real(n*(n - 1)*(n - 2), dp)/6.0_dp*(-xhi)**(n - 3) &
288  *exp_rab2*dsr_int(3)
289  dtemp(5) = real(n*(n - 1)*(n - 2)*(n - 3), dp)/24.0_dp*(-xhi)**(n - 4) &
290  *exp_rab2*dsr_int(4)
291  dtemp(6) = real(n*(n - 1)*(n - 2)*(n - 3)*(n - 4), dp)*(-xhi)**(n - 5) &
292  *k*exp_rab2
293  s(ipgfa, jpgfb, kpgfb, il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) &
294  + dtemp(4) + dtemp(5) + dtemp(6)
295  END DO
296  CASE DEFAULT
297  !*** general formula; factor 1.5-2 slower than explicit expressions
298  prefac = exp_rab2/sqrt_zet**(2*il + 3)
299  sr_int = 0.0_dp
300  DO i = 0, il
301  sr_int = sr_int + coeff_srs(i + 1, 1, il + 1)*(pfac)**i*rab2**i
302  END DO
303  s(ipgfa, jpgfb, kpgfb, il + 1, 1) = prefac*sr_int
304  DO ids = 2, nds
305  n = ids - 1
306  nfac = 1
307  dfsr_int = (-xhi)**n*sr_int
308  DO j = 1, il
309  temp = 0.0_dp
310  DO i = j, il
311  temp = temp + coeff_srs(i + 1, j + 1, il + 1)*(pfac)**i*rab2**(i - j)
312  END DO
313  nfac = nfac*(n - j + 1)
314  dfsr_int = dfsr_int + temp*real(nfac, dp)/fac(j)*(-xhi)**(n - j)
315  END DO
316  s(ipgfa, jpgfb, kpgfb, il + 1, ids) = prefac*dfsr_int
317  END DO
318 
319  END SELECT
320 
321  END DO
322 
323  END DO
324  END DO
325  END DO
326 
327  DEALLOCATE (dtemp, dsr_int)
328  DEALLOCATE (coeff_srs)
329 
330  END SUBROUTINE s_overlap_abb
331 
332 ! **************************************************************************************************
333 !> \brief calculates the uncontracted, not normalized [0a|ra^(2m)|0b] two-center integral,
334 !> where ra = r-Ra and Ra center of a
335 !> \param la_max maximal l quantum number on a
336 !> \param npgfa number of primitive Gaussian on a
337 !> \param zeta set of exponents on a
338 !> \param lb_max maximal l quantum number on b
339 !> \param npgfb number of primitive Gaussian on a
340 !> \param zetb set of exponents on a
341 !> \param m exponent of the ra operator
342 !> \param rab distance vector between a and b
343 !> \param s ...
344 !> \param calculate_forces ...
345 ! **************************************************************************************************
346  SUBROUTINE s_ra2m_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, m, rab, s, calculate_forces)
347 
348  INTEGER, INTENT(IN) :: la_max, npgfa
349  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta
350  INTEGER, INTENT(IN) :: lb_max, npgfb
351  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb
352  INTEGER, INTENT(IN) :: m
353  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
354  REAL(kind=dp), DIMENSION(:, :, :, :), &
355  INTENT(INOUT) :: s
356  LOGICAL, INTENT(IN) :: calculate_forces
357 
358  INTEGER :: i, ids, il, ipgfa, j, jpgfb, n, ndev, &
359  nds, nfac
360  REAL(kind=dp) :: a, b, dfsr_int, exp_rab2, k, pfac, &
361  prefac, rab2, sqrt_pi3, sqrt_zet, &
362  sr_int, temp, xhi, zet
363  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dsr_int, dtemp
364  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: coeff_srs
365 
366  ! Distance of the centers a and b
367  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
368  ndev = 0
369  IF (calculate_forces) ndev = 1
370 
371  ALLOCATE (dtemp(m + 1), dsr_int(m + 1))
372  ALLOCATE (coeff_srs(m + 1, m + 1, m + 1))
373  CALL get_prefac_sabb(m, coeff_srs)
374  !IF(m > 5) CALL get_prefac_sabb(m, coeff_srs)
375  sqrt_pi3 = sqrt(pi**3)
376 
377  ! Loops over all pairs of primitive Gaussian-type functions
378  DO ipgfa = 1, npgfa
379  DO jpgfb = 1, npgfb
380 
381  ! Calculate some prefactors
382  a = zeta(ipgfa)
383  b = zetb(jpgfb)
384  zet = a + b
385  xhi = a*b/zet
386  exp_rab2 = exp(-xhi*rab2)
387 
388  sqrt_zet = sqrt(zet)
389  pfac = b**2/zet
390 
391  nds = la_max + lb_max + ndev + 1
392  DO il = 0, m
393  SELECT CASE (il)
394  CASE (0)
395  ! [s|s] integral
396  s(ipgfa, jpgfb, m - il + 1, 1) = (pi/zet)**(1.5_dp)*exp_rab2
397  DO ids = 2, nds
398  n = ids - 1
399  s(ipgfa, jpgfb, m - il + 1, ids) = (-xhi)**n*s(ipgfa, jpgfb, m - il + 1, 1)
400  END DO
401  CASE (1)
402  ![s|r^2|s] integral
403  sr_int = sqrt_pi3/sqrt_zet**5*(3.0_dp + 2.0_dp*pfac*rab2)/2.0_dp
404  s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int
405  k = sqrt_pi3*b**2/sqrt_zet**7
406  DO ids = 2, nds
407  n = ids - 1
408  s(ipgfa, jpgfb, m - il + 1, ids) = (-xhi)**n*exp_rab2*sr_int &
409  + n*(-xhi)**(n - 1)*k*exp_rab2
410  END DO
411  CASE (2)
412  ![s|r^4|s] integral
413  prefac = sqrt_pi3/4.0_dp/sqrt_zet**7
414  temp = 15.0_dp + 20.0_dp*pfac*rab2 + 4.0_dp*(pfac*rab2)**2
415  sr_int = prefac*temp
416  s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int
417  !** derivatives
418  k = sqrt_pi3*b**4/sqrt_zet**11
419  dsr_int(1) = prefac*(20.0_dp*pfac + 8.0_dp*pfac**2*rab2)
420  DO ids = 2, nds
421  n = ids - 1
422  dtemp(1) = (-xhi)**n*exp_rab2*sr_int
423  dtemp(2) = n*(-xhi)**(n - 1)*exp_rab2*dsr_int(1)
424  dtemp(3) = (n**2 - n)*(-xhi)**(n - 2)*k*exp_rab2
425  s(ipgfa, jpgfb, m - il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3)
426  END DO
427  CASE (3)
428  ![s|r^6|s] integral
429  prefac = sqrt_pi3/8.0_dp/sqrt_zet**9
430  temp = 105.0_dp + 210.0_dp*pfac*rab2
431  temp = temp + 84.0_dp*(pfac*rab2)**2 + 8.0_dp*(pfac*rab2)**3
432  sr_int = prefac*temp
433  s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int
434  !** derivatives
435  k = sqrt_pi3*b**6/sqrt_zet**15
436  dsr_int(1) = prefac*(210.0_dp*pfac + 168.0_dp*pfac**2*rab2 &
437  + 24.0_dp*pfac**3*rab2**2)
438  dsr_int(2) = prefac*(168.0_dp*pfac**2 + 48.0_dp*pfac**3*rab2)
439  DO ids = 2, nds
440  n = ids - 1
441  dtemp(1) = (-xhi)**n*exp_rab2*sr_int
442  dtemp(2) = real(n, dp)*(-xhi)**(n - 1)*exp_rab2*dsr_int(1)
443  dtemp(3) = real(n**2 - n, dp)/2.0_dp*(-xhi)**(n - 2) &
444  *exp_rab2*dsr_int(2)
445  dtemp(4) = real(n*(n - 1)*(n - 2), dp)*(-xhi)**(n - 3)*k*exp_rab2
446  s(ipgfa, jpgfb, m - il + 1, ids) = dtemp(1) + dtemp(2) &
447  + dtemp(3) + dtemp(4)
448  END DO
449  CASE (4)
450  ![s|r^8|s] integral
451  prefac = sqrt_pi3/16.0_dp/sqrt_zet**11
452  temp = 945.0_dp + 2520.0_dp*pfac*rab2 + 1512.0_dp*(pfac*rab2)**2
453  temp = temp + 288.0_dp*(pfac*rab2)**3 + 16.0_dp*(pfac*rab2)**4
454  sr_int = prefac*temp
455  s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int
456  !** derivatives
457  k = sqrt_pi3*b**8/sqrt_zet**19
458  dsr_int(1) = 2520.0_dp*pfac + 3024.0_dp*pfac**2*rab2
459  dsr_int(1) = dsr_int(1) + 864.0_dp*pfac**3*rab2**2 &
460  + 64.0_dp*pfac**4*rab2**3
461  dsr_int(1) = prefac*dsr_int(1)
462  dsr_int(2) = 3024.0_dp*pfac**2 + 1728.0_dp*pfac**3*rab2
463  dsr_int(2) = dsr_int(2) + 192.0_dp*pfac**4*rab2**2
464  dsr_int(2) = prefac*dsr_int(2)
465  dsr_int(3) = 1728.0_dp*pfac**3 + 384.0_dp*pfac**4*rab2
466  dsr_int(3) = prefac*dsr_int(3)
467  DO ids = 2, nds
468  n = ids - 1
469  dtemp(1) = (-xhi)**n*exp_rab2*sr_int
470  dtemp(2) = real(n, dp)*(-xhi)**(n - 1)*exp_rab2*dsr_int(1)
471  dtemp(3) = real(n**2 - n, dp)/2.0_dp*(-xhi)**(n - 2) &
472  *exp_rab2*dsr_int(2)
473  dtemp(4) = real(n*(n - 1)*(n - 2), dp)/6.0_dp*(-xhi)**(n - 3) &
474  *exp_rab2*dsr_int(3)
475  dtemp(5) = real(n*(n - 1)*(n - 2)*(n - 3), dp)*(-xhi)**(n - 4) &
476  *k*exp_rab2
477  s(ipgfa, jpgfb, m - il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) &
478  + dtemp(4) + dtemp(5)
479  END DO
480  CASE (5)
481  ![s|r^10|s] integral
482  prefac = sqrt_pi3/32.0_dp/sqrt_zet**13
483  temp = 10395.0_dp + 34650.0_dp*pfac*rab2
484  temp = temp + 27720.0_dp*(pfac*rab2)**2 + 7920.0_dp*(pfac*rab2)**3
485  temp = temp + 880.0_dp*(pfac*rab2)**4 + 32.0_dp*(pfac*rab2)**5
486  sr_int = prefac*temp
487  s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int
488  !** derivatives
489  k = sqrt_pi3*b**10/sqrt_zet**23
490  dsr_int(1) = 34650.0_dp*pfac + 55440.0_dp*pfac**2*rab2
491  dsr_int(1) = dsr_int(1) + 23760.0_dp*pfac**3*rab2**2
492  dsr_int(1) = dsr_int(1) + 3520.0_dp*pfac**4*rab2**3
493  dsr_int(1) = dsr_int(1) + 160.0_dp*pfac**5*rab2**4
494  dsr_int(1) = prefac*dsr_int(1)
495  dsr_int(2) = 55440.0_dp*pfac**2 + 47520.0_dp*pfac**3*rab2
496  dsr_int(2) = dsr_int(2) + 10560.0_dp*pfac**4*rab2**2
497  dsr_int(2) = dsr_int(2) + 640.0_dp*pfac**5*rab2**3
498  dsr_int(2) = prefac*dsr_int(2)
499  dsr_int(3) = 47520.0_dp*pfac**3 + 21120.0_dp*pfac**4*rab2
500  dsr_int(3) = dsr_int(3) + 1920.0_dp*pfac**5*rab2**2
501  dsr_int(3) = prefac*dsr_int(3)
502  dsr_int(4) = 21120.0_dp*pfac**4 + 3840.0_dp*pfac**5*rab2
503  dsr_int(4) = prefac*dsr_int(4)
504  DO ids = 2, nds
505  n = ids - 1
506  dtemp(1) = (-xhi)**n*exp_rab2*sr_int
507  dtemp(2) = real(n, dp)*(-xhi)**(n - 1)*exp_rab2*dsr_int(1)
508  dtemp(3) = real(n**2 - n, dp)/2.0_dp*(-xhi)**(n - 2) &
509  *exp_rab2*dsr_int(2)
510  dtemp(4) = real(n*(n - 1)*(n - 2), dp)/6.0_dp*(-xhi)**(n - 3) &
511  *exp_rab2*dsr_int(3)
512  dtemp(5) = real(n*(n - 1)*(n - 2)*(n - 3), dp)/24.0_dp*(-xhi)**(n - 4) &
513  *exp_rab2*dsr_int(4)
514  dtemp(6) = real(n*(n - 1)*(n - 2)*(n - 3)*(n - 4), dp)*(-xhi)**(n - 5) &
515  *k*exp_rab2
516  s(ipgfa, jpgfb, m - il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) &
517  + dtemp(4) + dtemp(5) + dtemp(6)
518  END DO
519  CASE DEFAULT
520  prefac = exp_rab2/sqrt_zet**(2*il + 3)
521  sr_int = 0.0_dp
522  DO i = 0, il
523  sr_int = sr_int + coeff_srs(i + 1, 1, il + 1)*(pfac)**i*rab2**i
524  END DO
525  s(ipgfa, jpgfb, m - il + 1, 1) = prefac*sr_int
526  DO ids = 2, nds
527  n = ids - 1
528  nfac = 1
529  dfsr_int = (-xhi)**n*sr_int
530  DO j = 1, il
531  temp = 0.0_dp
532  DO i = j, il
533  temp = temp + coeff_srs(i + 1, j + 1, il + 1)*(pfac)**i*rab2**(i - j)
534  END DO
535  nfac = nfac*(n - j + 1)
536  dfsr_int = dfsr_int + temp*real(nfac, dp)/fac(j)*(-xhi)**(n - j)
537  END DO
538  s(ipgfa, jpgfb, m - il + 1, ids) = prefac*dfsr_int
539  END DO
540  END SELECT
541  END DO
542  END DO
543  END DO
544 
545  DEALLOCATE (coeff_srs)
546  DEALLOCATE (dtemp, dsr_int)
547 
548  END SUBROUTINE s_ra2m_ab
549 
550 ! **************************************************************************************************
551 !> \brief prefactor for the general formula to calculate the (0a|0b|0b) overlap integrals
552 !> \param nl ...
553 !> \param prefac ...
554 ! **************************************************************************************************
555  SUBROUTINE get_prefac_sabb(nl, prefac)
556  INTEGER, INTENT(IN) :: nl
557  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: prefac
558 
559  INTEGER :: il, j, k
560  REAL(kind=dp) :: sqrt_pi3, temp
561 
562  sqrt_pi3 = sqrt(pi**3)
563 
564  DO il = 0, nl
565  temp = dfac(2*il + 1)*sqrt_pi3*fac(il)/2.0_dp**il
566  DO j = 0, il
567  DO k = j, il
568  prefac(k + 1, j + 1, il + 1) = temp*2.0_dp**k/dfac(2*k + 1)/fac(il - k)/fac(k - j)
569  END DO
570  END DO
571  END DO
572 
573  END SUBROUTINE get_prefac_sabb
574 
575 ! **************************************************************************************************
576 !> \brief calculates the uncontracted, not normalized [s|1/r12|s] two-center coulomb integral
577 !> \param la_max maximal l quantum number on a
578 !> \param npgfa number of primitive Gaussian on a
579 !> \param zeta set of exponents on a
580 !> \param lb_max maximal l quantum number on b
581 !> \param npgfb number of primitive Gaussian on a
582 !> \param zetb set of exponents on a
583 !> \param omega parameter not needed, but given for the sake of the abstract interface
584 !> \param rab distance vector between a and b
585 !> \param v uncontracted coulomb integral of s functions
586 !> \param calculate_forces ...
587 ! **************************************************************************************************
588  SUBROUTINE s_coulomb_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
589 
590  INTEGER, INTENT(IN) :: la_max, npgfa
591  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta
592  INTEGER, INTENT(IN) :: lb_max, npgfb
593  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb
594  REAL(kind=dp), INTENT(IN) :: omega
595  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
596  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
597  LOGICAL, INTENT(IN) :: calculate_forces
598 
599  INTEGER :: ids, ipgfa, jpgfb, n, ndev, nmax
600  REAL(kind=dp) :: a, b, dummy, prefac, rab2, t, xhi, zet
601  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: f
602 
603  dummy = omega
604 
605  ! Distance of the centers a and b
606  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
607  ndev = 0
608  IF (calculate_forces) ndev = 1
609  nmax = la_max + lb_max + ndev + 1
610  ALLOCATE (f(0:nmax))
611  ! Loops over all pairs of primitive Gaussian-type functions
612  DO ipgfa = 1, npgfa
613  DO jpgfb = 1, npgfb
614 
615  ! Calculate some prefactors
616  a = zeta(ipgfa)
617  b = zetb(jpgfb)
618  zet = a + b
619  xhi = a*b/zet
620  prefac = 2.0_dp*sqrt(pi**5)/(a*b)/sqrt(zet)
621  t = xhi*rab2
622  CALL fgamma(nmax - 1, t, f)
623 
624  DO ids = 1, nmax
625  n = ids - 1
626  v(ipgfa, jpgfb, ids) = prefac*(-xhi)**n*f(n)
627  END DO
628 
629  END DO
630  END DO
631  DEALLOCATE (f)
632 
633  END SUBROUTINE s_coulomb_ab
634 
635 ! **************************************************************************************************
636 !> \brief calculates the uncontracted, not normalized [s|1/erf(omega*r12)/r12|s] two-center integral
637 !> \param la_max maximal l quantum number on a
638 !> \param npgfa number of primitive Gaussian on a
639 !> \param zeta set of exponents on a
640 !> \param lb_max maximal l quantum number on b
641 !> \param npgfb number of primitive Gaussian on a
642 !> \param zetb set of exponents on a
643 !> \param omega parameter in the operator
644 !> \param rab distance vector between a and b
645 !> \param v uncontracted erf(r)/r integral of s functions
646 !> \param calculate_forces ...
647 ! **************************************************************************************************
648  SUBROUTINE s_verf_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
649 
650  INTEGER, INTENT(IN) :: la_max, npgfa
651  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta
652  INTEGER, INTENT(IN) :: lb_max, npgfb
653  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb
654  REAL(kind=dp), INTENT(IN) :: omega
655  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
656  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
657  LOGICAL, INTENT(IN) :: calculate_forces
658 
659  INTEGER :: ids, ipgfa, jpgfb, n, ndev, nmax
660  REAL(kind=dp) :: a, arg, b, comega, prefac, rab2, t, xhi, &
661  zet
662  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: f
663 
664  ! Distance of the centers a and b
665  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
666  ndev = 0
667  IF (calculate_forces) ndev = 1
668  nmax = la_max + lb_max + ndev + 1
669  ALLOCATE (f(0:nmax))
670  ! Loops over all pairs of primitive Gaussian-type functions
671  DO ipgfa = 1, npgfa
672  DO jpgfb = 1, npgfb
673 
674  ! Calculate some prefactors
675  a = zeta(ipgfa)
676  b = zetb(jpgfb)
677  zet = a + b
678  xhi = a*b/zet
679  comega = omega**2/(omega**2 + xhi)
680  prefac = 2.0_dp*sqrt(pi**5)*sqrt(comega)/(a*b)/sqrt(zet)
681  t = xhi*rab2
682  arg = comega*t
683  CALL fgamma(nmax - 1, arg, f)
684 
685  DO ids = 1, nmax
686  n = ids - 1
687  v(ipgfa, jpgfb, ids) = prefac*(-xhi*comega)**n*f(n)
688  END DO
689 
690  END DO
691  END DO
692  DEALLOCATE (f)
693 
694  END SUBROUTINE s_verf_ab
695 
696 ! **************************************************************************************************
697 !> \brief calculates the uncontracted, not normalized [s|1/erfc(omega*r12)/r12|s] two-center
698 !> integral
699 !> \param la_max maximal l quantum number on a
700 !> \param npgfa number of primitive Gaussian on a
701 !> \param zeta set of exponents on a
702 !> \param lb_max maximal l quantum number on b
703 !> \param npgfb number of primitive Gaussian on a
704 !> \param zetb set of exponents on a
705 !> \param omega parameter in the operator
706 !> \param rab distance vector between a and b
707 !> \param v uncontracted erf(r)/r integral of s functions
708 !> \param calculate_forces ...
709 ! **************************************************************************************************
710  SUBROUTINE s_verfc_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
711 
712  INTEGER, INTENT(IN) :: la_max, npgfa
713  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta
714  INTEGER, INTENT(IN) :: lb_max, npgfb
715  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb
716  REAL(kind=dp), INTENT(IN) :: omega
717  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
718  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
719  LOGICAL, INTENT(IN) :: calculate_forces
720 
721  INTEGER :: ids, ipgfa, jpgfb, n, ndev, nmax
722  REAL(kind=dp) :: a, b, comega, comegat, prefac, rab2, t, &
723  xhi, zet
724  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: fv, fverf
725 
726  ! Distance of the centers a and b
727  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
728  ndev = 0
729  IF (calculate_forces) ndev = 1
730  nmax = la_max + lb_max + ndev + 1
731  ALLOCATE (fv(0:nmax), fverf(0:nmax))
732  ! Loops over all pairs of primitive Gaussian-type functions
733  DO ipgfa = 1, npgfa
734  DO jpgfb = 1, npgfb
735 
736  ! Calculate some prefactors
737  a = zeta(ipgfa)
738  b = zetb(jpgfb)
739  zet = a + b
740  xhi = a*b/zet
741  comega = omega**2/(omega**2 + xhi)
742  prefac = 2.0_dp*sqrt(pi**5)/(a*b)/sqrt(zet)
743  t = xhi*rab2
744  comegat = comega*t
745  CALL fgamma(nmax - 1, t, fv)
746  CALL fgamma(nmax - 1, comegat, fverf)
747 
748  DO ids = 1, nmax
749  n = ids - 1
750  v(ipgfa, jpgfb, ids) = prefac*(-xhi)**n*(fv(n) - sqrt(comega)*comega**n*fverf(n))
751  END DO
752 
753  END DO
754  END DO
755  DEALLOCATE (fv, fverf)
756 
757  END SUBROUTINE s_verfc_ab
758 
759 ! **************************************************************************************************
760 !> \brief calculates the uncontracted, not normalized [s|exp(-omega*r12**2)/r12|s] two-center
761 !> integral
762 !> \param la_max maximal l quantum number on a
763 !> \param npgfa number of primitive Gaussian on a
764 !> \param zeta set of exponents on a
765 !> \param lb_max maximal l quantum number on b
766 !> \param npgfb number of primitive Gaussian on a
767 !> \param zetb set of exponents on a
768 !> \param omega parameter in the operator
769 !> \param rab distance vector between a and b
770 !> \param v uncontracted exp(-omega*r**2)/r integral of s functions
771 !> \param calculate_forces ...
772 ! **************************************************************************************************
773  SUBROUTINE s_vgauss_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
774 
775  INTEGER, INTENT(IN) :: la_max, npgfa
776  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta
777  INTEGER, INTENT(IN) :: lb_max, npgfb
778  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb
779  REAL(kind=dp), INTENT(IN) :: omega
780  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
781  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
782  LOGICAL, INTENT(IN) :: calculate_forces
783 
784  INTEGER :: ids, ipgfa, j, jpgfb, n, ndev, nmax
785  REAL(kind=dp) :: a, b, eta, etat, expt, oeta, prefac, &
786  rab2, t, xeta, xhi, zet
787  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: f
788 
789  ! Distance of the centers a and b
790  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
791  ndev = 0
792  IF (calculate_forces) ndev = 1
793  nmax = la_max + lb_max + ndev + 1
794  ALLOCATE (f(0:nmax))
795  ! Loops over all pairs of primitive Gaussian-type functions
796  v = 0.0_dp
797  DO ipgfa = 1, npgfa
798  DO jpgfb = 1, npgfb
799 
800  ! Calculate some prefactors
801  a = zeta(ipgfa)
802  b = zetb(jpgfb)
803  zet = a + b
804  xhi = a*b/zet
805  eta = xhi/(xhi + omega)
806  oeta = omega*eta
807  xeta = xhi*eta
808  t = xhi*rab2
809  expt = exp(-omega/(omega + xhi)*t)
810  prefac = 2.0_dp*sqrt(pi**5/zet**3)/(xhi + omega)*expt
811  etat = eta*t
812  CALL fgamma(nmax - 1, etat, f)
813 
814  DO ids = 1, nmax
815  n = ids - 1
816  DO j = 0, n
817  v(ipgfa, jpgfb, ids) = v(ipgfa, jpgfb, ids) &
818  + prefac*fac(n)/fac(j)/fac(n - j)*(-oeta)**(n - j)*(-xeta)**j*f(j)
819  END DO
820  END DO
821 
822  END DO
823  END DO
824  DEALLOCATE (f)
825 
826  END SUBROUTINE s_vgauss_ab
827 
828 ! **************************************************************************************************
829 !> \brief calculates the uncontracted, not normalized [s|exp(-omega*r12**2)|s] two-center
830 !> integral
831 !> \param la_max maximal l quantum number on a
832 !> \param npgfa number of primitive Gaussian on a
833 !> \param zeta set of exponents on a
834 !> \param lb_max maximal l quantum number on b
835 !> \param npgfb number of primitive Gaussian on a
836 !> \param zetb set of exponents on a
837 !> \param omega parameter in the operator
838 !> \param rab distance vector between a and b
839 !> \param v uncontracted exp(-omega*r**2) integral of s functions
840 !> \param calculate_forces ...
841 ! **************************************************************************************************
842  SUBROUTINE s_gauss_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
843 
844  INTEGER, INTENT(IN) :: la_max, npgfa
845  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta
846  INTEGER, INTENT(IN) :: lb_max, npgfb
847  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb
848  REAL(kind=dp), INTENT(IN) :: omega
849  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
850  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
851  LOGICAL, INTENT(IN) :: calculate_forces
852 
853  INTEGER :: ids, ipgfa, jpgfb, n, ndev, nmax
854  REAL(kind=dp) :: a, b, eta, expt, oeta, prefac, rab2, t, &
855  xhi, zet
856  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: f
857 
858  ! Distance of the centers a and b
859  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
860  ndev = 0
861  IF (calculate_forces) ndev = 1
862  nmax = la_max + lb_max + ndev + 1
863  ALLOCATE (f(0:nmax))
864  ! Loops over all pairs of primitive Gaussian-type functions
865  DO ipgfa = 1, npgfa
866  DO jpgfb = 1, npgfb
867 
868  ! Calculate some prefactors
869  a = zeta(ipgfa)
870  b = zetb(jpgfb)
871  zet = a + b
872  xhi = a*b/zet
873  eta = xhi/(xhi + omega)
874  oeta = omega*eta
875  t = xhi*rab2
876  expt = exp(-omega/(omega + xhi)*t)
877  prefac = pi**3/sqrt(zet**3)/sqrt((xhi + omega)**3)*expt
878 
879  DO ids = 1, nmax
880  n = ids - 1
881  v(ipgfa, jpgfb, ids) = prefac*(-oeta)**n
882  END DO
883 
884  END DO
885  END DO
886  DEALLOCATE (f)
887 
888  END SUBROUTINE s_gauss_ab
889 
890 ! **************************************************************************************************
891 !> \brief Contraction and normalization of the [s|O(r12)|s] integrals and their scalar derivatives;
892 !> this routine is more efficient for uncontracted basis sets (clow), e.g. for ri basis sets
893 !> \param la set of l quantum numbers for a
894 !> \param npgfa number of primitive Gaussian on a
895 !> \param nshella number of shells for a
896 !> \param scona_shg SHG contraction/normalization matrix for a
897 !> \param lb set of l quantum numbers for b
898 !> \param npgfb number of primitive Gaussian on b
899 !> \param nshellb number of shells for b
900 !> \param sconb_shg SHG contraction/normalization matrix for b
901 !> \param swork matrix storing the uncontracted and unnormalized s-integrals and derivatives
902 !> \param swork_cont matrix storing the contracted and normalized s-integrals and derivatives
903 !> \param calculate_forces ...
904 ! **************************************************************************************************
905  SUBROUTINE contract_sint_ab_clow(la, npgfa, nshella, scona_shg, lb, npgfb, nshellb, sconb_shg, &
906  swork, swork_cont, calculate_forces)
907 
908  INTEGER, DIMENSION(:), INTENT(IN) :: la
909  INTEGER, INTENT(IN) :: npgfa, nshella
910  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: scona_shg
911  INTEGER, DIMENSION(:), INTENT(IN) :: lb
912  INTEGER, INTENT(IN) :: npgfb, nshellb
913  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: sconb_shg
914  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: swork
915  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: swork_cont
916  LOGICAL, INTENT(IN) :: calculate_forces
917 
918  INTEGER :: ids, ids_start, ipgfa, ishella, jpgfb, &
919  jshellb, lai, lbj, ndev, nds
920 
921  ndev = 0
922  IF (calculate_forces) ndev = 1
923 
924  swork_cont = 0.0_dp
925  DO ishella = 1, nshella
926  lai = la(ishella)
927  DO jshellb = 1, nshellb
928  lbj = lb(jshellb)
929  nds = lai + lbj + 1
930  ids_start = nds - min(lai, lbj)
931  DO ipgfa = 1, npgfa
932  DO jpgfb = 1, npgfb
933  DO ids = ids_start, nds + ndev
934  swork_cont(ids, ishella, jshellb) = swork_cont(ids, ishella, jshellb) &
935  + scona_shg(ipgfa, ishella) &
936  *sconb_shg(jpgfb, jshellb) &
937  *swork(ipgfa, jpgfb, ids)
938  END DO
939  END DO
940  END DO
941  END DO
942  END DO
943 
944  END SUBROUTINE contract_sint_ab_clow
945 
946 ! **************************************************************************************************
947 !> \brief Contraction and normalization of the [s|O(r12)|s] integrals and their scalar derivatives;
948 !> this routine is more efficient for highly contracted basis sets (chigh)
949 !> \param npgfa number of primitive Gaussian on a
950 !> \param nshella number of shells for a
951 !> \param scona SHG contraction/normalization matrix for a
952 !> \param npgfb number of primitive Gaussian on b
953 !> \param nshellb number of shells for b
954 !> \param sconb SHG contraction/normalization matrix for b
955 !> \param nds maximal derivative of [s|O(r12)|s] with respect to rab2
956 !> \param swork matrix storing the uncontracted and unnormalized s-integrals and derivatives
957 !> \param swork_cont matrix storing the contracted and normalized s-integrals and derivatives
958 ! **************************************************************************************************
959  SUBROUTINE contract_sint_ab_chigh(npgfa, nshella, scona, &
960  npgfb, nshellb, sconb, &
961  nds, swork, swork_cont)
962 
963  INTEGER, INTENT(IN) :: npgfa, nshella
964  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: scona
965  INTEGER, INTENT(IN) :: npgfb, nshellb
966  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: sconb
967  INTEGER, INTENT(IN) :: nds
968  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: swork
969  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: swork_cont
970 
971  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: work_pc
972 
973  swork_cont = 0.0_dp
974  ALLOCATE (work_pc(npgfb, nds, nshella))
975 
976  CALL dgemm("T", "N", npgfb*nds, nshella, npgfa, 1._dp, swork, npgfa, scona, npgfa, &
977  0.0_dp, work_pc, npgfb*nds)
978  CALL dgemm("T", "N", nds*nshella, nshellb, npgfb, 1._dp, work_pc, npgfb, sconb, npgfb, &
979  0.0_dp, swork_cont, nds*nshella)
980 
981  DEALLOCATE (work_pc)
982 
983  END SUBROUTINE contract_sint_ab_chigh
984 
985 ! **************************************************************************************************
986 !> \brief Contraction, normalization and combinatiorial combination of the [0a|(r-Ra)^(2m)|0b]
987 !> integrals and their scalar derivatives
988 !> \param npgfa number of primitive Gaussian on a
989 !> \param nshella number of shells for a
990 !> \param scon_ra2m contraction matrix on a containg the combinatorial factors
991 !> \param npgfb number of primitive Gaussian on b
992 !> \param nshellb number of shells for b
993 !> \param sconb SHG contraction/normalization matrix for b
994 !> \param swork matrix storing the uncontracted and unnormalized s-integrals and derivatives
995 !> \param swork_cont matrix storing the contracted and normalized s-integrals and derivatives
996 !> \param m exponent in operator (r-Ra)^(2m)
997 !> \param nds_max maximal derivative with respect to rab2
998 ! **************************************************************************************************
999  SUBROUTINE contract_s_ra2m_ab(npgfa, nshella, scon_ra2m, npgfb, nshellb, sconb, swork, swork_cont, &
1000  m, nds_max)
1001 
1002  INTEGER, INTENT(IN) :: npgfa, nshella
1003  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: scon_ra2m
1004  INTEGER, INTENT(IN) :: npgfb, nshellb
1005  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: sconb
1006  REAL(kind=dp), DIMENSION(:, :, :, :), &
1007  INTENT(INOUT) :: swork
1008  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: swork_cont
1009  INTEGER, INTENT(IN) :: m, nds_max
1010 
1011  INTEGER :: i, my_m
1012  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: work_pc
1013 
1014  my_m = m + 1
1015  ALLOCATE (work_pc(npgfb, nds_max, nshella))
1016  work_pc = 0.0_dp
1017  swork_cont = 0.0_dp
1018  DO i = 1, my_m
1019  CALL dgemm("T", "N", npgfb*nds_max, nshella, npgfa, 1.0_dp, swork(1:npgfa, 1:npgfb, i, 1:nds_max), npgfa, &
1020  scon_ra2m(1:npgfa, i, 1:nshella), npgfa, 1.0_dp, work_pc, npgfb*nds_max)
1021  END DO
1022  CALL dgemm("T", "N", nds_max*nshella, nshellb, npgfb, 1.0_dp, work_pc, npgfb, sconb, npgfb, 0.0_dp, &
1023  swork_cont, nds_max*nshella)
1024 
1025  DEALLOCATE (work_pc)
1026 
1027  END SUBROUTINE contract_s_ra2m_ab
1028 
1029 ! **************************************************************************************************
1030 !> \brief Contraction and normalization of the (abb) overlap
1031 !> \param la set of l quantum numbers on a
1032 !> \param npgfa number of primitive Gaussians functions on a; orbital basis
1033 !> \param nshella number of shells for a; orbital basis
1034 !> \param scona_shg SHG contraction/normalization matrix for a; orbital basis
1035 !> \param lb set of l quantum numbers on b; orbital basis
1036 !> \param npgfb number of primitive Gaussians functions on b; orbital basis
1037 !> \param nshellb number of shells for b; orbital basis
1038 !> \param lcb set of l quantum numbers on b; ri basis
1039 !> \param npgfcb number of primitive Gaussians functions on b; ri basis
1040 !> \param nshellcb number of shells for b; ri basis
1041 !> \param orbb_index index for orbital basis at B for sconb_mix
1042 !> \param rib_index index for ri basis at B for sconb_mix
1043 !> \param sconb_mix mixed contraction matrix for orbital and ri basis at B
1044 !> \param nl_max related to the parameter m in (a|rb^(2m)|b)
1045 !> \param nds_max derivative with respect to rab**2
1046 !> \param swork matrix of storing the uncontracted and unnormalized s-integrals and derivatives
1047 !> \param swork_cont matrix storing the contracted and normalized s-integrals and derivatives
1048 !> \param calculate_forces ...
1049 ! **************************************************************************************************
1050  SUBROUTINE contract_s_overlap_abb(la, npgfa, nshella, scona_shg, lb, npgfb, nshellb, &
1051  lcb, npgfcb, nshellcb, orbb_index, rib_index, sconb_mix, &
1052  nl_max, nds_max, swork, swork_cont, calculate_forces)
1053 
1054  INTEGER, DIMENSION(:), INTENT(IN) :: la
1055  INTEGER, INTENT(IN) :: npgfa, nshella
1056  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: scona_shg
1057  INTEGER, DIMENSION(:), INTENT(IN) :: lb
1058  INTEGER, INTENT(IN) :: npgfb, nshellb
1059  INTEGER, DIMENSION(:), INTENT(IN) :: lcb
1060  INTEGER, INTENT(IN) :: npgfcb, nshellcb
1061  INTEGER, DIMENSION(:, :), INTENT(IN) :: orbb_index, rib_index
1062  REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: sconb_mix
1063  INTEGER, INTENT(IN) :: nl_max, nds_max
1064  REAL(kind=dp), DIMENSION(:, :, :, :, :), &
1065  INTENT(IN) :: swork
1066  REAL(kind=dp), DIMENSION(:, 0:, :, :, :), &
1067  INTENT(INOUT) :: swork_cont
1068  LOGICAL, INTENT(IN) :: calculate_forces
1069 
1070  INTEGER :: forb, fri, ids, ids_start, iil, il, &
1071  ishella, jpgfb, jshellb, kpgfb, &
1072  kshellb, lai, lbb, lbj, lbk, ndev, &
1073  nds, nl
1074  REAL(kind=dp), ALLOCATABLE, &
1075  DIMENSION(:, :, :, :, :) :: work_ppc
1076 
1077  ndev = 0
1078  IF (calculate_forces) ndev = 1
1079 
1080  ALLOCATE (work_ppc(npgfb, npgfcb, nl_max, nds_max, nshella))
1081  work_ppc = 0.0_dp
1082  CALL dgemm("T", "N", npgfb*npgfcb*nl_max*nds_max, nshella, npgfa, 1._dp, swork, npgfa, &
1083  scona_shg, npgfa, 0.0_dp, work_ppc, npgfb*npgfcb*nl_max*nds_max)
1084  swork_cont = 0.0_dp
1085  DO kshellb = 1, nshellcb
1086  lbk = lcb(kshellb)
1087  DO jshellb = 1, nshellb
1088  lbj = lb(jshellb)
1089  nl = int((lbj + lbk)/2)
1090  IF (lbj == 0 .OR. lbk == 0) nl = 0
1091  DO ishella = 1, nshella
1092  lai = la(ishella)
1093  DO il = 0, nl
1094  lbb = lbj + lbk - 2*il
1095  nds = lai + lbb + 1
1096  ids_start = nds - min(lai, lbb)
1097  DO jpgfb = 1, npgfb
1098  forb = orbb_index(jpgfb, jshellb)
1099  DO kpgfb = 1, npgfcb
1100  fri = rib_index(kpgfb, kshellb)
1101  DO ids = ids_start, nds + ndev
1102  DO iil = 0, il
1103  swork_cont(ids, il, ishella, jshellb, kshellb) = &
1104  swork_cont(ids, il, ishella, jshellb, kshellb) &
1105  + sconb_mix(iil + 1, fri, forb, il + 1)*work_ppc(jpgfb, kpgfb, il - iil + 1, ids, ishella)
1106  END DO
1107  END DO
1108  END DO
1109  END DO
1110  END DO
1111  END DO
1112  END DO
1113  END DO
1114  DEALLOCATE (work_ppc)
1115 
1116  END SUBROUTINE contract_s_overlap_abb
1117 
1118 ! **************************************************************************************************
1119 !> \brief Contraction and normalization of the (aba) overlap
1120 !> \param la set of l quantum numbers on a; orbital basis
1121 !> \param npgfa number of primitive Gaussians functions on a; orbital basis
1122 !> \param nshella number of shells for a; orbital basis
1123 !> \param lb set of l quantum numbers on b; orbital basis
1124 !> \param npgfb number of primitive Gaussians functions on b; orbital basis
1125 !> \param nshellb number of shells for b; orbital basis
1126 !> \param sconb_shg SHG contraction/normalization matrix for b; orbital basis
1127 !> \param lca set of l quantum numbers on a; ri basis
1128 !> \param npgfca number of primitive Gaussians functions on a; ri basis
1129 !> \param nshellca number of shells for a; ri basis
1130 !> \param orba_index index for orbital basis at A for scona_mix
1131 !> \param ria_index index for ri basis at A for scona_mix
1132 !> \param scona_mix mixed contraction matrix for orbital and ri basis at A
1133 !> \param nl_max related to the parameter m in (a|ra^(2m)|b)
1134 !> \param nds_max maximal derivative with respect to rab**2
1135 !> \param swork matrix of storing the uncontracted and unnormalized s-integrals and derivatives
1136 !> \param swork_cont matrix storing the contracted and normalized s-integrals and derivatives
1137 !> \param calculate_forces ...
1138 ! **************************************************************************************************
1139  SUBROUTINE contract_s_overlap_aba(la, npgfa, nshella, lb, npgfb, nshellb, sconb_shg, &
1140  lca, npgfca, nshellca, orba_index, ria_index, scona_mix, &
1141  nl_max, nds_max, swork, swork_cont, calculate_forces)
1142 
1143  INTEGER, DIMENSION(:), INTENT(IN) :: la
1144  INTEGER, INTENT(IN) :: npgfa, nshella
1145  INTEGER, DIMENSION(:), INTENT(IN) :: lb
1146  INTEGER, INTENT(IN) :: npgfb, nshellb
1147  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: sconb_shg
1148  INTEGER, DIMENSION(:), INTENT(IN) :: lca
1149  INTEGER, INTENT(IN) :: npgfca, nshellca
1150  INTEGER, DIMENSION(:, :), INTENT(IN) :: orba_index, ria_index
1151  REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: scona_mix
1152  INTEGER, INTENT(IN) :: nl_max, nds_max
1153  REAL(kind=dp), DIMENSION(:, :, :, :, :), &
1154  INTENT(IN) :: swork
1155  REAL(kind=dp), DIMENSION(:, 0:, :, :, :), &
1156  INTENT(INOUT) :: swork_cont
1157  LOGICAL, INTENT(IN) :: calculate_forces
1158 
1159  INTEGER :: forb, fri, ids, ids_start, iil, il, &
1160  ipgfa, ishella, jshellb, kpgfa, &
1161  kshella, laa, lai, lak, lbj, ndev, &
1162  nds, nl
1163  REAL(kind=dp), ALLOCATABLE, &
1164  DIMENSION(:, :, :, :, :) :: work_ppc
1165 
1166  ndev = 0
1167  IF (calculate_forces) ndev = 1
1168 
1169  ALLOCATE (work_ppc(npgfa, npgfca, nl_max, nds_max, nshellb))
1170  work_ppc = 0.0_dp
1171  CALL dgemm("T", "N", npgfa*npgfca*nl_max*nds_max, nshellb, npgfb, 1._dp, swork, npgfb, &
1172  sconb_shg, npgfb, 0.0_dp, work_ppc, npgfa*npgfca*nl_max*nds_max)
1173  swork_cont = 0.0_dp
1174  DO kshella = 1, nshellca
1175  lak = lca(kshella)
1176  DO jshellb = 1, nshellb
1177  lbj = lb(jshellb)
1178  DO ishella = 1, nshella
1179  lai = la(ishella)
1180  nl = int((lai + lak)/2)
1181  IF (lai == 0 .OR. lak == 0) nl = 0
1182  DO il = 0, nl
1183  laa = lai + lak - 2*il
1184  nds = laa + lbj + 1
1185  ids_start = nds - min(laa, lbj)
1186  DO kpgfa = 1, npgfca
1187  fri = ria_index(kpgfa, kshella)
1188  DO ipgfa = 1, npgfa
1189  forb = orba_index(ipgfa, ishella)
1190  DO ids = ids_start, nds + ndev
1191  DO iil = 0, il
1192  swork_cont(ids, il, ishella, jshellb, kshella) = &
1193  swork_cont(ids, il, ishella, jshellb, kshella) &
1194  + scona_mix(iil + 1, fri, forb, il + 1)*work_ppc(ipgfa, kpgfa, il - iil + 1, ids, jshellb)
1195  END DO
1196  END DO
1197  END DO
1198  END DO
1199  END DO
1200  END DO
1201  END DO
1202  END DO
1203 
1204  DEALLOCATE (work_ppc)
1205  END SUBROUTINE contract_s_overlap_aba
1206 
1207 END MODULE s_contract_shg
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition: gamma.F:15
subroutine, public fgamma_0(nmax, t, f)
Calculation of the incomplete Gamma function F(t) for multicenter integrals over Gaussian functions....
Definition: gamma.F:154
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public pi
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
Definition: mathconstants.F:61
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition: mathconstants.F:37
Routines for calculating the s-integrals and their scalar derivatives with respect to rab2 over solid...
subroutine, public 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)
Contraction and normalization of the (aba) overlap.
subroutine, public s_verf_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
calculates the uncontracted, not normalized [s|1/erf(omega*r12)/r12|s] two-center integral
subroutine, public s_gauss_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
calculates the uncontracted, not normalized [s|exp(-omega*r12**2)|s] two-center integral
subroutine, public s_overlap_abb(la_max, npgfa, zeta, lb_max, npgfb, zetb, lcb_max, npgfcb, zetcb, rab, s, calculate_forces)
calculates [s|ra^n|s] integrals for [aba] and the [s|rb^n|s] integrals for [abb]
subroutine, public s_vgauss_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
calculates the uncontracted, not normalized [s|exp(-omega*r12**2)/r12|s] two-center integral
subroutine, public contract_s_ra2m_ab(npgfa, nshella, scon_ra2m, npgfb, nshellb, sconb, swork, swork_cont, m, nds_max)
Contraction, normalization and combinatiorial combination of the [0a|(r-Ra)^(2m)|0b] integrals and th...
subroutine, public contract_sint_ab_chigh(npgfa, nshella, scona, npgfb, nshellb, sconb, nds, swork, swork_cont)
Contraction and normalization of the [s|O(r12)|s] integrals and their scalar derivatives; this routin...
subroutine, public s_coulomb_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
calculates the uncontracted, not normalized [s|1/r12|s] two-center coulomb integral
subroutine, public 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)
Contraction and normalization of the (abb) overlap.
subroutine, public s_verfc_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
calculates the uncontracted, not normalized [s|1/erfc(omega*r12)/r12|s] two-center integral
subroutine, public contract_sint_ab_clow(la, npgfa, nshella, scona_shg, lb, npgfb, nshellb, sconb_shg, swork, swork_cont, calculate_forces)
Contraction and normalization of the [s|O(r12)|s] integrals and their scalar derivatives; this routin...
subroutine, public s_overlap_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, rab, s, calculate_forces)
calculates the uncontracted, not normalized [s|s] overlap
subroutine, public s_ra2m_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, m, rab, s, calculate_forces)
calculates the uncontracted, not normalized [0a|ra^(2m)|0b] two-center integral, where ra = r-Ra and ...