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