LCOV - code coverage report
Current view: top level - src - qs_dftb_utils.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 95.6 % 363 347
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 18 18

            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 Working with the DFTB parameter types.
      10              : !> \author JGH (24.02.2007)
      11              : ! **************************************************************************************************
      12              : MODULE qs_dftb_utils
      13              : 
      14              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      15              :                                               cp_logger_type
      16              :    USE cp_output_handling,              ONLY: cp_p_file,&
      17              :                                               cp_print_key_finished_output,&
      18              :                                               cp_print_key_should_output,&
      19              :                                               cp_print_key_unit_nr
      20              :    USE input_section_types,             ONLY: section_vals_type
      21              :    USE kinds,                           ONLY: default_string_length,&
      22              :                                               dp
      23              :    USE qs_dftb_types,                   ONLY: qs_dftb_atom_type
      24              : #include "./base/base_uses.f90"
      25              : 
      26              :    IMPLICIT NONE
      27              : 
      28              :    PRIVATE
      29              : 
      30              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dftb_utils'
      31              : 
      32              :    ! Maximum number of points used for interpolation
      33              :    INTEGER, PARAMETER                     :: max_inter = 5
      34              :    ! Maximum number of points used for extrapolation
      35              :    INTEGER, PARAMETER                     :: max_extra = 9
      36              :    ! see also qs_dftb_parameters
      37              :    REAL(dp), PARAMETER                    :: slako_d0 = 1._dp
      38              :    ! pointer to skab
      39              :    INTEGER, DIMENSION(0:3, 0:3, 0:3, 0:3, 0:3):: iptr
      40              :    ! small real number
      41              :    REAL(dp), PARAMETER                    :: rtiny = 1.e-10_dp
      42              :    ! eta(0) for mm atoms and non-scc qm atoms
      43              :    REAL(dp), PARAMETER                    :: eta_mm = 0.47_dp
      44              :    ! step size for qmmm finite difference
      45              :    REAL(dp), PARAMETER                    :: ddrmm = 0.0001_dp
      46              : 
      47              :    PUBLIC :: allocate_dftb_atom_param, &
      48              :              deallocate_dftb_atom_param, &
      49              :              get_dftb_atom_param, &
      50              :              set_dftb_atom_param, &
      51              :              write_dftb_atom_param
      52              :    PUBLIC :: compute_block_sk, &
      53              :              urep_egr, iptr
      54              : 
      55              : CONTAINS
      56              : 
      57              : ! **************************************************************************************************
      58              : !> \brief ...
      59              : !> \param dftb_parameter ...
      60              : ! **************************************************************************************************
      61          480 :    SUBROUTINE allocate_dftb_atom_param(dftb_parameter)
      62              : 
      63              :       TYPE(qs_dftb_atom_type), POINTER                   :: dftb_parameter
      64              : 
      65          480 :       IF (ASSOCIATED(dftb_parameter)) &
      66            0 :          CALL deallocate_dftb_atom_param(dftb_parameter)
      67              : 
      68         7200 :       ALLOCATE (dftb_parameter)
      69              : 
      70              :       dftb_parameter%defined = .FALSE.
      71          480 :       dftb_parameter%name = ""
      72          480 :       dftb_parameter%typ = "NONE"
      73              :       dftb_parameter%z = -1
      74              :       dftb_parameter%zeff = -1.0_dp
      75          480 :       dftb_parameter%natorb = 0
      76              :       dftb_parameter%lmax = -1
      77         2400 :       dftb_parameter%skself = 0.0_dp
      78         2400 :       dftb_parameter%occupation = 0.0_dp
      79         2400 :       dftb_parameter%eta = 0.0_dp
      80          480 :       dftb_parameter%energy = 0.0_dp
      81          480 :       dftb_parameter%xi = 0.0_dp
      82          480 :       dftb_parameter%di = 0.0_dp
      83          480 :       dftb_parameter%rcdisp = 0.0_dp
      84          480 :       dftb_parameter%dudq = 0.0_dp
      85              : 
      86          480 :    END SUBROUTINE allocate_dftb_atom_param
      87              : 
      88              : ! **************************************************************************************************
      89              : !> \brief ...
      90              : !> \param dftb_parameter ...
      91              : ! **************************************************************************************************
      92          480 :    SUBROUTINE deallocate_dftb_atom_param(dftb_parameter)
      93              : 
      94              :       TYPE(qs_dftb_atom_type), POINTER                   :: dftb_parameter
      95              : 
      96          480 :       CPASSERT(ASSOCIATED(dftb_parameter))
      97          480 :       DEALLOCATE (dftb_parameter)
      98              : 
      99          480 :    END SUBROUTINE deallocate_dftb_atom_param
     100              : 
     101              : ! **************************************************************************************************
     102              : !> \brief ...
     103              : !> \param dftb_parameter ...
     104              : !> \param name ...
     105              : !> \param typ ...
     106              : !> \param defined ...
     107              : !> \param z ...
     108              : !> \param zeff ...
     109              : !> \param natorb ...
     110              : !> \param lmax ...
     111              : !> \param skself ...
     112              : !> \param occupation ...
     113              : !> \param eta ...
     114              : !> \param energy ...
     115              : !> \param cutoff ...
     116              : !> \param xi ...
     117              : !> \param di ...
     118              : !> \param rcdisp ...
     119              : !> \param dudq ...
     120              : ! **************************************************************************************************
     121      6770901 :    SUBROUTINE get_dftb_atom_param(dftb_parameter, name, typ, defined, z, zeff, natorb, &
     122              :                                   lmax, skself, occupation, eta, energy, cutoff, xi, di, rcdisp, dudq)
     123              : 
     124              :       TYPE(qs_dftb_atom_type), POINTER                   :: dftb_parameter
     125              :       CHARACTER(LEN=default_string_length), &
     126              :          INTENT(OUT), OPTIONAL                           :: name, typ
     127              :       LOGICAL, INTENT(OUT), OPTIONAL                     :: defined
     128              :       INTEGER, INTENT(OUT), OPTIONAL                     :: z
     129              :       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: zeff
     130              :       INTEGER, INTENT(OUT), OPTIONAL                     :: natorb, lmax
     131              :       REAL(KIND=dp), DIMENSION(0:3), OPTIONAL            :: skself, occupation, eta
     132              :       REAL(KIND=dp), OPTIONAL                            :: energy, cutoff, xi, di, rcdisp, dudq
     133              : 
     134      6770901 :       CPASSERT(ASSOCIATED(dftb_parameter))
     135              : 
     136      6770901 :       IF (PRESENT(name)) name = dftb_parameter%name
     137      6770901 :       IF (PRESENT(typ)) typ = dftb_parameter%typ
     138      6770901 :       IF (PRESENT(defined)) defined = dftb_parameter%defined
     139      6770901 :       IF (PRESENT(z)) z = dftb_parameter%z
     140      6770901 :       IF (PRESENT(zeff)) zeff = dftb_parameter%zeff
     141      6770901 :       IF (PRESENT(natorb)) natorb = dftb_parameter%natorb
     142      6770901 :       IF (PRESENT(lmax)) lmax = dftb_parameter%lmax
     143      9677401 :       IF (PRESENT(skself)) skself = dftb_parameter%skself
     144     32126085 :       IF (PRESENT(eta)) eta = dftb_parameter%eta
     145      6770901 :       IF (PRESENT(energy)) energy = dftb_parameter%energy
     146      6770901 :       IF (PRESENT(cutoff)) cutoff = dftb_parameter%cutoff
     147      6773253 :       IF (PRESENT(occupation)) occupation = dftb_parameter%occupation
     148      6770901 :       IF (PRESENT(xi)) xi = dftb_parameter%xi
     149      6770901 :       IF (PRESENT(di)) di = dftb_parameter%di
     150      6770901 :       IF (PRESENT(rcdisp)) rcdisp = dftb_parameter%rcdisp
     151      6770901 :       IF (PRESENT(dudq)) dudq = dftb_parameter%dudq
     152              : 
     153      6770901 :    END SUBROUTINE get_dftb_atom_param
     154              : 
     155              : ! **************************************************************************************************
     156              : !> \brief ...
     157              : !> \param dftb_parameter ...
     158              : !> \param name ...
     159              : !> \param typ ...
     160              : !> \param defined ...
     161              : !> \param z ...
     162              : !> \param zeff ...
     163              : !> \param natorb ...
     164              : !> \param lmax ...
     165              : !> \param skself ...
     166              : !> \param occupation ...
     167              : !> \param eta ...
     168              : !> \param energy ...
     169              : !> \param cutoff ...
     170              : !> \param xi ...
     171              : !> \param di ...
     172              : !> \param rcdisp ...
     173              : !> \param dudq ...
     174              : ! **************************************************************************************************
     175         1798 :    SUBROUTINE set_dftb_atom_param(dftb_parameter, name, typ, defined, z, zeff, natorb, &
     176              :                                   lmax, skself, occupation, eta, energy, cutoff, xi, di, rcdisp, dudq)
     177              : 
     178              :       TYPE(qs_dftb_atom_type), POINTER                   :: dftb_parameter
     179              :       CHARACTER(LEN=default_string_length), INTENT(IN), &
     180              :          OPTIONAL                                        :: name, typ
     181              :       LOGICAL, INTENT(IN), OPTIONAL                      :: defined
     182              :       INTEGER, INTENT(IN), OPTIONAL                      :: z
     183              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: zeff
     184              :       INTEGER, INTENT(IN), OPTIONAL                      :: natorb, lmax
     185              :       REAL(KIND=dp), DIMENSION(0:3), OPTIONAL            :: skself, occupation, eta
     186              :       REAL(KIND=dp), OPTIONAL                            :: energy, cutoff, xi, di, rcdisp, dudq
     187              : 
     188         1798 :       CPASSERT(ASSOCIATED(dftb_parameter))
     189              : 
     190         1798 :       IF (PRESENT(name)) dftb_parameter%name = name
     191         1798 :       IF (PRESENT(typ)) dftb_parameter%typ = typ
     192         1798 :       IF (PRESENT(defined)) dftb_parameter%defined = defined
     193         1798 :       IF (PRESENT(z)) dftb_parameter%z = z
     194         1798 :       IF (PRESENT(zeff)) dftb_parameter%zeff = zeff
     195         1798 :       IF (PRESENT(natorb)) dftb_parameter%natorb = natorb
     196         1798 :       IF (PRESENT(lmax)) dftb_parameter%lmax = lmax
     197         3718 :       IF (PRESENT(skself)) dftb_parameter%skself = skself
     198         3718 :       IF (PRESENT(eta)) dftb_parameter%eta = eta
     199         3718 :       IF (PRESENT(occupation)) dftb_parameter%occupation = occupation
     200         1798 :       IF (PRESENT(energy)) dftb_parameter%energy = energy
     201         1798 :       IF (PRESENT(cutoff)) dftb_parameter%cutoff = cutoff
     202         1798 :       IF (PRESENT(xi)) dftb_parameter%xi = xi
     203         1798 :       IF (PRESENT(di)) dftb_parameter%di = di
     204         1798 :       IF (PRESENT(rcdisp)) dftb_parameter%rcdisp = rcdisp
     205         1798 :       IF (PRESENT(dudq)) dftb_parameter%dudq = dudq
     206              : 
     207         1798 :    END SUBROUTINE set_dftb_atom_param
     208              : 
     209              : ! **************************************************************************************************
     210              : !> \brief ...
     211              : !> \param dftb_parameter ...
     212              : !> \param subsys_section ...
     213              : ! **************************************************************************************************
     214          480 :    SUBROUTINE write_dftb_atom_param(dftb_parameter, subsys_section)
     215              : 
     216              :       TYPE(qs_dftb_atom_type), POINTER                   :: dftb_parameter
     217              :       TYPE(section_vals_type), POINTER                   :: subsys_section
     218              : 
     219              :       CHARACTER(LEN=default_string_length)               :: name, typ
     220              :       INTEGER                                            :: lmax, natorb, output_unit, z
     221              :       LOGICAL                                            :: defined
     222              :       REAL(dp)                                           :: zeff
     223              :       TYPE(cp_logger_type), POINTER                      :: logger
     224              : 
     225          480 :       NULLIFY (logger)
     226          480 :       logger => cp_get_default_logger()
     227          480 :       IF (ASSOCIATED(dftb_parameter) .AND. &
     228              :           BTEST(cp_print_key_should_output(logger%iter_info, subsys_section, &
     229              :                                            "PRINT%KINDS/POTENTIAL"), cp_p_file)) THEN
     230              : 
     231              :          output_unit = cp_print_key_unit_nr(logger, subsys_section, "PRINT%KINDS", &
     232            0 :                                             extension=".Log")
     233              : 
     234            0 :          IF (output_unit > 0) THEN
     235              :             CALL get_dftb_atom_param(dftb_parameter, name=name, typ=typ, defined=defined, &
     236            0 :                                      z=z, zeff=zeff, natorb=natorb, lmax=lmax)
     237              : 
     238              :             WRITE (UNIT=output_unit, FMT="(/,A,T67,A14)") &
     239            0 :                " DFTB  parameters: ", TRIM(name)
     240            0 :             IF (defined) THEN
     241              :                WRITE (UNIT=output_unit, FMT="(T16,A,T71,F10.2)") &
     242            0 :                   "Effective core charge:", zeff
     243              :                WRITE (UNIT=output_unit, FMT="(T16,A,T71,I10)") &
     244            0 :                   "Number of orbitals:", natorb
     245              :             ELSE
     246              :                WRITE (UNIT=output_unit, FMT="(T55,A)") &
     247            0 :                   "Parameters are not defined"
     248              :             END IF
     249              :          END IF
     250              :          CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
     251            0 :                                            "PRINT%KINDS")
     252              :       END IF
     253              : 
     254          480 :    END SUBROUTINE write_dftb_atom_param
     255              : 
     256              : ! **************************************************************************************************
     257              : !> \brief ...
     258              : !> \param block ...
     259              : !> \param smatij ...
     260              : !> \param smatji ...
     261              : !> \param rij ...
     262              : !> \param ngrd ...
     263              : !> \param ngrdcut ...
     264              : !> \param dgrd ...
     265              : !> \param llm ...
     266              : !> \param lmaxi ...
     267              : !> \param lmaxj ...
     268              : !> \param irow ...
     269              : !> \param iatom ...
     270              : ! **************************************************************************************************
     271      6036961 :    SUBROUTINE compute_block_sk(block, smatij, smatji, rij, ngrd, ngrdcut, dgrd, &
     272              :                                llm, lmaxi, lmaxj, irow, iatom)
     273              :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: block, smatij, smatji
     274              :       REAL(KIND=dp), DIMENSION(3)                        :: rij
     275              :       INTEGER                                            :: ngrd, ngrdcut
     276              :       REAL(KIND=dp)                                      :: dgrd
     277              :       INTEGER                                            :: llm, lmaxi, lmaxj, irow, iatom
     278              : 
     279              :       REAL(KIND=dp)                                      :: dr
     280              :       REAL(KIND=dp), DIMENSION(20)                       :: skabij, skabji
     281              : 
     282     24147844 :       dr = SQRT(SUM(rij(:)**2))
     283      6036961 :       CALL getskz(smatij, skabij, dr, ngrd, ngrdcut, dgrd, llm)
     284      6036961 :       CALL getskz(smatji, skabji, dr, ngrd, ngrdcut, dgrd, llm)
     285      6036961 :       IF (irow == iatom) THEN
     286      3500481 :          CALL turnsk(block, skabji, skabij, rij, dr, lmaxi, lmaxj)
     287              :       ELSE
     288     10145920 :          CALL turnsk(block, skabij, skabji, -rij, dr, lmaxj, lmaxi)
     289              :       END IF
     290              : 
     291      6036961 :    END SUBROUTINE compute_block_sk
     292              : 
     293              : ! **************************************************************************************************
     294              : !> \brief Gets matrix elements on z axis, as they are stored in the tables
     295              : !> \param slakotab ...
     296              : !> \param skpar ...
     297              : !> \param dx ...
     298              : !> \param ngrd ...
     299              : !> \param ngrdcut ...
     300              : !> \param dgrd ...
     301              : !> \param llm ...
     302              : !> \author 07. Feb. 2004, TH
     303              : ! **************************************************************************************************
     304     12073922 :    SUBROUTINE getskz(slakotab, skpar, dx, ngrd, ngrdcut, dgrd, llm)
     305              :       REAL(dp), INTENT(in)                               :: slakotab(:, :), dx
     306              :       INTEGER, INTENT(in)                                :: ngrd, ngrdcut
     307              :       REAL(dp), INTENT(in)                               :: dgrd
     308              :       INTEGER, INTENT(in)                                :: llm
     309              :       REAL(dp), INTENT(out)                              :: skpar(llm)
     310              : 
     311              :       INTEGER                                            :: clgp
     312              : 
     313     34051382 :       skpar = 0._dp
     314              :       !
     315              :       ! Determine closest grid point
     316              :       !
     317     12073922 :       clgp = NINT(dx/dgrd)
     318              :       !
     319              :       ! Screen elements which are too far away
     320              :       !
     321     12073922 :       IF (clgp > ngrdcut) RETURN
     322              :       !
     323              :       ! The grid point is either contained in the table --> matrix element
     324              :       ! can be interpolated, or it is outside the table --> matrix element
     325              :       ! needs to be extrapolated.
     326              :       !
     327     12072948 :       IF (clgp > ngrd) THEN
     328              :          !
     329              :          ! Extrapolate external matrix elements if table does not finish with zero
     330              :          !
     331      2922364 :          CALL extrapol(slakotab, skpar, dx, ngrd, dgrd, llm)
     332              :       ELSE
     333              :          !
     334              :          ! Interpolate tabulated matrix elements
     335              :          !
     336      9150584 :          CALL interpol(slakotab, skpar, dx, ngrd, dgrd, llm, clgp)
     337              :       END IF
     338              :    END SUBROUTINE getskz
     339              : 
     340              : ! **************************************************************************************************
     341              : !> \brief ...
     342              : !> \param slakotab ...
     343              : !> \param skpar ...
     344              : !> \param dx ...
     345              : !> \param ngrd ...
     346              : !> \param dgrd ...
     347              : !> \param llm ...
     348              : !> \param clgp ...
     349              : ! **************************************************************************************************
     350      9150584 :    SUBROUTINE interpol(slakotab, skpar, dx, ngrd, dgrd, llm, clgp)
     351              :       REAL(dp), INTENT(in)                               :: slakotab(:, :), dx
     352              :       INTEGER, INTENT(in)                                :: ngrd
     353              :       REAL(dp), INTENT(in)                               :: dgrd
     354              :       INTEGER, INTENT(in)                                :: llm
     355              :       REAL(dp), INTENT(out)                              :: skpar(llm)
     356              :       INTEGER, INTENT(in)                                :: clgp
     357              : 
     358              :       INTEGER                                            :: fgpm, k, l, lgpm
     359              :       REAL(dp)                                           :: error, xa(max_inter), ya(max_inter)
     360              : 
     361      9150584 :       lgpm = MIN(clgp + INT(max_inter/2.0), ngrd)
     362      9150584 :       fgpm = lgpm - max_inter + 1
     363     54903504 :       DO k = 0, max_inter - 1
     364     54903504 :          xa(k + 1) = (fgpm + k)*dgrd
     365              :       END DO
     366              :       !
     367              :       ! Interpolate matrix elements for all orbitals
     368              :       !
     369     25890490 :       DO l = 1, llm
     370              :          !
     371              :          ! Read SK parameters from table
     372              :          !
     373    100439436 :          ya(1:max_inter) = slakotab(fgpm:lgpm, l)
     374     25890490 :          CALL polint(xa, ya, max_inter, dx, skpar(l), error)
     375              :       END DO
     376      9150584 :    END SUBROUTINE interpol
     377              : 
     378              : ! **************************************************************************************************
     379              : !> \brief ...
     380              : !> \param slakotab ...
     381              : !> \param skpar ...
     382              : !> \param dx ...
     383              : !> \param ngrd ...
     384              : !> \param dgrd ...
     385              : !> \param llm ...
     386              : ! **************************************************************************************************
     387      2922364 :    SUBROUTINE extrapol(slakotab, skpar, dx, ngrd, dgrd, llm)
     388              :       REAL(dp), INTENT(in)                               :: slakotab(:, :), dx
     389              :       INTEGER, INTENT(in)                                :: ngrd
     390              :       REAL(dp), INTENT(in)                               :: dgrd
     391              :       INTEGER, INTENT(in)                                :: llm
     392              :       REAL(dp), INTENT(out)                              :: skpar(llm)
     393              : 
     394              :       INTEGER                                            :: fgp, k, l, lgp, ntable, nzero
     395              :       REAL(dp)                                           :: error, xa(max_extra), ya(max_extra)
     396              : 
     397      2922364 :       nzero = max_extra/3
     398      2922364 :       ntable = max_extra - nzero
     399              :       !
     400              :       ! Get the three last distances from the table
     401              :       !
     402     20456548 :       DO k = 1, ntable
     403     20456548 :          xa(k) = (ngrd - (max_extra - 3) + k)*dgrd
     404              :       END DO
     405     11689456 :       DO k = 1, nzero
     406      8767092 :          xa(ntable + k) = (ngrd + k - 1)*dgrd + slako_d0
     407     11689456 :          ya(ntable + k) = 0.0
     408              :       END DO
     409              :       !
     410              :       ! Extrapolate matrix elements for all orbitals
     411              :       !
     412      8158262 :       DO l = 1, llm
     413              :          !
     414              :          ! Read SK parameters from table
     415              :          !
     416      5235898 :          fgp = ngrd + 1 - (max_extra - 3)
     417      5235898 :          lgp = ngrd
     418     36651286 :          ya(1:max_extra - 3) = slakotab(fgp:lgp, l)
     419      8158262 :          CALL polint(xa, ya, max_extra, dx, skpar(l), error)
     420              :       END DO
     421      2922364 :    END SUBROUTINE extrapol
     422              : 
     423              : ! **************************************************************************************************
     424              : !> \brief   Turn matrix element from z-axis to orientation of dxv
     425              : !> \param mat ...
     426              : !> \param skab1 ...
     427              : !> \param skab2 ...
     428              : !> \param dxv ...
     429              : !> \param dx ...
     430              : !> \param lmaxa ...
     431              : !> \param lmaxb ...
     432              : !> \date    13. Jan 2004
     433              : !> \par Notes
     434              : !>          These routines are taken from an old TB code (unknown to TH).
     435              : !>          They are highly optimised and taken because they are time critical.
     436              : !>          They are explicit, so not recursive, and work up to d functions.
     437              : !>
     438              : !>          Set variables necessary for rotation of matrix elements
     439              : !>
     440              : !>          r_i^2/r, replicated in rr2(4:6) for index convenience later
     441              : !>          r_i/r, direction vector, rr(4:6) are replicated from 1:3
     442              : !>          lmax of A and B
     443              : !> \author  TH
     444              : !> \version 1.0
     445              : ! **************************************************************************************************
     446      6036961 :    SUBROUTINE turnsk(mat, skab1, skab2, dxv, dx, lmaxa, lmaxb)
     447              :       REAL(dp), INTENT(inout)                  :: mat(:, :)
     448              :       REAL(dp), INTENT(in)                     :: skab1(:), skab2(:), dxv(3), dx
     449              :       INTEGER, INTENT(in)                      :: lmaxa, lmaxb
     450              : 
     451              :       INTEGER                                  :: lmaxab, minlmaxab
     452              :       REAL(dp)                                 :: rinv, rr(6), rr2(6)
     453              : 
     454      6036961 :       lmaxab = MAX(lmaxa, lmaxb)
     455              :       ! Determine l quantum limits.
     456      6036961 :       IF (lmaxab .GT. 2) CPABORT('lmax=2')
     457      6036961 :       minlmaxab = MIN(lmaxa, lmaxb)
     458              :       !
     459              :       ! s-s interaction
     460              :       !
     461      6036961 :       CALL skss(skab1, mat)
     462              :       !
     463      6036961 :       IF (lmaxab .LE. 0) RETURN
     464              :       !
     465     13648156 :       rr2(1:3) = dxv(1:3)**2
     466     13648156 :       rr(1:3) = dxv(1:3)
     467      3412039 :       rinv = 1.0_dp/dx
     468              :       !
     469     13648156 :       rr(1:3) = rr(1:3)*rinv
     470     13648156 :       rr(4:6) = rr(1:3)
     471     13648156 :       rr2(1:3) = rr2(1:3)*rinv**2
     472     13648156 :       rr2(4:6) = rr2(1:3)
     473              :       !
     474              :       ! s-p, p-s and p-p interaction
     475              :       !
     476      3412039 :       IF (minlmaxab .GE. 1) THEN
     477       714089 :          CALL skpp(skab1, mat, iptr(:, :, :, lmaxa, lmaxb))
     478       714089 :          CALL sksp(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
     479       714089 :          CALL sksp(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.)
     480              :       ELSE
     481      2697950 :          IF (lmaxb .GE. 1) THEN
     482      1301518 :             CALL sksp(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
     483              :          ELSE
     484      1396432 :             CALL sksp(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.)
     485              :          END IF
     486              :       END IF
     487              :       !
     488              :       ! If there is only s-p interaction we have finished
     489              :       !
     490      3412039 :       IF (lmaxab .LE. 1) RETURN
     491              :       !
     492              :       ! at least one atom has d functions
     493              :       !
     494        18608 :       IF (minlmaxab .EQ. 2) THEN
     495              :          !
     496              :          ! in case both atoms have d functions
     497              :          !
     498        18576 :          CALL skdd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb))
     499        18576 :          CALL sksd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
     500        18576 :          CALL sksd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.)
     501        18576 :          CALL skpd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
     502        18576 :          CALL skpd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.)
     503              :       ELSE
     504              :          !
     505              :          ! One atom has d functions, the other has s or s and p functions
     506              :          !
     507           32 :          IF (lmaxa .EQ. 0) THEN
     508              :             !
     509              :             ! atom b has d, the atom a only s functions
     510              :             !
     511            0 :             CALL sksd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
     512           32 :          ELSE IF (lmaxa .EQ. 1) THEN
     513              :             !
     514              :             ! atom b has d, the atom a s and p functions
     515              :             !
     516            0 :             CALL sksd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
     517            0 :             CALL skpd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.)
     518              :          ELSE
     519              :             !
     520              :             ! atom a has d functions
     521              :             !
     522           32 :             IF (lmaxb .EQ. 0) THEN
     523              :                !
     524              :                ! atom a has d, atom b has only s functions
     525              :                !
     526            0 :                CALL sksd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.)
     527              :             ELSE
     528              :                !
     529              :                ! atom a has d, atom b has s and p functions
     530              :                !
     531           32 :                CALL sksd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.)
     532           32 :                CALL skpd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.)
     533              :             END IF
     534              :          END IF
     535              :       END IF
     536              :       !
     537              :    CONTAINS
     538              :       !
     539              :       ! The subroutines to turn the matrix elements are taken as internal subroutines
     540              :       ! as it is beneficial to inline them.
     541              :       !
     542              :       ! They are both turning the matrix elements and placing them appropriately
     543              :       ! into the matrix block
     544              :       !
     545              : ! **************************************************************************************************
     546              : !> \brief   s-s interaction (no rotation necessary)
     547              : !> \param skpar ...
     548              : !> \param mat ...
     549              : !> \version 1.0
     550              : ! **************************************************************************************************
     551      6036961 :       SUBROUTINE skss(skpar, mat)
     552              :       REAL(dp), INTENT(in)                               :: skpar(:)
     553              :       REAL(dp), INTENT(inout)                            :: mat(:, :)
     554              : 
     555      6036961 :          mat(1, 1) = mat(1, 1) + skpar(1)
     556              :          !
     557      6036961 :       END SUBROUTINE skss
     558              : 
     559              : ! **************************************************************************************************
     560              : !> \brief  s-p interaction (simple rotation)
     561              : !> \param skpar ...
     562              : !> \param mat ...
     563              : !> \param ind ...
     564              : !> \param transposed ...
     565              : !> \version 1.0
     566              : ! **************************************************************************************************
     567      4126128 :       SUBROUTINE sksp(skpar, mat, ind, transposed)
     568              :       REAL(dp), INTENT(in)                               :: skpar(:)
     569              :       REAL(dp), INTENT(inout)                            :: mat(:, :)
     570              :       INTEGER, INTENT(in)                                :: ind(0:, 0:, 0:)
     571              :       LOGICAL, INTENT(in)                                :: transposed
     572              : 
     573              :       INTEGER                                            :: l
     574              :       REAL(dp)                                           :: skp
     575              : 
     576      4126128 :          skp = skpar(ind(1, 0, 0))
     577      4126128 :          IF (transposed) THEN
     578      8062428 :             DO l = 1, 3
     579      8062428 :                mat(1, l + 1) = mat(1, l + 1) + rr(l)*skp
     580              :             END DO
     581              :          ELSE
     582      8442084 :             DO l = 1, 3
     583      8442084 :                mat(l + 1, 1) = mat(l + 1, 1) - rr(l)*skp
     584              :             END DO
     585              :          END IF
     586              :          !
     587      4126128 :       END SUBROUTINE sksp
     588              : 
     589              : ! **************************************************************************************************
     590              : !> \brief ...
     591              : !> \param skpar ...
     592              : !> \param mat ...
     593              : !> \param ind ...
     594              : ! **************************************************************************************************
     595       714089 :       SUBROUTINE skpp(skpar, mat, ind)
     596              :       REAL(dp), INTENT(in)                               :: skpar(:)
     597              :       REAL(dp), INTENT(inout)                            :: mat(:, :)
     598              :       INTEGER, INTENT(in)                                :: ind(0:, 0:, 0:)
     599              : 
     600              :       INTEGER                                            :: ii, ir, is, k, l
     601              :       REAL(dp)                                           :: epp(6), matel(6), skppp, skpps
     602              : 
     603      2856356 :          epp(1:3) = rr2(1:3)
     604      2856356 :          DO l = 1, 3
     605      2856356 :             epp(l + 3) = rr(l)*rr(l + 1)
     606              :          END DO
     607       714089 :          skppp = skpar(ind(1, 1, 1))
     608       714089 :          skpps = skpar(ind(1, 1, 0))
     609              :          !
     610      2856356 :          DO l = 1, 3
     611      2856356 :             matel(l) = epp(l)*skpps + (1._dp - epp(l))*skppp
     612              :          END DO
     613      2856356 :          DO l = 4, 6
     614      2856356 :             matel(l) = epp(l)*(skpps - skppp)
     615              :          END DO
     616              :          !
     617      2856356 :          DO ir = 1, 3
     618      4284534 :             DO is = 1, ir - 1
     619      2142267 :                ii = ir - is
     620      2142267 :                k = 3*ii - (ii*(ii - 1))/2 + is
     621      2142267 :                mat(is + 1, ir + 1) = mat(is + 1, ir + 1) + matel(k)
     622      4284534 :                mat(ir + 1, is + 1) = mat(ir + 1, is + 1) + matel(k)
     623              :             END DO
     624      2856356 :             mat(ir + 1, ir + 1) = mat(ir + 1, ir + 1) + matel(ir)
     625              :          END DO
     626       714089 :       END SUBROUTINE skpp
     627              : 
     628              : ! **************************************************************************************************
     629              : !> \brief ...
     630              : !> \param skpar ...
     631              : !> \param mat ...
     632              : !> \param ind ...
     633              : !> \param transposed ...
     634              : ! **************************************************************************************************
     635        37184 :       SUBROUTINE sksd(skpar, mat, ind, transposed)
     636              :       REAL(dp), INTENT(in)                               :: skpar(:)
     637              :       REAL(dp), INTENT(inout)                            :: mat(:, :)
     638              :       INTEGER, INTENT(in)                                :: ind(0:, 0:, 0:)
     639              :       LOGICAL, INTENT(in)                                :: transposed
     640              : 
     641              :       INTEGER                                            :: l
     642              :       REAL(dp)                                           :: d4, d5, es(5), r3, sksds
     643              : 
     644        37184 :          sksds = skpar(ind(2, 0, 0))
     645        37184 :          r3 = SQRT(3._dp)
     646        37184 :          d4 = rr2(3) - 0.5_dp*(rr2(1) + rr2(2))
     647        37184 :          d5 = rr2(1) - rr2(2)
     648              :          !
     649       148736 :          DO l = 1, 3
     650       148736 :             es(l) = r3*rr(l)*rr(l + 1)
     651              :          END DO
     652        37184 :          es(4) = 0.5_dp*r3*d5
     653        37184 :          es(5) = d4
     654              :          !
     655        37184 :          IF (transposed) THEN
     656       111456 :             DO l = 1, 5
     657       111456 :                mat(1, l + 4) = mat(1, l + 4) + es(l)*sksds
     658              :             END DO
     659              :          ELSE
     660       111648 :             DO l = 1, 5
     661       111648 :                mat(l + 4, 1) = mat(l + 4, 1) + es(l)*sksds
     662              :             END DO
     663              :          END IF
     664        37184 :       END SUBROUTINE sksd
     665              : 
     666              : ! **************************************************************************************************
     667              : !> \brief ...
     668              : !> \param skpar ...
     669              : !> \param mat ...
     670              : !> \param ind ...
     671              : !> \param transposed ...
     672              : ! **************************************************************************************************
     673        37184 :       SUBROUTINE skpd(skpar, mat, ind, transposed)
     674              :       REAL(dp), INTENT(in)                               :: skpar(:)
     675              :       REAL(dp), INTENT(inout)                            :: mat(:, :)
     676              :       INTEGER, INTENT(in)                                :: ind(0:, 0:, 0:)
     677              :       LOGICAL, INTENT(in)                                :: transposed
     678              : 
     679              :       INTEGER                                            :: ir, is, k, l, m
     680              :       REAL(dp)                                           :: d3, d4, d5, d6, dm(15), epd(13, 2), r3, &
     681              :                                                             sktmp
     682              : 
     683        37184 :          r3 = SQRT(3.0_dp)
     684        37184 :          d3 = rr2(1) + rr2(2)
     685        37184 :          d4 = rr2(3) - 0.5_dp*d3
     686        37184 :          d5 = rr2(1) - rr2(2)
     687        37184 :          d6 = rr(1)*rr(2)*rr(3)
     688       148736 :          DO l = 1, 3
     689       111552 :             epd(l, 1) = r3*rr2(l)*rr(l + 1)
     690       111552 :             epd(l, 2) = rr(l + 1)*(1.0_dp - 2._dp*rr2(l))
     691       111552 :             epd(l + 4, 1) = r3*rr2(l)*rr(l + 2)
     692       111552 :             epd(l + 4, 2) = rr(l + 2)*(1.0_dp - 2*rr2(l))
     693       111552 :             epd(l + 7, 1) = 0.5_dp*r3*rr(l)*d5
     694       148736 :             epd(l + 10, 1) = rr(l)*d4
     695              :          END DO
     696              :          !
     697        37184 :          epd(4, 1) = r3*d6
     698        37184 :          epd(4, 2) = -2._dp*d6
     699        37184 :          epd(8, 2) = rr(1)*(1.0_dp - d5)
     700        37184 :          epd(9, 2) = -rr(2)*(1.0_dp + d5)
     701        37184 :          epd(10, 2) = -rr(3)*d5
     702        37184 :          epd(11, 2) = -r3*rr(1)*rr2(3)
     703        37184 :          epd(12, 2) = -r3*rr(2)*rr2(3)
     704        37184 :          epd(13, 2) = r3*rr(3)*d3
     705              :          !
     706        37184 :          dm(1:15) = 0.0_dp
     707              :          !
     708       111552 :          DO m = 1, 2
     709        74368 :             sktmp = skpar(ind(2, 1, m - 1))
     710        74368 :             dm(1) = dm(1) + epd(1, m)*sktmp
     711        74368 :             dm(2) = dm(2) + epd(6, m)*sktmp
     712        74368 :             dm(3) = dm(3) + epd(4, m)*sktmp
     713        74368 :             dm(5) = dm(5) + epd(2, m)*sktmp
     714        74368 :             dm(6) = dm(6) + epd(7, m)*sktmp
     715        74368 :             dm(7) = dm(7) + epd(5, m)*sktmp
     716        74368 :             dm(9) = dm(9) + epd(3, m)*sktmp
     717       557760 :             DO l = 8, 13
     718       520576 :                dm(l + 2) = dm(l + 2) + epd(l, m)*sktmp
     719              :             END DO
     720              :          END DO
     721              :          !
     722        37184 :          dm(4) = dm(3)
     723        37184 :          dm(8) = dm(3)
     724              :          !
     725        37184 :          IF (transposed) THEN
     726       111456 :             DO ir = 1, 5
     727       390096 :                DO is = 1, 3
     728       278640 :                   k = 3*(ir - 1) + is
     729       371520 :                   mat(is + 1, ir + 4) = mat(is + 1, ir + 4) + dm(k)
     730              :                END DO
     731              :             END DO
     732              :          ELSE
     733       111648 :             DO ir = 1, 5
     734       390768 :                DO is = 1, 3
     735       279120 :                   k = 3*(ir - 1) + is
     736       372160 :                   mat(ir + 4, is + 1) = mat(ir + 4, is + 1) - dm(k)
     737              :                END DO
     738              :             END DO
     739              :          END IF
     740              :          !
     741        37184 :       END SUBROUTINE skpd
     742              : 
     743              : ! **************************************************************************************************
     744              : !> \brief ...
     745              : !> \param skpar ...
     746              : !> \param mat ...
     747              : !> \param ind ...
     748              : ! **************************************************************************************************
     749        18576 :       SUBROUTINE skdd(skpar, mat, ind)
     750              :       REAL(dp), INTENT(in)                               :: skpar(:)
     751              :       REAL(dp), INTENT(inout)                            :: mat(:, :)
     752              :       INTEGER, INTENT(in)                                :: ind(0:, 0:, 0:)
     753              : 
     754              :       INTEGER                                            :: ii, ir, is, k, l, m
     755              :       REAL(dp)                                           :: d3, d4, d5, dd(3), dm(15), e(15, 3), r3
     756              : 
     757        18576 :          r3 = SQRT(3._dp)
     758        18576 :          d3 = rr2(1) + rr2(2)
     759        18576 :          d4 = rr2(3) - 0.5_dp*d3
     760        18576 :          d5 = rr2(1) - rr2(2)
     761        74304 :          DO l = 1, 3
     762        55728 :             e(l, 1) = rr2(l)*rr2(l + 1)
     763        55728 :             e(l, 2) = rr2(l) + rr2(l + 1) - 4._dp*e(l, 1)
     764        55728 :             e(l, 3) = rr2(l + 2) + e(l, 1)
     765        74304 :             e(l, 1) = 3._dp*e(l, 1)
     766              :          END DO
     767        18576 :          e(4, 1) = d5**2
     768        18576 :          e(4, 2) = d3 - e(4, 1)
     769        18576 :          e(4, 3) = rr2(3) + 0.25_dp*e(4, 1)
     770        18576 :          e(4, 1) = 0.75_dp*e(4, 1)
     771        18576 :          e(5, 1) = d4**2
     772        18576 :          e(5, 2) = 3._dp*rr2(3)*d3
     773        18576 :          e(5, 3) = 0.75_dp*d3**2
     774        18576 :          dd(1) = rr(1)*rr(3)
     775        18576 :          dd(2) = rr(2)*rr(1)
     776        18576 :          dd(3) = rr(3)*rr(2)
     777        55728 :          DO l = 1, 2
     778        37152 :             e(l + 5, 1) = 3._dp*rr2(l + 1)*dd(l)
     779        37152 :             e(l + 5, 2) = dd(l)*(1._dp - 4._dp*rr2(l + 1))
     780        55728 :             e(l + 5, 3) = dd(l)*(rr2(l + 1) - 1._dp)
     781              :          END DO
     782        18576 :          e(8, 1) = dd(1)*d5*1.5_dp
     783        18576 :          e(8, 2) = dd(1)*(1.0_dp - 2.0_dp*d5)
     784        18576 :          e(8, 3) = dd(1)*(0.5_dp*d5 - 1.0_dp)
     785        18576 :          e(9, 1) = d5*0.5_dp*d4*r3
     786        18576 :          e(9, 2) = -d5*rr2(3)*r3
     787        18576 :          e(9, 3) = d5*0.25_dp*(1.0_dp + rr2(3))*r3
     788        18576 :          e(10, 1) = rr2(1)*dd(3)*3.0_dp
     789        18576 :          e(10, 2) = (0.25_dp - rr2(1))*dd(3)*4.0_dp
     790        18576 :          e(10, 3) = dd(3)*(rr2(1) - 1.0_dp)
     791        18576 :          e(11, 1) = 1.5_dp*dd(3)*d5
     792        18576 :          e(11, 2) = -dd(3)*(1.0_dp + 2.0_dp*d5)
     793        18576 :          e(11, 3) = dd(3)*(1.0_dp + 0.5_dp*d5)
     794        18576 :          e(13, 3) = 0.5_dp*d5*dd(2)
     795        18576 :          e(13, 2) = -2.0_dp*dd(2)*d5
     796        18576 :          e(13, 1) = e(13, 3)*3.0_dp
     797        18576 :          e(12, 1) = d4*dd(1)*r3
     798        18576 :          e(14, 1) = d4*dd(3)*r3
     799        18576 :          e(15, 1) = d4*dd(2)*r3
     800        18576 :          e(15, 2) = -2.0_dp*r3*dd(2)*rr2(3)
     801        18576 :          e(15, 3) = 0.5_dp*r3*(1.0_dp + rr2(3))*dd(2)
     802        18576 :          e(14, 2) = r3*dd(3)*(d3 - rr2(3))
     803        18576 :          e(14, 3) = -r3*0.5_dp*dd(3)*d3
     804        18576 :          e(12, 2) = r3*dd(1)*(d3 - rr2(3))
     805        18576 :          e(12, 3) = -r3*0.5_dp*dd(1)*d3
     806              :          !
     807        18576 :          dm(1:15) = 0._dp
     808       297216 :          DO l = 1, 15
     809      1133136 :             DO m = 1, 3
     810      1114560 :                dm(l) = dm(l) + e(l, m)*skpar(ind(2, 2, m - 1))
     811              :             END DO
     812              :          END DO
     813              :          !
     814       111456 :          DO ir = 1, 5
     815       278640 :             DO is = 1, ir - 1
     816       185760 :                ii = ir - is
     817       185760 :                k = 5*ii - (ii*(ii - 1))/2 + is
     818       185760 :                mat(ir + 4, is + 4) = mat(ir + 4, is + 4) + dm(k)
     819       278640 :                mat(is + 4, ir + 4) = mat(is + 4, ir + 4) + dm(k)
     820              :             END DO
     821       111456 :             mat(ir + 4, ir + 4) = mat(ir + 4, ir + 4) + dm(ir)
     822              :          END DO
     823        18576 :       END SUBROUTINE skdd
     824              :       !
     825              :    END SUBROUTINE turnsk
     826              : 
     827              : ! **************************************************************************************************
     828              : !> \brief ...
     829              : !> \param xa ...
     830              : !> \param ya ...
     831              : !> \param n ...
     832              : !> \param x ...
     833              : !> \param y ...
     834              : !> \param dy ...
     835              : ! **************************************************************************************************
     836     21975804 :    SUBROUTINE polint(xa, ya, n, x, y, dy)
     837              :       INTEGER, INTENT(in)                                :: n
     838              :       REAL(dp), INTENT(in)                               :: ya(n), xa(n), x
     839              :       REAL(dp), INTENT(out)                              :: y, dy
     840              : 
     841              :       INTEGER                                            :: i, m, ns
     842     43951608 :       REAL(dp)                                           :: c(n), d(n), den, dif, dift, ho, hp, w
     843              : 
     844              : !
     845              : !
     846              : 
     847     21975804 :       ns = 1
     848              : 
     849     21975804 :       dif = ABS(x - xa(1))
     850    152798416 :       DO i = 1, n
     851    130822612 :          dift = ABS(x - xa(i))
     852    130822612 :          IF (dift .LT. dif) THEN
     853     62888978 :             ns = i
     854     62888978 :             dif = dift
     855              :          END IF
     856    130822612 :          c(i) = ya(i)
     857    152798416 :          d(i) = ya(i)
     858              :       END DO
     859              :       !
     860     21975804 :       y = ya(ns)
     861     21975804 :       ns = ns - 1
     862    130822612 :       DO m = 1, n - 1
     863    464738196 :          DO i = 1, n - m
     864    355891388 :             ho = xa(i) - x
     865    355891388 :             hp = xa(i + m) - x
     866    355891388 :             w = c(i + 1) - d(i)
     867    355891388 :             den = ho - hp
     868    355891388 :             CPASSERT(den /= 0.0_dp)
     869    355891388 :             den = w/den
     870    355891388 :             d(i) = hp*den
     871    464738196 :             c(i) = ho*den
     872              :          END DO
     873    108846808 :          IF (2*ns .LT. n - m) THEN
     874     45957830 :             dy = c(ns + 1)
     875              :          ELSE
     876     62888978 :             dy = d(ns)
     877     62888978 :             ns = ns - 1
     878              :          END IF
     879    130822612 :          y = y + dy
     880              :       END DO
     881              : !
     882     21975804 :       RETURN
     883              :    END SUBROUTINE polint
     884              : 
     885              : ! **************************************************************************************************
     886              : !> \brief ...
     887              : !> \param rv ...
     888              : !> \param r ...
     889              : !> \param erep ...
     890              : !> \param derep ...
     891              : !> \param n_urpoly ...
     892              : !> \param urep ...
     893              : !> \param spdim ...
     894              : !> \param s_cut ...
     895              : !> \param srep ...
     896              : !> \param spxr ...
     897              : !> \param scoeff ...
     898              : !> \param surr ...
     899              : !> \param dograd ...
     900              : ! **************************************************************************************************
     901       542297 :    SUBROUTINE urep_egr(rv, r, erep, derep, &
     902       542297 :                        n_urpoly, urep, spdim, s_cut, srep, spxr, scoeff, surr, dograd)
     903              : 
     904              :       REAL(dp), INTENT(in)                               :: rv(3), r
     905              :       REAL(dp), INTENT(inout)                            :: erep, derep(3)
     906              :       INTEGER, INTENT(in)                                :: n_urpoly
     907              :       REAL(dp), INTENT(in)                               :: urep(:)
     908              :       INTEGER, INTENT(in)                                :: spdim
     909              :       REAL(dp), INTENT(in)                               :: s_cut, srep(3)
     910              :       REAL(dp), POINTER                                  :: spxr(:, :), scoeff(:, :)
     911              :       REAL(dp), INTENT(in)                               :: surr(2)
     912              :       LOGICAL, INTENT(in)                                :: dograd
     913              : 
     914              :       INTEGER                                            :: ic, isp, jsp, nsp
     915              :       REAL(dp)                                           :: de_z, rz
     916              : 
     917       542297 :       derep = 0._dp
     918       542297 :       de_z = 0._dp
     919       542297 :       IF (n_urpoly > 0) THEN
     920              :          !
     921              :          ! polynomial part
     922              :          !
     923         6975 :          rz = urep(1) - r
     924         6975 :          IF (rz <= rtiny) RETURN
     925        50562 :          DO ic = 2, n_urpoly
     926        50562 :             erep = erep + urep(ic)*rz**(ic)
     927              :          END DO
     928         6975 :          IF (dograd) THEN
     929        13889 :             DO ic = 2, n_urpoly
     930        13889 :                de_z = de_z - ic*urep(ic)*rz**(ic - 1)
     931              :             END DO
     932              :          END IF
     933       535322 :       ELSE IF (spdim > 0) THEN
     934              :          !
     935              :          ! spline part
     936              :          !
     937              :          ! This part is kind of proprietary Paderborn code and I won't reverse-engineer
     938              :          ! everything in detail. What is obvious is documented.
     939              :          !
     940              :          ! This part has 4 regions:
     941              :          ! a) very long range is screened
     942              :          ! b) short-range is extrapolated with e-functions
     943              :          ! ca) normal range is approximated with a spline
     944              :          ! cb) longer range is extrapolated with an higher degree spline
     945              :          !
     946       535322 :          IF (r > s_cut) RETURN ! screening (condition a)
     947              :          !
     948        15117 :          IF (r < spxr(1, 1)) THEN
     949              :             ! a) short range
     950            0 :             erep = erep + EXP(-srep(1)*r + srep(2)) + srep(3)
     951            0 :             IF (dograd) de_z = de_z - srep(1)*EXP(-srep(1)*r + srep(2))
     952              :          ELSE
     953              :             !
     954              :             ! condition c). First determine between which places the spline is located:
     955              :             !
     956       235681 :             ispg: DO isp = 1, spdim ! condition ca)
     957       235681 :                IF (r < spxr(isp, 1)) CYCLE ispg ! distance is smaller than this spline range
     958       235681 :                IF (r >= spxr(isp, 2)) CYCLE ispg ! distance is larger than this spline range
     959              :                ! at this point we have found the correct spline interval
     960        15117 :                rz = r - spxr(isp, 1)
     961        15117 :                IF (isp /= spdim) THEN
     962              :                   nsp = 3 ! condition ca
     963        59265 :                   DO jsp = 0, nsp
     964        59265 :                      erep = erep + scoeff(isp, jsp + 1)*rz**(jsp)
     965              :                   END DO
     966        11853 :                   IF (dograd) THEN
     967        19344 :                      DO jsp = 1, nsp
     968        19344 :                         de_z = de_z + jsp*scoeff(isp, jsp + 1)*rz**(jsp - 1)
     969              :                      END DO
     970              :                   END IF
     971              :                ELSE
     972              :                   nsp = 5 ! condition cb
     973        22848 :                   DO jsp = 0, nsp
     974        22848 :                      IF (jsp <= 3) THEN
     975        13056 :                         erep = erep + scoeff(isp, jsp + 1)*rz**(jsp)
     976              :                      ELSE
     977         6528 :                         erep = erep + surr(jsp - 3)*rz**(jsp)
     978              :                      END IF
     979              :                   END DO
     980         3264 :                   IF (dograd) THEN
     981         9792 :                      DO jsp = 1, nsp
     982         9792 :                         IF (jsp <= 3) THEN
     983         4896 :                            de_z = de_z + jsp*scoeff(isp, jsp + 1)*rz**(jsp - 1)
     984              :                         ELSE
     985         3264 :                            de_z = de_z + jsp*surr(jsp - 3)*rz**(jsp - 1)
     986              :                         END IF
     987              :                      END DO
     988              :                   END IF
     989              :                END IF
     990       220564 :                EXIT ispg
     991              :             END DO ispg
     992              :          END IF
     993              :       END IF
     994              :       !
     995        22092 :       IF (dograd) THEN
     996        33408 :          IF (r > 1.e-12_dp) derep(1:3) = (de_z/r)*rv(1:3)
     997              :       END IF
     998              : 
     999              :    END SUBROUTINE urep_egr
    1000              : 
    1001              : END MODULE qs_dftb_utils
    1002              : 
        

Generated by: LCOV version 2.0-1