Line data Source code
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 : ! **************************************************************************************************
18 : MODULE se_fock_matrix_integrals
19 :
20 : USE kinds, ONLY: dp
21 : USE semi_empirical_int_arrays, ONLY: se_orbital_pointer
22 : USE semi_empirical_integrals, ONLY: drotint, &
23 : drotnuc, &
24 : rotint, &
25 : rotnuc
26 : USE semi_empirical_store_int_types, ONLY: semi_empirical_si_type
27 : USE semi_empirical_types, ONLY: se_int_control_type, &
28 : se_taper_type, &
29 : semi_empirical_type
30 : #include "./base/base_uses.f90"
31 :
32 : #:include "ewalds_multipole_sr.fypp"
33 :
34 : IMPLICIT NONE
35 : PRIVATE
36 :
37 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'se_fock_matrix_integrals'
38 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
39 :
40 : PUBLIC :: fock2_1el, dfock2_1el, fock1_2el, fock2_1el_ew, fock2C_ew, &
41 : fock2C, dfock2C, fock2E, dfock2E, fock2_1el_r3, dfock2_1el_r3, &
42 : fock2C_r3, dfock2C_r3, se_coulomb_ij_interaction
43 :
44 : CONTAINS
45 :
46 : ! **************************************************************************************************
47 : !> \brief Construction of 2-center 1-electron Fock Matrix
48 : !> \param sepi ...
49 : !> \param sepj ...
50 : !> \param rij ...
51 : !> \param ksi_block DIMENSION(sepi%natorb, sepi%natorb)
52 : !> \param ksj_block DIMENSION(sepi%natorb, sepi%natorb)
53 : !> \param pi_block ...
54 : !> \param pj_block ...
55 : !> \param ecore ...
56 : !> \param itype ...
57 : !> \param anag ...
58 : !> \param se_int_control ...
59 : !> \param se_taper ...
60 : !> \param store_int_env ...
61 : !> \date 04.2008 [tlaino]
62 : !> \author Teodoro Laino [tlaino] - University of Zurich
63 : ! **************************************************************************************************
64 12353504 : SUBROUTINE fock2_1el(sepi, sepj, rij, ksi_block, ksj_block, pi_block, pj_block, &
65 : ecore, itype, anag, se_int_control, se_taper, store_int_env)
66 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
67 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rij
68 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ksi_block, ksj_block
69 : REAL(KIND=dp), &
70 : DIMENSION(sepi%natorb, sepi%natorb), INTENT(IN) :: pi_block
71 : REAL(KIND=dp), &
72 : DIMENSION(sepj%natorb, sepj%natorb), INTENT(IN) :: pj_block
73 : REAL(KIND=dp), DIMENSION(2), INTENT(INOUT) :: ecore
74 : INTEGER, INTENT(IN) :: itype
75 : LOGICAL, INTENT(IN) :: anag
76 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
77 : TYPE(se_taper_type), POINTER :: se_taper
78 : TYPE(semi_empirical_si_type), POINTER :: store_int_env
79 :
80 : INTEGER :: i1, i1L, i2, j1, j1L
81 : REAL(KIND=dp), DIMENSION(45) :: e1b, e2a
82 :
83 : ! Compute integrals
84 :
85 : CALL rotnuc(sepi, sepj, rij, e1b=e1b, e2a=e2a, itype=itype, anag=anag, &
86 12353504 : se_int_control=se_int_control, se_taper=se_taper, store_int_env=store_int_env)
87 : !
88 : ! Add the electron-nuclear attraction term for atom sepi
89 : !
90 12353504 : i2 = 0
91 49380712 : DO i1L = 1, sepi%natorb
92 37027208 : i1 = se_orbital_pointer(i1L)
93 88776416 : DO j1L = 1, i1L - 1
94 51749208 : j1 = se_orbital_pointer(j1L)
95 51749208 : i2 = i2 + 1
96 51749208 : ksi_block(i1, j1) = ksi_block(i1, j1) + e1b(i2)
97 51749208 : ksi_block(j1, i1) = ksi_block(i1, j1)
98 88776416 : ecore(1) = ecore(1) + 2.0_dp*e1b(i2)*pi_block(i1, j1)
99 : END DO
100 37027208 : j1 = se_orbital_pointer(j1L)
101 37027208 : i2 = i2 + 1
102 37027208 : ksi_block(i1, j1) = ksi_block(i1, j1) + e1b(i2)
103 49380712 : ecore(1) = ecore(1) + e1b(i2)*pi_block(i1, j1)
104 : END DO
105 : !
106 : ! Add the electron-nuclear attraction term for atom sepj
107 : !
108 12353504 : i2 = 0
109 49370282 : DO i1L = 1, sepj%natorb
110 37016778 : i1 = se_orbital_pointer(i1L)
111 88668066 : DO j1L = 1, i1L - 1
112 51651288 : j1 = se_orbital_pointer(j1L)
113 51651288 : i2 = i2 + 1
114 51651288 : ksj_block(i1, j1) = ksj_block(i1, j1) + e2a(i2)
115 51651288 : ksj_block(j1, i1) = ksj_block(i1, j1)
116 88668066 : ecore(2) = ecore(2) + 2.0_dp*e2a(i2)*pj_block(i1, j1)
117 : END DO
118 37016778 : j1 = se_orbital_pointer(j1L)
119 37016778 : i2 = i2 + 1
120 37016778 : ksj_block(i1, j1) = ksj_block(i1, j1) + e2a(i2)
121 49370282 : ecore(2) = ecore(2) + e2a(i2)*pj_block(i1, j1)
122 : END DO
123 :
124 12353504 : END SUBROUTINE fock2_1el
125 :
126 : ! **************************************************************************************************
127 : !> \brief Derivatives of 2-center 1-electron Fock Matrix
128 : !> \param sepi ...
129 : !> \param sepj ...
130 : !> \param rij ...
131 : !> \param pi_block ...
132 : !> \param pj_block ...
133 : !> \param itype ...
134 : !> \param anag ...
135 : !> \param se_int_control ...
136 : !> \param se_taper ...
137 : !> \param force ...
138 : !> \param delta ...
139 : !> \date 04.2008 [tlaino]
140 : !> \author Teodoro Laino [tlaino] - University of Zurich
141 : ! **************************************************************************************************
142 346191 : SUBROUTINE dfock2_1el(sepi, sepj, rij, pi_block, pj_block, itype, anag, &
143 : se_int_control, se_taper, force, delta)
144 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
145 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rij
146 : REAL(KIND=dp), &
147 : DIMENSION(sepi%natorb, sepi%natorb), INTENT(IN) :: pi_block
148 : REAL(KIND=dp), &
149 : DIMENSION(sepj%natorb, sepj%natorb), INTENT(IN) :: pj_block
150 : INTEGER, INTENT(IN) :: itype
151 : LOGICAL, INTENT(IN) :: anag
152 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
153 : TYPE(se_taper_type), POINTER :: se_taper
154 : REAL(KIND=dp), DIMENSION(3), INTENT(INOUT) :: force
155 : REAL(KIND=dp), INTENT(IN) :: delta
156 :
157 : INTEGER :: i1, i1L, i2, j1, j1L
158 : REAL(KIND=dp) :: tmp
159 : REAL(KIND=dp), DIMENSION(3, 45) :: de1b, de2a
160 :
161 : ! Compute integrals
162 :
163 : CALL drotnuc(sepi, sepj, rij, de1b=de1b, de2a=de2a, itype=itype, anag=anag, &
164 346191 : se_int_control=se_int_control, se_taper=se_taper, delta=delta)
165 : !
166 : ! Add the electron-nuclear attraction term for atom sepi
167 : !
168 346191 : i2 = 0
169 1291557 : DO i1L = 1, sepi%natorb
170 945366 : i1 = se_orbital_pointer(i1L)
171 2347296 : DO j1L = 1, i1L - 1
172 1401930 : j1 = se_orbital_pointer(j1L)
173 1401930 : i2 = i2 + 1
174 1401930 : tmp = 2.0_dp*pi_block(i1, j1)
175 1401930 : force(1) = force(1) + de1b(1, i2)*tmp
176 1401930 : force(2) = force(2) + de1b(2, i2)*tmp
177 2347296 : force(3) = force(3) + de1b(3, i2)*tmp
178 : END DO
179 945366 : j1 = se_orbital_pointer(j1L)
180 945366 : i2 = i2 + 1
181 945366 : force(1) = force(1) + de1b(1, i2)*pi_block(i1, j1)
182 945366 : force(2) = force(2) + de1b(2, i2)*pi_block(i1, j1)
183 1291557 : force(3) = force(3) + de1b(3, i2)*pi_block(i1, j1)
184 : END DO
185 : !
186 : ! Add the electron-nuclear attraction term for atom sepj
187 : !
188 346191 : i2 = 0
189 1288971 : DO i1L = 1, sepj%natorb
190 942780 : i1 = se_orbital_pointer(i1L)
191 2333658 : DO j1L = 1, i1L - 1
192 1390878 : j1 = se_orbital_pointer(j1L)
193 1390878 : i2 = i2 + 1
194 1390878 : tmp = 2.0_dp*pj_block(i1, j1)
195 1390878 : force(1) = force(1) + de2a(1, i2)*tmp
196 1390878 : force(2) = force(2) + de2a(2, i2)*tmp
197 2333658 : force(3) = force(3) + de2a(3, i2)*tmp
198 : END DO
199 942780 : j1 = se_orbital_pointer(j1L)
200 942780 : i2 = i2 + 1
201 942780 : force(1) = force(1) + de2a(1, i2)*pj_block(i1, j1)
202 942780 : force(2) = force(2) + de2a(2, i2)*pj_block(i1, j1)
203 1288971 : force(3) = force(3) + de2a(3, i2)*pj_block(i1, j1)
204 : END DO
205 :
206 346191 : END SUBROUTINE dfock2_1el
207 :
208 : ! **************************************************************************************************
209 : !> \brief Construction of 1-center 2-electron Fock Matrix
210 : !> \param sep ...
211 : !> \param p_tot ...
212 : !> \param p_mat ...
213 : !> \param f_mat DIMENSION(sep%natorb, sep%natorb)
214 : !> \param factor ...
215 : !> \date 04.2008 [tlaino]
216 : !> \author Teodoro Laino [tlaino] - University of Zurich
217 : ! **************************************************************************************************
218 197959 : SUBROUTINE fock1_2el(sep, p_tot, p_mat, f_mat, factor)
219 : TYPE(semi_empirical_type), POINTER :: sep
220 : REAL(KIND=dp), DIMENSION(45, 45), INTENT(IN) :: p_tot
221 : REAL(KIND=dp), DIMENSION(sep%natorb, sep%natorb), &
222 : INTENT(IN) :: p_mat
223 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: f_mat
224 : REAL(KIND=dp), INTENT(IN) :: factor
225 :
226 : INTEGER :: i, ijw, ikw, iL, im, j, jL, jlw, jm, k, &
227 : kL, klw, l, lL
228 : REAL(KIND=dp) :: sum
229 :
230 : ! One-center coulomb and exchange terms for semiempirical_type sep
231 : !
232 : ! F(i,j)=F(i,j)+sum(k,l)((PA(k,l)+PB(k,l))*<i,j|k,l>
233 : ! -(PA(k,l) )*<i,k|j,l>), k,l on type sep.
234 : !
235 :
236 747884 : DO iL = 1, sep%natorb
237 549925 : i = se_orbital_pointer(iL)
238 2157201 : DO jL = 1, iL
239 1409317 : j = se_orbital_pointer(jL)
240 :
241 : ! `J' Address IJ in W
242 1409317 : ijw = (iL*(iL - 1))/2 + jL
243 1409317 : sum = 0.0_dp
244 8514734 : DO kL = 1, sep%natorb
245 7105417 : k = se_orbital_pointer(kL)
246 52395951 : DO lL = 1, sep%natorb
247 43881217 : l = se_orbital_pointer(lL)
248 :
249 : ! `J' Address KL in W
250 43881217 : im = MAX(kL, lL)
251 43881217 : jm = MIN(kL, lL)
252 43881217 : klw = (im*(im - 1))/2 + jm
253 :
254 : ! `K' Address IK in W
255 43881217 : im = MAX(kL, jL)
256 43881217 : jm = MIN(kL, jL)
257 43881217 : ikw = (im*(im - 1))/2 + jm
258 :
259 : ! `K' Address JL in W
260 43881217 : im = MAX(lL, iL)
261 43881217 : jm = MIN(lL, iL)
262 43881217 : jlw = (im*(im - 1))/2 + jm
263 :
264 50986634 : sum = sum + p_tot(k, l)*sep%w(ijw, klw) - p_mat(k, l)*sep%w(ikw, jlw)
265 : END DO
266 : END DO
267 1409317 : f_mat(i, j) = f_mat(i, j) + factor*sum
268 1959242 : f_mat(j, i) = f_mat(i, j)
269 : END DO
270 : END DO
271 197959 : END SUBROUTINE fock1_2el
272 :
273 : ! **************************************************************************************************
274 : !> \brief Construction of 2-center 1-electron Fock Matrix (Ewald self term)
275 : !> \param sep ...
276 : !> \param rij ...
277 : !> \param ks_block DIMENSION(sep%natorb, sep%natorb)
278 : !> \param p_block ...
279 : !> \param ecore ...
280 : !> \param itype ...
281 : !> \param anag ...
282 : !> \param se_int_control ...
283 : !> \param se_taper ...
284 : !> \param store_int_env ...
285 : !> \date 04.2009 [jgh]
286 : !> \author jgh - University of Zurich
287 : ! **************************************************************************************************
288 159 : SUBROUTINE fock2_1el_ew(sep, rij, ks_block, p_block, ecore, itype, anag, &
289 : se_int_control, se_taper, store_int_env)
290 : TYPE(semi_empirical_type), POINTER :: sep
291 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rij
292 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ks_block
293 : REAL(KIND=dp), DIMENSION(sep%natorb, sep%natorb), &
294 : INTENT(IN) :: p_block
295 : REAL(KIND=dp), INTENT(INOUT) :: ecore
296 : INTEGER, INTENT(IN) :: itype
297 : LOGICAL, INTENT(IN) :: anag
298 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
299 : TYPE(se_taper_type), POINTER :: se_taper
300 : TYPE(semi_empirical_si_type), POINTER :: store_int_env
301 :
302 : INTEGER :: i1, i1L, i2, j1, j1L, n
303 : REAL(KIND=dp), DIMENSION(45) :: e1b, e2a
304 :
305 : ! Compute integrals
306 :
307 : CALL rotnuc(sep, sep, rij, e1b=e1b, e2a=e2a, itype=itype, anag=anag, &
308 159 : se_int_control=se_int_control, se_taper=se_taper, store_int_env=store_int_env)
309 : !
310 : ! Add the electron-nuclear attraction term for atom sep
311 : ! e1b == e2a
312 : !
313 159 : n = (sep%natorb*(sep%natorb + 1))/2
314 159 : i2 = 0
315 795 : DO i1L = 1, sep%natorb
316 636 : i1 = se_orbital_pointer(i1L)
317 1590 : DO j1L = 1, i1L - 1
318 954 : j1 = se_orbital_pointer(j1L)
319 954 : i2 = i2 + 1
320 954 : ks_block(i1, j1) = ks_block(i1, j1) + e1b(i2)
321 954 : ks_block(j1, i1) = ks_block(i1, j1)
322 1590 : ecore = ecore + 2._dp*e1b(i2)*p_block(i1, j1)
323 : END DO
324 : ! i1L == j1L
325 636 : i2 = i2 + 1
326 636 : ks_block(i1, i1) = ks_block(i1, i1) + e1b(i2)
327 795 : ecore = ecore + e1b(i2)*p_block(i1, i1)
328 : END DO
329 :
330 159 : END SUBROUTINE fock2_1el_ew
331 :
332 : ! **************************************************************************************************
333 : !> \brief Construction of 2-center Fock Matrix - Coulomb Self Terms (Ewald)
334 : !> \param sep ...
335 : !> \param rij ...
336 : !> \param p_tot ...
337 : !> \param f_mat DIMENSION(sep%natorb, sep%natorb)
338 : !> \param factor ...
339 : !> \param anag ...
340 : !> \param se_int_control ...
341 : !> \param se_taper ...
342 : !> \param store_int_env ...
343 : !> \date 04.2009 [jgh]
344 : !> \author jgh - University of Zurich
345 : ! **************************************************************************************************
346 159 : SUBROUTINE fock2C_ew(sep, rij, p_tot, f_mat, factor, anag, se_int_control, &
347 : se_taper, store_int_env)
348 : TYPE(semi_empirical_type), POINTER :: sep
349 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rij
350 : REAL(KIND=dp), DIMENSION(45, 45), INTENT(IN) :: p_tot
351 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: f_mat
352 : REAL(KIND=dp), INTENT(IN) :: factor
353 : LOGICAL, INTENT(IN) :: anag
354 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
355 : TYPE(se_taper_type), POINTER :: se_taper
356 : TYPE(semi_empirical_si_type), POINTER :: store_int_env
357 :
358 : INTEGER :: i, iL, j, jL, k, kL, kr, l, lL, natorb
359 : REAL(KIND=dp) :: a, aa, bb
360 : REAL(KIND=dp), DIMENSION(2025) :: w
361 :
362 : ! Evaluate integrals
363 :
364 : CALL rotint(sep, sep, rij, w, anag=anag, se_int_control=se_int_control, &
365 159 : se_taper=se_taper, store_int_env=store_int_env)
366 159 : kr = 0
367 159 : natorb = sep%natorb
368 795 : DO iL = 1, natorb
369 636 : i = se_orbital_pointer(iL)
370 636 : aa = 2.0_dp
371 2385 : DO jL = 1, iL
372 1590 : j = se_orbital_pointer(jL)
373 1590 : IF (i == j) THEN
374 636 : aa = 1.0_dp
375 : END IF
376 8586 : DO kL = 1, natorb
377 6360 : k = se_orbital_pointer(kL)
378 6360 : bb = 2.0_dp
379 23850 : DO lL = 1, kL
380 15900 : l = se_orbital_pointer(lL)
381 15900 : IF (k == l) THEN
382 6360 : bb = 1.0_dp
383 : END IF
384 15900 : kr = kr + 1
385 15900 : a = 0.5_dp*w(kr)*factor
386 : ! Coulomb
387 15900 : f_mat(i, j) = f_mat(i, j) + bb*a*p_tot(k, l)
388 15900 : f_mat(k, l) = f_mat(k, l) + aa*a*p_tot(i, j)
389 15900 : f_mat(j, i) = f_mat(i, j)
390 22260 : f_mat(l, k) = f_mat(k, l)
391 : END DO
392 : END DO
393 : END DO
394 : END DO
395 :
396 159 : END SUBROUTINE fock2C_ew
397 :
398 : ! **************************************************************************************************
399 : !> \brief Construction of 2-center Fock Matrix - Coulomb Terms
400 : !> \param sepi ...
401 : !> \param sepj ...
402 : !> \param rij ...
403 : !> \param switch ...
404 : !> \param pi_tot ...
405 : !> \param fi_mat DIMENSION(sepi%natorb, sepi%natorb)
406 : !> \param pj_tot DIMENSION(sepj%natorb, sepj%natorb)
407 : !> \param fj_mat ...
408 : !> \param factor ...
409 : !> \param anag ...
410 : !> \param se_int_control ...
411 : !> \param se_taper ...
412 : !> \param store_int_env ...
413 : !> \date 04.2008 [tlaino]
414 : !> \author Teodoro Laino [tlaino] - University of Zurich
415 : ! **************************************************************************************************
416 12353504 : SUBROUTINE fock2C(sepi, sepj, rij, switch, pi_tot, fi_mat, pj_tot, fj_mat, &
417 : factor, anag, se_int_control, se_taper, store_int_env)
418 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
419 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rij
420 : LOGICAL, INTENT(IN) :: switch
421 : REAL(KIND=dp), DIMENSION(45, 45), INTENT(IN) :: pi_tot
422 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: fi_mat
423 : REAL(KIND=dp), DIMENSION(45, 45), INTENT(IN) :: pj_tot
424 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: fj_mat
425 : REAL(KIND=dp), INTENT(IN) :: factor
426 : LOGICAL, INTENT(IN) :: anag
427 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
428 : TYPE(se_taper_type), POINTER :: se_taper
429 : TYPE(semi_empirical_si_type), POINTER :: store_int_env
430 :
431 : INTEGER :: i, iL, j, jL, k, kL, kr, l, lL, natorb(2)
432 : REAL(KIND=dp) :: a, aa, bb, irij(3)
433 : REAL(KIND=dp), DIMENSION(2025) :: w
434 :
435 : ! Evaluate integrals
436 :
437 12353504 : IF (.NOT. switch) THEN
438 : CALL rotint(sepi, sepj, rij, w, anag=anag, se_int_control=se_int_control, &
439 6533574 : se_taper=se_taper, store_int_env=store_int_env)
440 : ELSE
441 23279720 : irij = -rij
442 : CALL rotint(sepj, sepi, irij, w, anag=anag, se_int_control=se_int_control, &
443 5819930 : se_taper=se_taper, store_int_env=store_int_env)
444 : END IF
445 12353504 : kr = 0
446 12353504 : natorb(1) = sepi%natorb
447 12353504 : natorb(2) = sepj%natorb
448 12353504 : IF (switch) THEN
449 5819930 : natorb(1) = sepj%natorb
450 5819930 : natorb(2) = sepi%natorb
451 : END IF
452 51482370 : DO iL = 1, natorb(1)
453 39128866 : i = se_orbital_pointer(iL)
454 39128866 : aa = 2.0_dp
455 147156460 : DO jL = 1, iL
456 95674090 : j = se_orbital_pointer(jL)
457 95674090 : IF (i == j) THEN
458 39128866 : aa = 1.0_dp
459 : END IF
460 434812247 : DO kL = 1, natorb(2)
461 300009291 : k = se_orbital_pointer(kL)
462 300009291 : bb = 2.0_dp
463 1126318234 : DO lL = 1, kL
464 730634853 : l = se_orbital_pointer(lL)
465 730634853 : IF (k == l) THEN
466 300009291 : bb = 1.0_dp
467 : END IF
468 730634853 : kr = kr + 1
469 730634853 : a = w(kr)*factor
470 : ! Coulomb
471 1030644144 : IF (.NOT. switch) THEN
472 395654418 : fi_mat(i, j) = fi_mat(i, j) + bb*a*pj_tot(k, l)
473 395654418 : fj_mat(k, l) = fj_mat(k, l) + aa*a*pi_tot(i, j)
474 395654418 : fi_mat(j, i) = fi_mat(i, j)
475 395654418 : fj_mat(l, k) = fj_mat(k, l)
476 : ELSE
477 334980435 : fj_mat(i, j) = fj_mat(i, j) + bb*a*pi_tot(k, l)
478 334980435 : fi_mat(k, l) = fi_mat(k, l) + aa*a*pj_tot(i, j)
479 334980435 : fj_mat(j, i) = fj_mat(i, j)
480 334980435 : fi_mat(l, k) = fi_mat(k, l)
481 : END IF
482 : END DO
483 : END DO
484 : END DO
485 : END DO
486 :
487 12353504 : END SUBROUTINE fock2C
488 :
489 : ! **************************************************************************************************
490 : !> \brief Derivatives of 2-center Fock Matrix - Coulomb Terms
491 : !> \param sepi ...
492 : !> \param sepj ...
493 : !> \param rij ...
494 : !> \param switch ...
495 : !> \param pi_tot ...
496 : !> \param pj_tot ...
497 : !> \param factor ...
498 : !> \param anag ...
499 : !> \param se_int_control ...
500 : !> \param se_taper ...
501 : !> \param force ...
502 : !> \param delta ...
503 : !> \date 04.2008 [tlaino]
504 : !> \author Teodoro Laino [tlaino] - University of Zurich
505 : ! **************************************************************************************************
506 346191 : SUBROUTINE dfock2C(sepi, sepj, rij, switch, pi_tot, pj_tot, factor, anag, &
507 : se_int_control, se_taper, force, delta)
508 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
509 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rij
510 : LOGICAL, INTENT(IN) :: switch
511 : REAL(KIND=dp), DIMENSION(45, 45), INTENT(IN) :: pi_tot, pj_tot
512 : REAL(KIND=dp), INTENT(IN) :: factor
513 : LOGICAL, INTENT(IN) :: anag
514 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
515 : TYPE(se_taper_type), POINTER :: se_taper
516 : REAL(KIND=dp), DIMENSION(3), INTENT(INOUT) :: force
517 : REAL(KIND=dp), INTENT(IN) :: delta
518 :
519 : INTEGER :: i, iL, j, jL, k, kL, kr, l, lL, natorb(2)
520 : REAL(KIND=dp) :: aa, bb, tmp
521 : REAL(KIND=dp), DIMENSION(3) :: a, irij
522 : REAL(KIND=dp), DIMENSION(3, 2025) :: dw
523 :
524 : ! Evaluate integrals' derivatives
525 :
526 346191 : IF (.NOT. switch) THEN
527 : CALL drotint(sepi, sepj, rij, dw, delta, anag=anag, se_int_control=se_int_control, &
528 175380 : se_taper=se_taper)
529 : ELSE
530 683244 : irij = -rij
531 : CALL drotint(sepj, sepi, irij, dw, delta, anag=anag, se_int_control=se_int_control, &
532 170811 : se_taper=se_taper)
533 : END IF
534 :
535 346191 : kr = 0
536 346191 : natorb(1) = sepi%natorb
537 346191 : natorb(2) = sepj%natorb
538 346191 : IF (switch) THEN
539 170811 : natorb(1) = sepj%natorb
540 170811 : natorb(2) = sepi%natorb
541 : END IF
542 1379193 : DO iL = 1, natorb(1)
543 1033002 : i = se_orbital_pointer(iL)
544 1033002 : aa = 2.0_dp
545 3998217 : DO jL = 1, iL
546 2619024 : j = se_orbital_pointer(jL)
547 2619024 : IF (i == j) THEN
548 1033002 : aa = 1.0_dp
549 : END IF
550 11576268 : DO kL = 1, natorb(2)
551 7924242 : k = se_orbital_pointer(kL)
552 7924242 : bb = 2.0_dp
553 32601384 : DO lL = 1, kL
554 22058118 : l = se_orbital_pointer(lL)
555 22058118 : IF (k == l) THEN
556 7924242 : bb = 1.0_dp
557 : END IF
558 22058118 : kr = kr + 1
559 22058118 : a(1) = dw(1, kr)*factor
560 22058118 : a(2) = dw(2, kr)*factor
561 22058118 : a(3) = dw(3, kr)*factor
562 : ! Coulomb
563 22058118 : IF (.NOT. switch) THEN
564 11248759 : tmp = bb*aa*pj_tot(k, l)*pi_tot(i, j)
565 : ELSE
566 10809359 : tmp = bb*aa*pi_tot(k, l)*pj_tot(i, j)
567 : END IF
568 22058118 : force(1) = force(1) + a(1)*tmp
569 22058118 : force(2) = force(2) + a(2)*tmp
570 29982360 : force(3) = force(3) + a(3)*tmp
571 : END DO
572 : END DO
573 : END DO
574 : END DO
575 346191 : END SUBROUTINE dfock2C
576 :
577 : ! **************************************************************************************************
578 : !> \brief Construction of 2-center Fock Matrix - General Driver
579 : !> \param sepi ...
580 : !> \param sepj ...
581 : !> \param rij ...
582 : !> \param switch ...
583 : !> \param isize ...
584 : !> \param pi_mat ...
585 : !> \param fi_mat DIMENSION(isize(1), isize(2))
586 : !> \param factor ...
587 : !> \param anag ...
588 : !> \param se_int_control ...
589 : !> \param se_taper ...
590 : !> \param store_int_env ...
591 : !> \date 04.2008 [tlaino]
592 : !> \author Teodoro Laino [tlaino] - University of Zurich
593 : ! **************************************************************************************************
594 3614293 : SUBROUTINE fock2E(sepi, sepj, rij, switch, isize, pi_mat, fi_mat, factor, &
595 : anag, se_int_control, se_taper, store_int_env)
596 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
597 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rij
598 : LOGICAL, INTENT(IN) :: switch
599 : INTEGER, DIMENSION(2), INTENT(IN) :: isize
600 : REAL(KIND=dp), DIMENSION(isize(1), isize(2)), &
601 : INTENT(IN) :: pi_mat
602 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: fi_mat
603 : REAL(KIND=dp), INTENT(IN) :: factor
604 : LOGICAL, INTENT(IN) :: anag
605 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
606 : TYPE(se_taper_type), POINTER :: se_taper
607 : TYPE(semi_empirical_si_type), POINTER :: store_int_env
608 :
609 : INTEGER :: i, iL, j, jL, k, kL, kr, l, lL, natorb(2)
610 : REAL(KIND=dp) :: a, aa, bb, irij(3)
611 : REAL(KIND=dp), DIMENSION(2025) :: w
612 :
613 : ! Evaluate integrals
614 :
615 3614293 : IF (.NOT. switch) THEN
616 : CALL rotint(sepi, sepj, rij, w, anag=anag, se_int_control=se_int_control, &
617 1760733 : se_taper=se_taper, store_int_env=store_int_env)
618 : ELSE
619 7414240 : irij = -rij
620 : CALL rotint(sepj, sepi, irij, w, anag=anag, se_int_control=se_int_control, &
621 1853560 : se_taper=se_taper, store_int_env=store_int_env)
622 : END IF
623 3614293 : kr = 0
624 3614293 : natorb(1) = sepi%natorb
625 3614293 : natorb(2) = sepj%natorb
626 3614293 : IF (switch) THEN
627 1853560 : natorb(1) = sepj%natorb
628 1853560 : natorb(2) = sepi%natorb
629 : END IF
630 13798951 : DO iL = 1, natorb(1)
631 10184658 : i = se_orbital_pointer(iL)
632 10184658 : aa = 2.0_dp
633 38990839 : DO jL = 1, iL
634 25191888 : j = se_orbital_pointer(jL)
635 25191888 : IF (i == j) THEN
636 10184658 : aa = 1.0_dp
637 : END IF
638 104278015 : DO kL = 1, natorb(2)
639 68901469 : k = se_orbital_pointer(kL)
640 68901469 : bb = 2.0_dp
641 261419708 : DO lL = 1, kL
642 167326351 : l = se_orbital_pointer(lL)
643 167326351 : IF (k == l) THEN
644 68901469 : bb = 1.0_dp
645 : END IF
646 167326351 : kr = kr + 1
647 167326351 : a = w(kr)*factor
648 : ! Exchange
649 167326351 : a = a*aa*bb*0.25_dp
650 167326351 : fi_mat(i, k) = fi_mat(i, k) - a*pi_mat(j, l)
651 167326351 : fi_mat(i, l) = fi_mat(i, l) - a*pi_mat(j, k)
652 167326351 : fi_mat(j, k) = fi_mat(j, k) - a*pi_mat(i, l)
653 236227820 : fi_mat(j, l) = fi_mat(j, l) - a*pi_mat(i, k)
654 : END DO
655 : END DO
656 : END DO
657 : END DO
658 :
659 3614293 : END SUBROUTINE fock2E
660 :
661 : ! **************************************************************************************************
662 : !> \brief Derivatives of 2-center Fock Matrix - General Driver
663 : !> \param sepi ...
664 : !> \param sepj ...
665 : !> \param rij ...
666 : !> \param switch ...
667 : !> \param isize ...
668 : !> \param pi_mat ...
669 : !> \param factor ...
670 : !> \param anag ...
671 : !> \param se_int_control ...
672 : !> \param se_taper ...
673 : !> \param force ...
674 : !> \param delta ...
675 : !> \date 04.2008 [tlaino]
676 : !> \author Teodoro Laino [tlaino] - University of Zurich
677 : ! **************************************************************************************************
678 172433 : SUBROUTINE dfock2E(sepi, sepj, rij, switch, isize, pi_mat, factor, anag, &
679 : se_int_control, se_taper, force, delta)
680 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
681 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rij
682 : LOGICAL, INTENT(IN) :: switch
683 : INTEGER, DIMENSION(2), INTENT(IN) :: isize
684 : REAL(KIND=dp), DIMENSION(isize(1), isize(2)), &
685 : INTENT(IN) :: pi_mat
686 : REAL(KIND=dp), INTENT(IN) :: factor
687 : LOGICAL, INTENT(IN) :: anag
688 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
689 : TYPE(se_taper_type), POINTER :: se_taper
690 : REAL(KIND=dp), DIMENSION(3), INTENT(INOUT) :: force
691 : REAL(KIND=dp), INTENT(IN) :: delta
692 :
693 : INTEGER :: i, iL, j, jL, k, kL, kr, l, lL, natorb(2)
694 : REAL(KIND=dp) :: aa, bb, tmp, tmp1, tmp2, tmp3, tmp4
695 : REAL(KIND=dp), DIMENSION(3) :: a, irij
696 : REAL(KIND=dp), DIMENSION(3, 2025) :: dw
697 :
698 : ! Evaluate integrals' derivatives
699 :
700 172433 : IF (.NOT. switch) THEN
701 : CALL drotint(sepi, sepj, rij, dw, delta, anag=anag, se_int_control=se_int_control, &
702 83540 : se_taper=se_taper)
703 : ELSE
704 355572 : irij = -rij
705 : CALL drotint(sepj, sepi, irij, dw, delta, anag=anag, se_int_control=se_int_control, &
706 88893 : se_taper=se_taper)
707 : END IF
708 :
709 172433 : kr = 0
710 172433 : natorb(1) = sepi%natorb
711 172433 : natorb(2) = sepj%natorb
712 172433 : IF (switch) THEN
713 88893 : natorb(1) = sepj%natorb
714 88893 : natorb(2) = sepi%natorb
715 : END IF
716 668784 : DO iL = 1, natorb(1)
717 496351 : i = se_orbital_pointer(iL)
718 496351 : aa = 2.0_dp
719 1906051 : DO jL = 1, iL
720 1237267 : j = se_orbital_pointer(jL)
721 1237267 : IF (i == j) THEN
722 496351 : aa = 1.0_dp
723 : END IF
724 5189554 : DO kL = 1, natorb(2)
725 3455936 : k = se_orbital_pointer(kL)
726 3455936 : bb = 2.0_dp
727 13265477 : DO lL = 1, kL
728 8572274 : l = se_orbital_pointer(lL)
729 8572274 : IF (k == l) THEN
730 3455936 : bb = 1.0_dp
731 : END IF
732 8572274 : kr = kr + 1
733 8572274 : tmp = factor*aa*bb*0.25_dp
734 8572274 : a(1) = dw(1, kr)*tmp
735 8572274 : a(2) = dw(2, kr)*tmp
736 8572274 : a(3) = dw(3, kr)*tmp
737 : ! Exchange
738 8572274 : tmp1 = pi_mat(j, l)*pi_mat(i, k)
739 8572274 : tmp2 = pi_mat(j, k)*pi_mat(i, l)
740 8572274 : tmp3 = pi_mat(i, l)*pi_mat(j, k)
741 8572274 : tmp4 = pi_mat(i, k)*pi_mat(j, l)
742 :
743 8572274 : force(1) = force(1) - a(1)*tmp1
744 8572274 : force(1) = force(1) - a(1)*tmp2
745 8572274 : force(1) = force(1) - a(1)*tmp3
746 8572274 : force(1) = force(1) - a(1)*tmp4
747 :
748 8572274 : force(2) = force(2) - a(2)*tmp1
749 8572274 : force(2) = force(2) - a(2)*tmp2
750 8572274 : force(2) = force(2) - a(2)*tmp3
751 8572274 : force(2) = force(2) - a(2)*tmp4
752 :
753 8572274 : force(3) = force(3) - a(3)*tmp1
754 8572274 : force(3) = force(3) - a(3)*tmp2
755 8572274 : force(3) = force(3) - a(3)*tmp3
756 12028210 : force(3) = force(3) - a(3)*tmp4
757 : END DO
758 : END DO
759 : END DO
760 : END DO
761 172433 : END SUBROUTINE dfock2E
762 :
763 : ! **************************************************************************************************
764 : !> \brief Construction of 2-center 1-electron Fock Matrix for the residual
765 : !> (1/R^3) integral part
766 : !> \param sepi ...
767 : !> \param sepj ...
768 : !> \param ksi_block DIMENSION(sepi%natorb, sepi%natorb)
769 : !> \param ksj_block DIMENSION(sepj%natorb, sepj%natorb)
770 : !> \param pi_block ...
771 : !> \param pj_block ...
772 : !> \param e1b ...
773 : !> \param e2a ...
774 : !> \param ecore ...
775 : !> \param rp ...
776 : !> \date 12.2008 [tlaino]
777 : !> \author Teodoro Laino [tlaino]
778 : ! **************************************************************************************************
779 0 : SUBROUTINE fock2_1el_r3(sepi, sepj, ksi_block, ksj_block, pi_block, pj_block, &
780 0 : e1b, e2a, ecore, rp)
781 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
782 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ksi_block, ksj_block
783 : REAL(KIND=dp), &
784 : DIMENSION(sepi%natorb, sepi%natorb), INTENT(IN) :: pi_block
785 : REAL(KIND=dp), &
786 : DIMENSION(sepj%natorb, sepj%natorb), INTENT(IN) :: pj_block
787 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: e1b, e2a
788 : REAL(KIND=dp), DIMENSION(2), INTENT(INOUT) :: ecore
789 : REAL(KIND=dp), INTENT(IN) :: rp
790 :
791 : INTEGER :: i1, i1L, i2
792 :
793 : !
794 : ! Add the electron-nuclear residual attraction term for atom sepi
795 : !
796 :
797 0 : i2 = 0
798 0 : DO i1L = 1, sepi%natorb
799 0 : i2 = i2 + 1
800 0 : i1 = se_orbital_pointer(i1L)
801 0 : ksi_block(i1, i1) = ksi_block(i1, i1) + e1b(i2)*rp
802 0 : ecore(1) = ecore(1) + e1b(i2)*rp*pi_block(i1, i1)
803 : END DO
804 : !
805 : ! Add the electron-nuclear residual attraction term for atom sepj
806 : !
807 0 : i2 = 0
808 0 : DO i1L = 1, sepj%natorb
809 0 : i2 = i2 + 1
810 0 : i1 = se_orbital_pointer(i1L)
811 0 : ksj_block(i1, i1) = ksj_block(i1, i1) + e2a(i2)*rp
812 0 : ecore(2) = ecore(2) + e2a(i2)*rp*pj_block(i1, i1)
813 : END DO
814 :
815 0 : END SUBROUTINE fock2_1el_r3
816 :
817 : ! **************************************************************************************************
818 : !> \brief Derivatives of 2-center 1-electron Fock Matrix residual (1/R^3)
819 : !> integral part
820 : !> \param sepi ...
821 : !> \param sepj ...
822 : !> \param drp ...
823 : !> \param pi_block ...
824 : !> \param pj_block ...
825 : !> \param force ...
826 : !> \param e1b ...
827 : !> \param e2a ...
828 : !> \date 12.2008 [tlaino]
829 : !> \author Teodoro Laino [tlaino]
830 : ! **************************************************************************************************
831 0 : SUBROUTINE dfock2_1el_r3(sepi, sepj, drp, pi_block, pj_block, force, e1b, e2a)
832 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
833 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: drp
834 : REAL(KIND=dp), &
835 : DIMENSION(sepi%natorb, sepi%natorb), INTENT(IN) :: pi_block
836 : REAL(KIND=dp), &
837 : DIMENSION(sepj%natorb, sepj%natorb), INTENT(IN) :: pj_block
838 : REAL(KIND=dp), DIMENSION(3), INTENT(INOUT) :: force
839 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: e1b, e2a
840 :
841 : INTEGER :: i1, i1L, i2
842 : REAL(KIND=dp) :: tmp
843 :
844 : !
845 : ! Add the electron-nuclear residual attraction term for atom sepi
846 : !
847 :
848 0 : i2 = 0
849 0 : DO i1L = 1, sepi%natorb
850 0 : i1 = se_orbital_pointer(i1L)
851 0 : i2 = i2 + 1
852 0 : tmp = e1b(i2)*pi_block(i1, i1)
853 0 : force(1) = force(1) + tmp*drp(1)
854 0 : force(2) = force(2) + tmp*drp(2)
855 0 : force(3) = force(3) + tmp*drp(3)
856 : END DO
857 : !
858 : ! Add the electron-nuclear attraction term for atom sepj
859 : !
860 : i2 = 0
861 0 : DO i1L = 1, sepj%natorb
862 0 : i1 = se_orbital_pointer(i1L)
863 0 : i2 = i2 + 1
864 0 : tmp = e2a(i2)*pj_block(i1, i1)
865 0 : force(1) = force(1) + tmp*drp(1)
866 0 : force(2) = force(2) + tmp*drp(2)
867 0 : force(3) = force(3) + tmp*drp(3)
868 : END DO
869 :
870 0 : END SUBROUTINE dfock2_1el_r3
871 :
872 : ! **************************************************************************************************
873 : !> \brief Construction of 2-center Fock Matrix - Coulomb Terms for the residual
874 : !> (1/R^3) integral part
875 : !> \param sepi ...
876 : !> \param sepj ...
877 : !> \param switch ...
878 : !> \param pi_tot ...
879 : !> \param fi_mat DIMENSION(sepi%natorb, sepi%natorb)
880 : !> \param pj_tot ...
881 : !> \param fj_mat DIMENSION(sepj%natorb, sepj%natorb)
882 : !> \param factor ...
883 : !> \param w ...
884 : !> \param rp ...
885 : !> \date 12.2008 [tlaino]
886 : !> \author Teodoro Laino [tlaino]
887 : ! **************************************************************************************************
888 0 : SUBROUTINE fock2C_r3(sepi, sepj, switch, pi_tot, fi_mat, pj_tot, fj_mat, &
889 : factor, w, rp)
890 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
891 : LOGICAL, INTENT(IN) :: switch
892 : REAL(KIND=dp), DIMENSION(45, 45), INTENT(IN) :: pi_tot
893 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: fi_mat
894 : REAL(KIND=dp), DIMENSION(45, 45), INTENT(IN) :: pj_tot
895 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: fj_mat
896 : REAL(KIND=dp), INTENT(IN) :: factor
897 : REAL(KIND=dp), DIMENSION(81), INTENT(IN) :: w
898 : REAL(KIND=dp), INTENT(IN) :: rp
899 :
900 : INTEGER :: i, iL, ind, j, k, kL, kr, natorb(2)
901 : REAL(KIND=dp) :: a, w_l(81)
902 :
903 0 : natorb(1) = sepi%natorb
904 0 : natorb(2) = sepj%natorb
905 0 : IF (switch) THEN
906 0 : natorb(1) = sepj%natorb
907 0 : natorb(2) = sepi%natorb
908 : ! Reshuffle the integral array (natural storage order is sepi/sepj)
909 : kr = 0
910 0 : DO i = 1, sepj%natorb
911 0 : DO j = 1, sepi%natorb
912 0 : kr = kr + 1
913 0 : ind = (j - 1)*sepj%natorb + i
914 0 : w_l(kr) = w(ind)
915 : END DO
916 : END DO
917 : ELSE
918 0 : w_l = w
919 : END IF
920 :
921 : ! Modify the Fock Matrix
922 0 : kr = 0
923 0 : DO iL = 1, natorb(1)
924 0 : i = se_orbital_pointer(iL)
925 0 : DO kL = 1, natorb(2)
926 0 : k = se_orbital_pointer(kL)
927 0 : kr = kr + 1
928 0 : a = w_l(kr)*factor*rp
929 : ! Coulomb
930 0 : IF (.NOT. switch) THEN
931 0 : fi_mat(i, i) = fi_mat(i, i) + a*pj_tot(k, k)
932 0 : fj_mat(k, k) = fj_mat(k, k) + a*pi_tot(i, i)
933 : ELSE
934 0 : fj_mat(i, i) = fj_mat(i, i) + a*pi_tot(k, k)
935 0 : fi_mat(k, k) = fi_mat(k, k) + a*pj_tot(i, i)
936 : END IF
937 : END DO
938 : END DO
939 :
940 0 : END SUBROUTINE fock2C_r3
941 :
942 : ! **************************************************************************************************
943 : !> \brief Derivatives of 2-center Fock Matrix - Coulomb Terms for the residual
944 : !> (1/R^3) integral part
945 : !> \param sepi ...
946 : !> \param sepj ...
947 : !> \param switch ...
948 : !> \param pi_tot ...
949 : !> \param pj_tot ...
950 : !> \param factor ...
951 : !> \param w ...
952 : !> \param drp ...
953 : !> \param force ...
954 : !> \date 12.2008 [tlaino]
955 : !> \author Teodoro Laino [tlaino]
956 : ! **************************************************************************************************
957 0 : SUBROUTINE dfock2C_r3(sepi, sepj, switch, pi_tot, pj_tot, factor, w, drp, &
958 : force)
959 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
960 : LOGICAL, INTENT(IN) :: switch
961 : REAL(KIND=dp), DIMENSION(45, 45), INTENT(IN) :: pi_tot, pj_tot
962 : REAL(KIND=dp), INTENT(IN) :: factor
963 : REAL(KIND=dp), DIMENSION(81), INTENT(IN) :: w
964 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: drp
965 : REAL(KIND=dp), DIMENSION(3), INTENT(INOUT) :: force
966 :
967 : INTEGER :: i, iL, ind, j, k, kL, kr, natorb(2)
968 : REAL(KIND=dp) :: a(3), tmp, w_l(81)
969 :
970 0 : natorb(1) = sepi%natorb
971 0 : natorb(2) = sepj%natorb
972 0 : IF (switch) THEN
973 0 : natorb(1) = sepj%natorb
974 0 : natorb(2) = sepi%natorb
975 : ! Reshuffle the integral array (natural storage order is sepi/sepj)
976 : kr = 0
977 0 : DO i = 1, sepj%natorb
978 0 : DO j = 1, sepi%natorb
979 0 : kr = kr + 1
980 0 : ind = (j - 1)*sepj%natorb + i
981 0 : w_l(kr) = w(ind)
982 : END DO
983 : END DO
984 : ELSE
985 0 : w_l = w
986 : END IF
987 :
988 : ! Modify the Fock Matrix
989 0 : kr = 0
990 0 : DO iL = 1, natorb(1)
991 0 : i = se_orbital_pointer(iL)
992 0 : DO kL = 1, natorb(2)
993 0 : k = se_orbital_pointer(kL)
994 0 : kr = kr + 1
995 0 : tmp = w_l(kr)*factor
996 0 : a(1) = tmp*drp(1)
997 0 : a(2) = tmp*drp(2)
998 0 : a(3) = tmp*drp(3)
999 : ! Coulomb
1000 0 : IF (.NOT. switch) THEN
1001 0 : tmp = pj_tot(k, k)*pi_tot(i, i)
1002 : ELSE
1003 0 : tmp = pi_tot(k, k)*pj_tot(i, i)
1004 : END IF
1005 0 : force(1) = force(1) + a(1)*tmp
1006 0 : force(2) = force(2) + a(2)*tmp
1007 0 : force(3) = force(3) + a(3)*tmp
1008 : END DO
1009 : END DO
1010 :
1011 0 : END SUBROUTINE dfock2C_r3
1012 :
1013 : ! **************************************************************************************************
1014 : !> \brief Coulomb interaction multipolar correction
1015 : !> \param atom_a ...
1016 : !> \param atom_b ...
1017 : !> \param my_task ...
1018 : !> \param do_forces ...
1019 : !> \param do_efield ...
1020 : !> \param do_stress ...
1021 : !> \param charges ...
1022 : !> \param dipoles ...
1023 : !> \param quadrupoles ...
1024 : !> \param force_ab ...
1025 : !> \param efield0 ...
1026 : !> \param efield1 ...
1027 : !> \param efield2 ...
1028 : !> \param rab2 ...
1029 : !> \param rab ...
1030 : !> \param integral_value ...
1031 : !> \param ptens11 ...
1032 : !> \param ptens12 ...
1033 : !> \param ptens13 ...
1034 : !> \param ptens21 ...
1035 : !> \param ptens22 ...
1036 : !> \param ptens23 ...
1037 : !> \param ptens31 ...
1038 : !> \param ptens32 ...
1039 : !> \param ptens33 ...
1040 : !> \date 05.2009 [tlaino]
1041 : !> \author Teodoro Laino [tlaino]
1042 : ! **************************************************************************************************
1043 1523155 : SUBROUTINE se_coulomb_ij_interaction(atom_a, atom_b, my_task, do_forces, do_efield, &
1044 : do_stress, charges, dipoles, quadrupoles, force_ab, efield0, efield1, efield2, &
1045 : rab2, rab, integral_value, ptens11, ptens12, ptens13, ptens21, ptens22, ptens23, &
1046 : ptens31, ptens32, ptens33)
1047 : INTEGER, INTENT(IN) :: atom_a, atom_b
1048 : LOGICAL, DIMENSION(3) :: my_task
1049 : LOGICAL, INTENT(IN) :: do_forces, do_efield, do_stress
1050 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges
1051 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dipoles
1052 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: quadrupoles
1053 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: force_ab
1054 : REAL(KIND=dp), DIMENSION(:), POINTER :: efield0
1055 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: efield1, efield2
1056 : REAL(KIND=dp), INTENT(IN) :: rab2
1057 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
1058 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: integral_value
1059 : REAL(KIND=dp), INTENT(INOUT) :: ptens11, ptens12, ptens13, ptens21, &
1060 : ptens22, ptens23, ptens31, ptens32, &
1061 : ptens33
1062 :
1063 : INTEGER :: a, b, c, d, e, i, j, k
1064 : LOGICAL :: do_efield0, do_efield1, do_efield2, &
1065 : force_eval
1066 : LOGICAL, DIMENSION(3) :: do_task
1067 : LOGICAL, DIMENSION(3, 3) :: task
1068 : REAL(KIND=dp) :: ch_i, ch_j, ef0_i, ef0_j, eloc, energy, fac, fac_ij, ir, irab2, r, tij, &
1069 : tmp, tmp1, tmp11, tmp12, tmp13, tmp2, tmp21, tmp22, tmp23, tmp31, tmp32, tmp33, tmp_ij, &
1070 : tmp_ji
1071 : REAL(KIND=dp), DIMENSION(0:5) :: f
1072 : REAL(KIND=dp), DIMENSION(3) :: dp_i, dp_j, ef1_i, ef1_j, fr, tij_a
1073 : REAL(KIND=dp), DIMENSION(3, 3) :: ef2_i, ef2_j, qp_i, qp_j, tij_ab
1074 : REAL(KIND=dp), DIMENSION(3, 3, 3) :: tij_abc
1075 : REAL(KIND=dp), DIMENSION(3, 3, 3, 3) :: tij_abcd
1076 : REAL(KIND=dp), DIMENSION(3, 3, 3, 3, 3) :: tij_abcde
1077 :
1078 1523155 : do_task = my_task
1079 1523155 : energy = 0.0_dp
1080 6092620 : DO i = 1, 3
1081 6092620 : IF (do_task(i)) THEN
1082 908068 : SELECT CASE (i)
1083 : CASE (1)
1084 943309 : do_task(1) = (charges(atom_a) /= 0.0_dp) .OR. (charges(atom_b) /= 0.0_dp)
1085 : CASE (2)
1086 3438523 : do_task(2) = (ANY(dipoles(:, atom_a) /= 0.0_dp)) .OR. (ANY(dipoles(:, atom_b) /= 0.0_dp))
1087 : CASE (3)
1088 12846024 : do_task(3) = (ANY(quadrupoles(:, :, atom_a) /= 0.0_dp)) .OR. (ANY(quadrupoles(:, :, atom_b) /= 0.0_dp))
1089 : END SELECT
1090 : END IF
1091 : END DO
1092 6092620 : DO i = 1, 3
1093 15231550 : DO j = i, 3
1094 9138930 : task(j, i) = do_task(i) .AND. do_task(j)
1095 13708395 : task(i, j) = task(j, i)
1096 : END DO
1097 : END DO
1098 1523155 : do_efield0 = do_efield .AND. ASSOCIATED(efield0)
1099 1523155 : do_efield1 = do_efield .AND. ASSOCIATED(efield1)
1100 1523155 : do_efield2 = do_efield .AND. ASSOCIATED(efield2)
1101 :
1102 1523155 : fac_ij = 1.0_dp
1103 1523155 : IF (atom_a == atom_b) fac_ij = 0.5_dp
1104 :
1105 990370874 : $: ewalds_multipole_sr_macro(mode="PURE_COULOMB")
1106 :
1107 1523155 : IF (PRESENT(integral_value)) THEN
1108 0 : integral_value = eloc
1109 : END IF
1110 1523155 : IF (do_forces) THEN
1111 43876 : force_ab = fr
1112 : END IF
1113 1523155 : END SUBROUTINE se_coulomb_ij_interaction
1114 :
1115 : END MODULE se_fock_matrix_integrals
1116 :
|
