(git:3919fb7)
Loading...
Searching...
No Matches
ai_moments.F
Go to the documentation of this file.
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! **************************************************************************************************
17MODULE 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
49CONTAINS
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 SUBROUTINE contract_cossin(cos_block, sin_block, &
74 iatom, ncoa, nsgfa, sgfa, sphi_a, ldsa, &
75 jatom, ncob, nsgfb, sgfb, sphi_b, ldsb, &
76 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 0.0_dp, work(1, 1), ldwork)
95
96 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 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 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 0.0_dp, work(1, 1), ldwork)
115
116 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 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 SIZE(sin_block, 1))
128 END IF
129
130 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 SUBROUTINE cossin(la_max_set, npgfa, zeta, rpgfa, la_min_set, &
153 lb_max, npgfb, zetb, rpgfb, lb_min, &
154 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), & ncoset(lb_max), 3) :: dscos, dssin
173 REAL(kind=dp), &
174 DIMENSION(ncoset(la_max_set+1), ncoset(lb_max)) :: sc, ss
175
176 rab = rbc - rac
177 rab2 = sum(rab**2)
178 dab = sqrt(rab2)
179 k2 = kvec(1)*kvec(1) + kvec(2)*kvec(2) + kvec(3)*kvec(3)
180
181 IF (PRESENT(dcosab)) THEN
182 da_max = 1
183 la_max = la_max_set + 1
184 la_min = max(0, la_min_set - 1)
185 dscos = 0.0_dp
186 dssin = 0.0_dp
187 ELSE
188 da_max = 0
189 la_max = la_max_set
190 la_min = la_min_set
191 END IF
192
193 ! initialize all matrix elements to zero
194 IF (PRESENT(dcosab)) THEN
195 na = ncoset(la_max - 1)*npgfa
196 ELSE
197 na = ncoset(la_max)*npgfa
198 END IF
199 nb = ncoset(lb_max)*npgfb
200 cosab(1:na, 1:nb) = 0.0_dp
201 sinab(1:na, 1:nb) = 0.0_dp
202 IF (PRESENT(dcosab)) THEN
203 dcosab(1:na, 1:nb, :) = 0.0_dp
204 dsinab(1:na, 1:nb, :) = 0.0_dp
205 END IF
206! *** Loop over all pairs of primitive Gaussian-type functions ***
207
208 na = 0
209 DO ipgf = 1, npgfa
210
211 nb = 0
212
213 DO jpgf = 1, npgfb
214
215 ss = 0.0_dp
216 sc = 0.0_dp
217
218! *** Screening ***
219 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
220 nb = nb + ncoset(lb_max)
221 cycle
222 END IF
223
224! *** Calculate some prefactors ***
225
226 zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
227
228 f0 = (pi*zetp)**1.5_dp
229 f1 = zetb(jpgf)*zetp
230 f2 = 0.5_dp*zetp
231
232 kdp = zetp*dot_product(kvec, zeta(ipgf)*rac + zetb(jpgf)*rbc)
233
234! *** Calculate the basic two-center cos/sin integral [s|cos/sin|s] ***
235
236 s = f0*exp(-zeta(ipgf)*f1*rab2)*exp(-0.25_dp*k2*zetp)
237 sc(1, 1) = s*cos(kdp)
238 ss(1, 1) = s*sin(kdp)
239
240! *** Recurrence steps: [s|O|s] -> [a|O|b] ***
241
242 IF (la_max > 0) THEN
243
244! *** Vertical recurrence steps: [s|O|s] -> [a|O|s] ***
245
246 rap(:) = f1*rab(:)
247
248! *** [p|O|s] = (Pi - Ai)*[s|O|s] +[s|dO|s] (i = x,y,z) ***
249
250 sc(2, 1) = rap(1)*sc(1, 1) - f2*kvec(1)*ss(1, 1)
251 sc(3, 1) = rap(2)*sc(1, 1) - f2*kvec(2)*ss(1, 1)
252 sc(4, 1) = rap(3)*sc(1, 1) - f2*kvec(3)*ss(1, 1)
253 ss(2, 1) = rap(1)*ss(1, 1) + f2*kvec(1)*sc(1, 1)
254 ss(3, 1) = rap(2)*ss(1, 1) + f2*kvec(2)*sc(1, 1)
255 ss(4, 1) = rap(3)*ss(1, 1) + f2*kvec(3)*sc(1, 1)
256
257! *** [a|O|s] = (Pi - Ai)*[a-1i|O|s] + f2*Ni(a-1i)*[a-2i|s] ***
258! *** + [a-1i|dO|s] ***
259
260 DO la = 2, la_max
261
262! *** Increase the angular momentum component z of function a ***
263
264 sc(coset(0, 0, la), 1) = rap(3)*sc(coset(0, 0, la - 1), 1) + &
265 f2*real(la - 1, dp)*sc(coset(0, 0, la - 2), 1) - &
266 f2*kvec(3)*ss(coset(0, 0, la - 1), 1)
267 ss(coset(0, 0, la), 1) = rap(3)*ss(coset(0, 0, la - 1), 1) + &
268 f2*real(la - 1, dp)*ss(coset(0, 0, la - 2), 1) + &
269 f2*kvec(3)*sc(coset(0, 0, la - 1), 1)
270
271! *** Increase the angular momentum component y of function a ***
272
273 az = la - 1
274 sc(coset(0, 1, az), 1) = rap(2)*sc(coset(0, 0, az), 1) - &
275 f2*kvec(2)*ss(coset(0, 0, az), 1)
276 ss(coset(0, 1, az), 1) = rap(2)*ss(coset(0, 0, az), 1) + &
277 f2*kvec(2)*sc(coset(0, 0, az), 1)
278
279 DO ay = 2, la
280 az = la - ay
281 sc(coset(0, ay, az), 1) = rap(2)*sc(coset(0, ay - 1, az), 1) + &
282 f2*real(ay - 1, dp)*sc(coset(0, ay - 2, az), 1) - &
283 f2*kvec(2)*ss(coset(0, ay - 1, az), 1)
284 ss(coset(0, ay, az), 1) = rap(2)*ss(coset(0, ay - 1, az), 1) + &
285 f2*real(ay - 1, dp)*ss(coset(0, ay - 2, az), 1) + &
286 f2*kvec(2)*sc(coset(0, ay - 1, az), 1)
287 END DO
288
289! *** Increase the angular momentum component x of function a ***
290
291 DO ay = 0, la - 1
292 az = la - 1 - ay
293 sc(coset(1, ay, az), 1) = rap(1)*sc(coset(0, ay, az), 1) - &
294 f2*kvec(1)*ss(coset(0, ay, az), 1)
295 ss(coset(1, ay, az), 1) = rap(1)*ss(coset(0, ay, az), 1) + &
296 f2*kvec(1)*sc(coset(0, ay, az), 1)
297 END DO
298
299 DO ax = 2, la
300 f3 = f2*real(ax - 1, dp)
301 DO ay = 0, la - ax
302 az = la - ax - ay
303 sc(coset(ax, ay, az), 1) = rap(1)*sc(coset(ax - 1, ay, az), 1) + &
304 f3*sc(coset(ax - 2, ay, az), 1) - &
305 f2*kvec(1)*ss(coset(ax - 1, ay, az), 1)
306 ss(coset(ax, ay, az), 1) = rap(1)*ss(coset(ax - 1, ay, az), 1) + &
307 f3*ss(coset(ax - 2, ay, az), 1) + &
308 f2*kvec(1)*sc(coset(ax - 1, ay, az), 1)
309 END DO
310 END DO
311
312 END DO
313
314! *** Recurrence steps: [a|O|s] -> [a|O|b] ***
315
316 IF (lb_max > 0) THEN
317
318 DO j = 2, ncoset(lb_max)
319 DO i = 1, ncoset(la_max)
320 sc(i, j) = 0.0_dp
321 ss(i, j) = 0.0_dp
322 END DO
323 END DO
324
325! *** Horizontal recurrence steps ***
326
327 rbp(:) = rap(:) - rab(:)
328
329! *** [a|O|p] = [a+1i|O|s] - (Bi - Ai)*[a|O|s] ***
330
331 IF (lb_max == 1) THEN
332 la_start = la_min
333 ELSE
334 la_start = max(0, la_min - 1)
335 END IF
336
337 DO la = la_start, la_max - 1
338 DO ax = 0, la
339 DO ay = 0, la - ax
340 az = la - ax - ay
341 sc(coset(ax, ay, az), 2) = sc(coset(ax + 1, ay, az), 1) - &
342 rab(1)*sc(coset(ax, ay, az), 1)
343 sc(coset(ax, ay, az), 3) = sc(coset(ax, ay + 1, az), 1) - &
344 rab(2)*sc(coset(ax, ay, az), 1)
345 sc(coset(ax, ay, az), 4) = sc(coset(ax, ay, az + 1), 1) - &
346 rab(3)*sc(coset(ax, ay, az), 1)
347 ss(coset(ax, ay, az), 2) = ss(coset(ax + 1, ay, az), 1) - &
348 rab(1)*ss(coset(ax, ay, az), 1)
349 ss(coset(ax, ay, az), 3) = ss(coset(ax, ay + 1, az), 1) - &
350 rab(2)*ss(coset(ax, ay, az), 1)
351 ss(coset(ax, ay, az), 4) = ss(coset(ax, ay, az + 1), 1) - &
352 rab(3)*ss(coset(ax, ay, az), 1)
353 END DO
354 END DO
355 END DO
356
357! *** Vertical recurrence step ***
358
359! *** [a|O|p] = (Pi - Bi)*[a|O|s] + f2*Ni(a)*[a-1i|O|s] ***
360! *** + [a|dO|s] ***
361
362 DO ax = 0, la_max
363 fx = f2*real(ax, dp)
364 DO ay = 0, la_max - ax
365 fy = f2*real(ay, dp)
366 az = la_max - ax - ay
367 fz = f2*real(az, dp)
368 IF (ax == 0) THEN
369 sc(coset(ax, ay, az), 2) = rbp(1)*sc(coset(ax, ay, az), 1) - &
370 f2*kvec(1)*ss(coset(ax, ay, az), 1)
371 ss(coset(ax, ay, az), 2) = rbp(1)*ss(coset(ax, ay, az), 1) + &
372 f2*kvec(1)*sc(coset(ax, ay, az), 1)
373 ELSE
374 sc(coset(ax, ay, az), 2) = rbp(1)*sc(coset(ax, ay, az), 1) + &
375 fx*sc(coset(ax - 1, ay, az), 1) - &
376 f2*kvec(1)*ss(coset(ax, ay, az), 1)
377 ss(coset(ax, ay, az), 2) = rbp(1)*ss(coset(ax, ay, az), 1) + &
378 fx*ss(coset(ax - 1, ay, az), 1) + &
379 f2*kvec(1)*sc(coset(ax, ay, az), 1)
380 END IF
381 IF (ay == 0) THEN
382 sc(coset(ax, ay, az), 3) = rbp(2)*sc(coset(ax, ay, az), 1) - &
383 f2*kvec(2)*ss(coset(ax, ay, az), 1)
384 ss(coset(ax, ay, az), 3) = rbp(2)*ss(coset(ax, ay, az), 1) + &
385 f2*kvec(2)*sc(coset(ax, ay, az), 1)
386 ELSE
387 sc(coset(ax, ay, az), 3) = rbp(2)*sc(coset(ax, ay, az), 1) + &
388 fy*sc(coset(ax, ay - 1, az), 1) - &
389 f2*kvec(2)*ss(coset(ax, ay, az), 1)
390 ss(coset(ax, ay, az), 3) = rbp(2)*ss(coset(ax, ay, az), 1) + &
391 fy*ss(coset(ax, ay - 1, az), 1) + &
392 f2*kvec(2)*sc(coset(ax, ay, az), 1)
393 END IF
394 IF (az == 0) THEN
395 sc(coset(ax, ay, az), 4) = rbp(3)*sc(coset(ax, ay, az), 1) - &
396 f2*kvec(3)*ss(coset(ax, ay, az), 1)
397 ss(coset(ax, ay, az), 4) = rbp(3)*ss(coset(ax, ay, az), 1) + &
398 f2*kvec(3)*sc(coset(ax, ay, az), 1)
399 ELSE
400 sc(coset(ax, ay, az), 4) = rbp(3)*sc(coset(ax, ay, az), 1) + &
401 fz*sc(coset(ax, ay, az - 1), 1) - &
402 f2*kvec(3)*ss(coset(ax, ay, az), 1)
403 ss(coset(ax, ay, az), 4) = rbp(3)*ss(coset(ax, ay, az), 1) + &
404 fz*ss(coset(ax, ay, az - 1), 1) + &
405 f2*kvec(3)*sc(coset(ax, ay, az), 1)
406 END IF
407 END DO
408 END DO
409
410! *** Recurrence steps: [a|O|p] -> [a|O|b] ***
411
412 DO lb = 2, lb_max
413
414! *** Horizontal recurrence steps ***
415
416! *** [a|O|b] = [a+1i|O|b-1i] - (Bi - Ai)*[a|O|b-1i] ***
417
418 IF (lb == lb_max) THEN
419 la_start = la_min
420 ELSE
421 la_start = max(0, la_min - 1)
422 END IF
423
424 DO la = la_start, la_max - 1
425 DO ax = 0, la
426 DO ay = 0, la - ax
427 az = la - ax - ay
428
429! *** Shift of angular momentum component z from a to b ***
430
431 sc(coset(ax, ay, az), coset(0, 0, lb)) = &
432 sc(coset(ax, ay, az + 1), coset(0, 0, lb - 1)) - &
433 rab(3)*sc(coset(ax, ay, az), coset(0, 0, lb - 1))
434 ss(coset(ax, ay, az), coset(0, 0, lb)) = &
435 ss(coset(ax, ay, az + 1), coset(0, 0, lb - 1)) - &
436 rab(3)*ss(coset(ax, ay, az), coset(0, 0, lb - 1))
437
438! *** Shift of angular momentum component y from a to b ***
439
440 DO by = 1, lb
441 bz = lb - by
442 sc(coset(ax, ay, az), coset(0, by, bz)) = &
443 sc(coset(ax, ay + 1, az), coset(0, by - 1, bz)) - &
444 rab(2)*sc(coset(ax, ay, az), coset(0, by - 1, bz))
445 ss(coset(ax, ay, az), coset(0, by, bz)) = &
446 ss(coset(ax, ay + 1, az), coset(0, by - 1, bz)) - &
447 rab(2)*ss(coset(ax, ay, az), coset(0, by - 1, bz))
448 END DO
449
450! *** Shift of angular momentum component x from a to b ***
451
452 DO bx = 1, lb
453 DO by = 0, lb - bx
454 bz = lb - bx - by
455 sc(coset(ax, ay, az), coset(bx, by, bz)) = &
456 sc(coset(ax + 1, ay, az), coset(bx - 1, by, bz)) - &
457 rab(1)*sc(coset(ax, ay, az), coset(bx - 1, by, bz))
458 ss(coset(ax, ay, az), coset(bx, by, bz)) = &
459 ss(coset(ax + 1, ay, az), coset(bx - 1, by, bz)) - &
460 rab(1)*ss(coset(ax, ay, az), coset(bx - 1, by, bz))
461 END DO
462 END DO
463
464 END DO
465 END DO
466 END DO
467
468! *** Vertical recurrence step ***
469
470! *** [a|O|b] = (Pi - Bi)*[a|O|b-1i] + f2*Ni(a)*[a-1i|O|b-1i] + ***
471! *** f2*Ni(b-1i)*[a|O|b-2i] + [a|dO|b-1i] ***
472
473 DO ax = 0, la_max
474 fx = f2*real(ax, dp)
475 DO ay = 0, la_max - ax
476 fy = f2*real(ay, dp)
477 az = la_max - ax - ay
478 fz = f2*real(az, dp)
479
480! *** Increase the angular momentum component z of function b ***
481
482 f3 = f2*real(lb - 1, dp)
483
484 IF (az == 0) THEN
485 sc(coset(ax, ay, az), coset(0, 0, lb)) = &
486 rbp(3)*sc(coset(ax, ay, az), coset(0, 0, lb - 1)) + &
487 f3*sc(coset(ax, ay, az), coset(0, 0, lb - 2)) - &
488 f2*kvec(3)*ss(coset(ax, ay, az), coset(0, 0, lb - 1))
489 ss(coset(ax, ay, az), coset(0, 0, lb)) = &
490 rbp(3)*ss(coset(ax, ay, az), coset(0, 0, lb - 1)) + &
491 f3*ss(coset(ax, ay, az), coset(0, 0, lb - 2)) + &
492 f2*kvec(3)*sc(coset(ax, ay, az), coset(0, 0, lb - 1))
493 ELSE
494 sc(coset(ax, ay, az), coset(0, 0, lb)) = &
495 rbp(3)*sc(coset(ax, ay, az), coset(0, 0, lb - 1)) + &
496 fz*sc(coset(ax, ay, az - 1), coset(0, 0, lb - 1)) + &
497 f3*sc(coset(ax, ay, az), coset(0, 0, lb - 2)) - &
498 f2*kvec(3)*ss(coset(ax, ay, az), coset(0, 0, lb - 1))
499 ss(coset(ax, ay, az), coset(0, 0, lb)) = &
500 rbp(3)*ss(coset(ax, ay, az), coset(0, 0, lb - 1)) + &
501 fz*ss(coset(ax, ay, az - 1), coset(0, 0, lb - 1)) + &
502 f3*ss(coset(ax, ay, az), coset(0, 0, lb - 2)) + &
503 f2*kvec(3)*sc(coset(ax, ay, az), coset(0, 0, lb - 1))
504 END IF
505
506! *** Increase the angular momentum component y of function b ***
507
508 IF (ay == 0) THEN
509 bz = lb - 1
510 sc(coset(ax, ay, az), coset(0, 1, bz)) = &
511 rbp(2)*sc(coset(ax, ay, az), coset(0, 0, bz)) - &
512 f2*kvec(2)*ss(coset(ax, ay, az), coset(0, 0, bz))
513 ss(coset(ax, ay, az), coset(0, 1, bz)) = &
514 rbp(2)*ss(coset(ax, ay, az), coset(0, 0, bz)) + &
515 f2*kvec(2)*sc(coset(ax, ay, az), coset(0, 0, bz))
516 DO by = 2, lb
517 bz = lb - by
518 f3 = f2*real(by - 1, dp)
519 sc(coset(ax, ay, az), coset(0, by, bz)) = &
520 rbp(2)*sc(coset(ax, ay, az), coset(0, by - 1, bz)) + &
521 f3*sc(coset(ax, ay, az), coset(0, by - 2, bz)) - &
522 f2*kvec(2)*ss(coset(ax, ay, az), coset(0, by - 1, bz))
523 ss(coset(ax, ay, az), coset(0, by, bz)) = &
524 rbp(2)*ss(coset(ax, ay, az), coset(0, by - 1, bz)) + &
525 f3*ss(coset(ax, ay, az), coset(0, by - 2, bz)) + &
526 f2*kvec(2)*sc(coset(ax, ay, az), coset(0, by - 1, bz))
527 END DO
528 ELSE
529 bz = lb - 1
530 sc(coset(ax, ay, az), coset(0, 1, bz)) = &
531 rbp(2)*sc(coset(ax, ay, az), coset(0, 0, bz)) + &
532 fy*sc(coset(ax, ay - 1, az), coset(0, 0, bz)) - &
533 f2*kvec(2)*ss(coset(ax, ay, az), coset(0, 0, bz))
534 ss(coset(ax, ay, az), coset(0, 1, bz)) = &
535 rbp(2)*ss(coset(ax, ay, az), coset(0, 0, bz)) + &
536 fy*ss(coset(ax, ay - 1, az), coset(0, 0, bz)) + &
537 f2*kvec(2)*sc(coset(ax, ay, az), coset(0, 0, bz))
538 DO by = 2, lb
539 bz = lb - by
540 f3 = f2*real(by - 1, dp)
541 sc(coset(ax, ay, az), coset(0, by, bz)) = &
542 rbp(2)*sc(coset(ax, ay, az), coset(0, by - 1, bz)) + &
543 fy*sc(coset(ax, ay - 1, az), coset(0, by - 1, bz)) + &
544 f3*sc(coset(ax, ay, az), coset(0, by - 2, bz)) - &
545 f2*kvec(2)*ss(coset(ax, ay, az), coset(0, by - 1, bz))
546 ss(coset(ax, ay, az), coset(0, by, bz)) = &
547 rbp(2)*ss(coset(ax, ay, az), coset(0, by - 1, bz)) + &
548 fy*ss(coset(ax, ay - 1, az), coset(0, by - 1, bz)) + &
549 f3*ss(coset(ax, ay, az), coset(0, by - 2, bz)) + &
550 f2*kvec(2)*sc(coset(ax, ay, az), coset(0, by - 1, bz))
551 END DO
552 END IF
553
554! *** Increase the angular momentum component x of function b ***
555
556 IF (ax == 0) THEN
557 DO by = 0, lb - 1
558 bz = lb - 1 - by
559 sc(coset(ax, ay, az), coset(1, by, bz)) = &
560 rbp(1)*sc(coset(ax, ay, az), coset(0, by, bz)) - &
561 f2*kvec(1)*ss(coset(ax, ay, az), coset(0, by, bz))
562 ss(coset(ax, ay, az), coset(1, by, bz)) = &
563 rbp(1)*ss(coset(ax, ay, az), coset(0, by, bz)) + &
564 f2*kvec(1)*sc(coset(ax, ay, az), coset(0, by, bz))
565 END DO
566 DO bx = 2, lb
567 f3 = f2*real(bx - 1, dp)
568 DO by = 0, lb - bx
569 bz = lb - bx - by
570 sc(coset(ax, ay, az), coset(bx, by, bz)) = &
571 rbp(1)*sc(coset(ax, ay, az), &
572 coset(bx - 1, by, bz)) + &
573 f3*sc(coset(ax, ay, az), coset(bx - 2, by, bz)) - &
574 f2*kvec(1)*ss(coset(ax, ay, az), coset(bx - 1, by, bz))
575 ss(coset(ax, ay, az), coset(bx, by, bz)) = &
576 rbp(1)*ss(coset(ax, ay, az), &
577 coset(bx - 1, by, bz)) + &
578 f3*ss(coset(ax, ay, az), coset(bx - 2, by, bz)) + &
579 f2*kvec(1)*sc(coset(ax, ay, az), coset(bx - 1, by, bz))
580 END DO
581 END DO
582 ELSE
583 DO by = 0, lb - 1
584 bz = lb - 1 - by
585 sc(coset(ax, ay, az), coset(1, by, bz)) = &
586 rbp(1)*sc(coset(ax, ay, az), coset(0, by, bz)) + &
587 fx*sc(coset(ax - 1, ay, az), coset(0, by, bz)) - &
588 f2*kvec(1)*ss(coset(ax, ay, az), coset(0, by, bz))
589 ss(coset(ax, ay, az), coset(1, by, bz)) = &
590 rbp(1)*ss(coset(ax, ay, az), coset(0, by, bz)) + &
591 fx*ss(coset(ax - 1, ay, az), coset(0, by, bz)) + &
592 f2*kvec(1)*sc(coset(ax, ay, az), coset(0, by, bz))
593 END DO
594 DO bx = 2, lb
595 f3 = f2*real(bx - 1, dp)
596 DO by = 0, lb - bx
597 bz = lb - bx - by
598 sc(coset(ax, ay, az), coset(bx, by, bz)) = &
599 rbp(1)*sc(coset(ax, ay, az), &
600 coset(bx - 1, by, bz)) + &
601 fx*sc(coset(ax - 1, ay, az), coset(bx - 1, by, bz)) + &
602 f3*sc(coset(ax, ay, az), coset(bx - 2, by, bz)) - &
603 f2*kvec(1)*ss(coset(ax, ay, az), coset(bx - 1, by, bz))
604 ss(coset(ax, ay, az), coset(bx, by, bz)) = &
605 rbp(1)*ss(coset(ax, ay, az), &
606 coset(bx - 1, by, bz)) + &
607 fx*ss(coset(ax - 1, ay, az), coset(bx - 1, by, bz)) + &
608 f3*ss(coset(ax, ay, az), coset(bx - 2, by, bz)) + &
609 f2*kvec(1)*sc(coset(ax, ay, az), coset(bx - 1, by, bz))
610 END DO
611 END DO
612 END IF
613
614 END DO
615 END DO
616
617 END DO
618
619 END IF
620
621 ELSE
622
623 IF (lb_max > 0) THEN
624
625! *** Vertical recurrence steps: [s|O|s] -> [s|O|b] ***
626
627 rbp(:) = (f1 - 1.0_dp)*rab(:)
628
629! *** [s|O|p] = (Pi - Bi)*[s|O|s] + [s|dO|s] ***
630
631 sc(1, 2) = rbp(1)*sc(1, 1) - f2*kvec(1)*ss(1, 1)
632 sc(1, 3) = rbp(2)*sc(1, 1) - f2*kvec(2)*ss(1, 1)
633 sc(1, 4) = rbp(3)*sc(1, 1) - f2*kvec(3)*ss(1, 1)
634 ss(1, 2) = rbp(1)*ss(1, 1) + f2*kvec(1)*sc(1, 1)
635 ss(1, 3) = rbp(2)*ss(1, 1) + f2*kvec(2)*sc(1, 1)
636 ss(1, 4) = rbp(3)*ss(1, 1) + f2*kvec(3)*sc(1, 1)
637
638! *** [s|O|b] = (Pi - Bi)*[s|O|b-1i] + f2*Ni(b-1i)*[s|O|b-2i] ***
639! *** + [s|dO|b-1i] ***
640
641 DO lb = 2, lb_max
642
643! *** Increase the angular momentum component z of function b ***
644
645 sc(1, coset(0, 0, lb)) = rbp(3)*sc(1, coset(0, 0, lb - 1)) + &
646 f2*real(lb - 1, dp)*sc(1, coset(0, 0, lb - 2)) - &
647 f2*kvec(3)*ss(1, coset(0, 0, lb - 1))
648 ss(1, coset(0, 0, lb)) = rbp(3)*ss(1, coset(0, 0, lb - 1)) + &
649 f2*real(lb - 1, dp)*ss(1, coset(0, 0, lb - 2)) + &
650 f2*kvec(3)*sc(1, coset(0, 0, lb - 1))
651
652! *** Increase the angular momentum component y of function b ***
653
654 bz = lb - 1
655 sc(1, coset(0, 1, bz)) = rbp(2)*sc(1, coset(0, 0, bz)) - &
656 f2*kvec(2)*ss(1, coset(0, 0, bz))
657 ss(1, coset(0, 1, bz)) = rbp(2)*ss(1, coset(0, 0, bz)) + &
658 f2*kvec(2)*sc(1, coset(0, 0, bz))
659
660 DO by = 2, lb
661 bz = lb - by
662 sc(1, coset(0, by, bz)) = rbp(2)*sc(1, coset(0, by - 1, bz)) + &
663 f2*real(by - 1, dp)*sc(1, coset(0, by - 2, bz)) - &
664 f2*kvec(2)*ss(1, coset(0, by - 1, bz))
665 ss(1, coset(0, by, bz)) = rbp(2)*ss(1, coset(0, by - 1, bz)) + &
666 f2*real(by - 1, dp)*ss(1, coset(0, by - 2, bz)) + &
667 f2*kvec(2)*sc(1, coset(0, by - 1, bz))
668 END DO
669
670! *** Increase the angular momentum component x of function b ***
671
672 DO by = 0, lb - 1
673 bz = lb - 1 - by
674 sc(1, coset(1, by, bz)) = rbp(1)*sc(1, coset(0, by, bz)) - &
675 f2*kvec(1)*ss(1, coset(0, by, bz))
676 ss(1, coset(1, by, bz)) = rbp(1)*ss(1, coset(0, by, bz)) + &
677 f2*kvec(1)*sc(1, coset(0, by, bz))
678 END DO
679
680 DO bx = 2, lb
681 f3 = f2*real(bx - 1, dp)
682 DO by = 0, lb - bx
683 bz = lb - bx - by
684 sc(1, coset(bx, by, bz)) = rbp(1)*sc(1, coset(bx - 1, by, bz)) + &
685 f3*sc(1, coset(bx - 2, by, bz)) - &
686 f2*kvec(1)*ss(1, coset(bx - 1, by, bz))
687 ss(1, coset(bx, by, bz)) = rbp(1)*ss(1, coset(bx - 1, by, bz)) + &
688 f3*ss(1, coset(bx - 2, by, bz)) + &
689 f2*kvec(1)*sc(1, coset(bx - 1, by, bz))
690 END DO
691 END DO
692
693 END DO
694
695 END IF
696
697 END IF
698
699 DO j = ncoset(lb_min - 1) + 1, ncoset(lb_max)
700 DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
701 cosab(na + i, nb + j) = sc(i, j)
702 sinab(na + i, nb + j) = ss(i, j)
703 END DO
704 END DO
705
706 IF (PRESENT(dcosab)) THEN
707 la_start = 0
708 lb_start = 0
709 ELSE
710 la_start = la_min
711 lb_start = lb_min
712 END IF
713
714 DO da = 0, da_max - 1
715 ftz = 2.0_dp*zeta(ipgf)
716 DO dax = 0, da
717 DO day = 0, da - dax
718 daz = da - dax - day
719 cda = coset(dax, day, daz) - 1
720 cdax = coset(dax + 1, day, daz) - 1
721 cday = coset(dax, day + 1, daz) - 1
722 cdaz = coset(dax, day, daz + 1) - 1
723 !*** [da/dAi|O|b] = 2*zeta*[a+1i|O|b] - Ni(a)[a-1i|O|b] ***
724
725 DO la = la_start, la_max - da - 1
726 DO ax = 0, la
727 fax = real(ax, dp)
728 DO ay = 0, la - ax
729 fay = real(ay, dp)
730 az = la - ax - ay
731 faz = real(az, dp)
732 coa = coset(ax, ay, az)
733 coamx = coset(ax - 1, ay, az)
734 coamy = coset(ax, ay - 1, az)
735 coamz = coset(ax, ay, az - 1)
736 coapx = coset(ax + 1, ay, az)
737 coapy = coset(ax, ay + 1, az)
738 coapz = coset(ax, ay, az + 1)
739 DO lb = lb_start, lb_max
740 DO bx = 0, lb
741 DO by = 0, lb - bx
742 bz = lb - bx - by
743 cob = coset(bx, by, bz)
744 dscos(coa, cob, cdax) = ftz*sc(coapx, cob) - fax*sc(coamx, cob)
745 dscos(coa, cob, cday) = ftz*sc(coapy, cob) - fay*sc(coamy, cob)
746 dscos(coa, cob, cdaz) = ftz*sc(coapz, cob) - faz*sc(coamz, cob)
747 dssin(coa, cob, cdax) = ftz*ss(coapx, cob) - fax*ss(coamx, cob)
748 dssin(coa, cob, cday) = ftz*ss(coapy, cob) - fay*ss(coamy, cob)
749 dssin(coa, cob, cdaz) = ftz*ss(coapz, cob) - faz*ss(coamz, cob)
750 END DO
751 END DO
752 END DO
753 END DO
754 END DO
755 END DO
756
757 END DO
758 END DO
759 END DO
760
761 IF (PRESENT(dcosab)) THEN
762 DO k = 1, 3
763 DO j = 1, ncoset(lb_max)
764 DO i = 1, ncoset(la_max_set)
765 dcosab(na + i, nb + j, k) = dscos(i, j, k)
766 dsinab(na + i, nb + j, k) = dssin(i, j, k)
767 END DO
768 END DO
769 END DO
770 END IF
771
772 nb = nb + ncoset(lb_max)
773
774 END DO
775
776 na = na + ncoset(la_max_set)
777
778 END DO
779
780 END SUBROUTINE cossin
781
782! **************************************************************************************************
783!> \brief ...
784!> \param la_max ...
785!> \param npgfa ...
786!> \param zeta ...
787!> \param rpgfa ...
788!> \param la_min ...
789!> \param lb_max ...
790!> \param npgfb ...
791!> \param zetb ...
792!> \param rpgfb ...
793!> \param lc_max ...
794!> \param rac ...
795!> \param rbc ...
796!> \param mab ...
797! **************************************************************************************************
798 SUBROUTINE moment(la_max, npgfa, zeta, rpgfa, la_min, &
799 lb_max, npgfb, zetb, rpgfb, &
800 lc_max, rac, rbc, mab)
801
802 INTEGER, INTENT(IN) :: la_max, npgfa
803 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
804 INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
805 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
806 INTEGER, INTENT(IN) :: lc_max
807 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
808 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: mab
809
810 INTEGER :: ax, ay, az, bx, by, bz, i, ipgf, j, &
811 jpgf, k, l, l1, l2, la, la_start, lb, &
812 lx, lx1, ly, ly1, lz, lz1, na, nb, ni
813 REAL(kind=dp) :: dab, f0, f1, f2, f2x, f2y, f2z, f3, fx, &
814 fy, fz, rab2, zetp
815 REAL(kind=dp), DIMENSION(3) :: rab, rap, rbp, rpc
816 REAL(kind=dp), DIMENSION(ncoset(la_max), ncoset(& lb_max), ncoset(lc_max)) :: s
817
818 rab = rbc - rac
819 rab2 = sum(rab**2)
820 dab = sqrt(rab2)
821
822! *** Loop over all pairs of primitive Gaussian-type functions ***
823
824 na = 0
825
826 DO ipgf = 1, npgfa
827
828 nb = 0
829
830 DO jpgf = 1, npgfb
831
832 s = 0.0_dp
833! *** Screening ***
834
835 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
836 DO k = 1, ncoset(lc_max) - 1
837 DO j = nb + 1, nb + ncoset(lb_max)
838 DO i = na + 1, na + ncoset(la_max)
839 mab(i, j, k) = 0.0_dp
840 END DO
841 END DO
842 END DO
843 nb = nb + ncoset(lb_max)
844 cycle
845 END IF
846
847! *** Calculate some prefactors ***
848
849 zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
850
851 f0 = (pi*zetp)**1.5_dp
852 f1 = zetb(jpgf)*zetp
853 f2 = 0.5_dp*zetp
854
855! *** Calculate the basic two-center moment integral [s|M|s] ***
856
857 rpc = zetp*(zeta(ipgf)*rac + zetb(jpgf)*rbc)
858 s(1, 1, 1) = f0*exp(-zeta(ipgf)*f1*rab2)
859 DO l = 2, ncoset(lc_max)
860 lx = indco(1, l)
861 ly = indco(2, l)
862 lz = indco(3, l)
863 l2 = 0
864 IF (lz > 0) THEN
865 l1 = coset(lx, ly, lz - 1)
866 IF (lz > 1) l2 = coset(lx, ly, lz - 2)
867 ni = lz - 1
868 i = 3
869 ELSE IF (ly > 0) THEN
870 l1 = coset(lx, ly - 1, lz)
871 IF (ly > 1) l2 = coset(lx, ly - 2, lz)
872 ni = ly - 1
873 i = 2
874 ELSE IF (lx > 0) THEN
875 l1 = coset(lx - 1, ly, lz)
876 IF (lx > 1) l2 = coset(lx - 2, ly, lz)
877 ni = lx - 1
878 i = 1
879 END IF
880 s(1, 1, l) = rpc(i)*s(1, 1, l1)
881 IF (l2 > 0) s(1, 1, l) = s(1, 1, l) + f2*real(ni, dp)*s(1, 1, l2)
882 END DO
883
884! *** Recurrence steps: [s|M|s] -> [a|M|b] ***
885
886 DO l = 1, ncoset(lc_max)
887
888 lx = indco(1, l)
889 ly = indco(2, l)
890 lz = indco(3, l)
891 IF (lx > 0) THEN
892 lx1 = coset(lx - 1, ly, lz)
893 ELSE
894 lx1 = -1
895 END IF
896 IF (ly > 0) THEN
897 ly1 = coset(lx, ly - 1, lz)
898 ELSE
899 ly1 = -1
900 END IF
901 IF (lz > 0) THEN
902 lz1 = coset(lx, ly, lz - 1)
903 ELSE
904 lz1 = -1
905 END IF
906 f2x = f2*real(lx, dp)
907 f2y = f2*real(ly, dp)
908 f2z = f2*real(lz, dp)
909
910 IF (la_max > 0) THEN
911
912! *** Vertical recurrence steps: [s|M|s] -> [a|M|s] ***
913
914 rap(:) = f1*rab(:)
915
916! *** [p|M|s] = (Pi - Ai)*[s|M|s] + f2*Ni(m-1i)[s|M-1i|s] ***
917
918 s(2, 1, l) = rap(1)*s(1, 1, l)
919 s(3, 1, l) = rap(2)*s(1, 1, l)
920 s(4, 1, l) = rap(3)*s(1, 1, l)
921 IF (lx1 > 0) s(2, 1, l) = s(2, 1, l) + f2x*s(1, 1, lx1)
922 IF (ly1 > 0) s(3, 1, l) = s(3, 1, l) + f2y*s(1, 1, ly1)
923 IF (lz1 > 0) s(4, 1, l) = s(4, 1, l) + f2z*s(1, 1, lz1)
924
925! *** [a|M|s] = (Pi - Ai)*[a-1i|M|s] + f2*Ni(a-1i)*[a-2i|M|s] ***
926! *** + f2*Ni(m-1i)*[a-1i|M-1i|s] ***
927
928 DO la = 2, la_max
929
930! *** Increase the angular momentum component z of function a ***
931
932 s(coset(0, 0, la), 1, l) = rap(3)*s(coset(0, 0, la - 1), 1, l) + &
933 f2*real(la - 1, dp)*s(coset(0, 0, la - 2), 1, l)
934 IF (lz1 > 0) s(coset(0, 0, la), 1, l) = s(coset(0, 0, la), 1, l) + &
935 f2z*s(coset(0, 0, la - 1), 1, lz1)
936
937! *** Increase the angular momentum component y of function a ***
938
939 az = la - 1
940 s(coset(0, 1, az), 1, l) = rap(2)*s(coset(0, 0, az), 1, l)
941 IF (ly1 > 0) s(coset(0, 1, az), 1, l) = s(coset(0, 1, az), 1, l) + &
942 f2y*s(coset(0, 0, az), 1, ly1)
943
944 DO ay = 2, la
945 az = la - ay
946 s(coset(0, ay, az), 1, l) = rap(2)*s(coset(0, ay - 1, az), 1, l) + &
947 f2*real(ay - 1, dp)*s(coset(0, ay - 2, az), 1, l)
948 IF (ly1 > 0) s(coset(0, ay, az), 1, l) = s(coset(0, ay, az), 1, l) + &
949 f2y*s(coset(0, ay - 1, az), 1, ly1)
950 END DO
951
952! *** Increase the angular momentum component x of function a ***
953
954 DO ay = 0, la - 1
955 az = la - 1 - ay
956 s(coset(1, ay, az), 1, l) = rap(1)*s(coset(0, ay, az), 1, l)
957 IF (lx1 > 0) s(coset(1, ay, az), 1, l) = s(coset(1, ay, az), 1, l) + &
958 f2x*s(coset(0, ay, az), 1, lx1)
959 END DO
960
961 DO ax = 2, la
962 f3 = f2*real(ax - 1, dp)
963 DO ay = 0, la - ax
964 az = la - ax - ay
965 s(coset(ax, ay, az), 1, l) = rap(1)*s(coset(ax - 1, ay, az), 1, l) + &
966 f3*s(coset(ax - 2, ay, az), 1, l)
967 IF (lx1 > 0) s(coset(ax, ay, az), 1, l) = s(coset(ax, ay, az), 1, l) + &
968 f2x*s(coset(ax - 1, ay, az), 1, lx1)
969 END DO
970 END DO
971
972 END DO
973
974! *** Recurrence steps: [a|M|s] -> [a|M|b] ***
975
976 IF (lb_max > 0) THEN
977
978 DO j = 2, ncoset(lb_max)
979 DO i = 1, ncoset(la_max)
980 s(i, j, l) = 0.0_dp
981 END DO
982 END DO
983
984! *** Horizontal recurrence steps ***
985
986 rbp(:) = rap(:) - rab(:)
987
988! *** [a|M|p] = [a+1i|M|s] - (Bi - Ai)*[a|M|s] ***
989
990 IF (lb_max == 1) THEN
991 la_start = la_min
992 ELSE
993 la_start = max(0, la_min - 1)
994 END IF
995
996 DO la = la_start, la_max - 1
997 DO ax = 0, la
998 DO ay = 0, la - ax
999 az = la - ax - ay
1000 s(coset(ax, ay, az), 2, l) = s(coset(ax + 1, ay, az), 1, l) - &
1001 rab(1)*s(coset(ax, ay, az), 1, l)
1002 s(coset(ax, ay, az), 3, l) = s(coset(ax, ay + 1, az), 1, l) - &
1003 rab(2)*s(coset(ax, ay, az), 1, l)
1004 s(coset(ax, ay, az), 4, l) = s(coset(ax, ay, az + 1), 1, l) - &
1005 rab(3)*s(coset(ax, ay, az), 1, l)
1006 END DO
1007 END DO
1008 END DO
1009
1010! *** Vertical recurrence step ***
1011
1012! *** [a|M|p] = (Pi - Bi)*[a|M|s] + f2*Ni(a)*[a-1i|M|s] ***
1013! *** + f2*Ni(m)*[a|M-1i|s] ***
1014
1015 DO ax = 0, la_max
1016 fx = f2*real(ax, dp)
1017 DO ay = 0, la_max - ax
1018 fy = f2*real(ay, dp)
1019 az = la_max - ax - ay
1020 fz = f2*real(az, dp)
1021 IF (ax == 0) THEN
1022 s(coset(ax, ay, az), 2, l) = rbp(1)*s(coset(ax, ay, az), 1, l)
1023 ELSE
1024 s(coset(ax, ay, az), 2, l) = rbp(1)*s(coset(ax, ay, az), 1, l) + &
1025 fx*s(coset(ax - 1, ay, az), 1, l)
1026 END IF
1027 IF (lx1 > 0) s(coset(ax, ay, az), 2, l) = s(coset(ax, ay, az), 2, l) + &
1028 f2x*s(coset(ax, ay, az), 1, lx1)
1029 IF (ay == 0) THEN
1030 s(coset(ax, ay, az), 3, l) = rbp(2)*s(coset(ax, ay, az), 1, l)
1031 ELSE
1032 s(coset(ax, ay, az), 3, l) = rbp(2)*s(coset(ax, ay, az), 1, l) + &
1033 fy*s(coset(ax, ay - 1, az), 1, l)
1034 END IF
1035 IF (ly1 > 0) s(coset(ax, ay, az), 3, l) = s(coset(ax, ay, az), 3, l) + &
1036 f2y*s(coset(ax, ay, az), 1, ly1)
1037 IF (az == 0) THEN
1038 s(coset(ax, ay, az), 4, l) = rbp(3)*s(coset(ax, ay, az), 1, l)
1039 ELSE
1040 s(coset(ax, ay, az), 4, l) = rbp(3)*s(coset(ax, ay, az), 1, l) + &
1041 fz*s(coset(ax, ay, az - 1), 1, l)
1042 END IF
1043 IF (lz1 > 0) s(coset(ax, ay, az), 4, l) = s(coset(ax, ay, az), 4, l) + &
1044 f2z*s(coset(ax, ay, az), 1, lz1)
1045 END DO
1046 END DO
1047
1048! *** Recurrence steps: [a|M|p] -> [a|M|b] ***
1049
1050 DO lb = 2, lb_max
1051
1052! *** Horizontal recurrence steps ***
1053
1054! *** [a|M|b] = [a+1i|M|b-1i] - (Bi - Ai)*[a|M|b-1i] ***
1055
1056 IF (lb == lb_max) THEN
1057 la_start = la_min
1058 ELSE
1059 la_start = max(0, la_min - 1)
1060 END IF
1061
1062 DO la = la_start, la_max - 1
1063 DO ax = 0, la
1064 DO ay = 0, la - ax
1065 az = la - ax - ay
1066
1067! *** Shift of angular momentum component z from a to b ***
1068
1069 s(coset(ax, ay, az), coset(0, 0, lb), l) = &
1070 s(coset(ax, ay, az + 1), coset(0, 0, lb - 1), l) - &
1071 rab(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1), l)
1072
1073! *** Shift of angular momentum component y from a to b ***
1074
1075 DO by = 1, lb
1076 bz = lb - by
1077 s(coset(ax, ay, az), coset(0, by, bz), l) = &
1078 s(coset(ax, ay + 1, az), coset(0, by - 1, bz), l) - &
1079 rab(2)*s(coset(ax, ay, az), coset(0, by - 1, bz), l)
1080 END DO
1081
1082! *** Shift of angular momentum component x from a to b ***
1083
1084 DO bx = 1, lb
1085 DO by = 0, lb - bx
1086 bz = lb - bx - by
1087 s(coset(ax, ay, az), coset(bx, by, bz), l) = &
1088 s(coset(ax + 1, ay, az), coset(bx - 1, by, bz), l) - &
1089 rab(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz), l)
1090 END DO
1091 END DO
1092
1093 END DO
1094 END DO
1095 END DO
1096
1097! *** Vertical recurrence step ***
1098
1099! *** [a|M|b] = (Pi - Bi)*[a|M|b-1i] + f2*Ni(a)*[a-1i|M|b-1i] + ***
1100! *** f2*Ni(b-1i)*[a|M|b-2i] + f2*Ni(m)[a|M-1i|b-1i] ***
1101
1102 DO ax = 0, la_max
1103 fx = f2*real(ax, dp)
1104 DO ay = 0, la_max - ax
1105 fy = f2*real(ay, dp)
1106 az = la_max - ax - ay
1107 fz = f2*real(az, dp)
1108
1109! *** Shift of angular momentum component z from a to b ***
1110
1111 f3 = f2*real(lb - 1, dp)
1112
1113 IF (az == 0) THEN
1114 s(coset(ax, ay, az), coset(0, 0, lb), l) = &
1115 rbp(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1), l) + &
1116 f3*s(coset(ax, ay, az), coset(0, 0, lb - 2), l)
1117 ELSE
1118 s(coset(ax, ay, az), coset(0, 0, lb), l) = &
1119 rbp(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1), l) + &
1120 fz*s(coset(ax, ay, az - 1), coset(0, 0, lb - 1), l) + &
1121 f3*s(coset(ax, ay, az), coset(0, 0, lb - 2), l)
1122 END IF
1123 IF (lz1 > 0) s(coset(ax, ay, az), coset(0, 0, lb), l) = &
1124 s(coset(ax, ay, az), coset(0, 0, lb), l) + &
1125 f2z*s(coset(ax, ay, az), coset(0, 0, lb - 1), lz1)
1126
1127! *** Shift of angular momentum component y from a to b ***
1128
1129 IF (ay == 0) THEN
1130 bz = lb - 1
1131 s(coset(ax, ay, az), coset(0, 1, bz), l) = &
1132 rbp(2)*s(coset(ax, ay, az), coset(0, 0, bz), l)
1133 IF (ly1 > 0) s(coset(ax, ay, az), coset(0, 1, bz), l) = &
1134 s(coset(ax, ay, az), coset(0, 1, bz), l) + &
1135 f2y*s(coset(ax, ay, az), coset(0, 0, bz), ly1)
1136 DO by = 2, lb
1137 bz = lb - by
1138 f3 = f2*real(by - 1, dp)
1139 s(coset(ax, ay, az), coset(0, by, bz), l) = &
1140 rbp(2)*s(coset(ax, ay, az), coset(0, by - 1, bz), l) + &
1141 f3*s(coset(ax, ay, az), coset(0, by - 2, bz), l)
1142 IF (ly1 > 0) s(coset(ax, ay, az), coset(0, by, bz), l) = &
1143 s(coset(ax, ay, az), coset(0, by, bz), l) + &
1144 f2y*s(coset(ax, ay, az), coset(0, by - 1, bz), ly1)
1145 END DO
1146 ELSE
1147 bz = lb - 1
1148 s(coset(ax, ay, az), coset(0, 1, bz), l) = &
1149 rbp(2)*s(coset(ax, ay, az), coset(0, 0, bz), l) + &
1150 fy*s(coset(ax, ay - 1, az), coset(0, 0, bz), l)
1151 IF (ly1 > 0) s(coset(ax, ay, az), coset(0, 1, bz), l) = &
1152 s(coset(ax, ay, az), coset(0, 1, bz), l) + &
1153 f2y*s(coset(ax, ay, az), coset(0, 0, bz), ly1)
1154 DO by = 2, lb
1155 bz = lb - by
1156 f3 = f2*real(by - 1, dp)
1157 s(coset(ax, ay, az), coset(0, by, bz), l) = &
1158 rbp(2)*s(coset(ax, ay, az), coset(0, by - 1, bz), l) + &
1159 fy*s(coset(ax, ay - 1, az), coset(0, by - 1, bz), l) + &
1160 f3*s(coset(ax, ay, az), coset(0, by - 2, bz), l)
1161 IF (ly1 > 0) s(coset(ax, ay, az), coset(0, by, bz), l) = &
1162 s(coset(ax, ay, az), coset(0, by, bz), l) + &
1163 f2y*s(coset(ax, ay, az), coset(0, by - 1, bz), ly1)
1164 END DO
1165 END IF
1166
1167! *** Shift of angular momentum component x from a to b ***
1168
1169 IF (ax == 0) THEN
1170 DO by = 0, lb - 1
1171 bz = lb - 1 - by
1172 s(coset(ax, ay, az), coset(1, by, bz), l) = &
1173 rbp(1)*s(coset(ax, ay, az), coset(0, by, bz), l)
1174 IF (lx1 > 0) s(coset(ax, ay, az), coset(1, by, bz), l) = &
1175 s(coset(ax, ay, az), coset(1, by, bz), l) + &
1176 f2x*s(coset(ax, ay, az), coset(0, by, bz), lx1)
1177 END DO
1178 DO bx = 2, lb
1179 f3 = f2*real(bx - 1, dp)
1180 DO by = 0, lb - bx
1181 bz = lb - bx - by
1182 s(coset(ax, ay, az), coset(bx, by, bz), l) = &
1183 rbp(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz), l) + &
1184 f3*s(coset(ax, ay, az), coset(bx - 2, by, bz), l)
1185 IF (lx1 > 0) s(coset(ax, ay, az), coset(bx, by, bz), l) = &
1186 s(coset(ax, ay, az), coset(bx, by, bz), l) + &
1187 f2x*s(coset(ax, ay, az), coset(bx - 1, by, bz), lx1)
1188 END DO
1189 END DO
1190 ELSE
1191 DO by = 0, lb - 1
1192 bz = lb - 1 - by
1193 s(coset(ax, ay, az), coset(1, by, bz), l) = &
1194 rbp(1)*s(coset(ax, ay, az), coset(0, by, bz), l) + &
1195 fx*s(coset(ax - 1, ay, az), coset(0, by, bz), l)
1196 IF (lx1 > 0) s(coset(ax, ay, az), coset(1, by, bz), l) = &
1197 s(coset(ax, ay, az), coset(1, by, bz), l) + &
1198 f2x*s(coset(ax, ay, az), coset(0, by, bz), lx1)
1199 END DO
1200 DO bx = 2, lb
1201 f3 = f2*real(bx - 1, dp)
1202 DO by = 0, lb - bx
1203 bz = lb - bx - by
1204 s(coset(ax, ay, az), coset(bx, by, bz), l) = &
1205 rbp(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz), l) + &
1206 fx*s(coset(ax - 1, ay, az), coset(bx - 1, by, bz), l) + &
1207 f3*s(coset(ax, ay, az), coset(bx - 2, by, bz), l)
1208 IF (lx1 > 0) s(coset(ax, ay, az), coset(bx, by, bz), l) = &
1209 s(coset(ax, ay, az), coset(bx, by, bz), l) + &
1210 f2x*s(coset(ax, ay, az), coset(bx - 1, by, bz), lx1)
1211 END DO
1212 END DO
1213 END IF
1214
1215 END DO
1216 END DO
1217
1218 END DO
1219
1220 END IF
1221
1222 ELSE
1223
1224 IF (lb_max > 0) THEN
1225
1226! *** Vertical recurrence steps: [s|M|s] -> [s|M|b] ***
1227
1228 rbp(:) = (f1 - 1.0_dp)*rab(:)
1229
1230! *** [s|M|p] = (Pi - Bi)*[s|M|s] + f2*Ni(m)*[s|M-1i|s] ***
1231
1232 s(1, 2, l) = rbp(1)*s(1, 1, l)
1233 s(1, 3, l) = rbp(2)*s(1, 1, l)
1234 s(1, 4, l) = rbp(3)*s(1, 1, l)
1235 IF (lx1 > 0) s(1, 2, l) = s(1, 2, l) + f2x*s(1, 1, lx1)
1236 IF (ly1 > 0) s(1, 3, l) = s(1, 3, l) + f2y*s(1, 1, ly1)
1237 IF (lz1 > 0) s(1, 4, l) = s(1, 4, l) + f2z*s(1, 1, lz1)
1238
1239! *** [s|M|b] = (Pi - Bi)*[s|M|b-1i] + f2*Ni(b-1i)*[s|M|b-2i] ***
1240! *** + f2*Ni(m)*[s|M-1i|b-1i] ***
1241
1242 DO lb = 2, lb_max
1243
1244! *** Increase the angular momentum component z of function b ***
1245
1246 s(1, coset(0, 0, lb), l) = rbp(3)*s(1, coset(0, 0, lb - 1), l) + &
1247 f2*real(lb - 1, dp)*s(1, coset(0, 0, lb - 2), l)
1248 IF (lz1 > 0) s(1, coset(0, 0, lb), l) = s(1, coset(0, 0, lb), l) + &
1249 f2z*s(1, coset(0, 0, lb - 1), lz1)
1250
1251! *** Increase the angular momentum component y of function b ***
1252
1253 bz = lb - 1
1254 s(1, coset(0, 1, bz), l) = rbp(2)*s(1, coset(0, 0, bz), l)
1255 IF (ly1 > 0) s(1, coset(0, 1, bz), l) = s(1, coset(0, 1, bz), l) + &
1256 f2y*s(1, coset(0, 0, bz), ly1)
1257
1258 DO by = 2, lb
1259 bz = lb - by
1260 s(1, coset(0, by, bz), l) = rbp(2)*s(1, coset(0, by - 1, bz), l) + &
1261 f2*real(by - 1, dp)*s(1, coset(0, by - 2, bz), l)
1262 IF (ly1 > 0) s(1, coset(0, by, bz), l) = s(1, coset(0, by, bz), l) + &
1263 f2y*s(1, coset(0, by - 1, bz), ly1)
1264 END DO
1265
1266! *** Increase the angular momentum component x of function b ***
1267
1268 DO by = 0, lb - 1
1269 bz = lb - 1 - by
1270 s(1, coset(1, by, bz), l) = rbp(1)*s(1, coset(0, by, bz), l)
1271 IF (lx1 > 0) s(1, coset(1, by, bz), l) = s(1, coset(1, by, bz), l) + &
1272 f2x*s(1, coset(0, by, bz), lx1)
1273 END DO
1274
1275 DO bx = 2, lb
1276 f3 = f2*real(bx - 1, dp)
1277 DO by = 0, lb - bx
1278 bz = lb - bx - by
1279 s(1, coset(bx, by, bz), l) = rbp(1)*s(1, coset(bx - 1, by, bz), l) + &
1280 f3*s(1, coset(bx - 2, by, bz), l)
1281 IF (lx1 > 0) s(1, coset(bx, by, bz), l) = s(1, coset(bx, by, bz), l) + &
1282 f2x*s(1, coset(bx - 1, by, bz), lx1)
1283 END DO
1284 END DO
1285
1286 END DO
1287
1288 END IF
1289
1290 END IF
1291
1292 END DO
1293
1294 DO k = 2, ncoset(lc_max)
1295 DO j = 1, ncoset(lb_max)
1296 DO i = 1, ncoset(la_max)
1297 mab(na + i, nb + j, k - 1) = s(i, j, k)
1298 END DO
1299 END DO
1300 END DO
1301
1302 nb = nb + ncoset(lb_max)
1303
1304 END DO
1305
1306 na = na + ncoset(la_max)
1307
1308 END DO
1309
1310 END SUBROUTINE moment
1311
1312! **************************************************************************************************
1313!> \brief This returns the derivative of the moment integrals [a|\mu|b].
1314!> By default, it differentiates the primitive on the right:
1315!> [a|\mu|d/dR_bi] = 2*zetb*[a|\mu|b+1i] - Ni(b)[a|\mu|b-1i]
1316!> A weighted derivative combines the left and right primitive derivatives
1317!> using deltaR for the corresponding atom centers.
1318!> order indicates the max order of the moment operator to be calculated
1319!> 1: dipole
1320!> 2: quadrupole
1321!> ...
1322!> \param la_max ...
1323!> \param npgfa ...
1324!> \param zeta ...
1325!> \param rpgfa ...
1326!> \param la_min ...
1327!> \param lb_max ...
1328!> \param npgfb ...
1329!> \param zetb ...
1330!> \param rpgfb ...
1331!> \param lb_min ...
1332!> \param order ...
1333!> \param rac ...
1334!> \param rbc ...
1335!> \param difmab ...
1336!> \param mab_ext ...
1337!> \param deltaR optional weights for the left and right primitive derivatives
1338!> \param lambda optional atom selector for the factor (iatom == lambda) - (jatom == lambda)
1339!> \param iatom atom associated with the left basis function
1340!> \param jatom atom associated with the right basis function
1341!> \note
1342! **************************************************************************************************
1343 SUBROUTINE diff_momop(la_max, npgfa, zeta, rpgfa, la_min, &
1344 lb_max, npgfb, zetb, rpgfb, lb_min, &
1345 order, rac, rbc, difmab, mab_ext, deltaR, lambda, iatom, jatom)
1346
1347 INTEGER, INTENT(IN) :: la_max, npgfa
1348 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
1349 INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
1350 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
1351 INTEGER, INTENT(IN) :: lb_min, order
1352 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
1353 REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(OUT) :: difmab
1354 REAL(kind=dp), DIMENSION(:, :, :), OPTIONAL, &
1355 POINTER :: mab_ext
1356 REAL(kind=dp), DIMENSION(:, :), OPTIONAL, POINTER :: deltar
1357 INTEGER, INTENT(IN), OPTIONAL :: lambda, iatom, jatom
1358
1359 INTEGER :: ider, imom, lda, lda_min, ldb, ldb_min
1360 REAL(kind=dp) :: dab, rab(3), rab2
1361 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: difmab_tmp
1362 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: mab
1363
1364 rab = rbc - rac
1365 rab2 = sum(rab**2)
1366 dab = sqrt(rab2)
1367
1368 lda_min = max(0, la_min - 1)
1369 ldb_min = max(0, lb_min - 1)
1370 lda = ncoset(la_max)*npgfa
1371 ldb = ncoset(lb_max)*npgfb
1372 ALLOCATE (difmab_tmp(lda, ldb, 3))
1373
1374 IF (PRESENT(mab_ext)) THEN
1375 mab => mab_ext
1376 ELSE
1377 ALLOCATE (mab(npgfa*ncoset(la_max + 1), npgfb*ncoset(lb_max + 1), &
1378 ncoset(order) - 1))
1379 mab = 0.0_dp
1380! *** Calculate the primitive overlap integrals ***
1381 CALL moment(la_max + 1, npgfa, zeta, rpgfa, lda_min, &
1382 lb_max + 1, npgfb, zetb, rpgfb, &
1383 order, rac, rbc, mab)
1384
1385 END IF
1386 DO imom = 1, ncoset(order) - 1
1387 CALL adbdr(la_max, npgfa, rpgfa, la_min, &
1388 lb_max, npgfb, zetb, rpgfb, lb_min, &
1389 dab, mab(:, :, imom), difmab_tmp(:, :, 1), &
1390 difmab_tmp(:, :, 2), difmab_tmp(:, :, 3))
1391 IF (PRESENT(deltar)) THEN
1392 cpassert(ASSOCIATED(deltar))
1393 cpassert(PRESENT(iatom) .AND. PRESENT(jatom))
1394 DO ider = 1, 3
1395 difmab(1:lda, 1:ldb, imom, ider) = &
1396 difmab_tmp(1:lda, 1:ldb, ider)*deltar(ider, jatom)
1397 END DO
1398
1399 CALL dabdr(la_max, npgfa, zeta, rpgfa, la_min, &
1400 lb_max, npgfb, rpgfb, lb_min, &
1401 dab, mab(:, :, imom), difmab_tmp(:, :, 1), &
1402 difmab_tmp(:, :, 2), difmab_tmp(:, :, 3))
1403 DO ider = 1, 3
1404 difmab(1:lda, 1:ldb, imom, ider) = difmab(1:lda, 1:ldb, imom, ider) &
1405 + difmab_tmp(1:lda, 1:ldb, ider)*deltar(ider, iatom)
1406 END DO
1407 ELSE
1408 difmab(1:lda, 1:ldb, imom, :) = difmab_tmp(1:lda, 1:ldb, :)
1409 END IF
1410 END DO
1411
1412 IF (PRESENT(lambda)) THEN
1413 cpassert(.NOT. PRESENT(deltar))
1414 cpassert(PRESENT(iatom) .AND. PRESENT(jatom))
1415 IF (iatom == lambda .AND. jatom == lambda) THEN
1416 difmab = 0.0_dp
1417 ELSE IF (iatom == lambda) THEN
1418 ! The right-hand derivative is selected with a positive sign.
1419 ELSE IF (jatom == lambda) THEN
1420 difmab = -difmab
1421 ELSE
1422 difmab = 0.0_dp
1423 END IF
1424 END IF
1425
1426 IF (PRESENT(mab_ext)) THEN
1427 NULLIFY (mab)
1428 ELSE
1429 DEALLOCATE (mab)
1430 END IF
1431 DEALLOCATE (difmab_tmp)
1432
1433 END SUBROUTINE diff_momop
1434
1435! **************************************************************************************************
1436!> \brief This returns the derivative of the dipole integrals [a|x|b], with respect
1437!> to the position of the primitive on the left and right, i.e.
1438!> [da/dR_ai|\mu|b] = 2*zeta*[a+1i|\mu|b] - Ni(a)[a-1i|\mu|b]
1439!> \param la_max ...
1440!> \param npgfa ...
1441!> \param zeta ...
1442!> \param rpgfa ...
1443!> \param la_min ...
1444!> \param lb_max ...
1445!> \param npgfb ...
1446!> \param zetb ...
1447!> \param rpgfb ...
1448!> \param lb_min ...
1449!> \param order ...
1450!> \param rac ...
1451!> \param rbc ...
1452!> \param pab ...
1453!> \param forcea ...
1454!> \param forceb ...
1455!> \note
1456! **************************************************************************************************
1457 SUBROUTINE dipole_force(la_max, npgfa, zeta, rpgfa, la_min, &
1458 lb_max, npgfb, zetb, rpgfb, lb_min, &
1459 order, rac, rbc, pab, forcea, forceb)
1460
1461 INTEGER, INTENT(IN) :: la_max, npgfa
1462 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
1463 INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
1464 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
1465 INTEGER, INTENT(IN) :: lb_min, order
1466 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
1467 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: pab
1468 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: forcea, forceb
1469
1470 INTEGER :: i, imom, ipgf, j, jpgf, lda, lda_min, &
1471 ldb, ldb_min, na, nb
1472 REAL(kind=dp) :: dab, rab(3), rab2
1473 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: difmab, mab
1474
1475 cpassert(order == 1)
1476 mark_used(order)
1477
1478 rab = rbc - rac
1479 rab2 = sum(rab**2)
1480 dab = sqrt(rab2)
1481
1482 lda_min = max(0, la_min - 1)
1483 ldb_min = max(0, lb_min - 1)
1484 lda = ncoset(la_max)*npgfa
1485 ldb = ncoset(lb_max)*npgfb
1486 ALLOCATE (difmab(lda, ldb, 3))
1487 ALLOCATE (mab(npgfa*ncoset(la_max + 1), npgfb*ncoset(lb_max + 1), 3))
1488 mab = 0.0_dp
1489 CALL moment(la_max + 1, npgfa, zeta, rpgfa, lda_min, &
1490 lb_max + 1, npgfb, zetb, rpgfb, 1, rac, rbc, mab)
1491
1492 DO imom = 1, 3
1493 difmab = 0.0_dp
1494 CALL adbdr(la_max, npgfa, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, &
1495 dab, mab(:, :, imom), difmab(:, :, 1), difmab(:, :, 2), difmab(:, :, 3))
1496 na = 0
1497 DO ipgf = 1, npgfa
1498 nb = 0
1499 DO jpgf = 1, npgfb
1500 DO j = nb + ncoset(lb_min - 1) + 1, nb + ncoset(lb_max)
1501 DO i = na + ncoset(la_min - 1) + 1, na + ncoset(la_max)
1502 forceb(imom, 1) = forceb(imom, 1) + pab(i, j)*difmab(i, j, 1)
1503 forceb(imom, 2) = forceb(imom, 2) + pab(i, j)*difmab(i, j, 2)
1504 forceb(imom, 3) = forceb(imom, 3) + pab(i, j)*difmab(i, j, 3)
1505 END DO
1506 END DO
1507 nb = nb + ncoset(lb_max)
1508 END DO
1509 na = na + ncoset(la_max)
1510 END DO
1511
1512 difmab = 0.0_dp
1513 CALL dabdr(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, rpgfb, lb_min, &
1514 dab, mab(:, :, imom), difmab(:, :, 1), difmab(:, :, 2), difmab(:, :, 3))
1515 na = 0
1516 DO ipgf = 1, npgfa
1517 nb = 0
1518 DO jpgf = 1, npgfb
1519 DO j = nb + ncoset(lb_min - 1) + 1, nb + ncoset(lb_max)
1520 DO i = na + ncoset(la_min - 1) + 1, na + ncoset(la_max)
1521 forcea(imom, 1) = forcea(imom, 1) + pab(i, j)*difmab(i, j, 1)
1522 forcea(imom, 2) = forcea(imom, 2) + pab(i, j)*difmab(i, j, 2)
1523 forcea(imom, 3) = forcea(imom, 3) + pab(i, j)*difmab(i, j, 3)
1524 END DO
1525 END DO
1526 nb = nb + ncoset(lb_max)
1527 END DO
1528 na = na + ncoset(la_max)
1529 END DO
1530 END DO
1531
1532 DEALLOCATE (mab, difmab)
1533
1534 END SUBROUTINE dipole_force
1535
1536END MODULE ai_moments
1537
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:225
ai_derivatives
Calculate the first derivative of an integral block.
Definition ai_derivatives.F:34
ai_derivatives::adbdr
subroutine, public adbdr(la_max, npgfa, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, dab, ab, adbdx, adbdy, adbdz)
Calculate the first derivative of an integral block. This takes the derivative with respect to the at...
Definition ai_derivatives.F:231
ai_derivatives::dabdr
subroutine, public dabdr(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, rpgfb, lb_min, dab, ab, dabdx, dabdy, dabdz)
Calculate the first derivative of an integral block. This takes the derivative with respect to the at...
Definition ai_derivatives.F:82
ai_moments
Calculation of the moment integrals over Cartesian Gaussian-type functions.
Definition ai_moments.F:17
ai_moments::dipole_force
subroutine, public dipole_force(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, order, rac, rbc, pab, forcea, forceb)
This returns the derivative of the dipole integrals [a|x|b], with respect to the position of the prim...
Definition ai_moments.F:1462
ai_moments::diff_momop
subroutine, public diff_momop(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, order, rac, rbc, difmab, mab_ext, deltar, lambda, iatom, jatom)
This returns the derivative of the moment integrals [a|\mu|b]. By default, it differentiates the prim...
Definition ai_moments.F:1348
ai_moments::contract_cossin
subroutine, public contract_cossin(cos_block, sin_block, iatom, ncoa, nsgfa, sgfa, sphi_a, ldsa, jatom, ncob, nsgfb, sgfb, sphi_b, ldsb, cosab, sinab, ldab, work, ldwork)
...
Definition ai_moments.F:77
ai_moments::cossin
subroutine, public cossin(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max, npgfb, zetb, rpgfb, lb_min, rac, rbc, kvec, cosab, sinab, dcosab, dsinab)
...
Definition ai_moments.F:155
ai_moments::moment
subroutine, public moment(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lc_max, rac, rbc, mab)
...
Definition ai_moments.F:802
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
orbital_pointers::coset
integer, dimension(:, :, :), allocatable, public coset
Definition orbital_pointers.F:45
orbital_pointers::indco
integer, dimension(:, :), allocatable, public indco
Definition orbital_pointers.F:43