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 general three-center integrals over Cartesian
10 : !> Gaussian-type functions and a spherical operator centered at position C
11 : !>
12 : !> <a|V(local)|b> = <a|F(|r-C|)|b>
13 : !> \par Literature
14 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
15 : !> \par History
16 : !> - Based in part on code by MK
17 : !> \par Parameters
18 : !> - ax,ay,az : Angular momentum index numbers of orbital a.
19 : !> - bx,by,bz : Angular momentum index numbers of orbital b.
20 : !> - coset : Cartesian orbital set pointer.
21 : !> - dab : Distance between the atomic centers a and b.
22 : !> - dac : Distance between the atomic centers a and c.
23 : !> - dbc : Distance between the atomic centers b and c.
24 : !> - l{a,b} : Angular momentum quantum number of shell a or b.
25 : !> - l{a,b}_max : Maximum angular momentum quantum number of shell a or b.
26 : !> - ncoset : Number of Cartesian orbitals up to l.
27 : !> - rab : Distance vector between the atomic centers a and b.
28 : !> - rac : Distance vector between the atomic centers a and c.
29 : !> - rbc : Distance vector between the atomic centers b and c.
30 : !> - rpgf{a,b,c}: Radius of the primitive Gaussian-type function a or b.
31 : !> - zet{a,b} : Exponents of the Gaussian-type functions a or b.
32 : !> \author jhu (05.2011)
33 : ! **************************************************************************************************
34 : MODULE ai_oneelectron
35 :
36 : USE ai_cartesian_operators, ONLY: cartesian_ladder,&
37 : cartesian_ladder_map
38 : USE kinds, ONLY: dp
39 : USE orbital_pointers, ONLY: coset,&
40 : ncoset
41 : #include "../base/base_uses.f90"
42 :
43 : IMPLICIT NONE
44 :
45 : PRIVATE
46 :
47 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_oneelectron'
48 :
49 : ! *** Public subroutines ***
50 :
51 : PUBLIC :: os_3center, os_2center
52 :
53 : CONTAINS
54 :
55 : ! **************************************************************************************************
56 : !> \brief Calculation of three-center integrals <a|c|b> over
57 : !> Cartesian Gaussian functions and a spherical potential
58 : !>
59 : !> \param la_max_set ...
60 : !> \param la_min_set ...
61 : !> \param npgfa ...
62 : !> \param rpgfa ...
63 : !> \param zeta ...
64 : !> \param lb_max_set ...
65 : !> \param lb_min_set ...
66 : !> \param npgfb ...
67 : !> \param rpgfb ...
68 : !> \param zetb ...
69 : !> \param auxint ...
70 : !> \param rpgfc ...
71 : !> \param rab ...
72 : !> \param dab ...
73 : !> \param rac ...
74 : !> \param dac ...
75 : !> \param rbc ...
76 : !> \param dbc ...
77 : !> \param vab ...
78 : !> \param s ...
79 : !> \param pab ...
80 : !> \param force_a ...
81 : !> \param force_b ...
82 : !> \param fs ...
83 : !> \param vab2 The derivative of the 3-center integrals according to the weighting factors.
84 : !> \param vab2_work ...
85 : !> \param deltaR DIMENSION(3, natoms), weighting factors of the derivatives for each atom and direction
86 : !> \param iatom ...
87 : !> \param jatom ...
88 : !> \param katom ...
89 : !> \date May 2011
90 : !> \author Juerg Hutter
91 : !> \version 1.0
92 : !> \note Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
93 : ! **************************************************************************************************
94 51941909 : SUBROUTINE os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
95 51941909 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, auxint, rpgfc, &
96 103883818 : rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
97 103883818 : vab2, vab2_work, deltaR, iatom, jatom, katom)
98 : INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
99 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
100 : INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
101 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
102 : REAL(KIND=dp), DIMENSION(0:, :), INTENT(IN) :: auxint
103 : REAL(KIND=dp), INTENT(IN) :: rpgfc
104 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
105 : REAL(KIND=dp), INTENT(IN) :: dab
106 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
107 : REAL(KIND=dp), INTENT(IN) :: dac
108 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rbc
109 : REAL(KIND=dp), INTENT(IN) :: dbc
110 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
111 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: s
112 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
113 : OPTIONAL :: pab
114 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: force_a, force_b
115 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
116 : OPTIONAL :: fs, vab2, vab2_work
117 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
118 : OPTIONAL :: deltaR
119 : INTEGER, INTENT(IN), OPTIONAL :: iatom, jatom, katom
120 :
121 : INTEGER :: a_first, a_last, b_first, b_last, da_max, db_max, i, ia, iap, iax, iay, iaz, ib, &
122 : ibm, ibx, iby, ibz, idir, ii(3), iim(3), ij, ipgf, ir, ir1, ir2, irm(3), irr(3), irx, &
123 : iry, irz, ix, ixx(1), j, jj(3), jjp(3), jpgf, la_max, la_min, lb_max, lb_min, llr, m, ma, &
124 : mb, mmax, na, nb
125 51941909 : INTEGER, ALLOCATABLE :: map_a(:, :, :), map_b(:, :, :)
126 51941909 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: iiap
127 : LOGICAL :: calculate_force_a, calculate_force_b
128 : REAL(KIND=dp) :: aai, abx, orho, rho, s1, s2
129 : REAL(KIND=dp), DIMENSION(3) :: pai, pbi, pci
130 :
131 51941909 : IF (PRESENT(pab)) THEN
132 14731926 : CPASSERT(PRESENT(fs))
133 14731926 : IF (PRESENT(force_a)) THEN
134 : calculate_force_a = .TRUE.
135 : ELSE
136 0 : calculate_force_a = .FALSE.
137 : END IF
138 14731926 : IF (PRESENT(force_b)) THEN
139 : calculate_force_b = .TRUE.
140 : ELSE
141 0 : calculate_force_b = .FALSE.
142 : END IF
143 : ELSE
144 : calculate_force_a = .FALSE.
145 : calculate_force_b = .FALSE.
146 : END IF
147 :
148 14731926 : IF (calculate_force_a) THEN
149 14731926 : da_max = 1
150 14731926 : force_a = 0.0_dp
151 : ELSE
152 : da_max = 0
153 : END IF
154 :
155 51941909 : IF (calculate_force_b) THEN
156 14731926 : db_max = 1
157 14731926 : force_b = 0.0_dp
158 : ELSE
159 : db_max = 0
160 : END IF
161 :
162 51941909 : IF (PRESENT(vab2)) THEN
163 870 : da_max = 1
164 870 : db_max = 1
165 : END IF
166 :
167 51941909 : la_max = la_max_set + da_max
168 51941909 : la_min = MAX(0, la_min_set - da_max)
169 :
170 51941909 : lb_max = lb_max_set + db_max
171 51941909 : lb_min = MAX(0, lb_min_set - db_max)
172 :
173 51941909 : mmax = la_max + lb_max
174 :
175 51941909 : a_first = ncoset(la_min_set - 1) + 1
176 51941909 : a_last = ncoset(la_max_set)
177 51941909 : b_first = ncoset(lb_min_set - 1) + 1
178 51941909 : b_last = ncoset(lb_max_set)
179 51941909 : IF (da_max == 1) THEN
180 58931184 : ALLOCATE (map_a(4, a_last - a_first + 1, 3))
181 58931184 : DO idir = 1, 3
182 58931184 : CALL cartesian_ladder_map(la_max_set, la_min_set, idir, map_a(:, :, idir), source_lmin=0)
183 : END DO
184 : END IF
185 51941909 : IF (db_max == 1) THEN
186 58931184 : ALLOCATE (map_b(4, b_last - b_first + 1, 3))
187 58931184 : DO idir = 1, 3
188 58931184 : CALL cartesian_ladder_map(lb_max_set, lb_min_set, idir, map_b(:, :, idir), source_lmin=0)
189 : END DO
190 : END IF
191 :
192 : ! precalculate indices for horizontal recursion
193 155825727 : ALLOCATE (iiap(ncoset(mmax), 3))
194 240872336 : DO ma = 0, mmax
195 727755535 : DO iax = 0, ma
196 1728354831 : DO iay = 0, ma - iax
197 1052541205 : iaz = ma - iax - iay
198 1052541205 : ia = coset(iax, iay, iaz)
199 1052541205 : jj(1) = iax; jj(2) = iay; jj(3) = iaz
200 1052541205 : jjp = jj
201 1052541205 : jjp(1) = jjp(1) + 1
202 1052541205 : iap = coset(jjp(1), jjp(2), jjp(3))
203 1052541205 : iiap(ia, 1) = iap
204 1052541205 : jjp = jj
205 1052541205 : jjp(2) = jjp(2) + 1
206 1052541205 : iap = coset(jjp(1), jjp(2), jjp(3))
207 1052541205 : iiap(ia, 2) = iap
208 1052541205 : jjp = jj
209 1052541205 : jjp(3) = jjp(3) + 1
210 1052541205 : iap = coset(jjp(1), jjp(2), jjp(3))
211 1539424404 : iiap(ia, 3) = iap
212 : END DO
213 : END DO
214 : END DO
215 :
216 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
217 :
218 51941909 : na = 0
219 :
220 241692029 : DO ipgf = 1, npgfa
221 :
222 : ! *** Screening ***
223 189750120 : IF (rpgfa(ipgf) + rpgfc < dac) THEN
224 122068983 : na = na + ncoset(la_max_set)
225 122068983 : CYCLE
226 : END IF
227 :
228 67681137 : nb = 0
229 :
230 322539250 : DO jpgf = 1, npgfb
231 :
232 : ! *** Screening ***
233 254858113 : IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
234 : (rpgfa(ipgf) + rpgfb(jpgf) < dab)) THEN
235 165216536 : nb = nb + ncoset(lb_max_set)
236 165216536 : CYCLE
237 : END IF
238 :
239 : ! *** Calculate some prefactors ***
240 89641577 : rho = zeta(ipgf) + zetb(jpgf)
241 358566308 : pai(:) = zetb(jpgf)/rho*rab(:)
242 358566308 : pbi(:) = -zeta(ipgf)/rho*rab(:)
243 358566308 : pci(:) = -(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))/rho
244 89641577 : orho = 0.5_dp/rho
245 :
246 89641577 : ij = (ipgf - 1)*npgfb + jpgf
247 418339710 : s(1, 1, 1:mmax + 1) = auxint(0:mmax, ij)
248 :
249 89641577 : IF (la_max > 0) THEN
250 : ! *** Recurrence steps: [s|c|s] -> [a|c|s] ***
251 : ! *** [a|c|s](m) = (Pi - Ai)*[a-1i|c|s](m) - ***
252 : ! *** (Pi - Ci)*[a-1i|c|s](m+1)) + ***
253 : ! *** Ni(a-1i)/2(a+b)*[a-2i|c|s](m) - ***
254 : ! *** Ni(a-1i)/2(a+b)*[a-2i|c|s](m+1) ***
255 317634750 : DO llr = 1, mmax
256 317634750 : IF (llr == 1) THEN
257 317634750 : DO m = 0, mmax - llr
258 236491181 : s1 = s(1, 1, m + 1)
259 236491181 : s2 = s(1, 1, m + 2)
260 236491181 : s(2, 1, m + 1) = pai(1)*s1 - pci(1)*s2 ! [px|o|s]
261 236491181 : s(3, 1, m + 1) = pai(2)*s1 - pci(2)*s2 ! [py|o|s]
262 317634750 : s(4, 1, m + 1) = pai(3)*s1 - pci(3)*s2 ! [pz|o|s]
263 : END DO
264 155347612 : ELSE IF (llr == 2) THEN
265 234424026 : DO m = 0, mmax - llr
266 155347612 : s1 = s(1, 1, m + 1) - s(1, 1, m + 2)
267 155347612 : s(5, 1, m + 1) = pai(1)*s(2, 1, m + 1) - pci(1)*s(2, 1, m + 2) + orho*s1 ! [dx2|o|s]
268 155347612 : s(6, 1, m + 1) = pai(1)*s(3, 1, m + 1) - pci(1)*s(3, 1, m + 2) ! [dxy|o|s]
269 155347612 : s(7, 1, m + 1) = pai(1)*s(4, 1, m + 1) - pci(1)*s(4, 1, m + 2) ! [dxz|o|s]
270 155347612 : s(8, 1, m + 1) = pai(2)*s(3, 1, m + 1) - pci(2)*s(3, 1, m + 2) + orho*s1 ! [dy2|o|s]
271 155347612 : s(9, 1, m + 1) = pai(2)*s(4, 1, m + 1) - pci(2)*s(4, 1, m + 2) ! [dyz|o|s]
272 234424026 : s(10, 1, m + 1) = pai(3)*s(4, 1, m + 1) - pci(3)*s(4, 1, m + 2) + orho*s1 ! [dz2|o|s]
273 : END DO
274 76271198 : ELSE IF (llr == 3) THEN
275 114733183 : DO m = 0, mmax - llr
276 : s(11, 1, m + 1) = pai(1)*s(5, 1, m + 1) - pci(1)*s(5, 1, m + 2) & ! [fx3 |o|s]
277 76271198 : + 2._dp*orho*(s(2, 1, m + 1) - s(2, 1, m + 2))
278 : s(12, 1, m + 1) = pai(1)*s(6, 1, m + 1) - pci(1)*s(6, 1, m + 2) & ! [fx2y|o|s]
279 76271198 : + orho*(s(3, 1, m + 1) - s(3, 1, m + 2))
280 : s(13, 1, m + 1) = pai(1)*s(7, 1, m + 1) - pci(1)*s(7, 1, m + 2) & ! [fx2z|o|s]
281 76271198 : + orho*(s(4, 1, m + 1) - s(4, 1, m + 2))
282 : s(14, 1, m + 1) = pai(2)*s(6, 1, m + 1) - pci(2)*s(6, 1, m + 2) & ! [fxy2|o|s]
283 76271198 : + orho*(s(2, 1, m + 1) - s(2, 1, m + 2))
284 76271198 : s(15, 1, m + 1) = pai(1)*s(9, 1, m + 1) - pci(1)*s(9, 1, m + 2) ! [fxyz|o|s]
285 : s(16, 1, m + 1) = pai(3)*s(7, 1, m + 1) - pci(3)*s(7, 1, m + 2) & ! [fxz2|o|s]
286 76271198 : + orho*(s(2, 1, m + 1) - s(2, 1, m + 2))
287 : s(17, 1, m + 1) = pai(2)*s(8, 1, m + 1) - pci(2)*s(8, 1, m + 2) & ! [fy3 |o|s]
288 76271198 : + 2._dp*orho*(s(3, 1, m + 1) - s(3, 1, m + 2))
289 : s(18, 1, m + 1) = pai(2)*s(9, 1, m + 1) - pci(2)*s(9, 1, m + 2) & ! [fy2z|o|s]
290 76271198 : + orho*(s(4, 1, m + 1) - s(4, 1, m + 2))
291 : s(19, 1, m + 1) = pai(3)*s(9, 1, m + 1) - pci(3)*s(9, 1, m + 2) & ! [fyz2|o|s]
292 76271198 : + orho*(s(3, 1, m + 1) - s(3, 1, m + 2))
293 : s(20, 1, m + 1) = pai(3)*s(10, 1, m + 1) - pci(3)*s(10, 1, m + 2) & ! [fz3 |o|s]
294 114733183 : + 2._dp*orho*(s(4, 1, m + 1) - s(4, 1, m + 2))
295 : END DO
296 37809213 : ELSE IF (llr == 4) THEN
297 67259332 : DO m = 0, mmax - llr
298 : s(21, 1, m + 1) = pai(1)*s(11, 1, m + 1) - pci(1)*s(11, 1, m + 2) & ! [gx4 |s|s]
299 37809213 : + 3._dp*orho*(s(5, 1, m + 1) - s(5, 1, m + 2))
300 : s(22, 1, m + 1) = pai(1)*s(12, 1, m + 1) - pci(1)*s(12, 1, m + 2) & ! [gx3y |s|s]
301 37809213 : + 2._dp*orho*(s(6, 1, m + 1) - s(6, 1, m + 2))
302 : s(23, 1, m + 1) = pai(1)*s(13, 1, m + 1) - pci(1)*s(13, 1, m + 2) & ! [gx3z |s|s]
303 37809213 : + 2._dp*orho*(s(7, 1, m + 1) - s(7, 1, m + 2))
304 : s(24, 1, m + 1) = pai(1)*s(14, 1, m + 1) - pci(1)*s(14, 1, m + 2) & ! [gx2y2|s|s]
305 37809213 : + orho*(s(8, 1, m + 1) - s(8, 1, m + 2))
306 : s(25, 1, m + 1) = pai(1)*s(15, 1, m + 1) - pci(1)*s(15, 1, m + 2) & ! [gx2yz|s|s]
307 37809213 : + orho*(s(9, 1, m + 1) - s(9, 1, m + 2))
308 : s(26, 1, m + 1) = pai(1)*s(16, 1, m + 1) - pci(1)*s(16, 1, m + 2) & ! [gx2z2|s|s]
309 37809213 : + orho*(s(10, 1, m + 1) - s(10, 1, m + 2))
310 37809213 : s(27, 1, m + 1) = pai(1)*s(17, 1, m + 1) - pci(1)*s(17, 1, m + 2) ! [gxy3 |s|s]
311 37809213 : s(28, 1, m + 1) = pai(1)*s(18, 1, m + 1) - pci(1)*s(18, 1, m + 2) ! [gxy2z|s|s]
312 37809213 : s(29, 1, m + 1) = pai(1)*s(19, 1, m + 1) - pci(1)*s(19, 1, m + 2) ! [gxyz2|s|s]
313 37809213 : s(30, 1, m + 1) = pai(1)*s(20, 1, m + 1) - pci(1)*s(20, 1, m + 2) ! [gxz3 |s|s]
314 : s(31, 1, m + 1) = pai(2)*s(17, 1, m + 1) - pci(2)*s(17, 1, m + 2) & ! [gy4 |s|s]
315 37809213 : + 3._dp*orho*(s(8, 1, m + 1) - s(8, 1, m + 2))
316 : s(32, 1, m + 1) = pai(2)*s(18, 1, m + 1) - pci(2)*s(18, 1, m + 2) & ! [gy3z |s|s]
317 37809213 : + 2._dp*orho*(s(9, 1, m + 1) - s(9, 1, m + 2))
318 : s(33, 1, m + 1) = pai(2)*s(19, 1, m + 1) - pci(2)*s(19, 1, m + 2) & ! [gy2z2|s|s]
319 37809213 : + orho*(s(10, 1, m + 1) - s(10, 1, m + 2))
320 37809213 : s(34, 1, m + 1) = pai(2)*s(20, 1, m + 1) - pci(2)*s(20, 1, m + 2) ! [gyz3 |s|s]
321 : s(35, 1, m + 1) = pai(3)*s(20, 1, m + 1) - pci(3)*s(20, 1, m + 2) & ! [gz4 |s|s]
322 67259332 : + 3._dp*orho*(s(10, 1, m + 1) - s(10, 1, m + 2))
323 : END DO
324 : ELSE
325 61060505 : DO irx = 0, llr
326 254433254 : DO iry = 0, llr - irx
327 193372749 : irz = llr - irx - iry
328 193372749 : irr(1) = irx; irr(2) = iry; irr(3) = irz
329 966863745 : ixx = MAXLOC(irr)
330 193372749 : ix = ixx(1)
331 193372749 : ir = coset(irx, iry, irz)
332 193372749 : irm = irr
333 193372749 : irm(ix) = irm(ix) - 1
334 193372749 : aai = REAL(MAX(irm(ix), 0), dp)*orho
335 193372749 : ir1 = coset(irm(1), irm(2), irm(3))
336 193372749 : irm(ix) = irm(ix) - 1
337 193372749 : ir2 = coset(irm(1), irm(2), irm(3))
338 492959997 : DO m = 0, mmax - llr
339 : s(ir, 1, m + 1) = pai(ix)*s(ir1, 1, m + 1) - pci(ix)*s(ir1, 1, m + 2) &
340 440258586 : + aai*(s(ir2, 1, m + 1) - s(ir2, 1, m + 2))
341 : END DO
342 : END DO
343 : END DO
344 : END IF
345 : END DO
346 :
347 : ! *** Horizontal recurrence steps ***
348 : ! *** [a|c|b+1i] = [a+1i|c|b] + (Ai - Bi)*[a|c|b] ***
349 :
350 198153929 : DO mb = 1, lb_max
351 474580115 : DO ibx = 0, mb
352 875931178 : DO iby = 0, mb - ibx
353 482494632 : ibz = mb - ibx - iby
354 482494632 : ib = coset(ibx, iby, ibz)
355 482494632 : ii(1) = ibx; ii(2) = iby; ii(3) = ibz
356 2412473160 : ixx = MAXLOC(ii)
357 482494632 : ix = ixx(1)
358 482494632 : abx = -rab(ix)
359 482494632 : iim = ii
360 482494632 : iim(ix) = iim(ix) - 1
361 482494632 : ibm = coset(iim(1), iim(2), iim(3))
362 6828312246 : DO ia = 1, ncoset(mmax - mb)
363 6069391428 : iap = iiap(ia, ix)
364 6551886060 : s(ia, ib, 1) = s(iap, ibm, 1) + abx*s(ia, ibm, 1)
365 : END DO
366 : END DO
367 : END DO
368 : END DO
369 :
370 8498008 : ELSE IF (lb_max > 0) THEN
371 :
372 : ! *** Recurrence steps: [s|c|s] -> [s|c|b] ***
373 : ! *** [s|c|b](m) = (Pi - Bi)*[s|c|b-1i](m) - ***
374 : ! *** (Pi - Ci)*[s|c|b-1i](m+1)) + ***
375 : ! *** Ni(b-1i)/2(a+b)*[s|c|b-2i](m) - ***
376 : ! *** Ni(b-1i)/2(a+b)*[s|c|b-2i](m+1) ***
377 4891065 : DO llr = 1, lb_max
378 4891065 : IF (llr == 1) THEN
379 4891065 : DO m = 0, lb_max - llr
380 2565375 : s1 = s(1, 1, m + 1)
381 2565375 : s2 = s(1, 1, m + 2)
382 2565375 : s(1, 2, m + 1) = pbi(1)*s1 - pci(1)*s2 ! [px|o|s]
383 2565375 : s(1, 3, m + 1) = pbi(2)*s1 - pci(2)*s2 ! [py|o|s]
384 4891065 : s(1, 4, m + 1) = pbi(3)*s1 - pci(3)*s2 ! [pz|o|s]
385 : END DO
386 239685 : ELSE IF (llr == 2) THEN
387 467728 : DO m = 0, lb_max - llr
388 239685 : s1 = s(1, 1, m + 1) - s(1, 1, m + 2)
389 239685 : s(1, 5, m + 1) = pbi(1)*s(1, 2, m + 1) - pci(1)*s(1, 2, m + 2) + orho*s1 ! [dx2|o|s]
390 239685 : s(1, 6, m + 1) = pbi(1)*s(1, 3, m + 1) - pci(1)*s(1, 3, m + 2) ! [dxy|o|s]
391 239685 : s(1, 7, m + 1) = pbi(1)*s(1, 4, m + 1) - pci(1)*s(1, 4, m + 2) ! [dxz|o|s]
392 239685 : s(1, 8, m + 1) = pbi(2)*s(1, 3, m + 1) - pci(2)*s(1, 3, m + 2) + orho*s1 ! [dy2|o|s]
393 239685 : s(1, 9, m + 1) = pbi(2)*s(1, 4, m + 1) - pci(2)*s(1, 4, m + 2) ! [dyz|o|s]
394 467728 : s(1, 10, m + 1) = pbi(3)*s(1, 4, m + 1) - pci(3)*s(1, 4, m + 2) + orho*s1 ! [dz2|o|s]
395 : END DO
396 11642 : ELSE IF (llr == 3) THEN
397 23202 : DO m = 0, lb_max - llr
398 : s(1, 11, m + 1) = pbi(1)*s(1, 5, m + 1) - pci(1)*s(1, 5, m + 2) & ! [fx3 |o|s]
399 11642 : + 2._dp*orho*(s(1, 2, m + 1) - s(1, 2, m + 2))
400 : s(1, 12, m + 1) = pbi(1)*s(1, 6, m + 1) - pci(1)*s(1, 6, m + 2) & ! [fx2y|o|s]
401 11642 : + orho*(s(1, 3, m + 1) - s(1, 3, m + 2))
402 : s(1, 13, m + 1) = pbi(1)*s(1, 7, m + 1) - pci(1)*s(1, 7, m + 2) & ! [fx2z|o|s]
403 11642 : + orho*(s(1, 4, m + 1) - s(1, 4, m + 2))
404 : s(1, 14, m + 1) = pbi(2)*s(1, 6, m + 1) - pci(2)*s(1, 6, m + 2) & ! [fxy2|o|s]
405 11642 : + orho*(s(1, 2, m + 1) - s(1, 2, m + 2))
406 11642 : s(1, 15, m + 1) = pbi(1)*s(1, 9, m + 1) - pci(1)*s(1, 9, m + 2) ! [fxyz|o|s]
407 : s(1, 16, m + 1) = pbi(3)*s(1, 7, m + 1) - pci(3)*s(1, 7, m + 2) & ! [fxz2|o|s]
408 11642 : + orho*(s(1, 2, m + 1) - s(1, 2, m + 2))
409 : s(1, 17, m + 1) = pbi(2)*s(1, 8, m + 1) - pci(2)*s(1, 8, m + 2) & ! [fy3 |o|s]
410 11642 : + 2._dp*orho*(s(1, 3, m + 1) - s(1, 3, m + 2))
411 : s(1, 18, m + 1) = pbi(2)*s(1, 9, m + 1) - pci(2)*s(1, 9, m + 2) & ! [fy2z|o|s]
412 11642 : + orho*(s(1, 4, m + 1) - s(1, 4, m + 2))
413 : s(1, 19, m + 1) = pbi(3)*s(1, 9, m + 1) - pci(3)*s(1, 9, m + 2) & ! [fyz2|o|s]
414 11642 : + orho*(s(1, 3, m + 1) - s(1, 3, m + 2))
415 : s(1, 20, m + 1) = pbi(3)*s(1, 10, m + 1) - pci(3)*s(1, 10, m + 2) & ! [fz3 |o|s]
416 23202 : + 2._dp*orho*(s(1, 4, m + 1) - s(1, 4, m + 2))
417 : END DO
418 : ELSE
419 492 : DO irx = 0, llr
420 1722 : DO iry = 0, llr - irx
421 1230 : irz = llr - irx - iry
422 1230 : irr(1) = irx; irr(2) = iry; irr(3) = irz
423 6150 : ixx = MAXLOC(irr)
424 1230 : ix = ixx(1)
425 1230 : ir = coset(irx, iry, irz)
426 1230 : irm = irr
427 1230 : irm(ix) = irm(ix) - 1
428 1230 : aai = REAL(MAX(irm(ix), 0), dp)
429 1230 : ir1 = coset(irm(1), irm(2), irm(3))
430 1230 : irm(ix) = irm(ix) - 1
431 1230 : ir2 = coset(irm(1), irm(2), irm(3))
432 2870 : DO m = 0, lb_max - llr
433 : s(1, ir, m + 1) = pbi(ix)*s(1, ir1, m + 1) - pci(ix)*s(1, ir1, m + 2) &
434 2460 : + aai*orho*(s(1, ir2, m + 1) - s(1, ir2, m + 2))
435 : END DO
436 : END DO
437 : END DO
438 : END IF
439 : END DO
440 :
441 : END IF
442 :
443 : ! *** Store the primitive three-center overlap integrals ***
444 475339111 : DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
445 2468789004 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
446 2379147427 : vab(na + i, nb + j) = vab(na + i, nb + j) + s(i, j, 1)
447 : END DO
448 : END DO
449 :
450 : ! *** Calculate the requested derivatives with respect ***
451 : ! *** to the nuclear coordinates of the atomic center a ***
452 :
453 89641577 : IF (da_max == 1) THEN
454 103916784 : DO idir = 1, 3
455 : CALL cartesian_ladder(s(:, b_first:b_last, 1), fs(a_first:a_last, b_first:b_last, idir + 1), &
456 103916784 : map_a(:, :, idir), .TRUE., zeta(ipgf), 0.0_dp, .FALSE.)
457 : END DO
458 : END IF
459 :
460 : ! DFPT for APTs
461 89641577 : IF (PRESENT(vab2_work)) THEN
462 28512 : DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
463 72090 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
464 43578 : vab2_work(na + i, nb + j, 1) = vab2_work(na + i, nb + j, 1) + fs(i, j, 2)
465 43578 : vab2_work(na + i, nb + j, 2) = vab2_work(na + i, nb + j, 2) + fs(i, j, 3)
466 62676 : vab2_work(na + i, nb + j, 3) = vab2_work(na + i, nb + j, 3) + fs(i, j, 4)
467 : END DO
468 : END DO
469 : END IF
470 :
471 : ! *** Calculate the force contribution for the atomic center a ***
472 :
473 89641577 : IF (calculate_force_a) THEN
474 142242561 : DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
475 749372952 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
476 607130391 : force_a(1) = force_a(1) + pab(na + i, nb + j)*fs(i, j, 2)
477 607130391 : force_a(2) = force_a(2) + pab(na + i, nb + j)*fs(i, j, 3)
478 723403170 : force_a(3) = force_a(3) + pab(na + i, nb + j)*fs(i, j, 4)
479 : END DO
480 : END DO
481 : END IF
482 :
483 : ! *** Calculate the requested derivatives with respect ***
484 : ! *** to the nuclear coordinates of the atomic center b ***
485 :
486 89641577 : IF (db_max == 1) THEN
487 103916784 : DO idir = 1, 3
488 : CALL cartesian_ladder(s(a_first:a_last, :, 1), fs(a_first:a_last, b_first:b_last, idir + 1), &
489 103916784 : map_b(:, :, idir), .TRUE., zetb(jpgf), 0.0_dp, .TRUE.)
490 : END DO
491 : END IF
492 :
493 : ! DFPT for APTs
494 89641577 : IF (PRESENT(vab2_work)) THEN
495 28512 : DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
496 72090 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
497 43578 : vab2_work(na + i, nb + j, 4) = vab2_work(na + i, nb + j, 4) + fs(i, j, 2)
498 43578 : vab2_work(na + i, nb + j, 5) = vab2_work(na + i, nb + j, 5) + fs(i, j, 3)
499 62676 : vab2_work(na + i, nb + j, 6) = vab2_work(na + i, nb + j, 6) + fs(i, j, 4)
500 : END DO
501 : END DO
502 :
503 9414 : CPASSERT(PRESENT(iatom) .AND. PRESENT(jatom) .AND. PRESENT(katom))
504 9414 : CPASSERT(PRESENT(deltaR))
505 37656 : DO idir = 1, 3
506 94950 : DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
507 216270 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
508 : vab2(na + i, nb + j, idir) = vab2(na + i, nb + j, idir) &
509 : + vab2_work(na + i, nb + j, idir)*deltaR(idir, iatom) &
510 : - vab2_work(na + i, nb + j, idir)*deltaR(idir, katom) &
511 : + vab2_work(na + i, nb + j, idir + 3)*deltaR(idir, jatom) &
512 188028 : - vab2_work(na + i, nb + j, idir + 3)*deltaR(idir, katom)
513 : END DO
514 : END DO
515 : END DO
516 : END IF
517 :
518 : ! *** Calculate the force contribution for the atomic center b ***
519 :
520 89641577 : IF (calculate_force_b) THEN
521 142242561 : DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
522 749372952 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
523 607130391 : force_b(1) = force_b(1) + pab(na + i, nb + j)*fs(i, j, 2)
524 607130391 : force_b(2) = force_b(2) + pab(na + i, nb + j)*fs(i, j, 3)
525 723403170 : force_b(3) = force_b(3) + pab(na + i, nb + j)*fs(i, j, 4)
526 : END DO
527 : END DO
528 : END IF
529 :
530 157322714 : nb = nb + ncoset(lb_max_set)
531 :
532 : END DO
533 :
534 119623046 : na = na + ncoset(la_max_set)
535 :
536 : END DO
537 :
538 51941909 : DEALLOCATE (iiap)
539 :
540 51941909 : END SUBROUTINE os_3center
541 : ! **************************************************************************************************
542 : !> \brief Calculation of two-center integrals <a|c> over
543 : !> Cartesian Gaussian functions and a spherical potential
544 : !>
545 : !> \param la_max_set ...
546 : !> \param la_min_set ...
547 : !> \param npgfa ...
548 : !> \param rpgfa ...
549 : !> \param zeta ...
550 : !> \param auxint ...
551 : !> \param rpgfc ...
552 : !> \param rac ...
553 : !> \param dac ...
554 : !> \param va ...
555 : !> \param dva ...
556 : !> \date December 2017
557 : !> \author Juerg Hutter
558 : !> \version 1.0
559 : ! **************************************************************************************************
560 360 : SUBROUTINE os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
561 720 : auxint, rpgfc, rac, dac, va, dva)
562 : INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
563 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
564 : REAL(KIND=dp), DIMENSION(0:, :), INTENT(IN) :: auxint
565 : REAL(KIND=dp), INTENT(IN) :: rpgfc
566 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
567 : REAL(KIND=dp), INTENT(IN) :: dac
568 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: va
569 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
570 : OPTIONAL :: dva
571 :
572 : INTEGER :: a_first, a_last, da_max, i, idir, ipgf, &
573 : ir, ir1, ir2, irm(3), irr(3), irx, &
574 : iry, irz, ix, ixx(1), la_max, llr, m, &
575 : mmax, na
576 360 : INTEGER, ALLOCATABLE :: map_a(:, :, :)
577 : REAL(KIND=dp) :: aai, orho, s1
578 360 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: derivative, s
579 :
580 360 : IF (PRESENT(dva)) THEN
581 : da_max = 1
582 : ELSE
583 164 : da_max = 0
584 : END IF
585 :
586 360 : la_max = la_max_set + da_max
587 :
588 360 : mmax = la_max
589 :
590 1440 : ALLOCATE (s(ncoset(mmax), mmax + 1))
591 360 : a_first = ncoset(la_min_set - 1) + 1
592 360 : a_last = ncoset(la_max_set)
593 360 : IF (PRESENT(dva)) THEN
594 1176 : ALLOCATE (map_a(4, a_last - a_first + 1, 3), derivative(a_last - a_first + 1, 1))
595 784 : DO idir = 1, 3
596 784 : CALL cartesian_ladder_map(la_max_set, la_min_set, idir, map_a(:, :, idir), source_lmin=0)
597 : END DO
598 : END IF
599 360 : na = 0
600 752 : DO ipgf = 1, npgfa
601 392 : IF (rpgfa(ipgf) + rpgfc < dac) THEN
602 0 : na = na + ncoset(la_max_set)
603 0 : CYCLE
604 : END IF
605 1708 : s(1, 1:mmax + 1) = auxint(0:mmax, ipgf)
606 392 : IF (la_max > 0) THEN
607 : ! Recurrence steps: [s|c] -> [a|c]
608 : ! [a|c](m) = (Ci - Ai)*[a-1i|c](m+1) +
609 : ! Ni(a-1i)/2a*[a-2i|c](m) -
610 : ! Ni(a-1i)/2a*[a-2i|c](m+1)
611 : !
612 :
613 332 : orho = 0.5_dp/zeta(ipgf)
614 :
615 1256 : DO llr = 1, mmax
616 1256 : IF (llr == 1) THEN
617 1256 : DO m = 0, mmax - llr
618 924 : s1 = s(1, m + 2)
619 924 : s(2, m + 1) = -rac(1)*s1 ! [px|o]
620 924 : s(3, m + 1) = -rac(2)*s1 ! [py|o]
621 1256 : s(4, m + 1) = -rac(3)*s1 ! [pz|o]
622 : END DO
623 592 : ELSE IF (llr == 2) THEN
624 864 : DO m = 0, mmax - llr
625 592 : s1 = s(1, m + 1) - s(1, m + 2)
626 592 : s(5, m + 1) = -rac(1)*s(2, m + 2) + orho*s1 ! [dx2|o]
627 592 : s(6, m + 1) = -rac(1)*s(3, m + 2) ! [dxy|o]
628 592 : s(7, m + 1) = -rac(1)*s(4, m + 2) ! [dxz|o]
629 592 : s(8, m + 1) = -rac(2)*s(3, m + 2) + orho*s1 ! [dy2|o]
630 592 : s(9, m + 1) = -rac(2)*s(4, m + 2) ! [dyz|o]
631 864 : s(10, m + 1) = -rac(3)*s(4, m + 2) + orho*s1 ! [dz2|o]
632 : END DO
633 320 : ELSE IF (llr == 3) THEN
634 500 : DO m = 0, mmax - llr
635 320 : s(11, m + 1) = -rac(1)*s(5, m + 2) + 2._dp*orho*(s(2, m + 1) - s(2, m + 2)) ! [fx3 |o]
636 320 : s(12, m + 1) = -rac(1)*s(6, m + 2) + orho*(s(3, m + 1) - s(3, m + 2)) ! [fx2y|o]
637 320 : s(13, m + 1) = -rac(1)*s(7, m + 2) + orho*(s(4, m + 1) - s(4, m + 2)) ! [fx2z|o]
638 320 : s(14, m + 1) = -rac(2)*s(6, m + 2) + orho*(s(2, m + 1) - s(2, m + 2)) ! [fxy2|o]
639 320 : s(15, m + 1) = -rac(1)*s(9, m + 2) ! [fxyz|o]
640 320 : s(16, m + 1) = -rac(3)*s(7, m + 2) + orho*(s(2, m + 1) - s(2, m + 2)) ! [fxz2|o]
641 320 : s(17, m + 1) = -rac(2)*s(8, m + 2) + 2._dp*orho*(s(3, m + 1) - s(3, m + 2)) ! [fy3 |o]
642 320 : s(18, m + 1) = -rac(2)*s(9, m + 2) + orho*(s(4, m + 1) - s(4, m + 2)) ! [fy2z|o]
643 320 : s(19, m + 1) = -rac(3)*s(9, m + 2) + orho*(s(3, m + 1) - s(3, m + 2)) ! [fyz2|o]
644 500 : s(20, m + 1) = -rac(3)*s(10, m + 2) + 2._dp*orho*(s(4, m + 1) - s(4, m + 2)) ! [fz3 |o]
645 : END DO
646 140 : ELSE IF (llr == 4) THEN
647 216 : DO m = 0, mmax - llr
648 140 : s(21, m + 1) = -rac(1)*s(11, m + 2) + 3._dp*orho*(s(5, m + 1) - s(5, m + 2)) ! [gx4 |s]
649 140 : s(22, m + 1) = -rac(1)*s(12, m + 2) + 2._dp*orho*(s(6, m + 1) - s(6, m + 2)) ! [gx3y |s]
650 140 : s(23, m + 1) = -rac(1)*s(13, m + 2) + 2._dp*orho*(s(7, m + 1) - s(7, m + 2)) ! [gx3z |s]
651 140 : s(24, m + 1) = -rac(1)*s(14, m + 2) + orho*(s(8, m + 1) - s(8, m + 2)) ! [gx2y2|s]
652 140 : s(25, m + 1) = -rac(1)*s(15, m + 2) + orho*(s(9, m + 1) - s(9, m + 2)) ! [gx2yz|s]
653 140 : s(26, m + 1) = -rac(1)*s(16, m + 2) + orho*(s(10, m + 1) - s(10, m + 2)) ! [gx2z2|s]
654 140 : s(27, m + 1) = -rac(1)*s(17, m + 2) ! [gxy3 |s]
655 140 : s(28, m + 1) = -rac(1)*s(18, m + 2) ! [gxy2z|s]
656 140 : s(29, m + 1) = -rac(1)*s(19, m + 2) ! [gxyz2|s]
657 140 : s(30, m + 1) = -rac(1)*s(20, m + 2) ! [gxz3 |s]
658 140 : s(31, m + 1) = -rac(2)*s(17, m + 2) + 3._dp*orho*(s(8, m + 1) - s(8, m + 2)) ! [gy4 |s]
659 140 : s(32, m + 1) = -rac(2)*s(18, m + 2) + 2._dp*orho*(s(9, m + 1) - s(9, m + 2)) ! [gy3z |s]
660 140 : s(33, m + 1) = -rac(2)*s(19, m + 2) + orho*(s(10, m + 1) - s(10, m + 2)) ! [gy2z2|s]
661 140 : s(34, m + 1) = -rac(2)*s(20, m + 2) ! [gyz3 |s]
662 216 : s(35, m + 1) = -rac(3)*s(20, m + 2) + 3._dp*orho*(s(10, m + 1) - s(10, m + 2)) ! [gz4 |s]
663 : END DO
664 : ELSE
665 448 : DO irx = 0, llr
666 1792 : DO iry = 0, llr - irx
667 1344 : irz = llr - irx - iry
668 1344 : irr(1) = irx; irr(2) = iry; irr(3) = irz
669 6720 : ixx = MAXLOC(irr)
670 1344 : ix = ixx(1)
671 1344 : ir = coset(irx, iry, irz)
672 1344 : irm = irr
673 1344 : irm(ix) = irm(ix) - 1
674 1344 : aai = REAL(MAX(irm(ix), 0), dp)*orho
675 1344 : ir1 = coset(irm(1), irm(2), irm(3))
676 1344 : irm(ix) = irm(ix) - 1
677 1344 : ir2 = coset(irm(1), irm(2), irm(3))
678 3072 : DO m = 0, mmax - llr
679 2688 : s(ir, m + 1) = -rac(ix)*s(ir1, m + 2) + aai*(s(ir2, m + 1) - s(ir2, m + 2))
680 : END DO
681 : END DO
682 : END DO
683 : END IF
684 : END DO
685 :
686 : END IF
687 :
688 : ! Store the primitive three-center overlap integrals
689 4708 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
690 4708 : va(na + i) = va(na + i) + s(i, 1)
691 : END DO
692 :
693 : ! Calculate the requested derivatives with respect ***
694 : ! to the nuclear coordinates of the atomic center a ***
695 : ! [da/dAi|c] = 2*zeta*[a+1i|c] - Ni(a)[a-1i|c] ***
696 392 : IF (PRESENT(dva)) THEN
697 912 : DO idir = 1, 3
698 684 : CALL cartesian_ladder(s(:, 1:1), derivative, map_a(:, :, idir), .TRUE., zeta(ipgf), 0.0_dp, .FALSE.)
699 10746 : dva(na + a_first:na + a_last, idir) = dva(na + a_first:na + a_last, idir) + derivative(:, 1)
700 : END DO
701 : END IF
702 :
703 752 : na = na + ncoset(la_max_set)
704 :
705 : END DO
706 :
707 360 : DEALLOCATE (s)
708 :
709 360 : END SUBROUTINE os_2center
710 : ! **************************************************************************************************
711 :
712 : END MODULE ai_oneelectron
|
