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