(git:374b731)
Loading...
Searching...
No Matches
ai_derivatives.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 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! **************************************************************************************************
34MODULE 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
48CONTAINS
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 SUBROUTINE dabdr(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, rpgfb, lb_min, &
81 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 na = 0
100 nda = 0
101
102 dabdx = 0.0_dp
103 dabdy = 0.0_dp
104 dabdz = 0.0_dp
105
106 DO ipgf = 1, npgfa
107
108 fa = 2.0_dp*zeta(ipgf)
109
110 nb = 0
111 ndb = 0
112
113 DO jpgf = 1, npgfb
114
115 ! *** Screening ***
116
117 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
118 DO j = nb + ncoset(lb_min - 1) + 1, nb + ncoset(lb_max)
119 DO i = na + ncoset(la_min - 1) + 1, na + ncoset(la_max)
120 dabdx(i, j) = 0.0_dp
121 dabdy(i, j) = 0.0_dp
122 dabdz(i, j) = 0.0_dp
123 END DO
124 END DO
125 nb = nb + ncoset(lb_max)
126 ndb = ndb + ncoset(lb_max + 1)
127 cycle
128 END IF
129
130 ! *** [da/dRi|O|b] = 2*zeta*[a+1i|O|b] - Ni(a)[a-1i|O|b] ***
131
132 DO la = 0, la_max !MAX(0,la_min),la_max
133
134 IF (la == 0) THEN
135
136 coa = na + 1
137 coapx = nda + 2
138 coapy = nda + 3
139 coapz = nda + 4
140
141 DO lb = 0, lb_max !lb_min,lb_max
142 DO bx = 0, lb
143 DO by = 0, lb - bx
144 bz = lb - bx - by
145 cob = nb + coset(bx, by, bz)
146 codb = ndb + coset(bx, by, bz)
147 dabdx(coa, cob) = fa*ab(coapx, codb)
148 dabdy(coa, cob) = fa*ab(coapy, codb)
149 dabdz(coa, cob) = fa*ab(coapz, codb)
150 END DO
151 END DO
152 END DO
153
154 ELSE
155
156 DO ax = 0, la
157 DO ay = 0, la - ax
158 az = la - ax - ay
159
160 coa = na + coset(ax, ay, az)
161 coamx = nda + coset(max(0, ax - 1), ay, az)
162 coamy = nda + coset(ax, max(0, ay - 1), az)
163 coamz = nda + coset(ax, ay, max(0, az - 1))
164 coapx = nda + coset(ax + 1, ay, az)
165 coapy = nda + coset(ax, ay + 1, az)
166 coapz = nda + coset(ax, ay, az + 1)
167
168 fx = real(ax, dp)
169 fy = real(ay, dp)
170 fz = real(az, dp)
171
172 DO lb = 0, lb_max !lb_min,lb_max
173 DO bx = 0, lb
174 DO by = 0, lb - bx
175 bz = lb - bx - by
176 cob = nb + coset(bx, by, bz)
177 codb = ndb + coset(bx, by, bz)
178 dabdx(coa, cob) = fa*ab(coapx, codb) - fx*ab(coamx, codb)
179 dabdy(coa, cob) = fa*ab(coapy, codb) - fy*ab(coamy, codb)
180 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 nb = nb + ncoset(lb_max)
193 ndb = ndb + ncoset(lb_max + 1)
194
195 END DO
196
197 na = na + ncoset(la_max)
198 nda = nda + ncoset(la_max + 1)
199
200 END DO
201
202 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 SUBROUTINE adbdr(la_max, npgfa, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, &
230 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 na = 0
247 nda = 0
248
249 adbdx = 0.0_dp
250 adbdy = 0.0_dp
251 adbdz = 0.0_dp
252 ! *** Loop over all pairs of primitive Gaussian-type functions ***
253 DO ipgf = 1, npgfa
254
255 nb = 0
256 ndb = 0
257
258 DO jpgf = 1, npgfb
259
260 fb = 2.0_dp*zetb(jpgf)
261
262 ! *** Screening ***
263
264 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
265 DO j = nb + ncoset(lb_min - 1) + 1, nb + ncoset(lb_max)
266 DO i = na + ncoset(la_min - 1) + 1, na + ncoset(la_max)
267 adbdx(i, j) = 0.0_dp
268 adbdy(i, j) = 0.0_dp
269 adbdz(i, j) = 0.0_dp
270 END DO
271 END DO
272 nb = nb + ncoset(lb_max)
273 ndb = ndb + ncoset(lb_max + 1)
274 cycle
275 END IF
276
277 ! *** [a|O|db/dRi] = 2*zetb*[a|O|b+1i] - Ni(b)[a|O|b-1i] ***
278
279 DO lb = 0, lb_max
280
281 IF (lb == 0) THEN
282
283 cob = nb + 1
284 cobpx = ndb + 2
285 cobpy = ndb + 3
286 cobpz = ndb + 4
287
288 DO la = 0, la_max !la_min,la_max
289 DO ax = 0, la
290 DO ay = 0, la - ax
291 az = la - ax - ay
292 coa = na + coset(ax, ay, az)
293 coda = nda + coset(ax, ay, az)
294 adbdx(coa, cob) = fb*ab(coda, cobpx)
295 adbdy(coa, cob) = fb*ab(coda, cobpy)
296 adbdz(coa, cob) = fb*ab(coda, cobpz)
297 END DO
298 END DO
299 END DO
300 ELSE
301
302 DO bx = 0, lb
303 DO by = 0, lb - bx
304 bz = lb - bx - by
305
306 cob = nb + coset(bx, by, bz)
307 cobmx = ndb + coset(max(0, bx - 1), by, bz)
308 cobmy = ndb + coset(bx, max(0, by - 1), bz)
309 cobmz = ndb + coset(bx, by, max(0, bz - 1))
310 cobpx = ndb + coset(bx + 1, by, bz)
311 cobpy = ndb + coset(bx, by + 1, bz)
312 cobpz = ndb + coset(bx, by, bz + 1)
313
314 fx = real(bx, dp)
315 fy = real(by, dp)
316 fz = real(bz, dp)
317
318 DO la = 0, la_max !la_min,la_max
319 DO ax = 0, la
320 DO ay = 0, la - ax
321 az = la - ax - ay
322 coa = na + coset(ax, ay, az)
323 coda = nda + coset(ax, ay, az)
324 adbdx(coa, cob) = fb*ab(coda, cobpx) - fx*ab(coda, cobmx)
325 adbdy(coa, cob) = fb*ab(coda, cobpy) - fy*ab(coda, cobmy)
326 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 nb = nb + ncoset(lb_max)
339 ndb = ndb + ncoset(lb_max + 1)
340
341 END DO
342
343 na = na + ncoset(la_max)
344 nda = nda + ncoset(la_max + 1)
345
346 END DO
347
348 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 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 na = 0
381 nda = 0
382
383 dabdx = 0.0_dp
384 dabdy = 0.0_dp
385 dabdz = 0.0_dp
386
387 DO ipgf = 1, npgfa
388
389 fa = 2.0_dp*zeta(ipgf)
390
391 nb = 0
392 ndb = 0
393
394 DO jpgf = 1, npgfb
395
396 !*** [da/dRi|O|b] = 2*zeta*[a+1i|O|b] - Ni(a)[a-1i|O|b] ***
397
398 DO la = 0, la_max !MAX(0,la_min),la_max
399
400 IF (la == 0) THEN
401
402 coa = na + 1
403 coapx = nda + 2
404 coapy = nda + 3
405 coapz = nda + 4
406
407 DO lb = 0, lb_max !lb_min,lb_max
408 DO bx = 0, lb
409 DO by = 0, lb - bx
410 bz = lb - bx - by
411 cob = nb + coset(bx, by, bz)
412 codb = ndb + coset(bx, by, bz)
413 dabdx(coa, cob) = fa*ab(coapx, codb)
414 dabdy(coa, cob) = fa*ab(coapy, codb)
415 dabdz(coa, cob) = fa*ab(coapz, codb)
416 END DO
417 END DO
418 END DO
419
420 ELSE
421
422 DO ax = 0, la
423 DO ay = 0, la - ax
424 az = la - ax - ay
425
426 coa = na + coset(ax, ay, az)
427 coamx = nda + coset(max(0, ax - 1), ay, az)
428 coamy = nda + coset(ax, max(0, ay - 1), az)
429 coamz = nda + coset(ax, ay, max(0, az - 1))
430 coapx = nda + coset(ax + 1, ay, az)
431 coapy = nda + coset(ax, ay + 1, az)
432 coapz = nda + coset(ax, ay, az + 1)
433
434 fx = real(ax, dp)
435 fy = real(ay, dp)
436 fz = real(az, dp)
437
438 DO lb = 0, lb_max !lb_min,lb_max
439 DO bx = 0, lb
440 DO by = 0, lb - bx
441 bz = lb - bx - by
442 cob = nb + coset(bx, by, bz)
443 codb = ndb + coset(bx, by, bz)
444 dabdx(coa, cob) = fa*ab(coapx, codb) - fx*ab(coamx, codb)
445 dabdy(coa, cob) = fa*ab(coapy, codb) - fy*ab(coamy, codb)
446 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 nb = nb + ncoset(lb_max)
459 ndb = ndb + ncoset(lb_max + 1)
460
461 END DO
462
463 na = na + ncoset(la_max)
464 nda = nda + ncoset(la_max + 1)
465
466 END DO
467
468 END SUBROUTINE dabdr_noscreen
469
470END MODULE ai_derivatives
471
ai_derivatives
Calculate the first derivative of an integral block.
Definition ai_derivatives.F:34
ai_derivatives::dabdr_noscreen
subroutine, public dabdr_noscreen(la_max, npgfa, zeta, lb_max, npgfb, ab, dabdx, dabdy, dabdz)
Calculate the first derivative of an integral block. This takes the derivative with respect to the at...
Definition ai_derivatives.F:366
ai_derivatives::adbdr
subroutine, public adbdr(la_max, npgfa, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, dab, ab, adbdx, adbdy, adbdz)
Calculate the first derivative of an integral block. This takes the derivative with respect to the at...
Definition ai_derivatives.F:231
ai_derivatives::dabdr
subroutine, public dabdr(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, rpgfb, lb_min, dab, ab, dabdx, dabdy, dabdz)
Calculate the first derivative of an integral block. This takes the derivative with respect to the at...
Definition ai_derivatives.F:82
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