LCOV - code coverage report
Current view: top level - src - semi_empirical_par_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 425 434 97.9 %
Date: 2024-03-28 07:31:50 Functions: 17 17 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Utilities to post-process semi-empirical parameters
      10             : !> \par History
      11             : !>         [tlaino] 03.2008 - Splitting from semi_empirical_parameters and
      12             : !>                            keeping there only the setting of the SE params
      13             : !> \author Teodoro Laino [tlaino] - University of Zurich
      14             : !> \date   03.2008 [tlaino]
      15             : ! **************************************************************************************************
      16             : MODULE semi_empirical_par_utils
      17             : 
      18             :    USE kinds,                           ONLY: dp
      19             :    USE mathconstants,                   ONLY: fac
      20             :    USE mathlib,                         ONLY: binomial
      21             :    USE physcon,                         ONLY: bohr,&
      22             :                                               evolt
      23             :    USE semi_empirical_int_arrays,       ONLY: int_ij,&
      24             :                                               int_kl,&
      25             :                                               int_onec2el
      26             :    USE semi_empirical_types,            ONLY: get_se_param,&
      27             :                                               semi_empirical_type
      28             : #include "./base/base_uses.f90"
      29             : 
      30             :    IMPLICIT NONE
      31             : 
      32             :    PRIVATE
      33             : 
      34             :    INTEGER, PARAMETER, PRIVATE :: nelem = 106
      35             : 
      36             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_par_utils'
      37             : 
      38             : !                STANDARD MOPAC PARAMETERS USED FOR AM1, RM1, MNDO, PM3, PM6,
      39             : !                PM6-FM
      40             : !
      41             : !      H                                                                      He
      42             : !      Li Be                                                 B  C  N  O  F    Ne
      43             : !      Na Mg                                                 Al Si P  S  Cl   Ar
      44             : !      K  Ca Sc                Ti V  Cr Mn Fe Co Ni Cu Zn    Ga Ge As Se Br   Kr
      45             : !      Rb Sr Y                 Zr Nb Mo Tc Ru Rh Pd Ag Cd    In Sn Sb Te I    Xe
      46             : !      Cs Ba La Ce-Lu          Hf Ta W  Re Os Ir Pt Au Hg    Tl Pb Bi Po At   Rn
      47             : !      Fr Ra Ac Th Pa U        Np Pu Am Cm Bk Cf Es Fm Md    No Lr Rf Ha 106
      48             : 
      49             : !                                      "s" shell
      50             :    INTEGER, DIMENSION(0:nelem), PRIVATE :: Nos = (/-1, & !    0
      51             :                                                    1, 2, & !    2
      52             :                                                    1, 2, 2, 2, 2, 2, 2, 0, & !   10
      53             :                                                    1, 2, 2, 2, 2, 2, 2, 0, & !   18
      54             :                                                    1, 2, 2, 2, 2, 1, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 0, & !   36
      55             :                                                    1, 2, 2, 2, 1, 1, 2, 1, 1, 0, 1, 2, 2, 2, 2, 2, 2, 0, & !   54
      56             :                                                    1, 2, 2, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, &
      57             :                                                    2, 2, 1, 2, 2, 2, 1, 1, 2, 2, 2, 2, 2, 2, 0, & !   86
      58             :                                                    1, 1, 2, 4, 2, 2, 2, 2, 2, 2, 2, 2, 0, 3, -3, 1, 2, 1, -2, -1/)
      59             : 
      60             : !                                      "p" shell
      61             :    INTEGER, DIMENSION(0:nelem), PRIVATE :: Nop = (/-1, & !    0
      62             :                                                    0, 0, & !    2
      63             :                                                    0, 0, 1, 2, 3, 4, 5, 6, & !   10
      64             :                                                    0, 0, 1, 2, 3, 4, 5, 6, & !   18
      65             :                                                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 4, 5, 6, & !   36
      66             :                                                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 4, 5, 6, & !   54
      67             :                                                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
      68             :                                                    0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 4, 5, 6, & !   86
      69             :                                                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/)
      70             : 
      71             : !                                      "d" shell
      72             :    INTEGER, DIMENSION(0:nelem), PRIVATE :: Nod = (/-1, & !    0
      73             :                                                    0, 0, & !    2
      74             :                                                    0, 0, 0, 0, 0, 0, 0, 0, & !   10
      75             :                                                    0, 0, 0, 0, 0, 0, 0, 0, & !   18
      76             :                                                    0, 0, 1, 2, 3, 5, 5, 6, 7, 8, 10, 0, 0, 0, 0, 0, 0, 0, & !   36
      77             :                                                    0, 0, 1, 2, 4, 5, 5, 7, 8, 10, 10, 0, 0, 0, 0, 0, 0, 0, & !   54
      78             :                                                    0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, &
      79             :                                                    2, 3, 5, 5, 6, 7, 9, 10, 0, 0, 0, 0, 0, 0, 0, & !   86
      80             :                                                    0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/)
      81             : 
      82             : !      H          <Quantum Numbers for s, p, d and f orbitals>                He
      83             : !      Li Be                                                 B  C  N  O  F    Ne
      84             : !      Na Mg                                                 Al Si P  S  Cl   Ar
      85             : !      K  Ca Sc                Ti V  Cr Mn Fe Co Ni Cu Zn    Ga Ge As Se Br   Kr
      86             : !      Rb Sr Y                 Zr Nb Mo Tc Ru Rh Pd Ag Cd    In Sn Sb Te I    Xe
      87             : !      Cs Ba La Ce-Lu          Hf Ta W  Re Os Ir Pt Au Hg    Tl Pb Bi Po At   Rn
      88             : !      Fr Ra Ac Th Pa U        Np Pu Am Cm Bk Cf Es Fm Md    No Lr Rf Ha 106
      89             : 
      90             :    INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE ::   nqs = (/-1, & !    0
      91             :                                                                 1, 1, & !    2
      92             :                                                                 2, 2, 2, 2, 2, 2, 2, 2, & !   10
      93             :                                                                 3, 3, 3, 3, 3, 3, 3, 3, & !   18
      94             :                                                                 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, & !   36
      95             :                                                                 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & !   54
      96             :                                                                 6, 6, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 6, &
      97             :                                                                 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, -1, -1, -1, & !   86
      98             :                                                     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1/)
      99             : 
     100             :    INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE ::   nqp = (/-1, & !    0
     101             :                                                                 -1, -1, & !    2
     102             :                                                                 2, 2, 2, 2, 2, 2, 2, 2, & !   10
     103             :                                                                 3, 3, 3, 3, 3, 3, 3, 3, & !   18
     104             :                                                                 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, & !   36
     105             :                                                                 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & !   54
     106             :                                                                 6, 6, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 6, &
     107             :                                                                 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, -1, -1, -1, & !   86
     108             :                                                     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1/)
     109             : 
     110             :    INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE ::   nqd = (/-1, & !    0
     111             :                                                                 -1, -1, & !    2
     112             :                                                                 -1, -1, -1, -1, -1, -1, -1, -1, & !   10
     113             :                                                                 -1, -1, 3, 3, 3, 3, 3, -1, & !   18
     114             :                                                                 -1, -1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, -1, & !   36
     115             :                                                                 -1, -1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, -1, & !   54
     116             :                                                                 -1, -1, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 5, &
     117             :                                                                 5, 5, 5, 5, 5, 5, 5, 5, 5, -1, -1, -1, -1, -1, -1, & !   86
     118             :                                                     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1/)
     119             : 
     120             :    INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE ::   nqf = (/-1, & !    0
     121             :                                                                 -1, -1, & !    2
     122             :                                                                 -1, -1, -1, -1, -1, -1, -1, -1, & !   10
     123             :                                                                 -1, -1, -1, -1, -1, -1, -1, -1, & !   18
     124             :                                                     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, & !   36
     125             :                                                     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, & !   54
     126             :                                                                -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, &
     127             :                                                                 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, & !   86
     128             :                                                     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1/)
     129             : 
     130             :    ! Element Valence
     131             :    INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: zval = (/-1, & !    0
     132             :                                                                1, 2, & !    2
     133             :                                                                1, 2, 3, 4, 5, 6, 7, 8, & !   10
     134             :                                                                1, 2, 3, 4, 5, 6, 7, 8, & !   18
     135             :                                                                1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 2, 3, 4, 5, 6, 7, 8, & !   36
     136             :                                                                1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 2, 3, 4, 5, 6, 7, 8, & !   54
     137             :                                                                1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 3, &
     138             :                                                                4, 5, 6, 7, 8, 9, 10, 11, 2, 3, 4, 5, 6, 7, -1, & !   86
     139             :                                                       -1, -1, -1, 4, -1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1/)
     140             : 
     141             :    ! Number of 1 center 2 electron integrals involving partially filled d shells
     142             :    ! r016:  <SS|DD>
     143             :    INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: ir016 = (/0, & !    0
     144             :                                                                 0, 0, & !    2
     145             :                                                                 0, 0, 0, 0, 0, 0, 0, 0, & !   10
     146             :                                                                 0, 0, 0, 0, 0, 0, 0, 0, & !   18
     147             :                                                                 0, 0, 2, 4, 6, 5, 10, 12, 14, 16, 10, 0, 0, 0, 0, 0, 0, 0, & !   36
     148             :                                                                 0, 0, 4, 4, 4, 5, 10, 7, 8, 0, 10, 0, 0, 0, 0, 0, 0, 0, & !   54
     149             :                                                                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
     150             :                                                                 4, 6, 8, 10, 12, 14, 9, 10, 0, 0, 0, 0, 0, 0, 0, & !   86
     151             :                                                                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/)
     152             : 
     153             :    ! r066:  <DD|DD> "0" term
     154             :    INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: ir066 = (/0, & !    0
     155             :                                                                 0, 0, & !    2
     156             :                                                                 0, 0, 0, 0, 0, 0, 0, 0, & !   10
     157             :                                                                 0, 0, 0, 0, 0, 0, 0, 0, & !   18
     158             :                                                                 0, 0, 0, 1, 3, 10, 10, 15, 21, 28, 45, 0, 0, 0, 0, 0, 0, 0, & !   36
     159             :                                                                 0, 0, 0, 1, 6, 10, 10, 21, 28, 45, 45, 0, 0, 0, 0, 0, 0, 0, & !   54
     160             :                                                                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
     161             :                                                                 1, 3, 6, 10, 15, 21, 36, 45, 0, 0, 0, 0, 0, 0, 0, & !   86
     162             :                                                                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/)
     163             : 
     164             :    ! r244:  <SD|SD>
     165             :    INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: ir244 = (/0, & !    0
     166             :                                                                 0, 0, & !    2
     167             :                                                                 0, 0, 0, 0, 0, 0, 0, 0, & !   10
     168             :                                                                 0, 0, 0, 0, 0, 0, 0, 0, & !   18
     169             :                                                                 0, 0, 1, 2, 3, 5, 5, 6, 7, 8, 5, 0, 0, 0, 0, 0, 0, 0, & !   36
     170             :                                                                 0, 0, 1, 2, 4, 5, 5, 5, 5, 0, 5, 0, 0, 0, 0, 0, 0, 0, & !   54
     171             :                                                                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
     172             :                                                                 2, 3, 4, 5, 6, 7, 5, 5, 0, 0, 0, 0, 0, 0, 0, & !   86
     173             :                                                                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/)
     174             : 
     175             :    ! r266:  <DD|DD> "2" term
     176             :    INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: ir266 = (/0, & !    0
     177             :                                                                 0, 0, & !    2
     178             :                                                                 0, 0, 0, 0, 0, 0, 0, 0, & !   10
     179             :                                                                 0, 0, 0, 0, 0, 0, 0, 0, & !   18
     180             :                                                                 0, 0, 0, 8, 15, 35, 35, 35, 43, 50, 70, 0, 0, 0, 0, 0, 0, 0, & !   36
     181             :                                                                 0, 0, 0, 8, 21, 35, 35, 43, 50, 70, 70, 0, 0, 0, 0, 0, 0, 0, & !   54
     182             :                                                                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
     183             :                                                                 8, 15, 21, 35, 35, 43, 56, 70, 0, 0, 0, 0, 0, 0, 0, & !   86
     184             :                                                                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/)
     185             : 
     186             :    ! r466:  <DD|DD> "4" term
     187             :    INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: ir466 = (/0, & !    0
     188             :                                                                 0, 0, & !    2
     189             :                                                                 0, 0, 0, 0, 0, 0, 0, 0, & !   10
     190             :                                                                 0, 0, 0, 0, 0, 0, 0, 0, & !   18
     191             :                                                                 0, 0, 0, 1, 8, 35, 35, 35, 36, 43, 70, 0, 0, 0, 0, 0, 0, 0, & !   36
     192             :                                                                 0, 0, 0, 1, 21, 35, 35, 36, 43, 70, 70, 0, 0, 0, 0, 0, 0, 0, & !   54
     193             :                                                                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
     194             :                                                                 1, 8, 21, 35, 35, 36, 56, 70, 0, 0, 0, 0, 0, 0, 0, & !   86
     195             :                                                                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/)
     196             : 
     197             :    INTERFACE amn_l
     198             :       MODULE PROCEDURE amn_l1, amn_l2
     199             :    END INTERFACE
     200             : 
     201             :    PUBLIC :: convert_param_to_cp2k, calpar, valence_electrons, get_se_basis, &
     202             :              setup_1c_2el_int, amn_l
     203             : 
     204             : CONTAINS
     205             : 
     206             : ! **************************************************************************************************
     207             : !> \brief  Gives back the number of valence electrons for element z and also the
     208             : !>         number of atomic orbitals for that specific element
     209             : !> \param sep ...
     210             : !> \param extended_basis_set ...
     211             : ! **************************************************************************************************
     212        3964 :    SUBROUTINE valence_electrons(sep, extended_basis_set)
     213             :       TYPE(semi_empirical_type), POINTER                 :: sep
     214             :       LOGICAL, INTENT(IN)                                :: extended_basis_set
     215             : 
     216             :       INTEGER                                            :: natorb, z
     217             :       LOGICAL                                            :: check, use_p_orbitals
     218             :       REAL(KIND=dp)                                      :: zeff
     219             : 
     220        3964 :       use_p_orbitals = .TRUE.
     221        3964 :       z = sep%z
     222        3964 :       CPASSERT(z >= 0)
     223             :       ! Special case for Hydrogen.. If requested allow p-orbitals on it..
     224        1768 :       SELECT CASE (z)
     225             :       CASE (0, 2)
     226        1768 :          use_p_orbitals = .FALSE.
     227             :       CASE (1)
     228        3964 :          use_p_orbitals = sep%p_orbitals_on_h
     229             :       CASE DEFAULT
     230             :          ! Nothing to do..
     231             :       END SELECT
     232             :       ! Determine the number of atomic orbitals
     233        3964 :       natorb = 0
     234        3964 :       IF (nqs(z) > 0) natorb = natorb + 1
     235        3964 :       IF ((nqp(z) > 0) .OR. use_p_orbitals) natorb = natorb + 3
     236        3964 :       IF (extended_basis_set .AND. element_has_d(sep)) natorb = natorb + 5
     237        3964 :       IF (extended_basis_set .AND. element_has_f(sep)) natorb = natorb + 7
     238             :       ! Check and assignment
     239        3964 :       check = (natorb <= 4) .OR. (extended_basis_set)
     240           0 :       CPASSERT(check)
     241        3964 :       sep%natorb = natorb
     242        3964 :       sep%extended_basis_set = extended_basis_set
     243             :       ! Determine the Z eff
     244        3964 :       zeff = REAL(zval(z), KIND=dp)
     245        3964 :       sep%zeff = zeff
     246        3964 :    END SUBROUTINE valence_electrons
     247             : 
     248             : ! **************************************************************************************************
     249             : !> \brief  Gives back the number of basis function for each l
     250             : !> \param sep ...
     251             : !> \param l ...
     252             : !> \return ...
     253             : ! **************************************************************************************************
     254        4750 :    FUNCTION get_se_basis(sep, l) RESULT(n)
     255             :       TYPE(semi_empirical_type), POINTER                 :: sep
     256             :       INTEGER, INTENT(IN)                                :: l
     257             :       INTEGER                                            :: n
     258             : 
     259        4750 :       IF (sep%z < 0 .OR. sep%z > nelem) THEN
     260           0 :          CPABORT("Invalid atomic number !")
     261             :       ELSE
     262        4750 :          IF (l == 0) THEN
     263        2240 :             n = nqs(sep%z)
     264        2510 :          ELSEIF (l == 1) THEN
     265             :             ! Special case for Hydrogen.. If requested allow p-orbitals on it..
     266        1772 :             IF ((sep%z == 1) .AND. sep%p_orbitals_on_h) THEN
     267             :                n = 1
     268             :             ELSE
     269        1758 :                n = nqp(sep%z)
     270             :             END IF
     271         738 :          ELSEIF (l == 2) THEN
     272         738 :             n = nqd(sep%z)
     273           0 :          ELSEIF (l == 3) THEN
     274           0 :             n = nqf(sep%z)
     275             :          ELSE
     276           0 :             CPABORT("Invalid l quantum number !")
     277             :          END IF
     278        4736 :          IF (n < 0) THEN
     279           0 :             CPABORT("Invalid n quantum number !")
     280             :          END IF
     281             :       END IF
     282        4750 :    END FUNCTION get_se_basis
     283             : 
     284             : ! **************************************************************************************************
     285             : !> \brief  Converts parameter units to internal
     286             : !> \param sep ...
     287             : !> \date   03.2008 [tlaino]
     288             : !> \author Teodoro Laino [tlaino] - University of Zurich
     289             : ! **************************************************************************************************
     290        2240 :    SUBROUTINE convert_param_to_cp2k(sep)
     291             :       TYPE(semi_empirical_type), POINTER                 :: sep
     292             : 
     293       11200 :       sep%beta = sep%beta/evolt
     294        2240 :       sep%uss = sep%uss/evolt
     295        2240 :       sep%upp = sep%upp/evolt
     296        2240 :       sep%udd = sep%udd/evolt
     297        2240 :       sep%alp = sep%alp/bohr
     298        2240 :       sep%eisol = sep%eisol/evolt
     299        2240 :       sep%gss = sep%gss/evolt
     300        2240 :       sep%gsp = sep%gsp/evolt
     301        2240 :       sep%gpp = sep%gpp/evolt
     302        2240 :       sep%gp2 = sep%gp2/evolt
     303        2240 :       sep%gsd = sep%gsd/evolt
     304        2240 :       sep%gpd = sep%gpd/evolt
     305        2240 :       sep%gdd = sep%gdd/evolt
     306        2240 :       sep%hsp = sep%hsp/evolt
     307       11200 :       sep%fn1 = sep%fn1*bohr/evolt
     308       11200 :       sep%fn2 = sep%fn2/bohr/bohr
     309       11200 :       sep%fn3 = sep%fn3*bohr
     310       47040 :       sep%bfn1 = sep%bfn1*bohr/evolt
     311       47040 :       sep%bfn2 = sep%bfn2/bohr/bohr
     312       47040 :       sep%bfn3 = sep%bfn3*bohr
     313        2240 :       sep%f0sd = sep%f0sd
     314             :       sep%g2sd = sep%g2sd
     315        2240 :       sep%a = sep%a*bohr/evolt
     316        2240 :       sep%b = sep%b/bohr/bohr
     317        2240 :       sep%c = sep%c*bohr
     318        6720 :       sep%pre = sep%pre/evolt
     319        6720 :       sep%d = sep%d/bohr
     320             : 
     321        2240 :    END SUBROUTINE convert_param_to_cp2k
     322             : 
     323             : ! **************************************************************************************************
     324             : !> \brief  Calculates missing parameters
     325             : !> \param z ...
     326             : !> \param sep ...
     327             : !> \date   03.2008 [tlaino]
     328             : !> \author Teodoro Laino [tlaino] - University of Zurich
     329             : ! **************************************************************************************************
     330        2240 :    SUBROUTINE calpar(z, sep)
     331             :       INTEGER                                            :: z
     332             :       TYPE(semi_empirical_type), POINTER                 :: sep
     333             : 
     334             :       INTEGER                                            :: iod, iop, ios, j, jmax, k, l
     335             :       REAL(KIND=dp) :: ad, am, aq, d1, d2, d3, dd, df, eisol, gdd1, gp2, gp2c, gpp, gppc, gqq, &
     336             :          gsp, gspc, gss, gssc, hpp, hpp1, hpp2, hsp, hsp1, hsp2, hspc, p, p4, q1, q2, q3, qf, qn, &
     337             :          qq, udd, upp, uss, zp, zs
     338             : 
     339        2240 :       IF (.NOT. sep%defined) RETURN
     340        2240 :       uss = sep%uss
     341        2240 :       upp = sep%upp
     342        2240 :       udd = sep%udd
     343        2240 :       gss = sep%gss
     344        2240 :       gpp = sep%gpp
     345        2240 :       gsp = sep%gsp
     346        2240 :       gp2 = sep%gp2
     347        2240 :       hsp = sep%hsp
     348        2240 :       zs = sep%sto_exponents(0)
     349        2240 :       zp = sep%sto_exponents(1)
     350        2240 :       ios = Nos(z)
     351        2240 :       iop = Nop(z)
     352        2240 :       iod = Nod(z)
     353             : 
     354        2240 :       p = 2.0_dp
     355        2240 :       p4 = p**4
     356             :       !  GSSC is the number of two-electron terms of type <SS|SS>
     357        2240 :       gssc = REAL(MAX(ios - 1, 0), KIND=dp)
     358        2240 :       k = iop
     359             :       !  GSPC is the number of two-electron terms of type <SS|PP>
     360        2240 :       gspc = REAL(ios*k, KIND=dp)
     361        2240 :       l = MIN(k, 6 - k)
     362             :       !  GP2C is the number of two-electron terms of type <PP|PP>
     363             :       !       plus 0.5 of the number of HPP integrals.
     364             :       !  (HPP is not used; instead it is replaced by 0.5(GPP-GP2))
     365        2240 :       gp2c = REAL((k*(k - 1))/2, KIND=dp) + 0.5_dp*REAL((l*(l - 1))/2, KIND=dp)
     366             :       !  GPPC is minus 0.5 times the number of HPP integrals.
     367        2240 :       gppc = -0.5_dp*REAL((l*(l - 1))/2, KIND=dp)
     368             :       !  HSPC is the number of two-electron terms of type <SP|SP>.
     369             :       !       (S and P must have the same spin.  In all cases, if
     370             :       !  P is non-zero, there are two S electrons)
     371        2240 :       hspc = REAL(-k, KIND=dp)
     372             :       !  Constraint the value of the STO exponent
     373        2240 :       zp = MAX(0.3_dp, zp)
     374             :       !  Take into account constraints on the values of the integrals
     375        2240 :       hpp = 0.5_dp*(gpp - gp2)
     376        2240 :       hpp = MAX(0.1_dp, hpp)
     377        2240 :       hsp = MAX(1.E-7_dp, hsp)
     378             : 
     379             :       ! Evaluation of EISOL
     380        2240 :       eisol = uss*ios + upp*iop + udd*iod + gss*gssc + gpp*gppc + gsp*gspc + gp2*gp2c + hsp*hspc
     381             : 
     382             :       ! Principal quantum number
     383        2240 :       qn = REAL(nqs(z), KIND=dp)
     384        2240 :       CPASSERT(qn > 0)
     385             : 
     386             :       ! Charge separation evaluation
     387        2240 :       dd = (2.0_dp*qn + 1)*(4.0_dp*zs*zp)**(qn + 0.5_dp)/(zs + zp)**(2.0_dp*qn + 2)/SQRT(3.0_dp)
     388        2240 :       qq = SQRT((4.0_dp*qn*qn + 6.0_dp*qn + 2.0_dp)/20.0_dp)/zp
     389             : 
     390             :       ! Calculation of the additive terms in atomic units
     391        2240 :       jmax = 5
     392        2240 :       gdd1 = (hsp/(evolt*dd**2))**(1.0_dp/3.0_dp)
     393        2240 :       d1 = gdd1
     394        2240 :       d2 = gdd1 + 0.04_dp
     395       12126 :       DO j = 1, jmax
     396       10324 :          df = d2 - d1
     397       10324 :          hsp1 = 0.5_dp*d1 - 0.5_dp/SQRT(4.0_dp*dd**2 + 1.0_dp/d1**2)
     398       10324 :          hsp2 = 0.5_dp*d2 - 0.5_dp/SQRT(4.0_dp*dd**2 + 1.0_dp/d2**2)
     399       10324 :          IF (ABS(hsp2 - hsp1) < EPSILON(0.0_dp)) EXIT
     400        9886 :          d3 = d1 + df*(hsp/evolt - hsp1)/(hsp2 - hsp1)
     401        9886 :          d1 = d2
     402       12126 :          d2 = d3
     403             :       END DO
     404        2240 :       gqq = (p4*hpp/(evolt*48.0_dp*qq**4))**0.2_dp
     405        2240 :       q1 = gqq
     406        2240 :       q2 = gqq + 0.04_dp
     407       13440 :       DO j = 1, jmax
     408       11200 :          qf = q2 - q1
     409       11200 :          hpp1 = 0.25_dp*q1 - 0.5_dp/SQRT(4.0_dp*qq**2 + 1.0_dp/q1**2) + 0.25_dp/SQRT(8.0_dp*qq**2 + 1.0_dp/q1**2)
     410       11200 :          hpp2 = 0.25_dp*q2 - 0.5_dp/SQRT(4.0_dp*qq**2 + 1.0_dp/q2**2) + 0.25_dp/SQRT(8.0_dp*qq**2 + 1.0_dp/q2**2)
     411       11200 :          IF (ABS(hpp2 - hpp1) < EPSILON(0.0_dp)) EXIT
     412       11200 :          q3 = q1 + qf*(hpp/evolt - hpp1)/(hpp2 - hpp1)
     413       11200 :          q1 = q2
     414       13440 :          q2 = q3
     415             :       END DO
     416        2240 :       am = gss/evolt
     417        2240 :       ad = d2
     418        2240 :       aq = q2
     419        2240 :       IF (z == 1) THEN
     420         438 :          ad = am
     421         438 :          aq = am
     422         438 :          dd = 0.0_dp
     423         438 :          qq = 0.0_dp
     424             :       END IF
     425             :       ! Overwrite these parameters if they were undefined.. otherwise keep the defined
     426             :       ! value
     427        2240 :       IF (ABS(sep%eisol) < EPSILON(0.0_dp)) sep%eisol = eisol
     428        2240 :       IF (ABS(sep%dd) < EPSILON(0.0_dp)) sep%dd = dd
     429        2240 :       IF (ABS(sep%qq) < EPSILON(0.0_dp)) sep%qq = qq
     430        2240 :       IF (ABS(sep%am) < EPSILON(0.0_dp)) sep%am = am
     431        2240 :       IF (ABS(sep%ad) < EPSILON(0.0_dp)) sep%ad = ad
     432        2240 :       IF (ABS(sep%aq) < EPSILON(0.0_dp)) sep%aq = aq
     433             :       ! Proceed with d-orbitals and fill the Kolpman-Ohno and Charge Separation
     434             :       ! arrays
     435        2240 :       CALL calpar_d(sep)
     436             :    END SUBROUTINE calpar
     437             : 
     438             : ! **************************************************************************************************
     439             : !> \brief  Finalize the initialization of parameters, defining additional
     440             : !>         parameters for d-orbitals
     441             : !>
     442             : !> \param sep ...
     443             : !> \date   03.2008 [tlaino]
     444             : !> \author Teodoro Laino [tlaino] - University of Zurich
     445             : ! **************************************************************************************************
     446        2240 :    SUBROUTINE calpar_d(sep)
     447             :       TYPE(semi_empirical_type), POINTER                 :: sep
     448             : 
     449             :       REAL(KIND=dp), DIMENSION(6)                        :: amn
     450             : 
     451             : ! Determine if this element owns d-orbitals (only if the parametrization
     452             : ! supports the d-orbitals)
     453             : 
     454        2240 :       IF (sep%extended_basis_set) sep%dorb = element_has_d(sep)
     455        2240 :       IF (sep%dorb) THEN
     456         738 :          CALL amn_l(sep, amn)
     457         738 :          CALL eval_1c_2el_spd(sep)
     458         738 :          CALL eval_cs_ko(sep, amn)
     459             :       END IF
     460        2240 :       IF (.NOT. sep%dorb) THEN
     461             :          ! Use the old integral module
     462        1502 :          IF (ABS(sep%am) > EPSILON(0.0_dp)) THEN
     463        1502 :             sep%ko(1) = 0.5_dp/sep%am
     464             :          END IF
     465        1502 :          IF (ABS(sep%ad) > EPSILON(0.0_dp) .AND. (sep%z /= 1)) THEN
     466        1064 :             sep%ko(2) = 0.5_dp/sep%ad
     467             :          END IF
     468        1502 :          IF (ABS(sep%aq) > EPSILON(0.0_dp) .AND. (sep%z /= 1)) THEN
     469        1064 :             sep%ko(3) = 0.5_dp/sep%aq
     470             :          END IF
     471        1502 :          sep%ko(7) = sep%ko(1)
     472        1502 :          sep%ko(9) = sep%ko(1)
     473        1502 :          sep%cs(2) = sep%dd
     474        1502 :          sep%cs(3) = sep%qq*SQRT(2.0_dp)
     475             :       ELSE
     476             :          ! Use the new integral module
     477         738 :          sep%ko(9) = sep%ko(1)
     478         738 :          sep%aq = 0.5_dp/sep%ko(3)
     479             :       END IF
     480             :       ! In case the Klopman-Ohno CORE therm is provided let's overwrite the
     481             :       ! computed one
     482        2240 :       IF (ABS(sep%rho) > EPSILON(0.0_dp)) THEN
     483         130 :          sep%ko(9) = sep%rho
     484             :       END IF
     485        2240 :    END SUBROUTINE calpar_d
     486             : 
     487             : ! **************************************************************************************************
     488             : !> \brief  Determines if the elements has d-orbitals
     489             : !>
     490             : !> \param sep ...
     491             : !> \return ...
     492             : !> \date   05.2008 [tlaino]
     493             : !> \author Teodoro Laino [tlaino] - University of Zurich
     494             : ! **************************************************************************************************
     495        3464 :    FUNCTION element_has_d(sep) RESULT(res)
     496             :       TYPE(semi_empirical_type), POINTER                 :: sep
     497             :       LOGICAL                                            :: res
     498             : 
     499        3464 :       res = (nqd(sep%z) > 0 .AND. sep%sto_exponents(2) > EPSILON(0.0_dp))
     500        3464 :    END FUNCTION element_has_d
     501             : 
     502             : ! **************************************************************************************************
     503             : !> \brief  Determines if the elements has f-orbitals
     504             : !>
     505             : !> \param sep ...
     506             : !> \return ...
     507             : !> \date   05.2008 [tlaino]
     508             : !> \author Teodoro Laino [tlaino] - University of Zurich
     509             : ! **************************************************************************************************
     510        1732 :    FUNCTION element_has_f(sep) RESULT(res)
     511             :       TYPE(semi_empirical_type), POINTER                 :: sep
     512             :       LOGICAL                                            :: res
     513             : 
     514        1732 :       res = (nqf(sep%z) > 0 .AND. sep%sto_exponents(3) > EPSILON(0.0_dp))
     515        1732 :    END FUNCTION element_has_f
     516             : 
     517             : ! **************************************************************************************************
     518             : !> \brief  Computes the A^{\mu \nu}_l values for the evaluation of the two-center
     519             : !>          two-electron integrals. The term is the one reported in Eq.(7) of TCA
     520             : !>
     521             : !> \param sep ...
     522             : !> \param amn ...
     523             : !> \date   03.2008 [tlaino]
     524             : !> \par    Notation Index: 1 (SS), 2 (SP), 3 (SD), 4 (PP), 5 (PD), 6 (DD)
     525             : !>
     526             : !> \author Teodoro Laino [tlaino] - University of Zurich
     527             : ! **************************************************************************************************
     528         738 :    SUBROUTINE amn_l1(sep, amn)
     529             :       TYPE(semi_empirical_type), POINTER                 :: sep
     530             :       REAL(KIND=dp), DIMENSION(6), INTENT(OUT)           :: amn
     531             : 
     532             :       INTEGER                                            :: nd, nsp
     533             :       REAL(KIND=dp)                                      :: z1, z2, z3
     534             : 
     535         738 :       z1 = sep%sto_exponents(0)
     536         738 :       z2 = sep%sto_exponents(1)
     537         738 :       z3 = sep%sto_exponents(2)
     538         738 :       IF (z1 <= 0.0_dp) &
     539             :          CALL cp_abort(__LOCATION__, &
     540             :                        "Trying to use s-orbitals, but the STO exponents is set to 0. "// &
     541           0 :                        "Please check if your parameterization supports the usage of s orbitals! ")
     542         738 :       amn = 0.0_dp
     543         738 :       nsp = nqs(sep%z)
     544         738 :       IF (sep%natorb >= 4) THEN
     545         738 :          IF (z2 <= 0.0_dp) &
     546             :             CALL cp_abort(__LOCATION__, &
     547             :                           "Trying to use p-orbitals, but the STO exponents is set to 0. "// &
     548           0 :                           "Please check if your parameterization supports the usage of p orbitals! ")
     549         738 :          amn(2) = amn_l_low(z1, z2, nsp, nsp, 1)
     550         738 :          amn(3) = amn_l_low(z2, z2, nsp, nsp, 2)
     551         738 :          IF (sep%dorb) THEN
     552         738 :             IF (z3 <= 0.0_dp) &
     553             :                CALL cp_abort(__LOCATION__, &
     554             :                              "Trying to use d-orbitals, but the STO exponents is set to 0. "// &
     555           0 :                              "Please check if your parameterization supports the usage of d orbitals! ")
     556         738 :             nd = nqd(sep%z)
     557         738 :             amn(4) = amn_l_low(z1, z3, nsp, nd, 2)
     558         738 :             amn(5) = amn_l_low(z2, z3, nsp, nd, 1)
     559         738 :             amn(6) = amn_l_low(z3, z3, nd, nd, 2)
     560             :          END IF
     561             :       END IF
     562         738 :    END SUBROUTINE amn_l1
     563             : 
     564             : ! **************************************************************************************************
     565             : !> \brief  Computes the A^{\mu \nu}_l values for the evaluation of the two-center
     566             : !>          two-electron integrals. The term is the one reported in Eq.(7) of TCA
     567             : !>
     568             : !> \param sep ...
     569             : !> \param amn ...
     570             : !> \date   09.2008 [tlaino]
     571             : !> \par
     572             : !>
     573             : !> \author Teodoro Laino [tlaino] - University of Zurich
     574             : ! **************************************************************************************************
     575        2240 :    SUBROUTINE amn_l2(sep, amn)
     576             :       TYPE(semi_empirical_type), POINTER                 :: sep
     577             :       REAL(KIND=dp), DIMENSION(6, 0:2), INTENT(OUT)      :: amn
     578             : 
     579             :       INTEGER                                            :: nd, nsp
     580             :       REAL(KIND=dp)                                      :: z1, z2, z3
     581             : 
     582        2240 :       z1 = sep%sto_exponents(0)
     583        2240 :       z2 = sep%sto_exponents(1)
     584        2240 :       z3 = sep%sto_exponents(2)
     585        2240 :       CPASSERT(z1 > 0.0_dp)
     586        2240 :       amn = 0.0_dp
     587        2240 :       nsp = nqs(sep%z)
     588        2240 :       amn(1, 0) = amn_l_low(z1, z1, nsp, nsp, 0)
     589        2240 :       IF (sep%natorb >= 4) THEN
     590        1772 :          CPASSERT(z2 > 0.0_dp)
     591        1772 :          amn(2, 1) = amn_l_low(z1, z2, nsp, nsp, 1)
     592        1772 :          amn(3, 0) = amn_l_low(z2, z2, nsp, nsp, 0)
     593        1772 :          amn(3, 2) = amn_l_low(z2, z2, nsp, nsp, 2)
     594        1772 :          IF (sep%dorb) THEN
     595         738 :             CPASSERT(z3 > 0.0_dp)
     596         738 :             nd = nqd(sep%z)
     597         738 :             amn(4, 2) = amn_l_low(z1, z3, nsp, nd, 2)
     598         738 :             amn(5, 1) = amn_l_low(z2, z3, nsp, nd, 1)
     599         738 :             amn(6, 0) = amn_l_low(z3, z3, nd, nd, 0)
     600         738 :             amn(6, 2) = amn_l_low(z3, z3, nd, nd, 2)
     601             :          END IF
     602             :       END IF
     603        2240 :    END SUBROUTINE amn_l2
     604             : 
     605             : ! **************************************************************************************************
     606             : !> \brief  Low level for computing Eq.(7) of TCA
     607             : !> \param z1 ...
     608             : !> \param z2 ...
     609             : !> \param n1 ...
     610             : !> \param n2 ...
     611             : !> \param l ...
     612             : !> \return ...
     613             : !> \date   03.2008 [tlaino]
     614             : !> \author Teodoro Laino [tlaino] - University of Zurich
     615             : ! **************************************************************************************************
     616       14198 :    FUNCTION amn_l_low(z1, z2, n1, n2, l) RESULT(amnl)
     617             :       REAL(KIND=dp), INTENT(IN)                          :: z1, z2
     618             :       INTEGER, INTENT(IN)                                :: n1, n2, l
     619             :       REAL(KIND=dp)                                      :: amnl
     620             : 
     621             :       amnl = fac(n1 + n2 + l)/SQRT(fac(2*n1)*fac(2*n2))*(2.0_dp*z1/(z1 + z2))**n1* &
     622       14198 :              (2.0_dp*z2/(z1 + z2))**n2*2.0_dp*SQRT(z1*z2)/(z1 + z2)**(l + 1)
     623             : 
     624       14198 :    END FUNCTION amn_l_low
     625             : 
     626             : ! **************************************************************************************************
     627             : !> \brief  Calculation of chare separations and additive terms used for computing
     628             : !>         the two-center two-electron integrals with d-orbitals
     629             : !> \param sep ...
     630             : !> \param amn ...
     631             : !> \date   03.2008 [tlaino]
     632             : !> \par    Notation
     633             : !>         -) Charge separations [sep%cs(1:6)]  [see equations (12)-(16) of TCA]
     634             : !>         -) Additive terms of Klopman-Ohno terms [sep%ko(1:9)] [see equations
     635             : !>            (19)-(26) of TCA]
     636             : !>         -) Atomic core additive term stored in sep%ko(9): used in the calculation
     637             : !>            of the core-electron attractions and core-core repulsions
     638             : !> \author Teodoro Laino [tlaino] - University of Zurich
     639             : ! **************************************************************************************************
     640         738 :    SUBROUTINE eval_cs_ko(sep, amn)
     641             :       TYPE(semi_empirical_type), POINTER                 :: sep
     642             :       REAL(KIND=dp), DIMENSION(6), INTENT(IN)            :: amn
     643             : 
     644             :       REAL(KIND=dp)                                      :: d, fg
     645             : 
     646             : ! SS term
     647             : 
     648         738 :       fg = sep%gss
     649         738 :       sep%ko(1) = ko_ij(0, 1.0_dp, fg)
     650         738 :       IF (sep%natorb >= 4) THEN
     651             :          ! Other terms for SP basis
     652             :          ! SP
     653         738 :          d = amn(2)/SQRT(3.0_dp)
     654         738 :          fg = sep%hsp
     655         738 :          sep%cs(2) = d
     656         738 :          sep%ko(2) = ko_ij(1, d, fg)
     657             :          ! PP
     658         738 :          sep%ko(7) = sep%ko(1)
     659         738 :          d = SQRT(amn(3)*2.0_dp/5.0_dp)
     660         738 :          fg = 0.5_dp*(sep%gpp - sep%gp2)
     661         738 :          sep%cs(3) = d
     662         738 :          sep%ko(3) = ko_ij(2, d, fg)
     663             :          ! Terms involving d-orbitals
     664         738 :          IF (sep%dorb) THEN
     665             :             ! SD
     666         738 :             d = SQRT(amn(4)*2.0_dp/SQRT(15.0_dp))
     667         738 :             fg = sep%onec2el(19)
     668         738 :             sep%cs(4) = d
     669         738 :             sep%ko(4) = ko_ij(2, d, fg)
     670             :             ! PD
     671         738 :             d = amn(5)/SQRT(5.0_dp)
     672         738 :             fg = sep%onec2el(23) - 1.8_dp*sep%onec2el(35)
     673         738 :             sep%cs(5) = d
     674         738 :             sep%ko(5) = ko_ij(1, d, fg)
     675             :             ! DD
     676         738 :             fg = 0.2_dp*(sep%onec2el(29) + 2.0_dp*sep%onec2el(30) + 2.0_dp*sep%onec2el(31))
     677         738 :             sep%ko(8) = ko_ij(0, 1.0_dp, fg)
     678         738 :             d = SQRT(amn(6)*2.0_dp/7.0_dp)
     679         738 :             fg = sep%onec2el(44) - (20.0_dp/35.0_dp)*sep%onec2el(52)
     680         738 :             sep%cs(6) = d
     681         738 :             sep%ko(6) = ko_ij(2, d, fg)
     682             :          END IF
     683             :       END IF
     684         738 :    END SUBROUTINE eval_cs_ko
     685             : 
     686             : ! **************************************************************************************************
     687             : !> \brief  Computes the 1 center two-electrons integrals for a SPD basis
     688             : !> \param sep ...
     689             : !> \date   03.2008 [tlaino]
     690             : !> \author Teodoro Laino [tlaino] - University of Zurich
     691             : ! **************************************************************************************************
     692         738 :    SUBROUTINE eval_1c_2el_spd(sep)
     693             :       TYPE(semi_empirical_type), POINTER                 :: sep
     694             : 
     695             :       REAL(KIND=dp)                                      :: r016, r036, r066, r125, r155, r234, &
     696             :                                                             r236, r244, r246, r266, r355, r466, &
     697             :                                                             s15, s3, s5
     698             : 
     699         738 :       IF (sep%dorb) THEN
     700         738 :          s3 = SQRT(3.0_dp)
     701         738 :          s5 = SQRT(5.0_dp)
     702         738 :          s15 = SQRT(15.0_dp)
     703             : 
     704             :          ! We evaluate now the Slater-Condon parameters (Rlij)
     705             :          CALL sc_param(sep, r066, r266, r466, r016, r244, r036, r236, r155, r355, r125, &
     706         738 :                        r234, r246)
     707             : 
     708         738 :          IF (ABS(sep%f0sd) > EPSILON(0.0_dp)) THEN
     709         436 :             r016 = sep%f0sd
     710             :          END IF
     711         738 :          IF (ABS(sep%g2sd) > EPSILON(0.0_dp)) THEN
     712         436 :             r244 = sep%g2sd
     713             :          END IF
     714         738 :          CALL eisol_corr(sep, r016, r066, r244, r266, r466)
     715         738 :          sep%onec2el(1) = r016
     716         738 :          sep%onec2el(2) = 2.0_dp/(3.0_dp*s5)*r125
     717         738 :          sep%onec2el(3) = 1.0_dp/s15*r125
     718         738 :          sep%onec2el(4) = 2.0_dp/(5.0_dp*s5)*r234
     719         738 :          sep%onec2el(5) = r036 + 4.0_dp/35.0_dp*r236
     720         738 :          sep%onec2el(6) = r036 + 2.0_dp/35.0_dp*r236
     721         738 :          sep%onec2el(7) = r036 - 4.0_dp/35.0_dp*r236
     722         738 :          sep%onec2el(8) = -1.0_dp/(3.0_dp*s5)*r125
     723         738 :          sep%onec2el(9) = SQRT(3.0_dp/125.0_dp)*r234
     724         738 :          sep%onec2el(10) = s3/35.0_dp*r236
     725         738 :          sep%onec2el(11) = 3.0_dp/35.0_dp*r236
     726         738 :          sep%onec2el(12) = -1.0_dp/(5.0_dp*s5)*r234
     727         738 :          sep%onec2el(13) = r036 - 2.0_dp/35.0_dp*r236
     728         738 :          sep%onec2el(14) = -2.0_dp*s3/35.0_dp*r236
     729         738 :          sep%onec2el(15) = -sep%onec2el(3)
     730         738 :          sep%onec2el(16) = -sep%onec2el(11)
     731         738 :          sep%onec2el(17) = -sep%onec2el(9)
     732         738 :          sep%onec2el(18) = -sep%onec2el(14)
     733         738 :          sep%onec2el(19) = 1.0_dp/5.0_dp*r244
     734         738 :          sep%onec2el(20) = 2.0_dp/(7.0_dp*s5)*r246
     735         738 :          sep%onec2el(21) = sep%onec2el(20)/2.0_dp
     736         738 :          sep%onec2el(22) = -sep%onec2el(20)
     737         738 :          sep%onec2el(23) = 4.0_dp/15.0_dp*r155 + 27.0_dp/245.0_dp*r355
     738         738 :          sep%onec2el(24) = 2.0_dp*s3/15.0_dp*r155 - 9.0_dp*s3/245.0_dp*r355
     739         738 :          sep%onec2el(25) = 1.0_dp/15.0_dp*r155 + 18.0_dp/245.0_dp*r355
     740         738 :          sep%onec2el(26) = -s3/15.0_dp*r155 + 12.0_dp*s3/245.0_dp*r355
     741         738 :          sep%onec2el(27) = -s3/15.0_dp*r155 - 3.0_dp*s3/245.0_dp*r355
     742         738 :          sep%onec2el(28) = -sep%onec2el(27)
     743         738 :          sep%onec2el(29) = r066 + 4.0_dp/49.0_dp*r266 + 4.0_dp/49.0_dp*r466
     744         738 :          sep%onec2el(30) = r066 + 2.0_dp/49.0_dp*r266 - 24.0_dp/441.0_dp*r466
     745         738 :          sep%onec2el(31) = r066 - 4.0_dp/49.0_dp*r266 + 6.0_dp/441.0_dp*r466
     746         738 :          sep%onec2el(32) = SQRT(3.0_dp/245.0_dp)*r246
     747         738 :          sep%onec2el(33) = 1.0_dp/5.0_dp*r155 + 24.0_dp/245.0_dp*r355
     748         738 :          sep%onec2el(34) = 1.0_dp/5.0_dp*r155 - 6.0_dp/245.0_dp*r355
     749         738 :          sep%onec2el(35) = 3.0_dp/49.0_dp*r355
     750         738 :          sep%onec2el(36) = 1.0_dp/49.0_dp*r266 + 30.0_dp/441.0_dp*r466
     751         738 :          sep%onec2el(37) = s3/49.0_dp*r266 - 5.0_dp*s3/441.0_dp*r466
     752         738 :          sep%onec2el(38) = r066 - 2.0_dp/49.0_dp*r266 - 4.0_dp/441.0_dp*r466
     753         738 :          sep%onec2el(39) = -2.0_dp*s3/49.0_dp*r266 + 10.0_dp*s3/441.0_dp*r466
     754         738 :          sep%onec2el(40) = -sep%onec2el(32)
     755         738 :          sep%onec2el(41) = -sep%onec2el(34)
     756         738 :          sep%onec2el(42) = -sep%onec2el(35)
     757         738 :          sep%onec2el(43) = -sep%onec2el(37)
     758         738 :          sep%onec2el(44) = 3.0_dp/49.0_dp*r266 + 20.0_dp/441.0_dp*r466
     759         738 :          sep%onec2el(45) = -sep%onec2el(39)
     760         738 :          sep%onec2el(46) = 1.0_dp/5.0_dp*r155 - 3.0_dp/35.0_dp*r355
     761         738 :          sep%onec2el(47) = -sep%onec2el(46)
     762         738 :          sep%onec2el(48) = 4.0_dp/49.0_dp*r266 + 15.0_dp/441.0_dp*r466
     763         738 :          sep%onec2el(49) = 3.0_dp/49.0_dp*r266 - 5.0_dp/147.0_dp*r466
     764         738 :          sep%onec2el(50) = -sep%onec2el(49)
     765         738 :          sep%onec2el(51) = r066 + 4.0_dp/49.0_dp*r266 - 34.0_dp/441.0_dp*r466
     766         738 :          sep%onec2el(52) = 35.0_dp/441.0_dp*r466
     767         738 :          sep%f0dd = r066
     768         738 :          sep%f2dd = r266
     769         738 :          sep%f4dd = r466
     770         738 :          sep%f0sd = r016
     771         738 :          sep%g2sd = r244
     772         738 :          sep%f0pd = r036
     773         738 :          sep%f2pd = r236
     774         738 :          sep%g1pd = r155
     775         738 :          sep%g3pd = r355
     776             :       END IF
     777         738 :    END SUBROUTINE eval_1c_2el_spd
     778             : 
     779             : ! **************************************************************************************************
     780             : !> \brief  Slater-Condon parameters for 1 center 2 electrons integrals
     781             : !> \param sep ...
     782             : !> \param r066 ...
     783             : !> \param r266 ...
     784             : !> \param r466 ...
     785             : !> \param r016 ...
     786             : !> \param r244 ...
     787             : !> \param r036 ...
     788             : !> \param r236 ...
     789             : !> \param r155 ...
     790             : !> \param r355 ...
     791             : !> \param r125 ...
     792             : !> \param r234 ...
     793             : !> \param r246 ...
     794             : !> \date   03.2008 [tlaino]
     795             : !> \author Teodoro Laino [tlaino] - University of Zurich
     796             : ! **************************************************************************************************
     797         738 :    SUBROUTINE sc_param(sep, r066, r266, r466, r016, r244, r036, r236, r155, r355, &
     798             :                        r125, r234, r246)
     799             :       TYPE(semi_empirical_type), POINTER                 :: sep
     800             :       REAL(KIND=dp), INTENT(out)                         :: r066, r266, r466, r016, r244, r036, &
     801             :                                                             r236, r155, r355, r125, r234, r246
     802             : 
     803             :       INTEGER                                            :: nd, ns
     804             :       REAL(KIND=dp)                                      :: ed, ep, es
     805             : 
     806         738 :       ns = nqs(sep%z)
     807         738 :       nd = nqd(sep%z)
     808         738 :       es = sep%zn(0)
     809         738 :       ep = sep%zn(1)
     810         738 :       ed = sep%zn(2)
     811         738 :       r016 = sc_param_low(0, ns, es, ns, es, nd, ed, nd, ed)
     812         738 :       r036 = sc_param_low(0, ns, ep, ns, ep, nd, ed, nd, ed)
     813         738 :       r066 = sc_param_low(0, nd, ed, nd, ed, nd, ed, nd, ed)
     814         738 :       r155 = sc_param_low(1, ns, ep, nd, ed, ns, ep, nd, ed)
     815         738 :       r125 = sc_param_low(1, ns, es, ns, ep, ns, ep, nd, ed)
     816         738 :       r244 = sc_param_low(2, ns, es, nd, ed, ns, es, nd, ed)
     817         738 :       r236 = sc_param_low(2, ns, ep, ns, ep, nd, ed, nd, ed)
     818         738 :       r266 = sc_param_low(2, nd, ed, nd, ed, nd, ed, nd, ed)
     819         738 :       r234 = sc_param_low(2, ns, ep, ns, ep, ns, es, nd, ed)
     820         738 :       r246 = sc_param_low(2, ns, es, nd, ed, nd, ed, nd, ed)
     821         738 :       r355 = sc_param_low(3, ns, ep, nd, ed, ns, ep, nd, ed)
     822         738 :       r466 = sc_param_low(4, nd, ed, nd, ed, nd, ed, nd, ed)
     823         738 :    END SUBROUTINE sc_param
     824             : 
     825             : ! **************************************************************************************************
     826             : !> \brief  Slater-Condon parameters for 1 center 2 electrons integrals - Low level
     827             : !> \param k ...
     828             : !> \param na ...
     829             : !> \param ea ...
     830             : !> \param nb ...
     831             : !> \param eb ...
     832             : !> \param nc ...
     833             : !> \param ec ...
     834             : !> \param nd ...
     835             : !> \param ed ...
     836             : !> \return ...
     837             : !> \date   03.2008 [tlaino]
     838             : !> \par    Notation
     839             : !>         -) k:      Type of integral
     840             : !>         -) na,na:  Principle Quantum Number of AO,corresponding to electron 1
     841             : !>         -) ea,eb:  Exponents of AO,corresponding to electron 1
     842             : !>         -) nb,nc:  Principle Quantum Number of AO,corresponding to electron 2
     843             : !>         -) ec,ed:  Exponents of AO,corresponding to electron 2
     844             : !> \author Teodoro Laino [tlaino] - University of Zurich
     845             : ! **************************************************************************************************
     846        8856 :    FUNCTION sc_param_low(k, na, ea, nb, eb, nc, ec, nd, ed) RESULT(res)
     847             :       INTEGER, INTENT(in)                                :: k, na
     848             :       REAL(KIND=dp), INTENT(in)                          :: ea
     849             :       INTEGER, INTENT(in)                                :: nb
     850             :       REAL(KIND=dp), INTENT(in)                          :: eb
     851             :       INTEGER, INTENT(in)                                :: nc
     852             :       REAL(KIND=dp), INTENT(in)                          :: ec
     853             :       INTEGER, INTENT(in)                                :: nd
     854             :       REAL(KIND=dp), INTENT(in)                          :: ed
     855             :       REAL(KIND=dp)                                      :: res
     856             : 
     857             :       INTEGER                                            :: i, m, m1, m2, n, nab, ncd
     858             :       REAL(KIND=dp)                                      :: a2, aab, acd, ae, aea, aeb, aec, aed, c, &
     859             :                                                             e, eab, ecd, ff, s0, s1, s2, s3, tmp
     860             : 
     861        8856 :       CPASSERT(ea > 0.0_dp)
     862        8856 :       CPASSERT(eb > 0.0_dp)
     863        8856 :       CPASSERT(ec > 0.0_dp)
     864        8856 :       CPASSERT(ed > 0.0_dp)
     865        8856 :       aea = LOG(ea)
     866        8856 :       aeb = LOG(eb)
     867        8856 :       aec = LOG(ec)
     868        8856 :       aed = LOG(ed)
     869        8856 :       nab = na + nb
     870        8856 :       ncd = nc + nd
     871        8856 :       ecd = ec + ed
     872        8856 :       eab = ea + eb
     873        8856 :       e = ecd + eab
     874        8856 :       n = nab + ncd
     875        8856 :       ae = LOG(e)
     876        8856 :       a2 = LOG(2.0_dp)
     877        8856 :       acd = LOG(ecd)
     878        8856 :       aab = LOG(eab)
     879        8856 :       ff = fac(n - 1)/SQRT(fac(2*na)*fac(2*nb)*fac(2*nc)*fac(2*nd))
     880        8856 :       tmp = na*aea + nb*aeb + nc*aec + nd*aed + 0.5_dp*(aea + aeb + aec + aed) + a2*(n + 2) - ae*n
     881        8856 :       c = evolt*ff*EXP(tmp)
     882        8856 :       s0 = 1.0_dp/e
     883        8856 :       s1 = 0.0_dp
     884        8856 :       s2 = 0.0_dp
     885        8856 :       m = ncd - k
     886       57302 :       DO i = 1, m
     887       48446 :          s0 = s0*e/ecd
     888       57302 :          s1 = s1 + s0*(binomial(ncd - k - 1, i - 1) - binomial(ncd + k, i - 1))/binomial(n - 1, i - 1)
     889             :       END DO
     890        8856 :       m1 = m
     891        8856 :       m2 = ncd + k
     892       45756 :       DO i = m1, m2
     893       36900 :          s0 = s0*e/ecd
     894       45756 :          s2 = s2 + s0*binomial(m2, i)/binomial(n - 1, i)
     895             :       END DO
     896        8856 :       s3 = EXP(ae*n - acd*(m2 + 1) - aab*(nab - k))/binomial(n - 1, m2)
     897        8856 :       res = c*(s1 - s2 + s3)
     898        8856 :    END FUNCTION sc_param_low
     899             : 
     900             : ! **************************************************************************************************
     901             : !> \brief  Corrects the EISOL fo the one-center terms coming from those atoms
     902             : !>         that have partially filled "d" shells
     903             : !> \param sep ...
     904             : !> \param r016 ...
     905             : !> \param r066 ...
     906             : !> \param r244 ...
     907             : !> \param r266 ...
     908             : !> \param r466 ...
     909             : !> \date   03.2008 [tlaino]
     910             : !> \par    Notation
     911             : !>         r016:  <SS|DD>
     912             : !>         r066:  <DD|DD> "0" term
     913             : !>         r244:  <SD|SD>
     914             : !>         r266:  <DD|DD> "2" term
     915             : !>         r466:  <DD|DD> "4" term
     916             : !>
     917             : !> \author Teodoro Laino [tlaino] - University of Zurich
     918             : ! **************************************************************************************************
     919         738 :    SUBROUTINE eisol_corr(sep, r016, r066, r244, r266, r466)
     920             :       TYPE(semi_empirical_type), POINTER                 :: sep
     921             :       REAL(KIND=dp), INTENT(in)                          :: r016, r066, r244, r266, r466
     922             : 
     923             :       sep%eisol = sep%eisol + ir016(sep%z)*r016 + &
     924             :                   ir066(sep%z)*r066 - &
     925             :                   ir244(sep%z)*r244/5.0_dp - &
     926             :                   ir266(sep%z)*r266/49.0_dp - &
     927         738 :                   ir466(sep%z)*r466/49.0_dp
     928         738 :    END SUBROUTINE eisol_corr
     929             : 
     930             : ! **************************************************************************************************
     931             : !> \brief  Computes the Klopman-Ohno additive terms for 2-center 2-electron
     932             : !>         integrals requiring that the corresponding 1-center 2-electron integral
     933             : !>         is reproduced from the 2-center one for r->0
     934             : !> \param l ...
     935             : !> \param d ...
     936             : !> \param fg ...
     937             : !> \return ...
     938             : !> \date   03.2008 [tlaino]
     939             : !> \author Teodoro Laino [tlaino] - University of Zurich
     940             : ! **************************************************************************************************
     941        5166 :    FUNCTION ko_ij(l, d, fg) RESULT(res)
     942             :       INTEGER, INTENT(in)                                :: l
     943             :       REAL(KIND=dp), INTENT(in)                          :: d, fg
     944             :       REAL(KIND=dp)                                      :: res
     945             : 
     946             :       INTEGER, PARAMETER                                 :: niter = 100
     947             :       REAL(KIND=dp), PARAMETER                           :: epsil = 1.0E-08_dp, g1 = 0.382_dp, &
     948             :                                                             g2 = 0.618_dp
     949             : 
     950             :       INTEGER                                            :: i
     951             :       REAL(KIND=dp)                                      :: a1, a2, delta, dsq, ev4, ev8, f1, f2, &
     952             :                                                             y1, y2
     953             : 
     954        5166 :       CPASSERT(fg /= 0.0_dp)
     955             :       ! Term for SS
     956        5166 :       IF (l == 0) THEN
     957        1476 :          res = 0.5_dp*evolt/fg
     958        1476 :          RETURN
     959             :       END IF
     960             :       ! Term for Higher angular momentum
     961        3690 :       dsq = d*d
     962        3690 :       ev4 = evolt*0.25_dp
     963        3690 :       ev8 = evolt/8.0_dp
     964        3690 :       a1 = 0.1_dp
     965        3690 :       a2 = 5.0_dp
     966      158670 :       DO i = 1, niter
     967      158670 :          delta = a2 - a1
     968      158670 :          IF (delta < epsil) EXIT
     969      154980 :          y1 = a1 + delta*g1
     970      154980 :          y2 = a1 + delta*g2
     971      154980 :          IF (l == 1) THEN
     972       61992 :             f1 = (ev4*(1/y1 - 1/SQRT(y1**2 + dsq)) - fg)**2
     973       61992 :             f2 = (ev4*(1/y2 - 1/SQRT(y2**2 + dsq)) - fg)**2
     974       92988 :          ELSE IF (l == 2) THEN
     975       92988 :             f1 = (ev8*(1.0_dp/y1 - 2.0_dp/SQRT(y1**2 + dsq*0.5_dp) + 1.0_dp/SQRT(y1**2 + dsq)) - fg)**2
     976       92988 :             f2 = (ev8*(1/y2 - 2.0_dp/SQRT(y2**2 + dsq*0.5_dp) + 1.0_dp/SQRT(y2**2 + dsq)) - fg)**2
     977             :          END IF
     978      158670 :          IF (f1 < f2) THEN
     979             :             a2 = y2
     980             :          ELSE
     981       74844 :             a1 = y1
     982             :          END IF
     983             :       END DO
     984             :       ! Convergence reached.. define additive terms
     985        3690 :       IF (f1 >= f2) THEN
     986             :          res = a2
     987             :       ELSE
     988        1832 :          res = a1
     989             :       END IF
     990             :    END FUNCTION ko_ij
     991             : 
     992             : ! **************************************************************************************************
     993             : !> \brief  Fills the 1 center 2 electron integrals for the construction of the
     994             : !>         one-electron fock matrix
     995             : !> \param sep ...
     996             : !> \date   04.2008 [tlaino]
     997             : !> \author Teodoro Laino [tlaino] - University of Zurich
     998             : ! **************************************************************************************************
     999        3964 :    SUBROUTINE setup_1c_2el_int(sep)
    1000             :       TYPE(semi_empirical_type), POINTER                 :: sep
    1001             : 
    1002             :       INTEGER                                            :: i, ij, ij0, ind, ip, ipx, ipy, ipz, &
    1003             :                                                             isize, kl, natorb
    1004             :       LOGICAL                                            :: defined
    1005             :       REAL(KIND=dp)                                      :: gp2, gpp, gsp, gss, hsp
    1006             : 
    1007             :       CALL get_se_param(sep, defined=defined, natorb=natorb, &
    1008        3964 :                         gss=gss, gsp=gsp, gpp=gpp, gp2=gp2, hsp=hsp)
    1009        3964 :       CPASSERT(defined)
    1010             : 
    1011        3964 :       isize = natorb*(natorb + 1)/2
    1012       12408 :       ALLOCATE (sep%w(isize, isize))
    1013             :       ! Initialize array
    1014     1646300 :       sep%w = 0.0_dp
    1015             :       ! Fill the array
    1016        3964 :       IF (natorb > 0) THEN
    1017        2240 :          ip = 1
    1018        2240 :          sep%w(ip, ip) = gss
    1019        2240 :          IF (natorb > 2) THEN
    1020        1772 :             ipx = ip + 2
    1021        1772 :             ipy = ip + 5
    1022        1772 :             ipz = ip + 9
    1023        1772 :             sep%w(ipx, ip) = gsp
    1024        1772 :             sep%w(ipy, ip) = gsp
    1025        1772 :             sep%w(ipz, ip) = gsp
    1026        1772 :             sep%w(ip, ipx) = gsp
    1027        1772 :             sep%w(ip, ipy) = gsp
    1028        1772 :             sep%w(ip, ipz) = gsp
    1029        1772 :             sep%w(ipx, ipx) = gpp
    1030        1772 :             sep%w(ipy, ipy) = gpp
    1031        1772 :             sep%w(ipz, ipz) = gpp
    1032        1772 :             sep%w(ipy, ipx) = gp2
    1033        1772 :             sep%w(ipz, ipx) = gp2
    1034        1772 :             sep%w(ipz, ipy) = gp2
    1035        1772 :             sep%w(ipx, ipy) = gp2
    1036        1772 :             sep%w(ipx, ipz) = gp2
    1037        1772 :             sep%w(ipy, ipz) = gp2
    1038        1772 :             sep%w(ip + 1, ip + 1) = hsp
    1039        1772 :             sep%w(ip + 3, ip + 3) = hsp
    1040        1772 :             sep%w(ip + 6, ip + 6) = hsp
    1041        1772 :             sep%w(ip + 4, ip + 4) = 0.5_dp*(gpp - gp2)
    1042        1772 :             sep%w(ip + 7, ip + 7) = 0.5_dp*(gpp - gp2)
    1043        1772 :             sep%w(ip + 8, ip + 8) = 0.5_dp*(gpp - gp2)
    1044        1772 :             IF (sep%dorb) THEN
    1045             :                ij0 = ip - 1
    1046      180072 :                DO i = 1, 243
    1047      179334 :                   ij = int_ij(i)
    1048      179334 :                   kl = int_kl(i)
    1049      179334 :                   ind = int_onec2el(i)
    1050      180072 :                   sep%w(ij + ij0, kl + ij0) = sep%onec2el(ind)/evolt
    1051             :                END DO
    1052             :             END IF
    1053             :          END IF
    1054             :       END IF
    1055        3964 :    END SUBROUTINE setup_1c_2el_int
    1056             : 
    1057             : END MODULE semi_empirical_par_utils
    1058             : 

Generated by: LCOV version 1.15