(git:ed6f26b)
Loading...
Searching...
No Matches
ai_overlap_ppl.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Calculation of three-center overlap integrals over Cartesian
10!> Gaussian-type functions for the second term V(ppl) of the local
11!> part of the Goedecker pseudopotential (GTH):
12!>
13!> <a|V(local)|b> = <a|V(erf) + V(ppl)|b>
14!> = <a|V(erf)|b> + <a|V(ppl)|b>
15!> = <a|-Z(eff)*erf(SQRT(2)*alpha*r)/r +
16!> (C1 + C2*(alpha*r)**2 + C3*(alpha*r)**4 +
17!> C4*(alpha*r)**6)*exp(-(alpha*r)**2/2))|b>
18!> \par Literature
19!> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
20!> S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B 54, 1703 (1996)
21!> C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B 58, 3641 (1998)
22!> \par History
23!> - Derivatives added (17.05.2002,MK)
24!> - Complete refactoring (05.2011,jhu)
25!> \author Matthias Krack (04.10.2000)
26! **************************************************************************************************
27MODULE ai_overlap_ppl
28 USE ai_oneelectron, ONLY: os_2center,&
29 os_3center
30 USE ai_overlap_debug, ONLY: init_os_overlap2,&
31 os_overlap2
32 USE gamma, ONLY: fgamma => fgamma_0
33 USE gfun, ONLY: gfun_values
34 USE kinds, ONLY: dp
35 USE mathconstants, ONLY: pi
36 USE mathlib, ONLY: binomial
37 USE orbital_pointers, ONLY: indco,&
38 ncoset
39#include "../base/base_uses.f90"
40
41 IMPLICIT NONE
42
43 PRIVATE
44
45 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_overlap_ppl'
46
47! *** Public subroutines ***
48
49 PUBLIC :: ecploc_integral, ppl_integral, ppl_integral_ri
50
51CONTAINS
52
53! **************************************************************************************************
54!> \brief Calculation of three-center overlap integrals <a|c|b> over
55!> Cartesian Gaussian functions for the local part of the Goedecker
56!> pseudopotential (GTH). c is a primitive Gaussian-type function
57!> with a set of even angular momentum indices.
58!>
59!> <a|V(ppl)|b> = <a| (C1 + C2*(alpha*r)**2 + C3*(alpha*r)**4 +
60!> C4*(alpha*r)**6)*exp(-(alpha*r)**2/2))|b>
61!> zetc = alpha**2/2
62!>
63!> \param la_max_set ...
64!> \param la_min_set ...
65!> \param npgfa ...
66!> \param rpgfa ...
67!> \param zeta ...
68!> \param lb_max_set ...
69!> \param lb_min_set ...
70!> \param npgfb ...
71!> \param rpgfb ...
72!> \param zetb ...
73!> \param nexp_ppl ...
74!> \param alpha_ppl ...
75!> \param nct_ppl ...
76!> \param cexp_ppl ...
77!> \param rpgfc ...
78!> \param rab ...
79!> \param dab ...
80!> \param rac ...
81!> \param dac ...
82!> \param rbc ...
83!> \param dbc ...
84!> \param vab ...
85!> \param s ...
86!> \param pab ...
87!> \param force_a ...
88!> \param force_b ...
89!> \param fs ...
90!> \param hab2 The derivative of the ppl integrals according to the weighting factors deltaR
91!> \param hab2_work ...
92!> \param deltaR Weighting factors for the derivatives wrt. nuclear positions
93!> \param iatom ...
94!> \param jatom ...
95!> \param katom ...
96!> \date May 2011
97!> \author Juerg Hutter
98!> \version 1.0
99!> \note Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
100! **************************************************************************************************
101 SUBROUTINE ppl_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
102 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
103 rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
104 hab2, hab2_work, deltaR, iatom, jatom, katom)
105 INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
106 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
107 INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
108 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
109 INTEGER, INTENT(IN) :: nexp_ppl
110 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: alpha_ppl
111 INTEGER, DIMENSION(:), INTENT(IN) :: nct_ppl
112 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: cexp_ppl
113 REAL(kind=dp), INTENT(IN) :: rpgfc
114 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
115 REAL(kind=dp), INTENT(IN) :: dab
116 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
117 REAL(kind=dp), INTENT(IN) :: dac
118 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rbc
119 REAL(kind=dp), INTENT(IN) :: dbc
120 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
121 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: s
122 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
123 OPTIONAL :: pab
124 REAL(kind=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: force_a, force_b
125 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
126 OPTIONAL :: fs, hab2, hab2_work
127 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
128 OPTIONAL :: deltar
129 INTEGER, INTENT(IN), OPTIONAL :: iatom, jatom, katom
130
131 INTEGER :: iexp, ij, ipgf, jpgf, mmax, nexp
132 REAL(kind=dp) :: rho, sab, t, zetc
133 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: auxint
134 REAL(kind=dp), DIMENSION(3) :: pci
135
136 IF (PRESENT(pab)) THEN
137 cpassert(PRESENT(force_a))
138 cpassert(PRESENT(force_b))
139 cpassert(PRESENT(fs))
140 mmax = la_max_set + lb_max_set + 2
141 force_a(:) = 0.0_dp
142 force_b(:) = 0.0_dp
143 ELSE IF (PRESENT(hab2)) THEN
144 mmax = la_max_set + lb_max_set + 2
145 ELSE
146 mmax = la_max_set + lb_max_set
147 END IF
148
149 ALLOCATE (auxint(0:mmax, npgfa*npgfb))
150 auxint = 0._dp
151
152 ! *** Calculate auxiliary integrals ***
153
154 DO ipgf = 1, npgfa
155 ! *** Screening ***
156 IF (rpgfa(ipgf) + rpgfc < dac) cycle
157 DO jpgf = 1, npgfb
158 ! *** Screening ***
159 IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
160 (rpgfa(ipgf) + rpgfb(jpgf) < dab)) cycle
161 ij = (ipgf - 1)*npgfb + jpgf
162 rho = zeta(ipgf) + zetb(jpgf)
163 pci(:) = -(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))/rho
164 sab = exp(-(zeta(ipgf)*zetb(jpgf)/rho*dab*dab))
165 t = rho*sum(pci(:)*pci(:))
166
167 DO iexp = 1, nexp_ppl
168 nexp = nct_ppl(iexp)
169 zetc = alpha_ppl(iexp)
170 CALL ppl_aux(auxint(0:mmax, ij), mmax, t, rho, nexp, cexp_ppl(:, iexp), zetc)
171 END DO
172
173 auxint(0:mmax, ij) = sab*auxint(0:mmax, ij)
174
175 END DO
176 END DO
177
178 CALL os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
179 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, auxint, rpgfc, &
180 rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
181 vab2=hab2, vab2_work=hab2_work, &
182 deltar=deltar, iatom=iatom, jatom=jatom, katom=katom)
183
184 DEALLOCATE (auxint)
185
186 END SUBROUTINE ppl_integral
187
188! **************************************************************************************************
189!> \brief Calculation of three-center potential integrals <a|V(r)|b> over
190!> Cartesian Gaussian functions for the local part of ECP
191!> pseudopotential. Multiple terms C1-4 are possible.
192!>
193!> <a|V(ecploc)|b> = <a| C1/r*exp(-a1*r**2) + C2*exp(-a2*r**2) + C3*r*exp(-a3*r**2) +
194!> C4*r**2*exp(-a4*r**2)|b>
195!>
196!> \param la_max_set ...
197!> \param la_min_set ...
198!> \param npgfa ...
199!> \param rpgfa ...
200!> \param zeta ...
201!> \param lb_max_set ...
202!> \param lb_min_set ...
203!> \param npgfb ...
204!> \param rpgfb ...
205!> \param zetb ...
206!> \param nexp_ppl ...
207!> \param alpha_ppl ...
208!> \param nct_ppl ...
209!> \param cexp_ppl ...
210!> \param rpgfc ...
211!> \param rab ...
212!> \param dab ...
213!> \param rac ...
214!> \param dac ...
215!> \param rbc ...
216!> \param dbc ...
217!> \param vab ...
218!> \param s ...
219!> \param pab ...
220!> \param force_a ...
221!> \param force_b ...
222!> \param fs ...
223!> \param hab2 The derivative of the ppl integrals according to the weighting factors deltaR
224!> \param hab2_work ...
225!> \param deltaR Weighting factors for the derivatives wrt. nuclear positions
226!> \param iatom ...
227!> \param jatom ...
228!> \param katom ...
229!> \date 2025
230!> \author Juerg Hutter
231!> \version 1.0
232! **************************************************************************************************
233 SUBROUTINE ecploc_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
234 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
235 nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
236 rab, dab, rac, dac, rbc, dbc, vab, s, pab, &
237 force_a, force_b, fs, hab2, hab2_work, &
238 deltaR, iatom, jatom, katom)
239 INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
240 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
241 INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
242 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
243 INTEGER, INTENT(IN) :: nexp_ppl
244 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: alpha_ppl
245 INTEGER, DIMENSION(:), INTENT(IN) :: nct_ppl
246 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: cexp_ppl
247 REAL(kind=dp), INTENT(IN) :: rpgfc
248 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
249 REAL(kind=dp), INTENT(IN) :: dab
250 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
251 REAL(kind=dp), INTENT(IN) :: dac
252 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rbc
253 REAL(kind=dp), INTENT(IN) :: dbc
254 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
255 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: s
256 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
257 OPTIONAL :: pab
258 REAL(kind=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: force_a, force_b
259 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
260 OPTIONAL :: fs, hab2, hab2_work
261 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
262 OPTIONAL :: deltar
263 INTEGER, INTENT(IN), OPTIONAL :: iatom, jatom, katom
264
265 INTEGER :: iexp, ij, ipgf, jpgf, mmax, nexp
266 REAL(kind=dp) :: rho, sab, t, zetc
267 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: auxint
268 REAL(kind=dp), DIMENSION(3) :: pci
269
270 IF (PRESENT(pab)) THEN
271 cpassert(PRESENT(force_a))
272 cpassert(PRESENT(force_b))
273 cpassert(PRESENT(fs))
274 mmax = la_max_set + lb_max_set + 2
275 force_a(:) = 0.0_dp
276 force_b(:) = 0.0_dp
277 ELSE IF (PRESENT(hab2)) THEN
278 mmax = la_max_set + lb_max_set + 2
279 ELSE
280 mmax = la_max_set + lb_max_set
281 END IF
282
283 ALLOCATE (auxint(0:mmax, npgfa*npgfb))
284 auxint = 0._dp
285
286 ! *** Calculate auxiliary integrals ***
287
288 DO ipgf = 1, npgfa
289 ! *** Screening ***
290 IF (rpgfa(ipgf) + rpgfc < dac) cycle
291 DO jpgf = 1, npgfb
292 ! *** Screening ***
293 IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
294 (rpgfa(ipgf) + rpgfb(jpgf) < dab)) cycle
295 ij = (ipgf - 1)*npgfb + jpgf
296 rho = zeta(ipgf) + zetb(jpgf)
297 pci(:) = -(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))/rho
298 sab = exp(-(zeta(ipgf)*zetb(jpgf)/rho*dab*dab))
299 t = rho*sum(pci(:)*pci(:))
300
301 DO iexp = 1, nexp_ppl
302 nexp = nct_ppl(iexp)
303 zetc = alpha_ppl(iexp)
304 CALL ecploc_aux(auxint(0:mmax, ij), mmax, t, rho, nexp, cexp_ppl(1, iexp), zetc)
305 END DO
306
307 auxint(0:mmax, ij) = sab*auxint(0:mmax, ij)
308
309 END DO
310 END DO
311
312 CALL os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
313 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, auxint, rpgfc, &
314 rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
315 vab2=hab2, vab2_work=hab2_work, &
316 deltar=deltar, iatom=iatom, jatom=jatom, katom=katom)
317
318 DEALLOCATE (auxint)
319
320 END SUBROUTINE ecploc_integral
321! **************************************************************************************************
322!> \brief Calculation of two-center overlap integrals <a|c> over
323!> Cartesian Gaussian functions for the local part of the Goedecker
324!> pseudopotential (GTH). c is a primitive Gaussian-type function
325!> with a set of even angular momentum indices.
326!>
327!> <a|V(ppl)|b> = <a| (C1 + C2*(alpha*r)**2 + C3*(alpha*r)**4 +
328!> C4*(alpha*r)**6)*exp(-(alpha*r)**2/2))|b>
329!> zetc = alpha**2/2
330!>
331!> \param la_max_set ...
332!> \param la_min_set ...
333!> \param npgfa ...
334!> \param rpgfa ...
335!> \param zeta ...
336!> \param nexp_ppl ...
337!> \param alpha_ppl ...
338!> \param nct_ppl ...
339!> \param cexp_ppl ...
340!> \param rpgfc ...
341!> \param rac ...
342!> \param dac ...
343!> \param va ...
344!> \param dva ...
345!> \date December 2017
346!> \author Juerg Hutter
347!> \version 1.0
348! **************************************************************************************************
349 SUBROUTINE ppl_integral_ri(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
350 nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
351 rac, dac, va, dva)
352 INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
353 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
354 INTEGER, INTENT(IN) :: nexp_ppl
355 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: alpha_ppl
356 INTEGER, DIMENSION(:), INTENT(IN) :: nct_ppl
357 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: cexp_ppl
358 REAL(kind=dp), INTENT(IN) :: rpgfc
359 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
360 REAL(kind=dp), INTENT(IN) :: dac
361 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: va
362 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT), &
363 OPTIONAL :: dva
364
365 INTEGER :: i, iexp, ipgf, iw, mmax, na, nexp
366 INTEGER, DIMENSION(3) :: ani, anm, anp
367 LOGICAL :: debug
368 REAL(kind=dp) :: oint, oref, rho, t, zetc
369 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: auxint
370 REAL(kind=dp), DIMENSION(3) :: doint, doref
371
372 debug = .false.
373
374 IF (PRESENT(dva)) THEN
375 mmax = la_max_set + 1
376 ELSE
377 mmax = la_max_set
378 END IF
379
380 ALLOCATE (auxint(0:mmax, npgfa))
381 auxint = 0._dp
382
383 ! *** Calculate auxiliary integrals ***
384 DO ipgf = 1, npgfa
385 IF (rpgfa(ipgf) + rpgfc < dac) cycle
386 rho = zeta(ipgf)
387 t = rho*dac*dac
388
389 DO iexp = 1, nexp_ppl
390 nexp = nct_ppl(iexp)
391 zetc = alpha_ppl(iexp)
392 CALL ppl_aux(auxint(0:mmax, ipgf), mmax, t, rho, nexp, cexp_ppl(:, iexp), zetc)
393 END DO
394
395 END DO
396
397 IF (PRESENT(dva)) THEN
398 CALL os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
399 auxint, rpgfc, rac, dac, va, dva)
400 ELSE
401 CALL os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
402 auxint, rpgfc, rac, dac, va)
403 END IF
404
405 DEALLOCATE (auxint)
406
407 IF (debug) THEN
408 iw = 6
409 na = 0
410 DO ipgf = 1, npgfa
411 IF (rpgfa(ipgf) + rpgfc < dac) THEN
412 na = na + ncoset(la_max_set)
413 cycle
414 END IF
415 rho = zeta(ipgf)
416 DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
417 oref = va(na + i)
418 ani(1:3) = indco(1:3, i)
419 oint = ppl_ri_test(rho, ani, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
420 ! test
421 IF (abs(oint - oref) > 1.0e-12_dp) THEN
422 WRITE (iw, '(A,3i2,i5,F10.4,2G24.12)') "PPL int error ", ani, la_max_set, dac, oint, oref
423 END IF
424 IF (PRESENT(dva)) THEN
425 anp = ani + (/1, 0, 0/)
426 anm = ani - (/1, 0, 0/)
427 doint(1) = 2._dp*rho*ppl_ri_test(rho, anp, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl) &
428 - ani(1)*ppl_ri_test(rho, anm, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
429 anp = ani + (/0, 1, 0/)
430 anm = ani - (/0, 1, 0/)
431 doint(2) = 2._dp*rho*ppl_ri_test(rho, anp, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl) &
432 - ani(2)*ppl_ri_test(rho, anm, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
433 anp = ani + (/0, 0, 1/)
434 anm = ani - (/0, 0, 1/)
435 doint(3) = 2._dp*rho*ppl_ri_test(rho, anp, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl) &
436 - ani(3)*ppl_ri_test(rho, anm, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl)
437 doref(1:3) = dva(na + i, 1:3)
438 IF (any(abs(doint - doref) > 1.0e-6_dp)) THEN
439 WRITE (iw, '(A,3i2,i5,F10.4,2G24.12)') " PPL dint error ", &
440 ani, la_max_set, dac, sum(abs(doint)), sum(abs(doref))
441 END IF
442 END IF
443 END DO
444 na = na + ncoset(la_max_set)
445 END DO
446 END IF
447
448 END SUBROUTINE ppl_integral_ri
449
450! **************************************************************************************************
451!> \brief ...
452!> \param rho ...
453!> \param ani ...
454!> \param rac ...
455!> \param nexp_ppl ...
456!> \param nct_ppl ...
457!> \param alpha_ppl ...
458!> \param cexp_ppl ...
459!> \return ...
460! **************************************************************************************************
461 FUNCTION ppl_ri_test(rho, ani, rac, nexp_ppl, nct_ppl, alpha_ppl, cexp_ppl) RESULT(oint)
462 REAL(kind=dp), INTENT(IN) :: rho
463 INTEGER, DIMENSION(3), INTENT(IN) :: ani
464 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
465 INTEGER, INTENT(IN) :: nexp_ppl
466 INTEGER, DIMENSION(:), INTENT(IN) :: nct_ppl
467 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: alpha_ppl
468 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: cexp_ppl
469 REAL(kind=dp) :: oint
470
471 INTEGER :: iexp, nexp, ni
472 REAL(kind=dp) :: cn, zetc
473 REAL(kind=dp), DIMENSION(3) :: ra
474
475 oint = 0.0_dp
476 ra = 0.0_dp
477 DO iexp = 1, nexp_ppl
478 nexp = nct_ppl(iexp)
479 zetc = alpha_ppl(iexp)
480 CALL init_os_overlap2(rho, zetc, ra, -rac)
481 DO ni = 1, nexp
482 cn = cexp_ppl(ni, iexp)
483 SELECT CASE (ni)
484 CASE (1)
485 oint = oint + cn*os_overlap2(ani, (/0, 0, 0/))
486 CASE (2)
487 oint = oint + cn*os_overlap2(ani, (/2, 0, 0/))
488 oint = oint + cn*os_overlap2(ani, (/0, 2, 0/))
489 oint = oint + cn*os_overlap2(ani, (/0, 0, 2/))
490 CASE (3)
491 oint = oint + cn*os_overlap2(ani, (/4, 0, 0/))
492 oint = oint + cn*os_overlap2(ani, (/0, 4, 0/))
493 oint = oint + cn*os_overlap2(ani, (/0, 0, 4/))
494 oint = oint + 2.0_dp*cn*os_overlap2(ani, (/2, 2, 0/))
495 oint = oint + 2.0_dp*cn*os_overlap2(ani, (/0, 2, 2/))
496 oint = oint + 2.0_dp*cn*os_overlap2(ani, (/2, 0, 2/))
497 CASE (4)
498 oint = oint + cn*os_overlap2(ani, (/6, 0, 0/))
499 oint = oint + cn*os_overlap2(ani, (/0, 6, 0/))
500 oint = oint + cn*os_overlap2(ani, (/0, 0, 6/))
501 oint = oint + 3.0_dp*cn*os_overlap2(ani, (/4, 2, 0/))
502 oint = oint + 3.0_dp*cn*os_overlap2(ani, (/4, 0, 2/))
503 oint = oint + 3.0_dp*cn*os_overlap2(ani, (/2, 4, 0/))
504 oint = oint + 3.0_dp*cn*os_overlap2(ani, (/0, 4, 2/))
505 oint = oint + 3.0_dp*cn*os_overlap2(ani, (/2, 0, 4/))
506 oint = oint + 3.0_dp*cn*os_overlap2(ani, (/0, 2, 4/))
507 oint = oint + 6.0_dp*cn*os_overlap2(ani, (/2, 2, 2/))
508 CASE DEFAULT
509 cpabort("OVERLAP_PPL")
510 END SELECT
511 END DO
512 END DO
513
514 END FUNCTION ppl_ri_test
515
516! **************************************************************************************************
517!> \brief ...
518!> \param auxint ...
519!> \param mmax ...
520!> \param t ...
521!> \param rho ...
522!> \param nexp_ppl ...
523!> \param cexp_ppl ...
524!> \param zetc ...
525! **************************************************************************************************
526 SUBROUTINE ppl_aux(auxint, mmax, t, rho, nexp_ppl, cexp_ppl, zetc)
527 INTEGER, INTENT(IN) :: mmax
528 REAL(kind=dp), DIMENSION(0:mmax) :: auxint
529 REAL(kind=dp), INTENT(IN) :: t, rho
530 INTEGER, INTENT(IN) :: nexp_ppl
531 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: cexp_ppl
532 REAL(kind=dp), INTENT(IN) :: zetc
533
534 INTEGER :: i, j, ke, kp, pmax
535 REAL(kind=dp) :: a2, a3, a4, cc, f, q, q2, q4, q6, rho2, &
536 rho3, t2, t3
537 REAL(kind=dp), DIMENSION(0:6) :: polder
538 REAL(kind=dp), DIMENSION(0:mmax) :: expder
539
540 cpassert(nexp_ppl > 0)
541 q = rho + zetc
542 polder = 0._dp
543 pmax = 0
544 IF (nexp_ppl > 0) THEN
545 polder(0) = polder(0) + cexp_ppl(1)
546 pmax = 0
547 END IF
548 IF (nexp_ppl > 1) THEN
549 q2 = q*q
550 a2 = 0.5_dp/q2*cexp_ppl(2)
551 polder(0) = polder(0) + a2*(2._dp*rho*t + 3._dp*q)
552 polder(1) = polder(1) - a2*2._dp*rho
553 pmax = 1
554 END IF
555 IF (nexp_ppl > 2) THEN
556 q4 = q2*q2
557 rho2 = rho*rho
558 t2 = t*t
559 a3 = 0.25_dp/q4*cexp_ppl(3)
560 polder(0) = polder(0) + a3*(4._dp*rho2*t2 + 20._dp*rho*t*q + 15._dp*q2)
561 polder(1) = polder(1) - a3*(8._dp*rho2*t + 20._dp*rho*q)
562 polder(2) = polder(2) + a3*8._dp*rho2
563 pmax = 2
564 END IF
565 IF (nexp_ppl > 3) THEN
566 q6 = q4*q2
567 rho3 = rho2*rho
568 t3 = t2*t
569 a4 = 0.125_dp/q6*cexp_ppl(4)
570 polder(0) = polder(0) + a4*(8._dp*rho3*t3 + 84._dp*rho2*t2*q + 210._dp*rho*t*q2 + 105._dp*q*q2)
571 polder(1) = polder(1) - a4*(24._dp*rho3*t2 + 168._dp*rho2*t*q + 210._dp*rho*q2)
572 polder(2) = polder(2) + a4*(48._dp*rho3*t + 168._dp*rho2*q)
573 polder(3) = polder(3) - a4*48_dp*rho3
574 pmax = 3
575 END IF
576 IF (nexp_ppl > 4) THEN
577 cpabort("nexp_ppl > 4")
578 END IF
579
580 f = zetc/q
581 cc = (pi/q)**1.5_dp*exp(-t*f)
582
583 IF (mmax >= 0) expder(0) = cc
584 DO i = 1, mmax
585 expder(i) = f*expder(i - 1)
586 END DO
587
588 DO i = 0, mmax
589 DO j = 0, min(i, pmax)
590 kp = j
591 ke = i - j
592 auxint(i) = auxint(i) + expder(ke)*polder(kp)*binomial(i, j)
593 END DO
594 END DO
595
596 END SUBROUTINE ppl_aux
597! **************************************************************************************************
598!> \brief ...
599!> \param auxint ...
600!> \param mmax ...
601!> \param t ...
602!> \param rho ...
603!> \param nexp ...
604!> \param cexp ...
605!> \param zetc ...
606! **************************************************************************************************
607 SUBROUTINE ecploc_aux(auxint, mmax, t, rho, nexp, cexp, zetc)
608 INTEGER, INTENT(IN) :: mmax
609 REAL(kind=dp), DIMENSION(0:mmax) :: auxint
610 REAL(kind=dp), INTENT(IN) :: t, rho
611 INTEGER, INTENT(IN) :: nexp
612 REAL(kind=dp), INTENT(IN) :: cexp, zetc
613
614 INTEGER :: i, j, ke, kf
615 REAL(kind=dp) :: c0, c1, cc, cval, fa, fr, q, ts
616 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: expder, fdiff, funder, gfund
617
618 q = rho + zetc
619 fa = zetc/q
620 fr = rho/q
621 !
622 ALLOCATE (expder(0:mmax), funder(0:mmax + 1))
623 !
624 SELECT CASE (nexp)
625 CASE (0)
626 cval = 2.0_dp*cexp/sqrt(q)*pi**1.5_dp*exp(-t*fa)
627 expder(0) = cval
628 DO i = 1, mmax
629 expder(i) = fa*expder(i - 1)
630 END DO
631 ts = fr*t
632 ALLOCATE (gfund(0:mmax))
633 CALL gfun_values(mmax, ts, gfund)
634
635 funder(0) = gfund(0)
636 DO i = 1, mmax
637 funder(i) = 0.0_dp
638 DO j = 0, i
639 funder(i) = funder(i) + (-1)**j*binomial(i, j)*gfund(j)
640 END DO
641 END DO
642
643 DEALLOCATE (gfund)
644 DO i = 1, mmax
645 funder(i) = fr**i*funder(i)
646 END DO
647 DO i = 0, mmax
648 DO j = 0, i
649 kf = j
650 ke = i - j
651 auxint(i) = auxint(i) + expder(ke)*funder(kf)*binomial(i, j)
652 END DO
653 END DO
654 CASE (1)
655 cval = cexp*2._dp*pi/q*exp(-t*fa)
656 expder(0) = cval
657 DO i = 1, mmax
658 expder(i) = fa*expder(i - 1)
659 END DO
660 ts = fr*t
661 CALL fgamma(mmax, ts, funder)
662 DO i = 1, mmax
663 funder(i) = fr**i*funder(i)
664 END DO
665 DO i = 0, mmax
666 DO j = 0, i
667 kf = j
668 ke = i - j
669 auxint(i) = auxint(i) + expder(ke)*funder(kf)*binomial(i, j)
670 END DO
671 END DO
672 CASE (2)
673 cval = cexp*(pi/q)**1.5_dp*exp(-t*fa)
674 expder(0) = cval
675 DO i = 1, mmax
676 expder(i) = fa*expder(i - 1)
677 END DO
678 auxint(0:mmax) = auxint(0:mmax) + expder(0:mmax)
679 CASE (3)
680 cval = 2.*pi*cexp/q**2*exp(-t*fa)
681 expder(0) = cval
682 DO i = 1, mmax
683 expder(i) = fa*expder(i - 1)
684 END DO
685 ts = fr*t
686 CALL fgamma(mmax + 1, ts, funder)
687 ALLOCATE (fdiff(0:mmax))
688 fdiff(0) = (1.0_dp + ts)*funder(0) - ts*funder(1)
689 DO i = 1, mmax
690 fdiff(i) = fr**i*(-i*funder(i - 1) + (1.0_dp + ts)*funder(i) &
691 + i*funder(i) - ts*funder(i + 1))
692 END DO
693 DO i = 0, mmax
694 DO j = 0, i
695 kf = j
696 ke = i - j
697 auxint(i) = auxint(i) + expder(ke)*fdiff(kf)*binomial(i, j)
698 END DO
699 END DO
700 DEALLOCATE (fdiff)
701 CASE (4)
702 cval = cexp/(4._dp*q**2)*(pi/q)**1.5_dp*exp(-t*fa)
703 expder(0) = cval
704 DO i = 1, mmax
705 expder(i) = fa*expder(i - 1)
706 END DO
707 c0 = 4._dp*rho/fa
708 c1 = 6._dp*q + 4._dp*rho*t
709 DO i = 0, mmax
710 cc = -i*c0 + c1
711 expder(i) = cc*expder(i)
712 END DO
713 auxint(0:mmax) = auxint(0:mmax) + expder(0:mmax)
714 CASE DEFAULT
715 cpabort("nexp out of range [1..4]")
716 END SELECT
717 !
718 DEALLOCATE (expder, funder)
719
720 END SUBROUTINE ecploc_aux
721! **************************************************************************************************
722
723END MODULE ai_overlap_ppl
ai_oneelectron
Calculation of general three-center integrals over Cartesian Gaussian-type functions and a spherical ...
Definition ai_oneelectron.F:34
ai_oneelectron::os_3center
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...
Definition ai_oneelectron.F:96
ai_oneelectron::os_2center
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...
Definition ai_oneelectron.F:607
ai_overlap_debug
Two-center overlap integrals over Cartesian Gaussian-type functions.
Definition ai_overlap_debug.F:16
ai_overlap_debug::init_os_overlap2
subroutine, public init_os_overlap2(ya, yb, ra, rb)
Calculation of overlap integrals over Cartesian Gaussian-type functions.
Definition ai_overlap_debug.F:47
ai_overlap_debug::os_overlap2
recursive real(dp) function, public os_overlap2(an, bn)
...
Definition ai_overlap_debug.F:73
ai_overlap_debug::zeta
real(dp) zeta
Definition ai_overlap_debug.F:31
ai_overlap_ppl
Calculation of three-center overlap integrals over Cartesian Gaussian-type functions for the second t...
Definition ai_overlap_ppl.F:27
ai_overlap_ppl::ppl_integral
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...
Definition ai_overlap_ppl.F:105
ai_overlap_ppl::ecploc_integral
subroutine, public ecploc_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 potential integrals <a|V(r)|b> over Cartesian Gaussian functions for the ...
Definition ai_overlap_ppl.F:239
ai_overlap_ppl::ppl_integral_ri
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...
Definition ai_overlap_ppl.F:352
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
gamma::fgamma_0
subroutine, public fgamma_0(nmax, t, f)
Calculation of the incomplete Gamma function F(t) for multicenter integrals over Gaussian functions....
Definition gamma.F:154
gfun
Calculation of the G function G_n(t) for 1/R^2 operators.
Definition gfun.F:29
gfun::gfun_values
subroutine, public gfun_values(nmax, t, g)
...
Definition gfun.F:53
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
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::binomial
elemental real(kind=dp) function, public binomial(n, k)
The binomial coefficient n over k for 0 <= k <= n is calculated, otherwise zero is returned.
Definition mathlib.F:206
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::indco
integer, dimension(:, :), allocatable, public indco
Definition orbital_pointers.F:43