LCOV - code coverage report
Current view: top level - src - semi_empirical_int_num.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:58e3e09) Lines: 340 400 85.0 %
Date: 2024-03-29 07:50:05 Functions: 13 15 86.7 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Integrals for semi-empiric methods
      10             : !> \par History
      11             : !>      JGH (27.10.2004) : separate routine for nuclear attraction integrals
      12             : !>      JGH (13.03.2006) : tapering function
      13             : !>      Teodoro Laino (03.2008) [tlaino] - University of Zurich : new driver
      14             : !>                 for computing integrals
      15             : !> \author JGH (11.10.2004)
      16             : ! **************************************************************************************************
      17             : MODULE semi_empirical_int_num
      18             : 
      19             :    USE input_constants,                 ONLY: do_method_am1,&
      20             :                                               do_method_pchg,&
      21             :                                               do_method_pdg,&
      22             :                                               do_method_pm3,&
      23             :                                               do_method_pm6,&
      24             :                                               do_method_pm6fm,&
      25             :                                               do_method_undef,&
      26             :                                               do_se_IS_kdso_d
      27             :    USE kinds,                           ONLY: dp
      28             :    USE multipole_types,                 ONLY: do_multipole_none
      29             :    USE physcon,                         ONLY: angstrom,&
      30             :                                               evolt
      31             :    USE semi_empirical_int_arrays,       ONLY: &
      32             :         ijkl_ind, ijkl_sym, inddd, inddp, indexa, indexb, indpp, int2c_type, l_index, rij_threshold
      33             :    USE semi_empirical_int_utils,        ONLY: ijkl_d,&
      34             :                                               ijkl_sp,&
      35             :                                               rot_2el_2c_first,&
      36             :                                               rotmat,&
      37             :                                               store_2el_2c_diag
      38             :    USE semi_empirical_types,            ONLY: rotmat_create,&
      39             :                                               rotmat_release,&
      40             :                                               rotmat_type,&
      41             :                                               se_int_control_type,&
      42             :                                               se_int_screen_type,&
      43             :                                               se_taper_type,&
      44             :                                               semi_empirical_type,&
      45             :                                               setup_se_int_control_type
      46             :    USE taper_types,                     ONLY: taper_eval
      47             : #include "./base/base_uses.f90"
      48             : 
      49             :    IMPLICIT NONE
      50             : 
      51             :    PRIVATE
      52             : 
      53             :    LOGICAL, PARAMETER, PRIVATE          :: debug_this_module = .FALSE.
      54             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_int_num'
      55             :    PUBLIC :: rotint_num, rotnuc_num, corecore_num, &
      56             :              drotint_num, drotnuc_num, dcorecore_num, &
      57             :              ssss_nucint_num, core_nucint_num, terep_num, &
      58             :              nucint_sp_num, terep_sp_num, terep_d_num, &
      59             :              nucint_d_num, corecore_el_num, dcorecore_el_num
      60             : 
      61             : CONTAINS
      62             : 
      63             : ! **************************************************************************************************
      64             : !> \brief Computes the two particle interactions in the lab frame
      65             : !> \param sepi Atomic parameters of first atom
      66             : !> \param sepj Atomic parameters of second atom
      67             : !> \param rijv Coordinate vector i -> j
      68             : !> \param w Array of two-electron repulsion integrals.
      69             : !> \param se_int_control input parameters that control the calculation of SE
      70             : !>                           integrals (shortrange, R3 residual, screening type)
      71             : !> \param se_taper ...
      72             : !> \note routine adapted from mopac7 (rotate)
      73             : !>       written by Ernest R. Davidson, Indiana University.
      74             : !>       Teodoro Laino [tlaino] - University of Zurich 04.2008 : major rewriting
      75             : ! **************************************************************************************************
      76       51626 :    SUBROUTINE rotint_num(sepi, sepj, rijv, w, se_int_control, se_taper)
      77             :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
      78             :       REAL(dp), DIMENSION(3), INTENT(IN)                 :: rijv
      79             :       REAL(dp), DIMENSION(2025), INTENT(OUT)             :: w
      80             :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
      81             :       TYPE(se_taper_type), POINTER                       :: se_taper
      82             : 
      83             :       INTEGER                                            :: i, i1, ii, ij, ij1, iminus, istep, &
      84             :                                                             iw_loc, j, j1, jj, k, kk, kl, l, &
      85             :                                                             limij, limkl, mm
      86             :       LOGICAL, DIMENSION(45, 45)                         :: logv
      87             :       REAL(dp)                                           :: rij
      88             :       REAL(KIND=dp)                                      :: cc, wrepp
      89             :       REAL(KIND=dp), DIMENSION(2025)                     :: ww
      90             :       REAL(KIND=dp), DIMENSION(45, 45)                   :: v
      91             :       REAL(KIND=dp), DIMENSION(491)                      :: rep
      92             :       TYPE(rotmat_type), POINTER                         :: ij_matrix
      93             : 
      94       51626 :       NULLIFY (ij_matrix)
      95      206504 :       rij = DOT_PRODUCT(rijv, rijv)
      96       51626 :       IF (rij > rij_threshold) THEN
      97             :          ! The repulsion integrals over molecular frame (w) are stored in the
      98             :          ! order in which they will later be used.  ie.  (i,j/k,l) where
      99             :          ! j.le.i  and  l.le.k     and l varies most rapidly and i least
     100             :          ! rapidly.  (anti-normal computer storage)
     101       51626 :          rij = SQRT(rij)
     102             : 
     103             :          ! Create the rotation matrix
     104       51626 :          CALL rotmat_create(ij_matrix)
     105       51626 :          CALL rotmat(sepi, sepj, rijv, rij, ij_matrix, do_derivatives=.FALSE.)
     106             : 
     107             :          ! Compute Integrals in Diatomic Frame
     108       51626 :          CALL terep_num(sepi, sepj, rij, rep, se_taper=se_taper, se_int_control=se_int_control)
     109             : 
     110             :          ! Rotate Integrals
     111       51626 :          ii = sepi%natorb
     112       51626 :          kk = sepj%natorb
     113       51626 :          IF (ii*kk > 0) THEN
     114       51626 :             limij = sepi%atm_int_size
     115       51626 :             limkl = sepj%atm_int_size
     116       51626 :             istep = limkl*limij
     117     1644994 :             DO i1 = 1, istep
     118     1644994 :                ww(i1) = 0.0_dp
     119             :             END DO
     120             : 
     121             :             ! First step in rotation of integrals
     122             :             CALL rot_2el_2c_first(sepi, sepj, rijv, se_int_control, se_taper, .FALSE., ii, kk, rep, &
     123       51626 :                                   logv, ij_matrix, v, lgrad=.FALSE.)
     124             : 
     125             :             ! Second step in rotation of integrals
     126      209880 :             DO i1 = 1, ii
     127      582790 :                DO j1 = 1, i1
     128      372910 :                   ij = indexa(i1, j1)
     129      372910 :                   jj = indexb(i1, j1)
     130      372910 :                   mm = int2c_type(jj)
     131     1306600 :                   DO k = 1, kk
     132     2741714 :                      DO l = 1, k
     133     1593368 :                         kl = indexb(k, l)
     134     2368804 :                         IF (logv(ij, kl)) THEN
     135     1270248 :                            wrepp = v(ij, kl)
     136      202078 :                            SELECT CASE (mm)
     137             :                            CASE (1)
     138             :                               ! (SS/)
     139      202078 :                               i = 1
     140      202078 :                               j = 1
     141      202078 :                               iw_loc = (indexb(i, j) - 1)*limkl + kl
     142      202078 :                               ww(iw_loc) = wrepp
     143             :                            CASE (2)
     144             :                               ! (SP/)
     145     1397200 :                               j = 1
     146     1397200 :                               DO i = 1, 3
     147     1047900 :                                  iw_loc = (indexb(i + 1, j) - 1)*limkl + kl
     148     1397200 :                                  ww(iw_loc) = ww(iw_loc) + ij_matrix%sp(i1 - 1, i)*wrepp
     149             :                               END DO
     150             :                            CASE (3)
     151             :                               ! (PP/)
     152     2837896 :                               DO i = 1, 3
     153     2128422 :                                  cc = ij_matrix%pp(i, i1 - 1, j1 - 1)
     154     2128422 :                                  iw_loc = (indexb(i + 1, i + 1) - 1)*limkl + kl
     155     2128422 :                                  ww(iw_loc) = ww(iw_loc) + cc*wrepp
     156     2128422 :                                  iminus = i - 1
     157     2837896 :                                  IF (iminus /= 0) THEN
     158     3547370 :                                     DO j = 1, iminus
     159     2128422 :                                        cc = ij_matrix%pp(1 + i + j, i1 - 1, j1 - 1)
     160     2128422 :                                        iw_loc = (indexb(i + 1, j + 1) - 1)*limkl + kl
     161     3547370 :                                        ww(iw_loc) = ww(iw_loc) + cc*wrepp
     162             :                                     END DO
     163             :                                  END IF
     164             :                               END DO
     165             :                            CASE (4)
     166             :                               ! (SD/)
     167        8136 :                               j = 1
     168        8136 :                               DO i = 1, 5
     169        6780 :                                  iw_loc = (indexb(i + 4, j) - 1)*limkl + kl
     170        8136 :                                  ww(iw_loc) = ww(iw_loc) + ij_matrix%sd(i1 - 4, i)*wrepp
     171             :                               END DO
     172             :                            CASE (5)
     173             :                               ! (DP/)
     174       21528 :                               DO i = 1, 5
     175       75348 :                                  DO j = 1, 3
     176       53820 :                                     iw_loc = (indexb(i + 4, j + 1) - 1)*limkl + kl
     177       53820 :                                     ij1 = 3*(i - 1) + j
     178       71760 :                                     ww(iw_loc) = ww(iw_loc) + ij_matrix%pd(ij1, i1 - 4, j1 - 1)*wrepp
     179             :                                  END DO
     180             :                               END DO
     181             :                            CASE (6)
     182             :                               ! (DD/)
     183     1292508 :                               DO i = 1, 5
     184       22260 :                                  cc = ij_matrix%dd(i, i1 - 4, j1 - 4)
     185       22260 :                                  iw_loc = (indexb(i + 4, i + 4) - 1)*limkl + kl
     186       22260 :                                  ww(iw_loc) = ww(iw_loc) + cc*wrepp
     187       22260 :                                  iminus = i - 1
     188       26712 :                                  IF (iminus /= 0) THEN
     189       62328 :                                     DO j = 1, iminus
     190       44520 :                                        ij1 = inddd(i, j)
     191       44520 :                                        cc = ij_matrix%dd(ij1, i1 - 4, j1 - 4)
     192       44520 :                                        iw_loc = (indexb(i + 4, j + 4) - 1)*limkl + kl
     193       62328 :                                        ww(iw_loc) = ww(iw_loc) + cc*wrepp
     194             :                                     END DO
     195             :                                  END IF
     196             :                               END DO
     197             :                            END SELECT
     198             :                         END IF
     199             :                      END DO
     200             :                   END DO
     201             :                END DO
     202             :             END DO
     203             :          END IF
     204       51626 :          CALL rotmat_release(ij_matrix)
     205             : 
     206             :          ! Store two electron integrals in the triangular format
     207       51626 :          CALL store_2el_2c_diag(limij, limkl, ww, w)
     208             :       END IF
     209       51626 :    END SUBROUTINE rotint_num
     210             : 
     211             : ! **************************************************************************************************
     212             : !> \brief Calculates the derivative pf two-electron repulsion integrals and the
     213             : !>      nuclear attraction integrals w.r.t. |r|
     214             : !> \param sepi parameters of atom i
     215             : !> \param sepj parameters of atom j
     216             : !> \param rij interatomic distance
     217             : !> \param rep array of two-electron repulsion integrals
     218             : !> \param se_taper ...
     219             : !> \param se_int_control input parameters that control the calculation of SE
     220             : !>                         integrals (shortrange, R3 residual, screening type)
     221             : !> \par History
     222             : !>      03.2008 created [tlaino]
     223             : !> \author Teodoro Laino [tlaino] - Zurich University
     224             : ! **************************************************************************************************
     225       51626 :    SUBROUTINE terep_num(sepi, sepj, rij, rep, se_taper, se_int_control)
     226             :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
     227             :       REAL(dp), INTENT(IN)                               :: rij
     228             :       REAL(dp), DIMENSION(491), INTENT(OUT)              :: rep
     229             :       TYPE(se_taper_type), POINTER                       :: se_taper
     230             :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
     231             : 
     232             :       REAL(KIND=dp)                                      :: ft
     233             :       TYPE(se_int_screen_type)                           :: se_int_screen
     234             : 
     235       51626 :       ft = taper_eval(se_taper%taper, rij)
     236             :       ! In case of dumped integrals compute an additional taper term
     237       51626 :       IF (se_int_control%integral_screening == do_se_IS_kdso_d) THEN
     238           0 :          se_int_screen%ft = taper_eval(se_taper%taper_add, rij)
     239             :       END IF
     240             : 
     241             :       ! Contribution from sp shells
     242       51626 :       CALL terep_sp_num(sepi, sepj, rij, rep, se_int_control, se_int_screen, ft)
     243             : 
     244       51626 :       IF (sepi%dorb .OR. sepj%dorb) THEN
     245             :          ! Compute the contribution from d shells
     246             :          CALL terep_d_num(sepi, sepj, rij, rep, se_int_control, se_int_screen, &
     247          84 :                           ft)
     248             :       END IF
     249             : 
     250       51626 :    END SUBROUTINE terep_num
     251             : 
     252             : ! **************************************************************************************************
     253             : !> \brief Calculates the  two-electron repulsion integrals - sp shell only
     254             : !> \param sepi parameters of atom i
     255             : !> \param sepj parameters of atom j
     256             : !> \param rij ...
     257             : !> \param rep array of two-electron repulsion integrals
     258             : !> \param se_int_control input parameters that control the calculation of SE
     259             : !>                         integrals (shortrange, R3 residual, screening type)
     260             : !> \param se_int_screen contains information for computing the screened
     261             : !>                         integrals KDSO-D
     262             : !> \param ft ...
     263             : !> \par History
     264             : !>      Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
     265             : !>                 for computing integrals
     266             : ! **************************************************************************************************
     267     1197751 :    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     1197751 :       lasti = sepi%natorb
     281     1197751 :       lastj = sepj%natorb
     282     4525853 :       DO i = 1, MIN(lasti, 4)
     283     3328102 :          li = l_index(i)
     284    12114657 :          DO j = 1, i
     285     7588804 :             lj = l_index(j)
     286     7588804 :             ij = indexa(i, j)
     287    30167460 :             DO k = 1, MIN(lastj, 4)
     288    19250554 :                lk = l_index(k)
     289    69413412 :                DO l = 1, k
     290    42574054 :                   ll = l_index(l)
     291    42574054 :                   kl = indexa(k, l)
     292             : 
     293    42574054 :                   numb = ijkl_ind(ij, kl)
     294    61824608 :                   IF (numb > 0) THEN
     295    15499949 :                      nold = ijkl_sym(numb)
     296    15499949 :                      IF (nold > 0) THEN
     297     4968937 :                         rep(numb) = rep(nold)
     298    10531012 :                      ELSE IF (nold < 0) THEN
     299           0 :                         rep(numb) = -rep(-nold)
     300    10531012 :                      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    10531012 :                               *ft
     304    10531012 :                         rep(numb) = tmp
     305             :                      END IF
     306             :                   END IF
     307             :                END DO
     308             :             END DO
     309             :          END DO
     310             :       END DO
     311     1197751 :    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       55088 :    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       55088 :       lasti = sepi%natorb
     342       55088 :       lastj = sepj%natorb
     343      424574 :       DO i = 1, lasti
     344      369486 :          li = l_index(i)
     345     2085296 :          DO j = 1, i
     346     1660722 :             lj = l_index(j)
     347     1660722 :             ij = indexa(i, j)
     348     9657636 :             DO k = 1, lastj
     349     7627428 :                lk = l_index(k)
     350    37545630 :                DO l = 1, k
     351    28257480 :                   ll = l_index(l)
     352    28257480 :                   kl = indexa(k, l)
     353             : 
     354    28257480 :                   numb = ijkl_ind(ij, kl)
     355             :                   ! From 1 to 34 we store integrals involving sp shells
     356    35884908 :                   IF (numb > 34) THEN
     357     5922352 :                      nold = ijkl_sym(numb)
     358     5922352 :                      IF (nold > 34) THEN
     359     2762072 :                         rep(numb) = rep(nold)
     360     3160280 :                      ELSE IF (nold < -34) THEN
     361      538328 :                         rep(numb) = -rep(-nold)
     362     2621952 :                      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     2621952 :                               *ft
     366     2621952 :                         rep(numb) = tmp
     367             :                      END IF
     368             :                   END IF
     369             :                END DO
     370             :             END DO
     371             :          END DO
     372             :       END DO
     373             : 
     374       55088 :    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    42683706 :    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    42683706 :       l_core = PRESENT(core)
     852    42683706 :       l_ssss = PRESENT(ssss)
     853    42683706 :       IF (.NOT. (l_core .OR. l_ssss)) RETURN
     854    42683706 :       si = (sepi%natorb > 1)
     855    42683706 :       sj = (sepj%natorb > 1)
     856             : 
     857    42683706 :       ij = indexa(1, 1)
     858    42683706 :       IF (l_ssss) THEN
     859             :          ! Store the value for <S  S  | S  S  > (Used for computing the core-core interactions)
     860    26513442 :          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    42683706 :       IF (l_core) THEN
     864             :          !     <S  S  | S  S  >
     865    16170264 :          kl = indexa(1, 1)
     866    16170264 :          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    16170264 :          IF (si) THEN
     868             :             !  <S  P  | S  S  >
     869     7252842 :             kl = indexa(2, 1)
     870     7252842 :             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     7252842 :             kl = indexa(2, 2)
     873     7252842 :             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     7252842 :             kl = indexa(3, 3)
     876     7252842 :             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    16170264 :          kl = indexa(1, 1)
     881    16170264 :          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    16170264 :          IF (sj) THEN
     883             :             !  <S  S  | S  P  >
     884      423151 :             kl = indexa(2, 1)
     885      423151 :             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      423151 :             kl = indexa(2, 2)
     888      423151 :             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      423151 :             kl = indexa(3, 3)
     891      423151 :             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       38216 :    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       38216 :       IF (sepi%dorb .OR. sepj%dorb) THEN
     932       38216 :          ij = indexa(1, 1)
     933       38216 :          IF (sepj%dorb) THEN
     934             :             !  <S S | D S>
     935       21465 :             kl = indexa(5, 1)
     936       21465 :             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       21465 :             kl = indexa(5, 2)
     939       21465 :             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       21465 :             kl = indexa(5, 5)
     942       21465 :             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       21465 :             kl = indexa(6, 3)
     945       21465 :             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       21465 :             kl = indexa(6, 6)
     948       21465 :             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       21465 :             kl = indexa(8, 8)
     951       21465 :             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       38216 :          IF (sepi%dorb) THEN
     954             :             !  <D S | S S>
     955       22110 :             kl = indexa(5, 1)
     956       22110 :             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       22110 :             kl = indexa(5, 2)
     959       22110 :             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       22110 :             kl = indexa(5, 5)
     962       22110 :             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       22110 :             kl = indexa(6, 3)
     965       22110 :             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       22110 :             kl = indexa(6, 6)
     968       22110 :             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       22110 :             kl = indexa(8, 8)
     971       22110 :             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       38216 :    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 1.15