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