(git:b195825)
semi_empirical_int_ana.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 Analytical derivatives of Integrals for semi-empirical methods
10 !> \author Teodoro Laino - Zurich University 04.2007 [tlaino]
11 !> \par History
12 !> 23.11.2007 jhu short range version of integrals
13 !> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
14 !> for computing integrals
15 !> Teodoro Laino (05.2008) [tlaino] - University of Zurich : analytical
16 !> derivatives for d-orbitals
17 ! **************************************************************************************************
18 MODULE semi_empirical_int_ana
19 
20  USE input_constants, ONLY: do_method_am1, &
21  do_method_pchg, &
22  do_method_pdg, &
23  do_method_pm3, &
24  do_method_pm6, &
25  do_method_pm6fm, &
26  do_method_undef, &
27  do_se_is_kdso_d
28  USE kinds, ONLY: dp
29  USE multipole_types, ONLY: do_multipole_none
30  USE physcon, ONLY: angstrom, &
31  evolt
32  USE semi_empirical_int_arrays, ONLY: &
33  fac_x_to_z, ijkl_ind, ijkl_sym, inddd, inddp, indexa, indexb, indpp, int2c_type, l_index, &
34  map_x_to_z, rij_threshold
35  USE semi_empirical_int_num, ONLY: nucint_d_num, &
36  nucint_sp_num, &
37  terep_d_num, &
38  terep_sp_num
39  USE semi_empirical_int_utils, ONLY: d_ijkl_d, &
40  d_ijkl_sp, &
41  rot_2el_2c_first, &
42  rotmat, &
43  store_2el_2c_diag
44  USE semi_empirical_types, ONLY: rotmat_create, &
45  rotmat_release, &
46  rotmat_type, &
47  se_int_control_type, &
48  se_int_screen_type, &
49  se_taper_type, &
50  semi_empirical_type, &
51  setup_se_int_control_type
52  USE taper_types, ONLY: dtaper_eval, &
53  taper_eval
54 #include "./base/base_uses.f90"
55 
56  IMPLICIT NONE
57  PRIVATE
58 
59 
60 ! **************************************************************************************************
61 !> \brief Debug the derivatives of the the rotational matrices
62 !>
63 !> \author Teodoro Laino [tlaino] - University of Zurich
64 !> \date 04.2008 [tlaino]
65 ! **************************************************************************************************
66 INTERFACE
67  SUBROUTINE check_rotmat_der(sepi, sepj, rjiv, ij_matrix, do_invert)
68  USE kinds, ONLY: dp
69  USE semi_empirical_types, ONLY: rotmat_type, &
70  semi_empirical_type
71  TYPE(semi_empirical_type), POINTER :: sepi, sepj
72  REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rjiv
73  TYPE(rotmat_type), POINTER :: ij_matrix
74  LOGICAL, INTENT(IN) :: do_invert
75 
76  END SUBROUTINE check_rotmat_der
77 END INTERFACE
78 
79 ! **************************************************************************************************
80 !> \brief Check Numerical Vs Analytical NUCINT ssss
81 !> \note
82 !> Debug routine
83 !> \par History
84 !> 04.2008 created [tlaino]
85 !> \author Teodoro Laino - Zurich University
86 ! **************************************************************************************************
87 INTERFACE
88  SUBROUTINE check_dssss_nucint_ana(sepi, sepj, r, dssss, itype, se_int_control, &
89  se_taper)
90  USE kinds, ONLY: dp
91  USE semi_empirical_types, ONLY: semi_empirical_type, &
92  se_int_control_type, &
93  se_taper_type
94  TYPE(semi_empirical_type), POINTER :: sepi, sepj
95  REAL(dp), INTENT(IN) :: r, dssss
96  INTEGER, INTENT(IN) :: itype
97  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
98  TYPE(se_taper_type), POINTER :: se_taper
99 
100  END SUBROUTINE check_dssss_nucint_ana
101 END INTERFACE
102 
103 ! **************************************************************************************************
104 !> \brief Check Numerical Vs Analytical NUCINT core
105 !> \note
106 !> Debug routine
107 !> \par History
108 !> 04.2008 created [tlaino]
109 !> \author Teodoro Laino - Zurich University
110 ! **************************************************************************************************
111 INTERFACE
112  SUBROUTINE check_dcore_nucint_ana(sepi, sepj, r, dcore, itype, se_int_control, &
113  se_taper)
114  USE kinds, ONLY: dp
115  USE semi_empirical_types, ONLY: semi_empirical_type, &
116  se_int_control_type, &
117  se_taper_type
118  TYPE(semi_empirical_type), POINTER :: sepi, sepj
119  REAL(dp), INTENT(IN) :: r
120  REAL(dp), DIMENSION(10, 2), INTENT(IN) :: dcore
121  INTEGER, INTENT(IN) :: itype
122  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
123  TYPE(se_taper_type), POINTER :: se_taper
124 
125  END SUBROUTINE check_dcore_nucint_ana
126 END INTERFACE
127 
128 ! **************************************************************************************************
129 !> \brief Check Numerical Vs Analytical ROTNUC
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_drotnuc_ana(sepi, sepj, rijv, itype, se_int_control, se_taper, &
138  e1b, e2a, de1b, de2a)
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), DIMENSION(45), INTENT(IN), OPTIONAL :: e1b, e2a
149  REAL(dp), DIMENSION(45, 3), INTENT(IN), OPTIONAL :: de1b, de2a
150 
151  END SUBROUTINE check_drotnuc_ana
152 END INTERFACE
153 
154 ! **************************************************************************************************
155 !> \brief Check Numerical Vs Analytical CORECORE
156 !> \note
157 !> Debug routine
158 !> \par History
159 !> 04.2008 created [tlaino]
160 !> \author Teodoro Laino - Zurich University
161 ! **************************************************************************************************
162 INTERFACE
163  SUBROUTINE check_dcorecore_ana(sepi, sepj, rijv, itype, se_int_control, &
164  se_taper, enuc, denuc)
165  USE kinds, ONLY: dp
166  USE semi_empirical_types, ONLY: semi_empirical_type, &
167  se_int_control_type, &
168  se_taper_type
169  TYPE(semi_empirical_type), POINTER :: sepi, sepj
170  REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
171  INTEGER, INTENT(IN) :: itype
172  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
173  TYPE(se_taper_type), POINTER :: se_taper
174  REAL(dp), INTENT(IN), OPTIONAL :: enuc
175  REAL(dp), DIMENSION(3), INTENT(IN), OPTIONAL :: denuc
176 
177  END SUBROUTINE check_dcorecore_ana
178 
179 END INTERFACE
180 
181 ! **************************************************************************************************
182 !> \brief Check Numerical Vs Analytical rot_2el_2c_first
183 !> \note
184 !> Debug routine
185 !> \par History
186 !> 04.2008 created [tlaino]
187 !> \author Teodoro Laino - Zurich University
188 ! **************************************************************************************************
189 INTERFACE
190  SUBROUTINE rot_2el_2c_first_debug(sepi, sepj, rijv, se_int_control, se_taper, &
191  invert, ii, kk, v_d)
192  USE kinds, ONLY: dp
193  USE semi_empirical_types, ONLY: semi_empirical_type, &
194  se_int_control_type, &
195  se_taper_type
196  TYPE(semi_empirical_type), POINTER :: sepi, sepj
197  REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rijv
198  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
199  TYPE(se_taper_type), POINTER :: se_taper
200  LOGICAL, INTENT(IN) :: invert
201  INTEGER, INTENT(IN) :: ii, kk
202  REAL(KIND=dp), DIMENSION(45, 45, 3), INTENT(IN) :: v_d
203 
204  END SUBROUTINE rot_2el_2c_first_debug
205 END INTERFACE
206 
207 ! **************************************************************************************************
208 !> \brief Check Numerical Vs Analytical check_dterep_ana
209 !> \note
210 !> Debug routine
211 !> \par History
212 !> 04.2008 created [tlaino]
213 !> \author Teodoro Laino - Zurich University
214 ! **************************************************************************************************
215 INTERFACE
216  SUBROUTINE check_dterep_ana(sepi, sepj, r, ri, dri, se_int_control, se_taper, lgrad)
217  USE kinds, ONLY: dp
218  USE semi_empirical_types, ONLY: semi_empirical_type, &
219  se_int_control_type, &
220  se_taper_type
221  TYPE(semi_empirical_type), POINTER :: sepi, sepj
222  REAL(dp), INTENT(IN) :: r
223  REAL(dp), DIMENSION(491), INTENT(IN) :: ri, dri
224  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
225  TYPE(se_taper_type), POINTER :: se_taper
226  LOGICAL, INTENT(IN) :: lgrad
227 
228  END SUBROUTINE check_dterep_ana
229 END INTERFACE
230 
231 ! **************************************************************************************************
232 !> \brief Check Numerical Vs Analytical check_rotint_ana
233 !> \note
234 !> Debug routine
235 !> \par History
236 !> 04.2008 created [tlaino]
237 !> \author Teodoro Laino - Zurich University
238 ! **************************************************************************************************
239 INTERFACE
240  SUBROUTINE check_rotint_ana(sepi, sepj, rijv, w, dw, se_int_control, se_taper)
241  USE kinds, ONLY: dp
242  USE semi_empirical_types, ONLY: semi_empirical_type, &
243  se_int_control_type, &
244  se_taper_type
245  TYPE(semi_empirical_type), POINTER :: sepi, sepj
246  REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
247  REAL(dp), DIMENSION(2025), INTENT(IN), OPTIONAL :: w
248  REAL(dp), DIMENSION(2025, 3), INTENT(IN), OPTIONAL :: dw
249  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
250  TYPE(se_taper_type), POINTER :: se_taper
251 
252  END SUBROUTINE check_rotint_ana
253 END INTERFACE
254 
255  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_int_ana'
256  LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
257  PUBLIC :: rotnuc_ana, rotint_ana, corecore_ana, corecore_el_ana
258 
259 CONTAINS
260 
261 ! **************************************************************************************************
262 !> \brief Computes analytical gradients for semiempirical integrals
263 !> \param sepi Atomic parameters of first atom
264 !> \param sepj Atomic parameters of second atom
265 !> \param rijv Coordinate vector i -> j
266 !> \param itype ...
267 !> \param e1b Array of electron-nuclear attraction integrals, Electron on atom ni attracting nucleus of nj.
268 !> \param e2a Array of electron-nuclear attraction integrals, Electron on atom nj attracting nucleus of ni.
269 !> \param de1b derivative of e1b term
270 !> \param de2a derivative of e2a term
271 !> \param se_int_control input parameters that control the calculation of SE
272 !> integrals (shortrange, R3 residual, screening type)
273 !> \param se_taper ...
274 !> \par History
275 !> 04.2007 created [tlaino]
276 !> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
277 !> for computing integrals
278 !> Teodoro Laino [tlaino] - University of Zurich 04.2008 : removed the core-core part
279 !> \author Teodoro Laino [tlaino] - Zurich University
280 !> \note
281 !> Analytical version of the MOPAC rotnuc routine
282 ! **************************************************************************************************
283  RECURSIVE SUBROUTINE rotnuc_ana(sepi, sepj, rijv, itype, e1b, e2a, de1b, de2a, &
284  se_int_control, se_taper)
285  TYPE(semi_empirical_type), POINTER :: sepi, sepj
286  REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
287  INTEGER, INTENT(IN) :: itype
288  REAL(dp), DIMENSION(45), INTENT(OUT), OPTIONAL :: e1b, e2a
289  REAL(dp), DIMENSION(3, 45), INTENT(OUT), OPTIONAL :: de1b, de2a
290  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
291  TYPE(se_taper_type), POINTER :: se_taper
292 
293  INTEGER :: i, idd, idp, ind1, ind2, ipp, j, &
294  last_orbital(2), m, n
295  LOGICAL :: invert, l_de1b, l_de2a, l_e1b, l_e2a, &
296  lgrad, task(2)
297  REAL(kind=dp) :: rij, xtmp
298  REAL(kind=dp), DIMENSION(10, 2) :: core, dcore
299  REAL(kind=dp), DIMENSION(3) :: drij
300  REAL(kind=dp), DIMENSION(3, 45) :: tmp_d
301  REAL(kind=dp), DIMENSION(45) :: tmp
302  TYPE(rotmat_type), POINTER :: ij_matrix
303 
304  NULLIFY (ij_matrix)
305  rij = dot_product(rijv, rijv)
306  ! Initialization
307  l_e1b = PRESENT(e1b)
308  l_e2a = PRESENT(e2a)
309  l_de1b = PRESENT(de1b)
310  l_de2a = PRESENT(de2a)
311  lgrad = l_de1b .OR. l_de2a
312 
313  IF (rij > rij_threshold) THEN
314  ! Compute Integrals in diatomic frame opportunely inverted
315  rij = sqrt(rij)
316  ! Create the rotation matrix
317  CALL rotmat_create(ij_matrix)
318  CALL rotmat(sepi, sepj, rijv, rij, ij_matrix, do_derivatives=lgrad, do_invert=invert)
319 
320  IF (lgrad) THEN
321  drij(1) = rijv(1)/rij
322  drij(2) = rijv(2)/rij
323  drij(3) = rijv(3)/rij
324  ! Possibly Invert Frame
325  IF (invert) THEN
326  xtmp = drij(3)
327  drij(3) = drij(1)
328  drij(1) = xtmp
329  END IF
330  END IF
331 
332  CALL dcore_nucint_ana(sepi, sepj, rij, core=core, dcore=dcore, itype=itype, se_taper=se_taper, &
333  se_int_control=se_int_control, lgrad=lgrad)
334 
335  ! Copy parameters over to arrays for do loop.
336  last_orbital(1) = sepi%natorb
337  last_orbital(2) = sepj%natorb
338  task(1) = l_e1b
339  task(2) = l_e2a
340  DO n = 1, 2
341  IF (.NOT. task(n)) cycle
342  DO i = 1, last_orbital(n)
343  ind1 = i - 1
344  DO j = 1, i
345  ind2 = j - 1
346  m = (i*(i - 1))/2 + j
347  ! Perform Rotations ...
348  IF (ind2 == 0) THEN
349  IF (ind1 == 0) THEN
350  ! Type of Integral (SS/)
351  tmp(m) = core(1, n)
352  ELSE IF (ind1 < 4) THEN
353  ! Type of Integral (SP/)
354  tmp(m) = ij_matrix%sp(1, ind1)*core(2, n)
355  ELSE
356  ! Type of Integral (SD/)
357  tmp(m) = ij_matrix%sd(1, ind1 - 3)*core(5, n)
358  END IF
359  ELSE IF (ind2 < 4) THEN
360  IF (ind1 < 4) THEN
361  ! Type of Integral (PP/)
362  ipp = indpp(ind1, ind2)
363  tmp(m) = core(3, n)*ij_matrix%pp(ipp, 1, 1) + &
364  core(4, n)*(ij_matrix%pp(ipp, 2, 2) + ij_matrix%pp(ipp, 3, 3))
365  ELSE
366  ! Type of Integral (PD/)
367  idp = inddp(ind1 - 3, ind2)
368  tmp(m) = core(6, n)*ij_matrix%pd(idp, 1, 1) + &
369  core(8, n)*(ij_matrix%pd(idp, 2, 2) + ij_matrix%pd(idp, 3, 3))
370  END IF
371  ELSE
372  ! Type of Integral (DD/)
373  idd = inddd(ind1 - 3, ind2 - 3)
374  tmp(m) = core(7, n)*ij_matrix%dd(idd, 1, 1) + &
375  core(9, n)*(ij_matrix%dd(idd, 2, 2) + ij_matrix%dd(idd, 3, 3)) + &
376  core(10, n)*(ij_matrix%dd(idd, 4, 4) + ij_matrix%dd(idd, 5, 5))
377  END IF
378  END DO
379  END DO
380  IF (n == 1) THEN
381  DO i = 1, sepi%atm_int_size
382  e1b(i) = -tmp(i)
383  END DO
384  END IF
385  IF (n == 2) THEN
386  DO i = 1, sepj%atm_int_size
387  e2a(i) = -tmp(i)
388  END DO
389  END IF
390  END DO
391  IF (invert .AND. l_e1b) CALL invert_integral(sepi, sepi, int1el=e1b)
392  IF (invert .AND. l_e2a) CALL invert_integral(sepj, sepj, int1el=e2a)
393 
394  ! Possibly compute derivatives
395  task(1) = l_de1b
396  task(2) = l_de2a
397  DO n = 1, 2
398  IF (.NOT. task(n)) cycle
399  DO i = 1, last_orbital(n)
400  ind1 = i - 1
401  DO j = 1, i
402  ind2 = j - 1
403  m = (i*(i - 1))/2 + j
404  ! Perform Rotations ...
405  IF (ind2 == 0) THEN
406  IF (ind1 == 0) THEN
407  ! Type of Integral (SS/)
408  tmp_d(1, m) = dcore(1, n)*drij(1)
409  tmp_d(2, m) = dcore(1, n)*drij(2)
410  tmp_d(3, m) = dcore(1, n)*drij(3)
411  ELSE IF (ind1 < 4) THEN
412  ! Type of Integral (SP/)
413  tmp_d(1, m) = ij_matrix%sp_d(1, 1, ind1)*core(2, n) + &
414  ij_matrix%sp(1, ind1)*dcore(2, n)*drij(1)
415 
416  tmp_d(2, m) = ij_matrix%sp_d(2, 1, ind1)*core(2, n) + &
417  ij_matrix%sp(1, ind1)*dcore(2, n)*drij(2)
418 
419  tmp_d(3, m) = ij_matrix%sp_d(3, 1, ind1)*core(2, n) + &
420  ij_matrix%sp(1, ind1)*dcore(2, n)*drij(3)
421  ELSE
422  ! Type of Integral (SD/)
423  tmp_d(1, m) = ij_matrix%sd_d(1, 1, ind1 - 3)*core(5, n) + &
424  ij_matrix%sd(1, ind1 - 3)*dcore(5, n)*drij(1)
425 
426  tmp_d(2, m) = ij_matrix%sd_d(2, 1, ind1 - 3)*core(5, n) + &
427  ij_matrix%sd(1, ind1 - 3)*dcore(5, n)*drij(2)
428 
429  tmp_d(3, m) = ij_matrix%sd_d(3, 1, ind1 - 3)*core(5, n) + &
430  ij_matrix%sd(1, ind1 - 3)*dcore(5, n)*drij(3)
431  END IF
432  ELSE IF (ind2 < 4) THEN
433  IF (ind1 < 4) THEN
434  ! Type of Integral (PP/)
435  ipp = indpp(ind1, ind2)
436  tmp_d(1, m) = dcore(3, n)*drij(1)*ij_matrix%pp(ipp, 1, 1) + &
437  core(3, n)*ij_matrix%pp_d(1, ipp, 1, 1) + &
438  dcore(4, n)*drij(1)*(ij_matrix%pp(ipp, 2, 2) + ij_matrix%pp(ipp, 3, 3)) + &
439  core(4, n)*(ij_matrix%pp_d(1, ipp, 2, 2) + ij_matrix%pp_d(1, ipp, 3, 3))
440 
441  tmp_d(2, m) = dcore(3, n)*drij(2)*ij_matrix%pp(ipp, 1, 1) + &
442  core(3, n)*ij_matrix%pp_d(2, ipp, 1, 1) + &
443  dcore(4, n)*drij(2)*(ij_matrix%pp(ipp, 2, 2) + ij_matrix%pp(ipp, 3, 3)) + &
444  core(4, n)*(ij_matrix%pp_d(2, ipp, 2, 2) + ij_matrix%pp_d(2, ipp, 3, 3))
445 
446  tmp_d(3, m) = dcore(3, n)*drij(3)*ij_matrix%pp(ipp, 1, 1) + &
447  core(3, n)*ij_matrix%pp_d(3, ipp, 1, 1) + &
448  dcore(4, n)*drij(3)*(ij_matrix%pp(ipp, 2, 2) + ij_matrix%pp(ipp, 3, 3)) + &
449  core(4, n)*(ij_matrix%pp_d(3, ipp, 2, 2) + ij_matrix%pp_d(3, ipp, 3, 3))
450  ELSE
451  ! Type of Integral (PD/)
452  idp = inddp(ind1 - 3, ind2)
453  tmp_d(1, m) = dcore(6, n)*drij(1)*ij_matrix%pd(idp, 1, 1) + &
454  core(6, n)*ij_matrix%pd_d(1, idp, 1, 1) + &
455  dcore(8, n)*drij(1)*(ij_matrix%pd(idp, 2, 2) + ij_matrix%pd(idp, 3, 3)) + &
456  core(8, n)*(ij_matrix%pd_d(1, idp, 2, 2) + ij_matrix%pd_d(1, idp, 3, 3))
457 
458  tmp_d(2, m) = dcore(6, n)*drij(2)*ij_matrix%pd(idp, 1, 1) + &
459  core(6, n)*ij_matrix%pd_d(2, idp, 1, 1) + &
460  dcore(8, n)*drij(2)*(ij_matrix%pd(idp, 2, 2) + ij_matrix%pd(idp, 3, 3)) + &
461  core(8, n)*(ij_matrix%pd_d(2, idp, 2, 2) + ij_matrix%pd_d(2, idp, 3, 3))
462 
463  tmp_d(3, m) = dcore(6, n)*drij(3)*ij_matrix%pd(idp, 1, 1) + &
464  core(6, n)*ij_matrix%pd_d(3, idp, 1, 1) + &
465  dcore(8, n)*drij(3)*(ij_matrix%pd(idp, 2, 2) + ij_matrix%pd(idp, 3, 3)) + &
466  core(8, n)*(ij_matrix%pd_d(3, idp, 2, 2) + ij_matrix%pd_d(3, idp, 3, 3))
467  END IF
468  ELSE
469  ! Type of Integral (DD/)
470  idd = inddd(ind1 - 3, ind2 - 3)
471  tmp_d(1, m) = dcore(7, n)*drij(1)*ij_matrix%dd(idd, 1, 1) + &
472  core(7, n)*ij_matrix%dd_d(1, idd, 1, 1) + &
473  dcore(9, n)*drij(1)*(ij_matrix%dd(idd, 2, 2) + ij_matrix%dd(idd, 3, 3)) + &
474  core(9, n)*(ij_matrix%dd_d(1, idd, 2, 2) + ij_matrix%dd_d(1, idd, 3, 3)) + &
475  dcore(10, n)*drij(1)*(ij_matrix%dd(idd, 4, 4) + ij_matrix%dd(idd, 5, 5)) + &
476  core(10, n)*(ij_matrix%dd_d(1, idd, 4, 4) + ij_matrix%dd_d(1, idd, 5, 5))
477 
478  tmp_d(2, m) = dcore(7, n)*drij(2)*ij_matrix%dd(idd, 1, 1) + &
479  core(7, n)*ij_matrix%dd_d(2, idd, 1, 1) + &
480  dcore(9, n)*drij(2)*(ij_matrix%dd(idd, 2, 2) + ij_matrix%dd(idd, 3, 3)) + &
481  core(9, n)*(ij_matrix%dd_d(2, idd, 2, 2) + ij_matrix%dd_d(2, idd, 3, 3)) + &
482  dcore(10, n)*drij(2)*(ij_matrix%dd(idd, 4, 4) + ij_matrix%dd(idd, 5, 5)) + &
483  core(10, n)*(ij_matrix%dd_d(2, idd, 4, 4) + ij_matrix%dd_d(2, idd, 5, 5))
484 
485  tmp_d(3, m) = dcore(7, n)*drij(3)*ij_matrix%dd(idd, 1, 1) + &
486  core(7, n)*ij_matrix%dd_d(3, idd, 1, 1) + &
487  dcore(9, n)*drij(3)*(ij_matrix%dd(idd, 2, 2) + ij_matrix%dd(idd, 3, 3)) + &
488  core(9, n)*(ij_matrix%dd_d(3, idd, 2, 2) + ij_matrix%dd_d(3, idd, 3, 3)) + &
489  dcore(10, n)*drij(3)*(ij_matrix%dd(idd, 4, 4) + ij_matrix%dd(idd, 5, 5)) + &
490  core(10, n)*(ij_matrix%dd_d(3, idd, 4, 4) + ij_matrix%dd_d(3, idd, 5, 5))
491  END IF
492  END DO
493  END DO
494  IF (n == 1) THEN
495  DO i = 1, sepi%atm_int_size
496  de1b(1, i) = -tmp_d(1, i)
497  de1b(2, i) = -tmp_d(2, i)
498  de1b(3, i) = -tmp_d(3, i)
499  END DO
500  END IF
501  IF (n == 2) THEN
502  DO i = 1, sepj%atm_int_size
503  de2a(1, i) = -tmp_d(1, i)
504  de2a(2, i) = -tmp_d(2, i)
505  de2a(3, i) = -tmp_d(3, i)
506  END DO
507  END IF
508  END DO
509  IF (invert .AND. l_de1b) CALL invert_derivative(sepi, sepi, dint1el=de1b)
510  IF (invert .AND. l_de2a) CALL invert_derivative(sepj, sepj, dint1el=de2a)
511  CALL rotmat_release(ij_matrix)
512 
513  ! Possibly debug the analytical values versus the numerical ones
514  IF (debug_this_module) THEN
515  CALL check_drotnuc_ana(sepi, sepj, rijv, itype, se_int_control, se_taper, e1b, e2a, de1b, de2a)
516  END IF
517  END IF
518  END SUBROUTINE rotnuc_ana
519 
520 ! **************************************************************************************************
521 !> \brief Computes analytical gradients for semiempirical core-core interaction.
522 !> \param sepi Atomic parameters of first atom
523 !> \param sepj Atomic parameters of second atom
524 !> \param rijv Coordinate vector i -> j
525 !> \param itype ...
526 !> \param enuc nuclear-nuclear repulsion term.
527 !> \param denuc derivative of nuclear-nuclear repulsion term.
528 !> \param se_int_control input parameters that control the calculation of SE
529 !> integrals (shortrange, R3 residual, screening type)
530 !> \param se_taper ...
531 !> \par History
532 !> 04.2007 created [tlaino]
533 !> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
534 !> for computing integrals
535 !> Teodoro Laino [tlaino] - University of Zurich 04.2008 : removed the
536 !> core-core part
537 !> \author Teodoro Laino [tlaino] - Zurich University
538 !> \note
539 !> Analytical version of the MOPAC rotnuc routine
540 ! **************************************************************************************************
541  RECURSIVE SUBROUTINE corecore_ana(sepi, sepj, rijv, itype, enuc, denuc, se_int_control, se_taper)
542  TYPE(semi_empirical_type), POINTER :: sepi, sepj
543  REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
544  INTEGER, INTENT(IN) :: itype
545  REAL(dp), INTENT(OUT), OPTIONAL :: enuc
546  REAL(dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: denuc
547  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
548  TYPE(se_taper_type), POINTER :: se_taper
549 
550  INTEGER :: ig, nt
551  LOGICAL :: l_denuc, l_enuc
552  REAL(dp) :: aab, alpi, alpj, apdg, ax, dai, daj, dax, dbi, dbj, denuc_loc, dqcorr, drija, &
553  dscale, dssss, dssss_sr, dtmp, dzz, enuc_loc, pai, paj, pbi, pbj, qcorr, rij, rija, &
554  scale, ssss, ssss_sr, tmp, xab, xtmp, zaf, zbf, zz
555  REAL(dp), DIMENSION(3) :: drij
556  REAL(dp), DIMENSION(4) :: fni1, fni2, fni3, fnj1, fnj2, fnj3
557  TYPE(se_int_control_type) :: se_int_control_off
558 
559  rij = dot_product(rijv, rijv)
560  ! Initialization
561  l_enuc = PRESENT(enuc)
562  l_denuc = PRESENT(denuc)
563  IF ((rij > rij_threshold) .AND. (l_enuc .OR. l_denuc)) THEN
564  ! Compute Integrals in diatomic frame
565  rij = sqrt(rij)
566  CALL setup_se_int_control_type(se_int_control_off, shortrange=.false., do_ewald_r3=.false., &
567  do_ewald_gks=.false., integral_screening=se_int_control%integral_screening, &
568  max_multipole=do_multipole_none, pc_coulomb_int=.false.)
569  CALL dssss_nucint_ana(sepi, sepj, rij, ssss=ssss, dssss=dssss, itype=itype, se_taper=se_taper, &
570  se_int_control=se_int_control_off, lgrad=l_denuc)
571  ! In case let's compute the short-range part of the (ss|ss) integral
572  IF (se_int_control%shortrange) THEN
573  CALL dssss_nucint_ana(sepi, sepj, rij, ssss=ssss_sr, dssss=dssss_sr, itype=itype, &
574  se_taper=se_taper, se_int_control=se_int_control, lgrad=l_denuc)
575  ELSE
576  ssss_sr = ssss
577  dssss_sr = dssss
578  END IF
579  ! Zeroing local method dependent core-core corrections
580  enuc_loc = 0.0_dp
581  denuc_loc = 0.0_dp
582  qcorr = 0.0_dp
583  scale = 0.0_dp
584  dscale = 0.0_dp
585  dqcorr = 0.0_dp
586  zz = sepi%zeff*sepj%zeff
587  ! Core Core electrostatic contribution
588  IF (l_enuc) enuc_loc = zz*ssss_sr
589  IF (l_denuc) denuc_loc = zz*dssss_sr
590  ! Method dependent code
591  tmp = zz*ssss
592  IF (l_denuc) dtmp = zz*dssss
593  IF (itype /= do_method_pm6 .AND. itype /= do_method_pm6fm) THEN
594  alpi = sepi%alp
595  alpj = sepj%alp
596  scale = exp(-alpi*rij) + exp(-alpj*rij)
597  IF (l_denuc) THEN
598  dscale = -alpi*exp(-alpi*rij) - alpj*exp(-alpj*rij)
599  END IF
600  nt = sepi%z + sepj%z
601  IF (nt == 8 .OR. nt == 9) THEN
602  IF (sepi%z == 7 .OR. sepi%z == 8) THEN
603  scale = scale + (angstrom*rij - 1._dp)*exp(-alpi*rij)
604  IF (l_denuc) THEN
605  dscale = dscale + angstrom*exp(-alpi*rij) - (angstrom*rij - 1._dp)*alpi*exp(-alpi*rij)
606  END IF
607  END IF
608  IF (sepj%z == 7 .OR. sepj%z == 8) THEN
609  scale = scale + (angstrom*rij - 1._dp)*exp(-alpj*rij)
610  IF (l_denuc) THEN
611  dscale = dscale + angstrom*exp(-alpj*rij) - (angstrom*rij - 1._dp)*alpj*exp(-alpj*rij)
612  END IF
613  END IF
614  END IF
615  IF (l_denuc) THEN
616  dscale = sign(1.0_dp, scale*tmp)*(dscale*tmp + scale*dtmp)
617  dzz = -zz/rij**2
618  END IF
619  scale = abs(scale*tmp)
620  zz = zz/rij
621  IF (itype == do_method_am1 .OR. itype == do_method_pm3 .OR. itype == do_method_pdg) THEN
622  IF (itype == do_method_am1 .AND. sepi%z == 5) THEN
623  !special case AM1 Boron
624  SELECT CASE (sepj%z)
625  CASE DEFAULT
626  nt = 1
627  CASE (1)
628  nt = 2
629  CASE (6)
630  nt = 3
631  CASE (9, 17, 35, 53)
632  nt = 4
633  END SELECT
634  fni1(1) = sepi%bfn1(1, nt)
635  fni1(2) = sepi%bfn1(2, nt)
636  fni1(3) = sepi%bfn1(3, nt)
637  fni1(4) = sepi%bfn1(4, nt)
638  fni2(1) = sepi%bfn2(1, nt)
639  fni2(2) = sepi%bfn2(2, nt)
640  fni2(3) = sepi%bfn2(3, nt)
641  fni2(4) = sepi%bfn2(4, nt)
642  fni3(1) = sepi%bfn3(1, nt)
643  fni3(2) = sepi%bfn3(2, nt)
644  fni3(3) = sepi%bfn3(3, nt)
645  fni3(4) = sepi%bfn3(4, nt)
646  ELSE
647  fni1(1) = sepi%fn1(1)
648  fni1(2) = sepi%fn1(2)
649  fni1(3) = sepi%fn1(3)
650  fni1(4) = sepi%fn1(4)
651  fni2(1) = sepi%fn2(1)
652  fni2(2) = sepi%fn2(2)
653  fni2(3) = sepi%fn2(3)
654  fni2(4) = sepi%fn2(4)
655  fni3(1) = sepi%fn3(1)
656  fni3(2) = sepi%fn3(2)
657  fni3(3) = sepi%fn3(3)
658  fni3(4) = sepi%fn3(4)
659  END IF
660  IF (itype == do_method_am1 .AND. sepj%z == 5) THEN
661  !special case AM1 Boron
662  SELECT CASE (sepi%z)
663  CASE DEFAULT
664  nt = 1
665  CASE (1)
666  nt = 2
667  CASE (6)
668  nt = 3
669  CASE (9, 17, 35, 53)
670  nt = 4
671  END SELECT
672  fnj1(1) = sepj%bfn1(1, nt)
673  fnj1(2) = sepj%bfn1(2, nt)
674  fnj1(3) = sepj%bfn1(3, nt)
675  fnj1(4) = sepj%bfn1(4, nt)
676  fnj2(1) = sepj%bfn2(1, nt)
677  fnj2(2) = sepj%bfn2(2, nt)
678  fnj2(3) = sepj%bfn2(3, nt)
679  fnj2(4) = sepj%bfn2(4, nt)
680  fnj3(1) = sepj%bfn3(1, nt)
681  fnj3(2) = sepj%bfn3(2, nt)
682  fnj3(3) = sepj%bfn3(3, nt)
683  fnj3(4) = sepj%bfn3(4, nt)
684  ELSE
685  fnj1(1) = sepj%fn1(1)
686  fnj1(2) = sepj%fn1(2)
687  fnj1(3) = sepj%fn1(3)
688  fnj1(4) = sepj%fn1(4)
689  fnj2(1) = sepj%fn2(1)
690  fnj2(2) = sepj%fn2(2)
691  fnj2(3) = sepj%fn2(3)
692  fnj2(4) = sepj%fn2(4)
693  fnj3(1) = sepj%fn3(1)
694  fnj3(2) = sepj%fn3(2)
695  fnj3(3) = sepj%fn3(3)
696  fnj3(4) = sepj%fn3(4)
697  END IF
698  ! AM1/PM3/PDG correction to nuclear repulsion
699  DO ig = 1, SIZE(fni1)
700  IF (abs(fni1(ig)) > 0._dp) THEN
701  ax = fni2(ig)*(rij - fni3(ig))**2
702  IF (ax <= 25._dp) THEN
703  scale = scale + zz*fni1(ig)*exp(-ax)
704  IF (l_denuc) THEN
705  dax = fni2(ig)*2.0_dp*(rij - fni3(ig))
706  dscale = dscale + dzz*fni1(ig)*exp(-ax) - dax*zz*fni1(ig)*exp(-ax)
707  END IF
708  END IF
709  END IF
710  IF (abs(fnj1(ig)) > 0._dp) THEN
711  ax = fnj2(ig)*(rij - fnj3(ig))**2
712  IF (ax <= 25._dp) THEN
713  scale = scale + zz*fnj1(ig)*exp(-ax)
714  IF (l_denuc) THEN
715  dax = fnj2(ig)*2.0_dp*(rij - fnj3(ig))
716  dscale = dscale + dzz*fnj1(ig)*exp(-ax) - dax*zz*fnj1(ig)*exp(-ax)
717  END IF
718  END IF
719  END IF
720  END DO
721  END IF
722  IF (itype == do_method_pdg) THEN
723  ! PDDG function
724  zaf = sepi%zeff/nt
725  zbf = sepj%zeff/nt
726  pai = sepi%pre(1)
727  pbi = sepi%pre(2)
728  paj = sepj%pre(1)
729  pbj = sepj%pre(2)
730  dai = sepi%d(1)
731  dbi = sepi%d(2)
732  daj = sepj%d(1)
733  dbj = sepj%d(2)
734  apdg = 10._dp*angstrom**2
735  qcorr = (zaf*pai + zbf*paj)*exp(-apdg*(rij - dai - daj)**2) + &
736  (zaf*pai + zbf*pbj)*exp(-apdg*(rij - dai - dbj)**2) + &
737  (zaf*pbi + zbf*paj)*exp(-apdg*(rij - dbi - daj)**2) + &
738  (zaf*pbi + zbf*pbj)*exp(-apdg*(rij - dbi - dbj)**2)
739  IF (l_denuc) THEN
740  dqcorr = (zaf*pai + zbf*paj)*exp(-apdg*(rij - dai - daj)**2)*(-2.0_dp*apdg*(rij - dai - daj)) + &
741  (zaf*pai + zbf*pbj)*exp(-apdg*(rij - dai - dbj)**2)*(-2.0_dp*apdg*(rij - dai - dbj)) + &
742  (zaf*pbi + zbf*paj)*exp(-apdg*(rij - dbi - daj)**2)*(-2.0_dp*apdg*(rij - dbi - daj)) + &
743  (zaf*pbi + zbf*pbj)*exp(-apdg*(rij - dbi - dbj)**2)*(-2.0_dp*apdg*(rij - dbi - dbj))
744  END IF
745  ELSEIF (itype == do_method_pchg) THEN
746  qcorr = 0.0_dp
747  scale = 0.0_dp
748  dscale = 0.0_dp
749  dqcorr = 0.0_dp
750  ELSE
751  qcorr = 0.0_dp
752  dqcorr = 0.0_dp
753  END IF
754  ELSE
755  ! PM6 core-core terms
756  scale = tmp
757  IF (l_denuc) dscale = dtmp
758  drija = angstrom
759  rija = rij*drija
760  xab = sepi%xab(sepj%z)
761  aab = sepi%aab(sepj%z)
762  IF ((sepi%z == 1 .AND. (sepj%z == 6 .OR. sepj%z == 7 .OR. sepj%z == 8)) .OR. &
763  (sepj%z == 1 .AND. (sepi%z == 6 .OR. sepi%z == 7 .OR. sepi%z == 8))) THEN
764  ! Special Case O-H or N-H or C-H
765  IF (l_denuc) dscale = dscale*(2._dp*xab*exp(-aab*rija*rija)) - &
766  scale*2._dp*xab*exp(-aab*rija*rija)*(2.0_dp*aab*rija)*drija
767  IF (l_enuc) scale = scale*(2._dp*xab*exp(-aab*rija*rija))
768  ELSEIF (sepi%z == 6 .AND. sepj%z == 6) THEN
769  ! Special Case C-C
770  IF (l_denuc) dscale = &
771  dscale*(2._dp*xab*exp(-aab*(rija + 0.0003_dp*rija**6)) + 9.28_dp*exp(-5.98_dp*rija)) &
772  - scale*2._dp*xab*exp(-aab*(rija + 0.0003_dp*rija**6))*aab*(1.0_dp + 6.0_dp*0.0003_dp*rija**5)*drija &
773  - scale*9.28_dp*exp(-5.98_dp*rija)*5.98_dp*drija
774  IF (l_enuc) scale = scale*(2._dp*xab*exp(-aab*(rija + 0.0003_dp*rija**6)) + 9.28_dp*exp(-5.98_dp*rija))
775  ELSEIF ((sepi%z == 8 .AND. sepj%z == 14) .OR. &
776  (sepj%z == 8 .AND. sepi%z == 14)) THEN
777  ! Special Case Si-O
778  IF (l_denuc) dscale = &
779  dscale*(2._dp*xab*exp(-aab*(rija + 0.0003_dp*rija**6)) - 0.0007_dp*exp(-(rija - 2.9_dp)**2)) &
780  - scale*2._dp*xab*exp(-aab*(rija + 0.0003_dp*rija**6))*aab*(1.0_dp + 6.0_dp*0.0003_dp*rija**5)*drija + &
781  scale*0.0007_dp*exp(-(rija - 2.9_dp)**2)*(2.0_dp*(rija - 2.9_dp)*drija)
782  IF (l_enuc) scale = scale*(2._dp*xab*exp(-aab*(rija + 0.0003_dp*rija**6)) - 0.0007_dp*exp(-(rija - 2.9_dp)**2))
783  ELSE
784  ! General Case
785  ! Factor of 2 found by experiment
786  IF (l_denuc) dscale = dscale*(2._dp*xab*exp(-aab*(rija + 0.0003_dp*rija**6))) &
787  - scale*2._dp*xab*exp(-aab*(rija + 0.0003_dp*rija**6))*aab*(1.0_dp + 6.0_dp*0.0003_dp*rija**5)*drija
788  IF (l_enuc) scale = scale*(2._dp*xab*exp(-aab*(rija + 0.0003_dp*rija**6)))
789  END IF
790  ! General correction term a*exp(-b*(rij-c)^2)
791  xtmp = 1.e-8_dp/evolt*((real(sepi%z, dp)**(1._dp/3._dp) + real(sepj%z, dp)**(1._dp/3._dp))/rija)**12
792  IF (l_enuc) THEN
793  qcorr = (sepi%a*exp(-sepi%b*(rij - sepi%c)**2))*zz/rij + &
794  (sepj%a*exp(-sepj%b*(rij - sepj%c)**2))*zz/rij + &
795  ! Hard core repulsion
796  xtmp
797  END IF
798  IF (l_denuc) THEN
799  dqcorr = (sepi%a*exp(-sepi%b*(rij - sepi%c)**2)*(-2.0_dp*sepi%b*(rij - sepi%c)))*zz/rij - &
800  (sepi%a*exp(-sepi%b*(rij - sepi%c)**2))*zz/rij**2 + &
801  (sepj%a*exp(-sepj%b*(rij - sepj%c)**2)*(-2.0_dp*sepj%b*(rij - sepj%c)))*zz/rij - &
802  (sepj%a*exp(-sepj%b*(rij - sepj%c)**2))*zz/rij**2 + &
803  ! Hard core repulsion
804  (-12.0_dp*xtmp/rija*drija)
805  END IF
806  END IF
807 
808  ! Only at the very end let's sum-up the several contributions energy/derivatives
809  ! This assignment should be method independent
810  IF (l_enuc) THEN
811  enuc = enuc_loc + scale + qcorr
812  END IF
813  IF (l_denuc) THEN
814  drij(1) = rijv(1)/rij
815  drij(2) = rijv(2)/rij
816  drij(3) = rijv(3)/rij
817  denuc = (denuc_loc + dscale + dqcorr)*drij
818  END IF
819  ! Debug statement
820  IF (debug_this_module) THEN
821  CALL check_dcorecore_ana(sepi, sepj, rijv, itype, se_int_control, se_taper, enuc, denuc)
822  END IF
823  END IF
824  END SUBROUTINE corecore_ana
825 
826 ! **************************************************************************************************
827 !> \brief Computes analytical gradients for semiempirical core-core electrostatic
828 !> interaction only.
829 !> \param sepi Atomic parameters of first atom
830 !> \param sepj Atomic parameters of second atom
831 !> \param rijv Coordinate vector i -> j
832 !> \param itype ...
833 !> \param enuc nuclear-nuclear electrostatic repulsion term.
834 !> \param denuc derivative of nuclear-nuclear electrostatic
835 !> repulsion term.
836 !> \param se_int_control input parameters that control the calculation of SE
837 !> integrals (shortrange, R3 residual, screening type)
838 !> \param se_taper ...
839 !> \par History
840 !> 04.2007 created [tlaino]
841 !> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
842 !> for computing integrals
843 !> Teodoro Laino [tlaino] - University of Zurich 04.2008 : removed the
844 !> core-core part
845 !> \author Teodoro Laino [tlaino] - Zurich University
846 !> \note
847 !> Analytical version of the MOPAC rotnuc routine
848 ! **************************************************************************************************
849  RECURSIVE SUBROUTINE corecore_el_ana(sepi, sepj, rijv, itype, enuc, denuc, &
850  se_int_control, se_taper)
851  TYPE(semi_empirical_type), POINTER :: sepi, sepj
852  REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
853  INTEGER, INTENT(IN) :: itype
854  REAL(dp), INTENT(OUT), OPTIONAL :: enuc
855  REAL(dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: denuc
856  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
857  TYPE(se_taper_type), POINTER :: se_taper
858 
859  LOGICAL :: l_denuc, l_enuc
860  REAL(dp) :: drij(3), dssss, dssss_sr, rij, ssss, &
861  ssss_sr, tmp, zz
862  TYPE(se_int_control_type) :: se_int_control_off
863 
864  rij = dot_product(rijv, rijv)
865  ! Initialization
866  l_enuc = PRESENT(enuc)
867  l_denuc = PRESENT(denuc)
868  IF ((rij > rij_threshold) .AND. (l_enuc .OR. l_denuc)) THEN
869  ! Compute Integrals in diatomic frame
870  rij = sqrt(rij)
871  CALL setup_se_int_control_type(se_int_control_off, shortrange=.false., do_ewald_r3=.false., &
872  do_ewald_gks=.false., integral_screening=se_int_control%integral_screening, &
873  max_multipole=do_multipole_none, pc_coulomb_int=.false.)
874  CALL dssss_nucint_ana(sepi, sepj, rij, ssss=ssss, dssss=dssss, itype=itype, se_taper=se_taper, &
875  se_int_control=se_int_control_off, lgrad=l_denuc)
876  ! In case let's compute the short-range part of the (ss|ss) integral
877  IF (se_int_control%shortrange .OR. se_int_control%pc_coulomb_int) THEN
878  CALL dssss_nucint_ana(sepi, sepj, rij, ssss=ssss_sr, dssss=dssss_sr, itype=itype, &
879  se_taper=se_taper, se_int_control=se_int_control, lgrad=l_denuc)
880  ELSE
881  ssss_sr = ssss
882  dssss_sr = dssss
883  END IF
884  zz = sepi%zeff*sepj%zeff
885  ! Core Core electrostatic contribution
886  IF (l_enuc) enuc = zz*ssss_sr
887  IF (l_denuc) THEN
888  drij(1) = rijv(1)/rij
889  drij(2) = rijv(2)/rij
890  drij(3) = rijv(3)/rij
891  tmp = zz*dssss_sr
892  denuc = tmp*drij
893  END IF
894  END IF
895  END SUBROUTINE corecore_el_ana
896 
897 ! **************************************************************************************************
898 !> \brief Exploits inversion symmetry to avoid divergence
899 !> \param sepi ...
900 !> \param sepj ...
901 !> \param int1el ...
902 !> \param int2el ...
903 !> \par History
904 !> 04.2007 created [tlaino]
905 !> 05.2008 New driver for integral invertion (supports d-orbitals)
906 !> \author Teodoro Laino - Zurich University
907 ! **************************************************************************************************
908  SUBROUTINE invert_integral(sepi, sepj, int1el, int2el)
909  TYPE(semi_empirical_type), POINTER :: sepi, sepj
910  REAL(dp), DIMENSION(:), INTENT(INOUT), OPTIONAL :: int1el, int2el
911 
912  INTEGER :: fdim, gind, gknd, i, imap, ind, j, jmap, &
913  jnd, k, kmap, knd, l, lmap, lnd, ndim, &
914  sdim, tdim, tind
915  REAL(kind=dp) :: ifac, jfac, kfac, lfac
916  REAL(kind=dp), DIMENSION(2025) :: tmp2el
917  REAL(kind=dp), DIMENSION(45) :: tmp1el
918 
919 ! One-electron integral
920 
921  IF (PRESENT(int1el)) THEN
922  fdim = sepi%atm_int_size
923  ndim = 0
924  DO i = 1, fdim
925  tmp1el(i) = 0.0_dp
926  END DO
927  DO i = 1, sepi%natorb
928  DO j = 1, i
929  ndim = ndim + 1
930 
931  ! Get the integral in the original frame (along z)
932  DO ind = 1, 2
933  imap = map_x_to_z(ind, i)
934  IF (imap == 0) cycle
935  ifac = fac_x_to_z(ind, i)
936  DO jnd = 1, 2
937  jmap = map_x_to_z(jnd, j)
938  IF (jmap == 0) cycle
939  jfac = fac_x_to_z(jnd, j)
940  gind = indexb(imap, jmap)
941 
942  tmp1el(ndim) = tmp1el(ndim) + ifac*jfac*int1el(gind)
943  END DO
944  END DO
945  END DO
946  END DO
947  DO i = 1, fdim
948  int1el(i) = tmp1el(i)
949  END DO
950  END IF
951 
952  ! Two electron integrals
953  IF (PRESENT(int2el)) THEN
954  sdim = sepi%atm_int_size
955  tdim = sepj%atm_int_size
956  fdim = sdim*tdim
957  ndim = 0
958  DO i = 1, fdim
959  tmp2el(i) = 0.0_dp
960  END DO
961  DO i = 1, sepi%natorb
962  DO j = 1, i
963  DO k = 1, sepj%natorb
964  DO l = 1, k
965  ndim = ndim + 1
966 
967  ! Get the integral in the original frame (along z)
968  DO ind = 1, 2
969  imap = map_x_to_z(ind, i)
970  IF (imap == 0) cycle
971  ifac = fac_x_to_z(ind, i)
972  DO jnd = 1, 2
973  jmap = map_x_to_z(jnd, j)
974  IF (jmap == 0) cycle
975  jfac = fac_x_to_z(jnd, j)
976  gind = indexb(imap, jmap)
977 
978  ! Get the integral in the original frame (along z)
979  DO knd = 1, 2
980  kmap = map_x_to_z(knd, k)
981  IF (kmap == 0) cycle
982  kfac = fac_x_to_z(knd, k)
983  DO lnd = 1, 2
984  lmap = map_x_to_z(lnd, l)
985  IF (lmap == 0) cycle
986  lfac = fac_x_to_z(lnd, l)
987  gknd = indexb(kmap, lmap)
988 
989  tind = (gind - 1)*tdim + gknd
990  tmp2el(ndim) = tmp2el(ndim) + ifac*jfac*lfac*kfac*int2el(tind)
991  END DO
992  END DO
993 
994  END DO
995  END DO
996 
997  END DO
998  END DO
999  END DO
1000  END DO
1001  DO i = 1, fdim
1002  int2el(i) = tmp2el(i)
1003  END DO
1004  END IF
1005  END SUBROUTINE invert_integral
1006 
1007 ! **************************************************************************************************
1008 !> \brief Exploits inversion symmetry to avoid divergence
1009 !> \param sepi ...
1010 !> \param sepj ...
1011 !> \param dint1el ...
1012 !> \param dint2el ...
1013 !> \par History
1014 !> 04.2007 created [tlaino]
1015 !> \author Teodoro Laino - Zurich University
1016 ! **************************************************************************************************
1017  SUBROUTINE invert_derivative(sepi, sepj, dint1el, dint2el)
1018  TYPE(semi_empirical_type), POINTER :: sepi, sepj
1019  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT), &
1020  OPTIONAL :: dint1el, dint2el
1021 
1022  INTEGER :: i, m
1023  REAL(kind=dp) :: tmp
1024 
1025 ! Integral part
1026 
1027  DO i = 1, 3
1028  IF (PRESENT(dint1el)) THEN
1029  CALL invert_integral(sepi, sepj, int1el=dint1el(i, :))
1030  END IF
1031  IF (PRESENT(dint2el)) THEN
1032  CALL invert_integral(sepi, sepj, int2el=dint2el(i, :))
1033  END IF
1034  END DO
1035 
1036  ! Derivatives part
1037  IF (PRESENT(dint1el)) THEN
1038  DO m = 1, SIZE(dint1el, 2)
1039  tmp = dint1el(3, m)
1040  dint1el(3, m) = dint1el(1, m)
1041  dint1el(1, m) = tmp
1042  END DO
1043  END IF
1044  IF (PRESENT(dint2el)) THEN
1045  DO m = 1, SIZE(dint2el, 2)
1046  tmp = dint2el(3, m)
1047  dint2el(3, m) = dint2el(1, m)
1048  dint2el(1, m) = tmp
1049  END DO
1050  END IF
1051  END SUBROUTINE invert_derivative
1052 
1053 ! **************************************************************************************************
1054 !> \brief Calculates the ssss integral and analytical derivatives (main driver)
1055 !> \param sepi parameters of atom i
1056 !> \param sepj parameters of atom j
1057 !> \param rij interatomic distance
1058 !> \param ssss ...
1059 !> \param dssss derivative of (ssss) integral
1060 !> derivatives are intended w.r.t. rij
1061 !> \param itype ...
1062 !> \param se_taper ...
1063 !> \param se_int_control input parameters that control the calculation of SE
1064 !> integrals (shortrange, R3 residual, screening type)
1065 !> \param lgrad ...
1066 !> \par History
1067 !> 03.2007 created [tlaino]
1068 !> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
1069 !> for computing integrals
1070 !> \author Teodoro Laino - Zurich University
1071 !> \note
1072 !> Analytical version - Analytical evaluation of gradients
1073 !> Teodoro Laino - Zurich University 04.2007
1074 !>
1075 ! **************************************************************************************************
1076  SUBROUTINE dssss_nucint_ana(sepi, sepj, rij, ssss, dssss, itype, se_taper, se_int_control, &
1077  lgrad)
1078  TYPE(semi_empirical_type), POINTER :: sepi, sepj
1079  REAL(dp), INTENT(IN) :: rij
1080  REAL(dp), INTENT(OUT) :: ssss, dssss
1081  INTEGER, INTENT(IN) :: itype
1082  TYPE(se_taper_type), POINTER :: se_taper
1083  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
1084  LOGICAL, INTENT(IN) :: lgrad
1085 
1086  REAL(kind=dp) :: dft, ft
1087  TYPE(se_int_screen_type) :: se_int_screen
1088 
1089 ! Compute the Tapering function
1090 
1091  ft = 1.0_dp
1092  dft = 0.0_dp
1093  IF (itype /= do_method_pchg) THEN
1094  ft = taper_eval(se_taper%taper, rij)
1095  dft = dtaper_eval(se_taper%taper, rij)
1096  END IF
1097  ! Evaluate additional taper function for dumped integrals
1098  IF (se_int_control%integral_screening == do_se_is_kdso_d) THEN
1099  se_int_screen%ft = 1.0_dp
1100  se_int_screen%dft = 0.0_dp
1101  IF (itype /= do_method_pchg) THEN
1102  se_int_screen%ft = taper_eval(se_taper%taper_add, rij)
1103  se_int_screen%dft = dtaper_eval(se_taper%taper_add, rij)
1104  END IF
1105  END IF
1106 
1107  ! Value of the integrals for sp shell
1108  CALL nucint_sp_num(sepi, sepj, rij, ssss=ssss, itype=itype, se_int_control=se_int_control, &
1109  se_int_screen=se_int_screen)
1110 
1111  IF (lgrad) THEN
1112  ! Integrals derivatives for sp shell
1113  CALL dnucint_sp_ana(sepi, sepj, rij, dssss=dssss, itype=itype, se_int_control=se_int_control, &
1114  se_int_screen=se_int_screen)
1115  END IF
1116 
1117  ! Tapering the value of the integrals
1118  IF (lgrad) THEN
1119  dssss = ft*dssss + dft*ssss
1120  END IF
1121  ssss = ft*ssss
1122 
1123  ! Debug Procedure.. Check valifity of analytical gradients of nucint
1124  IF (debug_this_module .AND. lgrad) THEN
1125  CALL check_dssss_nucint_ana(sepi, sepj, rij, dssss, itype, se_int_control, se_taper=se_taper)
1126  END IF
1127  END SUBROUTINE dssss_nucint_ana
1128 
1129 ! **************************************************************************************************
1130 !> \brief Calculates the nuclear attraction integrals and analytical integrals (main driver)
1131 !> \param sepi parameters of atom i
1132 !> \param sepj parameters of atom j
1133 !> \param rij interatomic distance
1134 !> \param core ...
1135 !> \param dcore derivative of 4 X 2 array of electron-core attraction integrals
1136 !> derivatives are intended w.r.t. rij
1137 !> The storage of the nuclear attraction integrals core(kl/ij) iS
1138 !> (SS/)=1, (SO/)=2, (OO/)=3, (PP/)=4
1139 !> where ij=1 if the orbitals centred on atom i, =2 if on atom j.
1140 !> \param itype type of semi_empirical model
1141 !> extension to the original routine to compute qm/mm
1142 !> integrals
1143 !> \param se_taper ...
1144 !> \param se_int_control input parameters that control the calculation of SE
1145 !> integrals (shortrange, R3 residual, screening type)
1146 !> \param lgrad ...
1147 !> \par History
1148 !> 03.2007 created [tlaino]
1149 !> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
1150 !> for computing integrals
1151 !> \author Teodoro Laino - Zurich University
1152 !> \note
1153 !> Analytical version - Analytical evaluation of gradients
1154 !> Teodoro Laino - Zurich University 04.2007
1155 !>
1156 ! **************************************************************************************************
1157  SUBROUTINE dcore_nucint_ana(sepi, sepj, rij, core, dcore, itype, se_taper, &
1158  se_int_control, lgrad)
1159  TYPE(semi_empirical_type), POINTER :: sepi, sepj
1160  REAL(dp), INTENT(IN) :: rij
1161  REAL(dp), DIMENSION(10, 2), INTENT(OUT) :: core, dcore
1162  INTEGER, INTENT(IN) :: itype
1163  TYPE(se_taper_type), POINTER :: se_taper
1164  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
1165  LOGICAL, INTENT(IN) :: lgrad
1166 
1167  INTEGER :: i
1168  REAL(kind=dp) :: dft, ft
1169  TYPE(se_int_screen_type) :: se_int_screen
1170 
1171 ! Compute the Tapering function
1172 
1173  ft = 1.0_dp
1174  dft = 0.0_dp
1175  IF (itype /= do_method_pchg) THEN
1176  ft = taper_eval(se_taper%taper, rij)
1177  dft = dtaper_eval(se_taper%taper, rij)
1178  END IF
1179  ! Evaluate additional taper function for dumped integrals
1180  IF (se_int_control%integral_screening == do_se_is_kdso_d) THEN
1181  se_int_screen%ft = 1.0_dp
1182  se_int_screen%dft = 0.0_dp
1183  IF (itype /= do_method_pchg) THEN
1184  se_int_screen%ft = taper_eval(se_taper%taper_add, rij)
1185  se_int_screen%dft = dtaper_eval(se_taper%taper_add, rij)
1186  END IF
1187  END IF
1188 
1189  ! Value of the integrals for sp shell
1190  CALL nucint_sp_num(sepi, sepj, rij, core=core, itype=itype, &
1191  se_int_control=se_int_control, se_int_screen=se_int_screen)
1192 
1193  IF (sepi%dorb .OR. sepj%dorb) THEN
1194  ! Compute the contribution from d-orbitals
1195  CALL nucint_d_num(sepi, sepj, rij, core, itype, &
1196  se_int_control=se_int_control, se_int_screen=se_int_screen)
1197  END IF
1198 
1199  IF (lgrad) THEN
1200  ! Integrals derivatives for sp shell
1201  CALL dnucint_sp_ana(sepi, sepj, rij, dcore=dcore, itype=itype, &
1202  se_int_control=se_int_control, se_int_screen=se_int_screen)
1203 
1204  IF (sepi%dorb .OR. sepj%dorb) THEN
1205  ! Integral derivatives involving d-orbitals
1206  CALL dnucint_d_ana(sepi, sepj, rij, dcore=dcore, itype=itype, &
1207  se_int_control=se_int_control, se_int_screen=se_int_screen)
1208  END IF
1209  END IF
1210 
1211  ! Tapering the value of the integrals
1212  IF (lgrad) THEN
1213  DO i = 1, sepi%core_size
1214  dcore(i, 1) = ft*dcore(i, 1) + dft*core(i, 1)
1215  END DO
1216  DO i = 1, sepj%core_size
1217  dcore(i, 2) = ft*dcore(i, 2) + dft*core(i, 2)
1218  END DO
1219  END IF
1220  DO i = 1, sepi%core_size
1221  core(i, 1) = ft*core(i, 1)
1222  END DO
1223  DO i = 1, sepj%core_size
1224  core(i, 2) = ft*core(i, 2)
1225  END DO
1226 
1227  ! Debug Procedure.. Check valifity of analytical gradients of nucint
1228  IF (debug_this_module .AND. lgrad) THEN
1229  CALL check_dcore_nucint_ana(sepi, sepj, rij, dcore, itype, se_int_control, se_taper=se_taper)
1230  END IF
1231  END SUBROUTINE dcore_nucint_ana
1232 
1233 ! **************************************************************************************************
1234 !> \brief Calculates the nuclear attraction integrals and derivatives for sp basis
1235 !> \param sepi parameters of atom i
1236 !> \param sepj parameters of atom j
1237 !> \param rij interatomic distance
1238 !> \param dssss derivative of (ssss) integral
1239 !> derivatives are intended w.r.t. rij
1240 !> where ij=1 if the orbitals centred on atom i, =2 if on atom j.
1241 !> \param dcore derivative of 4 X 2 array of electron-core attraction integrals
1242 !> The storage of the nuclear attraction integrals core(kl/ij) iS
1243 !> (SS/)=1, (SP/)=2, (PP/)=3, (P+P+/)=4, (SD/)=5,
1244 !> (DP/)=6, (DD/)=7, (D+P+)=8, (D+D+/)=9, (D#D#)=10
1245 !> \param itype type of semi_empirical model
1246 !> extension to the original routine to compute qm/mm
1247 !> integrals
1248 !> \param se_int_control input parameters that control the calculation of SE
1249 !> integrals (shortrange, R3 residual, screening type)
1250 !> \param se_int_screen ...
1251 !> \par History
1252 !> 04.2007 created [tlaino]
1253 !> Teodoro Laino (03.2008) [tlaino] - University of Zurich : new driver
1254 !> for computing integrals
1255 !> 05.2008 Teodoro Laino [tlaino] - University of Zurich: major rewriting
1256 !> \author Teodoro Laino - Zurich University
1257 !> \note
1258 !> Analytical version - Analytical evaluation of gradients
1259 !> Teodoro Laino - Zurich University 04.2007
1260 !> routine adapted from mopac7 (repp)
1261 !> vector version written by Ernest R. Davidson, Indiana University
1262 ! **************************************************************************************************
1263  SUBROUTINE dnucint_sp_ana(sepi, sepj, rij, dssss, dcore, itype, se_int_control, &
1264  se_int_screen)
1265  TYPE(semi_empirical_type), POINTER :: sepi, sepj
1266  REAL(dp), INTENT(IN) :: rij
1267  REAL(dp), INTENT(INOUT), OPTIONAL :: dssss
1268  REAL(dp), DIMENSION(10, 2), INTENT(INOUT), &
1269  OPTIONAL :: dcore
1270  INTEGER, INTENT(IN) :: itype
1271  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
1272  TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen
1273 
1274  INTEGER :: ij, kl
1275  LOGICAL :: l_core, l_ssss, si, sj
1276 
1277  l_core = PRESENT(dcore)
1278  l_ssss = PRESENT(dssss)
1279  IF (.NOT. (l_core .OR. l_ssss)) RETURN
1280 
1281  si = (sepi%natorb > 1)
1282  sj = (sepj%natorb > 1)
1283 
1284  ij = indexa(1, 1)
1285  IF (l_ssss) THEN
1286  ! Store the value for the derivative of <S S | S S > (Used for computing the core-core interactions)
1287  dssss = d_ijkl_sp(sepi, sepj, ij, ij, 0, 0, 0, 0, -1, rij, se_int_control, se_int_screen, itype)
1288  END IF
1289 
1290  IF (l_core) THEN
1291  ! <S S | S S >
1292  kl = indexa(1, 1)
1293  dcore(1, 1) = d_ijkl_sp(sepi, sepj, kl, ij, 0, 0, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
1294  IF (si) THEN
1295  ! <S P | S S >
1296  kl = indexa(2, 1)
1297  dcore(2, 1) = d_ijkl_sp(sepi, sepj, kl, ij, 0, 1, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
1298  ! <P P | S S >
1299  kl = indexa(2, 2)
1300  dcore(3, 1) = d_ijkl_sp(sepi, sepj, kl, ij, 1, 1, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
1301  ! <P+ P+ | S S >
1302  kl = indexa(3, 3)
1303  dcore(4, 1) = d_ijkl_sp(sepi, sepj, kl, ij, 1, 1, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
1304  END IF
1305 
1306  ! <S S | S S >
1307  kl = indexa(1, 1)
1308  dcore(1, 2) = d_ijkl_sp(sepi, sepj, ij, kl, 0, 0, 0, 0, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
1309  IF (sj) THEN
1310  ! <S S | S P >
1311  kl = indexa(2, 1)
1312  dcore(2, 2) = d_ijkl_sp(sepi, sepj, ij, kl, 0, 0, 0, 1, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
1313  ! <S S | P P >
1314  kl = indexa(2, 2)
1315  dcore(3, 2) = d_ijkl_sp(sepi, sepj, ij, kl, 0, 0, 1, 1, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
1316  ! <S S | P+ P+ >
1317  kl = indexa(3, 3)
1318  dcore(4, 2) = d_ijkl_sp(sepi, sepj, ij, kl, 0, 0, 1, 1, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
1319  END IF
1320  END IF
1321  END SUBROUTINE dnucint_sp_ana
1322 
1323 ! **************************************************************************************************
1324 !> \brief Calculates the analytical derivative of the nuclear attraction
1325 !> integrals involving d orbitals
1326 !> \param sepi parameters of atom i
1327 !> \param sepj parameters of atom j
1328 !> \param rij interatomic distance
1329 !> \param dcore 4 X 2 array of electron-core attraction integrals
1330 !> The storage of the nuclear attraction integrals core(kl/ij) iS
1331 !> (SS/)=1, (SP/)=2, (PP/)=3, (P+P+/)=4, (SD/)=5,
1332 !> (DP/)=6, (DD/)=7, (D+P+)=8, (D+D+/)=9, (D#D#)=10
1333 !>
1334 !> where ij=1 if the orbitals centred on atom i, =2 if on atom j.
1335 !> \param itype type of semi_empirical model
1336 !> extension to the original routine to compute qm/mm
1337 !> integrals
1338 !> \param se_int_control input parameters that control the calculation of SE
1339 !> integrals (shortrange, R3 residual, screening type)
1340 !> \param se_int_screen ...
1341 !> \author
1342 !> Teodoro Laino (05.2008) [tlaino] - University of Zurich: created
1343 ! **************************************************************************************************
1344  SUBROUTINE dnucint_d_ana(sepi, sepj, rij, dcore, itype, se_int_control, &
1345  se_int_screen)
1346  TYPE(semi_empirical_type), POINTER :: sepi, sepj
1347  REAL(dp), INTENT(IN) :: rij
1348  REAL(dp), DIMENSION(10, 2), INTENT(INOUT) :: dcore
1349  INTEGER, INTENT(IN) :: itype
1350  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
1351  TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen
1352 
1353  INTEGER :: ij, kl
1354 
1355 ! Check if d-orbitals are present
1356 
1357  IF (sepi%dorb .OR. sepj%dorb) THEN
1358  ij = indexa(1, 1)
1359  IF (sepj%dorb) THEN
1360  ! <S S | D S>
1361  kl = indexa(5, 1)
1362  dcore(5, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 0, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
1363  ! <S S | D P >
1364  kl = indexa(5, 2)
1365  dcore(6, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 1, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
1366  ! <S S | D D >
1367  kl = indexa(5, 5)
1368  dcore(7, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
1369  ! <S S | D+P+>
1370  kl = indexa(6, 3)
1371  dcore(8, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 1, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
1372  ! <S S | D+D+>
1373  kl = indexa(6, 6)
1374  dcore(9, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
1375  ! <S S | D#D#>
1376  kl = indexa(8, 8)
1377  dcore(10, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
1378  END IF
1379  IF (sepi%dorb) THEN
1380  ! <D S | S S>
1381  kl = indexa(5, 1)
1382  dcore(5, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 0, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
1383  ! <D P | S S >
1384  kl = indexa(5, 2)
1385  dcore(6, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 1, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
1386  ! <D D | S S >
1387  kl = indexa(5, 5)
1388  dcore(7, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
1389  ! <D+P+| S S >
1390  kl = indexa(6, 3)
1391  dcore(8, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 1, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
1392  ! <D+D+| S S >
1393  kl = indexa(6, 6)
1394  dcore(9, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
1395  ! <D#D#| S S >
1396  kl = indexa(8, 8)
1397  dcore(10, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
1398  END IF
1399  END IF
1400  END SUBROUTINE dnucint_d_ana
1401 
1402 ! **************************************************************************************************
1403 !> \brief calculates the derivative of the two-particle interactions
1404 !> \param sepi Atomic parameters of first atom
1405 !> \param sepj Atomic parameters of second atom
1406 !> \param rijv Coordinate vector i -> j
1407 !> \param w Array of two-electron repulsion integrals.
1408 !> \param dw ...
1409 !> \param se_int_control ...
1410 !> \param se_taper ...
1411 !> \par History
1412 !> 04.2007 created [tlaino]
1413 !> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
1414 !> for computing integrals
1415 !> \author Teodoro Laino - Zurich University
1416 !> \note
1417 !> Analytical version - Analytical evaluation of gradients
1418 !> Teodoro Laino - Zurich University 04.2007
1419 !> routine adapted from mopac7 (repp)
1420 !> vector version written by Ernest R. Davidson, Indiana University
1421 ! **************************************************************************************************
1422  RECURSIVE SUBROUTINE rotint_ana(sepi, sepj, rijv, w, dw, se_int_control, se_taper)
1423  TYPE(semi_empirical_type), POINTER :: sepi, sepj
1424  REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
1425  REAL(dp), DIMENSION(2025), INTENT(OUT), OPTIONAL :: w
1426  REAL(dp), DIMENSION(3, 2025), INTENT(OUT), &
1427  OPTIONAL :: dw
1428  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
1429  TYPE(se_taper_type), POINTER :: se_taper
1430 
1431  INTEGER :: i, i1, ii, ij, ij1, iminus, istep, &
1432  iw_loc, j, j1, jj, k, kk, kl, l, &
1433  limij, limkl, mm
1434  LOGICAL :: invert, l_w, lgrad
1435  LOGICAL, DIMENSION(45, 45) :: logv, logv_d
1436  REAL(dp) :: rij, xtmp
1437  REAL(dp), DIMENSION(3) :: drij
1438  REAL(kind=dp) :: cc, cc_d(3), wrepp, wrepp_d(3)
1439  REAL(kind=dp), DIMENSION(2025) :: ww
1440  REAL(kind=dp), DIMENSION(3, 2025) :: ww_d
1441  REAL(kind=dp), DIMENSION(3, 45, 45) :: v_d
1442  REAL(kind=dp), DIMENSION(45, 45) :: v
1443  REAL(kind=dp), DIMENSION(491) :: rep, rep_d
1444  TYPE(rotmat_type), POINTER :: ij_matrix
1445 
1446  NULLIFY (ij_matrix)
1447  l_w = PRESENT(w)
1448  lgrad = PRESENT(dw)
1449  IF (.NOT. (l_w .OR. lgrad)) RETURN
1450 
1451  rij = dot_product(rijv, rijv)
1452  IF (rij > rij_threshold) THEN
1453  ! The repulsion integrals over molecular frame (w) are stored in the
1454  ! order in which they will later be used. ie. (i,j/k,l) where
1455  ! j.le.i and l.le.k and l varies most rapidly and i least
1456  ! rapidly. (anti-normal computer storage)
1457  rij = sqrt(rij)
1458 
1459  ! Create the rotation matrix
1460  CALL rotmat_create(ij_matrix)
1461  CALL rotmat(sepi, sepj, rijv, rij, ij_matrix, do_derivatives=lgrad, do_invert=invert)
1462 
1463  ! Compute integrals in diatomic frame as well their derivatives (if requested)
1464  CALL dterep_ana(sepi, sepj, rij, rep, rep_d, se_taper, se_int_control, lgrad=lgrad)
1465 
1466  IF (lgrad) THEN
1467  drij(1) = rijv(1)/rij
1468  drij(2) = rijv(2)/rij
1469  drij(3) = rijv(3)/rij
1470  ! Possibly Invert Frame
1471  IF (invert) THEN
1472  xtmp = drij(3)
1473  drij(3) = drij(1)
1474  drij(1) = xtmp
1475  END IF
1476  END IF
1477 
1478  ii = sepi%natorb
1479  kk = sepj%natorb
1480  ! First step in rotation of integrals
1481  CALL rot_2el_2c_first(sepi, sepj, rijv, se_int_control, se_taper, invert, ii, kk, rep, logv, ij_matrix, &
1482  v, lgrad, rep_d, v_d, logv_d, drij)
1483 
1484  ! Integrals if requested
1485  IF (l_w) THEN
1486  ! Rotate Integrals
1487  IF (ii*kk > 0) THEN
1488  limij = sepi%atm_int_size
1489  limkl = sepj%atm_int_size
1490  istep = limkl*limij
1491  DO i1 = 1, istep
1492  ww(i1) = 0.0_dp
1493  END DO
1494  ! Second step in rotation of integrals
1495  DO i1 = 1, ii
1496  DO j1 = 1, i1
1497  ij = indexa(i1, j1)
1498  jj = indexb(i1, j1)
1499  mm = int2c_type(jj)
1500  DO k = 1, kk
1501  DO l = 1, k
1502  kl = indexb(k, l)
1503  IF (logv(ij, kl)) THEN
1504  wrepp = v(ij, kl)
1505  SELECT CASE (mm)
1506  CASE (1)
1507  ! (SS/)
1508  i = 1
1509  j = 1
1510  iw_loc = (indexb(i, j) - 1)*limkl + kl
1511  ww(iw_loc) = wrepp
1512  CASE (2)
1513  ! (SP/)
1514  j = 1
1515  DO i = 1, 3
1516  iw_loc = (indexb(i + 1, j) - 1)*limkl + kl
1517  ww(iw_loc) = ww(iw_loc) + ij_matrix%sp(i1 - 1, i)*wrepp
1518  END DO
1519  CASE (3)
1520  ! (PP/)
1521  DO i = 1, 3
1522  cc = ij_matrix%pp(i, i1 - 1, j1 - 1)
1523  iw_loc = (indexb(i + 1, i + 1) - 1)*limkl + kl
1524  ww(iw_loc) = ww(iw_loc) + cc*wrepp
1525  iminus = i - 1
1526  IF (iminus /= 0) THEN
1527  DO j = 1, iminus
1528  cc = ij_matrix%pp(1 + i + j, i1 - 1, j1 - 1)
1529  iw_loc = (indexb(i + 1, j + 1) - 1)*limkl + kl
1530  ww(iw_loc) = ww(iw_loc) + cc*wrepp
1531  END DO
1532  END IF
1533  END DO
1534  CASE (4)
1535  ! (SD/)
1536  j = 1
1537  DO i = 1, 5
1538  iw_loc = (indexb(i + 4, j) - 1)*limkl + kl
1539  ww(iw_loc) = ww(iw_loc) + ij_matrix%sd(i1 - 4, i)*wrepp
1540  END DO
1541  CASE (5)
1542  ! (DP/)
1543  DO i = 1, 5
1544  DO j = 1, 3
1545  iw_loc = (indexb(i + 4, j + 1) - 1)*limkl + kl
1546  ij1 = 3*(i - 1) + j
1547  ww(iw_loc) = ww(iw_loc) + ij_matrix%pd(ij1, i1 - 4, j1 - 1)*wrepp
1548  END DO
1549  END DO
1550  CASE (6)
1551  ! (DD/)
1552  DO i = 1, 5
1553  cc = ij_matrix%dd(i, i1 - 4, j1 - 4)
1554  iw_loc = (indexb(i + 4, i + 4) - 1)*limkl + kl
1555  ww(iw_loc) = ww(iw_loc) + cc*wrepp
1556  iminus = i - 1
1557  IF (iminus /= 0) THEN
1558  DO j = 1, iminus
1559  ij1 = inddd(i, j)
1560  cc = ij_matrix%dd(ij1, i1 - 4, j1 - 4)
1561  iw_loc = (indexb(i + 4, j + 4) - 1)*limkl + kl
1562  ww(iw_loc) = ww(iw_loc) + cc*wrepp
1563  END DO
1564  END IF
1565  END DO
1566  END SELECT
1567  END IF
1568  END DO
1569  END DO
1570  END DO
1571  END DO
1572  ! Store two electron integrals in the triangular format
1573  CALL store_2el_2c_diag(limij, limkl, ww(1:istep), w)
1574  IF (invert) CALL invert_integral(sepi, sepj, int2el=w)
1575  END IF
1576 
1577  IF (debug_this_module) THEN
1578  ! Check value of integrals
1579  CALL check_rotint_ana(sepi, sepj, rijv, w, se_int_control=se_int_control, se_taper=se_taper)
1580  END IF
1581  END IF
1582 
1583  ! Gradients if requested
1584  IF (lgrad) THEN
1585  ! Rotate Integrals derivatives
1586  IF (ii*kk > 0) THEN
1587  limij = sepi%atm_int_size
1588  limkl = sepj%atm_int_size
1589  istep = limkl*limij
1590  DO i1 = 1, istep
1591  ww_d(1, i1) = 0.0_dp
1592  ww_d(2, i1) = 0.0_dp
1593  ww_d(3, i1) = 0.0_dp
1594  END DO
1595 
1596  ! Second step in rotation of integrals
1597  DO i1 = 1, ii
1598  DO j1 = 1, i1
1599  ij = indexa(i1, j1)
1600  jj = indexb(i1, j1)
1601  mm = int2c_type(jj)
1602  DO k = 1, kk
1603  DO l = 1, k
1604  kl = indexb(k, l)
1605  IF (logv_d(ij, kl)) THEN
1606  wrepp_d(1) = v_d(1, ij, kl)
1607  wrepp_d(2) = v_d(2, ij, kl)
1608  wrepp_d(3) = v_d(3, ij, kl)
1609  wrepp = v(ij, kl)
1610  SELECT CASE (mm)
1611  CASE (1)
1612  ! (SS/)
1613  i = 1
1614  j = 1
1615  iw_loc = (indexb(i, j) - 1)*limkl + kl
1616  ww_d(1, iw_loc) = wrepp_d(1)
1617  ww_d(2, iw_loc) = wrepp_d(2)
1618  ww_d(3, iw_loc) = wrepp_d(3)
1619  CASE (2)
1620  ! (SP/)
1621  j = 1
1622  DO i = 1, 3
1623  iw_loc = (indexb(i + 1, j) - 1)*limkl + kl
1624  ww_d(1, iw_loc) = ww_d(1, iw_loc) + ij_matrix%sp_d(1, i1 - 1, i)*wrepp + &
1625  ij_matrix%sp(i1 - 1, i)*wrepp_d(1)
1626 
1627  ww_d(2, iw_loc) = ww_d(2, iw_loc) + ij_matrix%sp_d(2, i1 - 1, i)*wrepp + &
1628  ij_matrix%sp(i1 - 1, i)*wrepp_d(2)
1629 
1630  ww_d(3, iw_loc) = ww_d(3, iw_loc) + ij_matrix%sp_d(3, i1 - 1, i)*wrepp + &
1631  ij_matrix%sp(i1 - 1, i)*wrepp_d(3)
1632  END DO
1633  CASE (3)
1634  ! (PP/)
1635  DO i = 1, 3
1636  cc = ij_matrix%pp(i, i1 - 1, j1 - 1)
1637  cc_d(1) = ij_matrix%pp_d(1, i, i1 - 1, j1 - 1)
1638  cc_d(2) = ij_matrix%pp_d(2, i, i1 - 1, j1 - 1)
1639  cc_d(3) = ij_matrix%pp_d(3, i, i1 - 1, j1 - 1)
1640  iw_loc = (indexb(i + 1, i + 1) - 1)*limkl + kl
1641  ww_d(1, iw_loc) = ww_d(1, iw_loc) + cc_d(1)*wrepp + cc*wrepp_d(1)
1642  ww_d(2, iw_loc) = ww_d(2, iw_loc) + cc_d(2)*wrepp + cc*wrepp_d(2)
1643  ww_d(3, iw_loc) = ww_d(3, iw_loc) + cc_d(3)*wrepp + cc*wrepp_d(3)
1644  iminus = i - 1
1645  IF (iminus /= 0) THEN
1646  DO j = 1, iminus
1647  cc = ij_matrix%pp(1 + i + j, i1 - 1, j1 - 1)
1648  cc_d(1) = ij_matrix%pp_d(1, 1 + i + j, i1 - 1, j1 - 1)
1649  cc_d(2) = ij_matrix%pp_d(2, 1 + i + j, i1 - 1, j1 - 1)
1650  cc_d(3) = ij_matrix%pp_d(3, 1 + i + j, i1 - 1, j1 - 1)
1651  iw_loc = (indexb(i + 1, j + 1) - 1)*limkl + kl
1652  ww_d(1, iw_loc) = ww_d(1, iw_loc) + cc_d(1)*wrepp + cc*wrepp_d(1)
1653  ww_d(2, iw_loc) = ww_d(2, iw_loc) + cc_d(2)*wrepp + cc*wrepp_d(2)
1654  ww_d(3, iw_loc) = ww_d(3, iw_loc) + cc_d(3)*wrepp + cc*wrepp_d(3)
1655  END DO
1656  END IF
1657  END DO
1658  CASE (4)
1659  ! (SD/)
1660  j = 1
1661  DO i = 1, 5
1662  iw_loc = (indexb(i + 4, j) - 1)*limkl + kl
1663  ww_d(1, iw_loc) = ww_d(1, iw_loc) + ij_matrix%sd_d(1, i1 - 4, i)*wrepp + &
1664  ij_matrix%sd(i1 - 4, i)*wrepp_d(1)
1665 
1666  ww_d(2, iw_loc) = ww_d(2, iw_loc) + ij_matrix%sd_d(2, i1 - 4, i)*wrepp + &
1667  ij_matrix%sd(i1 - 4, i)*wrepp_d(2)
1668 
1669  ww_d(3, iw_loc) = ww_d(3, iw_loc) + ij_matrix%sd_d(3, i1 - 4, i)*wrepp + &
1670  ij_matrix%sd(i1 - 4, i)*wrepp_d(3)
1671  END DO
1672  CASE (5)
1673  ! (DP/)
1674  DO i = 1, 5
1675  DO j = 1, 3
1676  iw_loc = (indexb(i + 4, j + 1) - 1)*limkl + kl
1677  ij1 = 3*(i - 1) + j
1678  ww_d(1, iw_loc) = ww_d(1, iw_loc) + ij_matrix%pd_d(1, ij1, i1 - 4, j1 - 1)*wrepp + &
1679  ij_matrix%pd(ij1, i1 - 4, j1 - 1)*wrepp_d(1)
1680 
1681  ww_d(2, iw_loc) = ww_d(2, iw_loc) + ij_matrix%pd_d(2, ij1, i1 - 4, j1 - 1)*wrepp + &
1682  ij_matrix%pd(ij1, i1 - 4, j1 - 1)*wrepp_d(2)
1683 
1684  ww_d(3, iw_loc) = ww_d(3, iw_loc) + ij_matrix%pd_d(3, ij1, i1 - 4, j1 - 1)*wrepp + &
1685  ij_matrix%pd(ij1, i1 - 4, j1 - 1)*wrepp_d(3)
1686  END DO
1687  END DO
1688  CASE (6)
1689  ! (DD/)
1690  DO i = 1, 5
1691  cc = ij_matrix%dd(i, i1 - 4, j1 - 4)
1692  cc_d = ij_matrix%dd_d(:, i, i1 - 4, j1 - 4)
1693  iw_loc = (indexb(i + 4, i + 4) - 1)*limkl + kl
1694  ww_d(1, iw_loc) = ww_d(1, iw_loc) + cc_d(1)*wrepp + cc*wrepp_d(1)
1695  ww_d(2, iw_loc) = ww_d(2, iw_loc) + cc_d(2)*wrepp + cc*wrepp_d(2)
1696  ww_d(3, iw_loc) = ww_d(3, iw_loc) + cc_d(3)*wrepp + cc*wrepp_d(3)
1697  iminus = i - 1
1698  IF (iminus /= 0) THEN
1699  DO j = 1, iminus
1700  ij1 = inddd(i, j)
1701  cc = ij_matrix%dd(ij1, i1 - 4, j1 - 4)
1702  cc_d(1) = ij_matrix%dd_d(1, ij1, i1 - 4, j1 - 4)
1703  cc_d(2) = ij_matrix%dd_d(2, ij1, i1 - 4, j1 - 4)
1704  cc_d(3) = ij_matrix%dd_d(3, ij1, i1 - 4, j1 - 4)
1705  iw_loc = (indexb(i + 4, j + 4) - 1)*limkl + kl
1706  ww_d(1, iw_loc) = ww_d(1, iw_loc) + cc_d(1)*wrepp + cc*wrepp_d(1)
1707  ww_d(2, iw_loc) = ww_d(2, iw_loc) + cc_d(2)*wrepp + cc*wrepp_d(2)
1708  ww_d(3, iw_loc) = ww_d(3, iw_loc) + cc_d(3)*wrepp + cc*wrepp_d(3)
1709  END DO
1710  END IF
1711  END DO
1712  END SELECT
1713  END IF
1714  END DO
1715  END DO
1716  END DO
1717  END DO
1718  ! Store two electron integrals in the triangular format
1719  CALL store_2el_2c_diag(limij, limkl, ww_dx=ww_d(1, 1:istep), ww_dy=ww_d(2, 1:istep), ww_dz=ww_d(3, 1:istep), &
1720  dw=dw)
1721  IF (invert) CALL invert_derivative(sepi, sepj, dint2el=dw)
1722  END IF
1723 
1724  IF (debug_this_module) THEN
1725  ! Check derivatives
1726  CALL check_rotint_ana(sepi, sepj, rijv, dw=dw, se_int_control=se_int_control, se_taper=se_taper)
1727  END IF
1728  END IF
1729  CALL rotmat_release(ij_matrix)
1730  END IF
1731  END SUBROUTINE rotint_ana
1732 
1733 ! **************************************************************************************************
1734 !> \brief Calculates the derivative and the value of two-electron repulsion
1735 !> integrals and the nuclear attraction integrals w.r.t. |r|
1736 !> \param sepi parameters of atom i
1737 !> \param sepj parameters of atom j
1738 !> \param rij interatomic distance
1739 !> \param rep rray of two-electron repulsion integrals
1740 !> \param rep_d array of two-electron repulsion integrals derivatives
1741 !> \param se_taper ...
1742 !> \param se_int_control input parameters that control the calculation of SE
1743 !> integrals (shortrange, R3 residual, screening type)
1744 !> \param lgrad ...
1745 !> \par History
1746 !> 03.2008 created [tlaino]
1747 !> \author Teodoro Laino [tlaino] - Zurich University
1748 ! **************************************************************************************************
1749  RECURSIVE SUBROUTINE dterep_ana(sepi, sepj, rij, rep, rep_d, se_taper, &
1750  se_int_control, lgrad)
1751  TYPE(semi_empirical_type), POINTER :: sepi, sepj
1752  REAL(kind=dp), INTENT(IN) :: rij
1753  REAL(kind=dp), DIMENSION(491), INTENT(OUT) :: rep, rep_d
1754  TYPE(se_taper_type), POINTER :: se_taper
1755  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
1756  LOGICAL, INTENT(IN) :: lgrad
1757 
1758  INTEGER :: i, ij, j, k, kl, l, lasti, lastj, li, &
1759  lj, lk, ll, numb
1760  REAL(kind=dp) :: dft, ft, ft1
1761  TYPE(se_int_screen_type) :: se_int_screen
1762 
1763 ! Compute the tapering function and its derivatives
1764 
1765  ft = taper_eval(se_taper%taper, rij)
1766  dft = 0.0_dp
1767  ft1 = ft
1768  IF (lgrad) THEN
1769  ft1 = 1.0_dp
1770  dft = dtaper_eval(se_taper%taper, rij)
1771  END IF
1772  ! Evaluate additional taper function for dumped integrals
1773  IF (se_int_control%integral_screening == do_se_is_kdso_d) THEN
1774  se_int_screen%ft = taper_eval(se_taper%taper_add, rij)
1775  IF (lgrad) &
1776  se_int_screen%dft = dtaper_eval(se_taper%taper_add, rij)
1777  END IF
1778 
1779  ! Integral Values for sp shells only
1780  CALL terep_sp_num(sepi, sepj, rij, rep, se_int_control=se_int_control, &
1781  se_int_screen=se_int_screen, ft=ft1)
1782 
1783  IF (sepi%dorb .OR. sepj%dorb) THEN
1784  ! Compute the contribution from d-orbitals
1785  CALL terep_d_num(sepi, sepj, rij, rep, se_int_control=se_int_control, &
1786  se_int_screen=se_int_screen, ft=ft1)
1787  END IF
1788 
1789  IF (lgrad) THEN
1790  ! Integral Derivatives
1791  CALL dterep_sp_ana(sepi, sepj, rij, rep_d, rep, se_int_control, &
1792  se_int_screen, ft, dft)
1793 
1794  IF (sepi%dorb .OR. sepj%dorb) THEN
1795  ! Compute the derivatives from d-orbitals
1796  CALL dterep_d_ana(sepi, sepj, rij, rep_d, rep, se_int_control, &
1797  se_int_screen, ft, dft)
1798  END IF
1799 
1800  ! Tapering Integral values
1801  lasti = sepi%natorb
1802  lastj = sepj%natorb
1803  DO i = 1, lasti
1804  li = l_index(i)
1805  DO j = 1, i
1806  lj = l_index(j)
1807  ij = indexa(i, j)
1808  DO k = 1, lastj
1809  lk = l_index(k)
1810  DO l = 1, k
1811  ll = l_index(l)
1812  kl = indexa(k, l)
1813  numb = ijkl_ind(ij, kl)
1814  IF (numb > 0) rep(numb) = rep(numb)*ft
1815  END DO
1816  END DO
1817  END DO
1818  END DO
1819  END IF
1820 
1821  ! Possibly debug 2el 2cent integrals and derivatives
1822  IF (debug_this_module) THEN
1823  CALL check_dterep_ana(sepi, sepj, rij, rep, rep_d, se_int_control, se_taper=se_taper, &
1824  lgrad=lgrad)
1825  END IF
1826  END SUBROUTINE dterep_ana
1827 
1828 ! **************************************************************************************************
1829 !> \brief Calculates the derivative and the value of two-electron repulsion
1830 !> integrals and the nuclear attraction integrals w.r.t. |r| - sp shells only
1831 !> \param sepi parameters of atom i
1832 !> \param sepj parameters of atom j
1833 !> \param rij interatomic distance
1834 !> \param drep array of derivatives of two-electron repulsion integrals
1835 !> \param rep array of two-electron repulsion integrals
1836 !> \param se_int_control input parameters that control the calculation of SE
1837 !> integrals (shortrange, R3 residual, screening type)
1838 !> \param se_int_screen ...
1839 !> \param ft ...
1840 !> \param dft ...
1841 !> \par History
1842 !> 04.2007 created [tlaino]
1843 !> Teodoro Laino (03.2008) [tlaino] - University of Zurich : new driver
1844 !> for computing integrals
1845 !> 05.2008 Teodoro Laino [tlaino] - University of Zurich: major rewriting
1846 !> \author Teodoro Laino - Zurich University
1847 !> \note
1848 !> Analytical version - Analytical evaluation of gradients
1849 !> Teodoro Laino - Zurich University 04.2007
1850 !> routine adapted from mopac7 (repp)
1851 !> vector version written by Ernest R. Davidson, Indiana University
1852 ! **************************************************************************************************
1853  SUBROUTINE dterep_sp_ana(sepi, sepj, rij, drep, rep, se_int_control, &
1854  se_int_screen, ft, dft)
1855  TYPE(semi_empirical_type), POINTER :: sepi, sepj
1856  REAL(dp), INTENT(IN) :: rij
1857  REAL(dp), DIMENSION(491), INTENT(OUT) :: drep
1858  REAL(dp), DIMENSION(491), INTENT(IN) :: rep
1859  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
1860  TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen
1861  REAL(dp), INTENT(IN) :: ft, dft
1862 
1863  INTEGER :: i, ij, j, k, kl, l, lasti, lastj, li, &
1864  lj, lk, ll, nold, numb
1865  REAL(kind=dp) :: tmp
1866 
1867  lasti = sepi%natorb
1868  lastj = sepj%natorb
1869  DO i = 1, min(lasti, 4)
1870  li = l_index(i)
1871  DO j = 1, i
1872  lj = l_index(j)
1873  ij = indexa(i, j)
1874  DO k = 1, min(lastj, 4)
1875  lk = l_index(k)
1876  DO l = 1, k
1877  ll = l_index(l)
1878  kl = indexa(k, l)
1879 
1880  numb = ijkl_ind(ij, kl)
1881  IF (numb > 0) THEN
1882  nold = ijkl_sym(numb)
1883  IF (nold > 0) THEN
1884  drep(numb) = drep(nold)
1885  ELSE IF (nold < 0) THEN
1886  drep(numb) = -drep(-nold)
1887  ELSE IF (nold == 0) THEN
1888  tmp = d_ijkl_sp(sepi, sepj, ij, kl, li, lj, lk, ll, 0, rij, &
1889  se_int_control, se_int_screen, do_method_undef)
1890  drep(numb) = dft*rep(numb) + ft*tmp
1891  END IF
1892  END IF
1893  END DO
1894  END DO
1895  END DO
1896  END DO
1897  END SUBROUTINE dterep_sp_ana
1898 
1899 ! **************************************************************************************************
1900 !> \brief Calculates the derivatives of the two-electron repulsion integrals - d shell only
1901 !> \param sepi parameters of atom i
1902 !> \param sepj parameters of atom j
1903 !> \param rij interatomic distance
1904 !> \param drep ...
1905 !> \param rep array of two-electron repulsion integrals
1906 !> \param se_int_control input parameters that control the calculation of
1907 !> integrals (shortrange, R3 residual, screening type)
1908 !> \param se_int_screen ...
1909 !> \param ft ...
1910 !> \param dft ...
1911 !> \par History
1912 !> Teodoro Laino (05.2008) [tlaino] - University of Zurich : new driver
1913 !> for computing integral derivatives for d-orbitals
1914 ! **************************************************************************************************
1915  SUBROUTINE dterep_d_ana(sepi, sepj, rij, drep, rep, se_int_control, &
1916  se_int_screen, ft, dft)
1917  TYPE(semi_empirical_type), POINTER :: sepi, sepj
1918  REAL(dp), INTENT(IN) :: rij
1919  REAL(dp), DIMENSION(491), INTENT(INOUT) :: drep
1920  REAL(dp), DIMENSION(491), INTENT(IN) :: rep
1921  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
1922  TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen
1923  REAL(dp), INTENT(IN) :: ft, dft
1924 
1925  INTEGER :: i, ij, j, k, kl, l, lasti, lastj, li, &
1926  lj, lk, ll, nold, numb
1927  REAL(kind=dp) :: tmp
1928 
1929  lasti = sepi%natorb
1930  lastj = sepj%natorb
1931  DO i = 1, lasti
1932  li = l_index(i)
1933  DO j = 1, i
1934  lj = l_index(j)
1935  ij = indexa(i, j)
1936  DO k = 1, lastj
1937  lk = l_index(k)
1938  DO l = 1, k
1939  ll = l_index(l)
1940  kl = indexa(k, l)
1941 
1942  numb = ijkl_ind(ij, kl)
1943  ! From 1 to 34 we store integrals involving sp shells
1944  IF (numb > 34) THEN
1945  nold = ijkl_sym(numb)
1946  IF (nold > 34) THEN
1947  drep(numb) = drep(nold)
1948  ELSE IF (nold < -34) THEN
1949  drep(numb) = -drep(-nold)
1950  ELSE IF (nold == 0) THEN
1951  tmp = d_ijkl_d(sepi, sepj, ij, kl, li, lj, lk, ll, 0, rij, &
1952  se_int_control, se_int_screen, do_method_undef)
1953  drep(numb) = dft*rep(numb) + ft*tmp
1954  END IF
1955  END IF
1956  END DO
1957  END DO
1958  END DO
1959  END DO
1960  END SUBROUTINE dterep_d_ana
1961 
1962 END MODULE semi_empirical_int_ana
semi_empirical_int_ana::check_dcore_nucint_ana
Check Numerical Vs Analytical NUCINT core.
Definition: semi_empirical_int_ana.F:112
semi_empirical_int_ana::check_dcorecore_ana
Check Numerical Vs Analytical CORECORE.
Definition: semi_empirical_int_ana.F:163
semi_empirical_int_ana::check_drotnuc_ana
Check Numerical Vs Analytical ROTNUC.
Definition: semi_empirical_int_ana.F:137
semi_empirical_int_ana::check_dssss_nucint_ana
Check Numerical Vs Analytical NUCINT ssss.
Definition: semi_empirical_int_ana.F:88
semi_empirical_int_ana::check_dterep_ana
Check Numerical Vs Analytical check_dterep_ana.
Definition: semi_empirical_int_ana.F:216
semi_empirical_int_ana::check_rotint_ana
Check Numerical Vs Analytical check_rotint_ana.
Definition: semi_empirical_int_ana.F:240
semi_empirical_int_ana::rot_2el_2c_first_debug
Check Numerical Vs Analytical rot_2el_2c_first.
Definition: semi_empirical_int_ana.F:190
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_method_pdg
integer, parameter, public do_method_pdg
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
input_constants::do_method_undef
integer, parameter, public do_method_undef
Definition: input_constants.F:250
input_constants::do_method_pm3
integer, parameter, public do_method_pm3
Definition: input_constants.F:250
input_constants::do_method_am1
integer, parameter, public do_method_am1
Definition: input_constants.F:250
input_constants::do_method_pm6fm
integer, parameter, public do_method_pm6fm
Definition: input_constants.F:250
input_constants::do_method_pm6
integer, parameter, public do_method_pm6
Definition: input_constants.F:250
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
multipole_types
Multipole structure: for multipole (fixed and induced) in FF based MD.
Definition: multipole_types.F:12
multipole_types::do_multipole_none
integer, parameter, public do_multipole_none
Definition: multipole_types.F:31
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition: physcon.F:183
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition: physcon.F:144
semi_empirical_int_ana
Analytical derivatives of Integrals for semi-empirical methods.
Definition: semi_empirical_int_ana.F:18
semi_empirical_int_ana::rotint_ana
recursive subroutine, public rotint_ana(sepi, sepj, rijv, w, dw, se_int_control, se_taper)
calculates the derivative of the two-particle interactions
Definition: semi_empirical_int_ana.F:1423
semi_empirical_int_ana::corecore_el_ana
recursive subroutine, public corecore_el_ana(sepi, sepj, rijv, itype, enuc, denuc, se_int_control, se_taper)
Computes analytical gradients for semiempirical core-core electrostatic interaction only.
Definition: semi_empirical_int_ana.F:851
semi_empirical_int_ana::corecore_ana
recursive subroutine, public corecore_ana(sepi, sepj, rijv, itype, enuc, denuc, se_int_control, se_taper)
Computes analytical gradients for semiempirical core-core interaction.
Definition: semi_empirical_int_ana.F:542
semi_empirical_int_ana::rotnuc_ana
recursive subroutine, public rotnuc_ana(sepi, sepj, rijv, itype, e1b, e2a, de1b, de2a, se_int_control, se_taper)
Computes analytical gradients for semiempirical integrals.
Definition: semi_empirical_int_ana.F:285
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::rij_threshold
real(kind=dp), parameter, public rij_threshold
Definition: semi_empirical_int_arrays.F:27
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::ijkl_sym
integer, dimension(491), public ijkl_sym
Definition: semi_empirical_int_arrays.F:87
semi_empirical_int_arrays::inddp
integer, dimension(5, 3), public inddp
Definition: semi_empirical_int_arrays.F:91
semi_empirical_int_arrays::indpp
integer, dimension(3, 3), public indpp
Definition: semi_empirical_int_arrays.F:90
semi_empirical_int_arrays::fac_x_to_z
real(kind=dp), dimension(2, 9), parameter, public fac_x_to_z
Definition: semi_empirical_int_arrays.F:46
semi_empirical_int_arrays::map_x_to_z
integer, dimension(2, 9), parameter, public map_x_to_z
Definition: semi_empirical_int_arrays.F:35
semi_empirical_int_arrays::l_index
integer, dimension(9), parameter, public l_index
Definition: semi_empirical_int_arrays.F:80
semi_empirical_int_arrays::indexa
integer, dimension(9, 9), public indexa
Definition: semi_empirical_int_arrays.F:71
semi_empirical_int_arrays::inddd
integer, dimension(5, 5), public inddd
Definition: semi_empirical_int_arrays.F:92
semi_empirical_int_arrays::ijkl_ind
integer, dimension(45, 45), public ijkl_ind
Definition: semi_empirical_int_arrays.F:84
semi_empirical_int_num
Integrals for semi-empiric methods.
Definition: semi_empirical_int_num.F:17
semi_empirical_int_num::terep_d_num
subroutine, public terep_d_num(sepi, sepj, rij, rep, se_int_control, se_int_screen, ft)
Calculates the two-electron repulsion integrals - d shell only.
Definition: semi_empirical_int_num.F:330
semi_empirical_int_num::terep_sp_num
subroutine, public terep_sp_num(sepi, sepj, rij, rep, se_int_control, se_int_screen, ft)
Calculates the two-electron repulsion integrals - sp shell only.
Definition: semi_empirical_int_num.F:269
semi_empirical_int_num::nucint_d_num
subroutine, public nucint_d_num(sepi, sepj, rij, core, itype, se_int_control, se_int_screen)
Calculates the nuclear attraction integrals involving d orbitals.
Definition: semi_empirical_int_num.F:920
semi_empirical_int_num::nucint_sp_num
subroutine, public nucint_sp_num(sepi, sepj, rij, ssss, core, itype, se_int_control, se_int_screen)
...
Definition: semi_empirical_int_num.F:839
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::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::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
semi_empirical_types::rotmat_release
subroutine, public rotmat_release(rotmat)
Releases rotmat type.
Definition: semi_empirical_types.F:534
semi_empirical_types::setup_se_int_control_type
subroutine, public setup_se_int_control_type(se_int_control, shortrange, do_ewald_r3, do_ewald_gks, integral_screening, max_multipole, pc_coulomb_int)
Setup the Semiempirical integral control type.
Definition: semi_empirical_types.F:555
semi_empirical_types::rotmat_create
subroutine, public rotmat_create(rotmat)
Creates rotmat type.
Definition: semi_empirical_types.F:522
taper_types
Definition of the semi empirical parameter types.
Definition: taper_types.F:12
taper_types::dtaper_eval
real(kind=dp) function, public dtaper_eval(taper, rij)
Analytical derivatives for taper function.
Definition: taper_types.F:99
taper_types::taper_eval
real(kind=dp) function, public taper_eval(taper, rij)
Taper functions.
Definition: taper_types.F:79
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