Line data Source code
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 : ! **************************************************************************************************
17 : MODULE 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 :
61 : CONTAINS
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 51626 : 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 51626 : NULLIFY (ij_matrix)
95 206504 : rij = DOT_PRODUCT(rijv, rijv)
96 51626 : 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 51626 : rij = SQRT(rij)
102 :
103 : ! Create the rotation matrix
104 51626 : CALL rotmat_create(ij_matrix)
105 51626 : CALL rotmat(sepi, sepj, rijv, rij, ij_matrix, do_derivatives=.FALSE.)
106 :
107 : ! Compute Integrals in Diatomic Frame
108 51626 : CALL terep_num(sepi, sepj, rij, rep, se_taper=se_taper, se_int_control=se_int_control)
109 :
110 : ! Rotate Integrals
111 51626 : ii = sepi%natorb
112 51626 : kk = sepj%natorb
113 51626 : IF (ii*kk > 0) THEN
114 51626 : limij = sepi%atm_int_size
115 51626 : limkl = sepj%atm_int_size
116 51626 : istep = limkl*limij
117 1644994 : DO i1 = 1, istep
118 1644994 : 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 51626 : logv, ij_matrix, v, lgrad=.FALSE.)
124 :
125 : ! Second step in rotation of integrals
126 209880 : DO i1 = 1, ii
127 582790 : DO j1 = 1, i1
128 372910 : ij = indexa(i1, j1)
129 372910 : jj = indexb(i1, j1)
130 372910 : mm = int2c_type(jj)
131 1306600 : DO k = 1, kk
132 2741714 : DO l = 1, k
133 1593368 : kl = indexb(k, l)
134 2368804 : IF (logv(ij, kl)) THEN
135 1270248 : wrepp = v(ij, kl)
136 202078 : SELECT CASE (mm)
137 : CASE (1)
138 : ! (SS/)
139 202078 : i = 1
140 202078 : j = 1
141 202078 : iw_loc = (indexb(i, j) - 1)*limkl + kl
142 202078 : ww(iw_loc) = wrepp
143 : CASE (2)
144 : ! (SP/)
145 1397200 : j = 1
146 1397200 : DO i = 1, 3
147 1047900 : iw_loc = (indexb(i + 1, j) - 1)*limkl + kl
148 1397200 : ww(iw_loc) = ww(iw_loc) + ij_matrix%sp(i1 - 1, i)*wrepp
149 : END DO
150 : CASE (3)
151 : ! (PP/)
152 2837896 : DO i = 1, 3
153 2128422 : cc = ij_matrix%pp(i, i1 - 1, j1 - 1)
154 2128422 : iw_loc = (indexb(i + 1, i + 1) - 1)*limkl + kl
155 2128422 : ww(iw_loc) = ww(iw_loc) + cc*wrepp
156 2128422 : iminus = i - 1
157 2837896 : IF (iminus /= 0) THEN
158 3547370 : DO j = 1, iminus
159 2128422 : cc = ij_matrix%pp(1 + i + j, i1 - 1, j1 - 1)
160 2128422 : iw_loc = (indexb(i + 1, j + 1) - 1)*limkl + kl
161 3547370 : ww(iw_loc) = ww(iw_loc) + cc*wrepp
162 : END DO
163 : END IF
164 : END DO
165 : CASE (4)
166 : ! (SD/)
167 8136 : j = 1
168 8136 : DO i = 1, 5
169 6780 : iw_loc = (indexb(i + 4, j) - 1)*limkl + kl
170 8136 : ww(iw_loc) = ww(iw_loc) + ij_matrix%sd(i1 - 4, i)*wrepp
171 : END DO
172 : CASE (5)
173 : ! (DP/)
174 21528 : DO i = 1, 5
175 75348 : DO j = 1, 3
176 53820 : iw_loc = (indexb(i + 4, j + 1) - 1)*limkl + kl
177 53820 : ij1 = 3*(i - 1) + j
178 71760 : 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 1292508 : DO i = 1, 5
184 22260 : cc = ij_matrix%dd(i, i1 - 4, j1 - 4)
185 22260 : iw_loc = (indexb(i + 4, i + 4) - 1)*limkl + kl
186 22260 : ww(iw_loc) = ww(iw_loc) + cc*wrepp
187 22260 : iminus = i - 1
188 26712 : IF (iminus /= 0) THEN
189 62328 : DO j = 1, iminus
190 44520 : ij1 = inddd(i, j)
191 44520 : cc = ij_matrix%dd(ij1, i1 - 4, j1 - 4)
192 44520 : iw_loc = (indexb(i + 4, j + 4) - 1)*limkl + kl
193 62328 : 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 51626 : CALL rotmat_release(ij_matrix)
205 :
206 : ! Store two electron integrals in the triangular format
207 51626 : CALL store_2el_2c_diag(limij, limkl, ww, w)
208 : END IF
209 51626 : 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 51626 : 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 51626 : ft = taper_eval(se_taper%taper, rij)
236 : ! In case of dumped integrals compute an additional taper term
237 51626 : IF (se_int_control%integral_screening == do_se_IS_kdso_d) THEN
238 0 : se_int_screen%ft = taper_eval(se_taper%taper_add, rij)
239 : END IF
240 :
241 : ! Contribution from sp shells
242 51626 : CALL terep_sp_num(sepi, sepj, rij, rep, se_int_control, se_int_screen, ft)
243 :
244 51626 : 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 84 : ft)
248 : END IF
249 :
250 51626 : 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 1195735 : 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 1195735 : lasti = sepi%natorb
281 1195735 : lastj = sepj%natorb
282 4519013 : DO i = 1, MIN(lasti, 4)
283 3323278 : li = l_index(i)
284 12097377 : DO j = 1, i
285 7578364 : lj = l_index(j)
286 7578364 : ij = indexa(i, j)
287 30134196 : DO k = 1, MIN(lastj, 4)
288 19232554 : lk = l_index(k)
289 69351852 : DO l = 1, k
290 42540934 : ll = l_index(l)
291 42540934 : kl = indexa(k, l)
292 :
293 42540934 : numb = ijkl_ind(ij, kl)
294 61773488 : IF (numb > 0) THEN
295 15486485 : nold = ijkl_sym(numb)
296 15486485 : IF (nold > 0) THEN
297 4965265 : rep(numb) = rep(nold)
298 10521220 : ELSE IF (nold < 0) THEN
299 0 : rep(numb) = -rep(-nold)
300 10521220 : 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 10521220 : *ft
304 10521220 : rep(numb) = tmp
305 : END IF
306 : END IF
307 : END DO
308 : END DO
309 : END DO
310 : END DO
311 1195735 : 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 54584 : 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 54584 : lasti = sepi%natorb
342 54584 : lastj = sepj%natorb
343 423134 : DO i = 1, lasti
344 368550 : li = l_index(i)
345 2082056 : DO j = 1, i
346 1658922 : lj = l_index(j)
347 1658922 : ij = indexa(i, j)
348 9638700 : DO k = 1, lastj
349 7611228 : lk = l_index(k)
350 37446630 : DO l = 1, k
351 28176480 : ll = l_index(l)
352 28176480 : kl = indexa(k, l)
353 :
354 28176480 : numb = ijkl_ind(ij, kl)
355 : ! From 1 to 34 we store integrals involving sp shells
356 35787708 : IF (numb > 34) THEN
357 5906296 : nold = ijkl_sym(numb)
358 5906296 : IF (nold > 34) THEN
359 2755232 : rep(numb) = rep(nold)
360 3151064 : ELSE IF (nold < -34) THEN
361 537320 : rep(numb) = -rep(-nold)
362 2613744 : 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 2613744 : *ft
366 2613744 : rep(numb) = tmp
367 : END IF
368 : END IF
369 : END DO
370 : END DO
371 : END DO
372 : END DO
373 :
374 54584 : 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 26956 : 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 26956 : NULLIFY (ij_matrix)
410 26956 : l_e1b = PRESENT(e1b)
411 26956 : l_e2a = PRESENT(e2a)
412 107824 : rij = DOT_PRODUCT(rijv, rijv)
413 :
414 26956 : IF (rij > rij_threshold) THEN
415 26956 : rij = SQRT(rij)
416 : ! Create the rotation matrix
417 26956 : CALL rotmat_create(ij_matrix)
418 26956 : 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 26956 : se_int_control=se_int_control)
423 :
424 : ! Copy parameters over to arrays for do loop.
425 26956 : last_orbital(1) = sepi%natorb
426 26956 : last_orbital(2) = sepj%natorb
427 26956 : task(1) = l_e1b
428 26956 : task(2) = l_e2a
429 80868 : DO n = 1, 2
430 53912 : IF (.NOT. task(n)) CYCLE
431 186279 : DO i = 1, last_orbital(n)
432 132538 : ind1 = i - 1
433 479335 : DO j = 1, i
434 293056 : ind2 = j - 1
435 293056 : m = (i*(i - 1))/2 + j
436 : ! Perform Rotations ...
437 425594 : IF (ind2 == 0) THEN
438 132538 : IF (ind1 == 0) THEN
439 : ! Type of Integral (SS/)
440 52769 : tmp(m) = core(1, n)
441 79769 : ELSE IF (ind1 < 4) THEN
442 : ! Type of Integral (SP/)
443 79524 : tmp(m) = ij_matrix%sp(1, ind1)*core(2, n)
444 : ELSE
445 : ! Type of Integral (SD/)
446 245 : tmp(m) = ij_matrix%sd(1, ind1 - 3)*core(5, n)
447 : END IF
448 160518 : ELSE IF (ind2 < 4) THEN
449 159783 : IF (ind1 < 4) THEN
450 : ! Type of Integral (PP/)
451 159048 : ipp = indpp(ind1, ind2)
452 : tmp(m) = core(3, n)*ij_matrix%pp(ipp, 1, 1) + &
453 159048 : core(4, n)*(ij_matrix%pp(ipp, 2, 2) + ij_matrix%pp(ipp, 3, 3))
454 : ELSE
455 : ! Type of Integral (PD/)
456 735 : idp = inddp(ind1 - 3, ind2)
457 : tmp(m) = core(6, n)*ij_matrix%pd(idp, 1, 1) + &
458 735 : 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 735 : 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 735 : 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 53741 : IF (n == 1) THEN
470 174507 : DO i = 1, sepi%atm_int_size
471 174507 : e1b(i) = -tmp(i)
472 : END DO
473 : END IF
474 80697 : IF (n == 2) THEN
475 172290 : DO i = 1, sepj%atm_int_size
476 172290 : e2a(i) = -tmp(i)
477 : END DO
478 : END IF
479 : END DO
480 26956 : CALL rotmat_release(ij_matrix)
481 : END IF
482 26956 : 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 29893 : 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 119572 : rij = DOT_PRODUCT(rijv, rijv)
516 29893 : enuc = 0.0_dp
517 29893 : IF (rij > rij_threshold) THEN
518 29893 : 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 29893 : 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 29893 : se_int_control=se_int_control_off)
524 : ! In case let's compute the short-range part of the (ss|ss) integral
525 29893 : IF (se_int_control%shortrange) THEN
526 : CALL ssss_nucint_num(sepi, sepj, rij, ssss=ssss_sr, itype=itype, se_taper=se_taper, &
527 0 : se_int_control=se_int_control)
528 : ELSE
529 29893 : ssss_sr = ssss
530 : END IF
531 29893 : zz = sepi%zeff*sepj%zeff
532 : ! Nuclear-Nuclear electrostatic contribution
533 29893 : enuc = zz*ssss_sr
534 : ! Method dependent correction terms
535 29893 : IF (itype /= do_method_pm6 .AND. itype /= do_method_pm6fm) THEN
536 28136 : alpi = sepi%alp
537 28136 : alpj = sepj%alp
538 28136 : scale = EXP(-alpi*rij) + EXP(-alpj*rij)
539 :
540 28136 : nt = sepi%z + sepj%z
541 28136 : IF (nt == 8 .OR. nt == 9) THEN
542 11319 : IF (sepi%z == 7 .OR. sepi%z == 8) scale = scale + (angstrom*rij - 1._dp)*EXP(-alpi*rij)
543 11319 : IF (sepj%z == 7 .OR. sepj%z == 8) scale = scale + (angstrom*rij - 1._dp)*EXP(-alpj*rij)
544 : END IF
545 28136 : scale = ABS(scale*zz*ssss)
546 28136 : zz = zz/rij
547 28136 : IF (itype == do_method_am1 .OR. itype == do_method_pm3 .OR. itype == do_method_pdg) THEN
548 16710 : IF (itype == do_method_am1 .AND. sepi%z == 5) THEN
549 : !special case AM1 Boron
550 : SELECT CASE (sepj%z)
551 : CASE DEFAULT
552 0 : nt = 1
553 : CASE (1)
554 0 : nt = 2
555 : CASE (6)
556 0 : nt = 3
557 : CASE (9, 17, 35, 53)
558 0 : nt = 4
559 : END SELECT
560 0 : fni1(:) = sepi%bfn1(:, nt)
561 0 : fni2(:) = sepi%bfn2(:, nt)
562 0 : fni3(:) = sepi%bfn3(:, nt)
563 : ELSE
564 83550 : fni1(:) = sepi%fn1(:)
565 83550 : fni2(:) = sepi%fn2(:)
566 83550 : fni3(:) = sepi%fn3(:)
567 : END IF
568 16710 : IF (itype == do_method_am1 .AND. sepj%z == 5) THEN
569 : !special case AM1 Boron
570 : SELECT CASE (sepi%z)
571 : CASE DEFAULT
572 0 : nt = 1
573 : CASE (1)
574 0 : nt = 2
575 : CASE (6)
576 0 : nt = 3
577 : CASE (9, 17, 35, 53)
578 0 : nt = 4
579 : END SELECT
580 0 : fnj1(:) = sepj%bfn1(:, nt)
581 0 : fnj2(:) = sepj%bfn2(:, nt)
582 0 : fnj3(:) = sepj%bfn3(:, nt)
583 : ELSE
584 83550 : fnj1(:) = sepj%fn1(:)
585 83550 : fnj2(:) = sepj%fn2(:)
586 83550 : fnj3(:) = sepj%fn3(:)
587 : END IF
588 : ! AM1/PM3/PDG correction to nuclear repulsion
589 83550 : DO ig = 1, SIZE(fni1)
590 66840 : IF (ABS(fni1(ig)) > 0._dp) THEN
591 34372 : ax = fni2(ig)*(rij - fni3(ig))**2
592 34372 : IF (ax <= 25._dp) THEN
593 19892 : scale = scale + zz*fni1(ig)*EXP(-ax)
594 : END IF
595 : END IF
596 83550 : IF (ABS(fnj1(ig)) > 0._dp) THEN
597 34281 : ax = fnj2(ig)*(rij - fnj3(ig))**2
598 34281 : IF (ax <= 25._dp) THEN
599 19985 : scale = scale + zz*fnj1(ig)*EXP(-ax)
600 : END IF
601 : END IF
602 : END DO
603 : END IF
604 28136 : IF (itype == do_method_pdg) THEN
605 : ! PDDG function
606 0 : zaf = sepi%zeff/nt
607 0 : zbf = sepj%zeff/nt
608 0 : pai = sepi%pre(1)
609 0 : pbi = sepi%pre(2)
610 0 : paj = sepj%pre(1)
611 0 : pbj = sepj%pre(2)
612 0 : dai = sepi%d(1)
613 0 : dbi = sepi%d(2)
614 0 : daj = sepj%d(1)
615 0 : dbj = sepj%d(2)
616 0 : 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 0 : (zaf*pbi + zbf*pbj)*EXP(-apdg*(rij - dbi - dbj)**2)
622 28136 : ELSEIF (itype == do_method_pchg) THEN
623 : qcorr = 0.0_dp
624 : scale = 0.0_dp
625 : ELSE
626 26993 : qcorr = 0.0_dp
627 : END IF
628 : ELSE
629 1757 : rija = rij*angstrom
630 1757 : scale = zz*ssss
631 : ! PM6 core-core terms
632 1757 : xab = sepi%xab(sepj%z)
633 1757 : aab = sepi%aab(sepj%z)
634 1757 : 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 246 : scale = scale*(2._dp*xab*EXP(-aab*rija*rija))
638 1511 : ELSEIF (sepi%z == 6 .AND. sepj%z == 6) THEN
639 : ! Special Case C-C
640 0 : scale = scale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6)) + 9.28_dp*EXP(-5.98_dp*rija))
641 1511 : ELSEIF ((sepi%z == 8 .AND. sepj%z == 14) .OR. &
642 : (sepj%z == 8 .AND. sepi%z == 14)) THEN
643 : ! Special Case Si-O
644 0 : 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 1511 : 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 1757 : (sepj%a*EXP(-sepj%b*(rij - sepj%c)**2))*sepi%zeff*sepj%zeff/rij
653 : ! Hard core repulsion
654 1757 : tmp = (REAL(sepi%z, dp)**(1._dp/3._dp) + REAL(sepj%z, dp)**(1._dp/3._dp))
655 1757 : qcorr = qcorr + 1.e-8_dp/evolt*(tmp/rija)**12
656 : END IF
657 29893 : enuc = enuc + scale + qcorr
658 : END IF
659 29893 : 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 0 : 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 0 : rij = DOT_PRODUCT(rijv, rijv)
685 0 : enuc = 0.0_dp
686 0 : IF (rij > rij_threshold) THEN
687 0 : 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 0 : 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 0 : se_int_control=se_int_control_off)
693 : ! In case let's compute the short-range part of the (ss|ss) integral
694 0 : 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 0 : se_int_control=se_int_control)
697 : ELSE
698 0 : ssss_sr = ssss
699 : END IF
700 0 : zz = sepi%zeff*sepj%zeff
701 : ! Nuclear-Nuclear electrostatic contribution
702 0 : enuc = zz*ssss_sr
703 : END IF
704 0 : 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 29893 : 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 29893 : ft = 1.0_dp
736 29893 : IF (itype /= do_method_pchg) THEN
737 28750 : ft = taper_eval(se_taper%taper, rij)
738 : END IF
739 :
740 : ! In case of dumped integrals compute an additional taper term
741 29893 : IF (se_int_control%integral_screening == do_se_IS_kdso_d) THEN
742 0 : se_int_screen%ft = 1.0_dp
743 0 : IF (itype /= do_method_pchg) THEN
744 0 : 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 29893 : se_int_control=se_int_control, se_int_screen=se_int_screen)
751 :
752 : ! Tapering the integrals
753 29893 : ssss = ft*ssss
754 :
755 29893 : 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 26956 : 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 26956 : ft = 1.0_dp
788 26956 : IF (itype /= do_method_pchg) THEN
789 25813 : ft = taper_eval(se_taper%taper, rij)
790 : END IF
791 :
792 : ! In case of dumped integrals compute an additional taper term
793 26956 : IF (se_int_control%integral_screening == do_se_IS_kdso_d) THEN
794 0 : se_int_screen%ft = 1.0_dp
795 0 : IF (itype /= do_method_pchg) THEN
796 0 : 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 26956 : se_int_control=se_int_control, se_int_screen=se_int_screen)
803 :
804 26956 : 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 42 : se_int_screen)
808 : END IF
809 :
810 : ! Tapering the integrals
811 94031 : DO i = 1, sepi%core_size
812 94031 : core(i, 1) = ft*core(i, 1)
813 : END DO
814 92468 : DO i = 1, sepj%core_size
815 92468 : core(i, 2) = ft*core(i, 2)
816 : END DO
817 :
818 26956 : 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 42675234 : 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 42675234 : l_core = PRESENT(core)
852 42675234 : l_ssss = PRESENT(ssss)
853 42675234 : IF (.NOT. (l_core .OR. l_ssss)) RETURN
854 42675234 : si = (sepi%natorb > 1)
855 42675234 : sj = (sepj%natorb > 1)
856 :
857 42675234 : ij = indexa(1, 1)
858 42675234 : IF (l_ssss) THEN
859 : ! Store the value for <S S | S S > (Used for computing the core-core interactions)
860 26505978 : 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 42675234 : IF (l_core) THEN
864 : ! <S S | S S >
865 16169256 : kl = indexa(1, 1)
866 16169256 : 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 16169256 : IF (si) THEN
868 : ! <S P | S S >
869 7252446 : kl = indexa(2, 1)
870 7252446 : 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 7252446 : kl = indexa(2, 2)
873 7252446 : 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 7252446 : kl = indexa(3, 3)
876 7252446 : 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 16169256 : kl = indexa(1, 1)
881 16169256 : 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 16169256 : IF (sj) THEN
883 : ! <S S | S P >
884 422791 : kl = indexa(2, 1)
885 422791 : 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 422791 : kl = indexa(2, 2)
888 422791 : 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 422791 : kl = indexa(3, 3)
891 422791 : 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 37964 : 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 37964 : IF (sepi%dorb .OR. sepj%dorb) THEN
932 37964 : ij = indexa(1, 1)
933 37964 : IF (sepj%dorb) THEN
934 : ! <S S | D S>
935 21357 : kl = indexa(5, 1)
936 21357 : 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 21357 : kl = indexa(5, 2)
939 21357 : 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 21357 : kl = indexa(5, 5)
942 21357 : 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 21357 : kl = indexa(6, 3)
945 21357 : 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 21357 : kl = indexa(6, 6)
948 21357 : 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 21357 : kl = indexa(8, 8)
951 21357 : 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 37964 : IF (sepi%dorb) THEN
954 : ! <D S | S S>
955 21966 : kl = indexa(5, 1)
956 21966 : 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 21966 : kl = indexa(5, 2)
959 21966 : 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 21966 : kl = indexa(5, 5)
962 21966 : 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 21966 : kl = indexa(6, 3)
965 21966 : 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 21966 : kl = indexa(6, 6)
968 21966 : 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 21966 : kl = indexa(8, 8)
971 21966 : 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 37964 : 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 7318 : 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 7318 : od = 0.5_dp/delta
1001 7318 : nsize = sepi%atm_int_size*sepj%atm_int_size
1002 29272 : DO i = 1, 3
1003 21954 : rr = r
1004 21954 : rr(i) = rr(i) + delta
1005 21954 : CALL rotint_num(sepi, sepj, rr, wp, se_int_control, se_taper=se_taper)
1006 21954 : rr(i) = rr(i) - 2._dp*delta
1007 21954 : CALL rotint_num(sepi, sepj, rr, wm, se_int_control, se_taper=se_taper)
1008 706372 : DO j = 1, nsize
1009 699054 : dw(i, j) = od*(wp(j) - wm(j))
1010 : END DO
1011 : END DO
1012 :
1013 7318 : 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 3821 : 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 3821 : l_de1b = PRESENT(de1b)
1043 3821 : l_de2a = PRESENT(de2a)
1044 3821 : IF (.NOT. (l_de1b .OR. l_de2a)) RETURN
1045 3821 : od = 0.5_dp/delta
1046 15284 : DO i = 1, 3
1047 11463 : rr = r
1048 11463 : rr(i) = rr(i) + delta
1049 11463 : CALL rotnuc_num(sepi, sepj, rr, e1p, e2p, itype, se_int_control, se_taper=se_taper)
1050 11463 : rr(i) = rr(i) - 2._dp*delta
1051 11463 : CALL rotnuc_num(sepi, sepj, rr, e1m, e2m, itype, se_int_control, se_taper=se_taper)
1052 11463 : IF (l_de1b) THEN
1053 74193 : DO j = 1, sepi%atm_int_size
1054 74193 : de1b(i, j) = od*(e1p(j) - e1m(j))
1055 : END DO
1056 : END IF
1057 15284 : IF (l_de2a) THEN
1058 72834 : DO j = 1, sepj%atm_int_size
1059 72834 : 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 3732 : 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 3732 : od = 0.5_dp/delta
1090 14928 : DO i = 1, 3
1091 11196 : rr = r
1092 11196 : rr(i) = rr(i) + delta
1093 11196 : CALL corecore_num(sepi, sepj, rr, enucp, itype, se_int_control, se_taper=se_taper)
1094 11196 : rr(i) = rr(i) - 2._dp*delta
1095 11196 : CALL corecore_num(sepi, sepj, rr, enucm, itype, se_int_control, se_taper=se_taper)
1096 14928 : denuc(i) = od*(enucp - enucm)
1097 : END DO
1098 3732 : 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 0 : 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 0 : od = 0.5_dp/delta
1125 0 : DO i = 1, 3
1126 0 : rr = r
1127 0 : rr(i) = rr(i) + delta
1128 0 : CALL corecore_el_num(sepi, sepj, rr, enucp, itype, se_int_control, se_taper=se_taper)
1129 0 : rr(i) = rr(i) - 2._dp*delta
1130 0 : CALL corecore_el_num(sepi, sepj, rr, enucm, itype, se_int_control, se_taper=se_taper)
1131 0 : denuc(i) = od*(enucp - enucm)
1132 : END DO
1133 0 : END SUBROUTINE dcorecore_el_num
1134 :
1135 : END MODULE semi_empirical_int_num
|
