Line data Source code
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 Utilities for Integrals for semi-empiric methods
10 : !> \author Teodoro Laino (03.2008) [tlaino] - University of Zurich
11 : ! **************************************************************************************************
12 : MODULE semi_empirical_int_utils
13 :
14 : USE input_constants, ONLY: do_method_pchg, &
15 : do_se_IS_kdso_d
16 : USE kinds, ONLY: dp
17 : USE semi_empirical_int3_utils, ONLY: charg_int_3, &
18 : dcharg_int_3, &
19 : ijkl_low_3
20 : USE semi_empirical_int_arrays, ONLY: &
21 : CLMp, CLMxx, CLMxy, CLMyy, CLMz, CLMzp, CLMzz, clm_d, clm_sp, ijkl_ind, indexa, indexb, &
22 : int2c_type
23 : USE semi_empirical_types, ONLY: rotmat_type, &
24 : se_int_control_type, &
25 : se_int_screen_type, &
26 : se_taper_type, &
27 : semi_empirical_type
28 : #include "./base/base_uses.f90"
29 :
30 : IMPLICIT NONE
31 : PRIVATE
32 :
33 : #:include 'semi_empirical_int_debug.fypp'
34 :
35 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
36 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_int_utils'
37 :
38 : PUBLIC :: ijkl_sp, ijkl_d, rotmat, rot_2el_2c_first, store_2el_2c_diag, &
39 : d_ijkl_sp, d_ijkl_d
40 :
41 : ABSTRACT INTERFACE
42 : ! **************************************************************************************************
43 : !> \brief ...
44 : !> \param r ...
45 : !> \param l1_i ...
46 : !> \param l2_i ...
47 : !> \param m1_i ...
48 : !> \param m2_i ...
49 : !> \param da_i ...
50 : !> \param db_i ...
51 : !> \param add0 ...
52 : !> \param fact_screen ...
53 : !> \return ...
54 : ! **************************************************************************************************
55 : FUNCTION eval_func_sp(r, l1_i, l2_i, m1_i, m2_i, da_i, db_i, add0, fact_screen) RESULT(charg)
56 : USE kinds, ONLY: dp
57 : REAL(KIND=dp), INTENT(IN) :: r
58 : INTEGER, INTENT(IN) :: l1_i, l2_i, m1_i, m2_i
59 : REAL(KIND=dp), INTENT(IN) :: da_i, db_i, add0, fact_screen
60 : REAL(KIND=dp) :: charg
61 :
62 : END FUNCTION eval_func_sp
63 : END INTERFACE
64 :
65 : ABSTRACT INTERFACE
66 : ! **************************************************************************************************
67 : !> \brief ...
68 : !> \param r ...
69 : !> \param l1_i ...
70 : !> \param l2_i ...
71 : !> \param m ...
72 : !> \param da_i ...
73 : !> \param db_i ...
74 : !> \param add0 ...
75 : !> \param fact_screen ...
76 : !> \return ...
77 : ! **************************************************************************************************
78 : FUNCTION eval_func_d(r, l1_i, l2_i, m, da_i, db_i, add0, fact_screen) RESULT(charg)
79 : USE kinds, ONLY: dp
80 : REAL(KIND=dp), INTENT(IN) :: r
81 : INTEGER, INTENT(IN) :: l1_i, l2_i, m
82 : REAL(KIND=dp), INTENT(IN) :: da_i, db_i, add0, fact_screen
83 : REAL(KIND=dp) :: charg
84 :
85 : END FUNCTION eval_func_d
86 : END INTERFACE
87 :
88 : CONTAINS
89 :
90 : ! **************************************************************************************************
91 : !> \brief General driver for computing semi-empirical integrals <ij|kl> with
92 : !> sp basis set. This code uses the old definitions of quadrupoles and
93 : !> therefore cannot be used for integrals involving d-orbitals (which
94 : !> require a definition of quadrupoles based on the rotational invariant
95 : !> property)
96 : !>
97 : !> \param sepi ...
98 : !> \param sepj ...
99 : !> \param ij ...
100 : !> \param kl ...
101 : !> \param li ...
102 : !> \param lj ...
103 : !> \param lk ...
104 : !> \param ll ...
105 : !> \param ic ...
106 : !> \param r ...
107 : !> \param se_int_control ...
108 : !> \param se_int_screen ...
109 : !> \param itype ...
110 : !> \return ...
111 : !> \date 04.2008 [tlaino]
112 : !> \author Teodoro Laino [tlaino] - University of Zurich
113 : ! **************************************************************************************************
114 92391421 : FUNCTION ijkl_sp(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_control, &
115 : se_int_screen, itype) RESULT(res)
116 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
117 : INTEGER, INTENT(IN) :: ij, kl, li, lj, lk, ll, ic
118 : REAL(KIND=dp), INTENT(IN) :: r
119 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
120 : TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen
121 : INTEGER, INTENT(IN) :: itype
122 : REAL(KIND=dp) :: res
123 :
124 : res = ijkl_sp_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, &
125 : se_int_control%integral_screening, se_int_control%shortrange, &
126 : se_int_control%pc_coulomb_int, se_int_control%max_multipole, &
127 92391421 : itype, charg_int_nri)
128 :
129 : ! If only the shortrange component is requested we can skip the rest
130 92391421 : IF ((.NOT. se_int_control%pc_coulomb_int) .AND. (itype /= do_method_pchg)) THEN
131 : ! Handle the 1/r^3 term, this term is ALWAYS false for KDSO-D integrals
132 23342132 : IF (se_int_control%shortrange .AND. se_int_control%do_ewald_r3) THEN
133 : res = res - ijkl_low_3(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, &
134 0 : itype, charg_int_3)
135 : END IF
136 : END IF
137 92391421 : END FUNCTION ijkl_sp
138 :
139 : ! **************************************************************************************************
140 : !> \brief General driver for computing derivatives of semi-empirical integrals
141 : !> <ij|kl> with sp basis set.
142 : !> This code uses the old definitions of quadrupoles and therefore
143 : !> cannot be used for integrals involving d-orbitals (which requires a
144 : !> definition of quadrupoles based on the rotational invariant property)
145 : !>
146 : !> \param sepi ...
147 : !> \param sepj ...
148 : !> \param ij ...
149 : !> \param kl ...
150 : !> \param li ...
151 : !> \param lj ...
152 : !> \param lk ...
153 : !> \param ll ...
154 : !> \param ic ...
155 : !> \param r ...
156 : !> \param se_int_control ...
157 : !> \param se_int_screen ...
158 : !> \param itype ...
159 : !> \return ...
160 : !> \date 05.2008 [tlaino]
161 : !> \author Teodoro Laino [tlaino] - University of Zurich
162 : ! **************************************************************************************************
163 40138828 : FUNCTION d_ijkl_sp(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_control, &
164 : se_int_screen, itype) RESULT(res)
165 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
166 : INTEGER, INTENT(IN) :: ij, kl, li, lj, lk, ll, ic
167 : REAL(KIND=dp), INTENT(IN) :: r
168 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
169 : TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen
170 : INTEGER, INTENT(IN) :: itype
171 : REAL(KIND=dp) :: res
172 :
173 : REAL(KIND=dp) :: dfs, srd
174 :
175 40138828 : IF (se_int_control%integral_screening == do_se_IS_kdso_d) THEN
176 : res = ijkl_sp_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, &
177 : se_int_control%integral_screening, .FALSE., &
178 : se_int_control%pc_coulomb_int, se_int_control%max_multipole, &
179 2232298 : itype, dcharg_int_nri)
180 :
181 2232298 : IF (.NOT. se_int_control%pc_coulomb_int) THEN
182 : dfs = ijkl_sp_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, &
183 : se_int_control%integral_screening, .FALSE., .FALSE., &
184 1662932 : se_int_control%max_multipole, itype, dcharg_int_nri_fs)
185 1662932 : res = res + dfs*se_int_screen%dft
186 :
187 : ! In case we need the shortrange part we have to evaluate an additional derivative
188 : ! to handle the derivative of the Tapering term
189 1662932 : IF (se_int_control%shortrange) THEN
190 : srd = ijkl_sp_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, &
191 : se_int_control%integral_screening, .FALSE., .TRUE., &
192 1438711 : se_int_control%max_multipole, itype, dcharg_int_nri)
193 1438711 : res = res - srd
194 : END IF
195 : END IF
196 : ELSE
197 : res = ijkl_sp_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, &
198 : se_int_control%integral_screening, se_int_control%shortrange, &
199 : se_int_control%pc_coulomb_int, se_int_control%max_multipole, &
200 37906530 : itype, dcharg_int_nri)
201 : END IF
202 :
203 : ! If only the shortrange component is requested we can skip the rest
204 40138828 : IF ((.NOT. se_int_control%pc_coulomb_int) .AND. (itype /= do_method_pchg)) THEN
205 : ! Handle the 1/r^3 term, this term is ALWAYS false for KDSO-D integrals
206 6310902 : IF (se_int_control%shortrange .AND. se_int_control%do_ewald_r3) THEN
207 : res = res - ijkl_low_3(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, &
208 0 : itype, dcharg_int_3)
209 : END IF
210 : END IF
211 :
212 40138828 : END FUNCTION d_ijkl_sp
213 :
214 : ! **************************************************************************************************
215 : !> \brief Low level general driver for computing semi-empirical integrals
216 : !> <ij|kl> and their derivatives with sp basis set only.
217 : !> This code uses the old definitions of quadrupoles and
218 : !> therefore cannot be used for integrals involving d-orbitals (which
219 : !> require a definition of quadrupoles based on the rotational invariant
220 : !> property)
221 : !>
222 : !> \param sepi ...
223 : !> \param sepj ...
224 : !> \param ij ...
225 : !> \param kl ...
226 : !> \param li ...
227 : !> \param lj ...
228 : !> \param lk ...
229 : !> \param ll ...
230 : !> \param ic ...
231 : !> \param r ...
232 : !> \param se_int_screen ...
233 : !> \param iscreen ...
234 : !> \param shortrange ...
235 : !> \param pc_coulomb_int ...
236 : !> \param max_multipole ...
237 : !> \param itype ...
238 : !> \param eval a function without explicit interface
239 : !> \return ...
240 : !> \date 05.2008 [tlaino]
241 : !> \author Teodoro Laino [tlaino] - University of Zurich
242 : ! **************************************************************************************************
243 135631892 : FUNCTION ijkl_sp_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, &
244 : iscreen, shortrange, pc_coulomb_int, max_multipole, itype, eval) RESULT(res)
245 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
246 : INTEGER, INTENT(IN) :: ij, kl, li, lj, lk, ll, ic
247 : REAL(KIND=dp), INTENT(IN) :: r
248 : TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen
249 : INTEGER, INTENT(IN) :: iscreen
250 : LOGICAL, INTENT(IN) :: shortrange, pc_coulomb_int
251 : INTEGER, INTENT(IN) :: max_multipole, itype
252 :
253 : PROCEDURE(eval_func_sp) :: eval
254 : REAL(KIND=dp) :: res
255 :
256 : INTEGER :: ccc, l1, l1max, l1min, l2, l2max, l2min, &
257 : lij, lkl, lmin, m
258 : REAL(KIND=dp) :: add, chrg, dij, dkl, fact_ij, fact_kl, &
259 : fact_screen, pij, pkl, s1, sum
260 :
261 135631892 : l1min = ABS(li - lj)
262 135631892 : l1max = li + lj
263 135631892 : lij = indexb(li + 1, lj + 1)
264 135631892 : l2min = ABS(lk - ll)
265 135631892 : l2max = lk + ll
266 135631892 : lkl = indexb(lk + 1, ll + 1)
267 :
268 135631892 : l1max = MIN(l1max, 2)
269 135631892 : l1min = MIN(l1min, 2)
270 135631892 : l2max = MIN(l2max, 2)
271 135631892 : l2min = MIN(l2min, 2)
272 135631892 : sum = 0.0_dp
273 135631892 : dij = 0.0_dp
274 135631892 : dkl = 0.0_dp
275 135631892 : fact_ij = 1.0_dp
276 135631892 : fact_kl = 1.0_dp
277 135631892 : fact_screen = 1.0_dp
278 135631892 : IF (lij == 3) fact_ij = SQRT(2.0_dp)
279 135631892 : IF (lkl == 3) fact_kl = SQRT(2.0_dp)
280 135631892 : IF (.NOT. pc_coulomb_int) THEN
281 131096923 : IF (iscreen == do_se_IS_kdso_d) fact_screen = se_int_screen%ft
282 : ! Standard value of the integral
283 320404514 : DO l1 = l1min, l1max
284 189307591 : IF (l1 == 0) THEN
285 116544256 : IF (lij == 1) THEN
286 87438922 : pij = sepi%ko(1)
287 87438922 : IF (ic == -1 .OR. ic == 1) THEN
288 59273813 : pij = sepi%ko(9)
289 : END IF
290 29105334 : ELSE IF (lij == 3) THEN
291 29105334 : pij = sepi%ko(7)
292 : END IF
293 : ELSE
294 72763335 : dij = sepi%cs(lij)*fact_ij
295 72763335 : pij = sepi%ko(lij)
296 : END IF
297 : !
298 541310017 : DO l2 = l2min, l2max
299 220905503 : IF (l2 == 0) THEN
300 181950934 : IF (lkl == 1) THEN
301 166151978 : pkl = sepj%ko(1)
302 166151978 : IF (ic == -1 .OR. ic == 2) THEN
303 133371675 : pkl = sepj%ko(9)
304 : END IF
305 15798956 : ELSE IF (lkl == 3) THEN
306 15798956 : pkl = sepj%ko(7)
307 : END IF
308 : ELSE
309 38954569 : dkl = sepj%cs(lkl)*fact_kl
310 38954569 : pkl = sepj%ko(lkl)
311 : END IF
312 220905503 : IF (itype == do_method_pchg) THEN
313 140755593 : add = 0.0_dp
314 : ELSE
315 80149910 : add = (pij + pkl)**2
316 : END IF
317 220905503 : lmin = MAX(l1, l2)
318 220905503 : s1 = 0.0_dp
319 753969752 : DO m = -lmin, lmin
320 533064249 : ccc = clm_sp(ij, l1, m)*clm_sp(kl, l2, m)
321 753969752 : IF (ABS(ccc) > EPSILON(0.0_dp)) THEN
322 167858898 : chrg = eval(r, l1, l2, clm_sp(ij, l1, m), clm_sp(kl, l2, m), dij, dkl, add, fact_screen)
323 167858898 : s1 = s1 + chrg
324 : END IF
325 : END DO
326 410213094 : sum = sum + s1
327 : END DO
328 : END DO
329 : res = sum
330 : END IF
331 : ! Shortrange: Possibly computes pure Coulomb and subtract from the original integral value
332 135631892 : IF (shortrange .OR. pc_coulomb_int) THEN
333 7651207 : sum = 0.0_dp
334 7651207 : dij = 0.0_dp
335 7651207 : dkl = 0.0_dp
336 7651207 : add = 0.0_dp
337 7651207 : fact_screen = 0.0_dp
338 20119018 : DO l1 = l1min, l1max
339 12467811 : IF (l1 > max_multipole) CYCLE
340 10111884 : IF (l1 /= 0) THEN
341 4993855 : dij = sepi%cs(lij)*fact_ij
342 : END IF
343 : !
344 34650287 : DO l2 = l2min, l2max
345 16887196 : IF (l2 > max_multipole) CYCLE
346 16887196 : IF (l2 /= 0) THEN
347 8352477 : dkl = sepj%cs(lkl)*fact_kl
348 : END IF
349 16887196 : lmin = MAX(l1, l2)
350 16887196 : s1 = 0.0_dp
351 71004126 : DO m = -lmin, lmin
352 54116930 : ccc = clm_sp(ij, l1, m)*clm_sp(kl, l2, m)
353 71004126 : IF (ABS(ccc) > EPSILON(0.0_dp)) THEN
354 9752053 : chrg = eval(r, l1, l2, clm_sp(ij, l1, m), clm_sp(kl, l2, m), dij, dkl, add, fact_screen)
355 9752053 : s1 = s1 + chrg
356 : END IF
357 : END DO
358 29355007 : sum = sum + s1
359 : END DO
360 : END DO
361 7651207 : IF (pc_coulomb_int) res = sum
362 7651207 : IF (shortrange) res = res - sum
363 : END IF
364 135631892 : END FUNCTION ijkl_sp_low
365 :
366 : ! **************************************************************************************************
367 : !> \brief Interaction function between two point-charge configurations NDDO sp-code
368 : !> Non-Rotational Invariant definition of quadrupoles
369 : !> r - Distance r12
370 : !> l1,m - Quantum numbers for multipole of configuration 1
371 : !> l2,m - Quantum numbers for multipole of configuration 2
372 : !> da - charge separation of configuration 1
373 : !> db - charge separation of configuration 2
374 : !> add - additive term
375 : !>
376 : !> \param r ...
377 : !> \param l1_i ...
378 : !> \param l2_i ...
379 : !> \param m1_i ...
380 : !> \param m2_i ...
381 : !> \param da_i ...
382 : !> \param db_i ...
383 : !> \param add0 ...
384 : !> \param fact_screen ...
385 : !> \return ...
386 : !> \date 04.2008 [tlaino]
387 : !> \author Teodoro Laino [tlaino] - University of Zurich
388 : ! **************************************************************************************************
389 120907701 : FUNCTION charg_int_nri(r, l1_i, l2_i, m1_i, m2_i, da_i, db_i, add0, fact_screen) RESULT(charg)
390 : REAL(KIND=dp), INTENT(in) :: r
391 : INTEGER, INTENT(in) :: l1_i, l2_i, m1_i, m2_i
392 : REAL(KIND=dp), INTENT(in) :: da_i, db_i, add0, fact_screen
393 : REAL(KIND=dp) :: charg
394 :
395 : INTEGER :: l1, l2, m1, m2
396 : REAL(KIND=dp) :: add, da, db, dxdx, dxqxz, dzdz, dzqxx, &
397 : dzqzz, fact, qqxx, qqzz, qxxqxx, &
398 : qxxqyy, qxzqxz, xyxy, zzzz
399 :
400 : ! Computing only Integral Values
401 :
402 120907701 : IF (l1_i < l2_i) THEN
403 7594517 : l1 = l1_i
404 7594517 : l2 = l2_i
405 7594517 : m1 = m1_i
406 7594517 : m2 = m2_i
407 7594517 : da = da_i
408 7594517 : db = db_i
409 7594517 : fact = 1.0_dp
410 113313184 : ELSE IF (l1_i > l2_i) THEN
411 28726373 : l1 = l2_i
412 28726373 : l2 = l1_i
413 28726373 : m1 = m2_i
414 28726373 : m2 = m1_i
415 28726373 : da = db_i
416 28726373 : db = da_i
417 28726373 : fact = (-1.0_dp)**(l1 + l2)
418 84586811 : ELSE IF (l1_i == l2_i) THEN
419 84586811 : l1 = l1_i
420 84586811 : l2 = l2_i
421 84586811 : IF (m1_i <= m2_i) THEN
422 84158349 : m1 = m1_i
423 84158349 : m2 = m2_i
424 84158349 : da = da_i
425 84158349 : db = db_i
426 : ELSE
427 428462 : m1 = m2_i
428 428462 : m2 = m1_i
429 428462 : da = db_i
430 428462 : db = da_i
431 : END IF
432 : fact = 1.0_dp
433 : END IF
434 120907701 : add = add0*fact_screen
435 120907701 : charg = 0.0_dp
436 : ! Q - Q.
437 120907701 : IF (l1 == 0 .AND. l2 == 0) THEN
438 80730653 : charg = fact/SQRT(r**2 + add)
439 80730653 : RETURN
440 : END IF
441 : ! Q - Z.
442 40177048 : IF (l1 == 0 .AND. l2 == 1 .AND. m2 == CLMz) THEN
443 10964398 : charg = 1.0_dp/SQRT((r + db)**2 + add) - 1.0_dp/SQRT((r - db)**2 + add)
444 10964398 : charg = charg*0.5_dp*fact
445 10964398 : RETURN
446 : END IF
447 : ! Z - Z.
448 29212650 : IF (l1 == 1 .AND. l2 == 1 .AND. m1 == CLMz .AND. m2 == CLMz) THEN
449 : dzdz = &
450 : +1.0_dp/SQRT((r + da - db)**2 + add) + 1.0_dp/SQRT((r - da + db)**2 + add) &
451 428462 : - 1.0_dp/SQRT((r - da - db)**2 + add) - 1.0_dp/SQRT((r + da + db)**2 + add)
452 428462 : charg = dzdz*0.25_dp*fact
453 428462 : RETURN
454 : END IF
455 : ! X - X
456 28784188 : IF (l1 == 1 .AND. l2 == 1 .AND. m1 == CLMp .AND. m2 == CLMp) THEN
457 428462 : dxdx = 2.0_dp/SQRT(r**2 + (da - db)**2 + add) - 2.0_dp/SQRT(r**2 + (da + db)**2 + add)
458 428462 : charg = dxdx*0.25_dp*fact
459 428462 : RETURN
460 : END IF
461 : ! Q - ZZ
462 28355726 : IF (l1 == 0 .AND. l2 == 2 .AND. m2 == CLMzz) THEN
463 10964398 : qqzz = 1.0_dp/SQRT((r - db)**2 + add) - 2.0_dp/SQRT(r**2 + add) + 1.0_dp/SQRT((r + db)**2 + add)
464 10964398 : charg = qqzz*0.25_dp*fact
465 10964398 : RETURN
466 : END IF
467 : ! Q - XX
468 17391328 : IF (l1 == 0 .AND. l2 == 2 .AND. (m2 == CLMyy .OR. m2 == CLMxx)) THEN
469 11821322 : qqxx = -1.0_dp/SQRT(r**2 + add) + 1.0_dp/SQRT(r**2 + add + db**2)
470 11821322 : charg = qqxx*0.5_dp*fact
471 11821322 : RETURN
472 : END IF
473 : ! Z - ZZ
474 5570006 : IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMz .AND. m2 == CLMzz) THEN
475 : dzqzz = &
476 : +1.0_dp/SQRT((r - da - db)**2 + add) - 2.0_dp/SQRT((r - da)**2 + add) &
477 : + 1.0_dp/SQRT((r - da + db)**2 + add) - 1.0_dp/SQRT((r + da - db)**2 + add) &
478 856924 : + 2.0_dp/SQRT((r + da)**2 + add) - 1.0_dp/SQRT((r + da + db)**2 + add)
479 856924 : charg = dzqzz*0.125_dp*fact
480 856924 : RETURN
481 : END IF
482 : ! Z - XX
483 4713082 : IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMz .AND. (m2 == CLMyy .OR. m2 == CLMxx)) THEN
484 : dzqxx = &
485 : +1.0_dp/SQRT((r + da)**2 + add) - 1.0_dp/SQRT((r + da)**2 + add + db**2) &
486 856924 : - 1.0_dp/SQRT((r - da)**2 + add) + 1.0_dp/SQRT((r - da)**2 + add + db**2)
487 856924 : charg = dzqxx*0.25_dp*fact
488 856924 : RETURN
489 : END IF
490 : ! ZZ - ZZ
491 3856158 : IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzz .AND. m2 == CLMzz) THEN
492 : zzzz = &
493 : +1.0_dp/SQRT((r - da - db)**2 + add) + 1.0_dp/SQRT((r + da + db)**2 + add) &
494 428462 : + 1.0_dp/SQRT((r - da + db)**2 + add) + 1.0_dp/SQRT((r + da - db)**2 + add)
495 : xyxy = &
496 : +1.0_dp/SQRT((r - da)**2 + add) + 1.0_dp/SQRT((r + da)**2 + add) &
497 : + 1.0_dp/SQRT((r - db)**2 + add) + 1.0_dp/SQRT((r + db)**2 + add) &
498 428462 : - 2.0_dp/SQRT(r**2 + add)
499 428462 : charg = (zzzz*0.0625_dp - xyxy*0.125_dp)*fact
500 428462 : RETURN
501 : END IF
502 : ! ZZ - XX
503 3427696 : IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzz .AND. (m2 == CLMxx .OR. m2 == CLMyy)) THEN
504 : zzzz = &
505 : -1.0_dp/SQRT((r + da)**2 + add) + 1.0_dp/SQRT((r + da)**2 + db**2 + add) &
506 856924 : - 1.0_dp/SQRT((r - da)**2 + add) + 1.0_dp/SQRT((r - da)**2 + db**2 + add)
507 : xyxy = &
508 856924 : +1.0_dp/SQRT(r**2 + db**2 + add) - 1.0_dp/SQRT(r**2 + add)
509 856924 : charg = (zzzz*0.125_dp - xyxy*0.25_dp)*fact
510 856924 : RETURN
511 : END IF
512 : ! X - ZX
513 2570772 : IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMp .AND. m2 == CLMzp) THEN
514 856924 : db = db/2.0_dp
515 : dxqxz = &
516 : -1.0_dp/SQRT((r - db)**2 + (da - db)**2 + add) + 1.0_dp/SQRT((r + db)**2 + (da - db)**2 + add) &
517 856924 : + 1.0_dp/SQRT((r - db)**2 + (da + db)**2 + add) - 1.0_dp/SQRT((r + db)**2 + (da + db)**2 + add)
518 856924 : charg = dxqxz*0.25_dp*fact
519 856924 : RETURN
520 : END IF
521 : ! ZX - ZX
522 1713848 : IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzp .AND. m2 == CLMzp) THEN
523 428462 : da = da/2.0_dp
524 428462 : db = db/2.0_dp
525 : qxzqxz = &
526 : +1.0_dp/SQRT((r + da - db)**2 + (da - db)**2 + add) - 1.0_dp/SQRT((r + da + db)**2 + (da - db)**2 + add) &
527 : - 1.0_dp/SQRT((r - da - db)**2 + (da - db)**2 + add) + 1.0_dp/SQRT((r - da + db)**2 + (da - db)**2 + add) &
528 : - 1.0_dp/SQRT((r + da - db)**2 + (da + db)**2 + add) + 1.0_dp/SQRT((r + da + db)**2 + (da + db)**2 + add) &
529 428462 : + 1.0_dp/SQRT((r - da - db)**2 + (da + db)**2 + add) - 1.0_dp/SQRT((r - da + db)**2 + (da + db)**2 + add)
530 428462 : charg = qxzqxz*0.125_dp*fact
531 428462 : RETURN
532 : END IF
533 : ! XX - XX
534 1285386 : IF (l1 == 2 .AND. l2 == 2 .AND. (((m1 == CLMyy) .AND. (m2 == CLMyy)) .OR. ((m1 == CLMxx) .AND. (m2 == CLMxx)))) THEN
535 : qxxqxx = &
536 : +2.0_dp/SQRT(r**2 + add) + 1.0_dp/SQRT(r**2 + (da - db)**2 + add) &
537 : + 1.0_dp/SQRT(r**2 + (da + db)**2 + add) - 2.0_dp/SQRT(r**2 + da**2 + add) &
538 428462 : - 2.0_dp/SQRT(r**2 + db**2 + add)
539 428462 : charg = qxxqxx*0.125_dp*fact
540 428462 : RETURN
541 : END IF
542 : ! XX - YY
543 856924 : IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMyy .AND. m2 == CLMxx) THEN
544 : qxxqyy = &
545 : +1.0_dp/SQRT(r**2 + add) - 1.0_dp/SQRT(r**2 + da**2 + add) &
546 428462 : - 1.0_dp/SQRT(r**2 + db**2 + add) + 1.0_dp/SQRT(r**2 + da**2 + db**2 + add)
547 428462 : charg = qxxqyy*0.25_dp*fact
548 428462 : RETURN
549 : END IF
550 : ! XY - XY
551 428462 : IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMxy .AND. m2 == CLMxy) THEN
552 : qxxqxx = &
553 : +2.0_dp/SQRT(r**2 + add) + 1.0_dp/SQRT(r**2 + (da - db)**2 + add) &
554 : + 1.0_dp/SQRT(r**2 + (da + db)**2 + add) - 2.0_dp/SQRT(r**2 + da**2 + add) &
555 428462 : - 2.0_dp/SQRT(r**2 + db**2 + add)
556 : qxxqyy = &
557 : +1.0_dp/SQRT(r**2 + add) - 1.0_dp/SQRT(r**2 + da**2 + add) &
558 428462 : - 1.0_dp/SQRT(r**2 + db**2 + add) + 1.0_dp/SQRT(r**2 + da**2 + db**2 + add)
559 428462 : charg = 0.5_dp*(qxxqxx*0.125_dp - qxxqyy*0.25_dp)*fact
560 428462 : RETURN
561 : END IF
562 : ! We should NEVER reach this point
563 0 : CPABORT("")
564 0 : END FUNCTION charg_int_nri
565 :
566 : ! **************************************************************************************************
567 : !> \brief Derivatives of interaction function between two point-charge
568 : !> configurations NDDO sp-code.
569 : !> Non-Rotational Invariant definition of quadrupoles
570 : !>
571 : !> r - Distance r12
572 : !> l1,m - Quantum numbers for multipole of configuration 1
573 : !> l2,m - Quantum numbers for multipole of configuration 2
574 : !> da - charge separation of configuration 1
575 : !> db - charge separation of configuration 2
576 : !> add - additive term
577 : !>
578 : !> \param r ...
579 : !> \param l1_i ...
580 : !> \param l2_i ...
581 : !> \param m1_i ...
582 : !> \param m2_i ...
583 : !> \param da_i ...
584 : !> \param db_i ...
585 : !> \param add0 ...
586 : !> \param fact_screen ...
587 : !> \return ...
588 : !> \date 04.2008 [tlaino]
589 : !> \author Teodoro Laino [tlaino] - University of Zurich
590 : ! **************************************************************************************************
591 53899489 : FUNCTION dcharg_int_nri(r, l1_i, l2_i, m1_i, m2_i, da_i, db_i, add0, fact_screen) RESULT(charg)
592 : REAL(KIND=dp), INTENT(in) :: r
593 : INTEGER, INTENT(in) :: l1_i, l2_i, m1_i, m2_i
594 : REAL(KIND=dp), INTENT(in) :: da_i, db_i, add0, fact_screen
595 : REAL(KIND=dp) :: charg
596 :
597 : INTEGER :: l1, l2, m1, m2
598 : REAL(KIND=dp) :: add, da, db, dxdx, dxqxz, dzdz, dzqxx, &
599 : dzqzz, fact, qqxx, qqzz, qxxqxx, &
600 : qxxqyy, qxzqxz, xyxy, zzzz
601 :
602 : ! Computing only Integral Derivatives
603 :
604 53899489 : IF (l1_i < l2_i) THEN
605 3385714 : l1 = l1_i
606 3385714 : l2 = l2_i
607 3385714 : m1 = m1_i
608 3385714 : m2 = m2_i
609 3385714 : da = da_i
610 3385714 : db = db_i
611 3385714 : fact = 1.0_dp
612 50513775 : ELSE IF (l1_i > l2_i) THEN
613 13898398 : l1 = l2_i
614 13898398 : l2 = l1_i
615 13898398 : m1 = m2_i
616 13898398 : m2 = m1_i
617 13898398 : da = db_i
618 13898398 : db = da_i
619 13898398 : fact = (-1.0_dp)**(l1 + l2)
620 36615377 : ELSE IF (l1_i == l2_i) THEN
621 36615377 : l1 = l1_i
622 36615377 : l2 = l2_i
623 36615377 : IF (m1_i <= m2_i) THEN
624 36422710 : m1 = m1_i
625 36422710 : m2 = m2_i
626 36422710 : da = da_i
627 36422710 : db = db_i
628 : ELSE
629 192667 : m1 = m2_i
630 192667 : m2 = m1_i
631 192667 : da = db_i
632 192667 : db = da_i
633 : END IF
634 : fact = 1.0_dp
635 : END IF
636 53899489 : charg = 0.0_dp
637 53899489 : add = add0*fact_screen
638 : ! Q - Q.
639 53899489 : IF (l1 == 0 .AND. l2 == 0) THEN
640 34881374 : charg = r/SQRT(r**2 + add)**3
641 34881374 : charg = -charg*fact
642 34881374 : RETURN
643 : END IF
644 : ! Q - Z.
645 19018115 : IF (l1 == 0 .AND. l2 == 1 .AND. m2 == CLMz) THEN
646 5247592 : charg = (r + db)/SQRT((r + db)**2 + add)**3 - (r - db)/SQRT((r - db)**2 + add)**3
647 5247592 : charg = -charg*0.5_dp*fact
648 5247592 : RETURN
649 : END IF
650 : ! Z - Z.
651 13770523 : IF (l1 == 1 .AND. l2 == 1 .AND. m1 == CLMz .AND. m2 == CLMz) THEN
652 : dzdz = &
653 : +(r + da - db)/SQRT((r + da - db)**2 + add)**3 + (r - da + db)/SQRT((r - da + db)**2 + add)**3 &
654 192667 : - (r - da - db)/SQRT((r - da - db)**2 + add)**3 - (r + da + db)/SQRT((r + da + db)**2 + add)**3
655 192667 : charg = -dzdz*0.25_dp*fact
656 192667 : RETURN
657 : END IF
658 : ! X - X
659 13577856 : IF (l1 == 1 .AND. l2 == 1 .AND. m1 == CLMp .AND. m2 == CLMp) THEN
660 192667 : dxdx = 2.0_dp*r/SQRT(r**2 + (da - db)**2 + add)**3 - 2.0_dp*r/SQRT(r**2 + (da + db)**2 + add)**3
661 192667 : charg = -dxdx*0.25_dp*fact
662 192667 : RETURN
663 : END IF
664 : ! Q - ZZ
665 13385189 : IF (l1 == 0 .AND. l2 == 2 .AND. m2 == CLMzz) THEN
666 5247592 : qqzz = (r - db)/SQRT((r - db)**2 + add)**3 - 2.0_dp*r/SQRT(r**2 + add)**3 + (r + db)/SQRT((r + db)**2 + add)**3
667 5247592 : charg = -qqzz*0.25_dp*fact
668 5247592 : RETURN
669 : END IF
670 : ! Q - XX
671 8137597 : IF (l1 == 0 .AND. l2 == 2 .AND. (m2 == CLMyy .OR. m2 == CLMxx)) THEN
672 5632926 : qqxx = -r/SQRT(r**2 + add)**3 + r/SQRT(r**2 + add + db**2)**3
673 5632926 : charg = -qqxx*0.5_dp*fact
674 5632926 : RETURN
675 : END IF
676 : ! Z - ZZ
677 2504671 : IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMz .AND. m2 == CLMzz) THEN
678 : dzqzz = &
679 : +(r - da - db)/SQRT((r - da - db)**2 + add)**3 - 2.0_dp*(r - da)/SQRT((r - da)**2 + add)**3 &
680 : + (r - da + db)/SQRT((r - da + db)**2 + add)**3 - (r + da - db)/SQRT((r + da - db)**2 + add)**3 &
681 385334 : + 2.0_dp*(r + da)/SQRT((r + da)**2 + add)**3 - (r + da + db)/SQRT((r + da + db)**2 + add)**3
682 385334 : charg = -dzqzz*0.125_dp*fact
683 385334 : RETURN
684 : END IF
685 : ! Z - XX
686 2119337 : IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMz .AND. (m2 == CLMyy .OR. m2 == CLMxx)) THEN
687 : dzqxx = &
688 : +(r + da)/SQRT((r + da)**2 + add)**3 - (r + da)/SQRT((r + da)**2 + add + db**2)**3 &
689 385334 : - (r - da)/SQRT((r - da)**2 + add)**3 + (r - da)/SQRT((r - da)**2 + add + db**2)**3
690 385334 : charg = -dzqxx*0.25_dp*fact
691 385334 : RETURN
692 : END IF
693 : ! ZZ - ZZ
694 1734003 : IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzz .AND. m2 == CLMzz) THEN
695 : zzzz = &
696 : +(r - da - db)/SQRT((r - da - db)**2 + add)**3 + (r + da + db)/SQRT((r + da + db)**2 + add)**3 &
697 192667 : + (r - da + db)/SQRT((r - da + db)**2 + add)**3 + (r + da - db)/SQRT((r + da - db)**2 + add)**3
698 : xyxy = &
699 : +(r - da)/SQRT((r - da)**2 + add)**3 + (r + da)/SQRT((r + da)**2 + add)**3 &
700 : + (r - db)/SQRT((r - db)**2 + add)**3 + (r + db)/SQRT((r + db)**2 + add)**3 &
701 192667 : - 2.0_dp*r/SQRT(r**2 + add)**3
702 192667 : charg = -(zzzz*0.0625_dp - xyxy*0.125_dp)*fact
703 192667 : RETURN
704 : END IF
705 : ! ZZ - XX
706 1541336 : IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzz .AND. (m2 == CLMxx .OR. m2 == CLMyy)) THEN
707 : zzzz = &
708 : -(r + da)/SQRT((r + da)**2 + add)**3 + (r + da)/SQRT((r + da)**2 + db**2 + add)**3 &
709 385334 : - (r - da)/SQRT((r - da)**2 + add)**3 + (r - da)/SQRT((r - da)**2 + db**2 + add)**3
710 385334 : xyxy = r/SQRT(r**2 + db**2 + add)**3 - r/SQRT(r**2 + add)**3
711 385334 : charg = -(zzzz*0.125_dp - xyxy*0.25_dp)*fact
712 385334 : RETURN
713 : END IF
714 : ! X - ZX
715 1156002 : IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMp .AND. m2 == CLMzp) THEN
716 385334 : db = db/2.0_dp
717 : dxqxz = &
718 : -(r - db)/SQRT((r - db)**2 + (da - db)**2 + add)**3 + (r + db)/SQRT((r + db)**2 + (da - db)**2 + add)**3 &
719 385334 : + (r - db)/SQRT((r - db)**2 + (da + db)**2 + add)**3 - (r + db)/SQRT((r + db)**2 + (da + db)**2 + add)**3
720 385334 : charg = -dxqxz*0.25_dp*fact
721 385334 : RETURN
722 : END IF
723 : ! ZX - ZX
724 770668 : IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzp .AND. m2 == CLMzp) THEN
725 192667 : da = da/2.0_dp
726 192667 : db = db/2.0_dp
727 : qxzqxz = &
728 : +(r + da - db)/SQRT((r + da - db)**2 + (da - db)**2 + add)**3 - (r + da + db)/SQRT((r + da + db)**2 + (da - db)**2 + add)**3 &
729 : - (r - da - db)/SQRT((r - da - db)**2 + (da - db)**2 + add)**3 + (r - da + db)/SQRT((r - da + db)**2 + (da - db)**2 + add)**3 &
730 : - (r + da - db)/SQRT((r + da - db)**2 + (da + db)**2 + add)**3 + (r + da + db)/SQRT((r + da + db)**2 + (da + db)**2 + add)**3 &
731 192667 : + (r - da - db)/SQRT((r - da - db)**2 + (da + db)**2 + add)**3 - (r - da + db)/SQRT((r - da + db)**2 + (da + db)**2 + add)**3
732 192667 : charg = -qxzqxz*0.125_dp*fact
733 192667 : RETURN
734 : END IF
735 : ! XX - XX
736 578001 : IF (l1 == 2 .AND. l2 == 2 .AND. (((m1 == CLMyy) .AND. (m2 == CLMyy)) .OR. ((m1 == CLMxx) .AND. (m2 == CLMxx)))) THEN
737 : qxxqxx = &
738 : +2.0_dp*r/SQRT(r**2 + add)**3 + r/SQRT(r**2 + (da - db)**2 + add)**3 &
739 : + r/SQRT(r**2 + (da + db)**2 + add)**3 - 2.0_dp*r/SQRT(r**2 + da**2 + add)**3 &
740 192667 : - 2.0_dp*r/SQRT(r**2 + db**2 + add)**3
741 192667 : charg = -qxxqxx*0.125_dp*fact
742 192667 : RETURN
743 : END IF
744 : ! XX - YY
745 385334 : IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMyy .AND. m2 == CLMxx) THEN
746 : qxxqyy = &
747 : +r/SQRT(r**2 + add)**3 - r/SQRT(r**2 + da**2 + add)**3 &
748 192667 : - r/SQRT(r**2 + db**2 + add)**3 + r/SQRT(r**2 + da**2 + db**2 + add)**3
749 192667 : charg = -qxxqyy*0.25_dp*fact
750 192667 : RETURN
751 : END IF
752 : ! XY - XY
753 192667 : IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMxy .AND. m2 == CLMxy) THEN
754 : qxxqxx = &
755 : +2.0_dp*r/SQRT(r**2 + add)**3 + r/SQRT(r**2 + (da - db)**2 + add)**3 &
756 : + r/SQRT(r**2 + (da + db)**2 + add)**3 - 2.0_dp*r/SQRT(r**2 + da**2 + add)**3 &
757 192667 : - 2.0_dp*r/SQRT(r**2 + db**2 + add)**3
758 : qxxqyy = &
759 : +r/SQRT(r**2 + add)**3 - r/SQRT(r**2 + da**2 + add)**3 &
760 192667 : - r/SQRT(r**2 + db**2 + add)**3 + r/SQRT(r**2 + da**2 + db**2 + add)**3
761 192667 : charg = -0.5_dp*(qxxqxx*0.125_dp - qxxqyy*0.25_dp)*fact
762 192667 : RETURN
763 : END IF
764 : ! We should NEVER reach this point
765 0 : CPABORT("")
766 0 : END FUNCTION dcharg_int_nri
767 :
768 : ! **************************************************************************************************
769 : !> \brief Derivatives of interaction function between two point-charge
770 : !> configurations NDDO sp-code. The derivative takes care of the screening
771 : !> term only.
772 : !> Non-Rotational Invariant definition of quadrupoles
773 : !>
774 : !> r - Distance r12
775 : !> l1,m - Quantum numbers for multipole of configuration 1
776 : !> l2,m - Quantum numbers for multipole of configuration 2
777 : !> da - charge separation of configuration 1
778 : !> db - charge separation of configuration 2
779 : !> add - additive term
780 : !>
781 : !> \param r ...
782 : !> \param l1_i ...
783 : !> \param l2_i ...
784 : !> \param m1_i ...
785 : !> \param m2_i ...
786 : !> \param da_i ...
787 : !> \param db_i ...
788 : !> \param add0 ...
789 : !> \param fact_screen ...
790 : !> \return ...
791 : !> \date 04.2008 [tlaino]
792 : !> \author Teodoro Laino [tlaino] - University of Zurich
793 : ! **************************************************************************************************
794 2803761 : FUNCTION dcharg_int_nri_fs(r, l1_i, l2_i, m1_i, m2_i, da_i, db_i, add0, fact_screen) RESULT(charg)
795 : REAL(KIND=dp), INTENT(in) :: r
796 : INTEGER, INTENT(in) :: l1_i, l2_i, m1_i, m2_i
797 : REAL(KIND=dp), INTENT(in) :: da_i, db_i, add0, fact_screen
798 : REAL(KIND=dp) :: charg
799 :
800 : INTEGER :: l1, l2, m1, m2
801 : REAL(KIND=dp) :: add, da, db, dxdx, dxqxz, dzdz, dzqxx, &
802 : dzqzz, fact, qqxx, qqzz, qxxqxx, &
803 : qxxqyy, qxzqxz, xyxy, zzzz
804 :
805 : ! Computing only Integral Derivatives
806 :
807 2803761 : IF (l1_i < l2_i) THEN
808 653103 : l1 = l1_i
809 653103 : l2 = l2_i
810 653103 : m1 = m1_i
811 653103 : m2 = m2_i
812 653103 : da = da_i
813 653103 : db = db_i
814 653103 : fact = 1.0_dp
815 2150658 : ELSE IF (l1_i > l2_i) THEN
816 732123 : l1 = l2_i
817 732123 : l2 = l1_i
818 732123 : m1 = m2_i
819 732123 : m2 = m1_i
820 732123 : da = db_i
821 732123 : db = da_i
822 732123 : fact = (-1.0_dp)**(l1 + l2)
823 1418535 : ELSE IF (l1_i == l2_i) THEN
824 1418535 : l1 = l1_i
825 1418535 : l2 = l2_i
826 1418535 : IF (m1_i <= m2_i) THEN
827 1380180 : m1 = m1_i
828 1380180 : m2 = m2_i
829 1380180 : da = da_i
830 1380180 : db = db_i
831 : ELSE
832 38355 : m1 = m2_i
833 38355 : m2 = m1_i
834 38355 : da = db_i
835 38355 : db = da_i
836 : END IF
837 : fact = 1.0_dp
838 : END IF
839 2803761 : charg = 0.0_dp
840 2803761 : add = add0*fact_screen
841 : ! The 0.5 factor handles the derivative of the SQRT
842 2803761 : fact = fact*0.5_dp
843 : ! Q - Q.
844 2803761 : IF (l1 == 0 .AND. l2 == 0) THEN
845 1073340 : charg = add0/SQRT(r**2 + add)**3
846 1073340 : charg = -charg*fact
847 1073340 : RETURN
848 : END IF
849 : ! Q - Z.
850 1730421 : IF (l1 == 0 .AND. l2 == 1 .AND. m2 == CLMz) THEN
851 359462 : charg = add0/SQRT((r + db)**2 + add)**3 - add0/SQRT((r - db)**2 + add)**3
852 359462 : charg = -charg*0.5_dp*fact
853 359462 : RETURN
854 : END IF
855 : ! Z - Z.
856 1370959 : IF (l1 == 1 .AND. l2 == 1 .AND. m1 == CLMz .AND. m2 == CLMz) THEN
857 : dzdz = &
858 : +add0/SQRT((r + da - db)**2 + add)**3 + add0/SQRT((r - da + db)**2 + add)**3 &
859 38355 : - add0/SQRT((r - da - db)**2 + add)**3 - add0/SQRT((r + da + db)**2 + add)**3
860 38355 : charg = -dzdz*0.25_dp*fact
861 38355 : RETURN
862 : END IF
863 : ! X - X
864 1332604 : IF (l1 == 1 .AND. l2 == 1 .AND. m1 == CLMp .AND. m2 == CLMp) THEN
865 38355 : dxdx = 2.0_dp*add0/SQRT(r**2 + (da - db)**2 + add)**3 - 2.0_dp*add0/SQRT(r**2 + (da + db)**2 + add)**3
866 38355 : charg = -dxdx*0.25_dp*fact
867 38355 : RETURN
868 : END IF
869 : ! Q - ZZ
870 1294249 : IF (l1 == 0 .AND. l2 == 2 .AND. m2 == CLMzz) THEN
871 359462 : qqzz = add0/SQRT((r - db)**2 + add)**3 - 2.0_dp*add0/SQRT(r**2 + add)**3 + add0/SQRT((r + db)**2 + add)**3
872 359462 : charg = -qqzz*0.25_dp*fact
873 359462 : RETURN
874 : END IF
875 : ! Q - XX
876 934787 : IF (l1 == 0 .AND. l2 == 2 .AND. (m2 == CLMyy .OR. m2 == CLMxx)) THEN
877 436172 : qqxx = -add0/SQRT(r**2 + add)**3 + add0/SQRT(r**2 + add + db**2)**3
878 436172 : charg = -qqxx*0.5_dp*fact
879 436172 : RETURN
880 : END IF
881 : ! Z - ZZ
882 498615 : IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMz .AND. m2 == CLMzz) THEN
883 : dzqzz = &
884 : +add0/SQRT((r - da - db)**2 + add)**3 - 2.0_dp*add0/SQRT((r - da)**2 + add)**3 &
885 : + add0/SQRT((r - da + db)**2 + add)**3 - add0/SQRT((r + da - db)**2 + add)**3 &
886 76710 : + 2.0_dp*add0/SQRT((r + da)**2 + add)**3 - add0/SQRT((r + da + db)**2 + add)**3
887 76710 : charg = -dzqzz*0.125_dp*fact
888 76710 : RETURN
889 : END IF
890 : ! Z - XX
891 421905 : IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMz .AND. (m2 == CLMyy .OR. m2 == CLMxx)) THEN
892 : dzqxx = &
893 : +add0/SQRT((r + da)**2 + add)**3 - add0/SQRT((r + da)**2 + add + db**2)**3 &
894 76710 : - add0/SQRT((r - da)**2 + add)**3 + add0/SQRT((r - da)**2 + add + db**2)**3
895 76710 : charg = -dzqxx*0.25_dp*fact
896 76710 : RETURN
897 : END IF
898 : ! ZZ - ZZ
899 345195 : IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzz .AND. m2 == CLMzz) THEN
900 : zzzz = &
901 : +add0/SQRT((r - da - db)**2 + add)**3 + add0/SQRT((r + da + db)**2 + add)**3 &
902 38355 : + add0/SQRT((r - da + db)**2 + add)**3 + add0/SQRT((r + da - db)**2 + add)**3
903 : xyxy = &
904 : +add0/SQRT((r - da)**2 + add)**3 + add0/SQRT((r + da)**2 + add)**3 &
905 : + add0/SQRT((r - db)**2 + add)**3 + add0/SQRT((r + db)**2 + add)**3 &
906 38355 : - 2.0_dp*add0/SQRT(r**2 + add)**3
907 38355 : charg = -(zzzz*0.0625_dp - xyxy*0.125_dp)*fact
908 38355 : RETURN
909 : END IF
910 : ! ZZ - XX
911 306840 : IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzz .AND. (m2 == CLMxx .OR. m2 == CLMyy)) THEN
912 : zzzz = &
913 : -add0/SQRT((r + da)**2 + add)**3 + add0/SQRT((r + da)**2 + db**2 + add)**3 &
914 76710 : - add0/SQRT((r - da)**2 + add)**3 + add0/SQRT((r - da)**2 + db**2 + add)**3
915 76710 : xyxy = add0/SQRT(r**2 + db**2 + add)**3 - add0/SQRT(r**2 + add)**3
916 76710 : charg = -(zzzz*0.125_dp - xyxy*0.25_dp)*fact
917 76710 : RETURN
918 : END IF
919 : ! X - ZX
920 230130 : IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMp .AND. m2 == CLMzp) THEN
921 76710 : db = db/2.0_dp
922 : dxqxz = &
923 : -add0/SQRT((r - db)**2 + (da - db)**2 + add)**3 + add0/SQRT((r + db)**2 + (da - db)**2 + add)**3 &
924 76710 : + add0/SQRT((r - db)**2 + (da + db)**2 + add)**3 - add0/SQRT((r + db)**2 + (da + db)**2 + add)**3
925 76710 : charg = -dxqxz*0.25_dp*fact
926 76710 : RETURN
927 : END IF
928 : ! ZX - ZX
929 153420 : IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzp .AND. m2 == CLMzp) THEN
930 38355 : da = da/2.0_dp
931 38355 : db = db/2.0_dp
932 : qxzqxz = &
933 : +add0/SQRT((r + da - db)**2 + (da - db)**2 + add)**3 - add0/SQRT((r + da + db)**2 + (da - db)**2 + add)**3 &
934 : - add0/SQRT((r - da - db)**2 + (da - db)**2 + add)**3 + add0/SQRT((r - da + db)**2 + (da - db)**2 + add)**3 &
935 : - add0/SQRT((r + da - db)**2 + (da + db)**2 + add)**3 + add0/SQRT((r + da + db)**2 + (da + db)**2 + add)**3 &
936 38355 : + add0/SQRT((r - da - db)**2 + (da + db)**2 + add)**3 - add0/SQRT((r - da + db)**2 + (da + db)**2 + add)**3
937 38355 : charg = -qxzqxz*0.125_dp*fact
938 38355 : RETURN
939 : END IF
940 : ! XX - XX
941 115065 : IF (l1 == 2 .AND. l2 == 2 .AND. (((m1 == CLMyy) .AND. (m2 == CLMyy)) .OR. ((m1 == CLMxx) .AND. (m2 == CLMxx)))) THEN
942 : qxxqxx = &
943 : +2.0_dp*add0/SQRT(r**2 + add)**3 + add0/SQRT(r**2 + (da - db)**2 + add)**3 &
944 : + add0/SQRT(r**2 + (da + db)**2 + add)**3 - 2.0_dp*add0/SQRT(r**2 + da**2 + add)**3 &
945 38355 : - 2.0_dp*add0/SQRT(r**2 + db**2 + add)**3
946 38355 : charg = -qxxqxx*0.125_dp*fact
947 38355 : RETURN
948 : END IF
949 : ! XX - YY
950 76710 : IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMyy .AND. m2 == CLMxx) THEN
951 : qxxqyy = &
952 : +add0/SQRT(r**2 + add)**3 - add0/SQRT(r**2 + da**2 + add)**3 &
953 38355 : - add0/SQRT(r**2 + db**2 + add)**3 + add0/SQRT(r**2 + da**2 + db**2 + add)**3
954 38355 : charg = -qxxqyy*0.25_dp*fact
955 38355 : RETURN
956 : END IF
957 : ! XY - XY
958 38355 : IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMxy .AND. m2 == CLMxy) THEN
959 : qxxqxx = &
960 : +2.0_dp*add0/SQRT(r**2 + add)**3 + add0/SQRT(r**2 + (da - db)**2 + add)**3 &
961 : + add0/SQRT(r**2 + (da + db)**2 + add)**3 - 2.0_dp*add0/SQRT(r**2 + da**2 + add)**3 &
962 38355 : - 2.0_dp*add0/SQRT(r**2 + db**2 + add)**3
963 : qxxqyy = &
964 : +add0/SQRT(r**2 + add)**3 - add0/SQRT(r**2 + da**2 + add)**3 &
965 38355 : - add0/SQRT(r**2 + db**2 + add)**3 + add0/SQRT(r**2 + da**2 + db**2 + add)**3
966 38355 : charg = -0.5_dp*(qxxqxx*0.125_dp - qxxqyy*0.25_dp)*fact
967 38355 : RETURN
968 : END IF
969 : ! We should NEVER reach this point
970 0 : CPABORT("")
971 0 : END FUNCTION dcharg_int_nri_fs
972 :
973 : ! **************************************************************************************************
974 : !> \brief General driver for computing semi-empirical integrals <ij|kl>
975 : !> involving d-orbitals.
976 : !> The choice of the linear quadrupole was REALLY unhappy
977 : !> in the first development of the NDDO codes. That choice makes
978 : !> impossible the merging of the integral code with or without d-orbitals
979 : !> unless a reparametrization of all NDDO codes for s and p orbitals will
980 : !> be performed.. more over the choice of the linear quadrupole does not make
981 : !> calculations rotational invariants (of course the rotational invariant
982 : !> can be forced). The definitions of quadrupoles for d-orbitals is the
983 : !> correct one in order to have the rotational invariant property by
984 : !> construction..
985 : !>
986 : !> \param sepi ...
987 : !> \param sepj ...
988 : !> \param ij ...
989 : !> \param kl ...
990 : !> \param li ...
991 : !> \param lj ...
992 : !> \param lk ...
993 : !> \param ll ...
994 : !> \param ic ...
995 : !> \param r ...
996 : !> \param se_int_control ...
997 : !> \param se_int_screen ...
998 : !> \param itype ...
999 : !> \return ...
1000 : !> \date 03.2008 [tlaino]
1001 : !> \author Teodoro Laino [tlaino] - University of Zurich
1002 : ! **************************************************************************************************
1003 2873682 : FUNCTION ijkl_d(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_control, &
1004 : se_int_screen, itype) RESULT(res)
1005 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
1006 : INTEGER, INTENT(IN) :: ij, kl, li, lj, lk, ll, ic
1007 : REAL(KIND=dp), INTENT(IN) :: r
1008 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
1009 : TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen
1010 : INTEGER, INTENT(IN) :: itype
1011 : REAL(KIND=dp) :: res
1012 :
1013 : res = ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, &
1014 : se_int_control%integral_screening, se_int_control%shortrange, &
1015 : se_int_control%pc_coulomb_int, se_int_control%max_multipole, &
1016 2873682 : itype, charg_int_ri)
1017 :
1018 : ! If only the shortrange component is requested we can skip the rest
1019 2873682 : IF ((.NOT. se_int_control%pc_coulomb_int) .AND. (itype /= do_method_pchg)) THEN
1020 : ! Handle the 1/r^3 term, this term is ALWAYS false for KDSO-D integrals
1021 2674266 : IF (se_int_control%shortrange .AND. se_int_control%do_ewald_r3) THEN
1022 : res = res - ijkl_low_3(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, &
1023 0 : itype, charg_int_3)
1024 : END IF
1025 : END IF
1026 2873682 : END FUNCTION ijkl_d
1027 :
1028 : ! **************************************************************************************************
1029 : !> \brief General driver for computing the derivatives of semi-empirical integrals <ij|kl>
1030 : !> involving d-orbitals.
1031 : !> The choice of the linear quadrupole was REALLY unhappy
1032 : !> in the first development of the NDDO codes. That choice makes
1033 : !> impossible the merging of the integral code with or without d-orbitals
1034 : !> unless a reparametrization of all NDDO codes for s and p orbitals will
1035 : !> be performed.. more over the choice of the linear quadrupole does not make
1036 : !> calculations rotational invariants (of course the rotational invariant
1037 : !> can be forced). The definitions of quadrupoles for d-orbitals is the
1038 : !> correct one in order to have the rotational invariant property by
1039 : !> construction..
1040 : !>
1041 : !> \param sepi ...
1042 : !> \param sepj ...
1043 : !> \param ij ...
1044 : !> \param kl ...
1045 : !> \param li ...
1046 : !> \param lj ...
1047 : !> \param lk ...
1048 : !> \param ll ...
1049 : !> \param ic ...
1050 : !> \param r ...
1051 : !> \param se_int_control ...
1052 : !> \param se_int_screen ...
1053 : !> \param itype ...
1054 : !> \return ...
1055 : !> \date 03.2008 [tlaino]
1056 : !> \author Teodoro Laino [tlaino] - University of Zurich
1057 : ! **************************************************************************************************
1058 1373222 : FUNCTION d_ijkl_d(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_control, &
1059 : se_int_screen, itype) RESULT(res)
1060 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
1061 : INTEGER, INTENT(IN) :: ij, kl, li, lj, lk, ll, ic
1062 : REAL(KIND=dp), INTENT(IN) :: r
1063 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
1064 : TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen
1065 : INTEGER, INTENT(IN) :: itype
1066 : REAL(KIND=dp) :: res
1067 :
1068 : REAL(KIND=dp) :: dfs, srd
1069 :
1070 1373222 : IF (se_int_control%integral_screening == do_se_IS_kdso_d) THEN
1071 : res = ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, &
1072 : se_int_control%integral_screening, .FALSE., &
1073 : se_int_control%pc_coulomb_int, se_int_control%max_multipole, &
1074 810708 : itype, dcharg_int_ri)
1075 :
1076 810708 : IF (.NOT. se_int_control%pc_coulomb_int) THEN
1077 : dfs = ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, &
1078 : se_int_control%integral_screening, .FALSE., .FALSE., &
1079 711000 : se_int_control%max_multipole, itype, dcharg_int_ri_fs)
1080 711000 : res = res + dfs*se_int_screen%dft
1081 :
1082 : ! In case we need the shortrange part we have to evaluate an additional derivative
1083 : ! to handle the derivative of the Tapering term
1084 711000 : IF (se_int_control%shortrange) THEN
1085 : srd = ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, &
1086 : se_int_control%integral_screening, .FALSE., .TRUE., &
1087 711000 : se_int_control%max_multipole, itype, dcharg_int_ri)
1088 711000 : res = res - srd
1089 : END IF
1090 : END IF
1091 : ELSE
1092 : res = ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, &
1093 : se_int_control%integral_screening, se_int_control%shortrange, &
1094 : se_int_control%pc_coulomb_int, se_int_control%max_multipole, &
1095 562514 : itype, dcharg_int_ri)
1096 : END IF
1097 :
1098 : ! If only the shortrange component is requested we can skip the rest
1099 1373222 : IF ((.NOT. se_int_control%pc_coulomb_int) .AND. (itype /= do_method_pchg)) THEN
1100 : ! Handle the 1/r^3 term, this term is ALWAYS false for KDSO-D integrals
1101 1273514 : IF (se_int_control%shortrange .AND. se_int_control%do_ewald_r3) THEN
1102 : res = res - ijkl_low_3(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, &
1103 0 : itype, dcharg_int_3)
1104 : END IF
1105 : END IF
1106 :
1107 1373222 : END FUNCTION d_ijkl_d
1108 :
1109 : ! **************************************************************************************************
1110 : !> \brief Low level general driver for computing semi-empirical integrals <ij|kl>
1111 : !> and their derivatives involving d-orbitals.
1112 : !> The choice of the linear quadrupole was REALLY unhappy
1113 : !> in the first development of the NDDO codes. That choice makes
1114 : !> impossible the merging of the integral code with or without d-orbitals
1115 : !> unless a reparametrization of all NDDO codes for s and p orbitals will
1116 : !> be performed.. more over the choice of the linear quadrupole does not make
1117 : !> calculations rotational invariants (of course the rotational invariant
1118 : !> can be forced). The definitions of quadrupoles for d-orbitals is the
1119 : !> correct one in order to have the rotational invariant property by
1120 : !> construction..
1121 : !>
1122 : !> \param sepi ...
1123 : !> \param sepj ...
1124 : !> \param ij ...
1125 : !> \param kl ...
1126 : !> \param li ...
1127 : !> \param lj ...
1128 : !> \param lk ...
1129 : !> \param ll ...
1130 : !> \param ic ...
1131 : !> \param r ...
1132 : !> \param se_int_screen ...
1133 : !> \param iscreen ...
1134 : !> \param shortrange ...
1135 : !> \param pc_coulomb_int ...
1136 : !> \param max_multipole ...
1137 : !> \param itype ...
1138 : !> \param eval a function without explicit interface
1139 : !> \return ...
1140 : !> \date 03.2008 [tlaino]
1141 : !> \author Teodoro Laino [tlaino] - University of Zurich
1142 : ! **************************************************************************************************
1143 5668904 : FUNCTION ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, &
1144 : iscreen, shortrange, pc_coulomb_int, max_multipole, itype, eval) RESULT(res)
1145 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
1146 : INTEGER, INTENT(IN) :: ij, kl, li, lj, lk, ll, ic
1147 : REAL(KIND=dp), INTENT(IN) :: r
1148 : TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen
1149 : INTEGER, INTENT(IN) :: iscreen
1150 : LOGICAL, INTENT(IN) :: shortrange, pc_coulomb_int
1151 : INTEGER, INTENT(IN) :: max_multipole, itype
1152 :
1153 : PROCEDURE(eval_func_d) :: eval
1154 : REAL(KIND=dp) :: res
1155 :
1156 : INTEGER :: l1, l1max, l1min, l2, l2max, l2min, lij, &
1157 : lkl, lmin, m, mm
1158 : REAL(KIND=dp) :: add, ccc, chrg, dij, dkl, fact_screen, &
1159 : pij, pkl, s1, sum
1160 :
1161 5668904 : l1min = ABS(li - lj)
1162 5668904 : l1max = li + lj
1163 5668904 : lij = indexb(li + 1, lj + 1)
1164 5668904 : l2min = ABS(lk - ll)
1165 5668904 : l2max = lk + ll
1166 5668904 : lkl = indexb(lk + 1, ll + 1)
1167 :
1168 5668904 : l1max = MIN(l1max, 2)
1169 5668904 : l1min = MIN(l1min, 2)
1170 5668904 : l2max = MIN(l2max, 2)
1171 5668904 : l2min = MIN(l2min, 2)
1172 5668904 : sum = 0.0_dp
1173 5668904 : dij = 0.0_dp
1174 5668904 : dkl = 0.0_dp
1175 5668904 : fact_screen = 1.0_dp
1176 5668904 : IF (.NOT. pc_coulomb_int) THEN
1177 4658780 : IF (iscreen == do_se_IS_kdso_d) fact_screen = se_int_screen%ft
1178 : ! Standard value of the integral
1179 14783386 : DO l1 = l1min, l1max
1180 10124606 : IF (l1 == 0) THEN
1181 2651891 : IF (lij == 1) THEN
1182 466356 : pij = sepi%ko(1)
1183 466356 : IF (ic == 1) THEN
1184 205998 : pij = sepi%ko(9)
1185 : END IF
1186 2185535 : ELSE IF (lij == 3) THEN
1187 866784 : pij = sepi%ko(7)
1188 1318751 : ELSE IF (lij == 6) THEN
1189 1318751 : pij = sepi%ko(8)
1190 : END IF
1191 : ELSE
1192 7472715 : dij = sepi%cs(lij)
1193 7472715 : pij = sepi%ko(lij)
1194 : END IF
1195 : !
1196 36584358 : DO l2 = l2min, l2max
1197 21800972 : IF (l2 == 0) THEN
1198 5970952 : IF (lkl == 1) THEN
1199 1214094 : pkl = sepj%ko(1)
1200 1214094 : IF (ic == 2) THEN
1201 493346 : pkl = sepj%ko(9)
1202 : END IF
1203 4756858 : ELSE IF (lkl == 3) THEN
1204 2141370 : pkl = sepj%ko(7)
1205 2615488 : ELSE IF (lkl == 6) THEN
1206 2615488 : pkl = sepj%ko(8)
1207 : END IF
1208 : ELSE
1209 15830020 : dkl = sepj%cs(lkl)
1210 15830020 : pkl = sepj%ko(lkl)
1211 : END IF
1212 21800972 : IF (itype == do_method_pchg) THEN
1213 0 : add = 0.0_dp
1214 : ELSE
1215 21800972 : add = (pij + pkl)**2
1216 : END IF
1217 21800972 : lmin = MIN(l1, l2)
1218 21800972 : s1 = 0.0_dp
1219 72210402 : DO m = -lmin, lmin
1220 50409430 : ccc = clm_d(ij, l1, m)*clm_d(kl, l2, m)
1221 72210402 : IF (ABS(ccc) > EPSILON(0.0_dp)) THEN
1222 8019068 : mm = ABS(m)
1223 8019068 : chrg = eval(r, l1, l2, mm, dij, dkl, add, fact_screen)
1224 8019068 : s1 = s1 + chrg*ccc
1225 : END IF
1226 : END DO
1227 31925578 : sum = sum + s1
1228 : END DO
1229 : END DO
1230 : res = sum
1231 : END IF
1232 : ! Shortrange: Possibly computes pure Coulomb and subtract from the original integral valeu
1233 5668904 : IF (shortrange .OR. pc_coulomb_int) THEN
1234 2432124 : sum = 0.0_dp
1235 2432124 : dij = 0.0_dp
1236 2432124 : dkl = 0.0_dp
1237 2432124 : add = 0.0_dp
1238 2432124 : fact_screen = 0.0_dp
1239 7751484 : DO l1 = l1min, l1max
1240 5319360 : IF (l1 > max_multipole) CYCLE
1241 5319360 : IF (l1 /= 0) THEN
1242 3961368 : dij = sepi%cs(lij)
1243 : END IF
1244 : !
1245 19373688 : DO l2 = l2min, l2max
1246 11622204 : IF (l2 > max_multipole) CYCLE
1247 11622204 : IF (l2 /= 0) THEN
1248 8607240 : dkl = sepj%cs(lkl)
1249 : END IF
1250 11622204 : lmin = MIN(l1, l2)
1251 11622204 : s1 = 0.0_dp
1252 39100932 : DO m = -lmin, lmin
1253 27478728 : ccc = clm_d(ij, l1, m)*clm_d(kl, l2, m)
1254 39100932 : IF (ABS(ccc) > EPSILON(0.0_dp)) THEN
1255 4133502 : mm = ABS(m)
1256 4133502 : chrg = eval(r, l1, l2, mm, dij, dkl, add, fact_screen)
1257 4133502 : s1 = s1 + chrg*ccc
1258 : END IF
1259 : END DO
1260 16941564 : sum = sum + s1
1261 : END DO
1262 : END DO
1263 2432124 : IF (pc_coulomb_int) res = sum
1264 2432124 : IF (shortrange) res = res - sum
1265 : END IF
1266 5668904 : END FUNCTION ijkl_d_low
1267 :
1268 : ! **************************************************************************************************
1269 : !> \brief Interaction function between two point-charge configurations (MNDO/d)
1270 : !> Rotational invariant property built-in in the quadrupole definition
1271 : !> r - Distance r12
1272 : !> l1,m - Quantum numbers for multipole of configuration 1
1273 : !> l2,m - Quantum numbers for multipole of configuration 2
1274 : !> da - charge separation of configuration 1
1275 : !> db - charge separation of configuration 2
1276 : !> add - additive term
1277 : !>
1278 : !> \param r ...
1279 : !> \param l1_i ...
1280 : !> \param l2_i ...
1281 : !> \param m ...
1282 : !> \param da_i ...
1283 : !> \param db_i ...
1284 : !> \param add0 ...
1285 : !> \param fact_screen ...
1286 : !> \return ...
1287 : !> \date 03.2008 [tlaino]
1288 : !> \author Teodoro Laino [tlaino] - University of Zurich
1289 : ! **************************************************************************************************
1290 7367875 : FUNCTION charg_int_ri(r, l1_i, l2_i, m, da_i, db_i, add0, fact_screen) RESULT(charg)
1291 : REAL(KIND=dp), INTENT(in) :: r
1292 : INTEGER, INTENT(in) :: l1_i, l2_i, m
1293 : REAL(KIND=dp), INTENT(in) :: da_i, db_i, add0, fact_screen
1294 : REAL(KIND=dp) :: charg
1295 :
1296 : INTEGER :: l1, l2
1297 : REAL(KIND=dp) :: aa, ab, add, da, db, dxdx, dxqxz, dzdz, &
1298 : dzqzz, fact, qqzz, qxzqxz, xyxy, zzzz
1299 :
1300 7367875 : IF (l1_i <= l2_i) THEN
1301 5187661 : l1 = l1_i
1302 5187661 : l2 = l2_i
1303 5187661 : da = da_i
1304 5187661 : db = db_i
1305 5187661 : fact = 1.0_dp
1306 : ELSE IF (l1_i > l2_i) THEN
1307 2180214 : l1 = l2_i
1308 2180214 : l2 = l1_i
1309 2180214 : da = db_i
1310 2180214 : db = da_i
1311 2180214 : fact = (-1.0_dp)**(l1 + l2)
1312 : END IF
1313 7367875 : charg = 0.0_dp
1314 7367875 : add = add0*fact_screen
1315 : ! Q - Q.
1316 7367875 : IF (l1 == 0 .AND. l2 == 0) THEN
1317 1020785 : charg = fact/SQRT(r**2 + add)
1318 1020785 : RETURN
1319 : END IF
1320 : ! Q - Z.
1321 6347090 : IF (l1 == 0 .AND. l2 == 1) THEN
1322 972950 : charg = 1.0_dp/SQRT((r + db)**2 + add) - 1.0_dp/SQRT((r - db)**2 + add)
1323 972950 : charg = charg*0.5_dp*fact
1324 972950 : RETURN
1325 : END IF
1326 : ! Z - Z.
1327 5374140 : IF (l1 == 1 .AND. l2 == 1 .AND. m == 0) THEN
1328 : dzdz = &
1329 : +1.0_dp/SQRT((r + da - db)**2 + add) + 1.0_dp/SQRT((r - da + db)**2 + add) &
1330 185584 : - 1.0_dp/SQRT((r - da - db)**2 + add) - 1.0_dp/SQRT((r + da + db)**2 + add)
1331 185584 : charg = dzdz*0.25_dp*fact
1332 185584 : RETURN
1333 : END IF
1334 : ! X - X
1335 5188556 : IF (l1 == 1 .AND. l2 == 1 .AND. m == 1) THEN
1336 303696 : dxdx = 2.0_dp/SQRT(r**2 + (da - db)**2 + add) - 2.0_dp/SQRT(r**2 + (da + db)**2 + add)
1337 303696 : charg = dxdx*0.25_dp*fact
1338 303696 : RETURN
1339 : END IF
1340 : ! Q - ZZ
1341 4884860 : IF (l1 == 0 .AND. l2 == 2) THEN
1342 1945900 : qqzz = 1.0_dp/SQRT((r - db)**2 + add) - 2.0_dp/SQRT(r**2 + db**2 + add) + 1.0_dp/SQRT((r + db)**2 + add)
1343 1945900 : charg = qqzz*0.25_dp*fact
1344 1945900 : RETURN
1345 : END IF
1346 : ! Z - ZZ
1347 2938960 : IF (l1 == 1 .AND. l2 == 2 .AND. m == 0) THEN
1348 : dzqzz = &
1349 : +1.0_dp/SQRT((r - da - db)**2 + add) - 2.0_dp/SQRT((r - da)**2 + db**2 + add) &
1350 : + 1.0_dp/SQRT((r + db - da)**2 + add) - 1.0_dp/SQRT((r - db + da)**2 + add) &
1351 789552 : + 2.0_dp/SQRT((r + da)**2 + db**2 + add) - 1.0_dp/SQRT((r + da + db)**2 + add)
1352 789552 : charg = dzqzz*0.125_dp*fact
1353 789552 : RETURN
1354 : END IF
1355 : ! ZZ - ZZ
1356 2149408 : IF (l1 == 2 .AND. l2 == 2 .AND. m == 0) THEN
1357 : zzzz = &
1358 : +1.0_dp/SQRT((r - da - db)**2 + add) + 1.0_dp/SQRT((r + da + db)**2 + add) &
1359 : + 1.0_dp/SQRT((r - da + db)**2 + add) + 1.0_dp/SQRT((r + da - db)**2 + add) &
1360 : - 2.0_dp/SQRT((r - da)**2 + db**2 + add) - 2.0_dp/SQRT((r - db)**2 + da**2 + add) &
1361 : - 2.0_dp/SQRT((r + da)**2 + db**2 + add) - 2.0_dp/SQRT((r + db)**2 + da**2 + add) &
1362 789552 : + 2.0_dp/SQRT(r**2 + (da - db)**2 + add) + 2.0_dp/SQRT(r**2 + (da + db)**2 + add)
1363 : xyxy = &
1364 : +4.0_dp/SQRT(r**2 + (da - db)**2 + add) + 4.0_dp/SQRT(r**2 + (da + db)**2 + add) &
1365 789552 : - 8.0_dp/SQRT(r**2 + da**2 + db**2 + add)
1366 789552 : charg = (zzzz*0.0625_dp - xyxy*0.015625_dp)*fact
1367 789552 : RETURN
1368 : END IF
1369 : ! X - ZX
1370 1359856 : IF (l1 == 1 .AND. l2 == 2 .AND. m == 1) THEN
1371 607392 : ab = db/SQRT(2.0_dp)
1372 : dxqxz = &
1373 : -2.0_dp/SQRT((r - ab)**2 + (da - ab)**2 + add) + 2.0_dp/SQRT((r + ab)**2 + (da - ab)**2 + add) &
1374 607392 : + 2.0_dp/SQRT((r - ab)**2 + (da + ab)**2 + add) - 2.0_dp/SQRT((r + ab)**2 + (da + ab)**2 + add)
1375 607392 : charg = dxqxz*0.125_dp*fact
1376 607392 : RETURN
1377 : END IF
1378 : ! ZX - ZX
1379 752464 : IF (l1 == 2 .AND. l2 == 2 .AND. m == 1) THEN
1380 303696 : aa = da/SQRT(2.0_dp)
1381 303696 : ab = db/SQRT(2.0_dp)
1382 : qxzqxz = &
1383 : +2.0_dp/SQRT((r + aa - ab)**2 + (aa - ab)**2 + add) - 2.0_dp/SQRT((r + aa + ab)**2 + (aa - ab)**2 + add) &
1384 : - 2.0_dp/SQRT((r - aa - ab)**2 + (aa - ab)**2 + add) + 2.0_dp/SQRT((r - aa + ab)**2 + (aa - ab)**2 + add) &
1385 : - 2.0_dp/SQRT((r + aa - ab)**2 + (aa + ab)**2 + add) + 2.0_dp/SQRT((r + aa + ab)**2 + (aa + ab)**2 + add) &
1386 303696 : + 2.0_dp/SQRT((r - aa - ab)**2 + (aa + ab)**2 + add) - 2.0_dp/SQRT((r - aa + ab)**2 + (aa + ab)**2 + add)
1387 303696 : charg = qxzqxz*0.0625_dp*fact
1388 303696 : RETURN
1389 : END IF
1390 : ! XX - XX
1391 448768 : IF (l1 == 2 .AND. l2 == 2 .AND. m == 2) THEN
1392 448768 : xyxy = 4.0_dp/SQRT(r**2 + (da - db)**2 + add) + 4.0_dp/SQRT(r**2 + (da + db)**2 + add) - 8.0_dp/SQRT(r**2 + da**2 + db**2 + add)
1393 448768 : charg = xyxy*0.0625_dp*fact
1394 448768 : RETURN
1395 : END IF
1396 : ! We should NEVER reach this point
1397 0 : CPABORT("")
1398 0 : END FUNCTION charg_int_ri
1399 :
1400 : ! **************************************************************************************************
1401 : !> \brief Derivatives of the interaction function between two point-charge
1402 : !> configurations (MNDO/d)
1403 : !> Rotational invariant property built-in in the quadrupole definition
1404 : !> r - Distance r12
1405 : !> l1,m - Quantum numbers for multipole of configuration 1
1406 : !> l2,m - Quantum numbers for multipole of configuration 2
1407 : !> da - charge separation of configuration 1
1408 : !> db - charge separation of configuration 2
1409 : !> add - additive term
1410 : !>
1411 : !> \param r ...
1412 : !> \param l1_i ...
1413 : !> \param l2_i ...
1414 : !> \param m ...
1415 : !> \param da_i ...
1416 : !> \param db_i ...
1417 : !> \param add0 ...
1418 : !> \param fact_screen ...
1419 : !> \return ...
1420 : !> \date 03.2008 [tlaino]
1421 : !> \author Teodoro Laino [tlaino] - University of Zurich
1422 : ! **************************************************************************************************
1423 3575779 : FUNCTION dcharg_int_ri(r, l1_i, l2_i, m, da_i, db_i, add0, fact_screen) RESULT(charg)
1424 : REAL(KIND=dp), INTENT(in) :: r
1425 : INTEGER, INTENT(in) :: l1_i, l2_i, m
1426 : REAL(KIND=dp), INTENT(in) :: da_i, db_i, add0, fact_screen
1427 : REAL(KIND=dp) :: charg
1428 :
1429 : INTEGER :: l1, l2
1430 : REAL(KIND=dp) :: aa, ab, add, da, db, dxdx, dxqxz, dzdz, &
1431 : dzqzz, fact, qqzz, qxzqxz, xyxy, zzzz
1432 :
1433 3575779 : IF (l1_i <= l2_i) THEN
1434 2520327 : l1 = l1_i
1435 2520327 : l2 = l2_i
1436 2520327 : da = da_i
1437 2520327 : db = db_i
1438 2520327 : fact = 1.0_dp
1439 : ELSE IF (l1_i > l2_i) THEN
1440 1055452 : l1 = l2_i
1441 1055452 : l2 = l1_i
1442 1055452 : da = db_i
1443 1055452 : db = da_i
1444 1055452 : fact = (-1.0_dp)**(l1 + l2)
1445 : END IF
1446 3575779 : charg = 0.0_dp
1447 3575779 : add = add0*fact_screen
1448 : ! Q - Q.
1449 3575779 : IF (l1 == 0 .AND. l2 == 0) THEN
1450 492107 : charg = r/SQRT(r**2 + add)**3
1451 492107 : charg = -fact*charg
1452 492107 : RETURN
1453 : END IF
1454 : ! Q - Z.
1455 3083672 : IF (l1 == 0 .AND. l2 == 1) THEN
1456 470618 : charg = (r + db)/SQRT((r + db)**2 + add)**3 - (r - db)/SQRT((r - db)**2 + add)**3
1457 470618 : charg = -charg*0.5_dp*fact
1458 470618 : RETURN
1459 : END IF
1460 : ! Z - Z.
1461 2613054 : IF (l1 == 1 .AND. l2 == 1 .AND. m == 0) THEN
1462 : dzdz = &
1463 : +(r + da - db)/SQRT((r + da - db)**2 + add)**3 + (r - da + db)/SQRT((r - da + db)**2 + add)**3 &
1464 90503 : - (r - da - db)/SQRT((r - da - db)**2 + add)**3 - (r + da + db)/SQRT((r + da + db)**2 + add)**3
1465 90503 : charg = -dzdz*0.25_dp*fact
1466 90503 : RETURN
1467 : END IF
1468 : ! X - X
1469 2522551 : IF (l1 == 1 .AND. l2 == 1 .AND. m == 1) THEN
1470 148192 : dxdx = 2.0_dp*r/SQRT(r**2 + (da - db)**2 + add)**3 - 2.0_dp*r/SQRT(r**2 + (da + db)**2 + add)**3
1471 148192 : charg = -dxdx*0.25_dp*fact
1472 148192 : RETURN
1473 : END IF
1474 : ! Q - ZZ
1475 2374359 : IF (l1 == 0 .AND. l2 == 2) THEN
1476 941236 : qqzz = (r - db)/SQRT((r - db)**2 + add)**3 - 2.0_dp*r/SQRT(r**2 + db**2 + add)**3 + (r + db)/SQRT((r + db)**2 + add)**3
1477 941236 : charg = -qqzz*0.25_dp*fact
1478 941236 : RETURN
1479 : END IF
1480 : ! Z - ZZ
1481 1433123 : IF (l1 == 1 .AND. l2 == 2 .AND. m == 0) THEN
1482 : dzqzz = &
1483 : +(r - da - db)/SQRT((r - da - db)**2 + add)**3 - 2.0_dp*(r - da)/SQRT((r - da)**2 + db**2 + add)**3 &
1484 : + (r + db - da)/SQRT((r + db - da)**2 + add)**3 - (r - db + da)/SQRT((r - db + da)**2 + add)**3 &
1485 384876 : + 2.0_dp*(r + da)/SQRT((r + da)**2 + db**2 + add)**3 - (r + da + db)/SQRT((r + da + db)**2 + add)**3
1486 384876 : charg = -dzqzz*0.125_dp*fact
1487 384876 : RETURN
1488 : END IF
1489 : ! ZZ - ZZ
1490 1048247 : IF (l1 == 2 .AND. l2 == 2 .AND. m == 0) THEN
1491 : zzzz = &
1492 : +(r - da - db)/SQRT((r - da - db)**2 + add)**3 + (r + da + db)/SQRT((r + da + db)**2 + add)**3 &
1493 : + (r - da + db)/SQRT((r - da + db)**2 + add)**3 + (r + da - db)/SQRT((r + da - db)**2 + add)**3 &
1494 : - 2.0_dp*(r - da)/SQRT((r - da)**2 + db**2 + add)**3 - 2.0_dp*(r - db)/SQRT((r - db)**2 + da**2 + add)**3 &
1495 : - 2.0_dp*(r + da)/SQRT((r + da)**2 + db**2 + add)**3 - 2.0_dp*(r + db)/SQRT((r + db)**2 + da**2 + add)**3 &
1496 384876 : + 2.0_dp*r/SQRT(r**2 + (da - db)**2 + add)**3 + 2.0_dp*r/SQRT(r**2 + (da + db)**2 + add)**3
1497 : xyxy = &
1498 : +4.0_dp*r/SQRT(r**2 + (da - db)**2 + add)**3 + 4.0_dp*r/SQRT(r**2 + (da + db)**2 + add)**3 &
1499 384876 : - 8.0_dp*r/SQRT(r**2 + da**2 + db**2 + add)**3
1500 384876 : charg = -(zzzz*0.0625_dp - xyxy*0.015625_dp)*fact
1501 384876 : RETURN
1502 : END IF
1503 : ! X - ZX
1504 663371 : IF (l1 == 1 .AND. l2 == 2 .AND. m == 1) THEN
1505 296384 : ab = db/SQRT(2.0_dp)
1506 : dxqxz = &
1507 : -2.0_dp*(r - ab)/SQRT((r - ab)**2 + (da - ab)**2 + add)**3 + 2.0_dp*(r + ab)/SQRT((r + ab)**2 + (da - ab)**2 + add)**3 &
1508 296384 : + 2.0_dp*(r - ab)/SQRT((r - ab)**2 + (da + ab)**2 + add)**3 - 2.0_dp*(r + ab)/SQRT((r + ab)**2 + (da + ab)**2 + add)**3
1509 296384 : charg = -dxqxz*0.125_dp*fact
1510 296384 : RETURN
1511 : END IF
1512 : ! ZX - ZX
1513 366987 : IF (l1 == 2 .AND. l2 == 2 .AND. m == 1) THEN
1514 148192 : aa = da/SQRT(2.0_dp)
1515 148192 : ab = db/SQRT(2.0_dp)
1516 : qxzqxz = &
1517 : +2.0_dp*(r+aa-ab)/SQRT((r+aa-ab)**2+(aa-ab)**2+add)**3-2.0_dp*(r+aa+ab)/SQRT((r+aa+ab)**2+(aa-ab)**2+add)**3 &
1518 : -2.0_dp*(r-aa-ab)/SQRT((r-aa-ab)**2+(aa-ab)**2+add)**3+2.0_dp*(r-aa+ab)/SQRT((r-aa+ab)**2+(aa-ab)**2+add)**3 &
1519 : -2.0_dp*(r+aa-ab)/SQRT((r+aa-ab)**2+(aa+ab)**2+add)**3+2.0_dp*(r+aa+ab)/SQRT((r+aa+ab)**2+(aa+ab)**2+add)**3 &
1520 148192 : +2.0_dp*(r-aa-ab)/SQRT((r-aa-ab)**2+(aa+ab)**2+add)**3-2.0_dp*(r-aa+ab)/SQRT((r-aa+ab)**2+(aa+ab)**2+add)**3
1521 148192 : charg = -qxzqxz*0.0625_dp*fact
1522 148192 : RETURN
1523 : END IF
1524 : ! XX - XX
1525 218795 : IF (l1 == 2 .AND. l2 == 2 .AND. m == 2) THEN
1526 218795 : xyxy = 4.0_dp*r/SQRT(r**2+(da-db)**2+add)**3+4.0_dp*r/SQRT(r**2+(da+db)**2+add)**3-8.0_dp*r/SQRT(r**2+da**2+db**2+add)**3
1527 218795 : charg = -xyxy*0.0625_dp*fact
1528 218795 : RETURN
1529 : END IF
1530 : ! We should NEVER reach this point
1531 0 : CPABORT("")
1532 0 : END FUNCTION dcharg_int_ri
1533 :
1534 : ! **************************************************************************************************
1535 : !> \brief Derivatives of the interaction function between two point-charge
1536 : !> configurations (MNDO/d). This derivative takes into account only the
1537 : !> tapering term
1538 : !> Rotational invariant property built-in in the quadrupole definition
1539 : !> r - Distance r12
1540 : !> l1,m - Quantum numbers for multipole of configuration 1
1541 : !> l2,m - Quantum numbers for multipole of configuration 2
1542 : !> da - charge separation of configuration 1
1543 : !> db - charge separation of configuration 2
1544 : !> add - additive term
1545 : !>
1546 : !> \param r ...
1547 : !> \param l1_i ...
1548 : !> \param l2_i ...
1549 : !> \param m ...
1550 : !> \param da_i ...
1551 : !> \param db_i ...
1552 : !> \param add0 ...
1553 : !> \param fact_screen ...
1554 : !> \return ...
1555 : !> \date 03.2008 [tlaino]
1556 : !> \author Teodoro Laino [tlaino] - University of Zurich
1557 : ! **************************************************************************************************
1558 1208916 : FUNCTION dcharg_int_ri_fs(r, l1_i, l2_i, m, da_i, db_i, add0, fact_screen) RESULT(charg)
1559 : REAL(KIND=dp), INTENT(in) :: r
1560 : INTEGER, INTENT(in) :: l1_i, l2_i, m
1561 : REAL(KIND=dp), INTENT(in) :: da_i, db_i, add0, fact_screen
1562 : REAL(KIND=dp) :: charg
1563 :
1564 : INTEGER :: l1, l2
1565 : REAL(KIND=dp) :: aa, ab, add, da, db, dxdx, dxqxz, dzdz, &
1566 : dzqzz, fact, qqzz, qxzqxz, xyxy, zzzz
1567 :
1568 1208916 : IF (l1_i <= l2_i) THEN
1569 856818 : l1 = l1_i
1570 856818 : l2 = l2_i
1571 856818 : da = da_i
1572 856818 : db = db_i
1573 856818 : fact = 1.0_dp
1574 : ELSE IF (l1_i > l2_i) THEN
1575 352098 : l1 = l2_i
1576 352098 : l2 = l1_i
1577 352098 : da = db_i
1578 352098 : db = da_i
1579 352098 : fact = (-1.0_dp)**(l1 + l2)
1580 : END IF
1581 1208916 : charg = 0.0_dp
1582 1208916 : add = add0*fact_screen
1583 : ! The 0.5 factor handles the derivative of the SQRT
1584 1208916 : fact = fact*0.5_dp
1585 : ! Q - Q.
1586 1208916 : IF (l1 == 0 .AND. l2 == 0) THEN
1587 159552 : charg = add0/SQRT(r**2 + add)**3
1588 159552 : charg = -fact*charg
1589 159552 : RETURN
1590 : END IF
1591 : ! Q - Z.
1592 1049364 : IF (l1 == 0 .AND. l2 == 1) THEN
1593 155448 : charg = add0/SQRT((r + db)**2 + add)**3 - add0/SQRT((r - db)**2 + add)**3
1594 155448 : charg = -charg*0.5_dp*fact
1595 155448 : RETURN
1596 : END IF
1597 : ! Z - Z.
1598 893916 : IF (l1 == 1 .AND. l2 == 1 .AND. m == 0) THEN
1599 : dzdz = &
1600 : +add0/SQRT((r + da - db)**2 + add)**3 + add0/SQRT((r - da + db)**2 + add)**3 &
1601 31572 : - add0/SQRT((r - da - db)**2 + add)**3 - add0/SQRT((r + da + db)**2 + add)**3
1602 31572 : charg = -dzdz*0.25_dp*fact
1603 31572 : RETURN
1604 : END IF
1605 : ! X - X
1606 862344 : IF (l1 == 1 .AND. l2 == 1 .AND. m == 1) THEN
1607 52668 : dxdx = 2.0_dp*add0/SQRT(r**2 + (da - db)**2 + add)**3 - 2.0_dp*add0/SQRT(r**2 + (da + db)**2 + add)**3
1608 52668 : charg = -dxdx*0.25_dp*fact
1609 52668 : RETURN
1610 : END IF
1611 : ! Q - ZZ
1612 809676 : IF (l1 == 0 .AND. l2 == 2) THEN
1613 310896 : qqzz = add0/SQRT((r - db)**2 + add)**3 - 2.0_dp*add0/SQRT(r**2 + db**2 + add)**3 + add0/SQRT((r + db)**2 + add)**3
1614 310896 : charg = -qqzz*0.25_dp*fact
1615 310896 : RETURN
1616 : END IF
1617 : ! Z - ZZ
1618 498780 : IF (l1 == 1 .AND. l2 == 2 .AND. m == 0) THEN
1619 : dzqzz = &
1620 : +add0/SQRT((r - da - db)**2 + add)**3 - 2.0_dp*add0/SQRT((r - da)**2 + db**2 + add)**3 &
1621 : + add0/SQRT((r + db - da)**2 + add)**3 - add0/SQRT((r - db + da)**2 + add)**3 &
1622 132516 : + 2.0_dp*add0/SQRT((r + da)**2 + db**2 + add)**3 - add0/SQRT((r + da + db)**2 + add)**3
1623 132516 : charg = -dzqzz*0.125_dp*fact
1624 132516 : RETURN
1625 : END IF
1626 : ! ZZ - ZZ
1627 366264 : IF (l1 == 2 .AND. l2 == 2 .AND. m == 0) THEN
1628 : zzzz = &
1629 : +add0/SQRT((r - da - db)**2 + add)**3 + add0/SQRT((r + da + db)**2 + add)**3 &
1630 : + add0/SQRT((r - da + db)**2 + add)**3 + add0/SQRT((r + da - db)**2 + add)**3 &
1631 : - 2.0_dp*add0/SQRT((r - da)**2 + db**2 + add)**3 - 2.0_dp*add0/SQRT((r - db)**2 + da**2 + add)**3 &
1632 : - 2.0_dp*add0/SQRT((r + da)**2 + db**2 + add)**3 - 2.0_dp*add0/SQRT((r + db)**2 + da**2 + add)**3 &
1633 132516 : + 2.0_dp*add0/SQRT(r**2 + (da - db)**2 + add)**3 + 2.0_dp*add0/SQRT(r**2 + (da + db)**2 + add)**3
1634 : xyxy = &
1635 : +4.0_dp*add0/SQRT(r**2 + (da - db)**2 + add)**3 + 4.0_dp*add0/SQRT(r**2 + (da + db)**2 + add)**3 &
1636 132516 : - 8.0_dp*add0/SQRT(r**2 + da**2 + db**2 + add)**3
1637 132516 : charg = -(zzzz*0.0625_dp - xyxy*0.015625_dp)*fact
1638 132516 : RETURN
1639 : END IF
1640 : ! X - ZX
1641 233748 : IF (l1 == 1 .AND. l2 == 2 .AND. m == 1) THEN
1642 105336 : ab = db/SQRT(2.0_dp)
1643 : dxqxz = &
1644 : -2.0_dp*add0/SQRT((r - ab)**2 + (da - ab)**2 + add)**3 + 2.0_dp*add0/SQRT((r + ab)**2 + (da - ab)**2 + add)**3 &
1645 105336 : + 2.0_dp*add0/SQRT((r - ab)**2 + (da + ab)**2 + add)**3 - 2.0_dp*add0/SQRT((r + ab)**2 + (da + ab)**2 + add)**3
1646 105336 : charg = -dxqxz*0.125_dp*fact
1647 105336 : RETURN
1648 : END IF
1649 : ! ZX - ZX
1650 128412 : IF (l1 == 2 .AND. l2 == 2 .AND. m == 1) THEN
1651 52668 : aa = da/SQRT(2.0_dp)
1652 52668 : ab = db/SQRT(2.0_dp)
1653 : qxzqxz = &
1654 : +2.0_dp*add0/SQRT((r + aa - ab)**2 + (aa - ab)**2 + add)**3 - 2.0_dp*add0/SQRT((r + aa + ab)**2 + (aa - ab)**2 + add)**3 &
1655 : - 2.0_dp*add0/SQRT((r - aa - ab)**2 + (aa - ab)**2 + add)**3 + 2.0_dp*add0/SQRT((r - aa + ab)**2 + (aa - ab)**2 + add)**3 &
1656 : - 2.0_dp*add0/SQRT((r + aa - ab)**2 + (aa + ab)**2 + add)**3 + 2.0_dp*add0/SQRT((r + aa + ab)**2 + (aa + ab)**2 + add)**3 &
1657 52668 : + 2.0_dp*add0/SQRT((r - aa - ab)**2 + (aa + ab)**2 + add)**3 - 2.0_dp*add0/SQRT((r - aa + ab)**2 + (aa + ab)**2 + add)**3
1658 52668 : charg = -qxzqxz*0.0625_dp*fact
1659 52668 : RETURN
1660 : END IF
1661 : ! XX - XX
1662 75744 : IF (l1 == 2 .AND. l2 == 2 .AND. m == 2) THEN
1663 : xyxy = 4.0_dp*add0/SQRT(r**2 + (da - db)**2 + add)**3 + 4.0_dp*add0/SQRT(r**2 + (da + db)**2 + add)**3 - &
1664 75744 : 8.0_dp*add0/SQRT(r**2 + da**2 + db**2 + add)**3
1665 75744 : charg = -xyxy*0.0625_dp*fact
1666 75744 : RETURN
1667 : END IF
1668 : ! We should NEVER reach this point
1669 0 : CPABORT("")
1670 0 : END FUNCTION dcharg_int_ri_fs
1671 :
1672 : ! **************************************************************************************************
1673 : !> \brief Computes the general rotation matrices for the integrals
1674 : !> between I and J (J>=I)
1675 : !>
1676 : !> \param sepi ...
1677 : !> \param sepj ...
1678 : !> \param rjiv ...
1679 : !> \param r ...
1680 : !> \param ij_matrix ...
1681 : !> \param do_derivatives ...
1682 : !> \param do_invert ...
1683 : !> \param debug_invert ...
1684 : !> \date 04.2008 [tlaino]
1685 : !> \author Teodoro Laino [tlaino] - University of Zurich
1686 : ! **************************************************************************************************
1687 17364991 : RECURSIVE SUBROUTINE rotmat(sepi, sepj, rjiv, r, ij_matrix, do_derivatives, &
1688 : do_invert, debug_invert)
1689 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
1690 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rjiv
1691 : REAL(KIND=dp), INTENT(IN) :: r
1692 : TYPE(rotmat_type), POINTER :: ij_matrix
1693 : LOGICAL, INTENT(IN) :: do_derivatives
1694 : LOGICAL, INTENT(OUT), OPTIONAL :: do_invert
1695 : LOGICAL, INTENT(IN), OPTIONAL :: debug_invert
1696 :
1697 : INTEGER :: imap(3), k, l
1698 : LOGICAL :: dbg_inv, eval
1699 : REAL(KIND=dp) :: b, c2a, c2b, ca, ca2, cb, cb2, check, &
1700 : pt5sq3, r2i, s2a, s2b, sa, sb, sb2, &
1701 : sqb, sqb2i, x11, x22, x33
1702 : REAL(KIND=dp), DIMENSION(3) :: b_d, c2a_d, c2b_d, ca2_d, ca_d, cb2_d, &
1703 : cb_d, r_d, s2a_d, s2b_d, sa_d, sb2_d, &
1704 : sb_d, sqb_d, x11_d, x22_d, x33_d
1705 : REAL(KIND=dp), DIMENSION(3, 3) :: p
1706 : REAL(KIND=dp), DIMENSION(3, 3, 3) :: p_d
1707 : REAL(KIND=dp), DIMENSION(3, 5, 5) :: d_d
1708 : REAL(KIND=dp), DIMENSION(5, 5) :: d
1709 :
1710 17364991 : CPASSERT(ASSOCIATED(ij_matrix))
1711 17364991 : IF (PRESENT(do_invert)) do_invert = .FALSE.
1712 17364991 : IF ((sepi%natorb > 1) .OR. (sepj%natorb > 1)) THEN
1713 : ! Compute Geomtric data and interatomic distance
1714 : ! CA = COS(PHI) , SA = SIN(PHI)
1715 : ! CB = COS(THETA) , SB = SIN(THETA)
1716 : ! C2A = COS(2*PHI) , S2A = SIN(2*PHI)
1717 : ! C2B = COS(2*THETA), S2B = SIN(2*PHI)
1718 8300795 : eval = .TRUE.
1719 8300795 : dbg_inv = .FALSE.
1720 8300795 : pt5sq3 = 0.5_dp*SQRT(3.0_dp)
1721 8300795 : imap(1) = 1
1722 8300795 : imap(2) = 2
1723 8300795 : imap(3) = 3
1724 8300795 : check = rjiv(3)/r
1725 8300795 : IF (PRESENT(debug_invert)) dbg_inv = debug_invert
1726 : ! Check for special case: both atoms lie on the Z Axis
1727 8300795 : IF (ABS(check) > 0.99999999_dp) THEN
1728 5682 : IF (PRESENT(do_invert)) THEN
1729 5389 : imap(1) = 3
1730 5389 : imap(2) = 2
1731 5389 : imap(3) = 1
1732 5389 : eval = .TRUE.
1733 5389 : do_invert = .TRUE.
1734 : ELSE
1735 293 : sa = 0.0_dp
1736 293 : sb = 0.0_dp
1737 293 : IF (check < 0.0_dp) THEN
1738 : ca = -1.0_dp
1739 : cb = -1.0_dp
1740 186 : ELSE IF (check > 0.0_dp) THEN
1741 : ca = 1.0_dp
1742 : cb = 1.0_dp
1743 : ELSE
1744 0 : ca = 0.0_dp
1745 0 : cb = 0.0_dp
1746 : END IF
1747 : eval = .FALSE.
1748 : END IF
1749 : END IF
1750 8300795 : IF (dbg_inv) THEN
1751 : ! When debugging the derivatives of the rotation matrix we must possibly force
1752 : ! the inversion of the reference frame
1753 0 : CPASSERT(.NOT. do_derivatives)
1754 : imap(1) = 3
1755 : imap(2) = 2
1756 : imap(3) = 1
1757 : eval = .TRUE.
1758 : END IF
1759 8300795 : IF (eval) THEN
1760 8300502 : x11 = rjiv(imap(1))
1761 8300502 : x22 = rjiv(imap(2))
1762 8300502 : x33 = rjiv(imap(3))
1763 8300502 : cb = x33/r
1764 8300502 : b = x11**2 + x22**2
1765 8300502 : sqb = SQRT(b)
1766 8300502 : ca = x11/sqb
1767 8300502 : sa = x22/sqb
1768 8300502 : sb = sqb/r
1769 : END IF
1770 : ! Calculate rotation matrix elements
1771 8300795 : p(1, 1) = ca*sb
1772 8300795 : p(2, 1) = ca*cb
1773 8300795 : p(3, 1) = -sa
1774 8300795 : p(1, 2) = sa*sb
1775 8300795 : p(2, 2) = sa*cb
1776 8300795 : p(3, 2) = ca
1777 8300795 : p(1, 3) = cb
1778 8300795 : p(2, 3) = -sb
1779 8300795 : p(3, 3) = 0.0_dp
1780 8300795 : IF (sepi%dorb .OR. sepj%dorb) THEN
1781 92548 : ca2 = ca**2
1782 92548 : cb2 = cb**2
1783 92548 : sb2 = sb**2
1784 92548 : c2a = 2.0_dp*ca2 - 1.0_dp
1785 92548 : c2b = 2.0_dp*cb2 - 1.0_dp
1786 92548 : s2a = 2.0_dp*sa*ca
1787 92548 : s2b = 2.0_dp*sb*cb
1788 92548 : d(1, 1) = pt5sq3*c2a*sb2
1789 92548 : d(2, 1) = 0.5_dp*c2a*s2b
1790 92548 : d(3, 1) = -s2a*sb
1791 92548 : d(4, 1) = c2a*(cb2 + 0.5_dp*sb2)
1792 92548 : d(5, 1) = -s2a*cb
1793 92548 : d(1, 2) = pt5sq3*ca*s2b
1794 92548 : d(2, 2) = ca*c2b
1795 92548 : d(3, 2) = -sa*cb
1796 92548 : d(4, 2) = -0.5_dp*ca*s2b
1797 92548 : d(5, 2) = sa*sb
1798 92548 : d(1, 3) = cb2 - 0.5_dp*sb2
1799 92548 : d(2, 3) = -pt5sq3*s2b
1800 92548 : d(3, 3) = 0.0_dp
1801 92548 : d(4, 3) = pt5sq3*sb2
1802 92548 : d(5, 3) = 0.0_dp
1803 92548 : d(1, 4) = pt5sq3*sa*s2b
1804 92548 : d(2, 4) = sa*c2b
1805 92548 : d(3, 4) = ca*cb
1806 92548 : d(4, 4) = -0.5_dp*sa*s2b
1807 92548 : d(5, 4) = -ca*sb
1808 92548 : d(1, 5) = pt5sq3*s2a*sb2
1809 92548 : d(2, 5) = 0.5_dp*s2a*s2b
1810 92548 : d(3, 5) = c2a*sb
1811 92548 : d(4, 5) = s2a*(cb2 + 0.5_dp*sb2)
1812 92548 : d(5, 5) = c2a*cb
1813 : END IF
1814 : ! Rotation Elements for : S-P
1815 33203180 : DO k = 1, 3
1816 107910335 : DO l = 1, 3
1817 99609540 : ij_matrix%sp(k, l) = p(k, l)
1818 : END DO
1819 : END DO
1820 : ! Rotation Elements for : P-P
1821 33203180 : DO k = 1, 3
1822 24902385 : ij_matrix%pp(1, k, k) = p(k, 1)*p(k, 1)
1823 24902385 : ij_matrix%pp(2, k, k) = p(k, 2)*p(k, 2)
1824 24902385 : ij_matrix%pp(3, k, k) = p(k, 3)*p(k, 3)
1825 24902385 : ij_matrix%pp(4, k, k) = p(k, 1)*p(k, 2)
1826 24902385 : ij_matrix%pp(5, k, k) = p(k, 1)*p(k, 3)
1827 24902385 : ij_matrix%pp(6, k, k) = p(k, 2)*p(k, 3)
1828 33203180 : IF (k /= 1) THEN
1829 41503975 : DO l = 1, k - 1
1830 24902385 : ij_matrix%pp(1, k, l) = 2.0_dp*p(k, 1)*p(l, 1)
1831 24902385 : ij_matrix%pp(2, k, l) = 2.0_dp*p(k, 2)*p(l, 2)
1832 24902385 : ij_matrix%pp(3, k, l) = 2.0_dp*p(k, 3)*p(l, 3)
1833 24902385 : ij_matrix%pp(4, k, l) = p(k, 1)*p(l, 2) + p(k, 2)*p(l, 1)
1834 24902385 : ij_matrix%pp(5, k, l) = p(k, 1)*p(l, 3) + p(k, 3)*p(l, 1)
1835 41503975 : ij_matrix%pp(6, k, l) = p(k, 2)*p(l, 3) + p(k, 3)*p(l, 2)
1836 : END DO
1837 : END IF
1838 : END DO
1839 8300795 : IF (sepi%dorb .OR. sepj%dorb) THEN
1840 : ! Rotation Elements for : S-D
1841 555288 : DO k = 1, 5
1842 2868988 : DO l = 1, 5
1843 2776440 : ij_matrix%sd(k, l) = d(k, l)
1844 : END DO
1845 : END DO
1846 : ! Rotation Elements for : D-P
1847 555288 : DO k = 1, 5
1848 1943508 : DO l = 1, 3
1849 1388220 : ij_matrix%pd(1, k, l) = d(k, 1)*p(l, 1)
1850 1388220 : ij_matrix%pd(2, k, l) = d(k, 1)*p(l, 2)
1851 1388220 : ij_matrix%pd(3, k, l) = d(k, 1)*p(l, 3)
1852 1388220 : ij_matrix%pd(4, k, l) = d(k, 2)*p(l, 1)
1853 1388220 : ij_matrix%pd(5, k, l) = d(k, 2)*p(l, 2)
1854 1388220 : ij_matrix%pd(6, k, l) = d(k, 2)*p(l, 3)
1855 1388220 : ij_matrix%pd(7, k, l) = d(k, 3)*p(l, 1)
1856 1388220 : ij_matrix%pd(8, k, l) = d(k, 3)*p(l, 2)
1857 1388220 : ij_matrix%pd(9, k, l) = d(k, 3)*p(l, 3)
1858 1388220 : ij_matrix%pd(10, k, l) = d(k, 4)*p(l, 1)
1859 1388220 : ij_matrix%pd(11, k, l) = d(k, 4)*p(l, 2)
1860 1388220 : ij_matrix%pd(12, k, l) = d(k, 4)*p(l, 3)
1861 1388220 : ij_matrix%pd(13, k, l) = d(k, 5)*p(l, 1)
1862 1388220 : ij_matrix%pd(14, k, l) = d(k, 5)*p(l, 2)
1863 1850960 : ij_matrix%pd(15, k, l) = d(k, 5)*p(l, 3)
1864 : END DO
1865 : END DO
1866 : ! Rotation Elements for : D-D
1867 555288 : DO k = 1, 5
1868 462740 : ij_matrix%dd(1, k, k) = d(k, 1)*d(k, 1)
1869 462740 : ij_matrix%dd(2, k, k) = d(k, 2)*d(k, 2)
1870 462740 : ij_matrix%dd(3, k, k) = d(k, 3)*d(k, 3)
1871 462740 : ij_matrix%dd(4, k, k) = d(k, 4)*d(k, 4)
1872 462740 : ij_matrix%dd(5, k, k) = d(k, 5)*d(k, 5)
1873 462740 : ij_matrix%dd(6, k, k) = d(k, 1)*d(k, 2)
1874 462740 : ij_matrix%dd(7, k, k) = d(k, 1)*d(k, 3)
1875 462740 : ij_matrix%dd(8, k, k) = d(k, 2)*d(k, 3)
1876 462740 : ij_matrix%dd(9, k, k) = d(k, 1)*d(k, 4)
1877 462740 : ij_matrix%dd(10, k, k) = d(k, 2)*d(k, 4)
1878 462740 : ij_matrix%dd(11, k, k) = d(k, 3)*d(k, 4)
1879 462740 : ij_matrix%dd(12, k, k) = d(k, 1)*d(k, 5)
1880 462740 : ij_matrix%dd(13, k, k) = d(k, 2)*d(k, 5)
1881 462740 : ij_matrix%dd(14, k, k) = d(k, 3)*d(k, 5)
1882 462740 : ij_matrix%dd(15, k, k) = d(k, 4)*d(k, 5)
1883 8763535 : IF (k /= 1) THEN
1884 1295672 : DO l = 1, k - 1
1885 925480 : ij_matrix%dd(1, k, l) = 2.0_dp*d(k, 1)*d(l, 1)
1886 925480 : ij_matrix%dd(2, k, l) = 2.0_dp*d(k, 2)*d(l, 2)
1887 925480 : ij_matrix%dd(3, k, l) = 2.0_dp*d(k, 3)*d(l, 3)
1888 925480 : ij_matrix%dd(4, k, l) = 2.0_dp*d(k, 4)*d(l, 4)
1889 925480 : ij_matrix%dd(5, k, l) = 2.0_dp*d(k, 5)*d(l, 5)
1890 925480 : ij_matrix%dd(6, k, l) = d(k, 1)*d(l, 2) + d(k, 2)*d(l, 1)
1891 925480 : ij_matrix%dd(7, k, l) = d(k, 1)*d(l, 3) + d(k, 3)*d(l, 1)
1892 925480 : ij_matrix%dd(8, k, l) = d(k, 2)*d(l, 3) + d(k, 3)*d(l, 2)
1893 925480 : ij_matrix%dd(9, k, l) = d(k, 1)*d(l, 4) + d(k, 4)*d(l, 1)
1894 925480 : ij_matrix%dd(10, k, l) = d(k, 2)*d(l, 4) + d(k, 4)*d(l, 2)
1895 925480 : ij_matrix%dd(11, k, l) = d(k, 3)*d(l, 4) + d(k, 4)*d(l, 3)
1896 925480 : ij_matrix%dd(12, k, l) = d(k, 1)*d(l, 5) + d(k, 5)*d(l, 1)
1897 925480 : ij_matrix%dd(13, k, l) = d(k, 2)*d(l, 5) + d(k, 5)*d(l, 2)
1898 925480 : ij_matrix%dd(14, k, l) = d(k, 3)*d(l, 5) + d(k, 5)*d(l, 3)
1899 1295672 : ij_matrix%dd(15, k, l) = d(k, 4)*d(l, 5) + d(k, 5)*d(l, 4)
1900 : END DO
1901 : END IF
1902 : END DO
1903 : END IF
1904 : ! Evaluate Derivatives if required
1905 8300795 : IF (do_derivatives) THEN
1906 : ! This condition is necessary because derivatives uses the invertion of the
1907 : ! axis for treating the divergence point along z-axis
1908 4044956 : CPASSERT(eval)
1909 4044956 : x11_d = 0.0_dp; x11_d(1) = 1.0_dp
1910 4044956 : x22_d = 0.0_dp; x22_d(2) = 1.0_dp
1911 4044956 : x33_d = 0.0_dp; x33_d(3) = 1.0_dp
1912 16179824 : r_d = (/x11, x22, x33/)/r
1913 16179824 : b_d = 2.0_dp*(x11*x11_d + x22*x22_d)
1914 16179824 : sqb_d = (0.5_dp/sqb)*b_d
1915 4044956 : r2i = 1.0_dp/r**2
1916 4044956 : sqb2i = 1.0_dp/sqb**2
1917 16179824 : cb_d = (x33_d*r - x33*r_d)*r2i
1918 16179824 : ca_d = (x11_d*sqb - x11*sqb_d)*sqb2i
1919 16179824 : sa_d = (x22_d*sqb - x22*sqb_d)*sqb2i
1920 16179824 : sb_d = (sqb_d*r - sqb*r_d)*r2i
1921 : ! Calculate derivatives of rotation matrix elements
1922 16179824 : p_d(:, 1, 1) = ca_d*sb + ca*sb_d
1923 16179824 : p_d(:, 2, 1) = ca_d*cb + ca*cb_d
1924 16179824 : p_d(:, 3, 1) = -sa_d
1925 16179824 : p_d(:, 1, 2) = sa_d*sb + sa*sb_d
1926 16179824 : p_d(:, 2, 2) = sa_d*cb + sa*cb_d
1927 16179824 : p_d(:, 3, 2) = ca_d
1928 16179824 : p_d(:, 1, 3) = cb_d
1929 16179824 : p_d(:, 2, 3) = -sb_d
1930 16179824 : p_d(:, 3, 3) = 0.0_dp
1931 4044956 : IF (sepi%dorb .OR. sepj%dorb) THEN
1932 171240 : ca2_d = 2.0_dp*ca*ca_d
1933 171240 : cb2_d = 2.0_dp*cb*cb_d
1934 171240 : sb2_d = 2.0_dp*sb*sb_d
1935 171240 : c2a_d = 2.0_dp*ca2_d
1936 171240 : c2b_d = 2.0_dp*cb2_d
1937 171240 : s2a_d = 2.0_dp*(sa_d*ca + sa*ca_d)
1938 171240 : s2b_d = 2.0_dp*(sb_d*cb + sb*cb_d)
1939 171240 : d_d(:, 1, 1) = pt5sq3*(c2a_d*sb2 + c2a*sb2_d)
1940 171240 : d_d(:, 2, 1) = 0.5_dp*(c2a_d*s2b + c2a*s2b_d)
1941 171240 : d_d(:, 3, 1) = -s2a_d*sb - s2a*sb_d
1942 171240 : d_d(:, 4, 1) = c2a_d*(cb2 + 0.5_dp*sb2) + c2a*(cb2_d + 0.5_dp*sb2_d)
1943 171240 : d_d(:, 5, 1) = -s2a_d*cb - s2a*cb_d
1944 171240 : d_d(:, 1, 2) = pt5sq3*(ca_d*s2b + ca*s2b_d)
1945 171240 : d_d(:, 2, 2) = ca_d*c2b + ca*c2b_d
1946 171240 : d_d(:, 3, 2) = -sa_d*cb - sa*cb_d
1947 171240 : d_d(:, 4, 2) = -0.5_dp*(ca_d*s2b + ca*s2b_d)
1948 171240 : d_d(:, 5, 2) = sa_d*sb + sa*sb_d
1949 171240 : d_d(:, 1, 3) = cb2_d - 0.5_dp*sb2_d
1950 171240 : d_d(:, 2, 3) = -pt5sq3*s2b_d
1951 171240 : d_d(:, 3, 3) = 0.0_dp
1952 171240 : d_d(:, 4, 3) = pt5sq3*sb2_d
1953 171240 : d_d(:, 5, 3) = 0.0_dp
1954 171240 : d_d(:, 1, 4) = pt5sq3*(sa_d*s2b + sa*s2b_d)
1955 171240 : d_d(:, 2, 4) = sa_d*c2b + sa*c2b_d
1956 171240 : d_d(:, 3, 4) = ca_d*cb + ca*cb_d
1957 171240 : d_d(:, 4, 4) = -0.5_dp*(sa_d*s2b + sa*s2b_d)
1958 171240 : d_d(:, 5, 4) = -ca_d*sb - ca*sb_d
1959 171240 : d_d(:, 1, 5) = pt5sq3*(s2a_d*sb2 + s2a*sb2_d)
1960 171240 : d_d(:, 2, 5) = 0.5_dp*(s2a_d*s2b + s2a*s2b_d)
1961 171240 : d_d(:, 3, 5) = c2a_d*sb + c2a*sb_d
1962 171240 : d_d(:, 4, 5) = s2a_d*(cb2 + 0.5_dp*sb2) + s2a*(cb2_d + 0.5_dp*sb2_d)
1963 4173386 : d_d(:, 5, 5) = c2a_d*cb + c2a*cb_d
1964 : END IF
1965 : ! Derivative for Rotation Elements for : S-P
1966 16179824 : DO k = 1, 3
1967 52584428 : DO l = 1, 3
1968 157753284 : ij_matrix%sp_d(:, k, l) = p_d(:, k, l)
1969 : END DO
1970 : END DO
1971 : ! Derivative for Rotation Elements for : P-P
1972 16179824 : DO k = 1, 3
1973 48539472 : ij_matrix%pp_d(:, 1, k, k) = p_d(:, k, 1)*p(k, 1) + p(k, 1)*p_d(:, k, 1)
1974 48539472 : ij_matrix%pp_d(:, 2, k, k) = p_d(:, k, 2)*p(k, 2) + p(k, 2)*p_d(:, k, 2)
1975 48539472 : ij_matrix%pp_d(:, 3, k, k) = p_d(:, k, 3)*p(k, 3) + p(k, 3)*p_d(:, k, 3)
1976 48539472 : ij_matrix%pp_d(:, 4, k, k) = p_d(:, k, 1)*p(k, 2) + p(k, 1)*p_d(:, k, 2)
1977 48539472 : ij_matrix%pp_d(:, 5, k, k) = p_d(:, k, 1)*p(k, 3) + p(k, 1)*p_d(:, k, 3)
1978 48539472 : ij_matrix%pp_d(:, 6, k, k) = p_d(:, k, 2)*p(k, 3) + p(k, 2)*p_d(:, k, 3)
1979 16179824 : IF (k /= 1) THEN
1980 20224780 : DO l = 1, k - 1
1981 48539472 : ij_matrix%pp_d(:, 1, k, l) = 2.0_dp*(p_d(:, k, 1)*p(l, 1) + p(k, 1)*p_d(:, l, 1))
1982 48539472 : ij_matrix%pp_d(:, 2, k, l) = 2.0_dp*(p_d(:, k, 2)*p(l, 2) + p(k, 2)*p_d(:, l, 2))
1983 48539472 : ij_matrix%pp_d(:, 3, k, l) = 2.0_dp*(p_d(:, k, 3)*p(l, 3) + p(k, 3)*p_d(:, l, 3))
1984 : ij_matrix%pp_d(:, 4, k, l) = (p_d(:, k, 1)*p(l, 2) + p(k, 1)*p_d(:, l, 2)) + &
1985 48539472 : (p_d(:, k, 2)*p(l, 1) + p(k, 2)*p_d(:, l, 1))
1986 : ij_matrix%pp_d(:, 5, k, l) = (p_d(:, k, 1)*p(l, 3) + p(k, 1)*p_d(:, l, 3)) + &
1987 48539472 : (p_d(:, k, 3)*p(l, 1) + p(k, 3)*p_d(:, l, 1))
1988 : ij_matrix%pp_d(:, 6, k, l) = (p_d(:, k, 2)*p(l, 3) + p(k, 2)*p_d(:, l, 3)) + &
1989 56629384 : (p_d(:, k, 3)*p(l, 2) + p(k, 3)*p_d(:, l, 2))
1990 : END DO
1991 : END IF
1992 : END DO
1993 4044956 : IF (sepi%dorb .OR. sepj%dorb) THEN
1994 : ! Rotation Elements for : S-D
1995 256860 : DO k = 1, 5
1996 1327110 : DO l = 1, 5
1997 4495050 : ij_matrix%sd_d(:, k, l) = d_d(:, k, l)
1998 : END DO
1999 : END DO
2000 : ! Rotation Elements for : D-P
2001 256860 : DO k = 1, 5
2002 899010 : DO l = 1, 3
2003 2568600 : ij_matrix%pd_d(:, 1, k, l) = (d_d(:, k, 1)*p(l, 1) + d(k, 1)*p_d(:, l, 1))
2004 2568600 : ij_matrix%pd_d(:, 2, k, l) = (d_d(:, k, 1)*p(l, 2) + d(k, 1)*p_d(:, l, 2))
2005 2568600 : ij_matrix%pd_d(:, 3, k, l) = (d_d(:, k, 1)*p(l, 3) + d(k, 1)*p_d(:, l, 3))
2006 2568600 : ij_matrix%pd_d(:, 4, k, l) = (d_d(:, k, 2)*p(l, 1) + d(k, 2)*p_d(:, l, 1))
2007 2568600 : ij_matrix%pd_d(:, 5, k, l) = (d_d(:, k, 2)*p(l, 2) + d(k, 2)*p_d(:, l, 2))
2008 2568600 : ij_matrix%pd_d(:, 6, k, l) = (d_d(:, k, 2)*p(l, 3) + d(k, 2)*p_d(:, l, 3))
2009 2568600 : ij_matrix%pd_d(:, 7, k, l) = (d_d(:, k, 3)*p(l, 1) + d(k, 3)*p_d(:, l, 1))
2010 2568600 : ij_matrix%pd_d(:, 8, k, l) = (d_d(:, k, 3)*p(l, 2) + d(k, 3)*p_d(:, l, 2))
2011 2568600 : ij_matrix%pd_d(:, 9, k, l) = (d_d(:, k, 3)*p(l, 3) + d(k, 3)*p_d(:, l, 3))
2012 2568600 : ij_matrix%pd_d(:, 10, k, l) = (d_d(:, k, 4)*p(l, 1) + d(k, 4)*p_d(:, l, 1))
2013 2568600 : ij_matrix%pd_d(:, 11, k, l) = (d_d(:, k, 4)*p(l, 2) + d(k, 4)*p_d(:, l, 2))
2014 2568600 : ij_matrix%pd_d(:, 12, k, l) = (d_d(:, k, 4)*p(l, 3) + d(k, 4)*p_d(:, l, 3))
2015 2568600 : ij_matrix%pd_d(:, 13, k, l) = (d_d(:, k, 5)*p(l, 1) + d(k, 5)*p_d(:, l, 1))
2016 2568600 : ij_matrix%pd_d(:, 14, k, l) = (d_d(:, k, 5)*p(l, 2) + d(k, 5)*p_d(:, l, 2))
2017 2782650 : ij_matrix%pd_d(:, 15, k, l) = (d_d(:, k, 5)*p(l, 3) + d(k, 5)*p_d(:, l, 3))
2018 : END DO
2019 : END DO
2020 : ! Rotation Elements for : D-D
2021 256860 : DO k = 1, 5
2022 856200 : ij_matrix%dd_d(:, 1, k, k) = (d_d(:, k, 1)*d(k, 1) + d(k, 1)*d_d(:, k, 1))
2023 856200 : ij_matrix%dd_d(:, 2, k, k) = (d_d(:, k, 2)*d(k, 2) + d(k, 2)*d_d(:, k, 2))
2024 856200 : ij_matrix%dd_d(:, 3, k, k) = (d_d(:, k, 3)*d(k, 3) + d(k, 3)*d_d(:, k, 3))
2025 856200 : ij_matrix%dd_d(:, 4, k, k) = (d_d(:, k, 4)*d(k, 4) + d(k, 4)*d_d(:, k, 4))
2026 856200 : ij_matrix%dd_d(:, 5, k, k) = (d_d(:, k, 5)*d(k, 5) + d(k, 5)*d_d(:, k, 5))
2027 856200 : ij_matrix%dd_d(:, 6, k, k) = (d_d(:, k, 1)*d(k, 2) + d(k, 1)*d_d(:, k, 2))
2028 856200 : ij_matrix%dd_d(:, 7, k, k) = (d_d(:, k, 1)*d(k, 3) + d(k, 1)*d_d(:, k, 3))
2029 856200 : ij_matrix%dd_d(:, 8, k, k) = (d_d(:, k, 2)*d(k, 3) + d(k, 2)*d_d(:, k, 3))
2030 856200 : ij_matrix%dd_d(:, 9, k, k) = (d_d(:, k, 1)*d(k, 4) + d(k, 1)*d_d(:, k, 4))
2031 856200 : ij_matrix%dd_d(:, 10, k, k) = (d_d(:, k, 2)*d(k, 4) + d(k, 2)*d_d(:, k, 4))
2032 856200 : ij_matrix%dd_d(:, 11, k, k) = (d_d(:, k, 3)*d(k, 4) + d(k, 3)*d_d(:, k, 4))
2033 856200 : ij_matrix%dd_d(:, 12, k, k) = (d_d(:, k, 1)*d(k, 5) + d(k, 1)*d_d(:, k, 5))
2034 856200 : ij_matrix%dd_d(:, 13, k, k) = (d_d(:, k, 2)*d(k, 5) + d(k, 2)*d_d(:, k, 5))
2035 856200 : ij_matrix%dd_d(:, 14, k, k) = (d_d(:, k, 3)*d(k, 5) + d(k, 3)*d_d(:, k, 5))
2036 856200 : ij_matrix%dd_d(:, 15, k, k) = (d_d(:, k, 4)*d(k, 5) + d(k, 4)*d_d(:, k, 5))
2037 4259006 : IF (k /= 1) THEN
2038 599340 : DO l = 1, k - 1
2039 1712400 : ij_matrix%dd_d(:, 1, k, l) = 2.0_dp*(d_d(:, k, 1)*d(l, 1) + d(k, 1)*d_d(:, l, 1))
2040 1712400 : ij_matrix%dd_d(:, 2, k, l) = 2.0_dp*(d_d(:, k, 2)*d(l, 2) + d(k, 2)*d_d(:, l, 2))
2041 1712400 : ij_matrix%dd_d(:, 3, k, l) = 2.0_dp*(d_d(:, k, 3)*d(l, 3) + d(k, 3)*d_d(:, l, 3))
2042 1712400 : ij_matrix%dd_d(:, 4, k, l) = 2.0_dp*(d_d(:, k, 4)*d(l, 4) + d(k, 4)*d_d(:, l, 4))
2043 1712400 : ij_matrix%dd_d(:, 5, k, l) = 2.0_dp*(d_d(:, k, 5)*d(l, 5) + d(k, 5)*d_d(:, l, 5))
2044 : ij_matrix%dd_d(:, 6, k, l) = (d_d(:, k, 1)*d(l, 2) + d(k, 1)*d_d(:, l, 2)) + &
2045 1712400 : (d_d(:, k, 2)*d(l, 1) + d(k, 2)*d_d(:, l, 1))
2046 : ij_matrix%dd_d(:, 7, k, l) = (d_d(:, k, 1)*d(l, 3) + d(k, 1)*d_d(:, l, 3)) + &
2047 1712400 : (d_d(:, k, 3)*d(l, 1) + d(k, 3)*d_d(:, l, 1))
2048 : ij_matrix%dd_d(:, 8, k, l) = (d_d(:, k, 2)*d(l, 3) + d(k, 2)*d_d(:, l, 3)) + &
2049 1712400 : (d_d(:, k, 3)*d(l, 2) + d(k, 3)*d_d(:, l, 2))
2050 : ij_matrix%dd_d(:, 9, k, l) = (d_d(:, k, 1)*d(l, 4) + d(k, 1)*d_d(:, l, 4)) + &
2051 1712400 : (d_d(:, k, 4)*d(l, 1) + d(k, 4)*d_d(:, l, 1))
2052 : ij_matrix%dd_d(:, 10, k, l) = (d_d(:, k, 2)*d(l, 4) + d(k, 2)*d_d(:, l, 4)) + &
2053 1712400 : (d_d(:, k, 4)*d(l, 2) + d(k, 4)*d_d(:, l, 2))
2054 : ij_matrix%dd_d(:, 11, k, l) = (d_d(:, k, 3)*d(l, 4) + d(k, 3)*d_d(:, l, 4)) + &
2055 1712400 : (d_d(:, k, 4)*d(l, 3) + d(k, 4)*d_d(:, l, 3))
2056 : ij_matrix%dd_d(:, 12, k, l) = (d_d(:, k, 1)*d(l, 5) + d(k, 1)*d_d(:, l, 5)) + &
2057 1712400 : (d_d(:, k, 5)*d(l, 1) + d(k, 5)*d_d(:, l, 1))
2058 : ij_matrix%dd_d(:, 13, k, l) = (d_d(:, k, 2)*d(l, 5) + d(k, 2)*d_d(:, l, 5)) + &
2059 1712400 : (d_d(:, k, 5)*d(l, 2) + d(k, 5)*d_d(:, l, 2))
2060 : ij_matrix%dd_d(:, 14, k, l) = (d_d(:, k, 3)*d(l, 5) + d(k, 3)*d_d(:, l, 5)) + &
2061 1712400 : (d_d(:, k, 5)*d(l, 3) + d(k, 5)*d_d(:, l, 3))
2062 : ij_matrix%dd_d(:, 15, k, l) = (d_d(:, k, 4)*d(l, 5) + d(k, 4)*d_d(:, l, 5)) + &
2063 1883640 : (d_d(:, k, 5)*d(l, 4) + d(k, 5)*d_d(:, l, 4))
2064 : END DO
2065 : END IF
2066 : END DO
2067 : END IF
2068 : IF (debug_this_module) THEN
2069 : CALL check_rotmat_der(sepi, sepj, rjiv, ij_matrix, do_invert=do_invert)
2070 : END IF
2071 : END IF
2072 : END IF
2073 17364991 : END SUBROUTINE rotmat
2074 :
2075 : ! **************************************************************************************************
2076 : !> \brief First Step of the rotation of the two-electron two-center integrals in
2077 : !> SPD basis
2078 : !>
2079 : !> \param sepi ...
2080 : !> \param sepj ...
2081 : !> \param rijv ...
2082 : !> \param se_int_control ...
2083 : !> \param se_taper ...
2084 : !> \param invert ...
2085 : !> \param ii ...
2086 : !> \param kk ...
2087 : !> \param rep ...
2088 : !> \param logv ...
2089 : !> \param ij_matrix ...
2090 : !> \param v ...
2091 : !> \param lgrad ...
2092 : !> \param rep_d ...
2093 : !> \param v_d ...
2094 : !> \param logv_d ...
2095 : !> \param drij ...
2096 : !> \date 04.2008 [tlaino]
2097 : !> \author Teodoro Laino [tlaino] - University of Zurich
2098 : ! **************************************************************************************************
2099 1195735 : RECURSIVE SUBROUTINE rot_2el_2c_first(sepi, sepj, rijv, se_int_control, se_taper, &
2100 : invert, ii, kk, rep, logv, ij_matrix, v, lgrad, rep_d, v_d, logv_d, drij)
2101 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
2102 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rijv
2103 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
2104 : TYPE(se_taper_type), POINTER :: se_taper
2105 : LOGICAL, INTENT(IN) :: invert
2106 : INTEGER, INTENT(IN) :: ii, kk
2107 : REAL(KIND=dp), DIMENSION(491), INTENT(IN) :: rep
2108 : LOGICAL, DIMENSION(45, 45), INTENT(OUT) :: logv
2109 : TYPE(rotmat_type), POINTER :: ij_matrix
2110 : REAL(KIND=dp), DIMENSION(45, 45), INTENT(OUT) :: v
2111 : LOGICAL, INTENT(IN) :: lgrad
2112 : REAL(KIND=dp), DIMENSION(491), INTENT(IN), &
2113 : OPTIONAL :: rep_d
2114 : REAL(KIND=dp), DIMENSION(3, 45, 45), INTENT(OUT), &
2115 : OPTIONAL :: v_d
2116 : LOGICAL, DIMENSION(45, 45), INTENT(OUT), OPTIONAL :: logv_d
2117 : REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: drij
2118 :
2119 : INTEGER :: i, i1, ij, j1, k, k1, kl, l, l1, limkl, &
2120 : ll, mm, nd
2121 : REAL(KIND=dp) :: wrepp, wrepp_d(3)
2122 :
2123 1195735 : IF (lgrad) THEN
2124 511306 : CPASSERT(PRESENT(rep_d))
2125 511306 : CPASSERT(PRESENT(v_d))
2126 511306 : CPASSERT(PRESENT(logv_d))
2127 511306 : CPASSERT(PRESENT(drij))
2128 : END IF
2129 1195735 : limkl = indexb(kk, kk)
2130 8231336 : DO k = 1, limkl
2131 324833381 : DO i = 1, 45
2132 316602045 : logv(i, k) = .FALSE.
2133 323637646 : v(i, k) = 0.0_dp
2134 : END DO
2135 : END DO
2136 : !
2137 4684623 : DO i1 = 1, ii
2138 13422257 : DO j1 = 1, i1
2139 8737634 : ij = indexa(i1, j1)
2140 : !
2141 33847168 : DO k1 = 1, kk
2142 : !
2143 100361142 : DO l1 = 1, k1
2144 66513974 : kl = indexa(k1, l1)
2145 66513974 : nd = ijkl_ind(ij, kl)
2146 91623508 : IF (nd /= 0) THEN
2147 : !
2148 21392781 : wrepp = rep(nd)
2149 21392781 : ll = indexb(k1, l1)
2150 21392781 : mm = int2c_type(ll)
2151 : !
2152 : IF (mm == 1) THEN
2153 4330557 : v(ij, 1) = wrepp
2154 : ELSE IF (mm == 2) THEN
2155 4121126 : k = k1 - 1
2156 4121126 : v(ij, 2) = v(ij, 2) + ij_matrix%sp(k, 1)*wrepp
2157 4121126 : v(ij, 4) = v(ij, 4) + ij_matrix%sp(k, 2)*wrepp
2158 4121126 : v(ij, 7) = v(ij, 7) + ij_matrix%sp(k, 3)*wrepp
2159 : ELSE IF (mm == 3) THEN
2160 9277700 : k = k1 - 1
2161 9277700 : l = l1 - 1
2162 9277700 : v(ij, 3) = v(ij, 3) + ij_matrix%pp(1, k, l)*wrepp
2163 9277700 : v(ij, 6) = v(ij, 6) + ij_matrix%pp(2, k, l)*wrepp
2164 9277700 : v(ij, 10) = v(ij, 10) + ij_matrix%pp(3, k, l)*wrepp
2165 9277700 : v(ij, 5) = v(ij, 5) + ij_matrix%pp(4, k, l)*wrepp
2166 9277700 : v(ij, 8) = v(ij, 8) + ij_matrix%pp(5, k, l)*wrepp
2167 9277700 : v(ij, 9) = v(ij, 9) + ij_matrix%pp(6, k, l)*wrepp
2168 : ELSE IF (mm == 4) THEN
2169 495092 : k = k1 - 4
2170 495092 : v(ij, 11) = v(ij, 11) + ij_matrix%sd(k, 1)*wrepp
2171 495092 : v(ij, 16) = v(ij, 16) + ij_matrix%sd(k, 2)*wrepp
2172 495092 : v(ij, 22) = v(ij, 22) + ij_matrix%sd(k, 3)*wrepp
2173 495092 : v(ij, 29) = v(ij, 29) + ij_matrix%sd(k, 4)*wrepp
2174 495092 : v(ij, 37) = v(ij, 37) + ij_matrix%sd(k, 5)*wrepp
2175 : ELSE IF (mm == 5) THEN
2176 1532316 : k = k1 - 4
2177 1532316 : l = l1 - 1
2178 1532316 : v(ij, 12) = v(ij, 12) + ij_matrix%pd(1, k, l)*wrepp
2179 1532316 : v(ij, 13) = v(ij, 13) + ij_matrix%pd(2, k, l)*wrepp
2180 1532316 : v(ij, 14) = v(ij, 14) + ij_matrix%pd(3, k, l)*wrepp
2181 1532316 : v(ij, 17) = v(ij, 17) + ij_matrix%pd(4, k, l)*wrepp
2182 1532316 : v(ij, 18) = v(ij, 18) + ij_matrix%pd(5, k, l)*wrepp
2183 1532316 : v(ij, 19) = v(ij, 19) + ij_matrix%pd(6, k, l)*wrepp
2184 1532316 : v(ij, 23) = v(ij, 23) + ij_matrix%pd(7, k, l)*wrepp
2185 1532316 : v(ij, 24) = v(ij, 24) + ij_matrix%pd(8, k, l)*wrepp
2186 1532316 : v(ij, 25) = v(ij, 25) + ij_matrix%pd(9, k, l)*wrepp
2187 1532316 : v(ij, 30) = v(ij, 30) + ij_matrix%pd(10, k, l)*wrepp
2188 1532316 : v(ij, 31) = v(ij, 31) + ij_matrix%pd(11, k, l)*wrepp
2189 1532316 : v(ij, 32) = v(ij, 32) + ij_matrix%pd(12, k, l)*wrepp
2190 1532316 : v(ij, 38) = v(ij, 38) + ij_matrix%pd(13, k, l)*wrepp
2191 1532316 : v(ij, 39) = v(ij, 39) + ij_matrix%pd(14, k, l)*wrepp
2192 1532316 : v(ij, 40) = v(ij, 40) + ij_matrix%pd(15, k, l)*wrepp
2193 : ELSE IF (mm == 6) THEN
2194 1635990 : k = k1 - 4
2195 1635990 : l = l1 - 4
2196 1635990 : v(ij, 15) = v(ij, 15) + ij_matrix%dd(1, k, l)*wrepp
2197 1635990 : v(ij, 21) = v(ij, 21) + ij_matrix%dd(2, k, l)*wrepp
2198 1635990 : v(ij, 28) = v(ij, 28) + ij_matrix%dd(3, k, l)*wrepp
2199 1635990 : v(ij, 36) = v(ij, 36) + ij_matrix%dd(4, k, l)*wrepp
2200 1635990 : v(ij, 45) = v(ij, 45) + ij_matrix%dd(5, k, l)*wrepp
2201 1635990 : v(ij, 20) = v(ij, 20) + ij_matrix%dd(6, k, l)*wrepp
2202 1635990 : v(ij, 26) = v(ij, 26) + ij_matrix%dd(7, k, l)*wrepp
2203 1635990 : v(ij, 27) = v(ij, 27) + ij_matrix%dd(8, k, l)*wrepp
2204 1635990 : v(ij, 33) = v(ij, 33) + ij_matrix%dd(9, k, l)*wrepp
2205 1635990 : v(ij, 34) = v(ij, 34) + ij_matrix%dd(10, k, l)*wrepp
2206 1635990 : v(ij, 35) = v(ij, 35) + ij_matrix%dd(11, k, l)*wrepp
2207 1635990 : v(ij, 41) = v(ij, 41) + ij_matrix%dd(12, k, l)*wrepp
2208 1635990 : v(ij, 42) = v(ij, 42) + ij_matrix%dd(13, k, l)*wrepp
2209 1635990 : v(ij, 43) = v(ij, 43) + ij_matrix%dd(14, k, l)*wrepp
2210 1635990 : v(ij, 44) = v(ij, 44) + ij_matrix%dd(15, k, l)*wrepp
2211 : END IF
2212 : END IF
2213 : END DO
2214 : END DO
2215 78740496 : DO kl = 1, limkl
2216 75251608 : logv(ij, kl) = (ABS(v(ij, kl)) > 0.0_dp)
2217 : END DO
2218 : END DO
2219 : END DO
2220 : ! Gradients
2221 1195735 : IF (lgrad) THEN
2222 3642985 : DO k = 1, limkl
2223 144568540 : DO i = 1, 45
2224 140925555 : logv_d(i, k) = .FALSE.
2225 140925555 : v_d(1, i, k) = 0.0_dp
2226 140925555 : v_d(2, i, k) = 0.0_dp
2227 144057234 : v_d(3, i, k) = 0.0_dp
2228 : END DO
2229 : END DO
2230 2018237 : DO i1 = 1, ii
2231 5821698 : DO j1 = 1, i1
2232 3803461 : ij = indexa(i1, j1)
2233 : !
2234 15073799 : DO k1 = 1, kk
2235 : !
2236 45478491 : DO l1 = 1, k1
2237 30404692 : kl = indexa(k1, l1)
2238 30404692 : nd = ijkl_ind(ij, kl)
2239 41675030 : IF (nd /= 0) THEN
2240 : !
2241 9688399 : wrepp = rep(nd)
2242 38753596 : wrepp_d = rep_d(nd)*drij
2243 9688399 : ll = indexb(k1, l1)
2244 9688399 : mm = int2c_type(ll)
2245 : !
2246 : IF (mm == 1) THEN
2247 1874422 : v_d(1, ij, 1) = wrepp_d(1)
2248 1874422 : v_d(2, ij, 1) = wrepp_d(2)
2249 1874422 : v_d(3, ij, 1) = wrepp_d(3)
2250 : ELSE IF (mm == 2) THEN
2251 1857556 : k = k1 - 1
2252 1857556 : v_d(1, ij, 2) = v_d(1, ij, 2) + ij_matrix%sp_d(1, k, 1)*wrepp + ij_matrix%sp(k, 1)*wrepp_d(1)
2253 1857556 : v_d(1, ij, 4) = v_d(1, ij, 4) + ij_matrix%sp_d(1, k, 2)*wrepp + ij_matrix%sp(k, 2)*wrepp_d(1)
2254 1857556 : v_d(1, ij, 7) = v_d(1, ij, 7) + ij_matrix%sp_d(1, k, 3)*wrepp + ij_matrix%sp(k, 3)*wrepp_d(1)
2255 :
2256 1857556 : v_d(2, ij, 2) = v_d(2, ij, 2) + ij_matrix%sp_d(2, k, 1)*wrepp + ij_matrix%sp(k, 1)*wrepp_d(2)
2257 1857556 : v_d(2, ij, 4) = v_d(2, ij, 4) + ij_matrix%sp_d(2, k, 2)*wrepp + ij_matrix%sp(k, 2)*wrepp_d(2)
2258 1857556 : v_d(2, ij, 7) = v_d(2, ij, 7) + ij_matrix%sp_d(2, k, 3)*wrepp + ij_matrix%sp(k, 3)*wrepp_d(2)
2259 :
2260 1857556 : v_d(3, ij, 2) = v_d(3, ij, 2) + ij_matrix%sp_d(3, k, 1)*wrepp + ij_matrix%sp(k, 1)*wrepp_d(3)
2261 1857556 : v_d(3, ij, 4) = v_d(3, ij, 4) + ij_matrix%sp_d(3, k, 2)*wrepp + ij_matrix%sp(k, 2)*wrepp_d(3)
2262 1857556 : v_d(3, ij, 7) = v_d(3, ij, 7) + ij_matrix%sp_d(3, k, 3)*wrepp + ij_matrix%sp(k, 3)*wrepp_d(3)
2263 : ELSE IF (mm == 3) THEN
2264 4179597 : k = k1 - 1
2265 4179597 : l = l1 - 1
2266 4179597 : v_d(1, ij, 3) = v_d(1, ij, 3) + ij_matrix%pp_d(1, 1, k, l)*wrepp + ij_matrix%pp(1, k, l)*wrepp_d(1)
2267 4179597 : v_d(1, ij, 6) = v_d(1, ij, 6) + ij_matrix%pp_d(1, 2, k, l)*wrepp + ij_matrix%pp(2, k, l)*wrepp_d(1)
2268 4179597 : v_d(1, ij, 10) = v_d(1, ij, 10) + ij_matrix%pp_d(1, 3, k, l)*wrepp + ij_matrix%pp(3, k, l)*wrepp_d(1)
2269 4179597 : v_d(1, ij, 5) = v_d(1, ij, 5) + ij_matrix%pp_d(1, 4, k, l)*wrepp + ij_matrix%pp(4, k, l)*wrepp_d(1)
2270 4179597 : v_d(1, ij, 8) = v_d(1, ij, 8) + ij_matrix%pp_d(1, 5, k, l)*wrepp + ij_matrix%pp(5, k, l)*wrepp_d(1)
2271 4179597 : v_d(1, ij, 9) = v_d(1, ij, 9) + ij_matrix%pp_d(1, 6, k, l)*wrepp + ij_matrix%pp(6, k, l)*wrepp_d(1)
2272 :
2273 4179597 : v_d(2, ij, 3) = v_d(2, ij, 3) + ij_matrix%pp_d(2, 1, k, l)*wrepp + ij_matrix%pp(1, k, l)*wrepp_d(2)
2274 4179597 : v_d(2, ij, 6) = v_d(2, ij, 6) + ij_matrix%pp_d(2, 2, k, l)*wrepp + ij_matrix%pp(2, k, l)*wrepp_d(2)
2275 4179597 : v_d(2, ij, 10) = v_d(2, ij, 10) + ij_matrix%pp_d(2, 3, k, l)*wrepp + ij_matrix%pp(3, k, l)*wrepp_d(2)
2276 4179597 : v_d(2, ij, 5) = v_d(2, ij, 5) + ij_matrix%pp_d(2, 4, k, l)*wrepp + ij_matrix%pp(4, k, l)*wrepp_d(2)
2277 4179597 : v_d(2, ij, 8) = v_d(2, ij, 8) + ij_matrix%pp_d(2, 5, k, l)*wrepp + ij_matrix%pp(5, k, l)*wrepp_d(2)
2278 4179597 : v_d(2, ij, 9) = v_d(2, ij, 9) + ij_matrix%pp_d(2, 6, k, l)*wrepp + ij_matrix%pp(6, k, l)*wrepp_d(2)
2279 :
2280 4179597 : v_d(3, ij, 3) = v_d(3, ij, 3) + ij_matrix%pp_d(3, 1, k, l)*wrepp + ij_matrix%pp(1, k, l)*wrepp_d(3)
2281 4179597 : v_d(3, ij, 6) = v_d(3, ij, 6) + ij_matrix%pp_d(3, 2, k, l)*wrepp + ij_matrix%pp(2, k, l)*wrepp_d(3)
2282 4179597 : v_d(3, ij, 10) = v_d(3, ij, 10) + ij_matrix%pp_d(3, 3, k, l)*wrepp + ij_matrix%pp(3, k, l)*wrepp_d(3)
2283 4179597 : v_d(3, ij, 5) = v_d(3, ij, 5) + ij_matrix%pp_d(3, 4, k, l)*wrepp + ij_matrix%pp(4, k, l)*wrepp_d(3)
2284 4179597 : v_d(3, ij, 8) = v_d(3, ij, 8) + ij_matrix%pp_d(3, 5, k, l)*wrepp + ij_matrix%pp(5, k, l)*wrepp_d(3)
2285 4179597 : v_d(3, ij, 9) = v_d(3, ij, 9) + ij_matrix%pp_d(3, 6, k, l)*wrepp + ij_matrix%pp(6, k, l)*wrepp_d(3)
2286 : ELSE IF (mm == 4) THEN
2287 240203 : k = k1 - 4
2288 240203 : v_d(1, ij, 11) = v_d(1, ij, 11) + ij_matrix%sd_d(1, k, 1)*wrepp + ij_matrix%sd(k, 1)*wrepp_d(1)
2289 240203 : v_d(1, ij, 16) = v_d(1, ij, 16) + ij_matrix%sd_d(1, k, 2)*wrepp + ij_matrix%sd(k, 2)*wrepp_d(1)
2290 240203 : v_d(1, ij, 22) = v_d(1, ij, 22) + ij_matrix%sd_d(1, k, 3)*wrepp + ij_matrix%sd(k, 3)*wrepp_d(1)
2291 240203 : v_d(1, ij, 29) = v_d(1, ij, 29) + ij_matrix%sd_d(1, k, 4)*wrepp + ij_matrix%sd(k, 4)*wrepp_d(1)
2292 240203 : v_d(1, ij, 37) = v_d(1, ij, 37) + ij_matrix%sd_d(1, k, 5)*wrepp + ij_matrix%sd(k, 5)*wrepp_d(1)
2293 :
2294 240203 : v_d(2, ij, 11) = v_d(2, ij, 11) + ij_matrix%sd_d(2, k, 1)*wrepp + ij_matrix%sd(k, 1)*wrepp_d(2)
2295 240203 : v_d(2, ij, 16) = v_d(2, ij, 16) + ij_matrix%sd_d(2, k, 2)*wrepp + ij_matrix%sd(k, 2)*wrepp_d(2)
2296 240203 : v_d(2, ij, 22) = v_d(2, ij, 22) + ij_matrix%sd_d(2, k, 3)*wrepp + ij_matrix%sd(k, 3)*wrepp_d(2)
2297 240203 : v_d(2, ij, 29) = v_d(2, ij, 29) + ij_matrix%sd_d(2, k, 4)*wrepp + ij_matrix%sd(k, 4)*wrepp_d(2)
2298 240203 : v_d(2, ij, 37) = v_d(2, ij, 37) + ij_matrix%sd_d(2, k, 5)*wrepp + ij_matrix%sd(k, 5)*wrepp_d(2)
2299 :
2300 240203 : v_d(3, ij, 11) = v_d(3, ij, 11) + ij_matrix%sd_d(3, k, 1)*wrepp + ij_matrix%sd(k, 1)*wrepp_d(3)
2301 240203 : v_d(3, ij, 16) = v_d(3, ij, 16) + ij_matrix%sd_d(3, k, 2)*wrepp + ij_matrix%sd(k, 2)*wrepp_d(3)
2302 240203 : v_d(3, ij, 22) = v_d(3, ij, 22) + ij_matrix%sd_d(3, k, 3)*wrepp + ij_matrix%sd(k, 3)*wrepp_d(3)
2303 240203 : v_d(3, ij, 29) = v_d(3, ij, 29) + ij_matrix%sd_d(3, k, 4)*wrepp + ij_matrix%sd(k, 4)*wrepp_d(3)
2304 240203 : v_d(3, ij, 37) = v_d(3, ij, 37) + ij_matrix%sd_d(3, k, 5)*wrepp + ij_matrix%sd(k, 5)*wrepp_d(3)
2305 : ELSE IF (mm == 5) THEN
2306 743417 : k = k1 - 4
2307 743417 : l = l1 - 1
2308 743417 : v_d(1, ij, 12) = v_d(1, ij, 12) + ij_matrix%pd_d(1, 1, k, l)*wrepp + ij_matrix%pd(1, k, l)*wrepp_d(1)
2309 743417 : v_d(1, ij, 13) = v_d(1, ij, 13) + ij_matrix%pd_d(1, 2, k, l)*wrepp + ij_matrix%pd(2, k, l)*wrepp_d(1)
2310 743417 : v_d(1, ij, 14) = v_d(1, ij, 14) + ij_matrix%pd_d(1, 3, k, l)*wrepp + ij_matrix%pd(3, k, l)*wrepp_d(1)
2311 743417 : v_d(1, ij, 17) = v_d(1, ij, 17) + ij_matrix%pd_d(1, 4, k, l)*wrepp + ij_matrix%pd(4, k, l)*wrepp_d(1)
2312 743417 : v_d(1, ij, 18) = v_d(1, ij, 18) + ij_matrix%pd_d(1, 5, k, l)*wrepp + ij_matrix%pd(5, k, l)*wrepp_d(1)
2313 743417 : v_d(1, ij, 19) = v_d(1, ij, 19) + ij_matrix%pd_d(1, 6, k, l)*wrepp + ij_matrix%pd(6, k, l)*wrepp_d(1)
2314 743417 : v_d(1, ij, 23) = v_d(1, ij, 23) + ij_matrix%pd_d(1, 7, k, l)*wrepp + ij_matrix%pd(7, k, l)*wrepp_d(1)
2315 743417 : v_d(1, ij, 24) = v_d(1, ij, 24) + ij_matrix%pd_d(1, 8, k, l)*wrepp + ij_matrix%pd(8, k, l)*wrepp_d(1)
2316 743417 : v_d(1, ij, 25) = v_d(1, ij, 25) + ij_matrix%pd_d(1, 9, k, l)*wrepp + ij_matrix%pd(9, k, l)*wrepp_d(1)
2317 743417 : v_d(1, ij, 30) = v_d(1, ij, 30) + ij_matrix%pd_d(1, 10, k, l)*wrepp + ij_matrix%pd(10, k, l)*wrepp_d(1)
2318 743417 : v_d(1, ij, 31) = v_d(1, ij, 31) + ij_matrix%pd_d(1, 11, k, l)*wrepp + ij_matrix%pd(11, k, l)*wrepp_d(1)
2319 743417 : v_d(1, ij, 32) = v_d(1, ij, 32) + ij_matrix%pd_d(1, 12, k, l)*wrepp + ij_matrix%pd(12, k, l)*wrepp_d(1)
2320 743417 : v_d(1, ij, 38) = v_d(1, ij, 38) + ij_matrix%pd_d(1, 13, k, l)*wrepp + ij_matrix%pd(13, k, l)*wrepp_d(1)
2321 743417 : v_d(1, ij, 39) = v_d(1, ij, 39) + ij_matrix%pd_d(1, 14, k, l)*wrepp + ij_matrix%pd(14, k, l)*wrepp_d(1)
2322 743417 : v_d(1, ij, 40) = v_d(1, ij, 40) + ij_matrix%pd_d(1, 15, k, l)*wrepp + ij_matrix%pd(15, k, l)*wrepp_d(1)
2323 :
2324 743417 : v_d(2, ij, 12) = v_d(2, ij, 12) + ij_matrix%pd_d(2, 1, k, l)*wrepp + ij_matrix%pd(1, k, l)*wrepp_d(2)
2325 743417 : v_d(2, ij, 13) = v_d(2, ij, 13) + ij_matrix%pd_d(2, 2, k, l)*wrepp + ij_matrix%pd(2, k, l)*wrepp_d(2)
2326 743417 : v_d(2, ij, 14) = v_d(2, ij, 14) + ij_matrix%pd_d(2, 3, k, l)*wrepp + ij_matrix%pd(3, k, l)*wrepp_d(2)
2327 743417 : v_d(2, ij, 17) = v_d(2, ij, 17) + ij_matrix%pd_d(2, 4, k, l)*wrepp + ij_matrix%pd(4, k, l)*wrepp_d(2)
2328 743417 : v_d(2, ij, 18) = v_d(2, ij, 18) + ij_matrix%pd_d(2, 5, k, l)*wrepp + ij_matrix%pd(5, k, l)*wrepp_d(2)
2329 743417 : v_d(2, ij, 19) = v_d(2, ij, 19) + ij_matrix%pd_d(2, 6, k, l)*wrepp + ij_matrix%pd(6, k, l)*wrepp_d(2)
2330 743417 : v_d(2, ij, 23) = v_d(2, ij, 23) + ij_matrix%pd_d(2, 7, k, l)*wrepp + ij_matrix%pd(7, k, l)*wrepp_d(2)
2331 743417 : v_d(2, ij, 24) = v_d(2, ij, 24) + ij_matrix%pd_d(2, 8, k, l)*wrepp + ij_matrix%pd(8, k, l)*wrepp_d(2)
2332 743417 : v_d(2, ij, 25) = v_d(2, ij, 25) + ij_matrix%pd_d(2, 9, k, l)*wrepp + ij_matrix%pd(9, k, l)*wrepp_d(2)
2333 743417 : v_d(2, ij, 30) = v_d(2, ij, 30) + ij_matrix%pd_d(2, 10, k, l)*wrepp + ij_matrix%pd(10, k, l)*wrepp_d(2)
2334 743417 : v_d(2, ij, 31) = v_d(2, ij, 31) + ij_matrix%pd_d(2, 11, k, l)*wrepp + ij_matrix%pd(11, k, l)*wrepp_d(2)
2335 743417 : v_d(2, ij, 32) = v_d(2, ij, 32) + ij_matrix%pd_d(2, 12, k, l)*wrepp + ij_matrix%pd(12, k, l)*wrepp_d(2)
2336 743417 : v_d(2, ij, 38) = v_d(2, ij, 38) + ij_matrix%pd_d(2, 13, k, l)*wrepp + ij_matrix%pd(13, k, l)*wrepp_d(2)
2337 743417 : v_d(2, ij, 39) = v_d(2, ij, 39) + ij_matrix%pd_d(2, 14, k, l)*wrepp + ij_matrix%pd(14, k, l)*wrepp_d(2)
2338 743417 : v_d(2, ij, 40) = v_d(2, ij, 40) + ij_matrix%pd_d(2, 15, k, l)*wrepp + ij_matrix%pd(15, k, l)*wrepp_d(2)
2339 :
2340 743417 : v_d(3, ij, 12) = v_d(3, ij, 12) + ij_matrix%pd_d(3, 1, k, l)*wrepp + ij_matrix%pd(1, k, l)*wrepp_d(3)
2341 743417 : v_d(3, ij, 13) = v_d(3, ij, 13) + ij_matrix%pd_d(3, 2, k, l)*wrepp + ij_matrix%pd(2, k, l)*wrepp_d(3)
2342 743417 : v_d(3, ij, 14) = v_d(3, ij, 14) + ij_matrix%pd_d(3, 3, k, l)*wrepp + ij_matrix%pd(3, k, l)*wrepp_d(3)
2343 743417 : v_d(3, ij, 17) = v_d(3, ij, 17) + ij_matrix%pd_d(3, 4, k, l)*wrepp + ij_matrix%pd(4, k, l)*wrepp_d(3)
2344 743417 : v_d(3, ij, 18) = v_d(3, ij, 18) + ij_matrix%pd_d(3, 5, k, l)*wrepp + ij_matrix%pd(5, k, l)*wrepp_d(3)
2345 743417 : v_d(3, ij, 19) = v_d(3, ij, 19) + ij_matrix%pd_d(3, 6, k, l)*wrepp + ij_matrix%pd(6, k, l)*wrepp_d(3)
2346 743417 : v_d(3, ij, 23) = v_d(3, ij, 23) + ij_matrix%pd_d(3, 7, k, l)*wrepp + ij_matrix%pd(7, k, l)*wrepp_d(3)
2347 743417 : v_d(3, ij, 24) = v_d(3, ij, 24) + ij_matrix%pd_d(3, 8, k, l)*wrepp + ij_matrix%pd(8, k, l)*wrepp_d(3)
2348 743417 : v_d(3, ij, 25) = v_d(3, ij, 25) + ij_matrix%pd_d(3, 9, k, l)*wrepp + ij_matrix%pd(9, k, l)*wrepp_d(3)
2349 743417 : v_d(3, ij, 30) = v_d(3, ij, 30) + ij_matrix%pd_d(3, 10, k, l)*wrepp + ij_matrix%pd(10, k, l)*wrepp_d(3)
2350 743417 : v_d(3, ij, 31) = v_d(3, ij, 31) + ij_matrix%pd_d(3, 11, k, l)*wrepp + ij_matrix%pd(11, k, l)*wrepp_d(3)
2351 743417 : v_d(3, ij, 32) = v_d(3, ij, 32) + ij_matrix%pd_d(3, 12, k, l)*wrepp + ij_matrix%pd(12, k, l)*wrepp_d(3)
2352 743417 : v_d(3, ij, 38) = v_d(3, ij, 38) + ij_matrix%pd_d(3, 13, k, l)*wrepp + ij_matrix%pd(13, k, l)*wrepp_d(3)
2353 743417 : v_d(3, ij, 39) = v_d(3, ij, 39) + ij_matrix%pd_d(3, 14, k, l)*wrepp + ij_matrix%pd(14, k, l)*wrepp_d(3)
2354 743417 : v_d(3, ij, 40) = v_d(3, ij, 40) + ij_matrix%pd_d(3, 15, k, l)*wrepp + ij_matrix%pd(15, k, l)*wrepp_d(3)
2355 : ELSE IF (mm == 6) THEN
2356 793204 : k = k1 - 4
2357 793204 : l = l1 - 4
2358 793204 : v_d(1, ij, 15) = v_d(1, ij, 15) + ij_matrix%dd_d(1, 1, k, l)*wrepp + ij_matrix%dd(1, k, l)*wrepp_d(1)
2359 793204 : v_d(1, ij, 21) = v_d(1, ij, 21) + ij_matrix%dd_d(1, 2, k, l)*wrepp + ij_matrix%dd(2, k, l)*wrepp_d(1)
2360 793204 : v_d(1, ij, 28) = v_d(1, ij, 28) + ij_matrix%dd_d(1, 3, k, l)*wrepp + ij_matrix%dd(3, k, l)*wrepp_d(1)
2361 793204 : v_d(1, ij, 36) = v_d(1, ij, 36) + ij_matrix%dd_d(1, 4, k, l)*wrepp + ij_matrix%dd(4, k, l)*wrepp_d(1)
2362 793204 : v_d(1, ij, 45) = v_d(1, ij, 45) + ij_matrix%dd_d(1, 5, k, l)*wrepp + ij_matrix%dd(5, k, l)*wrepp_d(1)
2363 793204 : v_d(1, ij, 20) = v_d(1, ij, 20) + ij_matrix%dd_d(1, 6, k, l)*wrepp + ij_matrix%dd(6, k, l)*wrepp_d(1)
2364 793204 : v_d(1, ij, 26) = v_d(1, ij, 26) + ij_matrix%dd_d(1, 7, k, l)*wrepp + ij_matrix%dd(7, k, l)*wrepp_d(1)
2365 793204 : v_d(1, ij, 27) = v_d(1, ij, 27) + ij_matrix%dd_d(1, 8, k, l)*wrepp + ij_matrix%dd(8, k, l)*wrepp_d(1)
2366 793204 : v_d(1, ij, 33) = v_d(1, ij, 33) + ij_matrix%dd_d(1, 9, k, l)*wrepp + ij_matrix%dd(9, k, l)*wrepp_d(1)
2367 793204 : v_d(1, ij, 34) = v_d(1, ij, 34) + ij_matrix%dd_d(1, 10, k, l)*wrepp + ij_matrix%dd(10, k, l)*wrepp_d(1)
2368 793204 : v_d(1, ij, 35) = v_d(1, ij, 35) + ij_matrix%dd_d(1, 11, k, l)*wrepp + ij_matrix%dd(11, k, l)*wrepp_d(1)
2369 793204 : v_d(1, ij, 41) = v_d(1, ij, 41) + ij_matrix%dd_d(1, 12, k, l)*wrepp + ij_matrix%dd(12, k, l)*wrepp_d(1)
2370 793204 : v_d(1, ij, 42) = v_d(1, ij, 42) + ij_matrix%dd_d(1, 13, k, l)*wrepp + ij_matrix%dd(13, k, l)*wrepp_d(1)
2371 793204 : v_d(1, ij, 43) = v_d(1, ij, 43) + ij_matrix%dd_d(1, 14, k, l)*wrepp + ij_matrix%dd(14, k, l)*wrepp_d(1)
2372 793204 : v_d(1, ij, 44) = v_d(1, ij, 44) + ij_matrix%dd_d(1, 15, k, l)*wrepp + ij_matrix%dd(15, k, l)*wrepp_d(1)
2373 :
2374 793204 : v_d(2, ij, 15) = v_d(2, ij, 15) + ij_matrix%dd_d(2, 1, k, l)*wrepp + ij_matrix%dd(1, k, l)*wrepp_d(2)
2375 793204 : v_d(2, ij, 21) = v_d(2, ij, 21) + ij_matrix%dd_d(2, 2, k, l)*wrepp + ij_matrix%dd(2, k, l)*wrepp_d(2)
2376 793204 : v_d(2, ij, 28) = v_d(2, ij, 28) + ij_matrix%dd_d(2, 3, k, l)*wrepp + ij_matrix%dd(3, k, l)*wrepp_d(2)
2377 793204 : v_d(2, ij, 36) = v_d(2, ij, 36) + ij_matrix%dd_d(2, 4, k, l)*wrepp + ij_matrix%dd(4, k, l)*wrepp_d(2)
2378 793204 : v_d(2, ij, 45) = v_d(2, ij, 45) + ij_matrix%dd_d(2, 5, k, l)*wrepp + ij_matrix%dd(5, k, l)*wrepp_d(2)
2379 793204 : v_d(2, ij, 20) = v_d(2, ij, 20) + ij_matrix%dd_d(2, 6, k, l)*wrepp + ij_matrix%dd(6, k, l)*wrepp_d(2)
2380 793204 : v_d(2, ij, 26) = v_d(2, ij, 26) + ij_matrix%dd_d(2, 7, k, l)*wrepp + ij_matrix%dd(7, k, l)*wrepp_d(2)
2381 793204 : v_d(2, ij, 27) = v_d(2, ij, 27) + ij_matrix%dd_d(2, 8, k, l)*wrepp + ij_matrix%dd(8, k, l)*wrepp_d(2)
2382 793204 : v_d(2, ij, 33) = v_d(2, ij, 33) + ij_matrix%dd_d(2, 9, k, l)*wrepp + ij_matrix%dd(9, k, l)*wrepp_d(2)
2383 793204 : v_d(2, ij, 34) = v_d(2, ij, 34) + ij_matrix%dd_d(2, 10, k, l)*wrepp + ij_matrix%dd(10, k, l)*wrepp_d(2)
2384 793204 : v_d(2, ij, 35) = v_d(2, ij, 35) + ij_matrix%dd_d(2, 11, k, l)*wrepp + ij_matrix%dd(11, k, l)*wrepp_d(2)
2385 793204 : v_d(2, ij, 41) = v_d(2, ij, 41) + ij_matrix%dd_d(2, 12, k, l)*wrepp + ij_matrix%dd(12, k, l)*wrepp_d(2)
2386 793204 : v_d(2, ij, 42) = v_d(2, ij, 42) + ij_matrix%dd_d(2, 13, k, l)*wrepp + ij_matrix%dd(13, k, l)*wrepp_d(2)
2387 793204 : v_d(2, ij, 43) = v_d(2, ij, 43) + ij_matrix%dd_d(2, 14, k, l)*wrepp + ij_matrix%dd(14, k, l)*wrepp_d(2)
2388 793204 : v_d(2, ij, 44) = v_d(2, ij, 44) + ij_matrix%dd_d(2, 15, k, l)*wrepp + ij_matrix%dd(15, k, l)*wrepp_d(2)
2389 :
2390 793204 : v_d(3, ij, 15) = v_d(3, ij, 15) + ij_matrix%dd_d(3, 1, k, l)*wrepp + ij_matrix%dd(1, k, l)*wrepp_d(3)
2391 793204 : v_d(3, ij, 21) = v_d(3, ij, 21) + ij_matrix%dd_d(3, 2, k, l)*wrepp + ij_matrix%dd(2, k, l)*wrepp_d(3)
2392 793204 : v_d(3, ij, 28) = v_d(3, ij, 28) + ij_matrix%dd_d(3, 3, k, l)*wrepp + ij_matrix%dd(3, k, l)*wrepp_d(3)
2393 793204 : v_d(3, ij, 36) = v_d(3, ij, 36) + ij_matrix%dd_d(3, 4, k, l)*wrepp + ij_matrix%dd(4, k, l)*wrepp_d(3)
2394 793204 : v_d(3, ij, 45) = v_d(3, ij, 45) + ij_matrix%dd_d(3, 5, k, l)*wrepp + ij_matrix%dd(5, k, l)*wrepp_d(3)
2395 793204 : v_d(3, ij, 20) = v_d(3, ij, 20) + ij_matrix%dd_d(3, 6, k, l)*wrepp + ij_matrix%dd(6, k, l)*wrepp_d(3)
2396 793204 : v_d(3, ij, 26) = v_d(3, ij, 26) + ij_matrix%dd_d(3, 7, k, l)*wrepp + ij_matrix%dd(7, k, l)*wrepp_d(3)
2397 793204 : v_d(3, ij, 27) = v_d(3, ij, 27) + ij_matrix%dd_d(3, 8, k, l)*wrepp + ij_matrix%dd(8, k, l)*wrepp_d(3)
2398 793204 : v_d(3, ij, 33) = v_d(3, ij, 33) + ij_matrix%dd_d(3, 9, k, l)*wrepp + ij_matrix%dd(9, k, l)*wrepp_d(3)
2399 793204 : v_d(3, ij, 34) = v_d(3, ij, 34) + ij_matrix%dd_d(3, 10, k, l)*wrepp + ij_matrix%dd(10, k, l)*wrepp_d(3)
2400 793204 : v_d(3, ij, 35) = v_d(3, ij, 35) + ij_matrix%dd_d(3, 11, k, l)*wrepp + ij_matrix%dd(11, k, l)*wrepp_d(3)
2401 793204 : v_d(3, ij, 41) = v_d(3, ij, 41) + ij_matrix%dd_d(3, 12, k, l)*wrepp + ij_matrix%dd(12, k, l)*wrepp_d(3)
2402 793204 : v_d(3, ij, 42) = v_d(3, ij, 42) + ij_matrix%dd_d(3, 13, k, l)*wrepp + ij_matrix%dd(13, k, l)*wrepp_d(3)
2403 793204 : v_d(3, ij, 43) = v_d(3, ij, 43) + ij_matrix%dd_d(3, 14, k, l)*wrepp + ij_matrix%dd(14, k, l)*wrepp_d(3)
2404 793204 : v_d(3, ij, 44) = v_d(3, ij, 44) + ij_matrix%dd_d(3, 15, k, l)*wrepp + ij_matrix%dd(15, k, l)*wrepp_d(3)
2405 : END IF
2406 : END IF
2407 : END DO
2408 : END DO
2409 35715084 : DO kl = 1, limkl
2410 69054396 : logv_d(ij, kl) = (ANY(ABS(v_d(1:3, ij, kl)) > 0.0_dp))
2411 : END DO
2412 : END DO
2413 : END DO
2414 : IF (debug_this_module) THEN
2415 : CALL rot_2el_2c_first_debug(sepi, sepj, rijv, se_int_control, se_taper, invert, ii, kk, v_d)
2416 : END IF
2417 : END IF
2418 1195735 : END SUBROUTINE rot_2el_2c_first
2419 :
2420 : ! **************************************************************************************************
2421 : !> \brief Store the two-electron two-center integrals in diagonal form
2422 : !>
2423 : !> \param limij ...
2424 : !> \param limkl ...
2425 : !> \param ww ...
2426 : !> \param w ...
2427 : !> \param ww_dx ...
2428 : !> \param ww_dy ...
2429 : !> \param ww_dz ...
2430 : !> \param dw ...
2431 : !> \date 04.2008 [tlaino]
2432 : !> \author Teodoro Laino [tlaino] - University of Zurich
2433 : ! **************************************************************************************************
2434 3414082 : SUBROUTINE store_2el_2c_diag(limij, limkl, ww, w, ww_dx, ww_dy, ww_dz, dw)
2435 :
2436 : INTEGER, INTENT(IN) :: limij, limkl
2437 : REAL(KIND=dp), DIMENSION(limkl, limij), &
2438 : INTENT(IN), OPTIONAL :: ww
2439 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
2440 : OPTIONAL :: w
2441 : REAL(KIND=dp), DIMENSION(limkl, limij), &
2442 : INTENT(IN), OPTIONAL :: ww_dx, ww_dy, ww_dz
2443 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
2444 : OPTIONAL :: dw
2445 :
2446 : INTEGER :: ij, kl, l
2447 :
2448 1195735 : l = 0
2449 1195735 : IF (PRESENT(ww) .AND. PRESENT(w)) THEN
2450 5618602 : DO ij = 1, limij
2451 41727884 : DO kl = 1, limkl
2452 36109282 : l = l + 1
2453 41043455 : w(l) = ww(kl, ij)
2454 : END DO
2455 : END DO
2456 511306 : ELSE IF (PRESENT(ww_dx) .AND. PRESENT(ww_dy) .AND. PRESENT(ww_dz) .AND. PRESENT(dw)) THEN
2457 4314767 : DO ij = 1, limij
2458 34719459 : DO kl = 1, limkl
2459 30404692 : l = l + 1
2460 30404692 : dw(1, l) = ww_dx(kl, ij)
2461 30404692 : dw(2, l) = ww_dy(kl, ij)
2462 34208153 : dw(3, l) = ww_dz(kl, ij)
2463 : END DO
2464 : END DO
2465 : ELSE
2466 0 : CPABORT("")
2467 : END IF
2468 :
2469 3414082 : END SUBROUTINE store_2el_2c_diag
2470 :
2471 : END MODULE semi_empirical_int_utils
|
