37 #include "../base/base_uses.f90"
43 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'ai_overlap_ppl'
99 SUBROUTINE ppl_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
100 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
101 rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
102 hab2, hab2_work, deltaR, iatom, jatom, katom)
103 INTEGER,
INTENT(IN) :: la_max_set, la_min_set, npgfa
104 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: rpgfa, zeta
105 INTEGER,
INTENT(IN) :: lb_max_set, lb_min_set, npgfb
106 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: rpgfb, zetb
107 INTEGER,
INTENT(IN) :: nexp_ppl
108 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: alpha_ppl
109 INTEGER,
DIMENSION(:),
INTENT(IN) :: nct_ppl
110 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN) :: cexp_ppl
111 REAL(kind=
dp),
INTENT(IN) :: rpgfc
112 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: rab
113 REAL(kind=
dp),
INTENT(IN) :: dab
114 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: rac
115 REAL(kind=
dp),
INTENT(IN) :: dac
116 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: rbc
117 REAL(kind=
dp),
INTENT(IN) :: dbc
118 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: vab
119 REAL(kind=
dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: s
120 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN), &
122 REAL(kind=
dp),
DIMENSION(3),
INTENT(OUT),
OPTIONAL :: force_a, force_b
123 REAL(kind=
dp),
DIMENSION(:, :, :),
INTENT(INOUT), &
124 OPTIONAL :: fs, hab2, hab2_work
125 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN), &
127 INTEGER,
INTENT(IN),
OPTIONAL :: iatom, jatom, katom
129 INTEGER :: iexp, ij, ipgf, jpgf, mmax, nexp
130 REAL(kind=
dp) :: rho, sab, t, zetc
131 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: auxint
132 REAL(kind=
dp),
DIMENSION(3) :: pci
134 IF (
PRESENT(pab))
THEN
135 cpassert(
PRESENT(force_a))
136 cpassert(
PRESENT(force_b))
137 cpassert(
PRESENT(fs))
138 mmax = la_max_set + lb_max_set + 2
141 ELSE IF (
PRESENT(hab2))
THEN
142 mmax = la_max_set + lb_max_set + 2
144 mmax = la_max_set + lb_max_set
147 ALLOCATE (auxint(0:mmax, npgfa*npgfb))
154 IF (rpgfa(ipgf) + rpgfc < dac) cycle
157 IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
158 (rpgfa(ipgf) + rpgfb(jpgf) < dab)) cycle
159 ij = (ipgf - 1)*npgfb + jpgf
160 rho = zeta(ipgf) + zetb(jpgf)
161 pci(:) = -(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))/rho
162 sab = exp(-(zeta(ipgf)*zetb(jpgf)/rho*dab*dab))
163 t = rho*sum(pci(:)*pci(:))
165 DO iexp = 1, nexp_ppl
167 zetc = alpha_ppl(iexp)
168 CALL ppl_aux(auxint(0:mmax, ij), mmax, t, rho, nexp, cexp_ppl(:, iexp), zetc)
171 auxint(0:mmax, ij) = sab*auxint(0:mmax, ij)
176 CALL os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
177 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, auxint, rpgfc, &
178 rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
179 vab2=hab2, vab2_work=hab2_work, &
180 deltar=deltar, iatom=iatom, jatom=jatom, katom=katom)
214 nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
216 INTEGER,
INTENT(IN) :: la_max_set, la_min_set, npgfa
217 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: rpgfa, zeta
218 INTEGER,
INTENT(IN) :: nexp_ppl
219 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: alpha_ppl
220 INTEGER,
DIMENSION(:),
INTENT(IN) :: nct_ppl
221 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN) :: cexp_ppl
222 REAL(kind=
dp),
INTENT(IN) :: rpgfc
223 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: rac
224 REAL(kind=
dp),
INTENT(IN) :: dac
225 REAL(kind=
dp),
DIMENSION(:),
INTENT(INOUT) :: va
226 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT), &
229 INTEGER :: i, iexp, ipgf, iw, mmax, na, nexp
230 INTEGER,
DIMENSION(3) :: ani, anm, anp
232 REAL(kind=
dp) :: oint, oref, rho, t, zetc
233 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: auxint
234 REAL(kind=
dp),
DIMENSION(3) :: doint, doref
238 IF (
PRESENT(dva))
THEN
239 mmax = la_max_set + 1
244 ALLOCATE (auxint(0:mmax, npgfa))
249 IF (rpgfa(ipgf) + rpgfc < dac) cycle
253 DO iexp = 1, nexp_ppl
255 zetc = alpha_ppl(iexp)
256 CALL ppl_aux(auxint(0:mmax, ipgf), mmax, t, rho, nexp, cexp_ppl(:, iexp), zetc)
261 IF (
PRESENT(dva))
THEN
262 CALL os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
263 auxint, rpgfc, rac, dac, va, dva)
265 CALL os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
266 auxint, rpgfc, rac, dac, va)
275 IF (rpgfa(ipgf) + rpgfc < dac)
THEN
276 na = na +
ncoset(la_max_set)
282 ani(1:3) =
indco(1:3, i)
283 oint = ppl_ri_test(rho, ani, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
285 IF (abs(oint - oref) > 1.0e-12_dp)
THEN
286 WRITE (iw,
'(A,3i2,i5,F10.4,2G24.12)')
"PPL int error ", ani, la_max_set, dac, oint, oref
288 IF (
PRESENT(dva))
THEN
289 anp = ani + (/1, 0, 0/)
290 anm = ani - (/1, 0, 0/)
291 doint(1) = 2._dp*rho*ppl_ri_test(rho, anp, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl) &
292 - ani(1)*ppl_ri_test(rho, anm, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
293 anp = ani + (/0, 1, 0/)
294 anm = ani - (/0, 1, 0/)
295 doint(2) = 2._dp*rho*ppl_ri_test(rho, anp, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl) &
296 - ani(2)*ppl_ri_test(rho, anm, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
297 anp = ani + (/0, 0, 1/)
298 anm = ani - (/0, 0, 1/)
299 doint(3) = 2._dp*rho*ppl_ri_test(rho, anp, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl) &
300 - ani(3)*ppl_ri_test(rho, anm, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
301 doref(1:3) = dva(na + i, 1:3)
302 IF (any(abs(doint - doref) > 1.0e-6_dp))
THEN
303 WRITE (iw,
'(A,3i2,i5,F10.4,2G24.12)')
" PPL dint error ", &
304 ani, la_max_set, dac, sum(abs(doint)), sum(abs(doref))
308 na = na +
ncoset(la_max_set)
325 FUNCTION ppl_ri_test(rho, ani, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
RESULT(oint)
326 REAL(kind=
dp),
INTENT(IN) :: rho
327 INTEGER,
DIMENSION(3),
INTENT(IN) :: ani
328 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: rac
329 INTEGER,
INTENT(IN) :: nexp_ppl
330 INTEGER,
DIMENSION(:),
INTENT(IN) :: nct_ppl
331 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: alpha_ppl
332 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN) :: cexp_ppl
333 REAL(kind=
dp) :: oint
335 INTEGER :: iexp, nexp, ni
336 REAL(kind=
dp) :: cn, zetc
337 REAL(kind=
dp),
DIMENSION(3) :: ra
341 DO iexp = 1, nexp_ppl
343 zetc = alpha_ppl(iexp)
346 cn = cexp_ppl(ni, iexp)
358 oint = oint + 2.0_dp*cn*
os_overlap2(ani, (/2, 2, 0/))
359 oint = oint + 2.0_dp*cn*
os_overlap2(ani, (/0, 2, 2/))
360 oint = oint + 2.0_dp*cn*
os_overlap2(ani, (/2, 0, 2/))
365 oint = oint + 3.0_dp*cn*
os_overlap2(ani, (/4, 2, 0/))
366 oint = oint + 3.0_dp*cn*
os_overlap2(ani, (/4, 0, 2/))
367 oint = oint + 3.0_dp*cn*
os_overlap2(ani, (/2, 4, 0/))
368 oint = oint + 3.0_dp*cn*
os_overlap2(ani, (/0, 4, 2/))
369 oint = oint + 3.0_dp*cn*
os_overlap2(ani, (/2, 0, 4/))
370 oint = oint + 3.0_dp*cn*
os_overlap2(ani, (/0, 2, 4/))
371 oint = oint + 6.0_dp*cn*
os_overlap2(ani, (/2, 2, 2/))
373 cpabort(
"OVERLAP_PPL")
378 END FUNCTION ppl_ri_test
390 SUBROUTINE ppl_aux(auxint, mmax, t, rho, nexp_ppl, cexp_ppl, zetc)
391 INTEGER,
INTENT(IN) :: mmax
392 REAL(kind=
dp),
DIMENSION(0:mmax) :: auxint
393 REAL(kind=
dp),
INTENT(IN) :: t, rho
394 INTEGER,
INTENT(IN) :: nexp_ppl
395 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: cexp_ppl
396 REAL(kind=
dp),
INTENT(IN) :: zetc
398 INTEGER :: i, j, ke, kp, pmax
399 REAL(kind=
dp) :: a2, a3, a4, cc, f, q, q2, q4, q6, rho2, &
401 REAL(kind=
dp),
DIMENSION(0:6) :: polder
402 REAL(kind=
dp),
DIMENSION(0:mmax) :: expder
404 cpassert(nexp_ppl > 0)
408 IF (nexp_ppl > 0)
THEN
409 polder(0) = polder(0) + cexp_ppl(1)
412 IF (nexp_ppl > 1)
THEN
414 a2 = 0.5_dp/q2*cexp_ppl(2)
415 polder(0) = polder(0) + a2*(2._dp*rho*t + 3._dp*q)
416 polder(1) = polder(1) - a2*2._dp*rho
419 IF (nexp_ppl > 2)
THEN
423 a3 = 0.25_dp/q4*cexp_ppl(3)
424 polder(0) = polder(0) + a3*(4._dp*rho2*t2 + 20._dp*rho*t*q + 15._dp*q2)
425 polder(1) = polder(1) - a3*(8._dp*rho2*t + 20._dp*rho*q)
426 polder(2) = polder(2) + a3*8._dp*rho2
429 IF (nexp_ppl > 3)
THEN
433 a4 = 0.125_dp/q6*cexp_ppl(4)
434 polder(0) = polder(0) + a4*(8._dp*rho3*t3 + 84._dp*rho2*t2*q + 210._dp*rho*t*q2 + 105._dp*q*q2)
435 polder(1) = polder(1) - a4*(24._dp*rho3*t2 + 168._dp*rho2*t*q + 210._dp*rho*q2)
436 polder(2) = polder(2) + a4*(48._dp*rho3*t + 168._dp*rho2*q)
437 polder(3) = polder(3) - a4*48_dp*rho3
440 IF (nexp_ppl > 4)
THEN
441 cpabort(
"nexp_ppl > 4")
445 cc = (
pi/q)**1.5_dp*exp(-t*f)
447 IF (mmax >= 0) expder(0) = cc
449 expder(i) = f*expder(i - 1)
453 DO j = 0, min(i, pmax)
456 auxint(i) = auxint(i) + expder(ke)*polder(kp)*choose(i, j)
460 END SUBROUTINE ppl_aux
467 FUNCTION choose(n, k)
469 INTEGER,
INTENT(IN) :: n, k
470 REAL(kind=
dp) :: choose
473 choose = real(nint(
fac(n)/(
fac(k)*
fac(n - k))), kind=
dp)
Calculation of general three-center integrals over Cartesian Gaussian-type functions and a spherical ...
subroutine, public os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, auxint, rpgfc, rac, dac, va, dva)
Calculation of two-center integrals <a|c> over Cartesian Gaussian functions and a spherical potential...
subroutine, public os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, auxint, rpgfc, rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, vab2, vab2_work, deltaR, iatom, jatom, katom)
Calculation of three-center integrals <a|c|b> over Cartesian Gaussian functions and a spherical poten...
Two-center overlap integrals over Cartesian Gaussian-type functions.
recursive real(dp) function, public os_overlap2(an, bn)
...
subroutine, public init_os_overlap2(ya, yb, rA, rB)
Calculation of overlap integrals over Cartesian Gaussian-type functions.
Calculation of three-center overlap integrals over Cartesian Gaussian-type functions for the second t...
subroutine, public ppl_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, hab2, hab2_work, deltaR, iatom, jatom, katom)
Calculation of three-center overlap integrals <a|c|b> over Cartesian Gaussian functions for the local...
subroutine, public ppl_integral_ri(la_max_set, la_min_set, npgfa, rpgfa, zeta, nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, rac, dac, va, dva)
Calculation of two-center overlap integrals <a|c> over Cartesian Gaussian functions for the local par...
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
real(kind=dp), dimension(0:maxfac), parameter, public fac
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public ncoset
integer, dimension(:, :), allocatable, public indco