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