(git:ed6f26b)
Loading...
Searching...
No Matches
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-2025 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! **************************************************************************************************
18MODULE 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
43CONTAINS
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 f = huge(0.0_dp)
1106 tij = huge(0.0_dp)
1107 tij_a = huge(0.0_dp)
1108 tij_ab = huge(0.0_dp)
1109 tij_abc = huge(0.0_dp)
1110 tij_abcd = huge(0.0_dp)
1111 tij_abcde = huge(0.0_dp)
1112 r = sqrt(rab2)
1113 irab2 = 1.0_dp/rab2
1114 ir = 1.0_dp/r
1115
1116 ! Compute the radial function
1117 ! code for point multipole without screening
1118 f(0) = ir
1119 DO i = 1, 5
1120 f(i) = irab2*f(i - 1)
1121 END DO
1122
1123
1124 ! Compute the Tensor components
1125 force_eval = do_stress
1126 IF (task(1, 1)) THEN
1127 tij = f(0)*fac_ij
1128 force_eval = do_forces .OR. do_efield1
1129 END IF
1130 IF (task(2, 2)) force_eval = force_eval .OR. do_efield0
1131 IF (task(1, 2) .OR. force_eval) THEN
1132 force_eval = do_stress
1133 tij_a = -rab*f(1)*fac_ij
1134 IF (task(1, 2)) force_eval = force_eval .OR. do_forces
1135 END IF
1136 IF (task(1, 1)) force_eval = force_eval .OR. do_efield2
1137 IF (task(3, 3)) force_eval = force_eval .OR. do_efield0
1138 IF (task(2, 2) .OR. task(3, 1) .OR. force_eval) THEN
1139 force_eval = do_stress
1140 DO b = 1, 3
1141 DO a = 1, 3
1142 tmp = rab(a)*rab(b)*fac_ij
1143 tij_ab(a, b) = 3.0_dp*tmp*f(2)
1144 IF (a == b) tij_ab(a, b) = tij_ab(a, b) - f(1)*fac_ij
1145 END DO
1146 END DO
1147 IF (task(2, 2) .OR. task(3, 1)) force_eval = force_eval .OR. do_forces
1148 END IF
1149 IF (task(2, 2)) force_eval = force_eval .OR. do_efield2
1150 IF (task(3, 3)) force_eval = force_eval .OR. do_efield1
1151 IF (task(3, 2) .OR. force_eval) THEN
1152 force_eval = do_stress
1153 DO c = 1, 3
1154 DO b = 1, 3
1155 DO a = 1, 3
1156 tmp = rab(a)*rab(b)*rab(c)*fac_ij
1157 tij_abc(a, b, c) = -15.0_dp*tmp*f(3)
1158 tmp = 3.0_dp*f(2)*fac_ij
1159 IF (a == b) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(c)
1160 IF (a == c) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(b)
1161 IF (b == c) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(a)
1162 END DO
1163 END DO
1164 END DO
1165 IF (task(3, 2)) force_eval = force_eval .OR. do_forces
1166 END IF
1167 IF (task(3, 3) .OR. force_eval) THEN
1168 force_eval = do_stress
1169 DO d = 1, 3
1170 DO c = 1, 3
1171 DO b = 1, 3
1172 DO a = 1, 3
1173 tmp = rab(a)*rab(b)*rab(c)*rab(d)*fac_ij
1174 tij_abcd(a, b, c, d) = 105.0_dp*tmp*f(4)
1175 tmp1 = 15.0_dp*f(3)*fac_ij
1176 tmp2 = 3.0_dp*f(2)*fac_ij
1177 IF (a == b) THEN
1178 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(c)*rab(d)
1179 IF (c == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
1180 END IF
1181 IF (a == c) THEN
1182 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(b)*rab(d)
1183 IF (b == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
1184 END IF
1185 IF (a == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(b)*rab(c)
1186 IF (b == c) THEN
1187 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(d)
1188 IF (a == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
1189 END IF
1190 IF (b == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(c)
1191 IF (c == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(b)
1192 END DO
1193 END DO
1194 END DO
1195 END DO
1196 IF (task(3, 3)) force_eval = force_eval .OR. do_forces
1197 END IF
1198 IF (force_eval) THEN
1199 force_eval = do_stress
1200 DO e = 1, 3
1201 DO d = 1, 3
1202 DO c = 1, 3
1203 DO b = 1, 3
1204 DO a = 1, 3
1205 tmp = rab(a)*rab(b)*rab(c)*rab(d)*rab(e)*fac_ij
1206 tij_abcde(a, b, c, d, e) = -945.0_dp*tmp*f(5)
1207 tmp1 = 105.0_dp*f(4)*fac_ij
1208 tmp2 = 15.0_dp*f(3)*fac_ij
1209 IF (a == b) THEN
1210 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(c)*rab(d)*rab(e)
1211 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
1212 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
1213 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
1214 END IF
1215 IF (a == c) THEN
1216 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(d)*rab(e)
1217 IF (b == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
1218 IF (b == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
1219 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
1220 END IF
1221 IF (a == d) THEN
1222 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(c)*rab(e)
1223 IF (b == c) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
1224 IF (b == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
1225 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
1226 END IF
1227 IF (a == e) THEN
1228 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(c)*rab(d)
1229 IF (b == c) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
1230 IF (b == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
1231 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
1232 END IF
1233 IF (b == c) THEN
1234 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(d)*rab(e)
1235 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
1236 END IF
1237 IF (b == d) THEN
1238 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(c)*rab(e)
1239 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
1240 END IF
1241 IF (b == e) THEN
1242 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(c)*rab(d)
1243 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
1244 END IF
1245 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(e)
1246 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(d)
1247 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(c)
1248 END DO
1249 END DO
1250 END DO
1251 END DO
1252 END DO
1253 END IF
1254 eloc = 0.0_dp
1255 fr = 0.0_dp
1256 ef0_i = 0.0_dp
1257 ef0_j = 0.0_dp
1258 ef1_j = 0.0_dp
1259 ef1_i = 0.0_dp
1260 ef2_j = 0.0_dp
1261 ef2_i = 0.0_dp
1262
1263
1264 ! Initialize the charge, dipole and quadrupole for atom A and B
1265 IF (debug_this_module) THEN
1266 ch_j = huge(0.0_dp)
1267 ch_i = huge(0.0_dp)
1268 dp_j = huge(0.0_dp)
1269 dp_i = huge(0.0_dp)
1270 qp_j = huge(0.0_dp)
1271 qp_i = huge(0.0_dp)
1272 END IF
1273 IF (any(task(1, :))) THEN
1274 ch_j = charges(atom_a)
1275 ch_i = charges(atom_b)
1276 END IF
1277 IF (any(task(2, :))) THEN
1278 dp_j = dipoles(:, atom_a)
1279 dp_i = dipoles(:, atom_b)
1280 END IF
1281 IF (any(task(3, :))) THEN
1282 qp_j = quadrupoles(:, :, atom_a)
1283 qp_i = quadrupoles(:, :, atom_b)
1284 END IF
1285 IF (task(1, 1)) THEN
1286 ! Charge - Charge
1287 eloc = eloc + ch_i*tij*ch_j
1288 ! Forces on particle i (locally b)
1289 IF (do_forces .OR. do_stress) THEN
1290 fr(1) = fr(1) - ch_j*tij_a(1)*ch_i
1291 fr(2) = fr(2) - ch_j*tij_a(2)*ch_i
1292 fr(3) = fr(3) - ch_j*tij_a(3)*ch_i
1293 END IF
1294 ! Electric fields
1295 IF (do_efield) THEN
1296 ! Potential
1297 IF (do_efield0) THEN
1298 ef0_i = ef0_i + tij*ch_j
1299
1300 ef0_j = ef0_j + tij*ch_i
1301 END IF
1302 ! Electric field
1303 IF (do_efield1) THEN
1304 ef1_i(1) = ef1_i(1) - tij_a(1)*ch_j
1305 ef1_i(2) = ef1_i(2) - tij_a(2)*ch_j
1306 ef1_i(3) = ef1_i(3) - tij_a(3)*ch_j
1307
1308 ef1_j(1) = ef1_j(1) + tij_a(1)*ch_i
1309 ef1_j(2) = ef1_j(2) + tij_a(2)*ch_i
1310 ef1_j(3) = ef1_j(3) + tij_a(3)*ch_i
1311
1312
1313 END IF
1314 ! Electric field gradient
1315 IF (do_efield2) THEN
1316 ef2_i(1, 1) = ef2_i(1, 1) - tij_ab(1, 1)*ch_j
1317 ef2_i(2, 1) = ef2_i(2, 1) - tij_ab(2, 1)*ch_j
1318 ef2_i(3, 1) = ef2_i(3, 1) - tij_ab(3, 1)*ch_j
1319 ef2_i(1, 2) = ef2_i(1, 2) - tij_ab(1, 2)*ch_j
1320 ef2_i(2, 2) = ef2_i(2, 2) - tij_ab(2, 2)*ch_j
1321 ef2_i(3, 2) = ef2_i(3, 2) - tij_ab(3, 2)*ch_j
1322 ef2_i(1, 3) = ef2_i(1, 3) - tij_ab(1, 3)*ch_j
1323 ef2_i(2, 3) = ef2_i(2, 3) - tij_ab(2, 3)*ch_j
1324 ef2_i(3, 3) = ef2_i(3, 3) - tij_ab(3, 3)*ch_j
1325
1326 ef2_j(1, 1) = ef2_j(1, 1) - tij_ab(1, 1)*ch_i
1327 ef2_j(2, 1) = ef2_j(2, 1) - tij_ab(2, 1)*ch_i
1328 ef2_j(3, 1) = ef2_j(3, 1) - tij_ab(3, 1)*ch_i
1329 ef2_j(1, 2) = ef2_j(1, 2) - tij_ab(1, 2)*ch_i
1330 ef2_j(2, 2) = ef2_j(2, 2) - tij_ab(2, 2)*ch_i
1331 ef2_j(3, 2) = ef2_j(3, 2) - tij_ab(3, 2)*ch_i
1332 ef2_j(1, 3) = ef2_j(1, 3) - tij_ab(1, 3)*ch_i
1333 ef2_j(2, 3) = ef2_j(2, 3) - tij_ab(2, 3)*ch_i
1334 ef2_j(3, 3) = ef2_j(3, 3) - tij_ab(3, 3)*ch_i
1335 END IF
1336 END IF
1337 END IF
1338 IF (task(2, 2)) THEN
1339 ! Dipole - Dipole
1340 tmp = -(dp_i(1)*(tij_ab(1, 1)*dp_j(1) + &
1341 tij_ab(2, 1)*dp_j(2) + &
1342 tij_ab(3, 1)*dp_j(3)) + &
1343 dp_i(2)*(tij_ab(1, 2)*dp_j(1) + &
1344 tij_ab(2, 2)*dp_j(2) + &
1345 tij_ab(3, 2)*dp_j(3)) + &
1346 dp_i(3)*(tij_ab(1, 3)*dp_j(1) + &
1347 tij_ab(2, 3)*dp_j(2) + &
1348 tij_ab(3, 3)*dp_j(3)))
1349 eloc = eloc + tmp
1350 ! Forces on particle i (locally b)
1351 IF (do_forces .OR. do_stress) THEN
1352 DO k = 1, 3
1353 fr(k) = fr(k) + dp_i(1)*(tij_abc(1, 1, k)*dp_j(1) + &
1354 tij_abc(2, 1, k)*dp_j(2) + &
1355 tij_abc(3, 1, k)*dp_j(3)) &
1356 + dp_i(2)*(tij_abc(1, 2, k)*dp_j(1) + &
1357 tij_abc(2, 2, k)*dp_j(2) + &
1358 tij_abc(3, 2, k)*dp_j(3)) &
1359 + dp_i(3)*(tij_abc(1, 3, k)*dp_j(1) + &
1360 tij_abc(2, 3, k)*dp_j(2) + &
1361 tij_abc(3, 3, k)*dp_j(3))
1362 END DO
1363 END IF
1364 ! Electric fields
1365 IF (do_efield) THEN
1366 ! Potential
1367 IF (do_efield0) THEN
1368 ef0_i = ef0_i - (tij_a(1)*dp_j(1) + &
1369 tij_a(2)*dp_j(2) + &
1370 tij_a(3)*dp_j(3))
1371
1372 ef0_j = ef0_j + (tij_a(1)*dp_i(1) + &
1373 tij_a(2)*dp_i(2) + &
1374 tij_a(3)*dp_i(3))
1375 END IF
1376 ! Electric field
1377 IF (do_efield1) THEN
1378 ef1_i(1) = ef1_i(1) + (tij_ab(1, 1)*dp_j(1) + &
1379 tij_ab(2, 1)*dp_j(2) + &
1380 tij_ab(3, 1)*dp_j(3))
1381 ef1_i(2) = ef1_i(2) + (tij_ab(1, 2)*dp_j(1) + &
1382 tij_ab(2, 2)*dp_j(2) + &
1383 tij_ab(3, 2)*dp_j(3))
1384 ef1_i(3) = ef1_i(3) + (tij_ab(1, 3)*dp_j(1) + &
1385 tij_ab(2, 3)*dp_j(2) + &
1386 tij_ab(3, 3)*dp_j(3))
1387
1388 ef1_j(1) = ef1_j(1) + (tij_ab(1, 1)*dp_i(1) + &
1389 tij_ab(2, 1)*dp_i(2) + &
1390 tij_ab(3, 1)*dp_i(3))
1391 ef1_j(2) = ef1_j(2) + (tij_ab(1, 2)*dp_i(1) + &
1392 tij_ab(2, 2)*dp_i(2) + &
1393 tij_ab(3, 2)*dp_i(3))
1394 ef1_j(3) = ef1_j(3) + (tij_ab(1, 3)*dp_i(1) + &
1395 tij_ab(2, 3)*dp_i(2) + &
1396 tij_ab(3, 3)*dp_i(3))
1397 END IF
1398 ! Electric field gradient
1399 IF (do_efield2) THEN
1400 ef2_i(1, 1) = ef2_i(1, 1) + (tij_abc(1, 1, 1)*dp_j(1) + &
1401 tij_abc(2, 1, 1)*dp_j(2) + &
1402 tij_abc(3, 1, 1)*dp_j(3))
1403 ef2_i(1, 2) = ef2_i(1, 2) + (tij_abc(1, 1, 2)*dp_j(1) + &
1404 tij_abc(2, 1, 2)*dp_j(2) + &
1405 tij_abc(3, 1, 2)*dp_j(3))
1406 ef2_i(1, 3) = ef2_i(1, 3) + (tij_abc(1, 1, 3)*dp_j(1) + &
1407 tij_abc(2, 1, 3)*dp_j(2) + &
1408 tij_abc(3, 1, 3)*dp_j(3))
1409 ef2_i(2, 1) = ef2_i(2, 1) + (tij_abc(1, 2, 1)*dp_j(1) + &
1410 tij_abc(2, 2, 1)*dp_j(2) + &
1411 tij_abc(3, 2, 1)*dp_j(3))
1412 ef2_i(2, 2) = ef2_i(2, 2) + (tij_abc(1, 2, 2)*dp_j(1) + &
1413 tij_abc(2, 2, 2)*dp_j(2) + &
1414 tij_abc(3, 2, 2)*dp_j(3))
1415 ef2_i(2, 3) = ef2_i(2, 3) + (tij_abc(1, 2, 3)*dp_j(1) + &
1416 tij_abc(2, 2, 3)*dp_j(2) + &
1417 tij_abc(3, 2, 3)*dp_j(3))
1418 ef2_i(3, 1) = ef2_i(3, 1) + (tij_abc(1, 3, 1)*dp_j(1) + &
1419 tij_abc(2, 3, 1)*dp_j(2) + &
1420 tij_abc(3, 3, 1)*dp_j(3))
1421 ef2_i(3, 2) = ef2_i(3, 2) + (tij_abc(1, 3, 2)*dp_j(1) + &
1422 tij_abc(2, 3, 2)*dp_j(2) + &
1423 tij_abc(3, 3, 2)*dp_j(3))
1424 ef2_i(3, 3) = ef2_i(3, 3) + (tij_abc(1, 3, 3)*dp_j(1) + &
1425 tij_abc(2, 3, 3)*dp_j(2) + &
1426 tij_abc(3, 3, 3)*dp_j(3))
1427
1428 ef2_j(1, 1) = ef2_j(1, 1) - (tij_abc(1, 1, 1)*dp_i(1) + &
1429 tij_abc(2, 1, 1)*dp_i(2) + &
1430 tij_abc(3, 1, 1)*dp_i(3))
1431 ef2_j(1, 2) = ef2_j(1, 2) - (tij_abc(1, 1, 2)*dp_i(1) + &
1432 tij_abc(2, 1, 2)*dp_i(2) + &
1433 tij_abc(3, 1, 2)*dp_i(3))
1434 ef2_j(1, 3) = ef2_j(1, 3) - (tij_abc(1, 1, 3)*dp_i(1) + &
1435 tij_abc(2, 1, 3)*dp_i(2) + &
1436 tij_abc(3, 1, 3)*dp_i(3))
1437 ef2_j(2, 1) = ef2_j(2, 1) - (tij_abc(1, 2, 1)*dp_i(1) + &
1438 tij_abc(2, 2, 1)*dp_i(2) + &
1439 tij_abc(3, 2, 1)*dp_i(3))
1440 ef2_j(2, 2) = ef2_j(2, 2) - (tij_abc(1, 2, 2)*dp_i(1) + &
1441 tij_abc(2, 2, 2)*dp_i(2) + &
1442 tij_abc(3, 2, 2)*dp_i(3))
1443 ef2_j(2, 3) = ef2_j(2, 3) - (tij_abc(1, 2, 3)*dp_i(1) + &
1444 tij_abc(2, 2, 3)*dp_i(2) + &
1445 tij_abc(3, 2, 3)*dp_i(3))
1446 ef2_j(3, 1) = ef2_j(3, 1) - (tij_abc(1, 3, 1)*dp_i(1) + &
1447 tij_abc(2, 3, 1)*dp_i(2) + &
1448 tij_abc(3, 3, 1)*dp_i(3))
1449 ef2_j(3, 2) = ef2_j(3, 2) - (tij_abc(1, 3, 2)*dp_i(1) + &
1450 tij_abc(2, 3, 2)*dp_i(2) + &
1451 tij_abc(3, 3, 2)*dp_i(3))
1452 ef2_j(3, 3) = ef2_j(3, 3) - (tij_abc(1, 3, 3)*dp_i(1) + &
1453 tij_abc(2, 3, 3)*dp_i(2) + &
1454 tij_abc(3, 3, 3)*dp_i(3))
1455 END IF
1456 END IF
1457 END IF
1458 IF (task(2, 1)) THEN
1459 ! Dipole - Charge
1460 tmp = ch_j*(tij_a(1)*dp_i(1) + &
1461 tij_a(2)*dp_i(2) + &
1462 tij_a(3)*dp_i(3)) &
1463 - ch_i*(tij_a(1)*dp_j(1) + &
1464 tij_a(2)*dp_j(2) + &
1465 tij_a(3)*dp_j(3))
1466 eloc = eloc + tmp
1467 ! Forces on particle i (locally b)
1468 IF (do_forces .OR. do_stress) THEN
1469 DO k = 1, 3
1470 fr(k) = fr(k) - ch_j*(tij_ab(1, k)*dp_i(1) + &
1471 tij_ab(2, k)*dp_i(2) + &
1472 tij_ab(3, k)*dp_i(3)) &
1473 + ch_i*(tij_ab(1, k)*dp_j(1) + &
1474 tij_ab(2, k)*dp_j(2) + &
1475 tij_ab(3, k)*dp_j(3))
1476 END DO
1477 END IF
1478 END IF
1479 IF (task(3, 3)) THEN
1480 ! Quadrupole - Quadrupole
1481 fac = 1.0_dp/9.0_dp
1482 tmp11 = qp_i(1, 1)*(tij_abcd(1, 1, 1, 1)*qp_j(1, 1) + &
1483 tij_abcd(2, 1, 1, 1)*qp_j(2, 1) + &
1484 tij_abcd(3, 1, 1, 1)*qp_j(3, 1) + &
1485 tij_abcd(1, 2, 1, 1)*qp_j(1, 2) + &
1486 tij_abcd(2, 2, 1, 1)*qp_j(2, 2) + &
1487 tij_abcd(3, 2, 1, 1)*qp_j(3, 2) + &
1488 tij_abcd(1, 3, 1, 1)*qp_j(1, 3) + &
1489 tij_abcd(2, 3, 1, 1)*qp_j(2, 3) + &
1490 tij_abcd(3, 3, 1, 1)*qp_j(3, 3))
1491 tmp21 = qp_i(2, 1)*(tij_abcd(1, 1, 1, 2)*qp_j(1, 1) + &
1492 tij_abcd(2, 1, 1, 2)*qp_j(2, 1) + &
1493 tij_abcd(3, 1, 1, 2)*qp_j(3, 1) + &
1494 tij_abcd(1, 2, 1, 2)*qp_j(1, 2) + &
1495 tij_abcd(2, 2, 1, 2)*qp_j(2, 2) + &
1496 tij_abcd(3, 2, 1, 2)*qp_j(3, 2) + &
1497 tij_abcd(1, 3, 1, 2)*qp_j(1, 3) + &
1498 tij_abcd(2, 3, 1, 2)*qp_j(2, 3) + &
1499 tij_abcd(3, 3, 1, 2)*qp_j(3, 3))
1500 tmp31 = qp_i(3, 1)*(tij_abcd(1, 1, 1, 3)*qp_j(1, 1) + &
1501 tij_abcd(2, 1, 1, 3)*qp_j(2, 1) + &
1502 tij_abcd(3, 1, 1, 3)*qp_j(3, 1) + &
1503 tij_abcd(1, 2, 1, 3)*qp_j(1, 2) + &
1504 tij_abcd(2, 2, 1, 3)*qp_j(2, 2) + &
1505 tij_abcd(3, 2, 1, 3)*qp_j(3, 2) + &
1506 tij_abcd(1, 3, 1, 3)*qp_j(1, 3) + &
1507 tij_abcd(2, 3, 1, 3)*qp_j(2, 3) + &
1508 tij_abcd(3, 3, 1, 3)*qp_j(3, 3))
1509 tmp22 = qp_i(2, 2)*(tij_abcd(1, 1, 2, 2)*qp_j(1, 1) + &
1510 tij_abcd(2, 1, 2, 2)*qp_j(2, 1) + &
1511 tij_abcd(3, 1, 2, 2)*qp_j(3, 1) + &
1512 tij_abcd(1, 2, 2, 2)*qp_j(1, 2) + &
1513 tij_abcd(2, 2, 2, 2)*qp_j(2, 2) + &
1514 tij_abcd(3, 2, 2, 2)*qp_j(3, 2) + &
1515 tij_abcd(1, 3, 2, 2)*qp_j(1, 3) + &
1516 tij_abcd(2, 3, 2, 2)*qp_j(2, 3) + &
1517 tij_abcd(3, 3, 2, 2)*qp_j(3, 3))
1518 tmp32 = qp_i(3, 2)*(tij_abcd(1, 1, 2, 3)*qp_j(1, 1) + &
1519 tij_abcd(2, 1, 2, 3)*qp_j(2, 1) + &
1520 tij_abcd(3, 1, 2, 3)*qp_j(3, 1) + &
1521 tij_abcd(1, 2, 2, 3)*qp_j(1, 2) + &
1522 tij_abcd(2, 2, 2, 3)*qp_j(2, 2) + &
1523 tij_abcd(3, 2, 2, 3)*qp_j(3, 2) + &
1524 tij_abcd(1, 3, 2, 3)*qp_j(1, 3) + &
1525 tij_abcd(2, 3, 2, 3)*qp_j(2, 3) + &
1526 tij_abcd(3, 3, 2, 3)*qp_j(3, 3))
1527 tmp33 = qp_i(3, 3)*(tij_abcd(1, 1, 3, 3)*qp_j(1, 1) + &
1528 tij_abcd(2, 1, 3, 3)*qp_j(2, 1) + &
1529 tij_abcd(3, 1, 3, 3)*qp_j(3, 1) + &
1530 tij_abcd(1, 2, 3, 3)*qp_j(1, 2) + &
1531 tij_abcd(2, 2, 3, 3)*qp_j(2, 2) + &
1532 tij_abcd(3, 2, 3, 3)*qp_j(3, 2) + &
1533 tij_abcd(1, 3, 3, 3)*qp_j(1, 3) + &
1534 tij_abcd(2, 3, 3, 3)*qp_j(2, 3) + &
1535 tij_abcd(3, 3, 3, 3)*qp_j(3, 3))
1536 tmp12 = tmp21
1537 tmp13 = tmp31
1538 tmp23 = tmp32
1539 tmp = tmp11 + tmp12 + tmp13 + &
1540 tmp21 + tmp22 + tmp23 + &
1541 tmp31 + tmp32 + tmp33
1542
1543 eloc = eloc + fac*tmp
1544 ! Forces on particle i (locally b)
1545 IF (do_forces .OR. do_stress) THEN
1546 DO k = 1, 3
1547 tmp11 = qp_i(1, 1)*(tij_abcde(1, 1, 1, 1, k)*qp_j(1, 1) + &
1548 tij_abcde(2, 1, 1, 1, k)*qp_j(2, 1) + &
1549 tij_abcde(3, 1, 1, 1, k)*qp_j(3, 1) + &
1550 tij_abcde(1, 2, 1, 1, k)*qp_j(1, 2) + &
1551 tij_abcde(2, 2, 1, 1, k)*qp_j(2, 2) + &
1552 tij_abcde(3, 2, 1, 1, k)*qp_j(3, 2) + &
1553 tij_abcde(1, 3, 1, 1, k)*qp_j(1, 3) + &
1554 tij_abcde(2, 3, 1, 1, k)*qp_j(2, 3) + &
1555 tij_abcde(3, 3, 1, 1, k)*qp_j(3, 3))
1556 tmp21 = qp_i(2, 1)*(tij_abcde(1, 1, 2, 1, k)*qp_j(1, 1) + &
1557 tij_abcde(2, 1, 2, 1, k)*qp_j(2, 1) + &
1558 tij_abcde(3, 1, 2, 1, k)*qp_j(3, 1) + &
1559 tij_abcde(1, 2, 2, 1, k)*qp_j(1, 2) + &
1560 tij_abcde(2, 2, 2, 1, k)*qp_j(2, 2) + &
1561 tij_abcde(3, 2, 2, 1, k)*qp_j(3, 2) + &
1562 tij_abcde(1, 3, 2, 1, k)*qp_j(1, 3) + &
1563 tij_abcde(2, 3, 2, 1, k)*qp_j(2, 3) + &
1564 tij_abcde(3, 3, 2, 1, k)*qp_j(3, 3))
1565 tmp31 = qp_i(3, 1)*(tij_abcde(1, 1, 3, 1, k)*qp_j(1, 1) + &
1566 tij_abcde(2, 1, 3, 1, k)*qp_j(2, 1) + &
1567 tij_abcde(3, 1, 3, 1, k)*qp_j(3, 1) + &
1568 tij_abcde(1, 2, 3, 1, k)*qp_j(1, 2) + &
1569 tij_abcde(2, 2, 3, 1, k)*qp_j(2, 2) + &
1570 tij_abcde(3, 2, 3, 1, k)*qp_j(3, 2) + &
1571 tij_abcde(1, 3, 3, 1, k)*qp_j(1, 3) + &
1572 tij_abcde(2, 3, 3, 1, k)*qp_j(2, 3) + &
1573 tij_abcde(3, 3, 3, 1, k)*qp_j(3, 3))
1574 tmp22 = qp_i(2, 2)*(tij_abcde(1, 1, 2, 2, k)*qp_j(1, 1) + &
1575 tij_abcde(2, 1, 2, 2, k)*qp_j(2, 1) + &
1576 tij_abcde(3, 1, 2, 2, k)*qp_j(3, 1) + &
1577 tij_abcde(1, 2, 2, 2, k)*qp_j(1, 2) + &
1578 tij_abcde(2, 2, 2, 2, k)*qp_j(2, 2) + &
1579 tij_abcde(3, 2, 2, 2, k)*qp_j(3, 2) + &
1580 tij_abcde(1, 3, 2, 2, k)*qp_j(1, 3) + &
1581 tij_abcde(2, 3, 2, 2, k)*qp_j(2, 3) + &
1582 tij_abcde(3, 3, 2, 2, k)*qp_j(3, 3))
1583 tmp32 = qp_i(3, 2)*(tij_abcde(1, 1, 3, 2, k)*qp_j(1, 1) + &
1584 tij_abcde(2, 1, 3, 2, k)*qp_j(2, 1) + &
1585 tij_abcde(3, 1, 3, 2, k)*qp_j(3, 1) + &
1586 tij_abcde(1, 2, 3, 2, k)*qp_j(1, 2) + &
1587 tij_abcde(2, 2, 3, 2, k)*qp_j(2, 2) + &
1588 tij_abcde(3, 2, 3, 2, k)*qp_j(3, 2) + &
1589 tij_abcde(1, 3, 3, 2, k)*qp_j(1, 3) + &
1590 tij_abcde(2, 3, 3, 2, k)*qp_j(2, 3) + &
1591 tij_abcde(3, 3, 3, 2, k)*qp_j(3, 3))
1592 tmp33 = qp_i(3, 3)*(tij_abcde(1, 1, 3, 3, k)*qp_j(1, 1) + &
1593 tij_abcde(2, 1, 3, 3, k)*qp_j(2, 1) + &
1594 tij_abcde(3, 1, 3, 3, k)*qp_j(3, 1) + &
1595 tij_abcde(1, 2, 3, 3, k)*qp_j(1, 2) + &
1596 tij_abcde(2, 2, 3, 3, k)*qp_j(2, 2) + &
1597 tij_abcde(3, 2, 3, 3, k)*qp_j(3, 2) + &
1598 tij_abcde(1, 3, 3, 3, k)*qp_j(1, 3) + &
1599 tij_abcde(2, 3, 3, 3, k)*qp_j(2, 3) + &
1600 tij_abcde(3, 3, 3, 3, k)*qp_j(3, 3))
1601 tmp12 = tmp21
1602 tmp13 = tmp31
1603 tmp23 = tmp32
1604 fr(k) = fr(k) - fac*(tmp11 + tmp12 + tmp13 + &
1605 tmp21 + tmp22 + tmp23 + &
1606 tmp31 + tmp32 + tmp33)
1607 END DO
1608 END IF
1609 ! Electric field
1610 IF (do_efield) THEN
1611 fac = 1.0_dp/3.0_dp
1612 ! Potential
1613 IF (do_efield0) THEN
1614 ef0_i = ef0_i + fac*(tij_ab(1, 1)*qp_j(1, 1) + &
1615 tij_ab(2, 1)*qp_j(2, 1) + &
1616 tij_ab(3, 1)*qp_j(3, 1) + &
1617 tij_ab(1, 2)*qp_j(1, 2) + &
1618 tij_ab(2, 2)*qp_j(2, 2) + &
1619 tij_ab(3, 2)*qp_j(3, 2) + &
1620 tij_ab(1, 3)*qp_j(1, 3) + &
1621 tij_ab(2, 3)*qp_j(2, 3) + &
1622 tij_ab(3, 3)*qp_j(3, 3))
1623
1624 ef0_j = ef0_j + fac*(tij_ab(1, 1)*qp_i(1, 1) + &
1625 tij_ab(2, 1)*qp_i(2, 1) + &
1626 tij_ab(3, 1)*qp_i(3, 1) + &
1627 tij_ab(1, 2)*qp_i(1, 2) + &
1628 tij_ab(2, 2)*qp_i(2, 2) + &
1629 tij_ab(3, 2)*qp_i(3, 2) + &
1630 tij_ab(1, 3)*qp_i(1, 3) + &
1631 tij_ab(2, 3)*qp_i(2, 3) + &
1632 tij_ab(3, 3)*qp_i(3, 3))
1633 END IF
1634 ! Electric field
1635 IF (do_efield1) THEN
1636 ef1_i(1) = ef1_i(1) - fac*(tij_abc(1, 1, 1)*qp_j(1, 1) + &
1637 tij_abc(2, 1, 1)*qp_j(2, 1) + &
1638 tij_abc(3, 1, 1)*qp_j(3, 1) + &
1639 tij_abc(1, 2, 1)*qp_j(1, 2) + &
1640 tij_abc(2, 2, 1)*qp_j(2, 2) + &
1641 tij_abc(3, 2, 1)*qp_j(3, 2) + &
1642 tij_abc(1, 3, 1)*qp_j(1, 3) + &
1643 tij_abc(2, 3, 1)*qp_j(2, 3) + &
1644 tij_abc(3, 3, 1)*qp_j(3, 3))
1645 ef1_i(2) = ef1_i(2) - fac*(tij_abc(1, 1, 2)*qp_j(1, 1) + &
1646 tij_abc(2, 1, 2)*qp_j(2, 1) + &
1647 tij_abc(3, 1, 2)*qp_j(3, 1) + &
1648 tij_abc(1, 2, 2)*qp_j(1, 2) + &
1649 tij_abc(2, 2, 2)*qp_j(2, 2) + &
1650 tij_abc(3, 2, 2)*qp_j(3, 2) + &
1651 tij_abc(1, 3, 2)*qp_j(1, 3) + &
1652 tij_abc(2, 3, 2)*qp_j(2, 3) + &
1653 tij_abc(3, 3, 2)*qp_j(3, 3))
1654 ef1_i(3) = ef1_i(3) - fac*(tij_abc(1, 1, 3)*qp_j(1, 1) + &
1655 tij_abc(2, 1, 3)*qp_j(2, 1) + &
1656 tij_abc(3, 1, 3)*qp_j(3, 1) + &
1657 tij_abc(1, 2, 3)*qp_j(1, 2) + &
1658 tij_abc(2, 2, 3)*qp_j(2, 2) + &
1659 tij_abc(3, 2, 3)*qp_j(3, 2) + &
1660 tij_abc(1, 3, 3)*qp_j(1, 3) + &
1661 tij_abc(2, 3, 3)*qp_j(2, 3) + &
1662 tij_abc(3, 3, 3)*qp_j(3, 3))
1663
1664 ef1_j(1) = ef1_j(1) + fac*(tij_abc(1, 1, 1)*qp_i(1, 1) + &
1665 tij_abc(2, 1, 1)*qp_i(2, 1) + &
1666 tij_abc(3, 1, 1)*qp_i(3, 1) + &
1667 tij_abc(1, 2, 1)*qp_i(1, 2) + &
1668 tij_abc(2, 2, 1)*qp_i(2, 2) + &
1669 tij_abc(3, 2, 1)*qp_i(3, 2) + &
1670 tij_abc(1, 3, 1)*qp_i(1, 3) + &
1671 tij_abc(2, 3, 1)*qp_i(2, 3) + &
1672 tij_abc(3, 3, 1)*qp_i(3, 3))
1673 ef1_j(2) = ef1_j(2) + fac*(tij_abc(1, 1, 2)*qp_i(1, 1) + &
1674 tij_abc(2, 1, 2)*qp_i(2, 1) + &
1675 tij_abc(3, 1, 2)*qp_i(3, 1) + &
1676 tij_abc(1, 2, 2)*qp_i(1, 2) + &
1677 tij_abc(2, 2, 2)*qp_i(2, 2) + &
1678 tij_abc(3, 2, 2)*qp_i(3, 2) + &
1679 tij_abc(1, 3, 2)*qp_i(1, 3) + &
1680 tij_abc(2, 3, 2)*qp_i(2, 3) + &
1681 tij_abc(3, 3, 2)*qp_i(3, 3))
1682 ef1_j(3) = ef1_j(3) + fac*(tij_abc(1, 1, 3)*qp_i(1, 1) + &
1683 tij_abc(2, 1, 3)*qp_i(2, 1) + &
1684 tij_abc(3, 1, 3)*qp_i(3, 1) + &
1685 tij_abc(1, 2, 3)*qp_i(1, 2) + &
1686 tij_abc(2, 2, 3)*qp_i(2, 2) + &
1687 tij_abc(3, 2, 3)*qp_i(3, 2) + &
1688 tij_abc(1, 3, 3)*qp_i(1, 3) + &
1689 tij_abc(2, 3, 3)*qp_i(2, 3) + &
1690 tij_abc(3, 3, 3)*qp_i(3, 3))
1691 END IF
1692 ! Electric field gradient
1693 IF (do_efield2) THEN
1694 tmp11 = fac*(tij_abcd(1, 1, 1, 1)*qp_j(1, 1) + &
1695 tij_abcd(2, 1, 1, 1)*qp_j(2, 1) + &
1696 tij_abcd(3, 1, 1, 1)*qp_j(3, 1) + &
1697 tij_abcd(1, 2, 1, 1)*qp_j(1, 2) + &
1698 tij_abcd(2, 2, 1, 1)*qp_j(2, 2) + &
1699 tij_abcd(3, 2, 1, 1)*qp_j(3, 2) + &
1700 tij_abcd(1, 3, 1, 1)*qp_j(1, 3) + &
1701 tij_abcd(2, 3, 1, 1)*qp_j(2, 3) + &
1702 tij_abcd(3, 3, 1, 1)*qp_j(3, 3))
1703 tmp12 = fac*(tij_abcd(1, 1, 1, 2)*qp_j(1, 1) + &
1704 tij_abcd(2, 1, 1, 2)*qp_j(2, 1) + &
1705 tij_abcd(3, 1, 1, 2)*qp_j(3, 1) + &
1706 tij_abcd(1, 2, 1, 2)*qp_j(1, 2) + &
1707 tij_abcd(2, 2, 1, 2)*qp_j(2, 2) + &
1708 tij_abcd(3, 2, 1, 2)*qp_j(3, 2) + &
1709 tij_abcd(1, 3, 1, 2)*qp_j(1, 3) + &
1710 tij_abcd(2, 3, 1, 2)*qp_j(2, 3) + &
1711 tij_abcd(3, 3, 1, 2)*qp_j(3, 3))
1712 tmp13 = fac*(tij_abcd(1, 1, 1, 3)*qp_j(1, 1) + &
1713 tij_abcd(2, 1, 1, 3)*qp_j(2, 1) + &
1714 tij_abcd(3, 1, 1, 3)*qp_j(3, 1) + &
1715 tij_abcd(1, 2, 1, 3)*qp_j(1, 2) + &
1716 tij_abcd(2, 2, 1, 3)*qp_j(2, 2) + &
1717 tij_abcd(3, 2, 1, 3)*qp_j(3, 2) + &
1718 tij_abcd(1, 3, 1, 3)*qp_j(1, 3) + &
1719 tij_abcd(2, 3, 1, 3)*qp_j(2, 3) + &
1720 tij_abcd(3, 3, 1, 3)*qp_j(3, 3))
1721 tmp22 = fac*(tij_abcd(1, 1, 2, 2)*qp_j(1, 1) + &
1722 tij_abcd(2, 1, 2, 2)*qp_j(2, 1) + &
1723 tij_abcd(3, 1, 2, 2)*qp_j(3, 1) + &
1724 tij_abcd(1, 2, 2, 2)*qp_j(1, 2) + &
1725 tij_abcd(2, 2, 2, 2)*qp_j(2, 2) + &
1726 tij_abcd(3, 2, 2, 2)*qp_j(3, 2) + &
1727 tij_abcd(1, 3, 2, 2)*qp_j(1, 3) + &
1728 tij_abcd(2, 3, 2, 2)*qp_j(2, 3) + &
1729 tij_abcd(3, 3, 2, 2)*qp_j(3, 3))
1730 tmp23 = fac*(tij_abcd(1, 1, 2, 3)*qp_j(1, 1) + &
1731 tij_abcd(2, 1, 2, 3)*qp_j(2, 1) + &
1732 tij_abcd(3, 1, 2, 3)*qp_j(3, 1) + &
1733 tij_abcd(1, 2, 2, 3)*qp_j(1, 2) + &
1734 tij_abcd(2, 2, 2, 3)*qp_j(2, 2) + &
1735 tij_abcd(3, 2, 2, 3)*qp_j(3, 2) + &
1736 tij_abcd(1, 3, 2, 3)*qp_j(1, 3) + &
1737 tij_abcd(2, 3, 2, 3)*qp_j(2, 3) + &
1738 tij_abcd(3, 3, 2, 3)*qp_j(3, 3))
1739 tmp33 = fac*(tij_abcd(1, 1, 3, 3)*qp_j(1, 1) + &
1740 tij_abcd(2, 1, 3, 3)*qp_j(2, 1) + &
1741 tij_abcd(3, 1, 3, 3)*qp_j(3, 1) + &
1742 tij_abcd(1, 2, 3, 3)*qp_j(1, 2) + &
1743 tij_abcd(2, 2, 3, 3)*qp_j(2, 2) + &
1744 tij_abcd(3, 2, 3, 3)*qp_j(3, 2) + &
1745 tij_abcd(1, 3, 3, 3)*qp_j(1, 3) + &
1746 tij_abcd(2, 3, 3, 3)*qp_j(2, 3) + &
1747 tij_abcd(3, 3, 3, 3)*qp_j(3, 3))
1748
1749 ef2_i(1, 1) = ef2_i(1, 1) - tmp11
1750 ef2_i(1, 2) = ef2_i(1, 2) - tmp12
1751 ef2_i(1, 3) = ef2_i(1, 3) - tmp13
1752 ef2_i(2, 1) = ef2_i(2, 1) - tmp12
1753 ef2_i(2, 2) = ef2_i(2, 2) - tmp22
1754 ef2_i(2, 3) = ef2_i(2, 3) - tmp23
1755 ef2_i(3, 1) = ef2_i(3, 1) - tmp13
1756 ef2_i(3, 2) = ef2_i(3, 2) - tmp23
1757 ef2_i(3, 3) = ef2_i(3, 3) - tmp33
1758
1759 tmp11 = fac*(tij_abcd(1, 1, 1, 1)*qp_i(1, 1) + &
1760 tij_abcd(2, 1, 1, 1)*qp_i(2, 1) + &
1761 tij_abcd(3, 1, 1, 1)*qp_i(3, 1) + &
1762 tij_abcd(1, 2, 1, 1)*qp_i(1, 2) + &
1763 tij_abcd(2, 2, 1, 1)*qp_i(2, 2) + &
1764 tij_abcd(3, 2, 1, 1)*qp_i(3, 2) + &
1765 tij_abcd(1, 3, 1, 1)*qp_i(1, 3) + &
1766 tij_abcd(2, 3, 1, 1)*qp_i(2, 3) + &
1767 tij_abcd(3, 3, 1, 1)*qp_i(3, 3))
1768 tmp12 = fac*(tij_abcd(1, 1, 1, 2)*qp_i(1, 1) + &
1769 tij_abcd(2, 1, 1, 2)*qp_i(2, 1) + &
1770 tij_abcd(3, 1, 1, 2)*qp_i(3, 1) + &
1771 tij_abcd(1, 2, 1, 2)*qp_i(1, 2) + &
1772 tij_abcd(2, 2, 1, 2)*qp_i(2, 2) + &
1773 tij_abcd(3, 2, 1, 2)*qp_i(3, 2) + &
1774 tij_abcd(1, 3, 1, 2)*qp_i(1, 3) + &
1775 tij_abcd(2, 3, 1, 2)*qp_i(2, 3) + &
1776 tij_abcd(3, 3, 1, 2)*qp_i(3, 3))
1777 tmp13 = fac*(tij_abcd(1, 1, 1, 3)*qp_i(1, 1) + &
1778 tij_abcd(2, 1, 1, 3)*qp_i(2, 1) + &
1779 tij_abcd(3, 1, 1, 3)*qp_i(3, 1) + &
1780 tij_abcd(1, 2, 1, 3)*qp_i(1, 2) + &
1781 tij_abcd(2, 2, 1, 3)*qp_i(2, 2) + &
1782 tij_abcd(3, 2, 1, 3)*qp_i(3, 2) + &
1783 tij_abcd(1, 3, 1, 3)*qp_i(1, 3) + &
1784 tij_abcd(2, 3, 1, 3)*qp_i(2, 3) + &
1785 tij_abcd(3, 3, 1, 3)*qp_i(3, 3))
1786 tmp22 = fac*(tij_abcd(1, 1, 2, 2)*qp_i(1, 1) + &
1787 tij_abcd(2, 1, 2, 2)*qp_i(2, 1) + &
1788 tij_abcd(3, 1, 2, 2)*qp_i(3, 1) + &
1789 tij_abcd(1, 2, 2, 2)*qp_i(1, 2) + &
1790 tij_abcd(2, 2, 2, 2)*qp_i(2, 2) + &
1791 tij_abcd(3, 2, 2, 2)*qp_i(3, 2) + &
1792 tij_abcd(1, 3, 2, 2)*qp_i(1, 3) + &
1793 tij_abcd(2, 3, 2, 2)*qp_i(2, 3) + &
1794 tij_abcd(3, 3, 2, 2)*qp_i(3, 3))
1795 tmp23 = fac*(tij_abcd(1, 1, 2, 3)*qp_i(1, 1) + &
1796 tij_abcd(2, 1, 2, 3)*qp_i(2, 1) + &
1797 tij_abcd(3, 1, 2, 3)*qp_i(3, 1) + &
1798 tij_abcd(1, 2, 2, 3)*qp_i(1, 2) + &
1799 tij_abcd(2, 2, 2, 3)*qp_i(2, 2) + &
1800 tij_abcd(3, 2, 2, 3)*qp_i(3, 2) + &
1801 tij_abcd(1, 3, 2, 3)*qp_i(1, 3) + &
1802 tij_abcd(2, 3, 2, 3)*qp_i(2, 3) + &
1803 tij_abcd(3, 3, 2, 3)*qp_i(3, 3))
1804 tmp33 = fac*(tij_abcd(1, 1, 3, 3)*qp_i(1, 1) + &
1805 tij_abcd(2, 1, 3, 3)*qp_i(2, 1) + &
1806 tij_abcd(3, 1, 3, 3)*qp_i(3, 1) + &
1807 tij_abcd(1, 2, 3, 3)*qp_i(1, 2) + &
1808 tij_abcd(2, 2, 3, 3)*qp_i(2, 2) + &
1809 tij_abcd(3, 2, 3, 3)*qp_i(3, 2) + &
1810 tij_abcd(1, 3, 3, 3)*qp_i(1, 3) + &
1811 tij_abcd(2, 3, 3, 3)*qp_i(2, 3) + &
1812 tij_abcd(3, 3, 3, 3)*qp_i(3, 3))
1813
1814 ef2_j(1, 1) = ef2_j(1, 1) - tmp11
1815 ef2_j(1, 2) = ef2_j(1, 2) - tmp12
1816 ef2_j(1, 3) = ef2_j(1, 3) - tmp13
1817 ef2_j(2, 1) = ef2_j(2, 1) - tmp12
1818 ef2_j(2, 2) = ef2_j(2, 2) - tmp22
1819 ef2_j(2, 3) = ef2_j(2, 3) - tmp23
1820 ef2_j(3, 1) = ef2_j(3, 1) - tmp13
1821 ef2_j(3, 2) = ef2_j(3, 2) - tmp23
1822 ef2_j(3, 3) = ef2_j(3, 3) - tmp33
1823 END IF
1824 END IF
1825 END IF
1826 IF (task(3, 2)) THEN
1827 ! Quadrupole - Dipole
1828 fac = 1.0_dp/3.0_dp
1829 ! Dipole i (locally B) - Quadrupole j (locally A)
1830 tmp_ij = dp_i(1)*(tij_abc(1, 1, 1)*qp_j(1, 1) + &
1831 tij_abc(2, 1, 1)*qp_j(2, 1) + &
1832 tij_abc(3, 1, 1)*qp_j(3, 1) + &
1833 tij_abc(1, 2, 1)*qp_j(1, 2) + &
1834 tij_abc(2, 2, 1)*qp_j(2, 2) + &
1835 tij_abc(3, 2, 1)*qp_j(3, 2) + &
1836 tij_abc(1, 3, 1)*qp_j(1, 3) + &
1837 tij_abc(2, 3, 1)*qp_j(2, 3) + &
1838 tij_abc(3, 3, 1)*qp_j(3, 3)) + &
1839 dp_i(2)*(tij_abc(1, 1, 2)*qp_j(1, 1) + &
1840 tij_abc(2, 1, 2)*qp_j(2, 1) + &
1841 tij_abc(3, 1, 2)*qp_j(3, 1) + &
1842 tij_abc(1, 2, 2)*qp_j(1, 2) + &
1843 tij_abc(2, 2, 2)*qp_j(2, 2) + &
1844 tij_abc(3, 2, 2)*qp_j(3, 2) + &
1845 tij_abc(1, 3, 2)*qp_j(1, 3) + &
1846 tij_abc(2, 3, 2)*qp_j(2, 3) + &
1847 tij_abc(3, 3, 2)*qp_j(3, 3)) + &
1848 dp_i(3)*(tij_abc(1, 1, 3)*qp_j(1, 1) + &
1849 tij_abc(2, 1, 3)*qp_j(2, 1) + &
1850 tij_abc(3, 1, 3)*qp_j(3, 1) + &
1851 tij_abc(1, 2, 3)*qp_j(1, 2) + &
1852 tij_abc(2, 2, 3)*qp_j(2, 2) + &
1853 tij_abc(3, 2, 3)*qp_j(3, 2) + &
1854 tij_abc(1, 3, 3)*qp_j(1, 3) + &
1855 tij_abc(2, 3, 3)*qp_j(2, 3) + &
1856 tij_abc(3, 3, 3)*qp_j(3, 3))
1857
1858 ! Dipole j (locally A) - Quadrupole i (locally B)
1859 tmp_ji = dp_j(1)*(tij_abc(1, 1, 1)*qp_i(1, 1) + &
1860 tij_abc(2, 1, 1)*qp_i(2, 1) + &
1861 tij_abc(3, 1, 1)*qp_i(3, 1) + &
1862 tij_abc(1, 2, 1)*qp_i(1, 2) + &
1863 tij_abc(2, 2, 1)*qp_i(2, 2) + &
1864 tij_abc(3, 2, 1)*qp_i(3, 2) + &
1865 tij_abc(1, 3, 1)*qp_i(1, 3) + &
1866 tij_abc(2, 3, 1)*qp_i(2, 3) + &
1867 tij_abc(3, 3, 1)*qp_i(3, 3)) + &
1868 dp_j(2)*(tij_abc(1, 1, 2)*qp_i(1, 1) + &
1869 tij_abc(2, 1, 2)*qp_i(2, 1) + &
1870 tij_abc(3, 1, 2)*qp_i(3, 1) + &
1871 tij_abc(1, 2, 2)*qp_i(1, 2) + &
1872 tij_abc(2, 2, 2)*qp_i(2, 2) + &
1873 tij_abc(3, 2, 2)*qp_i(3, 2) + &
1874 tij_abc(1, 3, 2)*qp_i(1, 3) + &
1875 tij_abc(2, 3, 2)*qp_i(2, 3) + &
1876 tij_abc(3, 3, 2)*qp_i(3, 3)) + &
1877 dp_j(3)*(tij_abc(1, 1, 3)*qp_i(1, 1) + &
1878 tij_abc(2, 1, 3)*qp_i(2, 1) + &
1879 tij_abc(3, 1, 3)*qp_i(3, 1) + &
1880 tij_abc(1, 2, 3)*qp_i(1, 2) + &
1881 tij_abc(2, 2, 3)*qp_i(2, 2) + &
1882 tij_abc(3, 2, 3)*qp_i(3, 2) + &
1883 tij_abc(1, 3, 3)*qp_i(1, 3) + &
1884 tij_abc(2, 3, 3)*qp_i(2, 3) + &
1885 tij_abc(3, 3, 3)*qp_i(3, 3))
1886
1887 tmp = fac*(tmp_ij - tmp_ji)
1888 eloc = eloc + tmp
1889 IF (do_forces .OR. do_stress) THEN
1890 DO k = 1, 3
1891 ! Dipole i (locally B) - Quadrupole j (locally A)
1892 tmp_ij = dp_i(1)*(tij_abcd(1, 1, 1, k)*qp_j(1, 1) + &
1893 tij_abcd(2, 1, 1, k)*qp_j(2, 1) + &
1894 tij_abcd(3, 1, 1, k)*qp_j(3, 1) + &
1895 tij_abcd(1, 2, 1, k)*qp_j(1, 2) + &
1896 tij_abcd(2, 2, 1, k)*qp_j(2, 2) + &
1897 tij_abcd(3, 2, 1, k)*qp_j(3, 2) + &
1898 tij_abcd(1, 3, 1, k)*qp_j(1, 3) + &
1899 tij_abcd(2, 3, 1, k)*qp_j(2, 3) + &
1900 tij_abcd(3, 3, 1, k)*qp_j(3, 3)) + &
1901 dp_i(2)*(tij_abcd(1, 1, 2, k)*qp_j(1, 1) + &
1902 tij_abcd(2, 1, 2, k)*qp_j(2, 1) + &
1903 tij_abcd(3, 1, 2, k)*qp_j(3, 1) + &
1904 tij_abcd(1, 2, 2, k)*qp_j(1, 2) + &
1905 tij_abcd(2, 2, 2, k)*qp_j(2, 2) + &
1906 tij_abcd(3, 2, 2, k)*qp_j(3, 2) + &
1907 tij_abcd(1, 3, 2, k)*qp_j(1, 3) + &
1908 tij_abcd(2, 3, 2, k)*qp_j(2, 3) + &
1909 tij_abcd(3, 3, 2, k)*qp_j(3, 3)) + &
1910 dp_i(3)*(tij_abcd(1, 1, 3, k)*qp_j(1, 1) + &
1911 tij_abcd(2, 1, 3, k)*qp_j(2, 1) + &
1912 tij_abcd(3, 1, 3, k)*qp_j(3, 1) + &
1913 tij_abcd(1, 2, 3, k)*qp_j(1, 2) + &
1914 tij_abcd(2, 2, 3, k)*qp_j(2, 2) + &
1915 tij_abcd(3, 2, 3, k)*qp_j(3, 2) + &
1916 tij_abcd(1, 3, 3, k)*qp_j(1, 3) + &
1917 tij_abcd(2, 3, 3, k)*qp_j(2, 3) + &
1918 tij_abcd(3, 3, 3, k)*qp_j(3, 3))
1919
1920 ! Dipole j (locally A) - Quadrupole i (locally B)
1921 tmp_ji = dp_j(1)*(tij_abcd(1, 1, 1, k)*qp_i(1, 1) + &
1922 tij_abcd(2, 1, 1, k)*qp_i(2, 1) + &
1923 tij_abcd(3, 1, 1, k)*qp_i(3, 1) + &
1924 tij_abcd(1, 2, 1, k)*qp_i(1, 2) + &
1925 tij_abcd(2, 2, 1, k)*qp_i(2, 2) + &
1926 tij_abcd(3, 2, 1, k)*qp_i(3, 2) + &
1927 tij_abcd(1, 3, 1, k)*qp_i(1, 3) + &
1928 tij_abcd(2, 3, 1, k)*qp_i(2, 3) + &
1929 tij_abcd(3, 3, 1, k)*qp_i(3, 3)) + &
1930 dp_j(2)*(tij_abcd(1, 1, 2, k)*qp_i(1, 1) + &
1931 tij_abcd(2, 1, 2, k)*qp_i(2, 1) + &
1932 tij_abcd(3, 1, 2, k)*qp_i(3, 1) + &
1933 tij_abcd(1, 2, 2, k)*qp_i(1, 2) + &
1934 tij_abcd(2, 2, 2, k)*qp_i(2, 2) + &
1935 tij_abcd(3, 2, 2, k)*qp_i(3, 2) + &
1936 tij_abcd(1, 3, 2, k)*qp_i(1, 3) + &
1937 tij_abcd(2, 3, 2, k)*qp_i(2, 3) + &
1938 tij_abcd(3, 3, 2, k)*qp_i(3, 3)) + &
1939 dp_j(3)*(tij_abcd(1, 1, 3, k)*qp_i(1, 1) + &
1940 tij_abcd(2, 1, 3, k)*qp_i(2, 1) + &
1941 tij_abcd(3, 1, 3, k)*qp_i(3, 1) + &
1942 tij_abcd(1, 2, 3, k)*qp_i(1, 2) + &
1943 tij_abcd(2, 2, 3, k)*qp_i(2, 2) + &
1944 tij_abcd(3, 2, 3, k)*qp_i(3, 2) + &
1945 tij_abcd(1, 3, 3, k)*qp_i(1, 3) + &
1946 tij_abcd(2, 3, 3, k)*qp_i(2, 3) + &
1947 tij_abcd(3, 3, 3, k)*qp_i(3, 3))
1948
1949 fr(k) = fr(k) - fac*(tmp_ij - tmp_ji)
1950 END DO
1951 END IF
1952 END IF
1953 IF (task(3, 1)) THEN
1954 ! Quadrupole - Charge
1955 fac = 1.0_dp/3.0_dp
1956
1957 ! Quadrupole j (locally A) - Charge j (locally B)
1958 tmp_ij = ch_i*(tij_ab(1, 1)*qp_j(1, 1) + &
1959 tij_ab(2, 1)*qp_j(2, 1) + &
1960 tij_ab(3, 1)*qp_j(3, 1) + &
1961 tij_ab(1, 2)*qp_j(1, 2) + &
1962 tij_ab(2, 2)*qp_j(2, 2) + &
1963 tij_ab(3, 2)*qp_j(3, 2) + &
1964 tij_ab(1, 3)*qp_j(1, 3) + &
1965 tij_ab(2, 3)*qp_j(2, 3) + &
1966 tij_ab(3, 3)*qp_j(3, 3))
1967
1968 ! Quadrupole i (locally B) - Charge j (locally A)
1969 tmp_ji = ch_j*(tij_ab(1, 1)*qp_i(1, 1) + &
1970 tij_ab(2, 1)*qp_i(2, 1) + &
1971 tij_ab(3, 1)*qp_i(3, 1) + &
1972 tij_ab(1, 2)*qp_i(1, 2) + &
1973 tij_ab(2, 2)*qp_i(2, 2) + &
1974 tij_ab(3, 2)*qp_i(3, 2) + &
1975 tij_ab(1, 3)*qp_i(1, 3) + &
1976 tij_ab(2, 3)*qp_i(2, 3) + &
1977 tij_ab(3, 3)*qp_i(3, 3))
1978
1979 eloc = eloc + fac*(tmp_ij + tmp_ji)
1980 IF (do_forces .OR. do_stress) THEN
1981 DO k = 1, 3
1982 ! Quadrupole j (locally A) - Charge i (locally B)
1983 tmp_ij = ch_i*(tij_abc(1, 1, k)*qp_j(1, 1) + &
1984 tij_abc(2, 1, k)*qp_j(2, 1) + &
1985 tij_abc(3, 1, k)*qp_j(3, 1) + &
1986 tij_abc(1, 2, k)*qp_j(1, 2) + &
1987 tij_abc(2, 2, k)*qp_j(2, 2) + &
1988 tij_abc(3, 2, k)*qp_j(3, 2) + &
1989 tij_abc(1, 3, k)*qp_j(1, 3) + &
1990 tij_abc(2, 3, k)*qp_j(2, 3) + &
1991 tij_abc(3, 3, k)*qp_j(3, 3))
1992
1993 ! Quadrupole i (locally B) - Charge j (locally A)
1994 tmp_ji = ch_j*(tij_abc(1, 1, k)*qp_i(1, 1) + &
1995 tij_abc(2, 1, k)*qp_i(2, 1) + &
1996 tij_abc(3, 1, k)*qp_i(3, 1) + &
1997 tij_abc(1, 2, k)*qp_i(1, 2) + &
1998 tij_abc(2, 2, k)*qp_i(2, 2) + &
1999 tij_abc(3, 2, k)*qp_i(3, 2) + &
2000 tij_abc(1, 3, k)*qp_i(1, 3) + &
2001 tij_abc(2, 3, k)*qp_i(2, 3) + &
2002 tij_abc(3, 3, k)*qp_i(3, 3))
2003
2004 fr(k) = fr(k) - fac*(tmp_ij + tmp_ji)
2005 END DO
2006 END IF
2007 END IF
2008 ! Electric fields
2009 IF (do_efield) THEN
2010 ! Potential
2011 IF (do_efield0) THEN
2012 efield0(atom_a) = efield0(atom_a) + ef0_j
2013
2014 efield0(atom_b) = efield0(atom_b) + ef0_i
2015 END IF
2016 ! Electric field
2017 IF (do_efield1) THEN
2018 efield1(1, atom_a) = efield1(1, atom_a) + ef1_j(1)
2019 efield1(2, atom_a) = efield1(2, atom_a) + ef1_j(2)
2020 efield1(3, atom_a) = efield1(3, atom_a) + ef1_j(3)
2021
2022 efield1(1, atom_b) = efield1(1, atom_b) + ef1_i(1)
2023 efield1(2, atom_b) = efield1(2, atom_b) + ef1_i(2)
2024 efield1(3, atom_b) = efield1(3, atom_b) + ef1_i(3)
2025 END IF
2026 ! Electric field gradient
2027 IF (do_efield2) THEN
2028 efield2(1, atom_a) = efield2(1, atom_a) + ef2_j(1, 1)
2029 efield2(2, atom_a) = efield2(2, atom_a) + ef2_j(1, 2)
2030 efield2(3, atom_a) = efield2(3, atom_a) + ef2_j(1, 3)
2031 efield2(4, atom_a) = efield2(4, atom_a) + ef2_j(2, 1)
2032 efield2(5, atom_a) = efield2(5, atom_a) + ef2_j(2, 2)
2033 efield2(6, atom_a) = efield2(6, atom_a) + ef2_j(2, 3)
2034 efield2(7, atom_a) = efield2(7, atom_a) + ef2_j(3, 1)
2035 efield2(8, atom_a) = efield2(8, atom_a) + ef2_j(3, 2)
2036 efield2(9, atom_a) = efield2(9, atom_a) + ef2_j(3, 3)
2037
2038 efield2(1, atom_b) = efield2(1, atom_b) + ef2_i(1, 1)
2039 efield2(2, atom_b) = efield2(2, atom_b) + ef2_i(1, 2)
2040 efield2(3, atom_b) = efield2(3, atom_b) + ef2_i(1, 3)
2041 efield2(4, atom_b) = efield2(4, atom_b) + ef2_i(2, 1)
2042 efield2(5, atom_b) = efield2(5, atom_b) + ef2_i(2, 2)
2043 efield2(6, atom_b) = efield2(6, atom_b) + ef2_i(2, 3)
2044 efield2(7, atom_b) = efield2(7, atom_b) + ef2_i(3, 1)
2045 efield2(8, atom_b) = efield2(8, atom_b) + ef2_i(3, 2)
2046 efield2(9, atom_b) = efield2(9, atom_b) + ef2_i(3, 3)
2047 END IF
2048 END IF
2049 IF (do_stress) THEN
2050 ptens11 = ptens11 + rab(1)*fr(1)
2051 ptens21 = ptens21 + rab(2)*fr(1)
2052 ptens31 = ptens31 + rab(3)*fr(1)
2053 ptens12 = ptens12 + rab(1)*fr(2)
2054 ptens22 = ptens22 + rab(2)*fr(2)
2055 ptens32 = ptens32 + rab(3)*fr(2)
2056 ptens13 = ptens13 + rab(1)*fr(3)
2057 ptens23 = ptens23 + rab(2)*fr(3)
2058 ptens33 = ptens33 + rab(3)*fr(3)
2059 END IF
2060
2061
2062 IF (PRESENT(integral_value)) THEN
2063 integral_value = eloc
2064 END IF
2065 IF (do_forces) THEN
2066 force_ab = fr
2067 END IF
2068 END SUBROUTINE se_coulomb_ij_interaction
2069
2070END MODULE se_fock_matrix_integrals
2071
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition grid_common.h:51
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
semi_empirical_store_int_types::semi_empirical_si_type
Semi-empirical store integrals type.
Definition semi_empirical_store_int_types.F:43
semi_empirical_types::se_int_control_type
Definition semi_empirical_types.F:137
semi_empirical_types::se_taper_type
Taper type use in semi-empirical calculations.
Definition semi_empirical_types.F:158
semi_empirical_types::semi_empirical_type
Semi-empirical type.
Definition semi_empirical_types.F:57