(git:1f285aa)
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 ! **************************************************************************************************
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  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 
470 END 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