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 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 3721031 : SUBROUTINE kinetic(la_max, la_min, npgfa, rpgfa, zeta, &
60 3721031 : lb_max, lb_min, npgfb, rpgfb, zetb, &
61 3721031 : 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 3721031 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr, tt
79 : REAL(KIND=dp), DIMENSION(3) :: rap, rbp
80 :
81 : ! Distance of the centers a and b
82 :
83 3721031 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
84 3721031 : tab = SQRT(rab2)
85 :
86 : ! Maximum l for auxiliary integrals
87 3721031 : IF (PRESENT(kab)) THEN
88 3721031 : lma = la_max + 1
89 3721031 : lmb = lb_max + 1
90 : END IF
91 3721031 : IF (PRESENT(dab)) THEN
92 662936 : lma = la_max + 2
93 662936 : lmb = lb_max + 1
94 662936 : idx = coset(1, 0, 0) - coset(0, 0, 0)
95 662936 : idy = coset(0, 1, 0) - coset(0, 0, 0)
96 662936 : idz = coset(0, 0, 1) - coset(0, 0, 0)
97 : END IF
98 3721031 : ldrr = MAX(lma, lmb) + 1
99 :
100 : ! Allocate space for auxiliary integrals
101 26047217 : ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3), tt(0:ldrr - 1, 0:ldrr - 1, 3))
102 :
103 : ! Number of integrals, check size of arrays
104 3721031 : ofa = ncoset(la_min - 1)
105 3721031 : ofb = ncoset(lb_min - 1)
106 3721031 : na = ncoset(la_max) - ofa
107 3721031 : nb = ncoset(lb_max) - ofb
108 3721031 : IF (PRESENT(kab)) THEN
109 3721031 : CPASSERT((SIZE(kab, 1) >= na*npgfa))
110 3721031 : CPASSERT((SIZE(kab, 2) >= nb*npgfb))
111 : END IF
112 3721031 : IF (PRESENT(dab)) THEN
113 662936 : CPASSERT((SIZE(dab, 1) >= na*npgfa))
114 662936 : CPASSERT((SIZE(dab, 2) >= nb*npgfb))
115 662936 : CPASSERT((SIZE(dab, 3) >= 3))
116 : END IF
117 :
118 : ! Loops over all pairs of primitive Gaussian-type functions
119 3721031 : ma = 0
120 13967394 : DO ipgf = 1, npgfa
121 10246363 : mb = 0
122 45086970 : DO jpgf = 1, npgfb
123 : ! Distance Screening
124 34840607 : IF (rpgfa(ipgf) + rpgfb(jpgf) < tab) THEN
125 650775257 : IF (PRESENT(kab)) kab(ma + 1:ma + na, mb + 1:mb + nb) = 0.0_dp
126 469761930 : IF (PRESENT(dab)) dab(ma + 1:ma + na, mb + 1:mb + nb, 1:3) = 0.0_dp
127 25177407 : mb = mb + nb
128 25177407 : CYCLE
129 : END IF
130 :
131 : ! Calculate some prefactors
132 9663200 : a = zeta(ipgf)
133 9663200 : b = zetb(jpgf)
134 9663200 : zet = a + b
135 9663200 : xhi = a*b/zet
136 38652800 : rap = b*rab/zet
137 38652800 : rbp = -a*rab/zet
138 :
139 : ! [s|s] integral
140 9663200 : f0 = 0.5_dp*(pi/zet)**(1.5_dp)*EXP(-xhi*rab2)
141 :
142 : ! Calculate the recurrence relation, overlap
143 9663200 : CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
144 :
145 : ! kinetic energy auxiliary integrals, overlap of [da/dx|db/dx]
146 29160395 : DO la = 0, lma - 1
147 67279612 : DO lb = 0, lmb - 1
148 38119217 : tt(la, lb, 1) = 4.0_dp*a*b*rr(la + 1, lb + 1, 1)
149 38119217 : tt(la, lb, 2) = 4.0_dp*a*b*rr(la + 1, lb + 1, 2)
150 38119217 : tt(la, lb, 3) = 4.0_dp*a*b*rr(la + 1, lb + 1, 3)
151 38119217 : IF (la > 0 .AND. lb > 0) THEN
152 10740928 : tt(la, lb, 1) = tt(la, lb, 1) + REAL(la*lb, dp)*rr(la - 1, lb - 1, 1)
153 10740928 : tt(la, lb, 2) = tt(la, lb, 2) + REAL(la*lb, dp)*rr(la - 1, lb - 1, 2)
154 10740928 : tt(la, lb, 3) = tt(la, lb, 3) + REAL(la*lb, dp)*rr(la - 1, lb - 1, 3)
155 : END IF
156 38119217 : IF (la > 0) THEN
157 20574923 : tt(la, lb, 1) = tt(la, lb, 1) - 2.0_dp*REAL(la, dp)*b*rr(la - 1, lb + 1, 1)
158 20574923 : tt(la, lb, 2) = tt(la, lb, 2) - 2.0_dp*REAL(la, dp)*b*rr(la - 1, lb + 1, 2)
159 20574923 : tt(la, lb, 3) = tt(la, lb, 3) - 2.0_dp*REAL(la, dp)*b*rr(la - 1, lb + 1, 3)
160 : END IF
161 57616412 : IF (lb > 0) THEN
162 18622022 : tt(la, lb, 1) = tt(la, lb, 1) - 2.0_dp*REAL(lb, dp)*a*rr(la + 1, lb - 1, 1)
163 18622022 : tt(la, lb, 2) = tt(la, lb, 2) - 2.0_dp*REAL(lb, dp)*a*rr(la + 1, lb - 1, 2)
164 18622022 : tt(la, lb, 3) = tt(la, lb, 3) - 2.0_dp*REAL(lb, dp)*a*rr(la + 1, lb - 1, 3)
165 : END IF
166 : END DO
167 : END DO
168 :
169 24915745 : DO lb = lb_min, lb_max
170 48897158 : DO bx = 0, lb
171 73713398 : DO by = 0, lb - bx
172 34479440 : bz = lb - bx - by
173 34479440 : cob = coset(bx, by, bz) - ofb
174 34479440 : ib = mb + cob
175 125033484 : DO la = la_min, la_max
176 217064879 : DO ax = 0, la
177 362399143 : DO ay = 0, la - ax
178 179813704 : az = la - ax - ay
179 179813704 : coa = coset(ax, ay, az) - ofa
180 179813704 : ia = ma + coa
181 : ! integrals
182 179813704 : IF (PRESENT(kab)) THEN
183 : kab(ia, ib) = f0*(tt(ax, bx, 1)*rr(ay, by, 2)*rr(az, bz, 3) + &
184 : rr(ax, bx, 1)*tt(ay, by, 2)*rr(az, bz, 3) + &
185 179813704 : rr(ax, bx, 1)*rr(ay, by, 2)*tt(az, bz, 3))
186 : END IF
187 : ! first derivatives
188 295826512 : IF (PRESENT(dab)) THEN
189 : ! dx
190 42720886 : dsx = 2.0_dp*a*rr(ax + 1, bx, 1)
191 42720886 : IF (ax > 0) dsx = dsx - REAL(ax, dp)*rr(ax - 1, bx, 1)
192 42720886 : dtx = 2.0_dp*a*tt(ax + 1, bx, 1)
193 42720886 : IF (ax > 0) dtx = dtx - REAL(ax, dp)*tt(ax - 1, bx, 1)
194 : dab(ia, ib, idx) = dtx*rr(ay, by, 2)*rr(az, bz, 3) + &
195 42720886 : dsx*(tt(ay, by, 2)*rr(az, bz, 3) + rr(ay, by, 2)*tt(az, bz, 3))
196 : ! dy
197 42720886 : dsy = 2.0_dp*a*rr(ay + 1, by, 2)
198 42720886 : IF (ay > 0) dsy = dsy - REAL(ay, dp)*rr(ay - 1, by, 2)
199 42720886 : dty = 2.0_dp*a*tt(ay + 1, by, 2)
200 42720886 : IF (ay > 0) dty = dty - REAL(ay, dp)*tt(ay - 1, by, 2)
201 : dab(ia, ib, idy) = dty*rr(ax, bx, 1)*rr(az, bz, 3) + &
202 42720886 : dsy*(tt(ax, bx, 1)*rr(az, bz, 3) + rr(ax, bx, 1)*tt(az, bz, 3))
203 : ! dz
204 42720886 : dsz = 2.0_dp*a*rr(az + 1, bz, 3)
205 42720886 : IF (az > 0) dsz = dsz - REAL(az, dp)*rr(az - 1, bz, 3)
206 42720886 : dtz = 2.0_dp*a*tt(az + 1, bz, 3)
207 42720886 : IF (az > 0) dtz = dtz - REAL(az, dp)*tt(az - 1, bz, 3)
208 : dab(ia, ib, idz) = dtz*rr(ax, bx, 1)*rr(ay, by, 2) + &
209 42720886 : dsz*(tt(ax, bx, 1)*rr(ay, by, 2) + rr(ax, bx, 1)*tt(ay, by, 2))
210 : ! scale
211 170883544 : dab(ia, ib, 1:3) = f0*dab(ia, ib, 1:3)
212 : END IF
213 : !
214 : END DO
215 : END DO
216 : END DO !la
217 : END DO
218 : END DO
219 : END DO !lb
220 :
221 19909563 : mb = mb + nb
222 : END DO
223 13967394 : ma = ma + na
224 : END DO
225 :
226 3721031 : DEALLOCATE (rr, tt)
227 :
228 3721031 : END SUBROUTINE kinetic
229 :
230 : END MODULE ai_kinetic
|
