(git:374b731)
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-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! **************************************************************************************************
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 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
2072END 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
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:136
semi_empirical_types::se_taper_type
Taper type use in semi-empirical calculations.
Definition semi_empirical_types.F:157
semi_empirical_types::semi_empirical_type
Semi-empirical type.
Definition semi_empirical_types.F:57