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