Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of the kinetic energy integrals over Cartesian
10 : !> Gaussian-type functions.
11 : !>
12 : !> [a|T|b] = [a|-nabla**2/2|b]
13 : !> \par Literature
14 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
15 : !> \par History
16 : !> - Derivatives added (10.05.2002,MK)
17 : !> - Fully refactored (07.07.2014,JGH)
18 : !> \author Matthias Krack (31.07.2000)
19 : ! **************************************************************************************************
20 : MODULE ai_kinetic
21 : USE ai_os_rr, ONLY: os_rr_ovlp
22 : USE kinds, ONLY: dp
23 : USE mathconstants, ONLY: pi
24 : USE orbital_pointers, ONLY: coset,&
25 : ncoset
26 : #include "../base/base_uses.f90"
27 :
28 : IMPLICIT NONE
29 :
30 : PRIVATE
31 :
32 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_kinetic'
33 :
34 : ! *** Public subroutines ***
35 :
36 : PUBLIC :: kinetic
37 :
38 : CONTAINS
39 :
40 : ! **************************************************************************************************
41 : !> \brief Calculation of the two-center kinetic energy integrals [a|T|b] over
42 : !> Cartesian Gaussian-type functions.
43 : !> \param la_max Maximum L of basis on A
44 : !> \param la_min Minimum L of basis on A
45 : !> \param npgfa Number of primitive functions in set of basis on A
46 : !> \param rpgfa Range of functions on A (used for prescreening)
47 : !> \param zeta Exponents of basis on center A
48 : !> \param lb_max Maximum L of basis on A
49 : !> \param lb_min Minimum L of basis on A
50 : !> \param npgfb Number of primitive functions in set of basis on B
51 : !> \param rpgfb Range of functions on B (used for prescreening)
52 : !> \param zetb Exponents of basis on center B
53 : !> \param rab Distance vector between centers A and B
54 : !> \param kab Kinetic energy integrals, optional
55 : !> \param dab First derivatives of Kinetic energy integrals, optional
56 : !> \date 07.07.2014
57 : !> \author JGH
58 : ! **************************************************************************************************
59 4672817 : SUBROUTINE kinetic(la_max, la_min, npgfa, rpgfa, zeta, &
60 4672817 : lb_max, lb_min, npgfb, rpgfb, zetb, &
61 4672817 : rab, kab, dab)
62 : INTEGER, INTENT(IN) :: la_max, la_min, npgfa
63 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
64 : INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
65 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
66 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
67 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
68 : OPTIONAL :: kab
69 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
70 : OPTIONAL :: dab
71 :
72 : INTEGER :: ax, ay, az, bx, by, bz, coa, cob, ia, &
73 : ib, idx, idy, idz, ipgf, jpgf, la, lb, &
74 : ldrr, lma, lmb, ma, mb, na, nb, ofa, &
75 : ofb
76 : REAL(KIND=dp) :: a, b, dsx, dsy, dsz, dtx, dty, dtz, f0, &
77 : rab2, tab, xhi, zet
78 4672817 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr, tt
79 : REAL(KIND=dp), DIMENSION(3) :: rap, rbp
80 :
81 4672817 : CPASSERT(PRESENT(kab) .OR. PRESENT(dab))
82 :
83 : ! Distance of the centers a and b
84 :
85 4672817 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
86 4672817 : tab = SQRT(rab2)
87 :
88 : ! Maximum l for auxiliary integrals
89 4672817 : IF (PRESENT(kab)) THEN
90 4672817 : lma = la_max + 1
91 4672817 : lmb = lb_max + 1
92 : END IF
93 4672817 : IF (PRESENT(dab)) THEN
94 923077 : lma = la_max + 2
95 923077 : lmb = lb_max + 1
96 923077 : idx = coset(1, 0, 0) - coset(0, 0, 0)
97 923077 : idy = coset(0, 1, 0) - coset(0, 0, 0)
98 923077 : idz = coset(0, 0, 1) - coset(0, 0, 0)
99 : END IF
100 4672817 : ldrr = MAX(lma, lmb) + 1
101 :
102 : ! Allocate space for auxiliary integrals
103 32709719 : ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3), tt(0:ldrr - 1, 0:ldrr - 1, 3))
104 :
105 : ! Number of integrals, check size of arrays
106 4672817 : ofa = ncoset(la_min - 1)
107 4672817 : ofb = ncoset(lb_min - 1)
108 4672817 : na = ncoset(la_max) - ofa
109 4672817 : nb = ncoset(lb_max) - ofb
110 4672817 : IF (PRESENT(kab)) THEN
111 4672817 : CPASSERT((SIZE(kab, 1) >= na*npgfa))
112 4672817 : CPASSERT((SIZE(kab, 2) >= nb*npgfb))
113 : END IF
114 4672817 : IF (PRESENT(dab)) THEN
115 923077 : CPASSERT((SIZE(dab, 1) >= na*npgfa))
116 923077 : CPASSERT((SIZE(dab, 2) >= nb*npgfb))
117 923077 : CPASSERT((SIZE(dab, 3) >= 3))
118 : END IF
119 :
120 : ! Loops over all pairs of primitive Gaussian-type functions
121 4672817 : ma = 0
122 18270302 : DO ipgf = 1, npgfa
123 13597485 : mb = 0
124 61171065 : DO jpgf = 1, npgfb
125 : ! Distance Screening
126 47573580 : IF (rpgfa(ipgf) + rpgfb(jpgf) < tab) THEN
127 866682912 : IF (PRESENT(kab)) kab(ma + 1:ma + na, mb + 1:mb + nb) = 0.0_dp
128 650164644 : IF (PRESENT(dab)) dab(ma + 1:ma + na, mb + 1:mb + nb, 1:3) = 0.0_dp
129 34759668 : mb = mb + nb
130 34759668 : CYCLE
131 : END IF
132 :
133 : ! Calculate some prefactors
134 12813912 : a = zeta(ipgf)
135 12813912 : b = zetb(jpgf)
136 12813912 : zet = a + b
137 12813912 : xhi = a*b/zet
138 51255648 : rap = b*rab/zet
139 51255648 : rbp = -a*rab/zet
140 :
141 : ! [s|s] integral
142 12813912 : f0 = 0.5_dp*(pi/zet)**(1.5_dp)*EXP(-xhi*rab2)
143 :
144 : ! Calculate the recurrence relation, overlap
145 12813912 : CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
146 :
147 : ! kinetic energy auxiliary integrals, overlap of [da/dx|db/dx]
148 39022510 : DO la = 0, lma - 1
149 90152387 : DO lb = 0, lmb - 1
150 51129877 : tt(la, lb, 1) = 4.0_dp*a*b*rr(la + 1, lb + 1, 1)
151 51129877 : tt(la, lb, 2) = 4.0_dp*a*b*rr(la + 1, lb + 1, 2)
152 51129877 : tt(la, lb, 3) = 4.0_dp*a*b*rr(la + 1, lb + 1, 3)
153 51129877 : IF (la > 0 .AND. lb > 0) THEN
154 14230164 : tt(la, lb, 1) = tt(la, lb, 1) + REAL(la*lb, dp)*rr(la - 1, lb - 1, 1)
155 14230164 : tt(la, lb, 2) = tt(la, lb, 2) + REAL(la*lb, dp)*rr(la - 1, lb - 1, 2)
156 14230164 : tt(la, lb, 3) = tt(la, lb, 3) + REAL(la*lb, dp)*rr(la - 1, lb - 1, 3)
157 : END IF
158 51129877 : IF (la > 0) THEN
159 27624850 : tt(la, lb, 1) = tt(la, lb, 1) - 2.0_dp*REAL(la, dp)*b*rr(la - 1, lb + 1, 1)
160 27624850 : tt(la, lb, 2) = tt(la, lb, 2) - 2.0_dp*REAL(la, dp)*b*rr(la - 1, lb + 1, 2)
161 27624850 : tt(la, lb, 3) = tt(la, lb, 3) - 2.0_dp*REAL(la, dp)*b*rr(la - 1, lb + 1, 3)
162 : END IF
163 77338475 : IF (lb > 0) THEN
164 24921279 : tt(la, lb, 1) = tt(la, lb, 1) - 2.0_dp*REAL(lb, dp)*a*rr(la + 1, lb - 1, 1)
165 24921279 : tt(la, lb, 2) = tt(la, lb, 2) - 2.0_dp*REAL(lb, dp)*a*rr(la + 1, lb - 1, 2)
166 24921279 : tt(la, lb, 3) = tt(la, lb, 3) - 2.0_dp*REAL(lb, dp)*a*rr(la + 1, lb - 1, 3)
167 : END IF
168 : END DO
169 : END DO
170 :
171 33574323 : DO lb = lb_min, lb_max
172 65886392 : DO bx = 0, lb
173 98791858 : DO by = 0, lb - bx
174 45719378 : bz = lb - bx - by
175 45719378 : cob = coset(bx, by, bz) - ofb
176 45719378 : ib = mb + cob
177 164703039 : DO la = la_min, la_max
178 277249899 : DO ax = 0, la
179 446855785 : DO ay = 0, la - ax
180 215325264 : az = la - ax - ay
181 215325264 : coa = coset(ax, ay, az) - ofa
182 215325264 : ia = ma + coa
183 : ! integrals
184 215325264 : IF (PRESENT(kab)) THEN
185 : kab(ia, ib) = f0*(tt(ax, bx, 1)*rr(ay, by, 2)*rr(az, bz, 3) + &
186 : rr(ax, bx, 1)*tt(ay, by, 2)*rr(az, bz, 3) + &
187 215325264 : rr(ax, bx, 1)*rr(ay, by, 2)*tt(az, bz, 3))
188 : END IF
189 : ! first derivatives
190 360184193 : IF (PRESENT(dab)) THEN
191 : ! dx
192 49202146 : dsx = 2.0_dp*a*rr(ax + 1, bx, 1)
193 49202146 : IF (ax > 0) dsx = dsx - REAL(ax, dp)*rr(ax - 1, bx, 1)
194 49202146 : dtx = 2.0_dp*a*tt(ax + 1, bx, 1)
195 49202146 : IF (ax > 0) dtx = dtx - REAL(ax, dp)*tt(ax - 1, bx, 1)
196 : dab(ia, ib, idx) = dtx*rr(ay, by, 2)*rr(az, bz, 3) + &
197 49202146 : dsx*(tt(ay, by, 2)*rr(az, bz, 3) + rr(ay, by, 2)*tt(az, bz, 3))
198 : ! dy
199 49202146 : dsy = 2.0_dp*a*rr(ay + 1, by, 2)
200 49202146 : IF (ay > 0) dsy = dsy - REAL(ay, dp)*rr(ay - 1, by, 2)
201 49202146 : dty = 2.0_dp*a*tt(ay + 1, by, 2)
202 49202146 : IF (ay > 0) dty = dty - REAL(ay, dp)*tt(ay - 1, by, 2)
203 : dab(ia, ib, idy) = dty*rr(ax, bx, 1)*rr(az, bz, 3) + &
204 49202146 : dsy*(tt(ax, bx, 1)*rr(az, bz, 3) + rr(ax, bx, 1)*tt(az, bz, 3))
205 : ! dz
206 49202146 : dsz = 2.0_dp*a*rr(az + 1, bz, 3)
207 49202146 : IF (az > 0) dsz = dsz - REAL(az, dp)*rr(az - 1, bz, 3)
208 49202146 : dtz = 2.0_dp*a*tt(az + 1, bz, 3)
209 49202146 : IF (az > 0) dtz = dtz - REAL(az, dp)*tt(az - 1, bz, 3)
210 : dab(ia, ib, idz) = dtz*rr(ax, bx, 1)*rr(ay, by, 2) + &
211 49202146 : dsz*(tt(ax, bx, 1)*rr(ay, by, 2) + rr(ax, bx, 1)*tt(ay, by, 2))
212 : ! scale
213 196808584 : dab(ia, ib, 1:3) = f0*dab(ia, ib, 1:3)
214 : END IF
215 : !
216 : END DO
217 : END DO
218 : END DO !la
219 : END DO
220 : END DO
221 : END DO !lb
222 :
223 26411397 : mb = mb + nb
224 : END DO
225 18270302 : ma = ma + na
226 : END DO
227 :
228 4672817 : DEALLOCATE (rr, tt)
229 :
230 4672817 : END SUBROUTINE kinetic
231 :
232 : END MODULE ai_kinetic
|
