Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Integral GKS scheme: The order of the integrals in makeCoul reflects
10 : !> the standard order by MOPAC
11 : !> \par History
12 : !> Teodoro Laino [tlaino] - 04.2009 : Adapted size arrays to d-orbitals and
13 : !> get rid of the alternative ordering. Using the
14 : !> CP2K one.
15 : !> Teodoro Laino [tlaino] - 04.2009 : Skip nullification (speed-up)
16 : !> Teodoro Laino [tlaino] - 04.2009 : Speed-up due to fortran arrays order
17 : !> optimization and collection of common pieces of
18 : !> code
19 : ! **************************************************************************************************
20 : MODULE semi_empirical_int_gks
21 :
22 : USE dg_rho0_types, ONLY: dg_rho0_type
23 : USE dg_types, ONLY: dg_get,&
24 : dg_type
25 : USE kinds, ONLY: dp
26 : USE mathconstants, ONLY: fourpi,&
27 : oorootpi
28 : USE multipole_types, ONLY: do_multipole_none
29 : USE pw_grid_types, ONLY: pw_grid_type
30 : USE pw_pool_types, ONLY: pw_pool_type
31 : USE semi_empirical_int_arrays, ONLY: indexb,&
32 : rij_threshold
33 : USE semi_empirical_mpole_types, ONLY: semi_empirical_mpole_type
34 : USE semi_empirical_types, ONLY: se_int_control_type,&
35 : semi_empirical_type,&
36 : setup_se_int_control_type
37 : #include "./base/base_uses.f90"
38 :
39 : IMPLICIT NONE
40 : PRIVATE
41 :
42 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_int_gks'
43 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
44 :
45 : PUBLIC :: corecore_gks, rotnuc_gks, drotnuc_gks, rotint_gks, drotint_gks
46 :
47 : CONTAINS
48 :
49 : ! **************************************************************************************************
50 : !> \brief Computes the electron-nuclei integrals
51 : !> \param sepi ...
52 : !> \param sepj ...
53 : !> \param rij ...
54 : !> \param e1b ...
55 : !> \param e2a ...
56 : !> \param se_int_control ...
57 : ! **************************************************************************************************
58 3813 : SUBROUTINE rotnuc_gks(sepi, sepj, rij, e1b, e2a, se_int_control)
59 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
60 : REAL(dp), DIMENSION(3), INTENT(IN) :: rij
61 : REAL(dp), DIMENSION(45), INTENT(OUT), OPTIONAL :: e1b, e2a
62 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
63 :
64 : INTEGER :: i, mu, nu
65 : REAL(KIND=dp), DIMENSION(3) :: rab
66 : REAL(kind=dp), DIMENSION(45, 45) :: Coul
67 :
68 15252 : rab = -rij
69 :
70 3813 : IF (se_int_control%do_ewald_gks) THEN
71 24 : IF (DOT_PRODUCT(rij, rij) > rij_threshold) THEN
72 3 : CALL makeCoulE(rab, sepi, sepj, Coul, se_int_control)
73 : ELSE
74 3 : CALL makeCoulE0(sepi, Coul, se_int_control)
75 : END IF
76 : ELSE
77 3807 : CALL makeCoul(rab, sepi, sepj, Coul, se_int_control)
78 : END IF
79 :
80 3813 : i = 0
81 19065 : DO mu = 1, sepi%natorb
82 57195 : DO nu = 1, mu
83 38130 : i = i + 1
84 53382 : e1b(i) = -Coul(i, 1)*sepj%zeff
85 : END DO
86 : END DO
87 :
88 3813 : i = 0
89 19065 : DO mu = 1, sepj%natorb
90 57195 : DO nu = 1, mu
91 38130 : i = i + 1
92 53382 : e2a(i) = -Coul(1, i)*sepi%zeff
93 : END DO
94 : END DO
95 :
96 3813 : END SUBROUTINE rotnuc_gks
97 :
98 : ! **************************************************************************************************
99 : !> \brief Computes the electron-electron integrals
100 : !> \param sepi ...
101 : !> \param sepj ...
102 : !> \param rij ...
103 : !> \param w ...
104 : !> \param se_int_control ...
105 : ! **************************************************************************************************
106 4754 : SUBROUTINE rotint_gks(sepi, sepj, rij, w, se_int_control)
107 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
108 : REAL(dp), DIMENSION(3), INTENT(IN) :: rij
109 : REAL(dp), DIMENSION(2025), INTENT(OUT), OPTIONAL :: w
110 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
111 :
112 : INTEGER :: i, ind1, ind2, lam, mu, nu, sig
113 : REAL(KIND=dp), DIMENSION(3) :: rab
114 : REAL(kind=dp), DIMENSION(45, 45) :: Coul
115 :
116 19016 : rab = -rij
117 :
118 4754 : IF (se_int_control%do_ewald_gks) THEN
119 24 : IF (DOT_PRODUCT(rij, rij) > rij_threshold) THEN
120 3 : CALL makeCoulE(rab, sepi, sepj, Coul, se_int_control)
121 : ELSE
122 3 : CALL makeCoulE0(sepi, Coul, se_int_control)
123 : END IF
124 : ELSE
125 4748 : CALL makeCoul(rab, sepi, sepj, Coul, se_int_control)
126 : END IF
127 :
128 4754 : i = 0
129 4754 : ind1 = 0
130 23770 : DO mu = 1, sepi%natorb
131 71310 : DO nu = 1, mu
132 47540 : ind1 = ind1 + 1
133 47540 : ind2 = 0
134 256716 : DO lam = 1, sepj%natorb
135 713100 : DO sig = 1, lam
136 475400 : i = i + 1
137 475400 : ind2 = ind2 + 1
138 665560 : w(i) = Coul(ind1, ind2)
139 : END DO
140 : END DO
141 : END DO
142 : END DO
143 :
144 4754 : END SUBROUTINE rotint_gks
145 :
146 : ! **************************************************************************************************
147 : !> \brief Computes the derivatives of the electron-nuclei integrals
148 : !> \param sepi ...
149 : !> \param sepj ...
150 : !> \param rij ...
151 : !> \param de1b ...
152 : !> \param de2a ...
153 : !> \param se_int_control ...
154 : ! **************************************************************************************************
155 0 : SUBROUTINE drotnuc_gks(sepi, sepj, rij, de1b, de2a, se_int_control)
156 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
157 : REAL(dp), DIMENSION(3), INTENT(IN) :: rij
158 : REAL(dp), DIMENSION(3, 45), INTENT(OUT), OPTIONAL :: de1b, de2a
159 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
160 :
161 : INTEGER :: i, mu, nu
162 : REAL(KIND=dp), DIMENSION(3) :: rab
163 : REAL(kind=dp), DIMENSION(3, 45, 45) :: dCoul
164 :
165 0 : rab = -rij
166 :
167 0 : IF (se_int_control%do_ewald_gks) THEN
168 0 : CALL makedCoulE(rab, sepi, sepj, dCoul, se_int_control)
169 : ELSE
170 0 : CALL makedCoul(rab, sepi, sepj, dCoul, se_int_control)
171 : END IF
172 :
173 0 : i = 0
174 0 : DO mu = 1, sepi%natorb
175 0 : DO nu = 1, mu
176 0 : i = i + 1
177 0 : de1b(1, i) = dCoul(1, i, 1)*sepj%zeff
178 0 : de1b(2, i) = dCoul(2, i, 1)*sepj%zeff
179 0 : de1b(3, i) = dCoul(3, i, 1)*sepj%zeff
180 : END DO
181 : END DO
182 :
183 0 : i = 0
184 0 : DO mu = 1, sepj%natorb
185 0 : DO nu = 1, mu
186 0 : i = i + 1
187 0 : de2a(1, i) = dCoul(1, 1, i)*sepi%zeff
188 0 : de2a(2, i) = dCoul(2, 1, i)*sepi%zeff
189 0 : de2a(3, i) = dCoul(3, 1, i)*sepi%zeff
190 : END DO
191 : END DO
192 :
193 0 : END SUBROUTINE drotnuc_gks
194 :
195 : ! **************************************************************************************************
196 : !> \brief Computes the derivatives of the electron-electron integrals
197 : !> \param sepi ...
198 : !> \param sepj ...
199 : !> \param rij ...
200 : !> \param dw ...
201 : !> \param se_int_control ...
202 : ! **************************************************************************************************
203 0 : SUBROUTINE drotint_gks(sepi, sepj, rij, dw, se_int_control)
204 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
205 : REAL(dp), DIMENSION(3), INTENT(IN) :: rij
206 : REAL(dp), DIMENSION(3, 2025), INTENT(OUT) :: dw
207 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
208 :
209 : INTEGER :: i, ind1, ind2, lam, mu, nu, sig
210 : REAL(KIND=dp), DIMENSION(3) :: rab
211 : REAL(kind=dp), DIMENSION(3, 45, 45) :: dCoul
212 :
213 0 : rab = -rij
214 :
215 0 : IF (se_int_control%do_ewald_gks) THEN
216 0 : CALL makedCoulE(rab, sepi, sepj, dCoul, se_int_control)
217 : ELSE
218 0 : CALL makedCoul(rab, sepi, sepj, dCoul, se_int_control)
219 : END IF
220 :
221 0 : i = 0
222 0 : ind1 = 0
223 0 : DO mu = 1, sepi%natorb
224 0 : DO nu = 1, mu
225 0 : ind1 = ind1 + 1
226 0 : ind2 = 0
227 0 : DO lam = 1, sepj%natorb
228 0 : DO sig = 1, lam
229 0 : i = i + 1
230 0 : ind2 = ind2 + 1
231 0 : dw(1, i) = -dCoul(1, ind1, ind2)
232 0 : dw(2, i) = -dCoul(2, ind1, ind2)
233 0 : dw(3, i) = -dCoul(3, ind1, ind2)
234 : END DO
235 : END DO
236 : END DO
237 : END DO
238 :
239 0 : END SUBROUTINE drotint_gks
240 :
241 : ! **************************************************************************************************
242 : !> \brief Computes the primitives of the integrals (non-periodic case)
243 : !> \param RAB ...
244 : !> \param sepi ...
245 : !> \param sepj ...
246 : !> \param Coul ...
247 : !> \param se_int_control ...
248 : ! **************************************************************************************************
249 32338 : SUBROUTINE makeCoul(RAB, sepi, sepj, Coul, se_int_control)
250 : REAL(kind=dp), DIMENSION(3) :: RAB
251 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
252 : REAL(kind=dp), DIMENSION(45, 45), INTENT(OUT) :: Coul
253 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
254 :
255 : INTEGER :: iA, iB, imA, imB, jA, jB, k1, k2, k3, k4
256 : LOGICAL :: shortrange
257 : REAL(kind=dp) :: a2, ACOULA, ACOULB, d1, d1f(3), d2, &
258 : d2f(3, 3), d3, d3f(3, 3, 3), d4, &
259 : d4f(3, 3, 3, 3), f, rr, w, w0, w1, w2, &
260 : w3, w4, w5
261 : REAL(kind=dp), DIMENSION(3) :: v
262 : REAL(kind=dp), DIMENSION(3, 3, 45) :: M2A, M2B
263 : REAL(kind=dp), DIMENSION(3, 45) :: M1A, M1B
264 : REAL(kind=dp), DIMENSION(45) :: M0A, M0B
265 :
266 16169 : shortrange = se_int_control%shortrange
267 16169 : CALL get_se_slater_multipole(sepi, M0A, M1A, M2A, ACOULA)
268 16169 : CALL get_se_slater_multipole(sepj, M0B, M1B, M2B, ACOULB)
269 :
270 16169 : v(1) = RAB(1)
271 16169 : v(2) = RAB(2)
272 16169 : v(3) = RAB(3)
273 64676 : rr = SQRT(DOT_PRODUCT(v, v))
274 :
275 16169 : a2 = 0.5_dp*(1.0_dp/ACOULA + 1.0_dp/ACOULB)
276 16169 : w0 = a2*rr
277 16169 : w = EXP(-w0)
278 16169 : w1 = (1.0_dp + 0.5_dp*w0)
279 16169 : w2 = (w1 + 0.5_dp*w0 + 0.5_dp*w0**2)
280 16169 : w3 = (w2 + w0**3/6.0_dp)
281 16169 : w4 = (w3 + w0**4/30.0_dp)
282 16169 : w5 = (w3 + 8.0_dp*w0**4/210.0_dp + w0**5/210.0_dp)
283 :
284 16169 : IF (shortrange) THEN
285 11412 : f = (-w*w1)/rr
286 11412 : d1 = -1.0_dp*(-w*w2)/rr**3
287 11412 : d2 = 3.0_dp*(-w*w3)/rr**5
288 11412 : d3 = -15.0_dp*(-w*w4)/rr**7
289 11412 : d4 = 105.0_dp*(-w*w5)/rr**9
290 : ELSE
291 4757 : f = (1.0_dp - w*w1)/rr
292 4757 : d1 = -1.0_dp*(1.0_dp - w*w2)/rr**3
293 4757 : d2 = 3.0_dp*(1.0_dp - w*w3)/rr**5
294 4757 : d3 = -15.0_dp*(1.0_dp - w*w4)/rr**7
295 4757 : d4 = 105.0_dp*(1.0_dp - w*w5)/rr**9
296 : END IF
297 :
298 16169 : CALL build_d_tensor_gks(d1f, d2f, d3f, d4f, v=v, d1=d1, d2=d2, d3=d3, d4=d4)
299 :
300 16169 : imA = 0
301 80845 : DO iA = 1, sepi%natorb
302 242535 : DO jA = 1, iA
303 161690 : imA = imA + 1
304 :
305 161690 : imB = 0
306 873126 : DO iB = 1, sepj%natorb
307 2425350 : DO jB = 1, iB
308 1616900 : imB = imB + 1
309 :
310 1616900 : w = M0A(imA)*M0B(imB)*f
311 6467600 : DO k1 = 1, 3
312 6467600 : w = w + (M1A(k1, imA)*M0B(imB) - M0A(imA)*M1B(k1, imB))*d1f(k1)
313 : END DO
314 6467600 : DO k2 = 1, 3
315 21019700 : DO k1 = 1, 3
316 19402800 : w = w + (M2A(k1, k2, imA)*M0B(imB) - M1A(k1, imA)*M1B(k2, imB) + M0A(imA)*M2B(k1, k2, imB))*d2f(k1, k2)
317 : END DO
318 : END DO
319 6467600 : DO k3 = 1, 3
320 21019700 : DO k2 = 1, 3
321 63059100 : DO k1 = 1, 3
322 58208400 : w = w + (-M2A(k1, k2, imA)*M1B(k3, imB) + M1A(k1, imA)*M2B(k2, k3, imB))*d3f(k1, k2, k3)
323 : END DO
324 : END DO
325 : END DO
326 :
327 6467600 : DO k4 = 1, 3
328 21019700 : DO k3 = 1, 3
329 63059100 : DO k2 = 1, 3
330 189177300 : DO k1 = 1, 3
331 174625200 : w = w + M2A(k1, k2, imA)*M2B(k3, k4, imB)*d4f(k1, k2, k3, k4)
332 : END DO
333 : END DO
334 : END DO
335 : END DO
336 :
337 2263660 : Coul(imA, imB) = w
338 : END DO
339 : END DO
340 : END DO
341 : END DO
342 :
343 16169 : END SUBROUTINE makeCoul
344 :
345 : ! **************************************************************************************************
346 : !> \brief Computes the derivatives of the primitives of the integrals
347 : !> (non-periodic case)
348 : !> \param RAB ...
349 : !> \param sepi ...
350 : !> \param sepj ...
351 : !> \param dCoul ...
352 : !> \param se_int_control ...
353 : ! **************************************************************************************************
354 0 : SUBROUTINE makedCoul(RAB, sepi, sepj, dCoul, se_int_control)
355 : REAL(kind=dp), DIMENSION(3) :: RAB
356 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
357 : REAL(kind=dp), DIMENSION(3, 45, 45), INTENT(OUT) :: dCoul
358 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
359 :
360 : INTEGER :: iA, iB, imA, imB, jA, jB, k1, k2, k3, k4
361 : LOGICAL :: shortrange
362 : REAL(kind=dp) :: a2, ACOULA, ACOULB, d1, d1f(3), d2, d2f(3, 3), d3, d3f(3, 3, 3), d4, &
363 : d4f(3, 3, 3, 3), d5, d5f(3, 3, 3, 3, 3), f, rr, tmp, w, w0, w1, w2, w3, w4, w5, w6
364 : REAL(kind=dp), DIMENSION(3) :: v, wv
365 : REAL(kind=dp), DIMENSION(3, 3, 45) :: M2A, M2B
366 : REAL(kind=dp), DIMENSION(3, 45) :: M1A, M1B
367 : REAL(kind=dp), DIMENSION(45) :: M0A, M0B
368 :
369 0 : shortrange = se_int_control%shortrange
370 0 : CALL get_se_slater_multipole(sepi, M0A, M1A, M2A, ACOULA)
371 0 : CALL get_se_slater_multipole(sepj, M0B, M1B, M2B, ACOULB)
372 :
373 0 : v(1) = RAB(1)
374 0 : v(2) = RAB(2)
375 0 : v(3) = RAB(3)
376 0 : rr = SQRT(DOT_PRODUCT(v, v))
377 :
378 0 : a2 = 0.5_dp*(1.0_dp/ACOULA + 1.0_dp/ACOULB)
379 0 : w0 = a2*rr
380 0 : w = EXP(-w0)
381 0 : w1 = (1.0_dp + 0.5_dp*w0)
382 0 : w2 = (w1 + 0.5_dp*w0 + 0.5_dp*w0**2)
383 0 : w3 = (w2 + w0**3/6.0_dp)
384 0 : w4 = (w3 + w0**4/30.0_dp)
385 0 : w5 = (w3 + 4.0_dp*w0**4/105.0_dp + w0**5/210.0_dp)
386 0 : w6 = (w3 + 15.0_dp*w0**4/378.0_dp + 2.0_dp*w0**5/315.0_dp + w0**6/1890.0_dp)
387 :
388 0 : IF (shortrange) THEN
389 0 : f = (-w*w1)/rr
390 0 : d1 = -1.0_dp*(-w*w2)/rr**3
391 0 : d2 = 3.0_dp*(-w*w3)/rr**5
392 0 : d3 = -15.0_dp*(-w*w4)/rr**7
393 0 : d4 = 105.0_dp*(-w*w5)/rr**9
394 0 : d5 = -945.0_dp*(-w*w6)/rr**11
395 : ELSE
396 0 : f = (1.0_dp - w*w1)/rr
397 0 : d1 = -1.0_dp*(1.0_dp - w*w2)/rr**3
398 0 : d2 = 3.0_dp*(1.0_dp - w*w3)/rr**5
399 0 : d3 = -15.0_dp*(1.0_dp - w*w4)/rr**7
400 0 : d4 = 105.0_dp*(1.0_dp - w*w5)/rr**9
401 0 : d5 = -945.0_dp*(1.0_dp - w*w6)/rr**11
402 : END IF
403 :
404 0 : CALL build_d_tensor_gks(d1f, d2f, d3f, d4f, d5f, v, d1, d2, d3, d4, d5)
405 :
406 0 : imA = 0
407 0 : DO iA = 1, sepi%natorb
408 0 : DO jA = 1, iA
409 0 : imA = imA + 1
410 :
411 0 : imB = 0
412 0 : DO iB = 1, sepj%natorb
413 0 : DO jB = 1, iB
414 0 : imB = imB + 1
415 :
416 0 : tmp = M0A(imA)*M0B(imB)
417 0 : wv(1) = tmp*d1f(1)
418 0 : wv(2) = tmp*d1f(2)
419 0 : wv(3) = tmp*d1f(3)
420 0 : DO k1 = 1, 3
421 0 : tmp = M1A(k1, imA)*M0B(imB) - M0A(imA)*M1B(k1, imB)
422 0 : wv(1) = wv(1) + tmp*d2f(1, k1)
423 0 : wv(2) = wv(2) + tmp*d2f(2, k1)
424 0 : wv(3) = wv(3) + tmp*d2f(3, k1)
425 : END DO
426 0 : DO k2 = 1, 3
427 0 : DO k1 = 1, 3
428 0 : tmp = M2A(k1, k2, imA)*M0B(imB) - M1A(k1, imA)*M1B(k2, imB) + M0A(imA)*M2B(k1, k2, imB)
429 0 : wv(1) = wv(1) + tmp*d3f(1, k1, k2)
430 0 : wv(2) = wv(2) + tmp*d3f(2, k1, k2)
431 0 : wv(3) = wv(3) + tmp*d3f(3, k1, k2)
432 : END DO
433 : END DO
434 0 : DO k3 = 1, 3
435 0 : DO k2 = 1, 3
436 0 : DO k1 = 1, 3
437 0 : tmp = -M2A(k1, k2, imA)*M1B(k3, imB) + M1A(k1, imA)*M2B(k2, k3, imB)
438 0 : wv(1) = wv(1) + tmp*d4f(1, k1, k2, k3)
439 0 : wv(2) = wv(2) + tmp*d4f(2, k1, k2, k3)
440 0 : wv(3) = wv(3) + tmp*d4f(3, k1, k2, k3)
441 : END DO
442 : END DO
443 : END DO
444 :
445 0 : DO k4 = 1, 3
446 0 : DO k3 = 1, 3
447 0 : DO k2 = 1, 3
448 0 : DO k1 = 1, 3
449 0 : tmp = M2A(k1, k2, imA)*M2B(k3, k4, imB)
450 0 : wv(1) = wv(1) + tmp*d5f(1, k1, k2, k3, k4)
451 0 : wv(2) = wv(2) + tmp*d5f(2, k1, k2, k3, k4)
452 0 : wv(3) = wv(3) + tmp*d5f(3, k1, k2, k3, k4)
453 : END DO
454 : END DO
455 : END DO
456 : END DO
457 :
458 0 : dCoul(1, imA, imB) = wv(1)
459 0 : dCoul(2, imA, imB) = wv(2)
460 0 : dCoul(3, imA, imB) = wv(3)
461 : END DO
462 : END DO
463 : END DO
464 : END DO
465 :
466 0 : END SUBROUTINE makedCoul
467 :
468 : ! **************************************************************************************************
469 : !> \brief Computes nuclei-nuclei interactions
470 : !> \param sepi ...
471 : !> \param sepj ...
472 : !> \param rijv ...
473 : !> \param enuc ...
474 : !> \param denuc ...
475 : !> \param se_int_control ...
476 : ! **************************************************************************************************
477 3813 : SUBROUTINE corecore_gks(sepi, sepj, rijv, enuc, denuc, se_int_control)
478 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
479 : REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
480 : REAL(dp), INTENT(OUT), OPTIONAL :: enuc
481 : REAL(dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: denuc
482 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
483 :
484 : LOGICAL :: l_denuc, l_enuc
485 : REAL(dp) :: alpi, alpj, dscale, rij, scale, zz
486 : REAL(kind=dp), DIMENSION(3, 45, 45) :: dCoul, dCoulE
487 : REAL(kind=dp), DIMENSION(45, 45) :: Coul, CoulE
488 : TYPE(se_int_control_type) :: se_int_control_off
489 :
490 15252 : rij = DOT_PRODUCT(rijv, rijv)
491 :
492 3813 : l_enuc = PRESENT(enuc)
493 3813 : l_denuc = PRESENT(denuc)
494 3813 : IF ((rij > rij_threshold) .AND. (l_enuc .OR. l_denuc)) THEN
495 :
496 3810 : rij = SQRT(rij)
497 :
498 3810 : IF (se_int_control%shortrange) THEN
499 : CALL setup_se_int_control_type(se_int_control_off, shortrange=.FALSE., do_ewald_r3=.FALSE., &
500 : do_ewald_gks=.FALSE., integral_screening=se_int_control%integral_screening, &
501 3807 : max_multipole=do_multipole_none, pc_coulomb_int=.FALSE.)
502 3807 : CALL makeCoul(rijv, sepi, sepj, Coul, se_int_control_off)
503 3807 : IF (l_denuc) CALL makedCoul(rijv, sepi, sepj, dCoul, se_int_control_off)
504 3807 : IF (se_int_control%do_ewald_gks) THEN
505 3 : CALL makeCoulE(rijv, sepi, sepj, CoulE, se_int_control)
506 3 : IF (l_denuc) CALL makedCoulE(rijv, sepi, sepj, dCoulE, se_int_control)
507 : ELSE
508 3804 : CALL makeCoul(rijv, sepi, sepj, CoulE, se_int_control)
509 3804 : IF (l_denuc) CALL makedCoul(rijv, sepi, sepj, dCoulE, se_int_control)
510 : END IF
511 : ELSE
512 3 : CALL makeCoul(rijv, sepi, sepj, Coul, se_int_control)
513 3 : CoulE = Coul
514 3 : IF (l_denuc) CALL makedCoul(rijv, sepi, sepj, dCoul, se_int_control)
515 3 : IF (l_denuc) dCoulE = dCoul
516 : END IF
517 :
518 3810 : scale = 0.0_dp
519 3810 : dscale = 0.0_dp
520 3810 : zz = sepi%zeff*sepj%zeff
521 3810 : alpi = sepi%alp
522 3810 : alpj = sepj%alp
523 3810 : scale = EXP(-alpi*rij) + EXP(-alpj*rij)
524 3810 : IF (l_enuc) THEN
525 3810 : enuc = zz*CoulE(1, 1) + scale*zz*Coul(1, 1)
526 : END IF
527 3810 : IF (l_denuc) THEN
528 0 : dscale = -alpi*EXP(-alpi*rij) - alpj*EXP(-alpj*rij)
529 0 : denuc(1) = zz*dCoulE(1, 1, 1) + dscale*(rijv(1)/rij)*zz*Coul(1, 1) + scale*zz*dCoul(1, 1, 1)
530 0 : denuc(2) = zz*dCoulE(2, 1, 1) + dscale*(rijv(2)/rij)*zz*Coul(1, 1) + scale*zz*dCoul(2, 1, 1)
531 0 : denuc(3) = zz*dCoulE(3, 1, 1) + dscale*(rijv(3)/rij)*zz*Coul(1, 1) + scale*zz*dCoul(3, 1, 1)
532 : END IF
533 :
534 : ELSE
535 :
536 3 : IF (se_int_control%do_ewald_gks) THEN
537 3 : zz = sepi%zeff*sepi%zeff
538 3 : CALL makeCoulE0(sepi, CoulE, se_int_control)
539 3 : IF (l_enuc) THEN
540 3 : enuc = zz*CoulE(1, 1)
541 : END IF
542 3 : IF (l_denuc) THEN
543 0 : denuc = 0._dp
544 : END IF
545 : END IF
546 :
547 : END IF
548 3813 : END SUBROUTINE corecore_gks
549 :
550 : ! **************************************************************************************************
551 : !> \brief Computes the primitives of the integrals (periodic case)
552 : !> \param RAB ...
553 : !> \param sepi ...
554 : !> \param sepj ...
555 : !> \param Coul ...
556 : !> \param se_int_control ...
557 : ! **************************************************************************************************
558 27 : SUBROUTINE makeCoulE(RAB, sepi, sepj, Coul, se_int_control)
559 : REAL(KIND=dp), DIMENSION(3) :: RAB
560 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
561 : REAL(KIND=dp), DIMENSION(45, 45), INTENT(OUT) :: Coul
562 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
563 :
564 : INTEGER :: gpt, imA, imB, k1, k2, k3, k4, lp, mp, np
565 9 : INTEGER, DIMENSION(:, :), POINTER :: bds
566 : REAL(KIND=dp) :: a2, ACOULA, ACOULB, alpha, cc, d1, d1f(3), d2, d2f(3, 3), d3, d3f(3, 3, 3), &
567 : d4, d4f(3, 3, 3, 3), f, ff, kr, kr2, r1, r2, r3, r5, r7, r9, rr, ss, w, w0, w1, w2, w3, &
568 : w4, w5
569 : REAL(KIND=dp), DIMENSION(3) :: kk, v
570 : REAL(KIND=dp), DIMENSION(3, 3, 45) :: M2A, M2B
571 : REAL(KIND=dp), DIMENSION(3, 45) :: M1A, M1B
572 : REAL(KIND=dp), DIMENSION(45) :: M0A, M0B
573 9 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rho0
574 : TYPE(dg_rho0_type), POINTER :: dg_rho0
575 : TYPE(dg_type), POINTER :: dg
576 : TYPE(pw_grid_type), POINTER :: pw_grid
577 : TYPE(pw_pool_type), POINTER :: pw_pool
578 :
579 9 : alpha = se_int_control%ewald_gks%alpha
580 9 : pw_pool => se_int_control%ewald_gks%pw_pool
581 9 : dg => se_int_control%ewald_gks%dg
582 9 : CALL dg_get(dg, dg_rho0=dg_rho0)
583 9 : rho0 => dg_rho0%density%array
584 9 : pw_grid => pw_pool%pw_grid
585 9 : bds => pw_grid%bounds
586 :
587 9 : CALL get_se_slater_multipole(sepi, M0A, M1A, M2A, ACOULA)
588 9 : CALL get_se_slater_multipole(sepj, M0B, M1B, M2B, ACOULB)
589 :
590 9 : v(1) = RAB(1)
591 9 : v(2) = RAB(2)
592 9 : v(3) = RAB(3)
593 36 : rr = SQRT(DOT_PRODUCT(v, v))
594 :
595 9 : r1 = 1.0_dp/rr
596 9 : r2 = r1*r1
597 9 : r3 = r2*r1
598 9 : r5 = r3*r2
599 9 : r7 = r5*r2
600 9 : r9 = r7*r2
601 :
602 9 : a2 = 0.5_dp*(1.0_dp/ACOULA + 1.0_dp/ACOULB)
603 :
604 9 : w0 = a2*rr
605 9 : w = EXP(-w0)
606 9 : w1 = (1.0_dp + 0.5_dp*w0)
607 9 : w2 = (w1 + 0.5_dp*w0 + 0.5_dp*w0**2)
608 9 : w3 = (w2 + w0**3/6.0_dp)
609 9 : w4 = (w3 + w0**4/30.0_dp)
610 9 : w5 = (w3 + 8.0_dp*w0**4/210.0_dp + w0**5/210.0_dp)
611 :
612 9 : f = (1.0_dp - w*w1)*r1
613 9 : d1 = -1.0_dp*(1.0_dp - w*w2)*r3
614 9 : d2 = 3.0_dp*(1.0_dp - w*w3)*r5
615 9 : d3 = -15.0_dp*(1.0_dp - w*w4)*r7
616 9 : d4 = 105.0_dp*(1.0_dp - w*w5)*r9
617 :
618 9 : kr = alpha*rr
619 9 : kr2 = kr*kr
620 9 : w0 = 1.0_dp - erfc(kr)
621 9 : w1 = 2.0_dp*oorootpi*EXP(-kr2)
622 9 : w2 = w1*kr
623 :
624 9 : f = f - w0*r1
625 9 : d1 = d1 + (-w2 + w0)*r3
626 9 : d2 = d2 + (w2*(3.0_dp + kr2*2.0_dp) - 3.0_dp*w0)*r5
627 9 : d3 = d3 + (-w2*(15.0_dp + kr2*(10.0_dp + kr2*4.0_dp)) + 15.0_dp*w0)*r7
628 9 : d4 = d4 + (w2*(105.0_dp + kr2*(70.0_dp + kr2*(28.0_dp + kr2*8.0_dp))) - 105.0_dp*w0)*r9
629 :
630 9 : CALL build_d_tensor_gks(d1f, d2f, d3f, d4f, v=v, d1=d1, d2=d2, d3=d3, d4=d4)
631 :
632 99 : DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2
633 999 : DO imB = 1, (sepj%natorb*(sepj%natorb + 1))/2
634 :
635 900 : w = M0A(imA)*M0B(imB)*f
636 3600 : DO k1 = 1, 3
637 3600 : w = w + (M1A(k1, imA)*M0B(imB) - M0A(imA)*M1B(k1, imB))*d1f(k1)
638 : END DO
639 3600 : DO k2 = 1, 3
640 11700 : DO k1 = 1, 3
641 10800 : w = w + (M2A(k1, k2, imA)*M0B(imB) - M1A(k1, imA)*M1B(k2, imB) + M0A(imA)*M2B(k1, k2, imB))*d2f(k1, k2)
642 : END DO
643 : END DO
644 3600 : DO k3 = 1, 3
645 11700 : DO k2 = 1, 3
646 35100 : DO k1 = 1, 3
647 32400 : w = w + (-M2A(k1, k2, imA)*M1B(k3, imB) + M1A(k1, imA)*M2B(k2, k3, imB))*d3f(k1, k2, k3)
648 : END DO
649 : END DO
650 : END DO
651 :
652 3600 : DO k4 = 1, 3
653 11700 : DO k3 = 1, 3
654 35100 : DO k2 = 1, 3
655 105300 : DO k1 = 1, 3
656 97200 : w = w + M2A(k1, k2, imA)*M2B(k3, k4, imB)*d4f(k1, k2, k3, k4)
657 : END DO
658 : END DO
659 : END DO
660 : END DO
661 :
662 990 : Coul(imA, imB) = w
663 : END DO
664 : END DO
665 :
666 : v(1) = RAB(1)
667 : v(2) = RAB(2)
668 : v(3) = RAB(3)
669 :
670 9 : f = 0.0_dp
671 9 : d1f = 0.0_dp
672 9 : d2f = 0.0_dp
673 9 : d3f = 0.0_dp
674 9 : d4f = 0.0_dp
675 :
676 150183 : DO gpt = 1, pw_grid%ngpts_cut_local
677 150174 : lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
678 150174 : mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
679 150174 : np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
680 :
681 150174 : lp = lp + bds(1, 1)
682 150174 : mp = mp + bds(1, 2)
683 150174 : np = np + bds(1, 3)
684 :
685 150174 : IF (pw_grid%gsq(gpt) == 0.0_dp) CYCLE
686 600660 : kk(:) = pw_grid%g(:, gpt)
687 150165 : ff = 2.0_dp*fourpi*rho0(lp, mp, np)**2*pw_grid%vol/pw_grid%gsq(gpt)
688 :
689 600660 : kr = DOT_PRODUCT(kk, v)
690 150165 : cc = COS(kr)
691 150165 : ss = SIN(kr)
692 :
693 150165 : f = f + cc*ff
694 600660 : DO k1 = 1, 3
695 600660 : d1f(k1) = d1f(k1) - kk(k1)*ss*ff
696 : END DO
697 600660 : DO k2 = 1, 3
698 1952145 : DO k1 = 1, 3
699 1801980 : d2f(k1, k2) = d2f(k1, k2) - kk(k1)*kk(k2)*cc*ff
700 : END DO
701 : END DO
702 600660 : DO k3 = 1, 3
703 1952145 : DO k2 = 1, 3
704 5856435 : DO k1 = 1, 3
705 5405940 : d3f(k1, k2, k3) = d3f(k1, k2, k3) + kk(k1)*kk(k2)*kk(k3)*ss*ff
706 : END DO
707 : END DO
708 : END DO
709 600669 : DO k4 = 1, 3
710 1952154 : DO k3 = 1, 3
711 5856435 : DO k2 = 1, 3
712 17569305 : DO k1 = 1, 3
713 16217820 : d4f(k1, k2, k3, k4) = d4f(k1, k2, k3, k4) + kk(k1)*kk(k2)*kk(k3)*kk(k4)*cc*ff
714 : END DO
715 : END DO
716 : END DO
717 : END DO
718 :
719 : END DO
720 :
721 99 : DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2
722 999 : DO imB = 1, (sepj%natorb*(sepj%natorb + 1))/2
723 :
724 900 : w = M0A(imA)*M0B(imB)*f
725 3600 : DO k1 = 1, 3
726 3600 : w = w + (M1A(k1, imA)*M0B(imB) - M0A(imA)*M1B(k1, imB))*d1f(k1)
727 : END DO
728 3600 : DO k2 = 1, 3
729 11700 : DO k1 = 1, 3
730 10800 : w = w + (M2A(k1, k2, imA)*M0B(imB) - M1A(k1, imA)*M1B(k2, imB) + M0A(imA)*M2B(k1, k2, imB))*d2f(k1, k2)
731 : END DO
732 : END DO
733 3600 : DO k3 = 1, 3
734 11700 : DO k2 = 1, 3
735 35100 : DO k1 = 1, 3
736 32400 : w = w + (-M2A(k1, k2, imA)*M1B(k3, imB) + M1A(k1, imA)*M2B(k2, k3, imB))*d3f(k1, k2, k3)
737 : END DO
738 : END DO
739 : END DO
740 :
741 3600 : DO k4 = 1, 3
742 11700 : DO k3 = 1, 3
743 35100 : DO k2 = 1, 3
744 105300 : DO k1 = 1, 3
745 97200 : w = w + M2A(k1, k2, imA)*M2B(k3, k4, imB)*d4f(k1, k2, k3, k4)
746 : END DO
747 : END DO
748 : END DO
749 : END DO
750 :
751 990 : Coul(imA, imB) = Coul(imA, imB) + w
752 :
753 : END DO
754 : END DO
755 :
756 99 : DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2
757 999 : DO imB = 1, (sepj%natorb*(sepj%natorb + 1))/2
758 900 : w = -M0A(imA)*M0B(imB)*0.25_dp*fourpi/(pw_grid%vol*alpha**2)
759 990 : Coul(imA, imB) = Coul(imA, imB) + w
760 : END DO
761 : END DO
762 :
763 9 : END SUBROUTINE makeCoulE
764 :
765 : ! **************************************************************************************************
766 : !> \brief Computes the derivatives of the primitives of the integrals
767 : !> (periodic case)
768 : !> \param RAB ...
769 : !> \param sepi ...
770 : !> \param sepj ...
771 : !> \param dCoul ...
772 : !> \param se_int_control ...
773 : ! **************************************************************************************************
774 0 : SUBROUTINE makedCoulE(RAB, sepi, sepj, dCoul, se_int_control)
775 : REAL(KIND=dp), DIMENSION(3) :: RAB
776 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
777 : REAL(KIND=dp), DIMENSION(3, 45, 45), INTENT(OUT) :: dCoul
778 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
779 :
780 : INTEGER :: gpt, imA, imB, k1, k2, k3, k4, k5, lp, &
781 : mp, np
782 0 : INTEGER, DIMENSION(:, :), POINTER :: bds
783 : REAL(KIND=dp) :: a2, ACOULA, ACOULB, alpha, cc, d1, d1f(3), d2, d2f(3, 3), d3, d3f(3, 3, 3), &
784 : d4, d4f(3, 3, 3, 3), d5, d5f(3, 3, 3, 3, 3), f, ff, kr, kr2, r1, r11, r2, r3, r5, r7, r9, &
785 : rr, ss, tmp, w, w0, w1, w2, w3, w4, w5, w6
786 : REAL(KIND=dp), DIMENSION(3) :: kk, v, wv
787 : REAL(kind=dp), DIMENSION(3, 3, 45) :: M2A, M2B
788 : REAL(kind=dp), DIMENSION(3, 45) :: M1A, M1B
789 : REAL(kind=dp), DIMENSION(45) :: M0A, M0B
790 0 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rho0
791 : TYPE(dg_rho0_type), POINTER :: dg_rho0
792 : TYPE(dg_type), POINTER :: dg
793 : TYPE(pw_grid_type), POINTER :: pw_grid
794 : TYPE(pw_pool_type), POINTER :: pw_pool
795 :
796 0 : alpha = se_int_control%ewald_gks%alpha
797 0 : pw_pool => se_int_control%ewald_gks%pw_pool
798 0 : dg => se_int_control%ewald_gks%dg
799 0 : CALL dg_get(dg, dg_rho0=dg_rho0)
800 0 : rho0 => dg_rho0%density%array
801 0 : pw_grid => pw_pool%pw_grid
802 0 : bds => pw_grid%bounds
803 :
804 0 : CALL get_se_slater_multipole(sepi, M0A, M1A, M2A, ACOULA)
805 0 : CALL get_se_slater_multipole(sepj, M0B, M1B, M2B, ACOULB)
806 :
807 0 : v(1) = RAB(1)
808 0 : v(2) = RAB(2)
809 0 : v(3) = RAB(3)
810 0 : rr = SQRT(DOT_PRODUCT(v, v))
811 :
812 0 : a2 = 0.5_dp*(1.0_dp/ACOULA + 1.0_dp/ACOULB)
813 :
814 0 : r1 = 1.0_dp/rr
815 0 : r2 = r1*r1
816 0 : r3 = r2*r1
817 0 : r5 = r3*r2
818 0 : r7 = r5*r2
819 0 : r9 = r7*r2
820 0 : r11 = r9*r2
821 :
822 0 : w0 = a2*rr
823 0 : w = EXP(-w0)
824 0 : w1 = (1.0_dp + 0.5_dp*w0)
825 0 : w2 = (w1 + 0.5_dp*w0 + 0.5_dp*w0**2)
826 0 : w3 = (w2 + w0**3/6.0_dp)
827 0 : w4 = (w3 + w0**4/30.0_dp)
828 0 : w5 = (w3 + 8.0_dp*w0**4/210.0_dp + w0**5/210.0_dp)
829 0 : w6 = (w3 + 5.0_dp*w0**4/126.0_dp + 2.0_dp*w0**5/315.0_dp + w0**6/1890.0_dp)
830 :
831 0 : f = (1.0_dp - w*w1)*r1
832 0 : d1 = -1.0_dp*(1.0_dp - w*w2)*r3
833 0 : d2 = 3.0_dp*(1.0_dp - w*w3)*r5
834 0 : d3 = -15.0_dp*(1.0_dp - w*w4)*r7
835 0 : d4 = 105.0_dp*(1.0_dp - w*w5)*r9
836 0 : d5 = -945.0_dp*(1.0_dp - w*w6)*r11
837 :
838 0 : kr = alpha*rr
839 0 : kr2 = kr*kr
840 0 : w0 = 1.0_dp - erfc(kr)
841 0 : w1 = 2.0_dp*oorootpi*EXP(-kr2)
842 0 : w2 = w1*kr
843 :
844 0 : f = f - w0*r1
845 0 : d1 = d1 + (-w2 + w0)*r3
846 0 : d2 = d2 + (w2*(3.0_dp + kr2*2.0_dp) - 3.0_dp*w0)*r5
847 0 : d3 = d3 + (-w2*(15.0_dp + kr2*(10.0_dp + kr2*4.0_dp)) + 15.0_dp*w0)*r7
848 0 : d4 = d4 + (w2*(105.0_dp + kr2*(70.0_dp + kr2*(28.0_dp + kr2*8.0_dp))) - 105.0_dp*w0)*r9
849 0 : d5 = d5 + (-w2*(945.0_dp + kr2*(630.0_dp + kr2*(252.0_dp + kr2*(72.0_dp + kr2*16.0_dp)))) + 945.0_dp*w0)*r11
850 :
851 0 : CALL build_d_tensor_gks(d1f, d2f, d3f, d4f, d5f, v, d1, d2, d3, d4, d5)
852 :
853 0 : DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2
854 0 : DO imB = 1, (sepj%natorb*(sepj%natorb + 1))/2
855 :
856 0 : tmp = M0A(imA)*M0B(imB)
857 0 : wv(1) = tmp*d1f(1)
858 0 : wv(2) = tmp*d1f(2)
859 0 : wv(3) = tmp*d1f(3)
860 :
861 0 : DO k1 = 1, 3
862 0 : tmp = M1A(k1, imA)*M0B(imB) - M0A(imA)*M1B(k1, imB)
863 0 : wv(1) = wv(1) + tmp*d2f(1, k1)
864 0 : wv(2) = wv(2) + tmp*d2f(2, k1)
865 0 : wv(3) = wv(3) + tmp*d2f(3, k1)
866 : END DO
867 0 : DO k2 = 1, 3
868 0 : DO k1 = 1, 3
869 0 : tmp = M2A(k1, k2, imA)*M0B(imB) - M1A(k1, imA)*M1B(k2, imB) + M0A(imA)*M2B(k1, k2, imB)
870 0 : wv(1) = wv(1) + tmp*d3f(1, k1, k2)
871 0 : wv(2) = wv(2) + tmp*d3f(2, k1, k2)
872 0 : wv(3) = wv(3) + tmp*d3f(3, k1, k2)
873 : END DO
874 : END DO
875 0 : DO k3 = 1, 3
876 0 : DO k2 = 1, 3
877 0 : DO k1 = 1, 3
878 0 : tmp = -M2A(k1, k2, imA)*M1B(k3, imB) + M1A(k1, imA)*M2B(k2, k3, imB)
879 0 : wv(1) = wv(1) + tmp*d4f(1, k1, k2, k3)
880 0 : wv(2) = wv(2) + tmp*d4f(2, k1, k2, k3)
881 0 : wv(3) = wv(3) + tmp*d4f(3, k1, k2, k3)
882 : END DO
883 : END DO
884 : END DO
885 :
886 0 : DO k4 = 1, 3
887 0 : DO k3 = 1, 3
888 0 : DO k2 = 1, 3
889 0 : DO k1 = 1, 3
890 0 : tmp = M2A(k1, k2, imA)*M2B(k3, k4, imB)
891 0 : wv(1) = wv(1) + tmp*d5f(1, k1, k2, k3, k4)
892 0 : wv(2) = wv(2) + tmp*d5f(2, k1, k2, k3, k4)
893 0 : wv(3) = wv(3) + tmp*d5f(3, k1, k2, k3, k4)
894 : END DO
895 : END DO
896 : END DO
897 : END DO
898 :
899 0 : dCoul(1, imA, imB) = wv(1)
900 0 : dCoul(2, imA, imB) = wv(2)
901 0 : dCoul(3, imA, imB) = wv(3)
902 : END DO
903 : END DO
904 :
905 : v(1) = RAB(1)
906 : v(2) = RAB(2)
907 : v(3) = RAB(3)
908 :
909 0 : f = 0.0_dp
910 0 : d1f = 0.0_dp
911 0 : d2f = 0.0_dp
912 0 : d3f = 0.0_dp
913 0 : d4f = 0.0_dp
914 0 : d5f = 0.0_dp
915 :
916 0 : DO gpt = 1, pw_grid%ngpts_cut_local
917 0 : lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
918 0 : mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
919 0 : np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
920 :
921 0 : lp = lp + bds(1, 1)
922 0 : mp = mp + bds(1, 2)
923 0 : np = np + bds(1, 3)
924 :
925 0 : IF (pw_grid%gsq(gpt) == 0.0_dp) CYCLE
926 0 : kk(:) = pw_grid%g(:, gpt)
927 0 : ff = 2.0_dp*fourpi*rho0(lp, mp, np)**2*pw_grid%vol/pw_grid%gsq(gpt)
928 :
929 0 : kr = DOT_PRODUCT(kk, v)
930 0 : cc = COS(kr)
931 0 : ss = SIN(kr)
932 :
933 0 : f = f + cc*ff
934 0 : DO k1 = 1, 3
935 0 : d1f(k1) = d1f(k1) - kk(k1)*ss*ff
936 : END DO
937 0 : DO k2 = 1, 3
938 0 : DO k1 = 1, 3
939 0 : d2f(k1, k2) = d2f(k1, k2) - kk(k1)*kk(k2)*cc*ff
940 : END DO
941 : END DO
942 0 : DO k3 = 1, 3
943 0 : DO k2 = 1, 3
944 0 : DO k1 = 1, 3
945 0 : d3f(k1, k2, k3) = d3f(k1, k2, k3) + kk(k1)*kk(k2)*kk(k3)*ss*ff
946 : END DO
947 : END DO
948 : END DO
949 0 : DO k4 = 1, 3
950 0 : DO k3 = 1, 3
951 0 : DO k2 = 1, 3
952 0 : DO k1 = 1, 3
953 0 : d4f(k1, k2, k3, k4) = d4f(k1, k2, k3, k4) + kk(k1)*kk(k2)*kk(k3)*kk(k4)*cc*ff
954 : END DO
955 : END DO
956 : END DO
957 : END DO
958 0 : DO k5 = 1, 3
959 0 : DO k4 = 1, 3
960 0 : DO k3 = 1, 3
961 0 : DO k2 = 1, 3
962 0 : DO k1 = 1, 3
963 0 : d5f(k1, k2, k3, k4, k5) = d5f(k1, k2, k3, k4, k5) - kk(k1)*kk(k2)*kk(k3)*kk(k4)*kk(k5)*ss*ff
964 : END DO
965 : END DO
966 : END DO
967 : END DO
968 : END DO
969 : END DO
970 :
971 0 : DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2
972 0 : DO imB = 1, (sepj%natorb*(sepj%natorb + 1))/2
973 0 : tmp = M0A(imA)*M0B(imB)
974 0 : wv(1) = tmp*d1f(1)
975 0 : wv(2) = tmp*d1f(2)
976 0 : wv(3) = tmp*d1f(3)
977 0 : DO k1 = 1, 3
978 0 : tmp = M1A(k1, imA)*M0B(imB) - M0A(imA)*M1B(k1, imB)
979 0 : wv(1) = wv(1) + tmp*d2f(1, k1)
980 0 : wv(2) = wv(2) + tmp*d2f(2, k1)
981 0 : wv(3) = wv(3) + tmp*d2f(3, k1)
982 : END DO
983 0 : DO k2 = 1, 3
984 0 : DO k1 = 1, 3
985 0 : tmp = M2A(k1, k2, imA)*M0B(imB) - M1A(k1, imA)*M1B(k2, imB) + M0A(imA)*M2B(k1, k2, imB)
986 0 : wv(1) = wv(1) + tmp*d3f(1, k1, k2)
987 0 : wv(2) = wv(2) + tmp*d3f(2, k1, k2)
988 0 : wv(3) = wv(3) + tmp*d3f(3, k1, k2)
989 : END DO
990 : END DO
991 0 : DO k3 = 1, 3
992 0 : DO k2 = 1, 3
993 0 : DO k1 = 1, 3
994 0 : tmp = -M2A(k1, k2, imA)*M1B(k3, imB) + M1A(k1, imA)*M2B(k2, k3, imB)
995 0 : wv(1) = wv(1) + tmp*d4f(1, k1, k2, k3)
996 0 : wv(2) = wv(2) + tmp*d4f(2, k1, k2, k3)
997 0 : wv(3) = wv(3) + tmp*d4f(3, k1, k2, k3)
998 : END DO
999 : END DO
1000 : END DO
1001 :
1002 0 : DO k4 = 1, 3
1003 0 : DO k3 = 1, 3
1004 0 : DO k2 = 1, 3
1005 0 : DO k1 = 1, 3
1006 0 : tmp = M2A(k1, k2, imA)*M2B(k3, k4, imB)
1007 0 : wv(1) = wv(1) + tmp*d5f(1, k1, k2, k3, k4)
1008 0 : wv(2) = wv(2) + tmp*d5f(2, k1, k2, k3, k4)
1009 0 : wv(3) = wv(3) + tmp*d5f(3, k1, k2, k3, k4)
1010 : END DO
1011 : END DO
1012 : END DO
1013 : END DO
1014 :
1015 0 : dCoul(1, imA, imB) = dCoul(1, imA, imB) + wv(1)
1016 0 : dCoul(2, imA, imB) = dCoul(2, imA, imB) + wv(2)
1017 0 : dCoul(3, imA, imB) = dCoul(3, imA, imB) + wv(3)
1018 : END DO
1019 : END DO
1020 :
1021 0 : END SUBROUTINE makedCoulE
1022 :
1023 : ! **************************************************************************************************
1024 : !> \brief Builds the tensor for the evaluation of the integrals with the
1025 : !> cartesian multipoles
1026 : !> \param d1f ...
1027 : !> \param d2f ...
1028 : !> \param d3f ...
1029 : !> \param d4f ...
1030 : !> \param d5f ...
1031 : !> \param v ...
1032 : !> \param d1 ...
1033 : !> \param d2 ...
1034 : !> \param d3 ...
1035 : !> \param d4 ...
1036 : !> \param d5 ...
1037 : ! **************************************************************************************************
1038 16178 : SUBROUTINE build_d_tensor_gks(d1f, d2f, d3f, d4f, d5f, v, d1, d2, d3, d4, d5)
1039 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: d1f
1040 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: d2f
1041 : REAL(KIND=dp), DIMENSION(3, 3, 3), INTENT(OUT) :: d3f
1042 : REAL(KIND=dp), DIMENSION(3, 3, 3, 3), INTENT(OUT) :: d4f
1043 : REAL(KIND=dp), DIMENSION(3, 3, 3, 3, 3), &
1044 : INTENT(OUT), OPTIONAL :: d5f
1045 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: v
1046 : REAL(KIND=dp), INTENT(IN) :: d1, d2, d3, d4
1047 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: d5
1048 :
1049 : INTEGER :: k1, k2, k3, k4, k5
1050 : REAL(KIND=dp) :: w
1051 :
1052 16178 : d1f = 0.0_dp
1053 16178 : d2f = 0.0_dp
1054 16178 : d3f = 0.0_dp
1055 16178 : d4f = 0.0_dp
1056 64712 : DO k1 = 1, 3
1057 64712 : d1f(k1) = d1f(k1) + v(k1)*d1
1058 : END DO
1059 64712 : DO k1 = 1, 3
1060 194136 : DO k2 = 1, 3
1061 194136 : d2f(k2, k1) = d2f(k2, k1) + v(k1)*v(k2)*d2
1062 : END DO
1063 64712 : d2f(k1, k1) = d2f(k1, k1) + d1
1064 : END DO
1065 64712 : DO k1 = 1, 3
1066 210314 : DO k2 = 1, 3
1067 582408 : DO k3 = 1, 3
1068 582408 : d3f(k3, k2, k1) = d3f(k3, k2, k1) + v(k1)*v(k2)*v(k3)*d3
1069 : END DO
1070 145602 : w = v(k1)*d2
1071 145602 : d3f(k1, k2, k2) = d3f(k1, k2, k2) + w
1072 145602 : d3f(k2, k1, k2) = d3f(k2, k1, k2) + w
1073 194136 : d3f(k2, k2, k1) = d3f(k2, k2, k1) + w
1074 : END DO
1075 : END DO
1076 64712 : DO k1 = 1, 3
1077 210314 : DO k2 = 1, 3
1078 582408 : DO k3 = 1, 3
1079 1747224 : DO k4 = 1, 3
1080 1747224 : d4f(k4, k3, k2, k1) = d4f(k4, k3, k2, k1) + v(k1)*v(k2)*v(k3)*v(k4)*d4
1081 : END DO
1082 436806 : w = v(k1)*v(k2)*d3
1083 436806 : d4f(k1, k2, k3, k3) = d4f(k1, k2, k3, k3) + w
1084 436806 : d4f(k1, k3, k2, k3) = d4f(k1, k3, k2, k3) + w
1085 436806 : d4f(k3, k1, k2, k3) = d4f(k3, k1, k2, k3) + w
1086 436806 : d4f(k1, k3, k3, k2) = d4f(k1, k3, k3, k2) + w
1087 436806 : d4f(k3, k1, k3, k2) = d4f(k3, k1, k3, k2) + w
1088 582408 : d4f(k3, k3, k1, k2) = d4f(k3, k3, k1, k2) + w
1089 : END DO
1090 145602 : d4f(k1, k1, k2, k2) = d4f(k1, k1, k2, k2) + d2
1091 145602 : d4f(k1, k2, k1, k2) = d4f(k1, k2, k1, k2) + d2
1092 194136 : d4f(k1, k2, k2, k1) = d4f(k1, k2, k2, k1) + d2
1093 : END DO
1094 : END DO
1095 16178 : IF (PRESENT(d5f) .AND. PRESENT(d5)) THEN
1096 0 : d5f = 0.0_dp
1097 :
1098 0 : DO k1 = 1, 3
1099 0 : DO k2 = 1, 3
1100 0 : DO k3 = 1, 3
1101 0 : DO k4 = 1, 3
1102 0 : DO k5 = 1, 3
1103 0 : d5f(k5, k4, k3, k2, k1) = d5f(k5, k4, k3, k2, k1) + v(k1)*v(k2)*v(k3)*v(k4)*v(k5)*d5
1104 : END DO
1105 0 : w = v(k1)*v(k2)*v(k3)*d4
1106 0 : d5f(k1, k2, k3, k4, k4) = d5f(k1, k2, k3, k4, k4) + w
1107 0 : d5f(k1, k2, k4, k3, k4) = d5f(k1, k2, k4, k3, k4) + w
1108 0 : d5f(k1, k4, k2, k3, k4) = d5f(k1, k4, k2, k3, k4) + w
1109 0 : d5f(k4, k1, k2, k3, k4) = d5f(k4, k1, k2, k3, k4) + w
1110 0 : d5f(k1, k2, k4, k4, k3) = d5f(k1, k2, k4, k4, k3) + w
1111 0 : d5f(k1, k4, k2, k4, k3) = d5f(k1, k4, k2, k4, k3) + w
1112 0 : d5f(k4, k1, k2, k4, k3) = d5f(k4, k1, k2, k4, k3) + w
1113 0 : d5f(k1, k4, k4, k2, k3) = d5f(k1, k4, k4, k2, k3) + w
1114 0 : d5f(k4, k1, k4, k2, k3) = d5f(k4, k1, k4, k2, k3) + w
1115 0 : d5f(k4, k4, k1, k2, k3) = d5f(k4, k4, k1, k2, k3) + w
1116 : END DO
1117 0 : w = v(k1)*d3
1118 0 : d5f(k1, k2, k2, k3, k3) = d5f(k1, k2, k2, k3, k3) + w
1119 0 : d5f(k1, k2, k3, k2, k3) = d5f(k1, k2, k3, k2, k3) + w
1120 0 : d5f(k1, k2, k3, k3, k2) = d5f(k1, k2, k3, k3, k2) + w
1121 0 : d5f(k2, k1, k2, k3, k3) = d5f(k2, k1, k2, k3, k3) + w
1122 0 : d5f(k2, k1, k3, k2, k3) = d5f(k2, k1, k3, k2, k3) + w
1123 0 : d5f(k2, k1, k3, k3, k2) = d5f(k2, k1, k3, k3, k2) + w
1124 0 : d5f(k2, k2, k1, k3, k3) = d5f(k2, k2, k1, k3, k3) + w
1125 0 : d5f(k2, k3, k1, k2, k3) = d5f(k2, k3, k1, k2, k3) + w
1126 0 : d5f(k2, k3, k1, k3, k2) = d5f(k2, k3, k1, k3, k2) + w
1127 0 : d5f(k2, k2, k3, k1, k3) = d5f(k2, k2, k3, k1, k3) + w
1128 0 : d5f(k2, k3, k2, k1, k3) = d5f(k2, k3, k2, k1, k3) + w
1129 0 : d5f(k2, k3, k3, k1, k2) = d5f(k2, k3, k3, k1, k2) + w
1130 0 : d5f(k2, k2, k3, k3, k1) = d5f(k2, k2, k3, k3, k1) + w
1131 0 : d5f(k2, k3, k2, k3, k1) = d5f(k2, k3, k2, k3, k1) + w
1132 0 : d5f(k2, k3, k3, k2, k1) = d5f(k2, k3, k3, k2, k1) + w
1133 : END DO
1134 : END DO
1135 : END DO
1136 : END IF
1137 16178 : END SUBROUTINE build_d_tensor_gks
1138 :
1139 : ! **************************************************************************************************
1140 : !> \brief ...
1141 : !> \param sepi ...
1142 : !> \param Coul ...
1143 : !> \param se_int_control ...
1144 : ! **************************************************************************************************
1145 9 : SUBROUTINE makeCoulE0(sepi, Coul, se_int_control)
1146 : TYPE(semi_empirical_type), POINTER :: sepi
1147 : REAL(KIND=dp), DIMENSION(45, 45), INTENT(OUT) :: Coul
1148 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
1149 :
1150 : INTEGER :: gpt, imA, imB, k1, k2, k3, k4, lp, mp, np
1151 9 : INTEGER, DIMENSION(:, :), POINTER :: bds
1152 : REAL(KIND=dp) :: alpha, d2f(3, 3), d4f(3, 3, 3, 3), f, &
1153 : ff, w
1154 : REAL(KIND=dp), DIMENSION(3) :: kk
1155 : REAL(KIND=dp), DIMENSION(3, 3, 45) :: M2A
1156 : REAL(KIND=dp), DIMENSION(3, 45) :: M1A
1157 : REAL(KIND=dp), DIMENSION(45) :: M0A
1158 9 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rho0
1159 : TYPE(dg_rho0_type), POINTER :: dg_rho0
1160 : TYPE(dg_type), POINTER :: dg
1161 : TYPE(pw_grid_type), POINTER :: pw_grid
1162 : TYPE(pw_pool_type), POINTER :: pw_pool
1163 :
1164 9 : alpha = se_int_control%ewald_gks%alpha
1165 9 : pw_pool => se_int_control%ewald_gks%pw_pool
1166 9 : dg => se_int_control%ewald_gks%dg
1167 9 : CALL dg_get(dg, dg_rho0=dg_rho0)
1168 9 : rho0 => dg_rho0%density%array
1169 9 : pw_grid => pw_pool%pw_grid
1170 9 : bds => pw_grid%bounds
1171 :
1172 9 : CALL get_se_slater_multipole(sepi, M0A, M1A, M2A)
1173 :
1174 9 : f = 0.0_dp
1175 9 : d2f = 0.0_dp
1176 9 : d4f = 0.0_dp
1177 :
1178 150183 : DO gpt = 1, pw_grid%ngpts_cut_local
1179 150174 : lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
1180 150174 : mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
1181 150174 : np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
1182 :
1183 150174 : lp = lp + bds(1, 1)
1184 150174 : mp = mp + bds(1, 2)
1185 150174 : np = np + bds(1, 3)
1186 :
1187 150174 : IF (pw_grid%gsq(gpt) == 0.0_dp) CYCLE
1188 600660 : kk(:) = pw_grid%g(:, gpt)
1189 150165 : ff = 2.0_dp*fourpi*rho0(lp, mp, np)**2*pw_grid%vol/pw_grid%gsq(gpt)
1190 :
1191 150165 : f = f + ff
1192 600660 : DO k2 = 1, 3
1193 1952145 : DO k1 = 1, 3
1194 1801980 : d2f(k1, k2) = d2f(k1, k2) - kk(k1)*kk(k2)*ff
1195 : END DO
1196 : END DO
1197 600669 : DO k4 = 1, 3
1198 1952154 : DO k3 = 1, 3
1199 5856435 : DO k2 = 1, 3
1200 17569305 : DO k1 = 1, 3
1201 16217820 : d4f(k1, k2, k3, k4) = d4f(k1, k2, k3, k4) + kk(k1)*kk(k2)*kk(k3)*kk(k4)*ff
1202 : END DO
1203 : END DO
1204 : END DO
1205 : END DO
1206 :
1207 : END DO
1208 :
1209 99 : DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2
1210 999 : DO imB = 1, (sepi%natorb*(sepi%natorb + 1))/2
1211 :
1212 900 : w = M0A(imA)*M0A(imB)*f
1213 3600 : DO k2 = 1, 3
1214 11700 : DO k1 = 1, 3
1215 10800 : w = w + (M2A(k1, k2, imA)*M0A(imB) - M1A(k1, imA)*M1A(k2, imB) + M0A(imA)*M2A(k1, k2, imB))*d2f(k1, k2)
1216 : END DO
1217 : END DO
1218 :
1219 3600 : DO k4 = 1, 3
1220 11700 : DO k3 = 1, 3
1221 35100 : DO k2 = 1, 3
1222 105300 : DO k1 = 1, 3
1223 97200 : w = w + M2A(k1, k2, imA)*M2A(k3, k4, imB)*d4f(k1, k2, k3, k4)
1224 : END DO
1225 : END DO
1226 : END DO
1227 : END DO
1228 :
1229 990 : Coul(imA, imB) = w
1230 :
1231 : END DO
1232 : END DO
1233 :
1234 99 : DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2
1235 999 : DO imB = 1, (sepi%natorb*(sepi%natorb + 1))/2
1236 900 : w = -M0A(imA)*M0A(imB)*0.25_dp*fourpi/(pw_grid%vol*alpha**2)
1237 990 : Coul(imA, imB) = Coul(imA, imB) + w
1238 : END DO
1239 : END DO
1240 :
1241 99 : DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2
1242 999 : DO imB = 1, (sepi%natorb*(sepi%natorb + 1))/2
1243 :
1244 900 : w = M0A(imA)*M0A(imB)
1245 900 : Coul(imA, imB) = Coul(imA, imB) - 2.0_dp*alpha*oorootpi*w
1246 :
1247 900 : w = 0.0_dp
1248 3600 : DO k1 = 1, 3
1249 2700 : w = w + M1A(k1, imA)*M1A(k1, imB)
1250 2700 : w = w - M0A(imA)*M2A(k1, k1, imB)
1251 3600 : w = w - M2A(k1, k1, imA)*M0A(imB)
1252 : END DO
1253 900 : Coul(imA, imB) = Coul(imA, imB) - 4.0_dp*alpha**3*oorootpi*w/3.0_dp
1254 :
1255 900 : w = 0.0_dp
1256 3600 : DO k2 = 1, 3
1257 11700 : DO k1 = 1, 3
1258 8100 : w = w + 2.0_dp*M2A(k1, k2, imA)*M2A(k1, k2, imB)
1259 10800 : w = w + M2A(k1, k1, imA)*M2A(k2, k2, imB)
1260 : END DO
1261 : END DO
1262 990 : Coul(imA, imB) = Coul(imA, imB) - 8.0_dp*alpha**5*oorootpi*w/5.0_dp
1263 : END DO
1264 : END DO
1265 9 : END SUBROUTINE makeCoulE0
1266 :
1267 : ! **************************************************************************************************
1268 : !> \brief Retrieves the multipole for the Slater integral evaluation
1269 : !> \param sepi ...
1270 : !> \param M0 ...
1271 : !> \param M1 ...
1272 : !> \param M2 ...
1273 : !> \param ACOUL ...
1274 : ! **************************************************************************************************
1275 32365 : SUBROUTINE get_se_slater_multipole(sepi, M0, M1, M2, ACOUL)
1276 : TYPE(semi_empirical_type), POINTER :: sepi
1277 : REAL(kind=dp), DIMENSION(45), INTENT(OUT) :: M0
1278 : REAL(kind=dp), DIMENSION(3, 45), INTENT(OUT) :: M1
1279 : REAL(kind=dp), DIMENSION(3, 3, 45), INTENT(OUT) :: M2
1280 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: ACOUL
1281 :
1282 : INTEGER :: i, j, jint, size_1c_int
1283 : TYPE(semi_empirical_mpole_type), POINTER :: mpole
1284 :
1285 32365 : NULLIFY (mpole)
1286 32365 : size_1c_int = SIZE(sepi%w_mpole)
1287 356015 : DO jint = 1, size_1c_int
1288 323650 : mpole => sepi%w_mpole(jint)%mpole
1289 323650 : i = mpole%indi
1290 323650 : j = mpole%indj
1291 323650 : M0(indexb(i, j)) = -mpole%cs
1292 :
1293 323650 : M1(1, indexb(i, j)) = -mpole%ds(1)
1294 323650 : M1(2, indexb(i, j)) = -mpole%ds(2)
1295 323650 : M1(3, indexb(i, j)) = -mpole%ds(3)
1296 :
1297 323650 : M2(1, 1, indexb(i, j)) = -mpole%qq(1, 1)/3.0_dp
1298 323650 : M2(2, 1, indexb(i, j)) = -mpole%qq(2, 1)/3.0_dp
1299 323650 : M2(3, 1, indexb(i, j)) = -mpole%qq(3, 1)/3.0_dp
1300 :
1301 323650 : M2(1, 2, indexb(i, j)) = -mpole%qq(1, 2)/3.0_dp
1302 323650 : M2(2, 2, indexb(i, j)) = -mpole%qq(2, 2)/3.0_dp
1303 323650 : M2(3, 2, indexb(i, j)) = -mpole%qq(3, 2)/3.0_dp
1304 :
1305 323650 : M2(1, 3, indexb(i, j)) = -mpole%qq(1, 3)/3.0_dp
1306 323650 : M2(2, 3, indexb(i, j)) = -mpole%qq(2, 3)/3.0_dp
1307 356015 : M2(3, 3, indexb(i, j)) = -mpole%qq(3, 3)/3.0_dp
1308 : END DO
1309 32365 : IF (PRESENT(ACOUL)) ACOUL = sepi%acoul
1310 32365 : END SUBROUTINE get_se_slater_multipole
1311 :
1312 : END MODULE semi_empirical_int_gks
|
