(git:f2099e5)
Loading...
Searching...
No Matches
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-2026 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! **************************************************************************************************
14MODULE debug_os_integrals
15
16 USE ai_overlap3_debug, ONLY: init_os_overlap3,&
17 os_overlap3
18 USE ai_overlap_debug, ONLY: init_os_overlap2,&
19 os_overlap2
20 USE kinds, ONLY: dp
21 USE orbital_pointers, ONLY: coset,&
22 indco,&
23 ncoset
24#include "./base/base_uses.f90"
25
26 IMPLICIT NONE
27
28 PRIVATE
29
30! **************************************************************************************************
31
32 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'debug_os_integrals'
33
34 PUBLIC :: overlap_ab_test, overlap_abc_test, overlap_aabb_test
35
36! **************************************************************************************************
37
38CONTAINS
39
40! ***************************************************************************************************
41!> \brief recursive test routines for integral (a,b)
42!> \param la_max ...
43!> \param la_min ...
44!> \param npgfa ...
45!> \param zeta ...
46!> \param lb_max ...
47!> \param lb_min ...
48!> \param npgfb ...
49!> \param zetb ...
50!> \param ra ...
51!> \param rb ...
52!> \param sab ...
53!> \param dmax ...
54! **************************************************************************************************
55 SUBROUTINE overlap_ab_test(la_max, la_min, npgfa, zeta, lb_max, lb_min, npgfb, zetb, &
56 ra, rb, sab, dmax)
57
58 INTEGER, INTENT(IN) :: la_max, la_min, npgfa
59 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta
60 INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
61 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb
62 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: ra, rb
63 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: sab
64 REAL(kind=dp), INTENT(INOUT) :: dmax
65
66 INTEGER :: coa, cob, ia1, iax, iay, iaz, ib1, ibx, &
67 iby, ibz, ipgf, jpgf, ma, mb
68 INTEGER, DIMENSION(3) :: na, nb
69 REAL(kind=dp) :: res1, res2, xa, xb
70 REAL(kind=dp), DIMENSION(3) :: a, b
71
72 coa = 0
73 DO ipgf = 1, npgfa
74 cob = 0
75 DO jpgf = 1, npgfb
76 xa = zeta(ipgf) !exponents
77 xb = zetb(jpgf)
78 a = ra !positions
79 b = rb
80 CALL init_os_overlap2(xa, xb, a, b)
81 DO ma = la_min, la_max
82 DO mb = lb_min, lb_max
83 DO iax = 0, ma
84 DO iay = 0, ma - iax
85 iaz = ma - iax - iay
86 na(1) = iax; na(2) = iay; na(3) = iaz
87 ia1 = coset(iax, iay, iaz)
88 DO ibx = 0, mb
89 DO iby = 0, mb - ibx
90 ibz = mb - ibx - iby
91 nb(1) = ibx; nb(2) = iby; nb(3) = ibz
92 ib1 = coset(ibx, iby, ibz)
93 res1 = os_overlap2(na, nb)
94 res2 = sab(coa + ia1, cob + ib1)
95 dmax = max(dmax, abs(res1 - res2))
96 END DO
97 END DO
98 END DO
99 END DO
100 END DO
101 END DO
102 cob = cob + ncoset(lb_max)
103 END DO
104 coa = coa + ncoset(la_max)
105 END DO
106 !WRITE(*,*) "dmax overlap_ab_test", dmax
107
108 END SUBROUTINE overlap_ab_test
109
110! ***************************************************************************************************
111!> \brief recursive test routines for integral (a,b,c)
112!> \param la_max ...
113!> \param npgfa ...
114!> \param zeta ...
115!> \param la_min ...
116!> \param lb_max ...
117!> \param npgfb ...
118!> \param zetb ...
119!> \param lb_min ...
120!> \param lc_max ...
121!> \param npgfc ...
122!> \param zetc ...
123!> \param lc_min ...
124!> \param ra ...
125!> \param rb ...
126!> \param rc ...
127!> \param sabc ...
128!> \param dmax ...
129! **************************************************************************************************
130 SUBROUTINE overlap_abc_test(la_max, npgfa, zeta, la_min, &
131 lb_max, npgfb, zetb, lb_min, &
132 lc_max, npgfc, zetc, lc_min, &
133 ra, rb, rc, sabc, dmax)
134
135 INTEGER, INTENT(IN) :: la_max, npgfa
136 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta
137 INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
138 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb
139 INTEGER, INTENT(IN) :: lb_min, lc_max, npgfc
140 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetc
141 INTEGER, INTENT(IN) :: lc_min
142 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: ra, rb, rc
143 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: sabc
144 REAL(kind=dp), INTENT(INOUT) :: dmax
145
146 INTEGER :: coa, cob, coc, ia1, iax, iay, iaz, ib1, &
147 ibx, iby, ibz, ic1, icx, icy, icz, &
148 ipgf, jpgf, kpgf, ma, mb, mc
149 INTEGER, DIMENSION(3) :: na, nb, nc
150 REAL(kind=dp) :: res1, res2, xa, xb, xc
151 REAL(kind=dp), DIMENSION(3) :: a, b, c
152
153 coa = 0
154 DO ipgf = 1, npgfa
155 cob = 0
156 DO jpgf = 1, npgfb
157 coc = 0
158 DO kpgf = 1, npgfc
159
160 xa = zeta(ipgf) ! exponents
161 xb = zetb(jpgf)
162 xc = zetc(kpgf)
163
164 a = ra !positions
165 b = rb
166 c = rc
167
168 CALL init_os_overlap3(xa, xb, xc, a, b, c)
169
170 DO ma = la_min, la_max
171 DO mc = lc_min, lc_max
172 DO mb = lb_min, lb_max
173 DO iax = 0, ma
174 DO iay = 0, ma - iax
175 iaz = ma - iax - iay
176 na(1) = iax; na(2) = iay; na(3) = iaz
177 ia1 = coset(iax, iay, iaz)
178 DO icx = 0, mc
179 DO icy = 0, mc - icx
180 icz = mc - icx - icy
181 nc(1) = icx; nc(2) = icy; nc(3) = icz
182 ic1 = coset(icx, icy, icz)
183 DO ibx = 0, mb
184 DO iby = 0, mb - ibx
185 ibz = mb - ibx - iby
186 nb(1) = ibx; nb(2) = iby; nb(3) = ibz
187 ib1 = coset(ibx, iby, ibz)
188 res1 = os_overlap3(na, nc, nb)
189 res2 = sabc(coa + ia1, cob + ib1, coc + ic1)
190 dmax = max(dmax, abs(res1 - res2))
191 !IF(dmax > 1.E-10) WRITE(*,*) "dmax in loop", dmax
192 END DO
193 END DO
194 END DO
195 END DO
196 END DO
197 END DO
198 END DO
199 END DO
200 END DO
201 coc = coc + ncoset(lc_max)
202 END DO
203 cob = cob + ncoset(lb_max)
204 END DO
205 coa = coa + ncoset(la_max)
206 END DO
207 !WRITE(*,*) "dmax abc", dmax
208
209 END SUBROUTINE overlap_abc_test
210
211! ***************************************************************************************************
212!> \brief recursive test routines for integral (aa,bb)
213!> \param la_max1 ...
214!> \param la_min1 ...
215!> \param npgfa1 ...
216!> \param zeta1 ...
217!> \param la_max2 ...
218!> \param la_min2 ...
219!> \param npgfa2 ...
220!> \param zeta2 ...
221!> \param lb_max1 ...
222!> \param lb_min1 ...
223!> \param npgfb1 ...
224!> \param zetb1 ...
225!> \param lb_max2 ...
226!> \param lb_min2 ...
227!> \param npgfb2 ...
228!> \param zetb2 ...
229!> \param ra ...
230!> \param rb ...
231!> \param saabb ...
232!> \param dmax ...
233! **************************************************************************************************
234 SUBROUTINE overlap_aabb_test(la_max1, la_min1, npgfa1, zeta1, &
235 la_max2, la_min2, npgfa2, zeta2, &
236 lb_max1, lb_min1, npgfb1, zetb1, &
237 lb_max2, lb_min2, npgfb2, zetb2, &
238 ra, rb, saabb, dmax)
239
240 INTEGER, INTENT(IN) :: la_max1, la_min1, npgfa1
241 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta1
242 INTEGER, INTENT(IN) :: la_max2, la_min2, npgfa2
243 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta2
244 INTEGER, INTENT(IN) :: lb_max1, lb_min1, npgfb1
245 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb1
246 INTEGER, INTENT(IN) :: lb_max2, lb_min2, npgfb2
247 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb2
248 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: ra, rb
249 REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: saabb
250 REAL(kind=dp), INTENT(INOUT) :: dmax
251
252 INTEGER :: coa1, coa2, cob1, cob2, i, iax, iay, &
253 iaz, ibx, iby, ibz, ipgf, j, jpgf, k, &
254 kpgf, l, la_max, la_min, lb_max, &
255 lb_min, lpgf, ma, mb
256 INTEGER, DIMENSION(3) :: na, naa, nb, nbb
257 REAL(kind=dp) :: res1, xa, xb
258 REAL(kind=dp), DIMENSION(3) :: a, b
259
260 coa1 = 0
261 DO ipgf = 1, npgfa1
262 coa2 = 0
263 DO jpgf = 1, npgfa2
264 cob1 = 0
265 DO kpgf = 1, npgfb1
266 cob2 = 0
267 DO lpgf = 1, npgfb2
268
269 xa = zeta1(ipgf) + zeta2(jpgf) ! exponents
270 xb = zetb1(kpgf) + zetb2(lpgf) ! exponents
271 la_max = la_max1 + la_max2
272 lb_max = lb_max1 + lb_max2
273 la_min = la_min1 + la_min2
274 lb_min = lb_min1 + lb_min2
275
276 a = ra !positions
277 b = rb
278
279 CALL init_os_overlap2(xa, xb, a, b)
280
281 DO ma = la_min, la_max
282 DO mb = lb_min, lb_max
283 DO iax = 0, ma
284 DO iay = 0, ma - iax
285 iaz = ma - iax - iay
286 na(1) = iax; na(2) = iay; na(3) = iaz
287 DO ibx = 0, mb
288 DO iby = 0, mb - ibx
289 ibz = mb - ibx - iby
290 nb(1) = ibx; nb(2) = iby; nb(3) = ibz
291 res1 = os_overlap2(na, nb)
292 DO i = ncoset(la_min1 - 1) + 1, ncoset(la_max1)
293 DO j = ncoset(la_min2 - 1) + 1, ncoset(la_max2)
294 naa = indco(1:3, i) + indco(1:3, j)
295 DO k = ncoset(lb_min1 - 1) + 1, ncoset(lb_max1)
296 DO l = ncoset(lb_min2 - 1) + 1, ncoset(lb_max2)
297 nbb = indco(1:3, k) + indco(1:3, l)
298 IF (all(na == naa) .AND. all(nb == nbb)) THEN
299 dmax = max(dmax, abs(res1 - saabb(coa1 + i, coa2 + j, cob1 + k, cob2 + l)))
300 END IF
301 END DO
302 END DO
303 END DO
304 END DO
305 END DO
306 END DO
307 END DO
308 END DO
309 END DO
310 END DO
311 cob2 = cob2 + ncoset(lb_max2)
312 END DO
313 cob1 = cob1 + ncoset(lb_max1)
314 END DO
315 coa2 = coa2 + ncoset(la_max2)
316 END DO
317 coa1 = coa1 + ncoset(la_max1)
318 END DO
319
320 !WRITE(*,*) "dmax aabb", dmax
321
322 END SUBROUTINE overlap_aabb_test
323
324END MODULE debug_os_integrals
ai_overlap3_debug
Three-center integrals over Cartesian Gaussian-type functions.
Definition ai_overlap3_debug.F:16
ai_overlap3_debug::a
real(dp), dimension(3) a
Definition ai_overlap3_debug.F:29
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::c
real(dp), dimension(3) c
Definition ai_overlap3_debug.F:29
ai_overlap3_debug::b
real(dp), dimension(3) b
Definition ai_overlap3_debug.F:29
ai_overlap3_debug::xb
real(dp) xb
Definition ai_overlap3_debug.F:28
ai_overlap3_debug::os_overlap3
recursive real(dp) function, public os_overlap3(an, cn, bn)
...
Definition ai_overlap3_debug.F:84
ai_overlap3_debug::xa
real(dp) xa
Definition ai_overlap3_debug.F:28
ai_overlap3_debug::zeta
real(dp) zeta
Definition ai_overlap3_debug.F:31
ai_overlap_debug
Two-center overlap integrals over Cartesian Gaussian-type functions.
Definition ai_overlap_debug.F:16
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_debug::os_overlap2
recursive real(dp) function, public os_overlap2(an, bn)
...
Definition ai_overlap_debug.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:57
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:134
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:239
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