Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 kinds, ONLY: dp
35 : USE mathconstants, ONLY: pi
36 : USE orbital_pointers, ONLY: coset,&
37 : indco,&
38 : ncoset
39 : #include "../base/base_uses.f90"
40 :
41 : IMPLICIT NONE
42 :
43 : PRIVATE
44 :
45 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_moments'
46 :
47 : PUBLIC :: cossin, moment, diff_momop, contract_cossin, dipole_force
48 : PUBLIC :: diff_momop2, diff_momop_velocity
49 :
50 : CONTAINS
51 :
52 : ! *****************************************************************************
53 : !> \brief This returns the derivative of the moment integrals [a|\mu|b], with respect
54 : !> to the primitive on the right
55 : !> difmab(:, :, beta, alpha) = < a | r_beta | ∂_alpha b > * (iatom - jatom)
56 : !> \param la_max ...
57 : !> \param npgfa ...
58 : !> \param zeta ...
59 : !> \param rpgfa ...
60 : !> \param la_min ...
61 : !> \param lb_max ...
62 : !> \param npgfb ...
63 : !> \param zetb ...
64 : !> \param rpgfb ...
65 : !> \param lb_min ...
66 : !> \param order ...
67 : !> \param rac ...
68 : !> \param rbc ...
69 : !> \param difmab ...
70 : !> \param lambda The atom on which we take the derivative
71 : !> \param iatom ...
72 : !> \param jatom ...
73 : !> \author Edward Ditler
74 : ! **************************************************************************************************
75 27 : SUBROUTINE diff_momop_velocity(la_max, npgfa, zeta, rpgfa, la_min, &
76 27 : lb_max, npgfb, zetb, rpgfb, lb_min, &
77 27 : order, rac, rbc, difmab, lambda, iatom, jatom)
78 :
79 : INTEGER, INTENT(IN) :: la_max, npgfa
80 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
81 : INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
82 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
83 : INTEGER, INTENT(IN) :: lb_min, order
84 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
85 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(OUT) :: difmab
86 : INTEGER, INTENT(IN) :: lambda
87 : INTEGER, INTENT(IN), OPTIONAL :: iatom, jatom
88 :
89 : INTEGER :: alpha, beta, lda, lda_min, ldb, ldb_min
90 : REAL(KIND=dp) :: dab, rab(3)
91 27 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: difmab_tmp, mab
92 :
93 108 : rab = rbc - rac
94 108 : dab = SQRT(SUM(rab**2))
95 :
96 27 : lda_min = MAX(0, la_min - 1)
97 27 : ldb_min = MAX(0, lb_min - 1)
98 27 : lda = ncoset(la_max)*npgfa
99 27 : ldb = ncoset(lb_max)*npgfb
100 135 : ALLOCATE (difmab_tmp(lda, ldb, 3))
101 :
102 135 : ALLOCATE (mab(npgfa*ncoset(la_max + 1), npgfb*ncoset(lb_max + 1), ncoset(order) - 1))
103 : ! *** Calculate the primitive overlap integrals ***
104 : ! mab(1:3) = < a | r_beta - RC_beta | b >
105 27 : mab = 0.0_dp
106 : CALL moment(la_max + 1, npgfa, zeta, rpgfa, lda_min, &
107 : lb_max + 1, npgfb, zetb, rpgfb, &
108 27 : order, rac, rbc, mab)
109 :
110 17847 : difmab = 0.0_dp
111 108 : DO beta = 1, ncoset(order) - 1 ! beta was imom
112 :
113 81 : difmab_tmp = 0.0_dp
114 : CALL adbdr(la_max, npgfa, rpgfa, la_min, &
115 : lb_max, npgfb, zetb, rpgfb, lb_min, &
116 : dab, mab(:, :, beta), difmab_tmp(:, :, 1), &
117 81 : difmab_tmp(:, :, 2), difmab_tmp(:, :, 3))
118 :
119 : ! difmab(beta, alpha) = < a | r_beta - RC_beta | ∂_alpha b > * [(a==lambda) - (b==lambda)]
120 351 : DO alpha = 1, 3
121 6075 : IF (iatom == lambda) difmab(:, :, beta, alpha) = difmab(:, :, beta, alpha) + difmab_tmp(:, :, alpha)
122 6156 : IF (jatom == lambda) difmab(:, :, beta, alpha) = difmab(:, :, beta, alpha) - difmab_tmp(:, :, alpha)
123 : END DO
124 : END DO
125 :
126 27 : DEALLOCATE (mab)
127 27 : DEALLOCATE (difmab_tmp)
128 27 : END SUBROUTINE diff_momop_velocity
129 :
130 : ! *****************************************************************************
131 : !> \brief This returns the derivative of the moment integrals [a|\mu|b], with respect
132 : !> to the position of the primitive on the left and right, i.e.
133 : !> [da/dR_ai|\mu|b] + [a|\mu|d/dR_bi]
134 : !> [da/dR_ai|\mu|b] = 2*zeta*[a+1i|\mu|b] - Ni(a)[a-1i|\mu|b]
135 : !> [a|\mu|d/dR_bi] = 2*zetb*[a|\mu|b+1i] - Ni(b)[a|\mu|b-1i]
136 : !> order indicates the max order of the moment operator to be calculated
137 : !> 1: dipole
138 : !> 2: quadrupole
139 : !> ...
140 : !> \param la_max ...
141 : !> \param npgfa ...
142 : !> \param zeta ...
143 : !> \param rpgfa ...
144 : !> \param la_min ...
145 : !> \param lb_max ...
146 : !> \param npgfb ...
147 : !> \param zetb ...
148 : !> \param rpgfb ...
149 : !> \param lb_min ...
150 : !> \param order ...
151 : !> \param rac ...
152 : !> \param rbc ...
153 : !> \param difmab ...
154 : !> \param mab_ext ...
155 : !> \param deltaR needed for weighted derivative
156 : !> \param iatom ...
157 : !> \param jatom ...
158 : !> SL August 2015, ED 2021
159 : ! **************************************************************************************************
160 2268 : SUBROUTINE diff_momop2(la_max, npgfa, zeta, rpgfa, la_min, &
161 2268 : lb_max, npgfb, zetb, rpgfb, lb_min, &
162 2268 : order, rac, rbc, difmab, mab_ext, deltaR, iatom, jatom)
163 :
164 : INTEGER, INTENT(IN) :: la_max, npgfa
165 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
166 : INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
167 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
168 : INTEGER, INTENT(IN) :: lb_min, order
169 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
170 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(OUT) :: difmab
171 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
172 : POINTER :: mab_ext
173 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
174 : OPTIONAL, POINTER :: deltaR
175 : INTEGER, INTENT(IN), OPTIONAL :: iatom, jatom
176 :
177 : INTEGER :: imom, lda, lda_min, ldb, ldb_min
178 : REAL(KIND=dp) :: dab, rab(3)
179 2268 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: difmab_tmp
180 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mab
181 :
182 9072 : rab = rbc - rac
183 9072 : dab = SQRT(SUM(rab**2))
184 :
185 2268 : lda_min = MAX(0, la_min - 1)
186 2268 : ldb_min = MAX(0, lb_min - 1)
187 2268 : lda = ncoset(la_max)*npgfa
188 2268 : ldb = ncoset(lb_max)*npgfb
189 11340 : ALLOCATE (difmab_tmp(lda, ldb, 3))
190 :
191 2268 : IF (PRESENT(mab_ext)) THEN
192 0 : mab => mab_ext
193 : ELSE
194 : ALLOCATE (mab(npgfa*ncoset(la_max + 1), npgfb*ncoset(lb_max + 1), &
195 11340 : ncoset(order) - 1))
196 4091472 : mab = 0.0_dp
197 : ! *** Calculate the primitive moment integrals ***
198 : CALL moment(la_max + 1, npgfa, zeta, rpgfa, lda_min, &
199 : lb_max + 1, npgfb, zetb, rpgfb, &
200 2268 : order, rac, rbc, mab)
201 : END IF
202 9072 : DO imom = 1, ncoset(order) - 1
203 1496880 : difmab(:, :, imom, :) = 0.0_dp
204 :
205 6804 : difmab_tmp = 0.0_dp
206 : CALL adbdr(la_max, npgfa, rpgfa, la_min, &
207 : lb_max, npgfb, zetb, rpgfb, lb_min, &
208 : dab, mab(:, :, imom), difmab_tmp(:, :, 1), &
209 6804 : difmab_tmp(:, :, 2), difmab_tmp(:, :, 3))
210 :
211 496692 : difmab(:, :, imom, 1) = difmab_tmp(:, :, 1)*deltaR(1, jatom)
212 496692 : difmab(:, :, imom, 2) = difmab_tmp(:, :, 2)*deltaR(2, jatom)
213 496692 : difmab(:, :, imom, 3) = difmab_tmp(:, :, 3)*deltaR(3, jatom)
214 :
215 6804 : difmab_tmp = 0.0_dp
216 : CALL dabdr(la_max, npgfa, zeta, rpgfa, la_min, &
217 : lb_max, npgfb, rpgfb, lb_min, &
218 : dab, mab(:, :, imom), difmab_tmp(:, :, 1), &
219 6804 : difmab_tmp(:, :, 2), difmab_tmp(:, :, 3))
220 :
221 496692 : difmab(:, :, imom, 1) = difmab(:, :, imom, 1) + difmab_tmp(:, :, 1)*deltaR(1, iatom)
222 496692 : difmab(:, :, imom, 2) = difmab(:, :, imom, 2) + difmab_tmp(:, :, 2)*deltaR(2, iatom)
223 498960 : difmab(:, :, imom, 3) = difmab(:, :, imom, 3) + difmab_tmp(:, :, 3)*deltaR(3, iatom)
224 : END DO
225 :
226 2268 : IF (PRESENT(mab_ext)) THEN
227 : NULLIFY (mab)
228 : ELSE
229 2268 : DEALLOCATE (mab)
230 : END IF
231 2268 : DEALLOCATE (difmab_tmp)
232 2268 : END SUBROUTINE diff_momop2
233 :
234 : ! **************************************************************************************************
235 : !> \brief ...
236 : !> \param cos_block ...
237 : !> \param sin_block ...
238 : !> \param iatom ...
239 : !> \param ncoa ...
240 : !> \param nsgfa ...
241 : !> \param sgfa ...
242 : !> \param sphi_a ...
243 : !> \param ldsa ...
244 : !> \param jatom ...
245 : !> \param ncob ...
246 : !> \param nsgfb ...
247 : !> \param sgfb ...
248 : !> \param sphi_b ...
249 : !> \param ldsb ...
250 : !> \param cosab ...
251 : !> \param sinab ...
252 : !> \param ldab ...
253 : !> \param work ...
254 : !> \param ldwork ...
255 : ! **************************************************************************************************
256 1510708 : SUBROUTINE contract_cossin(cos_block, sin_block, &
257 3021416 : iatom, ncoa, nsgfa, sgfa, sphi_a, ldsa, &
258 3021416 : jatom, ncob, nsgfb, sgfb, sphi_b, ldsb, &
259 1510708 : cosab, sinab, ldab, work, ldwork)
260 :
261 : REAL(dp), DIMENSION(:, :), POINTER :: cos_block, sin_block
262 : INTEGER, INTENT(IN) :: iatom, ncoa, nsgfa, sgfa
263 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: sphi_a
264 : INTEGER, INTENT(IN) :: ldsa, jatom, ncob, nsgfb, sgfb
265 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: sphi_b
266 : INTEGER, INTENT(IN) :: ldsb
267 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: cosab, sinab
268 : INTEGER, INTENT(IN) :: ldab
269 : REAL(dp), DIMENSION(:, :) :: work
270 : INTEGER, INTENT(IN) :: ldwork
271 :
272 : ! Calculate cosine
273 :
274 : CALL dgemm("N", "N", ncoa, nsgfb, ncob, &
275 : 1.0_dp, cosab(1, 1), ldab, &
276 : sphi_b(1, sgfb), ldsb, &
277 1510708 : 0.0_dp, work(1, 1), ldwork)
278 :
279 1510708 : IF (iatom <= jatom) THEN
280 : CALL dgemm("T", "N", nsgfa, nsgfb, ncoa, &
281 : 1.0_dp, sphi_a(1, sgfa), ldsa, &
282 : work(1, 1), ldwork, &
283 : 1.0_dp, cos_block(sgfa, sgfb), &
284 935814 : SIZE(cos_block, 1))
285 : ELSE
286 : CALL dgemm("T", "N", nsgfb, nsgfa, ncoa, &
287 : 1.0_dp, work(1, 1), ldwork, &
288 : sphi_a(1, sgfa), ldsa, &
289 : 1.0_dp, cos_block(sgfb, sgfa), &
290 574894 : SIZE(cos_block, 1))
291 : END IF
292 :
293 : ! Calculate sine
294 : CALL dgemm("N", "N", ncoa, nsgfb, ncob, &
295 : 1.0_dp, sinab(1, 1), ldab, &
296 : sphi_b(1, sgfb), ldsb, &
297 1510708 : 0.0_dp, work(1, 1), ldwork)
298 :
299 1510708 : IF (iatom <= jatom) THEN
300 : CALL dgemm("T", "N", nsgfa, nsgfb, ncoa, &
301 : 1.0_dp, sphi_a(1, sgfa), ldsa, &
302 : work(1, 1), ldwork, &
303 : 1.0_dp, sin_block(sgfa, sgfb), &
304 935814 : SIZE(sin_block, 1))
305 : ELSE
306 : CALL dgemm("T", "N", nsgfb, nsgfa, ncoa, &
307 : 1.0_dp, work(1, 1), ldwork, &
308 : sphi_a(1, sgfa), ldsa, &
309 : 1.0_dp, sin_block(sgfb, sgfa), &
310 574894 : SIZE(sin_block, 1))
311 : END IF
312 :
313 1510708 : END SUBROUTINE contract_cossin
314 :
315 : ! **************************************************************************************************
316 : !> \brief ...
317 : !> \param la_max_set ...
318 : !> \param npgfa ...
319 : !> \param zeta ...
320 : !> \param rpgfa ...
321 : !> \param la_min_set ...
322 : !> \param lb_max ...
323 : !> \param npgfb ...
324 : !> \param zetb ...
325 : !> \param rpgfb ...
326 : !> \param lb_min ...
327 : !> \param rac ...
328 : !> \param rbc ...
329 : !> \param kvec ...
330 : !> \param cosab ...
331 : !> \param sinab ...
332 : !> \param dcosab ...
333 : !> \param dsinab ...
334 : ! **************************************************************************************************
335 1537826 : SUBROUTINE cossin(la_max_set, npgfa, zeta, rpgfa, la_min_set, &
336 1537826 : lb_max, npgfb, zetb, rpgfb, lb_min, &
337 1537826 : rac, rbc, kvec, cosab, sinab, dcosab, dsinab)
338 :
339 : INTEGER, INTENT(IN) :: la_max_set, npgfa
340 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
341 : INTEGER, INTENT(IN) :: la_min_set, lb_max, npgfb
342 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
343 : INTEGER, INTENT(IN) :: lb_min
344 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac, rbc, kvec
345 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: cosab, sinab
346 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
347 : OPTIONAL :: dcosab, dsinab
348 :
349 : INTEGER :: ax, ay, az, bx, by, bz, cda, cdax, cday, cdaz, coa, coamx, coamy, coamz, coapx, &
350 : coapy, coapz, cob, da, da_max, dax, day, daz, i, ipgf, j, jpgf, k, la, la_max, la_min, &
351 : la_start, lb, lb_start, na, nb
352 : REAL(KIND=dp) :: dab, f0, f1, f2, f3, fax, fay, faz, ftz, &
353 : fx, fy, fz, k2, kdp, rab2, s, zetp
354 : REAL(KIND=dp), DIMENSION(3) :: rab, rap, rbp
355 : REAL(KIND=dp), DIMENSION(ncoset(la_max_set), &
356 3075652 : ncoset(lb_max), 3) :: dscos, dssin
357 : REAL(KIND=dp), &
358 1537826 : DIMENSION(ncoset(la_max_set+1), ncoset(lb_max)) :: sc, ss
359 :
360 6151304 : rab = rbc - rac
361 6151304 : rab2 = SUM(rab**2)
362 1537826 : dab = SQRT(rab2)
363 1537826 : k2 = kvec(1)*kvec(1) + kvec(2)*kvec(2) + kvec(3)*kvec(3)
364 :
365 1537826 : IF (PRESENT(dcosab)) THEN
366 24916 : da_max = 1
367 24916 : la_max = la_max_set + 1
368 24916 : la_min = MAX(0, la_min_set - 1)
369 1041304 : dscos = 0.0_dp
370 1041304 : dssin = 0.0_dp
371 : ELSE
372 1512910 : da_max = 0
373 1512910 : la_max = la_max_set
374 1512910 : la_min = la_min_set
375 : END IF
376 :
377 : ! initialize all matrix elements to zero
378 1537826 : IF (PRESENT(dcosab)) THEN
379 24916 : na = ncoset(la_max - 1)*npgfa
380 : ELSE
381 1512910 : na = ncoset(la_max)*npgfa
382 : END IF
383 1537826 : nb = ncoset(lb_max)*npgfb
384 221889821 : cosab(1:na, 1:nb) = 0.0_dp
385 221889821 : sinab(1:na, 1:nb) = 0.0_dp
386 1537826 : IF (PRESENT(dcosab)) THEN
387 5627902 : dcosab(1:na, 1:nb, :) = 0.0_dp
388 5627902 : dsinab(1:na, 1:nb, :) = 0.0_dp
389 : END IF
390 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
391 :
392 1537826 : na = 0
393 5118996 : DO ipgf = 1, npgfa
394 :
395 : nb = 0
396 :
397 14971886 : DO jpgf = 1, npgfb
398 :
399 540530589 : ss = 0.0_dp
400 540530589 : sc = 0.0_dp
401 :
402 : ! *** Screening ***
403 11390716 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
404 7203732 : nb = nb + ncoset(lb_max)
405 7203732 : CYCLE
406 : END IF
407 :
408 : ! *** Calculate some prefactors ***
409 :
410 4186984 : zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
411 :
412 4186984 : f0 = (pi*zetp)**1.5_dp
413 4186984 : f1 = zetb(jpgf)*zetp
414 4186984 : f2 = 0.5_dp*zetp
415 :
416 16747936 : 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 4186984 : s = f0*EXP(-zeta(ipgf)*f1*rab2)*EXP(-0.25_dp*k2*zetp)
421 4186984 : sc(1, 1) = s*COS(kdp)
422 4186984 : ss(1, 1) = s*SIN(kdp)
423 :
424 : ! *** Recurrence steps: [s|O|s] -> [a|O|b] ***
425 :
426 4186984 : IF (la_max > 0) THEN
427 :
428 : ! *** Vertical recurrence steps: [s|O|s] -> [a|O|s] ***
429 :
430 10636492 : rap(:) = f1*rab(:)
431 :
432 : ! *** [p|O|s] = (Pi - Ai)*[s|O|s] +[s|dO|s] (i = x,y,z) ***
433 :
434 2659123 : sc(2, 1) = rap(1)*sc(1, 1) - f2*kvec(1)*ss(1, 1)
435 2659123 : sc(3, 1) = rap(2)*sc(1, 1) - f2*kvec(2)*ss(1, 1)
436 2659123 : sc(4, 1) = rap(3)*sc(1, 1) - f2*kvec(3)*ss(1, 1)
437 2659123 : ss(2, 1) = rap(1)*ss(1, 1) + f2*kvec(1)*sc(1, 1)
438 2659123 : ss(3, 1) = rap(2)*ss(1, 1) + f2*kvec(2)*sc(1, 1)
439 2659123 : 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 3181113 : 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 521990 : 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 521990 : f2*kvec(3)*sc(coset(0, 0, la - 1), 1)
454 :
455 : ! *** Increase the angular momentum component y of function a ***
456 :
457 521990 : az = la - 1
458 : sc(coset(0, 1, az), 1) = rap(2)*sc(coset(0, 0, az), 1) - &
459 521990 : 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 521990 : f2*kvec(2)*sc(coset(0, 0, az), 1)
462 :
463 1056700 : DO ay = 2, la
464 534710 : 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 534710 : 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 1056700 : 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 1578690 : DO ay = 0, la - 1
476 1056700 : az = la - 1 - ay
477 : sc(coset(1, ay, az), 1) = rap(1)*sc(coset(0, ay, az), 1) - &
478 1056700 : 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 1578690 : f2*kvec(1)*sc(coset(0, ay, az), 1)
481 : END DO
482 :
483 3715823 : DO ax = 2, la
484 534710 : f3 = f2*REAL(ax - 1, dp)
485 1604391 : DO ay = 0, la - ax
486 547691 : 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 547691 : 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 1082401 : 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 2659123 : IF (lb_max > 0) THEN
501 :
502 10419048 : DO j = 2, ncoset(lb_max)
503 58825461 : DO i = 1, ncoset(la_max)
504 48406413 : sc(i, j) = 0.0_dp
505 56800650 : ss(i, j) = 0.0_dp
506 : END DO
507 : END DO
508 :
509 : ! *** Horizontal recurrence steps ***
510 :
511 8099244 : rbp(:) = rap(:) - rab(:)
512 :
513 : ! *** [a|O|p] = [a+1i|O|s] - (Bi - Ai)*[a|O|s] ***
514 :
515 2024811 : IF (lb_max == 1) THEN
516 : la_start = la_min
517 : ELSE
518 378054 : la_start = MAX(0, la_min - 1)
519 : END IF
520 :
521 4022225 : DO la = la_start, la_max - 1
522 6289549 : DO ax = 0, la
523 6805140 : DO ay = 0, la - ax
524 2540402 : az = la - ax - ay
525 : sc(coset(ax, ay, az), 2) = sc(coset(ax + 1, ay, az), 1) - &
526 2540402 : rab(1)*sc(coset(ax, ay, az), 1)
527 : sc(coset(ax, ay, az), 3) = sc(coset(ax, ay + 1, az), 1) - &
528 2540402 : rab(2)*sc(coset(ax, ay, az), 1)
529 : sc(coset(ax, ay, az), 4) = sc(coset(ax, ay, az + 1), 1) - &
530 2540402 : rab(3)*sc(coset(ax, ay, az), 1)
531 : ss(coset(ax, ay, az), 2) = ss(coset(ax + 1, ay, az), 1) - &
532 2540402 : rab(1)*ss(coset(ax, ay, az), 1)
533 : ss(coset(ax, ay, az), 3) = ss(coset(ax, ay + 1, az), 1) - &
534 2540402 : rab(2)*ss(coset(ax, ay, az), 1)
535 : ss(coset(ax, ay, az), 4) = ss(coset(ax, ay, az + 1), 1) - &
536 4807726 : 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 6471254 : DO ax = 0, la_max
547 4446443 : fx = f2*REAL(ax, dp)
548 13743188 : DO ay = 0, la_max - ax
549 7271934 : fy = f2*REAL(ay, dp)
550 7271934 : az = la_max - ax - ay
551 7271934 : fz = f2*REAL(az, dp)
552 7271934 : IF (ax == 0) THEN
553 : sc(coset(ax, ay, az), 2) = rbp(1)*sc(coset(ax, ay, az), 1) - &
554 4446443 : 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 4446443 : 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 2825491 : 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 2825491 : f2*kvec(1)*sc(coset(ax, ay, az), 1)
564 : END IF
565 7271934 : IF (ay == 0) THEN
566 : sc(coset(ax, ay, az), 3) = rbp(2)*sc(coset(ax, ay, az), 1) - &
567 4446443 : 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 4446443 : 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 2825491 : 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 2825491 : f2*kvec(2)*sc(coset(ax, ay, az), 1)
577 : END IF
578 11718377 : IF (az == 0) THEN
579 : sc(coset(ax, ay, az), 4) = rbp(3)*sc(coset(ax, ay, az), 1) - &
580 4446443 : 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 4446443 : 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 2825491 : 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 2825491 : 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 2407950 : 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 383139 : IF (lb == lb_max) THEN
603 : la_start = la_min
604 : ELSE
605 5085 : la_start = MAX(0, la_min - 1)
606 : END IF
607 :
608 828453 : DO la = la_start, la_max - 1
609 1430735 : DO ax = 0, la
610 1807446 : DO ay = 0, la - ax
611 759850 : 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 759850 : 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 759850 : 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 2282229 : DO by = 1, lb
625 1522379 : 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 1522379 : 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 2282229 : 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 2884511 : DO bx = 1, lb
637 4569816 : DO by = 0, lb - bx
638 2287587 : 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 2287587 : 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 3809966 : 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 3382617 : DO ax = 0, la_max
658 974667 : fx = f2*REAL(ax, dp)
659 3134352 : DO ay = 0, la_max - ax
660 1776546 : fy = f2*REAL(ay, dp)
661 1776546 : az = la_max - ax - ay
662 1776546 : fz = f2*REAL(az, dp)
663 :
664 : ! *** Increase the angular momentum component z of function b ***
665 :
666 1776546 : f3 = f2*REAL(lb - 1, dp)
667 :
668 1776546 : 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 974667 : 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 974667 : 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 801879 : 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 801879 : 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 1776546 : IF (ay == 0) THEN
693 974667 : 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 974667 : 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 974667 : f2*kvec(2)*sc(coset(ax, ay, az), coset(0, 0, bz))
700 1961568 : DO by = 2, lb
701 986901 : bz = lb - by
702 986901 : 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 986901 : 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 1961568 : f2*kvec(2)*sc(coset(ax, ay, az), coset(0, by - 1, bz))
711 : END DO
712 : ELSE
713 801879 : 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 801879 : 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 801879 : f2*kvec(2)*sc(coset(ax, ay, az), coset(0, 0, bz))
722 1613082 : DO by = 2, lb
723 811203 : bz = lb - by
724 811203 : 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 811203 : 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 1613082 : 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 2751213 : IF (ax == 0) THEN
741 2936235 : DO by = 0, lb - 1
742 1961568 : 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 1961568 : 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 2936235 : f2*kvec(1)*sc(coset(ax, ay, az), coset(0, by, bz))
749 : END DO
750 1961568 : DO bx = 2, lb
751 986901 : f3 = f2*REAL(bx - 1, dp)
752 2961045 : DO by = 0, lb - bx
753 999477 : 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 999477 : 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 1986378 : f2*kvec(1)*sc(coset(ax, ay, az), coset(bx - 1, by, bz))
764 : END DO
765 : END DO
766 : ELSE
767 2414961 : DO by = 0, lb - 1
768 1613082 : 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 1613082 : 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 2414961 : f2*kvec(1)*sc(coset(ax, ay, az), coset(0, by, bz))
777 : END DO
778 1613082 : DO bx = 2, lb
779 811203 : f3 = f2*REAL(bx - 1, dp)
780 2433960 : DO by = 0, lb - bx
781 820878 : 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 820878 : 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 1632081 : 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 1527861 : IF (lb_max > 0) THEN
808 :
809 : ! *** Vertical recurrence steps: [s|O|s] -> [s|O|b] ***
810 :
811 2317160 : rbp(:) = (f1 - 1.0_dp)*rab(:)
812 :
813 : ! *** [s|O|p] = (Pi - Bi)*[s|O|s] + [s|dO|s] ***
814 :
815 579290 : sc(1, 2) = rbp(1)*sc(1, 1) - f2*kvec(1)*ss(1, 1)
816 579290 : sc(1, 3) = rbp(2)*sc(1, 1) - f2*kvec(2)*ss(1, 1)
817 579290 : sc(1, 4) = rbp(3)*sc(1, 1) - f2*kvec(3)*ss(1, 1)
818 579290 : ss(1, 2) = rbp(1)*ss(1, 1) + f2*kvec(1)*sc(1, 1)
819 579290 : ss(1, 3) = rbp(2)*ss(1, 1) + f2*kvec(2)*sc(1, 1)
820 579290 : 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 684369 : 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 105079 : 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 105079 : f2*kvec(3)*sc(1, coset(0, 0, lb - 1))
835 :
836 : ! *** Increase the angular momentum component y of function b ***
837 :
838 105079 : bz = lb - 1
839 : sc(1, coset(0, 1, bz)) = rbp(2)*sc(1, coset(0, 0, bz)) - &
840 105079 : 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 105079 : f2*kvec(2)*sc(1, coset(0, 0, bz))
843 :
844 214340 : DO by = 2, lb
845 109261 : 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 109261 : 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 214340 : 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 319419 : DO by = 0, lb - 1
857 214340 : bz = lb - 1 - by
858 : sc(1, coset(1, by, bz)) = rbp(1)*sc(1, coset(0, by, bz)) - &
859 214340 : 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 319419 : f2*kvec(1)*sc(1, coset(0, by, bz))
862 : END DO
863 :
864 793630 : DO bx = 2, lb
865 109261 : f3 = f2*REAL(bx - 1, dp)
866 327918 : DO by = 0, lb - bx
867 113578 : 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 113578 : 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 222839 : 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 17599063 : DO j = ncoset(lb_min - 1) + 1, ncoset(lb_max)
884 72570951 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
885 54971888 : cosab(na + i, nb + j) = sc(i, j)
886 68383967 : sinab(na + i, nb + j) = ss(i, j)
887 : END DO
888 : END DO
889 :
890 4186984 : IF (PRESENT(dcosab)) THEN
891 : la_start = 0
892 : lb_start = 0
893 : ELSE
894 4125740 : la_start = la_min
895 4125740 : lb_start = lb_min
896 : END IF
897 :
898 4248228 : DO da = 0, da_max - 1
899 61244 : ftz = 2.0_dp*zeta(ipgf)
900 4309472 : DO dax = 0, da
901 183732 : DO day = 0, da - dax
902 61244 : daz = da - dax - day
903 61244 : cda = coset(dax, day, daz) - 1
904 61244 : cdax = coset(dax + 1, day, daz) - 1
905 61244 : cday = coset(dax, day + 1, daz) - 1
906 61244 : 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 213755 : DO la = la_start, la_max - da - 1
910 276858 : DO ax = 0, la
911 124347 : fax = REAL(ax, dp)
912 376098 : DO ay = 0, la - ax
913 160484 : fay = REAL(ay, dp)
914 160484 : az = la - ax - ay
915 160484 : faz = REAL(az, dp)
916 160484 : coa = coset(ax, ay, az)
917 160484 : coamx = coset(ax - 1, ay, az)
918 160484 : coamy = coset(ax, ay - 1, az)
919 160484 : coamz = coset(ax, ay, az - 1)
920 160484 : coapx = coset(ax + 1, ay, az)
921 160484 : coapy = coset(ax, ay + 1, az)
922 160484 : coapz = coset(ax, ay, az + 1)
923 533290 : DO lb = lb_start, lb_max
924 754566 : DO bx = 0, lb
925 1046058 : DO by = 0, lb - bx
926 451976 : bz = lb - bx - by
927 451976 : cob = coset(bx, by, bz)
928 451976 : dscos(coa, cob, cdax) = ftz*sc(coapx, cob) - fax*sc(coamx, cob)
929 451976 : dscos(coa, cob, cday) = ftz*sc(coapy, cob) - fay*sc(coamy, cob)
930 451976 : dscos(coa, cob, cdaz) = ftz*sc(coapz, cob) - faz*sc(coamz, cob)
931 451976 : dssin(coa, cob, cdax) = ftz*ss(coapx, cob) - fax*ss(coamx, cob)
932 451976 : dssin(coa, cob, cday) = ftz*ss(coapy, cob) - fay*ss(coamy, cob)
933 797599 : 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 4186984 : IF (PRESENT(dcosab)) THEN
946 244976 : DO k = 1, 3
947 715556 : DO j = 1, ncoset(lb_max)
948 2010240 : DO i = 1, ncoset(la_max_set)
949 1355928 : dcosab(na + i, nb + j, k) = dscos(i, j, k)
950 1826508 : dsinab(na + i, nb + j, k) = dssin(i, j, k)
951 : END DO
952 : END DO
953 : END DO
954 : END IF
955 :
956 7768154 : nb = nb + ncoset(lb_max)
957 :
958 : END DO
959 :
960 5118996 : na = na + ncoset(la_max_set)
961 :
962 : END DO
963 :
964 1537826 : 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 648766 : SUBROUTINE moment(la_max, npgfa, zeta, rpgfa, la_min, &
983 1297532 : lb_max, npgfb, zetb, rpgfb, &
984 648766 : 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(&
1001 648766 : lb_max), ncoset(lc_max)) :: s
1002 :
1003 2595064 : rab = rbc - rac
1004 2595064 : rab2 = SUM(rab**2)
1005 648766 : dab = SQRT(rab2)
1006 :
1007 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
1008 :
1009 648766 : na = 0
1010 :
1011 2118970 : DO ipgf = 1, npgfa
1012 :
1013 1470204 : nb = 0
1014 :
1015 5448299 : DO jpgf = 1, npgfb
1016 :
1017 2903713515 : s = 0.0_dp
1018 : ! *** Screening ***
1019 :
1020 3978095 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
1021 24081817 : DO k = 1, ncoset(lc_max) - 1
1022 195504314 : DO j = nb + 1, nb + ncoset(lb_max)
1023 1713647151 : DO i = na + 1, na + ncoset(la_max)
1024 1692033863 : mab(i, j, k) = 0.0_dp
1025 : END DO
1026 : END DO
1027 : END DO
1028 2468529 : nb = nb + ncoset(lb_max)
1029 2468529 : CYCLE
1030 : END IF
1031 :
1032 : ! *** Calculate some prefactors ***
1033 :
1034 1509566 : zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
1035 :
1036 1509566 : f0 = (pi*zetp)**1.5_dp
1037 1509566 : f1 = zetb(jpgf)*zetp
1038 1509566 : f2 = 0.5_dp*zetp
1039 :
1040 : ! *** Calculate the basic two-center moment integral [s|M|s] ***
1041 :
1042 6038264 : rpc = zetp*(zeta(ipgf)*rac + zetb(jpgf)*rbc)
1043 1509566 : s(1, 1, 1) = f0*EXP(-zeta(ipgf)*f1*rab2)
1044 12898021 : DO l = 2, ncoset(lc_max)
1045 11388455 : lx = indco(1, l)
1046 11388455 : ly = indco(2, l)
1047 11388455 : lz = indco(3, l)
1048 11388455 : l2 = 0
1049 11388455 : IF (lz > 0) THEN
1050 4939966 : l1 = coset(lx, ly, lz - 1)
1051 4939966 : IF (lz > 1) l2 = coset(lx, ly, lz - 2)
1052 : ni = lz - 1
1053 : i = 3
1054 6448489 : ELSE IF (ly > 0) THEN
1055 3795953 : l1 = coset(lx, ly - 1, lz)
1056 3795953 : IF (ly > 1) l2 = coset(lx, ly - 2, lz)
1057 : ni = ly - 1
1058 : i = 2
1059 2652536 : ELSE IF (lx > 0) THEN
1060 2652536 : l1 = coset(lx - 1, ly, lz)
1061 2652536 : IF (lx > 1) l2 = coset(lx - 2, ly, lz)
1062 : ni = lx - 1
1063 : i = 1
1064 : END IF
1065 11388455 : s(1, 1, l) = rpc(i)*s(1, 1, l1)
1066 12898021 : IF (l2 > 0) s(1, 1, l) = s(1, 1, l) + f2*REAL(ni, dp)*s(1, 1, l2)
1067 : END DO
1068 :
1069 : ! *** Recurrence steps: [s|M|s] -> [a|M|b] ***
1070 :
1071 14407587 : DO l = 1, ncoset(lc_max)
1072 :
1073 12898021 : lx = indco(1, l)
1074 12898021 : ly = indco(2, l)
1075 12898021 : lz = indco(3, l)
1076 12898021 : IF (lx > 0) THEN
1077 4939966 : lx1 = coset(lx - 1, ly, lz)
1078 : ELSE
1079 : lx1 = -1
1080 : END IF
1081 12898021 : IF (ly > 0) THEN
1082 4939966 : ly1 = coset(lx, ly - 1, lz)
1083 : ELSE
1084 : ly1 = -1
1085 : END IF
1086 12898021 : IF (lz > 0) THEN
1087 4939966 : lz1 = coset(lx, ly, lz - 1)
1088 : ELSE
1089 : lz1 = -1
1090 : END IF
1091 12898021 : f2x = f2*REAL(lx, dp)
1092 12898021 : f2y = f2*REAL(ly, dp)
1093 12898021 : f2z = f2*REAL(lz, dp)
1094 :
1095 14407587 : IF (la_max > 0) THEN
1096 :
1097 : ! *** Vertical recurrence steps: [s|M|s] -> [a|M|s] ***
1098 :
1099 49121872 : rap(:) = f1*rab(:)
1100 :
1101 : ! *** [p|M|s] = (Pi - Ai)*[s|M|s] + f2*Ni(m-1i)[s|M-1i|s] ***
1102 :
1103 12280468 : s(2, 1, l) = rap(1)*s(1, 1, l)
1104 12280468 : s(3, 1, l) = rap(2)*s(1, 1, l)
1105 12280468 : s(4, 1, l) = rap(3)*s(1, 1, l)
1106 12280468 : IF (lx1 > 0) s(2, 1, l) = s(2, 1, l) + f2x*s(1, 1, lx1)
1107 12280468 : IF (ly1 > 0) s(3, 1, l) = s(3, 1, l) + f2y*s(1, 1, ly1)
1108 12280468 : IF (lz1 > 0) s(4, 1, l) = s(4, 1, l) + f2z*s(1, 1, lz1)
1109 :
1110 : ! *** [a|M|s] = (Pi - Ai)*[a-1i|M|s] + f2*Ni(a-1i)*[a-2i|M|s] ***
1111 : ! *** + f2*Ni(m-1i)*[a-1i|M-1i|s] ***
1112 :
1113 19967260 : DO la = 2, la_max
1114 :
1115 : ! *** Increase the angular momentum component z of function a ***
1116 :
1117 : s(coset(0, 0, la), 1, l) = rap(3)*s(coset(0, 0, la - 1), 1, l) + &
1118 7686792 : f2*REAL(la - 1, dp)*s(coset(0, 0, la - 2), 1, l)
1119 7686792 : IF (lz1 > 0) s(coset(0, 0, la), 1, l) = s(coset(0, 0, la), 1, l) + &
1120 3044433 : f2z*s(coset(0, 0, la - 1), 1, lz1)
1121 :
1122 : ! *** Increase the angular momentum component y of function a ***
1123 :
1124 7686792 : az = la - 1
1125 7686792 : s(coset(0, 1, az), 1, l) = rap(2)*s(coset(0, 0, az), 1, l)
1126 7686792 : IF (ly1 > 0) s(coset(0, 1, az), 1, l) = s(coset(0, 1, az), 1, l) + &
1127 3044433 : f2y*s(coset(0, 0, az), 1, ly1)
1128 :
1129 16205598 : DO ay = 2, la
1130 8518806 : az = la - ay
1131 : s(coset(0, ay, az), 1, l) = rap(2)*s(coset(0, ay - 1, az), 1, l) + &
1132 8518806 : f2*REAL(ay - 1, dp)*s(coset(0, ay - 2, az), 1, l)
1133 8518806 : IF (ly1 > 0) s(coset(0, ay, az), 1, l) = s(coset(0, ay, az), 1, l) + &
1134 11061912 : f2y*s(coset(0, ay - 1, az), 1, ly1)
1135 : END DO
1136 :
1137 : ! *** Increase the angular momentum component x of function a ***
1138 :
1139 23892390 : DO ay = 0, la - 1
1140 16205598 : az = la - 1 - ay
1141 16205598 : s(coset(1, ay, az), 1, l) = rap(1)*s(coset(0, ay, az), 1, l)
1142 16205598 : IF (lx1 > 0) s(coset(1, ay, az), 1, l) = s(coset(1, ay, az), 1, l) + &
1143 14106345 : f2x*s(coset(0, ay, az), 1, lx1)
1144 : END DO
1145 :
1146 28486066 : DO ax = 2, la
1147 8518806 : f3 = f2*REAL(ax - 1, dp)
1148 25556418 : DO ay = 0, la - ax
1149 9350820 : az = la - ax - ay
1150 : s(coset(ax, ay, az), 1, l) = rap(1)*s(coset(ax - 1, ay, az), 1, l) + &
1151 9350820 : f3*s(coset(ax - 2, ay, az), 1, l)
1152 9350820 : IF (lx1 > 0) s(coset(ax, ay, az), 1, l) = s(coset(ax, ay, az), 1, l) + &
1153 12224613 : f2x*s(coset(ax - 1, ay, az), 1, lx1)
1154 : END DO
1155 : END DO
1156 :
1157 : END DO
1158 :
1159 : ! *** Recurrence steps: [a|M|s] -> [a|M|b] ***
1160 :
1161 12280468 : IF (lb_max > 0) THEN
1162 :
1163 97107992 : DO j = 2, ncoset(lb_max)
1164 877822304 : DO i = 1, ncoset(la_max)
1165 865859864 : s(i, j, l) = 0.0_dp
1166 : END DO
1167 : END DO
1168 :
1169 : ! *** Horizontal recurrence steps ***
1170 :
1171 47849760 : rbp(:) = rap(:) - rab(:)
1172 :
1173 : ! *** [a|M|p] = [a+1i|M|s] - (Bi - Ai)*[a|M|s] ***
1174 :
1175 11962440 : IF (lb_max == 1) THEN
1176 5138078 : la_start = la_min
1177 : ELSE
1178 6824362 : la_start = MAX(0, la_min - 1)
1179 : END IF
1180 :
1181 31396232 : DO la = la_start, la_max - 1
1182 59288495 : DO ax = 0, la
1183 84507987 : DO ay = 0, la - ax
1184 37181932 : az = la - ax - ay
1185 : s(coset(ax, ay, az), 2, l) = s(coset(ax + 1, ay, az), 1, l) - &
1186 37181932 : rab(1)*s(coset(ax, ay, az), 1, l)
1187 : s(coset(ax, ay, az), 3, l) = s(coset(ax, ay + 1, az), 1, l) - &
1188 37181932 : rab(2)*s(coset(ax, ay, az), 1, l)
1189 : s(coset(ax, ay, az), 4, l) = s(coset(ax, ay, az + 1), 1, l) - &
1190 65074195 : rab(3)*s(coset(ax, ay, az), 1, l)
1191 : END DO
1192 : END DO
1193 : END DO
1194 :
1195 : ! *** Vertical recurrence step ***
1196 :
1197 : ! *** [a|M|p] = (Pi - Bi)*[a|M|s] + f2*Ni(a)*[a-1i|M|s] ***
1198 : ! *** + f2*Ni(m)*[a|M-1i|s] ***
1199 :
1200 43546082 : DO ax = 0, la_max
1201 31583642 : fx = f2*REAL(ax, dp)
1202 103241174 : DO ay = 0, la_max - ax
1203 59695092 : fy = f2*REAL(ay, dp)
1204 59695092 : az = la_max - ax - ay
1205 59695092 : fz = f2*REAL(az, dp)
1206 59695092 : IF (ax == 0) THEN
1207 31583642 : s(coset(ax, ay, az), 2, l) = rbp(1)*s(coset(ax, ay, az), 1, l)
1208 : ELSE
1209 : s(coset(ax, ay, az), 2, l) = rbp(1)*s(coset(ax, ay, az), 1, l) + &
1210 28111450 : fx*s(coset(ax - 1, ay, az), 1, l)
1211 : END IF
1212 59695092 : IF (lx1 > 0) s(coset(ax, ay, az), 2, l) = s(coset(ax, ay, az), 2, l) + &
1213 23483727 : f2x*s(coset(ax, ay, az), 1, lx1)
1214 59695092 : IF (ay == 0) THEN
1215 31583642 : s(coset(ax, ay, az), 3, l) = rbp(2)*s(coset(ax, ay, az), 1, l)
1216 : ELSE
1217 : s(coset(ax, ay, az), 3, l) = rbp(2)*s(coset(ax, ay, az), 1, l) + &
1218 28111450 : fy*s(coset(ax, ay - 1, az), 1, l)
1219 : END IF
1220 59695092 : IF (ly1 > 0) s(coset(ax, ay, az), 3, l) = s(coset(ax, ay, az), 3, l) + &
1221 23483727 : f2y*s(coset(ax, ay, az), 1, ly1)
1222 59695092 : IF (az == 0) THEN
1223 31583642 : s(coset(ax, ay, az), 4, l) = rbp(3)*s(coset(ax, ay, az), 1, l)
1224 : ELSE
1225 : s(coset(ax, ay, az), 4, l) = rbp(3)*s(coset(ax, ay, az), 1, l) + &
1226 28111450 : fz*s(coset(ax, ay, az - 1), 1, l)
1227 : END IF
1228 59695092 : IF (lz1 > 0) s(coset(ax, ay, az), 4, l) = s(coset(ax, ay, az), 4, l) + &
1229 55067369 : f2z*s(coset(ax, ay, az), 1, lz1)
1230 : END DO
1231 : END DO
1232 :
1233 : ! *** Recurrence steps: [a|M|p] -> [a|M|b] ***
1234 :
1235 19618008 : DO lb = 2, lb_max
1236 :
1237 : ! *** Horizontal recurrence steps ***
1238 :
1239 : ! *** [a|M|b] = [a+1i|M|b-1i] - (Bi - Ai)*[a|M|b-1i] ***
1240 :
1241 7655568 : IF (lb == lb_max) THEN
1242 6824362 : la_start = la_min
1243 : ELSE
1244 831206 : la_start = MAX(0, la_min - 1)
1245 : END IF
1246 :
1247 21558402 : DO la = la_start, la_max - 1
1248 43273052 : DO ax = 0, la
1249 65957684 : DO ay = 0, la - ax
1250 30340200 : az = la - ax - ay
1251 :
1252 : ! *** Shift of angular momentum component z from a to b ***
1253 :
1254 : s(coset(ax, ay, az), coset(0, 0, lb), l) = &
1255 : s(coset(ax, ay, az + 1), coset(0, 0, lb - 1), l) - &
1256 30340200 : rab(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1), l)
1257 :
1258 : ! *** Shift of angular momentum component y from a to b ***
1259 :
1260 94427978 : DO by = 1, lb
1261 64087778 : bz = lb - by
1262 : s(coset(ax, ay, az), coset(0, by, bz), l) = &
1263 : s(coset(ax, ay + 1, az), coset(0, by - 1, bz), l) - &
1264 94427978 : rab(2)*s(coset(ax, ay, az), coset(0, by - 1, bz), l)
1265 : END DO
1266 :
1267 : ! *** Shift of angular momentum component x from a to b ***
1268 :
1269 116142628 : DO bx = 1, lb
1270 195670712 : DO by = 0, lb - bx
1271 101242734 : bz = lb - bx - by
1272 : s(coset(ax, ay, az), coset(bx, by, bz), l) = &
1273 : s(coset(ax + 1, ay, az), coset(bx - 1, by, bz), l) - &
1274 165330512 : rab(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz), l)
1275 : END DO
1276 : END DO
1277 :
1278 : END DO
1279 : END DO
1280 : END DO
1281 :
1282 : ! *** Vertical recurrence step ***
1283 :
1284 : ! *** [a|M|b] = (Pi - Bi)*[a|M|b-1i] + f2*Ni(a)*[a-1i|M|b-1i] + ***
1285 : ! *** f2*Ni(b-1i)*[a|M|b-2i] + f2*Ni(m)[a|M-1i|b-1i] ***
1286 :
1287 41934331 : DO ax = 0, la_max
1288 22316323 : fx = f2*REAL(ax, dp)
1289 74768034 : DO ay = 0, la_max - ax
1290 44796143 : fy = f2*REAL(ay, dp)
1291 44796143 : az = la_max - ax - ay
1292 44796143 : fz = f2*REAL(az, dp)
1293 :
1294 : ! *** Shift of angular momentum component z from a to b ***
1295 :
1296 44796143 : f3 = f2*REAL(lb - 1, dp)
1297 :
1298 44796143 : IF (az == 0) THEN
1299 : s(coset(ax, ay, az), coset(0, 0, lb), l) = &
1300 : rbp(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1), l) + &
1301 22316323 : f3*s(coset(ax, ay, az), coset(0, 0, lb - 2), l)
1302 : ELSE
1303 : s(coset(ax, ay, az), coset(0, 0, lb), l) = &
1304 : rbp(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1), l) + &
1305 : fz*s(coset(ax, ay, az - 1), coset(0, 0, lb - 1), l) + &
1306 22479820 : f3*s(coset(ax, ay, az), coset(0, 0, lb - 2), l)
1307 : END IF
1308 44796143 : IF (lz1 > 0) s(coset(ax, ay, az), coset(0, 0, lb), l) = &
1309 : s(coset(ax, ay, az), coset(0, 0, lb), l) + &
1310 17793194 : f2z*s(coset(ax, ay, az), coset(0, 0, lb - 1), lz1)
1311 :
1312 : ! *** Shift of angular momentum component y from a to b ***
1313 :
1314 44796143 : IF (ay == 0) THEN
1315 22316323 : bz = lb - 1
1316 : s(coset(ax, ay, az), coset(0, 1, bz), l) = &
1317 22316323 : rbp(2)*s(coset(ax, ay, az), coset(0, 0, bz), l)
1318 22316323 : IF (ly1 > 0) s(coset(ax, ay, az), coset(0, 1, bz), l) = &
1319 : s(coset(ax, ay, az), coset(0, 1, bz), l) + &
1320 8859553 : f2y*s(coset(ax, ay, az), coset(0, 0, bz), ly1)
1321 47108844 : DO by = 2, lb
1322 24792521 : bz = lb - by
1323 24792521 : f3 = f2*REAL(by - 1, dp)
1324 : s(coset(ax, ay, az), coset(0, by, bz), l) = &
1325 : rbp(2)*s(coset(ax, ay, az), coset(0, by - 1, bz), l) + &
1326 24792521 : f3*s(coset(ax, ay, az), coset(0, by - 2, bz), l)
1327 24792521 : IF (ly1 > 0) s(coset(ax, ay, az), coset(0, by, bz), l) = &
1328 : s(coset(ax, ay, az), coset(0, by, bz), l) + &
1329 32160024 : f2y*s(coset(ax, ay, az), coset(0, by - 1, bz), ly1)
1330 : END DO
1331 : ELSE
1332 22479820 : bz = lb - 1
1333 : s(coset(ax, ay, az), coset(0, 1, bz), l) = &
1334 : rbp(2)*s(coset(ax, ay, az), coset(0, 0, bz), l) + &
1335 22479820 : fy*s(coset(ax, ay - 1, az), coset(0, 0, bz), l)
1336 22479820 : IF (ly1 > 0) s(coset(ax, ay, az), coset(0, 1, bz), l) = &
1337 : s(coset(ax, ay, az), coset(0, 1, bz), l) + &
1338 8933641 : f2y*s(coset(ax, ay, az), coset(0, 0, bz), ly1)
1339 47483786 : DO by = 2, lb
1340 25003966 : bz = lb - by
1341 25003966 : f3 = f2*REAL(by - 1, dp)
1342 : s(coset(ax, ay, az), coset(0, by, bz), l) = &
1343 : rbp(2)*s(coset(ax, ay, az), coset(0, by - 1, bz), l) + &
1344 : fy*s(coset(ax, ay - 1, az), coset(0, by - 1, bz), l) + &
1345 25003966 : f3*s(coset(ax, ay, az), coset(0, by - 2, bz), l)
1346 25003966 : IF (ly1 > 0) s(coset(ax, ay, az), coset(0, by, bz), l) = &
1347 : s(coset(ax, ay, az), coset(0, by, bz), l) + &
1348 32415815 : f2y*s(coset(ax, ay, az), coset(0, by - 1, bz), ly1)
1349 : END DO
1350 : END IF
1351 :
1352 : ! *** Shift of angular momentum component x from a to b ***
1353 :
1354 67112466 : IF (ax == 0) THEN
1355 69425167 : DO by = 0, lb - 1
1356 47108844 : bz = lb - 1 - by
1357 : s(coset(ax, ay, az), coset(1, by, bz), l) = &
1358 47108844 : rbp(1)*s(coset(ax, ay, az), coset(0, by, bz), l)
1359 47108844 : IF (lx1 > 0) s(coset(ax, ay, az), coset(1, by, bz), l) = &
1360 : s(coset(ax, ay, az), coset(1, by, bz), l) + &
1361 41019577 : f2x*s(coset(ax, ay, az), coset(0, by, bz), lx1)
1362 : END DO
1363 47108844 : DO bx = 2, lb
1364 24792521 : f3 = f2*REAL(bx - 1, dp)
1365 74377563 : DO by = 0, lb - bx
1366 27268719 : bz = lb - bx - by
1367 : s(coset(ax, ay, az), coset(bx, by, bz), l) = &
1368 : rbp(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz), l) + &
1369 27268719 : f3*s(coset(ax, ay, az), coset(bx - 2, by, bz), l)
1370 27268719 : IF (lx1 > 0) s(coset(ax, ay, az), coset(bx, by, bz), l) = &
1371 : s(coset(ax, ay, az), coset(bx, by, bz), l) + &
1372 35620370 : f2x*s(coset(ax, ay, az), coset(bx - 1, by, bz), lx1)
1373 : END DO
1374 : END DO
1375 : ELSE
1376 69963606 : DO by = 0, lb - 1
1377 47483786 : bz = lb - 1 - by
1378 : s(coset(ax, ay, az), coset(1, by, bz), l) = &
1379 : rbp(1)*s(coset(ax, ay, az), coset(0, by, bz), l) + &
1380 47483786 : fx*s(coset(ax - 1, ay, az), coset(0, by, bz), l)
1381 47483786 : IF (lx1 > 0) s(coset(ax, ay, az), coset(1, by, bz), l) = &
1382 : s(coset(ax, ay, az), coset(1, by, bz), l) + &
1383 41349456 : f2x*s(coset(ax, ay, az), coset(0, by, bz), lx1)
1384 : END DO
1385 47483786 : DO bx = 2, lb
1386 25003966 : f3 = f2*REAL(bx - 1, dp)
1387 75011898 : DO by = 0, lb - bx
1388 27528112 : bz = lb - bx - by
1389 : s(coset(ax, ay, az), coset(bx, by, bz), l) = &
1390 : rbp(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz), l) + &
1391 : fx*s(coset(ax - 1, ay, az), coset(bx - 1, by, bz), l) + &
1392 27528112 : f3*s(coset(ax, ay, az), coset(bx - 2, by, bz), l)
1393 27528112 : IF (lx1 > 0) s(coset(ax, ay, az), coset(bx, by, bz), l) = &
1394 : s(coset(ax, ay, az), coset(bx, by, bz), l) + &
1395 35942315 : f2x*s(coset(ax, ay, az), coset(bx - 1, by, bz), lx1)
1396 : END DO
1397 : END DO
1398 : END IF
1399 :
1400 : END DO
1401 : END DO
1402 :
1403 : END DO
1404 :
1405 : END IF
1406 :
1407 : ELSE
1408 :
1409 617553 : IF (lb_max > 0) THEN
1410 :
1411 : ! *** Vertical recurrence steps: [s|M|s] -> [s|M|b] ***
1412 :
1413 842448 : rbp(:) = (f1 - 1.0_dp)*rab(:)
1414 :
1415 : ! *** [s|M|p] = (Pi - Bi)*[s|M|s] + f2*Ni(m)*[s|M-1i|s] ***
1416 :
1417 210612 : s(1, 2, l) = rbp(1)*s(1, 1, l)
1418 210612 : s(1, 3, l) = rbp(2)*s(1, 1, l)
1419 210612 : s(1, 4, l) = rbp(3)*s(1, 1, l)
1420 210612 : IF (lx1 > 0) s(1, 2, l) = s(1, 2, l) + f2x*s(1, 1, lx1)
1421 210612 : IF (ly1 > 0) s(1, 3, l) = s(1, 3, l) + f2y*s(1, 1, ly1)
1422 210612 : IF (lz1 > 0) s(1, 4, l) = s(1, 4, l) + f2z*s(1, 1, lz1)
1423 :
1424 : ! *** [s|M|b] = (Pi - Bi)*[s|M|b-1i] + f2*Ni(b-1i)*[s|M|b-2i] ***
1425 : ! *** + f2*Ni(m)*[s|M-1i|b-1i] ***
1426 :
1427 237224 : DO lb = 2, lb_max
1428 :
1429 : ! *** Increase the angular momentum component z of function b ***
1430 :
1431 : s(1, coset(0, 0, lb), l) = rbp(3)*s(1, coset(0, 0, lb - 1), l) + &
1432 26612 : f2*REAL(lb - 1, dp)*s(1, coset(0, 0, lb - 2), l)
1433 26612 : IF (lz1 > 0) s(1, coset(0, 0, lb), l) = s(1, coset(0, 0, lb), l) + &
1434 6809 : f2z*s(1, coset(0, 0, lb - 1), lz1)
1435 :
1436 : ! *** Increase the angular momentum component y of function b ***
1437 :
1438 26612 : bz = lb - 1
1439 26612 : s(1, coset(0, 1, bz), l) = rbp(2)*s(1, coset(0, 0, bz), l)
1440 26612 : IF (ly1 > 0) s(1, coset(0, 1, bz), l) = s(1, coset(0, 1, bz), l) + &
1441 6809 : f2y*s(1, coset(0, 0, bz), ly1)
1442 :
1443 53752 : DO by = 2, lb
1444 27140 : bz = lb - by
1445 : s(1, coset(0, by, bz), l) = rbp(2)*s(1, coset(0, by - 1, bz), l) + &
1446 27140 : f2*REAL(by - 1, dp)*s(1, coset(0, by - 2, bz), l)
1447 27140 : IF (ly1 > 0) s(1, coset(0, by, bz), l) = s(1, coset(0, by, bz), l) + &
1448 33553 : f2y*s(1, coset(0, by - 1, bz), ly1)
1449 : END DO
1450 :
1451 : ! *** Increase the angular momentum component x of function b ***
1452 :
1453 80364 : DO by = 0, lb - 1
1454 53752 : bz = lb - 1 - by
1455 53752 : s(1, coset(1, by, bz), l) = rbp(1)*s(1, coset(0, by, bz), l)
1456 53752 : IF (lx1 > 0) s(1, coset(1, by, bz), l) = s(1, coset(1, by, bz), l) + &
1457 40362 : f2x*s(1, coset(0, by, bz), lx1)
1458 : END DO
1459 :
1460 264364 : DO bx = 2, lb
1461 27140 : f3 = f2*REAL(bx - 1, dp)
1462 81420 : DO by = 0, lb - bx
1463 27668 : bz = lb - bx - by
1464 : s(1, coset(bx, by, bz), l) = rbp(1)*s(1, coset(bx - 1, by, bz), l) + &
1465 27668 : f3*s(1, coset(bx - 2, by, bz), l)
1466 27668 : IF (lx1 > 0) s(1, coset(bx, by, bz), l) = s(1, coset(bx, by, bz), l) + &
1467 34213 : f2x*s(1, coset(bx - 1, by, bz), lx1)
1468 : END DO
1469 : END DO
1470 :
1471 : END DO
1472 :
1473 : END IF
1474 :
1475 : END IF
1476 :
1477 : END DO
1478 :
1479 12898021 : DO k = 2, ncoset(lc_max)
1480 101039220 : DO j = 1, ncoset(lb_max)
1481 888151017 : DO i = 1, ncoset(la_max)
1482 876762562 : mab(na + i, nb + j, k - 1) = s(i, j, k)
1483 : END DO
1484 : END DO
1485 : END DO
1486 :
1487 2979770 : nb = nb + ncoset(lb_max)
1488 :
1489 : END DO
1490 :
1491 2118970 : na = na + ncoset(la_max)
1492 :
1493 : END DO
1494 :
1495 648766 : END SUBROUTINE moment
1496 :
1497 : ! **************************************************************************************************
1498 : !> \brief This returns the derivative of the moment integrals [a|\mu|b], with respect
1499 : !> to the position of the primitive on the left, i.e.
1500 : !> [da/dR_ai|\mu|b] = 2*zeta*[a+1i|\mu|b] - Ni(a)[a-1i|\mu|b]
1501 : !> order indicates the max order of the moment operator to be calculated
1502 : !> 1: dipole
1503 : !> 2: quadrupole
1504 : !> ...
1505 : !> \param la_max ...
1506 : !> \param npgfa ...
1507 : !> \param zeta ...
1508 : !> \param rpgfa ...
1509 : !> \param la_min ...
1510 : !> \param lb_max ...
1511 : !> \param npgfb ...
1512 : !> \param zetb ...
1513 : !> \param rpgfb ...
1514 : !> \param lb_min ...
1515 : !> \param order ...
1516 : !> \param rac ...
1517 : !> \param rbc ...
1518 : !> \param difmab ...
1519 : !> \param mab_ext ...
1520 : !> \note
1521 : ! **************************************************************************************************
1522 597936 : SUBROUTINE diff_momop(la_max, npgfa, zeta, rpgfa, la_min, &
1523 597936 : lb_max, npgfb, zetb, rpgfb, lb_min, &
1524 597936 : order, rac, rbc, difmab, mab_ext)
1525 :
1526 : INTEGER, INTENT(IN) :: la_max, npgfa
1527 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
1528 : INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
1529 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
1530 : INTEGER, INTENT(IN) :: lb_min, order
1531 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
1532 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(OUT) :: difmab
1533 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
1534 : POINTER :: mab_ext
1535 :
1536 : INTEGER :: imom, lda, lda_min, ldb, ldb_min
1537 : REAL(KIND=dp) :: dab, rab(3), rab2
1538 597936 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: difmab_tmp
1539 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mab
1540 :
1541 2391744 : rab = rbc - rac
1542 2391744 : rab2 = SUM(rab**2)
1543 597936 : dab = SQRT(rab2)
1544 :
1545 597936 : lda_min = MAX(0, la_min - 1)
1546 597936 : ldb_min = MAX(0, lb_min - 1)
1547 597936 : lda = ncoset(la_max)*npgfa
1548 597936 : ldb = ncoset(lb_max)*npgfb
1549 2962848 : ALLOCATE (difmab_tmp(lda, ldb, 3))
1550 :
1551 597936 : IF (PRESENT(mab_ext)) THEN
1552 597936 : mab => mab_ext
1553 : ELSE
1554 : ALLOCATE (mab(npgfa*ncoset(la_max + 1), npgfb*ncoset(lb_max + 1), &
1555 0 : ncoset(order) - 1))
1556 0 : mab = 0.0_dp
1557 : ! *** Calculate the primitive overlap integrals ***
1558 : CALL moment(la_max + 1, npgfa, zeta, rpgfa, lda_min, &
1559 : lb_max + 1, npgfb, zetb, rpgfb, &
1560 0 : order, rac, rbc, mab)
1561 :
1562 : END IF
1563 5973672 : DO imom = 1, ncoset(order) - 1
1564 5375736 : difmab_tmp = 0.0_dp
1565 : CALL adbdr(la_max, npgfa, rpgfa, la_min, &
1566 : lb_max, npgfb, zetb, rpgfb, lb_min, &
1567 : dab, mab(:, :, imom), difmab_tmp(:, :, 1), &
1568 5375736 : difmab_tmp(:, :, 2), difmab_tmp(:, :, 3))
1569 400516668 : difmab(1:lda, 1:ldb, imom, 1) = difmab_tmp(1:lda, 1:ldb, 1)
1570 400516668 : difmab(1:lda, 1:ldb, imom, 2) = difmab_tmp(1:lda, 1:ldb, 2)
1571 401114604 : difmab(1:lda, 1:ldb, imom, 3) = difmab_tmp(1:lda, 1:ldb, 3)
1572 : END DO
1573 :
1574 597936 : IF (PRESENT(mab_ext)) THEN
1575 : NULLIFY (mab)
1576 : ELSE
1577 0 : DEALLOCATE (mab)
1578 : END IF
1579 597936 : DEALLOCATE (difmab_tmp)
1580 :
1581 597936 : END SUBROUTINE diff_momop
1582 :
1583 : ! **************************************************************************************************
1584 : !> \brief This returns the derivative of the dipole integrals [a|x|b], with respect
1585 : !> to the position of the primitive on the left and right, i.e.
1586 : !> [da/dR_ai|\mu|b] = 2*zeta*[a+1i|\mu|b] - Ni(a)[a-1i|\mu|b]
1587 : !> \param la_max ...
1588 : !> \param npgfa ...
1589 : !> \param zeta ...
1590 : !> \param rpgfa ...
1591 : !> \param la_min ...
1592 : !> \param lb_max ...
1593 : !> \param npgfb ...
1594 : !> \param zetb ...
1595 : !> \param rpgfb ...
1596 : !> \param lb_min ...
1597 : !> \param order ...
1598 : !> \param rac ...
1599 : !> \param rbc ...
1600 : !> \param pab ...
1601 : !> \param forcea ...
1602 : !> \param forceb ...
1603 : !> \note
1604 : ! **************************************************************************************************
1605 2124 : SUBROUTINE dipole_force(la_max, npgfa, zeta, rpgfa, la_min, &
1606 2124 : lb_max, npgfb, zetb, rpgfb, lb_min, &
1607 2124 : order, rac, rbc, pab, forcea, forceb)
1608 :
1609 : INTEGER, INTENT(IN) :: la_max, npgfa
1610 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
1611 : INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
1612 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
1613 : INTEGER, INTENT(IN) :: lb_min, order
1614 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
1615 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: pab
1616 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: forcea, forceb
1617 :
1618 : INTEGER :: i, imom, ipgf, j, jpgf, lda, lda_min, &
1619 : ldb, ldb_min, na, nb
1620 : REAL(KIND=dp) :: dab, rab(3), rab2
1621 2124 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: difmab, mab
1622 :
1623 2124 : CPASSERT(order == 1)
1624 : MARK_USED(order)
1625 :
1626 8496 : rab = rbc - rac
1627 8496 : rab2 = SUM(rab**2)
1628 2124 : dab = SQRT(rab2)
1629 :
1630 2124 : lda_min = MAX(0, la_min - 1)
1631 2124 : ldb_min = MAX(0, lb_min - 1)
1632 2124 : lda = ncoset(la_max)*npgfa
1633 2124 : ldb = ncoset(lb_max)*npgfb
1634 10620 : ALLOCATE (difmab(lda, ldb, 3))
1635 10620 : ALLOCATE (mab(npgfa*ncoset(la_max + 1), npgfb*ncoset(lb_max + 1), 3))
1636 2124 : mab = 0.0_dp
1637 : CALL moment(la_max + 1, npgfa, zeta, rpgfa, lda_min, &
1638 2124 : lb_max + 1, npgfb, zetb, rpgfb, 1, rac, rbc, mab)
1639 :
1640 8496 : DO imom = 1, 3
1641 6372 : difmab = 0.0_dp
1642 : CALL adbdr(la_max, npgfa, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, &
1643 6372 : dab, mab(:, :, imom), difmab(:, :, 1), difmab(:, :, 2), difmab(:, :, 3))
1644 6372 : na = 0
1645 24360 : DO ipgf = 1, npgfa
1646 : nb = 0
1647 69429 : DO jpgf = 1, npgfb
1648 171363 : DO j = nb + ncoset(lb_min - 1) + 1, nb + ncoset(lb_max)
1649 518880 : DO i = na + ncoset(la_min - 1) + 1, na + ncoset(la_max)
1650 347517 : forceb(imom, 1) = forceb(imom, 1) + pab(i, j)*difmab(i, j, 1)
1651 347517 : forceb(imom, 2) = forceb(imom, 2) + pab(i, j)*difmab(i, j, 2)
1652 467439 : forceb(imom, 3) = forceb(imom, 3) + pab(i, j)*difmab(i, j, 3)
1653 : END DO
1654 : END DO
1655 69429 : nb = nb + ncoset(lb_max)
1656 : END DO
1657 24360 : na = na + ncoset(la_max)
1658 : END DO
1659 :
1660 6372 : difmab = 0.0_dp
1661 : CALL dabdr(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, rpgfb, lb_min, &
1662 6372 : dab, mab(:, :, imom), difmab(:, :, 1), difmab(:, :, 2), difmab(:, :, 3))
1663 6372 : na = 0
1664 26484 : DO ipgf = 1, npgfa
1665 : nb = 0
1666 69429 : DO jpgf = 1, npgfb
1667 171363 : DO j = nb + ncoset(lb_min - 1) + 1, nb + ncoset(lb_max)
1668 518880 : DO i = na + ncoset(la_min - 1) + 1, na + ncoset(la_max)
1669 347517 : forcea(imom, 1) = forcea(imom, 1) + pab(i, j)*difmab(i, j, 1)
1670 347517 : forcea(imom, 2) = forcea(imom, 2) + pab(i, j)*difmab(i, j, 2)
1671 467439 : forcea(imom, 3) = forcea(imom, 3) + pab(i, j)*difmab(i, j, 3)
1672 : END DO
1673 : END DO
1674 69429 : nb = nb + ncoset(lb_max)
1675 : END DO
1676 24360 : na = na + ncoset(la_max)
1677 : END DO
1678 : END DO
1679 :
1680 2124 : DEALLOCATE (mab, difmab)
1681 :
1682 2124 : END SUBROUTINE dipole_force
1683 :
1684 : END MODULE ai_moments
|
