(git:58e3e09)
ai_moments.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 Calculation of the moment integrals over Cartesian Gaussian-type
10 !> functions.
11 !> \par Literature
12 !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
13 !> \par History
14 !> none
15 !> \author J. Hutter (16.02.2005)
16 ! **************************************************************************************************
17 MODULE ai_moments
18 
19 ! ax,ay,az : Angular momentum index numbers of orbital a.
20 ! bx,by,bz : Angular momentum index numbers of orbital b.
21 ! coset : Cartesian orbital set pointer.
22 ! dab : Distance between the atomic centers a and b.
23 ! l{a,b} : Angular momentum quantum number of shell a or b.
24 ! l{a,b}_max: Maximum angular momentum quantum number of shell a or b.
25 ! l{a,b}_min: Minimum angular momentum quantum number of shell a or b.
26 ! rac : Distance vector between the atomic center a and reference point c.
27 ! rbc : Distance vector between the atomic center b and reference point c.
28 ! rpgf{a,b} : Radius of the primitive Gaussian-type function a or b.
29 ! zet{a,b} : Exponents of the Gaussian-type functions a or b.
30 ! zetp : Reciprocal of the sum of the exponents of orbital a and b.
31 
32  USE ai_derivatives, ONLY: adbdr,&
33  dabdr
34  USE ai_overlap, ONLY: overlap
35  USE kinds, ONLY: dp
36  USE mathconstants, ONLY: pi
37  USE orbital_pointers, ONLY: coset,&
38  indco,&
39  ncoset
40 #include "../base/base_uses.f90"
41 
42  IMPLICIT NONE
43 
44  PRIVATE
45 
46  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_moments'
47 
50 
51 CONTAINS
52 
53 ! *****************************************************************************
54 !> \brief This returns the derivative of the moment integrals [a|\mu|b], with respect
55 !> to the primitive on the right
56 !> difmab(:, :, beta, alpha) = < a | r_beta | ∂_alpha b > * (iatom - jatom)
57 !> \param la_max ...
58 !> \param npgfa ...
59 !> \param zeta ...
60 !> \param rpgfa ...
61 !> \param la_min ...
62 !> \param lb_max ...
63 !> \param npgfb ...
64 !> \param zetb ...
65 !> \param rpgfb ...
66 !> \param lb_min ...
67 !> \param order ...
68 !> \param rac ...
69 !> \param rbc ...
70 !> \param difmab ...
71 !> \param lambda The atom on which we take the derivative
72 !> \param iatom ...
73 !> \param jatom ...
74 !> \author Edward Ditler
75 ! **************************************************************************************************
76  SUBROUTINE diff_momop_velocity(la_max, npgfa, zeta, rpgfa, la_min, &
77  lb_max, npgfb, zetb, rpgfb, lb_min, &
78  order, rac, rbc, difmab, lambda, iatom, jatom)
79 
80  INTEGER, INTENT(IN) :: la_max, npgfa
81  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
82  INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
83  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
84  INTEGER, INTENT(IN) :: lb_min, order
85  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
86  REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(OUT) :: difmab
87  INTEGER, INTENT(IN) :: lambda
88  INTEGER, INTENT(IN), OPTIONAL :: iatom, jatom
89 
90  INTEGER :: alpha, beta, lda, lda_min, ldb, ldb_min
91  REAL(kind=dp) :: dab, rab(3)
92  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: difmab_tmp, mab
93 
94  rab = rbc - rac
95  dab = sqrt(sum(rab**2))
96 
97  lda_min = max(0, la_min - 1)
98  ldb_min = max(0, lb_min - 1)
99  lda = ncoset(la_max)*npgfa
100  ldb = ncoset(lb_max)*npgfb
101  ALLOCATE (difmab_tmp(lda, ldb, 3))
102 
103  ALLOCATE (mab(npgfa*ncoset(la_max + 1), npgfb*ncoset(lb_max + 1), ncoset(order) - 1))
104  ! *** Calculate the primitive overlap integrals ***
105  ! mab(1:3) = < a | r_beta - RC_beta | b >
106  mab = 0.0_dp
107  CALL moment(la_max + 1, npgfa, zeta, rpgfa, lda_min, &
108  lb_max + 1, npgfb, zetb, rpgfb, &
109  order, rac, rbc, mab)
110 
111  difmab = 0.0_dp
112  DO beta = 1, ncoset(order) - 1 ! beta was imom
113 
114  difmab_tmp = 0.0_dp
115  CALL adbdr(la_max, npgfa, rpgfa, la_min, &
116  lb_max, npgfb, zetb, rpgfb, lb_min, &
117  dab, mab(:, :, beta), difmab_tmp(:, :, 1), &
118  difmab_tmp(:, :, 2), difmab_tmp(:, :, 3))
119 
120  ! difmab(beta, alpha) = < a | r_beta - RC_beta | ∂_alpha b > * [(a==lambda) - (b==lambda)]
121  DO alpha = 1, 3
122  IF (iatom == lambda) difmab(:, :, beta, alpha) = difmab(:, :, beta, alpha) + difmab_tmp(:, :, alpha)
123  IF (jatom == lambda) difmab(:, :, beta, alpha) = difmab(:, :, beta, alpha) - difmab_tmp(:, :, alpha)
124  END DO
125  END DO
126 
127  DEALLOCATE (mab)
128  DEALLOCATE (difmab_tmp)
129  END SUBROUTINE diff_momop_velocity
130 
131 ! *****************************************************************************
132 !> \brief This returns the derivative of the moment integrals [a|\mu|b], with respect
133 !> to the position of the primitive on the left and right, i.e.
134 !> [da/dR_ai|\mu|b] + [a|\mu|d/dR_bi]
135 !> [da/dR_ai|\mu|b] = 2*zeta*[a+1i|\mu|b] - Ni(a)[a-1i|\mu|b]
136 !> [a|\mu|d/dR_bi] = 2*zetb*[a|\mu|b+1i] - Ni(b)[a|\mu|b-1i]
137 !> order indicates the max order of the moment operator to be calculated
138 !> 1: dipole
139 !> 2: quadrupole
140 !> ...
141 !> \param la_max ...
142 !> \param npgfa ...
143 !> \param zeta ...
144 !> \param rpgfa ...
145 !> \param la_min ...
146 !> \param lb_max ...
147 !> \param npgfb ...
148 !> \param zetb ...
149 !> \param rpgfb ...
150 !> \param lb_min ...
151 !> \param order ...
152 !> \param rac ...
153 !> \param rbc ...
154 !> \param difmab ...
155 !> \param mab_ext ...
156 !> \param deltaR needed for weighted derivative
157 !> \param iatom ...
158 !> \param jatom ...
159 !> SL August 2015, ED 2021
160 ! **************************************************************************************************
161  SUBROUTINE diff_momop2(la_max, npgfa, zeta, rpgfa, la_min, &
162  lb_max, npgfb, zetb, rpgfb, lb_min, &
163  order, rac, rbc, difmab, mab_ext, deltaR, iatom, jatom)
164 
165  INTEGER, INTENT(IN) :: la_max, npgfa
166  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
167  INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
168  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
169  INTEGER, INTENT(IN) :: lb_min, order
170  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
171  REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(OUT) :: difmab
172  REAL(kind=dp), DIMENSION(:, :, :), OPTIONAL, &
173  POINTER :: mab_ext
174  REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
175  OPTIONAL, POINTER :: deltar
176  INTEGER, INTENT(IN), OPTIONAL :: iatom, jatom
177 
178  INTEGER :: imom, lda, lda_min, ldb, ldb_min
179  REAL(kind=dp) :: dab, rab(3)
180  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: difmab_tmp
181  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: mab
182 
183  rab = rbc - rac
184  dab = sqrt(sum(rab**2))
185 
186  lda_min = max(0, la_min - 1)
187  ldb_min = max(0, lb_min - 1)
188  lda = ncoset(la_max)*npgfa
189  ldb = ncoset(lb_max)*npgfb
190  ALLOCATE (difmab_tmp(lda, ldb, 3))
191 
192  IF (PRESENT(mab_ext)) THEN
193  mab => mab_ext
194  ELSE
195  ALLOCATE (mab(npgfa*ncoset(la_max + 1), npgfb*ncoset(lb_max + 1), &
196  ncoset(order) - 1))
197  mab = 0.0_dp
198 ! *** Calculate the primitive moment integrals ***
199  CALL moment(la_max + 1, npgfa, zeta, rpgfa, lda_min, &
200  lb_max + 1, npgfb, zetb, rpgfb, &
201  order, rac, rbc, mab)
202  END IF
203  DO imom = 1, ncoset(order) - 1
204  difmab(:, :, imom, :) = 0.0_dp
205 
206  difmab_tmp = 0.0_dp
207  CALL adbdr(la_max, npgfa, rpgfa, la_min, &
208  lb_max, npgfb, zetb, rpgfb, lb_min, &
209  dab, mab(:, :, imom), difmab_tmp(:, :, 1), &
210  difmab_tmp(:, :, 2), difmab_tmp(:, :, 3))
211 
212  difmab(:, :, imom, 1) = difmab_tmp(:, :, 1)*deltar(1, jatom)
213  difmab(:, :, imom, 2) = difmab_tmp(:, :, 2)*deltar(2, jatom)
214  difmab(:, :, imom, 3) = difmab_tmp(:, :, 3)*deltar(3, jatom)
215 
216  difmab_tmp = 0.0_dp
217  CALL dabdr(la_max, npgfa, zeta, rpgfa, la_min, &
218  lb_max, npgfb, rpgfb, lb_min, &
219  dab, mab(:, :, imom), difmab_tmp(:, :, 1), &
220  difmab_tmp(:, :, 2), difmab_tmp(:, :, 3))
221 
222  difmab(:, :, imom, 1) = difmab(:, :, imom, 1) + difmab_tmp(:, :, 1)*deltar(1, iatom)
223  difmab(:, :, imom, 2) = difmab(:, :, imom, 2) + difmab_tmp(:, :, 2)*deltar(2, iatom)
224  difmab(:, :, imom, 3) = difmab(:, :, imom, 3) + difmab_tmp(:, :, 3)*deltar(3, iatom)
225  END DO
226 
227  IF (PRESENT(mab_ext)) THEN
228  NULLIFY (mab)
229  ELSE
230  DEALLOCATE (mab)
231  END IF
232  DEALLOCATE (difmab_tmp)
233  END SUBROUTINE diff_momop2
234 
235 ! **************************************************************************************************
236 !> \brief ...
237 !> \param cos_block ...
238 !> \param sin_block ...
239 !> \param iatom ...
240 !> \param ncoa ...
241 !> \param nsgfa ...
242 !> \param sgfa ...
243 !> \param sphi_a ...
244 !> \param ldsa ...
245 !> \param jatom ...
246 !> \param ncob ...
247 !> \param nsgfb ...
248 !> \param sgfb ...
249 !> \param sphi_b ...
250 !> \param ldsb ...
251 !> \param cosab ...
252 !> \param sinab ...
253 !> \param ldab ...
254 !> \param work ...
255 !> \param ldwork ...
256 ! **************************************************************************************************
257  SUBROUTINE contract_cossin(cos_block, sin_block, &
258  iatom, ncoa, nsgfa, sgfa, sphi_a, ldsa, &
259  jatom, ncob, nsgfb, sgfb, sphi_b, ldsb, &
260  cosab, sinab, ldab, work, ldwork)
261 
262  REAL(dp), DIMENSION(:, :), POINTER :: cos_block, sin_block
263  INTEGER, INTENT(IN) :: iatom, ncoa, nsgfa, sgfa
264  REAL(dp), DIMENSION(:, :), INTENT(IN) :: sphi_a
265  INTEGER, INTENT(IN) :: ldsa, jatom, ncob, nsgfb, sgfb
266  REAL(dp), DIMENSION(:, :), INTENT(IN) :: sphi_b
267  INTEGER, INTENT(IN) :: ldsb
268  REAL(dp), DIMENSION(:, :), INTENT(IN) :: cosab, sinab
269  INTEGER, INTENT(IN) :: ldab
270  REAL(dp), DIMENSION(:, :) :: work
271  INTEGER, INTENT(IN) :: ldwork
272 
273 ! Calculate cosine
274 
275  CALL dgemm("N", "N", ncoa, nsgfb, ncob, &
276  1.0_dp, cosab(1, 1), ldab, &
277  sphi_b(1, sgfb), ldsb, &
278  0.0_dp, work(1, 1), ldwork)
279 
280  IF (iatom <= jatom) THEN
281  CALL dgemm("T", "N", nsgfa, nsgfb, ncoa, &
282  1.0_dp, sphi_a(1, sgfa), ldsa, &
283  work(1, 1), ldwork, &
284  1.0_dp, cos_block(sgfa, sgfb), &
285  SIZE(cos_block, 1))
286  ELSE
287  CALL dgemm("T", "N", nsgfb, nsgfa, ncoa, &
288  1.0_dp, work(1, 1), ldwork, &
289  sphi_a(1, sgfa), ldsa, &
290  1.0_dp, cos_block(sgfb, sgfa), &
291  SIZE(cos_block, 1))
292  END IF
293 
294  ! Calculate sine
295  CALL dgemm("N", "N", ncoa, nsgfb, ncob, &
296  1.0_dp, sinab(1, 1), ldab, &
297  sphi_b(1, sgfb), ldsb, &
298  0.0_dp, work(1, 1), ldwork)
299 
300  IF (iatom <= jatom) THEN
301  CALL dgemm("T", "N", nsgfa, nsgfb, ncoa, &
302  1.0_dp, sphi_a(1, sgfa), ldsa, &
303  work(1, 1), ldwork, &
304  1.0_dp, sin_block(sgfa, sgfb), &
305  SIZE(sin_block, 1))
306  ELSE
307  CALL dgemm("T", "N", nsgfb, nsgfa, ncoa, &
308  1.0_dp, work(1, 1), ldwork, &
309  sphi_a(1, sgfa), ldsa, &
310  1.0_dp, sin_block(sgfb, sgfa), &
311  SIZE(sin_block, 1))
312  END IF
313 
314  END SUBROUTINE contract_cossin
315 
316 ! **************************************************************************************************
317 !> \brief ...
318 !> \param la_max_set ...
319 !> \param npgfa ...
320 !> \param zeta ...
321 !> \param rpgfa ...
322 !> \param la_min_set ...
323 !> \param lb_max ...
324 !> \param npgfb ...
325 !> \param zetb ...
326 !> \param rpgfb ...
327 !> \param lb_min ...
328 !> \param rac ...
329 !> \param rbc ...
330 !> \param kvec ...
331 !> \param cosab ...
332 !> \param sinab ...
333 !> \param dcosab ...
334 !> \param dsinab ...
335 ! **************************************************************************************************
336  SUBROUTINE cossin(la_max_set, npgfa, zeta, rpgfa, la_min_set, &
337  lb_max, npgfb, zetb, rpgfb, lb_min, &
338  rac, rbc, kvec, cosab, sinab, dcosab, dsinab)
339 
340  INTEGER, INTENT(IN) :: la_max_set, npgfa
341  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
342  INTEGER, INTENT(IN) :: la_min_set, lb_max, npgfb
343  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
344  INTEGER, INTENT(IN) :: lb_min
345  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac, rbc, kvec
346  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: cosab, sinab
347  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
348  OPTIONAL :: dcosab, dsinab
349 
350  INTEGER :: ax, ay, az, bx, by, bz, cda, cdax, cday, cdaz, coa, coamx, coamy, coamz, coapx, &
351  coapy, coapz, cob, da, da_max, dax, day, daz, i, ipgf, j, jpgf, k, la, la_max, la_min, &
352  la_start, lb, lb_start, na, nb
353  REAL(kind=dp) :: dab, f0, f1, f2, f3, fax, fay, faz, ftz, &
354  fx, fy, fz, k2, kdp, rab2, s, zetp
355  REAL(kind=dp), DIMENSION(3) :: rab, rap, rbp
356  REAL(kind=dp), DIMENSION(ncoset(la_max_set), & ncoset(lb_max), 3) :: dscos, dssin
357  REAL(kind=dp), &
358  DIMENSION(ncoset(la_max_set+1), ncoset(lb_max)) :: sc, ss
359 
360  rab = rbc - rac
361  rab2 = sum(rab**2)
362  dab = sqrt(rab2)
363  k2 = kvec(1)*kvec(1) + kvec(2)*kvec(2) + kvec(3)*kvec(3)
364 
365  IF (PRESENT(dcosab)) THEN
366  da_max = 1
367  la_max = la_max_set + 1
368  la_min = max(0, la_min_set - 1)
369  dscos = 0.0_dp
370  dssin = 0.0_dp
371  ELSE
372  da_max = 0
373  la_max = la_max_set
374  la_min = la_min_set
375  END IF
376 
377  ! initialize all matrix elements to zero
378  IF (PRESENT(dcosab)) THEN
379  na = ncoset(la_max - 1)*npgfa
380  ELSE
381  na = ncoset(la_max)*npgfa
382  END IF
383  nb = ncoset(lb_max)*npgfb
384  cosab(1:na, 1:nb) = 0.0_dp
385  sinab(1:na, 1:nb) = 0.0_dp
386  IF (PRESENT(dcosab)) THEN
387  dcosab(1:na, 1:nb, :) = 0.0_dp
388  dsinab(1:na, 1:nb, :) = 0.0_dp
389  END IF
390 ! *** Loop over all pairs of primitive Gaussian-type functions ***
391 
392  na = 0
393  DO ipgf = 1, npgfa
394 
395  nb = 0
396 
397  DO jpgf = 1, npgfb
398 
399  ss = 0.0_dp
400  sc = 0.0_dp
401 
402 ! *** Screening ***
403  IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
404  nb = nb + ncoset(lb_max)
405  cycle
406  END IF
407 
408 ! *** Calculate some prefactors ***
409 
410  zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
411 
412  f0 = (pi*zetp)**1.5_dp
413  f1 = zetb(jpgf)*zetp
414  f2 = 0.5_dp*zetp
415 
416  kdp = zetp*dot_product(kvec, zeta(ipgf)*rac + zetb(jpgf)*rbc)
417 
418 ! *** Calculate the basic two-center cos/sin integral [s|cos/sin|s] ***
419 
420  s = f0*exp(-zeta(ipgf)*f1*rab2)*exp(-0.25_dp*k2*zetp)
421  sc(1, 1) = s*cos(kdp)
422  ss(1, 1) = s*sin(kdp)
423 
424 ! *** Recurrence steps: [s|O|s] -> [a|O|b] ***
425 
426  IF (la_max > 0) THEN
427 
428 ! *** Vertical recurrence steps: [s|O|s] -> [a|O|s] ***
429 
430  rap(:) = f1*rab(:)
431 
432 ! *** [p|O|s] = (Pi - Ai)*[s|O|s] +[s|dO|s] (i = x,y,z) ***
433 
434  sc(2, 1) = rap(1)*sc(1, 1) - f2*kvec(1)*ss(1, 1)
435  sc(3, 1) = rap(2)*sc(1, 1) - f2*kvec(2)*ss(1, 1)
436  sc(4, 1) = rap(3)*sc(1, 1) - f2*kvec(3)*ss(1, 1)
437  ss(2, 1) = rap(1)*ss(1, 1) + f2*kvec(1)*sc(1, 1)
438  ss(3, 1) = rap(2)*ss(1, 1) + f2*kvec(2)*sc(1, 1)
439  ss(4, 1) = rap(3)*ss(1, 1) + f2*kvec(3)*sc(1, 1)
440 
441 ! *** [a|O|s] = (Pi - Ai)*[a-1i|O|s] + f2*Ni(a-1i)*[a-2i|s] ***
442 ! *** + [a-1i|dO|s] ***
443 
444  DO la = 2, la_max
445 
446 ! *** Increase the angular momentum component z of function a ***
447 
448  sc(coset(0, 0, la), 1) = rap(3)*sc(coset(0, 0, la - 1), 1) + &
449  f2*real(la - 1, dp)*sc(coset(0, 0, la - 2), 1) - &
450  f2*kvec(3)*ss(coset(0, 0, la - 1), 1)
451  ss(coset(0, 0, la), 1) = rap(3)*ss(coset(0, 0, la - 1), 1) + &
452  f2*real(la - 1, dp)*ss(coset(0, 0, la - 2), 1) + &
453  f2*kvec(3)*sc(coset(0, 0, la - 1), 1)
454 
455 ! *** Increase the angular momentum component y of function a ***
456 
457  az = la - 1
458  sc(coset(0, 1, az), 1) = rap(2)*sc(coset(0, 0, az), 1) - &
459  f2*kvec(2)*ss(coset(0, 0, az), 1)
460  ss(coset(0, 1, az), 1) = rap(2)*ss(coset(0, 0, az), 1) + &
461  f2*kvec(2)*sc(coset(0, 0, az), 1)
462 
463  DO ay = 2, la
464  az = la - ay
465  sc(coset(0, ay, az), 1) = rap(2)*sc(coset(0, ay - 1, az), 1) + &
466  f2*real(ay - 1, dp)*sc(coset(0, ay - 2, az), 1) - &
467  f2*kvec(2)*ss(coset(0, ay - 1, az), 1)
468  ss(coset(0, ay, az), 1) = rap(2)*ss(coset(0, ay - 1, az), 1) + &
469  f2*real(ay - 1, dp)*ss(coset(0, ay - 2, az), 1) + &
470  f2*kvec(2)*sc(coset(0, ay - 1, az), 1)
471  END DO
472 
473 ! *** Increase the angular momentum component x of function a ***
474 
475  DO ay = 0, la - 1
476  az = la - 1 - ay
477  sc(coset(1, ay, az), 1) = rap(1)*sc(coset(0, ay, az), 1) - &
478  f2*kvec(1)*ss(coset(0, ay, az), 1)
479  ss(coset(1, ay, az), 1) = rap(1)*ss(coset(0, ay, az), 1) + &
480  f2*kvec(1)*sc(coset(0, ay, az), 1)
481  END DO
482 
483  DO ax = 2, la
484  f3 = f2*real(ax - 1, dp)
485  DO ay = 0, la - ax
486  az = la - ax - ay
487  sc(coset(ax, ay, az), 1) = rap(1)*sc(coset(ax - 1, ay, az), 1) + &
488  f3*sc(coset(ax - 2, ay, az), 1) - &
489  f2*kvec(1)*ss(coset(ax - 1, ay, az), 1)
490  ss(coset(ax, ay, az), 1) = rap(1)*ss(coset(ax - 1, ay, az), 1) + &
491  f3*ss(coset(ax - 2, ay, az), 1) + &
492  f2*kvec(1)*sc(coset(ax - 1, ay, az), 1)
493  END DO
494  END DO
495 
496  END DO
497 
498 ! *** Recurrence steps: [a|O|s] -> [a|O|b] ***
499 
500  IF (lb_max > 0) THEN
501 
502  DO j = 2, ncoset(lb_max)
503  DO i = 1, ncoset(la_max)
504  sc(i, j) = 0.0_dp
505  ss(i, j) = 0.0_dp
506  END DO
507  END DO
508 
509 ! *** Horizontal recurrence steps ***
510 
511  rbp(:) = rap(:) - rab(:)
512 
513 ! *** [a|O|p] = [a+1i|O|s] - (Bi - Ai)*[a|O|s] ***
514 
515  IF (lb_max == 1) THEN
516  la_start = la_min
517  ELSE
518  la_start = max(0, la_min - 1)
519  END IF
520 
521  DO la = la_start, la_max - 1
522  DO ax = 0, la
523  DO ay = 0, la - ax
524  az = la - ax - ay
525  sc(coset(ax, ay, az), 2) = sc(coset(ax + 1, ay, az), 1) - &
526  rab(1)*sc(coset(ax, ay, az), 1)
527  sc(coset(ax, ay, az), 3) = sc(coset(ax, ay + 1, az), 1) - &
528  rab(2)*sc(coset(ax, ay, az), 1)
529  sc(coset(ax, ay, az), 4) = sc(coset(ax, ay, az + 1), 1) - &
530  rab(3)*sc(coset(ax, ay, az), 1)
531  ss(coset(ax, ay, az), 2) = ss(coset(ax + 1, ay, az), 1) - &
532  rab(1)*ss(coset(ax, ay, az), 1)
533  ss(coset(ax, ay, az), 3) = ss(coset(ax, ay + 1, az), 1) - &
534  rab(2)*ss(coset(ax, ay, az), 1)
535  ss(coset(ax, ay, az), 4) = ss(coset(ax, ay, az + 1), 1) - &
536  rab(3)*ss(coset(ax, ay, az), 1)
537  END DO
538  END DO
539  END DO
540 
541 ! *** Vertical recurrence step ***
542 
543 ! *** [a|O|p] = (Pi - Bi)*[a|O|s] + f2*Ni(a)*[a-1i|O|s] ***
544 ! *** + [a|dO|s] ***
545 
546  DO ax = 0, la_max
547  fx = f2*real(ax, dp)
548  DO ay = 0, la_max - ax
549  fy = f2*real(ay, dp)
550  az = la_max - ax - ay
551  fz = f2*real(az, dp)
552  IF (ax == 0) THEN
553  sc(coset(ax, ay, az), 2) = rbp(1)*sc(coset(ax, ay, az), 1) - &
554  f2*kvec(1)*ss(coset(ax, ay, az), 1)
555  ss(coset(ax, ay, az), 2) = rbp(1)*ss(coset(ax, ay, az), 1) + &
556  f2*kvec(1)*sc(coset(ax, ay, az), 1)
557  ELSE
558  sc(coset(ax, ay, az), 2) = rbp(1)*sc(coset(ax, ay, az), 1) + &
559  fx*sc(coset(ax - 1, ay, az), 1) - &
560  f2*kvec(1)*ss(coset(ax, ay, az), 1)
561  ss(coset(ax, ay, az), 2) = rbp(1)*ss(coset(ax, ay, az), 1) + &
562  fx*ss(coset(ax - 1, ay, az), 1) + &
563  f2*kvec(1)*sc(coset(ax, ay, az), 1)
564  END IF
565  IF (ay == 0) THEN
566  sc(coset(ax, ay, az), 3) = rbp(2)*sc(coset(ax, ay, az), 1) - &
567  f2*kvec(2)*ss(coset(ax, ay, az), 1)
568  ss(coset(ax, ay, az), 3) = rbp(2)*ss(coset(ax, ay, az), 1) + &
569  f2*kvec(2)*sc(coset(ax, ay, az), 1)
570  ELSE
571  sc(coset(ax, ay, az), 3) = rbp(2)*sc(coset(ax, ay, az), 1) + &
572  fy*sc(coset(ax, ay - 1, az), 1) - &
573  f2*kvec(2)*ss(coset(ax, ay, az), 1)
574  ss(coset(ax, ay, az), 3) = rbp(2)*ss(coset(ax, ay, az), 1) + &
575  fy*ss(coset(ax, ay - 1, az), 1) + &
576  f2*kvec(2)*sc(coset(ax, ay, az), 1)
577  END IF
578  IF (az == 0) THEN
579  sc(coset(ax, ay, az), 4) = rbp(3)*sc(coset(ax, ay, az), 1) - &
580  f2*kvec(3)*ss(coset(ax, ay, az), 1)
581  ss(coset(ax, ay, az), 4) = rbp(3)*ss(coset(ax, ay, az), 1) + &
582  f2*kvec(3)*sc(coset(ax, ay, az), 1)
583  ELSE
584  sc(coset(ax, ay, az), 4) = rbp(3)*sc(coset(ax, ay, az), 1) + &
585  fz*sc(coset(ax, ay, az - 1), 1) - &
586  f2*kvec(3)*ss(coset(ax, ay, az), 1)
587  ss(coset(ax, ay, az), 4) = rbp(3)*ss(coset(ax, ay, az), 1) + &
588  fz*ss(coset(ax, ay, az - 1), 1) + &
589  f2*kvec(3)*sc(coset(ax, ay, az), 1)
590  END IF
591  END DO
592  END DO
593 
594 ! *** Recurrence steps: [a|O|p] -> [a|O|b] ***
595 
596  DO lb = 2, lb_max
597 
598 ! *** Horizontal recurrence steps ***
599 
600 ! *** [a|O|b] = [a+1i|O|b-1i] - (Bi - Ai)*[a|O|b-1i] ***
601 
602  IF (lb == lb_max) THEN
603  la_start = la_min
604  ELSE
605  la_start = max(0, la_min - 1)
606  END IF
607 
608  DO la = la_start, la_max - 1
609  DO ax = 0, la
610  DO ay = 0, la - ax
611  az = la - ax - ay
612 
613 ! *** Shift of angular momentum component z from a to b ***
614 
615  sc(coset(ax, ay, az), coset(0, 0, lb)) = &
616  sc(coset(ax, ay, az + 1), coset(0, 0, lb - 1)) - &
617  rab(3)*sc(coset(ax, ay, az), coset(0, 0, lb - 1))
618  ss(coset(ax, ay, az), coset(0, 0, lb)) = &
619  ss(coset(ax, ay, az + 1), coset(0, 0, lb - 1)) - &
620  rab(3)*ss(coset(ax, ay, az), coset(0, 0, lb - 1))
621 
622 ! *** Shift of angular momentum component y from a to b ***
623 
624  DO by = 1, lb
625  bz = lb - by
626  sc(coset(ax, ay, az), coset(0, by, bz)) = &
627  sc(coset(ax, ay + 1, az), coset(0, by - 1, bz)) - &
628  rab(2)*sc(coset(ax, ay, az), coset(0, by - 1, bz))
629  ss(coset(ax, ay, az), coset(0, by, bz)) = &
630  ss(coset(ax, ay + 1, az), coset(0, by - 1, bz)) - &
631  rab(2)*ss(coset(ax, ay, az), coset(0, by - 1, bz))
632  END DO
633 
634 ! *** Shift of angular momentum component x from a to b ***
635 
636  DO bx = 1, lb
637  DO by = 0, lb - bx
638  bz = lb - bx - by
639  sc(coset(ax, ay, az), coset(bx, by, bz)) = &
640  sc(coset(ax + 1, ay, az), coset(bx - 1, by, bz)) - &
641  rab(1)*sc(coset(ax, ay, az), coset(bx - 1, by, bz))
642  ss(coset(ax, ay, az), coset(bx, by, bz)) = &
643  ss(coset(ax + 1, ay, az), coset(bx - 1, by, bz)) - &
644  rab(1)*ss(coset(ax, ay, az), coset(bx - 1, by, bz))
645  END DO
646  END DO
647 
648  END DO
649  END DO
650  END DO
651 
652 ! *** Vertical recurrence step ***
653 
654 ! *** [a|O|b] = (Pi - Bi)*[a|O|b-1i] + f2*Ni(a)*[a-1i|O|b-1i] + ***
655 ! *** f2*Ni(b-1i)*[a|O|b-2i] + [a|dO|b-1i] ***
656 
657  DO ax = 0, la_max
658  fx = f2*real(ax, dp)
659  DO ay = 0, la_max - ax
660  fy = f2*real(ay, dp)
661  az = la_max - ax - ay
662  fz = f2*real(az, dp)
663 
664 ! *** Increase the angular momentum component z of function b ***
665 
666  f3 = f2*real(lb - 1, dp)
667 
668  IF (az == 0) THEN
669  sc(coset(ax, ay, az), coset(0, 0, lb)) = &
670  rbp(3)*sc(coset(ax, ay, az), coset(0, 0, lb - 1)) + &
671  f3*sc(coset(ax, ay, az), coset(0, 0, lb - 2)) - &
672  f2*kvec(3)*ss(coset(ax, ay, az), coset(0, 0, lb - 1))
673  ss(coset(ax, ay, az), coset(0, 0, lb)) = &
674  rbp(3)*ss(coset(ax, ay, az), coset(0, 0, lb - 1)) + &
675  f3*ss(coset(ax, ay, az), coset(0, 0, lb - 2)) + &
676  f2*kvec(3)*sc(coset(ax, ay, az), coset(0, 0, lb - 1))
677  ELSE
678  sc(coset(ax, ay, az), coset(0, 0, lb)) = &
679  rbp(3)*sc(coset(ax, ay, az), coset(0, 0, lb - 1)) + &
680  fz*sc(coset(ax, ay, az - 1), coset(0, 0, lb - 1)) + &
681  f3*sc(coset(ax, ay, az), coset(0, 0, lb - 2)) - &
682  f2*kvec(3)*ss(coset(ax, ay, az), coset(0, 0, lb - 1))
683  ss(coset(ax, ay, az), coset(0, 0, lb)) = &
684  rbp(3)*ss(coset(ax, ay, az), coset(0, 0, lb - 1)) + &
685  fz*ss(coset(ax, ay, az - 1), coset(0, 0, lb - 1)) + &
686  f3*ss(coset(ax, ay, az), coset(0, 0, lb - 2)) + &
687  f2*kvec(3)*sc(coset(ax, ay, az), coset(0, 0, lb - 1))
688  END IF
689 
690 ! *** Increase the angular momentum component y of function b ***
691 
692  IF (ay == 0) THEN
693  bz = lb - 1
694  sc(coset(ax, ay, az), coset(0, 1, bz)) = &
695  rbp(2)*sc(coset(ax, ay, az), coset(0, 0, bz)) - &
696  f2*kvec(2)*ss(coset(ax, ay, az), coset(0, 0, bz))
697  ss(coset(ax, ay, az), coset(0, 1, bz)) = &
698  rbp(2)*ss(coset(ax, ay, az), coset(0, 0, bz)) + &
699  f2*kvec(2)*sc(coset(ax, ay, az), coset(0, 0, bz))
700  DO by = 2, lb
701  bz = lb - by
702  f3 = f2*real(by - 1, dp)
703  sc(coset(ax, ay, az), coset(0, by, bz)) = &
704  rbp(2)*sc(coset(ax, ay, az), coset(0, by - 1, bz)) + &
705  f3*sc(coset(ax, ay, az), coset(0, by - 2, bz)) - &
706  f2*kvec(2)*ss(coset(ax, ay, az), coset(0, by - 1, bz))
707  ss(coset(ax, ay, az), coset(0, by, bz)) = &
708  rbp(2)*ss(coset(ax, ay, az), coset(0, by - 1, bz)) + &
709  f3*ss(coset(ax, ay, az), coset(0, by - 2, bz)) + &
710  f2*kvec(2)*sc(coset(ax, ay, az), coset(0, by - 1, bz))
711  END DO
712  ELSE
713  bz = lb - 1
714  sc(coset(ax, ay, az), coset(0, 1, bz)) = &
715  rbp(2)*sc(coset(ax, ay, az), coset(0, 0, bz)) + &
716  fy*sc(coset(ax, ay - 1, az), coset(0, 0, bz)) - &
717  f2*kvec(2)*ss(coset(ax, ay, az), coset(0, 0, bz))
718  ss(coset(ax, ay, az), coset(0, 1, bz)) = &
719  rbp(2)*ss(coset(ax, ay, az), coset(0, 0, bz)) + &
720  fy*ss(coset(ax, ay - 1, az), coset(0, 0, bz)) + &
721  f2*kvec(2)*sc(coset(ax, ay, az), coset(0, 0, bz))
722  DO by = 2, lb
723  bz = lb - by
724  f3 = f2*real(by - 1, dp)
725  sc(coset(ax, ay, az), coset(0, by, bz)) = &
726  rbp(2)*sc(coset(ax, ay, az), coset(0, by - 1, bz)) + &
727  fy*sc(coset(ax, ay - 1, az), coset(0, by - 1, bz)) + &
728  f3*sc(coset(ax, ay, az), coset(0, by - 2, bz)) - &
729  f2*kvec(2)*ss(coset(ax, ay, az), coset(0, by - 1, bz))
730  ss(coset(ax, ay, az), coset(0, by, bz)) = &
731  rbp(2)*ss(coset(ax, ay, az), coset(0, by - 1, bz)) + &
732  fy*ss(coset(ax, ay - 1, az), coset(0, by - 1, bz)) + &
733  f3*ss(coset(ax, ay, az), coset(0, by - 2, bz)) + &
734  f2*kvec(2)*sc(coset(ax, ay, az), coset(0, by - 1, bz))
735  END DO
736  END IF
737 
738 ! *** Increase the angular momentum component x of function b ***
739 
740  IF (ax == 0) THEN
741  DO by = 0, lb - 1
742  bz = lb - 1 - by
743  sc(coset(ax, ay, az), coset(1, by, bz)) = &
744  rbp(1)*sc(coset(ax, ay, az), coset(0, by, bz)) - &
745  f2*kvec(1)*ss(coset(ax, ay, az), coset(0, by, bz))
746  ss(coset(ax, ay, az), coset(1, by, bz)) = &
747  rbp(1)*ss(coset(ax, ay, az), coset(0, by, bz)) + &
748  f2*kvec(1)*sc(coset(ax, ay, az), coset(0, by, bz))
749  END DO
750  DO bx = 2, lb
751  f3 = f2*real(bx - 1, dp)
752  DO by = 0, lb - bx
753  bz = lb - bx - by
754  sc(coset(ax, ay, az), coset(bx, by, bz)) = &
755  rbp(1)*sc(coset(ax, ay, az), &
756  coset(bx - 1, by, bz)) + &
757  f3*sc(coset(ax, ay, az), coset(bx - 2, by, bz)) - &
758  f2*kvec(1)*ss(coset(ax, ay, az), coset(bx - 1, by, bz))
759  ss(coset(ax, ay, az), coset(bx, by, bz)) = &
760  rbp(1)*ss(coset(ax, ay, az), &
761  coset(bx - 1, by, bz)) + &
762  f3*ss(coset(ax, ay, az), coset(bx - 2, by, bz)) + &
763  f2*kvec(1)*sc(coset(ax, ay, az), coset(bx - 1, by, bz))
764  END DO
765  END DO
766  ELSE
767  DO by = 0, lb - 1
768  bz = lb - 1 - by
769  sc(coset(ax, ay, az), coset(1, by, bz)) = &
770  rbp(1)*sc(coset(ax, ay, az), coset(0, by, bz)) + &
771  fx*sc(coset(ax - 1, ay, az), coset(0, by, bz)) - &
772  f2*kvec(1)*ss(coset(ax, ay, az), coset(0, by, bz))
773  ss(coset(ax, ay, az), coset(1, by, bz)) = &
774  rbp(1)*ss(coset(ax, ay, az), coset(0, by, bz)) + &
775  fx*ss(coset(ax - 1, ay, az), coset(0, by, bz)) + &
776  f2*kvec(1)*sc(coset(ax, ay, az), coset(0, by, bz))
777  END DO
778  DO bx = 2, lb
779  f3 = f2*real(bx - 1, dp)
780  DO by = 0, lb - bx
781  bz = lb - bx - by
782  sc(coset(ax, ay, az), coset(bx, by, bz)) = &
783  rbp(1)*sc(coset(ax, ay, az), &
784  coset(bx - 1, by, bz)) + &
785  fx*sc(coset(ax - 1, ay, az), coset(bx - 1, by, bz)) + &
786  f3*sc(coset(ax, ay, az), coset(bx - 2, by, bz)) - &
787  f2*kvec(1)*ss(coset(ax, ay, az), coset(bx - 1, by, bz))
788  ss(coset(ax, ay, az), coset(bx, by, bz)) = &
789  rbp(1)*ss(coset(ax, ay, az), &
790  coset(bx - 1, by, bz)) + &
791  fx*ss(coset(ax - 1, ay, az), coset(bx - 1, by, bz)) + &
792  f3*ss(coset(ax, ay, az), coset(bx - 2, by, bz)) + &
793  f2*kvec(1)*sc(coset(ax, ay, az), coset(bx - 1, by, bz))
794  END DO
795  END DO
796  END IF
797 
798  END DO
799  END DO
800 
801  END DO
802 
803  END IF
804 
805  ELSE
806 
807  IF (lb_max > 0) THEN
808 
809 ! *** Vertical recurrence steps: [s|O|s] -> [s|O|b] ***
810 
811  rbp(:) = (f1 - 1.0_dp)*rab(:)
812 
813 ! *** [s|O|p] = (Pi - Bi)*[s|O|s] + [s|dO|s] ***
814 
815  sc(1, 2) = rbp(1)*sc(1, 1) - f2*kvec(1)*ss(1, 1)
816  sc(1, 3) = rbp(2)*sc(1, 1) - f2*kvec(2)*ss(1, 1)
817  sc(1, 4) = rbp(3)*sc(1, 1) - f2*kvec(3)*ss(1, 1)
818  ss(1, 2) = rbp(1)*ss(1, 1) + f2*kvec(1)*sc(1, 1)
819  ss(1, 3) = rbp(2)*ss(1, 1) + f2*kvec(2)*sc(1, 1)
820  ss(1, 4) = rbp(3)*ss(1, 1) + f2*kvec(3)*sc(1, 1)
821 
822 ! *** [s|O|b] = (Pi - Bi)*[s|O|b-1i] + f2*Ni(b-1i)*[s|O|b-2i] ***
823 ! *** + [s|dO|b-1i] ***
824 
825  DO lb = 2, lb_max
826 
827 ! *** Increase the angular momentum component z of function b ***
828 
829  sc(1, coset(0, 0, lb)) = rbp(3)*sc(1, coset(0, 0, lb - 1)) + &
830  f2*real(lb - 1, dp)*sc(1, coset(0, 0, lb - 2)) - &
831  f2*kvec(3)*ss(1, coset(0, 0, lb - 1))
832  ss(1, coset(0, 0, lb)) = rbp(3)*ss(1, coset(0, 0, lb - 1)) + &
833  f2*real(lb - 1, dp)*ss(1, coset(0, 0, lb - 2)) + &
834  f2*kvec(3)*sc(1, coset(0, 0, lb - 1))
835 
836 ! *** Increase the angular momentum component y of function b ***
837 
838  bz = lb - 1
839  sc(1, coset(0, 1, bz)) = rbp(2)*sc(1, coset(0, 0, bz)) - &
840  f2*kvec(2)*ss(1, coset(0, 0, bz))
841  ss(1, coset(0, 1, bz)) = rbp(2)*ss(1, coset(0, 0, bz)) + &
842  f2*kvec(2)*sc(1, coset(0, 0, bz))
843 
844  DO by = 2, lb
845  bz = lb - by
846  sc(1, coset(0, by, bz)) = rbp(2)*sc(1, coset(0, by - 1, bz)) + &
847  f2*real(by - 1, dp)*sc(1, coset(0, by - 2, bz)) - &
848  f2*kvec(2)*ss(1, coset(0, by - 1, bz))
849  ss(1, coset(0, by, bz)) = rbp(2)*ss(1, coset(0, by - 1, bz)) + &
850  f2*real(by - 1, dp)*ss(1, coset(0, by - 2, bz)) + &
851  f2*kvec(2)*sc(1, coset(0, by - 1, bz))
852  END DO
853 
854 ! *** Increase the angular momentum component x of function b ***
855 
856  DO by = 0, lb - 1
857  bz = lb - 1 - by
858  sc(1, coset(1, by, bz)) = rbp(1)*sc(1, coset(0, by, bz)) - &
859  f2*kvec(1)*ss(1, coset(0, by, bz))
860  ss(1, coset(1, by, bz)) = rbp(1)*ss(1, coset(0, by, bz)) + &
861  f2*kvec(1)*sc(1, coset(0, by, bz))
862  END DO
863 
864  DO bx = 2, lb
865  f3 = f2*real(bx - 1, dp)
866  DO by = 0, lb - bx
867  bz = lb - bx - by
868  sc(1, coset(bx, by, bz)) = rbp(1)*sc(1, coset(bx - 1, by, bz)) + &
869  f3*sc(1, coset(bx - 2, by, bz)) - &
870  f2*kvec(1)*ss(1, coset(bx - 1, by, bz))
871  ss(1, coset(bx, by, bz)) = rbp(1)*ss(1, coset(bx - 1, by, bz)) + &
872  f3*ss(1, coset(bx - 2, by, bz)) + &
873  f2*kvec(1)*sc(1, coset(bx - 1, by, bz))
874  END DO
875  END DO
876 
877  END DO
878 
879  END IF
880 
881  END IF
882 
883  DO j = ncoset(lb_min - 1) + 1, ncoset(lb_max)
884  DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
885  cosab(na + i, nb + j) = sc(i, j)
886  sinab(na + i, nb + j) = ss(i, j)
887  END DO
888  END DO
889 
890  IF (PRESENT(dcosab)) THEN
891  la_start = 0
892  lb_start = 0
893  ELSE
894  la_start = la_min
895  lb_start = lb_min
896  END IF
897 
898  DO da = 0, da_max - 1
899  ftz = 2.0_dp*zeta(ipgf)
900  DO dax = 0, da
901  DO day = 0, da - dax
902  daz = da - dax - day
903  cda = coset(dax, day, daz) - 1
904  cdax = coset(dax + 1, day, daz) - 1
905  cday = coset(dax, day + 1, daz) - 1
906  cdaz = coset(dax, day, daz + 1) - 1
907  !*** [da/dAi|O|b] = 2*zeta*[a+1i|O|b] - Ni(a)[a-1i|O|b] ***
908 
909  DO la = la_start, la_max - da - 1
910  DO ax = 0, la
911  fax = real(ax, dp)
912  DO ay = 0, la - ax
913  fay = real(ay, dp)
914  az = la - ax - ay
915  faz = real(az, dp)
916  coa = coset(ax, ay, az)
917  coamx = coset(ax - 1, ay, az)
918  coamy = coset(ax, ay - 1, az)
919  coamz = coset(ax, ay, az - 1)
920  coapx = coset(ax + 1, ay, az)
921  coapy = coset(ax, ay + 1, az)
922  coapz = coset(ax, ay, az + 1)
923  DO lb = lb_start, lb_max
924  DO bx = 0, lb
925  DO by = 0, lb - bx
926  bz = lb - bx - by
927  cob = coset(bx, by, bz)
928  dscos(coa, cob, cdax) = ftz*sc(coapx, cob) - fax*sc(coamx, cob)
929  dscos(coa, cob, cday) = ftz*sc(coapy, cob) - fay*sc(coamy, cob)
930  dscos(coa, cob, cdaz) = ftz*sc(coapz, cob) - faz*sc(coamz, cob)
931  dssin(coa, cob, cdax) = ftz*ss(coapx, cob) - fax*ss(coamx, cob)
932  dssin(coa, cob, cday) = ftz*ss(coapy, cob) - fay*ss(coamy, cob)
933  dssin(coa, cob, cdaz) = ftz*ss(coapz, cob) - faz*ss(coamz, cob)
934  END DO
935  END DO
936  END DO
937  END DO
938  END DO
939  END DO
940 
941  END DO
942  END DO
943  END DO
944 
945  IF (PRESENT(dcosab)) THEN
946  DO k = 1, 3
947  DO j = 1, ncoset(lb_max)
948  DO i = 1, ncoset(la_max_set)
949  dcosab(na + i, nb + j, k) = dscos(i, j, k)
950  dsinab(na + i, nb + j, k) = dssin(i, j, k)
951  END DO
952  END DO
953  END DO
954  END IF
955 
956  nb = nb + ncoset(lb_max)
957 
958  END DO
959 
960  na = na + ncoset(la_max_set)
961 
962  END DO
963 
964  END SUBROUTINE cossin
965 
966 ! **************************************************************************************************
967 !> \brief ...
968 !> \param la_max ...
969 !> \param npgfa ...
970 !> \param zeta ...
971 !> \param rpgfa ...
972 !> \param la_min ...
973 !> \param lb_max ...
974 !> \param npgfb ...
975 !> \param zetb ...
976 !> \param rpgfb ...
977 !> \param lc_max ...
978 !> \param rac ...
979 !> \param rbc ...
980 !> \param mab ...
981 ! **************************************************************************************************
982  SUBROUTINE moment(la_max, npgfa, zeta, rpgfa, la_min, &
983  lb_max, npgfb, zetb, rpgfb, &
984  lc_max, rac, rbc, mab)
985 
986  INTEGER, INTENT(IN) :: la_max, npgfa
987  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
988  INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
989  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
990  INTEGER, INTENT(IN) :: lc_max
991  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
992  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: mab
993 
994  INTEGER :: ax, ay, az, bx, by, bz, i, ipgf, j, &
995  jpgf, k, l, l1, l2, la, la_start, lb, &
996  lx, lx1, ly, ly1, lz, lz1, na, nb, ni
997  REAL(kind=dp) :: dab, f0, f1, f2, f2x, f2y, f2z, f3, fx, &
998  fy, fz, rab2, zetp
999  REAL(kind=dp), DIMENSION(3) :: rab, rap, rbp, rpc
1000  REAL(kind=dp), DIMENSION(ncoset(la_max), ncoset(& lb_max), ncoset(lc_max)) :: s
1001 
1002  rab = rbc - rac
1003  rab2 = sum(rab**2)
1004  dab = sqrt(rab2)
1005 
1006 ! *** Loop over all pairs of primitive Gaussian-type functions ***
1007 
1008  na = 0
1009 
1010  DO ipgf = 1, npgfa
1011 
1012  nb = 0
1013 
1014  DO jpgf = 1, npgfb
1015 
1016  s = 0.0_dp
1017 ! *** Screening ***
1018 
1019  IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
1020  DO k = 1, ncoset(lc_max) - 1
1021  DO j = nb + 1, nb + ncoset(lb_max)
1022  DO i = na + 1, na + ncoset(la_max)
1023  mab(i, j, k) = 0.0_dp
1024  END DO
1025  END DO
1026  END DO
1027  nb = nb + ncoset(lb_max)
1028  cycle
1029  END IF
1030 
1031 ! *** Calculate some prefactors ***
1032 
1033  zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
1034 
1035  f0 = (pi*zetp)**1.5_dp
1036  f1 = zetb(jpgf)*zetp
1037  f2 = 0.5_dp*zetp
1038 
1039 ! *** Calculate the basic two-center moment integral [s|M|s] ***
1040 
1041  rpc = zetp*(zeta(ipgf)*rac + zetb(jpgf)*rbc)
1042  s(1, 1, 1) = f0*exp(-zeta(ipgf)*f1*rab2)
1043  DO l = 2, ncoset(lc_max)
1044  lx = indco(1, l)
1045  ly = indco(2, l)
1046  lz = indco(3, l)
1047  l2 = 0
1048  IF (lz > 0) THEN
1049  l1 = coset(lx, ly, lz - 1)
1050  IF (lz > 1) l2 = coset(lx, ly, lz - 2)
1051  ni = lz - 1
1052  i = 3
1053  ELSE IF (ly > 0) THEN
1054  l1 = coset(lx, ly - 1, lz)
1055  IF (ly > 1) l2 = coset(lx, ly - 2, lz)
1056  ni = ly - 1
1057  i = 2
1058  ELSE IF (lx > 0) THEN
1059  l1 = coset(lx - 1, ly, lz)
1060  IF (lx > 1) l2 = coset(lx - 2, ly, lz)
1061  ni = lx - 1
1062  i = 1
1063  END IF
1064  s(1, 1, l) = rpc(i)*s(1, 1, l1)
1065  IF (l2 > 0) s(1, 1, l) = s(1, 1, l) + f2*real(ni, dp)*s(1, 1, l2)
1066  END DO
1067 
1068 ! *** Recurrence steps: [s|M|s] -> [a|M|b] ***
1069 
1070  DO l = 1, ncoset(lc_max)
1071 
1072  lx = indco(1, l)
1073  ly = indco(2, l)
1074  lz = indco(3, l)
1075  IF (lx > 0) THEN
1076  lx1 = coset(lx - 1, ly, lz)
1077  ELSE
1078  lx1 = -1
1079  END IF
1080  IF (ly > 0) THEN
1081  ly1 = coset(lx, ly - 1, lz)
1082  ELSE
1083  ly1 = -1
1084  END IF
1085  IF (lz > 0) THEN
1086  lz1 = coset(lx, ly, lz - 1)
1087  ELSE
1088  lz1 = -1
1089  END IF
1090  f2x = f2*real(lx, dp)
1091  f2y = f2*real(ly, dp)
1092  f2z = f2*real(lz, dp)
1093 
1094  IF (la_max > 0) THEN
1095 
1096 ! *** Vertical recurrence steps: [s|M|s] -> [a|M|s] ***
1097 
1098  rap(:) = f1*rab(:)
1099 
1100 ! *** [p|M|s] = (Pi - Ai)*[s|M|s] + f2*Ni(m-1i)[s|M-1i|s] ***
1101 
1102  s(2, 1, l) = rap(1)*s(1, 1, l)
1103  s(3, 1, l) = rap(2)*s(1, 1, l)
1104  s(4, 1, l) = rap(3)*s(1, 1, l)
1105  IF (lx1 > 0) s(2, 1, l) = s(2, 1, l) + f2x*s(1, 1, lx1)
1106  IF (ly1 > 0) s(3, 1, l) = s(3, 1, l) + f2y*s(1, 1, ly1)
1107  IF (lz1 > 0) s(4, 1, l) = s(4, 1, l) + f2z*s(1, 1, lz1)
1108 
1109 ! *** [a|M|s] = (Pi - Ai)*[a-1i|M|s] + f2*Ni(a-1i)*[a-2i|M|s] ***
1110 ! *** + f2*Ni(m-1i)*[a-1i|M-1i|s] ***
1111 
1112  DO la = 2, la_max
1113 
1114 ! *** Increase the angular momentum component z of function a ***
1115 
1116  s(coset(0, 0, la), 1, l) = rap(3)*s(coset(0, 0, la - 1), 1, l) + &
1117  f2*real(la - 1, dp)*s(coset(0, 0, la - 2), 1, l)
1118  IF (lz1 > 0) s(coset(0, 0, la), 1, l) = s(coset(0, 0, la), 1, l) + &
1119  f2z*s(coset(0, 0, la - 1), 1, lz1)
1120 
1121 ! *** Increase the angular momentum component y of function a ***
1122 
1123  az = la - 1
1124  s(coset(0, 1, az), 1, l) = rap(2)*s(coset(0, 0, az), 1, l)
1125  IF (ly1 > 0) s(coset(0, 1, az), 1, l) = s(coset(0, 1, az), 1, l) + &
1126  f2y*s(coset(0, 0, az), 1, ly1)
1127 
1128  DO ay = 2, la
1129  az = la - ay
1130  s(coset(0, ay, az), 1, l) = rap(2)*s(coset(0, ay - 1, az), 1, l) + &
1131  f2*real(ay - 1, dp)*s(coset(0, ay - 2, az), 1, l)
1132  IF (ly1 > 0) s(coset(0, ay, az), 1, l) = s(coset(0, ay, az), 1, l) + &
1133  f2y*s(coset(0, ay - 1, az), 1, ly1)
1134  END DO
1135 
1136 ! *** Increase the angular momentum component x of function a ***
1137 
1138  DO ay = 0, la - 1
1139  az = la - 1 - ay
1140  s(coset(1, ay, az), 1, l) = rap(1)*s(coset(0, ay, az), 1, l)
1141  IF (lx1 > 0) s(coset(1, ay, az), 1, l) = s(coset(1, ay, az), 1, l) + &
1142  f2x*s(coset(0, ay, az), 1, lx1)
1143  END DO
1144 
1145  DO ax = 2, la
1146  f3 = f2*real(ax - 1, dp)
1147  DO ay = 0, la - ax
1148  az = la - ax - ay
1149  s(coset(ax, ay, az), 1, l) = rap(1)*s(coset(ax - 1, ay, az), 1, l) + &
1150  f3*s(coset(ax - 2, ay, az), 1, l)
1151  IF (lx1 > 0) s(coset(ax, ay, az), 1, l) = s(coset(ax, ay, az), 1, l) + &
1152  f2x*s(coset(ax - 1, ay, az), 1, lx1)
1153  END DO
1154  END DO
1155 
1156  END DO
1157 
1158 ! *** Recurrence steps: [a|M|s] -> [a|M|b] ***
1159 
1160  IF (lb_max > 0) THEN
1161 
1162  DO j = 2, ncoset(lb_max)
1163  DO i = 1, ncoset(la_max)
1164  s(i, j, l) = 0.0_dp
1165  END DO
1166  END DO
1167 
1168 ! *** Horizontal recurrence steps ***
1169 
1170  rbp(:) = rap(:) - rab(:)
1171 
1172 ! *** [a|M|p] = [a+1i|M|s] - (Bi - Ai)*[a|M|s] ***
1173 
1174  IF (lb_max == 1) THEN
1175  la_start = la_min
1176  ELSE
1177  la_start = max(0, la_min - 1)
1178  END IF
1179 
1180  DO la = la_start, la_max - 1
1181  DO ax = 0, la
1182  DO ay = 0, la - ax
1183  az = la - ax - ay
1184  s(coset(ax, ay, az), 2, l) = s(coset(ax + 1, ay, az), 1, l) - &
1185  rab(1)*s(coset(ax, ay, az), 1, l)
1186  s(coset(ax, ay, az), 3, l) = s(coset(ax, ay + 1, az), 1, l) - &
1187  rab(2)*s(coset(ax, ay, az), 1, l)
1188  s(coset(ax, ay, az), 4, l) = s(coset(ax, ay, az + 1), 1, l) - &
1189  rab(3)*s(coset(ax, ay, az), 1, l)
1190  END DO
1191  END DO
1192  END DO
1193 
1194 ! *** Vertical recurrence step ***
1195 
1196 ! *** [a|M|p] = (Pi - Bi)*[a|M|s] + f2*Ni(a)*[a-1i|M|s] ***
1197 ! *** + f2*Ni(m)*[a|M-1i|s] ***
1198 
1199  DO ax = 0, la_max
1200  fx = f2*real(ax, dp)
1201  DO ay = 0, la_max - ax
1202  fy = f2*real(ay, dp)
1203  az = la_max - ax - ay
1204  fz = f2*real(az, dp)
1205  IF (ax == 0) THEN
1206  s(coset(ax, ay, az), 2, l) = rbp(1)*s(coset(ax, ay, az), 1, l)
1207  ELSE
1208  s(coset(ax, ay, az), 2, l) = rbp(1)*s(coset(ax, ay, az), 1, l) + &
1209  fx*s(coset(ax - 1, ay, az), 1, l)
1210  END IF
1211  IF (lx1 > 0) s(coset(ax, ay, az), 2, l) = s(coset(ax, ay, az), 2, l) + &
1212  f2x*s(coset(ax, ay, az), 1, lx1)
1213  IF (ay == 0) THEN
1214  s(coset(ax, ay, az), 3, l) = rbp(2)*s(coset(ax, ay, az), 1, l)
1215  ELSE
1216  s(coset(ax, ay, az), 3, l) = rbp(2)*s(coset(ax, ay, az), 1, l) + &
1217  fy*s(coset(ax, ay - 1, az), 1, l)
1218  END IF
1219  IF (ly1 > 0) s(coset(ax, ay, az), 3, l) = s(coset(ax, ay, az), 3, l) + &
1220  f2y*s(coset(ax, ay, az), 1, ly1)
1221  IF (az == 0) THEN
1222  s(coset(ax, ay, az), 4, l) = rbp(3)*s(coset(ax, ay, az), 1, l)
1223  ELSE
1224  s(coset(ax, ay, az), 4, l) = rbp(3)*s(coset(ax, ay, az), 1, l) + &
1225  fz*s(coset(ax, ay, az - 1), 1, l)
1226  END IF
1227  IF (lz1 > 0) s(coset(ax, ay, az), 4, l) = s(coset(ax, ay, az), 4, l) + &
1228  f2z*s(coset(ax, ay, az), 1, lz1)
1229  END DO
1230  END DO
1231 
1232 ! *** Recurrence steps: [a|M|p] -> [a|M|b] ***
1233 
1234  DO lb = 2, lb_max
1235 
1236 ! *** Horizontal recurrence steps ***
1237 
1238 ! *** [a|M|b] = [a+1i|M|b-1i] - (Bi - Ai)*[a|M|b-1i] ***
1239 
1240  IF (lb == lb_max) THEN
1241  la_start = la_min
1242  ELSE
1243  la_start = max(0, la_min - 1)
1244  END IF
1245 
1246  DO la = la_start, la_max - 1
1247  DO ax = 0, la
1248  DO ay = 0, la - ax
1249  az = la - ax - ay
1250 
1251 ! *** Shift of angular momentum component z from a to b ***
1252 
1253  s(coset(ax, ay, az), coset(0, 0, lb), l) = &
1254  s(coset(ax, ay, az + 1), coset(0, 0, lb - 1), l) - &
1255  rab(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1), l)
1256 
1257 ! *** Shift of angular momentum component y from a to b ***
1258 
1259  DO by = 1, lb
1260  bz = lb - by
1261  s(coset(ax, ay, az), coset(0, by, bz), l) = &
1262  s(coset(ax, ay + 1, az), coset(0, by - 1, bz), l) - &
1263  rab(2)*s(coset(ax, ay, az), coset(0, by - 1, bz), l)
1264  END DO
1265 
1266 ! *** Shift of angular momentum component x from a to b ***
1267 
1268  DO bx = 1, lb
1269  DO by = 0, lb - bx
1270  bz = lb - bx - by
1271  s(coset(ax, ay, az), coset(bx, by, bz), l) = &
1272  s(coset(ax + 1, ay, az), coset(bx - 1, by, bz), l) - &
1273  rab(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz), l)
1274  END DO
1275  END DO
1276 
1277  END DO
1278  END DO
1279  END DO
1280 
1281 ! *** Vertical recurrence step ***
1282 
1283 ! *** [a|M|b] = (Pi - Bi)*[a|M|b-1i] + f2*Ni(a)*[a-1i|M|b-1i] + ***
1284 ! *** f2*Ni(b-1i)*[a|M|b-2i] + f2*Ni(m)[a|M-1i|b-1i] ***
1285 
1286  DO ax = 0, la_max
1287  fx = f2*real(ax, dp)
1288  DO ay = 0, la_max - ax
1289  fy = f2*real(ay, dp)
1290  az = la_max - ax - ay
1291  fz = f2*real(az, dp)
1292 
1293 ! *** Shift of angular momentum component z from a to b ***
1294 
1295  f3 = f2*real(lb - 1, dp)
1296 
1297  IF (az == 0) THEN
1298  s(coset(ax, ay, az), coset(0, 0, lb), l) = &
1299  rbp(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1), l) + &
1300  f3*s(coset(ax, ay, az), coset(0, 0, lb - 2), l)
1301  ELSE
1302  s(coset(ax, ay, az), coset(0, 0, lb), l) = &
1303  rbp(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1), l) + &
1304  fz*s(coset(ax, ay, az - 1), coset(0, 0, lb - 1), l) + &
1305  f3*s(coset(ax, ay, az), coset(0, 0, lb - 2), l)
1306  END IF
1307  IF (lz1 > 0) s(coset(ax, ay, az), coset(0, 0, lb), l) = &
1308  s(coset(ax, ay, az), coset(0, 0, lb), l) + &
1309  f2z*s(coset(ax, ay, az), coset(0, 0, lb - 1), lz1)
1310 
1311 ! *** Shift of angular momentum component y from a to b ***
1312 
1313  IF (ay == 0) THEN
1314  bz = lb - 1
1315  s(coset(ax, ay, az), coset(0, 1, bz), l) = &
1316  rbp(2)*s(coset(ax, ay, az), coset(0, 0, bz), l)
1317  IF (ly1 > 0) s(coset(ax, ay, az), coset(0, 1, bz), l) = &
1318  s(coset(ax, ay, az), coset(0, 1, bz), l) + &
1319  f2y*s(coset(ax, ay, az), coset(0, 0, bz), ly1)
1320  DO by = 2, lb
1321  bz = lb - by
1322  f3 = f2*real(by - 1, dp)
1323  s(coset(ax, ay, az), coset(0, by, bz), l) = &
1324  rbp(2)*s(coset(ax, ay, az), coset(0, by - 1, bz), l) + &
1325  f3*s(coset(ax, ay, az), coset(0, by - 2, bz), l)
1326  IF (ly1 > 0) s(coset(ax, ay, az), coset(0, by, bz), l) = &
1327  s(coset(ax, ay, az), coset(0, by, bz), l) + &
1328  f2y*s(coset(ax, ay, az), coset(0, by - 1, bz), ly1)
1329  END DO
1330  ELSE
1331  bz = lb - 1
1332  s(coset(ax, ay, az), coset(0, 1, bz), l) = &
1333  rbp(2)*s(coset(ax, ay, az), coset(0, 0, bz), l) + &
1334  fy*s(coset(ax, ay - 1, az), coset(0, 0, bz), l)
1335  IF (ly1 > 0) s(coset(ax, ay, az), coset(0, 1, bz), l) = &
1336  s(coset(ax, ay, az), coset(0, 1, bz), l) + &
1337  f2y*s(coset(ax, ay, az), coset(0, 0, bz), ly1)
1338  DO by = 2, lb
1339  bz = lb - by
1340  f3 = f2*real(by - 1, dp)
1341  s(coset(ax, ay, az), coset(0, by, bz), l) = &
1342  rbp(2)*s(coset(ax, ay, az), coset(0, by - 1, bz), l) + &
1343  fy*s(coset(ax, ay - 1, az), coset(0, by - 1, bz), l) + &
1344  f3*s(coset(ax, ay, az), coset(0, by - 2, bz), l)
1345  IF (ly1 > 0) s(coset(ax, ay, az), coset(0, by, bz), l) = &
1346  s(coset(ax, ay, az), coset(0, by, bz), l) + &
1347  f2y*s(coset(ax, ay, az), coset(0, by - 1, bz), ly1)
1348  END DO
1349  END IF
1350 
1351 ! *** Shift of angular momentum component x from a to b ***
1352 
1353  IF (ax == 0) THEN
1354  DO by = 0, lb - 1
1355  bz = lb - 1 - by
1356  s(coset(ax, ay, az), coset(1, by, bz), l) = &
1357  rbp(1)*s(coset(ax, ay, az), coset(0, by, bz), l)
1358  IF (lx1 > 0) s(coset(ax, ay, az), coset(1, by, bz), l) = &
1359  s(coset(ax, ay, az), coset(1, by, bz), l) + &
1360  f2x*s(coset(ax, ay, az), coset(0, by, bz), lx1)
1361  END DO
1362  DO bx = 2, lb
1363  f3 = f2*real(bx - 1, dp)
1364  DO by = 0, lb - bx
1365  bz = lb - bx - by
1366  s(coset(ax, ay, az), coset(bx, by, bz), l) = &
1367  rbp(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz), l) + &
1368  f3*s(coset(ax, ay, az), coset(bx - 2, by, bz), l)
1369  IF (lx1 > 0) s(coset(ax, ay, az), coset(bx, by, bz), l) = &
1370  s(coset(ax, ay, az), coset(bx, by, bz), l) + &
1371  f2x*s(coset(ax, ay, az), coset(bx - 1, by, bz), lx1)
1372  END DO
1373  END DO
1374  ELSE
1375  DO by = 0, lb - 1
1376  bz = lb - 1 - by
1377  s(coset(ax, ay, az), coset(1, by, bz), l) = &
1378  rbp(1)*s(coset(ax, ay, az), coset(0, by, bz), l) + &
1379  fx*s(coset(ax - 1, ay, az), coset(0, by, bz), l)
1380  IF (lx1 > 0) s(coset(ax, ay, az), coset(1, by, bz), l) = &
1381  s(coset(ax, ay, az), coset(1, by, bz), l) + &
1382  f2x*s(coset(ax, ay, az), coset(0, by, bz), lx1)
1383  END DO
1384  DO bx = 2, lb
1385  f3 = f2*real(bx - 1, dp)
1386  DO by = 0, lb - bx
1387  bz = lb - bx - by
1388  s(coset(ax, ay, az), coset(bx, by, bz), l) = &
1389  rbp(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz), l) + &
1390  fx*s(coset(ax - 1, ay, az), coset(bx - 1, by, bz), l) + &
1391  f3*s(coset(ax, ay, az), coset(bx - 2, by, bz), l)
1392  IF (lx1 > 0) s(coset(ax, ay, az), coset(bx, by, bz), l) = &
1393  s(coset(ax, ay, az), coset(bx, by, bz), l) + &
1394  f2x*s(coset(ax, ay, az), coset(bx - 1, by, bz), lx1)
1395  END DO
1396  END DO
1397  END IF
1398 
1399  END DO
1400  END DO
1401 
1402  END DO
1403 
1404  END IF
1405 
1406  ELSE
1407 
1408  IF (lb_max > 0) THEN
1409 
1410 ! *** Vertical recurrence steps: [s|M|s] -> [s|M|b] ***
1411 
1412  rbp(:) = (f1 - 1.0_dp)*rab(:)
1413 
1414 ! *** [s|M|p] = (Pi - Bi)*[s|M|s] + f2*Ni(m)*[s|M-1i|s] ***
1415 
1416  s(1, 2, l) = rbp(1)*s(1, 1, l)
1417  s(1, 3, l) = rbp(2)*s(1, 1, l)
1418  s(1, 4, l) = rbp(3)*s(1, 1, l)
1419  IF (lx1 > 0) s(1, 2, l) = s(1, 2, l) + f2x*s(1, 1, lx1)
1420  IF (ly1 > 0) s(1, 3, l) = s(1, 3, l) + f2y*s(1, 1, ly1)
1421  IF (lz1 > 0) s(1, 4, l) = s(1, 4, l) + f2z*s(1, 1, lz1)
1422 
1423 ! *** [s|M|b] = (Pi - Bi)*[s|M|b-1i] + f2*Ni(b-1i)*[s|M|b-2i] ***
1424 ! *** + f2*Ni(m)*[s|M-1i|b-1i] ***
1425 
1426  DO lb = 2, lb_max
1427 
1428 ! *** Increase the angular momentum component z of function b ***
1429 
1430  s(1, coset(0, 0, lb), l) = rbp(3)*s(1, coset(0, 0, lb - 1), l) + &
1431  f2*real(lb - 1, dp)*s(1, coset(0, 0, lb - 2), l)
1432  IF (lz1 > 0) s(1, coset(0, 0, lb), l) = s(1, coset(0, 0, lb), l) + &
1433  f2z*s(1, coset(0, 0, lb - 1), lz1)
1434 
1435 ! *** Increase the angular momentum component y of function b ***
1436 
1437  bz = lb - 1
1438  s(1, coset(0, 1, bz), l) = rbp(2)*s(1, coset(0, 0, bz), l)
1439  IF (ly1 > 0) s(1, coset(0, 1, bz), l) = s(1, coset(0, 1, bz), l) + &
1440  f2y*s(1, coset(0, 0, bz), ly1)
1441 
1442  DO by = 2, lb
1443  bz = lb - by
1444  s(1, coset(0, by, bz), l) = rbp(2)*s(1, coset(0, by - 1, bz), l) + &
1445  f2*real(by - 1, dp)*s(1, coset(0, by - 2, bz), l)
1446  IF (ly1 > 0) s(1, coset(0, by, bz), l) = s(1, coset(0, by, bz), l) + &
1447  f2y*s(1, coset(0, by - 1, bz), ly1)
1448  END DO
1449 
1450 ! *** Increase the angular momentum component x of function b ***
1451 
1452  DO by = 0, lb - 1
1453  bz = lb - 1 - by
1454  s(1, coset(1, by, bz), l) = rbp(1)*s(1, coset(0, by, bz), l)
1455  IF (lx1 > 0) s(1, coset(1, by, bz), l) = s(1, coset(1, by, bz), l) + &
1456  f2x*s(1, coset(0, by, bz), lx1)
1457  END DO
1458 
1459  DO bx = 2, lb
1460  f3 = f2*real(bx - 1, dp)
1461  DO by = 0, lb - bx
1462  bz = lb - bx - by
1463  s(1, coset(bx, by, bz), l) = rbp(1)*s(1, coset(bx - 1, by, bz), l) + &
1464  f3*s(1, coset(bx - 2, by, bz), l)
1465  IF (lx1 > 0) s(1, coset(bx, by, bz), l) = s(1, coset(bx, by, bz), l) + &
1466  f2x*s(1, coset(bx - 1, by, bz), lx1)
1467  END DO
1468  END DO
1469 
1470  END DO
1471 
1472  END IF
1473 
1474  END IF
1475 
1476  END DO
1477 
1478  DO k = 2, ncoset(lc_max)
1479  DO j = 1, ncoset(lb_max)
1480  DO i = 1, ncoset(la_max)
1481  mab(na + i, nb + j, k - 1) = s(i, j, k)
1482  END DO
1483  END DO
1484  END DO
1485 
1486  nb = nb + ncoset(lb_max)
1487 
1488  END DO
1489 
1490  na = na + ncoset(la_max)
1491 
1492  END DO
1493 
1494  END SUBROUTINE moment
1495 
1496 ! **************************************************************************************************
1497 !> \brief This returns the derivative of the overlap integrals [a|b], with respect
1498 !> to the position of the primitive on the left, i.e.
1499 !> [da/dR_ai|b] = 2*zeta*[a+1i|b] - Ni(a)[a-1i|b]
1500 !> \param la_max ...
1501 !> \param npgfa ...
1502 !> \param zeta ...
1503 !> \param rpgfa ...
1504 !> \param la_min ...
1505 !> \param lb_max ...
1506 !> \param npgfb ...
1507 !> \param zetb ...
1508 !> \param rpgfb ...
1509 !> \param lb_min ...
1510 !> \param rab ...
1511 !> \param difab ...
1512 ! **************************************************************************************************
1513  SUBROUTINE diffop(la_max, npgfa, zeta, rpgfa, la_min, &
1514  lb_max, npgfb, zetb, rpgfb, lb_min, &
1515  rab, difab)
1516 
1517  INTEGER, INTENT(IN) :: la_max, npgfa
1518  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
1519  INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
1520  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
1521  INTEGER, INTENT(IN) :: lb_min
1522  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
1523  REAL(kind=dp), DIMENSION(:, :, :), INTENT(OUT) :: difab
1524 
1525  INTEGER :: lda_min, ldb_min, lds, lmax
1526  REAL(kind=dp) :: dab, rab2
1527  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: s
1528  REAL(kind=dp), DIMENSION(npgfa*ncoset(la_max+1), & npgfb*ncoset(lb_max+1)) :: sab
1529 
1530  rab2 = sum(rab**2)
1531  dab = sqrt(rab2)
1532 
1533  lda_min = max(0, la_min - 1)
1534  ldb_min = max(0, lb_min - 1)
1535  lmax = max(la_max + 1, lb_max + 1)
1536  lds = ncoset(lmax + 1)
1537  ALLOCATE (s(lds, lds, ncoset(1)))
1538  sab = 0.0_dp
1539  s = 0.0_dp
1540  CALL overlap(la_max + 1, lda_min, npgfa, rpgfa, zeta, &
1541  lb_max + 1, ldb_min, npgfb, rpgfb, zetb, &
1542  rab, dab, sab, 0, .false., s, lds)
1543 
1544  CALL dabdr(la_max, npgfa, zeta, rpgfa, la_min, &
1545  lb_max, npgfb, rpgfb, lb_min, &
1546  dab, sab, difab(:, :, 1), difab(:, :, 2), difab(:, :, 3))
1547 
1548  DEALLOCATE (s)
1549 
1550  END SUBROUTINE diffop
1551 
1552 ! **************************************************************************************************
1553 !> \brief This returns the derivative of the moment integrals [a|\mu|b], with respect
1554 !> to the position of the primitive on the left, i.e.
1555 !> [da/dR_ai|\mu|b] = 2*zeta*[a+1i|\mu|b] - Ni(a)[a-1i|\mu|b]
1556 !> order indicates the max order of the moment operator to be calculated
1557 !> 1: dipole
1558 !> 2: quadrupole
1559 !> ...
1560 !> \param la_max ...
1561 !> \param npgfa ...
1562 !> \param zeta ...
1563 !> \param rpgfa ...
1564 !> \param la_min ...
1565 !> \param lb_max ...
1566 !> \param npgfb ...
1567 !> \param zetb ...
1568 !> \param rpgfb ...
1569 !> \param lb_min ...
1570 !> \param order ...
1571 !> \param rac ...
1572 !> \param rbc ...
1573 !> \param difmab ...
1574 !> \param mab_ext ...
1575 !> \note
1576 ! **************************************************************************************************
1577  SUBROUTINE diff_momop(la_max, npgfa, zeta, rpgfa, la_min, &
1578  lb_max, npgfb, zetb, rpgfb, lb_min, &
1579  order, rac, rbc, difmab, mab_ext)
1581  INTEGER, INTENT(IN) :: la_max, npgfa
1582  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
1583  INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
1584  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
1585  INTEGER, INTENT(IN) :: lb_min, order
1586  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
1587  REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(OUT) :: difmab
1588  REAL(kind=dp), DIMENSION(:, :, :), OPTIONAL, &
1589  POINTER :: mab_ext
1590 
1591  INTEGER :: imom, lda, lda_min, ldb, ldb_min, lmax
1592  REAL(kind=dp) :: dab, rab(3), rab2
1593  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: difmab_tmp
1594  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: mab
1595 
1596  rab = rbc - rac
1597  rab2 = sum(rab**2)
1598  dab = sqrt(rab2)
1599 
1600  lda_min = max(0, la_min - 1)
1601  ldb_min = max(0, lb_min - 1)
1602  lmax = max(la_max + 1, lb_max + 1)
1603  lda = ncoset(la_max)*npgfa
1604  ldb = ncoset(lb_max)*npgfb
1605  ALLOCATE (difmab_tmp(lda, ldb, 3))
1606 
1607  IF (PRESENT(mab_ext)) THEN
1608  mab => mab_ext
1609  ELSE
1610  ALLOCATE (mab(npgfa*ncoset(la_max + 1), npgfb*ncoset(lb_max + 1), &
1611  ncoset(order) - 1))
1612  mab = 0.0_dp
1613 ! *** Calculate the primitive overlap integrals ***
1614  CALL moment(la_max + 1, npgfa, zeta, rpgfa, lda_min, &
1615  lb_max + 1, npgfb, zetb, rpgfb, &
1616  order, rac, rbc, mab)
1617 
1618  END IF
1619  DO imom = 1, ncoset(order) - 1
1620  difmab_tmp = 0.0_dp
1621  CALL adbdr(la_max, npgfa, rpgfa, la_min, &
1622  lb_max, npgfb, zetb, rpgfb, lb_min, &
1623  dab, mab(:, :, imom), difmab_tmp(:, :, 1), &
1624  difmab_tmp(:, :, 2), difmab_tmp(:, :, 3))
1625  difmab(1:lda, 1:ldb, imom, 1) = difmab_tmp(1:lda, 1:ldb, 1)
1626  difmab(1:lda, 1:ldb, imom, 2) = difmab_tmp(1:lda, 1:ldb, 2)
1627  difmab(1:lda, 1:ldb, imom, 3) = difmab_tmp(1:lda, 1:ldb, 3)
1628  END DO
1629 
1630  IF (PRESENT(mab_ext)) THEN
1631  NULLIFY (mab)
1632  ELSE
1633  DEALLOCATE (mab)
1634  END IF
1635  DEALLOCATE (difmab_tmp)
1636 
1637  END SUBROUTINE diff_momop
1638 
1639 ! **************************************************************************************************
1640 !> \brief This returns the derivative of the dipole integrals [a|x|b], with respect
1641 !> to the position of the primitive on the left and right, i.e.
1642 !> [da/dR_ai|\mu|b] = 2*zeta*[a+1i|\mu|b] - Ni(a)[a-1i|\mu|b]
1643 !> \param la_max ...
1644 !> \param npgfa ...
1645 !> \param zeta ...
1646 !> \param rpgfa ...
1647 !> \param la_min ...
1648 !> \param lb_max ...
1649 !> \param npgfb ...
1650 !> \param zetb ...
1651 !> \param rpgfb ...
1652 !> \param lb_min ...
1653 !> \param order ...
1654 !> \param rac ...
1655 !> \param rbc ...
1656 !> \param pab ...
1657 !> \param forcea ...
1658 !> \param forceb ...
1659 !> \note
1660 ! **************************************************************************************************
1661  SUBROUTINE dipole_force(la_max, npgfa, zeta, rpgfa, la_min, &
1662  lb_max, npgfb, zetb, rpgfb, lb_min, &
1663  order, rac, rbc, pab, forcea, forceb)
1665  INTEGER, INTENT(IN) :: la_max, npgfa
1666  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
1667  INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
1668  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
1669  INTEGER, INTENT(IN) :: lb_min, order
1670  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
1671  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: pab
1672  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: forcea, forceb
1673 
1674  INTEGER :: i, imom, ipgf, j, jpgf, lda, lda_min, &
1675  ldb, ldb_min, lmax, na, nb
1676  REAL(kind=dp) :: dab, rab(3), rab2
1677  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: difmab, mab
1678 
1679  cpassert(order == 1)
1680  mark_used(order)
1681 
1682  rab = rbc - rac
1683  rab2 = sum(rab**2)
1684  dab = sqrt(rab2)
1685 
1686  lda_min = max(0, la_min - 1)
1687  ldb_min = max(0, lb_min - 1)
1688  lmax = max(la_max + 1, lb_max + 1)
1689  lda = ncoset(la_max)*npgfa
1690  ldb = ncoset(lb_max)*npgfb
1691  ALLOCATE (difmab(lda, ldb, 3))
1692  ALLOCATE (mab(npgfa*ncoset(la_max + 1), npgfb*ncoset(lb_max + 1), 3))
1693  mab = 0.0_dp
1694  CALL moment(la_max + 1, npgfa, zeta, rpgfa, lda_min, &
1695  lb_max + 1, npgfb, zetb, rpgfb, 1, rac, rbc, mab)
1696 
1697  DO imom = 1, 3
1698  difmab = 0.0_dp
1699  CALL adbdr(la_max, npgfa, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, &
1700  dab, mab(:, :, imom), difmab(:, :, 1), difmab(:, :, 2), difmab(:, :, 3))
1701  na = 0
1702  DO ipgf = 1, npgfa
1703  nb = 0
1704  DO jpgf = 1, npgfb
1705  DO j = nb + ncoset(lb_min - 1) + 1, nb + ncoset(lb_max)
1706  DO i = na + ncoset(la_min - 1) + 1, na + ncoset(la_max)
1707  forceb(imom, 1) = forceb(imom, 1) + pab(i, j)*difmab(i, j, 1)
1708  forceb(imom, 2) = forceb(imom, 2) + pab(i, j)*difmab(i, j, 2)
1709  forceb(imom, 3) = forceb(imom, 3) + pab(i, j)*difmab(i, j, 3)
1710  END DO
1711  END DO
1712  nb = nb + ncoset(lb_max)
1713  END DO
1714  na = na + ncoset(la_max)
1715  END DO
1716 
1717  difmab = 0.0_dp
1718  CALL dabdr(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, rpgfb, lb_min, &
1719  dab, mab(:, :, imom), difmab(:, :, 1), difmab(:, :, 2), difmab(:, :, 3))
1720  na = 0
1721  DO ipgf = 1, npgfa
1722  nb = 0
1723  DO jpgf = 1, npgfb
1724  DO j = nb + ncoset(lb_min - 1) + 1, nb + ncoset(lb_max)
1725  DO i = na + ncoset(la_min - 1) + 1, na + ncoset(la_max)
1726  forcea(imom, 1) = forcea(imom, 1) + pab(i, j)*difmab(i, j, 1)
1727  forcea(imom, 2) = forcea(imom, 2) + pab(i, j)*difmab(i, j, 2)
1728  forcea(imom, 3) = forcea(imom, 3) + pab(i, j)*difmab(i, j, 3)
1729  END DO
1730  END DO
1731  nb = nb + ncoset(lb_max)
1732  END DO
1733  na = na + ncoset(la_max)
1734  END DO
1735  END DO
1736 
1737  DEALLOCATE (mab, difmab)
1738 
1739  END SUBROUTINE dipole_force
1740 
1741 END MODULE ai_moments
1742 
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Calculate the first derivative of an integral block.
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...
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...
Calculation of the moment integrals over Cartesian Gaussian-type functions.
Definition: ai_moments.F:17
subroutine, public diff_momop(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, order, rac, rbc, difmab, mab_ext)
This returns the derivative of the moment integrals [a|\mu|b], with respect to the position of the pr...
Definition: ai_moments.F:1583
subroutine, public diff_momop_velocity(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, order, rac, rbc, difmab, lambda, iatom, jatom)
This returns the derivative of the moment integrals [a|\mu|b], with respect to the primitive on the r...
Definition: ai_moments.F:79
subroutine, public dipole_force(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, order, rac, rbc, pab, forcea, forceb)
This returns the derivative of the dipole integrals [a|x|b], with respect to the position of the prim...
Definition: ai_moments.F:1667
subroutine, public contract_cossin(cos_block, sin_block, iatom, ncoa, nsgfa, sgfa, sphi_a, ldsa, jatom, ncob, nsgfb, sgfb, sphi_b, ldsb, cosab, sinab, ldab, work, ldwork)
...
Definition: ai_moments.F:261
subroutine, public diffop(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, rab, difab)
This returns the derivative of the overlap integrals [a|b], with respect to the position of the primi...
Definition: ai_moments.F:1518
subroutine, public cossin(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max, npgfb, zetb, rpgfb, lb_min, rac, rbc, kvec, cosab, sinab, dcosab, dsinab)
...
Definition: ai_moments.F:339
subroutine, public moment(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lc_max, rac, rbc, mab)
...
Definition: ai_moments.F:986
subroutine, public diff_momop2(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, order, rac, rbc, difmab, mab_ext, deltaR, iatom, jatom)
This returns the derivative of the moment integrals [a|\mu|b], with respect to the position of the pr...
Definition: ai_moments.F:164
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition: ai_overlap.F:18
subroutine, public overlap(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, rab, dab, sab, da_max_set, return_derivatives, s, lds, sdab, pab, force_a)
Purpose: Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions...
Definition: ai_overlap.F:73
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public pi
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public ncoset
integer, dimension(:, :, :), allocatable, public coset
integer, dimension(:, :), allocatable, public indco