(git:34ef472)
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 ! **************************************************************************************************
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  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 
1135 END 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