Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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_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 :
38 : CONTAINS
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 47 : SUBROUTINE overlap_ab_test(la_max, la_min, npgfa, zeta, lb_max, lb_min, npgfb, zetb, &
56 47 : 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 47 : coa = 0
73 100 : DO ipgf = 1, npgfa
74 53 : cob = 0
75 148 : DO jpgf = 1, npgfb
76 95 : xa = zeta(ipgf) !exponents
77 95 : xb = zetb(jpgf)
78 95 : A = ra !positions
79 95 : B = rb
80 95 : CALL init_os_overlap2(xa, xb, A, B)
81 314 : DO ma = la_min, la_max
82 827 : DO mb = lb_min, lb_max
83 1664 : DO iax = 0, ma
84 2939 : DO iay = 0, ma - iax
85 1494 : iaz = ma - iax - iay
86 1494 : na(1) = iax; na(2) = iay; na(3) = iaz
87 1494 : ia1 = coset(iax, iay, iaz)
88 5229 : DO ibx = 0, mb
89 8904 : DO iby = 0, mb - ibx
90 4607 : ibz = mb - ibx - iby
91 4607 : nb(1) = ibx; nb(2) = iby; nb(3) = ibz
92 4607 : ib1 = coset(ibx, iby, ibz)
93 4607 : res1 = os_overlap2(na, nb)
94 4607 : res2 = sab(coa + ia1, cob + ib1)
95 7410 : 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 148 : cob = cob + ncoset(lb_max)
103 : END DO
104 100 : coa = coa + ncoset(la_max)
105 : END DO
106 : !WRITE(*,*) "dmax overlap_ab_test", dmax
107 :
108 47 : 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 14 : SUBROUTINE overlap_abc_test(la_max, npgfa, zeta, la_min, &
131 14 : lb_max, npgfb, zetb, lb_min, &
132 14 : lc_max, npgfc, zetc, lc_min, &
133 14 : 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 14 : coa = 0
154 112 : DO ipgf = 1, npgfa
155 98 : cob = 0
156 784 : DO jpgf = 1, npgfb
157 686 : coc = 0
158 1372 : DO kpgf = 1, npgfc
159 :
160 686 : xa = zeta(ipgf) ! exponents
161 686 : xb = zetb(jpgf)
162 686 : xc = zetc(kpgf)
163 :
164 686 : A = Ra !positions
165 686 : B = Rb
166 686 : C = Rc
167 :
168 686 : CALL init_os_overlap3(xa, xb, xc, A, B, C)
169 :
170 2058 : DO ma = la_min, la_max
171 5586 : DO mc = lc_min, lc_max
172 11956 : DO mb = lb_min, lb_max
173 21168 : DO iax = 0, ma
174 31752 : DO iay = 0, ma - iax
175 14112 : iaz = ma - iax - iay
176 14112 : na(1) = iax; na(2) = iay; na(3) = iaz
177 14112 : ia1 = coset(iax, iay, iaz)
178 52920 : DO icx = 0, mc
179 90944 : DO icy = 0, mc - icx
180 48608 : icz = mc - icx - icy
181 48608 : nc(1) = icx; nc(2) = icy; nc(3) = icz
182 48608 : ic1 = coset(icx, icy, icz)
183 149744 : DO ibx = 0, mb
184 218736 : DO iby = 0, mb - ibx
185 97216 : ibz = mb - ibx - iby
186 97216 : nb(1) = ibx; nb(2) = iby; nb(3) = ibz
187 97216 : ib1 = coset(ibx, iby, ibz)
188 97216 : res1 = os_overlap3(na, nc, nb)
189 97216 : res2 = sabc(coa + ia1, cob + ib1, coc + ic1)
190 170128 : 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 1372 : coc = coc + ncoset(lc_max)
202 : END DO
203 784 : cob = cob + ncoset(lb_max)
204 : END DO
205 112 : coa = coa + ncoset(la_max)
206 : END DO
207 : !WRITE(*,*) "dmax abc", dmax
208 :
209 14 : 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 3 : SUBROUTINE overlap_aabb_test(la_max1, la_min1, npgfa1, zeta1, &
235 3 : la_max2, la_min2, npgfa2, zeta2, &
236 3 : lb_max1, lb_min1, npgfb1, zetb1, &
237 3 : lb_max2, lb_min2, npgfb2, zetb2, &
238 3 : 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 3 : coa1 = 0
261 24 : DO ipgf = 1, npgfa1
262 21 : coa2 = 0
263 168 : DO jpgf = 1, npgfa2
264 147 : cob1 = 0
265 1176 : DO kpgf = 1, npgfb1
266 1029 : cob2 = 0
267 8232 : DO lpgf = 1, npgfb2
268 :
269 7203 : xa = zeta1(ipgf) + zeta2(jpgf) ! exponents
270 7203 : xb = zetb1(kpgf) + zetb2(lpgf) ! exponents
271 7203 : la_max = la_max1 + la_max2
272 7203 : lb_max = lb_max1 + lb_max2
273 7203 : la_min = la_min1 + la_min2
274 7203 : lb_min = lb_min1 + lb_min2
275 :
276 7203 : A = ra !positions
277 7203 : B = rb
278 :
279 7203 : CALL init_os_overlap2(xa, xb, A, B)
280 :
281 28812 : DO ma = la_min, la_max
282 93639 : DO mb = lb_min, lb_max
283 216090 : DO iax = 0, ma
284 410571 : DO iay = 0, ma - iax
285 216090 : iaz = ma - iax - iay
286 216090 : na(1) = iax; na(2) = iay; na(3) = iaz
287 777924 : DO ibx = 0, mb
288 1368570 : DO iby = 0, mb - ibx
289 720300 : ibz = mb - ibx - iby
290 720300 : nb(1) = ibx; nb(2) = iby; nb(3) = ibz
291 720300 : res1 = os_overlap2(na, nb)
292 4033680 : DO i = ncoset(la_min1 - 1) + 1, ncoset(la_max1)
293 15126300 : DO j = ncoset(la_min2 - 1) + 1, ncoset(la_max2)
294 46099200 : naa = indco(1:3, i) + indco(1:3, j)
295 60505200 : DO k = ncoset(lb_min1 - 1) + 1, ncoset(lb_max1)
296 242020800 : DO l = ncoset(lb_min2 - 1) + 1, ncoset(lb_max2)
297 737587200 : nbb = indco(1:3, k) + indco(1:3, l)
298 693792960 : IF (ALL(na == naa) .AND. ALL(nb == nbb)) THEN
299 1843968 : 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 8232 : cob2 = cob2 + ncoset(lb_max2)
312 : END DO
313 1176 : cob1 = cob1 + ncoset(lb_max1)
314 : END DO
315 168 : coa2 = coa2 + ncoset(la_max2)
316 : END DO
317 24 : coa1 = coa1 + ncoset(la_max1)
318 : END DO
319 :
320 : !WRITE(*,*) "dmax aabb", dmax
321 :
322 3 : END SUBROUTINE overlap_aabb_test
323 :
324 : END MODULE debug_os_integrals
|
