(git:34ef472)
debug_os_integrals.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief 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  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  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  xa = 0.783300000000_dp ! exponents
59  xb = 1.239648746700_dp
60 
61  a = (/0.329309000000_dp, 0.28408240000_dp, 0.28408240000_dp/) !* bohr !positions
62  b = (/0.983983000000_dp, 0.00453720000_dp, 0.00432740000_dp/) !* bohr
63 
64  la_min = 0
65  lb_min = 0
66 
67  la_max = 3
68  lb_max = 4
69 
70  !---------------------------------------
71  ALLOCATE (sab(ncoset(la_max), ncoset(lb_max)))
72 
73  maxl = max(la_max, lb_max)
74  lds = ncoset(maxl)
75  ALLOCATE (swork(lds, lds, 1))
76  sab = 0._dp
77  rab(:) = b(:) - a(:)
78  dab = sqrt(dot_product(rab, rab))
79  xa_work(1) = xa
80  xb_work(1) = xb
81  rpgfa = 20._dp
82  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  rab=rab, dab=dab, sab=sab, da_max_set=0, return_derivatives=.false., s=swork, lds=lds)
86  !---------------------------------------
87 
88  CALL init_os_overlap2(xa, xb, a, b)
89 
90  dmax = 0._dp
91  DO ma = la_min, la_max
92  DO mb = lb_min, lb_max
93  DO iax = 0, ma
94  DO iay = 0, ma - iax
95  iaz = ma - iax - iay
96  na(1) = iax; na(2) = iay; na(3) = iaz
97  ia1 = coset(iax, iay, iaz)
98  DO ibx = 0, mb
99  DO iby = 0, mb - ibx
100  ibz = mb - ibx - iby
101  nb(1) = ibx; nb(2) = iby; nb(3) = ibz
102  ib1 = coset(ibx, iby, ibz)
103  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  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  DEALLOCATE (sab, swork)
115 
116  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  SUBROUTINE overlap_ab_test(la_max, la_min, npgfa, zeta, lb_max, lb_min, npgfb, zetb, &
134  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  coa = 0
151  DO ipgf = 1, npgfa
152  cob = 0
153  DO jpgf = 1, npgfb
154  xa = zeta(ipgf) !exponents
155  xb = zetb(jpgf)
156  a = ra !positions
157  b = rb
158  CALL init_os_overlap2(xa, xb, a, b)
159  DO ma = la_min, la_max
160  DO mb = lb_min, lb_max
161  DO iax = 0, ma
162  DO iay = 0, ma - iax
163  iaz = ma - iax - iay
164  na(1) = iax; na(2) = iay; na(3) = iaz
165  ia1 = coset(iax, iay, iaz)
166  DO ibx = 0, mb
167  DO iby = 0, mb - ibx
168  ibz = mb - ibx - iby
169  nb(1) = ibx; nb(2) = iby; nb(3) = ibz
170  ib1 = coset(ibx, iby, ibz)
171  res1 = os_overlap2(na, nb)
172  res2 = sab(coa + ia1, cob + ib1)
173  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  cob = cob + ncoset(lb_max)
181  END DO
182  coa = coa + ncoset(la_max)
183  END DO
184  !WRITE(*,*) "dmax overlap_ab_test", dmax
185 
186  END SUBROUTINE overlap_ab_test
187 
188 ! ***************************************************************************************************
189 !> \brief recursive test routines for integral (a,b,c) for only three exponents
190 ! **************************************************************************************************
191  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  xa = 0.783300000000_dp ! exponents
205  xb = 1.239648746700_dp
206  xc = 0.548370000000_dp
207 
208  a = (/0.329309000000_dp, 0.28408240000_dp, 0.28408240000_dp/) !* bohr !positions
209  b = (/0.983983000000_dp, 0.00453720000_dp, 0.00432740000_dp/) !* bohr
210  c = (/0.032380000000_dp, 1.23470000000_dp, 0.11137400000_dp/) !* bohr
211 
212  la_min = 0
213  lb_min = 0
214  lc_min = 0
215 
216  la_max = 0
217  lb_max = 0
218  lc_max = 1
219 
220  !---------------------------------------
221  rab(:) = b(:) - a(:)
222  dab = sqrt(dot_product(rab, rab))
223  rac(:) = c(:) - a(:)
224  dac = sqrt(dot_product(rac, rac))
225  rbc(:) = c(:) - b(:)
226  dbc = sqrt(dot_product(rbc, rbc))
227  ALLOCATE (sabc(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
228  xa_work(1) = xa
229  xb_work(1) = xb
230  xc_work(1) = xc
231  rpgfa = 20._dp
232  rpgfb = 20._dp
233  rpgfc = 20._dp
234  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  rab=rab, dab=dab, rac=rac, dac=dac, rbc=rbc, dbc=dbc, sabc=sabc)
240 
241  !---------------------------------------
242 
243  CALL init_os_overlap3(xa, xb, xc, a, b, c)
244 
245  dmax = 0._dp
246  DO ma = la_min, la_max
247  DO mc = lc_min, lc_max
248  DO mb = lb_min, lb_max
249  DO iax = 0, ma
250  DO iay = 0, ma - iax
251  iaz = ma - iax - iay
252  na(1) = iax; na(2) = iay; na(3) = iaz
253  ia1 = coset(iax, iay, iaz)
254  DO icx = 0, mc
255  DO icy = 0, mc - icx
256  icz = mc - icx - icy
257  nc(1) = icx; nc(2) = icy; nc(3) = icz
258  ic1 = coset(icx, icy, icz)
259  DO ibx = 0, mb
260  DO iby = 0, mb - ibx
261  ibz = mb - ibx - iby
262  nb(1) = ibx; nb(2) = iby; nb(3) = ibz
263  ib1 = coset(ibx, iby, ibz)
264  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  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  DEALLOCATE (sabc)
279 
280  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  SUBROUTINE overlap_abc_test(la_max, npgfa, zeta, la_min, &
303  lb_max, npgfb, zetb, lb_min, &
304  lc_max, npgfc, zetc, lc_min, &
305  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  coa = 0
326  DO ipgf = 1, npgfa
327  cob = 0
328  DO jpgf = 1, npgfb
329  coc = 0
330  DO kpgf = 1, npgfc
331 
332  xa = zeta(ipgf) ! exponents
333  xb = zetb(jpgf)
334  xc = zetc(kpgf)
335 
336  a = ra !positions
337  b = rb
338  c = rc
339 
340  CALL init_os_overlap3(xa, xb, xc, a, b, c)
341 
342  DO ma = la_min, la_max
343  DO mc = lc_min, lc_max
344  DO mb = lb_min, lb_max
345  DO iax = 0, ma
346  DO iay = 0, ma - iax
347  iaz = ma - iax - iay
348  na(1) = iax; na(2) = iay; na(3) = iaz
349  ia1 = coset(iax, iay, iaz)
350  DO icx = 0, mc
351  DO icy = 0, mc - icx
352  icz = mc - icx - icy
353  nc(1) = icx; nc(2) = icy; nc(3) = icz
354  ic1 = coset(icx, icy, icz)
355  DO ibx = 0, mb
356  DO iby = 0, mb - ibx
357  ibz = mb - ibx - iby
358  nb(1) = ibx; nb(2) = iby; nb(3) = ibz
359  ib1 = coset(ibx, iby, ibz)
360  res1 = os_overlap3(na, nc, nb)
361  res2 = sabc(coa + ia1, cob + ib1, coc + ic1)
362  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  coc = coc + ncoset(lc_max)
374  END DO
375  cob = cob + ncoset(lb_max)
376  END DO
377  coa = coa + ncoset(la_max)
378  END DO
379  !WRITE(*,*) "dmax abc", dmax
380 
381  END SUBROUTINE overlap_abc_test
382 
383 ! ***************************************************************************************************
384 !> \brief recursive test routines for integral (aa,bb) for only four exponents
385 ! **************************************************************************************************
386  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  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  xa1 = 0.783300000000_dp ! exponents
400  xb1 = 1.239648746700_dp
401  xa2 = 0.3400000000_dp ! exponents
402  xb2 = 2.76_dp
403 
404  a = (/0.329309000000_dp, 0.28408240000_dp, 0.28408240000_dp/) !* bohr !positions
405  b = (/0.983983000000_dp, 0.00453720000_dp, 0.00432740000_dp/) !* bohr
406 
407  la_min1 = 0
408  la_min2 = 0
409  lb_min1 = 3
410  lb_min2 = 1
411 
412  la_max1 = 1
413  la_max2 = 2
414  lb_max1 = 3
415  lb_max2 = 4
416 
417  !---------------------------------------
418  ALLOCATE (saabb(ncoset(la_max1), ncoset(la_max2), ncoset(lb_max1), ncoset(lb_max2)))
419 
420  maxl = max(la_max1 + la_max2, lb_max1 + lb_max2)
421  lds = ncoset(maxl)
422  ALLOCATE (swork(lds, lds))
423  saabb = 0._dp
424  rab(:) = b(:) - a(:)
425  dab = sqrt(dot_product(rab, rab))
426  xa_work1(1) = xa1
427  xa_work2(1) = xa2
428  xb_work1(1) = xb1
429  xb_work2(1) = xb2
430  rpgfa1 = 20._dp
431  rpgfa2 = 20._dp
432  rpgfb1 = 20._dp
433  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  asets_equal=.false., bsets_equal=.false., rab=rab, dab=dab, saabb=saabb, s=swork, lds=lds)
439  !---------------------------------------
440 
441  xa = xa1 + xa2
442  xb = xb1 + xb2
443  la_min = la_min1 + la_min2
444  la_max = la_max1 + la_max2
445  lb_min = lb_min1 + lb_min2
446  lb_max = lb_max1 + lb_max2
447 
448  CALL init_os_overlap2(xa, xb, a, b)
449 
450  dmax = 0._dp
451  DO ma = la_min, la_max
452  DO mb = lb_min, lb_max
453  DO iax = 0, ma
454  DO iay = 0, ma - iax
455  iaz = ma - iax - iay
456  na(1) = iax; na(2) = iay; na(3) = iaz
457  DO ibx = 0, mb
458  DO iby = 0, mb - ibx
459  ibz = mb - ibx - iby
460  nb(1) = ibx; nb(2) = iby; nb(3) = ibz
461  res1 = os_overlap2(na, nb)
462  DO i = ncoset(la_min1 - 1) + 1, ncoset(la_max1)
463  DO j = ncoset(la_min2 - 1) + 1, ncoset(la_max2)
464  naa = indco(1:3, i) + indco(1:3, j)
465  DO k = ncoset(lb_min1 - 1) + 1, ncoset(lb_max1)
466  DO l = ncoset(lb_min2 - 1) + 1, ncoset(lb_max2)
467  nbb = indco(1:3, k) + indco(1:3, l)
468  IF (all(na == naa) .AND. all(nb == nbb)) THEN
469  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  DEALLOCATE (saabb, swork)
483 
484  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  SUBROUTINE overlap_aabb_test(la_max1, la_min1, npgfa1, zeta1, &
510  la_max2, la_min2, npgfa2, zeta2, &
511  lb_max1, lb_min1, npgfb1, zetb1, &
512  lb_max2, lb_min2, npgfb2, zetb2, &
513  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  coa1 = 0
536  DO ipgf = 1, npgfa1
537  coa2 = 0
538  DO jpgf = 1, npgfa2
539  cob1 = 0
540  DO kpgf = 1, npgfb1
541  cob2 = 0
542  DO lpgf = 1, npgfb2
543 
544  xa = zeta1(ipgf) + zeta2(jpgf) ! exponents
545  xb = zetb1(kpgf) + zetb2(lpgf) ! exponents
546  la_max = la_max1 + la_max2
547  lb_max = lb_max1 + lb_max2
548  la_min = la_min1 + la_min2
549  lb_min = lb_min1 + lb_min2
550 
551  a = ra !positions
552  b = rb
553 
554  CALL init_os_overlap2(xa, xb, a, b)
555 
556  DO ma = la_min, la_max
557  DO mb = lb_min, lb_max
558  DO iax = 0, ma
559  DO iay = 0, ma - iax
560  iaz = ma - iax - iay
561  na(1) = iax; na(2) = iay; na(3) = iaz
562  DO ibx = 0, mb
563  DO iby = 0, mb - ibx
564  ibz = mb - ibx - iby
565  nb(1) = ibx; nb(2) = iby; nb(3) = ibz
566  res1 = os_overlap2(na, nb)
567  DO i = ncoset(la_min1 - 1) + 1, ncoset(la_max1)
568  DO j = ncoset(la_min2 - 1) + 1, ncoset(la_max2)
569  naa = indco(1:3, i) + indco(1:3, j)
570  DO k = ncoset(lb_min1 - 1) + 1, ncoset(lb_max1)
571  DO l = ncoset(lb_min2 - 1) + 1, ncoset(lb_max2)
572  nbb = indco(1:3, k) + indco(1:3, l)
573  IF (all(na == naa) .AND. all(nb == nbb)) THEN
574  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  cob2 = cob2 + ncoset(lb_max2)
587  END DO
588  cob1 = cob1 + ncoset(lb_max1)
589  END DO
590  coa2 = coa2 + ncoset(la_max2)
591  END DO
592  coa1 = coa1 + ncoset(la_max1)
593  END DO
594 
595  !WRITE(*,*) "dmax aabb", dmax
596 
597  END SUBROUTINE overlap_aabb_test
598 
599 END MODULE debug_os_integrals
ai_overlap3_debug
Three-center integrals over Cartesian Gaussian-type functions.
Definition: ai_overlap3_debug.F:16
ai_overlap3_debug::init_os_overlap3
subroutine, public init_os_overlap3(ya, yb, yc, rA, rB, rC)
Calculation of three-center integrals over Cartesian Gaussian-type functions.
Definition: ai_overlap3_debug.F:49
ai_overlap3_debug::os_overlap3
recursive real(dp) function, public os_overlap3(an, cn, bn)
...
Definition: ai_overlap3_debug.F:84
ai_overlap3
Definition: ai_overlap3.F:24
ai_overlap3::overlap3
subroutine, public overlap3(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max_set, npgfb, zetb, rpgfb, lb_min_set, lc_max_set, npgfc, zetc, rpgfc, lc_min_set, rab, dab, rac, dac, rbc, dbc, sabc, sdabc, sabdc, int_abc_ext)
Calculation of three-center overlap integrals [a|b|c] over primitive Cartesian Gaussian functions.
Definition: ai_overlap3.F:107
ai_overlap_aabb
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition: ai_overlap_aabb.F:27
ai_overlap_aabb::overlap_aabb
subroutine, public overlap_aabb(la_max_set1, la_min_set1, npgfa1, rpgfa1, zeta1, la_max_set2, la_min_set2, npgfa2, rpgfa2, zeta2, lb_max_set1, lb_min_set1, npgfb1, rpgfb1, zetb1, lb_max_set2, lb_min_set2, npgfb2, rpgfb2, zetb2, asets_equal, bsets_equal, rab, dab, saabb, s, lds)
Purpose: Calculation of the two-center overlap integrals [aa|bb] over Cartesian Gaussian-type functio...
Definition: ai_overlap_aabb.F:85
ai_overlap_debug
Two-center overlap integrals over Cartesian Gaussian-type functions.
Definition: ai_overlap_debug.F:16
ai_overlap_debug::os_overlap2
recursive real(dp) function, public os_overlap2(an, bn)
...
Definition: ai_overlap_debug.F:73
ai_overlap_debug::init_os_overlap2
subroutine, public init_os_overlap2(ya, yb, rA, rB)
Calculation of overlap integrals over Cartesian Gaussian-type functions.
Definition: ai_overlap_debug.F:47
ai_overlap
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition: ai_overlap.F:18
ai_overlap::overlap
subroutine, public overlap(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, rab, dab, sab, da_max_set, return_derivatives, s, lds, sdab, pab, force_a)
Purpose: Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions...
Definition: ai_overlap.F:73
debug_os_integrals
Debugs Obara-Saika integral matrices.
Definition: debug_os_integrals.F:14
debug_os_integrals::overlap_ab_test
subroutine, public overlap_ab_test(la_max, la_min, npgfa, zeta, lb_max, lb_min, npgfb, zetb, ra, rb, sab, dmax)
recursive test routines for integral (a,b)
Definition: debug_os_integrals.F:135
debug_os_integrals::overlap_abc_test
subroutine, public overlap_abc_test(la_max, npgfa, zeta, la_min, lb_max, npgfb, zetb, lb_min, lc_max, npgfc, zetc, lc_min, ra, rb, rc, sabc, dmax)
recursive test routines for integral (a,b,c)
Definition: debug_os_integrals.F:306
debug_os_integrals::overlap_aabb_test
subroutine, public overlap_aabb_test(la_max1, la_min1, npgfa1, zeta1, la_max2, la_min2, npgfa2, zeta2, lb_max1, lb_min1, npgfb1, zetb1, lb_max2, lb_min2, npgfb2, zetb2, ra, rb, saabb, dmax)
recursive test routines for integral (aa,bb)
Definition: debug_os_integrals.F:514
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition: orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition: orbital_pointers.F:42
orbital_pointers::coset
integer, dimension(:, :, :), allocatable, public coset
Definition: orbital_pointers.F:45
orbital_pointers::indco
integer, dimension(:, :), allocatable, public indco
Definition: orbital_pointers.F:43
xc
Exchange and Correlation functional calculations.
Definition: xc.F:17