(git:b279b6b)
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 ! **************************************************************************************************
12 MODULE semi_empirical_int_utils
13 
14  USE input_constants, ONLY: do_method_pchg, &
15  do_se_is_kdso_d
16  USE kinds, ONLY: dp
17  USE semi_empirical_int3_utils, ONLY: charg_int_3, &
18  dcharg_int_3, &
19  ijkl_low_3
20  USE semi_empirical_int_arrays, ONLY: &
21  clmp, clmxx, clmxy, clmyy, clmz, clmzp, clmzz, clm_d, clm_sp, ijkl_ind, indexa, indexb, &
22  int2c_type
23  USE semi_empirical_types, ONLY: rotmat_type, &
24  se_int_control_type, &
25  se_int_screen_type, &
26  se_taper_type, &
27  semi_empirical_type
28 #include "./base/base_uses.f90"
29 
30  IMPLICIT NONE
31  PRIVATE
32 
33 
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 ! **************************************************************************************************
40 INTERFACE
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
51 END 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 ! **************************************************************************************************
61 INTERFACE
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
75 END 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 ! **************************************************************************************************
85 INTERFACE
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
100 END 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 ! **************************************************************************************************
110 INTERFACE
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
126 END 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 ! **************************************************************************************************
136 INTERFACE
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 
153 END 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 ! **************************************************************************************************
163 INTERFACE
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
179 END 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 ! **************************************************************************************************
189 INTERFACE
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
203 END 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 ! **************************************************************************************************
213 INTERFACE
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
227 END 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 
282 CONTAINS
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 
2665 END MODULE semi_empirical_int_utils
semi_empirical_int_utils::check_dcore_nucint_ana
Check Numerical Vs Analytical NUCINT core.
Definition: semi_empirical_int_utils.F:86
semi_empirical_int_utils::check_dcorecore_ana
Check Numerical Vs Analytical CORECORE.
Definition: semi_empirical_int_utils.F:137
semi_empirical_int_utils::check_drotnuc_ana
Check Numerical Vs Analytical ROTNUC.
Definition: semi_empirical_int_utils.F:111
semi_empirical_int_utils::check_dssss_nucint_ana
Check Numerical Vs Analytical NUCINT ssss.
Definition: semi_empirical_int_utils.F:62
semi_empirical_int_utils::check_dterep_ana
Check Numerical Vs Analytical check_dterep_ana.
Definition: semi_empirical_int_utils.F:190
semi_empirical_int_utils::check_rotint_ana
Check Numerical Vs Analytical check_rotint_ana.
Definition: semi_empirical_int_utils.F:214
semi_empirical_int_utils::eval_func_d
...
Definition: semi_empirical_int_utils.F:272
semi_empirical_int_utils::rot_2el_2c_first_debug
Check Numerical Vs Analytical rot_2el_2c_first.
Definition: semi_empirical_int_utils.F:164
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