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 :
49 : CONTAINS
50 :
51 : ! **************************************************************************************************
52 : !> \brief ...
53 : !> \param cos_block ...
54 : !> \param sin_block ...
55 : !> \param iatom ...
56 : !> \param ncoa ...
57 : !> \param nsgfa ...
58 : !> \param sgfa ...
59 : !> \param sphi_a ...
60 : !> \param ldsa ...
61 : !> \param jatom ...
62 : !> \param ncob ...
63 : !> \param nsgfb ...
64 : !> \param sgfb ...
65 : !> \param sphi_b ...
66 : !> \param ldsb ...
67 : !> \param cosab ...
68 : !> \param sinab ...
69 : !> \param ldab ...
70 : !> \param work ...
71 : !> \param ldwork ...
72 : !> \param ordered retain the ordered atom pair instead of mapping to an upper block
73 : ! **************************************************************************************************
74 1514470 : SUBROUTINE contract_cossin(cos_block, sin_block, &
75 3028940 : iatom, ncoa, nsgfa, sgfa, sphi_a, ldsa, &
76 3028940 : jatom, ncob, nsgfb, sgfb, sphi_b, ldsb, &
77 1514470 : cosab, sinab, ldab, work, ldwork, ordered)
78 :
79 : REAL(dp), DIMENSION(:, :), POINTER :: cos_block, sin_block
80 : INTEGER, INTENT(IN) :: iatom, ncoa, nsgfa, sgfa
81 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: sphi_a
82 : INTEGER, INTENT(IN) :: ldsa, jatom, ncob, nsgfb, sgfb
83 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: sphi_b
84 : INTEGER, INTENT(IN) :: ldsb
85 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: cosab, sinab
86 : INTEGER, INTENT(IN) :: ldab
87 : REAL(dp), DIMENSION(:, :) :: work
88 : INTEGER, INTENT(IN) :: ldwork
89 : LOGICAL, INTENT(IN), OPTIONAL :: ordered
90 :
91 : LOGICAL :: direct_block
92 :
93 1514470 : direct_block = iatom <= jatom
94 1514470 : IF (PRESENT(ordered)) direct_block = direct_block .OR. ordered
95 :
96 : ! Calculate cosine
97 :
98 : CALL dgemm("N", "N", ncoa, nsgfb, ncob, &
99 : 1.0_dp, cosab(1, 1), ldab, &
100 : sphi_b(1, sgfb), ldsb, &
101 1514470 : 0.0_dp, work(1, 1), ldwork)
102 :
103 1514470 : IF (direct_block) THEN
104 : CALL dgemm("T", "N", nsgfa, nsgfb, ncoa, &
105 : 1.0_dp, sphi_a(1, sgfa), ldsa, &
106 : work(1, 1), ldwork, &
107 : 1.0_dp, cos_block(sgfa, sgfb), &
108 939576 : SIZE(cos_block, 1))
109 : ELSE
110 : CALL dgemm("T", "N", nsgfb, nsgfa, ncoa, &
111 : 1.0_dp, work(1, 1), ldwork, &
112 : sphi_a(1, sgfa), ldsa, &
113 : 1.0_dp, cos_block(sgfb, sgfa), &
114 574894 : SIZE(cos_block, 1))
115 : END IF
116 :
117 : ! Calculate sine
118 : CALL dgemm("N", "N", ncoa, nsgfb, ncob, &
119 : 1.0_dp, sinab(1, 1), ldab, &
120 : sphi_b(1, sgfb), ldsb, &
121 1514470 : 0.0_dp, work(1, 1), ldwork)
122 :
123 1514470 : IF (direct_block) THEN
124 : CALL dgemm("T", "N", nsgfa, nsgfb, ncoa, &
125 : 1.0_dp, sphi_a(1, sgfa), ldsa, &
126 : work(1, 1), ldwork, &
127 : 1.0_dp, sin_block(sgfa, sgfb), &
128 939576 : SIZE(sin_block, 1))
129 : ELSE
130 : CALL dgemm("T", "N", nsgfb, nsgfa, ncoa, &
131 : 1.0_dp, work(1, 1), ldwork, &
132 : sphi_a(1, sgfa), ldsa, &
133 : 1.0_dp, sin_block(sgfb, sgfa), &
134 574894 : SIZE(sin_block, 1))
135 : END IF
136 :
137 1514470 : END SUBROUTINE contract_cossin
138 :
139 : ! **************************************************************************************************
140 : !> \brief ...
141 : !> \param la_max_set ...
142 : !> \param npgfa ...
143 : !> \param zeta ...
144 : !> \param rpgfa ...
145 : !> \param la_min_set ...
146 : !> \param lb_max ...
147 : !> \param npgfb ...
148 : !> \param zetb ...
149 : !> \param rpgfb ...
150 : !> \param lb_min ...
151 : !> \param rac ...
152 : !> \param rbc ...
153 : !> \param kvec ...
154 : !> \param cosab ...
155 : !> \param sinab ...
156 : !> \param dcosab ...
157 : !> \param dsinab ...
158 : ! **************************************************************************************************
159 1541596 : SUBROUTINE cossin(la_max_set, npgfa, zeta, rpgfa, la_min_set, &
160 1541596 : lb_max, npgfb, zetb, rpgfb, lb_min, &
161 1541596 : rac, rbc, kvec, cosab, sinab, dcosab, dsinab)
162 :
163 : INTEGER, INTENT(IN) :: la_max_set, npgfa
164 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
165 : INTEGER, INTENT(IN) :: la_min_set, lb_max, npgfb
166 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
167 : INTEGER, INTENT(IN) :: lb_min
168 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac, rbc, kvec
169 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: cosab, sinab
170 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
171 : OPTIONAL :: dcosab, dsinab
172 :
173 : INTEGER :: ax, ay, az, bx, by, bz, cda, cdax, cday, cdaz, coa, coamx, coamy, coamz, coapx, &
174 : coapy, coapz, cob, da, da_max, dax, day, daz, i, ipgf, j, jpgf, k, la, la_max, la_min, &
175 : la_start, lb, lb_start, na, nb
176 : REAL(KIND=dp) :: dab, f0, f1, f2, f3, fax, fay, faz, ftz, &
177 : fx, fy, fz, k2, kdp, rab2, s, zetp
178 : REAL(KIND=dp), DIMENSION(3) :: rab, rap, rbp
179 : REAL(KIND=dp), DIMENSION(ncoset(la_max_set), &
180 3083192 : ncoset(lb_max), 3) :: dscos, dssin
181 : REAL(KIND=dp), &
182 1541596 : DIMENSION(ncoset(la_max_set+1), ncoset(lb_max)) :: sc, ss
183 :
184 6166384 : rab = rbc - rac
185 6166384 : rab2 = SUM(rab**2)
186 1541596 : dab = SQRT(rab2)
187 1541596 : k2 = kvec(1)*kvec(1) + kvec(2)*kvec(2) + kvec(3)*kvec(3)
188 :
189 1541596 : IF (PRESENT(dcosab)) THEN
190 24916 : da_max = 1
191 24916 : la_max = la_max_set + 1
192 24916 : la_min = MAX(0, la_min_set - 1)
193 1041304 : dscos = 0.0_dp
194 1041304 : dssin = 0.0_dp
195 : ELSE
196 1516680 : da_max = 0
197 1516680 : la_max = la_max_set
198 1516680 : la_min = la_min_set
199 : END IF
200 :
201 : ! initialize all matrix elements to zero
202 1541596 : IF (PRESENT(dcosab)) THEN
203 24916 : na = ncoset(la_max - 1)*npgfa
204 : ELSE
205 1516680 : na = ncoset(la_max)*npgfa
206 : END IF
207 1541596 : nb = ncoset(lb_max)*npgfb
208 225011239 : cosab(1:na, 1:nb) = 0.0_dp
209 225011239 : sinab(1:na, 1:nb) = 0.0_dp
210 1541596 : IF (PRESENT(dcosab)) THEN
211 5627902 : dcosab(1:na, 1:nb, :) = 0.0_dp
212 5627902 : dsinab(1:na, 1:nb, :) = 0.0_dp
213 : END IF
214 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
215 :
216 1541596 : na = 0
217 5144838 : DO ipgf = 1, npgfa
218 :
219 : nb = 0
220 :
221 15128222 : DO jpgf = 1, npgfb
222 :
223 548012893 : ss = 0.0_dp
224 548012893 : sc = 0.0_dp
225 :
226 : ! *** Screening ***
227 11524980 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
228 7318286 : nb = nb + ncoset(lb_max)
229 7318286 : CYCLE
230 : END IF
231 :
232 : ! *** Calculate some prefactors ***
233 :
234 4206694 : zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
235 :
236 4206694 : f0 = (pi*zetp)**1.5_dp
237 4206694 : f1 = zetb(jpgf)*zetp
238 4206694 : f2 = 0.5_dp*zetp
239 :
240 16826776 : kdp = zetp*DOT_PRODUCT(kvec, zeta(ipgf)*rac + zetb(jpgf)*rbc)
241 :
242 : ! *** Calculate the basic two-center cos/sin integral [s|cos/sin|s] ***
243 :
244 4206694 : s = f0*EXP(-zeta(ipgf)*f1*rab2)*EXP(-0.25_dp*k2*zetp)
245 4206694 : sc(1, 1) = s*COS(kdp)
246 4206694 : ss(1, 1) = s*SIN(kdp)
247 :
248 : ! *** Recurrence steps: [s|O|s] -> [a|O|b] ***
249 :
250 4206694 : IF (la_max > 0) THEN
251 :
252 : ! *** Vertical recurrence steps: [s|O|s] -> [a|O|s] ***
253 :
254 10695380 : rap(:) = f1*rab(:)
255 :
256 : ! *** [p|O|s] = (Pi - Ai)*[s|O|s] +[s|dO|s] (i = x,y,z) ***
257 :
258 2673845 : sc(2, 1) = rap(1)*sc(1, 1) - f2*kvec(1)*ss(1, 1)
259 2673845 : sc(3, 1) = rap(2)*sc(1, 1) - f2*kvec(2)*ss(1, 1)
260 2673845 : sc(4, 1) = rap(3)*sc(1, 1) - f2*kvec(3)*ss(1, 1)
261 2673845 : ss(2, 1) = rap(1)*ss(1, 1) + f2*kvec(1)*sc(1, 1)
262 2673845 : ss(3, 1) = rap(2)*ss(1, 1) + f2*kvec(2)*sc(1, 1)
263 2673845 : ss(4, 1) = rap(3)*ss(1, 1) + f2*kvec(3)*sc(1, 1)
264 :
265 : ! *** [a|O|s] = (Pi - Ai)*[a-1i|O|s] + f2*Ni(a-1i)*[a-2i|s] ***
266 : ! *** + [a-1i|dO|s] ***
267 :
268 3200203 : DO la = 2, la_max
269 :
270 : ! *** Increase the angular momentum component z of function a ***
271 :
272 : sc(coset(0, 0, la), 1) = rap(3)*sc(coset(0, 0, la - 1), 1) + &
273 : f2*REAL(la - 1, dp)*sc(coset(0, 0, la - 2), 1) - &
274 526358 : f2*kvec(3)*ss(coset(0, 0, la - 1), 1)
275 : ss(coset(0, 0, la), 1) = rap(3)*ss(coset(0, 0, la - 1), 1) + &
276 : f2*REAL(la - 1, dp)*ss(coset(0, 0, la - 2), 1) + &
277 526358 : f2*kvec(3)*sc(coset(0, 0, la - 1), 1)
278 :
279 : ! *** Increase the angular momentum component y of function a ***
280 :
281 526358 : az = la - 1
282 : sc(coset(0, 1, az), 1) = rap(2)*sc(coset(0, 0, az), 1) - &
283 526358 : f2*kvec(2)*ss(coset(0, 0, az), 1)
284 : ss(coset(0, 1, az), 1) = rap(2)*ss(coset(0, 0, az), 1) + &
285 526358 : f2*kvec(2)*sc(coset(0, 0, az), 1)
286 :
287 1065436 : DO ay = 2, la
288 539078 : az = la - ay
289 : sc(coset(0, ay, az), 1) = rap(2)*sc(coset(0, ay - 1, az), 1) + &
290 : f2*REAL(ay - 1, dp)*sc(coset(0, ay - 2, az), 1) - &
291 539078 : f2*kvec(2)*ss(coset(0, ay - 1, az), 1)
292 : ss(coset(0, ay, az), 1) = rap(2)*ss(coset(0, ay - 1, az), 1) + &
293 : f2*REAL(ay - 1, dp)*ss(coset(0, ay - 2, az), 1) + &
294 1065436 : f2*kvec(2)*sc(coset(0, ay - 1, az), 1)
295 : END DO
296 :
297 : ! *** Increase the angular momentum component x of function a ***
298 :
299 1591794 : DO ay = 0, la - 1
300 1065436 : az = la - 1 - ay
301 : sc(coset(1, ay, az), 1) = rap(1)*sc(coset(0, ay, az), 1) - &
302 1065436 : f2*kvec(1)*ss(coset(0, ay, az), 1)
303 : ss(coset(1, ay, az), 1) = rap(1)*ss(coset(0, ay, az), 1) + &
304 1591794 : f2*kvec(1)*sc(coset(0, ay, az), 1)
305 : END DO
306 :
307 3739281 : DO ax = 2, la
308 539078 : f3 = f2*REAL(ax - 1, dp)
309 1617495 : DO ay = 0, la - ax
310 552059 : az = la - ax - ay
311 : sc(coset(ax, ay, az), 1) = rap(1)*sc(coset(ax - 1, ay, az), 1) + &
312 : f3*sc(coset(ax - 2, ay, az), 1) - &
313 552059 : f2*kvec(1)*ss(coset(ax - 1, ay, az), 1)
314 : ss(coset(ax, ay, az), 1) = rap(1)*ss(coset(ax - 1, ay, az), 1) + &
315 : f3*ss(coset(ax - 2, ay, az), 1) + &
316 1091137 : f2*kvec(1)*sc(coset(ax - 1, ay, az), 1)
317 : END DO
318 : END DO
319 :
320 : END DO
321 :
322 : ! *** Recurrence steps: [a|O|s] -> [a|O|b] ***
323 :
324 2673845 : IF (lb_max > 0) THEN
325 :
326 10494944 : DO j = 2, ncoset(lb_max)
327 59350517 : DO i = 1, ncoset(la_max)
328 48855573 : sc(i, j) = 0.0_dp
329 57313284 : ss(i, j) = 0.0_dp
330 : END DO
331 : END DO
332 :
333 : ! *** Horizontal recurrence steps ***
334 :
335 8148932 : rbp(:) = rap(:) - rab(:)
336 :
337 : ! *** [a|O|p] = [a+1i|O|s] - (Bi - Ai)*[a|O|s] ***
338 :
339 2037233 : IF (lb_max == 1) THEN
340 : la_start = la_min
341 : ELSE
342 382422 : la_start = MAX(0, la_min - 1)
343 : END IF
344 :
345 4047295 : DO la = la_start, la_max - 1
346 6330507 : DO ax = 0, la
347 6852804 : DO ay = 0, la - ax
348 2559530 : az = la - ax - ay
349 : sc(coset(ax, ay, az), 2) = sc(coset(ax + 1, ay, az), 1) - &
350 2559530 : rab(1)*sc(coset(ax, ay, az), 1)
351 : sc(coset(ax, ay, az), 3) = sc(coset(ax, ay + 1, az), 1) - &
352 2559530 : rab(2)*sc(coset(ax, ay, az), 1)
353 : sc(coset(ax, ay, az), 4) = sc(coset(ax, ay, az + 1), 1) - &
354 2559530 : rab(3)*sc(coset(ax, ay, az), 1)
355 : ss(coset(ax, ay, az), 2) = ss(coset(ax + 1, ay, az), 1) - &
356 2559530 : rab(1)*ss(coset(ax, ay, az), 1)
357 : ss(coset(ax, ay, az), 3) = ss(coset(ax, ay + 1, az), 1) - &
358 2559530 : rab(2)*ss(coset(ax, ay, az), 1)
359 : ss(coset(ax, ay, az), 4) = ss(coset(ax, ay, az + 1), 1) - &
360 4842742 : rab(3)*ss(coset(ax, ay, az), 1)
361 : END DO
362 : END DO
363 : END DO
364 :
365 : ! *** Vertical recurrence step ***
366 :
367 : ! *** [a|O|p] = (Pi - Bi)*[a|O|s] + f2*Ni(a)*[a-1i|O|s] ***
368 : ! *** + [a|dO|s] ***
369 :
370 6512888 : DO ax = 0, la_max
371 4475655 : fx = f2*REAL(ax, dp)
372 13835192 : DO ay = 0, la_max - ax
373 7322304 : fy = f2*REAL(ay, dp)
374 7322304 : az = la_max - ax - ay
375 7322304 : fz = f2*REAL(az, dp)
376 7322304 : IF (ax == 0) THEN
377 : sc(coset(ax, ay, az), 2) = rbp(1)*sc(coset(ax, ay, az), 1) - &
378 4475655 : f2*kvec(1)*ss(coset(ax, ay, az), 1)
379 : ss(coset(ax, ay, az), 2) = rbp(1)*ss(coset(ax, ay, az), 1) + &
380 4475655 : f2*kvec(1)*sc(coset(ax, ay, az), 1)
381 : ELSE
382 : sc(coset(ax, ay, az), 2) = rbp(1)*sc(coset(ax, ay, az), 1) + &
383 : fx*sc(coset(ax - 1, ay, az), 1) - &
384 2846649 : f2*kvec(1)*ss(coset(ax, ay, az), 1)
385 : ss(coset(ax, ay, az), 2) = rbp(1)*ss(coset(ax, ay, az), 1) + &
386 : fx*ss(coset(ax - 1, ay, az), 1) + &
387 2846649 : f2*kvec(1)*sc(coset(ax, ay, az), 1)
388 : END IF
389 7322304 : IF (ay == 0) THEN
390 : sc(coset(ax, ay, az), 3) = rbp(2)*sc(coset(ax, ay, az), 1) - &
391 4475655 : f2*kvec(2)*ss(coset(ax, ay, az), 1)
392 : ss(coset(ax, ay, az), 3) = rbp(2)*ss(coset(ax, ay, az), 1) + &
393 4475655 : f2*kvec(2)*sc(coset(ax, ay, az), 1)
394 : ELSE
395 : sc(coset(ax, ay, az), 3) = rbp(2)*sc(coset(ax, ay, az), 1) + &
396 : fy*sc(coset(ax, ay - 1, az), 1) - &
397 2846649 : f2*kvec(2)*ss(coset(ax, ay, az), 1)
398 : ss(coset(ax, ay, az), 3) = rbp(2)*ss(coset(ax, ay, az), 1) + &
399 : fy*ss(coset(ax, ay - 1, az), 1) + &
400 2846649 : f2*kvec(2)*sc(coset(ax, ay, az), 1)
401 : END IF
402 11797959 : IF (az == 0) THEN
403 : sc(coset(ax, ay, az), 4) = rbp(3)*sc(coset(ax, ay, az), 1) - &
404 4475655 : f2*kvec(3)*ss(coset(ax, ay, az), 1)
405 : ss(coset(ax, ay, az), 4) = rbp(3)*ss(coset(ax, ay, az), 1) + &
406 4475655 : f2*kvec(3)*sc(coset(ax, ay, az), 1)
407 : ELSE
408 : sc(coset(ax, ay, az), 4) = rbp(3)*sc(coset(ax, ay, az), 1) + &
409 : fz*sc(coset(ax, ay, az - 1), 1) - &
410 2846649 : f2*kvec(3)*ss(coset(ax, ay, az), 1)
411 : ss(coset(ax, ay, az), 4) = rbp(3)*ss(coset(ax, ay, az), 1) + &
412 : fz*ss(coset(ax, ay, az - 1), 1) + &
413 2846649 : f2*kvec(3)*sc(coset(ax, ay, az), 1)
414 : END IF
415 : END DO
416 : END DO
417 :
418 : ! *** Recurrence steps: [a|O|p] -> [a|O|b] ***
419 :
420 2424740 : DO lb = 2, lb_max
421 :
422 : ! *** Horizontal recurrence steps ***
423 :
424 : ! *** [a|O|b] = [a+1i|O|b-1i] - (Bi - Ai)*[a|O|b-1i] ***
425 :
426 387507 : IF (lb == lb_max) THEN
427 : la_start = la_min
428 : ELSE
429 5085 : la_start = MAX(0, la_min - 1)
430 : END IF
431 :
432 840157 : DO la = la_start, la_max - 1
433 1452879 : DO ax = 0, la
434 1838766 : DO ay = 0, la - ax
435 773394 : az = la - ax - ay
436 :
437 : ! *** Shift of angular momentum component z from a to b ***
438 :
439 : sc(coset(ax, ay, az), coset(0, 0, lb)) = &
440 : sc(coset(ax, ay, az + 1), coset(0, 0, lb - 1)) - &
441 773394 : rab(3)*sc(coset(ax, ay, az), coset(0, 0, lb - 1))
442 : ss(coset(ax, ay, az), coset(0, 0, lb)) = &
443 : ss(coset(ax, ay, az + 1), coset(0, 0, lb - 1)) - &
444 773394 : rab(3)*ss(coset(ax, ay, az), coset(0, 0, lb - 1))
445 :
446 : ! *** Shift of angular momentum component y from a to b ***
447 :
448 2322861 : DO by = 1, lb
449 1549467 : bz = lb - by
450 : sc(coset(ax, ay, az), coset(0, by, bz)) = &
451 : sc(coset(ax, ay + 1, az), coset(0, by - 1, bz)) - &
452 1549467 : rab(2)*sc(coset(ax, ay, az), coset(0, by - 1, bz))
453 : ss(coset(ax, ay, az), coset(0, by, bz)) = &
454 : ss(coset(ax, ay + 1, az), coset(0, by - 1, bz)) - &
455 2322861 : rab(2)*ss(coset(ax, ay, az), coset(0, by - 1, bz))
456 : END DO
457 :
458 : ! *** Shift of angular momentum component x from a to b ***
459 :
460 2935583 : DO bx = 1, lb
461 4651080 : DO by = 0, lb - bx
462 2328219 : bz = lb - bx - by
463 : sc(coset(ax, ay, az), coset(bx, by, bz)) = &
464 : sc(coset(ax + 1, ay, az), coset(bx - 1, by, bz)) - &
465 2328219 : rab(1)*sc(coset(ax, ay, az), coset(bx - 1, by, bz))
466 : ss(coset(ax, ay, az), coset(bx, by, bz)) = &
467 : ss(coset(ax + 1, ay, az), coset(bx - 1, by, bz)) - &
468 3877686 : rab(1)*ss(coset(ax, ay, az), coset(bx - 1, by, bz))
469 : END DO
470 : END DO
471 :
472 : END DO
473 : END DO
474 : END DO
475 :
476 : ! *** Vertical recurrence step ***
477 :
478 : ! *** [a|O|b] = (Pi - Bi)*[a|O|b-1i] + f2*Ni(a)*[a-1i|O|b-1i] + ***
479 : ! *** f2*Ni(b-1i)*[a|O|b-2i] + [a|dO|b-1i] ***
480 :
481 3411383 : DO ax = 0, la_max
482 986643 : fx = f2*REAL(ax, dp)
483 3173520 : DO ay = 0, la_max - ax
484 1799370 : fy = f2*REAL(ay, dp)
485 1799370 : az = la_max - ax - ay
486 1799370 : fz = f2*REAL(az, dp)
487 :
488 : ! *** Increase the angular momentum component z of function b ***
489 :
490 1799370 : f3 = f2*REAL(lb - 1, dp)
491 :
492 1799370 : IF (az == 0) THEN
493 : sc(coset(ax, ay, az), coset(0, 0, lb)) = &
494 : rbp(3)*sc(coset(ax, ay, az), coset(0, 0, lb - 1)) + &
495 : f3*sc(coset(ax, ay, az), coset(0, 0, lb - 2)) - &
496 986643 : f2*kvec(3)*ss(coset(ax, ay, az), coset(0, 0, lb - 1))
497 : ss(coset(ax, ay, az), coset(0, 0, lb)) = &
498 : rbp(3)*ss(coset(ax, ay, az), coset(0, 0, lb - 1)) + &
499 : f3*ss(coset(ax, ay, az), coset(0, 0, lb - 2)) + &
500 986643 : f2*kvec(3)*sc(coset(ax, ay, az), coset(0, 0, lb - 1))
501 : ELSE
502 : sc(coset(ax, ay, az), coset(0, 0, lb)) = &
503 : rbp(3)*sc(coset(ax, ay, az), coset(0, 0, lb - 1)) + &
504 : fz*sc(coset(ax, ay, az - 1), coset(0, 0, lb - 1)) + &
505 : f3*sc(coset(ax, ay, az), coset(0, 0, lb - 2)) - &
506 812727 : f2*kvec(3)*ss(coset(ax, ay, az), coset(0, 0, lb - 1))
507 : ss(coset(ax, ay, az), coset(0, 0, lb)) = &
508 : rbp(3)*ss(coset(ax, ay, az), coset(0, 0, lb - 1)) + &
509 : fz*ss(coset(ax, ay, az - 1), coset(0, 0, lb - 1)) + &
510 : f3*ss(coset(ax, ay, az), coset(0, 0, lb - 2)) + &
511 812727 : f2*kvec(3)*sc(coset(ax, ay, az), coset(0, 0, lb - 1))
512 : END IF
513 :
514 : ! *** Increase the angular momentum component y of function b ***
515 :
516 1799370 : IF (ay == 0) THEN
517 986643 : bz = lb - 1
518 : sc(coset(ax, ay, az), coset(0, 1, bz)) = &
519 : rbp(2)*sc(coset(ax, ay, az), coset(0, 0, bz)) - &
520 986643 : f2*kvec(2)*ss(coset(ax, ay, az), coset(0, 0, bz))
521 : ss(coset(ax, ay, az), coset(0, 1, bz)) = &
522 : rbp(2)*ss(coset(ax, ay, az), coset(0, 0, bz)) + &
523 986643 : f2*kvec(2)*sc(coset(ax, ay, az), coset(0, 0, bz))
524 1985520 : DO by = 2, lb
525 998877 : bz = lb - by
526 998877 : f3 = f2*REAL(by - 1, dp)
527 : sc(coset(ax, ay, az), coset(0, by, bz)) = &
528 : rbp(2)*sc(coset(ax, ay, az), coset(0, by - 1, bz)) + &
529 : f3*sc(coset(ax, ay, az), coset(0, by - 2, bz)) - &
530 998877 : f2*kvec(2)*ss(coset(ax, ay, az), coset(0, by - 1, bz))
531 : ss(coset(ax, ay, az), coset(0, by, bz)) = &
532 : rbp(2)*ss(coset(ax, ay, az), coset(0, by - 1, bz)) + &
533 : f3*ss(coset(ax, ay, az), coset(0, by - 2, bz)) + &
534 1985520 : f2*kvec(2)*sc(coset(ax, ay, az), coset(0, by - 1, bz))
535 : END DO
536 : ELSE
537 812727 : bz = lb - 1
538 : sc(coset(ax, ay, az), coset(0, 1, bz)) = &
539 : rbp(2)*sc(coset(ax, ay, az), coset(0, 0, bz)) + &
540 : fy*sc(coset(ax, ay - 1, az), coset(0, 0, bz)) - &
541 812727 : f2*kvec(2)*ss(coset(ax, ay, az), coset(0, 0, bz))
542 : ss(coset(ax, ay, az), coset(0, 1, bz)) = &
543 : rbp(2)*ss(coset(ax, ay, az), coset(0, 0, bz)) + &
544 : fy*ss(coset(ax, ay - 1, az), coset(0, 0, bz)) + &
545 812727 : f2*kvec(2)*sc(coset(ax, ay, az), coset(0, 0, bz))
546 1634778 : DO by = 2, lb
547 822051 : bz = lb - by
548 822051 : f3 = f2*REAL(by - 1, dp)
549 : sc(coset(ax, ay, az), coset(0, by, bz)) = &
550 : rbp(2)*sc(coset(ax, ay, az), coset(0, by - 1, bz)) + &
551 : fy*sc(coset(ax, ay - 1, az), coset(0, by - 1, bz)) + &
552 : f3*sc(coset(ax, ay, az), coset(0, by - 2, bz)) - &
553 822051 : f2*kvec(2)*ss(coset(ax, ay, az), coset(0, by - 1, bz))
554 : ss(coset(ax, ay, az), coset(0, by, bz)) = &
555 : rbp(2)*ss(coset(ax, ay, az), coset(0, by - 1, bz)) + &
556 : fy*ss(coset(ax, ay - 1, az), coset(0, by - 1, bz)) + &
557 : f3*ss(coset(ax, ay, az), coset(0, by - 2, bz)) + &
558 1634778 : f2*kvec(2)*sc(coset(ax, ay, az), coset(0, by - 1, bz))
559 : END DO
560 : END IF
561 :
562 : ! *** Increase the angular momentum component x of function b ***
563 :
564 2786013 : IF (ax == 0) THEN
565 2972163 : DO by = 0, lb - 1
566 1985520 : bz = lb - 1 - by
567 : sc(coset(ax, ay, az), coset(1, by, bz)) = &
568 : rbp(1)*sc(coset(ax, ay, az), coset(0, by, bz)) - &
569 1985520 : f2*kvec(1)*ss(coset(ax, ay, az), coset(0, by, bz))
570 : ss(coset(ax, ay, az), coset(1, by, bz)) = &
571 : rbp(1)*ss(coset(ax, ay, az), coset(0, by, bz)) + &
572 2972163 : f2*kvec(1)*sc(coset(ax, ay, az), coset(0, by, bz))
573 : END DO
574 1985520 : DO bx = 2, lb
575 998877 : f3 = f2*REAL(bx - 1, dp)
576 2996973 : DO by = 0, lb - bx
577 1011453 : bz = lb - bx - by
578 : sc(coset(ax, ay, az), coset(bx, by, bz)) = &
579 : rbp(1)*sc(coset(ax, ay, az), &
580 : coset(bx - 1, by, bz)) + &
581 : f3*sc(coset(ax, ay, az), coset(bx - 2, by, bz)) - &
582 1011453 : f2*kvec(1)*ss(coset(ax, ay, az), coset(bx - 1, by, bz))
583 : ss(coset(ax, ay, az), coset(bx, by, bz)) = &
584 : rbp(1)*ss(coset(ax, ay, az), &
585 : coset(bx - 1, by, bz)) + &
586 : f3*ss(coset(ax, ay, az), coset(bx - 2, by, bz)) + &
587 2010330 : f2*kvec(1)*sc(coset(ax, ay, az), coset(bx - 1, by, bz))
588 : END DO
589 : END DO
590 : ELSE
591 2447505 : DO by = 0, lb - 1
592 1634778 : bz = lb - 1 - by
593 : sc(coset(ax, ay, az), coset(1, by, bz)) = &
594 : rbp(1)*sc(coset(ax, ay, az), coset(0, by, bz)) + &
595 : fx*sc(coset(ax - 1, ay, az), coset(0, by, bz)) - &
596 1634778 : f2*kvec(1)*ss(coset(ax, ay, az), coset(0, by, bz))
597 : ss(coset(ax, ay, az), coset(1, by, bz)) = &
598 : rbp(1)*ss(coset(ax, ay, az), coset(0, by, bz)) + &
599 : fx*ss(coset(ax - 1, ay, az), coset(0, by, bz)) + &
600 2447505 : f2*kvec(1)*sc(coset(ax, ay, az), coset(0, by, bz))
601 : END DO
602 1634778 : DO bx = 2, lb
603 822051 : f3 = f2*REAL(bx - 1, dp)
604 2466504 : DO by = 0, lb - bx
605 831726 : bz = lb - bx - by
606 : sc(coset(ax, ay, az), coset(bx, by, bz)) = &
607 : rbp(1)*sc(coset(ax, ay, az), &
608 : coset(bx - 1, by, bz)) + &
609 : fx*sc(coset(ax - 1, ay, az), coset(bx - 1, by, bz)) + &
610 : f3*sc(coset(ax, ay, az), coset(bx - 2, by, bz)) - &
611 831726 : f2*kvec(1)*ss(coset(ax, ay, az), coset(bx - 1, by, bz))
612 : ss(coset(ax, ay, az), coset(bx, by, bz)) = &
613 : rbp(1)*ss(coset(ax, ay, az), &
614 : coset(bx - 1, by, bz)) + &
615 : fx*ss(coset(ax - 1, ay, az), coset(bx - 1, by, bz)) + &
616 : f3*ss(coset(ax, ay, az), coset(bx - 2, by, bz)) + &
617 1653777 : f2*kvec(1)*sc(coset(ax, ay, az), coset(bx - 1, by, bz))
618 : END DO
619 : END DO
620 : END IF
621 :
622 : END DO
623 : END DO
624 :
625 : END DO
626 :
627 : END IF
628 :
629 : ELSE
630 :
631 1532849 : IF (lb_max > 0) THEN
632 :
633 : ! *** Vertical recurrence steps: [s|O|s] -> [s|O|b] ***
634 :
635 2326360 : rbp(:) = (f1 - 1.0_dp)*rab(:)
636 :
637 : ! *** [s|O|p] = (Pi - Bi)*[s|O|s] + [s|dO|s] ***
638 :
639 581590 : sc(1, 2) = rbp(1)*sc(1, 1) - f2*kvec(1)*ss(1, 1)
640 581590 : sc(1, 3) = rbp(2)*sc(1, 1) - f2*kvec(2)*ss(1, 1)
641 581590 : sc(1, 4) = rbp(3)*sc(1, 1) - f2*kvec(3)*ss(1, 1)
642 581590 : ss(1, 2) = rbp(1)*ss(1, 1) + f2*kvec(1)*sc(1, 1)
643 581590 : ss(1, 3) = rbp(2)*ss(1, 1) + f2*kvec(2)*sc(1, 1)
644 581590 : ss(1, 4) = rbp(3)*ss(1, 1) + f2*kvec(3)*sc(1, 1)
645 :
646 : ! *** [s|O|b] = (Pi - Bi)*[s|O|b-1i] + f2*Ni(b-1i)*[s|O|b-2i] ***
647 : ! *** + [s|dO|b-1i] ***
648 :
649 686669 : DO lb = 2, lb_max
650 :
651 : ! *** Increase the angular momentum component z of function b ***
652 :
653 : sc(1, coset(0, 0, lb)) = rbp(3)*sc(1, coset(0, 0, lb - 1)) + &
654 : f2*REAL(lb - 1, dp)*sc(1, coset(0, 0, lb - 2)) - &
655 105079 : f2*kvec(3)*ss(1, coset(0, 0, lb - 1))
656 : ss(1, coset(0, 0, lb)) = rbp(3)*ss(1, coset(0, 0, lb - 1)) + &
657 : f2*REAL(lb - 1, dp)*ss(1, coset(0, 0, lb - 2)) + &
658 105079 : f2*kvec(3)*sc(1, coset(0, 0, lb - 1))
659 :
660 : ! *** Increase the angular momentum component y of function b ***
661 :
662 105079 : bz = lb - 1
663 : sc(1, coset(0, 1, bz)) = rbp(2)*sc(1, coset(0, 0, bz)) - &
664 105079 : f2*kvec(2)*ss(1, coset(0, 0, bz))
665 : ss(1, coset(0, 1, bz)) = rbp(2)*ss(1, coset(0, 0, bz)) + &
666 105079 : f2*kvec(2)*sc(1, coset(0, 0, bz))
667 :
668 214340 : DO by = 2, lb
669 109261 : bz = lb - by
670 : sc(1, coset(0, by, bz)) = rbp(2)*sc(1, coset(0, by - 1, bz)) + &
671 : f2*REAL(by - 1, dp)*sc(1, coset(0, by - 2, bz)) - &
672 109261 : f2*kvec(2)*ss(1, coset(0, by - 1, bz))
673 : ss(1, coset(0, by, bz)) = rbp(2)*ss(1, coset(0, by - 1, bz)) + &
674 : f2*REAL(by - 1, dp)*ss(1, coset(0, by - 2, bz)) + &
675 214340 : f2*kvec(2)*sc(1, coset(0, by - 1, bz))
676 : END DO
677 :
678 : ! *** Increase the angular momentum component x of function b ***
679 :
680 319419 : DO by = 0, lb - 1
681 214340 : bz = lb - 1 - by
682 : sc(1, coset(1, by, bz)) = rbp(1)*sc(1, coset(0, by, bz)) - &
683 214340 : f2*kvec(1)*ss(1, coset(0, by, bz))
684 : ss(1, coset(1, by, bz)) = rbp(1)*ss(1, coset(0, by, bz)) + &
685 319419 : f2*kvec(1)*sc(1, coset(0, by, bz))
686 : END DO
687 :
688 795930 : DO bx = 2, lb
689 109261 : f3 = f2*REAL(bx - 1, dp)
690 327918 : DO by = 0, lb - bx
691 113578 : bz = lb - bx - by
692 : sc(1, coset(bx, by, bz)) = rbp(1)*sc(1, coset(bx - 1, by, bz)) + &
693 : f3*sc(1, coset(bx - 2, by, bz)) - &
694 113578 : f2*kvec(1)*ss(1, coset(bx - 1, by, bz))
695 : ss(1, coset(bx, by, bz)) = rbp(1)*ss(1, coset(bx - 1, by, bz)) + &
696 : f3*ss(1, coset(bx - 2, by, bz)) + &
697 222839 : f2*kvec(1)*sc(1, coset(bx - 1, by, bz))
698 : END DO
699 : END DO
700 :
701 : END DO
702 :
703 : END IF
704 :
705 : END IF
706 :
707 17699751 : DO j = ncoset(lb_min - 1) + 1, ncoset(lb_max)
708 73156133 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
709 55456382 : cosab(na + i, nb + j) = sc(i, j)
710 68949439 : sinab(na + i, nb + j) = ss(i, j)
711 : END DO
712 : END DO
713 :
714 4206694 : IF (PRESENT(dcosab)) THEN
715 : la_start = 0
716 : lb_start = 0
717 : ELSE
718 4145450 : la_start = la_min
719 4145450 : lb_start = lb_min
720 : END IF
721 :
722 4267938 : DO da = 0, da_max - 1
723 61244 : ftz = 2.0_dp*zeta(ipgf)
724 4329182 : DO dax = 0, da
725 183732 : DO day = 0, da - dax
726 61244 : daz = da - dax - day
727 61244 : cda = coset(dax, day, daz) - 1
728 61244 : cdax = coset(dax + 1, day, daz) - 1
729 61244 : cday = coset(dax, day + 1, daz) - 1
730 61244 : cdaz = coset(dax, day, daz + 1) - 1
731 : !*** [da/dAi|O|b] = 2*zeta*[a+1i|O|b] - Ni(a)[a-1i|O|b] ***
732 :
733 213755 : DO la = la_start, la_max - da - 1
734 276858 : DO ax = 0, la
735 124347 : fax = REAL(ax, dp)
736 376098 : DO ay = 0, la - ax
737 160484 : fay = REAL(ay, dp)
738 160484 : az = la - ax - ay
739 160484 : faz = REAL(az, dp)
740 160484 : coa = coset(ax, ay, az)
741 160484 : coamx = coset(ax - 1, ay, az)
742 160484 : coamy = coset(ax, ay - 1, az)
743 160484 : coamz = coset(ax, ay, az - 1)
744 160484 : coapx = coset(ax + 1, ay, az)
745 160484 : coapy = coset(ax, ay + 1, az)
746 160484 : coapz = coset(ax, ay, az + 1)
747 533290 : DO lb = lb_start, lb_max
748 754566 : DO bx = 0, lb
749 1046058 : DO by = 0, lb - bx
750 451976 : bz = lb - bx - by
751 451976 : cob = coset(bx, by, bz)
752 451976 : dscos(coa, cob, cdax) = ftz*sc(coapx, cob) - fax*sc(coamx, cob)
753 451976 : dscos(coa, cob, cday) = ftz*sc(coapy, cob) - fay*sc(coamy, cob)
754 451976 : dscos(coa, cob, cdaz) = ftz*sc(coapz, cob) - faz*sc(coamz, cob)
755 451976 : dssin(coa, cob, cdax) = ftz*ss(coapx, cob) - fax*ss(coamx, cob)
756 451976 : dssin(coa, cob, cday) = ftz*ss(coapy, cob) - fay*ss(coamy, cob)
757 797599 : dssin(coa, cob, cdaz) = ftz*ss(coapz, cob) - faz*ss(coamz, cob)
758 : END DO
759 : END DO
760 : END DO
761 : END DO
762 : END DO
763 : END DO
764 :
765 : END DO
766 : END DO
767 : END DO
768 :
769 4206694 : IF (PRESENT(dcosab)) THEN
770 244976 : DO k = 1, 3
771 715556 : DO j = 1, ncoset(lb_max)
772 2010240 : DO i = 1, ncoset(la_max_set)
773 1355928 : dcosab(na + i, nb + j, k) = dscos(i, j, k)
774 1826508 : dsinab(na + i, nb + j, k) = dssin(i, j, k)
775 : END DO
776 : END DO
777 : END DO
778 : END IF
779 :
780 7809936 : nb = nb + ncoset(lb_max)
781 :
782 : END DO
783 :
784 5144838 : na = na + ncoset(la_max_set)
785 :
786 : END DO
787 :
788 1541596 : END SUBROUTINE cossin
789 :
790 : ! **************************************************************************************************
791 : !> \brief ...
792 : !> \param la_max ...
793 : !> \param npgfa ...
794 : !> \param zeta ...
795 : !> \param rpgfa ...
796 : !> \param la_min ...
797 : !> \param lb_max ...
798 : !> \param npgfb ...
799 : !> \param zetb ...
800 : !> \param rpgfb ...
801 : !> \param lc_max ...
802 : !> \param rac ...
803 : !> \param rbc ...
804 : !> \param mab ...
805 : ! **************************************************************************************************
806 648802 : SUBROUTINE moment(la_max, npgfa, zeta, rpgfa, la_min, &
807 1297604 : lb_max, npgfb, zetb, rpgfb, &
808 648802 : lc_max, rac, rbc, mab)
809 :
810 : INTEGER, INTENT(IN) :: la_max, npgfa
811 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
812 : INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
813 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
814 : INTEGER, INTENT(IN) :: lc_max
815 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
816 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: mab
817 :
818 : INTEGER :: ax, ay, az, bx, by, bz, i, ipgf, j, &
819 : jpgf, k, l, l1, l2, la, la_start, lb, &
820 : lx, lx1, ly, ly1, lz, lz1, na, nb, ni
821 : REAL(KIND=dp) :: dab, f0, f1, f2, f2x, f2y, f2z, f3, fx, &
822 : fy, fz, rab2, zetp
823 : REAL(KIND=dp), DIMENSION(3) :: rab, rap, rbp, rpc
824 : REAL(KIND=dp), DIMENSION(ncoset(la_max), ncoset(&
825 648802 : lb_max), ncoset(lc_max)) :: s
826 :
827 2595208 : rab = rbc - rac
828 2595208 : rab2 = SUM(rab**2)
829 648802 : dab = SQRT(rab2)
830 :
831 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
832 :
833 648802 : na = 0
834 :
835 2118772 : DO ipgf = 1, npgfa
836 :
837 1469970 : nb = 0
838 :
839 5448011 : DO jpgf = 1, npgfb
840 :
841 2903726493 : s = 0.0_dp
842 : ! *** Screening ***
843 :
844 3978041 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
845 24081829 : DO k = 1, ncoset(lc_max) - 1
846 195504335 : DO j = nb + 1, nb + ncoset(lb_max)
847 1713647178 : DO i = na + 1, na + ncoset(la_max)
848 1692033881 : mab(i, j, k) = 0.0_dp
849 : END DO
850 : END DO
851 : END DO
852 2468532 : nb = nb + ncoset(lb_max)
853 2468532 : CYCLE
854 : END IF
855 :
856 : ! *** Calculate some prefactors ***
857 :
858 1509509 : zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
859 :
860 1509509 : f0 = (pi*zetp)**1.5_dp
861 1509509 : f1 = zetb(jpgf)*zetp
862 1509509 : f2 = 0.5_dp*zetp
863 :
864 : ! *** Calculate the basic two-center moment integral [s|M|s] ***
865 :
866 6038036 : rpc = zetp*(zeta(ipgf)*rac + zetb(jpgf)*rbc)
867 1509509 : s(1, 1, 1) = f0*EXP(-zeta(ipgf)*f1*rab2)
868 12897505 : DO l = 2, ncoset(lc_max)
869 11387996 : lx = indco(1, l)
870 11387996 : ly = indco(2, l)
871 11387996 : lz = indco(3, l)
872 11387996 : l2 = 0
873 11387996 : IF (lz > 0) THEN
874 4939765 : l1 = coset(lx, ly, lz - 1)
875 4939765 : IF (lz > 1) l2 = coset(lx, ly, lz - 2)
876 : ni = lz - 1
877 : i = 3
878 6448231 : ELSE IF (ly > 0) THEN
879 3795800 : l1 = coset(lx, ly - 1, lz)
880 3795800 : IF (ly > 1) l2 = coset(lx, ly - 2, lz)
881 : ni = ly - 1
882 : i = 2
883 2652431 : ELSE IF (lx > 0) THEN
884 2652431 : l1 = coset(lx - 1, ly, lz)
885 2652431 : IF (lx > 1) l2 = coset(lx - 2, ly, lz)
886 : ni = lx - 1
887 : i = 1
888 : END IF
889 11387996 : s(1, 1, l) = rpc(i)*s(1, 1, l1)
890 12897505 : IF (l2 > 0) s(1, 1, l) = s(1, 1, l) + f2*REAL(ni, dp)*s(1, 1, l2)
891 : END DO
892 :
893 : ! *** Recurrence steps: [s|M|s] -> [a|M|b] ***
894 :
895 14407014 : DO l = 1, ncoset(lc_max)
896 :
897 12897505 : lx = indco(1, l)
898 12897505 : ly = indco(2, l)
899 12897505 : lz = indco(3, l)
900 12897505 : IF (lx > 0) THEN
901 4939765 : lx1 = coset(lx - 1, ly, lz)
902 : ELSE
903 : lx1 = -1
904 : END IF
905 12897505 : IF (ly > 0) THEN
906 4939765 : ly1 = coset(lx, ly - 1, lz)
907 : ELSE
908 : ly1 = -1
909 : END IF
910 12897505 : IF (lz > 0) THEN
911 4939765 : lz1 = coset(lx, ly, lz - 1)
912 : ELSE
913 : lz1 = -1
914 : END IF
915 12897505 : f2x = f2*REAL(lx, dp)
916 12897505 : f2y = f2*REAL(ly, dp)
917 12897505 : f2z = f2*REAL(lz, dp)
918 :
919 14407014 : IF (la_max > 0) THEN
920 :
921 : ! *** Vertical recurrence steps: [s|M|s] -> [a|M|s] ***
922 :
923 49122160 : rap(:) = f1*rab(:)
924 :
925 : ! *** [p|M|s] = (Pi - Ai)*[s|M|s] + f2*Ni(m-1i)[s|M-1i|s] ***
926 :
927 12280540 : s(2, 1, l) = rap(1)*s(1, 1, l)
928 12280540 : s(3, 1, l) = rap(2)*s(1, 1, l)
929 12280540 : s(4, 1, l) = rap(3)*s(1, 1, l)
930 12280540 : IF (lx1 > 0) s(2, 1, l) = s(2, 1, l) + f2x*s(1, 1, lx1)
931 12280540 : IF (ly1 > 0) s(3, 1, l) = s(3, 1, l) + f2y*s(1, 1, ly1)
932 12280540 : IF (lz1 > 0) s(4, 1, l) = s(4, 1, l) + f2z*s(1, 1, lz1)
933 :
934 : ! *** [a|M|s] = (Pi - Ai)*[a-1i|M|s] + f2*Ni(a-1i)*[a-2i|M|s] ***
935 : ! *** + f2*Ni(m-1i)*[a-1i|M-1i|s] ***
936 :
937 19967452 : DO la = 2, la_max
938 :
939 : ! *** Increase the angular momentum component z of function a ***
940 :
941 : s(coset(0, 0, la), 1, l) = rap(3)*s(coset(0, 0, la - 1), 1, l) + &
942 7686912 : f2*REAL(la - 1, dp)*s(coset(0, 0, la - 2), 1, l)
943 7686912 : IF (lz1 > 0) s(coset(0, 0, la), 1, l) = s(coset(0, 0, la), 1, l) + &
944 3044463 : f2z*s(coset(0, 0, la - 1), 1, lz1)
945 :
946 : ! *** Increase the angular momentum component y of function a ***
947 :
948 7686912 : az = la - 1
949 7686912 : s(coset(0, 1, az), 1, l) = rap(2)*s(coset(0, 0, az), 1, l)
950 7686912 : IF (ly1 > 0) s(coset(0, 1, az), 1, l) = s(coset(0, 1, az), 1, l) + &
951 3044463 : f2y*s(coset(0, 0, az), 1, ly1)
952 :
953 16205838 : DO ay = 2, la
954 8518926 : az = la - ay
955 : s(coset(0, ay, az), 1, l) = rap(2)*s(coset(0, ay - 1, az), 1, l) + &
956 8518926 : f2*REAL(ay - 1, dp)*s(coset(0, ay - 2, az), 1, l)
957 8518926 : IF (ly1 > 0) s(coset(0, ay, az), 1, l) = s(coset(0, ay, az), 1, l) + &
958 11062062 : f2y*s(coset(0, ay - 1, az), 1, ly1)
959 : END DO
960 :
961 : ! *** Increase the angular momentum component x of function a ***
962 :
963 23892750 : DO ay = 0, la - 1
964 16205838 : az = la - 1 - ay
965 16205838 : s(coset(1, ay, az), 1, l) = rap(1)*s(coset(0, ay, az), 1, l)
966 16205838 : IF (lx1 > 0) s(coset(1, ay, az), 1, l) = s(coset(1, ay, az), 1, l) + &
967 14106525 : f2x*s(coset(0, ay, az), 1, lx1)
968 : END DO
969 :
970 28486378 : DO ax = 2, la
971 8518926 : f3 = f2*REAL(ax - 1, dp)
972 25556778 : DO ay = 0, la - ax
973 9350940 : az = la - ax - ay
974 : s(coset(ax, ay, az), 1, l) = rap(1)*s(coset(ax - 1, ay, az), 1, l) + &
975 9350940 : f3*s(coset(ax - 2, ay, az), 1, l)
976 9350940 : IF (lx1 > 0) s(coset(ax, ay, az), 1, l) = s(coset(ax, ay, az), 1, l) + &
977 12224763 : f2x*s(coset(ax - 1, ay, az), 1, lx1)
978 : END DO
979 : END DO
980 :
981 : END DO
982 :
983 : ! *** Recurrence steps: [a|M|s] -> [a|M|b] ***
984 :
985 12280540 : IF (lb_max > 0) THEN
986 :
987 97110248 : DO j = 2, ncoset(lb_max)
988 877833488 : DO i = 1, ncoset(la_max)
989 865870592 : s(i, j, l) = 0.0_dp
990 : END DO
991 : END DO
992 :
993 : ! *** Horizontal recurrence steps ***
994 :
995 47851584 : rbp(:) = rap(:) - rab(:)
996 :
997 : ! *** [a|M|p] = [a+1i|M|s] - (Bi - Ai)*[a|M|s] ***
998 :
999 11962896 : IF (lb_max == 1) THEN
1000 5138462 : la_start = la_min
1001 : ELSE
1002 6824434 : la_start = MAX(0, la_min - 1)
1003 : END IF
1004 :
1005 31398210 : DO la = la_start, la_max - 1
1006 59292007 : DO ax = 0, la
1007 84512589 : DO ay = 0, la - ax
1008 37183478 : az = la - ax - ay
1009 : s(coset(ax, ay, az), 2, l) = s(coset(ax + 1, ay, az), 1, l) - &
1010 37183478 : rab(1)*s(coset(ax, ay, az), 1, l)
1011 : s(coset(ax, ay, az), 3, l) = s(coset(ax, ay + 1, az), 1, l) - &
1012 37183478 : rab(2)*s(coset(ax, ay, az), 1, l)
1013 : s(coset(ax, ay, az), 4, l) = s(coset(ax, ay, az + 1), 1, l) - &
1014 65077275 : rab(3)*s(coset(ax, ay, az), 1, l)
1015 : END DO
1016 : END DO
1017 : END DO
1018 :
1019 : ! *** Vertical recurrence step ***
1020 :
1021 : ! *** [a|M|p] = (Pi - Bi)*[a|M|s] + f2*Ni(a)*[a-1i|M|s] ***
1022 : ! *** + f2*Ni(m)*[a|M-1i|s] ***
1023 :
1024 43547522 : DO ax = 0, la_max
1025 31584626 : fx = f2*REAL(ax, dp)
1026 103244198 : DO ay = 0, la_max - ax
1027 59696676 : fy = f2*REAL(ay, dp)
1028 59696676 : az = la_max - ax - ay
1029 59696676 : fz = f2*REAL(az, dp)
1030 59696676 : IF (ax == 0) THEN
1031 31584626 : s(coset(ax, ay, az), 2, l) = rbp(1)*s(coset(ax, ay, az), 1, l)
1032 : ELSE
1033 : s(coset(ax, ay, az), 2, l) = rbp(1)*s(coset(ax, ay, az), 1, l) + &
1034 28112050 : fx*s(coset(ax - 1, ay, az), 1, l)
1035 : END IF
1036 59696676 : IF (lx1 > 0) s(coset(ax, ay, az), 2, l) = s(coset(ax, ay, az), 2, l) + &
1037 23484123 : f2x*s(coset(ax, ay, az), 1, lx1)
1038 59696676 : IF (ay == 0) THEN
1039 31584626 : s(coset(ax, ay, az), 3, l) = rbp(2)*s(coset(ax, ay, az), 1, l)
1040 : ELSE
1041 : s(coset(ax, ay, az), 3, l) = rbp(2)*s(coset(ax, ay, az), 1, l) + &
1042 28112050 : fy*s(coset(ax, ay - 1, az), 1, l)
1043 : END IF
1044 59696676 : IF (ly1 > 0) s(coset(ax, ay, az), 3, l) = s(coset(ax, ay, az), 3, l) + &
1045 23484123 : f2y*s(coset(ax, ay, az), 1, ly1)
1046 59696676 : IF (az == 0) THEN
1047 31584626 : s(coset(ax, ay, az), 4, l) = rbp(3)*s(coset(ax, ay, az), 1, l)
1048 : ELSE
1049 : s(coset(ax, ay, az), 4, l) = rbp(3)*s(coset(ax, ay, az), 1, l) + &
1050 28112050 : fz*s(coset(ax, ay, az - 1), 1, l)
1051 : END IF
1052 59696676 : IF (lz1 > 0) s(coset(ax, ay, az), 4, l) = s(coset(ax, ay, az), 4, l) + &
1053 55068749 : f2z*s(coset(ax, ay, az), 1, lz1)
1054 : END DO
1055 : END DO
1056 :
1057 : ! *** Recurrence steps: [a|M|p] -> [a|M|b] ***
1058 :
1059 19618536 : DO lb = 2, lb_max
1060 :
1061 : ! *** Horizontal recurrence steps ***
1062 :
1063 : ! *** [a|M|b] = [a+1i|M|b-1i] - (Bi - Ai)*[a|M|b-1i] ***
1064 :
1065 7655640 : IF (lb == lb_max) THEN
1066 6824434 : la_start = la_min
1067 : ELSE
1068 831206 : la_start = MAX(0, la_min - 1)
1069 : END IF
1070 :
1071 21558852 : DO la = la_start, la_max - 1
1072 43273880 : DO ax = 0, la
1073 65958818 : DO ay = 0, la - ax
1074 30340578 : az = la - ax - ay
1075 :
1076 : ! *** Shift of angular momentum component z from a to b ***
1077 :
1078 : s(coset(ax, ay, az), coset(0, 0, lb), l) = &
1079 : s(coset(ax, ay, az + 1), coset(0, 0, lb - 1), l) - &
1080 30340578 : rab(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1), l)
1081 :
1082 : ! *** Shift of angular momentum component y from a to b ***
1083 :
1084 94429112 : DO by = 1, lb
1085 64088534 : bz = lb - by
1086 : s(coset(ax, ay, az), coset(0, by, bz), l) = &
1087 : s(coset(ax, ay + 1, az), coset(0, by - 1, bz), l) - &
1088 94429112 : rab(2)*s(coset(ax, ay, az), coset(0, by - 1, bz), l)
1089 : END DO
1090 :
1091 : ! *** Shift of angular momentum component x from a to b ***
1092 :
1093 116144140 : DO bx = 1, lb
1094 195672980 : DO by = 0, lb - bx
1095 101243868 : bz = lb - bx - by
1096 : s(coset(ax, ay, az), coset(bx, by, bz), l) = &
1097 : s(coset(ax + 1, ay, az), coset(bx - 1, by, bz), l) - &
1098 165332402 : rab(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz), l)
1099 : END DO
1100 : END DO
1101 :
1102 : END DO
1103 : END DO
1104 : END DO
1105 :
1106 : ! *** Vertical recurrence step ***
1107 :
1108 : ! *** [a|M|b] = (Pi - Bi)*[a|M|b-1i] + f2*Ni(a)*[a-1i|M|b-1i] + ***
1109 : ! *** f2*Ni(b-1i)*[a|M|b-2i] + f2*Ni(m)[a|M-1i|b-1i] ***
1110 :
1111 41935015 : DO ax = 0, la_max
1112 22316479 : fx = f2*REAL(ax, dp)
1113 74768514 : DO ay = 0, la_max - ax
1114 44796395 : fy = f2*REAL(ay, dp)
1115 44796395 : az = la_max - ax - ay
1116 44796395 : fz = f2*REAL(az, dp)
1117 :
1118 : ! *** Shift of angular momentum component z from a to b ***
1119 :
1120 44796395 : f3 = f2*REAL(lb - 1, dp)
1121 :
1122 44796395 : IF (az == 0) THEN
1123 : s(coset(ax, ay, az), coset(0, 0, lb), l) = &
1124 : rbp(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1), l) + &
1125 22316479 : f3*s(coset(ax, ay, az), coset(0, 0, lb - 2), l)
1126 : ELSE
1127 : s(coset(ax, ay, az), coset(0, 0, lb), l) = &
1128 : rbp(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1), l) + &
1129 : fz*s(coset(ax, ay, az - 1), coset(0, 0, lb - 1), l) + &
1130 22479916 : f3*s(coset(ax, ay, az), coset(0, 0, lb - 2), l)
1131 : END IF
1132 44796395 : IF (lz1 > 0) s(coset(ax, ay, az), coset(0, 0, lb), l) = &
1133 : s(coset(ax, ay, az), coset(0, 0, lb), l) + &
1134 17793257 : f2z*s(coset(ax, ay, az), coset(0, 0, lb - 1), lz1)
1135 :
1136 : ! *** Shift of angular momentum component y from a to b ***
1137 :
1138 44796395 : IF (ay == 0) THEN
1139 22316479 : bz = lb - 1
1140 : s(coset(ax, ay, az), coset(0, 1, bz), l) = &
1141 22316479 : rbp(2)*s(coset(ax, ay, az), coset(0, 0, bz), l)
1142 22316479 : IF (ly1 > 0) s(coset(ax, ay, az), coset(0, 1, bz), l) = &
1143 : s(coset(ax, ay, az), coset(0, 1, bz), l) + &
1144 8859592 : f2y*s(coset(ax, ay, az), coset(0, 0, bz), ly1)
1145 47109156 : DO by = 2, lb
1146 24792677 : bz = lb - by
1147 24792677 : f3 = f2*REAL(by - 1, dp)
1148 : s(coset(ax, ay, az), coset(0, by, bz), l) = &
1149 : rbp(2)*s(coset(ax, ay, az), coset(0, by - 1, bz), l) + &
1150 24792677 : f3*s(coset(ax, ay, az), coset(0, by - 2, bz), l)
1151 24792677 : IF (ly1 > 0) s(coset(ax, ay, az), coset(0, by, bz), l) = &
1152 : s(coset(ax, ay, az), coset(0, by, bz), l) + &
1153 32160219 : f2y*s(coset(ax, ay, az), coset(0, by - 1, bz), ly1)
1154 : END DO
1155 : ELSE
1156 22479916 : bz = lb - 1
1157 : s(coset(ax, ay, az), coset(0, 1, bz), l) = &
1158 : rbp(2)*s(coset(ax, ay, az), coset(0, 0, bz), l) + &
1159 22479916 : fy*s(coset(ax, ay - 1, az), coset(0, 0, bz), l)
1160 22479916 : IF (ly1 > 0) s(coset(ax, ay, az), coset(0, 1, bz), l) = &
1161 : s(coset(ax, ay, az), coset(0, 1, bz), l) + &
1162 8933665 : f2y*s(coset(ax, ay, az), coset(0, 0, bz), ly1)
1163 47483978 : DO by = 2, lb
1164 25004062 : bz = lb - by
1165 25004062 : f3 = f2*REAL(by - 1, dp)
1166 : s(coset(ax, ay, az), coset(0, by, bz), l) = &
1167 : rbp(2)*s(coset(ax, ay, az), coset(0, by - 1, bz), l) + &
1168 : fy*s(coset(ax, ay - 1, az), coset(0, by - 1, bz), l) + &
1169 25004062 : f3*s(coset(ax, ay, az), coset(0, by - 2, bz), l)
1170 25004062 : IF (ly1 > 0) s(coset(ax, ay, az), coset(0, by, bz), l) = &
1171 : s(coset(ax, ay, az), coset(0, by, bz), l) + &
1172 32415935 : f2y*s(coset(ax, ay, az), coset(0, by - 1, bz), ly1)
1173 : END DO
1174 : END IF
1175 :
1176 : ! *** Shift of angular momentum component x from a to b ***
1177 :
1178 67112874 : IF (ax == 0) THEN
1179 69425635 : DO by = 0, lb - 1
1180 47109156 : bz = lb - 1 - by
1181 : s(coset(ax, ay, az), coset(1, by, bz), l) = &
1182 47109156 : rbp(1)*s(coset(ax, ay, az), coset(0, by, bz), l)
1183 47109156 : IF (lx1 > 0) s(coset(ax, ay, az), coset(1, by, bz), l) = &
1184 : s(coset(ax, ay, az), coset(1, by, bz), l) + &
1185 41019811 : f2x*s(coset(ax, ay, az), coset(0, by, bz), lx1)
1186 : END DO
1187 47109156 : DO bx = 2, lb
1188 24792677 : f3 = f2*REAL(bx - 1, dp)
1189 74378031 : DO by = 0, lb - bx
1190 27268875 : bz = lb - bx - by
1191 : s(coset(ax, ay, az), coset(bx, by, bz), l) = &
1192 : rbp(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz), l) + &
1193 27268875 : f3*s(coset(ax, ay, az), coset(bx - 2, by, bz), l)
1194 27268875 : IF (lx1 > 0) s(coset(ax, ay, az), coset(bx, by, bz), l) = &
1195 : s(coset(ax, ay, az), coset(bx, by, bz), l) + &
1196 35620565 : f2x*s(coset(ax, ay, az), coset(bx - 1, by, bz), lx1)
1197 : END DO
1198 : END DO
1199 : ELSE
1200 69963894 : DO by = 0, lb - 1
1201 47483978 : bz = lb - 1 - by
1202 : s(coset(ax, ay, az), coset(1, by, bz), l) = &
1203 : rbp(1)*s(coset(ax, ay, az), coset(0, by, bz), l) + &
1204 47483978 : fx*s(coset(ax - 1, ay, az), coset(0, by, bz), l)
1205 47483978 : IF (lx1 > 0) s(coset(ax, ay, az), coset(1, by, bz), l) = &
1206 : s(coset(ax, ay, az), coset(1, by, bz), l) + &
1207 41349600 : f2x*s(coset(ax, ay, az), coset(0, by, bz), lx1)
1208 : END DO
1209 47483978 : DO bx = 2, lb
1210 25004062 : f3 = f2*REAL(bx - 1, dp)
1211 75012186 : DO by = 0, lb - bx
1212 27528208 : bz = lb - bx - by
1213 : s(coset(ax, ay, az), coset(bx, by, bz), l) = &
1214 : rbp(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz), l) + &
1215 : fx*s(coset(ax - 1, ay, az), coset(bx - 1, by, bz), l) + &
1216 27528208 : f3*s(coset(ax, ay, az), coset(bx - 2, by, bz), l)
1217 27528208 : IF (lx1 > 0) s(coset(ax, ay, az), coset(bx, by, bz), l) = &
1218 : s(coset(ax, ay, az), coset(bx, by, bz), l) + &
1219 35942435 : f2x*s(coset(ax, ay, az), coset(bx - 1, by, bz), lx1)
1220 : END DO
1221 : END DO
1222 : END IF
1223 :
1224 : END DO
1225 : END DO
1226 :
1227 : END DO
1228 :
1229 : END IF
1230 :
1231 : ELSE
1232 :
1233 616965 : IF (lb_max > 0) THEN
1234 :
1235 : ! *** Vertical recurrence steps: [s|M|s] -> [s|M|b] ***
1236 :
1237 843984 : rbp(:) = (f1 - 1.0_dp)*rab(:)
1238 :
1239 : ! *** [s|M|p] = (Pi - Bi)*[s|M|s] + f2*Ni(m)*[s|M-1i|s] ***
1240 :
1241 210996 : s(1, 2, l) = rbp(1)*s(1, 1, l)
1242 210996 : s(1, 3, l) = rbp(2)*s(1, 1, l)
1243 210996 : s(1, 4, l) = rbp(3)*s(1, 1, l)
1244 210996 : IF (lx1 > 0) s(1, 2, l) = s(1, 2, l) + f2x*s(1, 1, lx1)
1245 210996 : IF (ly1 > 0) s(1, 3, l) = s(1, 3, l) + f2y*s(1, 1, ly1)
1246 210996 : IF (lz1 > 0) s(1, 4, l) = s(1, 4, l) + f2z*s(1, 1, lz1)
1247 :
1248 : ! *** [s|M|b] = (Pi - Bi)*[s|M|b-1i] + f2*Ni(b-1i)*[s|M|b-2i] ***
1249 : ! *** + f2*Ni(m)*[s|M-1i|b-1i] ***
1250 :
1251 237656 : DO lb = 2, lb_max
1252 :
1253 : ! *** Increase the angular momentum component z of function b ***
1254 :
1255 : s(1, coset(0, 0, lb), l) = rbp(3)*s(1, coset(0, 0, lb - 1), l) + &
1256 26660 : f2*REAL(lb - 1, dp)*s(1, coset(0, 0, lb - 2), l)
1257 26660 : IF (lz1 > 0) s(1, coset(0, 0, lb), l) = s(1, coset(0, 0, lb), l) + &
1258 6821 : f2z*s(1, coset(0, 0, lb - 1), lz1)
1259 :
1260 : ! *** Increase the angular momentum component y of function b ***
1261 :
1262 26660 : bz = lb - 1
1263 26660 : s(1, coset(0, 1, bz), l) = rbp(2)*s(1, coset(0, 0, bz), l)
1264 26660 : IF (ly1 > 0) s(1, coset(0, 1, bz), l) = s(1, coset(0, 1, bz), l) + &
1265 6821 : f2y*s(1, coset(0, 0, bz), ly1)
1266 :
1267 53848 : DO by = 2, lb
1268 27188 : bz = lb - by
1269 : s(1, coset(0, by, bz), l) = rbp(2)*s(1, coset(0, by - 1, bz), l) + &
1270 27188 : f2*REAL(by - 1, dp)*s(1, coset(0, by - 2, bz), l)
1271 27188 : IF (ly1 > 0) s(1, coset(0, by, bz), l) = s(1, coset(0, by, bz), l) + &
1272 33613 : f2y*s(1, coset(0, by - 1, bz), ly1)
1273 : END DO
1274 :
1275 : ! *** Increase the angular momentum component x of function b ***
1276 :
1277 80508 : DO by = 0, lb - 1
1278 53848 : bz = lb - 1 - by
1279 53848 : s(1, coset(1, by, bz), l) = rbp(1)*s(1, coset(0, by, bz), l)
1280 53848 : IF (lx1 > 0) s(1, coset(1, by, bz), l) = s(1, coset(1, by, bz), l) + &
1281 40434 : f2x*s(1, coset(0, by, bz), lx1)
1282 : END DO
1283 :
1284 264844 : DO bx = 2, lb
1285 27188 : f3 = f2*REAL(bx - 1, dp)
1286 81564 : DO by = 0, lb - bx
1287 27716 : bz = lb - bx - by
1288 : s(1, coset(bx, by, bz), l) = rbp(1)*s(1, coset(bx - 1, by, bz), l) + &
1289 27716 : f3*s(1, coset(bx - 2, by, bz), l)
1290 27716 : IF (lx1 > 0) s(1, coset(bx, by, bz), l) = s(1, coset(bx, by, bz), l) + &
1291 34273 : f2x*s(1, coset(bx - 1, by, bz), lx1)
1292 : END DO
1293 : END DO
1294 :
1295 : END DO
1296 :
1297 : END IF
1298 :
1299 : END IF
1300 :
1301 : END DO
1302 :
1303 12897505 : DO k = 2, ncoset(lc_max)
1304 101040675 : DO j = 1, ncoset(lb_max)
1305 888160476 : DO i = 1, ncoset(la_max)
1306 876772480 : mab(na + i, nb + j, k - 1) = s(i, j, k)
1307 : END DO
1308 : END DO
1309 : END DO
1310 :
1311 2979479 : nb = nb + ncoset(lb_max)
1312 :
1313 : END DO
1314 :
1315 2118772 : na = na + ncoset(la_max)
1316 :
1317 : END DO
1318 :
1319 648802 : END SUBROUTINE moment
1320 :
1321 : ! **************************************************************************************************
1322 : !> \brief This returns the derivative of the moment integrals [a|\mu|b].
1323 : !> By default, it differentiates the primitive on the right:
1324 : !> [a|\mu|d/dR_bi] = 2*zetb*[a|\mu|b+1i] - Ni(b)[a|\mu|b-1i]
1325 : !> A weighted derivative combines the left and right primitive derivatives
1326 : !> using deltaR for the corresponding atom centers.
1327 : !> order indicates the max order of the moment operator to be calculated
1328 : !> 1: dipole
1329 : !> 2: quadrupole
1330 : !> ...
1331 : !> \param la_max ...
1332 : !> \param npgfa ...
1333 : !> \param zeta ...
1334 : !> \param rpgfa ...
1335 : !> \param la_min ...
1336 : !> \param lb_max ...
1337 : !> \param npgfb ...
1338 : !> \param zetb ...
1339 : !> \param rpgfb ...
1340 : !> \param lb_min ...
1341 : !> \param order ...
1342 : !> \param rac ...
1343 : !> \param rbc ...
1344 : !> \param difmab ...
1345 : !> \param mab_ext ...
1346 : !> \param deltaR optional weights for the left and right primitive derivatives
1347 : !> \param lambda optional atom selector for the factor (iatom == lambda) - (jatom == lambda)
1348 : !> \param iatom atom associated with the left basis function
1349 : !> \param jatom atom associated with the right basis function
1350 : !> \note
1351 : ! **************************************************************************************************
1352 600231 : SUBROUTINE diff_momop(la_max, npgfa, zeta, rpgfa, la_min, &
1353 600231 : lb_max, npgfb, zetb, rpgfb, lb_min, &
1354 600231 : order, rac, rbc, difmab, mab_ext, deltaR, lambda, iatom, jatom)
1355 :
1356 : INTEGER, INTENT(IN) :: la_max, npgfa
1357 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
1358 : INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
1359 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
1360 : INTEGER, INTENT(IN) :: lb_min, order
1361 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
1362 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(OUT) :: difmab
1363 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
1364 : POINTER :: mab_ext
1365 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: deltaR
1366 : INTEGER, INTENT(IN), OPTIONAL :: lambda, iatom, jatom
1367 :
1368 : INTEGER :: ider, imom, lda, lda_min, ldb, ldb_min
1369 : REAL(KIND=dp) :: dab, rab(3), rab2
1370 600231 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: difmab_tmp
1371 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mab
1372 :
1373 2400924 : rab = rbc - rac
1374 2400924 : rab2 = SUM(rab**2)
1375 600231 : dab = SQRT(rab2)
1376 :
1377 600231 : lda_min = MAX(0, la_min - 1)
1378 600231 : ldb_min = MAX(0, lb_min - 1)
1379 600231 : lda = ncoset(la_max)*npgfa
1380 600231 : ldb = ncoset(lb_max)*npgfb
1381 2974323 : ALLOCATE (difmab_tmp(lda, ldb, 3))
1382 :
1383 600231 : IF (PRESENT(mab_ext)) THEN
1384 600231 : mab => mab_ext
1385 : ELSE
1386 : ALLOCATE (mab(npgfa*ncoset(la_max + 1), npgfb*ncoset(lb_max + 1), &
1387 0 : ncoset(order) - 1))
1388 0 : mab = 0.0_dp
1389 : ! *** Calculate the primitive overlap integrals ***
1390 : CALL moment(la_max + 1, npgfa, zeta, rpgfa, lda_min, &
1391 : lb_max + 1, npgfb, zetb, rpgfb, &
1392 0 : order, rac, rbc, mab)
1393 :
1394 : END IF
1395 5982852 : DO imom = 1, ncoset(order) - 1
1396 : CALL adbdr(la_max, npgfa, rpgfa, la_min, &
1397 : lb_max, npgfb, zetb, rpgfb, lb_min, &
1398 : dab, mab(:, :, imom), difmab_tmp(:, :, 1), &
1399 5382621 : difmab_tmp(:, :, 2), difmab_tmp(:, :, 3))
1400 5982852 : IF (PRESENT(deltaR)) THEN
1401 6804 : CPASSERT(ASSOCIATED(deltaR))
1402 6804 : CPASSERT(PRESENT(iatom) .AND. PRESENT(jatom))
1403 27216 : DO ider = 1, 3
1404 : difmab(1:lda, 1:ldb, imom, ider) = &
1405 1496880 : difmab_tmp(1:lda, 1:ldb, ider)*deltaR(ider, jatom)
1406 : END DO
1407 :
1408 : CALL dabdr(la_max, npgfa, zeta, rpgfa, la_min, &
1409 : lb_max, npgfb, rpgfb, lb_min, &
1410 : dab, mab(:, :, imom), difmab_tmp(:, :, 1), &
1411 6804 : difmab_tmp(:, :, 2), difmab_tmp(:, :, 3))
1412 27216 : DO ider = 1, 3
1413 : difmab(1:lda, 1:ldb, imom, ider) = difmab(1:lda, 1:ldb, imom, ider) &
1414 1496880 : + difmab_tmp(1:lda, 1:ldb, ider)*deltaR(ider, iatom)
1415 : END DO
1416 : ELSE
1417 1206943560 : difmab(1:lda, 1:ldb, imom, :) = difmab_tmp(1:lda, 1:ldb, :)
1418 : END IF
1419 : END DO
1420 :
1421 600231 : IF (PRESENT(lambda)) THEN
1422 27 : CPASSERT(.NOT. PRESENT(deltaR))
1423 27 : CPASSERT(PRESENT(iatom) .AND. PRESENT(jatom))
1424 27 : IF (iatom == lambda .AND. jatom == lambda) THEN
1425 7383 : difmab = 0.0_dp
1426 24 : ELSE IF (iatom == lambda) THEN
1427 : ! The right-hand derivative is selected with a positive sign.
1428 18 : ELSE IF (jatom == lambda) THEN
1429 14766 : difmab = -difmab
1430 : ELSE
1431 29532 : difmab = 0.0_dp
1432 : END IF
1433 : END IF
1434 :
1435 600231 : IF (PRESENT(mab_ext)) THEN
1436 : NULLIFY (mab)
1437 : ELSE
1438 0 : DEALLOCATE (mab)
1439 : END IF
1440 600231 : DEALLOCATE (difmab_tmp)
1441 :
1442 600231 : END SUBROUTINE diff_momop
1443 :
1444 : ! **************************************************************************************************
1445 : !> \brief This returns the derivative of the dipole integrals [a|x|b], with respect
1446 : !> to the position of the primitive on the left and right, i.e.
1447 : !> [da/dR_ai|\mu|b] = 2*zeta*[a+1i|\mu|b] - Ni(a)[a-1i|\mu|b]
1448 : !> \param la_max ...
1449 : !> \param npgfa ...
1450 : !> \param zeta ...
1451 : !> \param rpgfa ...
1452 : !> \param la_min ...
1453 : !> \param lb_max ...
1454 : !> \param npgfb ...
1455 : !> \param zetb ...
1456 : !> \param rpgfb ...
1457 : !> \param lb_min ...
1458 : !> \param order ...
1459 : !> \param rac ...
1460 : !> \param rbc ...
1461 : !> \param pab ...
1462 : !> \param forcea ...
1463 : !> \param forceb ...
1464 : !> \note
1465 : ! **************************************************************************************************
1466 2124 : SUBROUTINE dipole_force(la_max, npgfa, zeta, rpgfa, la_min, &
1467 2124 : lb_max, npgfb, zetb, rpgfb, lb_min, &
1468 2124 : order, rac, rbc, pab, forcea, forceb)
1469 :
1470 : INTEGER, INTENT(IN) :: la_max, npgfa
1471 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
1472 : INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
1473 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
1474 : INTEGER, INTENT(IN) :: lb_min, order
1475 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
1476 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: pab
1477 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: forcea, forceb
1478 :
1479 : INTEGER :: i, imom, ipgf, j, jpgf, lda, lda_min, &
1480 : ldb, ldb_min, na, nb
1481 : REAL(KIND=dp) :: dab, rab(3), rab2
1482 2124 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: difmab, mab
1483 :
1484 2124 : CPASSERT(order == 1)
1485 : MARK_USED(order)
1486 :
1487 8496 : rab = rbc - rac
1488 8496 : rab2 = SUM(rab**2)
1489 2124 : dab = SQRT(rab2)
1490 :
1491 2124 : lda_min = MAX(0, la_min - 1)
1492 2124 : ldb_min = MAX(0, lb_min - 1)
1493 2124 : lda = ncoset(la_max)*npgfa
1494 2124 : ldb = ncoset(lb_max)*npgfb
1495 10620 : ALLOCATE (difmab(lda, ldb, 3))
1496 10620 : ALLOCATE (mab(npgfa*ncoset(la_max + 1), npgfb*ncoset(lb_max + 1), 3))
1497 2124 : mab = 0.0_dp
1498 : CALL moment(la_max + 1, npgfa, zeta, rpgfa, lda_min, &
1499 2124 : lb_max + 1, npgfb, zetb, rpgfb, 1, rac, rbc, mab)
1500 :
1501 8496 : DO imom = 1, 3
1502 6372 : difmab = 0.0_dp
1503 : CALL adbdr(la_max, npgfa, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, &
1504 6372 : dab, mab(:, :, imom), difmab(:, :, 1), difmab(:, :, 2), difmab(:, :, 3))
1505 6372 : na = 0
1506 24360 : DO ipgf = 1, npgfa
1507 : nb = 0
1508 69429 : DO jpgf = 1, npgfb
1509 171363 : DO j = nb + ncoset(lb_min - 1) + 1, nb + ncoset(lb_max)
1510 518880 : DO i = na + ncoset(la_min - 1) + 1, na + ncoset(la_max)
1511 347517 : forceb(imom, 1) = forceb(imom, 1) + pab(i, j)*difmab(i, j, 1)
1512 347517 : forceb(imom, 2) = forceb(imom, 2) + pab(i, j)*difmab(i, j, 2)
1513 467439 : forceb(imom, 3) = forceb(imom, 3) + pab(i, j)*difmab(i, j, 3)
1514 : END DO
1515 : END DO
1516 69429 : nb = nb + ncoset(lb_max)
1517 : END DO
1518 24360 : na = na + ncoset(la_max)
1519 : END DO
1520 :
1521 6372 : difmab = 0.0_dp
1522 : CALL dabdr(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, rpgfb, lb_min, &
1523 6372 : dab, mab(:, :, imom), difmab(:, :, 1), difmab(:, :, 2), difmab(:, :, 3))
1524 6372 : na = 0
1525 26484 : DO ipgf = 1, npgfa
1526 : nb = 0
1527 69429 : DO jpgf = 1, npgfb
1528 171363 : DO j = nb + ncoset(lb_min - 1) + 1, nb + ncoset(lb_max)
1529 518880 : DO i = na + ncoset(la_min - 1) + 1, na + ncoset(la_max)
1530 347517 : forcea(imom, 1) = forcea(imom, 1) + pab(i, j)*difmab(i, j, 1)
1531 347517 : forcea(imom, 2) = forcea(imom, 2) + pab(i, j)*difmab(i, j, 2)
1532 467439 : forcea(imom, 3) = forcea(imom, 3) + pab(i, j)*difmab(i, j, 3)
1533 : END DO
1534 : END DO
1535 69429 : nb = nb + ncoset(lb_max)
1536 : END DO
1537 24360 : na = na + ncoset(la_max)
1538 : END DO
1539 : END DO
1540 :
1541 2124 : DEALLOCATE (mab, difmab)
1542 :
1543 2124 : END SUBROUTINE dipole_force
1544 :
1545 : END MODULE ai_moments
|
