Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of contracted, spherical Gaussian integrals using the (OS) integral
10 : !> scheme. Routines for the following two-center integrals:
11 : !> i) (a|O(r12)|b) where O(r12) is the overlap, coulomb operator etc.
12 : !> ii) (aba) and (abb) s-overlaps
13 : !> \par History
14 : !> created [06.2015]
15 : !> 05.2019: Added truncated coulomb operator (A. Bussy)
16 : !> \author Dorothea Golze
17 : ! **************************************************************************************************
18 : MODULE generic_os_integrals
19 : USE ai_cartesian_operators, ONLY: apply_cartesian_operator
20 : USE ai_contraction_sphi, ONLY: ab_contract,&
21 : abc_contract,&
22 : abcd_contract
23 : USE ai_operator_ra2m, ONLY: operator_ra2m
24 : USE ai_operators_r12, ONLY: ab_sint_os,&
25 : cps_coulomb2,&
26 : cps_gauss2,&
27 : cps_truncated2,&
28 : cps_verf2,&
29 : cps_verfc2,&
30 : cps_vgauss2,&
31 : operator2
32 : USE ai_overlap, ONLY: overlap_aab,&
33 : overlap_ab,&
34 : overlap_abb
35 : USE ai_overlap_aabb, ONLY: overlap_aabb
36 : USE basis_set_types, ONLY: get_gto_basis_set,&
37 : gto_basis_set_type
38 : USE constants_operator, ONLY: operator_coulomb,&
39 : operator_gauss,&
40 : operator_truncated,&
41 : operator_verf,&
42 : operator_verfc,&
43 : operator_vgauss
44 : USE debug_os_integrals, ONLY: overlap_aabb_test,&
45 : overlap_ab_test,&
46 : overlap_abc_test
47 : USE kinds, ONLY: dp,&
48 : int_8
49 : USE orbital_pointers, ONLY: ncoset
50 : #include "./base/base_uses.f90"
51 :
52 : IMPLICIT NONE
53 :
54 : PRIVATE
55 :
56 : ! **************************************************************************************************
57 :
58 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'generic_os_integrals'
59 :
60 : PUBLIC :: int_operators_r12_ab_os, int_overlap_ab_os, int_ra2m_ab_os, int_overlap_aba_os, &
61 : int_overlap_abb_os, int_overlap_aabb_os
62 :
63 : CONTAINS
64 :
65 : ! **************************************************************************************************
66 : !> \brief Calcululates the two-center integrals of the type (a|O(r12)|b) using the OS scheme
67 : !> \param r12_operator the integral operator, which depends on r12=|r1-r2|
68 : !> \param vab integral matrix of spherical contracted Gaussian functions
69 : !> \param dvab derivative of the integrals
70 : !> \param rab distance vector between center A and B
71 : !> \param fba basis at center A
72 : !> \param fbb basis at center B
73 : !> \param omega parameter in the operator
74 : !> \param r_cutoff the cutoff in case of truncated coulomb operator
75 : !> \param calculate_forces ...
76 : ! **************************************************************************************************
77 284 : SUBROUTINE int_operators_r12_ab_os(r12_operator, vab, dvab, rab, fba, fbb, omega, r_cutoff, &
78 : calculate_forces)
79 :
80 : INTEGER, INTENT(IN) :: r12_operator
81 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
82 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
83 : OPTIONAL :: dvab
84 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
85 : TYPE(gto_basis_set_type), POINTER :: fba, fbb
86 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: omega, r_cutoff
87 : LOGICAL, INTENT(IN) :: calculate_forces
88 :
89 : CHARACTER(LEN=*), PARAMETER :: routineN = 'int_operators_r12_ab_os'
90 :
91 : INTEGER :: handle
92 : REAL(KIND=dp) :: my_omega, my_r_cutoff
93 :
94 : PROCEDURE(ab_sint_os), POINTER :: cps_operator2
95 :
96 284 : NULLIFY (cps_operator2)
97 284 : CALL timeset(routineN, handle)
98 :
99 284 : my_omega = 1.0_dp
100 284 : my_r_cutoff = 1.0_dp
101 :
102 401 : SELECT CASE (r12_operator)
103 : CASE (operator_coulomb)
104 117 : cps_operator2 => cps_coulomb2
105 : CASE (operator_verf)
106 32 : cps_operator2 => cps_verf2
107 32 : IF (PRESENT(omega)) my_omega = omega
108 : CASE (operator_verfc)
109 32 : cps_operator2 => cps_verfc2
110 32 : IF (PRESENT(omega)) my_omega = omega
111 : CASE (operator_vgauss)
112 32 : cps_operator2 => cps_vgauss2
113 32 : IF (PRESENT(omega)) my_omega = omega
114 : CASE (operator_gauss)
115 32 : cps_operator2 => cps_gauss2
116 32 : IF (PRESENT(omega)) my_omega = omega
117 : CASE (operator_truncated)
118 39 : cps_operator2 => cps_truncated2
119 39 : IF (PRESENT(r_cutoff)) my_r_cutoff = r_cutoff
120 : CASE DEFAULT
121 284 : CPABORT("Operator not available")
122 : END SELECT
123 :
124 : CALL int_operator_ab_os_low(cps_operator2, vab, dvab, rab, fba, fbb, my_omega, my_r_cutoff, &
125 408 : calculate_forces)
126 :
127 284 : CALL timestop(handle)
128 :
129 284 : END SUBROUTINE int_operators_r12_ab_os
130 :
131 : ! **************************************************************************************************
132 : !> \brief calculate integrals (a|O(r12)|b)
133 : !> \param cps_operator2 procedure pointer for the respective operator.
134 : !> \param vab integral matrix of spherical contracted Gaussian functions
135 : !> \param dvab derivative of the integrals
136 : !> \param rab distance vector between center A and B
137 : !> \param fba basis at center A
138 : !> \param fbb basis at center B
139 : !> \param omega parameter in the operator
140 : !> \param r_cutoff ...
141 : !> \param calculate_forces ...
142 : ! **************************************************************************************************
143 284 : SUBROUTINE int_operator_ab_os_low(cps_operator2, vab, dvab, rab, fba, fbb, omega, r_cutoff, &
144 : calculate_forces)
145 :
146 : PROCEDURE(ab_sint_os), POINTER :: cps_operator2
147 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
148 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
149 : INTENT(INOUT) :: dvab
150 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
151 : TYPE(gto_basis_set_type), POINTER :: fba, fbb
152 : REAL(KIND=dp), INTENT(IN) :: omega, r_cutoff
153 : LOGICAL, INTENT(IN) :: calculate_forces
154 :
155 : CHARACTER(LEN=*), PARAMETER :: routineN = 'int_operator_ab_os_low'
156 :
157 : INTEGER :: handle, i, iset, jset, lds, m1, m2, &
158 : maxco, maxcoa, maxcob, maxl, maxla, &
159 : maxlb, ncoa, ncoap, ncob, ncobp, &
160 : nseta, nsetb, sgfa, sgfb
161 284 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
162 284 : npgfb, nsgfa, nsgfb
163 284 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
164 : REAL(KIND=dp) :: dab, rab2
165 284 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vac, vac_plus
166 284 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: devab, vwork
167 284 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: sphi_a, sphi_b, zeta, zetb
168 :
169 284 : CALL timeset(routineN, handle)
170 284 : NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, &
171 284 : first_sgfa, first_sgfb, sphi_a, sphi_b, zeta, zetb)
172 :
173 : ! basis ikind
174 284 : first_sgfa => fba%first_sgf
175 284 : la_max => fba%lmax
176 284 : la_min => fba%lmin
177 284 : npgfa => fba%npgf
178 284 : nseta = fba%nset
179 284 : nsgfa => fba%nsgf_set
180 284 : sphi_a => fba%sphi
181 284 : zeta => fba%zet
182 : ! basis jkind
183 284 : first_sgfb => fbb%first_sgf
184 284 : lb_max => fbb%lmax
185 284 : lb_min => fbb%lmin
186 284 : npgfb => fbb%npgf
187 284 : nsetb = fbb%nset
188 284 : nsgfb => fbb%nsgf_set
189 284 : sphi_b => fbb%sphi
190 284 : zetb => fbb%zet
191 :
192 284 : CALL get_gto_basis_set(fba, maxco=maxcoa, maxl=maxla)
193 284 : CALL get_gto_basis_set(fbb, maxco=maxcob, maxl=maxlb)
194 284 : maxco = MAX(maxcoa, maxcob)
195 284 : IF (calculate_forces) THEN
196 : maxl = MAX(maxla + 1, maxlb)
197 : ELSE
198 284 : maxl = MAX(maxla, maxlb)
199 : END IF
200 284 : lds = ncoset(maxl)
201 :
202 1136 : rab2 = SUM(rab*rab)
203 : dab = SQRT(rab2)
204 :
205 5943564 : vab = 0.0_dp
206 14458524 : IF (calculate_forces) dvab = 0.0_dp
207 :
208 3274 : DO iset = 1, nseta
209 :
210 2990 : ncoa = npgfa(iset)*ncoset(la_max(iset))
211 2990 : ncoap = npgfa(iset)*ncoset(la_max(iset) + 1)
212 2990 : sgfa = first_sgfa(1, iset)
213 :
214 58118 : DO jset = 1, nsetb
215 :
216 54844 : ncob = npgfb(jset)*ncoset(lb_max(jset))
217 54844 : ncobp = npgfb(jset)*ncoset(lb_max(jset) + 1)
218 54844 : sgfb = first_sgfb(1, jset)
219 54844 : m1 = sgfa + nsgfa(iset) - 1
220 54844 : m2 = sgfb + nsgfb(jset) - 1
221 :
222 : ! calculate integrals
223 54844 : IF (calculate_forces) THEN
224 0 : ALLOCATE (vwork(ncoap, ncobp, la_max(iset) + lb_max(jset) + 3), &
225 360000 : vac(ncoa, ncob), vac_plus(ncoap, ncobp), devab(ncoa, ncob, 3))
226 24000 : devab = 0._dp
227 24000 : vwork = 0.0_dp
228 24000 : vac = 0.0_dp
229 : CALL operator2(cps_operator2, la_max(iset) + 1, npgfa(iset), zeta(:, iset), la_min(iset), &
230 : lb_max(jset) + 1, npgfb(jset), zetb(:, jset), lb_min(jset), &
231 24000 : omega, r_cutoff, rab, rab2, vac, vwork, maxder=1, vac_plus=vac_plus)
232 : CALL apply_cartesian_operator(la_max(iset), 0, npgfa(iset), lb_max(jset), 0, npgfb(jset), &
233 24000 : vac_plus, devab, right=.FALSE., derivative=.TRUE., exponents=zeta(:, iset))
234 96000 : DO i = 1, 3
235 : CALL ab_contract(dvab(sgfa:m1, sgfb:m2, i), devab(:, :, i), sphi_a(:, sgfa:), &
236 96000 : sphi_b(:, sgfb:), ncoa, ncob, nsgfa(iset), nsgfb(jset))
237 : END DO
238 :
239 : ELSE
240 0 : ALLOCATE (vwork(ncoa, ncob, la_max(iset) + lb_max(jset) + 1), &
241 462660 : vac(ncoa, ncob), vac_plus(ncoap, ncobp), devab(ncoa, ncob, 3))
242 30844 : vwork = 0.0_dp
243 30844 : vac = 0.0_dp
244 : CALL operator2(cps_operator2, la_max(iset), npgfa(iset), zeta(:, iset), la_min(iset), &
245 : lb_max(jset), npgfb(jset), zetb(:, jset), lb_min(jset), &
246 30844 : omega, r_cutoff, rab, rab2, vac, vwork)
247 : END IF
248 :
249 : CALL ab_contract(vab(sgfa:m1, sgfb:m2), vac(1:ncoa, 1:ncob), sphi_a(:, sgfa:), sphi_b(:, sgfb:), &
250 54844 : ncoa, ncob, nsgfa(iset), nsgfb(jset))
251 57834 : DEALLOCATE (vwork, vac, vac_plus, devab)
252 : END DO
253 : END DO
254 :
255 284 : CALL timestop(handle)
256 :
257 568 : END SUBROUTINE int_operator_ab_os_low
258 :
259 : ! **************************************************************************************************
260 : !> \brief calculate overlap integrals (a,b)
261 : !> \param sab integral (a,b)
262 : !> \param dsab derivative of sab with respect to A
263 : !> \param rab distance vector between center A and B
264 : !> \param fba basis at center A
265 : !> \param fbb basis at center B
266 : !> \param calculate_forces ...
267 : !> \param debug integrals are debugged by recursive routines if requested
268 : !> \param dmax maximal deviation between integrals when debugging
269 : ! **************************************************************************************************
270 2296 : SUBROUTINE int_overlap_ab_os(sab, dsab, rab, fba, fbb, calculate_forces, debug, dmax)
271 :
272 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: sab
273 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
274 : OPTIONAL :: dsab
275 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
276 : TYPE(gto_basis_set_type), POINTER :: fba, fbb
277 : LOGICAL, INTENT(IN) :: calculate_forces, debug
278 : REAL(KIND=dp), INTENT(INOUT) :: dmax
279 :
280 : CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_ab_os'
281 :
282 : INTEGER :: handle, i, iset, jset, m1, m2, maxco, &
283 : maxcoa, maxcob, maxl, maxla, maxlb, &
284 : ncoa, ncob, nseta, nsetb, sgfa, sgfb
285 2296 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
286 2296 : npgfb, nsgfa, nsgfb
287 2296 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
288 : LOGICAL :: failure
289 : REAL(KIND=dp) :: dab, ra(3), rb(3)
290 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sint
291 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dsint
292 2296 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
293 2296 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, scon_a, scon_b, zeta, zetb
294 :
295 2296 : failure = .FALSE.
296 2296 : CALL timeset(routineN, handle)
297 2296 : NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, &
298 2296 : first_sgfa, first_sgfb, scon_a, scon_b, zeta, zetb)
299 :
300 : ! basis ikind
301 2296 : first_sgfa => fba%first_sgf
302 2296 : la_max => fba%lmax
303 2296 : la_min => fba%lmin
304 2296 : npgfa => fba%npgf
305 2296 : nseta = fba%nset
306 2296 : nsgfa => fba%nsgf_set
307 2296 : rpgfa => fba%pgf_radius
308 2296 : set_radius_a => fba%set_radius
309 2296 : scon_a => fba%scon
310 2296 : zeta => fba%zet
311 : ! basis jkind
312 2296 : first_sgfb => fbb%first_sgf
313 2296 : lb_max => fbb%lmax
314 2296 : lb_min => fbb%lmin
315 2296 : npgfb => fbb%npgf
316 2296 : nsetb = fbb%nset
317 2296 : nsgfb => fbb%nsgf_set
318 2296 : rpgfb => fbb%pgf_radius
319 2296 : set_radius_b => fbb%set_radius
320 2296 : scon_b => fbb%scon
321 2296 : zetb => fbb%zet
322 :
323 2296 : CALL get_gto_basis_set(fba, maxco=maxcoa, maxl=maxla)
324 2296 : CALL get_gto_basis_set(fbb, maxco=maxcob, maxl=maxlb)
325 2296 : maxco = MAX(maxcoa, maxcob)
326 2296 : maxl = MAX(maxla, maxlb)
327 9184 : ALLOCATE (sint(maxco, maxco))
328 11480 : ALLOCATE (dsint(maxco, maxco, 3))
329 :
330 9184 : dab = SQRT(SUM(rab**2))
331 :
332 4523933 : sab = 0.0_dp
333 2296 : IF (calculate_forces) THEN
334 6857394 : IF (PRESENT(dsab)) dsab = 0.0_dp
335 : END IF
336 :
337 11826 : DO iset = 1, nseta
338 :
339 9530 : ncoa = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
340 9530 : sgfa = first_sgfa(1, iset)
341 :
342 75365 : DO jset = 1, nsetb
343 :
344 63539 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
345 :
346 22981 : ncob = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
347 22981 : sgfb = first_sgfb(1, jset)
348 22981 : m1 = sgfa + nsgfa(iset) - 1
349 22981 : m2 = sgfb + nsgfb(jset) - 1
350 :
351 22981 : IF (calculate_forces) THEN
352 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
353 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
354 11517 : rab, sint, dsint)
355 : ELSE
356 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
357 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
358 11464 : rab, sint)
359 :
360 : END IF
361 :
362 : ! debug if requested
363 22981 : IF (debug) THEN
364 47 : ra = 0.0_dp
365 47 : rb = rab
366 : CALL overlap_ab_test(la_max(iset), la_min(iset), npgfa(iset), zeta(:, iset), &
367 : lb_max(jset), lb_min(jset), npgfb(jset), zetb(:, jset), &
368 47 : ra, rb, sint, dmax)
369 : END IF
370 :
371 : CALL ab_contract(sab(sgfa:m1, sgfb:m2), sint(1:ncoa, 1:ncob), scon_a(:, sgfa:), scon_b(:, sgfb:), &
372 22981 : ncoa, ncob, nsgfa(iset), nsgfb(jset))
373 32511 : IF (calculate_forces) THEN
374 46068 : DO i = 1, 3
375 : CALL ab_contract(dsab(sgfa:m1, sgfb:m2, i), dsint(1:ncoa, 1:ncob, i), scon_a(:, sgfa:), &
376 46068 : scon_b(:, sgfb:), ncoa, ncob, nsgfa(iset), nsgfb(jset))
377 : END DO
378 : END IF
379 : END DO
380 : END DO
381 :
382 2296 : DEALLOCATE (sint, dsint)
383 :
384 2296 : CALL timestop(handle)
385 :
386 2296 : END SUBROUTINE int_overlap_ab_os
387 :
388 : ! **************************************************************************************************
389 : !> \brief calculate integrals (a|(r-Ra)^(2m)|b)
390 : !> \param sab integral (a|(r-Ra)^(2m)|b)
391 : !> \param dsab derivative of sab with respect to A
392 : !> \param rab distance vector between center A and B
393 : !> \param fba fba basis at center A
394 : !> \param fbb fbb basis at center B
395 : !> \param m exponent m in operator (r-Ra)^(2m)
396 : !> \param calculate_forces ...
397 : ! **************************************************************************************************
398 32 : SUBROUTINE int_ra2m_ab_os(sab, dsab, rab, fba, fbb, m, calculate_forces)
399 :
400 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: sab
401 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
402 : OPTIONAL :: dsab
403 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
404 : TYPE(gto_basis_set_type), POINTER :: fba, fbb
405 : INTEGER, INTENT(IN) :: m
406 : LOGICAL, INTENT(IN) :: calculate_forces
407 :
408 : CHARACTER(LEN=*), PARAMETER :: routineN = 'int_ra2m_ab_os'
409 :
410 : INTEGER :: handle, i, iset, jset, m1, m2, maxco, &
411 : maxcoa, maxcob, maxl, maxla, maxlb, &
412 : ncoa, ncob, nseta, nsetb, sgfa, sgfb
413 32 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
414 32 : npgfb, nsgfa, nsgfb
415 32 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
416 : LOGICAL :: failure
417 : REAL(KIND=dp) :: dab
418 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sint
419 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dsint
420 32 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: scon_a, scon_b, zeta, zetb
421 :
422 32 : failure = .FALSE.
423 32 : CALL timeset(routineN, handle)
424 32 : NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, &
425 32 : first_sgfa, first_sgfb, scon_a, scon_b, zeta, zetb)
426 :
427 : ! basis ikind
428 32 : first_sgfa => fba%first_sgf
429 32 : la_max => fba%lmax
430 32 : la_min => fba%lmin
431 32 : npgfa => fba%npgf
432 32 : nseta = fba%nset
433 32 : nsgfa => fba%nsgf_set
434 32 : scon_a => fba%scon
435 32 : zeta => fba%zet
436 : ! basis jkind
437 32 : first_sgfb => fbb%first_sgf
438 32 : lb_max => fbb%lmax
439 32 : lb_min => fbb%lmin
440 32 : npgfb => fbb%npgf
441 32 : nsetb = fbb%nset
442 32 : nsgfb => fbb%nsgf_set
443 32 : scon_b => fbb%scon
444 32 : zetb => fbb%zet
445 :
446 32 : CALL get_gto_basis_set(fba, maxco=maxcoa, maxl=maxla)
447 32 : CALL get_gto_basis_set(fbb, maxco=maxcob, maxl=maxlb)
448 32 : maxco = MAX(maxcoa, maxcob)
449 32 : maxl = MAX(maxla, maxlb)
450 128 : ALLOCATE (sint(maxco, maxco))
451 160 : ALLOCATE (dsint(maxco, maxco, 3))
452 :
453 128 : dab = SQRT(SUM(rab**2))
454 :
455 963872 : sab = 0.0_dp
456 32 : IF (calculate_forces) THEN
457 2891680 : IF (PRESENT(dsab)) dsab = 0.0_dp
458 : END IF
459 :
460 352 : DO iset = 1, nseta
461 :
462 320 : ncoa = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
463 320 : sgfa = first_sgfa(1, iset)
464 :
465 5152 : DO jset = 1, nsetb
466 :
467 4800 : ncob = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
468 4800 : sgfb = first_sgfb(1, jset)
469 4800 : m1 = sgfa + nsgfa(iset) - 1
470 4800 : m2 = sgfb + nsgfb(jset) - 1
471 :
472 : CALL operator_ra2m(la_max(iset), la_min(iset), npgfa(iset), zeta(:, iset), &
473 : lb_max(jset), lb_min(jset), npgfb(jset), zetb(:, jset), &
474 4800 : m, rab, sint, dsint, calculate_forces)
475 :
476 : CALL ab_contract(sab(sgfa:m1, sgfb:m2), sint, scon_a(:, sgfa:), scon_b(:, sgfb:), &
477 4800 : ncoa, ncob, nsgfa(iset), nsgfb(jset))
478 5120 : IF (calculate_forces) THEN
479 19200 : DO i = 1, 3
480 : CALL ab_contract(dsab(sgfa:m1, sgfb:m2, i), dsint(:, :, i), scon_a(:, sgfa:), &
481 19200 : scon_b(:, sgfb:), ncoa, ncob, nsgfa(iset), nsgfb(jset))
482 : END DO
483 : END IF
484 : END DO
485 : END DO
486 :
487 32 : DEALLOCATE (sint, dsint)
488 :
489 32 : CALL timestop(handle)
490 :
491 32 : END SUBROUTINE int_ra2m_ab_os
492 :
493 : ! **************************************************************************************************
494 : !> \brief calculate integrals (a,b,fa)
495 : !> \param abaint integral (a,b,fa)
496 : !> \param dabdaint derivative of abaint with respect to A
497 : !> \param rab distance vector between center A and B
498 : !> \param oba orbital basis at center A
499 : !> \param obb orbital basis at center B
500 : !> \param fba auxiliary basis set at center A
501 : !> \param calculate_forces ...
502 : !> \param debug integrals are debugged by recursive routines if requested
503 : !> \param dmax maximal deviation between integrals when debugging
504 : ! **************************************************************************************************
505 1148 : SUBROUTINE int_overlap_aba_os(abaint, dabdaint, rab, oba, obb, fba, &
506 : calculate_forces, debug, dmax)
507 :
508 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: abaint
509 : REAL(KIND=dp), ALLOCATABLE, &
510 : DIMENSION(:, :, :, :), OPTIONAL :: dabdaint
511 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
512 : TYPE(gto_basis_set_type), POINTER :: oba, obb, fba
513 : LOGICAL, INTENT(IN) :: calculate_forces, debug
514 : REAL(KIND=dp), INTENT(INOUT) :: dmax
515 :
516 : CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_aba_os'
517 :
518 : INTEGER :: handle, i, iset, jset, kaset, m1, m2, &
519 : m3, ncoa, ncob, ncoc, nseta, nsetb, &
520 : nsetca, sgfa, sgfb, sgfc
521 1148 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, lca_max, &
522 1148 : lca_min, npgfa, npgfb, npgfca, nsgfa, &
523 1148 : nsgfb, nsgfca
524 1148 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, first_sgfca
525 : REAL(KIND=dp) :: dab, dac, dbc
526 1148 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: saba
527 1148 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: sdaba
528 : REAL(KIND=dp), DIMENSION(3) :: ra, rac, rb, rbc
529 1148 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b, set_radius_ca
530 1148 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, rpgfca, scon_a, scon_b, &
531 1148 : scon_ca, zeta, zetb, zetca
532 :
533 1148 : CALL timeset(routineN, handle)
534 1148 : NULLIFY (la_max, la_min, lb_max, lb_min, lca_max, lca_min, npgfa, npgfb, &
535 1148 : npgfca, nsgfa, nsgfb, nsgfca)
536 1148 : NULLIFY (first_sgfa, first_sgfb, first_sgfca, set_radius_a, set_radius_b, &
537 1148 : set_radius_ca, rpgfa, rpgfb, rpgfca, scon_a, scon_b, scon_ca, &
538 1148 : zeta, zetb, zetca)
539 :
540 : ! basis ikind
541 1148 : first_sgfa => oba%first_sgf
542 1148 : la_max => oba%lmax
543 1148 : la_min => oba%lmin
544 1148 : npgfa => oba%npgf
545 1148 : nseta = oba%nset
546 1148 : nsgfa => oba%nsgf_set
547 1148 : rpgfa => oba%pgf_radius
548 1148 : set_radius_a => oba%set_radius
549 1148 : scon_a => oba%scon
550 1148 : zeta => oba%zet
551 : ! basis jkind
552 1148 : first_sgfb => obb%first_sgf
553 1148 : lb_max => obb%lmax
554 1148 : lb_min => obb%lmin
555 1148 : npgfb => obb%npgf
556 1148 : nsetb = obb%nset
557 1148 : nsgfb => obb%nsgf_set
558 1148 : rpgfb => obb%pgf_radius
559 1148 : set_radius_b => obb%set_radius
560 1148 : scon_b => obb%scon
561 1148 : zetb => obb%zet
562 :
563 : ! basis RI A
564 1148 : first_sgfca => fba%first_sgf
565 1148 : lca_max => fba%lmax
566 1148 : lca_min => fba%lmin
567 1148 : npgfca => fba%npgf
568 1148 : nsetca = fba%nset
569 1148 : nsgfca => fba%nsgf_set
570 1148 : rpgfca => fba%pgf_radius
571 1148 : set_radius_ca => fba%set_radius
572 1148 : scon_ca => fba%scon
573 1148 : zetca => fba%zet
574 :
575 4592 : dab = SQRT(SUM(rab**2))
576 :
577 1148 : abaint = 0.0_dp
578 1148 : IF (calculate_forces) THEN
579 462 : IF (PRESENT(dabdaint)) dabdaint = 0.0_dp
580 : END IF
581 :
582 2352 : DO iset = 1, nseta
583 :
584 1204 : ncoa = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
585 1204 : sgfa = first_sgfa(1, iset)
586 :
587 3836 : DO jset = 1, nsetb
588 :
589 1484 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
590 :
591 1484 : ncob = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
592 1484 : sgfb = first_sgfb(1, jset)
593 1484 : m1 = sgfa + nsgfa(iset) - 1
594 1484 : m2 = sgfb + nsgfb(jset) - 1
595 :
596 : ! calculate integrals abaint and derivative [d(a,b,a)/dA] dabdaint if requested
597 : rac = 0._dp
598 1484 : dac = 0._dp
599 5936 : rbc = -rab
600 : dbc = dab
601 :
602 14742 : DO kaset = 1, nsetca
603 :
604 12054 : IF (set_radius_b(jset) + set_radius_ca(kaset) < dab) CYCLE
605 :
606 7997 : ncoc = npgfca(kaset)*(ncoset(lca_max(kaset)) - ncoset(lca_min(kaset) - 1))
607 7997 : sgfc = first_sgfca(1, kaset)
608 7997 : m3 = sgfc + nsgfca(kaset) - 1
609 9481 : IF (ncoa*ncob*ncoc > 0) THEN
610 39985 : ALLOCATE (saba(ncoa, ncob, ncoc))
611 : ! integrals
612 7997 : IF (calculate_forces) THEN
613 20382 : ALLOCATE (sdaba(ncoa, ncob, ncoc, 3))
614 : CALL overlap_aab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
615 : lca_max(kaset), lca_min(kaset), npgfca(kaset), rpgfca(:, kaset), zetca(:, kaset), &
616 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
617 3397 : rab, saba=saba, daba=sdaba)
618 :
619 13588 : DO i = 1, 3
620 : CALL abc_contract(dabdaint(sgfa:m1, sgfb:m2, sgfc:m3, i), sdaba(1:ncoa, 1:ncob, 1:ncoc, i), &
621 : scon_a(:, sgfa:), scon_b(:, sgfb:), scon_ca(:, sgfc:), &
622 13588 : ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfca(kaset))
623 : END DO
624 :
625 3397 : DEALLOCATE (sdaba)
626 : ELSE
627 : CALL overlap_aab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
628 : lca_max(kaset), lca_min(kaset), npgfca(kaset), rpgfca(:, kaset), zetca(:, kaset), &
629 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
630 4600 : rab, saba=saba)
631 :
632 : END IF
633 : ! debug if requested
634 7997 : IF (debug) THEN
635 7 : ra = 0.0_dp
636 7 : rb = rab
637 : CALL overlap_abc_test(la_max(iset), npgfa(iset), zeta(:, iset), la_min(iset), &
638 : lb_max(jset), npgfb(jset), zetb(:, jset), lb_min(jset), &
639 : lca_max(kaset), npgfca(kaset), zetca(:, kaset), lca_min(kaset), &
640 7 : ra, rb, ra, saba, dmax)
641 : END IF
642 : CALL abc_contract(abaint(sgfa:m1, sgfb:m2, sgfc:m3), saba(1:ncoa, 1:ncob, 1:ncoc), &
643 : scon_a(:, sgfa:), scon_b(:, sgfb:), scon_ca(:, sgfc:), &
644 7997 : ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfca(kaset))
645 7997 : DEALLOCATE (saba)
646 : END IF
647 : END DO
648 : END DO
649 : END DO
650 :
651 1148 : CALL timestop(handle)
652 :
653 2296 : END SUBROUTINE int_overlap_aba_os
654 :
655 : ! **************************************************************************************************
656 : !> \brief calculate integrals (a,b,fb)
657 : !> \param abbint integral (a,b,fb)
658 : !> \param dabbint derivative of abbint with respect to A
659 : !> \param rab distance vector between center A and B
660 : !> \param oba orbital basis at center A
661 : !> \param obb orbital basis at center B
662 : !> \param fbb auxiliary basis set at center B
663 : !> \param calculate_forces ...
664 : !> \param debug integrals are debugged by recursive routines if requested
665 : !> \param dmax maximal deviation between integrals when debugging
666 : ! **************************************************************************************************
667 1148 : SUBROUTINE int_overlap_abb_os(abbint, dabbint, rab, oba, obb, fbb, calculate_forces, debug, dmax)
668 :
669 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: abbint
670 : REAL(KIND=dp), ALLOCATABLE, &
671 : DIMENSION(:, :, :, :), OPTIONAL :: dabbint
672 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
673 : TYPE(gto_basis_set_type), POINTER :: oba, obb, fbb
674 : LOGICAL, INTENT(IN) :: calculate_forces, debug
675 : REAL(KIND=dp), INTENT(INOUT) :: dmax
676 :
677 : CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_abb_os'
678 :
679 : INTEGER :: handle, i, iset, jset, kbset, m1, m2, &
680 : m3, ncoa, ncob, ncoc, nseta, nsetb, &
681 : nsetcb, sgfa, sgfb, sgfc
682 1148 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, lcb_max, &
683 1148 : lcb_min, npgfa, npgfb, npgfcb, nsgfa, &
684 1148 : nsgfb, nsgfcb
685 1148 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, first_sgfcb
686 : REAL(KIND=dp) :: dab, dac, dbc
687 1148 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: sabb
688 1148 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: sdabb
689 : REAL(KIND=dp), DIMENSION(3) :: ra, rac, rb, rbc
690 1148 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b, set_radius_cb
691 1148 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, rpgfcb, scon_a, scon_b, &
692 1148 : scon_cb, zeta, zetb, zetcb
693 :
694 1148 : CALL timeset(routineN, handle)
695 1148 : NULLIFY (la_max, la_min, lb_max, lb_min, lcb_max, lcb_min, npgfa, npgfb, &
696 1148 : npgfcb, nsgfa, nsgfb, nsgfcb)
697 1148 : NULLIFY (first_sgfa, first_sgfb, first_sgfcb, set_radius_a, set_radius_b, &
698 1148 : set_radius_cb, rpgfa, rpgfb, rpgfcb, scon_a, scon_b, scon_cb, &
699 1148 : zeta, zetb, zetcb)
700 :
701 : ! basis ikind
702 1148 : first_sgfa => oba%first_sgf
703 1148 : la_max => oba%lmax
704 1148 : la_min => oba%lmin
705 1148 : npgfa => oba%npgf
706 1148 : nseta = oba%nset
707 1148 : nsgfa => oba%nsgf_set
708 1148 : rpgfa => oba%pgf_radius
709 1148 : set_radius_a => oba%set_radius
710 1148 : scon_a => oba%scon
711 1148 : zeta => oba%zet
712 : ! basis jkind
713 1148 : first_sgfb => obb%first_sgf
714 1148 : lb_max => obb%lmax
715 1148 : lb_min => obb%lmin
716 1148 : npgfb => obb%npgf
717 1148 : nsetb = obb%nset
718 1148 : nsgfb => obb%nsgf_set
719 1148 : rpgfb => obb%pgf_radius
720 1148 : set_radius_b => obb%set_radius
721 1148 : scon_b => obb%scon
722 1148 : zetb => obb%zet
723 :
724 : ! basis RI B
725 1148 : first_sgfcb => fbb%first_sgf
726 1148 : lcb_max => fbb%lmax
727 1148 : lcb_min => fbb%lmin
728 1148 : npgfcb => fbb%npgf
729 1148 : nsetcb = fbb%nset
730 1148 : nsgfcb => fbb%nsgf_set
731 1148 : rpgfcb => fbb%pgf_radius
732 1148 : set_radius_cb => fbb%set_radius
733 1148 : scon_cb => fbb%scon
734 1148 : zetcb => fbb%zet
735 :
736 4592 : dab = SQRT(SUM(rab**2))
737 :
738 1148 : abbint = 0.0_dp
739 1148 : IF (calculate_forces) THEN
740 462 : IF (PRESENT(dabbint)) dabbint = 0.0_dp
741 : END IF
742 :
743 2352 : DO iset = 1, nseta
744 :
745 1204 : ncoa = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
746 1204 : sgfa = first_sgfa(1, iset)
747 :
748 3836 : DO jset = 1, nsetb
749 :
750 1484 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
751 :
752 1484 : ncob = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
753 1484 : sgfb = first_sgfb(1, jset)
754 1484 : m1 = sgfa + nsgfa(iset) - 1
755 1484 : m2 = sgfb + nsgfb(jset) - 1
756 :
757 : ! calculate integrals abbint and derivative [d(a,b,b)/dA] dabbint if requested
758 1484 : rac = rab
759 1484 : dac = dab
760 : rbc = 0._dp
761 1484 : dbc = 0._dp
762 :
763 14742 : DO kbset = 1, nsetcb
764 :
765 12054 : IF (set_radius_a(iset) + set_radius_cb(kbset) < dab) CYCLE
766 :
767 7997 : ncoc = npgfcb(kbset)*(ncoset(lcb_max(kbset)) - ncoset(lcb_min(kbset) - 1))
768 7997 : sgfc = first_sgfcb(1, kbset)
769 7997 : m3 = sgfc + nsgfcb(kbset) - 1
770 9481 : IF (ncoa*ncob*ncoc > 0) THEN
771 39985 : ALLOCATE (sabb(ncoa, ncob, ncoc))
772 7997 : IF (calculate_forces) THEN
773 20382 : ALLOCATE (sdabb(ncoa, ncob, ncoc, 3))
774 : CALL overlap_abb(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
775 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
776 : lcb_max(kbset), lcb_min(kbset), npgfcb(kbset), rpgfcb(:, kbset), zetcb(:, kbset), &
777 3397 : rab, sabb, sdabb)
778 :
779 13588 : DO i = 1, 3
780 : CALL abc_contract(dabbint(sgfa:m1, sgfb:m2, sgfc:m3, i), sdabb(1:ncoa, 1:ncob, 1:ncoc, i), &
781 : scon_a(:, sgfa:), scon_b(:, sgfb:), scon_cb(:, sgfc:), &
782 13588 : ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfcb(kbset))
783 : END DO
784 3397 : DEALLOCATE (sdabb)
785 : ELSE
786 : CALL overlap_abb(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
787 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
788 : lcb_max(kbset), lcb_min(kbset), npgfcb(kbset), rpgfcb(:, kbset), zetcb(:, kbset), &
789 4600 : rab, sabb)
790 : END IF
791 : ! debug if requested
792 7997 : IF (debug) THEN
793 7 : ra = 0.0_dp
794 7 : rb = rab
795 : CALL overlap_abc_test(la_max(iset), npgfa(iset), zeta(:, iset), la_min(iset), &
796 : lb_max(jset), npgfb(jset), zetb(:, jset), lb_min(jset), &
797 : lcb_max(kbset), npgfcb(kbset), zetcb(:, kbset), lcb_min(kbset), &
798 7 : ra, rb, rb, sabb, dmax)
799 : END IF
800 :
801 : CALL abc_contract(abbint(sgfa:m1, sgfb:m2, sgfc:m3), sabb(1:ncoa, 1:ncob, 1:ncoc), &
802 : scon_a(:, sgfa:), scon_b(:, sgfb:), scon_cb(:, sgfc:), &
803 7997 : ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfcb(kbset))
804 7997 : DEALLOCATE (sabb)
805 : END IF
806 : END DO
807 :
808 : END DO
809 : END DO
810 :
811 1148 : CALL timestop(handle)
812 :
813 2296 : END SUBROUTINE int_overlap_abb_os
814 :
815 : ! **************************************************************************************************
816 : !> \brief calculate overlap integrals (aa,bb)
817 : !> \param saabb integral (aa,bb)
818 : !> \param oba orbital basis at center A
819 : !> \param obb orbital basis at center B
820 : !> \param rab ...
821 : !> \param debug integrals are debugged by recursive routines if requested
822 : !> \param dmax maximal deviation between integrals when debugging
823 : ! **************************************************************************************************
824 9 : SUBROUTINE int_overlap_aabb_os(saabb, oba, obb, rab, debug, dmax)
825 :
826 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: saabb
827 : TYPE(gto_basis_set_type), POINTER :: oba, obb
828 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
829 : LOGICAL, INTENT(IN) :: debug
830 : REAL(KIND=dp), INTENT(INOUT) :: dmax
831 :
832 : CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_aabb_os'
833 :
834 : INTEGER :: handle, iset, isgfa1, jset, jsgfa2, kset, ksgfb1, lset, lsgfb2, m1, m2, m3, m4, &
835 : maxcoa, maxcob, ncoa1, ncoa2, ncob1, ncob2, nseta, nsetb, sgfa1, sgfa2, sgfb1, sgfb2
836 9 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
837 9 : npgfb, nsgfa, nsgfb
838 9 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
839 : REAL(KIND=dp) :: dab, ra(3), rb(3)
840 9 : REAL(KIND=dp), ALLOCATABLE, TARGET :: buffer(:)
841 9 : REAL(KIND=dp), CONTIGUOUS, POINTER :: sint(:, :, :, :)
842 9 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
843 :
844 9 : CALL timeset(routineN, handle)
845 9 : NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, &
846 9 : first_sgfa, first_sgfb, rpgfa, rpgfb, &
847 9 : sphi_a, sphi_b, zeta, zetb)
848 :
849 : ! basis ikind
850 9 : first_sgfa => oba%first_sgf
851 9 : la_max => oba%lmax
852 9 : la_min => oba%lmin
853 9 : npgfa => oba%npgf
854 9 : nseta = oba%nset
855 9 : nsgfa => oba%nsgf_set
856 9 : rpgfa => oba%pgf_radius
857 9 : sphi_a => oba%sphi
858 9 : zeta => oba%zet
859 : ! basis jkind
860 9 : first_sgfb => obb%first_sgf
861 9 : lb_max => obb%lmax
862 9 : lb_min => obb%lmin
863 9 : npgfb => obb%npgf
864 9 : nsetb = obb%nset
865 9 : nsgfb => obb%nsgf_set
866 9 : rpgfb => obb%pgf_radius
867 9 : sphi_b => obb%sphi
868 9 : zetb => obb%zet
869 :
870 : ! Reserve the largest quartet once; each view uses its own contiguous layout.
871 9 : CALL get_gto_basis_set(oba, maxco=maxcoa)
872 9 : CALL get_gto_basis_set(obb, maxco=maxcob)
873 27 : ALLOCATE (buffer(INT(maxcoa, int_8)**2*INT(maxcob, int_8)**2))
874 :
875 36 : dab = SQRT(SUM(rab**2))
876 :
877 18 : DO iset = 1, nseta
878 :
879 9 : ncoa1 = npgfa(iset)*ncoset(la_max(iset))
880 9 : sgfa1 = first_sgfa(1, iset)
881 9 : m1 = sgfa1 + nsgfa(iset) - 1
882 :
883 27 : DO jset = iset, nseta
884 :
885 9 : ncoa2 = npgfa(jset)*ncoset(la_max(jset))
886 9 : sgfa2 = first_sgfa(1, jset)
887 9 : m2 = sgfa2 + nsgfa(jset) - 1
888 :
889 27 : DO kset = 1, nsetb
890 :
891 9 : ncob1 = npgfb(kset)*ncoset(lb_max(kset))
892 9 : sgfb1 = first_sgfb(1, kset)
893 9 : m3 = sgfb1 + nsgfb(kset) - 1
894 :
895 27 : DO lset = kset, nsetb
896 :
897 9 : ncob2 = npgfb(lset)*ncoset(lb_max(lset))
898 9 : sgfb2 = first_sgfb(1, lset)
899 9 : m4 = sgfb2 + nsgfb(lset) - 1
900 :
901 : sint(1:ncoa1, 1:ncoa2, 1:ncob1, 1:ncob2) => &
902 9 : buffer(1:INT(ncoa1, int_8)*ncoa2*ncob1*ncob2)
903 : ! calculate integrals
904 : CALL overlap_aabb(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
905 : la_max(jset), la_min(jset), npgfa(jset), rpgfa(:, jset), zeta(:, jset), &
906 : lb_max(kset), lb_min(kset), npgfb(kset), rpgfb(:, kset), zetb(:, kset), &
907 : lb_max(lset), lb_min(lset), npgfb(lset), rpgfb(:, lset), zetb(:, lset), &
908 9 : iset == jset, kset == lset, rab, dab, sint)
909 : ! debug if requested
910 9 : IF (debug) THEN
911 3 : ra = 0.0_dp
912 3 : rb = rab
913 : CALL overlap_aabb_test(la_max(iset), la_min(iset), npgfa(iset), zeta(:, iset), &
914 : la_max(jset), la_min(jset), npgfa(jset), zeta(:, jset), &
915 : lb_max(kset), lb_min(kset), npgfb(kset), zetb(:, kset), &
916 : lb_max(lset), lb_min(lset), npgfb(lset), zetb(:, lset), &
917 3 : ra, rb, sint, dmax)
918 : END IF
919 :
920 : CALL abcd_contract(saabb(sgfa1:m1, sgfa2:m2, sgfb1:m3, sgfb2:m4), sint, sphi_a(:, sgfa1:), &
921 : sphi_a(:, sgfa2:), sphi_b(:, sgfb1:), sphi_b(:, sgfb2:), ncoa1, ncoa2, &
922 9 : ncob1, ncob2, nsgfa(iset), nsgfa(jset), nsgfb(kset), nsgfb(lset))
923 :
924 : ! account for the fact that some integrals are alike
925 54 : DO isgfa1 = sgfa1, m1
926 189 : DO jsgfa2 = sgfa2, m2
927 756 : DO ksgfb1 = sgfb1, m3
928 3024 : DO lsgfb2 = sgfb2, m4
929 2304 : saabb(jsgfa2, isgfa1, ksgfb1, lsgfb2) = saabb(isgfa1, jsgfa2, ksgfb1, lsgfb2)
930 2304 : saabb(isgfa1, jsgfa2, lsgfb2, ksgfb1) = saabb(isgfa1, jsgfa2, ksgfb1, lsgfb2)
931 2880 : saabb(jsgfa2, isgfa1, lsgfb2, ksgfb1) = saabb(isgfa1, jsgfa2, ksgfb1, lsgfb2)
932 : END DO
933 : END DO
934 : END DO
935 : END DO
936 :
937 : END DO
938 : END DO
939 : END DO
940 : END DO
941 :
942 9 : DEALLOCATE (buffer)
943 9 : CALL timestop(handle)
944 :
945 9 : END SUBROUTINE int_overlap_aabb_os
946 :
947 : END MODULE generic_os_integrals
|
