Line data Source code
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 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 37050 : SUBROUTINE dabdr(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, rpgfb, lb_min, &
81 74100 : 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 37050 : na = 0
100 37050 : nda = 0
101 :
102 4814739 : dabdx = 0.0_dp
103 4814739 : dabdy = 0.0_dp
104 4814739 : dabdz = 0.0_dp
105 :
106 118074 : DO ipgf = 1, npgfa
107 :
108 81024 : fa = 2.0_dp*zeta(ipgf)
109 :
110 81024 : nb = 0
111 81024 : ndb = 0
112 :
113 291730 : DO jpgf = 1, npgfb
114 :
115 : ! *** Screening ***
116 :
117 210706 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
118 74097 : DO j = nb + ncoset(lb_min - 1) + 1, nb + ncoset(lb_max)
119 226720 : DO i = na + ncoset(la_min - 1) + 1, na + ncoset(la_max)
120 152623 : dabdx(i, j) = 0.0_dp
121 152623 : dabdy(i, j) = 0.0_dp
122 202862 : dabdz(i, j) = 0.0_dp
123 : END DO
124 : END DO
125 23858 : nb = nb + ncoset(lb_max)
126 23858 : ndb = ndb + ncoset(lb_max + 1)
127 23858 : CYCLE
128 : END IF
129 :
130 : ! *** [da/dRi|O|b] = 2*zeta*[a+1i|O|b] - Ni(a)[a-1i|O|b] ***
131 :
132 455689 : DO la = 0, la_max !MAX(0,la_min),la_max
133 :
134 455689 : IF (la == 0) THEN
135 :
136 186848 : coa = na + 1
137 186848 : coapx = nda + 2
138 186848 : coapy = nda + 3
139 186848 : coapz = nda + 4
140 :
141 454852 : DO lb = 0, lb_max !lb_min,lb_max
142 811395 : DO bx = 0, lb
143 1077626 : DO by = 0, lb - bx
144 453079 : bz = lb - bx - by
145 453079 : cob = nb + coset(bx, by, bz)
146 453079 : codb = ndb + coset(bx, by, bz)
147 453079 : dabdx(coa, cob) = fa*ab(coapx, codb)
148 453079 : dabdy(coa, cob) = fa*ab(coapy, codb)
149 809622 : dabdz(coa, cob) = fa*ab(coapz, codb)
150 : END DO
151 : END DO
152 : END DO
153 :
154 : ELSE
155 :
156 253476 : DO ax = 0, la
157 522560 : DO ay = 0, la - ax
158 269084 : az = la - ax - ay
159 :
160 269084 : coa = na + coset(ax, ay, az)
161 269084 : coamx = nda + coset(MAX(0, ax - 1), ay, az)
162 269084 : coamy = nda + coset(ax, MAX(0, ay - 1), az)
163 269084 : coamz = nda + coset(ax, ay, MAX(0, az - 1))
164 269084 : coapx = nda + coset(ax + 1, ay, az)
165 269084 : coapy = nda + coset(ax, ay + 1, az)
166 269084 : coapz = nda + coset(ax, ay, az + 1)
167 :
168 269084 : fx = REAL(ax, dp)
169 269084 : fy = REAL(ay, dp)
170 269084 : fz = REAL(az, dp)
171 :
172 858591 : DO lb = 0, lb_max !lb_min,lb_max
173 1275857 : DO bx = 0, lb
174 1790450 : DO by = 0, lb - bx
175 783677 : bz = lb - bx - by
176 783677 : cob = nb + coset(bx, by, bz)
177 783677 : codb = ndb + coset(bx, by, bz)
178 783677 : dabdx(coa, cob) = fa*ab(coapx, codb) - fx*ab(coamx, codb)
179 783677 : dabdy(coa, cob) = fa*ab(coapy, codb) - fy*ab(coamy, codb)
180 1372426 : 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 186848 : nb = nb + ncoset(lb_max)
193 267872 : ndb = ndb + ncoset(lb_max + 1)
194 :
195 : END DO
196 :
197 81024 : na = na + ncoset(la_max)
198 118074 : nda = nda + ncoset(la_max + 1)
199 :
200 : END DO
201 :
202 37050 : 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 5387625 : SUBROUTINE adbdr(la_max, npgfa, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, &
230 10775250 : 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 5387625 : na = 0
247 5387625 : nda = 0
248 :
249 401013270 : adbdx = 0.0_dp
250 401013270 : adbdy = 0.0_dp
251 401013270 : adbdz = 0.0_dp
252 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
253 17159812 : DO ipgf = 1, npgfa
254 :
255 11772187 : nb = 0
256 11772187 : ndb = 0
257 :
258 43068334 : DO jpgf = 1, npgfb
259 :
260 31296147 : fb = 2.0_dp*zetb(jpgf)
261 :
262 : ! *** Screening ***
263 :
264 31296147 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
265 81028821 : DO j = nb + ncoset(lb_min - 1) + 1, nb + ncoset(lb_max)
266 284437047 : DO i = na + ncoset(la_min - 1) + 1, na + ncoset(la_max)
267 203408226 : adbdx(i, j) = 0.0_dp
268 203408226 : adbdy(i, j) = 0.0_dp
269 263301969 : adbdz(i, j) = 0.0_dp
270 : END DO
271 : END DO
272 21135078 : nb = nb + ncoset(lb_max)
273 21135078 : ndb = ndb + ncoset(lb_max + 1)
274 21135078 : CYCLE
275 : END IF
276 :
277 : ! *** [a|O|db/dRi] = 2*zetb*[a|O|b+1i] - Ni(b)[a|O|b-1i] ***
278 :
279 27047088 : DO lb = 0, lb_max
280 :
281 27047088 : IF (lb == 0) THEN
282 :
283 10161069 : cob = nb + 1
284 10161069 : cobpx = ndb + 2
285 10161069 : cobpy = ndb + 3
286 10161069 : cobpz = ndb + 4
287 :
288 27048006 : DO la = 0, la_max !la_min,la_max
289 51398562 : DO ax = 0, la
290 73789419 : DO ay = 0, la - ax
291 32551926 : az = la - ax - ay
292 32551926 : coa = na + coset(ax, ay, az)
293 32551926 : coda = nda + coset(ax, ay, az)
294 32551926 : adbdx(coa, cob) = fb*ab(coda, cobpx)
295 32551926 : adbdy(coa, cob) = fb*ab(coda, cobpy)
296 56902482 : adbdz(coa, cob) = fb*ab(coda, cobpz)
297 : END DO
298 : END DO
299 : END DO
300 : ELSE
301 :
302 20912499 : DO bx = 0, lb
303 43300296 : DO by = 0, lb - bx
304 22387797 : bz = lb - bx - by
305 :
306 22387797 : cob = nb + coset(bx, by, bz)
307 22387797 : cobmx = ndb + coset(MAX(0, bx - 1), by, bz)
308 22387797 : cobmy = ndb + coset(bx, MAX(0, by - 1), bz)
309 22387797 : cobmz = ndb + coset(bx, by, MAX(0, bz - 1))
310 22387797 : cobpx = ndb + coset(bx + 1, by, bz)
311 22387797 : cobpy = ndb + coset(bx, by + 1, bz)
312 22387797 : cobpz = ndb + coset(bx, by, bz + 1)
313 :
314 22387797 : fx = REAL(bx, dp)
315 22387797 : fy = REAL(by, dp)
316 22387797 : fz = REAL(bz, dp)
317 :
318 79697718 : DO la = 0, la_max !la_min,la_max
319 131685135 : DO ax = 0, la
320 200842917 : DO ay = 0, la - ax
321 91545579 : az = la - ax - ay
322 91545579 : coa = na + coset(ax, ay, az)
323 91545579 : coda = nda + coset(ax, ay, az)
324 91545579 : adbdx(coa, cob) = fb*ab(coda, cobpx) - fx*ab(coda, cobmx)
325 91545579 : adbdy(coa, cob) = fb*ab(coda, cobpy) - fy*ab(coda, cobmy)
326 157720545 : 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 10161069 : nb = nb + ncoset(lb_max)
339 21933256 : ndb = ndb + ncoset(lb_max + 1)
340 :
341 : END DO
342 :
343 11772187 : na = na + ncoset(la_max)
344 17159812 : nda = nda + ncoset(la_max + 1)
345 :
346 : END DO
347 :
348 5387625 : 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 :
|
