Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Debugs Obara-Saika integral matrices
10 : !> \par History
11 : !> created [07.2014]
12 : !> \authors Dorothea Golze
13 : ! **************************************************************************************************
14 : MODULE debug_os_integrals
15 :
16 : USE ai_overlap, ONLY: overlap
17 : USE ai_overlap3, ONLY: overlap3
18 : USE ai_overlap3_debug, ONLY: init_os_overlap3,&
19 : os_overlap3
20 : USE ai_overlap_aabb, ONLY: overlap_aabb
21 : USE ai_overlap_debug, ONLY: init_os_overlap2,&
22 : os_overlap2
23 : USE kinds, ONLY: dp
24 : USE orbital_pointers, ONLY: coset,&
25 : indco,&
26 : ncoset
27 : #include "./base/base_uses.f90"
28 :
29 : IMPLICIT NONE
30 :
31 : PRIVATE
32 :
33 : ! **************************************************************************************************
34 :
35 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'debug_os_integrals'
36 :
37 : PUBLIC :: overlap_ab_test, overlap_abc_test, overlap_aabb_test
38 :
39 : ! **************************************************************************************************
40 :
41 : CONTAINS
42 :
43 : ! ***************************************************************************************************
44 : !> \brief recursive test routines for integral (a,b) for only two exponents
45 : ! **************************************************************************************************
46 0 : SUBROUTINE overlap_ab_test_simple()
47 :
48 : INTEGER :: ia1, iax, iay, iaz, ib1, ibx, iby, ibz, &
49 : la_max, la_min, lb_max, lb_min, lds, &
50 : ma, maxl, mb
51 : INTEGER, DIMENSION(3) :: na, nb
52 : REAL(KIND=dp) :: dab, dmax, res1, xa, xb
53 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sab
54 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: swork
55 : REAL(KIND=dp), DIMENSION(1) :: rpgfa, rpgfb, xa_work, xb_work
56 : REAL(KIND=dp), DIMENSION(3) :: A, B, rab
57 :
58 0 : xa = 0.783300000000_dp ! exponents
59 0 : xb = 1.239648746700_dp
60 :
61 0 : A = (/0.329309000000_dp, 0.28408240000_dp, 0.28408240000_dp/) !* bohr !positions
62 0 : B = (/0.983983000000_dp, 0.00453720000_dp, 0.00432740000_dp/) !* bohr
63 :
64 0 : la_min = 0
65 0 : lb_min = 0
66 :
67 0 : la_max = 3
68 0 : lb_max = 4
69 :
70 : !---------------------------------------
71 0 : ALLOCATE (sab(ncoset(la_max), ncoset(lb_max)))
72 :
73 0 : maxl = MAX(la_max, lb_max)
74 0 : lds = ncoset(maxl)
75 0 : ALLOCATE (swork(lds, lds, 1))
76 0 : sab = 0._dp
77 0 : rab(:) = B(:) - A(:)
78 0 : dab = SQRT(DOT_PRODUCT(rab, rab))
79 0 : xa_work(1) = xa
80 0 : xb_work(1) = xb
81 0 : rpgfa = 20._dp
82 0 : rpgfb = 20._dp
83 : CALL overlap(la_max_set=la_max, la_min_set=la_min, npgfa=1, rpgfa=rpgfa, zeta=xa_work, &
84 : lb_max_set=lb_max, lb_min_set=lb_min, npgfb=1, rpgfb=rpgfb, zetb=xb_work, &
85 0 : rab=rab, dab=dab, sab=sab, da_max_set=0, return_derivatives=.FALSE., s=swork, lds=lds)
86 : !---------------------------------------
87 :
88 0 : CALL init_os_overlap2(xa, xb, A, B)
89 :
90 0 : dmax = 0._dp
91 0 : DO ma = la_min, la_max
92 0 : DO mb = lb_min, lb_max
93 0 : DO iax = 0, ma
94 0 : DO iay = 0, ma - iax
95 0 : iaz = ma - iax - iay
96 0 : na(1) = iax; na(2) = iay; na(3) = iaz
97 0 : ia1 = coset(iax, iay, iaz)
98 0 : DO ibx = 0, mb
99 0 : DO iby = 0, mb - ibx
100 0 : ibz = mb - ibx - iby
101 0 : nb(1) = ibx; nb(2) = iby; nb(3) = ibz
102 0 : ib1 = coset(ibx, iby, ibz)
103 0 : res1 = os_overlap2(na, nb)
104 : ! write(*,*) "la, lb,na, nb, res1", ma, mb, na, nb, res1
105 : ! write(*,*) "sab ia1, ib1", ia1, ib1, sab(ia1,ib1)
106 0 : dmax = MAX(dmax, ABS(res1 - sab(ia1, ib1)))
107 : END DO
108 : END DO
109 : END DO
110 : END DO
111 : END DO
112 : END DO
113 :
114 0 : DEALLOCATE (sab, swork)
115 :
116 0 : END SUBROUTINE overlap_ab_test_simple
117 :
118 : ! ***************************************************************************************************
119 : !> \brief recursive test routines for integral (a,b)
120 : !> \param la_max ...
121 : !> \param la_min ...
122 : !> \param npgfa ...
123 : !> \param zeta ...
124 : !> \param lb_max ...
125 : !> \param lb_min ...
126 : !> \param npgfb ...
127 : !> \param zetb ...
128 : !> \param ra ...
129 : !> \param rb ...
130 : !> \param sab ...
131 : !> \param dmax ...
132 : ! **************************************************************************************************
133 47 : SUBROUTINE overlap_ab_test(la_max, la_min, npgfa, zeta, lb_max, lb_min, npgfb, zetb, &
134 47 : ra, rb, sab, dmax)
135 :
136 : INTEGER, INTENT(IN) :: la_max, la_min, npgfa
137 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta
138 : INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
139 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb
140 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rb
141 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sab
142 : REAL(KIND=dp), INTENT(INOUT) :: dmax
143 :
144 : INTEGER :: coa, cob, ia1, iax, iay, iaz, ib1, ibx, &
145 : iby, ibz, ipgf, jpgf, ma, mb
146 : INTEGER, DIMENSION(3) :: na, nb
147 : REAL(KIND=dp) :: res1, res2, xa, xb
148 : REAL(KIND=dp), DIMENSION(3) :: A, B
149 :
150 47 : coa = 0
151 100 : DO ipgf = 1, npgfa
152 53 : cob = 0
153 148 : DO jpgf = 1, npgfb
154 95 : xa = zeta(ipgf) !exponents
155 95 : xb = zetb(jpgf)
156 95 : A = ra !positions
157 95 : B = rb
158 95 : CALL init_os_overlap2(xa, xb, A, B)
159 314 : DO ma = la_min, la_max
160 827 : DO mb = lb_min, lb_max
161 1664 : DO iax = 0, ma
162 2939 : DO iay = 0, ma - iax
163 1494 : iaz = ma - iax - iay
164 1494 : na(1) = iax; na(2) = iay; na(3) = iaz
165 1494 : ia1 = coset(iax, iay, iaz)
166 5229 : DO ibx = 0, mb
167 8904 : DO iby = 0, mb - ibx
168 4607 : ibz = mb - ibx - iby
169 4607 : nb(1) = ibx; nb(2) = iby; nb(3) = ibz
170 4607 : ib1 = coset(ibx, iby, ibz)
171 4607 : res1 = os_overlap2(na, nb)
172 4607 : res2 = sab(coa + ia1, cob + ib1)
173 7410 : dmax = MAX(dmax, ABS(res1 - res2))
174 : END DO
175 : END DO
176 : END DO
177 : END DO
178 : END DO
179 : END DO
180 148 : cob = cob + ncoset(lb_max)
181 : END DO
182 100 : coa = coa + ncoset(la_max)
183 : END DO
184 : !WRITE(*,*) "dmax overlap_ab_test", dmax
185 :
186 47 : END SUBROUTINE overlap_ab_test
187 :
188 : ! ***************************************************************************************************
189 : !> \brief recursive test routines for integral (a,b,c) for only three exponents
190 : ! **************************************************************************************************
191 0 : SUBROUTINE overlap_abc_test_simple()
192 :
193 : INTEGER :: ia1, iax, iay, iaz, ib1, ibx, iby, ibz, &
194 : ic1, icx, icy, icz, la_max, la_min, &
195 : lb_max, lb_min, lc_max, lc_min, ma, &
196 : mb, mc
197 : INTEGER, DIMENSION(3) :: na, nb, nc
198 : REAL(KIND=dp) :: dab, dac, dbc, dmax, res1, xa, xb, xc
199 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: sabc
200 : REAL(KIND=dp), DIMENSION(1) :: rpgfa, rpgfb, rpgfc, xa_work, xb_work, &
201 : xc_work
202 : REAL(KIND=dp), DIMENSION(3) :: A, B, C, rab, rac, rbc
203 :
204 0 : xa = 0.783300000000_dp ! exponents
205 0 : xb = 1.239648746700_dp
206 0 : xc = 0.548370000000_dp
207 :
208 0 : A = (/0.329309000000_dp, 0.28408240000_dp, 0.28408240000_dp/) !* bohr !positions
209 0 : B = (/0.983983000000_dp, 0.00453720000_dp, 0.00432740000_dp/) !* bohr
210 0 : C = (/0.032380000000_dp, 1.23470000000_dp, 0.11137400000_dp/) !* bohr
211 :
212 0 : la_min = 0
213 0 : lb_min = 0
214 0 : lc_min = 0
215 :
216 0 : la_max = 0
217 0 : lb_max = 0
218 0 : lc_max = 1
219 :
220 : !---------------------------------------
221 0 : rab(:) = B(:) - A(:)
222 0 : dab = SQRT(DOT_PRODUCT(rab, rab))
223 0 : rac(:) = C(:) - A(:)
224 0 : dac = SQRT(DOT_PRODUCT(rac, rac))
225 0 : rbc(:) = C(:) - B(:)
226 0 : dbc = SQRT(DOT_PRODUCT(rbc, rbc))
227 0 : ALLOCATE (sabc(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
228 0 : xa_work(1) = xa
229 0 : xb_work(1) = xb
230 0 : xc_work(1) = xc
231 0 : rpgfa = 20._dp
232 0 : rpgfb = 20._dp
233 0 : rpgfc = 20._dp
234 0 : sabc = 0._dp
235 :
236 : CALL overlap3(la_max_set=la_max, npgfa=1, zeta=xa_work, rpgfa=rpgfa, la_min_set=la_min, &
237 : lb_max_set=lb_max, npgfb=1, zetb=xb_work, rpgfb=rpgfb, lb_min_set=lb_min, &
238 : lc_max_set=lc_max, npgfc=1, zetc=xc_work, rpgfc=rpgfc, lc_min_set=lc_min, &
239 0 : rab=rab, dab=dab, rac=rac, dac=dac, rbc=rbc, dbc=dbc, sabc=sabc)
240 :
241 : !---------------------------------------
242 :
243 0 : CALL init_os_overlap3(xa, xb, xc, A, B, C)
244 :
245 0 : dmax = 0._dp
246 0 : DO ma = la_min, la_max
247 0 : DO mc = lc_min, lc_max
248 0 : DO mb = lb_min, lb_max
249 0 : DO iax = 0, ma
250 0 : DO iay = 0, ma - iax
251 0 : iaz = ma - iax - iay
252 0 : na(1) = iax; na(2) = iay; na(3) = iaz
253 0 : ia1 = coset(iax, iay, iaz)
254 0 : DO icx = 0, mc
255 0 : DO icy = 0, mc - icx
256 0 : icz = mc - icx - icy
257 0 : nc(1) = icx; nc(2) = icy; nc(3) = icz
258 0 : ic1 = coset(icx, icy, icz)
259 0 : DO ibx = 0, mb
260 0 : DO iby = 0, mb - ibx
261 0 : ibz = mb - ibx - iby
262 0 : nb(1) = ibx; nb(2) = iby; nb(3) = ibz
263 0 : ib1 = coset(ibx, iby, ibz)
264 0 : res1 = os_overlap3(na, nc, nb)
265 : !write(*,*) "la, lc, lb,na,nc, nb, res1", ma, mc, mb, na, nc, nb, res1
266 : !write(*,*) "sabc ia1, ib1, ic1", ia1, ib1, ic1, sabc(ia1,ib1,ic1)
267 0 : dmax = MAX(dmax, ABS(res1 - sabc(ia1, ib1, ic1)))
268 : END DO
269 : END DO
270 : END DO
271 : END DO
272 : END DO
273 : END DO
274 : END DO
275 : END DO
276 : END DO
277 :
278 0 : DEALLOCATE (sabc)
279 :
280 0 : END SUBROUTINE overlap_abc_test_simple
281 :
282 : ! ***************************************************************************************************
283 : !> \brief recursive test routines for integral (a,b,c)
284 : !> \param la_max ...
285 : !> \param npgfa ...
286 : !> \param zeta ...
287 : !> \param la_min ...
288 : !> \param lb_max ...
289 : !> \param npgfb ...
290 : !> \param zetb ...
291 : !> \param lb_min ...
292 : !> \param lc_max ...
293 : !> \param npgfc ...
294 : !> \param zetc ...
295 : !> \param lc_min ...
296 : !> \param ra ...
297 : !> \param rb ...
298 : !> \param rc ...
299 : !> \param sabc ...
300 : !> \param dmax ...
301 : ! **************************************************************************************************
302 14 : SUBROUTINE overlap_abc_test(la_max, npgfa, zeta, la_min, &
303 14 : lb_max, npgfb, zetb, lb_min, &
304 14 : lc_max, npgfc, zetc, lc_min, &
305 14 : ra, rb, rc, sabc, dmax)
306 :
307 : INTEGER, INTENT(IN) :: la_max, npgfa
308 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta
309 : INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
310 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb
311 : INTEGER, INTENT(IN) :: lb_min, lc_max, npgfc
312 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetc
313 : INTEGER, INTENT(IN) :: lc_min
314 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rb, rc
315 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: sabc
316 : REAL(KIND=dp), INTENT(INOUT) :: dmax
317 :
318 : INTEGER :: coa, cob, coc, ia1, iax, iay, iaz, ib1, &
319 : ibx, iby, ibz, ic1, icx, icy, icz, &
320 : ipgf, jpgf, kpgf, ma, mb, mc
321 : INTEGER, DIMENSION(3) :: na, nb, nc
322 : REAL(KIND=dp) :: res1, res2, xa, xb, xc
323 : REAL(KIND=dp), DIMENSION(3) :: A, B, C
324 :
325 14 : coa = 0
326 112 : DO ipgf = 1, npgfa
327 98 : cob = 0
328 784 : DO jpgf = 1, npgfb
329 686 : coc = 0
330 1372 : DO kpgf = 1, npgfc
331 :
332 686 : xa = zeta(ipgf) ! exponents
333 686 : xb = zetb(jpgf)
334 686 : xc = zetc(kpgf)
335 :
336 686 : A = Ra !positions
337 686 : B = Rb
338 686 : C = Rc
339 :
340 686 : CALL init_os_overlap3(xa, xb, xc, A, B, C)
341 :
342 2058 : DO ma = la_min, la_max
343 5586 : DO mc = lc_min, lc_max
344 11956 : DO mb = lb_min, lb_max
345 21168 : DO iax = 0, ma
346 31752 : DO iay = 0, ma - iax
347 14112 : iaz = ma - iax - iay
348 14112 : na(1) = iax; na(2) = iay; na(3) = iaz
349 14112 : ia1 = coset(iax, iay, iaz)
350 52920 : DO icx = 0, mc
351 90944 : DO icy = 0, mc - icx
352 48608 : icz = mc - icx - icy
353 48608 : nc(1) = icx; nc(2) = icy; nc(3) = icz
354 48608 : ic1 = coset(icx, icy, icz)
355 149744 : DO ibx = 0, mb
356 218736 : DO iby = 0, mb - ibx
357 97216 : ibz = mb - ibx - iby
358 97216 : nb(1) = ibx; nb(2) = iby; nb(3) = ibz
359 97216 : ib1 = coset(ibx, iby, ibz)
360 97216 : res1 = os_overlap3(na, nc, nb)
361 97216 : res2 = sabc(coa + ia1, cob + ib1, coc + ic1)
362 170128 : dmax = MAX(dmax, ABS(res1 - res2))
363 : !IF(dmax > 1.E-10) WRITE(*,*) "dmax in loop", dmax
364 : END DO
365 : END DO
366 : END DO
367 : END DO
368 : END DO
369 : END DO
370 : END DO
371 : END DO
372 : END DO
373 1372 : coc = coc + ncoset(lc_max)
374 : END DO
375 784 : cob = cob + ncoset(lb_max)
376 : END DO
377 112 : coa = coa + ncoset(la_max)
378 : END DO
379 : !WRITE(*,*) "dmax abc", dmax
380 :
381 14 : END SUBROUTINE overlap_abc_test
382 :
383 : ! ***************************************************************************************************
384 : !> \brief recursive test routines for integral (aa,bb) for only four exponents
385 : ! **************************************************************************************************
386 0 : SUBROUTINE overlap_aabb_test_simple()
387 :
388 : INTEGER :: i, iax, iay, iaz, ibx, iby, ibz, j, k, l, la_max, la_max1, la_max2, la_min, &
389 : la_min1, la_min2, lb_max, lb_max1, lb_max2, lb_min, lb_min1, lb_min2, lds, ma, maxl, mb
390 : INTEGER, DIMENSION(3) :: na, naa, nb, nbb
391 : REAL(KIND=dp) :: dab, dmax, res1, xa, xa1, xa2, xb, xb1, &
392 : xb2
393 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: swork
394 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: saabb
395 : REAL(KIND=dp), DIMENSION(1) :: rpgfa1, rpgfa2, rpgfb1, rpgfb2, &
396 : xa_work1, xa_work2, xb_work1, xb_work2
397 : REAL(KIND=dp), DIMENSION(3) :: A, B, rab
398 :
399 0 : xa1 = 0.783300000000_dp ! exponents
400 0 : xb1 = 1.239648746700_dp
401 0 : xa2 = 0.3400000000_dp ! exponents
402 0 : xb2 = 2.76_dp
403 :
404 0 : A = (/0.329309000000_dp, 0.28408240000_dp, 0.28408240000_dp/) !* bohr !positions
405 0 : B = (/0.983983000000_dp, 0.00453720000_dp, 0.00432740000_dp/) !* bohr
406 :
407 0 : la_min1 = 0
408 0 : la_min2 = 0
409 0 : lb_min1 = 3
410 0 : lb_min2 = 1
411 :
412 0 : la_max1 = 1
413 0 : la_max2 = 2
414 0 : lb_max1 = 3
415 0 : lb_max2 = 4
416 :
417 : !---------------------------------------
418 0 : ALLOCATE (saabb(ncoset(la_max1), ncoset(la_max2), ncoset(lb_max1), ncoset(lb_max2)))
419 :
420 0 : maxl = MAX(la_max1 + la_max2, lb_max1 + lb_max2)
421 0 : lds = ncoset(maxl)
422 0 : ALLOCATE (swork(lds, lds))
423 0 : saabb = 0._dp
424 0 : rab(:) = B(:) - A(:)
425 0 : dab = SQRT(DOT_PRODUCT(rab, rab))
426 0 : xa_work1(1) = xa1
427 0 : xa_work2(1) = xa2
428 0 : xb_work1(1) = xb1
429 0 : xb_work2(1) = xb2
430 0 : rpgfa1 = 20._dp
431 0 : rpgfa2 = 20._dp
432 0 : rpgfb1 = 20._dp
433 0 : rpgfb2 = 20._dp
434 : CALL overlap_aabb(la_max_set1=la_max1, la_min_set1=la_min1, npgfa1=1, rpgfa1=rpgfa1, zeta1=xa_work1, &
435 : la_max_set2=la_max2, la_min_set2=la_min2, npgfa2=1, rpgfa2=rpgfa2, zeta2=xa_work2, &
436 : lb_max_set1=lb_max1, lb_min_set1=lb_min1, npgfb1=1, rpgfb1=rpgfb1, zetb1=xb_work1, &
437 : lb_max_set2=lb_max2, lb_min_set2=lb_min2, npgfb2=1, rpgfb2=rpgfb2, zetb2=xb_work2, &
438 0 : asets_equal=.FALSE., bsets_equal=.FALSE., rab=rab, dab=dab, saabb=saabb, s=swork, lds=lds)
439 : !---------------------------------------
440 :
441 0 : xa = xa1 + xa2
442 0 : xb = xb1 + xb2
443 0 : la_min = la_min1 + la_min2
444 0 : la_max = la_max1 + la_max2
445 0 : lb_min = lb_min1 + lb_min2
446 0 : lb_max = lb_max1 + lb_max2
447 :
448 0 : CALL init_os_overlap2(xa, xb, A, B)
449 :
450 0 : dmax = 0._dp
451 0 : DO ma = la_min, la_max
452 0 : DO mb = lb_min, lb_max
453 0 : DO iax = 0, ma
454 0 : DO iay = 0, ma - iax
455 0 : iaz = ma - iax - iay
456 0 : na(1) = iax; na(2) = iay; na(3) = iaz
457 0 : DO ibx = 0, mb
458 0 : DO iby = 0, mb - ibx
459 0 : ibz = mb - ibx - iby
460 0 : nb(1) = ibx; nb(2) = iby; nb(3) = ibz
461 0 : res1 = os_overlap2(na, nb)
462 0 : DO i = ncoset(la_min1 - 1) + 1, ncoset(la_max1)
463 0 : DO j = ncoset(la_min2 - 1) + 1, ncoset(la_max2)
464 0 : naa = indco(1:3, i) + indco(1:3, j)
465 0 : DO k = ncoset(lb_min1 - 1) + 1, ncoset(lb_max1)
466 0 : DO l = ncoset(lb_min2 - 1) + 1, ncoset(lb_max2)
467 0 : nbb = indco(1:3, k) + indco(1:3, l)
468 0 : IF (ALL(na == naa) .AND. ALL(nb == nbb)) THEN
469 0 : dmax = MAX(dmax, ABS(res1 - saabb(i, j, k, l)))
470 : END IF
471 : END DO
472 : END DO
473 : END DO
474 : END DO
475 : END DO
476 : END DO
477 : END DO
478 : END DO
479 : END DO
480 : END DO
481 :
482 0 : DEALLOCATE (saabb, swork)
483 :
484 0 : END SUBROUTINE overlap_aabb_test_simple
485 :
486 : ! ***************************************************************************************************
487 : !> \brief recursive test routines for integral (aa,bb)
488 : !> \param la_max1 ...
489 : !> \param la_min1 ...
490 : !> \param npgfa1 ...
491 : !> \param zeta1 ...
492 : !> \param la_max2 ...
493 : !> \param la_min2 ...
494 : !> \param npgfa2 ...
495 : !> \param zeta2 ...
496 : !> \param lb_max1 ...
497 : !> \param lb_min1 ...
498 : !> \param npgfb1 ...
499 : !> \param zetb1 ...
500 : !> \param lb_max2 ...
501 : !> \param lb_min2 ...
502 : !> \param npgfb2 ...
503 : !> \param zetb2 ...
504 : !> \param ra ...
505 : !> \param rb ...
506 : !> \param saabb ...
507 : !> \param dmax ...
508 : ! **************************************************************************************************
509 3 : SUBROUTINE overlap_aabb_test(la_max1, la_min1, npgfa1, zeta1, &
510 3 : la_max2, la_min2, npgfa2, zeta2, &
511 3 : lb_max1, lb_min1, npgfb1, zetb1, &
512 3 : lb_max2, lb_min2, npgfb2, zetb2, &
513 3 : ra, rb, saabb, dmax)
514 :
515 : INTEGER, INTENT(IN) :: la_max1, la_min1, npgfa1
516 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta1
517 : INTEGER, INTENT(IN) :: la_max2, la_min2, npgfa2
518 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta2
519 : INTEGER, INTENT(IN) :: lb_max1, lb_min1, npgfb1
520 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb1
521 : INTEGER, INTENT(IN) :: lb_max2, lb_min2, npgfb2
522 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb2
523 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rb
524 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: saabb
525 : REAL(KIND=dp), INTENT(INOUT) :: dmax
526 :
527 : INTEGER :: coa1, coa2, cob1, cob2, i, iax, iay, &
528 : iaz, ibx, iby, ibz, ipgf, j, jpgf, k, &
529 : kpgf, l, la_max, la_min, lb_max, &
530 : lb_min, lpgf, ma, mb
531 : INTEGER, DIMENSION(3) :: na, naa, nb, nbb
532 : REAL(KIND=dp) :: res1, xa, xb
533 : REAL(KIND=dp), DIMENSION(3) :: A, B
534 :
535 3 : coa1 = 0
536 24 : DO ipgf = 1, npgfa1
537 21 : coa2 = 0
538 168 : DO jpgf = 1, npgfa2
539 147 : cob1 = 0
540 1176 : DO kpgf = 1, npgfb1
541 1029 : cob2 = 0
542 8232 : DO lpgf = 1, npgfb2
543 :
544 7203 : xa = zeta1(ipgf) + zeta2(jpgf) ! exponents
545 7203 : xb = zetb1(kpgf) + zetb2(lpgf) ! exponents
546 7203 : la_max = la_max1 + la_max2
547 7203 : lb_max = lb_max1 + lb_max2
548 7203 : la_min = la_min1 + la_min2
549 7203 : lb_min = lb_min1 + lb_min2
550 :
551 7203 : A = ra !positions
552 7203 : B = rb
553 :
554 7203 : CALL init_os_overlap2(xa, xb, A, B)
555 :
556 28812 : DO ma = la_min, la_max
557 93639 : DO mb = lb_min, lb_max
558 216090 : DO iax = 0, ma
559 410571 : DO iay = 0, ma - iax
560 216090 : iaz = ma - iax - iay
561 216090 : na(1) = iax; na(2) = iay; na(3) = iaz
562 777924 : DO ibx = 0, mb
563 1368570 : DO iby = 0, mb - ibx
564 720300 : ibz = mb - ibx - iby
565 720300 : nb(1) = ibx; nb(2) = iby; nb(3) = ibz
566 720300 : res1 = os_overlap2(na, nb)
567 4033680 : DO i = ncoset(la_min1 - 1) + 1, ncoset(la_max1)
568 15126300 : DO j = ncoset(la_min2 - 1) + 1, ncoset(la_max2)
569 46099200 : naa = indco(1:3, i) + indco(1:3, j)
570 60505200 : DO k = ncoset(lb_min1 - 1) + 1, ncoset(lb_max1)
571 242020800 : DO l = ncoset(lb_min2 - 1) + 1, ncoset(lb_max2)
572 737587200 : nbb = indco(1:3, k) + indco(1:3, l)
573 693792960 : IF (ALL(na == naa) .AND. ALL(nb == nbb)) THEN
574 1843968 : dmax = MAX(dmax, ABS(res1 - saabb(coa1 + i, coa2 + j, cob1 + k, cob2 + l)))
575 : END IF
576 : END DO
577 : END DO
578 : END DO
579 : END DO
580 : END DO
581 : END DO
582 : END DO
583 : END DO
584 : END DO
585 : END DO
586 8232 : cob2 = cob2 + ncoset(lb_max2)
587 : END DO
588 1176 : cob1 = cob1 + ncoset(lb_max1)
589 : END DO
590 168 : coa2 = coa2 + ncoset(la_max2)
591 : END DO
592 24 : coa1 = coa1 + ncoset(la_max1)
593 : END DO
594 :
595 : !WRITE(*,*) "dmax aabb", dmax
596 :
597 3 : END SUBROUTINE overlap_aabb_test
598 :
599 : END MODULE debug_os_integrals
|
