(git:ccc2433)
se_fock_matrix_integrals.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 Provides the low level routines to build both the exchange and the
10 !> Coulomb Fock matrices.. This routines support d-orbitals and should
11 !> be changed only if one knows exactly what he is doing..
12 !> \author Teodoro Laino [tlaino] (05.2009) - Split and module reorganization
13 !> \par History
14 !> Teodoro Laino (04.2008) [tlaino] - University of Zurich : d-orbitals
15 !> Teodoro Laino (09.2008) [tlaino] - University of Zurich : Speed-up
16 !> Teodoro Laino (09.2008) [tlaino] - University of Zurich : Periodic SE
17 ! **************************************************************************************************
18 MODULE se_fock_matrix_integrals
19 
20  USE kinds, ONLY: dp
21  USE semi_empirical_int_arrays, ONLY: se_orbital_pointer
22  USE semi_empirical_integrals, ONLY: drotint, &
23  drotnuc, &
24  rotint, &
25  rotnuc
26  USE semi_empirical_store_int_types, ONLY: semi_empirical_si_type
27  USE semi_empirical_types, ONLY: se_int_control_type, &
28  se_taper_type, &
29  semi_empirical_type
30 #include "./base/base_uses.f90"
31 
32 
33  IMPLICIT NONE
34  PRIVATE
35 
36  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'se_fock_matrix_integrals'
37  LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
38 
39  PUBLIC :: fock2_1el, dfock2_1el, fock1_2el, fock2_1el_ew, fock2c_ew, &
40  fock2c, dfock2c, fock2e, dfock2e, fock2_1el_r3, dfock2_1el_r3, &
41  fock2c_r3, dfock2c_r3, se_coulomb_ij_interaction
42 
43 CONTAINS
44 
45 ! **************************************************************************************************
46 !> \brief Construction of 2-center 1-electron Fock Matrix
47 !> \param sepi ...
48 !> \param sepj ...
49 !> \param rij ...
50 !> \param ksi_block DIMENSION(sepi%natorb, sepi%natorb)
51 !> \param ksj_block DIMENSION(sepi%natorb, sepi%natorb)
52 !> \param pi_block ...
53 !> \param pj_block ...
54 !> \param ecore ...
55 !> \param itype ...
56 !> \param anag ...
57 !> \param se_int_control ...
58 !> \param se_taper ...
59 !> \param store_int_env ...
60 !> \date 04.2008 [tlaino]
61 !> \author Teodoro Laino [tlaino] - University of Zurich
62 ! **************************************************************************************************
63  SUBROUTINE fock2_1el(sepi, sepj, rij, ksi_block, ksj_block, pi_block, pj_block, &
64  ecore, itype, anag, se_int_control, se_taper, store_int_env)
65  TYPE(semi_empirical_type), POINTER :: sepi, sepj
66  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rij
67  REAL(kind=dp), DIMENSION(:, :), POINTER :: ksi_block, ksj_block
68  REAL(kind=dp), &
69  DIMENSION(sepi%natorb, sepi%natorb), INTENT(IN) :: pi_block
70  REAL(kind=dp), &
71  DIMENSION(sepj%natorb, sepj%natorb), INTENT(IN) :: pj_block
72  REAL(kind=dp), DIMENSION(2), INTENT(INOUT) :: ecore
73  INTEGER, INTENT(IN) :: itype
74  LOGICAL, INTENT(IN) :: anag
75  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
76  TYPE(se_taper_type), POINTER :: se_taper
77  TYPE(semi_empirical_si_type), POINTER :: store_int_env
78 
79  INTEGER :: i1, i1l, i2, j1, j1l
80  REAL(kind=dp), DIMENSION(45) :: e1b, e2a
81 
82 ! Compute integrals
83 
84  CALL rotnuc(sepi, sepj, rij, e1b=e1b, e2a=e2a, itype=itype, anag=anag, &
85  se_int_control=se_int_control, se_taper=se_taper, store_int_env=store_int_env)
86  !
87  ! Add the electron-nuclear attraction term for atom sepi
88  !
89  i2 = 0
90  DO i1l = 1, sepi%natorb
91  i1 = se_orbital_pointer(i1l)
92  DO j1l = 1, i1l - 1
93  j1 = se_orbital_pointer(j1l)
94  i2 = i2 + 1
95  ksi_block(i1, j1) = ksi_block(i1, j1) + e1b(i2)
96  ksi_block(j1, i1) = ksi_block(i1, j1)
97  ecore(1) = ecore(1) + 2.0_dp*e1b(i2)*pi_block(i1, j1)
98  END DO
99  j1 = se_orbital_pointer(j1l)
100  i2 = i2 + 1
101  ksi_block(i1, j1) = ksi_block(i1, j1) + e1b(i2)
102  ecore(1) = ecore(1) + e1b(i2)*pi_block(i1, j1)
103  END DO
104  !
105  ! Add the electron-nuclear attraction term for atom sepj
106  !
107  i2 = 0
108  DO i1l = 1, sepj%natorb
109  i1 = se_orbital_pointer(i1l)
110  DO j1l = 1, i1l - 1
111  j1 = se_orbital_pointer(j1l)
112  i2 = i2 + 1
113  ksj_block(i1, j1) = ksj_block(i1, j1) + e2a(i2)
114  ksj_block(j1, i1) = ksj_block(i1, j1)
115  ecore(2) = ecore(2) + 2.0_dp*e2a(i2)*pj_block(i1, j1)
116  END DO
117  j1 = se_orbital_pointer(j1l)
118  i2 = i2 + 1
119  ksj_block(i1, j1) = ksj_block(i1, j1) + e2a(i2)
120  ecore(2) = ecore(2) + e2a(i2)*pj_block(i1, j1)
121  END DO
122 
123  END SUBROUTINE fock2_1el
124 
125 ! **************************************************************************************************
126 !> \brief Derivatives of 2-center 1-electron Fock Matrix
127 !> \param sepi ...
128 !> \param sepj ...
129 !> \param rij ...
130 !> \param pi_block ...
131 !> \param pj_block ...
132 !> \param itype ...
133 !> \param anag ...
134 !> \param se_int_control ...
135 !> \param se_taper ...
136 !> \param force ...
137 !> \param delta ...
138 !> \date 04.2008 [tlaino]
139 !> \author Teodoro Laino [tlaino] - University of Zurich
140 ! **************************************************************************************************
141  SUBROUTINE dfock2_1el(sepi, sepj, rij, pi_block, pj_block, itype, anag, &
142  se_int_control, se_taper, force, delta)
143  TYPE(semi_empirical_type), POINTER :: sepi, sepj
144  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rij
145  REAL(kind=dp), &
146  DIMENSION(sepi%natorb, sepi%natorb), INTENT(IN) :: pi_block
147  REAL(kind=dp), &
148  DIMENSION(sepj%natorb, sepj%natorb), INTENT(IN) :: pj_block
149  INTEGER, INTENT(IN) :: itype
150  LOGICAL, INTENT(IN) :: anag
151  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
152  TYPE(se_taper_type), POINTER :: se_taper
153  REAL(kind=dp), DIMENSION(3), INTENT(INOUT) :: force
154  REAL(kind=dp), INTENT(IN) :: delta
155 
156  INTEGER :: i1, i1l, i2, j1, j1l
157  REAL(kind=dp) :: tmp
158  REAL(kind=dp), DIMENSION(3, 45) :: de1b, de2a
159 
160 ! Compute integrals
161 
162  CALL drotnuc(sepi, sepj, rij, de1b=de1b, de2a=de2a, itype=itype, anag=anag, &
163  se_int_control=se_int_control, se_taper=se_taper, delta=delta)
164  !
165  ! Add the electron-nuclear attraction term for atom sepi
166  !
167  i2 = 0
168  DO i1l = 1, sepi%natorb
169  i1 = se_orbital_pointer(i1l)
170  DO j1l = 1, i1l - 1
171  j1 = se_orbital_pointer(j1l)
172  i2 = i2 + 1
173  tmp = 2.0_dp*pi_block(i1, j1)
174  force(1) = force(1) + de1b(1, i2)*tmp
175  force(2) = force(2) + de1b(2, i2)*tmp
176  force(3) = force(3) + de1b(3, i2)*tmp
177  END DO
178  j1 = se_orbital_pointer(j1l)
179  i2 = i2 + 1
180  force(1) = force(1) + de1b(1, i2)*pi_block(i1, j1)
181  force(2) = force(2) + de1b(2, i2)*pi_block(i1, j1)
182  force(3) = force(3) + de1b(3, i2)*pi_block(i1, j1)
183  END DO
184  !
185  ! Add the electron-nuclear attraction term for atom sepj
186  !
187  i2 = 0
188  DO i1l = 1, sepj%natorb
189  i1 = se_orbital_pointer(i1l)
190  DO j1l = 1, i1l - 1
191  j1 = se_orbital_pointer(j1l)
192  i2 = i2 + 1
193  tmp = 2.0_dp*pj_block(i1, j1)
194  force(1) = force(1) + de2a(1, i2)*tmp
195  force(2) = force(2) + de2a(2, i2)*tmp
196  force(3) = force(3) + de2a(3, i2)*tmp
197  END DO
198  j1 = se_orbital_pointer(j1l)
199  i2 = i2 + 1
200  force(1) = force(1) + de2a(1, i2)*pj_block(i1, j1)
201  force(2) = force(2) + de2a(2, i2)*pj_block(i1, j1)
202  force(3) = force(3) + de2a(3, i2)*pj_block(i1, j1)
203  END DO
204 
205  END SUBROUTINE dfock2_1el
206 
207 ! **************************************************************************************************
208 !> \brief Construction of 1-center 2-electron Fock Matrix
209 !> \param sep ...
210 !> \param p_tot ...
211 !> \param p_mat ...
212 !> \param f_mat DIMENSION(sep%natorb, sep%natorb)
213 !> \param factor ...
214 !> \date 04.2008 [tlaino]
215 !> \author Teodoro Laino [tlaino] - University of Zurich
216 ! **************************************************************************************************
217  SUBROUTINE fock1_2el(sep, p_tot, p_mat, f_mat, factor)
218  TYPE(semi_empirical_type), POINTER :: sep
219  REAL(kind=dp), DIMENSION(45, 45), INTENT(IN) :: p_tot
220  REAL(kind=dp), DIMENSION(sep%natorb, sep%natorb), &
221  INTENT(IN) :: p_mat
222  REAL(kind=dp), DIMENSION(:, :), POINTER :: f_mat
223  REAL(kind=dp), INTENT(IN) :: factor
224 
225  INTEGER :: i, ijw, ikw, il, im, j, jl, jlw, jm, k, &
226  kl, klw, l, ll
227  REAL(kind=dp) :: sum
228 
229 ! One-center coulomb and exchange terms for semiempirical_type sep
230 !
231 ! F(i,j)=F(i,j)+sum(k,l)((PA(k,l)+PB(k,l))*<i,j|k,l>
232 ! -(PA(k,l) )*<i,k|j,l>), k,l on type sep.
233 !
234 
235  DO il = 1, sep%natorb
236  i = se_orbital_pointer(il)
237  DO jl = 1, il
238  j = se_orbital_pointer(jl)
239 
240  ! `J' Address IJ in W
241  ijw = (il*(il - 1))/2 + jl
242  sum = 0.0_dp
243  DO kl = 1, sep%natorb
244  k = se_orbital_pointer(kl)
245  DO ll = 1, sep%natorb
246  l = se_orbital_pointer(ll)
247 
248  ! `J' Address KL in W
249  im = max(kl, ll)
250  jm = min(kl, ll)
251  klw = (im*(im - 1))/2 + jm
252 
253  ! `K' Address IK in W
254  im = max(kl, jl)
255  jm = min(kl, jl)
256  ikw = (im*(im - 1))/2 + jm
257 
258  ! `K' Address JL in W
259  im = max(ll, il)
260  jm = min(ll, il)
261  jlw = (im*(im - 1))/2 + jm
262 
263  sum = sum + p_tot(k, l)*sep%w(ijw, klw) - p_mat(k, l)*sep%w(ikw, jlw)
264  END DO
265  END DO
266  f_mat(i, j) = f_mat(i, j) + factor*sum
267  f_mat(j, i) = f_mat(i, j)
268  END DO
269  END DO
270  END SUBROUTINE fock1_2el
271 
272 ! **************************************************************************************************
273 !> \brief Construction of 2-center 1-electron Fock Matrix (Ewald self term)
274 !> \param sep ...
275 !> \param rij ...
276 !> \param ks_block DIMENSION(sep%natorb, sep%natorb)
277 !> \param p_block ...
278 !> \param ecore ...
279 !> \param itype ...
280 !> \param anag ...
281 !> \param se_int_control ...
282 !> \param se_taper ...
283 !> \param store_int_env ...
284 !> \date 04.2009 [jgh]
285 !> \author jgh - University of Zurich
286 ! **************************************************************************************************
287  SUBROUTINE fock2_1el_ew(sep, rij, ks_block, p_block, ecore, itype, anag, &
288  se_int_control, se_taper, store_int_env)
289  TYPE(semi_empirical_type), POINTER :: sep
290  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rij
291  REAL(kind=dp), DIMENSION(:, :), POINTER :: ks_block
292  REAL(kind=dp), DIMENSION(sep%natorb, sep%natorb), &
293  INTENT(IN) :: p_block
294  REAL(kind=dp), INTENT(INOUT) :: ecore
295  INTEGER, INTENT(IN) :: itype
296  LOGICAL, INTENT(IN) :: anag
297  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
298  TYPE(se_taper_type), POINTER :: se_taper
299  TYPE(semi_empirical_si_type), POINTER :: store_int_env
300 
301  INTEGER :: i1, i1l, i2, j1, j1l, n
302  REAL(kind=dp), DIMENSION(45) :: e1b, e2a
303 
304 ! Compute integrals
305 
306  CALL rotnuc(sep, sep, rij, e1b=e1b, e2a=e2a, itype=itype, anag=anag, &
307  se_int_control=se_int_control, se_taper=se_taper, store_int_env=store_int_env)
308  !
309  ! Add the electron-nuclear attraction term for atom sep
310  ! e1b == e2a
311  !
312  n = (sep%natorb*(sep%natorb + 1))/2
313  i2 = 0
314  DO i1l = 1, sep%natorb
315  i1 = se_orbital_pointer(i1l)
316  DO j1l = 1, i1l - 1
317  j1 = se_orbital_pointer(j1l)
318  i2 = i2 + 1
319  ks_block(i1, j1) = ks_block(i1, j1) + e1b(i2)
320  ks_block(j1, i1) = ks_block(i1, j1)
321  ecore = ecore + 2._dp*e1b(i2)*p_block(i1, j1)
322  END DO
323  ! i1L == j1L
324  i2 = i2 + 1
325  ks_block(i1, i1) = ks_block(i1, i1) + e1b(i2)
326  ecore = ecore + e1b(i2)*p_block(i1, i1)
327  END DO
328 
329  END SUBROUTINE fock2_1el_ew
330 
331 ! **************************************************************************************************
332 !> \brief Construction of 2-center Fock Matrix - Coulomb Self Terms (Ewald)
333 !> \param sep ...
334 !> \param rij ...
335 !> \param p_tot ...
336 !> \param f_mat DIMENSION(sep%natorb, sep%natorb)
337 !> \param factor ...
338 !> \param anag ...
339 !> \param se_int_control ...
340 !> \param se_taper ...
341 !> \param store_int_env ...
342 !> \date 04.2009 [jgh]
343 !> \author jgh - University of Zurich
344 ! **************************************************************************************************
345  SUBROUTINE fock2c_ew(sep, rij, p_tot, f_mat, factor, anag, se_int_control, &
346  se_taper, store_int_env)
347  TYPE(semi_empirical_type), POINTER :: sep
348  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rij
349  REAL(kind=dp), DIMENSION(45, 45), INTENT(IN) :: p_tot
350  REAL(kind=dp), DIMENSION(:, :), POINTER :: f_mat
351  REAL(kind=dp), INTENT(IN) :: factor
352  LOGICAL, INTENT(IN) :: anag
353  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
354  TYPE(se_taper_type), POINTER :: se_taper
355  TYPE(semi_empirical_si_type), POINTER :: store_int_env
356 
357  INTEGER :: i, il, j, jl, k, kl, kr, l, ll, natorb
358  REAL(kind=dp) :: a, aa, bb
359  REAL(kind=dp), DIMENSION(2025) :: w
360 
361 ! Evaluate integrals
362 
363  CALL rotint(sep, sep, rij, w, anag=anag, se_int_control=se_int_control, &
364  se_taper=se_taper, store_int_env=store_int_env)
365  kr = 0
366  natorb = sep%natorb
367  DO il = 1, natorb
368  i = se_orbital_pointer(il)
369  aa = 2.0_dp
370  DO jl = 1, il
371  j = se_orbital_pointer(jl)
372  IF (i == j) THEN
373  aa = 1.0_dp
374  END IF
375  DO kl = 1, natorb
376  k = se_orbital_pointer(kl)
377  bb = 2.0_dp
378  DO ll = 1, kl
379  l = se_orbital_pointer(ll)
380  IF (k == l) THEN
381  bb = 1.0_dp
382  END IF
383  kr = kr + 1
384  a = 0.5_dp*w(kr)*factor
385  ! Coulomb
386  f_mat(i, j) = f_mat(i, j) + bb*a*p_tot(k, l)
387  f_mat(k, l) = f_mat(k, l) + aa*a*p_tot(i, j)
388  f_mat(j, i) = f_mat(i, j)
389  f_mat(l, k) = f_mat(k, l)
390  END DO
391  END DO
392  END DO
393  END DO
394 
395  END SUBROUTINE fock2c_ew
396 
397 ! **************************************************************************************************
398 !> \brief Construction of 2-center Fock Matrix - Coulomb Terms
399 !> \param sepi ...
400 !> \param sepj ...
401 !> \param rij ...
402 !> \param switch ...
403 !> \param pi_tot ...
404 !> \param fi_mat DIMENSION(sepi%natorb, sepi%natorb)
405 !> \param pj_tot DIMENSION(sepj%natorb, sepj%natorb)
406 !> \param fj_mat ...
407 !> \param factor ...
408 !> \param anag ...
409 !> \param se_int_control ...
410 !> \param se_taper ...
411 !> \param store_int_env ...
412 !> \date 04.2008 [tlaino]
413 !> \author Teodoro Laino [tlaino] - University of Zurich
414 ! **************************************************************************************************
415  SUBROUTINE fock2c(sepi, sepj, rij, switch, pi_tot, fi_mat, pj_tot, fj_mat, &
416  factor, anag, se_int_control, se_taper, store_int_env)
417  TYPE(semi_empirical_type), POINTER :: sepi, sepj
418  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rij
419  LOGICAL, INTENT(IN) :: switch
420  REAL(kind=dp), DIMENSION(45, 45), INTENT(IN) :: pi_tot
421  REAL(kind=dp), DIMENSION(:, :), POINTER :: fi_mat
422  REAL(kind=dp), DIMENSION(45, 45), INTENT(IN) :: pj_tot
423  REAL(kind=dp), DIMENSION(:, :), POINTER :: fj_mat
424  REAL(kind=dp), INTENT(IN) :: factor
425  LOGICAL, INTENT(IN) :: anag
426  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
427  TYPE(se_taper_type), POINTER :: se_taper
428  TYPE(semi_empirical_si_type), POINTER :: store_int_env
429 
430  INTEGER :: i, il, j, jl, k, kl, kr, l, ll, natorb(2)
431  REAL(kind=dp) :: a, aa, bb, irij(3)
432  REAL(kind=dp), DIMENSION(2025) :: w
433 
434 ! Evaluate integrals
435 
436  IF (.NOT. switch) THEN
437  CALL rotint(sepi, sepj, rij, w, anag=anag, se_int_control=se_int_control, &
438  se_taper=se_taper, store_int_env=store_int_env)
439  ELSE
440  irij = -rij
441  CALL rotint(sepj, sepi, irij, w, anag=anag, se_int_control=se_int_control, &
442  se_taper=se_taper, store_int_env=store_int_env)
443  END IF
444  kr = 0
445  natorb(1) = sepi%natorb
446  natorb(2) = sepj%natorb
447  IF (switch) THEN
448  natorb(1) = sepj%natorb
449  natorb(2) = sepi%natorb
450  END IF
451  DO il = 1, natorb(1)
452  i = se_orbital_pointer(il)
453  aa = 2.0_dp
454  DO jl = 1, il
455  j = se_orbital_pointer(jl)
456  IF (i == j) THEN
457  aa = 1.0_dp
458  END IF
459  DO kl = 1, natorb(2)
460  k = se_orbital_pointer(kl)
461  bb = 2.0_dp
462  DO ll = 1, kl
463  l = se_orbital_pointer(ll)
464  IF (k == l) THEN
465  bb = 1.0_dp
466  END IF
467  kr = kr + 1
468  a = w(kr)*factor
469  ! Coulomb
470  IF (.NOT. switch) THEN
471  fi_mat(i, j) = fi_mat(i, j) + bb*a*pj_tot(k, l)
472  fj_mat(k, l) = fj_mat(k, l) + aa*a*pi_tot(i, j)
473  fi_mat(j, i) = fi_mat(i, j)
474  fj_mat(l, k) = fj_mat(k, l)
475  ELSE
476  fj_mat(i, j) = fj_mat(i, j) + bb*a*pi_tot(k, l)
477  fi_mat(k, l) = fi_mat(k, l) + aa*a*pj_tot(i, j)
478  fj_mat(j, i) = fj_mat(i, j)
479  fi_mat(l, k) = fi_mat(k, l)
480  END IF
481  END DO
482  END DO
483  END DO
484  END DO
485 
486  END SUBROUTINE fock2c
487 
488 ! **************************************************************************************************
489 !> \brief Derivatives of 2-center Fock Matrix - Coulomb Terms
490 !> \param sepi ...
491 !> \param sepj ...
492 !> \param rij ...
493 !> \param switch ...
494 !> \param pi_tot ...
495 !> \param pj_tot ...
496 !> \param factor ...
497 !> \param anag ...
498 !> \param se_int_control ...
499 !> \param se_taper ...
500 !> \param force ...
501 !> \param delta ...
502 !> \date 04.2008 [tlaino]
503 !> \author Teodoro Laino [tlaino] - University of Zurich
504 ! **************************************************************************************************
505  SUBROUTINE dfock2c(sepi, sepj, rij, switch, pi_tot, pj_tot, factor, anag, &
506  se_int_control, se_taper, force, delta)
507  TYPE(semi_empirical_type), POINTER :: sepi, sepj
508  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rij
509  LOGICAL, INTENT(IN) :: switch
510  REAL(kind=dp), DIMENSION(45, 45), INTENT(IN) :: pi_tot, pj_tot
511  REAL(kind=dp), INTENT(IN) :: factor
512  LOGICAL, INTENT(IN) :: anag
513  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
514  TYPE(se_taper_type), POINTER :: se_taper
515  REAL(kind=dp), DIMENSION(3), INTENT(INOUT) :: force
516  REAL(kind=dp), INTENT(IN) :: delta
517 
518  INTEGER :: i, il, j, jl, k, kl, kr, l, ll, natorb(2)
519  REAL(kind=dp) :: aa, bb, tmp
520  REAL(kind=dp), DIMENSION(3) :: a, irij
521  REAL(kind=dp), DIMENSION(3, 2025) :: dw
522 
523 ! Evaluate integrals' derivatives
524 
525  IF (.NOT. switch) THEN
526  CALL drotint(sepi, sepj, rij, dw, delta, anag=anag, se_int_control=se_int_control, &
527  se_taper=se_taper)
528  ELSE
529  irij = -rij
530  CALL drotint(sepj, sepi, irij, dw, delta, anag=anag, se_int_control=se_int_control, &
531  se_taper=se_taper)
532  END IF
533 
534  kr = 0
535  natorb(1) = sepi%natorb
536  natorb(2) = sepj%natorb
537  IF (switch) THEN
538  natorb(1) = sepj%natorb
539  natorb(2) = sepi%natorb
540  END IF
541  DO il = 1, natorb(1)
542  i = se_orbital_pointer(il)
543  aa = 2.0_dp
544  DO jl = 1, il
545  j = se_orbital_pointer(jl)
546  IF (i == j) THEN
547  aa = 1.0_dp
548  END IF
549  DO kl = 1, natorb(2)
550  k = se_orbital_pointer(kl)
551  bb = 2.0_dp
552  DO ll = 1, kl
553  l = se_orbital_pointer(ll)
554  IF (k == l) THEN
555  bb = 1.0_dp
556  END IF
557  kr = kr + 1
558  a(1) = dw(1, kr)*factor
559  a(2) = dw(2, kr)*factor
560  a(3) = dw(3, kr)*factor
561  ! Coulomb
562  IF (.NOT. switch) THEN
563  tmp = bb*aa*pj_tot(k, l)*pi_tot(i, j)
564  ELSE
565  tmp = bb*aa*pi_tot(k, l)*pj_tot(i, j)
566  END IF
567  force(1) = force(1) + a(1)*tmp
568  force(2) = force(2) + a(2)*tmp
569  force(3) = force(3) + a(3)*tmp
570  END DO
571  END DO
572  END DO
573  END DO
574  END SUBROUTINE dfock2c
575 
576 ! **************************************************************************************************
577 !> \brief Construction of 2-center Fock Matrix - General Driver
578 !> \param sepi ...
579 !> \param sepj ...
580 !> \param rij ...
581 !> \param switch ...
582 !> \param isize ...
583 !> \param pi_mat ...
584 !> \param fi_mat DIMENSION(isize(1), isize(2))
585 !> \param factor ...
586 !> \param anag ...
587 !> \param se_int_control ...
588 !> \param se_taper ...
589 !> \param store_int_env ...
590 !> \date 04.2008 [tlaino]
591 !> \author Teodoro Laino [tlaino] - University of Zurich
592 ! **************************************************************************************************
593  SUBROUTINE fock2e(sepi, sepj, rij, switch, isize, pi_mat, fi_mat, factor, &
594  anag, se_int_control, se_taper, store_int_env)
595  TYPE(semi_empirical_type), POINTER :: sepi, sepj
596  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rij
597  LOGICAL, INTENT(IN) :: switch
598  INTEGER, DIMENSION(2), INTENT(IN) :: isize
599  REAL(kind=dp), DIMENSION(isize(1), isize(2)), &
600  INTENT(IN) :: pi_mat
601  REAL(kind=dp), DIMENSION(:, :), POINTER :: fi_mat
602  REAL(kind=dp), INTENT(IN) :: factor
603  LOGICAL, INTENT(IN) :: anag
604  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
605  TYPE(se_taper_type), POINTER :: se_taper
606  TYPE(semi_empirical_si_type), POINTER :: store_int_env
607 
608  INTEGER :: i, il, j, jl, k, kl, kr, l, ll, natorb(2)
609  REAL(kind=dp) :: a, aa, bb, irij(3)
610  REAL(kind=dp), DIMENSION(2025) :: w
611 
612 ! Evaluate integrals
613 
614  IF (.NOT. switch) THEN
615  CALL rotint(sepi, sepj, rij, w, anag=anag, se_int_control=se_int_control, &
616  se_taper=se_taper, store_int_env=store_int_env)
617  ELSE
618  irij = -rij
619  CALL rotint(sepj, sepi, irij, w, anag=anag, se_int_control=se_int_control, &
620  se_taper=se_taper, store_int_env=store_int_env)
621  END IF
622  kr = 0
623  natorb(1) = sepi%natorb
624  natorb(2) = sepj%natorb
625  IF (switch) THEN
626  natorb(1) = sepj%natorb
627  natorb(2) = sepi%natorb
628  END IF
629  DO il = 1, natorb(1)
630  i = se_orbital_pointer(il)
631  aa = 2.0_dp
632  DO jl = 1, il
633  j = se_orbital_pointer(jl)
634  IF (i == j) THEN
635  aa = 1.0_dp
636  END IF
637  DO kl = 1, natorb(2)
638  k = se_orbital_pointer(kl)
639  bb = 2.0_dp
640  DO ll = 1, kl
641  l = se_orbital_pointer(ll)
642  IF (k == l) THEN
643  bb = 1.0_dp
644  END IF
645  kr = kr + 1
646  a = w(kr)*factor
647  ! Exchange
648  a = a*aa*bb*0.25_dp
649  fi_mat(i, k) = fi_mat(i, k) - a*pi_mat(j, l)
650  fi_mat(i, l) = fi_mat(i, l) - a*pi_mat(j, k)
651  fi_mat(j, k) = fi_mat(j, k) - a*pi_mat(i, l)
652  fi_mat(j, l) = fi_mat(j, l) - a*pi_mat(i, k)
653  END DO
654  END DO
655  END DO
656  END DO
657 
658  END SUBROUTINE fock2e
659 
660 ! **************************************************************************************************
661 !> \brief Derivatives of 2-center Fock Matrix - General Driver
662 !> \param sepi ...
663 !> \param sepj ...
664 !> \param rij ...
665 !> \param switch ...
666 !> \param isize ...
667 !> \param pi_mat ...
668 !> \param factor ...
669 !> \param anag ...
670 !> \param se_int_control ...
671 !> \param se_taper ...
672 !> \param force ...
673 !> \param delta ...
674 !> \date 04.2008 [tlaino]
675 !> \author Teodoro Laino [tlaino] - University of Zurich
676 ! **************************************************************************************************
677  SUBROUTINE dfock2e(sepi, sepj, rij, switch, isize, pi_mat, factor, anag, &
678  se_int_control, se_taper, force, delta)
679  TYPE(semi_empirical_type), POINTER :: sepi, sepj
680  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rij
681  LOGICAL, INTENT(IN) :: switch
682  INTEGER, DIMENSION(2), INTENT(IN) :: isize
683  REAL(kind=dp), DIMENSION(isize(1), isize(2)), &
684  INTENT(IN) :: pi_mat
685  REAL(kind=dp), INTENT(IN) :: factor
686  LOGICAL, INTENT(IN) :: anag
687  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
688  TYPE(se_taper_type), POINTER :: se_taper
689  REAL(kind=dp), DIMENSION(3), INTENT(INOUT) :: force
690  REAL(kind=dp), INTENT(IN) :: delta
691 
692  INTEGER :: i, il, j, jl, k, kl, kr, l, ll, natorb(2)
693  REAL(kind=dp) :: aa, bb, tmp, tmp1, tmp2, tmp3, tmp4
694  REAL(kind=dp), DIMENSION(3) :: a, irij
695  REAL(kind=dp), DIMENSION(3, 2025) :: dw
696 
697 ! Evaluate integrals' derivatives
698 
699  IF (.NOT. switch) THEN
700  CALL drotint(sepi, sepj, rij, dw, delta, anag=anag, se_int_control=se_int_control, &
701  se_taper=se_taper)
702  ELSE
703  irij = -rij
704  CALL drotint(sepj, sepi, irij, dw, delta, anag=anag, se_int_control=se_int_control, &
705  se_taper=se_taper)
706  END IF
707 
708  kr = 0
709  natorb(1) = sepi%natorb
710  natorb(2) = sepj%natorb
711  IF (switch) THEN
712  natorb(1) = sepj%natorb
713  natorb(2) = sepi%natorb
714  END IF
715  DO il = 1, natorb(1)
716  i = se_orbital_pointer(il)
717  aa = 2.0_dp
718  DO jl = 1, il
719  j = se_orbital_pointer(jl)
720  IF (i == j) THEN
721  aa = 1.0_dp
722  END IF
723  DO kl = 1, natorb(2)
724  k = se_orbital_pointer(kl)
725  bb = 2.0_dp
726  DO ll = 1, kl
727  l = se_orbital_pointer(ll)
728  IF (k == l) THEN
729  bb = 1.0_dp
730  END IF
731  kr = kr + 1
732  tmp = factor*aa*bb*0.25_dp
733  a(1) = dw(1, kr)*tmp
734  a(2) = dw(2, kr)*tmp
735  a(3) = dw(3, kr)*tmp
736  ! Exchange
737  tmp1 = pi_mat(j, l)*pi_mat(i, k)
738  tmp2 = pi_mat(j, k)*pi_mat(i, l)
739  tmp3 = pi_mat(i, l)*pi_mat(j, k)
740  tmp4 = pi_mat(i, k)*pi_mat(j, l)
741 
742  force(1) = force(1) - a(1)*tmp1
743  force(1) = force(1) - a(1)*tmp2
744  force(1) = force(1) - a(1)*tmp3
745  force(1) = force(1) - a(1)*tmp4
746 
747  force(2) = force(2) - a(2)*tmp1
748  force(2) = force(2) - a(2)*tmp2
749  force(2) = force(2) - a(2)*tmp3
750  force(2) = force(2) - a(2)*tmp4
751 
752  force(3) = force(3) - a(3)*tmp1
753  force(3) = force(3) - a(3)*tmp2
754  force(3) = force(3) - a(3)*tmp3
755  force(3) = force(3) - a(3)*tmp4
756  END DO
757  END DO
758  END DO
759  END DO
760  END SUBROUTINE dfock2e
761 
762 ! **************************************************************************************************
763 !> \brief Construction of 2-center 1-electron Fock Matrix for the residual
764 !> (1/R^3) integral part
765 !> \param sepi ...
766 !> \param sepj ...
767 !> \param ksi_block DIMENSION(sepi%natorb, sepi%natorb)
768 !> \param ksj_block DIMENSION(sepj%natorb, sepj%natorb)
769 !> \param pi_block ...
770 !> \param pj_block ...
771 !> \param e1b ...
772 !> \param e2a ...
773 !> \param ecore ...
774 !> \param rp ...
775 !> \date 12.2008 [tlaino]
776 !> \author Teodoro Laino [tlaino]
777 ! **************************************************************************************************
778  SUBROUTINE fock2_1el_r3(sepi, sepj, ksi_block, ksj_block, pi_block, pj_block, &
779  e1b, e2a, ecore, rp)
780  TYPE(semi_empirical_type), POINTER :: sepi, sepj
781  REAL(kind=dp), DIMENSION(:, :), POINTER :: ksi_block, ksj_block
782  REAL(kind=dp), &
783  DIMENSION(sepi%natorb, sepi%natorb), INTENT(IN) :: pi_block
784  REAL(kind=dp), &
785  DIMENSION(sepj%natorb, sepj%natorb), INTENT(IN) :: pj_block
786  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: e1b, e2a
787  REAL(kind=dp), DIMENSION(2), INTENT(INOUT) :: ecore
788  REAL(kind=dp), INTENT(IN) :: rp
789 
790  INTEGER :: i1, i1l, i2
791 
792 !
793 ! Add the electron-nuclear residual attraction term for atom sepi
794 !
795 
796  i2 = 0
797  DO i1l = 1, sepi%natorb
798  i2 = i2 + 1
799  i1 = se_orbital_pointer(i1l)
800  ksi_block(i1, i1) = ksi_block(i1, i1) + e1b(i2)*rp
801  ecore(1) = ecore(1) + e1b(i2)*rp*pi_block(i1, i1)
802  END DO
803  !
804  ! Add the electron-nuclear residual attraction term for atom sepj
805  !
806  i2 = 0
807  DO i1l = 1, sepj%natorb
808  i2 = i2 + 1
809  i1 = se_orbital_pointer(i1l)
810  ksj_block(i1, i1) = ksj_block(i1, i1) + e2a(i2)*rp
811  ecore(2) = ecore(2) + e2a(i2)*rp*pj_block(i1, i1)
812  END DO
813 
814  END SUBROUTINE fock2_1el_r3
815 
816 ! **************************************************************************************************
817 !> \brief Derivatives of 2-center 1-electron Fock Matrix residual (1/R^3)
818 !> integral part
819 !> \param sepi ...
820 !> \param sepj ...
821 !> \param drp ...
822 !> \param pi_block ...
823 !> \param pj_block ...
824 !> \param force ...
825 !> \param e1b ...
826 !> \param e2a ...
827 !> \date 12.2008 [tlaino]
828 !> \author Teodoro Laino [tlaino]
829 ! **************************************************************************************************
830  SUBROUTINE dfock2_1el_r3(sepi, sepj, drp, pi_block, pj_block, force, e1b, e2a)
831  TYPE(semi_empirical_type), POINTER :: sepi, sepj
832  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: drp
833  REAL(kind=dp), &
834  DIMENSION(sepi%natorb, sepi%natorb), INTENT(IN) :: pi_block
835  REAL(kind=dp), &
836  DIMENSION(sepj%natorb, sepj%natorb), INTENT(IN) :: pj_block
837  REAL(kind=dp), DIMENSION(3), INTENT(INOUT) :: force
838  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: e1b, e2a
839 
840  INTEGER :: i1, i1l, i2
841  REAL(kind=dp) :: tmp
842 
843 !
844 ! Add the electron-nuclear residual attraction term for atom sepi
845 !
846 
847  i2 = 0
848  DO i1l = 1, sepi%natorb
849  i1 = se_orbital_pointer(i1l)
850  i2 = i2 + 1
851  tmp = e1b(i2)*pi_block(i1, i1)
852  force(1) = force(1) + tmp*drp(1)
853  force(2) = force(2) + tmp*drp(2)
854  force(3) = force(3) + tmp*drp(3)
855  END DO
856  !
857  ! Add the electron-nuclear attraction term for atom sepj
858  !
859  i2 = 0
860  DO i1l = 1, sepj%natorb
861  i1 = se_orbital_pointer(i1l)
862  i2 = i2 + 1
863  tmp = e2a(i2)*pj_block(i1, i1)
864  force(1) = force(1) + tmp*drp(1)
865  force(2) = force(2) + tmp*drp(2)
866  force(3) = force(3) + tmp*drp(3)
867  END DO
868 
869  END SUBROUTINE dfock2_1el_r3
870 
871 ! **************************************************************************************************
872 !> \brief Construction of 2-center Fock Matrix - Coulomb Terms for the residual
873 !> (1/R^3) integral part
874 !> \param sepi ...
875 !> \param sepj ...
876 !> \param switch ...
877 !> \param pi_tot ...
878 !> \param fi_mat DIMENSION(sepi%natorb, sepi%natorb)
879 !> \param pj_tot ...
880 !> \param fj_mat DIMENSION(sepj%natorb, sepj%natorb)
881 !> \param factor ...
882 !> \param w ...
883 !> \param rp ...
884 !> \date 12.2008 [tlaino]
885 !> \author Teodoro Laino [tlaino]
886 ! **************************************************************************************************
887  SUBROUTINE fock2c_r3(sepi, sepj, switch, pi_tot, fi_mat, pj_tot, fj_mat, &
888  factor, w, rp)
889  TYPE(semi_empirical_type), POINTER :: sepi, sepj
890  LOGICAL, INTENT(IN) :: switch
891  REAL(kind=dp), DIMENSION(45, 45), INTENT(IN) :: pi_tot
892  REAL(kind=dp), DIMENSION(:, :), POINTER :: fi_mat
893  REAL(kind=dp), DIMENSION(45, 45), INTENT(IN) :: pj_tot
894  REAL(kind=dp), DIMENSION(:, :), POINTER :: fj_mat
895  REAL(kind=dp), INTENT(IN) :: factor
896  REAL(kind=dp), DIMENSION(81), INTENT(IN) :: w
897  REAL(kind=dp), INTENT(IN) :: rp
898 
899  INTEGER :: i, il, ind, j, k, kl, kr, natorb(2)
900  REAL(kind=dp) :: a, w_l(81)
901 
902  natorb(1) = sepi%natorb
903  natorb(2) = sepj%natorb
904  IF (switch) THEN
905  natorb(1) = sepj%natorb
906  natorb(2) = sepi%natorb
907  ! Reshuffle the integral array (natural storage order is sepi/sepj)
908  kr = 0
909  DO i = 1, sepj%natorb
910  DO j = 1, sepi%natorb
911  kr = kr + 1
912  ind = (j - 1)*sepj%natorb + i
913  w_l(kr) = w(ind)
914  END DO
915  END DO
916  ELSE
917  w_l = w
918  END IF
919 
920  ! Modify the Fock Matrix
921  kr = 0
922  DO il = 1, natorb(1)
923  i = se_orbital_pointer(il)
924  DO kl = 1, natorb(2)
925  k = se_orbital_pointer(kl)
926  kr = kr + 1
927  a = w_l(kr)*factor*rp
928  ! Coulomb
929  IF (.NOT. switch) THEN
930  fi_mat(i, i) = fi_mat(i, i) + a*pj_tot(k, k)
931  fj_mat(k, k) = fj_mat(k, k) + a*pi_tot(i, i)
932  ELSE
933  fj_mat(i, i) = fj_mat(i, i) + a*pi_tot(k, k)
934  fi_mat(k, k) = fi_mat(k, k) + a*pj_tot(i, i)
935  END IF
936  END DO
937  END DO
938 
939  END SUBROUTINE fock2c_r3
940 
941 ! **************************************************************************************************
942 !> \brief Derivatives of 2-center Fock Matrix - Coulomb Terms for the residual
943 !> (1/R^3) integral part
944 !> \param sepi ...
945 !> \param sepj ...
946 !> \param switch ...
947 !> \param pi_tot ...
948 !> \param pj_tot ...
949 !> \param factor ...
950 !> \param w ...
951 !> \param drp ...
952 !> \param force ...
953 !> \date 12.2008 [tlaino]
954 !> \author Teodoro Laino [tlaino]
955 ! **************************************************************************************************
956  SUBROUTINE dfock2c_r3(sepi, sepj, switch, pi_tot, pj_tot, factor, w, drp, &
957  force)
958  TYPE(semi_empirical_type), POINTER :: sepi, sepj
959  LOGICAL, INTENT(IN) :: switch
960  REAL(kind=dp), DIMENSION(45, 45), INTENT(IN) :: pi_tot, pj_tot
961  REAL(kind=dp), INTENT(IN) :: factor
962  REAL(kind=dp), DIMENSION(81), INTENT(IN) :: w
963  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: drp
964  REAL(kind=dp), DIMENSION(3), INTENT(INOUT) :: force
965 
966  INTEGER :: i, il, ind, j, k, kl, kr, natorb(2)
967  REAL(kind=dp) :: a(3), tmp, w_l(81)
968 
969  natorb(1) = sepi%natorb
970  natorb(2) = sepj%natorb
971  IF (switch) THEN
972  natorb(1) = sepj%natorb
973  natorb(2) = sepi%natorb
974  ! Reshuffle the integral array (natural storage order is sepi/sepj)
975  kr = 0
976  DO i = 1, sepj%natorb
977  DO j = 1, sepi%natorb
978  kr = kr + 1
979  ind = (j - 1)*sepj%natorb + i
980  w_l(kr) = w(ind)
981  END DO
982  END DO
983  ELSE
984  w_l = w
985  END IF
986 
987  ! Modify the Fock Matrix
988  kr = 0
989  DO il = 1, natorb(1)
990  i = se_orbital_pointer(il)
991  DO kl = 1, natorb(2)
992  k = se_orbital_pointer(kl)
993  kr = kr + 1
994  tmp = w_l(kr)*factor
995  a(1) = tmp*drp(1)
996  a(2) = tmp*drp(2)
997  a(3) = tmp*drp(3)
998  ! Coulomb
999  IF (.NOT. switch) THEN
1000  tmp = pj_tot(k, k)*pi_tot(i, i)
1001  ELSE
1002  tmp = pi_tot(k, k)*pj_tot(i, i)
1003  END IF
1004  force(1) = force(1) + a(1)*tmp
1005  force(2) = force(2) + a(2)*tmp
1006  force(3) = force(3) + a(3)*tmp
1007  END DO
1008  END DO
1009 
1010  END SUBROUTINE dfock2c_r3
1011 
1012 ! **************************************************************************************************
1013 !> \brief Coulomb interaction multipolar correction
1014 !> \param atom_a ...
1015 !> \param atom_b ...
1016 !> \param my_task ...
1017 !> \param do_forces ...
1018 !> \param do_efield ...
1019 !> \param do_stress ...
1020 !> \param charges ...
1021 !> \param dipoles ...
1022 !> \param quadrupoles ...
1023 !> \param force_ab ...
1024 !> \param efield0 ...
1025 !> \param efield1 ...
1026 !> \param efield2 ...
1027 !> \param rab2 ...
1028 !> \param rab ...
1029 !> \param integral_value ...
1030 !> \param ptens11 ...
1031 !> \param ptens12 ...
1032 !> \param ptens13 ...
1033 !> \param ptens21 ...
1034 !> \param ptens22 ...
1035 !> \param ptens23 ...
1036 !> \param ptens31 ...
1037 !> \param ptens32 ...
1038 !> \param ptens33 ...
1039 !> \date 05.2009 [tlaino]
1040 !> \author Teodoro Laino [tlaino]
1041 ! **************************************************************************************************
1042  SUBROUTINE se_coulomb_ij_interaction(atom_a, atom_b, my_task, do_forces, do_efield, &
1043  do_stress, charges, dipoles, quadrupoles, force_ab, efield0, efield1, efield2, &
1044  rab2, rab, integral_value, ptens11, ptens12, ptens13, ptens21, ptens22, ptens23, &
1045  ptens31, ptens32, ptens33)
1046  INTEGER, INTENT(IN) :: atom_a, atom_b
1047  LOGICAL, DIMENSION(3) :: my_task
1048  LOGICAL, INTENT(IN) :: do_forces, do_efield, do_stress
1049  REAL(kind=dp), DIMENSION(:), POINTER :: charges
1050  REAL(kind=dp), DIMENSION(:, :), POINTER :: dipoles
1051  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: quadrupoles
1052  REAL(kind=dp), DIMENSION(3), INTENT(OUT) :: force_ab
1053  REAL(kind=dp), DIMENSION(:), POINTER :: efield0
1054  REAL(kind=dp), DIMENSION(:, :), POINTER :: efield1, efield2
1055  REAL(kind=dp), INTENT(IN) :: rab2
1056  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
1057  REAL(kind=dp), INTENT(OUT), OPTIONAL :: integral_value
1058  REAL(kind=dp), INTENT(INOUT) :: ptens11, ptens12, ptens13, ptens21, &
1059  ptens22, ptens23, ptens31, ptens32, &
1060  ptens33
1061 
1062  INTEGER :: a, b, c, d, e, i, j, k
1063  LOGICAL :: do_efield0, do_efield1, do_efield2, &
1064  force_eval
1065  LOGICAL, DIMENSION(3) :: do_task
1066  LOGICAL, DIMENSION(3, 3) :: task
1067  REAL(kind=dp) :: ch_i, ch_j, ef0_i, ef0_j, eloc, energy, fac, fac_ij, ir, irab2, r, tij, &
1068  tmp, tmp1, tmp11, tmp12, tmp13, tmp2, tmp21, tmp22, tmp23, tmp31, tmp32, tmp33, tmp_ij, &
1069  tmp_ji
1070  REAL(kind=dp), DIMENSION(0:5) :: f
1071  REAL(kind=dp), DIMENSION(3) :: dp_i, dp_j, ef1_i, ef1_j, fr, tij_a
1072  REAL(kind=dp), DIMENSION(3, 3) :: ef2_i, ef2_j, qp_i, qp_j, tij_ab
1073  REAL(kind=dp), DIMENSION(3, 3, 3) :: tij_abc
1074  REAL(kind=dp), DIMENSION(3, 3, 3, 3) :: tij_abcd
1075  REAL(kind=dp), DIMENSION(3, 3, 3, 3, 3) :: tij_abcde
1076 
1077  do_task = my_task
1078  energy = 0.0_dp
1079  DO i = 1, 3
1080  IF (do_task(i)) THEN
1081  SELECT CASE (i)
1082  CASE (1)
1083  do_task(1) = (charges(atom_a) /= 0.0_dp) .OR. (charges(atom_b) /= 0.0_dp)
1084  CASE (2)
1085  do_task(2) = (any(dipoles(:, atom_a) /= 0.0_dp)) .OR. (any(dipoles(:, atom_b) /= 0.0_dp))
1086  CASE (3)
1087  do_task(3) = (any(quadrupoles(:, :, atom_a) /= 0.0_dp)) .OR. (any(quadrupoles(:, :, atom_b) /= 0.0_dp))
1088  END SELECT
1089  END IF
1090  END DO
1091  DO i = 1, 3
1092  DO j = i, 3
1093  task(j, i) = do_task(i) .AND. do_task(j)
1094  task(i, j) = task(j, i)
1095  END DO
1096  END DO
1097  do_efield0 = do_efield .AND. ASSOCIATED(efield0)
1098  do_efield1 = do_efield .AND. ASSOCIATED(efield1)
1099  do_efield2 = do_efield .AND. ASSOCIATED(efield2)
1100 
1101  fac_ij = 1.0_dp
1102  IF (atom_a == atom_b) fac_ij = 0.5_dp
1103 
1104  ! Compute the Short Range constribution according the task
1105  IF (debug_this_module) THEN
1106  f = huge(0.0_dp)
1107  tij = huge(0.0_dp)
1108  tij_a = huge(0.0_dp)
1109  tij_ab = huge(0.0_dp)
1110  tij_abc = huge(0.0_dp)
1111  tij_abcd = huge(0.0_dp)
1112  tij_abcde = huge(0.0_dp)
1113  END IF
1114  r = sqrt(rab2)
1115  irab2 = 1.0_dp/rab2
1116  ir = 1.0_dp/r
1117 
1118  ! Compute the radial function
1119  ! code for point multipole without screening
1120  f(0) = ir
1121  DO i = 1, 5
1122  f(i) = irab2*f(i - 1)
1123  END DO
1124 
1125 
1126  ! Compute the Tensor components
1127  force_eval = do_stress
1128  IF (task(1, 1)) THEN
1129  tij = f(0)*fac_ij
1130  force_eval = do_forces .OR. do_efield1
1131  END IF
1132  IF (task(2, 2)) force_eval = force_eval .OR. do_efield0
1133  IF (task(1, 2) .OR. force_eval) THEN
1134  force_eval = do_stress
1135  tij_a = -rab*f(1)*fac_ij
1136  IF (task(1, 2)) force_eval = force_eval .OR. do_forces
1137  END IF
1138  IF (task(1, 1)) force_eval = force_eval .OR. do_efield2
1139  IF (task(3, 3)) force_eval = force_eval .OR. do_efield0
1140  IF (task(2, 2) .OR. task(3, 1) .OR. force_eval) THEN
1141  force_eval = do_stress
1142  DO b = 1, 3
1143  DO a = 1, 3
1144  tmp = rab(a)*rab(b)*fac_ij
1145  tij_ab(a, b) = 3.0_dp*tmp*f(2)
1146  IF (a == b) tij_ab(a, b) = tij_ab(a, b) - f(1)*fac_ij
1147  END DO
1148  END DO
1149  IF (task(2, 2) .OR. task(3, 1)) force_eval = force_eval .OR. do_forces
1150  END IF
1151  IF (task(2, 2)) force_eval = force_eval .OR. do_efield2
1152  IF (task(3, 3)) force_eval = force_eval .OR. do_efield1
1153  IF (task(3, 2) .OR. force_eval) THEN
1154  force_eval = do_stress
1155  DO c = 1, 3
1156  DO b = 1, 3
1157  DO a = 1, 3
1158  tmp = rab(a)*rab(b)*rab(c)*fac_ij
1159  tij_abc(a, b, c) = -15.0_dp*tmp*f(3)
1160  tmp = 3.0_dp*f(2)*fac_ij
1161  IF (a == b) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(c)
1162  IF (a == c) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(b)
1163  IF (b == c) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(a)
1164  END DO
1165  END DO
1166  END DO
1167  IF (task(3, 2)) force_eval = force_eval .OR. do_forces
1168  END IF
1169  IF (task(3, 3) .OR. force_eval) THEN
1170  force_eval = do_stress
1171  DO d = 1, 3
1172  DO c = 1, 3
1173  DO b = 1, 3
1174  DO a = 1, 3
1175  tmp = rab(a)*rab(b)*rab(c)*rab(d)*fac_ij
1176  tij_abcd(a, b, c, d) = 105.0_dp*tmp*f(4)
1177  tmp1 = 15.0_dp*f(3)*fac_ij
1178  tmp2 = 3.0_dp*f(2)*fac_ij
1179  IF (a == b) THEN
1180  tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(c)*rab(d)
1181  IF (c == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
1182  END IF
1183  IF (a == c) THEN
1184  tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(b)*rab(d)
1185  IF (b == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
1186  END IF
1187  IF (a == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(b)*rab(c)
1188  IF (b == c) THEN
1189  tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(d)
1190  IF (a == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
1191  END IF
1192  IF (b == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(c)
1193  IF (c == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(b)
1194  END DO
1195  END DO
1196  END DO
1197  END DO
1198  IF (task(3, 3)) force_eval = force_eval .OR. do_forces
1199  END IF
1200  IF (force_eval) THEN
1201  force_eval = do_stress
1202  DO e = 1, 3
1203  DO d = 1, 3
1204  DO c = 1, 3
1205  DO b = 1, 3
1206  DO a = 1, 3
1207  tmp = rab(a)*rab(b)*rab(c)*rab(d)*rab(e)*fac_ij
1208  tij_abcde(a, b, c, d, e) = -945.0_dp*tmp*f(5)
1209  tmp1 = 105.0_dp*f(4)*fac_ij
1210  tmp2 = 15.0_dp*f(3)*fac_ij
1211  IF (a == b) THEN
1212  tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(c)*rab(d)*rab(e)
1213  IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
1214  IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
1215  IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
1216  END IF
1217  IF (a == c) THEN
1218  tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(d)*rab(e)
1219  IF (b == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
1220  IF (b == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
1221  IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
1222  END IF
1223  IF (a == d) THEN
1224  tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(c)*rab(e)
1225  IF (b == c) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
1226  IF (b == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
1227  IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
1228  END IF
1229  IF (a == e) THEN
1230  tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(c)*rab(d)
1231  IF (b == c) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
1232  IF (b == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
1233  IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
1234  END IF
1235  IF (b == c) THEN
1236  tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(d)*rab(e)
1237  IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
1238  END IF
1239  IF (b == d) THEN
1240  tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(c)*rab(e)
1241  IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
1242  END IF
1243  IF (b == e) THEN
1244  tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(c)*rab(d)
1245  IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
1246  END IF
1247  IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(e)
1248  IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(d)
1249  IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(c)
1250  END DO
1251  END DO
1252  END DO
1253  END DO
1254  END DO
1255  END IF
1256  eloc = 0.0_dp
1257  fr = 0.0_dp
1258  ef0_i = 0.0_dp
1259  ef0_j = 0.0_dp
1260  ef1_j = 0.0_dp
1261  ef1_i = 0.0_dp
1262  ef2_j = 0.0_dp
1263  ef2_i = 0.0_dp
1264 
1265 
1266  ! Initialize the charge, dipole and quadrupole for atom A and B
1267  IF (debug_this_module) THEN
1268  ch_j = huge(0.0_dp)
1269  ch_i = huge(0.0_dp)
1270  dp_j = huge(0.0_dp)
1271  dp_i = huge(0.0_dp)
1272  qp_j = huge(0.0_dp)
1273  qp_i = huge(0.0_dp)
1274  END IF
1275  IF (any(task(1, :))) THEN
1276  ch_j = charges(atom_a)
1277  ch_i = charges(atom_b)
1278  END IF
1279  IF (any(task(2, :))) THEN
1280  dp_j = dipoles(:, atom_a)
1281  dp_i = dipoles(:, atom_b)
1282  END IF
1283  IF (any(task(3, :))) THEN
1284  qp_j = quadrupoles(:, :, atom_a)
1285  qp_i = quadrupoles(:, :, atom_b)
1286  END IF
1287  IF (task(1, 1)) THEN
1288  ! Charge - Charge
1289  eloc = eloc + ch_i*tij*ch_j
1290  ! Forces on particle i (locally b)
1291  IF (do_forces .OR. do_stress) THEN
1292  fr(1) = fr(1) - ch_j*tij_a(1)*ch_i
1293  fr(2) = fr(2) - ch_j*tij_a(2)*ch_i
1294  fr(3) = fr(3) - ch_j*tij_a(3)*ch_i
1295  END IF
1296  ! Electric fields
1297  IF (do_efield) THEN
1298  ! Potential
1299  IF (do_efield0) THEN
1300  ef0_i = ef0_i + tij*ch_j
1301 
1302  ef0_j = ef0_j + tij*ch_i
1303  END IF
1304  ! Electric field
1305  IF (do_efield1) THEN
1306  ef1_i(1) = ef1_i(1) - tij_a(1)*ch_j
1307  ef1_i(2) = ef1_i(2) - tij_a(2)*ch_j
1308  ef1_i(3) = ef1_i(3) - tij_a(3)*ch_j
1309 
1310  ef1_j(1) = ef1_j(1) + tij_a(1)*ch_i
1311  ef1_j(2) = ef1_j(2) + tij_a(2)*ch_i
1312  ef1_j(3) = ef1_j(3) + tij_a(3)*ch_i
1313 
1314 
1315  END IF
1316  ! Electric field gradient
1317  IF (do_efield2) THEN
1318  ef2_i(1, 1) = ef2_i(1, 1) - tij_ab(1, 1)*ch_j
1319  ef2_i(2, 1) = ef2_i(2, 1) - tij_ab(2, 1)*ch_j
1320  ef2_i(3, 1) = ef2_i(3, 1) - tij_ab(3, 1)*ch_j
1321  ef2_i(1, 2) = ef2_i(1, 2) - tij_ab(1, 2)*ch_j
1322  ef2_i(2, 2) = ef2_i(2, 2) - tij_ab(2, 2)*ch_j
1323  ef2_i(3, 2) = ef2_i(3, 2) - tij_ab(3, 2)*ch_j
1324  ef2_i(1, 3) = ef2_i(1, 3) - tij_ab(1, 3)*ch_j
1325  ef2_i(2, 3) = ef2_i(2, 3) - tij_ab(2, 3)*ch_j
1326  ef2_i(3, 3) = ef2_i(3, 3) - tij_ab(3, 3)*ch_j
1327 
1328  ef2_j(1, 1) = ef2_j(1, 1) - tij_ab(1, 1)*ch_i
1329  ef2_j(2, 1) = ef2_j(2, 1) - tij_ab(2, 1)*ch_i
1330  ef2_j(3, 1) = ef2_j(3, 1) - tij_ab(3, 1)*ch_i
1331  ef2_j(1, 2) = ef2_j(1, 2) - tij_ab(1, 2)*ch_i
1332  ef2_j(2, 2) = ef2_j(2, 2) - tij_ab(2, 2)*ch_i
1333  ef2_j(3, 2) = ef2_j(3, 2) - tij_ab(3, 2)*ch_i
1334  ef2_j(1, 3) = ef2_j(1, 3) - tij_ab(1, 3)*ch_i
1335  ef2_j(2, 3) = ef2_j(2, 3) - tij_ab(2, 3)*ch_i
1336  ef2_j(3, 3) = ef2_j(3, 3) - tij_ab(3, 3)*ch_i
1337  END IF
1338  END IF
1339  END IF
1340  IF (task(2, 2)) THEN
1341  ! Dipole - Dipole
1342  tmp = -(dp_i(1)*(tij_ab(1, 1)*dp_j(1) + &
1343  tij_ab(2, 1)*dp_j(2) + &
1344  tij_ab(3, 1)*dp_j(3)) + &
1345  dp_i(2)*(tij_ab(1, 2)*dp_j(1) + &
1346  tij_ab(2, 2)*dp_j(2) + &
1347  tij_ab(3, 2)*dp_j(3)) + &
1348  dp_i(3)*(tij_ab(1, 3)*dp_j(1) + &
1349  tij_ab(2, 3)*dp_j(2) + &
1350  tij_ab(3, 3)*dp_j(3)))
1351  eloc = eloc + tmp
1352  ! Forces on particle i (locally b)
1353  IF (do_forces .OR. do_stress) THEN
1354  DO k = 1, 3
1355  fr(k) = fr(k) + dp_i(1)*(tij_abc(1, 1, k)*dp_j(1) + &
1356  tij_abc(2, 1, k)*dp_j(2) + &
1357  tij_abc(3, 1, k)*dp_j(3)) &
1358  + dp_i(2)*(tij_abc(1, 2, k)*dp_j(1) + &
1359  tij_abc(2, 2, k)*dp_j(2) + &
1360  tij_abc(3, 2, k)*dp_j(3)) &
1361  + dp_i(3)*(tij_abc(1, 3, k)*dp_j(1) + &
1362  tij_abc(2, 3, k)*dp_j(2) + &
1363  tij_abc(3, 3, k)*dp_j(3))
1364  END DO
1365  END IF
1366  ! Electric fields
1367  IF (do_efield) THEN
1368  ! Potential
1369  IF (do_efield0) THEN
1370  ef0_i = ef0_i - (tij_a(1)*dp_j(1) + &
1371  tij_a(2)*dp_j(2) + &
1372  tij_a(3)*dp_j(3))
1373 
1374  ef0_j = ef0_j + (tij_a(1)*dp_i(1) + &
1375  tij_a(2)*dp_i(2) + &
1376  tij_a(3)*dp_i(3))
1377  END IF
1378  ! Electric field
1379  IF (do_efield1) THEN
1380  ef1_i(1) = ef1_i(1) + (tij_ab(1, 1)*dp_j(1) + &
1381  tij_ab(2, 1)*dp_j(2) + &
1382  tij_ab(3, 1)*dp_j(3))
1383  ef1_i(2) = ef1_i(2) + (tij_ab(1, 2)*dp_j(1) + &
1384  tij_ab(2, 2)*dp_j(2) + &
1385  tij_ab(3, 2)*dp_j(3))
1386  ef1_i(3) = ef1_i(3) + (tij_ab(1, 3)*dp_j(1) + &
1387  tij_ab(2, 3)*dp_j(2) + &
1388  tij_ab(3, 3)*dp_j(3))
1389 
1390  ef1_j(1) = ef1_j(1) + (tij_ab(1, 1)*dp_i(1) + &
1391  tij_ab(2, 1)*dp_i(2) + &
1392  tij_ab(3, 1)*dp_i(3))
1393  ef1_j(2) = ef1_j(2) + (tij_ab(1, 2)*dp_i(1) + &
1394  tij_ab(2, 2)*dp_i(2) + &
1395  tij_ab(3, 2)*dp_i(3))
1396  ef1_j(3) = ef1_j(3) + (tij_ab(1, 3)*dp_i(1) + &
1397  tij_ab(2, 3)*dp_i(2) + &
1398  tij_ab(3, 3)*dp_i(3))
1399  END IF
1400  ! Electric field gradient
1401  IF (do_efield2) THEN
1402  ef2_i(1, 1) = ef2_i(1, 1) + (tij_abc(1, 1, 1)*dp_j(1) + &
1403  tij_abc(2, 1, 1)*dp_j(2) + &
1404  tij_abc(3, 1, 1)*dp_j(3))
1405  ef2_i(1, 2) = ef2_i(1, 2) + (tij_abc(1, 1, 2)*dp_j(1) + &
1406  tij_abc(2, 1, 2)*dp_j(2) + &
1407  tij_abc(3, 1, 2)*dp_j(3))
1408  ef2_i(1, 3) = ef2_i(1, 3) + (tij_abc(1, 1, 3)*dp_j(1) + &
1409  tij_abc(2, 1, 3)*dp_j(2) + &
1410  tij_abc(3, 1, 3)*dp_j(3))
1411  ef2_i(2, 1) = ef2_i(2, 1) + (tij_abc(1, 2, 1)*dp_j(1) + &
1412  tij_abc(2, 2, 1)*dp_j(2) + &
1413  tij_abc(3, 2, 1)*dp_j(3))
1414  ef2_i(2, 2) = ef2_i(2, 2) + (tij_abc(1, 2, 2)*dp_j(1) + &
1415  tij_abc(2, 2, 2)*dp_j(2) + &
1416  tij_abc(3, 2, 2)*dp_j(3))
1417  ef2_i(2, 3) = ef2_i(2, 3) + (tij_abc(1, 2, 3)*dp_j(1) + &
1418  tij_abc(2, 2, 3)*dp_j(2) + &
1419  tij_abc(3, 2, 3)*dp_j(3))
1420  ef2_i(3, 1) = ef2_i(3, 1) + (tij_abc(1, 3, 1)*dp_j(1) + &
1421  tij_abc(2, 3, 1)*dp_j(2) + &
1422  tij_abc(3, 3, 1)*dp_j(3))
1423  ef2_i(3, 2) = ef2_i(3, 2) + (tij_abc(1, 3, 2)*dp_j(1) + &
1424  tij_abc(2, 3, 2)*dp_j(2) + &
1425  tij_abc(3, 3, 2)*dp_j(3))
1426  ef2_i(3, 3) = ef2_i(3, 3) + (tij_abc(1, 3, 3)*dp_j(1) + &
1427  tij_abc(2, 3, 3)*dp_j(2) + &
1428  tij_abc(3, 3, 3)*dp_j(3))
1429 
1430  ef2_j(1, 1) = ef2_j(1, 1) - (tij_abc(1, 1, 1)*dp_i(1) + &
1431  tij_abc(2, 1, 1)*dp_i(2) + &
1432  tij_abc(3, 1, 1)*dp_i(3))
1433  ef2_j(1, 2) = ef2_j(1, 2) - (tij_abc(1, 1, 2)*dp_i(1) + &
1434  tij_abc(2, 1, 2)*dp_i(2) + &
1435  tij_abc(3, 1, 2)*dp_i(3))
1436  ef2_j(1, 3) = ef2_j(1, 3) - (tij_abc(1, 1, 3)*dp_i(1) + &
1437  tij_abc(2, 1, 3)*dp_i(2) + &
1438  tij_abc(3, 1, 3)*dp_i(3))
1439  ef2_j(2, 1) = ef2_j(2, 1) - (tij_abc(1, 2, 1)*dp_i(1) + &
1440  tij_abc(2, 2, 1)*dp_i(2) + &
1441  tij_abc(3, 2, 1)*dp_i(3))
1442  ef2_j(2, 2) = ef2_j(2, 2) - (tij_abc(1, 2, 2)*dp_i(1) + &
1443  tij_abc(2, 2, 2)*dp_i(2) + &
1444  tij_abc(3, 2, 2)*dp_i(3))
1445  ef2_j(2, 3) = ef2_j(2, 3) - (tij_abc(1, 2, 3)*dp_i(1) + &
1446  tij_abc(2, 2, 3)*dp_i(2) + &
1447  tij_abc(3, 2, 3)*dp_i(3))
1448  ef2_j(3, 1) = ef2_j(3, 1) - (tij_abc(1, 3, 1)*dp_i(1) + &
1449  tij_abc(2, 3, 1)*dp_i(2) + &
1450  tij_abc(3, 3, 1)*dp_i(3))
1451  ef2_j(3, 2) = ef2_j(3, 2) - (tij_abc(1, 3, 2)*dp_i(1) + &
1452  tij_abc(2, 3, 2)*dp_i(2) + &
1453  tij_abc(3, 3, 2)*dp_i(3))
1454  ef2_j(3, 3) = ef2_j(3, 3) - (tij_abc(1, 3, 3)*dp_i(1) + &
1455  tij_abc(2, 3, 3)*dp_i(2) + &
1456  tij_abc(3, 3, 3)*dp_i(3))
1457  END IF
1458  END IF
1459  END IF
1460  IF (task(2, 1)) THEN
1461  ! Dipole - Charge
1462  tmp = ch_j*(tij_a(1)*dp_i(1) + &
1463  tij_a(2)*dp_i(2) + &
1464  tij_a(3)*dp_i(3)) &
1465  - ch_i*(tij_a(1)*dp_j(1) + &
1466  tij_a(2)*dp_j(2) + &
1467  tij_a(3)*dp_j(3))
1468  eloc = eloc + tmp
1469  ! Forces on particle i (locally b)
1470  IF (do_forces .OR. do_stress) THEN
1471  DO k = 1, 3
1472  fr(k) = fr(k) - ch_j*(tij_ab(1, k)*dp_i(1) + &
1473  tij_ab(2, k)*dp_i(2) + &
1474  tij_ab(3, k)*dp_i(3)) &
1475  + ch_i*(tij_ab(1, k)*dp_j(1) + &
1476  tij_ab(2, k)*dp_j(2) + &
1477  tij_ab(3, k)*dp_j(3))
1478  END DO
1479  END IF
1480  END IF
1481  IF (task(3, 3)) THEN
1482  ! Quadrupole - Quadrupole
1483  fac = 1.0_dp/9.0_dp
1484  tmp11 = qp_i(1, 1)*(tij_abcd(1, 1, 1, 1)*qp_j(1, 1) + &
1485  tij_abcd(2, 1, 1, 1)*qp_j(2, 1) + &
1486  tij_abcd(3, 1, 1, 1)*qp_j(3, 1) + &
1487  tij_abcd(1, 2, 1, 1)*qp_j(1, 2) + &
1488  tij_abcd(2, 2, 1, 1)*qp_j(2, 2) + &
1489  tij_abcd(3, 2, 1, 1)*qp_j(3, 2) + &
1490  tij_abcd(1, 3, 1, 1)*qp_j(1, 3) + &
1491  tij_abcd(2, 3, 1, 1)*qp_j(2, 3) + &
1492  tij_abcd(3, 3, 1, 1)*qp_j(3, 3))
1493  tmp21 = qp_i(2, 1)*(tij_abcd(1, 1, 1, 2)*qp_j(1, 1) + &
1494  tij_abcd(2, 1, 1, 2)*qp_j(2, 1) + &
1495  tij_abcd(3, 1, 1, 2)*qp_j(3, 1) + &
1496  tij_abcd(1, 2, 1, 2)*qp_j(1, 2) + &
1497  tij_abcd(2, 2, 1, 2)*qp_j(2, 2) + &
1498  tij_abcd(3, 2, 1, 2)*qp_j(3, 2) + &
1499  tij_abcd(1, 3, 1, 2)*qp_j(1, 3) + &
1500  tij_abcd(2, 3, 1, 2)*qp_j(2, 3) + &
1501  tij_abcd(3, 3, 1, 2)*qp_j(3, 3))
1502  tmp31 = qp_i(3, 1)*(tij_abcd(1, 1, 1, 3)*qp_j(1, 1) + &
1503  tij_abcd(2, 1, 1, 3)*qp_j(2, 1) + &
1504  tij_abcd(3, 1, 1, 3)*qp_j(3, 1) + &
1505  tij_abcd(1, 2, 1, 3)*qp_j(1, 2) + &
1506  tij_abcd(2, 2, 1, 3)*qp_j(2, 2) + &
1507  tij_abcd(3, 2, 1, 3)*qp_j(3, 2) + &
1508  tij_abcd(1, 3, 1, 3)*qp_j(1, 3) + &
1509  tij_abcd(2, 3, 1, 3)*qp_j(2, 3) + &
1510  tij_abcd(3, 3, 1, 3)*qp_j(3, 3))
1511  tmp22 = qp_i(2, 2)*(tij_abcd(1, 1, 2, 2)*qp_j(1, 1) + &
1512  tij_abcd(2, 1, 2, 2)*qp_j(2, 1) + &
1513  tij_abcd(3, 1, 2, 2)*qp_j(3, 1) + &
1514  tij_abcd(1, 2, 2, 2)*qp_j(1, 2) + &
1515  tij_abcd(2, 2, 2, 2)*qp_j(2, 2) + &
1516  tij_abcd(3, 2, 2, 2)*qp_j(3, 2) + &
1517  tij_abcd(1, 3, 2, 2)*qp_j(1, 3) + &
1518  tij_abcd(2, 3, 2, 2)*qp_j(2, 3) + &
1519  tij_abcd(3, 3, 2, 2)*qp_j(3, 3))
1520  tmp32 = qp_i(3, 2)*(tij_abcd(1, 1, 2, 3)*qp_j(1, 1) + &
1521  tij_abcd(2, 1, 2, 3)*qp_j(2, 1) + &
1522  tij_abcd(3, 1, 2, 3)*qp_j(3, 1) + &
1523  tij_abcd(1, 2, 2, 3)*qp_j(1, 2) + &
1524  tij_abcd(2, 2, 2, 3)*qp_j(2, 2) + &
1525  tij_abcd(3, 2, 2, 3)*qp_j(3, 2) + &
1526  tij_abcd(1, 3, 2, 3)*qp_j(1, 3) + &
1527  tij_abcd(2, 3, 2, 3)*qp_j(2, 3) + &
1528  tij_abcd(3, 3, 2, 3)*qp_j(3, 3))
1529  tmp33 = qp_i(3, 3)*(tij_abcd(1, 1, 3, 3)*qp_j(1, 1) + &
1530  tij_abcd(2, 1, 3, 3)*qp_j(2, 1) + &
1531  tij_abcd(3, 1, 3, 3)*qp_j(3, 1) + &
1532  tij_abcd(1, 2, 3, 3)*qp_j(1, 2) + &
1533  tij_abcd(2, 2, 3, 3)*qp_j(2, 2) + &
1534  tij_abcd(3, 2, 3, 3)*qp_j(3, 2) + &
1535  tij_abcd(1, 3, 3, 3)*qp_j(1, 3) + &
1536  tij_abcd(2, 3, 3, 3)*qp_j(2, 3) + &
1537  tij_abcd(3, 3, 3, 3)*qp_j(3, 3))
1538  tmp12 = tmp21
1539  tmp13 = tmp31
1540  tmp23 = tmp32
1541  tmp = tmp11 + tmp12 + tmp13 + &
1542  tmp21 + tmp22 + tmp23 + &
1543  tmp31 + tmp32 + tmp33
1544 
1545  eloc = eloc + fac*tmp
1546  ! Forces on particle i (locally b)
1547  IF (do_forces .OR. do_stress) THEN
1548  DO k = 1, 3
1549  tmp11 = qp_i(1, 1)*(tij_abcde(1, 1, 1, 1, k)*qp_j(1, 1) + &
1550  tij_abcde(2, 1, 1, 1, k)*qp_j(2, 1) + &
1551  tij_abcde(3, 1, 1, 1, k)*qp_j(3, 1) + &
1552  tij_abcde(1, 2, 1, 1, k)*qp_j(1, 2) + &
1553  tij_abcde(2, 2, 1, 1, k)*qp_j(2, 2) + &
1554  tij_abcde(3, 2, 1, 1, k)*qp_j(3, 2) + &
1555  tij_abcde(1, 3, 1, 1, k)*qp_j(1, 3) + &
1556  tij_abcde(2, 3, 1, 1, k)*qp_j(2, 3) + &
1557  tij_abcde(3, 3, 1, 1, k)*qp_j(3, 3))
1558  tmp21 = qp_i(2, 1)*(tij_abcde(1, 1, 2, 1, k)*qp_j(1, 1) + &
1559  tij_abcde(2, 1, 2, 1, k)*qp_j(2, 1) + &
1560  tij_abcde(3, 1, 2, 1, k)*qp_j(3, 1) + &
1561  tij_abcde(1, 2, 2, 1, k)*qp_j(1, 2) + &
1562  tij_abcde(2, 2, 2, 1, k)*qp_j(2, 2) + &
1563  tij_abcde(3, 2, 2, 1, k)*qp_j(3, 2) + &
1564  tij_abcde(1, 3, 2, 1, k)*qp_j(1, 3) + &
1565  tij_abcde(2, 3, 2, 1, k)*qp_j(2, 3) + &
1566  tij_abcde(3, 3, 2, 1, k)*qp_j(3, 3))
1567  tmp31 = qp_i(3, 1)*(tij_abcde(1, 1, 3, 1, k)*qp_j(1, 1) + &
1568  tij_abcde(2, 1, 3, 1, k)*qp_j(2, 1) + &
1569  tij_abcde(3, 1, 3, 1, k)*qp_j(3, 1) + &
1570  tij_abcde(1, 2, 3, 1, k)*qp_j(1, 2) + &
1571  tij_abcde(2, 2, 3, 1, k)*qp_j(2, 2) + &
1572  tij_abcde(3, 2, 3, 1, k)*qp_j(3, 2) + &
1573  tij_abcde(1, 3, 3, 1, k)*qp_j(1, 3) + &
1574  tij_abcde(2, 3, 3, 1, k)*qp_j(2, 3) + &
1575  tij_abcde(3, 3, 3, 1, k)*qp_j(3, 3))
1576  tmp22 = qp_i(2, 2)*(tij_abcde(1, 1, 2, 2, k)*qp_j(1, 1) + &
1577  tij_abcde(2, 1, 2, 2, k)*qp_j(2, 1) + &
1578  tij_abcde(3, 1, 2, 2, k)*qp_j(3, 1) + &
1579  tij_abcde(1, 2, 2, 2, k)*qp_j(1, 2) + &
1580  tij_abcde(2, 2, 2, 2, k)*qp_j(2, 2) + &
1581  tij_abcde(3, 2, 2, 2, k)*qp_j(3, 2) + &
1582  tij_abcde(1, 3, 2, 2, k)*qp_j(1, 3) + &
1583  tij_abcde(2, 3, 2, 2, k)*qp_j(2, 3) + &
1584  tij_abcde(3, 3, 2, 2, k)*qp_j(3, 3))
1585  tmp32 = qp_i(3, 2)*(tij_abcde(1, 1, 3, 2, k)*qp_j(1, 1) + &
1586  tij_abcde(2, 1, 3, 2, k)*qp_j(2, 1) + &
1587  tij_abcde(3, 1, 3, 2, k)*qp_j(3, 1) + &
1588  tij_abcde(1, 2, 3, 2, k)*qp_j(1, 2) + &
1589  tij_abcde(2, 2, 3, 2, k)*qp_j(2, 2) + &
1590  tij_abcde(3, 2, 3, 2, k)*qp_j(3, 2) + &
1591  tij_abcde(1, 3, 3, 2, k)*qp_j(1, 3) + &
1592  tij_abcde(2, 3, 3, 2, k)*qp_j(2, 3) + &
1593  tij_abcde(3, 3, 3, 2, k)*qp_j(3, 3))
1594  tmp33 = qp_i(3, 3)*(tij_abcde(1, 1, 3, 3, k)*qp_j(1, 1) + &
1595  tij_abcde(2, 1, 3, 3, k)*qp_j(2, 1) + &
1596  tij_abcde(3, 1, 3, 3, k)*qp_j(3, 1) + &
1597  tij_abcde(1, 2, 3, 3, k)*qp_j(1, 2) + &
1598  tij_abcde(2, 2, 3, 3, k)*qp_j(2, 2) + &
1599  tij_abcde(3, 2, 3, 3, k)*qp_j(3, 2) + &
1600  tij_abcde(1, 3, 3, 3, k)*qp_j(1, 3) + &
1601  tij_abcde(2, 3, 3, 3, k)*qp_j(2, 3) + &
1602  tij_abcde(3, 3, 3, 3, k)*qp_j(3, 3))
1603  tmp12 = tmp21
1604  tmp13 = tmp31
1605  tmp23 = tmp32
1606  fr(k) = fr(k) - fac*(tmp11 + tmp12 + tmp13 + &
1607  tmp21 + tmp22 + tmp23 + &
1608  tmp31 + tmp32 + tmp33)
1609  END DO
1610  END IF
1611  ! Electric field
1612  IF (do_efield) THEN
1613  fac = 1.0_dp/3.0_dp
1614  ! Potential
1615  IF (do_efield0) THEN
1616  ef0_i = ef0_i + fac*(tij_ab(1, 1)*qp_j(1, 1) + &
1617  tij_ab(2, 1)*qp_j(2, 1) + &
1618  tij_ab(3, 1)*qp_j(3, 1) + &
1619  tij_ab(1, 2)*qp_j(1, 2) + &
1620  tij_ab(2, 2)*qp_j(2, 2) + &
1621  tij_ab(3, 2)*qp_j(3, 2) + &
1622  tij_ab(1, 3)*qp_j(1, 3) + &
1623  tij_ab(2, 3)*qp_j(2, 3) + &
1624  tij_ab(3, 3)*qp_j(3, 3))
1625 
1626  ef0_j = ef0_j + fac*(tij_ab(1, 1)*qp_i(1, 1) + &
1627  tij_ab(2, 1)*qp_i(2, 1) + &
1628  tij_ab(3, 1)*qp_i(3, 1) + &
1629  tij_ab(1, 2)*qp_i(1, 2) + &
1630  tij_ab(2, 2)*qp_i(2, 2) + &
1631  tij_ab(3, 2)*qp_i(3, 2) + &
1632  tij_ab(1, 3)*qp_i(1, 3) + &
1633  tij_ab(2, 3)*qp_i(2, 3) + &
1634  tij_ab(3, 3)*qp_i(3, 3))
1635  END IF
1636  ! Electric field
1637  IF (do_efield1) THEN
1638  ef1_i(1) = ef1_i(1) - fac*(tij_abc(1, 1, 1)*qp_j(1, 1) + &
1639  tij_abc(2, 1, 1)*qp_j(2, 1) + &
1640  tij_abc(3, 1, 1)*qp_j(3, 1) + &
1641  tij_abc(1, 2, 1)*qp_j(1, 2) + &
1642  tij_abc(2, 2, 1)*qp_j(2, 2) + &
1643  tij_abc(3, 2, 1)*qp_j(3, 2) + &
1644  tij_abc(1, 3, 1)*qp_j(1, 3) + &
1645  tij_abc(2, 3, 1)*qp_j(2, 3) + &
1646  tij_abc(3, 3, 1)*qp_j(3, 3))
1647  ef1_i(2) = ef1_i(2) - fac*(tij_abc(1, 1, 2)*qp_j(1, 1) + &
1648  tij_abc(2, 1, 2)*qp_j(2, 1) + &
1649  tij_abc(3, 1, 2)*qp_j(3, 1) + &
1650  tij_abc(1, 2, 2)*qp_j(1, 2) + &
1651  tij_abc(2, 2, 2)*qp_j(2, 2) + &
1652  tij_abc(3, 2, 2)*qp_j(3, 2) + &
1653  tij_abc(1, 3, 2)*qp_j(1, 3) + &
1654  tij_abc(2, 3, 2)*qp_j(2, 3) + &
1655  tij_abc(3, 3, 2)*qp_j(3, 3))
1656  ef1_i(3) = ef1_i(3) - fac*(tij_abc(1, 1, 3)*qp_j(1, 1) + &
1657  tij_abc(2, 1, 3)*qp_j(2, 1) + &
1658  tij_abc(3, 1, 3)*qp_j(3, 1) + &
1659  tij_abc(1, 2, 3)*qp_j(1, 2) + &
1660  tij_abc(2, 2, 3)*qp_j(2, 2) + &
1661  tij_abc(3, 2, 3)*qp_j(3, 2) + &
1662  tij_abc(1, 3, 3)*qp_j(1, 3) + &
1663  tij_abc(2, 3, 3)*qp_j(2, 3) + &
1664  tij_abc(3, 3, 3)*qp_j(3, 3))
1665 
1666  ef1_j(1) = ef1_j(1) + fac*(tij_abc(1, 1, 1)*qp_i(1, 1) + &
1667  tij_abc(2, 1, 1)*qp_i(2, 1) + &
1668  tij_abc(3, 1, 1)*qp_i(3, 1) + &
1669  tij_abc(1, 2, 1)*qp_i(1, 2) + &
1670  tij_abc(2, 2, 1)*qp_i(2, 2) + &
1671  tij_abc(3, 2, 1)*qp_i(3, 2) + &
1672  tij_abc(1, 3, 1)*qp_i(1, 3) + &
1673  tij_abc(2, 3, 1)*qp_i(2, 3) + &
1674  tij_abc(3, 3, 1)*qp_i(3, 3))
1675  ef1_j(2) = ef1_j(2) + fac*(tij_abc(1, 1, 2)*qp_i(1, 1) + &
1676  tij_abc(2, 1, 2)*qp_i(2, 1) + &
1677  tij_abc(3, 1, 2)*qp_i(3, 1) + &
1678  tij_abc(1, 2, 2)*qp_i(1, 2) + &
1679  tij_abc(2, 2, 2)*qp_i(2, 2) + &
1680  tij_abc(3, 2, 2)*qp_i(3, 2) + &
1681  tij_abc(1, 3, 2)*qp_i(1, 3) + &
1682  tij_abc(2, 3, 2)*qp_i(2, 3) + &
1683  tij_abc(3, 3, 2)*qp_i(3, 3))
1684  ef1_j(3) = ef1_j(3) + fac*(tij_abc(1, 1, 3)*qp_i(1, 1) + &
1685  tij_abc(2, 1, 3)*qp_i(2, 1) + &
1686  tij_abc(3, 1, 3)*qp_i(3, 1) + &
1687  tij_abc(1, 2, 3)*qp_i(1, 2) + &
1688  tij_abc(2, 2, 3)*qp_i(2, 2) + &
1689  tij_abc(3, 2, 3)*qp_i(3, 2) + &
1690  tij_abc(1, 3, 3)*qp_i(1, 3) + &
1691  tij_abc(2, 3, 3)*qp_i(2, 3) + &
1692  tij_abc(3, 3, 3)*qp_i(3, 3))
1693  END IF
1694  ! Electric field gradient
1695  IF (do_efield2) THEN
1696  tmp11 = fac*(tij_abcd(1, 1, 1, 1)*qp_j(1, 1) + &
1697  tij_abcd(2, 1, 1, 1)*qp_j(2, 1) + &
1698  tij_abcd(3, 1, 1, 1)*qp_j(3, 1) + &
1699  tij_abcd(1, 2, 1, 1)*qp_j(1, 2) + &
1700  tij_abcd(2, 2, 1, 1)*qp_j(2, 2) + &
1701  tij_abcd(3, 2, 1, 1)*qp_j(3, 2) + &
1702  tij_abcd(1, 3, 1, 1)*qp_j(1, 3) + &
1703  tij_abcd(2, 3, 1, 1)*qp_j(2, 3) + &
1704  tij_abcd(3, 3, 1, 1)*qp_j(3, 3))
1705  tmp12 = fac*(tij_abcd(1, 1, 1, 2)*qp_j(1, 1) + &
1706  tij_abcd(2, 1, 1, 2)*qp_j(2, 1) + &
1707  tij_abcd(3, 1, 1, 2)*qp_j(3, 1) + &
1708  tij_abcd(1, 2, 1, 2)*qp_j(1, 2) + &
1709  tij_abcd(2, 2, 1, 2)*qp_j(2, 2) + &
1710  tij_abcd(3, 2, 1, 2)*qp_j(3, 2) + &
1711  tij_abcd(1, 3, 1, 2)*qp_j(1, 3) + &
1712  tij_abcd(2, 3, 1, 2)*qp_j(2, 3) + &
1713  tij_abcd(3, 3, 1, 2)*qp_j(3, 3))
1714  tmp13 = fac*(tij_abcd(1, 1, 1, 3)*qp_j(1, 1) + &
1715  tij_abcd(2, 1, 1, 3)*qp_j(2, 1) + &
1716  tij_abcd(3, 1, 1, 3)*qp_j(3, 1) + &
1717  tij_abcd(1, 2, 1, 3)*qp_j(1, 2) + &
1718  tij_abcd(2, 2, 1, 3)*qp_j(2, 2) + &
1719  tij_abcd(3, 2, 1, 3)*qp_j(3, 2) + &
1720  tij_abcd(1, 3, 1, 3)*qp_j(1, 3) + &
1721  tij_abcd(2, 3, 1, 3)*qp_j(2, 3) + &
1722  tij_abcd(3, 3, 1, 3)*qp_j(3, 3))
1723  tmp22 = fac*(tij_abcd(1, 1, 2, 2)*qp_j(1, 1) + &
1724  tij_abcd(2, 1, 2, 2)*qp_j(2, 1) + &
1725  tij_abcd(3, 1, 2, 2)*qp_j(3, 1) + &
1726  tij_abcd(1, 2, 2, 2)*qp_j(1, 2) + &
1727  tij_abcd(2, 2, 2, 2)*qp_j(2, 2) + &
1728  tij_abcd(3, 2, 2, 2)*qp_j(3, 2) + &
1729  tij_abcd(1, 3, 2, 2)*qp_j(1, 3) + &
1730  tij_abcd(2, 3, 2, 2)*qp_j(2, 3) + &
1731  tij_abcd(3, 3, 2, 2)*qp_j(3, 3))
1732  tmp23 = fac*(tij_abcd(1, 1, 2, 3)*qp_j(1, 1) + &
1733  tij_abcd(2, 1, 2, 3)*qp_j(2, 1) + &
1734  tij_abcd(3, 1, 2, 3)*qp_j(3, 1) + &
1735  tij_abcd(1, 2, 2, 3)*qp_j(1, 2) + &
1736  tij_abcd(2, 2, 2, 3)*qp_j(2, 2) + &
1737  tij_abcd(3, 2, 2, 3)*qp_j(3, 2) + &
1738  tij_abcd(1, 3, 2, 3)*qp_j(1, 3) + &
1739  tij_abcd(2, 3, 2, 3)*qp_j(2, 3) + &
1740  tij_abcd(3, 3, 2, 3)*qp_j(3, 3))
1741  tmp33 = fac*(tij_abcd(1, 1, 3, 3)*qp_j(1, 1) + &
1742  tij_abcd(2, 1, 3, 3)*qp_j(2, 1) + &
1743  tij_abcd(3, 1, 3, 3)*qp_j(3, 1) + &
1744  tij_abcd(1, 2, 3, 3)*qp_j(1, 2) + &
1745  tij_abcd(2, 2, 3, 3)*qp_j(2, 2) + &
1746  tij_abcd(3, 2, 3, 3)*qp_j(3, 2) + &
1747  tij_abcd(1, 3, 3, 3)*qp_j(1, 3) + &
1748  tij_abcd(2, 3, 3, 3)*qp_j(2, 3) + &
1749  tij_abcd(3, 3, 3, 3)*qp_j(3, 3))
1750 
1751  ef2_i(1, 1) = ef2_i(1, 1) - tmp11
1752  ef2_i(1, 2) = ef2_i(1, 2) - tmp12
1753  ef2_i(1, 3) = ef2_i(1, 3) - tmp13
1754  ef2_i(2, 1) = ef2_i(2, 1) - tmp12
1755  ef2_i(2, 2) = ef2_i(2, 2) - tmp22
1756  ef2_i(2, 3) = ef2_i(2, 3) - tmp23
1757  ef2_i(3, 1) = ef2_i(3, 1) - tmp13
1758  ef2_i(3, 2) = ef2_i(3, 2) - tmp23
1759  ef2_i(3, 3) = ef2_i(3, 3) - tmp33
1760 
1761  tmp11 = fac*(tij_abcd(1, 1, 1, 1)*qp_i(1, 1) + &
1762  tij_abcd(2, 1, 1, 1)*qp_i(2, 1) + &
1763  tij_abcd(3, 1, 1, 1)*qp_i(3, 1) + &
1764  tij_abcd(1, 2, 1, 1)*qp_i(1, 2) + &
1765  tij_abcd(2, 2, 1, 1)*qp_i(2, 2) + &
1766  tij_abcd(3, 2, 1, 1)*qp_i(3, 2) + &
1767  tij_abcd(1, 3, 1, 1)*qp_i(1, 3) + &
1768  tij_abcd(2, 3, 1, 1)*qp_i(2, 3) + &
1769  tij_abcd(3, 3, 1, 1)*qp_i(3, 3))
1770  tmp12 = fac*(tij_abcd(1, 1, 1, 2)*qp_i(1, 1) + &
1771  tij_abcd(2, 1, 1, 2)*qp_i(2, 1) + &
1772  tij_abcd(3, 1, 1, 2)*qp_i(3, 1) + &
1773  tij_abcd(1, 2, 1, 2)*qp_i(1, 2) + &
1774  tij_abcd(2, 2, 1, 2)*qp_i(2, 2) + &
1775  tij_abcd(3, 2, 1, 2)*qp_i(3, 2) + &
1776  tij_abcd(1, 3, 1, 2)*qp_i(1, 3) + &
1777  tij_abcd(2, 3, 1, 2)*qp_i(2, 3) + &
1778  tij_abcd(3, 3, 1, 2)*qp_i(3, 3))
1779  tmp13 = fac*(tij_abcd(1, 1, 1, 3)*qp_i(1, 1) + &
1780  tij_abcd(2, 1, 1, 3)*qp_i(2, 1) + &
1781  tij_abcd(3, 1, 1, 3)*qp_i(3, 1) + &
1782  tij_abcd(1, 2, 1, 3)*qp_i(1, 2) + &
1783  tij_abcd(2, 2, 1, 3)*qp_i(2, 2) + &
1784  tij_abcd(3, 2, 1, 3)*qp_i(3, 2) + &
1785  tij_abcd(1, 3, 1, 3)*qp_i(1, 3) + &
1786  tij_abcd(2, 3, 1, 3)*qp_i(2, 3) + &
1787  tij_abcd(3, 3, 1, 3)*qp_i(3, 3))
1788  tmp22 = fac*(tij_abcd(1, 1, 2, 2)*qp_i(1, 1) + &
1789  tij_abcd(2, 1, 2, 2)*qp_i(2, 1) + &
1790  tij_abcd(3, 1, 2, 2)*qp_i(3, 1) + &
1791  tij_abcd(1, 2, 2, 2)*qp_i(1, 2) + &
1792  tij_abcd(2, 2, 2, 2)*qp_i(2, 2) + &
1793  tij_abcd(3, 2, 2, 2)*qp_i(3, 2) + &
1794  tij_abcd(1, 3, 2, 2)*qp_i(1, 3) + &
1795  tij_abcd(2, 3, 2, 2)*qp_i(2, 3) + &
1796  tij_abcd(3, 3, 2, 2)*qp_i(3, 3))
1797  tmp23 = fac*(tij_abcd(1, 1, 2, 3)*qp_i(1, 1) + &
1798  tij_abcd(2, 1, 2, 3)*qp_i(2, 1) + &
1799  tij_abcd(3, 1, 2, 3)*qp_i(3, 1) + &
1800  tij_abcd(1, 2, 2, 3)*qp_i(1, 2) + &
1801  tij_abcd(2, 2, 2, 3)*qp_i(2, 2) + &
1802  tij_abcd(3, 2, 2, 3)*qp_i(3, 2) + &
1803  tij_abcd(1, 3, 2, 3)*qp_i(1, 3) + &
1804  tij_abcd(2, 3, 2, 3)*qp_i(2, 3) + &
1805  tij_abcd(3, 3, 2, 3)*qp_i(3, 3))
1806  tmp33 = fac*(tij_abcd(1, 1, 3, 3)*qp_i(1, 1) + &
1807  tij_abcd(2, 1, 3, 3)*qp_i(2, 1) + &
1808  tij_abcd(3, 1, 3, 3)*qp_i(3, 1) + &
1809  tij_abcd(1, 2, 3, 3)*qp_i(1, 2) + &
1810  tij_abcd(2, 2, 3, 3)*qp_i(2, 2) + &
1811  tij_abcd(3, 2, 3, 3)*qp_i(3, 2) + &
1812  tij_abcd(1, 3, 3, 3)*qp_i(1, 3) + &
1813  tij_abcd(2, 3, 3, 3)*qp_i(2, 3) + &
1814  tij_abcd(3, 3, 3, 3)*qp_i(3, 3))
1815 
1816  ef2_j(1, 1) = ef2_j(1, 1) - tmp11
1817  ef2_j(1, 2) = ef2_j(1, 2) - tmp12
1818  ef2_j(1, 3) = ef2_j(1, 3) - tmp13
1819  ef2_j(2, 1) = ef2_j(2, 1) - tmp12
1820  ef2_j(2, 2) = ef2_j(2, 2) - tmp22
1821  ef2_j(2, 3) = ef2_j(2, 3) - tmp23
1822  ef2_j(3, 1) = ef2_j(3, 1) - tmp13
1823  ef2_j(3, 2) = ef2_j(3, 2) - tmp23
1824  ef2_j(3, 3) = ef2_j(3, 3) - tmp33
1825  END IF
1826  END IF
1827  END IF
1828  IF (task(3, 2)) THEN
1829  ! Quadrupole - Dipole
1830  fac = 1.0_dp/3.0_dp
1831  ! Dipole i (locally B) - Quadrupole j (locally A)
1832  tmp_ij = dp_i(1)*(tij_abc(1, 1, 1)*qp_j(1, 1) + &
1833  tij_abc(2, 1, 1)*qp_j(2, 1) + &
1834  tij_abc(3, 1, 1)*qp_j(3, 1) + &
1835  tij_abc(1, 2, 1)*qp_j(1, 2) + &
1836  tij_abc(2, 2, 1)*qp_j(2, 2) + &
1837  tij_abc(3, 2, 1)*qp_j(3, 2) + &
1838  tij_abc(1, 3, 1)*qp_j(1, 3) + &
1839  tij_abc(2, 3, 1)*qp_j(2, 3) + &
1840  tij_abc(3, 3, 1)*qp_j(3, 3)) + &
1841  dp_i(2)*(tij_abc(1, 1, 2)*qp_j(1, 1) + &
1842  tij_abc(2, 1, 2)*qp_j(2, 1) + &
1843  tij_abc(3, 1, 2)*qp_j(3, 1) + &
1844  tij_abc(1, 2, 2)*qp_j(1, 2) + &
1845  tij_abc(2, 2, 2)*qp_j(2, 2) + &
1846  tij_abc(3, 2, 2)*qp_j(3, 2) + &
1847  tij_abc(1, 3, 2)*qp_j(1, 3) + &
1848  tij_abc(2, 3, 2)*qp_j(2, 3) + &
1849  tij_abc(3, 3, 2)*qp_j(3, 3)) + &
1850  dp_i(3)*(tij_abc(1, 1, 3)*qp_j(1, 1) + &
1851  tij_abc(2, 1, 3)*qp_j(2, 1) + &
1852  tij_abc(3, 1, 3)*qp_j(3, 1) + &
1853  tij_abc(1, 2, 3)*qp_j(1, 2) + &
1854  tij_abc(2, 2, 3)*qp_j(2, 2) + &
1855  tij_abc(3, 2, 3)*qp_j(3, 2) + &
1856  tij_abc(1, 3, 3)*qp_j(1, 3) + &
1857  tij_abc(2, 3, 3)*qp_j(2, 3) + &
1858  tij_abc(3, 3, 3)*qp_j(3, 3))
1859 
1860  ! Dipole j (locally A) - Quadrupole i (locally B)
1861  tmp_ji = dp_j(1)*(tij_abc(1, 1, 1)*qp_i(1, 1) + &
1862  tij_abc(2, 1, 1)*qp_i(2, 1) + &
1863  tij_abc(3, 1, 1)*qp_i(3, 1) + &
1864  tij_abc(1, 2, 1)*qp_i(1, 2) + &
1865  tij_abc(2, 2, 1)*qp_i(2, 2) + &
1866  tij_abc(3, 2, 1)*qp_i(3, 2) + &
1867  tij_abc(1, 3, 1)*qp_i(1, 3) + &
1868  tij_abc(2, 3, 1)*qp_i(2, 3) + &
1869  tij_abc(3, 3, 1)*qp_i(3, 3)) + &
1870  dp_j(2)*(tij_abc(1, 1, 2)*qp_i(1, 1) + &
1871  tij_abc(2, 1, 2)*qp_i(2, 1) + &
1872  tij_abc(3, 1, 2)*qp_i(3, 1) + &
1873  tij_abc(1, 2, 2)*qp_i(1, 2) + &
1874  tij_abc(2, 2, 2)*qp_i(2, 2) + &
1875  tij_abc(3, 2, 2)*qp_i(3, 2) + &
1876  tij_abc(1, 3, 2)*qp_i(1, 3) + &
1877  tij_abc(2, 3, 2)*qp_i(2, 3) + &
1878  tij_abc(3, 3, 2)*qp_i(3, 3)) + &
1879  dp_j(3)*(tij_abc(1, 1, 3)*qp_i(1, 1) + &
1880  tij_abc(2, 1, 3)*qp_i(2, 1) + &
1881  tij_abc(3, 1, 3)*qp_i(3, 1) + &
1882  tij_abc(1, 2, 3)*qp_i(1, 2) + &
1883  tij_abc(2, 2, 3)*qp_i(2, 2) + &
1884  tij_abc(3, 2, 3)*qp_i(3, 2) + &
1885  tij_abc(1, 3, 3)*qp_i(1, 3) + &
1886  tij_abc(2, 3, 3)*qp_i(2, 3) + &
1887  tij_abc(3, 3, 3)*qp_i(3, 3))
1888 
1889  tmp = fac*(tmp_ij - tmp_ji)
1890  eloc = eloc + tmp
1891  IF (do_forces .OR. do_stress) THEN
1892  DO k = 1, 3
1893  ! Dipole i (locally B) - Quadrupole j (locally A)
1894  tmp_ij = dp_i(1)*(tij_abcd(1, 1, 1, k)*qp_j(1, 1) + &
1895  tij_abcd(2, 1, 1, k)*qp_j(2, 1) + &
1896  tij_abcd(3, 1, 1, k)*qp_j(3, 1) + &
1897  tij_abcd(1, 2, 1, k)*qp_j(1, 2) + &
1898  tij_abcd(2, 2, 1, k)*qp_j(2, 2) + &
1899  tij_abcd(3, 2, 1, k)*qp_j(3, 2) + &
1900  tij_abcd(1, 3, 1, k)*qp_j(1, 3) + &
1901  tij_abcd(2, 3, 1, k)*qp_j(2, 3) + &
1902  tij_abcd(3, 3, 1, k)*qp_j(3, 3)) + &
1903  dp_i(2)*(tij_abcd(1, 1, 2, k)*qp_j(1, 1) + &
1904  tij_abcd(2, 1, 2, k)*qp_j(2, 1) + &
1905  tij_abcd(3, 1, 2, k)*qp_j(3, 1) + &
1906  tij_abcd(1, 2, 2, k)*qp_j(1, 2) + &
1907  tij_abcd(2, 2, 2, k)*qp_j(2, 2) + &
1908  tij_abcd(3, 2, 2, k)*qp_j(3, 2) + &
1909  tij_abcd(1, 3, 2, k)*qp_j(1, 3) + &
1910  tij_abcd(2, 3, 2, k)*qp_j(2, 3) + &
1911  tij_abcd(3, 3, 2, k)*qp_j(3, 3)) + &
1912  dp_i(3)*(tij_abcd(1, 1, 3, k)*qp_j(1, 1) + &
1913  tij_abcd(2, 1, 3, k)*qp_j(2, 1) + &
1914  tij_abcd(3, 1, 3, k)*qp_j(3, 1) + &
1915  tij_abcd(1, 2, 3, k)*qp_j(1, 2) + &
1916  tij_abcd(2, 2, 3, k)*qp_j(2, 2) + &
1917  tij_abcd(3, 2, 3, k)*qp_j(3, 2) + &
1918  tij_abcd(1, 3, 3, k)*qp_j(1, 3) + &
1919  tij_abcd(2, 3, 3, k)*qp_j(2, 3) + &
1920  tij_abcd(3, 3, 3, k)*qp_j(3, 3))
1921 
1922  ! Dipole j (locally A) - Quadrupole i (locally B)
1923  tmp_ji = dp_j(1)*(tij_abcd(1, 1, 1, k)*qp_i(1, 1) + &
1924  tij_abcd(2, 1, 1, k)*qp_i(2, 1) + &
1925  tij_abcd(3, 1, 1, k)*qp_i(3, 1) + &
1926  tij_abcd(1, 2, 1, k)*qp_i(1, 2) + &
1927  tij_abcd(2, 2, 1, k)*qp_i(2, 2) + &
1928  tij_abcd(3, 2, 1, k)*qp_i(3, 2) + &
1929  tij_abcd(1, 3, 1, k)*qp_i(1, 3) + &
1930  tij_abcd(2, 3, 1, k)*qp_i(2, 3) + &
1931  tij_abcd(3, 3, 1, k)*qp_i(3, 3)) + &
1932  dp_j(2)*(tij_abcd(1, 1, 2, k)*qp_i(1, 1) + &
1933  tij_abcd(2, 1, 2, k)*qp_i(2, 1) + &
1934  tij_abcd(3, 1, 2, k)*qp_i(3, 1) + &
1935  tij_abcd(1, 2, 2, k)*qp_i(1, 2) + &
1936  tij_abcd(2, 2, 2, k)*qp_i(2, 2) + &
1937  tij_abcd(3, 2, 2, k)*qp_i(3, 2) + &
1938  tij_abcd(1, 3, 2, k)*qp_i(1, 3) + &
1939  tij_abcd(2, 3, 2, k)*qp_i(2, 3) + &
1940  tij_abcd(3, 3, 2, k)*qp_i(3, 3)) + &
1941  dp_j(3)*(tij_abcd(1, 1, 3, k)*qp_i(1, 1) + &
1942  tij_abcd(2, 1, 3, k)*qp_i(2, 1) + &
1943  tij_abcd(3, 1, 3, k)*qp_i(3, 1) + &
1944  tij_abcd(1, 2, 3, k)*qp_i(1, 2) + &
1945  tij_abcd(2, 2, 3, k)*qp_i(2, 2) + &
1946  tij_abcd(3, 2, 3, k)*qp_i(3, 2) + &
1947  tij_abcd(1, 3, 3, k)*qp_i(1, 3) + &
1948  tij_abcd(2, 3, 3, k)*qp_i(2, 3) + &
1949  tij_abcd(3, 3, 3, k)*qp_i(3, 3))
1950 
1951  fr(k) = fr(k) - fac*(tmp_ij - tmp_ji)
1952  END DO
1953  END IF
1954  END IF
1955  IF (task(3, 1)) THEN
1956  ! Quadrupole - Charge
1957  fac = 1.0_dp/3.0_dp
1958 
1959  ! Quadrupole j (locally A) - Charge j (locally B)
1960  tmp_ij = ch_i*(tij_ab(1, 1)*qp_j(1, 1) + &
1961  tij_ab(2, 1)*qp_j(2, 1) + &
1962  tij_ab(3, 1)*qp_j(3, 1) + &
1963  tij_ab(1, 2)*qp_j(1, 2) + &
1964  tij_ab(2, 2)*qp_j(2, 2) + &
1965  tij_ab(3, 2)*qp_j(3, 2) + &
1966  tij_ab(1, 3)*qp_j(1, 3) + &
1967  tij_ab(2, 3)*qp_j(2, 3) + &
1968  tij_ab(3, 3)*qp_j(3, 3))
1969 
1970  ! Quadrupole i (locally B) - Charge j (locally A)
1971  tmp_ji = ch_j*(tij_ab(1, 1)*qp_i(1, 1) + &
1972  tij_ab(2, 1)*qp_i(2, 1) + &
1973  tij_ab(3, 1)*qp_i(3, 1) + &
1974  tij_ab(1, 2)*qp_i(1, 2) + &
1975  tij_ab(2, 2)*qp_i(2, 2) + &
1976  tij_ab(3, 2)*qp_i(3, 2) + &
1977  tij_ab(1, 3)*qp_i(1, 3) + &
1978  tij_ab(2, 3)*qp_i(2, 3) + &
1979  tij_ab(3, 3)*qp_i(3, 3))
1980 
1981  eloc = eloc + fac*(tmp_ij + tmp_ji)
1982  IF (do_forces .OR. do_stress) THEN
1983  DO k = 1, 3
1984  ! Quadrupole j (locally A) - Charge i (locally B)
1985  tmp_ij = ch_i*(tij_abc(1, 1, k)*qp_j(1, 1) + &
1986  tij_abc(2, 1, k)*qp_j(2, 1) + &
1987  tij_abc(3, 1, k)*qp_j(3, 1) + &
1988  tij_abc(1, 2, k)*qp_j(1, 2) + &
1989  tij_abc(2, 2, k)*qp_j(2, 2) + &
1990  tij_abc(3, 2, k)*qp_j(3, 2) + &
1991  tij_abc(1, 3, k)*qp_j(1, 3) + &
1992  tij_abc(2, 3, k)*qp_j(2, 3) + &
1993  tij_abc(3, 3, k)*qp_j(3, 3))
1994 
1995  ! Quadrupole i (locally B) - Charge j (locally A)
1996  tmp_ji = ch_j*(tij_abc(1, 1, k)*qp_i(1, 1) + &
1997  tij_abc(2, 1, k)*qp_i(2, 1) + &
1998  tij_abc(3, 1, k)*qp_i(3, 1) + &
1999  tij_abc(1, 2, k)*qp_i(1, 2) + &
2000  tij_abc(2, 2, k)*qp_i(2, 2) + &
2001  tij_abc(3, 2, k)*qp_i(3, 2) + &
2002  tij_abc(1, 3, k)*qp_i(1, 3) + &
2003  tij_abc(2, 3, k)*qp_i(2, 3) + &
2004  tij_abc(3, 3, k)*qp_i(3, 3))
2005 
2006  fr(k) = fr(k) - fac*(tmp_ij + tmp_ji)
2007  END DO
2008  END IF
2009  END IF
2010  ! Electric fields
2011  IF (do_efield) THEN
2012  ! Potential
2013  IF (do_efield0) THEN
2014  efield0(atom_a) = efield0(atom_a) + ef0_j
2015 
2016  efield0(atom_b) = efield0(atom_b) + ef0_i
2017  END IF
2018  ! Electric field
2019  IF (do_efield1) THEN
2020  efield1(1, atom_a) = efield1(1, atom_a) + ef1_j(1)
2021  efield1(2, atom_a) = efield1(2, atom_a) + ef1_j(2)
2022  efield1(3, atom_a) = efield1(3, atom_a) + ef1_j(3)
2023 
2024  efield1(1, atom_b) = efield1(1, atom_b) + ef1_i(1)
2025  efield1(2, atom_b) = efield1(2, atom_b) + ef1_i(2)
2026  efield1(3, atom_b) = efield1(3, atom_b) + ef1_i(3)
2027  END IF
2028  ! Electric field gradient
2029  IF (do_efield2) THEN
2030  efield2(1, atom_a) = efield2(1, atom_a) + ef2_j(1, 1)
2031  efield2(2, atom_a) = efield2(2, atom_a) + ef2_j(1, 2)
2032  efield2(3, atom_a) = efield2(3, atom_a) + ef2_j(1, 3)
2033  efield2(4, atom_a) = efield2(4, atom_a) + ef2_j(2, 1)
2034  efield2(5, atom_a) = efield2(5, atom_a) + ef2_j(2, 2)
2035  efield2(6, atom_a) = efield2(6, atom_a) + ef2_j(2, 3)
2036  efield2(7, atom_a) = efield2(7, atom_a) + ef2_j(3, 1)
2037  efield2(8, atom_a) = efield2(8, atom_a) + ef2_j(3, 2)
2038  efield2(9, atom_a) = efield2(9, atom_a) + ef2_j(3, 3)
2039 
2040  efield2(1, atom_b) = efield2(1, atom_b) + ef2_i(1, 1)
2041  efield2(2, atom_b) = efield2(2, atom_b) + ef2_i(1, 2)
2042  efield2(3, atom_b) = efield2(3, atom_b) + ef2_i(1, 3)
2043  efield2(4, atom_b) = efield2(4, atom_b) + ef2_i(2, 1)
2044  efield2(5, atom_b) = efield2(5, atom_b) + ef2_i(2, 2)
2045  efield2(6, atom_b) = efield2(6, atom_b) + ef2_i(2, 3)
2046  efield2(7, atom_b) = efield2(7, atom_b) + ef2_i(3, 1)
2047  efield2(8, atom_b) = efield2(8, atom_b) + ef2_i(3, 2)
2048  efield2(9, atom_b) = efield2(9, atom_b) + ef2_i(3, 3)
2049  END IF
2050  END IF
2051  IF (do_stress) THEN
2052  ptens11 = ptens11 + rab(1)*fr(1)
2053  ptens21 = ptens21 + rab(2)*fr(1)
2054  ptens31 = ptens31 + rab(3)*fr(1)
2055  ptens12 = ptens12 + rab(1)*fr(2)
2056  ptens22 = ptens22 + rab(2)*fr(2)
2057  ptens32 = ptens32 + rab(3)*fr(2)
2058  ptens13 = ptens13 + rab(1)*fr(3)
2059  ptens23 = ptens23 + rab(2)*fr(3)
2060  ptens33 = ptens33 + rab(3)*fr(3)
2061  END IF
2062 
2063 
2064  IF (PRESENT(integral_value)) THEN
2065  integral_value = eloc
2066  END IF
2067  IF (do_forces) THEN
2068  force_ab = fr
2069  END IF
2070  END SUBROUTINE se_coulomb_ij_interaction
2071 
2072 END MODULE se_fock_matrix_integrals
2073 
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition: grid_common.h:48
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
se_fock_matrix_integrals
Provides the low level routines to build both the exchange and the Coulomb Fock matrices....
Definition: se_fock_matrix_integrals.F:18
se_fock_matrix_integrals::fock2c_r3
subroutine, public fock2c_r3(sepi, sepj, switch, pi_tot, fi_mat, pj_tot, fj_mat, factor, w, rp)
Construction of 2-center Fock Matrix - Coulomb Terms for the residual (1/R^3) integral part.
Definition: se_fock_matrix_integrals.F:889
se_fock_matrix_integrals::fock2_1el_ew
subroutine, public fock2_1el_ew(sep, rij, ks_block, p_block, ecore, itype, anag, se_int_control, se_taper, store_int_env)
Construction of 2-center 1-electron Fock Matrix (Ewald self term)
Definition: se_fock_matrix_integrals.F:289
se_fock_matrix_integrals::dfock2c
subroutine, public dfock2c(sepi, sepj, rij, switch, pi_tot, pj_tot, factor, anag, se_int_control, se_taper, force, delta)
Derivatives of 2-center Fock Matrix - Coulomb Terms.
Definition: se_fock_matrix_integrals.F:507
se_fock_matrix_integrals::fock2_1el_r3
subroutine, public fock2_1el_r3(sepi, sepj, ksi_block, ksj_block, pi_block, pj_block, e1b, e2a, ecore, rp)
Construction of 2-center 1-electron Fock Matrix for the residual (1/R^3) integral part.
Definition: se_fock_matrix_integrals.F:780
se_fock_matrix_integrals::dfock2c_r3
subroutine, public dfock2c_r3(sepi, sepj, switch, pi_tot, pj_tot, factor, w, drp, force)
Derivatives of 2-center Fock Matrix - Coulomb Terms for the residual (1/R^3) integral part.
Definition: se_fock_matrix_integrals.F:958
se_fock_matrix_integrals::dfock2e
subroutine, public dfock2e(sepi, sepj, rij, switch, isize, pi_mat, factor, anag, se_int_control, se_taper, force, delta)
Derivatives of 2-center Fock Matrix - General Driver.
Definition: se_fock_matrix_integrals.F:679
se_fock_matrix_integrals::fock2c
subroutine, public fock2c(sepi, sepj, rij, switch, pi_tot, fi_mat, pj_tot, fj_mat, factor, anag, se_int_control, se_taper, store_int_env)
Construction of 2-center Fock Matrix - Coulomb Terms.
Definition: se_fock_matrix_integrals.F:417
se_fock_matrix_integrals::fock1_2el
subroutine, public fock1_2el(sep, p_tot, p_mat, f_mat, factor)
Construction of 1-center 2-electron Fock Matrix.
Definition: se_fock_matrix_integrals.F:218
se_fock_matrix_integrals::se_coulomb_ij_interaction
subroutine, public se_coulomb_ij_interaction(atom_a, atom_b, my_task, do_forces, do_efield, do_stress, charges, dipoles, quadrupoles, force_ab, efield0, efield1, efield2, rab2, rab, integral_value, ptens11, ptens12, ptens13, ptens21, ptens22, ptens23, ptens31, ptens32, ptens33)
Coulomb interaction multipolar correction.
Definition: se_fock_matrix_integrals.F:1046
se_fock_matrix_integrals::dfock2_1el_r3
subroutine, public dfock2_1el_r3(sepi, sepj, drp, pi_block, pj_block, force, e1b, e2a)
Derivatives of 2-center 1-electron Fock Matrix residual (1/R^3) integral part.
Definition: se_fock_matrix_integrals.F:831
se_fock_matrix_integrals::fock2_1el
subroutine, public fock2_1el(sepi, sepj, rij, ksi_block, ksj_block, pi_block, pj_block, ecore, itype, anag, se_int_control, se_taper, store_int_env)
Construction of 2-center 1-electron Fock Matrix.
Definition: se_fock_matrix_integrals.F:65
se_fock_matrix_integrals::fock2e
subroutine, public fock2e(sepi, sepj, rij, switch, isize, pi_mat, fi_mat, factor, anag, se_int_control, se_taper, store_int_env)
Construction of 2-center Fock Matrix - General Driver.
Definition: se_fock_matrix_integrals.F:595
se_fock_matrix_integrals::dfock2_1el
subroutine, public dfock2_1el(sepi, sepj, rij, pi_block, pj_block, itype, anag, se_int_control, se_taper, force, delta)
Derivatives of 2-center 1-electron Fock Matrix.
Definition: se_fock_matrix_integrals.F:143
se_fock_matrix_integrals::fock2c_ew
subroutine, public fock2c_ew(sep, rij, p_tot, f_mat, factor, anag, se_int_control, se_taper, store_int_env)
Construction of 2-center Fock Matrix - Coulomb Self Terms (Ewald)
Definition: se_fock_matrix_integrals.F:347
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::se_orbital_pointer
integer, dimension(9), public se_orbital_pointer
Definition: semi_empirical_int_arrays.F:30
semi_empirical_integrals
Set of wrappers for semi-empirical analytical/numerical Integrals routines.
Definition: semi_empirical_integrals.F:16
semi_empirical_integrals::rotnuc
subroutine, public rotnuc(sepi, sepj, rij, e1b, e2a, itype, anag, se_int_control, se_taper, store_int_env)
wrapper for numerical/analytical 1 center 1 electron integrals
Definition: semi_empirical_integrals.F:222
semi_empirical_integrals::drotint
subroutine, public drotint(sepi, sepj, rij, dw, delta, anag, se_int_control, se_taper)
wrapper for numerical/analytical routines
Definition: semi_empirical_integrals.F:470
semi_empirical_integrals::rotint
subroutine, public rotint(sepi, sepj, rij, w, anag, se_int_control, se_taper, store_int_env)
wrapper for numerical/analytical 2 center 2 electrons integrals routines with possibility of incore s...
Definition: semi_empirical_integrals.F:72
semi_empirical_integrals::drotnuc
subroutine, public drotnuc(sepi, sepj, rij, de1b, de2a, itype, delta, anag, se_int_control, se_taper)
wrapper for numerical/analytical routines
Definition: semi_empirical_integrals.F:507
semi_empirical_store_int_types
Type to store integrals for semi-empirical calculations.
Definition: semi_empirical_store_int_types.F:13
semi_empirical_types
Definition of the semi empirical parameter types.
Definition: semi_empirical_types.F:12