(git:374b731)
Loading...
Searching...
No Matches
semi_empirical_int_num.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 Integrals for semi-empiric methods
10!> \par History
11!> JGH (27.10.2004) : separate routine for nuclear attraction integrals
12!> JGH (13.03.2006) : tapering function
13!> Teodoro Laino (03.2008) [tlaino] - University of Zurich : new driver
14!> for computing integrals
15!> \author JGH (11.10.2004)
16! **************************************************************************************************
17MODULE semi_empirical_int_num
18
19 USE input_constants, ONLY: do_method_am1,&
20 do_method_pchg,&
21 do_method_pdg,&
22 do_method_pm3,&
23 do_method_pm6,&
24 do_method_pm6fm,&
25 do_method_undef,&
26 do_se_is_kdso_d
27 USE kinds, ONLY: dp
28 USE multipole_types, ONLY: do_multipole_none
29 USE physcon, ONLY: angstrom,&
30 evolt
31 USE semi_empirical_int_arrays, ONLY: &
32 ijkl_ind, ijkl_sym, inddd, inddp, indexa, indexb, indpp, int2c_type, l_index, rij_threshold
33 USE semi_empirical_int_utils, ONLY: ijkl_d,&
34 ijkl_sp,&
35 rot_2el_2c_first,&
36 rotmat,&
37 store_2el_2c_diag
38 USE semi_empirical_types, ONLY: rotmat_create,&
39 rotmat_release,&
40 rotmat_type,&
41 se_int_control_type,&
42 se_int_screen_type,&
43 se_taper_type,&
44 semi_empirical_type,&
45 setup_se_int_control_type
46 USE taper_types, ONLY: taper_eval
47#include "./base/base_uses.f90"
48
49 IMPLICIT NONE
50
51 PRIVATE
52
53 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
54 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_int_num'
55 PUBLIC :: rotint_num, rotnuc_num, corecore_num, &
56 drotint_num, drotnuc_num, dcorecore_num, &
57 ssss_nucint_num, core_nucint_num, terep_num, &
58 nucint_sp_num, terep_sp_num, terep_d_num, &
59 nucint_d_num, corecore_el_num, dcorecore_el_num
60
61CONTAINS
62
63! **************************************************************************************************
64!> \brief Computes the two particle interactions in the lab frame
65!> \param sepi Atomic parameters of first atom
66!> \param sepj Atomic parameters of second atom
67!> \param rijv Coordinate vector i -> j
68!> \param w Array of two-electron repulsion integrals.
69!> \param se_int_control input parameters that control the calculation of SE
70!> integrals (shortrange, R3 residual, screening type)
71!> \param se_taper ...
72!> \note routine adapted from mopac7 (rotate)
73!> written by Ernest R. Davidson, Indiana University.
74!> Teodoro Laino [tlaino] - University of Zurich 04.2008 : major rewriting
75! **************************************************************************************************
76 SUBROUTINE rotint_num(sepi, sepj, rijv, w, se_int_control, se_taper)
77 TYPE(semi_empirical_type), POINTER :: sepi, sepj
78 REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
79 REAL(dp), DIMENSION(2025), INTENT(OUT) :: w
80 TYPE(se_int_control_type), INTENT(IN) :: se_int_control
81 TYPE(se_taper_type), POINTER :: se_taper
82
83 INTEGER :: i, i1, ii, ij, ij1, iminus, istep, &
84 iw_loc, j, j1, jj, k, kk, kl, l, &
85 limij, limkl, mm
86 LOGICAL, DIMENSION(45, 45) :: logv
87 REAL(dp) :: rij
88 REAL(kind=dp) :: cc, wrepp
89 REAL(kind=dp), DIMENSION(2025) :: ww
90 REAL(kind=dp), DIMENSION(45, 45) :: v
91 REAL(kind=dp), DIMENSION(491) :: rep
92 TYPE(rotmat_type), POINTER :: ij_matrix
93
94 NULLIFY (ij_matrix)
95 rij = dot_product(rijv, rijv)
96 IF (rij > rij_threshold) THEN
97 ! The repulsion integrals over molecular frame (w) are stored in the
98 ! order in which they will later be used. ie. (i,j/k,l) where
99 ! j.le.i and l.le.k and l varies most rapidly and i least
100 ! rapidly. (anti-normal computer storage)
101 rij = sqrt(rij)
102
103 ! Create the rotation matrix
104 CALL rotmat_create(ij_matrix)
105 CALL rotmat(sepi, sepj, rijv, rij, ij_matrix, do_derivatives=.false.)
106
107 ! Compute Integrals in Diatomic Frame
108 CALL terep_num(sepi, sepj, rij, rep, se_taper=se_taper, se_int_control=se_int_control)
109
110 ! Rotate Integrals
111 ii = sepi%natorb
112 kk = sepj%natorb
113 IF (ii*kk > 0) THEN
114 limij = sepi%atm_int_size
115 limkl = sepj%atm_int_size
116 istep = limkl*limij
117 DO i1 = 1, istep
118 ww(i1) = 0.0_dp
119 END DO
120
121 ! First step in rotation of integrals
122 CALL rot_2el_2c_first(sepi, sepj, rijv, se_int_control, se_taper, .false., ii, kk, rep, &
123 logv, ij_matrix, v, lgrad=.false.)
124
125 ! Second step in rotation of integrals
126 DO i1 = 1, ii
127 DO j1 = 1, i1
128 ij = indexa(i1, j1)
129 jj = indexb(i1, j1)
130 mm = int2c_type(jj)
131 DO k = 1, kk
132 DO l = 1, k
133 kl = indexb(k, l)
134 IF (logv(ij, kl)) THEN
135 wrepp = v(ij, kl)
136 SELECT CASE (mm)
137 CASE (1)
138 ! (SS/)
139 i = 1
140 j = 1
141 iw_loc = (indexb(i, j) - 1)*limkl + kl
142 ww(iw_loc) = wrepp
143 CASE (2)
144 ! (SP/)
145 j = 1
146 DO i = 1, 3
147 iw_loc = (indexb(i + 1, j) - 1)*limkl + kl
148 ww(iw_loc) = ww(iw_loc) + ij_matrix%sp(i1 - 1, i)*wrepp
149 END DO
150 CASE (3)
151 ! (PP/)
152 DO i = 1, 3
153 cc = ij_matrix%pp(i, i1 - 1, j1 - 1)
154 iw_loc = (indexb(i + 1, i + 1) - 1)*limkl + kl
155 ww(iw_loc) = ww(iw_loc) + cc*wrepp
156 iminus = i - 1
157 IF (iminus /= 0) THEN
158 DO j = 1, iminus
159 cc = ij_matrix%pp(1 + i + j, i1 - 1, j1 - 1)
160 iw_loc = (indexb(i + 1, j + 1) - 1)*limkl + kl
161 ww(iw_loc) = ww(iw_loc) + cc*wrepp
162 END DO
163 END IF
164 END DO
165 CASE (4)
166 ! (SD/)
167 j = 1
168 DO i = 1, 5
169 iw_loc = (indexb(i + 4, j) - 1)*limkl + kl
170 ww(iw_loc) = ww(iw_loc) + ij_matrix%sd(i1 - 4, i)*wrepp
171 END DO
172 CASE (5)
173 ! (DP/)
174 DO i = 1, 5
175 DO j = 1, 3
176 iw_loc = (indexb(i + 4, j + 1) - 1)*limkl + kl
177 ij1 = 3*(i - 1) + j
178 ww(iw_loc) = ww(iw_loc) + ij_matrix%pd(ij1, i1 - 4, j1 - 1)*wrepp
179 END DO
180 END DO
181 CASE (6)
182 ! (DD/)
183 DO i = 1, 5
184 cc = ij_matrix%dd(i, i1 - 4, j1 - 4)
185 iw_loc = (indexb(i + 4, i + 4) - 1)*limkl + kl
186 ww(iw_loc) = ww(iw_loc) + cc*wrepp
187 iminus = i - 1
188 IF (iminus /= 0) THEN
189 DO j = 1, iminus
190 ij1 = inddd(i, j)
191 cc = ij_matrix%dd(ij1, i1 - 4, j1 - 4)
192 iw_loc = (indexb(i + 4, j + 4) - 1)*limkl + kl
193 ww(iw_loc) = ww(iw_loc) + cc*wrepp
194 END DO
195 END IF
196 END DO
197 END SELECT
198 END IF
199 END DO
200 END DO
201 END DO
202 END DO
203 END IF
204 CALL rotmat_release(ij_matrix)
205
206 ! Store two electron integrals in the triangular format
207 CALL store_2el_2c_diag(limij, limkl, ww, w)
208 END IF
209 END SUBROUTINE rotint_num
210
211! **************************************************************************************************
212!> \brief Calculates the derivative pf two-electron repulsion integrals and the
213!> nuclear attraction integrals w.r.t. |r|
214!> \param sepi parameters of atom i
215!> \param sepj parameters of atom j
216!> \param rij interatomic distance
217!> \param rep array of two-electron repulsion integrals
218!> \param se_taper ...
219!> \param se_int_control input parameters that control the calculation of SE
220!> integrals (shortrange, R3 residual, screening type)
221!> \par History
222!> 03.2008 created [tlaino]
223!> \author Teodoro Laino [tlaino] - Zurich University
224! **************************************************************************************************
225 SUBROUTINE terep_num(sepi, sepj, rij, rep, se_taper, se_int_control)
226 TYPE(semi_empirical_type), POINTER :: sepi, sepj
227 REAL(dp), INTENT(IN) :: rij
228 REAL(dp), DIMENSION(491), INTENT(OUT) :: rep
229 TYPE(se_taper_type), POINTER :: se_taper
230 TYPE(se_int_control_type), INTENT(IN) :: se_int_control
231
232 REAL(kind=dp) :: ft
233 TYPE(se_int_screen_type) :: se_int_screen
234
235 ft = taper_eval(se_taper%taper, rij)
236 ! In case of dumped integrals compute an additional taper term
237 IF (se_int_control%integral_screening == do_se_is_kdso_d) THEN
238 se_int_screen%ft = taper_eval(se_taper%taper_add, rij)
239 END IF
240
241 ! Contribution from sp shells
242 CALL terep_sp_num(sepi, sepj, rij, rep, se_int_control, se_int_screen, ft)
243
244 IF (sepi%dorb .OR. sepj%dorb) THEN
245 ! Compute the contribution from d shells
246 CALL terep_d_num(sepi, sepj, rij, rep, se_int_control, se_int_screen, &
247 ft)
248 END IF
249
250 END SUBROUTINE terep_num
251
252! **************************************************************************************************
253!> \brief Calculates the two-electron repulsion integrals - sp shell only
254!> \param sepi parameters of atom i
255!> \param sepj parameters of atom j
256!> \param rij ...
257!> \param rep array of two-electron repulsion integrals
258!> \param se_int_control input parameters that control the calculation of SE
259!> integrals (shortrange, R3 residual, screening type)
260!> \param se_int_screen contains information for computing the screened
261!> integrals KDSO-D
262!> \param ft ...
263!> \par History
264!> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
265!> for computing integrals
266! **************************************************************************************************
267 SUBROUTINE terep_sp_num(sepi, sepj, rij, rep, se_int_control, se_int_screen, &
268 ft)
269 TYPE(semi_empirical_type), POINTER :: sepi, sepj
270 REAL(dp), INTENT(IN) :: rij
271 REAL(dp), DIMENSION(491), INTENT(OUT) :: rep
272 TYPE(se_int_control_type), INTENT(IN) :: se_int_control
273 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen
274 REAL(dp), INTENT(IN) :: ft
275
276 INTEGER :: i, ij, j, k, kl, l, lasti, lastj, li, &
277 lj, lk, ll, nold, numb
278 REAL(kind=dp) :: tmp
279
280 lasti = sepi%natorb
281 lastj = sepj%natorb
282 DO i = 1, min(lasti, 4)
283 li = l_index(i)
284 DO j = 1, i
285 lj = l_index(j)
286 ij = indexa(i, j)
287 DO k = 1, min(lastj, 4)
288 lk = l_index(k)
289 DO l = 1, k
290 ll = l_index(l)
291 kl = indexa(k, l)
292
293 numb = ijkl_ind(ij, kl)
294 IF (numb > 0) THEN
295 nold = ijkl_sym(numb)
296 IF (nold > 0) THEN
297 rep(numb) = rep(nold)
298 ELSE IF (nold < 0) THEN
299 rep(numb) = -rep(-nold)
300 ELSE IF (nold == 0) THEN
301 tmp = ijkl_sp(sepi, sepj, ij, kl, li, lj, lk, ll, 0, rij, &
302 se_int_control, se_int_screen, do_method_undef) &
303 *ft
304 rep(numb) = tmp
305 END IF
306 END IF
307 END DO
308 END DO
309 END DO
310 END DO
311 END SUBROUTINE terep_sp_num
312
313! **************************************************************************************************
314!> \brief Calculates the two-electron repulsion integrals - d shell only
315!> \param sepi ...
316!> \param sepj ...
317!> \param rij interatomic distance
318!> \param rep array of two-electron repulsion integrals
319!> \param se_int_control input parameters that control the calculation of SE
320!> integrals (shortrange, R3 residual, screening type)
321!> \param se_int_screen contains information for computing the screened
322!> integrals KDSO-D
323!> \param ft ...
324!> \par History
325!> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
326!> for computing integrals
327! **************************************************************************************************
328 SUBROUTINE terep_d_num(sepi, sepj, rij, rep, se_int_control, se_int_screen, &
329 ft)
330 TYPE(semi_empirical_type), POINTER :: sepi, sepj
331 REAL(dp), INTENT(IN) :: rij
332 REAL(dp), DIMENSION(491), INTENT(INOUT) :: rep
333 TYPE(se_int_control_type), INTENT(IN) :: se_int_control
334 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen
335 REAL(dp), INTENT(IN) :: ft
336
337 INTEGER :: i, ij, j, k, kl, l, lasti, lastj, li, &
338 lj, lk, ll, nold, numb
339 REAL(kind=dp) :: tmp
340
341 lasti = sepi%natorb
342 lastj = sepj%natorb
343 DO i = 1, lasti
344 li = l_index(i)
345 DO j = 1, i
346 lj = l_index(j)
347 ij = indexa(i, j)
348 DO k = 1, lastj
349 lk = l_index(k)
350 DO l = 1, k
351 ll = l_index(l)
352 kl = indexa(k, l)
353
354 numb = ijkl_ind(ij, kl)
355 ! From 1 to 34 we store integrals involving sp shells
356 IF (numb > 34) THEN
357 nold = ijkl_sym(numb)
358 IF (nold > 34) THEN
359 rep(numb) = rep(nold)
360 ELSE IF (nold < -34) THEN
361 rep(numb) = -rep(-nold)
362 ELSE IF (nold == 0) THEN
363 tmp = ijkl_d(sepi, sepj, ij, kl, li, lj, lk, ll, 0, rij, &
364 se_int_control, se_int_screen, do_method_undef) &
365 *ft
366 rep(numb) = tmp
367 END IF
368 END IF
369 END DO
370 END DO
371 END DO
372 END DO
373
374 END SUBROUTINE terep_d_num
375
376! **************************************************************************************************
377!> \brief Computes the two-particle interactions.
378!> \param sepi Atomic parameters of first atom
379!> \param sepj Atomic parameters of second atom
380!> \param rijv Coordinate vector i -> j
381!> \param e1b Array of electron-nuclear attraction integrals,
382!> e1b = Electron on atom ni attracting nucleus of nj.
383!> \param e2a Array of electron-nuclear attraction integrals,
384!> e2a = Electron on atom nj attracting nucleus of ni.
385!> \param itype ...
386!> \param se_int_control ...
387!> \param se_taper ...
388!> \note routine adapted from mopac7 (rotate)
389!> written by Ernest R. Davidson, Indiana University.
390!> Teodoro Laino [tlaino] - University of Zurich 04.2008 : major rewriting
391!> Teodoro Laino [tlaino] - University of Zurich 04.2008 : removed the core-core part
392! **************************************************************************************************
393 SUBROUTINE rotnuc_num(sepi, sepj, rijv, e1b, e2a, itype, se_int_control, se_taper)
394 TYPE(semi_empirical_type), POINTER :: sepi, sepj
395 REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
396 REAL(dp), DIMENSION(45), INTENT(OUT), OPTIONAL :: e1b, e2a
397 INTEGER, INTENT(IN) :: itype
398 TYPE(se_int_control_type), INTENT(IN) :: se_int_control
399 TYPE(se_taper_type), POINTER :: se_taper
400
401 INTEGER :: i, idd, idp, ind1, ind2, ipp, j, &
402 last_orbital(2), m, n
403 LOGICAL :: l_e1b, l_e2a, task(2)
404 REAL(dp) :: rij
405 REAL(dp), DIMENSION(10, 2) :: core
406 REAL(dp), DIMENSION(45) :: tmp
407 TYPE(rotmat_type), POINTER :: ij_matrix
408
409 NULLIFY (ij_matrix)
410 l_e1b = PRESENT(e1b)
411 l_e2a = PRESENT(e2a)
412 rij = dot_product(rijv, rijv)
413
414 IF (rij > rij_threshold) THEN
415 rij = sqrt(rij)
416 ! Create the rotation matrix
417 CALL rotmat_create(ij_matrix)
418 CALL rotmat(sepi, sepj, rijv, rij, ij_matrix, do_derivatives=.false.)
419
420 ! Compute Integrals in Diatomic Frame
421 CALL core_nucint_num(sepi, sepj, rij, core=core, itype=itype, se_taper=se_taper, &
422 se_int_control=se_int_control)
423
424 ! Copy parameters over to arrays for do loop.
425 last_orbital(1) = sepi%natorb
426 last_orbital(2) = sepj%natorb
427 task(1) = l_e1b
428 task(2) = l_e2a
429 DO n = 1, 2
430 IF (.NOT. task(n)) cycle
431 DO i = 1, last_orbital(n)
432 ind1 = i - 1
433 DO j = 1, i
434 ind2 = j - 1
435 m = (i*(i - 1))/2 + j
436 ! Perform Rotations ...
437 IF (ind2 == 0) THEN
438 IF (ind1 == 0) THEN
439 ! Type of Integral (SS/)
440 tmp(m) = core(1, n)
441 ELSE IF (ind1 < 4) THEN
442 ! Type of Integral (SP/)
443 tmp(m) = ij_matrix%sp(1, ind1)*core(2, n)
444 ELSE
445 ! Type of Integral (SD/)
446 tmp(m) = ij_matrix%sd(1, ind1 - 3)*core(5, n)
447 END IF
448 ELSE IF (ind2 < 4) THEN
449 IF (ind1 < 4) THEN
450 ! Type of Integral (PP/)
451 ipp = indpp(ind1, ind2)
452 tmp(m) = core(3, n)*ij_matrix%pp(ipp, 1, 1) + &
453 core(4, n)*(ij_matrix%pp(ipp, 2, 2) + ij_matrix%pp(ipp, 3, 3))
454 ELSE
455 ! Type of Integral (PD/)
456 idp = inddp(ind1 - 3, ind2)
457 tmp(m) = core(6, n)*ij_matrix%pd(idp, 1, 1) + &
458 core(8, n)*(ij_matrix%pd(idp, 2, 2) + ij_matrix%pd(idp, 3, 3))
459 END IF
460 ELSE
461 ! Type of Integral (DD/)
462 idd = inddd(ind1 - 3, ind2 - 3)
463 tmp(m) = core(7, n)*ij_matrix%dd(idd, 1, 1) + &
464 core(9, n)*(ij_matrix%dd(idd, 2, 2) + ij_matrix%dd(idd, 3, 3)) + &
465 core(10, n)*(ij_matrix%dd(idd, 4, 4) + ij_matrix%dd(idd, 5, 5))
466 END IF
467 END DO
468 END DO
469 IF (n == 1) THEN
470 DO i = 1, sepi%atm_int_size
471 e1b(i) = -tmp(i)
472 END DO
473 END IF
474 IF (n == 2) THEN
475 DO i = 1, sepj%atm_int_size
476 e2a(i) = -tmp(i)
477 END DO
478 END IF
479 END DO
480 CALL rotmat_release(ij_matrix)
481 END IF
482 END SUBROUTINE rotnuc_num
483
484! **************************************************************************************************
485!> \brief Computes the core-core interactions.
486!> \param sepi Atomic parameters of first atom
487!> \param sepj Atomic parameters of second atom
488!> \param rijv Coordinate vector i -> j
489!> \param enuc nuclear-nuclear repulsion term.
490!> \param itype ...
491!> \param se_int_control input parameters that control the calculation of SE
492!> integrals (shortrange, R3 residual, screening type)
493!> \param se_taper ...
494!> \note routine adapted from mopac7 (rotate)
495!> written by Ernest R. Davidson, Indiana University.
496!> Teodoro Laino [tlaino] - University of Zurich 04.2008 : major rewriting
497!> Teodoro Laino [tlaino] - University of Zurich 04.2008 : splitted from rotnuc
498! **************************************************************************************************
499 SUBROUTINE corecore_num(sepi, sepj, rijv, enuc, itype, se_int_control, se_taper)
500 TYPE(semi_empirical_type), POINTER :: sepi, sepj
501 REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
502 REAL(dp), INTENT(OUT) :: enuc
503 INTEGER, INTENT(IN) :: itype
504 TYPE(se_int_control_type), INTENT(IN) :: se_int_control
505 TYPE(se_taper_type), POINTER :: se_taper
506
507 INTEGER :: ig, nt
508 REAL(dp) :: aab, alpi, alpj, apdg, ax, dai, daj, &
509 dbi, dbj, pai, paj, pbi, pbj, qcorr, &
510 rij, rija, scale, ssss, ssss_sr, tmp, &
511 xab, zaf, zbf, zz
512 REAL(dp), DIMENSION(4) :: fni1, fni2, fni3, fnj1, fnj2, fnj3
513 TYPE(se_int_control_type) :: se_int_control_off
514
515 rij = dot_product(rijv, rijv)
516 enuc = 0.0_dp
517 IF (rij > rij_threshold) THEN
518 rij = sqrt(rij)
519 CALL setup_se_int_control_type(se_int_control_off, shortrange=.false., do_ewald_r3=.false., &
520 do_ewald_gks=.false., integral_screening=se_int_control%integral_screening, &
521 max_multipole=do_multipole_none, pc_coulomb_int=.false.)
522 CALL ssss_nucint_num(sepi, sepj, rij, ssss=ssss, itype=itype, se_taper=se_taper, &
523 se_int_control=se_int_control_off)
524 ! In case let's compute the short-range part of the (ss|ss) integral
525 IF (se_int_control%shortrange) THEN
526 CALL ssss_nucint_num(sepi, sepj, rij, ssss=ssss_sr, itype=itype, se_taper=se_taper, &
527 se_int_control=se_int_control)
528 ELSE
529 ssss_sr = ssss
530 END IF
531 zz = sepi%zeff*sepj%zeff
532 ! Nuclear-Nuclear electrostatic contribution
533 enuc = zz*ssss_sr
534 ! Method dependent correction terms
535 IF (itype /= do_method_pm6 .AND. itype /= do_method_pm6fm) THEN
536 alpi = sepi%alp
537 alpj = sepj%alp
538 scale = exp(-alpi*rij) + exp(-alpj*rij)
539
540 nt = sepi%z + sepj%z
541 IF (nt == 8 .OR. nt == 9) THEN
542 IF (sepi%z == 7 .OR. sepi%z == 8) scale = scale + (angstrom*rij - 1._dp)*exp(-alpi*rij)
543 IF (sepj%z == 7 .OR. sepj%z == 8) scale = scale + (angstrom*rij - 1._dp)*exp(-alpj*rij)
544 END IF
545 scale = abs(scale*zz*ssss)
546 zz = zz/rij
547 IF (itype == do_method_am1 .OR. itype == do_method_pm3 .OR. itype == do_method_pdg) THEN
548 IF (itype == do_method_am1 .AND. sepi%z == 5) THEN
549 !special case AM1 Boron
550 SELECT CASE (sepj%z)
551 CASE DEFAULT
552 nt = 1
553 CASE (1)
554 nt = 2
555 CASE (6)
556 nt = 3
557 CASE (9, 17, 35, 53)
558 nt = 4
559 END SELECT
560 fni1(:) = sepi%bfn1(:, nt)
561 fni2(:) = sepi%bfn2(:, nt)
562 fni3(:) = sepi%bfn3(:, nt)
563 ELSE
564 fni1(:) = sepi%fn1(:)
565 fni2(:) = sepi%fn2(:)
566 fni3(:) = sepi%fn3(:)
567 END IF
568 IF (itype == do_method_am1 .AND. sepj%z == 5) THEN
569 !special case AM1 Boron
570 SELECT CASE (sepi%z)
571 CASE DEFAULT
572 nt = 1
573 CASE (1)
574 nt = 2
575 CASE (6)
576 nt = 3
577 CASE (9, 17, 35, 53)
578 nt = 4
579 END SELECT
580 fnj1(:) = sepj%bfn1(:, nt)
581 fnj2(:) = sepj%bfn2(:, nt)
582 fnj3(:) = sepj%bfn3(:, nt)
583 ELSE
584 fnj1(:) = sepj%fn1(:)
585 fnj2(:) = sepj%fn2(:)
586 fnj3(:) = sepj%fn3(:)
587 END IF
588 ! AM1/PM3/PDG correction to nuclear repulsion
589 DO ig = 1, SIZE(fni1)
590 IF (abs(fni1(ig)) > 0._dp) THEN
591 ax = fni2(ig)*(rij - fni3(ig))**2
592 IF (ax <= 25._dp) THEN
593 scale = scale + zz*fni1(ig)*exp(-ax)
594 END IF
595 END IF
596 IF (abs(fnj1(ig)) > 0._dp) THEN
597 ax = fnj2(ig)*(rij - fnj3(ig))**2
598 IF (ax <= 25._dp) THEN
599 scale = scale + zz*fnj1(ig)*exp(-ax)
600 END IF
601 END IF
602 END DO
603 END IF
604 IF (itype == do_method_pdg) THEN
605 ! PDDG function
606 zaf = sepi%zeff/nt
607 zbf = sepj%zeff/nt
608 pai = sepi%pre(1)
609 pbi = sepi%pre(2)
610 paj = sepj%pre(1)
611 pbj = sepj%pre(2)
612 dai = sepi%d(1)
613 dbi = sepi%d(2)
614 daj = sepj%d(1)
615 dbj = sepj%d(2)
616 apdg = 10._dp*angstrom**2
617 qcorr = &
618 (zaf*pai + zbf*paj)*exp(-apdg*(rij - dai - daj)**2) + &
619 (zaf*pai + zbf*pbj)*exp(-apdg*(rij - dai - dbj)**2) + &
620 (zaf*pbi + zbf*paj)*exp(-apdg*(rij - dbi - daj)**2) + &
621 (zaf*pbi + zbf*pbj)*exp(-apdg*(rij - dbi - dbj)**2)
622 ELSEIF (itype == do_method_pchg) THEN
623 qcorr = 0.0_dp
624 scale = 0.0_dp
625 ELSE
626 qcorr = 0.0_dp
627 END IF
628 ELSE
629 rija = rij*angstrom
630 scale = zz*ssss
631 ! PM6 core-core terms
632 xab = sepi%xab(sepj%z)
633 aab = sepi%aab(sepj%z)
634 IF ((sepi%z == 1 .AND. (sepj%z == 6 .OR. sepj%z == 7 .OR. sepj%z == 8)) .OR. &
635 (sepj%z == 1 .AND. (sepi%z == 6 .OR. sepi%z == 7 .OR. sepi%z == 8))) THEN
636 ! Special Case O-H or N-H or C-H
637 scale = scale*(2._dp*xab*exp(-aab*rija*rija))
638 ELSEIF (sepi%z == 6 .AND. sepj%z == 6) THEN
639 ! Special Case C-C
640 scale = scale*(2._dp*xab*exp(-aab*(rija + 0.0003_dp*rija**6)) + 9.28_dp*exp(-5.98_dp*rija))
641 ELSEIF ((sepi%z == 8 .AND. sepj%z == 14) .OR. &
642 (sepj%z == 8 .AND. sepi%z == 14)) THEN
643 ! Special Case Si-O
644 scale = scale*(2._dp*xab*exp(-aab*(rija + 0.0003_dp*rija**6)) - 0.0007_dp*exp(-(rija - 2.9_dp)**2))
645 ELSE
646 ! General Case
647 ! Factor of 2 found by experiment
648 scale = scale*(2._dp*xab*exp(-aab*(rija + 0.0003_dp*rija**6)))
649 END IF
650 ! General correction term a*exp(-b*(rij-c)^2)
651 qcorr = (sepi%a*exp(-sepi%b*(rij - sepi%c)**2))*sepi%zeff*sepj%zeff/rij + &
652 (sepj%a*exp(-sepj%b*(rij - sepj%c)**2))*sepi%zeff*sepj%zeff/rij
653 ! Hard core repulsion
654 tmp = (real(sepi%z, dp)**(1._dp/3._dp) + real(sepj%z, dp)**(1._dp/3._dp))
655 qcorr = qcorr + 1.e-8_dp/evolt*(tmp/rija)**12
656 END IF
657 enuc = enuc + scale + qcorr
658 END IF
659 END SUBROUTINE corecore_num
660
661! **************************************************************************************************
662!> \brief Computes the electrostatic core-core interactions only.
663!> \param sepi Atomic parameters of first atom
664!> \param sepj Atomic parameters of second atom
665!> \param rijv Coordinate vector i -> j
666!> \param enuc nuclear-nuclear repulsion term.
667!> \param itype ...
668!> \param se_int_control input parameters that control the calculation of SE
669!> integrals (shortrange, R3 residual, screening type)
670!> \param se_taper ...
671!> \author Teodoro Laino [tlaino] - 05.2009
672! **************************************************************************************************
673 SUBROUTINE corecore_el_num(sepi, sepj, rijv, enuc, itype, se_int_control, se_taper)
674 TYPE(semi_empirical_type), POINTER :: sepi, sepj
675 REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
676 REAL(dp), INTENT(OUT) :: enuc
677 INTEGER, INTENT(IN) :: itype
678 TYPE(se_int_control_type), INTENT(IN) :: se_int_control
679 TYPE(se_taper_type), POINTER :: se_taper
680
681 REAL(dp) :: rij, ssss, ssss_sr, zz
682 TYPE(se_int_control_type) :: se_int_control_off
683
684 rij = dot_product(rijv, rijv)
685 enuc = 0.0_dp
686 IF (rij > rij_threshold) THEN
687 rij = sqrt(rij)
688 CALL setup_se_int_control_type(se_int_control_off, shortrange=.false., do_ewald_r3=.false., &
689 do_ewald_gks=.false., integral_screening=se_int_control%integral_screening, &
690 max_multipole=do_multipole_none, pc_coulomb_int=.false.)
691 CALL ssss_nucint_num(sepi, sepj, rij, ssss=ssss, itype=itype, se_taper=se_taper, &
692 se_int_control=se_int_control_off)
693 ! In case let's compute the short-range part of the (ss|ss) integral
694 IF (se_int_control%shortrange .OR. se_int_control%pc_coulomb_int) THEN
695 CALL ssss_nucint_num(sepi, sepj, rij, ssss=ssss_sr, itype=itype, se_taper=se_taper, &
696 se_int_control=se_int_control)
697 ELSE
698 ssss_sr = ssss
699 END IF
700 zz = sepi%zeff*sepj%zeff
701 ! Nuclear-Nuclear electrostatic contribution
702 enuc = zz*ssss_sr
703 END IF
704 END SUBROUTINE corecore_el_num
705
706! **************************************************************************************************
707!> \brief Calculates the SSSS integrals (main driver)
708!> \param sepi parameters of atom i
709!> \param sepj parameters of atom j
710!> \param rij interatomic distance
711!> \param ssss derivative of (ssss) integral
712!> derivatives are intended w.r.t. rij
713!> \param itype type of semi_empirical model
714!> extension to the original routine to compute qm/mm integrals
715!> \param se_taper ...
716!> \param se_int_control input parameters that control the calculation of SE
717!> integrals (shortrange, R3 residual, screening type)
718!> \par History
719!> 03.2008 created [tlaino]
720!> \author Teodoro Laino - Zurich University
721! **************************************************************************************************
722 SUBROUTINE ssss_nucint_num(sepi, sepj, rij, ssss, itype, se_taper, se_int_control)
723 TYPE(semi_empirical_type), POINTER :: sepi, sepj
724 REAL(dp), INTENT(IN) :: rij
725 REAL(dp), INTENT(OUT) :: ssss
726 INTEGER, INTENT(IN) :: itype
727 TYPE(se_taper_type), POINTER :: se_taper
728 TYPE(se_int_control_type), INTENT(IN) :: se_int_control
729
730 REAL(kind=dp) :: ft
731 TYPE(se_int_screen_type) :: se_int_screen
732
733! Computing Tapering function
734
735 ft = 1.0_dp
736 IF (itype /= do_method_pchg) THEN
737 ft = taper_eval(se_taper%taper, rij)
738 END IF
739
740 ! In case of dumped integrals compute an additional taper term
741 IF (se_int_control%integral_screening == do_se_is_kdso_d) THEN
742 se_int_screen%ft = 1.0_dp
743 IF (itype /= do_method_pchg) THEN
744 se_int_screen%ft = taper_eval(se_taper%taper_add, rij)
745 END IF
746 END IF
747
748 ! Contribution from the sp shells
749 CALL nucint_sp_num(sepi, sepj, rij, ssss=ssss, itype=itype, &
750 se_int_control=se_int_control, se_int_screen=se_int_screen)
751
752 ! Tapering the integrals
753 ssss = ft*ssss
754
755 END SUBROUTINE ssss_nucint_num
756
757! **************************************************************************************************
758!> \brief Calculates the nuclear attraction integrals CORE (main driver)
759!> \param sepi parameters of atom i
760!> \param sepj parameters of atom j
761!> \param rij interatomic distance
762!> \param core derivative of 4 X 2 array of electron-core attraction integrals
763!> derivatives are intended w.r.t. rij
764!> \param itype type of semi_empirical model
765!> extension to the original routine to compute qm/mm integrals
766!> \param se_taper ...
767!> \param se_int_control input parameters that control the calculation of SE
768!> integrals (shortrange, R3 residual, screening type)
769!> \par History
770!> 03.2008 created [tlaino]
771!> \author Teodoro Laino - Zurich University
772! **************************************************************************************************
773 SUBROUTINE core_nucint_num(sepi, sepj, rij, core, itype, se_taper, se_int_control)
774 TYPE(semi_empirical_type), POINTER :: sepi, sepj
775 REAL(dp), INTENT(IN) :: rij
776 REAL(dp), DIMENSION(10, 2), INTENT(OUT) :: core
777 INTEGER, INTENT(IN) :: itype
778 TYPE(se_taper_type), POINTER :: se_taper
779 TYPE(se_int_control_type), INTENT(IN) :: se_int_control
780
781 INTEGER :: i
782 REAL(kind=dp) :: ft
783 TYPE(se_int_screen_type) :: se_int_screen
784
785! Computing the Tapering function
786
787 ft = 1.0_dp
788 IF (itype /= do_method_pchg) THEN
789 ft = taper_eval(se_taper%taper, rij)
790 END IF
791
792 ! In case of dumped integrals compute an additional taper term
793 IF (se_int_control%integral_screening == do_se_is_kdso_d) THEN
794 se_int_screen%ft = 1.0_dp
795 IF (itype /= do_method_pchg) THEN
796 se_int_screen%ft = taper_eval(se_taper%taper_add, rij)
797 END IF
798 END IF
799
800 ! Contribution from the sp shells
801 CALL nucint_sp_num(sepi, sepj, rij, core=core, itype=itype, &
802 se_int_control=se_int_control, se_int_screen=se_int_screen)
803
804 IF (sepi%dorb .OR. sepj%dorb) THEN
805 ! Compute the contribution from d shells
806 CALL nucint_d_num(sepi, sepj, rij, core, itype, se_int_control, &
807 se_int_screen)
808 END IF
809
810 ! Tapering the integrals
811 DO i = 1, sepi%core_size
812 core(i, 1) = ft*core(i, 1)
813 END DO
814 DO i = 1, sepj%core_size
815 core(i, 2) = ft*core(i, 2)
816 END DO
817
818 END SUBROUTINE core_nucint_num
819
820! **************************************************************************************************
821!> \brief ...
822!> \param sepi ...
823!> \param sepj ...
824!> \param rij ...
825!> \param ssss ...
826!> \param core ...
827!> \param itype ...
828!> \param se_int_control ...
829!> \param se_int_screen ...
830!> \par History
831!> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
832!> for computing integrals
833!> Teodoro Laino (04.2008) [tlaino] - Totally rewritten: nothing to do with
834!> the old version
835! **************************************************************************************************
836
837 SUBROUTINE nucint_sp_num(sepi, sepj, rij, ssss, core, itype, se_int_control, &
838 se_int_screen)
839 TYPE(semi_empirical_type), POINTER :: sepi, sepj
840 REAL(dp), INTENT(IN) :: rij
841 REAL(dp), INTENT(INOUT), OPTIONAL :: ssss
842 REAL(dp), DIMENSION(10, 2), INTENT(INOUT), &
843 OPTIONAL :: core
844 INTEGER, INTENT(IN) :: itype
845 TYPE(se_int_control_type), INTENT(IN) :: se_int_control
846 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen
847
848 INTEGER :: ij, kl
849 LOGICAL :: l_core, l_ssss, si, sj
850
851 l_core = PRESENT(core)
852 l_ssss = PRESENT(ssss)
853 IF (.NOT. (l_core .OR. l_ssss)) RETURN
854 si = (sepi%natorb > 1)
855 sj = (sepj%natorb > 1)
856
857 ij = indexa(1, 1)
858 IF (l_ssss) THEN
859 ! Store the value for <S S | S S > (Used for computing the core-core interactions)
860 ssss = ijkl_sp(sepi, sepj, ij, ij, 0, 0, 0, 0, -1, rij, se_int_control, se_int_screen, itype)
861 END IF
862
863 IF (l_core) THEN
864 ! <S S | S S >
865 kl = indexa(1, 1)
866 core(1, 1) = ijkl_sp(sepi, sepj, kl, ij, 0, 0, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
867 IF (si) THEN
868 ! <S P | S S >
869 kl = indexa(2, 1)
870 core(2, 1) = ijkl_sp(sepi, sepj, kl, ij, 0, 1, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
871 ! <P P | S S >
872 kl = indexa(2, 2)
873 core(3, 1) = ijkl_sp(sepi, sepj, kl, ij, 1, 1, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
874 ! <P+ P+ | S S >
875 kl = indexa(3, 3)
876 core(4, 1) = ijkl_sp(sepi, sepj, kl, ij, 1, 1, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
877 END IF
878
879 ! <S S | S S >
880 kl = indexa(1, 1)
881 core(1, 2) = ijkl_sp(sepi, sepj, ij, kl, 0, 0, 0, 0, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
882 IF (sj) THEN
883 ! <S S | S P >
884 kl = indexa(2, 1)
885 core(2, 2) = ijkl_sp(sepi, sepj, ij, kl, 0, 0, 0, 1, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
886 ! <S S | P P >
887 kl = indexa(2, 2)
888 core(3, 2) = ijkl_sp(sepi, sepj, ij, kl, 0, 0, 1, 1, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
889 ! <S S | P+ P+ >
890 kl = indexa(3, 3)
891 core(4, 2) = ijkl_sp(sepi, sepj, ij, kl, 0, 0, 1, 1, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
892 END IF
893 END IF
894
895 END SUBROUTINE nucint_sp_num
896
897! **************************************************************************************************
898!> \brief Calculates the nuclear attraction integrals involving d orbitals
899!> \param sepi parameters of atom i
900!> \param sepj parameters of atom j
901!> \param rij interatomic distance
902!> \param core 4 X 2 array of electron-core attraction integrals
903!>
904!> The storage of the nuclear attraction integrals core(kl/ij) iS
905!> (SS/)=1, (SP/)=2, (PP/)=3, (P+P+/)=4, (SD/)=5,
906!> (DP/)=6, (DD/)=7, (D+P+)=8, (D+D+/)=9, (D#D#)=10
907!>
908!> where ij=1 if the orbitals centred on atom i, =2 if on atom j.
909!> \param itype type of semi_empirical model
910!> extension to the original routine to compute qm/mm integrals
911!> \param se_int_control input parameters that control the calculation of SE
912!> integrals (shortrange, R3 residual, screening type)
913!> \param se_int_screen contains information for computing the screened
914!> integrals KDSO-D
915!> \author
916!> Teodoro Laino (03.2008) [tlaino] - University of Zurich
917! **************************************************************************************************
918 SUBROUTINE nucint_d_num(sepi, sepj, rij, core, itype, se_int_control, &
919 se_int_screen)
920 TYPE(semi_empirical_type), POINTER :: sepi, sepj
921 REAL(dp), INTENT(IN) :: rij
922 REAL(dp), DIMENSION(10, 2), INTENT(INOUT) :: core
923 INTEGER, INTENT(IN) :: itype
924 TYPE(se_int_control_type), INTENT(IN) :: se_int_control
925 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen
926
927 INTEGER :: ij, kl
928
929! Check if d-orbitals are present
930
931 IF (sepi%dorb .OR. sepj%dorb) THEN
932 ij = indexa(1, 1)
933 IF (sepj%dorb) THEN
934 ! <S S | D S>
935 kl = indexa(5, 1)
936 core(5, 2) = ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 0, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
937 ! <S S | D P >
938 kl = indexa(5, 2)
939 core(6, 2) = ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 1, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
940 ! <S S | D D >
941 kl = indexa(5, 5)
942 core(7, 2) = ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
943 ! <S S | D+P+>
944 kl = indexa(6, 3)
945 core(8, 2) = ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 1, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
946 ! <S S | D+D+>
947 kl = indexa(6, 6)
948 core(9, 2) = ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
949 ! <S S | D#D#>
950 kl = indexa(8, 8)
951 core(10, 2) = ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
952 END IF
953 IF (sepi%dorb) THEN
954 ! <D S | S S>
955 kl = indexa(5, 1)
956 core(5, 1) = ijkl_d(sepi, sepj, kl, ij, 2, 0, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
957 ! <D P | S S >
958 kl = indexa(5, 2)
959 core(6, 1) = ijkl_d(sepi, sepj, kl, ij, 2, 1, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
960 ! <D D | S S >
961 kl = indexa(5, 5)
962 core(7, 1) = ijkl_d(sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
963 ! <D+P+| S S >
964 kl = indexa(6, 3)
965 core(8, 1) = ijkl_d(sepi, sepj, kl, ij, 2, 1, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
966 ! <D+D+| S S >
967 kl = indexa(6, 6)
968 core(9, 1) = ijkl_d(sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
969 ! <D#D#| S S >
970 kl = indexa(8, 8)
971 core(10, 1) = ijkl_d(sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
972 END IF
973
974 END IF
975 END SUBROUTINE nucint_d_num
976
977! **************************************************************************************************
978!> \brief Numerical Derivatives for rotint
979!> \param sepi ...
980!> \param sepj ...
981!> \param r ...
982!> \param dw ...
983!> \param delta ...
984!> \param se_int_control ...
985!> \param se_taper ...
986! **************************************************************************************************
987 SUBROUTINE drotint_num(sepi, sepj, r, dw, delta, se_int_control, se_taper)
988 TYPE(semi_empirical_type), POINTER :: sepi, sepj
989 REAL(dp), DIMENSION(3), INTENT(IN) :: r
990 REAL(dp), DIMENSION(3, 2025), INTENT(OUT) :: dw
991 REAL(dp), INTENT(IN) :: delta
992 TYPE(se_int_control_type), INTENT(IN) :: se_int_control
993 TYPE(se_taper_type), POINTER :: se_taper
994
995 INTEGER :: i, j, nsize
996 REAL(dp) :: od
997 REAL(dp), DIMENSION(2025) :: wm, wp
998 REAL(dp), DIMENSION(3) :: rr
999
1000 od = 0.5_dp/delta
1001 nsize = sepi%atm_int_size*sepj%atm_int_size
1002 DO i = 1, 3
1003 rr = r
1004 rr(i) = rr(i) + delta
1005 CALL rotint_num(sepi, sepj, rr, wp, se_int_control, se_taper=se_taper)
1006 rr(i) = rr(i) - 2._dp*delta
1007 CALL rotint_num(sepi, sepj, rr, wm, se_int_control, se_taper=se_taper)
1008 DO j = 1, nsize
1009 dw(i, j) = od*(wp(j) - wm(j))
1010 END DO
1011 END DO
1012
1013 END SUBROUTINE drotint_num
1014
1015! **************************************************************************************************
1016!> \brief Numerical Derivatives for rotnuc
1017!> \param sepi ...
1018!> \param sepj ...
1019!> \param r ...
1020!> \param de1b ...
1021!> \param de2a ...
1022!> \param itype ...
1023!> \param delta ...
1024!> \param se_int_control ...
1025!> \param se_taper ...
1026! **************************************************************************************************
1027 SUBROUTINE drotnuc_num(sepi, sepj, r, de1b, de2a, itype, delta, se_int_control, se_taper)
1028 TYPE(semi_empirical_type), POINTER :: sepi, sepj
1029 REAL(dp), DIMENSION(3), INTENT(IN) :: r
1030 REAL(dp), DIMENSION(3, 45), INTENT(OUT), OPTIONAL :: de1b, de2a
1031 INTEGER, INTENT(IN) :: itype
1032 REAL(dp), INTENT(IN) :: delta
1033 TYPE(se_int_control_type), INTENT(IN) :: se_int_control
1034 TYPE(se_taper_type), POINTER :: se_taper
1035
1036 INTEGER :: i, j
1037 LOGICAL :: l_de1b, l_de2a
1038 REAL(dp) :: od
1039 REAL(dp), DIMENSION(3) :: rr
1040 REAL(dp), DIMENSION(45) :: e1m, e1p, e2m, e2p
1041
1042 l_de1b = PRESENT(de1b)
1043 l_de2a = PRESENT(de2a)
1044 IF (.NOT. (l_de1b .OR. l_de2a)) RETURN
1045 od = 0.5_dp/delta
1046 DO i = 1, 3
1047 rr = r
1048 rr(i) = rr(i) + delta
1049 CALL rotnuc_num(sepi, sepj, rr, e1p, e2p, itype, se_int_control, se_taper=se_taper)
1050 rr(i) = rr(i) - 2._dp*delta
1051 CALL rotnuc_num(sepi, sepj, rr, e1m, e2m, itype, se_int_control, se_taper=se_taper)
1052 IF (l_de1b) THEN
1053 DO j = 1, sepi%atm_int_size
1054 de1b(i, j) = od*(e1p(j) - e1m(j))
1055 END DO
1056 END IF
1057 IF (l_de2a) THEN
1058 DO j = 1, sepj%atm_int_size
1059 de2a(i, j) = od*(e2p(j) - e2m(j))
1060 END DO
1061 END IF
1062 END DO
1063 END SUBROUTINE drotnuc_num
1064
1065! **************************************************************************************************
1066!> \brief Numerical Derivatives for corecore
1067!> \param sepi ...
1068!> \param sepj ...
1069!> \param r ...
1070!> \param denuc ...
1071!> \param itype ...
1072!> \param delta ...
1073!> \param se_int_control ...
1074!> \param se_taper ...
1075! **************************************************************************************************
1076 SUBROUTINE dcorecore_num(sepi, sepj, r, denuc, itype, delta, se_int_control, se_taper)
1077 TYPE(semi_empirical_type), POINTER :: sepi, sepj
1078 REAL(dp), DIMENSION(3), INTENT(IN) :: r
1079 REAL(dp), DIMENSION(3), INTENT(OUT) :: denuc
1080 INTEGER, INTENT(IN) :: itype
1081 REAL(dp), INTENT(IN) :: delta
1082 TYPE(se_int_control_type), INTENT(IN) :: se_int_control
1083 TYPE(se_taper_type), POINTER :: se_taper
1084
1085 INTEGER :: i
1086 REAL(dp) :: enucm, enucp, od
1087 REAL(dp), DIMENSION(3) :: rr
1088
1089 od = 0.5_dp/delta
1090 DO i = 1, 3
1091 rr = r
1092 rr(i) = rr(i) + delta
1093 CALL corecore_num(sepi, sepj, rr, enucp, itype, se_int_control, se_taper=se_taper)
1094 rr(i) = rr(i) - 2._dp*delta
1095 CALL corecore_num(sepi, sepj, rr, enucm, itype, se_int_control, se_taper=se_taper)
1096 denuc(i) = od*(enucp - enucm)
1097 END DO
1098 END SUBROUTINE dcorecore_num
1099
1100! **************************************************************************************************
1101!> \brief Numerical Derivatives for corecore
1102!> \param sepi ...
1103!> \param sepj ...
1104!> \param r ...
1105!> \param denuc ...
1106!> \param itype ...
1107!> \param delta ...
1108!> \param se_int_control ...
1109!> \param se_taper ...
1110! **************************************************************************************************
1111 SUBROUTINE dcorecore_el_num(sepi, sepj, r, denuc, itype, delta, se_int_control, se_taper)
1112 TYPE(semi_empirical_type), POINTER :: sepi, sepj
1113 REAL(dp), DIMENSION(3), INTENT(IN) :: r
1114 REAL(dp), DIMENSION(3), INTENT(OUT) :: denuc
1115 INTEGER, INTENT(IN) :: itype
1116 REAL(dp), INTENT(IN) :: delta
1117 TYPE(se_int_control_type), INTENT(IN) :: se_int_control
1118 TYPE(se_taper_type), POINTER :: se_taper
1119
1120 INTEGER :: i
1121 REAL(dp) :: enucm, enucp, od
1122 REAL(dp), DIMENSION(3) :: rr
1123
1124 od = 0.5_dp/delta
1125 DO i = 1, 3
1126 rr = r
1127 rr(i) = rr(i) + delta
1128 CALL corecore_el_num(sepi, sepj, rr, enucp, itype, se_int_control, se_taper=se_taper)
1129 rr(i) = rr(i) - 2._dp*delta
1130 CALL corecore_el_num(sepi, sepj, rr, enucm, itype, se_int_control, se_taper=se_taper)
1131 denuc(i) = od*(enucp - enucm)
1132 END DO
1133 END SUBROUTINE dcorecore_el_num
1134
1135END MODULE semi_empirical_int_num
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_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::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::rotnuc_num
subroutine, public rotnuc_num(sepi, sepj, rijv, e1b, e2a, itype, se_int_control, se_taper)
Computes the two-particle interactions.
Definition semi_empirical_int_num.F:394
semi_empirical_int_num::corecore_el_num
subroutine, public corecore_el_num(sepi, sepj, rijv, enuc, itype, se_int_control, se_taper)
Computes the electrostatic core-core interactions only.
Definition semi_empirical_int_num.F:674
semi_empirical_int_num::rotint_num
subroutine, public rotint_num(sepi, sepj, rijv, w, se_int_control, se_taper)
Computes the two particle interactions in the lab frame.
Definition semi_empirical_int_num.F:77
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::drotint_num
subroutine, public drotint_num(sepi, sepj, r, dw, delta, se_int_control, se_taper)
Numerical Derivatives for rotint.
Definition semi_empirical_int_num.F:988
semi_empirical_int_num::core_nucint_num
subroutine, public core_nucint_num(sepi, sepj, rij, core, itype, se_taper, se_int_control)
Calculates the nuclear attraction integrals CORE (main driver)
Definition semi_empirical_int_num.F:774
semi_empirical_int_num::ssss_nucint_num
subroutine, public ssss_nucint_num(sepi, sepj, rij, ssss, itype, se_taper, se_int_control)
Calculates the SSSS integrals (main driver)
Definition semi_empirical_int_num.F:723
semi_empirical_int_num::drotnuc_num
subroutine, public drotnuc_num(sepi, sepj, r, de1b, de2a, itype, delta, se_int_control, se_taper)
Numerical Derivatives for rotnuc.
Definition semi_empirical_int_num.F:1028
semi_empirical_int_num::dcorecore_num
subroutine, public dcorecore_num(sepi, sepj, r, denuc, itype, delta, se_int_control, se_taper)
Numerical Derivatives for corecore.
Definition semi_empirical_int_num.F:1077
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::dcorecore_el_num
subroutine, public dcorecore_el_num(sepi, sepj, r, denuc, itype, delta, se_int_control, se_taper)
Numerical Derivatives for corecore.
Definition semi_empirical_int_num.F:1112
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_num::corecore_num
subroutine, public corecore_num(sepi, sepj, rijv, enuc, itype, se_int_control, se_taper)
Computes the core-core interactions.
Definition semi_empirical_int_num.F:500
semi_empirical_int_num::terep_num
subroutine, public terep_num(sepi, sepj, rij, rep, se_taper, se_int_control)
Calculates the derivative pf two-electron repulsion integrals and the nuclear attraction integrals w....
Definition semi_empirical_int_num.F:226
semi_empirical_int_utils
Utilities for Integrals for semi-empiric methods.
Definition semi_empirical_int_utils.F:12
semi_empirical_int_utils::ijkl_sp
real(kind=dp) function, public ijkl_sp(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_control, se_int_screen, itype)
General driver for computing semi-empirical integrals <ij|kl> with sp basis set. This code uses the o...
Definition semi_empirical_int_utils.F:310
semi_empirical_int_utils::store_2el_2c_diag
subroutine, public store_2el_2c_diag(limij, limkl, ww, w, ww_dx, ww_dy, ww_dz, dw)
Store the two-electron two-center integrals in diagonal form.
Definition semi_empirical_int_utils.F:2629
semi_empirical_int_utils::rotmat
recursive subroutine, public rotmat(sepi, sepj, rjiv, r, ij_matrix, do_derivatives, do_invert, debug_invert)
Computes the general rotation matrices for the integrals between I and J (J>=I)
Definition semi_empirical_int_utils.F:1883
semi_empirical_int_utils::rot_2el_2c_first
recursive subroutine, public rot_2el_2c_first(sepi, sepj, rijv, se_int_control, se_taper, invert, ii, kk, rep, logv, ij_matrix, v, lgrad, rep_d, v_d, logv_d, drij)
First Step of the rotation of the two-electron two-center integrals in SPD basis.
Definition semi_empirical_int_utils.F:2295
semi_empirical_int_utils::ijkl_d
real(kind=dp) function, public ijkl_d(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_control, se_int_screen, itype)
General driver for computing semi-empirical integrals <ij|kl> involving d-orbitals....
Definition semi_empirical_int_utils.F:1199
semi_empirical_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::taper_eval
real(kind=dp) function, public taper_eval(taper, rij)
Taper functions.
Definition taper_types.F:79
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