(git:374b731)
Loading...
Searching...
No Matches
semi_empirical_int_gks.F
Go to the documentation of this file.
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! **************************************************************************************************
20MODULE 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
47CONTAINS
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 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 rab = -rij
69
70 IF (se_int_control%do_ewald_gks) THEN
71 IF (dot_product(rij, rij) > rij_threshold) THEN
72 CALL makecoule(rab, sepi, sepj, coul, se_int_control)
73 ELSE
74 CALL makecoule0(sepi, coul, se_int_control)
75 END IF
76 ELSE
77 CALL makecoul(rab, sepi, sepj, coul, se_int_control)
78 END IF
79
80 i = 0
81 DO mu = 1, sepi%natorb
82 DO nu = 1, mu
83 i = i + 1
84 e1b(i) = -coul(i, 1)*sepj%zeff
85 END DO
86 END DO
87
88 i = 0
89 DO mu = 1, sepj%natorb
90 DO nu = 1, mu
91 i = i + 1
92 e2a(i) = -coul(1, i)*sepi%zeff
93 END DO
94 END DO
95
96 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 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 rab = -rij
117
118 IF (se_int_control%do_ewald_gks) THEN
119 IF (dot_product(rij, rij) > rij_threshold) THEN
120 CALL makecoule(rab, sepi, sepj, coul, se_int_control)
121 ELSE
122 CALL makecoule0(sepi, coul, se_int_control)
123 END IF
124 ELSE
125 CALL makecoul(rab, sepi, sepj, coul, se_int_control)
126 END IF
127
128 i = 0
129 ind1 = 0
130 DO mu = 1, sepi%natorb
131 DO nu = 1, mu
132 ind1 = ind1 + 1
133 ind2 = 0
134 DO lam = 1, sepj%natorb
135 DO sig = 1, lam
136 i = i + 1
137 ind2 = ind2 + 1
138 w(i) = coul(ind1, ind2)
139 END DO
140 END DO
141 END DO
142 END DO
143
144 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 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 rab = -rij
166
167 IF (se_int_control%do_ewald_gks) THEN
168 CALL makedcoule(rab, sepi, sepj, dcoul, se_int_control)
169 ELSE
170 CALL makedcoul(rab, sepi, sepj, dcoul, se_int_control)
171 END IF
172
173 i = 0
174 DO mu = 1, sepi%natorb
175 DO nu = 1, mu
176 i = i + 1
177 de1b(1, i) = dcoul(1, i, 1)*sepj%zeff
178 de1b(2, i) = dcoul(2, i, 1)*sepj%zeff
179 de1b(3, i) = dcoul(3, i, 1)*sepj%zeff
180 END DO
181 END DO
182
183 i = 0
184 DO mu = 1, sepj%natorb
185 DO nu = 1, mu
186 i = i + 1
187 de2a(1, i) = dcoul(1, 1, i)*sepi%zeff
188 de2a(2, i) = dcoul(2, 1, i)*sepi%zeff
189 de2a(3, i) = dcoul(3, 1, i)*sepi%zeff
190 END DO
191 END DO
192
193 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 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 rab = -rij
214
215 IF (se_int_control%do_ewald_gks) THEN
216 CALL makedcoule(rab, sepi, sepj, dcoul, se_int_control)
217 ELSE
218 CALL makedcoul(rab, sepi, sepj, dcoul, se_int_control)
219 END IF
220
221 i = 0
222 ind1 = 0
223 DO mu = 1, sepi%natorb
224 DO nu = 1, mu
225 ind1 = ind1 + 1
226 ind2 = 0
227 DO lam = 1, sepj%natorb
228 DO sig = 1, lam
229 i = i + 1
230 ind2 = ind2 + 1
231 dw(1, i) = -dcoul(1, ind1, ind2)
232 dw(2, i) = -dcoul(2, ind1, ind2)
233 dw(3, i) = -dcoul(3, ind1, ind2)
234 END DO
235 END DO
236 END DO
237 END DO
238
239 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 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 shortrange = se_int_control%shortrange
267 CALL get_se_slater_multipole(sepi, m0a, m1a, m2a, acoula)
268 CALL get_se_slater_multipole(sepj, m0b, m1b, m2b, acoulb)
269
270 v(1) = rab(1)
271 v(2) = rab(2)
272 v(3) = rab(3)
273 rr = sqrt(dot_product(v, v))
274
275 a2 = 0.5_dp*(1.0_dp/acoula + 1.0_dp/acoulb)
276 w0 = a2*rr
277 w = exp(-w0)
278 w1 = (1.0_dp + 0.5_dp*w0)
279 w2 = (w1 + 0.5_dp*w0 + 0.5_dp*w0**2)
280 w3 = (w2 + w0**3/6.0_dp)
281 w4 = (w3 + w0**4/30.0_dp)
282 w5 = (w3 + 8.0_dp*w0**4/210.0_dp + w0**5/210.0_dp)
283
284 IF (shortrange) THEN
285 f = (-w*w1)/rr
286 d1 = -1.0_dp*(-w*w2)/rr**3
287 d2 = 3.0_dp*(-w*w3)/rr**5
288 d3 = -15.0_dp*(-w*w4)/rr**7
289 d4 = 105.0_dp*(-w*w5)/rr**9
290 ELSE
291 f = (1.0_dp - w*w1)/rr
292 d1 = -1.0_dp*(1.0_dp - w*w2)/rr**3
293 d2 = 3.0_dp*(1.0_dp - w*w3)/rr**5
294 d3 = -15.0_dp*(1.0_dp - w*w4)/rr**7
295 d4 = 105.0_dp*(1.0_dp - w*w5)/rr**9
296 END IF
297
298 CALL build_d_tensor_gks(d1f, d2f, d3f, d4f, v=v, d1=d1, d2=d2, d3=d3, d4=d4)
299
300 ima = 0
301 DO ia = 1, sepi%natorb
302 DO ja = 1, ia
303 ima = ima + 1
304
305 imb = 0
306 DO ib = 1, sepj%natorb
307 DO jb = 1, ib
308 imb = imb + 1
309
310 w = m0a(ima)*m0b(imb)*f
311 DO k1 = 1, 3
312 w = w + (m1a(k1, ima)*m0b(imb) - m0a(ima)*m1b(k1, imb))*d1f(k1)
313 END DO
314 DO k2 = 1, 3
315 DO k1 = 1, 3
316 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 DO k3 = 1, 3
320 DO k2 = 1, 3
321 DO k1 = 1, 3
322 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 DO k4 = 1, 3
328 DO k3 = 1, 3
329 DO k2 = 1, 3
330 DO k1 = 1, 3
331 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 coul(ima, imb) = w
338 END DO
339 END DO
340 END DO
341 END DO
342
343 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 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 shortrange = se_int_control%shortrange
370 CALL get_se_slater_multipole(sepi, m0a, m1a, m2a, acoula)
371 CALL get_se_slater_multipole(sepj, m0b, m1b, m2b, acoulb)
372
373 v(1) = rab(1)
374 v(2) = rab(2)
375 v(3) = rab(3)
376 rr = sqrt(dot_product(v, v))
377
378 a2 = 0.5_dp*(1.0_dp/acoula + 1.0_dp/acoulb)
379 w0 = a2*rr
380 w = exp(-w0)
381 w1 = (1.0_dp + 0.5_dp*w0)
382 w2 = (w1 + 0.5_dp*w0 + 0.5_dp*w0**2)
383 w3 = (w2 + w0**3/6.0_dp)
384 w4 = (w3 + w0**4/30.0_dp)
385 w5 = (w3 + 4.0_dp*w0**4/105.0_dp + w0**5/210.0_dp)
386 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 IF (shortrange) THEN
389 f = (-w*w1)/rr
390 d1 = -1.0_dp*(-w*w2)/rr**3
391 d2 = 3.0_dp*(-w*w3)/rr**5
392 d3 = -15.0_dp*(-w*w4)/rr**7
393 d4 = 105.0_dp*(-w*w5)/rr**9
394 d5 = -945.0_dp*(-w*w6)/rr**11
395 ELSE
396 f = (1.0_dp - w*w1)/rr
397 d1 = -1.0_dp*(1.0_dp - w*w2)/rr**3
398 d2 = 3.0_dp*(1.0_dp - w*w3)/rr**5
399 d3 = -15.0_dp*(1.0_dp - w*w4)/rr**7
400 d4 = 105.0_dp*(1.0_dp - w*w5)/rr**9
401 d5 = -945.0_dp*(1.0_dp - w*w6)/rr**11
402 END IF
403
404 CALL build_d_tensor_gks(d1f, d2f, d3f, d4f, d5f, v, d1, d2, d3, d4, d5)
405
406 ima = 0
407 DO ia = 1, sepi%natorb
408 DO ja = 1, ia
409 ima = ima + 1
410
411 imb = 0
412 DO ib = 1, sepj%natorb
413 DO jb = 1, ib
414 imb = imb + 1
415
416 tmp = m0a(ima)*m0b(imb)
417 wv(1) = tmp*d1f(1)
418 wv(2) = tmp*d1f(2)
419 wv(3) = tmp*d1f(3)
420 DO k1 = 1, 3
421 tmp = m1a(k1, ima)*m0b(imb) - m0a(ima)*m1b(k1, imb)
422 wv(1) = wv(1) + tmp*d2f(1, k1)
423 wv(2) = wv(2) + tmp*d2f(2, k1)
424 wv(3) = wv(3) + tmp*d2f(3, k1)
425 END DO
426 DO k2 = 1, 3
427 DO k1 = 1, 3
428 tmp = m2a(k1, k2, ima)*m0b(imb) - m1a(k1, ima)*m1b(k2, imb) + m0a(ima)*m2b(k1, k2, imb)
429 wv(1) = wv(1) + tmp*d3f(1, k1, k2)
430 wv(2) = wv(2) + tmp*d3f(2, k1, k2)
431 wv(3) = wv(3) + tmp*d3f(3, k1, k2)
432 END DO
433 END DO
434 DO k3 = 1, 3
435 DO k2 = 1, 3
436 DO k1 = 1, 3
437 tmp = -m2a(k1, k2, ima)*m1b(k3, imb) + m1a(k1, ima)*m2b(k2, k3, imb)
438 wv(1) = wv(1) + tmp*d4f(1, k1, k2, k3)
439 wv(2) = wv(2) + tmp*d4f(2, k1, k2, k3)
440 wv(3) = wv(3) + tmp*d4f(3, k1, k2, k3)
441 END DO
442 END DO
443 END DO
444
445 DO k4 = 1, 3
446 DO k3 = 1, 3
447 DO k2 = 1, 3
448 DO k1 = 1, 3
449 tmp = m2a(k1, k2, ima)*m2b(k3, k4, imb)
450 wv(1) = wv(1) + tmp*d5f(1, k1, k2, k3, k4)
451 wv(2) = wv(2) + tmp*d5f(2, k1, k2, k3, k4)
452 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 dcoul(1, ima, imb) = wv(1)
459 dcoul(2, ima, imb) = wv(2)
460 dcoul(3, ima, imb) = wv(3)
461 END DO
462 END DO
463 END DO
464 END DO
465
466 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 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 rij = dot_product(rijv, rijv)
491
492 l_enuc = PRESENT(enuc)
493 l_denuc = PRESENT(denuc)
494 IF ((rij > rij_threshold) .AND. (l_enuc .OR. l_denuc)) THEN
495
496 rij = sqrt(rij)
497
498 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 max_multipole=do_multipole_none, pc_coulomb_int=.false.)
502 CALL makecoul(rijv, sepi, sepj, coul, se_int_control_off)
503 IF (l_denuc) CALL makedcoul(rijv, sepi, sepj, dcoul, se_int_control_off)
504 IF (se_int_control%do_ewald_gks) THEN
505 CALL makecoule(rijv, sepi, sepj, coule, se_int_control)
506 IF (l_denuc) CALL makedcoule(rijv, sepi, sepj, dcoule, se_int_control)
507 ELSE
508 CALL makecoul(rijv, sepi, sepj, coule, se_int_control)
509 IF (l_denuc) CALL makedcoul(rijv, sepi, sepj, dcoule, se_int_control)
510 END IF
511 ELSE
512 CALL makecoul(rijv, sepi, sepj, coul, se_int_control)
513 coule = coul
514 IF (l_denuc) CALL makedcoul(rijv, sepi, sepj, dcoul, se_int_control)
515 IF (l_denuc) dcoule = dcoul
516 END IF
517
518 scale = 0.0_dp
519 dscale = 0.0_dp
520 zz = sepi%zeff*sepj%zeff
521 alpi = sepi%alp
522 alpj = sepj%alp
523 scale = exp(-alpi*rij) + exp(-alpj*rij)
524 IF (l_enuc) THEN
525 enuc = zz*coule(1, 1) + scale*zz*coul(1, 1)
526 END IF
527 IF (l_denuc) THEN
528 dscale = -alpi*exp(-alpi*rij) - alpj*exp(-alpj*rij)
529 denuc(1) = zz*dcoule(1, 1, 1) + dscale*(rijv(1)/rij)*zz*coul(1, 1) + scale*zz*dcoul(1, 1, 1)
530 denuc(2) = zz*dcoule(2, 1, 1) + dscale*(rijv(2)/rij)*zz*coul(1, 1) + scale*zz*dcoul(2, 1, 1)
531 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 IF (se_int_control%do_ewald_gks) THEN
537 zz = sepi%zeff*sepi%zeff
538 CALL makecoule0(sepi, coule, se_int_control)
539 IF (l_enuc) THEN
540 enuc = zz*coule(1, 1)
541 END IF
542 IF (l_denuc) THEN
543 denuc = 0._dp
544 END IF
545 END IF
546
547 END IF
548 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 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 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 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 alpha = se_int_control%ewald_gks%alpha
580 pw_pool => se_int_control%ewald_gks%pw_pool
581 dg => se_int_control%ewald_gks%dg
582 CALL dg_get(dg, dg_rho0=dg_rho0)
583 rho0 => dg_rho0%density%array
584 pw_grid => pw_pool%pw_grid
585 bds => pw_grid%bounds
586
587 CALL get_se_slater_multipole(sepi, m0a, m1a, m2a, acoula)
588 CALL get_se_slater_multipole(sepj, m0b, m1b, m2b, acoulb)
589
590 v(1) = rab(1)
591 v(2) = rab(2)
592 v(3) = rab(3)
593 rr = sqrt(dot_product(v, v))
594
595 r1 = 1.0_dp/rr
596 r2 = r1*r1
597 r3 = r2*r1
598 r5 = r3*r2
599 r7 = r5*r2
600 r9 = r7*r2
601
602 a2 = 0.5_dp*(1.0_dp/acoula + 1.0_dp/acoulb)
603
604 w0 = a2*rr
605 w = exp(-w0)
606 w1 = (1.0_dp + 0.5_dp*w0)
607 w2 = (w1 + 0.5_dp*w0 + 0.5_dp*w0**2)
608 w3 = (w2 + w0**3/6.0_dp)
609 w4 = (w3 + w0**4/30.0_dp)
610 w5 = (w3 + 8.0_dp*w0**4/210.0_dp + w0**5/210.0_dp)
611
612 f = (1.0_dp - w*w1)*r1
613 d1 = -1.0_dp*(1.0_dp - w*w2)*r3
614 d2 = 3.0_dp*(1.0_dp - w*w3)*r5
615 d3 = -15.0_dp*(1.0_dp - w*w4)*r7
616 d4 = 105.0_dp*(1.0_dp - w*w5)*r9
617
618 kr = alpha*rr
619 kr2 = kr*kr
620 w0 = 1.0_dp - erfc(kr)
621 w1 = 2.0_dp*oorootpi*exp(-kr2)
622 w2 = w1*kr
623
624 f = f - w0*r1
625 d1 = d1 + (-w2 + w0)*r3
626 d2 = d2 + (w2*(3.0_dp + kr2*2.0_dp) - 3.0_dp*w0)*r5
627 d3 = d3 + (-w2*(15.0_dp + kr2*(10.0_dp + kr2*4.0_dp)) + 15.0_dp*w0)*r7
628 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 CALL build_d_tensor_gks(d1f, d2f, d3f, d4f, v=v, d1=d1, d2=d2, d3=d3, d4=d4)
631
632 DO ima = 1, (sepi%natorb*(sepi%natorb + 1))/2
633 DO imb = 1, (sepj%natorb*(sepj%natorb + 1))/2
634
635 w = m0a(ima)*m0b(imb)*f
636 DO k1 = 1, 3
637 w = w + (m1a(k1, ima)*m0b(imb) - m0a(ima)*m1b(k1, imb))*d1f(k1)
638 END DO
639 DO k2 = 1, 3
640 DO k1 = 1, 3
641 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 DO k3 = 1, 3
645 DO k2 = 1, 3
646 DO k1 = 1, 3
647 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 DO k4 = 1, 3
653 DO k3 = 1, 3
654 DO k2 = 1, 3
655 DO k1 = 1, 3
656 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 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 f = 0.0_dp
671 d1f = 0.0_dp
672 d2f = 0.0_dp
673 d3f = 0.0_dp
674 d4f = 0.0_dp
675
676 DO gpt = 1, pw_grid%ngpts_cut_local
677 lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
678 mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
679 np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
680
681 lp = lp + bds(1, 1)
682 mp = mp + bds(1, 2)
683 np = np + bds(1, 3)
684
685 IF (pw_grid%gsq(gpt) == 0.0_dp) cycle
686 kk(:) = pw_grid%g(:, gpt)
687 ff = 2.0_dp*fourpi*rho0(lp, mp, np)**2*pw_grid%vol/pw_grid%gsq(gpt)
688
689 kr = dot_product(kk, v)
690 cc = cos(kr)
691 ss = sin(kr)
692
693 f = f + cc*ff
694 DO k1 = 1, 3
695 d1f(k1) = d1f(k1) - kk(k1)*ss*ff
696 END DO
697 DO k2 = 1, 3
698 DO k1 = 1, 3
699 d2f(k1, k2) = d2f(k1, k2) - kk(k1)*kk(k2)*cc*ff
700 END DO
701 END DO
702 DO k3 = 1, 3
703 DO k2 = 1, 3
704 DO k1 = 1, 3
705 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 DO k4 = 1, 3
710 DO k3 = 1, 3
711 DO k2 = 1, 3
712 DO k1 = 1, 3
713 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 DO ima = 1, (sepi%natorb*(sepi%natorb + 1))/2
722 DO imb = 1, (sepj%natorb*(sepj%natorb + 1))/2
723
724 w = m0a(ima)*m0b(imb)*f
725 DO k1 = 1, 3
726 w = w + (m1a(k1, ima)*m0b(imb) - m0a(ima)*m1b(k1, imb))*d1f(k1)
727 END DO
728 DO k2 = 1, 3
729 DO k1 = 1, 3
730 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 DO k3 = 1, 3
734 DO k2 = 1, 3
735 DO k1 = 1, 3
736 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 DO k4 = 1, 3
742 DO k3 = 1, 3
743 DO k2 = 1, 3
744 DO k1 = 1, 3
745 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 coul(ima, imb) = coul(ima, imb) + w
752
753 END DO
754 END DO
755
756 DO ima = 1, (sepi%natorb*(sepi%natorb + 1))/2
757 DO imb = 1, (sepj%natorb*(sepj%natorb + 1))/2
758 w = -m0a(ima)*m0b(imb)*0.25_dp*fourpi/(pw_grid%vol*alpha**2)
759 coul(ima, imb) = coul(ima, imb) + w
760 END DO
761 END DO
762
763 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 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 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 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 alpha = se_int_control%ewald_gks%alpha
797 pw_pool => se_int_control%ewald_gks%pw_pool
798 dg => se_int_control%ewald_gks%dg
799 CALL dg_get(dg, dg_rho0=dg_rho0)
800 rho0 => dg_rho0%density%array
801 pw_grid => pw_pool%pw_grid
802 bds => pw_grid%bounds
803
804 CALL get_se_slater_multipole(sepi, m0a, m1a, m2a, acoula)
805 CALL get_se_slater_multipole(sepj, m0b, m1b, m2b, acoulb)
806
807 v(1) = rab(1)
808 v(2) = rab(2)
809 v(3) = rab(3)
810 rr = sqrt(dot_product(v, v))
811
812 a2 = 0.5_dp*(1.0_dp/acoula + 1.0_dp/acoulb)
813
814 r1 = 1.0_dp/rr
815 r2 = r1*r1
816 r3 = r2*r1
817 r5 = r3*r2
818 r7 = r5*r2
819 r9 = r7*r2
820 r11 = r9*r2
821
822 w0 = a2*rr
823 w = exp(-w0)
824 w1 = (1.0_dp + 0.5_dp*w0)
825 w2 = (w1 + 0.5_dp*w0 + 0.5_dp*w0**2)
826 w3 = (w2 + w0**3/6.0_dp)
827 w4 = (w3 + w0**4/30.0_dp)
828 w5 = (w3 + 8.0_dp*w0**4/210.0_dp + w0**5/210.0_dp)
829 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 f = (1.0_dp - w*w1)*r1
832 d1 = -1.0_dp*(1.0_dp - w*w2)*r3
833 d2 = 3.0_dp*(1.0_dp - w*w3)*r5
834 d3 = -15.0_dp*(1.0_dp - w*w4)*r7
835 d4 = 105.0_dp*(1.0_dp - w*w5)*r9
836 d5 = -945.0_dp*(1.0_dp - w*w6)*r11
837
838 kr = alpha*rr
839 kr2 = kr*kr
840 w0 = 1.0_dp - erfc(kr)
841 w1 = 2.0_dp*oorootpi*exp(-kr2)
842 w2 = w1*kr
843
844 f = f - w0*r1
845 d1 = d1 + (-w2 + w0)*r3
846 d2 = d2 + (w2*(3.0_dp + kr2*2.0_dp) - 3.0_dp*w0)*r5
847 d3 = d3 + (-w2*(15.0_dp + kr2*(10.0_dp + kr2*4.0_dp)) + 15.0_dp*w0)*r7
848 d4 = d4 + (w2*(105.0_dp + kr2*(70.0_dp + kr2*(28.0_dp + kr2*8.0_dp))) - 105.0_dp*w0)*r9
849 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 CALL build_d_tensor_gks(d1f, d2f, d3f, d4f, d5f, v, d1, d2, d3, d4, d5)
852
853 DO ima = 1, (sepi%natorb*(sepi%natorb + 1))/2
854 DO imb = 1, (sepj%natorb*(sepj%natorb + 1))/2
855
856 tmp = m0a(ima)*m0b(imb)
857 wv(1) = tmp*d1f(1)
858 wv(2) = tmp*d1f(2)
859 wv(3) = tmp*d1f(3)
860
861 DO k1 = 1, 3
862 tmp = m1a(k1, ima)*m0b(imb) - m0a(ima)*m1b(k1, imb)
863 wv(1) = wv(1) + tmp*d2f(1, k1)
864 wv(2) = wv(2) + tmp*d2f(2, k1)
865 wv(3) = wv(3) + tmp*d2f(3, k1)
866 END DO
867 DO k2 = 1, 3
868 DO k1 = 1, 3
869 tmp = m2a(k1, k2, ima)*m0b(imb) - m1a(k1, ima)*m1b(k2, imb) + m0a(ima)*m2b(k1, k2, imb)
870 wv(1) = wv(1) + tmp*d3f(1, k1, k2)
871 wv(2) = wv(2) + tmp*d3f(2, k1, k2)
872 wv(3) = wv(3) + tmp*d3f(3, k1, k2)
873 END DO
874 END DO
875 DO k3 = 1, 3
876 DO k2 = 1, 3
877 DO k1 = 1, 3
878 tmp = -m2a(k1, k2, ima)*m1b(k3, imb) + m1a(k1, ima)*m2b(k2, k3, imb)
879 wv(1) = wv(1) + tmp*d4f(1, k1, k2, k3)
880 wv(2) = wv(2) + tmp*d4f(2, k1, k2, k3)
881 wv(3) = wv(3) + tmp*d4f(3, k1, k2, k3)
882 END DO
883 END DO
884 END DO
885
886 DO k4 = 1, 3
887 DO k3 = 1, 3
888 DO k2 = 1, 3
889 DO k1 = 1, 3
890 tmp = m2a(k1, k2, ima)*m2b(k3, k4, imb)
891 wv(1) = wv(1) + tmp*d5f(1, k1, k2, k3, k4)
892 wv(2) = wv(2) + tmp*d5f(2, k1, k2, k3, k4)
893 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 dcoul(1, ima, imb) = wv(1)
900 dcoul(2, ima, imb) = wv(2)
901 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 f = 0.0_dp
910 d1f = 0.0_dp
911 d2f = 0.0_dp
912 d3f = 0.0_dp
913 d4f = 0.0_dp
914 d5f = 0.0_dp
915
916 DO gpt = 1, pw_grid%ngpts_cut_local
917 lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
918 mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
919 np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
920
921 lp = lp + bds(1, 1)
922 mp = mp + bds(1, 2)
923 np = np + bds(1, 3)
924
925 IF (pw_grid%gsq(gpt) == 0.0_dp) cycle
926 kk(:) = pw_grid%g(:, gpt)
927 ff = 2.0_dp*fourpi*rho0(lp, mp, np)**2*pw_grid%vol/pw_grid%gsq(gpt)
928
929 kr = dot_product(kk, v)
930 cc = cos(kr)
931 ss = sin(kr)
932
933 f = f + cc*ff
934 DO k1 = 1, 3
935 d1f(k1) = d1f(k1) - kk(k1)*ss*ff
936 END DO
937 DO k2 = 1, 3
938 DO k1 = 1, 3
939 d2f(k1, k2) = d2f(k1, k2) - kk(k1)*kk(k2)*cc*ff
940 END DO
941 END DO
942 DO k3 = 1, 3
943 DO k2 = 1, 3
944 DO k1 = 1, 3
945 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 DO k4 = 1, 3
950 DO k3 = 1, 3
951 DO k2 = 1, 3
952 DO k1 = 1, 3
953 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 DO k5 = 1, 3
959 DO k4 = 1, 3
960 DO k3 = 1, 3
961 DO k2 = 1, 3
962 DO k1 = 1, 3
963 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 DO ima = 1, (sepi%natorb*(sepi%natorb + 1))/2
972 DO imb = 1, (sepj%natorb*(sepj%natorb + 1))/2
973 tmp = m0a(ima)*m0b(imb)
974 wv(1) = tmp*d1f(1)
975 wv(2) = tmp*d1f(2)
976 wv(3) = tmp*d1f(3)
977 DO k1 = 1, 3
978 tmp = m1a(k1, ima)*m0b(imb) - m0a(ima)*m1b(k1, imb)
979 wv(1) = wv(1) + tmp*d2f(1, k1)
980 wv(2) = wv(2) + tmp*d2f(2, k1)
981 wv(3) = wv(3) + tmp*d2f(3, k1)
982 END DO
983 DO k2 = 1, 3
984 DO k1 = 1, 3
985 tmp = m2a(k1, k2, ima)*m0b(imb) - m1a(k1, ima)*m1b(k2, imb) + m0a(ima)*m2b(k1, k2, imb)
986 wv(1) = wv(1) + tmp*d3f(1, k1, k2)
987 wv(2) = wv(2) + tmp*d3f(2, k1, k2)
988 wv(3) = wv(3) + tmp*d3f(3, k1, k2)
989 END DO
990 END DO
991 DO k3 = 1, 3
992 DO k2 = 1, 3
993 DO k1 = 1, 3
994 tmp = -m2a(k1, k2, ima)*m1b(k3, imb) + m1a(k1, ima)*m2b(k2, k3, imb)
995 wv(1) = wv(1) + tmp*d4f(1, k1, k2, k3)
996 wv(2) = wv(2) + tmp*d4f(2, k1, k2, k3)
997 wv(3) = wv(3) + tmp*d4f(3, k1, k2, k3)
998 END DO
999 END DO
1000 END DO
1001
1002 DO k4 = 1, 3
1003 DO k3 = 1, 3
1004 DO k2 = 1, 3
1005 DO k1 = 1, 3
1006 tmp = m2a(k1, k2, ima)*m2b(k3, k4, imb)
1007 wv(1) = wv(1) + tmp*d5f(1, k1, k2, k3, k4)
1008 wv(2) = wv(2) + tmp*d5f(2, k1, k2, k3, k4)
1009 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 dcoul(1, ima, imb) = dcoul(1, ima, imb) + wv(1)
1016 dcoul(2, ima, imb) = dcoul(2, ima, imb) + wv(2)
1017 dcoul(3, ima, imb) = dcoul(3, ima, imb) + wv(3)
1018 END DO
1019 END DO
1020
1021 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 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 d1f = 0.0_dp
1053 d2f = 0.0_dp
1054 d3f = 0.0_dp
1055 d4f = 0.0_dp
1056 DO k1 = 1, 3
1057 d1f(k1) = d1f(k1) + v(k1)*d1
1058 END DO
1059 DO k1 = 1, 3
1060 DO k2 = 1, 3
1061 d2f(k2, k1) = d2f(k2, k1) + v(k1)*v(k2)*d2
1062 END DO
1063 d2f(k1, k1) = d2f(k1, k1) + d1
1064 END DO
1065 DO k1 = 1, 3
1066 DO k2 = 1, 3
1067 DO k3 = 1, 3
1068 d3f(k3, k2, k1) = d3f(k3, k2, k1) + v(k1)*v(k2)*v(k3)*d3
1069 END DO
1070 w = v(k1)*d2
1071 d3f(k1, k2, k2) = d3f(k1, k2, k2) + w
1072 d3f(k2, k1, k2) = d3f(k2, k1, k2) + w
1073 d3f(k2, k2, k1) = d3f(k2, k2, k1) + w
1074 END DO
1075 END DO
1076 DO k1 = 1, 3
1077 DO k2 = 1, 3
1078 DO k3 = 1, 3
1079 DO k4 = 1, 3
1080 d4f(k4, k3, k2, k1) = d4f(k4, k3, k2, k1) + v(k1)*v(k2)*v(k3)*v(k4)*d4
1081 END DO
1082 w = v(k1)*v(k2)*d3
1083 d4f(k1, k2, k3, k3) = d4f(k1, k2, k3, k3) + w
1084 d4f(k1, k3, k2, k3) = d4f(k1, k3, k2, k3) + w
1085 d4f(k3, k1, k2, k3) = d4f(k3, k1, k2, k3) + w
1086 d4f(k1, k3, k3, k2) = d4f(k1, k3, k3, k2) + w
1087 d4f(k3, k1, k3, k2) = d4f(k3, k1, k3, k2) + w
1088 d4f(k3, k3, k1, k2) = d4f(k3, k3, k1, k2) + w
1089 END DO
1090 d4f(k1, k1, k2, k2) = d4f(k1, k1, k2, k2) + d2
1091 d4f(k1, k2, k1, k2) = d4f(k1, k2, k1, k2) + d2
1092 d4f(k1, k2, k2, k1) = d4f(k1, k2, k2, k1) + d2
1093 END DO
1094 END DO
1095 IF (PRESENT(d5f) .AND. PRESENT(d5)) THEN
1096 d5f = 0.0_dp
1097
1098 DO k1 = 1, 3
1099 DO k2 = 1, 3
1100 DO k3 = 1, 3
1101 DO k4 = 1, 3
1102 DO k5 = 1, 3
1103 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 w = v(k1)*v(k2)*v(k3)*d4
1106 d5f(k1, k2, k3, k4, k4) = d5f(k1, k2, k3, k4, k4) + w
1107 d5f(k1, k2, k4, k3, k4) = d5f(k1, k2, k4, k3, k4) + w
1108 d5f(k1, k4, k2, k3, k4) = d5f(k1, k4, k2, k3, k4) + w
1109 d5f(k4, k1, k2, k3, k4) = d5f(k4, k1, k2, k3, k4) + w
1110 d5f(k1, k2, k4, k4, k3) = d5f(k1, k2, k4, k4, k3) + w
1111 d5f(k1, k4, k2, k4, k3) = d5f(k1, k4, k2, k4, k3) + w
1112 d5f(k4, k1, k2, k4, k3) = d5f(k4, k1, k2, k4, k3) + w
1113 d5f(k1, k4, k4, k2, k3) = d5f(k1, k4, k4, k2, k3) + w
1114 d5f(k4, k1, k4, k2, k3) = d5f(k4, k1, k4, k2, k3) + w
1115 d5f(k4, k4, k1, k2, k3) = d5f(k4, k4, k1, k2, k3) + w
1116 END DO
1117 w = v(k1)*d3
1118 d5f(k1, k2, k2, k3, k3) = d5f(k1, k2, k2, k3, k3) + w
1119 d5f(k1, k2, k3, k2, k3) = d5f(k1, k2, k3, k2, k3) + w
1120 d5f(k1, k2, k3, k3, k2) = d5f(k1, k2, k3, k3, k2) + w
1121 d5f(k2, k1, k2, k3, k3) = d5f(k2, k1, k2, k3, k3) + w
1122 d5f(k2, k1, k3, k2, k3) = d5f(k2, k1, k3, k2, k3) + w
1123 d5f(k2, k1, k3, k3, k2) = d5f(k2, k1, k3, k3, k2) + w
1124 d5f(k2, k2, k1, k3, k3) = d5f(k2, k2, k1, k3, k3) + w
1125 d5f(k2, k3, k1, k2, k3) = d5f(k2, k3, k1, k2, k3) + w
1126 d5f(k2, k3, k1, k3, k2) = d5f(k2, k3, k1, k3, k2) + w
1127 d5f(k2, k2, k3, k1, k3) = d5f(k2, k2, k3, k1, k3) + w
1128 d5f(k2, k3, k2, k1, k3) = d5f(k2, k3, k2, k1, k3) + w
1129 d5f(k2, k3, k3, k1, k2) = d5f(k2, k3, k3, k1, k2) + w
1130 d5f(k2, k2, k3, k3, k1) = d5f(k2, k2, k3, k3, k1) + w
1131 d5f(k2, k3, k2, k3, k1) = d5f(k2, k3, k2, k3, k1) + w
1132 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 END SUBROUTINE build_d_tensor_gks
1138
1139! **************************************************************************************************
1140!> \brief ...
1141!> \param sepi ...
1142!> \param Coul ...
1143!> \param se_int_control ...
1144! **************************************************************************************************
1145 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 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 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 alpha = se_int_control%ewald_gks%alpha
1165 pw_pool => se_int_control%ewald_gks%pw_pool
1166 dg => se_int_control%ewald_gks%dg
1167 CALL dg_get(dg, dg_rho0=dg_rho0)
1168 rho0 => dg_rho0%density%array
1169 pw_grid => pw_pool%pw_grid
1170 bds => pw_grid%bounds
1171
1172 CALL get_se_slater_multipole(sepi, m0a, m1a, m2a)
1173
1174 f = 0.0_dp
1175 d2f = 0.0_dp
1176 d4f = 0.0_dp
1177
1178 DO gpt = 1, pw_grid%ngpts_cut_local
1179 lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
1180 mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
1181 np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
1182
1183 lp = lp + bds(1, 1)
1184 mp = mp + bds(1, 2)
1185 np = np + bds(1, 3)
1186
1187 IF (pw_grid%gsq(gpt) == 0.0_dp) cycle
1188 kk(:) = pw_grid%g(:, gpt)
1189 ff = 2.0_dp*fourpi*rho0(lp, mp, np)**2*pw_grid%vol/pw_grid%gsq(gpt)
1190
1191 f = f + ff
1192 DO k2 = 1, 3
1193 DO k1 = 1, 3
1194 d2f(k1, k2) = d2f(k1, k2) - kk(k1)*kk(k2)*ff
1195 END DO
1196 END DO
1197 DO k4 = 1, 3
1198 DO k3 = 1, 3
1199 DO k2 = 1, 3
1200 DO k1 = 1, 3
1201 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 DO ima = 1, (sepi%natorb*(sepi%natorb + 1))/2
1210 DO imb = 1, (sepi%natorb*(sepi%natorb + 1))/2
1211
1212 w = m0a(ima)*m0a(imb)*f
1213 DO k2 = 1, 3
1214 DO k1 = 1, 3
1215 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 DO k4 = 1, 3
1220 DO k3 = 1, 3
1221 DO k2 = 1, 3
1222 DO k1 = 1, 3
1223 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 coul(ima, imb) = w
1230
1231 END DO
1232 END DO
1233
1234 DO ima = 1, (sepi%natorb*(sepi%natorb + 1))/2
1235 DO imb = 1, (sepi%natorb*(sepi%natorb + 1))/2
1236 w = -m0a(ima)*m0a(imb)*0.25_dp*fourpi/(pw_grid%vol*alpha**2)
1237 coul(ima, imb) = coul(ima, imb) + w
1238 END DO
1239 END DO
1240
1241 DO ima = 1, (sepi%natorb*(sepi%natorb + 1))/2
1242 DO imb = 1, (sepi%natorb*(sepi%natorb + 1))/2
1243
1244 w = m0a(ima)*m0a(imb)
1245 coul(ima, imb) = coul(ima, imb) - 2.0_dp*alpha*oorootpi*w
1246
1247 w = 0.0_dp
1248 DO k1 = 1, 3
1249 w = w + m1a(k1, ima)*m1a(k1, imb)
1250 w = w - m0a(ima)*m2a(k1, k1, imb)
1251 w = w - m2a(k1, k1, ima)*m0a(imb)
1252 END DO
1253 coul(ima, imb) = coul(ima, imb) - 4.0_dp*alpha**3*oorootpi*w/3.0_dp
1254
1255 w = 0.0_dp
1256 DO k2 = 1, 3
1257 DO k1 = 1, 3
1258 w = w + 2.0_dp*m2a(k1, k2, ima)*m2a(k1, k2, imb)
1259 w = w + m2a(k1, k1, ima)*m2a(k2, k2, imb)
1260 END DO
1261 END DO
1262 coul(ima, imb) = coul(ima, imb) - 8.0_dp*alpha**5*oorootpi*w/5.0_dp
1263 END DO
1264 END DO
1265 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 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 NULLIFY (mpole)
1286 size_1c_int = SIZE(sepi%w_mpole)
1287 DO jint = 1, size_1c_int
1288 mpole => sepi%w_mpole(jint)%mpole
1289 i = mpole%indi
1290 j = mpole%indj
1291 m0(indexb(i, j)) = -mpole%cs
1292
1293 m1(1, indexb(i, j)) = -mpole%ds(1)
1294 m1(2, indexb(i, j)) = -mpole%ds(2)
1295 m1(3, indexb(i, j)) = -mpole%ds(3)
1296
1297 m2(1, 1, indexb(i, j)) = -mpole%qq(1, 1)/3.0_dp
1298 m2(2, 1, indexb(i, j)) = -mpole%qq(2, 1)/3.0_dp
1299 m2(3, 1, indexb(i, j)) = -mpole%qq(3, 1)/3.0_dp
1300
1301 m2(1, 2, indexb(i, j)) = -mpole%qq(1, 2)/3.0_dp
1302 m2(2, 2, indexb(i, j)) = -mpole%qq(2, 2)/3.0_dp
1303 m2(3, 2, indexb(i, j)) = -mpole%qq(3, 2)/3.0_dp
1304
1305 m2(1, 3, indexb(i, j)) = -mpole%qq(1, 3)/3.0_dp
1306 m2(2, 3, indexb(i, j)) = -mpole%qq(2, 3)/3.0_dp
1307 m2(3, 3, indexb(i, j)) = -mpole%qq(3, 3)/3.0_dp
1308 END DO
1309 IF (PRESENT(acoul)) acoul = sepi%acoul
1310 END SUBROUTINE get_se_slater_multipole
1311
1312END MODULE semi_empirical_int_gks
dg_rho0_types
Definition dg_rho0_types.F:12
dg_types
Definition dg_types.F:12
dg_types::dg_get
subroutine, public dg_get(dg, dg_rho0)
Get the dg_type.
Definition dg_types.F:44
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::oorootpi
real(kind=dp), parameter, public oorootpi
Definition mathconstants.F:115
mathconstants::fourpi
real(kind=dp), parameter, public fourpi
Definition mathconstants.F:113
multipole_types
Multipole structure: for multipole (fixed and induced) in FF based MD.
Definition multipole_types.F:12
multipole_types::do_multipole_none
integer, parameter, public do_multipole_none
Definition multipole_types.F:31
pw_grid_types
Definition pw_grid_types.F:13
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
semi_empirical_int_arrays
Arrays of parameters used in the semi-empirical calculations \References Everywhere in this module TC...
Definition semi_empirical_int_arrays.F:17
semi_empirical_int_arrays::rij_threshold
real(kind=dp), parameter, public rij_threshold
Definition semi_empirical_int_arrays.F:27
semi_empirical_int_arrays::indexb
integer, dimension(9, 9), public indexb
Definition semi_empirical_int_arrays.F:71
semi_empirical_int_gks
Integral GKS scheme: The order of the integrals in makeCoul reflects the standard order by MOPAC.
Definition semi_empirical_int_gks.F:20
semi_empirical_int_gks::drotnuc_gks
subroutine, public drotnuc_gks(sepi, sepj, rij, de1b, de2a, se_int_control)
Computes the derivatives of the electron-nuclei integrals.
Definition semi_empirical_int_gks.F:156
semi_empirical_int_gks::corecore_gks
subroutine, public corecore_gks(sepi, sepj, rijv, enuc, denuc, se_int_control)
Computes nuclei-nuclei interactions.
Definition semi_empirical_int_gks.F:478
semi_empirical_int_gks::drotint_gks
subroutine, public drotint_gks(sepi, sepj, rij, dw, se_int_control)
Computes the derivatives of the electron-electron integrals.
Definition semi_empirical_int_gks.F:204
semi_empirical_int_gks::rotnuc_gks
subroutine, public rotnuc_gks(sepi, sepj, rij, e1b, e2a, se_int_control)
Computes the electron-nuclei integrals.
Definition semi_empirical_int_gks.F:59
semi_empirical_int_gks::rotint_gks
subroutine, public rotint_gks(sepi, sepj, rij, w, se_int_control)
Computes the electron-electron integrals.
Definition semi_empirical_int_gks.F:107
semi_empirical_mpole_types
Definition of the semi empirical multipole integral expansions types.
Definition semi_empirical_mpole_types.F:12
semi_empirical_types
Definition of the semi empirical parameter types.
Definition semi_empirical_types.F:12
semi_empirical_types::setup_se_int_control_type
subroutine, public setup_se_int_control_type(se_int_control, shortrange, do_ewald_r3, do_ewald_gks, integral_screening, max_multipole, pc_coulomb_int)
Setup the Semiempirical integral control type.
Definition semi_empirical_types.F:555
dg_rho0_types::dg_rho0_type
Type for Gaussian Densities type = type of gaussian (PME) grid = grid number gcc = Gaussian contracti...
Definition dg_rho0_types.F:39
dg_types::dg_type
Definition dg_types.F:23
pw_grid_types::pw_grid_type
Definition pw_grid_types.F:62
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
semi_empirical_mpole_types::semi_empirical_mpole_type
Semi-empirical integral multipole expansion type.
Definition semi_empirical_mpole_types.F:29
semi_empirical_types::se_int_control_type
Definition semi_empirical_types.F:136
semi_empirical_types::semi_empirical_type
Semi-empirical type.
Definition semi_empirical_types.F:57