LCOV - code coverage report
Current view: top level - src - semi_empirical_int_num.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 85.0 % 400 340
Test Date: 2025-07-25 12:55:17 Functions: 86.7 % 15 13

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 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     42876054 :    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     42876054 :       l_core = PRESENT(core)
     852     42876054 :       l_ssss = PRESENT(ssss)
     853     42876054 :       IF (.NOT. (l_core .OR. l_ssss)) RETURN
     854     42876054 :       si = (sepi%natorb > 1)
     855     42876054 :       sj = (sepj%natorb > 1)
     856              : 
     857     42876054 :       ij = indexa(1, 1)
     858     42876054 :       IF (l_ssss) THEN
     859              :          ! Store the value for <S  S  | S  S  > (Used for computing the core-core interactions)
     860     26706798 :          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     42876054 :       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
        

Generated by: LCOV version 2.0-1