Line data Source code
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 Three-center integrals over Cartesian Gaussian-type functions
10 : !> \par Literature
11 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
12 : !> \par History
13 : !> none
14 : !> \author Dorothea Golze
15 : ! **************************************************************************************************
16 : MODULE ai_overlap3_debug
17 :
18 : USE kinds, ONLY: dp
19 : USE mathconstants, ONLY: pi
20 : #include "../base/base_uses.f90"
21 :
22 : IMPLICIT NONE
23 :
24 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_overlap3_debug'
25 :
26 : INTEGER, PARAMETER :: lmax = 5
27 :
28 : REAL(dp) :: xa, xb, xc
29 : REAL(dp), DIMENSION(3) :: A, B, C
30 : REAL(dp), DIMENSION(3) :: P, G
31 : REAL(dp) :: xsi, zeta, sss
32 :
33 : PRIVATE
34 : PUBLIC :: init_os_overlap3, os_overlap3
35 :
36 : CONTAINS
37 :
38 : ! **************************************************************************************************
39 : !> \brief Calculation of three-center integrals over
40 : !> Cartesian Gaussian-type functions.
41 : !> \param ya ...
42 : !> \param yb ...
43 : !> \param yc ...
44 : !> \param rA ...
45 : !> \param rB ...
46 : !> \param rC ...
47 : ! **************************************************************************************************
48 686 : SUBROUTINE init_os_overlap3(ya, yb, yc, rA, rB, rC)
49 : REAL(dp) :: ya, yb, yc
50 : REAL(dp), DIMENSION(3) :: rA, rB, rC
51 :
52 : REAL(dp) :: fpc, ss
53 :
54 686 : xa = ya
55 686 : xb = yb
56 686 : xc = yc
57 686 : A = rA
58 686 : B = rB
59 686 : C = rC
60 :
61 686 : xsi = xa + xb
62 686 : zeta = xa*xb/xsi
63 :
64 2744 : P = (xa*A + xb*B)/xsi
65 2744 : G = (xsi*P + xc*C)/(xsi + xc)
66 :
67 2744 : ss = (pi/xsi)**(3._dp/2._dp)*EXP(-zeta*SUM((A - B)**2))
68 :
69 2744 : fpc = EXP(-xsi*xc/(xsi + xc)*SUM((P - C)**2))
70 686 : sss = (xsi/(xsi + xc))**(3._dp/2._dp)*ss*fpc
71 :
72 686 : END SUBROUTINE init_os_overlap3
73 :
74 : ! **************************************************************************************************
75 :
76 : ! **************************************************************************************************
77 : !> \brief ...
78 : !> \param an ...
79 : !> \param cn ...
80 : !> \param bn ...
81 : !> \return ...
82 : ! **************************************************************************************************
83 6277096 : RECURSIVE FUNCTION os_overlap3(an, cn, bn) RESULT(IACB)
84 : INTEGER, DIMENSION(3) :: an, cn, bn
85 : REAL(dp) :: IACB
86 :
87 : INTEGER, DIMENSION(3), PARAMETER :: i1 = [1, 0, 0], i2 = [0, 1, 0], &
88 : i3 = [0, 0, 1]
89 :
90 6277096 : IACB = 0._dp
91 21779324 : IF (ANY(an < 0)) RETURN
92 18224472 : IF (ANY(bn < 0)) RETURN
93 18224472 : IF (ANY(cn < 0)) RETURN
94 :
95 18224472 : IF (SUM(an + cn + bn) == 0) THEN
96 1466178 : IACB = sss
97 1466178 : RETURN
98 : END IF
99 :
100 3089940 : IF (bn(1) > 0) THEN
101 243040 : IACB = os_overlap3(an, cn + i1, bn - i1) + (C(1) - B(1))*os_overlap3(an, cn, bn - i1)
102 3065636 : ELSEIF (bn(2) > 0) THEN
103 243040 : IACB = os_overlap3(an, cn + i2, bn - i2) + (C(2) - B(2))*os_overlap3(an, cn, bn - i2)
104 3041332 : ELSEIF (bn(3) > 0) THEN
105 243040 : IACB = os_overlap3(an, cn + i3, bn - i3) + (C(3) - B(3))*os_overlap3(an, cn, bn - i3)
106 : ELSE
107 3017028 : IF (cn(1) > 0) THEN
108 1685600 : IACB = os_overlap3(an + i1, cn - i1, bn) + (A(1) - C(1))*os_overlap3(an, cn - i1, bn)
109 2848468 : ELSEIF (cn(2) > 0) THEN
110 2469600 : IACB = os_overlap3(an + i2, cn - i2, bn) + (A(2) - C(2))*os_overlap3(an, cn - i2, bn)
111 2601508 : ELSEIF (cn(3) > 0) THEN
112 3433920 : IACB = os_overlap3(an + i3, cn - i3, bn) + (A(3) - C(3))*os_overlap3(an, cn - i3, bn)
113 : ELSE
114 2258116 : IF (an(1) > 0) THEN
115 : IACB = (G(1) - A(1))*os_overlap3(an - i1, cn, bn) + &
116 4819836 : 0.5_dp*(an(1) - 1)/(xsi + xc)*os_overlap3(an - i1 - i1, cn, bn)
117 1569568 : ELSEIF (an(2) > 0) THEN
118 : IACB = (G(2) - A(2))*os_overlap3(an - i2, cn, bn) + &
119 5267108 : 0.5_dp*(an(2) - 1)/(xsi + xc)*os_overlap3(an - i2 - i2, cn, bn)
120 817124 : ELSEIF (an(3) > 0) THEN
121 : IACB = (G(3) - A(3))*os_overlap3(an - i3, cn, bn) + &
122 5719868 : 0.5_dp*(an(3) - 1)/(xsi + xc)*os_overlap3(an - i3 - i3, cn, bn)
123 : ELSE
124 0 : CPABORT("I(0000)")
125 : END IF
126 : END IF
127 : END IF
128 :
129 : END FUNCTION os_overlap3
130 :
131 : ! **************************************************************************************************
132 :
133 : END MODULE ai_overlap3_debug
|
