LCOV - code coverage report
Current view: top level - src - semi_empirical_int_ana.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 95.2 % 724 689
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 13 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 Analytical derivatives of Integrals for semi-empirical methods
      10              : !> \author Teodoro Laino - Zurich University 04.2007 [tlaino]
      11              : !> \par History
      12              : !>      23.11.2007 jhu   short range version of integrals
      13              : !>      Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
      14              : !>                 for computing integrals
      15              : !>      Teodoro Laino (05.2008) [tlaino] - University of Zurich : analytical
      16              : !>                 derivatives for d-orbitals
      17              : ! **************************************************************************************************
      18              : MODULE semi_empirical_int_ana
      19              : 
      20              :    USE input_constants, ONLY: do_method_am1, &
      21              :                               do_method_pchg, &
      22              :                               do_method_pdg, &
      23              :                               do_method_pm3, &
      24              :                               do_method_pm6, &
      25              :                               do_method_pm6fm, &
      26              :                               do_method_undef, &
      27              :                               do_se_IS_kdso_d
      28              :    USE kinds, ONLY: dp
      29              :    USE multipole_types, ONLY: do_multipole_none
      30              :    USE physcon, ONLY: angstrom, &
      31              :                       evolt
      32              :    USE semi_empirical_int_arrays, ONLY: &
      33              :       fac_x_to_z, ijkl_ind, ijkl_sym, inddd, inddp, indexa, indexb, indpp, int2c_type, l_index, &
      34              :       map_x_to_z, rij_threshold
      35              :    USE semi_empirical_int_num, ONLY: nucint_d_num, &
      36              :                                      nucint_sp_num, &
      37              :                                      terep_d_num, &
      38              :                                      terep_sp_num
      39              :    USE semi_empirical_int_utils, ONLY: d_ijkl_d, &
      40              :                                        d_ijkl_sp, &
      41              :                                        rot_2el_2c_first, &
      42              :                                        rotmat, &
      43              :                                        store_2el_2c_diag
      44              :    USE semi_empirical_types, ONLY: rotmat_create, &
      45              :                                    rotmat_release, &
      46              :                                    rotmat_type, &
      47              :                                    se_int_control_type, &
      48              :                                    se_int_screen_type, &
      49              :                                    se_taper_type, &
      50              :                                    semi_empirical_type, &
      51              :                                    setup_se_int_control_type
      52              :    USE taper_types, ONLY: dtaper_eval, &
      53              :                           taper_eval
      54              : #include "./base/base_uses.f90"
      55              : 
      56              :    IMPLICIT NONE
      57              :    PRIVATE
      58              : 
      59              :    #:include 'semi_empirical_int_debug.fypp'
      60              : 
      61              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_int_ana'
      62              :    LOGICAL, PARAMETER, PRIVATE          :: debug_this_module = .FALSE.
      63              :    PUBLIC :: rotnuc_ana, rotint_ana, corecore_ana, corecore_el_ana
      64              : 
      65              : CONTAINS
      66              : 
      67              : ! **************************************************************************************************
      68              : !> \brief Computes analytical gradients for semiempirical integrals
      69              : !> \param sepi Atomic parameters of first atom
      70              : !> \param sepj Atomic parameters of second atom
      71              : !> \param rijv Coordinate vector i -> j
      72              : !> \param itype ...
      73              : !> \param e1b Array of electron-nuclear attraction integrals, Electron on atom ni attracting nucleus of nj.
      74              : !> \param e2a Array of electron-nuclear attraction integrals, Electron on atom nj attracting nucleus of ni.
      75              : !> \param de1b derivative of e1b term
      76              : !> \param de2a derivative of e2a term
      77              : !> \param se_int_control input parameters that control the calculation of SE
      78              : !>                           integrals (shortrange, R3 residual, screening type)
      79              : !> \param se_taper ...
      80              : !> \par History
      81              : !>      04.2007 created [tlaino]
      82              : !>      Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
      83              : !>                 for computing integrals
      84              : !>      Teodoro Laino [tlaino] - University of Zurich 04.2008 : removed the core-core part
      85              : !> \author Teodoro Laino [tlaino] - Zurich University
      86              : !> \note
      87              : !>      Analytical version of the MOPAC rotnuc routine
      88              : ! **************************************************************************************************
      89     16142300 :    RECURSIVE SUBROUTINE rotnuc_ana(sepi, sepj, rijv, itype, e1b, e2a, de1b, de2a, &
      90              :                                    se_int_control, se_taper)
      91              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
      92              :       REAL(dp), DIMENSION(3), INTENT(IN)                 :: rijv
      93              :       INTEGER, INTENT(IN)                                :: itype
      94              :       REAL(dp), DIMENSION(45), INTENT(OUT), OPTIONAL     :: e1b, e2a
      95              :       REAL(dp), DIMENSION(3, 45), INTENT(OUT), OPTIONAL  :: de1b, de2a
      96              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
      97              :       TYPE(se_taper_type), POINTER                       :: se_taper
      98              : 
      99              :       INTEGER                                            :: i, idd, idp, ind1, ind2, ipp, j, &
     100              :                                                             last_orbital(2), m, n
     101              :       LOGICAL                                            :: invert, l_de1b, l_de2a, l_e1b, l_e2a, &
     102              :                                                             lgrad, task(2)
     103              :       REAL(KIND=dp)                                      :: rij, xtmp
     104              :       REAL(KIND=dp), DIMENSION(10, 2)                    :: core, dcore
     105              :       REAL(KIND=dp), DIMENSION(3)                        :: drij
     106              :       REAL(KIND=dp), DIMENSION(3, 45)                    :: tmp_d
     107              :       REAL(KIND=dp), DIMENSION(45)                       :: tmp
     108              :       TYPE(rotmat_type), POINTER                         :: ij_matrix
     109              : 
     110     16142300 :       NULLIFY (ij_matrix)
     111     64569200 :       rij = DOT_PRODUCT(rijv, rijv)
     112              :       ! Initialization
     113     16142300 :       l_e1b = PRESENT(e1b)
     114     16142300 :       l_e2a = PRESENT(e2a)
     115     16142300 :       l_de1b = PRESENT(de1b)
     116     16142300 :       l_de2a = PRESENT(de2a)
     117     16142300 :       lgrad = l_de1b .OR. l_de2a
     118              : 
     119     16142300 :       IF (rij > rij_threshold) THEN
     120              :          ! Compute Integrals in diatomic frame opportunely inverted
     121     16142300 :          rij = SQRT(rij)
     122              :          ! Create the rotation matrix
     123     16142300 :          CALL rotmat_create(ij_matrix)
     124     16142300 :          CALL rotmat(sepi, sepj, rijv, rij, ij_matrix, do_derivatives=lgrad, do_invert=invert)
     125              : 
     126     16142300 :          IF (lgrad) THEN
     127      8014152 :             drij(1) = rijv(1)/rij
     128      8014152 :             drij(2) = rijv(2)/rij
     129      8014152 :             drij(3) = rijv(3)/rij
     130              :             ! Possibly Invert Frame
     131      8014152 :             IF (invert) THEN
     132          879 :                xtmp = drij(3)
     133          879 :                drij(3) = drij(1)
     134          879 :                drij(1) = xtmp
     135              :             END IF
     136              :          END IF
     137              : 
     138              :          CALL dcore_nucint_ana(sepi, sepj, rij, core=core, dcore=dcore, itype=itype, se_taper=se_taper, &
     139     16142300 :                                se_int_control=se_int_control, lgrad=lgrad)
     140              : 
     141              :          ! Copy parameters over to arrays for do loop.
     142     16142300 :          last_orbital(1) = sepi%natorb
     143     16142300 :          last_orbital(2) = sepj%natorb
     144     16142300 :          task(1) = l_e1b
     145     16142300 :          task(2) = l_e2a
     146     48426900 :          DO n = 1, 2
     147     32284600 :             IF (.NOT. task(n)) CYCLE
     148     28902913 :             DO i = 1, last_orbital(n)
     149     20318894 :                ind1 = i - 1
     150     73155897 :                DO j = 1, i
     151     44252984 :                   ind2 = j - 1
     152     44252984 :                   m = (i*(i - 1))/2 + j
     153              :                   ! Perform Rotations ...
     154     64571878 :                   IF (ind2 == 0) THEN
     155     20318894 :                      IF (ind1 == 0) THEN
     156              :                         ! Type of Integral (SS/)
     157      8584019 :                         tmp(m) = core(1, n)
     158     11734875 :                      ELSE IF (ind1 < 4) THEN
     159              :                         ! Type of Integral (SP/)
     160     11618790 :                         tmp(m) = ij_matrix%sp(1, ind1)*core(2, n)
     161              :                      ELSE
     162              :                         ! Type of Integral (SD/)
     163       116085 :                         tmp(m) = ij_matrix%sd(1, ind1 - 3)*core(5, n)
     164              :                      END IF
     165     23934090 :                   ELSE IF (ind2 < 4) THEN
     166     23585835 :                      IF (ind1 < 4) THEN
     167              :                         ! Type of Integral (PP/)
     168     23237580 :                         ipp = indpp(ind1, ind2)
     169              :                         tmp(m) = core(3, n)*ij_matrix%pp(ipp, 1, 1) + &
     170     23237580 :                                  core(4, n)*(ij_matrix%pp(ipp, 2, 2) + ij_matrix%pp(ipp, 3, 3))
     171              :                      ELSE
     172              :                         ! Type of Integral (PD/)
     173       348255 :                         idp = inddp(ind1 - 3, ind2)
     174              :                         tmp(m) = core(6, n)*ij_matrix%pd(idp, 1, 1) + &
     175       348255 :                                  core(8, n)*(ij_matrix%pd(idp, 2, 2) + ij_matrix%pd(idp, 3, 3))
     176              :                      END IF
     177              :                   ELSE
     178              :                      ! Type of Integral (DD/)
     179       348255 :                      idd = inddd(ind1 - 3, ind2 - 3)
     180              :                      tmp(m) = core(7, n)*ij_matrix%dd(idd, 1, 1) + &
     181              :                               core(9, n)*(ij_matrix%dd(idd, 2, 2) + ij_matrix%dd(idd, 3, 3)) + &
     182       348255 :                               core(10, n)*(ij_matrix%dd(idd, 4, 4) + ij_matrix%dd(idd, 5, 5))
     183              :                   END IF
     184              :                END DO
     185              :             END DO
     186      8584019 :             IF (n == 1) THEN
     187     49462846 :                DO i = 1, sepi%atm_int_size
     188     49462846 :                   e1b(i) = -tmp(i)
     189              :                END DO
     190              :             END IF
     191     24726319 :             IF (n == 2) THEN
     192      3374157 :                DO i = 1, sepj%atm_int_size
     193      3374157 :                   e2a(i) = -tmp(i)
     194              :                END DO
     195              :             END IF
     196              :          END DO
     197     16142300 :          IF (invert .AND. l_e1b) CALL invert_integral(sepi, sepi, int1el=e1b)
     198     16142300 :          IF (invert .AND. l_e2a) CALL invert_integral(sepj, sepj, int1el=e2a)
     199              : 
     200              :          ! Possibly compute derivatives
     201     16142300 :          task(1) = l_de1b
     202     16142300 :          task(2) = l_de2a
     203     48426900 :          DO n = 1, 2
     204     32284600 :             IF (.NOT. task(n)) CYCLE
     205     28141050 :             DO i = 1, last_orbital(n)
     206     19784366 :                ind1 = i - 1
     207     71181920 :                DO j = 1, i
     208     43040870 :                   ind2 = j - 1
     209     43040870 :                   m = (i*(i - 1))/2 + j
     210              :                   ! Perform Rotations ...
     211     62825236 :                   IF (ind2 == 0) THEN
     212     19784366 :                      IF (ind1 == 0) THEN
     213              :                         ! Type of Integral (SS/)
     214      8356684 :                         tmp_d(1, m) = dcore(1, n)*drij(1)
     215      8356684 :                         tmp_d(2, m) = dcore(1, n)*drij(2)
     216      8356684 :                         tmp_d(3, m) = dcore(1, n)*drij(3)
     217     11427682 :                      ELSE IF (ind1 < 4) THEN
     218              :                         ! Type of Integral (SP/)
     219              :                         tmp_d(1, m) = ij_matrix%sp_d(1, 1, ind1)*core(2, n) + &
     220     11327397 :                                       ij_matrix%sp(1, ind1)*dcore(2, n)*drij(1)
     221              : 
     222              :                         tmp_d(2, m) = ij_matrix%sp_d(2, 1, ind1)*core(2, n) + &
     223     11327397 :                                       ij_matrix%sp(1, ind1)*dcore(2, n)*drij(2)
     224              : 
     225              :                         tmp_d(3, m) = ij_matrix%sp_d(3, 1, ind1)*core(2, n) + &
     226     11327397 :                                       ij_matrix%sp(1, ind1)*dcore(2, n)*drij(3)
     227              :                      ELSE
     228              :                         ! Type of Integral (SD/)
     229              :                         tmp_d(1, m) = ij_matrix%sd_d(1, 1, ind1 - 3)*core(5, n) + &
     230       100285 :                                       ij_matrix%sd(1, ind1 - 3)*dcore(5, n)*drij(1)
     231              : 
     232              :                         tmp_d(2, m) = ij_matrix%sd_d(2, 1, ind1 - 3)*core(5, n) + &
     233       100285 :                                       ij_matrix%sd(1, ind1 - 3)*dcore(5, n)*drij(2)
     234              : 
     235              :                         tmp_d(3, m) = ij_matrix%sd_d(3, 1, ind1 - 3)*core(5, n) + &
     236       100285 :                                       ij_matrix%sd(1, ind1 - 3)*dcore(5, n)*drij(3)
     237              :                      END IF
     238     23256504 :                   ELSE IF (ind2 < 4) THEN
     239     22955649 :                      IF (ind1 < 4) THEN
     240              :                         ! Type of Integral (PP/)
     241     22654794 :                         ipp = indpp(ind1, ind2)
     242              :                         tmp_d(1, m) = dcore(3, n)*drij(1)*ij_matrix%pp(ipp, 1, 1) + &
     243              :                                       core(3, n)*ij_matrix%pp_d(1, ipp, 1, 1) + &
     244              :                                       dcore(4, n)*drij(1)*(ij_matrix%pp(ipp, 2, 2) + ij_matrix%pp(ipp, 3, 3)) + &
     245     22654794 :                                       core(4, n)*(ij_matrix%pp_d(1, ipp, 2, 2) + ij_matrix%pp_d(1, ipp, 3, 3))
     246              : 
     247              :                         tmp_d(2, m) = dcore(3, n)*drij(2)*ij_matrix%pp(ipp, 1, 1) + &
     248              :                                       core(3, n)*ij_matrix%pp_d(2, ipp, 1, 1) + &
     249              :                                       dcore(4, n)*drij(2)*(ij_matrix%pp(ipp, 2, 2) + ij_matrix%pp(ipp, 3, 3)) + &
     250     22654794 :                                       core(4, n)*(ij_matrix%pp_d(2, ipp, 2, 2) + ij_matrix%pp_d(2, ipp, 3, 3))
     251              : 
     252              :                         tmp_d(3, m) = dcore(3, n)*drij(3)*ij_matrix%pp(ipp, 1, 1) + &
     253              :                                       core(3, n)*ij_matrix%pp_d(3, ipp, 1, 1) + &
     254              :                                       dcore(4, n)*drij(3)*(ij_matrix%pp(ipp, 2, 2) + ij_matrix%pp(ipp, 3, 3)) + &
     255     22654794 :                                       core(4, n)*(ij_matrix%pp_d(3, ipp, 2, 2) + ij_matrix%pp_d(3, ipp, 3, 3))
     256              :                      ELSE
     257              :                         ! Type of Integral (PD/)
     258       300855 :                         idp = inddp(ind1 - 3, ind2)
     259              :                         tmp_d(1, m) = dcore(6, n)*drij(1)*ij_matrix%pd(idp, 1, 1) + &
     260              :                                       core(6, n)*ij_matrix%pd_d(1, idp, 1, 1) + &
     261              :                                       dcore(8, n)*drij(1)*(ij_matrix%pd(idp, 2, 2) + ij_matrix%pd(idp, 3, 3)) + &
     262       300855 :                                       core(8, n)*(ij_matrix%pd_d(1, idp, 2, 2) + ij_matrix%pd_d(1, idp, 3, 3))
     263              : 
     264              :                         tmp_d(2, m) = dcore(6, n)*drij(2)*ij_matrix%pd(idp, 1, 1) + &
     265              :                                       core(6, n)*ij_matrix%pd_d(2, idp, 1, 1) + &
     266              :                                       dcore(8, n)*drij(2)*(ij_matrix%pd(idp, 2, 2) + ij_matrix%pd(idp, 3, 3)) + &
     267       300855 :                                       core(8, n)*(ij_matrix%pd_d(2, idp, 2, 2) + ij_matrix%pd_d(2, idp, 3, 3))
     268              : 
     269              :                         tmp_d(3, m) = dcore(6, n)*drij(3)*ij_matrix%pd(idp, 1, 1) + &
     270              :                                       core(6, n)*ij_matrix%pd_d(3, idp, 1, 1) + &
     271              :                                       dcore(8, n)*drij(3)*(ij_matrix%pd(idp, 2, 2) + ij_matrix%pd(idp, 3, 3)) + &
     272       300855 :                                       core(8, n)*(ij_matrix%pd_d(3, idp, 2, 2) + ij_matrix%pd_d(3, idp, 3, 3))
     273              :                      END IF
     274              :                   ELSE
     275              :                      ! Type of Integral (DD/)
     276       300855 :                      idd = inddd(ind1 - 3, ind2 - 3)
     277              :                      tmp_d(1, m) = dcore(7, n)*drij(1)*ij_matrix%dd(idd, 1, 1) + &
     278              :                                    core(7, n)*ij_matrix%dd_d(1, idd, 1, 1) + &
     279              :                                    dcore(9, n)*drij(1)*(ij_matrix%dd(idd, 2, 2) + ij_matrix%dd(idd, 3, 3)) + &
     280              :                                    core(9, n)*(ij_matrix%dd_d(1, idd, 2, 2) + ij_matrix%dd_d(1, idd, 3, 3)) + &
     281              :                                    dcore(10, n)*drij(1)*(ij_matrix%dd(idd, 4, 4) + ij_matrix%dd(idd, 5, 5)) + &
     282       300855 :                                    core(10, n)*(ij_matrix%dd_d(1, idd, 4, 4) + ij_matrix%dd_d(1, idd, 5, 5))
     283              : 
     284              :                      tmp_d(2, m) = dcore(7, n)*drij(2)*ij_matrix%dd(idd, 1, 1) + &
     285              :                                    core(7, n)*ij_matrix%dd_d(2, idd, 1, 1) + &
     286              :                                    dcore(9, n)*drij(2)*(ij_matrix%dd(idd, 2, 2) + ij_matrix%dd(idd, 3, 3)) + &
     287              :                                    core(9, n)*(ij_matrix%dd_d(2, idd, 2, 2) + ij_matrix%dd_d(2, idd, 3, 3)) + &
     288              :                                    dcore(10, n)*drij(2)*(ij_matrix%dd(idd, 4, 4) + ij_matrix%dd(idd, 5, 5)) + &
     289       300855 :                                    core(10, n)*(ij_matrix%dd_d(2, idd, 4, 4) + ij_matrix%dd_d(2, idd, 5, 5))
     290              : 
     291              :                      tmp_d(3, m) = dcore(7, n)*drij(3)*ij_matrix%dd(idd, 1, 1) + &
     292              :                                    core(7, n)*ij_matrix%dd_d(3, idd, 1, 1) + &
     293              :                                    dcore(9, n)*drij(3)*(ij_matrix%dd(idd, 2, 2) + ij_matrix%dd(idd, 3, 3)) + &
     294              :                                    core(9, n)*(ij_matrix%dd_d(3, idd, 2, 2) + ij_matrix%dd_d(3, idd, 3, 3)) + &
     295              :                                    dcore(10, n)*drij(3)*(ij_matrix%dd(idd, 4, 4) + ij_matrix%dd(idd, 5, 5)) + &
     296       300855 :                                    core(10, n)*(ij_matrix%dd_d(3, idd, 4, 4) + ij_matrix%dd_d(3, idd, 5, 5))
     297              :                   END IF
     298              :                END DO
     299              :             END DO
     300      8356684 :             IF (n == 1) THEN
     301     48741983 :                DO i = 1, sepi%atm_int_size
     302     40727831 :                   de1b(1, i) = -tmp_d(1, i)
     303     40727831 :                   de1b(2, i) = -tmp_d(2, i)
     304     48741983 :                   de1b(3, i) = -tmp_d(3, i)
     305              :                END DO
     306              :             END IF
     307     24498984 :             IF (n == 2) THEN
     308      2655571 :                DO i = 1, sepj%atm_int_size
     309      2313039 :                   de2a(1, i) = -tmp_d(1, i)
     310      2313039 :                   de2a(2, i) = -tmp_d(2, i)
     311      2655571 :                   de2a(3, i) = -tmp_d(3, i)
     312              :                END DO
     313              :             END IF
     314              :          END DO
     315     16142300 :          IF (invert .AND. l_de1b) CALL invert_derivative(sepi, sepi, dint1el=de1b)
     316     16142300 :          IF (invert .AND. l_de2a) CALL invert_derivative(sepj, sepj, dint1el=de2a)
     317     16142300 :          CALL rotmat_release(ij_matrix)
     318              : 
     319              :          ! Possibly debug the analytical values versus the numerical ones
     320              :          IF (debug_this_module) THEN
     321              :             CALL check_drotnuc_ana(sepi, sepj, rijv, itype, se_int_control, se_taper, e1b, e2a, de1b, de2a)
     322              :          END IF
     323              :       END IF
     324     16142300 :    END SUBROUTINE rotnuc_ana
     325              : 
     326              : ! **************************************************************************************************
     327              : !> \brief Computes analytical gradients for semiempirical core-core interaction.
     328              : !> \param sepi Atomic parameters of first atom
     329              : !> \param sepj Atomic parameters of second atom
     330              : !> \param rijv Coordinate vector i -> j
     331              : !> \param itype ...
     332              : !> \param enuc nuclear-nuclear repulsion term.
     333              : !> \param denuc derivative of nuclear-nuclear repulsion term.
     334              : !> \param se_int_control input parameters that control the calculation of SE
     335              : !>                           integrals (shortrange, R3 residual, screening type)
     336              : !> \param se_taper ...
     337              : !> \par History
     338              : !>      04.2007 created [tlaino]
     339              : !>      Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
     340              : !>                 for computing integrals
     341              : !>       Teodoro Laino [tlaino] - University of Zurich 04.2008 : removed the
     342              : !>                 core-core part
     343              : !> \author Teodoro Laino [tlaino] - Zurich University
     344              : !> \note
     345              : !>      Analytical version of the MOPAC rotnuc routine
     346              : ! **************************************************************************************************
     347     23294450 :    RECURSIVE SUBROUTINE corecore_ana(sepi, sepj, rijv, itype, enuc, denuc, se_int_control, se_taper)
     348              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
     349              :       REAL(dp), DIMENSION(3), INTENT(IN)                 :: rijv
     350              :       INTEGER, INTENT(IN)                                :: itype
     351              :       REAL(dp), INTENT(OUT), OPTIONAL                    :: enuc
     352              :       REAL(dp), DIMENSION(3), INTENT(OUT), OPTIONAL      :: denuc
     353              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
     354              :       TYPE(se_taper_type), POINTER                       :: se_taper
     355              : 
     356              :       INTEGER                                            :: ig, nt
     357              :       LOGICAL                                            :: l_denuc, l_enuc
     358              :       REAL(dp) :: aab, alpi, alpj, apdg, ax, dai, daj, dax, dbi, dbj, denuc_loc, dqcorr, drija, &
     359              :                   dscale, dssss, dssss_sr, dtmp, dzz, enuc_loc, pai, paj, pbi, pbj, qcorr, rij, rija, &
     360              :                   scale, ssss, ssss_sr, tmp, xab, xtmp, zaf, zbf, zz
     361              :       REAL(dp), DIMENSION(3)                             :: drij
     362              :       REAL(dp), DIMENSION(4)                             :: fni1, fni2, fni3, fnj1, fnj2, fnj3
     363              :       TYPE(se_int_control_type)                          :: se_int_control_off
     364              : 
     365     93177800 :       rij = DOT_PRODUCT(rijv, rijv)
     366              :       ! Initialization
     367     23294450 :       l_enuc = PRESENT(enuc)
     368     23294450 :       l_denuc = PRESENT(denuc)
     369     23294450 :       IF ((rij > rij_threshold) .AND. (l_enuc .OR. l_denuc)) THEN
     370              :          ! Compute Integrals in diatomic frame
     371     23294450 :          rij = SQRT(rij)
     372              :          CALL setup_se_int_control_type(se_int_control_off, shortrange=.FALSE., do_ewald_r3=.FALSE., &
     373              :                                         do_ewald_gks=.FALSE., integral_screening=se_int_control%integral_screening, &
     374     23294450 :                                         max_multipole=do_multipole_none, pc_coulomb_int=.FALSE.)
     375              :          CALL dssss_nucint_ana(sepi, sepj, rij, ssss=ssss, dssss=dssss, itype=itype, se_taper=se_taper, &
     376     23294450 :                                se_int_control=se_int_control_off, lgrad=l_denuc)
     377              :          ! In case let's compute the short-range part of the (ss|ss) integral
     378     23294450 :          IF (se_int_control%shortrange) THEN
     379              :             CALL dssss_nucint_ana(sepi, sepj, rij, ssss=ssss_sr, dssss=dssss_sr, itype=itype, &
     380       248393 :                                   se_taper=se_taper, se_int_control=se_int_control, lgrad=l_denuc)
     381              :          ELSE
     382     23046057 :             ssss_sr = ssss
     383     23046057 :             dssss_sr = dssss
     384              :          END IF
     385              :          ! Zeroing local method dependent core-core corrections
     386     23294450 :          enuc_loc = 0.0_dp
     387     23294450 :          denuc_loc = 0.0_dp
     388     23294450 :          qcorr = 0.0_dp
     389     23294450 :          scale = 0.0_dp
     390     23294450 :          dscale = 0.0_dp
     391     23294450 :          dqcorr = 0.0_dp
     392     23294450 :          zz = sepi%zeff*sepj%zeff
     393              :          ! Core Core electrostatic contribution
     394     23294450 :          IF (l_enuc) enuc_loc = zz*ssss_sr
     395     23294450 :          IF (l_denuc) denuc_loc = zz*dssss_sr
     396              :          ! Method dependent code
     397     23294450 :          tmp = zz*ssss
     398     23294450 :          IF (l_denuc) dtmp = zz*dssss
     399     23294450 :          IF (itype /= do_method_pm6 .AND. itype /= do_method_pm6fm) THEN
     400     15726037 :             alpi = sepi%alp
     401     15726037 :             alpj = sepj%alp
     402     15726037 :             scale = EXP(-alpi*rij) + EXP(-alpj*rij)
     403     15726037 :             IF (l_denuc) THEN
     404      7796641 :                dscale = -alpi*EXP(-alpi*rij) - alpj*EXP(-alpj*rij)
     405              :             END IF
     406     15726037 :             nt = sepi%z + sepj%z
     407     15726037 :             IF (nt == 8 .OR. nt == 9) THEN
     408      2993026 :                IF (sepi%z == 7 .OR. sepi%z == 8) THEN
     409      2961656 :                   scale = scale + (angstrom*rij - 1._dp)*EXP(-alpi*rij)
     410      2961656 :                   IF (l_denuc) THEN
     411      1475382 :                      dscale = dscale + angstrom*EXP(-alpi*rij) - (angstrom*rij - 1._dp)*alpi*EXP(-alpi*rij)
     412              :                   END IF
     413              :                END IF
     414      2993026 :                IF (sepj%z == 7 .OR. sepj%z == 8) THEN
     415        31370 :                   scale = scale + (angstrom*rij - 1._dp)*EXP(-alpj*rij)
     416        31370 :                   IF (l_denuc) THEN
     417        10395 :                      dscale = dscale + angstrom*EXP(-alpj*rij) - (angstrom*rij - 1._dp)*alpj*EXP(-alpj*rij)
     418              :                   END IF
     419              :                END IF
     420              :             END IF
     421     15726037 :             IF (l_denuc) THEN
     422      7796641 :                dscale = SIGN(1.0_dp, scale*tmp)*(dscale*tmp + scale*dtmp)
     423      7796641 :                dzz = -zz/rij**2
     424              :             END IF
     425     15726037 :             scale = ABS(scale*tmp)
     426     15726037 :             zz = zz/rij
     427     15726037 :             IF (itype == do_method_am1 .OR. itype == do_method_pm3 .OR. itype == do_method_pdg) THEN
     428       302787 :                IF (itype == do_method_am1 .AND. sepi%z == 5) THEN
     429              :                   !special case AM1 Boron
     430              :                   SELECT CASE (sepj%z)
     431              :                   CASE DEFAULT
     432            0 :                      nt = 1
     433              :                   CASE (1)
     434            0 :                      nt = 2
     435              :                   CASE (6)
     436            0 :                      nt = 3
     437              :                   CASE (9, 17, 35, 53)
     438            0 :                      nt = 4
     439              :                   END SELECT
     440            0 :                   fni1(1) = sepi%bfn1(1, nt)
     441            0 :                   fni1(2) = sepi%bfn1(2, nt)
     442            0 :                   fni1(3) = sepi%bfn1(3, nt)
     443            0 :                   fni1(4) = sepi%bfn1(4, nt)
     444            0 :                   fni2(1) = sepi%bfn2(1, nt)
     445            0 :                   fni2(2) = sepi%bfn2(2, nt)
     446            0 :                   fni2(3) = sepi%bfn2(3, nt)
     447            0 :                   fni2(4) = sepi%bfn2(4, nt)
     448            0 :                   fni3(1) = sepi%bfn3(1, nt)
     449            0 :                   fni3(2) = sepi%bfn3(2, nt)
     450            0 :                   fni3(3) = sepi%bfn3(3, nt)
     451            0 :                   fni3(4) = sepi%bfn3(4, nt)
     452              :                ELSE
     453       302787 :                   fni1(1) = sepi%fn1(1)
     454       302787 :                   fni1(2) = sepi%fn1(2)
     455       302787 :                   fni1(3) = sepi%fn1(3)
     456       302787 :                   fni1(4) = sepi%fn1(4)
     457       302787 :                   fni2(1) = sepi%fn2(1)
     458       302787 :                   fni2(2) = sepi%fn2(2)
     459       302787 :                   fni2(3) = sepi%fn2(3)
     460       302787 :                   fni2(4) = sepi%fn2(4)
     461       302787 :                   fni3(1) = sepi%fn3(1)
     462       302787 :                   fni3(2) = sepi%fn3(2)
     463       302787 :                   fni3(3) = sepi%fn3(3)
     464       302787 :                   fni3(4) = sepi%fn3(4)
     465              :                END IF
     466       302787 :                IF (itype == do_method_am1 .AND. sepj%z == 5) THEN
     467              :                   !special case AM1 Boron
     468              :                   SELECT CASE (sepi%z)
     469              :                   CASE DEFAULT
     470            0 :                      nt = 1
     471              :                   CASE (1)
     472            0 :                      nt = 2
     473              :                   CASE (6)
     474            0 :                      nt = 3
     475              :                   CASE (9, 17, 35, 53)
     476            0 :                      nt = 4
     477              :                   END SELECT
     478            0 :                   fnj1(1) = sepj%bfn1(1, nt)
     479            0 :                   fnj1(2) = sepj%bfn1(2, nt)
     480            0 :                   fnj1(3) = sepj%bfn1(3, nt)
     481            0 :                   fnj1(4) = sepj%bfn1(4, nt)
     482            0 :                   fnj2(1) = sepj%bfn2(1, nt)
     483            0 :                   fnj2(2) = sepj%bfn2(2, nt)
     484            0 :                   fnj2(3) = sepj%bfn2(3, nt)
     485            0 :                   fnj2(4) = sepj%bfn2(4, nt)
     486            0 :                   fnj3(1) = sepj%bfn3(1, nt)
     487            0 :                   fnj3(2) = sepj%bfn3(2, nt)
     488            0 :                   fnj3(3) = sepj%bfn3(3, nt)
     489            0 :                   fnj3(4) = sepj%bfn3(4, nt)
     490              :                ELSE
     491       302787 :                   fnj1(1) = sepj%fn1(1)
     492       302787 :                   fnj1(2) = sepj%fn1(2)
     493       302787 :                   fnj1(3) = sepj%fn1(3)
     494       302787 :                   fnj1(4) = sepj%fn1(4)
     495       302787 :                   fnj2(1) = sepj%fn2(1)
     496       302787 :                   fnj2(2) = sepj%fn2(2)
     497       302787 :                   fnj2(3) = sepj%fn2(3)
     498       302787 :                   fnj2(4) = sepj%fn2(4)
     499       302787 :                   fnj3(1) = sepj%fn3(1)
     500       302787 :                   fnj3(2) = sepj%fn3(2)
     501       302787 :                   fnj3(3) = sepj%fn3(3)
     502       302787 :                   fnj3(4) = sepj%fn3(4)
     503              :                END IF
     504              :                ! AM1/PM3/PDG correction to nuclear repulsion
     505      1513935 :                DO ig = 1, SIZE(fni1)
     506      1211148 :                   IF (ABS(fni1(ig)) > 0._dp) THEN
     507       836741 :                      ax = fni2(ig)*(rij - fni3(ig))**2
     508       836741 :                      IF (ax <= 25._dp) THEN
     509       238319 :                         scale = scale + zz*fni1(ig)*EXP(-ax)
     510       238319 :                         IF (l_denuc) THEN
     511        79200 :                            dax = fni2(ig)*2.0_dp*(rij - fni3(ig))
     512        79200 :                            dscale = dscale + dzz*fni1(ig)*EXP(-ax) - dax*zz*fni1(ig)*EXP(-ax)
     513              :                         END IF
     514              :                      END IF
     515              :                   END IF
     516      1513935 :                   IF (ABS(fnj1(ig)) > 0._dp) THEN
     517       836729 :                      ax = fnj2(ig)*(rij - fnj3(ig))**2
     518       836729 :                      IF (ax <= 25._dp) THEN
     519       237956 :                         scale = scale + zz*fnj1(ig)*EXP(-ax)
     520       237956 :                         IF (l_denuc) THEN
     521        79084 :                            dax = fnj2(ig)*2.0_dp*(rij - fnj3(ig))
     522        79084 :                            dscale = dscale + dzz*fnj1(ig)*EXP(-ax) - dax*zz*fnj1(ig)*EXP(-ax)
     523              :                         END IF
     524              :                      END IF
     525              :                   END IF
     526              :                END DO
     527              :             END IF
     528     15726037 :             IF (itype == do_method_pdg) THEN
     529              :                ! PDDG function
     530            9 :                zaf = sepi%zeff/nt
     531            9 :                zbf = sepj%zeff/nt
     532            9 :                pai = sepi%pre(1)
     533            9 :                pbi = sepi%pre(2)
     534            9 :                paj = sepj%pre(1)
     535            9 :                pbj = sepj%pre(2)
     536            9 :                dai = sepi%d(1)
     537            9 :                dbi = sepi%d(2)
     538            9 :                daj = sepj%d(1)
     539            9 :                dbj = sepj%d(2)
     540            9 :                apdg = 10._dp*angstrom**2
     541              :                qcorr = (zaf*pai + zbf*paj)*EXP(-apdg*(rij - dai - daj)**2) + &
     542              :                        (zaf*pai + zbf*pbj)*EXP(-apdg*(rij - dai - dbj)**2) + &
     543              :                        (zaf*pbi + zbf*paj)*EXP(-apdg*(rij - dbi - daj)**2) + &
     544            9 :                        (zaf*pbi + zbf*pbj)*EXP(-apdg*(rij - dbi - dbj)**2)
     545            9 :                IF (l_denuc) THEN
     546              :                   dqcorr = (zaf*pai + zbf*paj)*EXP(-apdg*(rij - dai - daj)**2)*(-2.0_dp*apdg*(rij - dai - daj)) + &
     547              :                            (zaf*pai + zbf*pbj)*EXP(-apdg*(rij - dai - dbj)**2)*(-2.0_dp*apdg*(rij - dai - dbj)) + &
     548              :                            (zaf*pbi + zbf*paj)*EXP(-apdg*(rij - dbi - daj)**2)*(-2.0_dp*apdg*(rij - dbi - daj)) + &
     549            3 :                            (zaf*pbi + zbf*pbj)*EXP(-apdg*(rij - dbi - dbj)**2)*(-2.0_dp*apdg*(rij - dbi - dbj))
     550              :                END IF
     551     15726028 :             ELSEIF (itype == do_method_pchg) THEN
     552              :                qcorr = 0.0_dp
     553              :                scale = 0.0_dp
     554              :                dscale = 0.0_dp
     555              :                dqcorr = 0.0_dp
     556              :             ELSE
     557       382131 :                qcorr = 0.0_dp
     558       382131 :                dqcorr = 0.0_dp
     559              :             END IF
     560              :          ELSE
     561              :             ! PM6 core-core terms
     562      7568413 :             scale = tmp
     563      7568413 :             IF (l_denuc) dscale = dtmp
     564      7568413 :             drija = angstrom
     565      7568413 :             rija = rij*drija
     566      7568413 :             xab = sepi%xab(sepj%z)
     567      7568413 :             aab = sepi%aab(sepj%z)
     568      7568413 :             IF ((sepi%z == 1 .AND. (sepj%z == 6 .OR. sepj%z == 7 .OR. sepj%z == 8)) .OR. &
     569              :                 (sepj%z == 1 .AND. (sepi%z == 6 .OR. sepi%z == 7 .OR. sepi%z == 8))) THEN
     570              :                ! Special Case O-H or N-H or C-H
     571       253301 :                IF (l_denuc) dscale = dscale*(2._dp*xab*EXP(-aab*rija*rija)) - &
     572        55301 :                                      scale*2._dp*xab*EXP(-aab*rija*rija)*(2.0_dp*aab*rija)*drija
     573       253301 :                IF (l_enuc) scale = scale*(2._dp*xab*EXP(-aab*rija*rija))
     574      7315112 :             ELSEIF (sepi%z == 6 .AND. sepj%z == 6) THEN
     575              :                ! Special Case C-C
     576       352747 :                IF (l_denuc) dscale = &
     577              :                   dscale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6)) + 9.28_dp*EXP(-5.98_dp*rija)) &
     578              :                   - scale*2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6))*aab*(1.0_dp + 6.0_dp*0.0003_dp*rija**5)*drija &
     579        12485 :                   - scale*9.28_dp*EXP(-5.98_dp*rija)*5.98_dp*drija
     580       352747 :                IF (l_enuc) scale = scale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6)) + 9.28_dp*EXP(-5.98_dp*rija))
     581      6962365 :             ELSEIF ((sepi%z == 8 .AND. sepj%z == 14) .OR. &
     582              :                     (sepj%z == 8 .AND. sepi%z == 14)) THEN
     583              :                ! Special Case Si-O
     584        97448 :                IF (l_denuc) dscale = &
     585              :                   dscale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6)) - 0.0007_dp*EXP(-(rija - 2.9_dp)**2)) &
     586              :                   - scale*2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6))*aab*(1.0_dp + 6.0_dp*0.0003_dp*rija**5)*drija + &
     587            2 :                   scale*0.0007_dp*EXP(-(rija - 2.9_dp)**2)*(2.0_dp*(rija - 2.9_dp)*drija)
     588        97448 :                IF (l_enuc) scale = scale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6)) - 0.0007_dp*EXP(-(rija - 2.9_dp)**2))
     589              :             ELSE
     590              :                ! General Case
     591              :                ! Factor of 2 found by experiment
     592      6864917 :                IF (l_denuc) dscale = dscale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6))) &
     593       105775 :                                 - scale*2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6))*aab*(1.0_dp + 6.0_dp*0.0003_dp*rija**5)*drija
     594      6864917 :                IF (l_enuc) scale = scale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6)))
     595              :             END IF
     596              :             ! General correction term a*exp(-b*(rij-c)^2)
     597      7568413 :             xtmp = 1.e-8_dp/evolt*((REAL(sepi%z, dp)**(1._dp/3._dp) + REAL(sepj%z, dp)**(1._dp/3._dp))/rija)**12
     598      7568413 :             IF (l_enuc) THEN
     599              :                qcorr = (sepi%a*EXP(-sepi%b*(rij - sepi%c)**2))*zz/rij + &
     600              :                        (sepj%a*EXP(-sepj%b*(rij - sepj%c)**2))*zz/rij + &
     601              :                        ! Hard core repulsion
     602      7394850 :                        xtmp
     603              :             END IF
     604      7568413 :             IF (l_denuc) THEN
     605              :                dqcorr = (sepi%a*EXP(-sepi%b*(rij - sepi%c)**2)*(-2.0_dp*sepi%b*(rij - sepi%c)))*zz/rij - &
     606              :                         (sepi%a*EXP(-sepi%b*(rij - sepi%c)**2))*zz/rij**2 + &
     607              :                         (sepj%a*EXP(-sepj%b*(rij - sepj%c)**2)*(-2.0_dp*sepj%b*(rij - sepj%c)))*zz/rij - &
     608              :                         (sepj%a*EXP(-sepj%b*(rij - sepj%c)**2))*zz/rij**2 + &
     609              :                         ! Hard core repulsion
     610       173563 :                         (-12.0_dp*xtmp/rija*drija)
     611              :             END IF
     612              :          END IF
     613              : 
     614              :          ! Only at the very end let's sum-up the several contributions energy/derivatives
     615              :          ! This assignment should be method independent
     616     23294450 :          IF (l_enuc) THEN
     617     15324246 :             enuc = enuc_loc + scale + qcorr
     618              :          END IF
     619     23294450 :          IF (l_denuc) THEN
     620      7970204 :             drij(1) = rijv(1)/rij
     621      7970204 :             drij(2) = rijv(2)/rij
     622      7970204 :             drij(3) = rijv(3)/rij
     623     31880816 :             denuc = (denuc_loc + dscale + dqcorr)*drij
     624              :          END IF
     625              :          ! Debug statement
     626              :          IF (debug_this_module) THEN
     627              :             CALL check_dcorecore_ana(sepi, sepj, rijv, itype, se_int_control, se_taper, enuc, denuc)
     628              :          END IF
     629              :       END IF
     630     23294450 :    END SUBROUTINE corecore_ana
     631              : 
     632              : ! **************************************************************************************************
     633              : !> \brief Computes analytical gradients for semiempirical core-core electrostatic
     634              : !>        interaction only.
     635              : !> \param sepi Atomic parameters of first atom
     636              : !> \param sepj Atomic parameters of second atom
     637              : !> \param rijv Coordinate vector i -> j
     638              : !> \param itype ...
     639              : !> \param enuc nuclear-nuclear electrostatic repulsion term.
     640              : !> \param denuc derivative of nuclear-nuclear electrostatic
     641              : !>                             repulsion term.
     642              : !> \param se_int_control input parameters that control the calculation of SE
     643              : !>                           integrals (shortrange, R3 residual, screening type)
     644              : !> \param se_taper ...
     645              : !> \par History
     646              : !>      04.2007 created [tlaino]
     647              : !>      Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
     648              : !>                 for computing integrals
     649              : !>       Teodoro Laino [tlaino] - University of Zurich 04.2008 : removed the
     650              : !>                 core-core part
     651              : !> \author Teodoro Laino [tlaino] - Zurich University
     652              : !> \note
     653              : !>      Analytical version of the MOPAC rotnuc routine
     654              : ! **************************************************************************************************
     655      1567031 :    RECURSIVE SUBROUTINE corecore_el_ana(sepi, sepj, rijv, itype, enuc, denuc, &
     656              :                                         se_int_control, se_taper)
     657              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
     658              :       REAL(dp), DIMENSION(3), INTENT(IN)                 :: rijv
     659              :       INTEGER, INTENT(IN)                                :: itype
     660              :       REAL(dp), INTENT(OUT), OPTIONAL                    :: enuc
     661              :       REAL(dp), DIMENSION(3), INTENT(OUT), OPTIONAL      :: denuc
     662              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
     663              :       TYPE(se_taper_type), POINTER                       :: se_taper
     664              : 
     665              :       LOGICAL                                            :: l_denuc, l_enuc
     666              :       REAL(dp)                                           :: drij(3), dssss, dssss_sr, rij, ssss, &
     667              :                                                             ssss_sr, tmp, zz
     668              :       TYPE(se_int_control_type)                          :: se_int_control_off
     669              : 
     670      6268124 :       rij = DOT_PRODUCT(rijv, rijv)
     671              :       ! Initialization
     672      1567031 :       l_enuc = PRESENT(enuc)
     673      1567031 :       l_denuc = PRESENT(denuc)
     674      1567031 :       IF ((rij > rij_threshold) .AND. (l_enuc .OR. l_denuc)) THEN
     675              :          ! Compute Integrals in diatomic frame
     676      1567031 :          rij = SQRT(rij)
     677              :          CALL setup_se_int_control_type(se_int_control_off, shortrange=.FALSE., do_ewald_r3=.FALSE., &
     678              :                                         do_ewald_gks=.FALSE., integral_screening=se_int_control%integral_screening, &
     679      1567031 :                                         max_multipole=do_multipole_none, pc_coulomb_int=.FALSE.)
     680              :          CALL dssss_nucint_ana(sepi, sepj, rij, ssss=ssss, dssss=dssss, itype=itype, se_taper=se_taper, &
     681      1567031 :                                se_int_control=se_int_control_off, lgrad=l_denuc)
     682              :          ! In case let's compute the short-range part of the (ss|ss) integral
     683      1567031 :          IF (se_int_control%shortrange .OR. se_int_control%pc_coulomb_int) THEN
     684              :             CALL dssss_nucint_ana(sepi, sepj, rij, ssss=ssss_sr, dssss=dssss_sr, itype=itype, &
     685      1567031 :                                   se_taper=se_taper, se_int_control=se_int_control, lgrad=l_denuc)
     686              :          ELSE
     687            0 :             ssss_sr = ssss
     688            0 :             dssss_sr = dssss
     689              :          END IF
     690      1567031 :          zz = sepi%zeff*sepj%zeff
     691              :          ! Core Core electrostatic contribution
     692      1567031 :          IF (l_enuc) enuc = zz*ssss_sr
     693      1567031 :          IF (l_denuc) THEN
     694        43876 :             drij(1) = rijv(1)/rij
     695        43876 :             drij(2) = rijv(2)/rij
     696        43876 :             drij(3) = rijv(3)/rij
     697        43876 :             tmp = zz*dssss_sr
     698       175504 :             denuc = tmp*drij
     699              :          END IF
     700              :       END IF
     701      1567031 :    END SUBROUTINE corecore_el_ana
     702              : 
     703              : ! **************************************************************************************************
     704              : !> \brief Exploits inversion symmetry to avoid divergence
     705              : !> \param sepi ...
     706              : !> \param sepj ...
     707              : !> \param int1el ...
     708              : !> \param int2el ...
     709              : !> \par History
     710              : !>      04.2007 created [tlaino]
     711              : !>      05.2008 New driver for integral invertion (supports d-orbitals)
     712              : !> \author Teodoro Laino - Zurich University
     713              : ! **************************************************************************************************
     714        13358 :    SUBROUTINE invert_integral(sepi, sepj, int1el, int2el)
     715              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
     716              :       REAL(dp), DIMENSION(:), INTENT(INOUT), OPTIONAL    :: int1el, int2el
     717              : 
     718              :       INTEGER                                            :: fdim, gind, gknd, i, imap, ind, j, jmap, &
     719              :                                                             jnd, k, kmap, knd, l, lmap, lnd, ndim, &
     720              :                                                             sdim, tdim, tind
     721              :       REAL(KIND=dp)                                      :: ifac, jfac, kfac, lfac
     722              :       REAL(KIND=dp), DIMENSION(2025)                     :: tmp2el
     723              :       REAL(KIND=dp), DIMENSION(45)                       :: tmp1el
     724              : 
     725              : ! One-electron integral
     726              : 
     727        13358 :       IF (PRESENT(int1el)) THEN
     728         8726 :          fdim = sepi%atm_int_size
     729         8726 :          ndim = 0
     730       108031 :          DO i = 1, fdim
     731       108031 :             tmp1el(i) = 0.0_dp
     732              :          END DO
     733        44845 :          DO i = 1, sepi%natorb
     734       144150 :          DO j = 1, i
     735        99305 :             ndim = ndim + 1
     736              : 
     737              :             ! Get the integral in the original frame (along z)
     738       334034 :             DO ind = 1, 2
     739       198610 :                imap = map_x_to_z(ind, i)
     740       198610 :                IF (imap == 0) CYCLE
     741       104345 :                ifac = fac_x_to_z(ind, i)
     742       412340 :                DO jnd = 1, 2
     743       208690 :                   jmap = map_x_to_z(jnd, j)
     744       208690 :                   IF (jmap == 0) CYCLE
     745       108965 :                   jfac = fac_x_to_z(jnd, j)
     746       108965 :                   gind = indexb(imap, jmap)
     747              : 
     748       407300 :                   tmp1el(ndim) = tmp1el(ndim) + ifac*jfac*int1el(gind)
     749              :                END DO
     750              :             END DO
     751              :          END DO
     752              :          END DO
     753       108031 :          DO i = 1, fdim
     754       108031 :             int1el(i) = tmp1el(i)
     755              :          END DO
     756              :       END IF
     757              : 
     758              :       ! Two electron integrals
     759        13358 :       IF (PRESENT(int2el)) THEN
     760         4632 :          sdim = sepi%atm_int_size
     761         4632 :          tdim = sepj%atm_int_size
     762         4632 :          fdim = sdim*tdim
     763         4632 :          ndim = 0
     764       976572 :          DO i = 1, fdim
     765       976572 :             tmp2el(i) = 0.0_dp
     766              :          END DO
     767        24516 :          DO i = 1, sepi%natorb
     768        80424 :          DO j = 1, i
     769       353124 :             DO k = 1, sepj%natorb
     770      1305180 :             DO l = 1, k
     771       971940 :                ndim = ndim + 1
     772              : 
     773              :                ! Get the integral in the original frame (along z)
     774      3193152 :                DO ind = 1, 2
     775      1943880 :                   imap = map_x_to_z(ind, i)
     776      1943880 :                   IF (imap == 0) CYCLE
     777      1120980 :                   ifac = fac_x_to_z(ind, i)
     778      4334880 :                   DO jnd = 1, 2
     779      2241960 :                      jmap = map_x_to_z(jnd, j)
     780      2241960 :                      IF (jmap == 0) CYCLE
     781      1257600 :                      jfac = fac_x_to_z(jnd, j)
     782      1257600 :                      gind = indexb(imap, jmap)
     783              : 
     784              :                      ! Get the integral in the original frame (along z)
     785      5716680 :                      DO knd = 1, 2
     786      2515200 :                         kmap = map_x_to_z(knd, k)
     787      2515200 :                         IF (kmap == 0) CYCLE
     788      1484832 :                         kfac = fac_x_to_z(knd, k)
     789      6696456 :                         DO lnd = 1, 2
     790      2969664 :                            lmap = map_x_to_z(lnd, l)
     791      2969664 :                            IF (lmap == 0) CYCLE
     792      1693128 :                            lfac = fac_x_to_z(lnd, l)
     793      1693128 :                            gknd = indexb(kmap, lmap)
     794              : 
     795      1693128 :                            tind = (gind - 1)*tdim + gknd
     796      5484864 :                            tmp2el(ndim) = tmp2el(ndim) + ifac*jfac*lfac*kfac*int2el(tind)
     797              :                         END DO
     798              :                      END DO
     799              : 
     800              :                   END DO
     801              :                END DO
     802              : 
     803              :             END DO
     804              :             END DO
     805              :          END DO
     806              :          END DO
     807       976572 :          DO i = 1, fdim
     808       976572 :             int2el(i) = tmp2el(i)
     809              :          END DO
     810              :       END IF
     811        13358 :    END SUBROUTINE invert_integral
     812              : 
     813              : ! **************************************************************************************************
     814              : !> \brief Exploits inversion symmetry to avoid divergence
     815              : !> \param sepi ...
     816              : !> \param sepj ...
     817              : !> \param dint1el ...
     818              : !> \param dint2el ...
     819              : !> \par History
     820              : !>      04.2007 created [tlaino]
     821              : !> \author Teodoro Laino - Zurich University
     822              : ! **************************************************************************************************
     823         2700 :    SUBROUTINE invert_derivative(sepi, sepj, dint1el, dint2el)
     824              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
     825              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
     826              :          OPTIONAL                                        :: dint1el, dint2el
     827              : 
     828              :       INTEGER                                            :: i, m
     829              :       REAL(KIND=dp)                                      :: tmp
     830              : 
     831              : ! Integral part
     832              : 
     833        10800 :       DO i = 1, 3
     834         8100 :          IF (PRESENT(dint1el)) THEN
     835         5259 :             CALL invert_integral(sepi, sepj, int1el=dint1el(i, :))
     836              :          END IF
     837        10800 :          IF (PRESENT(dint2el)) THEN
     838         2841 :             CALL invert_integral(sepi, sepj, int2el=dint2el(i, :))
     839              :          END IF
     840              :       END DO
     841              : 
     842              :       ! Derivatives part
     843         2700 :       IF (PRESENT(dint1el)) THEN
     844        80638 :          DO m = 1, SIZE(dint1el, 2)
     845        78885 :             tmp = dint1el(3, m)
     846        78885 :             dint1el(3, m) = dint1el(1, m)
     847        80638 :             dint1el(1, m) = tmp
     848              :          END DO
     849              :       END IF
     850         2700 :       IF (PRESENT(dint2el)) THEN
     851      1918622 :          DO m = 1, SIZE(dint2el, 2)
     852      1917675 :             tmp = dint2el(3, m)
     853      1917675 :             dint2el(3, m) = dint2el(1, m)
     854      1918622 :             dint2el(1, m) = tmp
     855              :          END DO
     856              :       END IF
     857         2700 :    END SUBROUTINE invert_derivative
     858              : 
     859              : ! **************************************************************************************************
     860              : !> \brief Calculates the ssss integral and analytical derivatives (main driver)
     861              : !> \param sepi parameters of atom i
     862              : !> \param sepj parameters of atom j
     863              : !> \param rij interatomic distance
     864              : !> \param ssss ...
     865              : !> \param dssss derivative of (ssss) integral
     866              : !>                          derivatives are intended w.r.t. rij
     867              : !> \param itype ...
     868              : !> \param se_taper ...
     869              : !> \param se_int_control input parameters that control the calculation of SE
     870              : !>                          integrals (shortrange, R3 residual, screening type)
     871              : !> \param lgrad ...
     872              : !> \par History
     873              : !>      03.2007 created [tlaino]
     874              : !>      Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
     875              : !>                 for computing integrals
     876              : !> \author Teodoro Laino - Zurich University
     877              : !> \note
     878              : !>      Analytical version - Analytical evaluation of gradients
     879              : !>      Teodoro Laino - Zurich University 04.2007
     880              : !>
     881              : ! **************************************************************************************************
     882     26676905 :    SUBROUTINE dssss_nucint_ana(sepi, sepj, rij, ssss, dssss, itype, se_taper, se_int_control, &
     883              :                                lgrad)
     884              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
     885              :       REAL(dp), INTENT(IN)                               :: rij
     886              :       REAL(dp), INTENT(OUT)                              :: ssss, dssss
     887              :       INTEGER, INTENT(IN)                                :: itype
     888              :       TYPE(se_taper_type), POINTER                       :: se_taper
     889              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
     890              :       LOGICAL, INTENT(IN)                                :: lgrad
     891              : 
     892              :       REAL(KIND=dp)                                      :: dft, ft
     893              :       TYPE(se_int_screen_type)                           :: se_int_screen
     894              : 
     895              : ! Compute the Tapering function
     896              : 
     897     26676905 :       ft = 1.0_dp
     898     26676905 :       dft = 0.0_dp
     899     26676905 :       IF (itype /= do_method_pchg) THEN
     900     11333008 :          ft = taper_eval(se_taper%taper, rij)
     901     11333008 :          dft = dtaper_eval(se_taper%taper, rij)
     902              :       END IF
     903              :       ! Evaluate additional taper function for dumped integrals
     904     26676905 :       IF (se_int_control%integral_screening == do_se_IS_kdso_d) THEN
     905      3629867 :          se_int_screen%ft = 1.0_dp
     906      3629867 :          se_int_screen%dft = 0.0_dp
     907      3629867 :          IF (itype /= do_method_pchg) THEN
     908      3629867 :             se_int_screen%ft = taper_eval(se_taper%taper_add, rij)
     909      3629867 :             se_int_screen%dft = dtaper_eval(se_taper%taper_add, rij)
     910              :          END IF
     911              :       END IF
     912              : 
     913              :       ! Value of the integrals for sp shell
     914              :       CALL nucint_sp_num(sepi, sepj, rij, ssss=ssss, itype=itype, se_int_control=se_int_control, &
     915     26676905 :                          se_int_screen=se_int_screen)
     916              : 
     917     26676905 :       IF (lgrad) THEN
     918              :          ! Integrals derivatives for sp shell
     919              :          CALL dnucint_sp_ana(sepi, sepj, rij, dssss=dssss, itype=itype, se_int_control=se_int_control, &
     920      8139105 :                              se_int_screen=se_int_screen)
     921              :       END IF
     922              : 
     923              :       ! Tapering the value of the integrals
     924     26676905 :       IF (lgrad) THEN
     925      8139105 :          dssss = ft*dssss + dft*ssss
     926              :       END IF
     927     26676905 :       ssss = ft*ssss
     928              : 
     929              :       ! Debug Procedure.. Check valifity of analytical gradients of nucint
     930              :       IF (debug_this_module .AND. lgrad) THEN
     931              :          CALL check_dssss_nucint_ana(sepi, sepj, rij, dssss, itype, se_int_control, se_taper=se_taper)
     932              :       END IF
     933     26676905 :    END SUBROUTINE dssss_nucint_ana
     934              : 
     935              : ! **************************************************************************************************
     936              : !> \brief Calculates the nuclear attraction integrals and analytical integrals (main driver)
     937              : !> \param sepi parameters of atom i
     938              : !> \param sepj parameters of atom j
     939              : !> \param rij interatomic distance
     940              : !> \param core ...
     941              : !> \param dcore derivative of 4 X 2 array of electron-core attraction integrals
     942              : !>                          derivatives are intended w.r.t. rij
     943              : !>         The storage of the nuclear attraction integrals  core(kl/ij) iS
     944              : !>         (SS/)=1,   (SO/)=2,   (OO/)=3,   (PP/)=4
     945              : !>         where ij=1 if the orbitals centred on atom i,  =2 if on atom j.
     946              : !> \param itype type of semi_empirical model
     947              : !>                          extension to the original routine to compute qm/mm
     948              : !>                          integrals
     949              : !> \param se_taper ...
     950              : !> \param se_int_control input parameters that control the calculation of SE
     951              : !>                          integrals (shortrange, R3 residual, screening type)
     952              : !> \param lgrad ...
     953              : !> \par History
     954              : !>      03.2007 created [tlaino]
     955              : !>      Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
     956              : !>                 for computing integrals
     957              : !> \author Teodoro Laino - Zurich University
     958              : !> \note
     959              : !>      Analytical version - Analytical evaluation of gradients
     960              : !>      Teodoro Laino - Zurich University 04.2007
     961              : !>
     962              : ! **************************************************************************************************
     963     16142300 :    SUBROUTINE dcore_nucint_ana(sepi, sepj, rij, core, dcore, itype, se_taper, &
     964              :                                se_int_control, lgrad)
     965              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
     966              :       REAL(dp), INTENT(IN)                               :: rij
     967              :       REAL(dp), DIMENSION(10, 2), INTENT(OUT)            :: core, dcore
     968              :       INTEGER, INTENT(IN)                                :: itype
     969              :       TYPE(se_taper_type), POINTER                       :: se_taper
     970              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
     971              :       LOGICAL, INTENT(IN)                                :: lgrad
     972              : 
     973              :       INTEGER                                            :: i
     974              :       REAL(KIND=dp)                                      :: dft, ft
     975              :       TYPE(se_int_screen_type)                           :: se_int_screen
     976              : 
     977              : ! Compute the Tapering function
     978              : 
     979     16142300 :       ft = 1.0_dp
     980     16142300 :       dft = 0.0_dp
     981     16142300 :       IF (itype /= do_method_pchg) THEN
     982       798403 :          ft = taper_eval(se_taper%taper, rij)
     983       798403 :          dft = dtaper_eval(se_taper%taper, rij)
     984              :       END IF
     985              :       ! Evaluate additional taper function for dumped integrals
     986     16142300 :       IF (se_int_control%integral_screening == do_se_IS_kdso_d) THEN
     987       260198 :          se_int_screen%ft = 1.0_dp
     988       260198 :          se_int_screen%dft = 0.0_dp
     989       260198 :          IF (itype /= do_method_pchg) THEN
     990       260198 :             se_int_screen%ft = taper_eval(se_taper%taper_add, rij)
     991       260198 :             se_int_screen%dft = dtaper_eval(se_taper%taper_add, rij)
     992              :          END IF
     993              :       END IF
     994              : 
     995              :       ! Value of the integrals for sp shell
     996              :       CALL nucint_sp_num(sepi, sepj, rij, core=core, itype=itype, &
     997     16142300 :                          se_int_control=se_int_control, se_int_screen=se_int_screen)
     998              : 
     999     16142300 :       IF (sepi%dorb .OR. sepj%dorb) THEN
    1000              :          ! Compute the contribution from d-orbitals
    1001              :          CALL nucint_d_num(sepi, sepj, rij, core, itype, &
    1002        37922 :                            se_int_control=se_int_control, se_int_screen=se_int_screen)
    1003              :       END IF
    1004              : 
    1005     16142300 :       IF (lgrad) THEN
    1006              :          ! Integrals derivatives for sp shell
    1007              :          CALL dnucint_sp_ana(sepi, sepj, rij, dcore=dcore, itype=itype, &
    1008      8014152 :                              se_int_control=se_int_control, se_int_screen=se_int_screen)
    1009              : 
    1010      8014152 :          IF (sepi%dorb .OR. sepj%dorb) THEN
    1011              :             ! Integral derivatives involving d-orbitals
    1012              :             CALL dnucint_d_ana(sepi, sepj, rij, dcore=dcore, itype=itype, &
    1013        17442 :                                se_int_control=se_int_control, se_int_screen=se_int_screen)
    1014              :          END IF
    1015              :       END IF
    1016              : 
    1017              :       ! Tapering the value of the integrals
    1018     16142300 :       IF (lgrad) THEN
    1019     26875194 :          DO i = 1, sepi%core_size
    1020     26875194 :             dcore(i, 1) = ft*dcore(i, 1) + dft*core(i, 1)
    1021              :          END DO
    1022      8957533 :          DO i = 1, sepj%core_size
    1023      8957533 :             dcore(i, 2) = ft*dcore(i, 2) + dft*core(i, 2)
    1024              :          END DO
    1025              :       END IF
    1026     54133615 :       DO i = 1, sepi%core_size
    1027     54133615 :          core(i, 1) = ft*core(i, 1)
    1028              :       END DO
    1029     18297519 :       DO i = 1, sepj%core_size
    1030     18297519 :          core(i, 2) = ft*core(i, 2)
    1031              :       END DO
    1032              : 
    1033              :       ! Debug Procedure.. Check valifity of analytical gradients of nucint
    1034              :       IF (debug_this_module .AND. lgrad) THEN
    1035              :          CALL check_dcore_nucint_ana(sepi, sepj, rij, dcore, itype, se_int_control, se_taper=se_taper)
    1036              :       END IF
    1037     16142300 :    END SUBROUTINE dcore_nucint_ana
    1038              : 
    1039              : ! **************************************************************************************************
    1040              : !> \brief Calculates the nuclear attraction integrals and derivatives for sp basis
    1041              : !> \param sepi parameters of atom i
    1042              : !> \param sepj parameters of atom j
    1043              : !> \param rij interatomic distance
    1044              : !> \param dssss derivative of (ssss) integral
    1045              : !>                          derivatives are intended w.r.t. rij
    1046              : !>         where ij=1 if the orbitals centred on atom i,  =2 if on atom j.
    1047              : !> \param dcore derivative of 4 X 2 array of electron-core attraction integrals
    1048              : !>         The storage of the nuclear attraction integrals  core(kl/ij) iS
    1049              : !>         (SS/)=1,   (SP/)=2,   (PP/)=3,   (P+P+/)=4,   (SD/)=5,
    1050              : !>         (DP/)=6,   (DD/)=7,   (D+P+)=8,  (D+D+/)=9,   (D#D#)=10
    1051              : !> \param itype type of semi_empirical model
    1052              : !>                          extension to the original routine to compute qm/mm
    1053              : !>                          integrals
    1054              : !> \param se_int_control input parameters that control the calculation of SE
    1055              : !>                          integrals (shortrange, R3 residual, screening type)
    1056              : !> \param se_int_screen ...
    1057              : !> \par History
    1058              : !>      04.2007 created [tlaino]
    1059              : !>      Teodoro Laino (03.2008) [tlaino] - University of Zurich : new driver
    1060              : !>                 for computing integrals
    1061              : !>      05.2008 Teodoro Laino [tlaino] - University of Zurich: major rewriting
    1062              : !> \author Teodoro Laino - Zurich University
    1063              : !> \note
    1064              : !>      Analytical version - Analytical evaluation of gradients
    1065              : !>      Teodoro Laino - Zurich University 04.2007
    1066              : !>      routine adapted from mopac7 (repp)
    1067              : !>      vector version written by Ernest R. Davidson, Indiana University
    1068              : ! **************************************************************************************************
    1069     16153257 :    SUBROUTINE dnucint_sp_ana(sepi, sepj, rij, dssss, dcore, itype, se_int_control, &
    1070              :                              se_int_screen)
    1071              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
    1072              :       REAL(dp), INTENT(IN)                               :: rij
    1073              :       REAL(dp), INTENT(INOUT), OPTIONAL                  :: dssss
    1074              :       REAL(dp), DIMENSION(10, 2), INTENT(INOUT), &
    1075              :          OPTIONAL                                        :: dcore
    1076              :       INTEGER, INTENT(IN)                                :: itype
    1077              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
    1078              :       TYPE(se_int_screen_type), INTENT(IN)               :: se_int_screen
    1079              : 
    1080              :       INTEGER                                            :: ij, kl
    1081              :       LOGICAL                                            :: l_core, l_ssss, si, sj
    1082              : 
    1083     16153257 :       l_core = PRESENT(dcore)
    1084     16153257 :       l_ssss = PRESENT(dssss)
    1085     16153257 :       IF (.NOT. (l_core .OR. l_ssss)) RETURN
    1086              : 
    1087     16153257 :       si = (sepi%natorb > 1)
    1088     16153257 :       sj = (sepj%natorb > 1)
    1089              : 
    1090     16153257 :       ij = indexa(1, 1)
    1091     16153257 :       IF (l_ssss) THEN
    1092              :          ! Store the value for the derivative of <S  S  | S  S  > (Used for computing the core-core interactions)
    1093      8139105 :          dssss = d_ijkl_sp(sepi, sepj, ij, ij, 0, 0, 0, 0, -1, rij, se_int_control, se_int_screen, itype)
    1094              :       END IF
    1095              : 
    1096     16153257 :       IF (l_core) THEN
    1097              :          !     <S  S  | S  S  >
    1098      8014152 :          kl = indexa(1, 1)
    1099      8014152 :          dcore(1, 1) = d_ijkl_sp(sepi, sepj, kl, ij, 0, 0, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
    1100      8014152 :          IF (si) THEN
    1101              :             !  <S  P  | S  S  >
    1102      3595276 :             kl = indexa(2, 1)
    1103      3595276 :             dcore(2, 1) = d_ijkl_sp(sepi, sepj, kl, ij, 0, 1, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
    1104              :             !  <P  P  | S  S  >
    1105      3595276 :             kl = indexa(2, 2)
    1106      3595276 :             dcore(3, 1) = d_ijkl_sp(sepi, sepj, kl, ij, 1, 1, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
    1107              :             !  <P+ P+ | S  S  >
    1108      3595276 :             kl = indexa(3, 3)
    1109      3595276 :             dcore(4, 1) = d_ijkl_sp(sepi, sepj, kl, ij, 1, 1, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
    1110              :          END IF
    1111              : 
    1112              :          !     <S  S  | S  S  >
    1113      8014152 :          kl = indexa(1, 1)
    1114      8014152 :          dcore(1, 2) = d_ijkl_sp(sepi, sepj, ij, kl, 0, 0, 0, 0, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
    1115      8014152 :          IF (sj) THEN
    1116              :             !  <S  S  | S  P  >
    1117       180523 :             kl = indexa(2, 1)
    1118       180523 :             dcore(2, 2) = d_ijkl_sp(sepi, sepj, ij, kl, 0, 0, 0, 1, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
    1119              :             !  <S  S  | P  P  >
    1120       180523 :             kl = indexa(2, 2)
    1121       180523 :             dcore(3, 2) = d_ijkl_sp(sepi, sepj, ij, kl, 0, 0, 1, 1, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
    1122              :             !  <S  S  | P+ P+ >
    1123       180523 :             kl = indexa(3, 3)
    1124       180523 :             dcore(4, 2) = d_ijkl_sp(sepi, sepj, ij, kl, 0, 0, 1, 1, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
    1125              :          END IF
    1126              :       END IF
    1127              :    END SUBROUTINE dnucint_sp_ana
    1128              : 
    1129              : ! **************************************************************************************************
    1130              : !> \brief Calculates the analytical derivative of the nuclear attraction
    1131              : !>        integrals involving d orbitals
    1132              : !> \param sepi parameters of atom i
    1133              : !> \param sepj parameters of atom j
    1134              : !> \param rij interatomic distance
    1135              : !> \param dcore 4 X 2 array of electron-core attraction integrals
    1136              : !>         The storage of the nuclear attraction integrals  core(kl/ij) iS
    1137              : !>         (SS/)=1,   (SP/)=2,   (PP/)=3,   (P+P+/)=4,   (SD/)=5,
    1138              : !>         (DP/)=6,   (DD/)=7,   (D+P+)=8,  (D+D+/)=9,   (D#D#)=10
    1139              : !>
    1140              : !>         where ij=1 if the orbitals centred on atom i,  =2 if on atom j.
    1141              : !> \param itype type of semi_empirical model
    1142              : !>                         extension to the original routine to compute qm/mm
    1143              : !>                         integrals
    1144              : !> \param se_int_control input parameters that control the calculation of SE
    1145              : !>                         integrals (shortrange, R3 residual, screening type)
    1146              : !> \param se_int_screen ...
    1147              : !> \author
    1148              : !>      Teodoro Laino (05.2008) [tlaino] - University of Zurich: created
    1149              : ! **************************************************************************************************
    1150        17442 :    SUBROUTINE dnucint_d_ana(sepi, sepj, rij, dcore, itype, se_int_control, &
    1151              :                             se_int_screen)
    1152              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
    1153              :       REAL(dp), INTENT(IN)                               :: rij
    1154              :       REAL(dp), DIMENSION(10, 2), INTENT(INOUT)          :: dcore
    1155              :       INTEGER, INTENT(IN)                                :: itype
    1156              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
    1157              :       TYPE(se_int_screen_type), INTENT(IN)               :: se_int_screen
    1158              : 
    1159              :       INTEGER                                            :: ij, kl
    1160              : 
    1161              : ! Check if d-orbitals are present
    1162              : 
    1163        17442 :       IF (sepi%dorb .OR. sepj%dorb) THEN
    1164        17442 :          ij = indexa(1, 1)
    1165        17442 :          IF (sepj%dorb) THEN
    1166              :             !  <S S | D S>
    1167         9880 :             kl = indexa(5, 1)
    1168         9880 :             dcore(5, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 0, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
    1169              :             !  <S S | D P >
    1170         9880 :             kl = indexa(5, 2)
    1171         9880 :             dcore(6, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 1, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
    1172              :             !  <S S | D D >
    1173         9880 :             kl = indexa(5, 5)
    1174         9880 :             dcore(7, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
    1175              :             !  <S S | D+P+>
    1176         9880 :             kl = indexa(6, 3)
    1177         9880 :             dcore(8, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 1, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
    1178              :             !  <S S | D+D+>
    1179         9880 :             kl = indexa(6, 6)
    1180         9880 :             dcore(9, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
    1181              :             !  <S S | D#D#>
    1182         9880 :             kl = indexa(8, 8)
    1183         9880 :             dcore(10, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
    1184              :          END IF
    1185        17442 :          IF (sepi%dorb) THEN
    1186              :             !  <D S | S S>
    1187        10177 :             kl = indexa(5, 1)
    1188        10177 :             dcore(5, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 0, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
    1189              :             !  <D P | S S >
    1190        10177 :             kl = indexa(5, 2)
    1191        10177 :             dcore(6, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 1, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
    1192              :             !  <D D | S S >
    1193        10177 :             kl = indexa(5, 5)
    1194        10177 :             dcore(7, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
    1195              :             !  <D+P+| S S >
    1196        10177 :             kl = indexa(6, 3)
    1197        10177 :             dcore(8, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 1, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
    1198              :             !  <D+D+| S S >
    1199        10177 :             kl = indexa(6, 6)
    1200        10177 :             dcore(9, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
    1201              :             !  <D#D#| S S >
    1202        10177 :             kl = indexa(8, 8)
    1203        10177 :             dcore(10, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
    1204              :          END IF
    1205              :       END IF
    1206        17442 :    END SUBROUTINE dnucint_d_ana
    1207              : 
    1208              : ! **************************************************************************************************
    1209              : !> \brief calculates the derivative of the two-particle interactions
    1210              : !> \param sepi Atomic parameters of first atom
    1211              : !> \param sepj Atomic parameters of second atom
    1212              : !> \param rijv Coordinate vector i -> j
    1213              : !> \param w Array of two-electron repulsion integrals.
    1214              : !> \param dw ...
    1215              : !> \param se_int_control ...
    1216              : !> \param se_taper ...
    1217              : !> \par History
    1218              : !>      04.2007 created [tlaino]
    1219              : !>      Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
    1220              : !>                 for computing integrals
    1221              : !> \author Teodoro Laino - Zurich University
    1222              : !> \note
    1223              : !>      Analytical version - Analytical evaluation of gradients
    1224              : !>      Teodoro Laino - Zurich University 04.2007
    1225              : !>      routine adapted from mopac7 (repp)
    1226              : !>      vector version written by Ernest R. Davidson, Indiana University
    1227              : ! **************************************************************************************************
    1228      1144109 :    RECURSIVE SUBROUTINE rotint_ana(sepi, sepj, rijv, w, dw, se_int_control, se_taper)
    1229              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
    1230              :       REAL(dp), DIMENSION(3), INTENT(IN)                 :: rijv
    1231              :       REAL(dp), DIMENSION(2025), INTENT(OUT), OPTIONAL   :: w
    1232              :       REAL(dp), DIMENSION(3, 2025), INTENT(OUT), &
    1233              :          OPTIONAL                                        :: dw
    1234              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
    1235              :       TYPE(se_taper_type), POINTER                       :: se_taper
    1236              : 
    1237              :       INTEGER                                            :: i, i1, ii, ij, ij1, iminus, istep, &
    1238              :                                                             iw_loc, j, j1, jj, k, kk, kl, l, &
    1239              :                                                             limij, limkl, mm
    1240              :       LOGICAL                                            :: invert, l_w, lgrad
    1241              :       LOGICAL, DIMENSION(45, 45)                         :: logv, logv_d
    1242              :       REAL(dp)                                           :: rij, xtmp
    1243              :       REAL(dp), DIMENSION(3)                             :: drij
    1244              :       REAL(KIND=dp)                                      :: cc, cc_d(3), wrepp, wrepp_d(3)
    1245              :       REAL(KIND=dp), DIMENSION(2025)                     :: ww
    1246              :       REAL(KIND=dp), DIMENSION(3, 2025)                  :: ww_d
    1247              :       REAL(KIND=dp), DIMENSION(3, 45, 45)                :: v_d
    1248              :       REAL(KIND=dp), DIMENSION(45, 45)                   :: v
    1249              :       REAL(KIND=dp), DIMENSION(491)                      :: rep, rep_d
    1250              :       TYPE(rotmat_type), POINTER                         :: ij_matrix
    1251              : 
    1252      1144109 :       NULLIFY (ij_matrix)
    1253      1144109 :       l_w = PRESENT(w)
    1254      1144109 :       lgrad = PRESENT(dw)
    1255      1144109 :       IF (.NOT. (l_w .OR. lgrad)) RETURN
    1256              : 
    1257      4576436 :       rij = DOT_PRODUCT(rijv, rijv)
    1258      1144109 :       IF (rij > rij_threshold) THEN
    1259              :          ! The repulsion integrals over molecular frame (w) are stored in the
    1260              :          ! order in which they will later be used.  ie.  (i,j/k,l) where
    1261              :          ! j.le.i  and  l.le.k     and l varies most rapidly and i least
    1262              :          ! rapidly.  (anti-normal computer storage)
    1263      1144109 :          rij = SQRT(rij)
    1264              : 
    1265              :          ! Create the rotation matrix
    1266      1144109 :          CALL rotmat_create(ij_matrix)
    1267      1144109 :          CALL rotmat(sepi, sepj, rijv, rij, ij_matrix, do_derivatives=lgrad, do_invert=invert)
    1268              : 
    1269              :          ! Compute integrals in diatomic frame as well their derivatives (if requested)
    1270      1144109 :          CALL dterep_ana(sepi, sepj, rij, rep, rep_d, se_taper, se_int_control, lgrad=lgrad)
    1271              : 
    1272      1144109 :          IF (lgrad) THEN
    1273       511306 :             drij(1) = rijv(1)/rij
    1274       511306 :             drij(2) = rijv(2)/rij
    1275       511306 :             drij(3) = rijv(3)/rij
    1276              :             ! Possibly Invert Frame
    1277       511306 :             IF (invert) THEN
    1278          947 :                xtmp = drij(3)
    1279          947 :                drij(3) = drij(1)
    1280          947 :                drij(1) = xtmp
    1281              :             END IF
    1282              :          END IF
    1283              : 
    1284      1144109 :          ii = sepi%natorb
    1285      1144109 :          kk = sepj%natorb
    1286              :          ! First step in rotation of integrals
    1287              :          CALL rot_2el_2c_first(sepi, sepj, rijv, se_int_control, se_taper, invert, ii, kk, rep, logv, ij_matrix, &
    1288      1144109 :                                v, lgrad, rep_d, v_d, logv_d, drij)
    1289              : 
    1290              :          ! Integrals if requested
    1291      1144109 :          IF (l_w) THEN
    1292              :             ! Rotate Integrals
    1293       632803 :             IF (ii*kk > 0) THEN
    1294       632803 :                limij = sepi%atm_int_size
    1295       632803 :                limkl = sepj%atm_int_size
    1296       632803 :                istep = limkl*limij
    1297     35148717 :                DO i1 = 1, istep
    1298     35148717 :                   ww(i1) = 0.0_dp
    1299              :                END DO
    1300              :                ! Second step in rotation of integrals
    1301      2456506 :                DO i1 = 1, ii
    1302      7017769 :                   DO j1 = 1, i1
    1303      4561263 :                      ij = indexa(i1, j1)
    1304      4561263 :                      jj = indexb(i1, j1)
    1305      4561263 :                      mm = int2c_type(jj)
    1306     19448726 :                      DO k = 1, kk
    1307     52140937 :                         DO l = 1, k
    1308     34515914 :                            kl = indexb(k, l)
    1309     47579674 :                            IF (logv(ij, kl)) THEN
    1310     21412441 :                               wrepp = v(ij, kl)
    1311      3258722 :                               SELECT CASE (mm)
    1312              :                               CASE (1)
    1313              :                                  ! (SS/)
    1314      3258722 :                                  i = 1
    1315      3258722 :                                  j = 1
    1316      3258722 :                                  iw_loc = (indexb(i, j) - 1)*limkl + kl
    1317      3258722 :                                  ww(iw_loc) = wrepp
    1318              :                               CASE (2)
    1319              :                                  ! (SP/)
    1320     20421512 :                                  j = 1
    1321     20421512 :                                  DO i = 1, 3
    1322     15316134 :                                     iw_loc = (indexb(i + 1, j) - 1)*limkl + kl
    1323     20421512 :                                     ww(iw_loc) = ww(iw_loc) + ij_matrix%sp(i1 - 1, i)*wrepp
    1324              :                                  END DO
    1325              :                               CASE (3)
    1326              :                                  ! (PP/)
    1327     40688688 :                                  DO i = 1, 3
    1328     30516516 :                                     cc = ij_matrix%pp(i, i1 - 1, j1 - 1)
    1329     30516516 :                                     iw_loc = (indexb(i + 1, i + 1) - 1)*limkl + kl
    1330     30516516 :                                     ww(iw_loc) = ww(iw_loc) + cc*wrepp
    1331     30516516 :                                     iminus = i - 1
    1332     40688688 :                                     IF (iminus /= 0) THEN
    1333     50860860 :                                        DO j = 1, iminus
    1334     30516516 :                                           cc = ij_matrix%pp(1 + i + j, i1 - 1, j1 - 1)
    1335     30516516 :                                           iw_loc = (indexb(i + 1, j + 1) - 1)*limkl + kl
    1336     50860860 :                                           ww(iw_loc) = ww(iw_loc) + cc*wrepp
    1337              :                                        END DO
    1338              :                                     END IF
    1339              :                                  END DO
    1340              :                               CASE (4)
    1341              :                                  ! (SD/)
    1342      2577900 :                                  j = 1
    1343      2577900 :                                  DO i = 1, 5
    1344      2148250 :                                     iw_loc = (indexb(i + 4, j) - 1)*limkl + kl
    1345      2577900 :                                     ww(iw_loc) = ww(iw_loc) + ij_matrix%sd(i1 - 4, i)*wrepp
    1346              :                                  END DO
    1347              :                               CASE (5)
    1348              :                                  ! (DP/)
    1349      6637350 :                                  DO i = 1, 5
    1350     23230725 :                                     DO j = 1, 3
    1351     16593375 :                                        iw_loc = (indexb(i + 4, j + 1) - 1)*limkl + kl
    1352     16593375 :                                        ij1 = 3*(i - 1) + j
    1353     22124500 :                                        ww(iw_loc) = ww(iw_loc) + ij_matrix%pd(ij1, i1 - 4, j1 - 1)*wrepp
    1354              :                                     END DO
    1355              :                                  END DO
    1356              :                               CASE (6)
    1357              :                                  ! (DD/)
    1358     28113911 :                                  DO i = 1, 5
    1359      6701470 :                                     cc = ij_matrix%dd(i, i1 - 4, j1 - 4)
    1360      6701470 :                                     iw_loc = (indexb(i + 4, i + 4) - 1)*limkl + kl
    1361      6701470 :                                     ww(iw_loc) = ww(iw_loc) + cc*wrepp
    1362      6701470 :                                     iminus = i - 1
    1363      8041764 :                                     IF (iminus /= 0) THEN
    1364     18764116 :                                        DO j = 1, iminus
    1365     13402940 :                                           ij1 = inddd(i, j)
    1366     13402940 :                                           cc = ij_matrix%dd(ij1, i1 - 4, j1 - 4)
    1367     13402940 :                                           iw_loc = (indexb(i + 4, j + 4) - 1)*limkl + kl
    1368     18764116 :                                           ww(iw_loc) = ww(iw_loc) + cc*wrepp
    1369              :                                        END DO
    1370              :                                     END IF
    1371              :                                  END DO
    1372              :                               END SELECT
    1373              :                            END IF
    1374              :                         END DO
    1375              :                      END DO
    1376              :                   END DO
    1377              :                END DO
    1378              :                ! Store two electron integrals in the triangular format
    1379       632803 :                CALL store_2el_2c_diag(limij, limkl, ww(1:istep), w)
    1380       632803 :                IF (invert) CALL invert_integral(sepi, sepj, int2el=w)
    1381              :             END IF
    1382              : 
    1383              :             IF (debug_this_module) THEN
    1384              :                ! Check value of integrals
    1385              :                CALL check_rotint_ana(sepi, sepj, rijv, w, se_int_control=se_int_control, se_taper=se_taper)
    1386              :             END IF
    1387              :          END IF
    1388              : 
    1389              :          ! Gradients if requested
    1390      1144109 :          IF (lgrad) THEN
    1391              :             ! Rotate Integrals derivatives
    1392       511306 :             IF (ii*kk > 0) THEN
    1393       511306 :                limij = sepi%atm_int_size
    1394       511306 :                limkl = sepj%atm_int_size
    1395       511306 :                istep = limkl*limij
    1396     30915998 :                DO i1 = 1, istep
    1397     30404692 :                   ww_d(1, i1) = 0.0_dp
    1398     30404692 :                   ww_d(2, i1) = 0.0_dp
    1399     30915998 :                   ww_d(3, i1) = 0.0_dp
    1400              :                END DO
    1401              : 
    1402              :                ! Second step in rotation of integrals
    1403      2018237 :                DO i1 = 1, ii
    1404      5821698 :                   DO j1 = 1, i1
    1405      3803461 :                      ij = indexa(i1, j1)
    1406      3803461 :                      jj = indexb(i1, j1)
    1407      3803461 :                      mm = int2c_type(jj)
    1408     16580730 :                      DO k = 1, kk
    1409     45478491 :                         DO l = 1, k
    1410     30404692 :                            kl = indexb(k, l)
    1411     41675030 :                            IF (logv_d(ij, kl)) THEN
    1412     19007015 :                               wrepp_d(1) = v_d(1, ij, kl)
    1413     19007015 :                               wrepp_d(2) = v_d(2, ij, kl)
    1414     19007015 :                               wrepp_d(3) = v_d(3, ij, kl)
    1415     19007015 :                               wrepp = v(ij, kl)
    1416      2769321 :                               SELECT CASE (mm)
    1417              :                               CASE (1)
    1418              :                                  ! (SS/)
    1419      2769321 :                                  i = 1
    1420      2769321 :                                  j = 1
    1421      2769321 :                                  iw_loc = (indexb(i, j) - 1)*limkl + kl
    1422      2769321 :                                  ww_d(1, iw_loc) = wrepp_d(1)
    1423      2769321 :                                  ww_d(2, iw_loc) = wrepp_d(2)
    1424      2769321 :                                  ww_d(3, iw_loc) = wrepp_d(3)
    1425              :                               CASE (2)
    1426              :                                  ! (SP/)
    1427     18274700 :                                  j = 1
    1428     18274700 :                                  DO i = 1, 3
    1429     13706025 :                                     iw_loc = (indexb(i + 1, j) - 1)*limkl + kl
    1430              :                                     ww_d(1, iw_loc) = ww_d(1, iw_loc) + ij_matrix%sp_d(1, i1 - 1, i)*wrepp + &
    1431     13706025 :                                                       ij_matrix%sp(i1 - 1, i)*wrepp_d(1)
    1432              : 
    1433              :                                     ww_d(2, iw_loc) = ww_d(2, iw_loc) + ij_matrix%sp_d(2, i1 - 1, i)*wrepp + &
    1434     13706025 :                                                       ij_matrix%sp(i1 - 1, i)*wrepp_d(2)
    1435              : 
    1436              :                                     ww_d(3, iw_loc) = ww_d(3, iw_loc) + ij_matrix%sp_d(3, i1 - 1, i)*wrepp + &
    1437     18274700 :                                                       ij_matrix%sp(i1 - 1, i)*wrepp_d(3)
    1438              :                                  END DO
    1439              :                               CASE (3)
    1440              :                                  ! (PP/)
    1441     36201524 :                                  DO i = 1, 3
    1442     27151143 :                                     cc = ij_matrix%pp(i, i1 - 1, j1 - 1)
    1443     27151143 :                                     cc_d(1) = ij_matrix%pp_d(1, i, i1 - 1, j1 - 1)
    1444     27151143 :                                     cc_d(2) = ij_matrix%pp_d(2, i, i1 - 1, j1 - 1)
    1445     27151143 :                                     cc_d(3) = ij_matrix%pp_d(3, i, i1 - 1, j1 - 1)
    1446     27151143 :                                     iw_loc = (indexb(i + 1, i + 1) - 1)*limkl + kl
    1447     27151143 :                                     ww_d(1, iw_loc) = ww_d(1, iw_loc) + cc_d(1)*wrepp + cc*wrepp_d(1)
    1448     27151143 :                                     ww_d(2, iw_loc) = ww_d(2, iw_loc) + cc_d(2)*wrepp + cc*wrepp_d(2)
    1449     27151143 :                                     ww_d(3, iw_loc) = ww_d(3, iw_loc) + cc_d(3)*wrepp + cc*wrepp_d(3)
    1450     27151143 :                                     iminus = i - 1
    1451     36201524 :                                     IF (iminus /= 0) THEN
    1452     45251905 :                                        DO j = 1, iminus
    1453     27151143 :                                           cc = ij_matrix%pp(1 + i + j, i1 - 1, j1 - 1)
    1454     27151143 :                                           cc_d(1) = ij_matrix%pp_d(1, 1 + i + j, i1 - 1, j1 - 1)
    1455     27151143 :                                           cc_d(2) = ij_matrix%pp_d(2, 1 + i + j, i1 - 1, j1 - 1)
    1456     27151143 :                                           cc_d(3) = ij_matrix%pp_d(3, 1 + i + j, i1 - 1, j1 - 1)
    1457     27151143 :                                           iw_loc = (indexb(i + 1, j + 1) - 1)*limkl + kl
    1458     27151143 :                                           ww_d(1, iw_loc) = ww_d(1, iw_loc) + cc_d(1)*wrepp + cc*wrepp_d(1)
    1459     27151143 :                                           ww_d(2, iw_loc) = ww_d(2, iw_loc) + cc_d(2)*wrepp + cc*wrepp_d(2)
    1460     45251905 :                                           ww_d(3, iw_loc) = ww_d(3, iw_loc) + cc_d(3)*wrepp + cc*wrepp_d(3)
    1461              :                                        END DO
    1462              :                                     END IF
    1463              :                                  END DO
    1464              :                               CASE (4)
    1465              :                                  ! (SD/)
    1466      2351052 :                                  j = 1
    1467      2351052 :                                  DO i = 1, 5
    1468      1959210 :                                     iw_loc = (indexb(i + 4, j) - 1)*limkl + kl
    1469              :                                     ww_d(1, iw_loc) = ww_d(1, iw_loc) + ij_matrix%sd_d(1, i1 - 4, i)*wrepp + &
    1470      1959210 :                                                       ij_matrix%sd(i1 - 4, i)*wrepp_d(1)
    1471              : 
    1472              :                                     ww_d(2, iw_loc) = ww_d(2, iw_loc) + ij_matrix%sd_d(2, i1 - 4, i)*wrepp + &
    1473      1959210 :                                                       ij_matrix%sd(i1 - 4, i)*wrepp_d(2)
    1474              : 
    1475              :                                     ww_d(3, iw_loc) = ww_d(3, iw_loc) + ij_matrix%sd_d(3, i1 - 4, i)*wrepp + &
    1476      2351052 :                                                       ij_matrix%sd(i1 - 4, i)*wrepp_d(3)
    1477              :                                  END DO
    1478              :                               CASE (5)
    1479              :                                  ! (DP/)
    1480      6077604 :                                  DO i = 1, 5
    1481     21271614 :                                     DO j = 1, 3
    1482     15194010 :                                        iw_loc = (indexb(i + 4, j + 1) - 1)*limkl + kl
    1483     15194010 :                                        ij1 = 3*(i - 1) + j
    1484              :                                        ww_d(1, iw_loc) = ww_d(1, iw_loc) + ij_matrix%pd_d(1, ij1, i1 - 4, j1 - 1)*wrepp + &
    1485     15194010 :                                                          ij_matrix%pd(ij1, i1 - 4, j1 - 1)*wrepp_d(1)
    1486              : 
    1487              :                                        ww_d(2, iw_loc) = ww_d(2, iw_loc) + ij_matrix%pd_d(2, ij1, i1 - 4, j1 - 1)*wrepp + &
    1488     15194010 :                                                          ij_matrix%pd(ij1, i1 - 4, j1 - 1)*wrepp_d(2)
    1489              : 
    1490              :                                        ww_d(3, iw_loc) = ww_d(3, iw_loc) + ij_matrix%pd_d(3, ij1, i1 - 4, j1 - 1)*wrepp + &
    1491     20258680 :                                                          ij_matrix%pd(ij1, i1 - 4, j1 - 1)*wrepp_d(3)
    1492              :                                     END DO
    1493              :                                  END DO
    1494              :                               CASE (6)
    1495              :                                  ! (DD/)
    1496     25076325 :                                  DO i = 1, 5
    1497      6069310 :                                     cc = ij_matrix%dd(i, i1 - 4, j1 - 4)
    1498     24277240 :                                     cc_d = ij_matrix%dd_d(:, i, i1 - 4, j1 - 4)
    1499      6069310 :                                     iw_loc = (indexb(i + 4, i + 4) - 1)*limkl + kl
    1500      6069310 :                                     ww_d(1, iw_loc) = ww_d(1, iw_loc) + cc_d(1)*wrepp + cc*wrepp_d(1)
    1501      6069310 :                                     ww_d(2, iw_loc) = ww_d(2, iw_loc) + cc_d(2)*wrepp + cc*wrepp_d(2)
    1502      6069310 :                                     ww_d(3, iw_loc) = ww_d(3, iw_loc) + cc_d(3)*wrepp + cc*wrepp_d(3)
    1503      6069310 :                                     iminus = i - 1
    1504      7283172 :                                     IF (iminus /= 0) THEN
    1505     16994068 :                                        DO j = 1, iminus
    1506     12138620 :                                           ij1 = inddd(i, j)
    1507     12138620 :                                           cc = ij_matrix%dd(ij1, i1 - 4, j1 - 4)
    1508     12138620 :                                           cc_d(1) = ij_matrix%dd_d(1, ij1, i1 - 4, j1 - 4)
    1509     12138620 :                                           cc_d(2) = ij_matrix%dd_d(2, ij1, i1 - 4, j1 - 4)
    1510     12138620 :                                           cc_d(3) = ij_matrix%dd_d(3, ij1, i1 - 4, j1 - 4)
    1511     12138620 :                                           iw_loc = (indexb(i + 4, j + 4) - 1)*limkl + kl
    1512     12138620 :                                           ww_d(1, iw_loc) = ww_d(1, iw_loc) + cc_d(1)*wrepp + cc*wrepp_d(1)
    1513     12138620 :                                           ww_d(2, iw_loc) = ww_d(2, iw_loc) + cc_d(2)*wrepp + cc*wrepp_d(2)
    1514     16994068 :                                           ww_d(3, iw_loc) = ww_d(3, iw_loc) + cc_d(3)*wrepp + cc*wrepp_d(3)
    1515              :                                        END DO
    1516              :                                     END IF
    1517              :                                  END DO
    1518              :                               END SELECT
    1519              :                            END IF
    1520              :                         END DO
    1521              :                      END DO
    1522              :                   END DO
    1523              :                END DO
    1524              :                ! Store two electron integrals in the triangular format
    1525              :                CALL store_2el_2c_diag(limij, limkl, ww_dx=ww_d(1, 1:istep), ww_dy=ww_d(2, 1:istep), ww_dz=ww_d(3, 1:istep), &
    1526     91725382 :                                       dw=dw)
    1527       511306 :                IF (invert) CALL invert_derivative(sepi, sepj, dint2el=dw)
    1528              :             END IF
    1529              : 
    1530              :             IF (debug_this_module) THEN
    1531              :                ! Check derivatives
    1532              :                CALL check_rotint_ana(sepi, sepj, rijv, dw=dw, se_int_control=se_int_control, se_taper=se_taper)
    1533              :             END IF
    1534              :          END IF
    1535      1144109 :          CALL rotmat_release(ij_matrix)
    1536              :       END IF
    1537              :    END SUBROUTINE rotint_ana
    1538              : 
    1539              : ! **************************************************************************************************
    1540              : !> \brief Calculates the derivative and the value of two-electron repulsion
    1541              : !>      integrals and the nuclear attraction integrals w.r.t. |r|
    1542              : !> \param sepi parameters of atom i
    1543              : !> \param sepj parameters of atom j
    1544              : !> \param rij interatomic distance
    1545              : !> \param rep rray of two-electron repulsion integrals
    1546              : !> \param rep_d array of two-electron repulsion integrals derivatives
    1547              : !> \param se_taper ...
    1548              : !> \param se_int_control input parameters that control the calculation of SE
    1549              : !>                         integrals (shortrange, R3 residual, screening type)
    1550              : !> \param lgrad ...
    1551              : !> \par History
    1552              : !>      03.2008 created [tlaino]
    1553              : !> \author Teodoro Laino [tlaino] - Zurich University
    1554              : ! **************************************************************************************************
    1555      1144109 :    RECURSIVE SUBROUTINE dterep_ana(sepi, sepj, rij, rep, rep_d, se_taper, &
    1556              :                                    se_int_control, lgrad)
    1557              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
    1558              :       REAL(KIND=dp), INTENT(IN)                          :: rij
    1559              :       REAL(KIND=dp), DIMENSION(491), INTENT(OUT)         :: rep, rep_d
    1560              :       TYPE(se_taper_type), POINTER                       :: se_taper
    1561              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
    1562              :       LOGICAL, INTENT(IN)                                :: lgrad
    1563              : 
    1564              :       INTEGER                                            :: i, ij, j, k, kl, l, lasti, lastj, li, &
    1565              :                                                             lj, lk, ll, numb
    1566              :       REAL(KIND=dp)                                      :: dft, ft, ft1
    1567              :       TYPE(se_int_screen_type)                           :: se_int_screen
    1568              : 
    1569              : ! Compute the tapering function and its derivatives
    1570              : 
    1571      1144109 :       ft = taper_eval(se_taper%taper, rij)
    1572      1144109 :       dft = 0.0_dp
    1573      1144109 :       ft1 = ft
    1574      1144109 :       IF (lgrad) THEN
    1575       511306 :          ft1 = 1.0_dp
    1576       511306 :          dft = dtaper_eval(se_taper%taper, rij)
    1577              :       END IF
    1578              :       ! Evaluate additional taper function for dumped integrals
    1579      1144109 :       IF (se_int_control%integral_screening == do_se_IS_kdso_d) THEN
    1580       260198 :          se_int_screen%ft = taper_eval(se_taper%taper_add, rij)
    1581       260198 :          IF (lgrad) &
    1582       130072 :             se_int_screen%dft = dtaper_eval(se_taper%taper_add, rij)
    1583              :       END IF
    1584              : 
    1585              :       ! Integral Values for sp shells only
    1586              :       CALL terep_sp_num(sepi, sepj, rij, rep, se_int_control=se_int_control, &
    1587      1144109 :                         se_int_screen=se_int_screen, ft=ft1)
    1588              : 
    1589      1144109 :       IF (sepi%dorb .OR. sepj%dorb) THEN
    1590              :          ! Compute the contribution from d-orbitals
    1591              :          CALL terep_d_num(sepi, sepj, rij, rep, se_int_control=se_int_control, &
    1592        54500 :                           se_int_screen=se_int_screen, ft=ft1)
    1593              :       END IF
    1594              : 
    1595      1144109 :       IF (lgrad) THEN
    1596              :          ! Integral Derivatives
    1597              :          CALL dterep_sp_ana(sepi, sepj, rij, rep_d, rep, se_int_control, &
    1598       511306 :                             se_int_screen, ft, dft)
    1599              : 
    1600       511306 :          IF (sepi%dorb .OR. sepj%dorb) THEN
    1601              :             ! Compute the derivatives from d-orbitals
    1602              :             CALL dterep_d_ana(sepi, sepj, rij, rep_d, rep, se_int_control, &
    1603        25368 :                               se_int_screen, ft, dft)
    1604              :          END IF
    1605              : 
    1606              :          ! Tapering Integral values
    1607       511306 :          lasti = sepi%natorb
    1608       511306 :          lastj = sepj%natorb
    1609      2018237 :          DO i = 1, lasti
    1610      1506931 :             li = l_index(i)
    1611      5821698 :             DO j = 1, i
    1612      3803461 :                lj = l_index(j)
    1613      3803461 :                ij = indexa(i, j)
    1614     16580730 :                DO k = 1, lastj
    1615     11270338 :                   lk = l_index(k)
    1616     45478491 :                   DO l = 1, k
    1617     30404692 :                      ll = l_index(l)
    1618     30404692 :                      kl = indexa(k, l)
    1619     30404692 :                      numb = ijkl_ind(ij, kl)
    1620     41675030 :                      IF (numb > 0) rep(numb) = rep(numb)*ft
    1621              :                   END DO
    1622              :                END DO
    1623              :             END DO
    1624              :          END DO
    1625              :       END IF
    1626              : 
    1627              :       ! Possibly debug 2el 2cent integrals and derivatives
    1628              :       IF (debug_this_module) THEN
    1629              :          CALL check_dterep_ana(sepi, sepj, rij, rep, rep_d, se_int_control, se_taper=se_taper, &
    1630              :                                lgrad=lgrad)
    1631              :       END IF
    1632      1144109 :    END SUBROUTINE dterep_ana
    1633              : 
    1634              : ! **************************************************************************************************
    1635              : !> \brief Calculates the derivative and the value of two-electron repulsion
    1636              : !>      integrals and the nuclear attraction integrals w.r.t. |r| - sp shells only
    1637              : !> \param sepi parameters of atom i
    1638              : !> \param sepj parameters of atom j
    1639              : !> \param rij interatomic distance
    1640              : !> \param drep array of derivatives of two-electron repulsion integrals
    1641              : !> \param rep array of two-electron repulsion integrals
    1642              : !> \param se_int_control input parameters that control the calculation of SE
    1643              : !>                         integrals (shortrange, R3 residual, screening type)
    1644              : !> \param se_int_screen ...
    1645              : !> \param ft ...
    1646              : !> \param dft ...
    1647              : !> \par History
    1648              : !>      04.2007 created [tlaino]
    1649              : !>      Teodoro Laino (03.2008) [tlaino] - University of Zurich : new driver
    1650              : !>                 for computing integrals
    1651              : !>      05.2008 Teodoro Laino [tlaino] - University of Zurich: major rewriting
    1652              : !> \author Teodoro Laino - Zurich University
    1653              : !> \note
    1654              : !>      Analytical version - Analytical evaluation of gradients
    1655              : !>      Teodoro Laino - Zurich University 04.2007
    1656              : !>      routine adapted from mopac7 (repp)
    1657              : !>      vector version written by Ernest R. Davidson, Indiana University
    1658              : ! **************************************************************************************************
    1659       511306 :    SUBROUTINE dterep_sp_ana(sepi, sepj, rij, drep, rep, se_int_control, &
    1660              :                             se_int_screen, ft, dft)
    1661              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
    1662              :       REAL(dp), INTENT(IN)                               :: rij
    1663              :       REAL(dp), DIMENSION(491), INTENT(OUT)              :: drep
    1664              :       REAL(dp), DIMENSION(491), INTENT(IN)               :: rep
    1665              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
    1666              :       TYPE(se_int_screen_type), INTENT(IN)               :: se_int_screen
    1667              :       REAL(dp), INTENT(IN)                               :: ft, dft
    1668              : 
    1669              :       INTEGER                                            :: i, ij, j, k, kl, l, lasti, lastj, li, &
    1670              :                                                             lj, lk, ll, nold, numb
    1671              :       REAL(KIND=dp)                                      :: tmp
    1672              : 
    1673       511306 :       lasti = sepi%natorb
    1674       511306 :       lastj = sepj%natorb
    1675      1941917 :       DO i = 1, MIN(lasti, 4)
    1676      1430611 :          li = l_index(i)
    1677      5211138 :          DO j = 1, i
    1678      3269221 :             lj = l_index(j)
    1679      3269221 :             ij = indexa(i, j)
    1680     13175340 :             DO k = 1, MIN(lastj, 4)
    1681      8475508 :                lk = l_index(k)
    1682     30632811 :                DO l = 1, k
    1683     18888082 :                   ll = l_index(l)
    1684     18888082 :                   kl = indexa(k, l)
    1685              : 
    1686     18888082 :                   numb = ijkl_ind(ij, kl)
    1687     27363590 :                   IF (numb > 0) THEN
    1688      6851959 :                      nold = ijkl_sym(numb)
    1689      6851959 :                      IF (nold > 0) THEN
    1690      2207937 :                         drep(numb) = drep(nold)
    1691      4644022 :                      ELSE IF (nold < 0) THEN
    1692            0 :                         drep(numb) = -drep(-nold)
    1693      4644022 :                      ELSE IF (nold == 0) THEN
    1694              :                         tmp = d_ijkl_sp(sepi, sepj, ij, kl, li, lj, lk, ll, 0, rij, &
    1695      4644022 :                                         se_int_control, se_int_screen, do_method_undef)
    1696      4644022 :                         drep(numb) = dft*rep(numb) + ft*tmp
    1697              :                      END IF
    1698              :                   END IF
    1699              :                END DO
    1700              :             END DO
    1701              :          END DO
    1702              :       END DO
    1703       511306 :    END SUBROUTINE dterep_sp_ana
    1704              : 
    1705              : ! **************************************************************************************************
    1706              : !> \brief Calculates the derivatives of the two-electron repulsion integrals - d shell only
    1707              : !> \param sepi parameters of atom i
    1708              : !> \param sepj parameters of atom j
    1709              : !> \param rij interatomic distance
    1710              : !> \param drep ...
    1711              : !> \param rep array of two-electron repulsion integrals
    1712              : !> \param se_int_control input parameters that control the calculation of
    1713              : !>                         integrals (shortrange, R3 residual, screening type)
    1714              : !> \param se_int_screen ...
    1715              : !> \param ft ...
    1716              : !> \param dft ...
    1717              : !> \par History
    1718              : !>      Teodoro Laino (05.2008) [tlaino] - University of Zurich : new driver
    1719              : !>                 for computing integral derivatives for d-orbitals
    1720              : ! **************************************************************************************************
    1721        25368 :    SUBROUTINE dterep_d_ana(sepi, sepj, rij, drep, rep, se_int_control, &
    1722              :                            se_int_screen, ft, dft)
    1723              :       TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
    1724              :       REAL(dp), INTENT(IN)                               :: rij
    1725              :       REAL(dp), DIMENSION(491), INTENT(INOUT)            :: drep
    1726              :       REAL(dp), DIMENSION(491), INTENT(IN)               :: rep
    1727              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
    1728              :       TYPE(se_int_screen_type), INTENT(IN)               :: se_int_screen
    1729              :       REAL(dp), INTENT(IN)                               :: ft, dft
    1730              : 
    1731              :       INTEGER                                            :: i, ij, j, k, kl, l, lasti, lastj, li, &
    1732              :                                                             lj, lk, ll, nold, numb
    1733              :       REAL(KIND=dp)                                      :: tmp
    1734              : 
    1735        25368 :       lasti = sepi%natorb
    1736        25368 :       lastj = sepj%natorb
    1737       196725 :       DO i = 1, lasti
    1738       171357 :          li = l_index(i)
    1739       965340 :          DO j = 1, i
    1740       768615 :             lj = l_index(j)
    1741       768615 :             ij = indexa(i, j)
    1742      4560222 :             DO k = 1, lastj
    1743      3620250 :                lk = l_index(k)
    1744     17912985 :                DO l = 1, k
    1745     13524120 :                   ll = l_index(l)
    1746     13524120 :                   kl = indexa(k, l)
    1747              : 
    1748     13524120 :                   numb = ijkl_ind(ij, kl)
    1749              :                   ! From 1 to 34 we store integrals involving sp shells
    1750     17144370 :                   IF (numb > 34) THEN
    1751      2836440 :                      nold = ijkl_sym(numb)
    1752      2836440 :                      IF (nold > 34) THEN
    1753      1324558 :                         drep(numb) = drep(nold)
    1754      1511882 :                      ELSE IF (nold < -34) THEN
    1755       259002 :                         drep(numb) = -drep(-nold)
    1756      1252880 :                      ELSE IF (nold == 0) THEN
    1757              :                         tmp = d_ijkl_d(sepi, sepj, ij, kl, li, lj, lk, ll, 0, rij, &
    1758      1252880 :                                        se_int_control, se_int_screen, do_method_undef)
    1759      1252880 :                         drep(numb) = dft*rep(numb) + ft*tmp
    1760              :                      END IF
    1761              :                   END IF
    1762              :                END DO
    1763              :             END DO
    1764              :          END DO
    1765              :       END DO
    1766        25368 :    END SUBROUTINE dterep_d_ana
    1767              : 
    1768              : END MODULE semi_empirical_int_ana
        

Generated by: LCOV version 2.0-1