(git:e7e05ae)
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 ! **************************************************************************************************
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  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 
1312 END 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_types::dg_type
Definition: dg_types.F:23