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 Calculate the first derivative of an integral block.
10 : !> \author Matthias Krack
11 : !> \date 05.10.2000
12 : !> \version 1.0
13 : !> \par Literature
14 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
15 : !> \par Parameters
16 : !> - ax,ay,az : Angular momentum index numbers of orbital a.
17 : !> - bx,by,bz : Angular momentum index numbers of orbital b.
18 : !> - coset : Cartesian orbital set pointer.
19 : !> - l{a,b} : Angular momentum quantum number of shell a or b.
20 : !> - l{a,b}_max: Maximum angular momentum quantum number of shell a or b.
21 : !> - l{a,b}_min: Minimum angular momentum quantum number of shell a or b.
22 : !> - ncoset : Number of orbitals in a Cartesian orbital set.
23 : !> - npgf{a,b} : Degree of contraction of shell a or b.
24 : !> - rab : Distance vector between the atomic centers a and b.
25 : !> - rab2 : Square of the distance between the atomic centers a and b.
26 : !> - rac : Distance vector between the atomic centers a and c.
27 : !> - rac2 : Square of the distance between the atomic centers a and c.
28 : !> - rbc : Distance vector between the atomic centers b and c.
29 : !> - rbc2 : Square of the distance between the atomic centers b and c.
30 : !> - rpgf{a,b} : Radius of the primitive Gaussian-type function a or b.
31 : !> - zet{a,b} : Exponents of the Gaussian-type functions a or b.
32 : !> - zetp : Reciprocal of the sum of the exponents of orbital a and b.
33 : ! **************************************************************************************************
34 : MODULE ai_derivatives
35 :
36 : USE kinds, ONLY: dp
37 : USE orbital_pointers, ONLY: coset,&
38 : ncoset
39 : #include "../base/base_uses.f90"
40 :
41 : IMPLICIT NONE
42 :
43 : PRIVATE
44 :
45 : ! *** Public subroutines ***
46 : PUBLIC :: adbdr, dabdr, dabdr_noscreen
47 :
48 : CONTAINS
49 :
50 : ! **************************************************************************************************
51 : !> \brief Calculate the first derivative of an integral block.
52 : !> This takes the derivative with respect to
53 : !> the atomic position Ra, i.e. the center of the primitive on the left.
54 : !> To get the derivative of the left primitive with respect to r (orbital
55 : !> coordinate), take the opposite sign.
56 : !> To get the derivative with respect to the center of the primitive on
57 : !> the right Rb, take the opposite sign.
58 : !> To get the derivative of the right primitive with respect to r,
59 : !> do not change the sign.
60 : !> [da/dRi|O|b] = 2*zeta*[a+1i|O|b] - Ni(a)[a-1i|O|b]
61 : !>
62 : !> \param la_max ...
63 : !> \param npgfa ...
64 : !> \param zeta ...
65 : !> \param rpgfa ...
66 : !> \param la_min ...
67 : !> \param lb_max ...
68 : !> \param npgfb ...
69 : !> \param rpgfb ...
70 : !> \param lb_min ...
71 : !> \param dab ...
72 : !> \param ab ...
73 : !> \param dabdx ...
74 : !> \param dabdy ...
75 : !> \param dabdz ...
76 : !> \date 05.10.2000
77 : !> \author Matthias Krack
78 : !> \version 1.0
79 : ! **************************************************************************************************
80 39170 : SUBROUTINE dabdr(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, rpgfb, lb_min, &
81 78340 : dab, ab, dabdx, dabdy, dabdz)
82 : INTEGER, INTENT(IN) :: la_max, npgfa
83 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
84 : INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
85 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb
86 : INTEGER, INTENT(IN) :: lb_min
87 : REAL(KIND=dp), INTENT(IN) :: dab
88 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ab
89 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: dabdx, dabdy, dabdz
90 :
91 : INTEGER :: ax, ay, az, bx, by, bz, coa, coamx, &
92 : coamy, coamz, coapx, coapy, coapz, &
93 : cob, codb, i, ipgf, j, jpgf, la, lb, &
94 : na, nb, nda, ndb
95 : REAL(KIND=dp) :: fa, fx, fy, fz
96 :
97 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
98 :
99 39170 : na = 0
100 39170 : nda = 0
101 :
102 5053962 : dabdx = 0.0_dp
103 5053962 : dabdy = 0.0_dp
104 5053962 : dabdz = 0.0_dp
105 :
106 126182 : DO ipgf = 1, npgfa
107 :
108 87012 : fa = 2.0_dp*zeta(ipgf)
109 :
110 87012 : nb = 0
111 87012 : ndb = 0
112 :
113 315828 : DO jpgf = 1, npgfb
114 :
115 : ! *** Screening ***
116 :
117 228816 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
118 75108 : DO j = nb + ncoset(lb_min - 1) + 1, nb + ncoset(lb_max)
119 228469 : DO i = na + ncoset(la_min - 1) + 1, na + ncoset(la_max)
120 153361 : dabdx(i, j) = 0.0_dp
121 153361 : dabdy(i, j) = 0.0_dp
122 204191 : dabdz(i, j) = 0.0_dp
123 : END DO
124 : END DO
125 24278 : nb = nb + ncoset(lb_max)
126 24278 : ndb = ndb + ncoset(lb_max + 1)
127 24278 : CYCLE
128 : END IF
129 :
130 : ! *** [da/dRi|O|b] = 2*zeta*[a+1i|O|b] - Ni(a)[a-1i|O|b] ***
131 :
132 499463 : DO la = 0, la_max !MAX(0,la_min),la_max
133 :
134 499463 : IF (la == 0) THEN
135 :
136 204538 : coa = na + 1
137 204538 : coapx = nda + 2
138 204538 : coapy = nda + 3
139 204538 : coapz = nda + 4
140 :
141 498702 : DO lb = 0, lb_max !lb_min,lb_max
142 890811 : DO bx = 0, lb
143 1185274 : DO by = 0, lb - bx
144 499001 : bz = lb - bx - by
145 499001 : cob = nb + coset(bx, by, bz)
146 499001 : codb = ndb + coset(bx, by, bz)
147 499001 : dabdx(coa, cob) = fa*ab(coapx, codb)
148 499001 : dabdy(coa, cob) = fa*ab(coapy, codb)
149 891110 : dabdz(coa, cob) = fa*ab(coapz, codb)
150 : END DO
151 : END DO
152 : END DO
153 :
154 : ELSE
155 :
156 279634 : DO ax = 0, la
157 576842 : DO ay = 0, la - ax
158 297208 : az = la - ax - ay
159 :
160 297208 : coa = na + coset(ax, ay, az)
161 297208 : coamx = nda + coset(MAX(0, ax - 1), ay, az)
162 297208 : coamy = nda + coset(ax, MAX(0, ay - 1), az)
163 297208 : coamz = nda + coset(ax, ay, MAX(0, az - 1))
164 297208 : coapx = nda + coset(ax + 1, ay, az)
165 297208 : coapy = nda + coset(ax, ay + 1, az)
166 297208 : coapz = nda + coset(ax, ay, az + 1)
167 :
168 297208 : fx = REAL(ax, dp)
169 297208 : fy = REAL(ay, dp)
170 297208 : fz = REAL(az, dp)
171 :
172 949443 : DO lb = 0, lb_max !lb_min,lb_max
173 1413223 : DO bx = 0, lb
174 1985804 : DO by = 0, lb - bx
175 869789 : bz = lb - bx - by
176 869789 : cob = nb + coset(bx, by, bz)
177 869789 : codb = ndb + coset(bx, by, bz)
178 869789 : dabdx(coa, cob) = fa*ab(coapx, codb) - fx*ab(coamx, codb)
179 869789 : dabdy(coa, cob) = fa*ab(coapy, codb) - fy*ab(coamy, codb)
180 1522816 : dabdz(coa, cob) = fa*ab(coapz, codb) - fz*ab(coamz, codb)
181 : END DO
182 : END DO
183 : END DO
184 :
185 : END DO
186 : END DO
187 :
188 : END IF
189 :
190 : END DO
191 :
192 204538 : nb = nb + ncoset(lb_max)
193 291550 : ndb = ndb + ncoset(lb_max + 1)
194 :
195 : END DO
196 :
197 87012 : na = na + ncoset(la_max)
198 126182 : nda = nda + ncoset(la_max + 1)
199 :
200 : END DO
201 :
202 39170 : END SUBROUTINE dabdr
203 :
204 : ! **************************************************************************************************
205 : !> \brief Calculate the first derivative of an integral block.
206 : !> This takes the derivative with respect to
207 : !> the atomic position Rb, i.e. the center of the primitive on the right.
208 : !> To get the derivative of the left primitive with respect to r
209 : !> (orbital coordinate), take the opposite sign.
210 : !> [a|O|db/dRi] = 2*zetb*[a|O|b+1i] - Ni(b)[a|O|b-1i]
211 : !> \param la_max ...
212 : !> \param npgfa ...
213 : !> \param rpgfa ...
214 : !> \param la_min ...
215 : !> \param lb_max ...
216 : !> \param npgfb ...
217 : !> \param zetb ...
218 : !> \param rpgfb ...
219 : !> \param lb_min ...
220 : !> \param dab ...
221 : !> \param ab ...
222 : !> \param adbdx ...
223 : !> \param adbdy ...
224 : !> \param adbdz ...
225 : !> \date 05.10.2000
226 : !> \author Matthias Krack
227 : !> \version 1.0
228 : ! **************************************************************************************************
229 5389415 : SUBROUTINE adbdr(la_max, npgfa, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, &
230 10778830 : dab, ab, adbdx, adbdy, adbdz)
231 : INTEGER, INTENT(IN) :: la_max, npgfa
232 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa
233 : INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
234 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
235 : INTEGER, INTENT(IN) :: lb_min
236 : REAL(KIND=dp), INTENT(IN) :: dab
237 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ab
238 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: adbdx, adbdy, adbdz
239 :
240 : INTEGER :: ax, ay, az, bx, by, bz, coa, cob, cobmx, &
241 : cobmy, cobmz, cobpx, cobpy, cobpz, &
242 : coda, i, ipgf, j, jpgf, la, lb, na, &
243 : nb, nda, ndb
244 : REAL(KIND=dp) :: fb, fx, fy, fz
245 :
246 5389415 : na = 0
247 5389415 : nda = 0
248 :
249 401151339 : adbdx = 0.0_dp
250 401151339 : adbdy = 0.0_dp
251 401151339 : adbdz = 0.0_dp
252 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
253 17166893 : DO ipgf = 1, npgfa
254 :
255 11777478 : nb = 0
256 11777478 : ndb = 0
257 :
258 43090109 : DO jpgf = 1, npgfb
259 :
260 31312631 : fb = 2.0_dp*zetb(jpgf)
261 :
262 : ! *** Screening ***
263 :
264 31312631 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
265 81029553 : DO j = nb + ncoset(lb_min - 1) + 1, nb + ncoset(lb_max)
266 284438334 : DO i = na + ncoset(la_min - 1) + 1, na + ncoset(la_max)
267 203408781 : adbdx(i, j) = 0.0_dp
268 203408781 : adbdy(i, j) = 0.0_dp
269 263302962 : adbdz(i, j) = 0.0_dp
270 : END DO
271 : END DO
272 21135372 : nb = nb + ncoset(lb_max)
273 21135372 : ndb = ndb + ncoset(lb_max + 1)
274 21135372 : CYCLE
275 : END IF
276 :
277 : ! *** [a|O|db/dRi] = 2*zetb*[a|O|b+1i] - Ni(b)[a|O|b-1i] ***
278 :
279 27087325 : DO lb = 0, lb_max
280 :
281 27087325 : IF (lb == 0) THEN
282 :
283 10177259 : cob = nb + 1
284 10177259 : cobpx = ndb + 2
285 10177259 : cobpy = ndb + 3
286 10177259 : cobpz = ndb + 4
287 :
288 27088168 : DO la = 0, la_max !la_min,la_max
289 51471322 : DO ax = 0, la
290 73888057 : DO ay = 0, la - ax
291 32593994 : az = la - ax - ay
292 32593994 : coa = na + coset(ax, ay, az)
293 32593994 : coda = nda + coset(ax, ay, az)
294 32593994 : adbdx(coa, cob) = fb*ab(coda, cobpx)
295 32593994 : adbdy(coa, cob) = fb*ab(coda, cobpy)
296 56977148 : adbdz(coa, cob) = fb*ab(coda, cobpz)
297 : END DO
298 : END DO
299 : END DO
300 : ELSE
301 :
302 20936874 : DO bx = 0, lb
303 43350654 : DO by = 0, lb - bx
304 22413780 : bz = lb - bx - by
305 :
306 22413780 : cob = nb + coset(bx, by, bz)
307 22413780 : cobmx = ndb + coset(MAX(0, bx - 1), by, bz)
308 22413780 : cobmy = ndb + coset(bx, MAX(0, by - 1), bz)
309 22413780 : cobmz = ndb + coset(bx, by, MAX(0, bz - 1))
310 22413780 : cobpx = ndb + coset(bx + 1, by, bz)
311 22413780 : cobpy = ndb + coset(bx, by + 1, bz)
312 22413780 : cobpz = ndb + coset(bx, by, bz + 1)
313 :
314 22413780 : fx = REAL(bx, dp)
315 22413780 : fy = REAL(by, dp)
316 22413780 : fz = REAL(bz, dp)
317 :
318 79781763 : DO la = 0, la_max !la_min,la_max
319 131812059 : DO ax = 0, la
320 201023400 : DO ay = 0, la - ax
321 91625121 : az = la - ax - ay
322 91625121 : coa = na + coset(ax, ay, az)
323 91625121 : coda = nda + coset(ax, ay, az)
324 91625121 : adbdx(coa, cob) = fb*ab(coda, cobpx) - fx*ab(coda, cobmx)
325 91625121 : adbdy(coa, cob) = fb*ab(coda, cobpy) - fy*ab(coda, cobmy)
326 157859484 : adbdz(coa, cob) = fb*ab(coda, cobpz) - fz*ab(coda, cobmz)
327 : END DO
328 : END DO
329 : END DO
330 :
331 : END DO
332 : END DO
333 :
334 : END IF
335 :
336 : END DO
337 :
338 10177259 : nb = nb + ncoset(lb_max)
339 21954737 : ndb = ndb + ncoset(lb_max + 1)
340 :
341 : END DO
342 :
343 11777478 : na = na + ncoset(la_max)
344 17166893 : nda = nda + ncoset(la_max + 1)
345 :
346 : END DO
347 :
348 5389415 : END SUBROUTINE adbdr
349 : ! **************************************************************************************************
350 : !> \brief Calculate the first derivative of an integral block.
351 : !> This takes the derivative with respect to
352 : !> the atomic position Ra, i.e. the center of the primitive on the left.
353 : !> Difference to routine dabdr: employs no (!!) screening, which is relevant when
354 : !> calculating the derivatives of Coulomb integrals
355 : !> \param la_max ...
356 : !> \param npgfa ...
357 : !> \param zeta ...
358 : !> \param lb_max ...
359 : !> \param npgfb ...
360 : !> \param ab ...
361 : !> \param dabdx ...
362 : !> \param dabdy ...
363 : !> \param dabdz ...
364 : ! **************************************************************************************************
365 24000 : SUBROUTINE dabdr_noscreen(la_max, npgfa, zeta, lb_max, npgfb, ab, dabdx, dabdy, dabdz)
366 : INTEGER, INTENT(IN) :: la_max, npgfa
367 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta
368 : INTEGER, INTENT(IN) :: lb_max, npgfb
369 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ab
370 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: dabdx, dabdy, dabdz
371 :
372 : INTEGER :: ax, ay, az, bx, by, bz, coa, coamx, &
373 : coamy, coamz, coapx, coapy, coapz, &
374 : cob, codb, ipgf, jpgf, la, lb, na, nb, &
375 : nda, ndb
376 : REAL(KIND=dp) :: fa, fx, fy, fz
377 :
378 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
379 :
380 24000 : na = 0
381 24000 : nda = 0
382 :
383 8456000 : dabdx = 0.0_dp
384 8456000 : dabdy = 0.0_dp
385 8456000 : dabdz = 0.0_dp
386 :
387 48000 : DO ipgf = 1, npgfa
388 :
389 24000 : fa = 2.0_dp*zeta(ipgf)
390 :
391 24000 : nb = 0
392 24000 : ndb = 0
393 :
394 48000 : DO jpgf = 1, npgfb
395 :
396 : !*** [da/dRi|O|b] = 2*zeta*[a+1i|O|b] - Ni(a)[a-1i|O|b] ***
397 :
398 103200 : DO la = 0, la_max !MAX(0,la_min),la_max
399 :
400 103200 : IF (la == 0) THEN
401 :
402 24000 : coa = na + 1
403 24000 : coapx = nda + 2
404 24000 : coapy = nda + 3
405 24000 : coapz = nda + 4
406 :
407 115200 : DO lb = 0, lb_max !lb_min,lb_max
408 361600 : DO bx = 0, lb
409 881600 : DO by = 0, lb - bx
410 544000 : bz = lb - bx - by
411 544000 : cob = nb + coset(bx, by, bz)
412 544000 : codb = ndb + coset(bx, by, bz)
413 544000 : dabdx(coa, cob) = fa*ab(coapx, codb)
414 544000 : dabdy(coa, cob) = fa*ab(coapy, codb)
415 790400 : dabdz(coa, cob) = fa*ab(coapz, codb)
416 : END DO
417 : END DO
418 : END DO
419 :
420 : ELSE
421 :
422 213600 : DO ax = 0, la
423 537600 : DO ay = 0, la - ax
424 324000 : az = la - ax - ay
425 :
426 324000 : coa = na + coset(ax, ay, az)
427 324000 : coamx = nda + coset(MAX(0, ax - 1), ay, az)
428 324000 : coamy = nda + coset(ax, MAX(0, ay - 1), az)
429 324000 : coamz = nda + coset(ax, ay, MAX(0, az - 1))
430 324000 : coapx = nda + coset(ax + 1, ay, az)
431 324000 : coapy = nda + coset(ax, ay + 1, az)
432 324000 : coapz = nda + coset(ax, ay, az + 1)
433 :
434 324000 : fx = REAL(ax, dp)
435 324000 : fy = REAL(ay, dp)
436 324000 : fz = REAL(az, dp)
437 :
438 1713600 : DO lb = 0, lb_max !lb_min,lb_max
439 4881600 : DO bx = 0, lb
440 11901600 : DO by = 0, lb - bx
441 7344000 : bz = lb - bx - by
442 7344000 : cob = nb + coset(bx, by, bz)
443 7344000 : codb = ndb + coset(bx, by, bz)
444 7344000 : dabdx(coa, cob) = fa*ab(coapx, codb) - fx*ab(coamx, codb)
445 7344000 : dabdy(coa, cob) = fa*ab(coapy, codb) - fy*ab(coamy, codb)
446 10670400 : dabdz(coa, cob) = fa*ab(coapz, codb) - fz*ab(coamz, codb)
447 : END DO
448 : END DO
449 : END DO
450 :
451 : END DO
452 : END DO
453 :
454 : END IF
455 :
456 : END DO
457 :
458 24000 : nb = nb + ncoset(lb_max)
459 48000 : ndb = ndb + ncoset(lb_max + 1)
460 :
461 : END DO
462 :
463 24000 : na = na + ncoset(la_max)
464 48000 : nda = nda + ncoset(la_max + 1)
465 :
466 : END DO
467 :
468 24000 : END SUBROUTINE dabdr_noscreen
469 :
470 : END MODULE ai_derivatives
471 :
|
