(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
18MODULE 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! **************************************************************************************************
66INTERFACE
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
77END 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! **************************************************************************************************
87INTERFACE
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
101END 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! **************************************************************************************************
111INTERFACE
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
126END 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! **************************************************************************************************
136INTERFACE
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
152END 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! **************************************************************************************************
162INTERFACE
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
179END 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! **************************************************************************************************
189INTERFACE
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
205END 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! **************************************************************************************************
215INTERFACE
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
229END 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! **************************************************************************************************
239INTERFACE
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
253END 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
259CONTAINS
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
1962END MODULE semi_empirical_int_ana
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
rot_2el_2c_first_debug
subroutine rot_2el_2c_first_debug(sepi, sepj, rijv, se_int_control, se_taper, invert, ii, kk, v_d)
Check Numerical Vs Analytical.
Definition semi_empirical_int_debug.F:152
check_dssss_nucint_ana
subroutine check_dssss_nucint_ana(sepi, sepj, r, dssss, itype, se_int_control, se_taper)
Check Numerical Vs Analytical ssss.
Definition semi_empirical_int_debug.F:253
check_dcorecore_ana
subroutine check_dcorecore_ana(sepi, sepj, rijv, itype, se_int_control, se_taper, enuc, denuc)
Check Numerical Vs Analytical CORE-CORE.
Definition semi_empirical_int_debug.F:510
check_rotint_ana
subroutine check_rotint_ana(sepi, sepj, rijv, w, dw, se_int_control, se_taper)
Check Numerical Vs Analytical ROTINT_ANA.
Definition semi_empirical_int_debug.F:644
check_drotnuc_ana
subroutine check_drotnuc_ana(sepi, sepj, rijv, itype, se_int_control, se_taper, e1b, e2a, de1b, de2a)
Check Numerical Vs Analytical.
Definition semi_empirical_int_debug.F:405
check_dcore_nucint_ana
subroutine check_dcore_nucint_ana(sepi, sepj, r, dcore, itype, se_int_control, se_taper)
Check Numerical Vs Analytical core.
Definition semi_empirical_int_debug.F:305
check_dterep_ana
subroutine check_dterep_ana(sepi, sepj, r, ri, dri, se_int_control, se_taper, lgrad)
Check Numerical Vs Analytical DTEREP_ANA.
Definition semi_empirical_int_debug.F:574
semi_empirical_types::rotmat_type
Rotation Matrix Type.
Definition semi_empirical_types.F:111
semi_empirical_types::se_int_control_type
Definition semi_empirical_types.F:136
semi_empirical_types::se_int_screen_type
Store the value of the tapering function and possibly its derivative for screened integrals.
Definition semi_empirical_types.F:150
semi_empirical_types::se_taper_type
Taper type use in semi-empirical calculations.
Definition semi_empirical_types.F:157
semi_empirical_types::semi_empirical_type
Semi-empirical type.
Definition semi_empirical_types.F:57