40 #include "../base/base_uses.f90"
46 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'ai_moments'
77 lb_max, npgfb, zetb, rpgfb, lb_min, &
78 order, rac, rbc, difmab, lambda, iatom, jatom)
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
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
95 dab = sqrt(sum(rab**2))
97 lda_min = max(0, la_min - 1)
98 ldb_min = max(0, lb_min - 1)
100 ldb =
ncoset(lb_max)*npgfb
101 ALLOCATE (difmab_tmp(lda, ldb, 3))
107 CALL moment(la_max + 1, npgfa, zeta, rpgfa, lda_min, &
108 lb_max + 1, npgfb, zetb, rpgfb, &
109 order, rac, rbc, mab)
112 DO beta = 1,
ncoset(order) - 1
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))
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)
128 DEALLOCATE (difmab_tmp)
162 lb_max, npgfb, zetb, rpgfb, lb_min, &
163 order, rac, rbc, difmab, mab_ext, deltaR, iatom, jatom)
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, &
174 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN), &
175 OPTIONAL,
POINTER :: deltar
176 INTEGER,
INTENT(IN),
OPTIONAL :: iatom, jatom
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
184 dab = sqrt(sum(rab**2))
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))
192 IF (
PRESENT(mab_ext))
THEN
195 ALLOCATE (mab(npgfa*
ncoset(la_max + 1), npgfb*
ncoset(lb_max + 1), &
199 CALL moment(la_max + 1, npgfa, zeta, rpgfa, lda_min, &
200 lb_max + 1, npgfb, zetb, rpgfb, &
201 order, rac, rbc, mab)
203 DO imom = 1,
ncoset(order) - 1
204 difmab(:, :, imom, :) = 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))
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)
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))
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)
227 IF (
PRESENT(mab_ext))
THEN
232 DEALLOCATE (difmab_tmp)
258 iatom, ncoa, nsgfa, sgfa, sphi_a, ldsa, &
259 jatom, ncob, nsgfb, sgfb, sphi_b, ldsb, &
260 cosab, sinab, ldab, work, ldwork)
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
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)
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), &
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), &
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)
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), &
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), &
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)
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
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
358 DIMENSION(ncoset(la_max_set+1), ncoset(lb_max)) :: sc, ss
363 k2 = kvec(1)*kvec(1) + kvec(2)*kvec(2) + kvec(3)*kvec(3)
365 IF (
PRESENT(dcosab))
THEN
367 la_max = la_max_set + 1
368 la_min = max(0, la_min_set - 1)
378 IF (
PRESENT(dcosab))
THEN
379 na =
ncoset(la_max - 1)*npgfa
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
403 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab)
THEN
410 zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
412 f0 = (
pi*zetp)**1.5_dp
416 kdp = zetp*dot_product(kvec, zeta(ipgf)*rac + zetb(jpgf)*rbc)
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)
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)
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)
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)
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)
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)
484 f3 = f2*real(ax - 1,
dp)
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)
511 rbp(:) = rap(:) - rab(:)
515 IF (lb_max == 1)
THEN
518 la_start = max(0, la_min - 1)
521 DO la = la_start, la_max - 1
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)
548 DO ay = 0, la_max - ax
550 az = la_max - ax - ay
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)
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)
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)
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)
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)
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)
602 IF (lb == lb_max)
THEN
605 la_start = max(0, la_min - 1)
608 DO la = la_start, la_max - 1
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))
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))
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))
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))
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))
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))
659 DO ay = 0, la_max - ax
661 az = la_max - ax - ay
666 f3 = f2*real(lb - 1,
dp)
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))
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))
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))
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))
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))
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))
702 f3 = f2*real(by - 1,
dp)
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))
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))
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))
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))
724 f3 = f2*real(by - 1,
dp)
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))
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))
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))
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))
751 f3 = f2*real(bx - 1,
dp)
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))
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))
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))
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))
779 f3 = f2*real(bx - 1,
dp)
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))
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))
811 rbp(:) = (f1 - 1.0_dp)*rab(:)
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)
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))
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))
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))
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))
865 f3 = f2*real(bx - 1,
dp)
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))
885 cosab(na + i, nb + j) = sc(i, j)
886 sinab(na + i, nb + j) = ss(i, j)
890 IF (
PRESENT(dcosab))
THEN
898 DO da = 0, da_max - 1
899 ftz = 2.0_dp*zeta(ipgf)
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
909 DO la = la_start, la_max - da - 1
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
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)
945 IF (
PRESENT(dcosab))
THEN
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)
960 na = na +
ncoset(la_max_set)
982 SUBROUTINE moment(la_max, npgfa, zeta, rpgfa, la_min, &
983 lb_max, npgfb, zetb, rpgfb, &
984 lc_max, rac, rbc, mab)
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
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, &
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
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
1033 zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
1035 f0 = (
pi*zetp)**1.5_dp
1036 f1 = zetb(jpgf)*zetp
1041 rpc = zetp*(zeta(ipgf)*rac + zetb(jpgf)*rbc)
1042 s(1, 1, 1) = f0*exp(-zeta(ipgf)*f1*rab2)
1049 l1 =
coset(lx, ly, lz - 1)
1050 IF (lz > 1) l2 =
coset(lx, ly, lz - 2)
1053 ELSE IF (ly > 0)
THEN
1054 l1 =
coset(lx, ly - 1, lz)
1055 IF (ly > 1) l2 =
coset(lx, ly - 2, lz)
1058 ELSE IF (lx > 0)
THEN
1059 l1 =
coset(lx - 1, ly, lz)
1060 IF (lx > 1) l2 =
coset(lx - 2, ly, lz)
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)
1076 lx1 =
coset(lx - 1, ly, lz)
1081 ly1 =
coset(lx, ly - 1, lz)
1086 lz1 =
coset(lx, ly, lz - 1)
1090 f2x = f2*real(lx,
dp)
1091 f2y = f2*real(ly,
dp)
1092 f2z = f2*real(lz,
dp)
1094 IF (la_max > 0)
THEN
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)
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)
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)
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)
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)
1146 f3 = f2*real(ax - 1,
dp)
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)
1160 IF (lb_max > 0)
THEN
1170 rbp(:) = rap(:) - rab(:)
1174 IF (lb_max == 1)
THEN
1177 la_start = max(0, la_min - 1)
1180 DO la = la_start, la_max - 1
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)
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)
1206 s(
coset(ax, ay, az), 2, l) = rbp(1)*s(
coset(ax, ay, az), 1, l)
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)
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)
1214 s(
coset(ax, ay, az), 3, l) = rbp(2)*s(
coset(ax, ay, az), 1, l)
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)
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)
1222 s(
coset(ax, ay, az), 4, l) = rbp(3)*s(
coset(ax, ay, az), 1, l)
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)
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)
1240 IF (lb == lb_max)
THEN
1243 la_start = max(0, la_min - 1)
1246 DO la = la_start, la_max - 1
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)
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)
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)
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)
1295 f3 = f2*real(lb - 1,
dp)
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)
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)
1307 IF (lz1 > 0) s(
coset(ax, ay, az),
coset(0, 0, lb), l) = &
1309 f2z*s(
coset(ax, ay, az),
coset(0, 0, lb - 1), lz1)
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) = &
1319 f2y*s(
coset(ax, ay, az),
coset(0, 0, bz), ly1)
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)
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) = &
1337 f2y*s(
coset(ax, ay, az),
coset(0, 0, bz), ly1)
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)
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)
1363 f3 = f2*real(bx - 1,
dp)
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)
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)
1385 f3 = f2*real(bx - 1,
dp)
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)
1408 IF (lb_max > 0)
THEN
1412 rbp(:) = (f1 - 1.0_dp)*rab(:)
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)
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)
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)
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)
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)
1460 f3 = f2*real(bx - 1,
dp)
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)
1481 mab(na + i, nb + j, k - 1) = s(i, j, k)
1513 SUBROUTINE diffop(la_max, npgfa, zeta, rpgfa, la_min, &
1514 lb_max, npgfb, zetb, rpgfb, lb_min, &
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
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
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)
1537 ALLOCATE (s(lds, lds,
ncoset(1)))
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)
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))
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)
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, &
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
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))
1607 IF (
PRESENT(mab_ext))
THEN
1610 ALLOCATE (mab(npgfa*
ncoset(la_max + 1), npgfb*
ncoset(lb_max + 1), &
1614 CALL moment(la_max + 1, npgfa, zeta, rpgfa, lda_min, &
1615 lb_max + 1, npgfb, zetb, rpgfb, &
1616 order, rac, rbc, mab)
1619 DO imom = 1,
ncoset(order) - 1
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)
1630 IF (
PRESENT(mab_ext))
THEN
1635 DEALLOCATE (difmab_tmp)
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)
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
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
1679 cpassert(order == 1)
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))
1694 CALL moment(la_max + 1, npgfa, zeta, rpgfa, lda_min, &
1695 lb_max + 1, npgfb, zetb, rpgfb, 1, rac, rbc, mab)
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))
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)
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))
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)
1737 DEALLOCATE (mab, difmab)
1742 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.
Calculate the first derivative of an integral block.
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...
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...
Calculation of the moment integrals over Cartesian Gaussian-type functions.
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...
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...
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...
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)
...
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...
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)
...
subroutine, public moment(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lc_max, rac, rbc, mab)
...
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...
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
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...
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public ncoset
integer, dimension(:, :, :), allocatable, public coset
integer, dimension(:, :), allocatable, public indco